!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||


 module forcing_sfwf 7,14

!BOP
! !MODULE: forcing_sfwf
! !DESCRIPTION:
!  Contains routines and variables used for determining the
!  surface fresh water flux.
!
! !REVISION HISTORY:
!  SVN:$Id: forcing_sfwf.F90 12674 2008-10-31 22:21:32Z njn01 $
!
! !USES:

   use kinds_mod
   use blocks
   use distribution
   use domain
   use constants
   use io
   use grid
   use global_reductions
   use forcing_tools
   use forcing_shf
   use ice
   use time_management
   use prognostic
   use exit_mod

   implicit none
   private
   save

! !PUBLIC MEMBER FUNCTIONS:

   public :: init_sfwf,     &
             set_sfwf

! !PUBLIC DATA MEMBERS:

   real (r8), public, allocatable, dimension(:,:,:,:) :: &
      SFWF_COMP

   real (r8), public, allocatable, dimension(:,:,:,:,:) :: &
      TFW_COMP

   real (r8), public :: &! public for use in restart
      sfwf_interp_last   ! time when last interpolation was done

   !*** water balance factors for bulk-NCEP forcing

   real (r8), public :: &! public for use in restart
      sum_precip,       &! global precip for water balance
      precip_fact = c1, &! factor for adjusting precip for water balance
      ssh_initial        ! initial ssh

   real (r8), dimension(km), public :: &
      sal_initial

   logical (log_kind), public :: &
      lfw_as_salt_flx          ! treat fw flux as virtual salt flux
                               ! even with var.thickness sfc layer
   logical (log_kind), public :: &
      lsend_precip_fact        ! if T,send precip_fact to cpl for use in fw balance
                               ! (partially-coupled option)


!EOP
!BOC
!-----------------------------------------------------------------------
!
!  module variables
!
!-----------------------------------------------------------------------

   real (r8), allocatable, dimension(:,:,:,:,:) :: &
      SFWF_DATA    ! forcing data used to get SFWF

   real (r8), dimension(12) :: &
      sfwf_data_time  ! time (hours) corresponding to surface fresh
                      ! water fluxes

   real (r8), dimension(20) :: &
      sfwf_data_renorm         ! factors for converting to model units

   real (r8) ::             &
      sfwf_data_inc,        &! time increment between values of forcing data
      sfwf_data_next,       &! time to be used for next value of forcing data
      sfwf_data_update,     &! time new forcing data needs to be added to interpolation set
      sfwf_interp_inc,      &! time increment between interpolation
      sfwf_interp_next,     &! time when next interpolation will be done
      sfwf_restore_tau,     &! restoring time scale
      sfwf_restore_rtau,    &! reciprocal of restoring time scale
      sfwf_weak_restore,    &!
      sfwf_strong_restore,  &!
      sfwf_strong_restore_ms !

   integer (int_kind) :: &
      sfwf_interp_order,      &! order of temporal interpolation
      sfwf_data_time_min_loc, &! time index for first SFWF_DATA point
      sfwf_data_num_fields      

   integer (int_kind), public :: &
      sfwf_num_comps           

   character (char_len), dimension(:), allocatable :: &
      sfwf_data_names          ! short names for input data fields

   integer (int_kind), dimension(:), allocatable :: &
      sfwf_bndy_loc,          &! location and field types for ghost
      sfwf_bndy_type           !    cell update routines

   !*** integer addresses for various forcing data fields

   integer (int_kind) :: & ! restoring and partially-coupled options
      sfwf_data_sss

   integer (int_kind), public :: &! bulk-NCEP and partially-coupled (some) options
      sfwf_data_precip,          &
      sfwf_comp_precip,          &
      sfwf_comp_evap,            &
      sfwf_comp_wrest,           &
      sfwf_comp_srest

   real (r8) ::                &
      ann_avg_precip,          &!
      !sum_fw,                 &!
      !ann_avg_fw,             &!
      ssh_final

   real (r8), dimension (km) :: &
      sal_final

   logical (log_kind) ::   &
      ladjust_precip

    integer (int_kind),public ::  &! used with the partially-coupled option
       sfwf_comp_cpl,             &
       sfwf_data_flxio,           &
       sfwf_comp_flxio,           &
       tfw_num_comps,             &
       tfw_comp_cpl,              & 
       tfw_comp_flxio       
 
   real (r8), parameter :: &
      precip_mean = 3.4e-5_r8

   character (char_len) :: &
      sfwf_filename,       &! name of file conainting forcing data
      sfwf_file_fmt,       &! format (bin or netcdf) of forcing file
      sfwf_interp_freq,    &! keyword for period of temporal interpolation
      sfwf_interp_type,    &!
      sfwf_data_label

   character (char_len),public :: &
      sfwf_data_type,      &! keyword for period of forcing data
      sfwf_formulation     

   logical (log_kind), public ::   & 
      lms_balance           ! control balancing of P,E,M,R,S in marginal seas
                            ! .T. only with sfc_layer_oldfree option

!EOC
!***********************************************************************

 contains

!***********************************************************************
!BOP
! !IROUTINE: init_sfwf
! !INTERFACE:


 subroutine init_sfwf(STF) 1,109

! !DESCRIPTION:
!  Initializes surface fresh water flux forcing by either calculating
!  or reading in the surface fresh water flux.  Also does initial
!  book-keeping concerning when new data is needed for the temporal
!  interpolation and when the forcing will need to be updated.
!
! !REVISION HISTORY:
!  same as module

! !INPUT/OUTPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block,nt,max_blocks_clinic), &
      intent(inout) :: &
      STF    !  surface tracer fluxes at current timestep

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer(int_kind) :: &
      k, n,             &! dummy loop indices
      iblock,           &! block loop index
      nml_error          ! namelist error flag

   character (char_len) :: &
      forcing_filename     ! full filename for forcing input

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
      WORK               ! temporary work space

   real (r8), dimension(:,:,:,:,:), allocatable :: &
      TEMP_DATA          ! temporary array for reading monthly data

   type (block) :: &
      this_block   ! block info for local block

   type (datafile) ::  &
      forcing_file     ! data file structure for input forcing file

   type (io_field_desc) :: &
      io_sss,         &! io field descriptor for input sss field
      io_precip,      &! io field descriptor for input precip field
      io_flxio         ! io field descriptor for input io_flxio field

   type (io_dim) :: &
      i_dim, j_dim, &! dimension descriptors for horiz dimensions
      month_dim      ! dimension descriptor  for monthly data

   namelist /forcing_sfwf_nml/ sfwf_data_type,      sfwf_data_inc,    &
                               sfwf_interp_type,    sfwf_interp_freq, &
                               sfwf_interp_inc,     sfwf_restore_tau, &
                               sfwf_filename,       sfwf_file_fmt,    &
                               sfwf_data_renorm,    sfwf_formulation, &
                               ladjust_precip,      sfwf_weak_restore,&
                               sfwf_strong_restore, lfw_as_salt_flx,  &
                               sfwf_strong_restore_ms,                &
                               lsend_precip_fact,   lms_balance

!-----------------------------------------------------------------------
!
!  read surface fresh water flux namelist input after setting
!  default values.
!
!-----------------------------------------------------------------------

   sfwf_formulation       = 'restoring'
   sfwf_data_type         = 'analytic'
   sfwf_data_inc          = 1.e20_r8
   sfwf_interp_type       = 'nearest'
   sfwf_interp_freq       = 'never'
   sfwf_interp_inc        = 1.e20_r8
   sfwf_restore_tau       = 1.e20_r8
   sfwf_filename          = 'unknown-sfwf'
   sfwf_file_fmt          = 'bin'
   sfwf_data_renorm       = c1
   !sfwf_data_renorm      = 1.e-3_r8  !  convert from psu to msu
   ladjust_precip         = .false.
   lms_balance            = .false.
   sfwf_weak_restore      = 0.092_r8
   sfwf_strong_restore_ms = 0.6648_r8
   sfwf_strong_restore    = 0.6648_r8
   lfw_as_salt_flx        = .false.
   lsend_precip_fact      = .false.

   if (my_task == master_task) then
      open (nml_in, file=nml_filename, status='old',iostat=nml_error)
      if (nml_error /= 0) then
         nml_error = -1
      else
         nml_error =  1
      endif
      do while (nml_error > 0)
         read(nml_in, nml=forcing_sfwf_nml,iostat=nml_error)
      end do
      if (nml_error == 0) close(nml_in)
   endif


   call broadcast_scalar(nml_error, master_task)
   if (nml_error /= 0) then
     call exit_POP(sigAbort,'ERROR reading forcing_sfwf_nml')
   endif

   call broadcast_scalar(sfwf_data_type,         master_task)
   call broadcast_scalar(sfwf_data_inc,          master_task)
   call broadcast_scalar(sfwf_interp_type,       master_task)
   call broadcast_scalar(sfwf_interp_freq,       master_task)
   call broadcast_scalar(sfwf_interp_inc,        master_task)
   call broadcast_scalar(sfwf_restore_tau,       master_task)
   call broadcast_scalar(sfwf_filename,          master_task)
   call broadcast_scalar(sfwf_file_fmt,          master_task)
   call broadcast_scalar(sfwf_formulation,       master_task)
   call broadcast_array (sfwf_data_renorm,       master_task)
   call broadcast_scalar(ladjust_precip,         master_task)
   call broadcast_scalar(sfwf_weak_restore,      master_task)
   call broadcast_scalar(sfwf_strong_restore,    master_task)
   call broadcast_scalar(sfwf_strong_restore_ms, master_task)
   call broadcast_scalar(lfw_as_salt_flx,        master_task)
   call broadcast_scalar(lsend_precip_fact,      master_task)
   call broadcast_scalar(lms_balance,            master_task)

!-----------------------------------------------------------------------
!
!  convert data_type to 'monthly-calendar' if input is 'monthly'
!
!-----------------------------------------------------------------------

   if (sfwf_data_type == 'monthly') sfwf_data_type = 'monthly-calendar'

!-----------------------------------------------------------------------
!
!  set values based on sfwf_formulation
!
!-----------------------------------------------------------------------

   select case (sfwf_formulation)
   case ('restoring')
      allocate(sfwf_data_names(1), &
               sfwf_bndy_loc  (1), &
               sfwf_bndy_type (1))
      sfwf_data_num_fields = 1
      sfwf_data_sss = 1
      sfwf_data_names(sfwf_data_sss) = 'SSS'
      sfwf_bndy_loc  (sfwf_data_sss) = field_loc_center
      sfwf_bndy_type (sfwf_data_sss) = field_type_scalar

   case ('bulk-NCEP')
      sfwf_data_num_fields = 2
      sfwf_data_sss = 1
      sfwf_data_precip = 2

      allocate(sfwf_data_names(sfwf_data_num_fields), &
               sfwf_bndy_loc  (sfwf_data_num_fields), &
               sfwf_bndy_type (sfwf_data_num_fields))
      sfwf_data_names(sfwf_data_sss) = 'SSS'
      sfwf_bndy_loc  (sfwf_data_sss) = field_loc_center
      sfwf_bndy_type (sfwf_data_sss) = field_type_scalar
      sfwf_data_names(sfwf_data_precip) = 'PRECIPITATION'
      sfwf_bndy_loc  (sfwf_data_precip) = field_loc_center
      sfwf_bndy_type (sfwf_data_precip) = field_type_scalar

      sfwf_num_comps = 4
      sfwf_comp_precip = 1
      sfwf_comp_evap = 2
      sfwf_comp_wrest = 3
      sfwf_comp_srest = 4


   case ('partially-coupled')
      sfwf_data_num_fields = 2 
      sfwf_data_sss   = 1
      sfwf_data_flxio = 2

      allocate(sfwf_data_names(sfwf_data_num_fields), &
               sfwf_bndy_loc  (sfwf_data_num_fields), &
               sfwf_bndy_type (sfwf_data_num_fields))
      sfwf_data_names(sfwf_data_sss) = 'SSS'
      sfwf_bndy_loc  (sfwf_data_sss) = field_loc_center
      sfwf_bndy_type (sfwf_data_sss) = field_type_scalar
      sfwf_data_names(sfwf_data_flxio) = 'FLXIO'
      sfwf_bndy_loc  (sfwf_data_flxio) = field_loc_center
      sfwf_bndy_type (sfwf_data_flxio) = field_type_scalar

 
      sfwf_num_comps  = 4 
      sfwf_comp_wrest = 1
      sfwf_comp_srest = 2
      sfwf_comp_cpl   = 3
      sfwf_comp_flxio = 4

      tfw_num_comps  = 2
      tfw_comp_cpl   = 1
      tfw_comp_flxio = 2

   case default
      call exit_POP(sigAbort, &
                    'init_sfwf: Unknown value for sfwf_formulation')

   end select
   
   if ( sfwf_formulation == 'bulk-NCEP' .or.  &
        sfwf_formulation == 'partially-coupled' ) then
      !*** calculate initial salinity profile for ocean points that are
      !*** not marginal seas.

      !*** very first step of run
      if (ladjust_precip .and. nsteps_total == 0) then

         sum_precip = c0
         ssh_initial = c0
         !sum_fw = c0

         do k = 1,km
            !$OMP PARALLEL DO PRIVATE(iblock, this_block)
            do iblock=1,nblocks_clinic
               this_block = get_block(blocks_clinic(iblock),iblock)

               if (partial_bottom_cells) then
                  WORK(:,:,iblock) = &
                         merge(TRACER(:,:,k,2,curtime,iblock)*      &
                               TAREA(:,:,iblock)*DZT(:,:,k,iblock), &
                               c0, k <= KMT(:,:,iblock) .and.       &
                                   MASK_SR(:,:,iblock) > 0)
               else
                  WORK(:,:,iblock) = &
                         merge(TRACER(:,:,k,2,curtime,iblock)*  &
                               TAREA(:,:,iblock)*dz(k),         &
                               c0, k <= KMT(:,:,iblock) .and.   &
                                   MASK_SR(:,:,iblock) > 0)
               endif
            end do
            !$OMP END PARALLEL DO

            sal_initial(k) = global_sum(WORK,distrb_clinic,field_loc_center)/ &
                             (volume_t_k(k) - volume_t_marg_k(k)) ! in m^3
         enddo
      endif
   endif

!-----------------------------------------------------------------------
!
!  calculate inverse of restoring time scale and convert to seconds.
!
!-----------------------------------------------------------------------

   sfwf_restore_rtau = c1/(seconds_in_day*sfwf_restore_tau)

!-----------------------------------------------------------------------
!
!  convert interp_type to corresponding integer value.
!
!-----------------------------------------------------------------------

   select case (sfwf_interp_type)
   case ('nearest')
      sfwf_interp_order = 1

   case ('linear')
      sfwf_interp_order = 2

   case ('4point')
      sfwf_interp_order = 4

   case default
      call exit_POP(sigAbort, &
                    'init_sfwf: Unknown value for sfwf_interp_type')

   end select

!-----------------------------------------------------------------------
!
!  set values of the surface fresh water flux arrays (SFWF or
!  SFWF_DATA) depending on type of the surface fresh water flux
!  data.
!
!-----------------------------------------------------------------------

   select case (sfwf_data_type)

!-----------------------------------------------------------------------
!
!  no surface fresh water flux, therefore no interpolation in time
!  is needed (sfwf_interp_freq = 'none'), nor are there any new
!  values to be used (sfwf_data_next = sfwf_data_update = never).
!
!-----------------------------------------------------------------------

   case ('none')

     STF(:,:,2,:) = c0
     sfwf_data_next = never
     sfwf_data_update = never
     sfwf_interp_freq = 'never'

!-----------------------------------------------------------------------
!
!  simple analytic surface salinity that is constant in time,
!  therefore no new values will be needed.
!
!-----------------------------------------------------------------------

   case ('analytic')

      allocate(SFWF_DATA(nx_block,ny_block,max_blocks_clinic, &
                         sfwf_data_num_fields,1))
      SFWF_DATA = c0

      select case (sfwf_formulation)
      case ('restoring')
         SFWF_DATA(:,:,:,sfwf_data_sss,1) = 0.035_r8
      end select

      sfwf_data_next = never
      sfwf_data_update = never
      sfwf_interp_freq = 'never'

!-----------------------------------------------------------------------
!
!  annual mean climatological surface salinity (read in from a file)
!  that is constant in time, therefore no new values will be needed.
!
!-----------------------------------------------------------------------

   case ('annual')

      allocate(SFWF_DATA(nx_block,ny_block,max_blocks_clinic, &
                         sfwf_data_num_fields,1))
      SFWF_DATA = c0

      forcing_file = construct_file(sfwf_file_fmt,                 &
                                    full_name=trim(sfwf_filename), &
                                    record_length = rec_type_dbl,  &
                                    recl_words=nx_global*ny_global)

      call data_set(forcing_file,'open_read')

      i_dim = construct_io_dim('i',nx_global)
      j_dim = construct_io_dim('j',nx_global)

      select case (sfwf_formulation)
      case ('restoring')

         io_sss = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_sss)),       &
                       dim1=i_dim, dim2=j_dim,                               &
                       field_loc = sfwf_bndy_loc(sfwf_data_sss),   &
                       field_type = sfwf_bndy_type(sfwf_data_sss), &
                       d2d_array=SFWF_DATA(:,:,:,sfwf_data_sss,1))
         call data_set(forcing_file,'define',io_sss)
         call data_set(forcing_file,'read'  ,io_sss)
         call destroy_io_field(io_sss)

      case ('bulk-NCEP')

         io_sss = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_sss)),       &
                       dim1=i_dim, dim2=j_dim,                               &
                       field_loc = sfwf_bndy_loc(sfwf_data_sss),   &
                       field_type = sfwf_bndy_type(sfwf_data_sss), &
                       d2d_array=SFWF_DATA(:,:,:,sfwf_data_sss,1))
         io_precip = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_precip)),       &
                       dim1=i_dim, dim2=j_dim,                                  &
                       field_loc = sfwf_bndy_loc(sfwf_data_precip),   &
                       field_type = sfwf_bndy_type(sfwf_data_precip), &
                       d2d_array=SFWF_DATA(:,:,:,sfwf_data_precip,1))
         call data_set(forcing_file,'define',io_sss)
         call data_set(forcing_file,'define',io_precip)
         call data_set(forcing_file,'read'  ,io_sss)
         call data_set(forcing_file,'read'  ,io_precip)
         call destroy_io_field(io_sss)
         call destroy_io_field(io_precip)

         allocate( SFWF_COMP(nx_block,ny_block,max_blocks_clinic, &
                             sfwf_num_comps))
         SFWF_COMP = c0  ! initialize

      case ('partially-coupled')

         io_sss = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_sss)),       &
                       dim1=i_dim, dim2=j_dim,                               &
                       field_loc = sfwf_bndy_loc(sfwf_data_sss),   &
                       field_type = sfwf_bndy_type(sfwf_data_sss), &
                       d2d_array=SFWF_DATA(:,:,:,sfwf_data_sss,1))
         io_flxio = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_flxio)),       &
                       dim1=i_dim, dim2=j_dim,                                 &
                       field_loc = sfwf_bndy_loc(sfwf_data_flxio),   &
                       field_type = sfwf_bndy_type(sfwf_data_flxio), &
                       d2d_array=SFWF_DATA(:,:,:,sfwf_data_flxio,1))
         call data_set(forcing_file,'define',io_sss)
         call data_set(forcing_file,'define',io_flxio)
         call data_set(forcing_file,'read'  ,io_sss)
         call data_set(forcing_file,'read'  ,io_flxio)
         call destroy_io_field(io_sss)
         call destroy_io_field(io_flxio)

         allocate( SFWF_COMP(nx_block,ny_block,max_blocks_clinic, &
                             sfwf_num_comps))
         SFWF_COMP = c0  ! initialize

      end select

      call data_set(forcing_file,'close')
      call destroy_file(forcing_file)

      !*** renormalize values if necessary to compensate for
      !*** different units.

      do n = 1,sfwf_data_num_fields
         if (sfwf_data_renorm(n) /= c1) SFWF_DATA(:,:,:,n,:) = &
                    sfwf_data_renorm(n)*SFWF_DATA(:,:,:,n,:)
      enddo

      sfwf_data_next = never
      sfwf_data_update = never
      sfwf_interp_freq = 'never'

      if (my_task == master_task) then
         write(stdout,blank_fmt)
         write(stdout,'(a25,a)') ' SFWF Annual file read: ', &
                                 trim(sfwf_filename)
      endif

!-----------------------------------------------------------------------
!
!  monthly mean climatological surface fresh water flux. all
!  12 months are read in from a file. interpolation order
!  (sfwf_interp_order) may be specified with namelist input.
!
!-----------------------------------------------------------------------

   case ('monthly-equal','monthly-calendar')

      allocate(SFWF_DATA(nx_block,ny_block,max_blocks_clinic,    &
                                     sfwf_data_num_fields,0:12), &
               TEMP_DATA(nx_block,ny_block,12,max_blocks_clinic, &
                                           sfwf_data_num_fields) )
      SFWF_DATA = c0

      call find_forcing_times(sfwf_data_time,   sfwf_data_inc,  &
                              sfwf_interp_type, sfwf_data_next, &
                              sfwf_data_time_min_loc,           &
                              sfwf_data_update, sfwf_data_type)

      forcing_file = construct_file(sfwf_file_fmt,                 &
                                    full_name=trim(sfwf_filename), &
                                    record_length = rec_type_dbl,  &
                                    recl_words=nx_global*ny_global)
      call data_set(forcing_file,'open_read')

      i_dim = construct_io_dim('i',nx_global)
      j_dim = construct_io_dim('j',nx_global)
      month_dim = construct_io_dim('month',12)

      select case (sfwf_formulation)
      case ('restoring')
         io_sss = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_sss)),       &
                       dim1=i_dim, dim2=j_dim, dim3=month_dim,               &
                       field_loc = sfwf_bndy_loc(sfwf_data_sss),   &
                       field_type = sfwf_bndy_type(sfwf_data_sss), &
                       d3d_array=TEMP_DATA(:,:,:,:,sfwf_data_sss))
         call data_set(forcing_file,'define',io_sss)
         call data_set(forcing_file,'read'  ,io_sss)
         call destroy_io_field(io_sss)

      case ('bulk-NCEP')
         io_sss = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_sss)),       &
                       dim1=i_dim, dim2=j_dim, dim3=month_dim,               &
                       field_loc = sfwf_bndy_loc(sfwf_data_sss),   &
                       field_type = sfwf_bndy_type(sfwf_data_sss), &
                       d3d_array=TEMP_DATA(:,:,:,:,sfwf_data_sss))
         io_precip = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_precip)),       &
                       dim1=i_dim, dim2=j_dim, dim3=month_dim,                  &
                       field_loc = sfwf_bndy_loc(sfwf_data_precip),   &
                       field_type = sfwf_bndy_type(sfwf_data_precip), &
                       d3d_array=TEMP_DATA(:,:,:,:,sfwf_data_precip))
         call data_set(forcing_file,'define',io_sss)
         call data_set(forcing_file,'define',io_precip)
         call data_set(forcing_file,'read'  ,io_sss)
         call data_set(forcing_file,'read'  ,io_precip)
         call destroy_io_field(io_sss)
         call destroy_io_field(io_precip)

         allocate(SFWF_COMP(nx_block,ny_block,max_blocks_clinic, &
                                                 sfwf_num_comps))
         SFWF_COMP = c0  ! initialize

      case ('partially-coupled')
         io_sss = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_sss)),       &
                       dim1=i_dim, dim2=j_dim, dim3=month_dim,               &
                       field_loc = sfwf_bndy_loc(sfwf_data_sss),   &
                       field_type = sfwf_bndy_type(sfwf_data_sss), &
                       d3d_array=TEMP_DATA(:,:,:,:,sfwf_data_sss))
         io_flxio  = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_flxio )),       &
                       dim1=i_dim, dim2=j_dim, dim3=month_dim,                  &
                       field_loc = sfwf_bndy_loc(sfwf_data_flxio ),   &
                       field_type = sfwf_bndy_type(sfwf_data_flxio ), &
                       d3d_array=TEMP_DATA(:,:,:,:,sfwf_data_flxio ))
         call data_set(forcing_file,'define',io_sss)
         call data_set(forcing_file,'define',io_flxio )
         call data_set(forcing_file,'read'  ,io_sss)
         call data_set(forcing_file,'read'  ,io_flxio )
         call destroy_io_field(io_sss)
         call destroy_io_field(io_flxio )

         allocate(SFWF_COMP(nx_block,ny_block,max_blocks_clinic, &
                                                 sfwf_num_comps))
         SFWF_COMP = c0  ! initialize

      end select

      call data_set(forcing_file,'close')
      call destroy_file(forcing_file)

      !*** re-order data and renormalize values if necessary to
      !*** compensate for different units.

      !$OMP PARALLEL DO PRIVATE(iblock, k, n)
      do iblock=1,nblocks_clinic
         do k=1,sfwf_data_num_fields
            if (sfwf_data_renorm(k) /= c1) then
               do n=1,12
                  SFWF_DATA(:,:,iblock,k,n) = &
                  TEMP_DATA(:,:,n,iblock,k)*sfwf_data_renorm(k)
               end do
            else
               do n=1,12
                  SFWF_DATA(:,:,iblock,k,n) = TEMP_DATA(:,:,n,iblock,k)
               end do
            endif
         end do
      end do
      !$OMP END PARALLEL DO

      deallocate(TEMP_DATA)

      if (my_task == master_task) then
         write(stdout,blank_fmt)
         write(stdout,'(a25,a)') ' SFWF Monthly file read: ', &
                                 trim(sfwf_filename)
      endif

!-----------------------------------------------------------------------
!
!  surface salinity specified every n-hours, where the n-hour
!  increment should be specified with namelist input
!  (sfwf_data_inc). only as many times as are necessary based on
!  the order of the temporal interpolation scheme
!  (sfwf_interp_order) reside in memory at any given time.
!
!-----------------------------------------------------------------------

   case ('n-hour')

      allocate( SFWF_DATA(nx_block,ny_block,max_blocks_clinic, &
                          sfwf_data_num_fields,0:sfwf_interp_order))
      SFWF_DATA = c0

      call find_forcing_times(sfwf_data_time,   sfwf_data_inc,  &
                              sfwf_interp_type, sfwf_data_next, &
                              sfwf_data_time_min_loc,           &
                              sfwf_data_update, sfwf_data_type)

      do n=1,sfwf_interp_order
         call get_forcing_filename(forcing_filename,  sfwf_filename, &
                                   sfwf_data_time(n), sfwf_data_inc)

         forcing_file = construct_file(sfwf_file_fmt,                 &
                                       full_name=trim(sfwf_filename), &
                                       record_length = rec_type_dbl,  &
                                       recl_words=nx_global*ny_global)
         call data_set(forcing_file,'open_read')

         i_dim = construct_io_dim('i',nx_global)
         j_dim = construct_io_dim('j',nx_global)

         select case (sfwf_formulation)
         case ('restoring')
            io_sss = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_sss)),       &
                       dim1=i_dim, dim2=j_dim,                               &
                       field_loc = sfwf_bndy_loc(sfwf_data_sss),   &
                       field_type = sfwf_bndy_type(sfwf_data_sss), &
                       d2d_array=SFWF_DATA(:,:,:,sfwf_data_sss,n))
            call data_set(forcing_file,'define',io_sss)
            call data_set(forcing_file,'read'  ,io_sss)
            call destroy_io_field(io_sss)

         case ('bulk-NCEP')
            io_sss = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_sss)),       &
                       dim1=i_dim, dim2=j_dim,                               &
                       field_loc = sfwf_bndy_loc(sfwf_data_sss),   &
                       field_type = sfwf_bndy_type(sfwf_data_sss), &
                       d2d_array=SFWF_DATA(:,:,:,sfwf_data_sss,n))
            io_precip = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_precip)),       &
                       dim1=i_dim, dim2=j_dim,                                  &
                       field_loc = sfwf_bndy_loc(sfwf_data_precip),   &
                       field_type = sfwf_bndy_type(sfwf_data_precip), &
                       d2d_array=SFWF_DATA(:,:,:,sfwf_data_precip,n))
            call data_set(forcing_file,'define',io_sss)
            call data_set(forcing_file,'define',io_precip)
            call data_set(forcing_file,'read'  ,io_sss)
            call data_set(forcing_file,'read'  ,io_precip)
            call destroy_io_field(io_sss)
            call destroy_io_field(io_precip)

         case ('partially-coupled')
            io_sss = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_sss)),       &
                       dim1=i_dim, dim2=j_dim,                               &
                       field_loc = sfwf_bndy_loc(sfwf_data_sss),   &
                       field_type = sfwf_bndy_type(sfwf_data_sss), &
                       d2d_array=SFWF_DATA(:,:,:,sfwf_data_sss,n))
            io_flxio  = construct_io_field( &
                       trim(sfwf_data_names(sfwf_data_flxio )),       &
                       dim1=i_dim, dim2=j_dim,                                  &
                       field_loc = sfwf_bndy_loc(sfwf_data_flxio ),   &
                       field_type = sfwf_bndy_type(sfwf_data_flxio ), &
                       d2d_array=SFWF_DATA(:,:,:,sfwf_data_flxio ,n))
            call data_set(forcing_file,'define',io_sss)
            call data_set(forcing_file,'define',io_flxio )
            call data_set(forcing_file,'read'  ,io_sss)
            call data_set(forcing_file,'read'  ,io_flxio )
            call destroy_io_field(io_sss)
            call destroy_io_field(io_flxio )

         end select

         call data_set(forcing_file,'close')
         call destroy_file(forcing_file)

         if (my_task == master_task) then
            write(stdout,blank_fmt)
            write(stdout,'(a24,a)') ' SFWF n-hour file read: ', &
                                    trim(forcing_filename)
         endif
      enddo

      if (sfwf_formulation == 'bulk-NCEP' .or.  &
          sfwf_formulation == 'partially-coupled') then
         allocate(SFWF_COMP(nx_block,ny_block,max_blocks_clinic, &
                                                 sfwf_num_comps))
         SFWF_COMP = c0  ! initialize
      endif

      !*** renormalize values if necessary to compensate for different
      !*** units.

      do n = 1,sfwf_data_num_fields
         if (sfwf_data_renorm(n) /= c1) SFWF_DATA(:,:,:,n,:) = &
                    sfwf_data_renorm(n)*SFWF_DATA(:,:,:,n,:)
      enddo

   case default

      call exit_POP(sigAbort, &
                    'init_sfwf: Unknown value for sfwf_data_type')

   end select

   if ( sfwf_formulation == 'partially-coupled' ) then
     allocate ( TFW_COMP(nx_block,ny_block,nt,max_blocks_clinic, &
                                              tfw_num_comps))
     TFW_COMP = c0
   endif


!-----------------------------------------------------------------------
!
!  now check interpolation period (sfwf_interp_freq) to set the
!    time for the next temporal interpolation (sfwf_interp_next).
!
!  if no interpolation is to be done, set next interpolation time
!    to a large number so the surface fresh water flux update test
!    in routine set_surface_forcing will always be false.
!
!  if interpolation is to be done every n-hours, find the first
!    interpolation time greater than the current time.
!
!  if interpolation is to be done every timestep, set next interpolation
!    time to a large negative number so the surface fresh water flux
!    update test in routine set_surface_forcing will always be true.
!
!-----------------------------------------------------------------------

   select case (sfwf_interp_freq)

   case ('never')

      sfwf_interp_next = never
      sfwf_interp_last = never
      sfwf_interp_inc  = c0

   case ('n-hour')

      call find_interp_time(sfwf_interp_inc, sfwf_interp_next)

   case ('every-timestep')

      sfwf_interp_next = always
      sfwf_interp_inc  = c0

   case default

      call exit_POP(sigAbort, &
                    'init_sfwf: Unknown value for sfwf_interp_freq')

   end select

   if(nsteps_total == 0) sfwf_interp_last = thour00

!-----------------------------------------------------------------------
!
!  echo forcing options to stdout.
!
!-----------------------------------------------------------------------

   sfwf_data_label = 'Surface Fresh Water Flux'
   call echo_forcing_options(sfwf_data_type,                     &
                             sfwf_formulation, sfwf_data_inc,    &
                             sfwf_interp_freq, sfwf_interp_type, &
                             sfwf_interp_inc,  sfwf_data_label)

   if (my_task == master_task) then
      if (lfw_as_salt_flx .and. sfc_layer_type == sfc_layer_varthick) &
         write(stdout,'(a47)') &
         '    Fresh water flux input as virtual salt flux'
   endif

!-----------------------------------------------------------------------
!EOC

 call flushm (stdout)

 end subroutine init_sfwf

!***********************************************************************
!BOP
! !IROUTINE: set_sfwf
! !INTERFACE:


 subroutine set_sfwf(STF,FW,TFW) 1,14

! !DESCRIPTION:
!  Updates the current value of the surface fresh water flux arrays
!  by interpolating fields or computing fields at the current time.
!  If new data are necessary for the interpolation, the new data are
!  read from a file.
!
! !REVISION HISTORY:
!  same as module

! !INPUT/OUTPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block,nt,max_blocks_clinic), &
      intent(inout) :: &
      STF,    &! surface tracer fluxes at current timestep
      TFW      ! tracer concentration in fresh water flux

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic), &
      intent(inout) :: &
      FW       ! fresh water flux if using varthick sfc layer

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      iblock             ! local address for current block

   type (block) :: &
      this_block   ! block info for current block

!-----------------------------------------------------------------------
!
!  check if new data is necessary for interpolation.  if yes, then
!    shuffle indices in SFWF_DATA and sfwf_data_time arrays
!    and read in new data if necessary ('n-hour' case).  note
!    that no new data is necessary for 'analytic' and 'annual' cases.
!  then perform interpolation or computation of fluxes at current time
!    using updated forcing data.
!
!-----------------------------------------------------------------------

   select case(sfwf_data_type)

   case ('analytic')

      select case (sfwf_formulation)
      case ('restoring')
         !$OMP PARALLEL DO PRIVATE(iblock, this_block)
         do iblock=1,nblocks_clinic
            this_block = get_block(blocks_clinic(iblock),iblock)

            STF(:,:,2,iblock) = &
                              (SFWF_DATA(:,:,iblock,sfwf_data_sss,1) - &
                               TRACER(:,:,1,2,curtime,iblock))*        &
                              sfwf_restore_rtau*dz(1)
         end do
         !$OMP END PARALLEL DO

      end select

   case('annual')

      select case (sfwf_formulation)
      case ('restoring')
         !$OMP PARALLEL DO PRIVATE(iblock, this_block)
         do iblock=1,nblocks_clinic
            this_block = get_block(blocks_clinic(iblock),iblock)

            STF(:,:,2,iblock) = &
                              (SFWF_DATA(:,:,iblock,sfwf_data_sss,1) - &
                               TRACER(:,:,1,2,curtime,iblock))*        &
                              sfwf_restore_rtau*dz(1)
         end do
         !$OMP END PARALLEL DO

      case ('bulk-NCEP')
         call calc_sfwf_bulk_ncep(STF,FW,TFW,1)

      case ('partially-coupled')
         call calc_sfwf_partially_coupled(1)

      end select

   case ('monthly-equal','monthly-calendar')

      sfwf_data_label = 'SFWF Monthly'
      if (thour00 >= sfwf_data_update) then
         call update_forcing_data(                   sfwf_data_time,   &
                             sfwf_data_time_min_loc, sfwf_interp_type, &
                             sfwf_data_next,         sfwf_data_update, &
                             sfwf_data_type,         sfwf_data_inc,    &
                             SFWF_DATA(:,:,:,:,1:12),sfwf_data_renorm, &
                             sfwf_data_label,        sfwf_data_names,  &
                             sfwf_bndy_loc,          sfwf_bndy_type,   &
                             sfwf_filename,          sfwf_file_fmt)
      endif

      if (thour00 >= sfwf_interp_next .or. nsteps_run == 0) then
         call interpolate_forcing(SFWF_DATA(:,:,:,:,0),               &
                                  SFWF_DATA(:,:,:,:,1:12),            &
                            sfwf_data_time,         sfwf_interp_type, &
                            sfwf_data_time_min_loc, sfwf_interp_freq, &
                            sfwf_interp_inc,        sfwf_interp_next, &
                            sfwf_interp_last,       nsteps_run)
         if (nsteps_run /= 0) sfwf_interp_next = &
                              sfwf_interp_next + sfwf_interp_inc
      endif

      select case (sfwf_formulation)
      case ('restoring')

         !$OMP PARALLEL DO PRIVATE(iblock, this_block)
         do iblock=1,nblocks_clinic
            this_block = get_block(blocks_clinic(iblock),iblock)

            STF(:,:,2,iblock) = &
                              (SFWF_DATA(:,:,iblock,sfwf_data_sss,0) - &
                               TRACER(:,:,1,2,curtime,iblock))*        &
                              sfwf_restore_rtau*dz(1)
         end do
         !$OMP END PARALLEL DO

      case ('bulk-NCEP')

         call calc_sfwf_bulk_ncep(STF,FW,TFW,12)

      case ('partially-coupled')

         call calc_sfwf_partially_coupled(12)

      end select

   case('n-hour')

      sfwf_data_label = 'SFWF n-hour'
      if (thour00 >= sfwf_data_update) then
         call update_forcing_data(                   sfwf_data_time,   &
                             sfwf_data_time_min_loc, sfwf_interp_type, &
                             sfwf_data_next,         sfwf_data_update, &
                             sfwf_data_type,         sfwf_data_inc,    &
                             SFWF_DATA(:,:,:,:,1:sfwf_interp_order),   &
                             sfwf_data_renorm,                         &
                             sfwf_data_label, sfwf_data_names,         &
                             sfwf_bndy_loc,   sfwf_bndy_type,          &
                             sfwf_filename,   sfwf_file_fmt)
      endif

      if (thour00 >= sfwf_interp_next .or. nsteps_run == 0) then
         call interpolate_forcing(SFWF_DATA(:,:,:,:,0),               &
                            SFWF_DATA(:,:,:,:,1:sfwf_interp_order),   &
                            sfwf_data_time,         sfwf_interp_type, &
                            sfwf_data_time_min_loc, sfwf_interp_freq, &
                            sfwf_interp_inc,        sfwf_interp_next, &
                            sfwf_interp_last,       nsteps_run)
         if (nsteps_run /= 0) sfwf_interp_next = &
                              sfwf_interp_next + sfwf_interp_inc
      endif

      select case (sfwf_formulation)
      case ('restoring')

         !$OMP PARALLEL DO PRIVATE(iblock, this_block)
         do iblock=1,nblocks_clinic
            this_block = get_block(blocks_clinic(iblock),iblock)

            STF(:,:,2,iblock) = &
                              (SFWF_DATA(:,:,iblock,sfwf_data_sss,0) - &
                               TRACER(:,:,1,2,curtime,iblock))*        &
                              sfwf_restore_rtau*dz(1)
         end do
         !$OMP END PARALLEL DO

      case ('bulk-NCEP')

         call calc_sfwf_bulk_ncep(STF, FW, TFW, sfwf_interp_order)

      case ('partially-coupled')

         call calc_sfwf_partially_coupled(sfwf_interp_order)

      end select

   end select

!-----------------------------------------------------------------------
!EOC

 end subroutine set_sfwf

!***********************************************************************
!BOP
! !IROUTINE: calc_sfwf_bulk_ncep
! !INTERFACE:


 subroutine calc_sfwf_bulk_ncep(STF, FW, TFW, time_dim) 3,5

! !DESCRIPTION:
!  Calculates surface freshwater flux from a combination of
!  air-sea fluxes (precipitation,  evaporation based on
!  latent heat flux computed in calc\_shf\_bulk\_ncep),
!  and restoring terms (due to restoring fields of SSS).
!
!  Notes:
!    the forcing data (on t-grid) are computed and
!    stored in SFWF\_DATA(:,:,sfwf\_comp\_*,now) where:
!       sfwf\_data\_sss    is   restoring SSS       (psu)
!       sfwf\_data\_precip is   precipitation       (m/y)
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   integer (int_kind), intent(in) ::  &
      time_dim                ! number of time points for interpolation

! !INPUT/OUTPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block,nt,max_blocks_clinic), &
      intent(inout) :: &
      STF,     &! surface tracer fluxes for all tracers
      TFW       ! tracer concentration in fresh water flux

! !OUTPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic), &
      intent(out) :: &
      FW        !  fresh water flux if using varthick sfc layer

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      now,               &! index for location of interpolated data
      k, n,              &! dummy loop indices
      iblock              ! block loop index

   real (r8) ::           &
      dttmp,              &! temporary time step variable
      fres_hor_ave,       &! area-weighted mean of weak restoring
      fres_hor_area,      &! total area of weak restoring
      area_glob_m_marg,   &! total ocean area - marginal sea area
      vol_glob_m_marg,    &! total ocean volume - marginal sea volume
      weak_mean            ! mean weak restoring

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
      WORK1, WORK2   ! temporary work space

   type(block) :: &
      this_block  ! block info for current block

!-----------------------------------------------------------------------
!
!  sfwf_weak_restore= weak(non-ice) restoring h2o flux per msu  (kg/s/m^2/msu)
!  sfwf_strong_restore= strong  (ice) ..   ..  ..   ..  ..   ..
!
!  to calculate restoring factors, use mixed layer of 50m,
!  and restoring time constant tau (days):
!
!           F (kg/s/m^2/msu)
!   tau =   6  :     2.77
!   tau =  30  :     0.55
!   tau = 182.5:     0.092
!   tau = 365  :     0.046
!   tau = 730  :     0.023
!   tau = Inf  :     0.0
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  set location where interpolated data exists.
!
!-----------------------------------------------------------------------

   if (sfwf_data_type == 'annual') then
      now = 1
   else
      now = 0
   endif

!-----------------------------------------------------------------------
!
!  compute forcing terms for each block
!
!-----------------------------------------------------------------------

   !$OMP PARALLEL DO PRIVATE(iblock, this_block)
   do iblock=1,nblocks_clinic
      this_block = get_block(blocks_clinic(iblock),iblock)

      !***  compute evaporation from latent heat computed in shf forcing

      SFWF_COMP(:,:,iblock,sfwf_comp_evap) = &
       SHF_COMP(:,:,iblock,shf_comp_qlat)/latent_heat_vapor

      !*** precipitation (kg/m^2/s)

      SFWF_COMP(:,:,iblock,sfwf_comp_precip) = &
      SFWF_DATA(:,:,iblock,sfwf_data_precip,now)*precip_fact

      ! *c1000/seconds_in_year ! convert m/y to Kg/m^2/s if needed

      !*** weak salinity restoring term
      !*** (note: OCN_WGT = 0. at land points)
      !*** will be subtracting global mean later, so compute
      !*** necessary terms for global mean

      SFWF_COMP(:,:,iblock,sfwf_comp_wrest) = &
                         -sfwf_weak_restore*OCN_WGT(:,:,iblock)*    &
                          MASK_SR(:,:,iblock)*                      &
                         (SFWF_DATA(:,:,iblock,sfwf_data_sss,now) - &
                          TRACER(:,:,1,2,curtime,iblock))

      WORK1(:,:,iblock) = TAREA(:,:,iblock)*                        &
                          SFWF_COMP(:,:,iblock,sfwf_comp_wrest)
      WORK2(:,:,iblock) = TAREA(:,:,iblock)*OCN_WGT(:,:,iblock)*    &
                          MASK_SR(:,:,iblock)

      !*** strong salinity restoring term

      where (KMT(:,:,iblock) > 0)
         SFWF_COMP(:,:,iblock,sfwf_comp_srest) = &
                     -sfwf_strong_restore*(c1 - OCN_WGT(:,:,iblock))* &
                      (SFWF_DATA(:,:,iblock,sfwf_data_sss,now) -      &
                       TRACER(:,:,1,2,curtime,iblock))
      endwhere

      where (KMT(:,:,iblock) > 0 .and. MASK_SR(:,:,iblock) == 0) 
         SFWF_COMP(:,:,iblock,sfwf_comp_srest) = &
                     -sfwf_strong_restore_ms*    &
                      (SFWF_DATA(:,:,iblock,sfwf_data_sss,now) -      &
                       TRACER(:,:,1,2,curtime,iblock))
      endwhere

   end do ! block loop
   !$OMP END PARALLEL DO

!----------------------------------------------------------------------
!
!  compute global mean of weak restoring term
!
!----------------------------------------------------------------------

   fres_hor_ave  = global_sum(WORK1, distrb_clinic, field_loc_center)
   fres_hor_area = global_sum(WORK2, distrb_clinic, field_loc_center)
   weak_mean = fres_hor_ave/fres_hor_area

!-----------------------------------------------------------------------
!
!  finish computing forcing terms for each block
!
!-----------------------------------------------------------------------

   !$OMP PARALLEL DO PRIVATE(iblock, this_block)
   do iblock=1,nblocks_clinic
      this_block = get_block(blocks_clinic(iblock),iblock)

      !*** subtract mean from weak restoring term

      SFWF_COMP(:,:,iblock,sfwf_comp_wrest) = &
      SFWF_COMP(:,:,iblock,sfwf_comp_wrest) - OCN_WGT(:,:,iblock)* &
                                     MASK_SR(:,:,iblock)*weak_mean

      !  if variable thickness surface layer, compute net surface
      !  freshwater flux (kg/m^2/s) due to restoring terms only
      !  then compute freshwater flux due to P-E and convert to (m/s)
      !  then set the tracer content in the freshwater flux
      !  defaults are FW*SST for tracer 1 (temperature)
      !             0 for salinity (really fresh water)
      !             0 for all other tracers
      !
      !  IF DATA ARE AVAILABLE...
      !  IMPLEMENT SUM OVER FRESHWATER TYPES (E,P,R,F,M) HERE:
      !
      !  TFW(:,:,n) = FW_EVAP*TW_EVAP(:,:,n) + FW_PRECIP*TW_PRECIP(:,:,n)
      !             + FW_ROFF*TW_ROFF(:,:,n) + FW_MELT*TW_MELT(:,:,n)
      !             + FW_FREEZE*TW_FREEZE(:,:,n)
      !
      !     where, for example  FW_ROFF is the water flux from rivers
      !     and TW_ROFF(:,:,n) is the concentration of the nth tracer
      !     in the river water;  similarly for water fluxes due to
      !     evaporation, precipitation, ice freezing, and ice melting.

      if (sfc_layer_type == sfc_layer_varthick .and. &
          .not. lfw_as_salt_flx) then

         STF(:,:,2,iblock) = SFWF_COMP(:,:,iblock,sfwf_comp_wrest) + &
                             SFWF_COMP(:,:,iblock,sfwf_comp_srest)

         FW(:,:,iblock) = OCN_WGT(:,:,iblock)*MASK_SR(:,:,iblock)*  &
                          (SFWF_COMP(:,:,iblock,sfwf_comp_evap) +   &
                           SFWF_COMP(:,:,iblock,sfwf_comp_precip))* &
                          fwmass_to_fwflux

         !*** fw same temp as ocean and no tracers in FW input
         TFW(:,:,1,iblock) = FW(:,:,iblock)* &
                             TRACER(:,:,1,1,curtime,iblock)
         TFW(:,:,2:nt,iblock) = c0

      !*** if rigid lid or old free surface form, compute
      !*** net surface freshwater flux (kg/m^2/s)

      else

         STF(:,:,2,iblock) = OCN_WGT(:,:,iblock)*MASK_SR(:,:,iblock)* &
                          (SFWF_COMP(:,:,iblock,sfwf_comp_evap) +     &
                           SFWF_COMP(:,:,iblock,sfwf_comp_precip)) +  &
                           SFWF_COMP(:,:,iblock,sfwf_comp_wrest)+     &
                           SFWF_COMP(:,:,iblock,sfwf_comp_srest)

      endif

      !*** convert surface freshwater flux (kg/m^2/s) to
      !*** salinity flux (msu*cm/s)

      STF(:,:,2,iblock) = STF(:,:,2,iblock)*salinity_factor

      !*** compute fields for accumulating annual-mean precipitation
      !*** over ocean points that are not marginal seas.

      WORK1(:,:,iblock) = merge(SFWF_COMP(:,:,iblock,sfwf_comp_precip)*&
                                TAREA(:,:,iblock)*OCN_WGT(:,:,iblock), &
                                c0, MASK_SR(:,:,iblock) > 0)

      !WORK2 = merge(FW_OLD(:,:,iblock)*TAREA(:,:,iblock), &
      !              c0, MASK_SR(:,:,iblock) > 0)

   end do ! block loop
   !$OMP END PARALLEL DO

!-----------------------------------------------------------------------
!
!  accumulate annual-mean precipitation over ocean points that are
!  not marginal seas.
!
!-----------------------------------------------------------------------

   if (avg_ts .or. back_to_back) then
      dttmp = p5*dtt
   else
      dttmp = dtt
   endif

   area_glob_m_marg = area_t - area_t_marg

   sum_precip = sum_precip + dttmp*1.0e-4_r8*  &
                global_sum(WORK1,distrb_clinic,field_loc_center)/ &
                area_glob_m_marg

   !sum_fw    = sum_fw     + &
   !            dttmp*global_sum(WORK2,distrb_clinic,field_loc_center)/ &
   !                  area_glob_m_marg

!-----------------------------------------------------------------------
!
!  Perform end of year adjustment calculations
!
!-----------------------------------------------------------------------

   if (eoy) then

      !*** Compute the surface volume-averaged salinity and
      !*** average surface height (for variable thickness sfc layer)
      !*** note that it is evaluated at the current time level.

      !$OMP PARALLEL DO PRIVATE(iblock, this_block)
      do iblock=1,nblocks_clinic
         this_block = get_block(blocks_clinic(iblock),iblock)

         if (sfc_layer_type == sfc_layer_varthick) then
            WORK1(:,:,iblock) = merge( &
                  TRACER(:,:,1,2,curtime,iblock)*TAREA(:,:,iblock)* &
                  (dz(1) + PSURF(:,:,curtime,iblock)/grav),         &
                 c0, KMT(:,:,iblock) > 0 .and. MASK_SR(:,:,iblock) > 0)
            WORK2(:,:,iblock) = merge(PSURF(:,:,curtime,iblock)*  &
                                      TAREA(:,:,iblock)/grav, c0, &
                                      KMT(:,:,iblock) > 0 .and.   &
                                      MASK_SR(:,:,iblock) > 0)
         else
            WORK1(:,:,iblock) = merge(TRACER(:,:,1,2,curtime,iblock)* &
                                      TAREA(:,:,iblock)*dz(1),        &
                                      c0, KMT(:,:,iblock) > 0 .and.   &
                                      MASK_SR(:,:,iblock) > 0)
         endif
      end do
      !$OMP END PARALLEL DO

      vol_glob_m_marg = volume_t_k(1) - volume_t_marg_k(1)
      sal_final(1) = 1.0e-6_r8*   &
                     global_sum(WORK1, distrb_clinic, field_loc_center)/ &
                     vol_glob_m_marg

      if (sfc_layer_type == sfc_layer_varthick) then
         ssh_final = 1.0e-4_r8*   &
                     global_sum(WORK2, distrb_clinic, field_loc_center)/ &
                     area_glob_m_marg/seconds_in_year
         if (my_task == master_task) then
            write(stdout,blank_fmt)
            write(stdout,'(a22,1pe23.15)') &
               'annual change in SSH: ', ssh_final
         endif
         ssh_final = ssh_final*10.0_r8 ! convert (cm/s) -> kg/m^2/s
                                       ! (cm/s)x0.01(m/cm)x1000kg/m^3
      else
         ssh_final = c0
      endif

      !*** Compute the volume-averaged salinity for each level.

      do k=2,km

         !$OMP PARALLEL DO PRIVATE(iblock, this_block)
         do iblock=1,nblocks_clinic
            this_block = get_block(blocks_clinic(iblock),iblock)

            if (partial_bottom_cells) then
               WORK1(:,:,iblock) =                               &
                  merge(TRACER(:,:,k,2,curtime,iblock)*          &
                        TAREA(:,:,iblock)*DZT(:,:,k,iblock), c0, &
                        k <= KMT(:,:,iblock) .and.               &
                        MASK_SR(:,:,iblock) > 0)
            else
               WORK1(:,:,iblock) =                               & 
                  merge(TRACER(:,:,k,2,curtime,iblock)*          &
                        TAREA(:,:,iblock)*dz(k), c0,             &
                        k <= KMT(:,:,iblock) .and.               &
                        MASK_SR(:,:,iblock) > 0)
            endif
         end do
         !$OMP END PARALLEL DO

         vol_glob_m_marg = volume_t_k(k) - volume_t_marg_k(k)
         if (vol_glob_m_marg == 0) vol_glob_m_marg = 1.e+20_r8
         sal_final(k) = 1.0e-6_r8*  &
                        global_sum(WORK1, distrb_clinic, field_loc_center)/ &
                        vol_glob_m_marg
      enddo

      !*** find annual mean precip and reset annual counters

      !ann_avg_fw = sum_fw / seconds_in_year
      !if (my_task == master_task) then
      !   write(stdout,blank_fmt)
      !   write(stdout,'(a32,1pe22.15)') &
      !     'annual average freshwater flux: ', ann_avg_fw
      !endif
      !sum_fw = c0

      ann_avg_precip = sum_precip / seconds_in_year
      sum_precip = c0

      if (ladjust_precip) call precip_adjustment

      sal_initial = sal_final
      ssh_initial = ssh_final

   endif ! end of year calculations

!-----------------------------------------------------------------------
!EOC

 end subroutine calc_sfwf_bulk_ncep

!***********************************************************************
!BOP
! !IROUTINE: calc_sfwf_partially_coupled
! !INTERFACE:


 subroutine calc_sfwf_partially_coupled(time_dim) 3,6

! !DESCRIPTION:
!  Calculate ice-ocean flux, weak restoring, and strong restoring 
!  components of surface freshwater flux for partially-coupled formulation.
!  these components will later be used in
!  set\_surface\_forcing (forcing.F) to form the total surface freshwater
!  (salt) flux.
!
!  the forcing data (on t-grid) sets needed are 
!     sfwf\_data\_sss,    restoring SSS                       (msu)
!     sfwf\_data\_flxio,  diagnosed ("climatological")        (kg/m^2/s)
!                       ice-ocean freshwater flux 

!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   integer (int_kind), intent(in) ::  &
      time_dim                ! number of time points for interpolation



!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      now,               &! index for location of interpolated data
      k, n,              &! dummy loop indices
      iblock              ! block loop index

   real (r8) ::           &
      dttmp,              &! temporary time step variable
      fres_hor_ave,       &! area-weighted mean of weak restoring
      fres_hor_area,      &! total area of weak restoring
      area_glob_m_marg,   &! total ocean area - marginal sea area
      vol_glob_m_marg      ! total ocean volume - marginal sea volume

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
      WORK,WORK1, WORK2   ! temporary work space

   type(block) :: &
      this_block  ! block info for current block


!-----------------------------------------------------------------------
!
!  set location where interpolated data exists.
!
!-----------------------------------------------------------------------

   area_glob_m_marg = 1.0e-4*(area_t - area_t_marg)

   if (sfwf_data_type == 'annual') then
      now = 1
   else
      now = 0
   endif

!-----------------------------------------------------------------------
!
!  compute forcing terms for each block
!
!-----------------------------------------------------------------------

   !$OMP PARALLEL DO PRIVATE(iblock, this_block)
   do iblock=1,nblocks_clinic
      this_block = get_block(blocks_clinic(iblock),iblock)

      if (.not. luse_cpl_ifrac) then
        WORK(:,:,iblock) = OCN_WGT(:,:,iblock)*MASK_SR(:,:,iblock)
      else
        WORK(:,:,iblock) = MASK_SR(:,:,iblock)
      endif

      !*** weak salinity restoring term
      !*** (note: MASK_SR = 0. at land and marginal-sea points)
      !*** will be subtracting global mean later, so compute
      !*** necessary terms for global mean

      SFWF_COMP(:,:,iblock,sfwf_comp_wrest) = &
                         -sfwf_weak_restore*WORK(:,:,iblock)*    &
                         (SFWF_DATA(:,:,iblock,sfwf_data_sss,now) - &
                          TRACER(:,:,1,2,curtime,iblock))

      WORK1(:,:,iblock) = TAREA(:,:,iblock)* &
                          SFWF_COMP(:,:,iblock,sfwf_comp_wrest)
      WORK2(:,:,iblock) = TAREA(:,:,iblock)*WORK(:,:,iblock)


   end do ! block loop
   !$OMP END PARALLEL DO

!----------------------------------------------------------------------
!
!  compute global mean of weak restoring term
!
!----------------------------------------------------------------------

   fres_hor_ave  = global_sum(WORK1, distrb_clinic, field_loc_center)
   fres_hor_area = global_sum(WORK2, distrb_clinic, field_loc_center)


!-----------------------------------------------------------------------
!
!  finish computing forcing terms for each block
!
!-----------------------------------------------------------------------

   !$OMP PARALLEL DO PRIVATE(iblock, this_block)
   do iblock=1,nblocks_clinic
      this_block = get_block(blocks_clinic(iblock),iblock)

    ! subtract global mean from weak restoring term
      SFWF_COMP(:,:,iblock,sfwf_comp_wrest) =  &
      SFWF_COMP(:,:,iblock,sfwf_comp_wrest) - WORK(:,:,iblock)*  &
      fres_hor_ave/fres_hor_area

      where (KMT(:,:,iblock) > 0)
         SFWF_COMP(:,:,iblock,sfwf_comp_srest) = &
                     -sfwf_strong_restore*(c1 - OCN_WGT(:,:,iblock))* &
                      (SFWF_DATA(:,:,iblock,sfwf_data_sss,now) -      &
                       TRACER(:,:,1,2,curtime,iblock))
      endwhere

      where (KMT(:,:,iblock) > 0 .and. MASK_SR(:,:,iblock) == 0)
         SFWF_COMP(:,:,iblock,sfwf_comp_srest) = &
                     -sfwf_strong_restore_ms*                         &
                      (SFWF_DATA(:,:,iblock,sfwf_data_sss,now) -      &
                       TRACER(:,:,1,2,curtime,iblock))
      endwhere


      !*** ice-ocean climatological flux term
      if ( .not. lactive_ice ) then
        where (KMT(:,:,iblock) > 0)
          SFWF_COMP(:,:,iblock,sfwf_comp_flxio) =       &
                           (c1 - OCN_WGT(:,:,iblock))*  &
                            SFWF_DATA(:,:,iblock,sfwf_data_flxio,now)
        endwhere
        if ( .not. lms_balance )  &
          SFWF_COMP(:,:,iblock,sfwf_comp_flxio) = &
          SFWF_COMP(:,:,iblock,sfwf_comp_flxio)*MASK_SR(:,:,iblock)
      endif

      !*** convert surface freshwater flux components (kg/m^2/s) to
      !*** salinity flux components (msu*cm/s)

      SFWF_COMP(:,:,iblock,sfwf_comp_wrest) =  &
                           SFWF_COMP(:,:,iblock,sfwf_comp_wrest)*salinity_factor
      SFWF_COMP(:,:,iblock,sfwf_comp_srest) =  &
                           SFWF_COMP(:,:,iblock,sfwf_comp_srest)*salinity_factor

      if ( sfc_layer_type == sfc_layer_varthick .and.             &
                                                .not. lfw_as_salt_flx) then
        WORK(:,:,iblock) = fwmass_to_fwflux *                     &
                           SFWF_COMP(:,:,iblock,sfwf_comp_flxio)
        SFWF_COMP(:,:,iblock,sfwf_comp_flxio) = WORK(:,:,iblock)

        call tmelt(WORK1(:,:,iblock),TRACER(:,:,1,2,curtime,iblock))

        TFW_COMP(:,:,1,   iblock,tfw_comp_flxio) = WORK(:,:,iblock)*  &
                                                   WORK1(:,:,iblock)
        TFW_COMP(:,:,2:nt,iblock,tfw_comp_flxio) = c0
      else
        SFWF_COMP(:,:,iblock,sfwf_comp_flxio)    = salinity_factor*  &
                           SFWF_COMP(:,:,iblock,sfwf_comp_flxio)
      endif


   end do ! block loop
   !$OMP END PARALLEL DO


!-----------------------------------------------------------------------
!
!  Perform end of year adjustment calculations
!
!-----------------------------------------------------------------------

   if (eoy) then

      !*** Compute the surface volume-averaged salinity and
      !*** average surface height (for variable thickness sfc layer)
      !*** note that it is evaluated at the current time level.

      !$OMP PARALLEL DO PRIVATE(iblock, this_block)
      do iblock=1,nblocks_clinic
         this_block = get_block(blocks_clinic(iblock),iblock)

         if (sfc_layer_type == sfc_layer_varthick) then
            WORK1(:,:,iblock) = merge( &
                  TRACER(:,:,1,2,curtime,iblock)*TAREA(:,:,iblock)* &
                  (dz(1) + PSURF(:,:,curtime,iblock)/grav),         &
                 c0, KMT(:,:,iblock) > 0 .and. MASK_SR(:,:,iblock) > 0)
            WORK2(:,:,iblock) = merge(PSURF(:,:,curtime,iblock)*  &
                                      TAREA(:,:,iblock)/grav, c0, &
                                      KMT(:,:,iblock) > 0 .and.   &
                                      MASK_SR(:,:,iblock) > 0)
         else
            WORK1(:,:,iblock) = merge(TRACER(:,:,1,2,curtime,iblock)* &
                                      TAREA(:,:,iblock)*dz(1),        &
                                      c0, KMT(:,:,iblock) > 0 .and.   &
                                      MASK_SR(:,:,iblock) > 0)
         endif
      end do
      !$OMP END PARALLEL DO

      vol_glob_m_marg = 1.0e-6_r8*(volume_t_k(1) - volume_t_marg_k(1))
      sal_final(1) = 1.0e-6_r8*  &! convert to m^3
                     global_sum(WORK1,distrb_clinic,field_loc_center)/vol_glob_m_marg
      if (sfc_layer_type == sfc_layer_varthick) then
         ssh_final = 1.0e-4_r8* &! convert to m^3
                     global_sum(WORK2, distrb_clinic, field_loc_center)/ &
                     area_glob_m_marg/seconds_in_year
         if (my_task == master_task) then
            write(stdout,blank_fmt)
            write(stdout,'(a22,1pe23.15)') &
               'annual change in SSH: ', ssh_final
         endif
         ssh_final = ssh_final*10.0_r8 ! convert (cm/s) -> kg/m^2/s
                                       ! (cm/s)x0.01(m/cm)x1000kg/m^3
      else
         ssh_final = c0
      endif

      !*** Compute the volume-averaged salinity for each level.

      do k=2,km

         !$OMP PARALLEL DO PRIVATE(iblock, this_block)
         do iblock=1,nblocks_clinic
            this_block = get_block(blocks_clinic(iblock),iblock)

            if (partial_bottom_cells) then
               WORK1(:,:,iblock) =                               &
                  merge(TRACER(:,:,k,2,curtime,iblock)*          &
                        TAREA(:,:,iblock)*DZT(:,:,k,iblock), c0, &
                        k <= KMT(:,:,iblock) .and.               &
                        MASK_SR(:,:,iblock) > 0)
            else
               WORK1(:,:,iblock) =                               & 
                  merge(TRACER(:,:,k,2,curtime,iblock)*          &
                        TAREA(:,:,iblock)*dz(k), c0,             &
                        k <= KMT(:,:,iblock) .and.               &
                        MASK_SR(:,:,iblock) > 0)
            endif
         end do
         !$OMP END PARALLEL DO

         vol_glob_m_marg = 1.0e-6_r8*(volume_t_k(k) - volume_t_marg_k(k))
         if (vol_glob_m_marg == 0) vol_glob_m_marg = 1.e+20_r8
         sal_final(k) = 1.0e-6_r8* &! convert to m^3
                        global_sum(WORK1,distrb_clinic,field_loc_center)/vol_glob_m_marg
      enddo


      if (ladjust_precip) call precip_adjustment

      sal_initial = sal_final
      ssh_initial = ssh_final

   endif ! end of year calculations

!-----------------------------------------------------------------------
!EOC

 end subroutine calc_sfwf_partially_coupled
!***********************************************************************
!BOP
! !IROUTINE: precip_adjustment
! !INTERFACE:


 subroutine precip_adjustment 2

! !DESCRIPTION:
!  Computes a precipitation factor to multiply the fresh water flux
!  due to precipitation uniformly to insure a balance of fresh water
!  at the ocean surface.
!
! !REVISION HISTORY:
!  same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

   real (r8) ::  &
      sal_tendency, fw_tendency, precip_tav, &
      area_glob_m_marg, &! global ocean area - marginal sea area (cm^2)
      vol_glob_m_marg    ! global ocean vol  - marginal sea vol  (cm^3)

  integer (int_kind) :: k

!-----------------------------------------------------------------------
!
!  compute tendency of salinity for eack "k" layer, considering the
!  effects of depth acceleration
!
!-----------------------------------------------------------------------

   do k=1,km
      sal_initial(k) = (sal_final(k) - sal_initial(k))/ &
                       (dttxcel(k)*seconds_in_year)
   enddo

!-----------------------------------------------------------------------
!
!  form the global volume-averaged tendency to be used in "precip_fact"
!  computation
!
!-----------------------------------------------------------------------

   sal_tendency = c0
   do k=1,km
      vol_glob_m_marg = 1.0e-6_r8*(volume_t_k(k) - volume_t_marg_k(k))
      sal_tendency = sal_tendency + vol_glob_m_marg*sal_initial(k)
   enddo

   vol_glob_m_marg = 1.0e-6_r8*(volume_t - volume_t_marg)
   sal_tendency = sal_tendency/vol_glob_m_marg

   if (my_task == master_task) then
      write (stdout,'(a58,1pe22.15)') &
         ' precip_adjustment: volume-averaged salinity tendency = ', &
         sal_tendency
   endif

!-----------------------------------------------------------------------
!
!  convert "sal_tendency" from (msu/s) to -(kg/m^2/s). note that
!  areag in cm^2 and volgt in cm^3
!  assumes density of fresh water = 1000 kg/m**3
!
!-----------------------------------------------------------------------

   area_glob_m_marg = 1.0e-4*(area_t - area_t_marg)
   sal_tendency = - sal_tendency*vol_glob_m_marg*  &
                    1.0e6_r8/area_glob_m_marg/ocn_ref_salinity

!-----------------------------------------------------------------------
!
!  compute annual change in mass due to freshwater flux (kg/m^2/s)
!
!-----------------------------------------------------------------------

   fw_tendency = ssh_final - ssh_initial

   if (my_task == master_task) then
      write (stdout,'(a22)') ' precip_adjustment: '
      write (stdout,'(a28,1pe22.15)') '   sal_tendency (kg/m^2/s): ', &
                                      sal_tendency
      write (stdout,'(a28,1pe22.15)') '    fw_tendency (kg/m^2/s): ', &
                                      fw_tendency
   endif

!-----------------------------------------------------------------------
!
!  change "precip_fact" based on tendency of freshwater and previous
!  amount of precipitation 
!
!-----------------------------------------------------------------------
   if (sfwf_formulation == 'partially-coupled') then
      precip_tav = precip_mean
   else
      precip_tav = ann_avg_precip/precip_fact
   endif

   precip_fact = precip_fact - &
                 (sal_tendency + fw_tendency)/precip_tav

   if (my_task == master_task) then
      write (stdout,'(a33,e14.8)') ' Changed precipitation factor to ',&
                                  precip_fact
   endif

!-----------------------------------------------------------------------
!EOC

 end subroutine precip_adjustment

!***********************************************************************


 end module forcing_sfwf

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||