!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module forcing_shf 10,11
!BOP
! !MODULE: forcing_shf
! !DESCRIPTION:
! Contains routines and variables used for determining the surface
! heat flux.
!
! !REVISION HISTORY:
! SVN:$Id: forcing_shf.F90 14725 2009-03-04 22:50:06Z njn01 $
!
! !USES:
use kinds_mod
use blocks
use distribution
use domain
use constants
use io
use grid
use forcing_tools
use time_management
use prognostic
use exit_mod
implicit none
private
save
! !PUBLIC MEMBER FUNCTIONS:
public :: init_shf, &
set_shf
! !PUBLIC DATA MEMBERS:
real (r8), dimension(nx_block,ny_block,max_blocks_clinic), &
public, target :: &
SHF_QSW, & ! incoming short wave
SHF_QSW_RAW ! no masking, no diurnal cycle
logical (log_kind), public :: &
lsw_absorb ! true if short wave available as separate flux
! (use penetrative short wave)
!*** the following must be shared with sfwf forcing in
!*** bulk-NCEP option
real (r8), allocatable, dimension(:,:,:,:), public :: &
SHF_COMP
real (r8), allocatable, dimension(:,:,:), public :: &
OCN_WGT
integer (int_kind), allocatable, dimension(:,:,:), public :: &
MASK_SR ! strong restoring mask for marginal seas
integer (int_kind), public :: &
shf_data_tair, &
shf_data_qair, &
shf_data_cldfrac, &
shf_data_windspd, &
shf_comp_qsw, &
shf_comp_qlw, &
shf_comp_qlat, &
shf_comp_qsens, &
shf_comp_wrest, &
shf_comp_srest, &
shf_comp_cpl
!*** the following are needed by restart
real (r8), public :: &
shf_interp_last ! time when last interpolation was done
!EOP
!BOC
!-----------------------------------------------------------------------
!
! module variables
!
!-----------------------------------------------------------------------
real (r8), allocatable, dimension(:,:,:,:,:) :: &
SHF_DATA ! forcing data to use for computing SHF
real (r8), dimension(12) :: &
shf_data_time ! time (hours) corresponding to surface heat fluxes
real (r8), dimension(20) :: &
shf_data_renorm ! factors for converting to model units
real (r8), parameter, private :: &
T_strong_restore_limit = -1.8_r8, &
T_weak_restore_limit = -0.8_r8, &
dT_restore_limit = T_weak_restore_limit - T_strong_restore_limit
real (r8) :: &
shf_data_inc, &! time increment between values of forcing data
shf_data_next, &! time that will be used for the next value of forcing data that is needed
shf_data_update, &! time when the a new forcing value needs to be added to interpolation set
shf_interp_inc, &! time increment between interpolation
shf_interp_next, &! time when next interpolation will be done
shf_restore_tau, &
shf_restore_rtau, &
shf_weak_restore, &! heat flux weak restoring max time scale
shf_strong_restore,&! heat flux strong restoring max time scale
shf_strong_restore_ms
integer (int_kind) :: &
shf_interp_order, &! order of temporal interpolation
shf_data_time_min_loc, &! time index for first shf_data point
shf_data_num_fields
integer (int_kind), public :: &
shf_num_comps
character (char_len), dimension(:), allocatable :: &
shf_data_names ! short names for input data fields
integer (int_kind), dimension(:), allocatable :: &
shf_bndy_loc, &! location and field type for ghost
shf_bndy_type ! cell updates
! the following is necessary for sst restoring and partially-coupled
integer (int_kind) :: &
shf_data_sst
! the following are necessary for Barnier-restoring
integer (int_kind) :: &
shf_data_tstar, &
shf_data_tau, &
shf_data_ice, &
shf_data_qsw
character (char_len) :: &
shf_interp_freq, &! keyword for period of temporal interpolation
shf_filename, &! file containing forcing data
shf_file_fmt, &! format (bin or netcdf) of shf file
shf_interp_type, &
shf_data_label
character (char_len), public :: &
shf_data_type, &! keyword for period of forcing data
shf_formulation
! the following is necessary for partially-coupled
! luse_cpl_ifrac = .T. use fractional ice coverage
! sent by the coupler from the (dummy) ice,
! .F. use fractional ice coverage based on the
! STR SST climatology.
logical (log_kind), public :: &
luse_cpl_ifrac
!-----------------------------------------------------------------------
!
! the following are needed for long-wave heat flux
! with bulk-NCEP forcing
!
!-----------------------------------------------------------------------
real (r8), allocatable, dimension (:,:,:) :: &
CCINT
real (r8), dimension(21) :: &
cc = (/ 0.88_r8, 0.84_r8, 0.80_r8, &
0.76_r8, 0.72_r8, 0.68_r8, &
0.63_r8, 0.59_r8, 0.52_r8, &
0.50_r8, 0.50_r8, 0.50_r8, &
0.52_r8, 0.59_r8, 0.63_r8, &
0.68_r8, 0.72_r8, 0.76_r8, &
0.80_r8, 0.84_r8, 0.88_r8 /)
real (r8), dimension(21) :: &
clat = (/ -90.0_r8, -80.0_r8, -70.0_r8, &
-60.0_r8, -50.0_r8, -40.0_r8, &
-30.0_r8, -20.0_r8, -10.0_r8, &
-5.0_r8, 0.0_r8, 5.0_r8, &
10.0_r8, 20.0_r8, 30.0_r8, &
40.0_r8, 50.0_r8, 60.0_r8, &
70.0_r8, 80.0_r8, 90.0_r8 /)
!EOC
!***********************************************************************
contains
!***********************************************************************
!BOP
! !IROUTINE: init_shf
! !INTERFACE:
subroutine init_shf(STF) 1,185
! !DESCRIPTION:
! Initializes surface heat flux forcing by either calculating
! or reading in the surface heat flux. Also do 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
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,nt,max_blocks_clinic), &
intent(out) :: &
STF ! surface tracer flux - this routine only modifies
! the slice corresponding to temperature (tracer 1)
!EOP
!BOC
!----------------------------------------------------------------------
!
! local variables
!
!----------------------------------------------------------------------
integer (int_kind) :: &
i,j, k, n, iblock, &! loop indices
nml_error ! namelist error flag
character (char_len) :: &
forcing_filename ! temp for full filename of forcing file
logical (log_kind) :: &
no_region_mask ! flag for existence of region mask
real (r8), dimension(:,:,:,:,:), allocatable :: &
TEMP_DATA ! temporary data array for monthly forcing
type (datafile) :: &
forcing_file ! file containing forcing fields
type (io_field_desc) :: & ! io descriptors for various input fields
io_sst, &
io_tstar, &
io_tau, &
io_ice, &
io_qsw, &
io_tair, &
io_qair, &
io_cldfrac, &
io_windspd
type (io_dim) :: &
i_dim, j_dim, &! dimension descriptors for horiz dims
month_dim ! dimension descriptor for monthly data
namelist /forcing_shf_nml/ shf_data_type, shf_data_inc, &
shf_interp_type, shf_interp_freq, &
shf_interp_inc, shf_restore_tau, &
shf_filename, shf_file_fmt, &
shf_data_renorm, &
shf_formulation, &
shf_weak_restore, shf_strong_restore,&
shf_strong_restore_ms, &
luse_cpl_ifrac
!-----------------------------------------------------------------------
!
! read surface heat flux namelist input after setting default values.
!
!-----------------------------------------------------------------------
shf_formulation = 'restoring'
shf_data_type = 'analytic'
shf_data_inc = 1.e20_r8
shf_interp_type = 'nearest'
shf_interp_freq = 'never'
shf_interp_inc = 1.e20_r8
shf_restore_tau = 1.e20_r8
shf_filename = 'unknown-shf'
shf_file_fmt = 'bin'
shf_data_renorm = c1
shf_weak_restore = c0
shf_strong_restore = 92.64_r8
shf_strong_restore_ms = 92.64_r8
luse_cpl_ifrac = .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_shf_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_shf_nml')
endif
call broadcast_scalar(shf_formulation, master_task)
call broadcast_scalar(shf_data_type, master_task)
call broadcast_scalar(shf_data_inc, master_task)
call broadcast_scalar(shf_interp_type, master_task)
call broadcast_scalar(shf_interp_freq, master_task)
call broadcast_scalar(shf_interp_inc, master_task)
call broadcast_scalar(shf_restore_tau, master_task)
call broadcast_scalar(shf_filename, master_task)
call broadcast_scalar(shf_file_fmt, master_task)
call broadcast_array (shf_data_renorm, master_task)
call broadcast_scalar(shf_weak_restore, master_task)
call broadcast_scalar(shf_strong_restore, master_task)
call broadcast_scalar(shf_strong_restore_ms, master_task)
call broadcast_scalar(luse_cpl_ifrac, master_task)
!-----------------------------------------------------------------------
!
! convert data_type to 'monthly-calendar' if input is 'monthly'
!
!-----------------------------------------------------------------------
if (shf_data_type == 'monthly') shf_data_type = 'monthly-calendar'
!-----------------------------------------------------------------------
!
! set values based on shf_formulation
!
!-----------------------------------------------------------------------
select case (shf_formulation)
case ('restoring')
lsw_absorb = .false.
shf_data_num_fields = 1
allocate(shf_data_names(shf_data_num_fields), &
shf_bndy_loc (shf_data_num_fields), &
shf_bndy_type (shf_data_num_fields))
shf_data_sst = 1
shf_data_names(shf_data_sst) = 'SST'
shf_bndy_loc (shf_data_sst) = field_loc_center
shf_bndy_type (shf_data_sst) = field_type_scalar
case ('Barnier-restoring')
lsw_absorb = .true.
shf_data_num_fields = 4
allocate(shf_data_names(shf_data_num_fields), &
shf_bndy_loc (shf_data_num_fields), &
shf_bndy_type (shf_data_num_fields))
shf_data_tstar = 1
shf_data_tau = 2
shf_data_ice = 3
shf_data_qsw = 4
shf_data_names(shf_data_tstar) = 'TSTAR'
shf_bndy_loc (shf_data_tstar) = field_loc_center
shf_bndy_type (shf_data_tstar) = field_type_scalar
shf_data_names(shf_data_tau) = 'TAU'
shf_bndy_loc (shf_data_tau) = field_loc_center
shf_bndy_type (shf_data_tau) = field_type_scalar
shf_data_names(shf_data_ice) = 'ICE'
shf_bndy_loc (shf_data_ice) = field_loc_center
shf_bndy_type (shf_data_ice) = field_type_scalar
shf_data_names(shf_data_qsw) = 'QSW'
shf_bndy_loc (shf_data_qsw) = field_loc_center
shf_bndy_type (shf_data_qsw) = field_type_scalar
case ('bulk-NCEP')
lsw_absorb = .true.
shf_data_num_fields = 6
allocate(shf_data_names(shf_data_num_fields), &
shf_bndy_loc (shf_data_num_fields), &
shf_bndy_type (shf_data_num_fields))
shf_data_sst = 1
shf_data_tair = 2
shf_data_qair = 3
shf_data_qsw = 4
shf_data_cldfrac = 5
shf_data_windspd = 6
shf_data_names(shf_data_sst) = 'SST'
shf_bndy_loc (shf_data_sst) = field_loc_center
shf_bndy_type (shf_data_sst) = field_type_scalar
shf_data_names(shf_data_tair) = 'TAIR'
shf_bndy_loc (shf_data_tair) = field_loc_center
shf_bndy_type (shf_data_tair) = field_type_scalar
shf_data_names(shf_data_qair) = 'QAIR'
shf_bndy_loc (shf_data_qair) = field_loc_center
shf_bndy_type (shf_data_qair) = field_type_scalar
shf_data_names(shf_data_qsw) = 'QSW'
shf_bndy_loc (shf_data_qsw) = field_loc_center
shf_bndy_type (shf_data_qsw) = field_type_scalar
shf_data_names(shf_data_cldfrac) = 'CLDFRAC'
shf_bndy_loc (shf_data_cldfrac) = field_loc_center
shf_bndy_type (shf_data_cldfrac) = field_type_scalar
shf_data_names(shf_data_windspd) = 'WINDSPD'
shf_bndy_loc (shf_data_windspd) = field_loc_center
shf_bndy_type (shf_data_windspd) = field_type_scalar
shf_num_comps = 6
shf_comp_qsw = 1
shf_comp_qlw = 2
shf_comp_qlat = 3
shf_comp_qsens = 4
shf_comp_wrest = 5
shf_comp_srest = 6
!*** initialize CCINT (cloud factor used in long-wave heat flux
!*** with bulk-NCEP forcing).
allocate(CCINT(nx_block,ny_block,max_blocks_clinic))
!$OMP PARALLEL DO PRIVATE(iblock,i,j)
do iblock=1,nblocks_clinic
do j=1,ny_block
do i=1,20
where ((TLAT(:,j,iblock)*radian > clat(i )) .and. &
(TLAT(:,j,iblock)*radian <= clat(i+1)))
CCINT(:,j,iblock) = cc(i) + (cc(i+1)-cc(i))* &
(TLAT(:,j,iblock)*radian - clat(i))/ &
(clat(i+1)-clat(i))
endwhere
end do ! i
end do ! j
end do ! block loop
!$OMP END PARALLEL DO
case ('partially-coupled')
lsw_absorb = .false.
shf_data_num_fields = 1
allocate(shf_data_names(shf_data_num_fields), &
shf_bndy_loc (shf_data_num_fields), &
shf_bndy_type (shf_data_num_fields))
shf_data_sst = 1
shf_data_names(shf_data_sst) = 'SST'
shf_bndy_loc (shf_data_sst) = field_loc_center
shf_bndy_type (shf_data_sst) = field_type_scalar
shf_num_comps = 4
shf_comp_wrest = 1
shf_comp_srest = 2
shf_comp_cpl = 3
shf_comp_qsw = 4
case default
call exit_POP(sigAbort, &
'init_shf: Unknown value for shf_formulation')
end select
!-----------------------------------------------------------------------
!
! calculate inverse of restoring time scale and convert to seconds.
!
!-----------------------------------------------------------------------
shf_restore_tau = seconds_in_day*shf_restore_tau
shf_restore_rtau = c1/shf_restore_tau
!-----------------------------------------------------------------------
!
! initialize SHF_QSW in case a value is needed but not
! supplied by data: for example, with KPP and restoring.
!
!-----------------------------------------------------------------------
SHF_QSW = c0
SHF_QSW_RAW = c0
!-----------------------------------------------------------------------
!
! set strong restoring mask to 0 only at ocean points that are
! marginal seas and land.
!
!-----------------------------------------------------------------------
if (allocated(REGION_MASK)) then
allocate( MASK_SR(nx_block,ny_block,max_blocks_clinic))
no_region_mask = .false.
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock=1,nblocks_clinic
MASK_SR(:,:,iblock) = merge(0, 1, &
REGION_MASK(:,:,iblock) <= 0)
end do
!$OMP END PARALLEL DO
else
no_region_mask = .true.
endif
!-----------------------------------------------------------------------
!
! convert interp_type to corresponding integer value.
!
!-----------------------------------------------------------------------
select case (shf_interp_type)
case ('nearest')
shf_interp_order = 1
case ('linear')
shf_interp_order = 2
case ('4point')
shf_interp_order = 4
case default
call exit_POP(sigAbort, &
'init_shf: Unknown value for shf_interp_type')
end select
!-----------------------------------------------------------------------
!
! set values of the surface heat flux arrays (STF or SHF_DATA)
! depending on the type of the surface heat flux data.
!
!-----------------------------------------------------------------------
select case (shf_data_type)
!-----------------------------------------------------------------------
!
! no surface heat flux, therefore no interpolation in time
! needed, nor are there any new values to be used.
!
!-----------------------------------------------------------------------
case ('none')
STF(:,:,1,:) = c0
shf_data_next = never
shf_data_update = never
shf_interp_freq = 'never'
!-----------------------------------------------------------------------
!
! simple analytic surface temperature that is constant in
! time, therefore no new values will be needed.
!
!-----------------------------------------------------------------------
case ('analytic')
allocate( SHF_DATA(nx_block,ny_block,max_blocks_clinic, &
shf_data_num_fields,1))
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock=1,nblocks_clinic
select case (shf_formulation)
case ('restoring')
SHF_DATA(:,:,iblock,shf_data_sst,1) = &
28.0_r8*(c1 - sin(ULAT(:,:,iblock)))
end select
end do ! block loop
!$OMP END PARALLEL DO
shf_data_next = never
shf_data_update = never
shf_interp_freq = 'never'
!-----------------------------------------------------------------------
!
! annual mean climatological surface temperature (read in from file)
! that is constant in time, therefore no new values will be needed
! (shf_data_next = shf_data_update = never).
!
!-----------------------------------------------------------------------
case ('annual')
allocate( SHF_DATA(nx_block,ny_block,max_blocks_clinic, &
shf_data_num_fields,1))
SHF_DATA = c0
forcing_file = construct_file(shf_file_fmt, &
full_name=trim(shf_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',ny_global)
select case (shf_formulation)
case ('restoring')
io_sst = construct_io_field( &
trim(shf_data_names(shf_data_sst)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_sst), &
field_type = shf_bndy_type(shf_data_sst), &
d2d_array=SHF_DATA(:,:,:,shf_data_sst,1))
call data_set(forcing_file,'define',io_sst)
call data_set(forcing_file,'read' ,io_sst)
call destroy_io_field(io_sst)
case ('partially-coupled')
io_sst = construct_io_field( &
trim(shf_data_names(shf_data_sst)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_sst), &
field_type = shf_bndy_type(shf_data_sst), &
d2d_array=SHF_DATA(:,:,:,shf_data_sst,1))
call data_set(forcing_file,'define',io_sst)
call data_set(forcing_file,'read' ,io_sst)
call destroy_io_field(io_sst)
allocate(SHF_COMP(nx_block,ny_block,max_blocks_clinic, &
shf_num_comps), &
OCN_WGT (nx_block,ny_block,max_blocks_clinic))
SHF_COMP = c0 ! initialize
case ('Barnier-restoring')
io_tstar = construct_io_field( &
trim(shf_data_names(shf_data_tstar)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_tstar), &
field_type = shf_bndy_type(shf_data_tstar), &
d2d_array=SHF_DATA(:,:,:,shf_data_tstar,1))
io_tau = construct_io_field( &
trim(shf_data_names(shf_data_tau)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_tau), &
field_type = shf_bndy_type(shf_data_tau), &
d2d_array=SHF_DATA(:,:,:,shf_data_tau,1))
io_ice = construct_io_field( &
trim(shf_data_names(shf_data_ice)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_ice), &
field_type = shf_bndy_type(shf_data_ice), &
d2d_array=SHF_DATA(:,:,:,shf_data_ice,1))
io_qsw = construct_io_field( &
trim(shf_data_names(shf_data_qsw)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_qsw), &
field_type = shf_bndy_type(shf_data_qsw), &
d2d_array=SHF_DATA(:,:,:,shf_data_qsw,1))
call data_set(forcing_file,'define',io_tstar)
call data_set(forcing_file,'define',io_tau)
call data_set(forcing_file,'define',io_ice)
call data_set(forcing_file,'define',io_qsw)
call data_set(forcing_file,'read',io_tstar)
call data_set(forcing_file,'read',io_tau)
call data_set(forcing_file,'read',io_ice)
call data_set(forcing_file,'read',io_qsw)
call destroy_io_field(io_tstar)
call destroy_io_field(io_tau)
call destroy_io_field(io_ice)
call destroy_io_field(io_qsw)
SHF_DATA(:,:,:,shf_data_tau,1) = seconds_in_day* &
SHF_DATA(:,:,:,shf_data_tau,1)
case ('bulk-NCEP')
io_sst = construct_io_field( &
trim(shf_data_names(shf_data_sst)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_sst), &
field_type = shf_bndy_type(shf_data_sst), &
d2d_array=SHF_DATA(:,:,:,shf_data_sst,1))
io_tair = construct_io_field( &
trim(shf_data_names(shf_data_tair)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_tair), &
field_type = shf_bndy_type(shf_data_tair), &
d2d_array=SHF_DATA(:,:,:,shf_data_tair,1))
io_qair = construct_io_field( &
trim(shf_data_names(shf_data_qair)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_qair), &
field_type = shf_bndy_type(shf_data_qair), &
d2d_array=SHF_DATA(:,:,:,shf_data_qair,1))
io_qsw = construct_io_field( &
trim(shf_data_names(shf_data_qsw)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_qsw), &
field_type = shf_bndy_type(shf_data_qsw), &
d2d_array=SHF_DATA(:,:,:,shf_data_qsw,1))
io_cldfrac = construct_io_field( &
trim(shf_data_names(shf_data_cldfrac)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_cldfrac), &
field_type = shf_bndy_type(shf_data_cldfrac), &
d2d_array=SHF_DATA(:,:,:,shf_data_cldfrac,1))
io_windspd = construct_io_field( &
trim(shf_data_names(shf_data_windspd)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_windspd), &
field_type = shf_bndy_type(shf_data_windspd), &
d2d_array=SHF_DATA(:,:,:,shf_data_windspd,1))
call data_set(forcing_file, 'define', io_sst)
call data_set(forcing_file, 'define', io_tair)
call data_set(forcing_file, 'define', io_qair)
call data_set(forcing_file, 'define', io_qsw)
call data_set(forcing_file, 'define', io_cldfrac)
call data_set(forcing_file, 'define', io_windspd)
call data_set(forcing_file, 'read', io_sst)
call data_set(forcing_file, 'read', io_tair)
call data_set(forcing_file, 'read', io_qair)
call data_set(forcing_file, 'read', io_qsw)
call data_set(forcing_file, 'read', io_cldfrac)
call data_set(forcing_file, 'read', io_windspd)
call destroy_io_field(io_sst)
call destroy_io_field(io_tair)
call destroy_io_field(io_qair)
call destroy_io_field(io_qsw)
call destroy_io_field(io_cldfrac)
call destroy_io_field(io_windspd)
allocate(SHF_COMP(nx_block,ny_block,max_blocks_clinic, &
shf_num_comps), &
OCN_WGT (nx_block,ny_block,max_blocks_clinic))
SHF_COMP = c0 ! initialize
end select
call data_set(forcing_file,'close')
!*** renormalize values if necessary to compensate for different
!*** units
do n = 1,shf_data_num_fields
if (shf_data_renorm(n) /= c1) SHF_DATA(:,:,:,n,:) = &
shf_data_renorm(n)*SHF_DATA(:,:,:,n,:)
enddo
shf_data_next = never
shf_data_update = never
shf_interp_freq = 'never'
if (my_task == master_task) then
write(stdout,blank_fmt)
write(stdout,'(a23,a)') ' SHF Annual file read: ', &
trim(forcing_file%full_name)
endif
call destroy_file(forcing_file)
!-----------------------------------------------------------------------
! monthly mean climatological surface heat flux. all
! 12 months are read in from a file. interpolation order
! (shf_interp_order) may be specified with namelist input.
!-----------------------------------------------------------------------
case ('monthly-equal','monthly-calendar')
allocate(SHF_DATA(nx_block,ny_block,max_blocks_clinic, &
shf_data_num_fields,0:12), &
TEMP_DATA(nx_block,ny_block,12,max_blocks_clinic, &
shf_data_num_fields))
SHF_DATA = c0
call find_forcing_times(shf_data_time, shf_data_inc, &
shf_interp_type, shf_data_next, &
shf_data_time_min_loc, shf_data_update, &
shf_data_type)
forcing_file = construct_file(shf_file_fmt, &
full_name = trim(shf_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',ny_global)
month_dim = construct_io_dim('month',12)
select case (shf_formulation)
case ('restoring')
io_sst = construct_io_field( &
trim(shf_data_names(shf_data_sst)), &
dim1=i_dim, dim2=j_dim, dim3=month_dim, &
field_loc = shf_bndy_loc(shf_data_sst), &
field_type = shf_bndy_type(shf_data_sst), &
d3d_array=TEMP_DATA(:,:,:,:,shf_data_sst))
call data_set(forcing_file,'define',io_sst)
call data_set(forcing_file,'read' ,io_sst)
call destroy_io_field(io_sst)
!$OMP PARALLEL DO PRIVATE(iblock, n)
do iblock=1,nblocks_clinic
do n=1,12
SHF_DATA (:,:,iblock,shf_data_sst,n) = &
TEMP_DATA(:,:,n,iblock,shf_data_sst)
end do
end do
!$OMP END PARALLEL DO
case ('partially-coupled')
io_sst = construct_io_field( &
trim(shf_data_names(shf_data_sst)), &
dim1=i_dim, dim2=j_dim, dim3=month_dim, &
field_loc = shf_bndy_loc(shf_data_sst), &
field_type = shf_bndy_type(shf_data_sst), &
d3d_array=TEMP_DATA(:,:,:,:,shf_data_sst))
call data_set(forcing_file,'define',io_sst)
call data_set(forcing_file,'read' ,io_sst)
call destroy_io_field(io_sst)
!$OMP PARALLEL DO PRIVATE(iblock, n)
do iblock=1,nblocks_clinic
do n=1,12
SHF_DATA (:,:,iblock,shf_data_sst,n) = &
TEMP_DATA(:,:,n,iblock,shf_data_sst)
end do
end do
!$OMP END PARALLEL DO
allocate(SHF_COMP(nx_block,ny_block,max_blocks_clinic, &
shf_num_comps), &
OCN_WGT (nx_block,ny_block,max_blocks_clinic))
SHF_COMP = c0 ! initialize
case ('Barnier-restoring')
io_tstar = construct_io_field( &
trim(shf_data_names(shf_data_tstar)), &
dim1=i_dim, dim2=j_dim, dim3=month_dim, &
field_loc = shf_bndy_loc(shf_data_tstar), &
field_type = shf_bndy_type(shf_data_tstar), &
d3d_array=TEMP_DATA(:,:,:,:,shf_data_tstar))
io_tau = construct_io_field( &
trim(shf_data_names(shf_data_tau)), &
dim1=i_dim, dim2=j_dim, dim3=month_dim, &
field_loc = shf_bndy_loc(shf_data_tau), &
field_type = shf_bndy_type(shf_data_tau), &
d3d_array=TEMP_DATA(:,:,:,:,shf_data_tau))
io_ice = construct_io_field( &
trim(shf_data_names(shf_data_ice)), &
dim1=i_dim, dim2=j_dim, dim3=month_dim, &
field_loc = shf_bndy_loc(shf_data_ice), &
field_type = shf_bndy_type(shf_data_ice), &
d3d_array=TEMP_DATA(:,:,:,:,shf_data_ice))
io_qsw = construct_io_field( &
trim(shf_data_names(shf_data_qsw)), &
dim1=i_dim, dim2=j_dim, dim3=month_dim, &
field_loc = shf_bndy_loc(shf_data_qsw), &
field_type = shf_bndy_type(shf_data_qsw), &
d3d_array=TEMP_DATA(:,:,:,:,shf_data_qsw))
call data_set(forcing_file,'define',io_tstar)
call data_set(forcing_file,'define',io_tau)
call data_set(forcing_file,'define',io_ice)
call data_set(forcing_file,'define',io_qsw)
call data_set(forcing_file,'read',io_tstar)
call data_set(forcing_file,'read',io_tau)
call data_set(forcing_file,'read',io_ice)
call data_set(forcing_file,'read',io_qsw)
!$OMP PARALLEL DO PRIVATE(iblock,n)
do iblock=1,nblocks_clinic
do n=1,12
SHF_DATA (:,:,iblock,shf_data_tstar,n) = &
TEMP_DATA(:,:,n,iblock,shf_data_tstar)
SHF_DATA (:,:,iblock,shf_data_tau,n) = &
TEMP_DATA(:,:,n,iblock,shf_data_tau)*seconds_in_day
SHF_DATA (:,:,iblock,shf_data_ice,n) = &
TEMP_DATA(:,:,n,iblock,shf_data_ice)
SHF_DATA (:,:,iblock,shf_data_qsw,n) = &
TEMP_DATA(:,:,n,iblock,shf_data_qsw)
end do
end do
!$OMP END PARALLEL DO
call destroy_io_field(io_tstar)
call destroy_io_field(io_tau)
call destroy_io_field(io_ice)
call destroy_io_field(io_qsw)
case ('bulk-NCEP')
io_sst = construct_io_field( &
trim(shf_data_names(shf_data_sst)), &
dim1=i_dim, dim2=j_dim, dim3=month_dim, &
field_loc = shf_bndy_loc(shf_data_sst), &
field_type = shf_bndy_type(shf_data_sst), &
d3d_array=TEMP_DATA(:,:,:,:,shf_data_sst))
io_tair = construct_io_field( &
trim(shf_data_names(shf_data_tair)), &
dim1=i_dim, dim2=j_dim, dim3=month_dim, &
field_loc = shf_bndy_loc(shf_data_tair), &
field_type = shf_bndy_type(shf_data_tair), &
d3d_array=TEMP_DATA(:,:,:,:,shf_data_tair))
io_qair = construct_io_field( &
trim(shf_data_names(shf_data_qair)), &
dim1=i_dim, dim2=j_dim, dim3=month_dim, &
field_loc = shf_bndy_loc(shf_data_qair), &
field_type = shf_bndy_type(shf_data_qair), &
d3d_array=TEMP_DATA(:,:,:,:,shf_data_qair))
io_qsw = construct_io_field( &
trim(shf_data_names(shf_data_qsw)), &
dim1=i_dim, dim2=j_dim, dim3=month_dim, &
field_loc = shf_bndy_loc(shf_data_qsw), &
field_type = shf_bndy_type(shf_data_qsw), &
d3d_array=TEMP_DATA(:,:,:,:,shf_data_qsw))
io_cldfrac = construct_io_field( &
trim(shf_data_names(shf_data_cldfrac)), &
dim1=i_dim, dim2=j_dim, dim3=month_dim, &
field_loc = shf_bndy_loc(shf_data_cldfrac), &
field_type = shf_bndy_type(shf_data_cldfrac), &
d3d_array=TEMP_DATA(:,:,:,:,shf_data_cldfrac))
io_windspd = construct_io_field( &
trim(shf_data_names(shf_data_windspd)), &
dim1=i_dim, dim2=j_dim, dim3=month_dim, &
field_loc = shf_bndy_loc(shf_data_windspd), &
field_type = shf_bndy_type(shf_data_windspd), &
d3d_array=TEMP_DATA(:,:,:,:,shf_data_windspd))
call data_set(forcing_file, 'define', io_sst)
call data_set(forcing_file, 'define', io_tair)
call data_set(forcing_file, 'define', io_qair)
call data_set(forcing_file, 'define', io_qsw)
call data_set(forcing_file, 'define', io_cldfrac)
call data_set(forcing_file, 'define', io_windspd)
call data_set(forcing_file, 'read', io_sst)
call data_set(forcing_file, 'read', io_tair)
call data_set(forcing_file, 'read', io_qair)
call data_set(forcing_file, 'read', io_qsw)
call data_set(forcing_file, 'read', io_cldfrac)
call data_set(forcing_file, 'read', io_windspd)
!$OMP PARALLEL DO PRIVATE(iblock, n)
do iblock=1,nblocks_clinic
do n=1,12
SHF_DATA (:,:,iblock,shf_data_sst,n) = &
TEMP_DATA(:,:,n,iblock,shf_data_sst)
SHF_DATA (:,:,iblock,shf_data_tair,n) = &
TEMP_DATA(:,:,n,iblock,shf_data_tair)
SHF_DATA (:,:,iblock,shf_data_qair,n) = &
TEMP_DATA(:,:,n,iblock,shf_data_qair)
SHF_DATA (:,:,iblock,shf_data_qsw,n) = &
TEMP_DATA(:,:,n,iblock,shf_data_qsw)
SHF_DATA (:,:,iblock,shf_data_cldfrac,n) = &
TEMP_DATA(:,:,n,iblock,shf_data_cldfrac)
SHF_DATA (:,:,iblock,shf_data_windspd,n) = &
TEMP_DATA(:,:,n,iblock,shf_data_windspd)
end do
end do
!$OMP END PARALLEL DO
call destroy_io_field(io_sst)
call destroy_io_field(io_tair)
call destroy_io_field(io_qair)
call destroy_io_field(io_qsw)
call destroy_io_field(io_cldfrac)
call destroy_io_field(io_windspd)
allocate( SHF_COMP(nx_block,ny_block,max_blocks_clinic, &
shf_num_comps), &
OCN_WGT(nx_block,ny_block,max_blocks_clinic))
SHF_COMP = c0 ! initialize
end select
deallocate(TEMP_DATA)
call data_set(forcing_file,'close')
call destroy_file(forcing_file)
!*** renormalize values if necessary to compensate for different
!*** units.
do n = 1,shf_data_num_fields
if (shf_data_renorm(n) /= c1) SHF_DATA(:,:,:,n,:) = &
shf_data_renorm(n)*SHF_DATA(:,:,:,n,:)
enddo
if (my_task == master_task) then
write(stdout,blank_fmt)
write(stdout,'(a24,a)') ' SHF Monthly file read: ', &
trim(shf_filename)
endif
!-----------------------------------------------------------------------
!
! surface temperature specified every n-hours, where the n-hour
! increment should be specified with namelist input
! (shf_data_inc). only as many times as are necessary based on
! the order of the temporal interpolation scheme
! (shf_interp_order) reside in memory at any given time.
!
!-----------------------------------------------------------------------
case ('n-hour')
allocate(SHF_DATA(nx_block,ny_block,max_blocks_clinic, &
shf_data_num_fields,0:shf_interp_order))
SHF_DATA = c0
call find_forcing_times(shf_data_time, shf_data_inc, &
shf_interp_type, shf_data_next, &
shf_data_time_min_loc, shf_data_update, &
shf_data_type)
do n = 1, shf_interp_order
call get_forcing_filename(forcing_filename, shf_filename, &
shf_data_time(n), shf_data_inc)
forcing_file = construct_file(shf_file_fmt, &
full_name=trim(forcing_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',ny_global)
select case (shf_formulation)
case ('restoring','partially-coupled')
io_sst = construct_io_field( &
trim(shf_data_names(shf_data_sst)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_sst), &
field_type = shf_bndy_type(shf_data_sst), &
d2d_array=SHF_DATA(:,:,:,shf_data_sst,n))
call data_set(forcing_file,'define',io_sst)
call data_set(forcing_file,'read' ,io_sst)
call destroy_io_field(io_sst)
case ('Barnier-restoring')
io_tstar = construct_io_field( &
trim(shf_data_names(shf_data_tstar)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_tstar), &
field_type = shf_bndy_type(shf_data_tstar), &
d2d_array=SHF_DATA(:,:,:,shf_data_tstar,n))
io_tau = construct_io_field( &
trim(shf_data_names(shf_data_tau)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_tau), &
field_type = shf_bndy_type(shf_data_tau), &
d2d_array=SHF_DATA(:,:,:,shf_data_tau ,n))
io_ice = construct_io_field( &
trim(shf_data_names(shf_data_ice)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_ice), &
field_type = shf_bndy_type(shf_data_ice), &
d2d_array=SHF_DATA(:,:,:,shf_data_ice ,n))
io_qsw = construct_io_field( &
trim(shf_data_names(shf_data_qsw)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_qsw), &
field_type = shf_bndy_type(shf_data_qsw), &
d2d_array=SHF_DATA(:,:,:,shf_data_qsw ,n))
call data_set(forcing_file,'define',io_tstar)
call data_set(forcing_file,'define',io_tau)
call data_set(forcing_file,'define',io_ice)
call data_set(forcing_file,'define',io_qsw)
call data_set(forcing_file,'read',io_tstar)
call data_set(forcing_file,'read',io_tau)
call data_set(forcing_file,'read',io_ice)
call data_set(forcing_file,'read',io_qsw)
call destroy_io_field(io_tstar)
call destroy_io_field(io_tau)
call destroy_io_field(io_ice)
call destroy_io_field(io_qsw)
SHF_DATA(:,:,:,shf_data_tau ,n) = &
SHF_DATA(:,:,:,shf_data_tau ,n)*seconds_in_day
case ('bulk-NCEP')
io_sst = construct_io_field( &
trim(shf_data_names(shf_data_sst)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_sst), &
field_type = shf_bndy_type(shf_data_sst), &
d2d_array=SHF_DATA(:,:,:,shf_data_sst,n))
io_tair = construct_io_field( &
trim(shf_data_names(shf_data_tair)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_tair), &
field_type = shf_bndy_type(shf_data_tair), &
d2d_array=SHF_DATA(:,:,:,shf_data_tair,n))
io_qair = construct_io_field( &
trim(shf_data_names(shf_data_qair)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_qair), &
field_type = shf_bndy_type(shf_data_qair), &
d2d_array=SHF_DATA(:,:,:,shf_data_qair,n))
io_qsw = construct_io_field( &
trim(shf_data_names(shf_data_qsw)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_qsw), &
field_type = shf_bndy_type(shf_data_qsw), &
d2d_array=SHF_DATA(:,:,:,shf_data_qsw,n))
io_cldfrac = construct_io_field( &
trim(shf_data_names(shf_data_cldfrac)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_cldfrac), &
field_type = shf_bndy_type(shf_data_cldfrac), &
d2d_array=SHF_DATA(:,:,:,shf_data_cldfrac,n))
io_windspd = construct_io_field( &
trim(shf_data_names(shf_data_windspd)), &
dim1=i_dim, dim2=j_dim, &
field_loc = shf_bndy_loc(shf_data_windspd), &
field_type = shf_bndy_type(shf_data_windspd), &
d2d_array=SHF_DATA(:,:,:,shf_data_windspd,n))
call data_set(forcing_file, 'define', io_sst)
call data_set(forcing_file, 'define', io_tair)
call data_set(forcing_file, 'define', io_qair)
call data_set(forcing_file, 'define', io_qsw)
call data_set(forcing_file, 'define', io_cldfrac)
call data_set(forcing_file, 'define', io_windspd)
call data_set(forcing_file, 'read', io_sst)
call data_set(forcing_file, 'read', io_tair)
call data_set(forcing_file, 'read', io_qair)
call data_set(forcing_file, 'read', io_qsw)
call data_set(forcing_file, 'read', io_cldfrac)
call data_set(forcing_file, 'read', io_windspd)
call destroy_io_field(io_sst)
call destroy_io_field(io_tair)
call destroy_io_field(io_qair)
call destroy_io_field(io_qsw)
call destroy_io_field(io_cldfrac)
call destroy_io_field(io_windspd)
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,'(a23,a)') ' SHF n-hour file read: ', &
trim(forcing_filename)
endif
enddo
if (shf_formulation == 'bulk-NCEP' .or. &
shf_formulation == 'partially-coupled') then
allocate(SHF_COMP(nx_block,ny_block,max_blocks_clinic, &
shf_num_comps), &
OCN_WGT(nx_block,ny_block,max_blocks_clinic))
SHF_COMP = c0 ! initialize
endif
!*** renormalize values if necessary to compensate for different
!*** units.
do n = 1,shf_data_num_fields
if (shf_data_renorm(n) /= c1) SHF_DATA(:,:,:,n,:) = &
shf_data_renorm(n)*SHF_DATA(:,:,:,n,:)
enddo
case default
call exit_POP(sigAbort,'init_shf: Unknown value for shf_data_type')
end select
!-----------------------------------------------------------------------
!
! now check interpolation period (shf_interp_freq) to set the
! time for the next temporal interpolation (shf_interp_next).
!
! if no interpolation is to be done, set next interpolation time
! to a large number so the surface heat 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 heat flux
! update test in routine set_surface_forcing will always be true.
!
!-----------------------------------------------------------------------
select case (shf_interp_freq)
case ('never')
shf_interp_next = never
shf_interp_last = never
shf_interp_inc = c0
case ('n-hour')
call find_interp_time(shf_interp_inc, shf_interp_next)
case ('every-timestep')
shf_interp_next = always
shf_interp_inc = c0
case default
call exit_POP(sigAbort, &
'init_shf: Unknown value for shf_interp_freq')
end select
if (nsteps_total == 0) shf_interp_last = thour00
!-----------------------------------------------------------------------
!
! echo forcing options to stdout.
!
!-----------------------------------------------------------------------
shf_data_label = 'Surface Heat Flux'
call echo_forcing_options(shf_data_type, shf_formulation, &
shf_data_inc, shf_interp_freq, &
shf_interp_type, shf_interp_inc, &
shf_data_label)
!-----------------------------------------------------------------------
!EOC
call flushm (stdout)
end subroutine init_shf
!***********************************************************************
!BOP
! !IROUTINE: set_shf
! !INTERFACE:
subroutine set_shf(STF) 1,13
! !DESCRIPTION:
! Updates the current value of the surface heat flux array
! (shf) by interpolating to the current time or calculating
! fluxes based on states at current time. If new data are
! required for interpolation, new data are read.
!
! !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) :: &
iblock
!-----------------------------------------------------------------------
!
! check if new data is necessary for interpolation. if yes, then
! shuffle indices in SHF_DATA and shf_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 using updated shf data or compute fluxes
! based on current or interpolated state data.
!
!-----------------------------------------------------------------------
select case(shf_data_type)
case ('analytic')
select case (shf_formulation)
case ('restoring')
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock=1,nblocks_clinic
STF(:,:,1,iblock) = (SHF_DATA(:,:,iblock,shf_data_sst,1) - &
TRACER(:,:,1,1,curtime,iblock))* &
shf_restore_rtau*dz(1)
end do
!$OMP END PARALLEL DO
end select
case ('annual')
select case (shf_formulation)
case ('restoring')
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock=1,nblocks_clinic
STF(:,:,1,iblock) = (SHF_DATA(:,:,iblock,shf_data_sst,1) - &
TRACER(:,:,1,1,curtime,iblock))* &
shf_restore_rtau*dz(1)
end do
!$OMP END PARALLEL DO
case ('Barnier-restoring')
call calc_shf_barnier_restoring(STF,1)
case ('bulk-NCEP')
call calc_shf_bulk_ncep(STF,1)
case ('partially-coupled')
call calc_shf_partially_coupled(1)
end select
case ('monthly-equal','monthly-calendar')
shf_data_label = 'SHF Monthly'
if (thour00 >= shf_data_update) then
call update_forcing_data(shf_data_time, shf_data_time_min_loc,&
shf_interp_type, shf_data_next, &
shf_data_update, shf_data_type, &
shf_data_inc, SHF_DATA(:,:,:,:,1:12),&
shf_data_renorm, &
shf_data_label, shf_data_names, &
shf_bndy_loc, shf_bndy_type, &
shf_filename, shf_file_fmt)
endif
if (thour00 >= shf_interp_next .or. nsteps_run == 0) then
call interpolate_forcing(SHF_DATA(:,:,:,:,0), &
SHF_DATA(:,:,:,:,1:12), &
shf_data_time, shf_interp_type, &
shf_data_time_min_loc, shf_interp_freq, &
shf_interp_inc, shf_interp_next, &
shf_interp_last, nsteps_run)
if (nsteps_run /= 0) shf_interp_next = &
shf_interp_next + shf_interp_inc
endif
select case (shf_formulation)
case ('restoring')
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock=1,nblocks_clinic
STF(:,:,1,iblock) = (SHF_DATA(:,:,iblock,shf_data_sst,0) - &
TRACER(:,:,1,1,curtime,iblock))* &
shf_restore_rtau*dz(1)
end do
!$OMP END PARALLEL DO
case ('Barnier-restoring')
call calc_shf_barnier_restoring(STF,12)
case ('bulk-NCEP')
call calc_shf_bulk_ncep(STF,12)
case ('partially-coupled')
call calc_shf_partially_coupled(12)
end select
case('n-hour')
shf_data_label = 'SHF n-hour'
if (thour00 >= shf_data_update) then
call update_forcing_data(shf_data_time, shf_data_time_min_loc,&
shf_interp_type, shf_data_next, &
shf_data_update, shf_data_type, &
shf_data_inc, &
SHF_DATA(:,:,:,:,1:shf_interp_order),&
shf_data_renorm, &
shf_data_label, shf_data_names, &
shf_bndy_loc, shf_bndy_type, &
shf_filename, shf_file_fmt)
endif
if (thour00 >= shf_interp_next .or. nsteps_run == 0) then
call interpolate_forcing(SHF_DATA(:,:,:,:,0), &
SHF_DATA(:,:,:,:,1:shf_interp_order), &
shf_data_time, shf_interp_type, &
shf_data_time_min_loc, shf_interp_freq, &
shf_interp_inc, shf_interp_next, &
shf_interp_last, nsteps_run)
if (nsteps_run /= 0) shf_interp_next = &
shf_interp_next + shf_interp_inc
endif
select case (shf_formulation)
case ('restoring')
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock=1,nblocks_clinic
STF(:,:,1,iblock) = (SHF_DATA(:,:,iblock,shf_data_sst,0) - &
TRACER(:,:,1,1,curtime,iblock))* &
shf_restore_rtau*dz(1)
end do
!$OMP END PARALLEL DO
case ('Barnier-restoring')
call calc_shf_barnier_restoring(STF, shf_interp_order)
case ('partially-coupled')
call calc_shf_partially_coupled(shf_interp_order)
case ('bulk-NCEP')
call calc_shf_bulk_ncep(STF, shf_interp_order)
end select
end select ! shf_data_type
!-----------------------------------------------------------------------
!EOC
end subroutine set_shf
!***********************************************************************
!BOP
! !IROUTINE: calc_shf_barnier_restoring
! !INTERFACE:
subroutine calc_shf_barnier_restoring(STF, time_dim) 3
! !DESCRIPTION:
! calculates surface heat fluxes
!
! !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 heat flux at current timestep
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
nearest_data, now, &! indices for nearest,interpolated time slices
iblock ! local address of current block
real (r8) :: &
tcheck, ice_cutoff, ice_restore_temp
!-----------------------------------------------------------------------
!
! local parameters
!
!-----------------------------------------------------------------------
ice_cutoff = 0.9_r8
ice_restore_temp = -2.0_r8
!-----------------------------------------------------------------------
!
! if annual forcing, no interpolation to current time is necessary.
! otherwise, interpolated fields in index=0 slice of data array
!
!-----------------------------------------------------------------------
if (shf_data_type == 'annual') then
now = 1
nearest_data = 1
else
now = 0
!*** find nearest data time and use it for determining the ice
!*** mask in place of interpolated field.
!*** NOTE: this is for backward compatibility. perhaps
!*** interpolating and using a cut-off of .45 would be acceptable.
tcheck = (shf_data_update - thour00)/shf_data_inc
select case(shf_interp_type)
case ('nearest')
nearest_data = shf_data_time_min_loc
case ('linear')
if (tcheck > 0.5) then
nearest_data = shf_data_time_min_loc
else
nearest_data = shf_data_time_min_loc + 1
endif
case ('4point')
if (tcheck > 0.5) then
nearest_data = shf_data_time_min_loc + 1
else
nearest_data = shf_data_time_min_loc + 2
endif
end select
if ((nearest_data - time_dim) > 0 ) nearest_data = &
nearest_data - time_dim
endif
!-----------------------------------------------------------------------
!
! calculate forcing for each block
!
!-----------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock=1,nblocks_clinic
!-----------------------------------------------------------------------
!
! check for ice concentration >= ice_cutoff in the nearest month.
! if there is ice, set TAU to be constant and set TSTAR to
! ice_restore_temp.
!
!-----------------------------------------------------------------------
where (SHF_DATA(:,:,iblock,shf_data_ice,nearest_data) >= &
ice_cutoff)
SHF_DATA(:,:,iblock,shf_data_tau,now) = shf_restore_tau
SHF_DATA(:,:,iblock,shf_data_tstar,now) = ice_restore_temp
endwhere
!-----------------------------------------------------------------------
!
! apply restoring only where TAU is defined.
!
!-----------------------------------------------------------------------
where (SHF_DATA(:,:,iblock,shf_data_tau,now) > c0)
STF(:,:,1,iblock) =(SHF_DATA(:,:,iblock,shf_data_tstar,now) - &
TRACER(:,:,1,1,curtime,iblock))* &
dz(1)/SHF_DATA(:,:,iblock,shf_data_tau,now)
elsewhere
STF(:,:,1,iblock) = c0
end where
!-----------------------------------------------------------------------
!
! copy penetrative shortwave into its own array (SHF_QSW) and
! convert to T flux from W/m^2.
!
!-----------------------------------------------------------------------
SHF_QSW(:,:,iblock) = SHF_DATA(:,:,iblock,shf_data_qsw,now)* &
hflux_factor
SHF_QSW_RAW(:,:,iblock) = SHF_QSW(:,:,iblock)
end do
!$OMP END PARALLEL DO
!-----------------------------------------------------------------------
!EOC
end subroutine calc_shf_barnier_restoring
!***********************************************************************
!BOP
! !IROUTINE: calc_shf_bulk_ncep
! !INTERFACE:
subroutine calc_shf_bulk_ncep(STF, time_dim) 3,2
! !DESCRIPTION:
! Calculates surface heat flux from a combination of
! air-sea fluxes (based on air temperature, specific humidity,
! solar short wave flux, cloud fraction, and windspeed)
! and restoring terms (due to restoring fields of SST).
!
! Notes:
! the forcing data (on t-grid)
! are computed as SHF\_DATA(:,:,shf\_comp\_*,now) where:
!
! shf\_data\_sst, restoring SST (C)
! shf\_data\_tair, surface air temp. at tair\_height (K)
! shf\_data\_qair, specific humidity at qair\_height (kg/kg)
! shf\_data\_qsw, surface short wave flux ($W/m^2$)
! shf\_data\_cldfrac, cloud fraction (0.-1.)
! shf\_data\_windspd , windspeed at height windspd\_height (m/s)
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: &
time_dim
! !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) :: &
n, now, k, &
iblock ! local address of current block
real (r8), dimension(nx_block,ny_block) :: &
RTEA, &! work array
FRAC_CLOUD_COVER ! fractional cloud cover
real (r8), parameter :: &
windspd_height = 10.0_r8, &
tair_height = 2.0_r8, &
qair_height = 2.0_r8, &
qair_mod_fact = 0.94_r8, &! factor to modify humidity
sw_mod_fact = 0.875_r8, &! factor to modify short-wave flux
sw_mod_albedo = 0.93_r8 ! factor to modify albedo
!-----------------------------------------------------------------------
!
! shf_weak_restore= weak(non-ice) restoring heatflux per degree (W/m2/C)
! shf_strong_restore= strong (ice) .. .. .. .. .. ..
!
! to calculate restoring factors, use mixed layer of 50m,
! and restoring time constant tau (days):
!
! Q (W/m2/C)
! tau = 6 : 386.0
! tau = 30 : 77.2
! tau = 182.5: 12.0
! tau = 365 : 6.0
! tau = 730 : 3.0
! tau = Inf : 0.0
!
!---------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! set location of interpolated data
!
!-----------------------------------------------------------------------
if (shf_data_type == 'annual') then
now = 1
else
now = 0
endif
!----------------------------------------------------------------------
!
! compute ocean weights (fraction of ocean vs. ice) every timestep
!
!----------------------------------------------------------------------
call ocean_weights(now)
!----------------------------------------------------------------------
!
! do the rest of the computation for each block
!
!----------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblock,FRAC_CLOUD_COVER,RTEA)
do iblock=1,nblocks_clinic
!----------------------------------------------------------------------
!
! compute sensible and latent heat fluxes
!
!----------------------------------------------------------------------
call sen_lat_flux( &
SHF_DATA(:,:,iblock,shf_data_windspd,now), windspd_height, &
TRACER(:,:,1,1,curtime,iblock), &
SHF_DATA(:,:,iblock,shf_data_tair,now), tair_height, &
SHF_DATA(:,:,iblock,shf_data_qair,now), qair_height, &
T0_Kelvin, SHF_COMP(:,:,iblock,shf_comp_qsens), &
SHF_COMP(:,:,iblock,shf_comp_qlat))
!----------------------------------------------------------------------
!
! compute short wave and long wave fluxes
!
!----------------------------------------------------------------------
SHF_COMP(:,:,iblock,shf_comp_qsw) = sw_mod_albedo*sw_mod_fact* &
SHF_DATA(:,:,iblock,shf_data_qsw,now)
FRAC_CLOUD_COVER = c1 - CCINT(:,:,iblock)* &
SHF_DATA(:,:,iblock,shf_data_cldfrac,now)**2
RTEA = sqrt( c1000*SHF_DATA(:,:,iblock,shf_data_qair,now) &
/(0.622_r8 + 0.378_r8 &
*SHF_DATA(:,:,iblock,shf_data_qair,now)) + eps2 )
SHF_COMP(:,:,iblock,shf_comp_qlw) = -emissivity*stefan_boltzmann*&
SHF_DATA(:,:,iblock,shf_data_tair,now)**3* &
(SHF_DATA(:,:,iblock,shf_data_tair,now)* &
(0.39_r8-0.05_r8*RTEA)*FRAC_CLOUD_COVER + &
c4*(TRACER(:,:,1,1,curtime,iblock) + &
T0_Kelvin - &
SHF_DATA(:,:,iblock,shf_data_tair,now)) )
!----------------------------------------------------------------------
!
! weak temperature restoring term (note: OCN_WGT = 0 at land pts)
!
!----------------------------------------------------------------------
SHF_COMP(:,:,iblock,shf_comp_wrest) = shf_weak_restore* &
MASK_SR(:,:,iblock)*OCN_WGT(:,:,iblock)* &
(SHF_DATA(:,:,iblock,shf_data_sst,now) - &
TRACER(:,:,1,1,curtime,iblock))
!----------------------------------------------------------------------
!
! strong temperature restoring term
!
!----------------------------------------------------------------------
where (KMT(:,:,iblock) > 0)
SHF_COMP(:,:,iblock,shf_comp_srest) = shf_strong_restore* &
(c1-OCN_WGT(:,:,iblock))* &
(SHF_DATA(:,:,iblock,shf_data_sst,now) - &
TRACER(:,:,1,1,curtime,iblock))
endwhere
where (KMT(:,:,iblock) > 0 .and. MASK_SR(:,:,iblock) == 0)
SHF_COMP(:,:,iblock,shf_comp_srest) = shf_strong_restore_ms* &
(SHF_DATA(:,:,iblock,shf_data_sst,now) - &
TRACER(:,:,1,1,curtime,iblock))
endwhere
!----------------------------------------------------------------------
!
! net surface heat flux (W/m^2) (except penetrative shortwave flux)
! convert to model units
!
!----------------------------------------------------------------------
STF(:,:,1,iblock) = hflux_factor* &
(OCN_WGT(:,:,iblock)*MASK_SR(:,:,iblock)* &
(SHF_COMP(:,:,iblock,shf_comp_qsens) + &
SHF_COMP(:,:,iblock,shf_comp_qlat ) + &
SHF_COMP(:,:,iblock,shf_comp_qlw )) + &
SHF_COMP(:,:,iblock,shf_comp_wrest) + &
SHF_COMP(:,:,iblock,shf_comp_srest))
!----------------------------------------------------------------------
!
! copy penetrative shortwave flux into its own array (SHF_QSW) and
! convert it and SHF to model units.
!
!----------------------------------------------------------------------
SHF_QSW(:,:,iblock) = SHF_COMP(:,:,iblock,shf_comp_qsw)* &
OCN_WGT(:,:,iblock)*MASK_SR(:,:,iblock)* &
hflux_factor
SHF_QSW_RAW(:,:,iblock) = SHF_COMP(:,:,iblock,shf_comp_qsw)* &
hflux_factor
end do
!$OMP END PARALLEL DO
!----------------------------------------------------------------------
!EOC
end subroutine calc_shf_bulk_ncep
!***********************************************************************
!BOP
! !IROUTINE: calc_shf_partially_coupled
! !INTERFACE:
subroutine calc_shf_partially_coupled(time_dim) 3,1
! !DESCRIPTION:
! Calculates weak and strong restoring components of surface heat flux
! for partially-coupled formulation. These components will later be
! added to shf_comp_cpl component in set_coupled_forcing
! (forcing_coupled) to form the total surface heat flux.
!
! The only forcing dataset (on t-grid) is
! shf_data_sst, restoring SST
!
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: &
time_dim
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
n, now, k, &
iblock ! local address of current block
real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
WORK1 ! work array
!-----------------------------------------------------------------------
!
! set location of interpolated data
!
!-----------------------------------------------------------------------
if (shf_data_type == 'annual') then
now = 1
else
now = 0
endif
!----------------------------------------------------------------------
!
! compute ocean weights (fraction of ocean vs. ice) every timestep,
! if needed
!
!----------------------------------------------------------------------
if ( .not. luse_cpl_ifrac ) then
call ocean_weights (now)
WORK1 = OCN_WGT*MASK_SR
else
WORK1 = MASK_SR
endif
!----------------------------------------------------------------------
!
! do the rest of the computation for each block
!
!----------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock=1,nblocks_clinic
!----------------------------------------------------------------------
!
! weak temperature restoring term (note: MASK_SR = 0. at land and
! marginal sea points)
! note that weak restoring may be applied to every non-marginal-sea
! ocean point.
!
!----------------------------------------------------------------------
SHF_COMP(:,:,iblock,shf_comp_wrest) = shf_weak_restore* &
WORK1(:,:,iblock)* &
(SHF_DATA(:,:,iblock,shf_data_sst,now) - &
TRACER(:,:,1,1,curtime,iblock))
!----------------------------------------------------------------------
!
! strong temperature restoring term
! note that strong restoring may be applied only in marginal seas.
! in under-ice regions, the ice formation term may replace the
! strong-restoring term.
!
!----------------------------------------------------------------------
where (KMT(:,:,iblock) > 0)
SHF_COMP(:,:,iblock,shf_comp_srest) = shf_strong_restore* &
(c1-OCN_WGT(:,:,iblock))* &
(SHF_DATA(:,:,iblock,shf_data_sst,now) - &
TRACER(:,:,1,1,curtime,iblock))
endwhere
where (KMT(:,:,iblock) > 0 .and. MASK_SR(:,:,iblock) == 0)
SHF_COMP(:,:,iblock,shf_comp_srest) = shf_strong_restore_ms* &
(SHF_DATA(:,:,iblock,shf_data_sst,now) - &
TRACER(:,:,1,1,curtime,iblock))
endwhere
!----------------------------------------------------------------------
!
! convert to model units: (W/m^2) to (C*cm/s)
!
!----------------------------------------------------------------------
SHF_COMP(:,:,iblock,shf_comp_wrest) = &
SHF_COMP(:,:,iblock,shf_comp_wrest)*hflux_factor
SHF_COMP(:,:,iblock,shf_comp_srest) = &
SHF_COMP(:,:,iblock,shf_comp_srest)*hflux_factor
end do
!$OMP END PARALLEL DO
!----------------------------------------------------------------------
!EOC
end subroutine calc_shf_partially_coupled
!***********************************************************************
!BOP
! !IROUTINE: sen_lat_flux
! !INTERFACE:
subroutine sen_lat_flux(US,hu,SST,TH,ht,QH,hq,tk0,HS,HL) 1,3
! !DESCRIPTION:
! Computes latent and sensible heat fluxes following bulk formulae and
! coefficients in Large and Pond (1981; 1982)
!
! Assume 1) a neutral 10m drag coefficient = cdn =
! .0027/u10 + .000142 + .0000764 u10
! 2) a neutral 10m stanton number ctn= .0327 sqrt(cdn), unstable
! ctn= .0180 sqrt(cdn), stable
! 3) a neutral 10m dalton number cen= .0346 sqrt(cdn)
! 4) the saturation humidity of air at t(k) = qsat(t) ($kg/m^3$)
!
! note 1) here, tstar = <wt>/u*, and qstar = <wq>/u*.
! 2) wind speedx should all be above a minimum speed say 0.5 m/s
! 3) with optional interation loop, niter=3, should suffice
!
! *** this version is for analyses inputs with hu = 10m and ht = hq **
! *** also, SST enters in Celsius ***************************
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension (nx_block,ny_block), intent(in) :: &
US, &! mean wind speed (m/s) at height hu (m)
TH, &! mean air temperature (k) at height ht (m)
QH, &! mean air humidity (kg/kg) at height hq (m)
SST ! sea surface temperature (K)
real (r8), intent(in) :: &
hu, &! height (m) for mean wind speed
ht, &! height (m) for mean air temperature
hq, &! height (m) for mean air humidity
tk0 ! Celsius zero point
! !OUTPUT PARAMETERS:
real (r8), dimension (nx_block,ny_block), intent(out) :: &
HS, &! sensible heat flux (w/m^2), into ocean
HL ! latent heat flux (w/m^2), into ocean
!EOP
!BOC
!--------------------------------------------------------------------------
!
! local variables
!
!--------------------------------------------------------------------------
real (r8), dimension (nx_block,ny_block) :: &
SH,T0,DELP,DELQ,STABLETMP,RDN,RHN,USTARR,TSTARR,QSTARR,TAU, &
HUOL,HTOL,HQOL,SSHUM,PSIMH,PSIXH,RD,UZN,RH,RE,QSAT
real (r8) :: &
ren,umin,zolmin,vonk,lapse_rate,gravity_mks,f1,refhgt,aln,czol
!-----------------------------------------------------------------------
!
! constants
!
!-----------------------------------------------------------------------
umin = 0.5_r8 ! minimum wind speed
zolmin = -100._r8 ! minimum stability parameter
vonk = 0.4_r8 ! Von Karman''s constant
lapse_rate = 0.01_r8 ! abiabatic lapse rate deg/m
gravity_mks = grav/100.0_r8 ! gravity m/s/s
f1 = 0.606_r8
refhgt = 10.0_r8 ! reference height
aln = log(ht/refhgt)
czol = hu*vonk*gravity_mks
SH = max(US,umin)
!-----------------------------------------------------------------------
!
! initial guess z/l=0.0; hu=ht=hq=z
!
!-----------------------------------------------------------------------
T0 = TH * (c1 + f1 * QH) ! virtual temperature (k)
QSAT = 640380._r8 / exp(5107.4_r8/(SST+tk0))
SSHUM = 0.98_r8 * QSAT/rho_air
! sea surface humidity (kg/kg)
DELP = TH + lapse_rate*ht - SST - tk0 ! pot temperature diff (k)
DELQ = QH - SSHUM
STABLETMP = 0.5_r8 + sign(0.5_r8 , DELP)
RDN = sqrt(CDN(SH))
RHN = (c1-STABLETMP)* 0.0327_r8 + STABLETMP * 0.0180_r8
ren = 0.0346_r8
USTARR = RDN * SH
TSTARR = RHN * DELP
QSTARR = REN * DELQ
!-----------------------------------------------------------------------
!
! first iteration loop
!
!-----------------------------------------------------------------------
HUOL = czol * (TSTARR/T0 + QSTARR/(c1/f1+QH)) / USTARR**2
HUOL = max(HUOL,zolmin)
STABLETMP = 0.5_r8 + sign(0.5_r8 , HUOL)
HTOL = HUOL * ht / hu
HQOL = HUOL * hq / hu
!-----------------------------------------------------------------------
!
! evaluate all stability functions assuming hq = ht
!
!-----------------------------------------------------------------------
SSHUM = max(sqrt(abs(c1 - 16._r8*HUOL)),c1)
SSHUM = sqrt(SSHUM)
PSIMH = -5._r8 * HUOL * STABLETMP + (c1-STABLETMP) &
* log((c1+SSHUM*(c2+SSHUM))*(c1+SSHUM*SSHUM)/8._r8) &
- c2*atan(SSHUM)+1.571_r8
SSHUM = max(sqrt(abs(c1 - 16._r8*HTOL)),c1)
PSIXH = -5._r8*HTOL*STABLETMP + (c1-STABLETMP)*c2*log((c1+SSHUM)/c2)
!-----------------------------------------------------------------------
!
! shift wind speed using old coefficient
!
!-----------------------------------------------------------------------
RD = RDN / (c1-RDN/vonk*PSIMH )
UZN = max(SH * RD / RDN , umin)
!-----------------------------------------------------------------------
!
! update the transfer coefficients at 10 meters and neutral stability
!
!-----------------------------------------------------------------------
RDN = sqrt(CDN(UZN))
ren = 0.0346_r8
RHN = (c1-STABLETMP)*0.0327_r8 + STABLETMP *0.0180_r8
!-----------------------------------------------------------------------
!
! shift all coefficients to the measurement height and stability
!
!-----------------------------------------------------------------------
RD = RDN / (c1-RDN/vonk*PSIMH )
RH = RHN / (c1+RHN/vonk*( aln -PSIXH) )
RE = ren / (c1+ren/vonk*( aln -PSIXH) )
!-----------------------------------------------------------------------
!
! update USTARR, TSTARR, QSTARR using updated, shifted coefficients
!
!-----------------------------------------------------------------------
USTARR = RD * SH
QSTARR = RE * DELQ
TSTARR = RH * DELP
!-----------------------------------------------------------------------
!
! second iteration to converge on z/l and hence the fluxes
!
!-----------------------------------------------------------------------
HUOL= czol * (TSTARR/T0+QSTARR/(c1/f1+QH)) / USTARR**2
HUOL= max(HUOL,zolmin)
STABLETMP = 0.5_r8 + sign(0.5_r8 , HUOL)
HTOL = HUOL * ht / hu
HQOL = HUOL * hq / hu
!-----------------------------------------------------------------------
!
! evaluate all stability functions assuming hq = ht
!
!-----------------------------------------------------------------------
SSHUM = max(sqrt(abs(c1 - 16.*HUOL)),c1)
SSHUM = sqrt(SSHUM)
PSIMH = -5._r8 * HUOL * STABLETMP + (c1-STABLETMP) &
* log((c1+SSHUM*(c2+SSHUM))*(c1+SSHUM*SSHUM)/8._r8) &
- c2*atan(SSHUM)+1.571_r8
SSHUM = max(sqrt(abs(c1 - 16._r8*HTOL)),c1)
PSIXH = -5._r8*HTOL*STABLETMP + (c1-STABLETMP)*c2*log((c1+SSHUM)/c2)
!-----------------------------------------------------------------------
!
! shift wind speed using old coefficient
!
!-----------------------------------------------------------------------
RD = RDN / (c1-RDN/vonk*PSIMH )
UZN = max(SH * RD / RDN , umin)
!-----------------------------------------------------------------------
!
! update the transfer coefficients at 10 meters and neutral stability
!
!-----------------------------------------------------------------------
RDN = sqrt(CDN(UZN))
ren = 0.0346_r8
RHN = (c1-STABLETMP)*0.0327_r8 + STABLETMP*0.0180_r8
!-----------------------------------------------------------------------
!
! shift all coefficients to the measurement height and stability
!
!-----------------------------------------------------------------------
RD = RDN / (c1-RDN/vonk*PSIMH )
RH = RHN / (c1+RHN/vonk*( aln -PSIXH) )
RE = ren / (c1+ren/vonk*( aln -PSIXH) )
!-----------------------------------------------------------------------
!
! update USTARR, TSTARR, QSTARR using updated, shifted coefficients
!
!-----------------------------------------------------------------------
USTARR = RD * SH
QSTARR = RE * DELQ
TSTARR = RH * DELP
!-----------------------------------------------------------------------
!
! done >>>> compute the fluxes
!
!-----------------------------------------------------------------------
TAU = rho_air * USTARR**2
TAU = TAU * US / SH
HS = cp_air* TAU * TSTARR / USTARR
HL = latent_heat_vapor * TAU * QSTARR / USTARR
!-----------------------------------------------------------------------
!EOC
end subroutine sen_lat_flux
!***********************************************************************
!BOP
! !IROUTINE: CDN
! !INTERFACE:
function CDN(UMPS) 7
! !DESCRIPTION:
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(in) :: &
UMPS
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block) :: &
CDN
!EOP
!BOC
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
CDN = 0.0027_r8/UMPS + .000142_r8 + .0000764_r8*UMPS
!-----------------------------------------------------------------------
!EOC
end function CDN
!***********************************************************************
!BOP
! !IROUTINE: ocean_weights
! !INTERFACE:
subroutine ocean_weights(now) 2
! !DESCRIPTION:
! Compute ocean weights (fraction of ocean vs. ice) every timestep
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: &
now
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
iblock
!----------------------------------------------------------------------
!
! compute ocean weights (fraction of ocean vs. ice) every timestep
!
!----------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock=1,nblocks_clinic
where (SHF_DATA(:,:,iblock,shf_data_sst,now) <= &
T_strong_restore_limit)
OCN_WGT(:,:,iblock) = c0
elsewhere
OCN_WGT(:,:,iblock) =(SHF_DATA(:,:,iblock,shf_data_sst,now) - &
T_strong_restore_limit)/dT_restore_limit
endwhere
where (SHF_DATA(:,:,iblock,shf_data_sst,now) >= &
T_weak_restore_limit) OCN_WGT(:,:,iblock) = c1
!*** zero OCN_WGT at land pts
where (KMT(:,:,iblock) == 0) OCN_WGT(:,:,iblock) = c0
end do
!$OMP END PARALLEL DO
!-----------------------------------------------------------------------
!EOC
end subroutine ocean_weights
end module forcing_shf
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||