!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| module forcing_sfwf !BOP ! !MODULE: forcing_sfwf ! !DESCRIPTION: ! Contains routines and variables used for determining the ! surface fresh water flux. ! ! !REVISION HISTORY: ! SVN:$Id: forcing_sfwf.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 global_reductions use forcing_tools use forcing_shf use ice use time_management use prognostic use exit_mod implicit none private save ! !PUBLIC MEMBER FUNCTIONS: public :: init_sfwf, & set_sfwf ! !PUBLIC DATA MEMBERS: real (r8), public, allocatable, dimension(:,:,:,:) :: & SFWF_COMP real (r8), public, allocatable, dimension(:,:,:,:,:) :: & TFW_COMP real (r8), public :: &! public for use in restart sfwf_interp_last ! time when last interpolation was done !*** water balance factors for bulk-NCEP forcing real (r8), public :: &! public for use in restart sum_precip, &! global precip for water balance precip_fact = c1, &! factor for adjusting precip for water balance ssh_initial ! initial ssh real (r8), dimension(km), public :: & sal_initial logical (log_kind), public :: & lfw_as_salt_flx ! treat fw flux as virtual salt flux ! even with var.thickness sfc layer logical (log_kind), public :: & lsend_precip_fact ! if T,send precip_fact to cpl for use in fw balance ! (partially-coupled option) !EOP !BOC !----------------------------------------------------------------------- ! ! module variables ! !----------------------------------------------------------------------- real (r8), allocatable, dimension(:,:,:,:,:) :: & SFWF_DATA ! forcing data used to get SFWF real (r8), dimension(12) :: & sfwf_data_time ! time (hours) corresponding to surface fresh ! water fluxes real (r8), dimension(20) :: & sfwf_data_renorm ! factors for converting to model units real (r8) :: & sfwf_data_inc, &! time increment between values of forcing data sfwf_data_next, &! time to be used for next value of forcing data sfwf_data_update, &! time new forcing data needs to be added to interpolation set sfwf_interp_inc, &! time increment between interpolation sfwf_interp_next, &! time when next interpolation will be done sfwf_restore_tau, &! restoring time scale sfwf_restore_rtau, &! reciprocal of restoring time scale sfwf_weak_restore, &! sfwf_strong_restore, &! sfwf_strong_restore_ms ! integer (int_kind) :: & sfwf_interp_order, &! order of temporal interpolation sfwf_data_time_min_loc, &! time index for first SFWF_DATA point sfwf_data_num_fields integer (int_kind), public :: & sfwf_num_comps character (char_len), dimension(:), allocatable :: & sfwf_data_names ! short names for input data fields integer (int_kind), dimension(:), allocatable :: & sfwf_bndy_loc, &! location and field types for ghost sfwf_bndy_type ! cell update routines !*** integer addresses for various forcing data fields integer (int_kind) :: & ! restoring and partially-coupled options sfwf_data_sss integer (int_kind), public :: &! bulk-NCEP and partially-coupled (some) options sfwf_data_precip, & sfwf_comp_precip, & sfwf_comp_evap, & sfwf_comp_wrest, & sfwf_comp_srest real (r8) :: & ann_avg_precip, &! !sum_fw, &! !ann_avg_fw, &! ssh_final real (r8), dimension (km) :: & sal_final logical (log_kind) :: & ladjust_precip integer (int_kind),public :: &! used with the partially-coupled option sfwf_comp_cpl, & sfwf_data_flxio, & sfwf_comp_flxio, & tfw_num_comps, & tfw_comp_cpl, & tfw_comp_flxio real (r8), parameter :: & precip_mean = 3.4e-5_r8 character (char_len) :: & sfwf_filename, &! name of file conainting forcing data sfwf_file_fmt, &! format (bin or netcdf) of forcing file sfwf_interp_freq, &! keyword for period of temporal interpolation sfwf_interp_type, &! sfwf_data_label character (char_len),public :: & sfwf_data_type, &! keyword for period of forcing data sfwf_formulation logical (log_kind), public :: & lms_balance ! control balancing of P,E,M,R,S in marginal seas ! .T. only with sfc_layer_oldfree option !EOC !*********************************************************************** contains !*********************************************************************** !BOP ! !IROUTINE: init_sfwf ! !INTERFACE: subroutine init_sfwf(STF) ! !DESCRIPTION: ! Initializes surface fresh water flux forcing by either calculating ! or reading in the surface fresh water flux. Also does initial ! book-keeping concerning when new data is needed for the temporal ! interpolation and when the forcing will need to be updated. ! ! !REVISION HISTORY: ! same as module ! !INPUT/OUTPUT PARAMETERS: real (r8), dimension(nx_block,ny_block,nt,max_blocks_clinic), & intent(inout) :: & STF ! surface tracer fluxes at current timestep !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer(int_kind) :: & k, n, &! dummy loop indices iblock, &! block loop index nml_error ! namelist error flag character (char_len) :: & forcing_filename ! full filename for forcing input real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: & WORK ! temporary work space real (r8), dimension(:,:,:,:,:), allocatable :: & TEMP_DATA ! temporary array for reading monthly data type (block) :: & this_block ! block info for local block type (datafile) :: & forcing_file ! data file structure for input forcing file type (io_field_desc) :: & io_sss, &! io field descriptor for input sss field io_precip, &! io field descriptor for input precip field io_flxio ! io field descriptor for input io_flxio field type (io_dim) :: & i_dim, j_dim, &! dimension descriptors for horiz dimensions month_dim ! dimension descriptor for monthly data namelist /forcing_sfwf_nml/ sfwf_data_type, sfwf_data_inc, & sfwf_interp_type, sfwf_interp_freq, & sfwf_interp_inc, sfwf_restore_tau, & sfwf_filename, sfwf_file_fmt, & sfwf_data_renorm, sfwf_formulation, & ladjust_precip, sfwf_weak_restore,& sfwf_strong_restore, lfw_as_salt_flx, & sfwf_strong_restore_ms, & lsend_precip_fact, lms_balance !----------------------------------------------------------------------- ! ! read surface fresh water flux namelist input after setting ! default values. ! !----------------------------------------------------------------------- sfwf_formulation = 'restoring' sfwf_data_type = 'analytic' sfwf_data_inc = 1.e20_r8 sfwf_interp_type = 'nearest' sfwf_interp_freq = 'never' sfwf_interp_inc = 1.e20_r8 sfwf_restore_tau = 1.e20_r8 sfwf_filename = 'unknown-sfwf' sfwf_file_fmt = 'bin' sfwf_data_renorm = c1 !sfwf_data_renorm = 1.e-3_r8 ! convert from psu to msu ladjust_precip = .false. lms_balance = .false. sfwf_weak_restore = 0.092_r8 sfwf_strong_restore_ms = 0.6648_r8 sfwf_strong_restore = 0.6648_r8 lfw_as_salt_flx = .false. lsend_precip_fact = .false. if (my_task == master_task) then open (nml_in, file=nml_filename, status='old',iostat=nml_error) if (nml_error /= 0) then nml_error = -1 else nml_error = 1 endif do while (nml_error > 0) read(nml_in, nml=forcing_sfwf_nml,iostat=nml_error) end do if (nml_error == 0) close(nml_in) endif call broadcast_scalar(nml_error, master_task) if (nml_error /= 0) then call exit_POP(sigAbort,'ERROR reading forcing_sfwf_nml') endif call broadcast_scalar(sfwf_data_type, master_task) call broadcast_scalar(sfwf_data_inc, master_task) call broadcast_scalar(sfwf_interp_type, master_task) call broadcast_scalar(sfwf_interp_freq, master_task) call broadcast_scalar(sfwf_interp_inc, master_task) call broadcast_scalar(sfwf_restore_tau, master_task) call broadcast_scalar(sfwf_filename, master_task) call broadcast_scalar(sfwf_file_fmt, master_task) call broadcast_scalar(sfwf_formulation, master_task) call broadcast_array (sfwf_data_renorm, master_task) call broadcast_scalar(ladjust_precip, master_task) call broadcast_scalar(sfwf_weak_restore, master_task) call broadcast_scalar(sfwf_strong_restore, master_task) call broadcast_scalar(sfwf_strong_restore_ms, master_task) call broadcast_scalar(lfw_as_salt_flx, master_task) call broadcast_scalar(lsend_precip_fact, master_task) call broadcast_scalar(lms_balance, master_task) !----------------------------------------------------------------------- ! ! convert data_type to 'monthly-calendar' if input is 'monthly' ! !----------------------------------------------------------------------- if (sfwf_data_type == 'monthly') sfwf_data_type = 'monthly-calendar' !----------------------------------------------------------------------- ! ! set values based on sfwf_formulation ! !----------------------------------------------------------------------- select case (sfwf_formulation) case ('restoring') allocate(sfwf_data_names(1), & sfwf_bndy_loc (1), & sfwf_bndy_type (1)) sfwf_data_num_fields = 1 sfwf_data_sss = 1 sfwf_data_names(sfwf_data_sss) = 'SSS' sfwf_bndy_loc (sfwf_data_sss) = field_loc_center sfwf_bndy_type (sfwf_data_sss) = field_type_scalar case ('bulk-NCEP') sfwf_data_num_fields = 2 sfwf_data_sss = 1 sfwf_data_precip = 2 allocate(sfwf_data_names(sfwf_data_num_fields), & sfwf_bndy_loc (sfwf_data_num_fields), & sfwf_bndy_type (sfwf_data_num_fields)) sfwf_data_names(sfwf_data_sss) = 'SSS' sfwf_bndy_loc (sfwf_data_sss) = field_loc_center sfwf_bndy_type (sfwf_data_sss) = field_type_scalar sfwf_data_names(sfwf_data_precip) = 'PRECIPITATION' sfwf_bndy_loc (sfwf_data_precip) = field_loc_center sfwf_bndy_type (sfwf_data_precip) = field_type_scalar sfwf_num_comps = 4 sfwf_comp_precip = 1 sfwf_comp_evap = 2 sfwf_comp_wrest = 3 sfwf_comp_srest = 4 case ('partially-coupled') sfwf_data_num_fields = 2 sfwf_data_sss = 1 sfwf_data_flxio = 2 allocate(sfwf_data_names(sfwf_data_num_fields), & sfwf_bndy_loc (sfwf_data_num_fields), & sfwf_bndy_type (sfwf_data_num_fields)) sfwf_data_names(sfwf_data_sss) = 'SSS' sfwf_bndy_loc (sfwf_data_sss) = field_loc_center sfwf_bndy_type (sfwf_data_sss) = field_type_scalar sfwf_data_names(sfwf_data_flxio) = 'FLXIO' sfwf_bndy_loc (sfwf_data_flxio) = field_loc_center sfwf_bndy_type (sfwf_data_flxio) = field_type_scalar sfwf_num_comps = 4 sfwf_comp_wrest = 1 sfwf_comp_srest = 2 sfwf_comp_cpl = 3 sfwf_comp_flxio = 4 tfw_num_comps = 2 tfw_comp_cpl = 1 tfw_comp_flxio = 2 case default call exit_POP(sigAbort, & 'init_sfwf: Unknown value for sfwf_formulation') end select if ( sfwf_formulation == 'bulk-NCEP' .or. & sfwf_formulation == 'partially-coupled' ) then !*** calculate initial salinity profile for ocean points that are !*** not marginal seas. !*** very first step of run if (ladjust_precip .and. nsteps_total == 0) then sum_precip = c0 ssh_initial = c0 !sum_fw = c0 do k = 1,km !$OMP PARALLEL DO PRIVATE(iblock, this_block) do iblock=1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) if (partial_bottom_cells) then WORK(:,:,iblock) = & merge(TRACER(:,:,k,2,curtime,iblock)* & TAREA(:,:,iblock)*DZT(:,:,k,iblock), & c0, k <= KMT(:,:,iblock) .and. & MASK_SR(:,:,iblock) > 0) else WORK(:,:,iblock) = & merge(TRACER(:,:,k,2,curtime,iblock)* & TAREA(:,:,iblock)*dz(k), & c0, k <= KMT(:,:,iblock) .and. & MASK_SR(:,:,iblock) > 0) endif end do !$OMP END PARALLEL DO sal_initial(k) = global_sum(WORK,distrb_clinic,field_loc_center)/ & (volume_t_k(k) - volume_t_marg_k(k)) ! in m^3 enddo endif endif !----------------------------------------------------------------------- ! ! calculate inverse of restoring time scale and convert to seconds. ! !----------------------------------------------------------------------- sfwf_restore_rtau = c1/(seconds_in_day*sfwf_restore_tau) !----------------------------------------------------------------------- ! ! convert interp_type to corresponding integer value. ! !----------------------------------------------------------------------- select case (sfwf_interp_type) case ('nearest') sfwf_interp_order = 1 case ('linear') sfwf_interp_order = 2 case ('4point') sfwf_interp_order = 4 case default call exit_POP(sigAbort, & 'init_sfwf: Unknown value for sfwf_interp_type') end select !----------------------------------------------------------------------- ! ! set values of the surface fresh water flux arrays (SFWF or ! SFWF_DATA) depending on type of the surface fresh water flux ! data. ! !----------------------------------------------------------------------- select case (sfwf_data_type) !----------------------------------------------------------------------- ! ! no surface fresh water flux, therefore no interpolation in time ! is needed (sfwf_interp_freq = 'none'), nor are there any new ! values to be used (sfwf_data_next = sfwf_data_update = never). ! !----------------------------------------------------------------------- case ('none') STF(:,:,2,:) = c0 sfwf_data_next = never sfwf_data_update = never sfwf_interp_freq = 'never' !----------------------------------------------------------------------- ! ! simple analytic surface salinity that is constant in time, ! therefore no new values will be needed. ! !----------------------------------------------------------------------- case ('analytic') allocate(SFWF_DATA(nx_block,ny_block,max_blocks_clinic, & sfwf_data_num_fields,1)) SFWF_DATA = c0 select case (sfwf_formulation) case ('restoring') SFWF_DATA(:,:,:,sfwf_data_sss,1) = 0.035_r8 end select sfwf_data_next = never sfwf_data_update = never sfwf_interp_freq = 'never' !----------------------------------------------------------------------- ! ! annual mean climatological surface salinity (read in from a file) ! that is constant in time, therefore no new values will be needed. ! !----------------------------------------------------------------------- case ('annual') allocate(SFWF_DATA(nx_block,ny_block,max_blocks_clinic, & sfwf_data_num_fields,1)) SFWF_DATA = c0 forcing_file = construct_file(sfwf_file_fmt, & full_name=trim(sfwf_filename), & record_length = rec_type_dbl, & recl_words=nx_global*ny_global) call data_set(forcing_file,'open_read') i_dim = construct_io_dim('i',nx_global) j_dim = construct_io_dim('j',nx_global) select case (sfwf_formulation) case ('restoring') io_sss = construct_io_field( & trim(sfwf_data_names(sfwf_data_sss)), & dim1=i_dim, dim2=j_dim, & field_loc = sfwf_bndy_loc(sfwf_data_sss), & field_type = sfwf_bndy_type(sfwf_data_sss), & d2d_array=SFWF_DATA(:,:,:,sfwf_data_sss,1)) call data_set(forcing_file,'define',io_sss) call data_set(forcing_file,'read' ,io_sss) call destroy_io_field(io_sss) case ('bulk-NCEP') io_sss = construct_io_field( & trim(sfwf_data_names(sfwf_data_sss)), & dim1=i_dim, dim2=j_dim, & field_loc = sfwf_bndy_loc(sfwf_data_sss), & field_type = sfwf_bndy_type(sfwf_data_sss), & d2d_array=SFWF_DATA(:,:,:,sfwf_data_sss,1)) io_precip = construct_io_field( & trim(sfwf_data_names(sfwf_data_precip)), & dim1=i_dim, dim2=j_dim, & field_loc = sfwf_bndy_loc(sfwf_data_precip), & field_type = sfwf_bndy_type(sfwf_data_precip), & d2d_array=SFWF_DATA(:,:,:,sfwf_data_precip,1)) call data_set(forcing_file,'define',io_sss) call data_set(forcing_file,'define',io_precip) call data_set(forcing_file,'read' ,io_sss) call data_set(forcing_file,'read' ,io_precip) call destroy_io_field(io_sss) call destroy_io_field(io_precip) allocate( SFWF_COMP(nx_block,ny_block,max_blocks_clinic, & sfwf_num_comps)) SFWF_COMP = c0 ! initialize case ('partially-coupled') io_sss = construct_io_field( & trim(sfwf_data_names(sfwf_data_sss)), & dim1=i_dim, dim2=j_dim, & field_loc = sfwf_bndy_loc(sfwf_data_sss), & field_type = sfwf_bndy_type(sfwf_data_sss), & d2d_array=SFWF_DATA(:,:,:,sfwf_data_sss,1)) io_flxio = construct_io_field( & trim(sfwf_data_names(sfwf_data_flxio)), & dim1=i_dim, dim2=j_dim, & field_loc = sfwf_bndy_loc(sfwf_data_flxio), & field_type = sfwf_bndy_type(sfwf_data_flxio), & d2d_array=SFWF_DATA(:,:,:,sfwf_data_flxio,1)) call data_set(forcing_file,'define',io_sss) call data_set(forcing_file,'define',io_flxio) call data_set(forcing_file,'read' ,io_sss) call data_set(forcing_file,'read' ,io_flxio) call destroy_io_field(io_sss) call destroy_io_field(io_flxio) allocate( SFWF_COMP(nx_block,ny_block,max_blocks_clinic, & sfwf_num_comps)) SFWF_COMP = c0 ! initialize end select call data_set(forcing_file,'close') call destroy_file(forcing_file) !*** renormalize values if necessary to compensate for !*** different units. do n = 1,sfwf_data_num_fields if (sfwf_data_renorm(n) /= c1) SFWF_DATA(:,:,:,n,:) = & sfwf_data_renorm(n)*SFWF_DATA(:,:,:,n,:) enddo sfwf_data_next = never sfwf_data_update = never sfwf_interp_freq = 'never' if (my_task == master_task) then write(stdout,blank_fmt) write(stdout,'(a25,a)') ' SFWF Annual file read: ', & trim(sfwf_filename) endif !----------------------------------------------------------------------- ! ! monthly mean climatological surface fresh water flux. all ! 12 months are read in from a file. interpolation order ! (sfwf_interp_order) may be specified with namelist input. ! !----------------------------------------------------------------------- case ('monthly-equal','monthly-calendar') allocate(SFWF_DATA(nx_block,ny_block,max_blocks_clinic, & sfwf_data_num_fields,0:12), & TEMP_DATA(nx_block,ny_block,12,max_blocks_clinic, & sfwf_data_num_fields) ) SFWF_DATA = c0 call find_forcing_times(sfwf_data_time, sfwf_data_inc, & sfwf_interp_type, sfwf_data_next, & sfwf_data_time_min_loc, & sfwf_data_update, sfwf_data_type) forcing_file = construct_file(sfwf_file_fmt, & full_name=trim(sfwf_filename), & record_length = rec_type_dbl, & recl_words=nx_global*ny_global) call data_set(forcing_file,'open_read') i_dim = construct_io_dim('i',nx_global) j_dim = construct_io_dim('j',nx_global) month_dim = construct_io_dim('month',12) select case (sfwf_formulation) case ('restoring') io_sss = construct_io_field( & trim(sfwf_data_names(sfwf_data_sss)), & dim1=i_dim, dim2=j_dim, dim3=month_dim, & field_loc = sfwf_bndy_loc(sfwf_data_sss), & field_type = sfwf_bndy_type(sfwf_data_sss), & d3d_array=TEMP_DATA(:,:,:,:,sfwf_data_sss)) call data_set(forcing_file,'define',io_sss) call data_set(forcing_file,'read' ,io_sss) call destroy_io_field(io_sss) case ('bulk-NCEP') io_sss = construct_io_field( & trim(sfwf_data_names(sfwf_data_sss)), & dim1=i_dim, dim2=j_dim, dim3=month_dim, & field_loc = sfwf_bndy_loc(sfwf_data_sss), & field_type = sfwf_bndy_type(sfwf_data_sss), & d3d_array=TEMP_DATA(:,:,:,:,sfwf_data_sss)) io_precip = construct_io_field( & trim(sfwf_data_names(sfwf_data_precip)), & dim1=i_dim, dim2=j_dim, dim3=month_dim, & field_loc = sfwf_bndy_loc(sfwf_data_precip), & field_type = sfwf_bndy_type(sfwf_data_precip), & d3d_array=TEMP_DATA(:,:,:,:,sfwf_data_precip)) call data_set(forcing_file,'define',io_sss) call data_set(forcing_file,'define',io_precip) call data_set(forcing_file,'read' ,io_sss) call data_set(forcing_file,'read' ,io_precip) call destroy_io_field(io_sss) call destroy_io_field(io_precip) allocate(SFWF_COMP(nx_block,ny_block,max_blocks_clinic, & sfwf_num_comps)) SFWF_COMP = c0 ! initialize case ('partially-coupled') io_sss = construct_io_field( & trim(sfwf_data_names(sfwf_data_sss)), & dim1=i_dim, dim2=j_dim, dim3=month_dim, & field_loc = sfwf_bndy_loc(sfwf_data_sss), & field_type = sfwf_bndy_type(sfwf_data_sss), & d3d_array=TEMP_DATA(:,:,:,:,sfwf_data_sss)) io_flxio = construct_io_field( & trim(sfwf_data_names(sfwf_data_flxio )), & dim1=i_dim, dim2=j_dim, dim3=month_dim, & field_loc = sfwf_bndy_loc(sfwf_data_flxio ), & field_type = sfwf_bndy_type(sfwf_data_flxio ), & d3d_array=TEMP_DATA(:,:,:,:,sfwf_data_flxio )) call data_set(forcing_file,'define',io_sss) call data_set(forcing_file,'define',io_flxio ) call data_set(forcing_file,'read' ,io_sss) call data_set(forcing_file,'read' ,io_flxio ) call destroy_io_field(io_sss) call destroy_io_field(io_flxio ) allocate(SFWF_COMP(nx_block,ny_block,max_blocks_clinic, & sfwf_num_comps)) SFWF_COMP = c0 ! initialize end select call data_set(forcing_file,'close') call destroy_file(forcing_file) !*** re-order data and renormalize values if necessary to !*** compensate for different units. !$OMP PARALLEL DO PRIVATE(iblock, k, n) do iblock=1,nblocks_clinic do k=1,sfwf_data_num_fields if (sfwf_data_renorm(k) /= c1) then do n=1,12 SFWF_DATA(:,:,iblock,k,n) = & TEMP_DATA(:,:,n,iblock,k)*sfwf_data_renorm(k) end do else do n=1,12 SFWF_DATA(:,:,iblock,k,n) = TEMP_DATA(:,:,n,iblock,k) end do endif end do end do !$OMP END PARALLEL DO deallocate(TEMP_DATA) if (my_task == master_task) then write(stdout,blank_fmt) write(stdout,'(a25,a)') ' SFWF Monthly file read: ', & trim(sfwf_filename) endif !----------------------------------------------------------------------- ! ! surface salinity specified every n-hours, where the n-hour ! increment should be specified with namelist input ! (sfwf_data_inc). only as many times as are necessary based on ! the order of the temporal interpolation scheme ! (sfwf_interp_order) reside in memory at any given time. ! !----------------------------------------------------------------------- case ('n-hour') allocate( SFWF_DATA(nx_block,ny_block,max_blocks_clinic, & sfwf_data_num_fields,0:sfwf_interp_order)) SFWF_DATA = c0 call find_forcing_times(sfwf_data_time, sfwf_data_inc, & sfwf_interp_type, sfwf_data_next, & sfwf_data_time_min_loc, & sfwf_data_update, sfwf_data_type) do n=1,sfwf_interp_order call get_forcing_filename(forcing_filename, sfwf_filename, & sfwf_data_time(n), sfwf_data_inc) forcing_file = construct_file(sfwf_file_fmt, & full_name=trim(sfwf_filename), & record_length = rec_type_dbl, & recl_words=nx_global*ny_global) call data_set(forcing_file,'open_read') i_dim = construct_io_dim('i',nx_global) j_dim = construct_io_dim('j',nx_global) select case (sfwf_formulation) case ('restoring') io_sss = construct_io_field( & trim(sfwf_data_names(sfwf_data_sss)), & dim1=i_dim, dim2=j_dim, & field_loc = sfwf_bndy_loc(sfwf_data_sss), & field_type = sfwf_bndy_type(sfwf_data_sss), & d2d_array=SFWF_DATA(:,:,:,sfwf_data_sss,n)) call data_set(forcing_file,'define',io_sss) call data_set(forcing_file,'read' ,io_sss) call destroy_io_field(io_sss) case ('bulk-NCEP') io_sss = construct_io_field( & trim(sfwf_data_names(sfwf_data_sss)), & dim1=i_dim, dim2=j_dim, & field_loc = sfwf_bndy_loc(sfwf_data_sss), & field_type = sfwf_bndy_type(sfwf_data_sss), & d2d_array=SFWF_DATA(:,:,:,sfwf_data_sss,n)) io_precip = construct_io_field( & trim(sfwf_data_names(sfwf_data_precip)), & dim1=i_dim, dim2=j_dim, & field_loc = sfwf_bndy_loc(sfwf_data_precip), & field_type = sfwf_bndy_type(sfwf_data_precip), & d2d_array=SFWF_DATA(:,:,:,sfwf_data_precip,n)) call data_set(forcing_file,'define',io_sss) call data_set(forcing_file,'define',io_precip) call data_set(forcing_file,'read' ,io_sss) call data_set(forcing_file,'read' ,io_precip) call destroy_io_field(io_sss) call destroy_io_field(io_precip) case ('partially-coupled') io_sss = construct_io_field( & trim(sfwf_data_names(sfwf_data_sss)), & dim1=i_dim, dim2=j_dim, & field_loc = sfwf_bndy_loc(sfwf_data_sss), & field_type = sfwf_bndy_type(sfwf_data_sss), & d2d_array=SFWF_DATA(:,:,:,sfwf_data_sss,n)) io_flxio = construct_io_field( & trim(sfwf_data_names(sfwf_data_flxio )), & dim1=i_dim, dim2=j_dim, & field_loc = sfwf_bndy_loc(sfwf_data_flxio ), & field_type = sfwf_bndy_type(sfwf_data_flxio ), & d2d_array=SFWF_DATA(:,:,:,sfwf_data_flxio ,n)) call data_set(forcing_file,'define',io_sss) call data_set(forcing_file,'define',io_flxio ) call data_set(forcing_file,'read' ,io_sss) call data_set(forcing_file,'read' ,io_flxio ) call destroy_io_field(io_sss) call destroy_io_field(io_flxio ) end select call data_set(forcing_file,'close') call destroy_file(forcing_file) if (my_task == master_task) then write(stdout,blank_fmt) write(stdout,'(a24,a)') ' SFWF n-hour file read: ', & trim(forcing_filename) endif enddo if (sfwf_formulation == 'bulk-NCEP' .or. & sfwf_formulation == 'partially-coupled') then allocate(SFWF_COMP(nx_block,ny_block,max_blocks_clinic, & sfwf_num_comps)) SFWF_COMP = c0 ! initialize endif !*** renormalize values if necessary to compensate for different !*** units. do n = 1,sfwf_data_num_fields if (sfwf_data_renorm(n) /= c1) SFWF_DATA(:,:,:,n,:) = & sfwf_data_renorm(n)*SFWF_DATA(:,:,:,n,:) enddo case default call exit_POP(sigAbort, & 'init_sfwf: Unknown value for sfwf_data_type') end select if ( sfwf_formulation == 'partially-coupled' ) then allocate ( TFW_COMP(nx_block,ny_block,nt,max_blocks_clinic, & tfw_num_comps)) TFW_COMP = c0 endif !----------------------------------------------------------------------- ! ! now check interpolation period (sfwf_interp_freq) to set the ! time for the next temporal interpolation (sfwf_interp_next). ! ! if no interpolation is to be done, set next interpolation time ! to a large number so the surface fresh water flux update test ! in routine set_surface_forcing will always be false. ! ! if interpolation is to be done every n-hours, find the first ! interpolation time greater than the current time. ! ! if interpolation is to be done every timestep, set next interpolation ! time to a large negative number so the surface fresh water flux ! update test in routine set_surface_forcing will always be true. ! !----------------------------------------------------------------------- select case (sfwf_interp_freq) case ('never') sfwf_interp_next = never sfwf_interp_last = never sfwf_interp_inc = c0 case ('n-hour') call find_interp_time(sfwf_interp_inc, sfwf_interp_next) case ('every-timestep') sfwf_interp_next = always sfwf_interp_inc = c0 case default call exit_POP(sigAbort, & 'init_sfwf: Unknown value for sfwf_interp_freq') end select if(nsteps_total == 0) sfwf_interp_last = thour00 !----------------------------------------------------------------------- ! ! echo forcing options to stdout. ! !----------------------------------------------------------------------- sfwf_data_label = 'Surface Fresh Water Flux' call echo_forcing_options(sfwf_data_type, & sfwf_formulation, sfwf_data_inc, & sfwf_interp_freq, sfwf_interp_type, & sfwf_interp_inc, sfwf_data_label) if (my_task == master_task) then if (lfw_as_salt_flx .and. sfc_layer_type == sfc_layer_varthick) & write(stdout,'(a47)') & ' Fresh water flux input as virtual salt flux' endif !----------------------------------------------------------------------- !EOC call flushm (stdout) end subroutine init_sfwf !*********************************************************************** !BOP ! !IROUTINE: set_sfwf ! !INTERFACE: subroutine set_sfwf(STF,FW,TFW) ! !DESCRIPTION: ! Updates the current value of the surface fresh water flux arrays ! by interpolating fields or computing fields at the current time. ! If new data are necessary for the interpolation, the new data are ! read from a file. ! ! !REVISION HISTORY: ! same as module ! !INPUT/OUTPUT PARAMETERS: real (r8), dimension(nx_block,ny_block,nt,max_blocks_clinic), & intent(inout) :: & STF, &! surface tracer fluxes at current timestep TFW ! tracer concentration in fresh water flux real (r8), dimension(nx_block,ny_block,max_blocks_clinic), & intent(inout) :: & FW ! fresh water flux if using varthick sfc layer !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (int_kind) :: & iblock ! local address for current block type (block) :: & this_block ! block info for current block !----------------------------------------------------------------------- ! ! check if new data is necessary for interpolation. if yes, then ! shuffle indices in SFWF_DATA and sfwf_data_time arrays ! and read in new data if necessary ('n-hour' case). note ! that no new data is necessary for 'analytic' and 'annual' cases. ! then perform interpolation or computation of fluxes at current time ! using updated forcing data. ! !----------------------------------------------------------------------- select case(sfwf_data_type) case ('analytic') select case (sfwf_formulation) case ('restoring') !$OMP PARALLEL DO PRIVATE(iblock, this_block) do iblock=1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) STF(:,:,2,iblock) = & (SFWF_DATA(:,:,iblock,sfwf_data_sss,1) - & TRACER(:,:,1,2,curtime,iblock))* & sfwf_restore_rtau*dz(1) end do !$OMP END PARALLEL DO end select case('annual') select case (sfwf_formulation) case ('restoring') !$OMP PARALLEL DO PRIVATE(iblock, this_block) do iblock=1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) STF(:,:,2,iblock) = & (SFWF_DATA(:,:,iblock,sfwf_data_sss,1) - & TRACER(:,:,1,2,curtime,iblock))* & sfwf_restore_rtau*dz(1) end do !$OMP END PARALLEL DO case ('bulk-NCEP') call calc_sfwf_bulk_ncep(STF,FW,TFW,1) case ('partially-coupled') call calc_sfwf_partially_coupled(1) end select case ('monthly-equal','monthly-calendar') sfwf_data_label = 'SFWF Monthly' if (thour00 >= sfwf_data_update) then call update_forcing_data( sfwf_data_time, & sfwf_data_time_min_loc, sfwf_interp_type, & sfwf_data_next, sfwf_data_update, & sfwf_data_type, sfwf_data_inc, & SFWF_DATA(:,:,:,:,1:12),sfwf_data_renorm, & sfwf_data_label, sfwf_data_names, & sfwf_bndy_loc, sfwf_bndy_type, & sfwf_filename, sfwf_file_fmt) endif if (thour00 >= sfwf_interp_next .or. nsteps_run == 0) then call interpolate_forcing(SFWF_DATA(:,:,:,:,0), & SFWF_DATA(:,:,:,:,1:12), & sfwf_data_time, sfwf_interp_type, & sfwf_data_time_min_loc, sfwf_interp_freq, & sfwf_interp_inc, sfwf_interp_next, & sfwf_interp_last, nsteps_run) if (nsteps_run /= 0) sfwf_interp_next = & sfwf_interp_next + sfwf_interp_inc endif select case (sfwf_formulation) case ('restoring') !$OMP PARALLEL DO PRIVATE(iblock, this_block) do iblock=1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) STF(:,:,2,iblock) = & (SFWF_DATA(:,:,iblock,sfwf_data_sss,0) - & TRACER(:,:,1,2,curtime,iblock))* & sfwf_restore_rtau*dz(1) end do !$OMP END PARALLEL DO case ('bulk-NCEP') call calc_sfwf_bulk_ncep(STF,FW,TFW,12) case ('partially-coupled') call calc_sfwf_partially_coupled(12) end select case('n-hour') sfwf_data_label = 'SFWF n-hour' if (thour00 >= sfwf_data_update) then call update_forcing_data( sfwf_data_time, & sfwf_data_time_min_loc, sfwf_interp_type, & sfwf_data_next, sfwf_data_update, & sfwf_data_type, sfwf_data_inc, & SFWF_DATA(:,:,:,:,1:sfwf_interp_order), & sfwf_data_renorm, & sfwf_data_label, sfwf_data_names, & sfwf_bndy_loc, sfwf_bndy_type, & sfwf_filename, sfwf_file_fmt) endif if (thour00 >= sfwf_interp_next .or. nsteps_run == 0) then call interpolate_forcing(SFWF_DATA(:,:,:,:,0), & SFWF_DATA(:,:,:,:,1:sfwf_interp_order), & sfwf_data_time, sfwf_interp_type, & sfwf_data_time_min_loc, sfwf_interp_freq, & sfwf_interp_inc, sfwf_interp_next, & sfwf_interp_last, nsteps_run) if (nsteps_run /= 0) sfwf_interp_next = & sfwf_interp_next + sfwf_interp_inc endif select case (sfwf_formulation) case ('restoring') !$OMP PARALLEL DO PRIVATE(iblock, this_block) do iblock=1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) STF(:,:,2,iblock) = & (SFWF_DATA(:,:,iblock,sfwf_data_sss,0) - & TRACER(:,:,1,2,curtime,iblock))* & sfwf_restore_rtau*dz(1) end do !$OMP END PARALLEL DO case ('bulk-NCEP') call calc_sfwf_bulk_ncep(STF, FW, TFW, sfwf_interp_order) case ('partially-coupled') call calc_sfwf_partially_coupled(sfwf_interp_order) end select end select !----------------------------------------------------------------------- !EOC end subroutine set_sfwf !*********************************************************************** !BOP ! !IROUTINE: calc_sfwf_bulk_ncep ! !INTERFACE: subroutine calc_sfwf_bulk_ncep(STF, FW, TFW, time_dim) ! !DESCRIPTION: ! Calculates surface freshwater flux from a combination of ! air-sea fluxes (precipitation, evaporation based on ! latent heat flux computed in calc\_shf\_bulk\_ncep), ! and restoring terms (due to restoring fields of SSS). ! ! Notes: ! the forcing data (on t-grid) are computed and ! stored in SFWF\_DATA(:,:,sfwf\_comp\_*,now) where: ! sfwf\_data\_sss is restoring SSS (psu) ! sfwf\_data\_precip is precipitation (m/y) ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: integer (int_kind), intent(in) :: & time_dim ! number of time points for interpolation ! !INPUT/OUTPUT PARAMETERS: real (r8), dimension(nx_block,ny_block,nt,max_blocks_clinic), & intent(inout) :: & STF, &! surface tracer fluxes for all tracers TFW ! tracer concentration in fresh water flux ! !OUTPUT PARAMETERS: real (r8), dimension(nx_block,ny_block,max_blocks_clinic), & intent(out) :: & FW ! fresh water flux if using varthick sfc layer !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (int_kind) :: & now, &! index for location of interpolated data k, n, &! dummy loop indices iblock ! block loop index real (r8) :: & dttmp, &! temporary time step variable fres_hor_ave, &! area-weighted mean of weak restoring fres_hor_area, &! total area of weak restoring area_glob_m_marg, &! total ocean area - marginal sea area vol_glob_m_marg, &! total ocean volume - marginal sea volume weak_mean ! mean weak restoring real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: & WORK1, WORK2 ! temporary work space type(block) :: & this_block ! block info for current block !----------------------------------------------------------------------- ! ! sfwf_weak_restore= weak(non-ice) restoring h2o flux per msu (kg/s/m^2/msu) ! sfwf_strong_restore= strong (ice) .. .. .. .. .. .. ! ! to calculate restoring factors, use mixed layer of 50m, ! and restoring time constant tau (days): ! ! F (kg/s/m^2/msu) ! tau = 6 : 2.77 ! tau = 30 : 0.55 ! tau = 182.5: 0.092 ! tau = 365 : 0.046 ! tau = 730 : 0.023 ! tau = Inf : 0.0 ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! set location where interpolated data exists. ! !----------------------------------------------------------------------- if (sfwf_data_type == 'annual') then now = 1 else now = 0 endif !----------------------------------------------------------------------- ! ! compute forcing terms for each block ! !----------------------------------------------------------------------- !$OMP PARALLEL DO PRIVATE(iblock, this_block) do iblock=1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) !*** compute evaporation from latent heat computed in shf forcing SFWF_COMP(:,:,iblock,sfwf_comp_evap) = & SHF_COMP(:,:,iblock,shf_comp_qlat)/latent_heat_vapor_mks !*** precipitation (kg/m^2/s) SFWF_COMP(:,:,iblock,sfwf_comp_precip) = & SFWF_DATA(:,:,iblock,sfwf_data_precip,now)*precip_fact ! *c1000/seconds_in_year ! convert m/y to Kg/m^2/s if needed !*** weak salinity restoring term !*** (note: OCN_WGT = 0. at land points) !*** will be subtracting global mean later, so compute !*** necessary terms for global mean SFWF_COMP(:,:,iblock,sfwf_comp_wrest) = & -sfwf_weak_restore*OCN_WGT(:,:,iblock)* & MASK_SR(:,:,iblock)* & (SFWF_DATA(:,:,iblock,sfwf_data_sss,now) - & TRACER(:,:,1,2,curtime,iblock)) WORK1(:,:,iblock) = TAREA(:,:,iblock)* & SFWF_COMP(:,:,iblock,sfwf_comp_wrest) WORK2(:,:,iblock) = TAREA(:,:,iblock)*OCN_WGT(:,:,iblock)* & MASK_SR(:,:,iblock) !*** strong salinity restoring term where (KMT(:,:,iblock) > 0) SFWF_COMP(:,:,iblock,sfwf_comp_srest) = & -sfwf_strong_restore*(c1 - OCN_WGT(:,:,iblock))* & (SFWF_DATA(:,:,iblock,sfwf_data_sss,now) - & TRACER(:,:,1,2,curtime,iblock)) endwhere where (KMT(:,:,iblock) > 0 .and. MASK_SR(:,:,iblock) == 0) SFWF_COMP(:,:,iblock,sfwf_comp_srest) = & -sfwf_strong_restore_ms* & (SFWF_DATA(:,:,iblock,sfwf_data_sss,now) - & TRACER(:,:,1,2,curtime,iblock)) endwhere end do ! block loop !$OMP END PARALLEL DO !---------------------------------------------------------------------- ! ! compute global mean of weak restoring term ! !---------------------------------------------------------------------- fres_hor_ave = global_sum(WORK1, distrb_clinic, field_loc_center) fres_hor_area = global_sum(WORK2, distrb_clinic, field_loc_center) weak_mean = fres_hor_ave/fres_hor_area !----------------------------------------------------------------------- ! ! finish computing forcing terms for each block ! !----------------------------------------------------------------------- !$OMP PARALLEL DO PRIVATE(iblock, this_block) do iblock=1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) !*** subtract mean from weak restoring term SFWF_COMP(:,:,iblock,sfwf_comp_wrest) = & SFWF_COMP(:,:,iblock,sfwf_comp_wrest) - OCN_WGT(:,:,iblock)* & MASK_SR(:,:,iblock)*weak_mean ! if variable thickness surface layer, compute net surface ! freshwater flux (kg/m^2/s) due to restoring terms only ! then compute freshwater flux due to P-E and convert to (m/s) ! then set the tracer content in the freshwater flux ! defaults are FW*SST for tracer 1 (temperature) ! 0 for salinity (really fresh water) ! 0 for all other tracers ! ! IF DATA ARE AVAILABLE... ! IMPLEMENT SUM OVER FRESHWATER TYPES (E,P,R,F,M) HERE: ! ! TFW(:,:,n) = FW_EVAP*TW_EVAP(:,:,n) + FW_PRECIP*TW_PRECIP(:,:,n) ! + FW_ROFF*TW_ROFF(:,:,n) + FW_MELT*TW_MELT(:,:,n) ! + FW_FREEZE*TW_FREEZE(:,:,n) ! ! where, for example FW_ROFF is the water flux from rivers ! and TW_ROFF(:,:,n) is the concentration of the nth tracer ! in the river water; similarly for water fluxes due to ! evaporation, precipitation, ice freezing, and ice melting. if (sfc_layer_type == sfc_layer_varthick .and. & .not. lfw_as_salt_flx) then STF(:,:,2,iblock) = SFWF_COMP(:,:,iblock,sfwf_comp_wrest) + & SFWF_COMP(:,:,iblock,sfwf_comp_srest) FW(:,:,iblock) = OCN_WGT(:,:,iblock)*MASK_SR(:,:,iblock)* & (SFWF_COMP(:,:,iblock,sfwf_comp_evap) + & SFWF_COMP(:,:,iblock,sfwf_comp_precip))* & fwmass_to_fwflux !*** fw same temp as ocean and no tracers in FW input TFW(:,:,1,iblock) = FW(:,:,iblock)* & TRACER(:,:,1,1,curtime,iblock) TFW(:,:,2:nt,iblock) = c0 !*** if rigid lid or old free surface form, compute !*** net surface freshwater flux (kg/m^2/s) else STF(:,:,2,iblock) = OCN_WGT(:,:,iblock)*MASK_SR(:,:,iblock)* & (SFWF_COMP(:,:,iblock,sfwf_comp_evap) + & SFWF_COMP(:,:,iblock,sfwf_comp_precip)) + & SFWF_COMP(:,:,iblock,sfwf_comp_wrest)+ & SFWF_COMP(:,:,iblock,sfwf_comp_srest) endif !*** convert surface freshwater flux (kg/m^2/s) to !*** salinity flux (msu*cm/s) STF(:,:,2,iblock) = STF(:,:,2,iblock)*salinity_factor !*** compute fields for accumulating annual-mean precipitation !*** over ocean points that are not marginal seas. WORK1(:,:,iblock) = merge(SFWF_COMP(:,:,iblock,sfwf_comp_precip)*& TAREA(:,:,iblock)*OCN_WGT(:,:,iblock), & c0, MASK_SR(:,:,iblock) > 0) !WORK2 = merge(FW_OLD(:,:,iblock)*TAREA(:,:,iblock), & ! c0, MASK_SR(:,:,iblock) > 0) end do ! block loop !$OMP END PARALLEL DO !----------------------------------------------------------------------- ! ! accumulate annual-mean precipitation over ocean points that are ! not marginal seas. ! !----------------------------------------------------------------------- if (avg_ts .or. back_to_back) then dttmp = p5*dtt else dttmp = dtt endif area_glob_m_marg = area_t - area_t_marg sum_precip = sum_precip + dttmp*1.0e-4_r8* & global_sum(WORK1,distrb_clinic,field_loc_center)/ & area_glob_m_marg !sum_fw = sum_fw + & ! dttmp*global_sum(WORK2,distrb_clinic,field_loc_center)/ & ! area_glob_m_marg !----------------------------------------------------------------------- ! ! Perform end of year adjustment calculations ! !----------------------------------------------------------------------- if (eoy) then !*** Compute the surface volume-averaged salinity and !*** average surface height (for variable thickness sfc layer) !*** note that it is evaluated at the current time level. !$OMP PARALLEL DO PRIVATE(iblock, this_block) do iblock=1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) if (sfc_layer_type == sfc_layer_varthick) then WORK1(:,:,iblock) = merge( & TRACER(:,:,1,2,curtime,iblock)*TAREA(:,:,iblock)* & (dz(1) + PSURF(:,:,curtime,iblock)/grav), & c0, KMT(:,:,iblock) > 0 .and. MASK_SR(:,:,iblock) > 0) WORK2(:,:,iblock) = merge(PSURF(:,:,curtime,iblock)* & TAREA(:,:,iblock)/grav, c0, & KMT(:,:,iblock) > 0 .and. & MASK_SR(:,:,iblock) > 0) else WORK1(:,:,iblock) = merge(TRACER(:,:,1,2,curtime,iblock)* & TAREA(:,:,iblock)*dz(1), & c0, KMT(:,:,iblock) > 0 .and. & MASK_SR(:,:,iblock) > 0) endif end do !$OMP END PARALLEL DO vol_glob_m_marg = volume_t_k(1) - volume_t_marg_k(1) sal_final(1) = 1.0e-6_r8* & global_sum(WORK1, distrb_clinic, field_loc_center)/ & vol_glob_m_marg if (sfc_layer_type == sfc_layer_varthick) then ssh_final = 1.0e-4_r8* & global_sum(WORK2, distrb_clinic, field_loc_center)/ & area_glob_m_marg/seconds_in_year if (my_task == master_task) then write(stdout,blank_fmt) write(stdout,'(a22,1pe23.15)') & 'annual change in SSH: ', ssh_final endif ssh_final = ssh_final*10.0_r8 ! convert (cm/s) -> kg/m^2/s ! (cm/s)x0.01(m/cm)x1000kg/m^3 else ssh_final = c0 endif !*** Compute the volume-averaged salinity for each level. do k=2,km !$OMP PARALLEL DO PRIVATE(iblock, this_block) do iblock=1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) if (partial_bottom_cells) then WORK1(:,:,iblock) = & merge(TRACER(:,:,k,2,curtime,iblock)* & TAREA(:,:,iblock)*DZT(:,:,k,iblock), c0, & k <= KMT(:,:,iblock) .and. & MASK_SR(:,:,iblock) > 0) else WORK1(:,:,iblock) = & merge(TRACER(:,:,k,2,curtime,iblock)* & TAREA(:,:,iblock)*dz(k), c0, & k <= KMT(:,:,iblock) .and. & MASK_SR(:,:,iblock) > 0) endif end do !$OMP END PARALLEL DO vol_glob_m_marg = volume_t_k(k) - volume_t_marg_k(k) if (vol_glob_m_marg == 0) vol_glob_m_marg = 1.e+20_r8 sal_final(k) = 1.0e-6_r8* & global_sum(WORK1, distrb_clinic, field_loc_center)/ & vol_glob_m_marg enddo !*** find annual mean precip and reset annual counters !ann_avg_fw = sum_fw / seconds_in_year !if (my_task == master_task) then ! write(stdout,blank_fmt) ! write(stdout,'(a32,1pe22.15)') & ! 'annual average freshwater flux: ', ann_avg_fw !endif !sum_fw = c0 ann_avg_precip = sum_precip / seconds_in_year sum_precip = c0 if (ladjust_precip) call precip_adjustment sal_initial = sal_final ssh_initial = ssh_final endif ! end of year calculations !----------------------------------------------------------------------- !EOC end subroutine calc_sfwf_bulk_ncep !*********************************************************************** !BOP ! !IROUTINE: calc_sfwf_partially_coupled ! !INTERFACE: subroutine calc_sfwf_partially_coupled(time_dim) ! !DESCRIPTION: ! Calculate ice-ocean flux, weak restoring, and strong restoring ! components of surface freshwater flux for partially-coupled formulation. ! these components will later be used in ! set\_surface\_forcing (forcing.F) to form the total surface freshwater ! (salt) flux. ! ! the forcing data (on t-grid) sets needed are ! sfwf\_data\_sss, restoring SSS (msu) ! sfwf\_data\_flxio, diagnosed ("climatological") (kg/m^2/s) ! ice-ocean freshwater flux ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: integer (int_kind), intent(in) :: & time_dim ! number of time points for interpolation !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (int_kind) :: & now, &! index for location of interpolated data k, n, &! dummy loop indices iblock ! block loop index real (r8) :: & dttmp, &! temporary time step variable fres_hor_ave, &! area-weighted mean of weak restoring fres_hor_area, &! total area of weak restoring area_glob_m_marg, &! total ocean area - marginal sea area vol_glob_m_marg ! total ocean volume - marginal sea volume real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: & WORK,WORK1, WORK2 ! temporary work space type(block) :: & this_block ! block info for current block !----------------------------------------------------------------------- ! ! set location where interpolated data exists. ! !----------------------------------------------------------------------- area_glob_m_marg = 1.0e-4*(area_t - area_t_marg) if (sfwf_data_type == 'annual') then now = 1 else now = 0 endif !----------------------------------------------------------------------- ! ! compute forcing terms for each block ! !----------------------------------------------------------------------- !$OMP PARALLEL DO PRIVATE(iblock, this_block) do iblock=1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) if (.not. luse_cpl_ifrac) then WORK(:,:,iblock) = OCN_WGT(:,:,iblock)*MASK_SR(:,:,iblock) else WORK(:,:,iblock) = MASK_SR(:,:,iblock) endif !*** weak salinity restoring term !*** (note: MASK_SR = 0. at land and marginal-sea points) !*** will be subtracting global mean later, so compute !*** necessary terms for global mean SFWF_COMP(:,:,iblock,sfwf_comp_wrest) = & -sfwf_weak_restore*WORK(:,:,iblock)* & (SFWF_DATA(:,:,iblock,sfwf_data_sss,now) - & TRACER(:,:,1,2,curtime,iblock)) WORK1(:,:,iblock) = TAREA(:,:,iblock)* & SFWF_COMP(:,:,iblock,sfwf_comp_wrest) WORK2(:,:,iblock) = TAREA(:,:,iblock)*WORK(:,:,iblock) end do ! block loop !$OMP END PARALLEL DO !---------------------------------------------------------------------- ! ! compute global mean of weak restoring term ! !---------------------------------------------------------------------- fres_hor_ave = global_sum(WORK1, distrb_clinic, field_loc_center) fres_hor_area = global_sum(WORK2, distrb_clinic, field_loc_center) !----------------------------------------------------------------------- ! ! finish computing forcing terms for each block ! !----------------------------------------------------------------------- !$OMP PARALLEL DO PRIVATE(iblock, this_block) do iblock=1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) ! subtract global mean from weak restoring term SFWF_COMP(:,:,iblock,sfwf_comp_wrest) = & SFWF_COMP(:,:,iblock,sfwf_comp_wrest) - WORK(:,:,iblock)* & fres_hor_ave/fres_hor_area where (KMT(:,:,iblock) > 0) SFWF_COMP(:,:,iblock,sfwf_comp_srest) = & -sfwf_strong_restore*(c1 - OCN_WGT(:,:,iblock))* & (SFWF_DATA(:,:,iblock,sfwf_data_sss,now) - & TRACER(:,:,1,2,curtime,iblock)) endwhere where (KMT(:,:,iblock) > 0 .and. MASK_SR(:,:,iblock) == 0) SFWF_COMP(:,:,iblock,sfwf_comp_srest) = & -sfwf_strong_restore_ms* & (SFWF_DATA(:,:,iblock,sfwf_data_sss,now) - & TRACER(:,:,1,2,curtime,iblock)) endwhere !*** ice-ocean climatological flux term if ( .not. lactive_ice ) then where (KMT(:,:,iblock) > 0) SFWF_COMP(:,:,iblock,sfwf_comp_flxio) = & (c1 - OCN_WGT(:,:,iblock))* & SFWF_DATA(:,:,iblock,sfwf_data_flxio,now) endwhere if ( .not. lms_balance ) & SFWF_COMP(:,:,iblock,sfwf_comp_flxio) = & SFWF_COMP(:,:,iblock,sfwf_comp_flxio)*MASK_SR(:,:,iblock) endif !*** convert surface freshwater flux components (kg/m^2/s) to !*** salinity flux components (msu*cm/s) SFWF_COMP(:,:,iblock,sfwf_comp_wrest) = & SFWF_COMP(:,:,iblock,sfwf_comp_wrest)*salinity_factor SFWF_COMP(:,:,iblock,sfwf_comp_srest) = & SFWF_COMP(:,:,iblock,sfwf_comp_srest)*salinity_factor if ( sfc_layer_type == sfc_layer_varthick .and. & .not. lfw_as_salt_flx) then WORK(:,:,iblock) = fwmass_to_fwflux * & SFWF_COMP(:,:,iblock,sfwf_comp_flxio) SFWF_COMP(:,:,iblock,sfwf_comp_flxio) = WORK(:,:,iblock) call tmelt(WORK1(:,:,iblock),TRACER(:,:,1,2,curtime,iblock)) TFW_COMP(:,:,1, iblock,tfw_comp_flxio) = WORK(:,:,iblock)* & WORK1(:,:,iblock) TFW_COMP(:,:,2:nt,iblock,tfw_comp_flxio) = c0 else SFWF_COMP(:,:,iblock,sfwf_comp_flxio) = salinity_factor* & SFWF_COMP(:,:,iblock,sfwf_comp_flxio) endif end do ! block loop !$OMP END PARALLEL DO !----------------------------------------------------------------------- ! ! Perform end of year adjustment calculations ! !----------------------------------------------------------------------- if (eoy) then !*** Compute the surface volume-averaged salinity and !*** average surface height (for variable thickness sfc layer) !*** note that it is evaluated at the current time level. !$OMP PARALLEL DO PRIVATE(iblock, this_block) do iblock=1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) if (sfc_layer_type == sfc_layer_varthick) then WORK1(:,:,iblock) = merge( & TRACER(:,:,1,2,curtime,iblock)*TAREA(:,:,iblock)* & (dz(1) + PSURF(:,:,curtime,iblock)/grav), & c0, KMT(:,:,iblock) > 0 .and. MASK_SR(:,:,iblock) > 0) WORK2(:,:,iblock) = merge(PSURF(:,:,curtime,iblock)* & TAREA(:,:,iblock)/grav, c0, & KMT(:,:,iblock) > 0 .and. & MASK_SR(:,:,iblock) > 0) else WORK1(:,:,iblock) = merge(TRACER(:,:,1,2,curtime,iblock)* & TAREA(:,:,iblock)*dz(1), & c0, KMT(:,:,iblock) > 0 .and. & MASK_SR(:,:,iblock) > 0) endif end do !$OMP END PARALLEL DO vol_glob_m_marg = 1.0e-6_r8*(volume_t_k(1) - volume_t_marg_k(1)) sal_final(1) = 1.0e-6_r8* &! convert to m^3 global_sum(WORK1,distrb_clinic,field_loc_center)/vol_glob_m_marg if (sfc_layer_type == sfc_layer_varthick) then ssh_final = 1.0e-4_r8* &! convert to m^3 global_sum(WORK2, distrb_clinic, field_loc_center)/ & area_glob_m_marg/seconds_in_year if (my_task == master_task) then write(stdout,blank_fmt) write(stdout,'(a22,1pe23.15)') & 'annual change in SSH: ', ssh_final endif ssh_final = ssh_final*10.0_r8 ! convert (cm/s) -> kg/m^2/s ! (cm/s)x0.01(m/cm)x1000kg/m^3 else ssh_final = c0 endif !*** Compute the volume-averaged salinity for each level. do k=2,km !$OMP PARALLEL DO PRIVATE(iblock, this_block) do iblock=1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) if (partial_bottom_cells) then WORK1(:,:,iblock) = & merge(TRACER(:,:,k,2,curtime,iblock)* & TAREA(:,:,iblock)*DZT(:,:,k,iblock), c0, & k <= KMT(:,:,iblock) .and. & MASK_SR(:,:,iblock) > 0) else WORK1(:,:,iblock) = & merge(TRACER(:,:,k,2,curtime,iblock)* & TAREA(:,:,iblock)*dz(k), c0, & k <= KMT(:,:,iblock) .and. & MASK_SR(:,:,iblock) > 0) endif end do !$OMP END PARALLEL DO vol_glob_m_marg = 1.0e-6_r8*(volume_t_k(k) - volume_t_marg_k(k)) if (vol_glob_m_marg == 0) vol_glob_m_marg = 1.e+20_r8 sal_final(k) = 1.0e-6_r8* &! convert to m^3 global_sum(WORK1,distrb_clinic,field_loc_center)/vol_glob_m_marg enddo if (ladjust_precip) call precip_adjustment sal_initial = sal_final ssh_initial = ssh_final endif ! end of year calculations !----------------------------------------------------------------------- !EOC end subroutine calc_sfwf_partially_coupled !*********************************************************************** !BOP ! !IROUTINE: precip_adjustment ! !INTERFACE: subroutine precip_adjustment ! !DESCRIPTION: ! Computes a precipitation factor to multiply the fresh water flux ! due to precipitation uniformly to insure a balance of fresh water ! at the ocean surface. ! ! !REVISION HISTORY: ! same as module !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- real (r8) :: & sal_tendency, fw_tendency, precip_tav, & area_glob_m_marg, &! global ocean area - marginal sea area (cm^2) vol_glob_m_marg ! global ocean vol - marginal sea vol (cm^3) integer (int_kind) :: k !----------------------------------------------------------------------- ! ! compute tendency of salinity for eack "k" layer, considering the ! effects of depth acceleration ! !----------------------------------------------------------------------- do k=1,km sal_initial(k) = (sal_final(k) - sal_initial(k))/ & (dttxcel(k)*seconds_in_year) enddo !----------------------------------------------------------------------- ! ! form the global volume-averaged tendency to be used in "precip_fact" ! computation ! !----------------------------------------------------------------------- sal_tendency = c0 do k=1,km vol_glob_m_marg = 1.0e-6_r8*(volume_t_k(k) - volume_t_marg_k(k)) sal_tendency = sal_tendency + vol_glob_m_marg*sal_initial(k) enddo vol_glob_m_marg = 1.0e-6_r8*(volume_t - volume_t_marg) sal_tendency = sal_tendency/vol_glob_m_marg if (my_task == master_task) then write (stdout,'(a58,1pe22.15)') & ' precip_adjustment: volume-averaged salinity tendency = ', & sal_tendency endif !----------------------------------------------------------------------- ! ! convert "sal_tendency" from (msu/s) to -(kg/m^2/s). note that ! areag in cm^2 and volgt in cm^3 ! assumes density of fresh water = 1000 kg/m**3 ! !----------------------------------------------------------------------- area_glob_m_marg = 1.0e-4*(area_t - area_t_marg) sal_tendency = - sal_tendency*vol_glob_m_marg* & 1.0e6_r8/area_glob_m_marg/ocn_ref_salinity !----------------------------------------------------------------------- ! ! compute annual change in mass due to freshwater flux (kg/m^2/s) ! !----------------------------------------------------------------------- fw_tendency = ssh_final - ssh_initial if (my_task == master_task) then write (stdout,'(a22)') ' precip_adjustment: ' write (stdout,'(a28,1pe22.15)') ' sal_tendency (kg/m^2/s): ', & sal_tendency write (stdout,'(a28,1pe22.15)') ' fw_tendency (kg/m^2/s): ', & fw_tendency endif !----------------------------------------------------------------------- ! ! change "precip_fact" based on tendency of freshwater and previous ! amount of precipitation ! !----------------------------------------------------------------------- if (sfwf_formulation == 'partially-coupled') then precip_tav = precip_mean else precip_tav = ann_avg_precip/precip_fact endif precip_fact = precip_fact - & (sal_tendency + fw_tendency)/precip_tav if (my_task == master_task) then write (stdout,'(a33,e14.8)') ' Changed precipitation factor to ',& precip_fact endif !----------------------------------------------------------------------- !EOC end subroutine precip_adjustment !*********************************************************************** end module forcing_sfwf !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||