!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| module forcing_shf !BOP ! !MODULE: forcing_shf ! !DESCRIPTION: ! Contains routines and variables used for determining the surface ! heat flux. ! ! !REVISION HISTORY: ! SVN:$Id: forcing_shf.F90 24379 2010-08-13 19:54:51Z 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) ! !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) ! !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) ! !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) ! !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) ! !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) ! !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 = /u*, and qstar = /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_mks * TAU * QSTARR / USTARR !----------------------------------------------------------------------- !EOC end subroutine sen_lat_flux !*********************************************************************** !BOP ! !IROUTINE: CDN ! !INTERFACE: function CDN(UMPS) ! !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) ! !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 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||