From d4bd9a61cc37b35bfdfbb4b1b071a897194c2519 Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Wed, 5 Jul 2023 13:00:17 -0600 Subject: [PATCH] Add Leith+E This commit adds the 2D Leith+E closure, which uses a modified 2D Leith biharmonic viscosity paired with a harmonic backscatter. ('Modified' here is not used in the same sense as 'modified 2D Leith'; it just means that the biharmonic coefficient is modified to account for enstrophy backscatter.) Variables are often named 'leithy' to refer to Leith+E. The parameterization is controlled by three main entries in user_nl_mom: 1. USE_LEITHY = True 2. LEITH_CK = 1.0 3. LEITH_BI_CONST = 8.0 To use Leith+E you should have LAPLACIAN=True and BIHARMONIC=True. (It doesn't hurt to be explicit and also set LEITH_AH=False, along with any other viscous closures, but this is not required. If USE_LEITHY=True it will not use any of the other schemes.) LEITH_CK is the fraction of energy dissipated by the biharmonic term that gets backscattered by the harmonic term (it's a target; the backscatter rate is not exact.) Recommended values between 0 and 1. LEITH_BI_CONST is Upsilon^6 where Upsilon is the ratio between the dissipation scale for enstrophy and the grid scale. Values should probably be greater than or equal to 1. This commit also includes a merge with dev/ncar. Leith+E only changes MOM_hor_visc.F90; all other changes are associated with the merge. --- config_src/drivers/nuopc_cap/mom_cap.F90 | 64 ++- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 20 +- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 40 +- .../solo_driver/MOM_surface_forcing.F90 | 3 +- .../infra/FMS1/MOM_ensemble_manager_infra.F90 | 13 +- .../infra/FMS2/MOM_ensemble_manager_infra.F90 | 15 +- src/core/MOM.F90 | 8 +- src/core/MOM_forcing_type.F90 | 44 +- src/core/MOM_unit_tests.F90 | 2 +- src/core/MOM_variables.F90 | 15 +- src/diagnostics/MOM_obsolete_params.F90 | 2 + src/framework/MOM_file_parser.F90 | 9 +- src/framework/MOM_get_input.F90 | 2 +- .../lateral/MOM_hor_visc.F90 | 378 ++++++++++++++++-- src/tracer/MOM_CFC_cap.F90 | 334 +++++++++------- ...iffusion.F90 => MOM_hor_bnd_diffusion.F90} | 345 ++++++++++------ src/tracer/MOM_neutral_diffusion.F90 | 2 +- src/tracer/MOM_tracer_flow_control.F90 | 11 +- src/tracer/MOM_tracer_hor_diff.F90 | 79 ++-- src/tracer/MOM_tracer_registry.F90 | 74 ++-- src/tracer/MOM_tracer_types.F90 | 30 +- src/tracer/oil_tracer.F90 | 2 +- 22 files changed, 967 insertions(+), 525 deletions(-) mode change 100644 => 100755 src/parameterizations/lateral/MOM_hor_visc.F90 rename src/tracer/{MOM_lateral_boundary_diffusion.F90 => MOM_hor_bnd_diffusion.F90} (80%) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index b7d651bf55..2841c7196c 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -28,9 +28,11 @@ module MOM_cap_mod use MOM_cap_methods, only: mom_import, mom_export, mom_set_geomtype, mod2med_areacor use MOM_cap_methods, only: med2mod_areacor, state_diagnose use MOM_cap_methods, only: ChkErr +use MOM_ensemble_manager, only: ensemble_manager_init #ifdef CESMCOUPLED use shr_log_mod, only: shr_log_setLogUnit +use nuopc_shr_methods, only: get_component_instance #endif use time_utils_mod, only: esmf2fms_time @@ -146,7 +148,8 @@ module MOM_cap_mod logical :: cesm_coupled = .false. type(ESMF_GeomType_Flag) :: geomtype #endif -character(len=8) :: restart_mode = 'alarms' +character(len=8) :: restart_mode = 'alarms' +character(len=16) :: inst_suffix = '' contains @@ -422,6 +425,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ! (same as restartfile if single restart file) character(len=*), parameter :: subname='(MOM_cap:InitializeAdvertise)' character(len=32) :: calendar + character(len=:), allocatable :: rpointer_filename + integer :: inst_index !-------------------------------- rc = ESMF_SUCCESS @@ -451,6 +456,13 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) CALL ESMF_TimeIntervalGet(TINT, S=DT_OCEAN, RC=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return +#ifdef CESMCOUPLED + call get_component_instance(gcomp, inst_suffix, inst_index, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ensemble_manager_init(inst_suffix) + rpointer_filename = 'rpointer.ocn'//trim(inst_suffix) +#endif + ! reset shr logging to my log file if (localPet==0) then call NUOPC_CompAttributeGet(gcomp, name="diro", & @@ -460,11 +472,19 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) isPresent=isPresentLogfile, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return if (isPresentDiro .and. isPresentLogfile) then - call NUOPC_CompAttributeGet(gcomp, name="diro", value=diro, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call NUOPC_CompAttributeGet(gcomp, name="logfile", value=logfile, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - open(newunit=stdout,file=trim(diro)//"/"//trim(logfile)) + call NUOPC_CompAttributeGet(gcomp, name="diro", value=diro, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call NUOPC_CompAttributeGet(gcomp, name="logfile", value=logfile, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + if (cesm_coupled) then + ! Multiinstance logfile name needs a correction + if(logfile(4:4) == '_') then + logfile = logfile(1:3)//trim(inst_suffix)//logfile(9:) + endif + endif + + open(newunit=stdout,file=trim(diro)//"/"//trim(logfile)) else stdout = output_unit endif @@ -521,12 +541,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) time0 = set_date (YEAR,MONTH,DAY,HOUR,MINUTE,SECOND) - - ! rsd need to figure out how to get this without share code - !call shr_nuopc_get_component_instance(gcomp, inst_suffix, inst_index) - !inst_name = "OCN"//trim(inst_suffix) - - if (is_root_pe()) then write(stdout,*) subname//'start time: y,m,d-',year,month,day,'h,m,s=',hour,minute,second endif @@ -581,9 +595,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) if (localPet == 0) then ! this hard coded for rpointer.ocn right now - open(newunit=readunit, file='rpointer.ocn', form='formatted', status='old', iostat=iostat) + open(newunit=readunit, file=rpointer_filename, form='formatted', status='old', iostat=iostat) if (iostat /= 0) then - call ESMF_LogSetError(ESMF_RC_FILE_OPEN, msg=subname//' ERROR opening rpointer.ocn', & + call ESMF_LogSetError(ESMF_RC_FILE_OPEN, msg=subname//' ERROR opening '//rpointer_filename, & line=__LINE__, file=u_FILE_u, rcToReturn=rc) return endif @@ -593,7 +607,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) if (len(trim(restartfiles))>1 .and. iostat<0) then exit ! done reading restart files list. else - call ESMF_LogSetError(ESMF_RC_FILE_READ, msg=subname//' ERROR reading rpointer.ocn', & + call ESMF_LogSetError(ESMF_RC_FILE_READ, msg=subname//' ERROR reading '//rpointer_filename, & line=__LINE__, file=u_FILE_u, rcToReturn=rc) return endif @@ -616,7 +630,12 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) endif ocean_public%is_ocean_pe = .true. - call ocean_model_init(ocean_public, ocean_state, time0, time_start, input_restart_file=trim(adjustl(restartfiles))) + if (cesm_coupled .and. len_trim(inst_suffix)>0) then + call ocean_model_init(ocean_public, ocean_state, time0, time_start, & + input_restart_file=trim(adjustl(restartfiles)), inst_index=inst_index) + else + call ocean_model_init(ocean_public, ocean_state, time0, time_start, input_restart_file=trim(adjustl(restartfiles))) + endif ! GMM, this call is not needed in CESM. Check with EMC if it can be deleted. call ocean_model_flux_init(ocean_state) @@ -1489,6 +1508,7 @@ subroutine ModelAdvance(gcomp, rc) character(len=128) :: fldname character(len=*),parameter :: subname='(MOM_cap:ModelAdvance)' character(len=8) :: suffix + character(len=:), allocatable :: rpointer_filename integer :: num_rest_files rc = ESMF_SUCCESS @@ -1658,6 +1678,8 @@ subroutine ModelAdvance(gcomp, rc) call ESMF_VMGet(vm, localPet=localPet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + rpointer_filename = 'rpointer.ocn'//trim(inst_suffix) + write(restartname,'(A,".mom6.r.",I4.4,"-",I2.2,"-",I2.2,"-",I5.5)') & trim(casename), year, month, day, seconds call ESMF_LogWrite("MOM_cap: Writing restart : "//trim(restartname), ESMF_LOGMSG_INFO) @@ -1665,13 +1687,17 @@ subroutine ModelAdvance(gcomp, rc) call ocean_model_restart(ocean_state, restartname=restartname, num_rest_files=num_rest_files) if (localPet == 0) then ! Write name of restart file in the rpointer file - this is currently hard-coded for the ocean - open(newunit=writeunit, file='rpointer.ocn', form='formatted', status='unknown', iostat=iostat) + open(newunit=writeunit, file=rpointer_filename, form='formatted', status='unknown', iostat=iostat) if (iostat /= 0) then call ESMF_LogSetError(ESMF_RC_FILE_OPEN, & - msg=subname//' ERROR opening rpointer.ocn', line=__LINE__, file=u_FILE_u, rcToReturn=rc) + msg=subname//' ERROR opening '//rpointer_filename, line=__LINE__, file=u_FILE_u, rcToReturn=rc) return endif - write(writeunit,'(a)') trim(restartname)//'.nc' + if (len_trim(inst_suffix) == 0) then + write(writeunit,'(a)') trim(restartname)//'.nc' + else + write(writeunit,'(a)') trim(restartname)//'.'//trim(inst_suffix)//'.nc' + endif if (num_rest_files > 1) then ! append i.th restart file name to rpointer diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 1283b98ba0..04b60b0d37 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -41,7 +41,7 @@ module MOM_ocean_model_nuopc use MOM_time_manager, only : operator(/=), operator(<=), operator(>=) use MOM_time_manager, only : operator(<), real_to_time_type, time_type_to_real use MOM_interpolate, only : time_interp_external_init -use MOM_tracer_flow_control, only : call_tracer_flux_init +use MOM_tracer_flow_control, only : tracer_flow_control_CS, call_tracer_flux_init, call_tracer_set_forcing use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type @@ -210,6 +210,8 @@ module MOM_ocean_model_nuopc type(marine_ice_CS), pointer :: & marine_ice_CSp => NULL() !< A pointer to the control structure for the !! marine ice effects module. + type(tracer_flow_control_CS), pointer :: & + tracer_flow_CSp => NULL() !< A pointer to the tracer flow control structure type(wave_parameters_CS), pointer, public :: & Waves => NULL() !< A pointer to the surface wave control structure type(surface_forcing_CS), pointer :: & @@ -229,7 +231,7 @@ module MOM_ocean_model_nuopc !! This subroutine initializes both the ocean state and the ocean surface type. !! Because of the way that indicies and domains are handled, Ocean_sfc must have !! been used in a previous call to initialize_ocean_type. -subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, input_restart_file) +subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, input_restart_file, inst_index) type(ocean_public_type), target, & intent(inout) :: Ocean_sfc !< A structure containing various publicly !! visible ocean surface properties after initialization, @@ -246,6 +248,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i !! tracer fluxes, and can be used to spawn related !! internal variables in the ice model. character(len=*), optional, intent(in) :: input_restart_file !< If present, name of restart file to read + integer, optional :: inst_index !< Ensemble index provided by the cap (instead of FMS ensemble manager) ! Local variables real :: Rho0 ! The Boussinesq ocean density, in kg m-3. @@ -255,7 +258,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i !! min(HFrz, OBLD), where OBLD is the boundary layer depth. !! If HFrz <= 0 (default), melt potential will not be computed. logical :: use_melt_pot !< If true, allocate melt_potential array - logical :: use_CFC !< If true, allocated arrays for surface CFCs. ! This include declares and sets the variable "version". @@ -283,7 +285,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i call initialize_MOM(OS%Time, Time_init, param_file, OS%dirs, OS%MOM_CSp, & OS%restart_CSp, Time_in, offline_tracer_mode=OS%offline_tracer_mode, & input_restart_file=input_restart_file, & - diag_ptr=OS%diag, count_calls=.true., waves_CSp=OS%Waves) + diag_ptr=OS%diag, count_calls=.true., tracer_flow_CSp=OS%tracer_flow_CSp, & + waves_CSp=OS%Waves, ensemble_num=inst_index) call get_MOM_state_elements(OS%MOM_CSp, G=OS%grid, GV=OS%GV, US=OS%US, C_p=OS%C_p, & C_p_scaled=OS%fluxes%C_p, use_temp=use_temperature) @@ -376,8 +379,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i use_melt_pot=.false. endif - call get_param(param_file, mdl, "USE_CFC_CAP", use_CFC, & - default=.false., do_not_log=.true.) call get_param(param_file, mdl, "USE_WAVES", OS%Use_Waves, & "If true, enables surface wave modules.", default=.false.) @@ -385,7 +386,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i ! vertical integrals, since the related 3-d sums are not negligible in cost. call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, & do_integrals=.true., gas_fields_ocn=gas_fields_ocn, & - use_meltpot=use_melt_pot, use_cfcs=use_CFC) + use_meltpot=use_melt_pot) call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, & OS%forcing_CSp, OS%restore_salinity, OS%restore_temp, OS%use_waves) @@ -611,6 +612,11 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call finish_MOM_initialization(OS%Time, OS%dirs, OS%MOM_CSp, OS%restart_CSp) endif + if (do_thermo) & + call call_tracer_set_forcing(OS%sfc_state, OS%fluxes, OS%Time, & + real_to_time_type(dt_coupling), OS%grid, OS%US, OS%GV%Rho0, & + OS%tracer_flow_CSp) + call disable_averaging(OS%diag) Master_time = OS%Time ; Time1 = OS%Time diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index 2c8e3db8bd..b921f7355d 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -26,7 +26,6 @@ module MOM_surface_forcing_nuopc use MOM_grid, only : ocean_grid_type use MOM_interpolate, only : init_external_field, time_interp_external use MOM_interpolate, only : time_interp_external_init -use MOM_CFC_cap, only : CFC_cap_fluxes use MOM_io, only : slasher, write_version_number, MOM_read_data use MOM_io, only : stdout use MOM_restart, only : register_restart_field, restart_init, MOM_restart_CS @@ -129,7 +128,6 @@ module MOM_surface_forcing_nuopc type(diag_ctrl), pointer :: diag !< structure to regulate diagnostic output timing character(len=200) :: inputdir !< directory where NetCDF input files are - character(len=200) :: CFC_BC_file !< filename with cfc11 and cfc12 data character(len=200) :: salt_restore_file !< filename for salt restoring data character(len=30) :: salt_restore_var_name !< name of surface salinity in salt_restore_file logical :: mask_srestore !< if true, apply a 2-dimensional mask to the surface @@ -143,13 +141,9 @@ module MOM_surface_forcing_nuopc !! temperature restoring fluxes. The masking file should be !! in inputdir/temp_restore_mask.nc and the field should !! be named 'mask' - character(len=30) :: cfc11_var_name !< name of cfc11 in CFC_BC_file - character(len=30) :: cfc12_var_name !< name of cfc11 in CFC_BC_file real, pointer, dimension(:,:) :: trestore_mask => NULL() !< mask for SST restoring - integer :: id_srestore = -1 !< id number for time_interp_external. - integer :: id_trestore = -1 !< id number for time_interp_external. - integer :: id_cfc11_atm = -1 !< id number for time_interp_external. - integer :: id_cfc12_atm = -1 !< id number for time_interp_external. + integer :: id_srestore = -1 !< id number for time_interp_external. + integer :: id_trestore = -1 !< id number for time_interp_external. ! Diagnostics handles type(forcing_diags), public :: handles @@ -245,8 +239,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! local variables real, dimension(SZI_(G),SZJ_(G)) :: & - cfc11_atm, & !< CFC11 concentration in the atmopshere [???????] - cfc12_atm, & !< CFC11 concentration in the atmopshere [???????] data_restore, & !< The surface value toward which to restore [S ~> ppt] or [C ~> degC] PmE_adj, & !< The adjustment to PminusE that will cause the salinity !! to be restored toward its target value [kg/(m^2 * s)] @@ -594,11 +586,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, fluxes%accumulate_p_surf = .true. ! Multiple components may contribute to surface pressure. endif - ! CFCs - if (CS%use_CFC) then - call CFC_cap_fluxes(fluxes, sfc_state, G, US, CS%Rho0, Time, CS%id_cfc11_atm, CS%id_cfc11_atm) - endif - if (associated(IOB%salt_flux)) then do j=js,je ; do i=is,ie fluxes%salt_flux(i,j) = G%mask2dT(i,j)*(fluxes%salt_flux(i,j) + kg_m2_s_conversion*IOB%salt_flux(i-i0,j-j0)) @@ -1412,29 +1399,6 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, endif endif ; endif - ! Do not log these params here since they are logged in the CFC cap module - if (CS%use_CFC) then - call get_param(param_file, mdl, "CFC_BC_FILE", CS%CFC_BC_file, & - "The file in which the CFC-11 and CFC-12 atm concentrations can be "//& - "found (units must be parts per trillion), or an empty string for "//& - "internal BC generation (TODO).", default=" ", do_not_log=.true.) - if ((len_trim(CS%CFC_BC_file) > 0) .and. (scan(CS%CFC_BC_file,'/') == 0)) then - ! Add the directory if CFC_BC_file is not already a complete path. - CS%CFC_BC_file = trim(CS%inputdir) // trim(CS%CFC_BC_file) - endif - if (len_trim(CS%CFC_BC_file) > 0) then - call get_param(param_file, mdl, "CFC11_VARIABLE", CS%cfc11_var_name, & - "The name of the variable representing CFC-11 in "//& - "CFC_BC_FILE.", default="CFC_11", do_not_log=.true.) - call get_param(param_file, mdl, "CFC12_VARIABLE", CS%cfc12_var_name, & - "The name of the variable representing CFC-12 in "//& - "CFC_BC_FILE.", default="CFC_12", do_not_log=.true.) - - CS%id_cfc11_atm = init_external_field(CS%CFC_BC_file, CS%cfc11_var_name, domain=G%Domain%mpp_domain) - CS%id_cfc12_atm = init_external_field(CS%CFC_BC_file, CS%cfc12_var_name, domain=G%Domain%mpp_domain) - endif - endif - ! Set up any restart fields associated with the forcing. call restart_init(param_file, CS%restart_CSp, "MOM_forcing.res") call restart_init_end(CS%restart_CSp) diff --git a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 b/config_src/drivers/solo_driver/MOM_surface_forcing.F90 index 092bc9e513..522420e004 100644 --- a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 +++ b/config_src/drivers/solo_driver/MOM_surface_forcing.F90 @@ -346,7 +346,8 @@ subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, G, US endif if (associated(CS%tracer_flow_CSp)) then - call call_tracer_set_forcing(sfc_state, fluxes, day_start, day_interval, G, CS%tracer_flow_CSp) + call call_tracer_set_forcing(sfc_state, fluxes, day_start, day_interval, G, US, CS%Rho0, & + CS%tracer_flow_CSp) endif ! Allow for user-written code to alter the fluxes after all the above diff --git a/config_src/infra/FMS1/MOM_ensemble_manager_infra.F90 b/config_src/infra/FMS1/MOM_ensemble_manager_infra.F90 index 66bbb86e2f..3ab9d591da 100644 --- a/config_src/infra/FMS1/MOM_ensemble_manager_infra.F90 +++ b/config_src/infra/FMS1/MOM_ensemble_manager_infra.F90 @@ -9,6 +9,7 @@ module MOM_ensemble_manager_infra use ensemble_manager_mod, only : FMS_get_ensemble_size => get_ensemble_size use ensemble_manager_mod, only : FMS_get_ensemble_pelist => get_ensemble_pelist use ensemble_manager_mod, only : FMS_get_ensemble_filter_pelist => get_ensemble_filter_pelist +use fms_io_mod, only : fms_io_set_filename_appendix=>set_filename_appendix implicit none ; private @@ -20,9 +21,15 @@ module MOM_ensemble_manager_infra !> Initializes the ensemble manager which divides available resources !! in order to concurrently execute an ensemble of model realizations. -subroutine ensemble_manager_init() - - call FMS_ensemble_manager_init() +subroutine ensemble_manager_init(ensemble_suffix) + character(len=*), optional, intent(in) :: ensemble_suffix !> Ensemble suffix provided by the cap. This may be + !! provided to bypass FMS ensemble manager. + + if (present(ensemble_suffix)) then + call fms_io_set_filename_appendix(trim(ensemble_suffix)) + else + call FMS_ensemble_manager_init() + endif end subroutine ensemble_manager_init diff --git a/config_src/infra/FMS2/MOM_ensemble_manager_infra.F90 b/config_src/infra/FMS2/MOM_ensemble_manager_infra.F90 index 66bbb86e2f..c9eb067e54 100644 --- a/config_src/infra/FMS2/MOM_ensemble_manager_infra.F90 +++ b/config_src/infra/FMS2/MOM_ensemble_manager_infra.F90 @@ -9,6 +9,8 @@ module MOM_ensemble_manager_infra use ensemble_manager_mod, only : FMS_get_ensemble_size => get_ensemble_size use ensemble_manager_mod, only : FMS_get_ensemble_pelist => get_ensemble_pelist use ensemble_manager_mod, only : FMS_get_ensemble_filter_pelist => get_ensemble_filter_pelist +use fms2_io_mod, only : fms2_io_set_filename_appendix=>set_filename_appendix +use fms_io_mod, only : fms_io_set_filename_appendix=>set_filename_appendix implicit none ; private @@ -20,9 +22,16 @@ module MOM_ensemble_manager_infra !> Initializes the ensemble manager which divides available resources !! in order to concurrently execute an ensemble of model realizations. -subroutine ensemble_manager_init() - - call FMS_ensemble_manager_init() +subroutine ensemble_manager_init(ensemble_suffix) + character(len=*), optional, intent(in) :: ensemble_suffix !> Ensemble suffix provided by the cap. This may be + !! provided to bypass FMS ensemble manager. + + if (present(ensemble_suffix)) then + call fms2_io_set_filename_appendix(trim(ensemble_suffix)) + call fms_io_set_filename_appendix(trim(ensemble_suffix)) + else + call FMS_ensemble_manager_init() + endif end subroutine ensemble_manager_init diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 84eb5fc90a..ba7152ea30 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1862,7 +1862,7 @@ end subroutine step_offline !! initializing the ocean state variables, and initializing subsidiary modules subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & Time_in, offline_tracer_mode, input_restart_file, diag_ptr, & - count_calls, tracer_flow_CSp, ice_shelf_CSp, waves_CSp) + count_calls, tracer_flow_CSp, ice_shelf_CSp, waves_CSp, ensemble_num) type(time_type), target, intent(inout) :: Time !< model time, set in this routine type(time_type), intent(in) :: Time_init !< The start time for the coupled model's calendar type(param_file_type), intent(out) :: param_file !< structure indicating parameter file to parse @@ -1885,7 +1885,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & !! dynamics timesteps. type(ice_shelf_CS), optional, pointer :: ice_shelf_CSp !< A pointer to an ice shelf control structure type(Wave_parameters_CS), & - optional, pointer :: Waves_CSp !< An optional pointer to a wave property CS + optional, pointer :: Waves_CSp !< An optional pointer to a wave property CS + integer, optional :: ensemble_num !< Ensemble index provided by the cap (instead of FMS + !! ensemble manager) ! local variables type(ocean_grid_type), pointer :: G => NULL() ! A pointer to the metric grid use for the run type(ocean_grid_type), pointer :: G_in => NULL() ! Pointer to the input grid @@ -1994,7 +1996,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! Read paths and filenames from namelist and store in "dirs". ! Also open the parsed input parameter file(s) and setup param_file. - call get_MOM_input(param_file, dirs, default_input_filename=input_restart_file) + call get_MOM_input(param_file, dirs, default_input_filename=input_restart_file, ensemble_num=ensemble_num) verbosity = 2 ; call read_param(param_file, "VERBOSITY", verbosity) call MOM_set_verbosity(verbosity) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 9a5e1f48f5..a3b7d604dd 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -191,10 +191,8 @@ module MOM_forcing_type real :: C_p !< heat capacity of seawater [Q C-1 ~> J kg-1 degC-1]. !! C_p is is the same value as in thermovar_ptrs_type. - ! CFC-related arrays needed in the MOM_CFC_cap module + ! arrays needed in the some tracer modules, e.g., MOM_CFC_cap real, pointer, dimension(:,:) :: & - cfc11_flux => NULL(), & !< flux of cfc_11 into the ocean [CU R Z T-1 ~> mol m-2 s-1] - cfc12_flux => NULL(), & !< flux of cfc_12 into the ocean [CU R Z T-1 ~> mol m-2 s-1] ice_fraction => NULL(), & !< fraction of sea ice coverage at h-cells, from 0 to 1 [nondim]. u10_sqr => NULL() !< wind magnitude at 10 m squared [L2 T-2 ~> m2 s-2] @@ -364,9 +362,7 @@ module MOM_forcing_type integer :: id_TKE_tidal = -1 integer :: id_buoy = -1 - ! cfc-related diagnostics handles - integer :: id_cfc11 = -1 - integer :: id_cfc12 = -1 + ! tracer surface flux related diagnostics handles integer :: id_ice_fraction = -1 integer :: id_u10_sqr = -1 @@ -1129,10 +1125,6 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) call hchksum(fluxes%u10_sqr, mesg//" fluxes%u10_sqr", G%HI, haloshift=hshift, scale=US%L_to_m**2*US%s_to_T**2) if (associated(fluxes%ice_fraction)) & call hchksum(fluxes%ice_fraction, mesg//" fluxes%ice_fraction", G%HI, haloshift=hshift) - if (associated(fluxes%cfc11_flux)) & - call hchksum(fluxes%cfc11_flux, mesg//" fluxes%cfc11_flux", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) - if (associated(fluxes%cfc12_flux)) & - call hchksum(fluxes%cfc12_flux, mesg//" fluxes%cfc12_flux", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) if (associated(fluxes%salt_flux)) & call hchksum(fluxes%salt_flux, mesg//" fluxes%salt_flux", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) if (associated(fluxes%TKE_tidal)) & @@ -1340,26 +1332,9 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, endif endif - ! units for cfc11_flux and cfc12_flux are [Conc R Z T-1 ~> mol m-2 s-1] ! See: - ! http://clipc-services.ceda.ac.uk/dreq/u/0940cbee6105037e4b7aa5579004f124.html - ! http://clipc-services.ceda.ac.uk/dreq/u/e9e21426e4810d0bb2d3dddb24dbf4dc.html if (present(use_cfcs)) then if (use_cfcs) then - handles%id_cfc11 = register_diag_field('ocean_model', 'cfc11_flux', diag%axesT1, Time, & - 'Gas exchange flux of CFC11 into the ocean ', & - 'mol m-2 s-1', conversion=US%RZ_T_to_kg_m2s, & - cmor_field_name='fgcfc11', & - cmor_long_name='Surface Downward CFC11 Flux', & - cmor_standard_name='surface_downward_cfc11_flux') - - handles%id_cfc12 = register_diag_field('ocean_model', 'cfc12_flux', diag%axesT1, Time, & - 'Gas exchange flux of CFC12 into the ocean ', & - 'mol m-2 s-1', conversion=US%RZ_T_to_kg_m2s, & - cmor_field_name='fgcfc12', & - cmor_long_name='Surface Downward CFC12 Flux', & - cmor_standard_name='surface_downward_cfc12_flux') - handles%id_ice_fraction = register_diag_field('ocean_model', 'ice_fraction', diag%axesT1, Time, & 'Fraction of cell area covered by sea ice', 'm2 m-2') @@ -2921,13 +2896,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (handles%id_netFWGlobalScl > 0) & call post_data(handles%id_netFWGlobalScl, fluxes%netFWGlobalScl, diag) - ! post diagnostics related to cfcs ==================================== - - if ((handles%id_cfc11 > 0) .and. associated(fluxes%cfc11_flux)) & - call post_data(handles%id_cfc11, fluxes%cfc11_flux, diag) - - if ((handles%id_cfc11 > 0) .and. associated(fluxes%cfc12_flux)) & - call post_data(handles%id_cfc12, fluxes%cfc12_flux, diag) + ! post diagnostics related to tracer surface fluxes ======================== if ((handles%id_ice_fraction > 0) .and. associated(fluxes%ice_fraction)) & call post_data(handles%id_ice_fraction, fluxes%ice_fraction, diag) @@ -2989,7 +2958,8 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & logical, optional, intent(in) :: salt !< If present and true, allocate salt fluxes logical, optional, intent(in) :: fix_accum_bug !< If present and true, avoid using a bug in !! accumulation of ustar_gustless - logical, optional, intent(in) :: cfc !< If present and true, allocate cfc fluxes + logical, optional, intent(in) :: cfc !< If present and true, allocate fields needed + !! for cfc surface fluxes logical, optional, intent(in) :: waves !< If present and true, allocate wave fields logical, optional, intent(in) :: shelf_sfc_accumulation !< If present and true, and shelf is true, !! then allocate surface flux deposition from the atmosphere @@ -3064,8 +3034,6 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%mass_berg,isd,ied,jsd,jed, iceberg) !These fields should only on allocated when USE_CFC_CAP is activated. - call myAlloc(fluxes%cfc11_flux,isd,ied,jsd,jed, cfc) - call myAlloc(fluxes%cfc12_flux,isd,ied,jsd,jed, cfc) call myAlloc(fluxes%ice_fraction,isd,ied,jsd,jed, cfc) call myAlloc(fluxes%u10_sqr,isd,ied,jsd,jed, cfc) @@ -3322,8 +3290,6 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%mass_berg)) deallocate(fluxes%mass_berg) if (associated(fluxes%ice_fraction)) deallocate(fluxes%ice_fraction) if (associated(fluxes%u10_sqr)) deallocate(fluxes%u10_sqr) - if (associated(fluxes%cfc11_flux)) deallocate(fluxes%cfc11_flux) - if (associated(fluxes%cfc12_flux)) deallocate(fluxes%cfc12_flux) call coupler_type_destructor(fluxes%tr_fluxes) diff --git a/src/core/MOM_unit_tests.F90 b/src/core/MOM_unit_tests.F90 index 10782e8890..d13be05ffd 100644 --- a/src/core/MOM_unit_tests.F90 +++ b/src/core/MOM_unit_tests.F90 @@ -9,7 +9,7 @@ module MOM_unit_tests use MOM_remapping, only : remapping_unit_tests use MOM_neutral_diffusion, only : neutral_diffusion_unit_tests use MOM_random, only : random_unit_tests -use MOM_lateral_boundary_diffusion, only : near_boundary_unit_tests +use MOM_hor_bnd_diffusion, only : near_boundary_unit_tests use MOM_CFC_cap, only : CFC_cap_unit_tests implicit none ; private diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 6aa94f584f..bf4b33af11 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -44,8 +44,6 @@ module MOM_variables SST, & !< The sea surface temperature [C ~> degC]. SSS, & !< The sea surface salinity [S ~> psu or gSalt/kg]. sfc_density, & !< The mixed layer density [R ~> kg m-3]. - sfc_cfc11, & !< Sea surface concentration of CFC11 [mol kg-1]. - sfc_cfc12, & !< Sea surface concentration of CFC12 [mol kg-1]. Hml, & !< The mixed layer depth [Z ~> m]. u, & !< The mixed layer zonal velocity [L T-1 ~> m s-1]. v, & !< The mixed layer meridional velocity [L T-1 ~> m s-1]. @@ -332,7 +330,7 @@ module MOM_variables !! the ocean model. Unused fields are unallocated. subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & gas_fields_ocn, use_meltpot, use_iceshelves, & - omit_frazil, use_cfcs) + omit_frazil) type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(surface), intent(inout) :: sfc_state !< ocean surface state type to be allocated. logical, optional, intent(in) :: use_temperature !< If true, allocate the space for thermodynamic variables. @@ -345,14 +343,13 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & !! tracer fluxes, and can be used to spawn related !! internal variables in the ice model. logical, optional, intent(in) :: use_meltpot !< If true, allocate the space for melt potential - logical, optional, intent(in) :: use_cfcs !< If true, allocate the space for cfcs logical, optional, intent(in) :: use_iceshelves !< If true, allocate the space for the stresses !! under ice shelves. logical, optional, intent(in) :: omit_frazil !< If present and false, do not allocate the space to !! pass frazil fluxes to the coupler ! local variables - logical :: use_temp, alloc_integ, use_melt_potential, alloc_iceshelves, alloc_frazil, alloc_cfcs + logical :: use_temp, alloc_integ, use_melt_potential, alloc_iceshelves, alloc_frazil integer :: is, ie, js, je, isd, ied, jsd, jed integer :: isdB, iedB, jsdB, jedB @@ -363,7 +360,6 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & use_temp = .true. ; if (present(use_temperature)) use_temp = use_temperature alloc_integ = .true. ; if (present(do_integrals)) alloc_integ = do_integrals use_melt_potential = .false. ; if (present(use_meltpot)) use_melt_potential = use_meltpot - alloc_cfcs = .false. ; if (present(use_cfcs)) alloc_cfcs = use_cfcs alloc_iceshelves = .false. ; if (present(use_iceshelves)) alloc_iceshelves = use_iceshelves alloc_frazil = .true. ; if (present(omit_frazil)) alloc_frazil = .not.omit_frazil @@ -387,11 +383,6 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & allocate(sfc_state%melt_potential(isd:ied,jsd:jed), source=0.0) endif - if (alloc_cfcs) then - allocate(sfc_state%sfc_cfc11(isd:ied,jsd:jed), source=0.0) - allocate(sfc_state%sfc_cfc12(isd:ied,jsd:jed), source=0.0) - endif - if (alloc_integ) then ! Allocate structures for the vertically integrated ocean_mass, ocean_heat, and ocean_salt. allocate(sfc_state%ocean_mass(isd:ied,jsd:jed), source=0.0) @@ -431,8 +422,6 @@ subroutine deallocate_surface_state(sfc_state) if (allocated(sfc_state%ocean_mass)) deallocate(sfc_state%ocean_mass) if (allocated(sfc_state%ocean_heat)) deallocate(sfc_state%ocean_heat) if (allocated(sfc_state%ocean_salt)) deallocate(sfc_state%ocean_salt) - if (allocated(sfc_state%sfc_cfc11)) deallocate(sfc_state%sfc_cfc11) - if (allocated(sfc_state%sfc_cfc12)) deallocate(sfc_state%sfc_cfc12) call coupler_type_destructor(sfc_state%tr_fields) sfc_state%arrays_allocated = .false. diff --git a/src/diagnostics/MOM_obsolete_params.F90 b/src/diagnostics/MOM_obsolete_params.F90 index 7564137de8..1f1a8e0d36 100644 --- a/src/diagnostics/MOM_obsolete_params.F90 +++ b/src/diagnostics/MOM_obsolete_params.F90 @@ -94,6 +94,8 @@ subroutine find_obsolete_params(param_file) call obsolete_real(param_file, "FIRST_GUESS_SURFACE_LAYER_DEPTH") call obsolete_logical(param_file, "CORRECT_SURFACE_LAYER_AVERAGE") call obsolete_int(param_file, "SEAMOUNT_LENGTH_SCALE", hint="Use SEAMOUNT_X_LENGTH_SCALE instead.") + call obsolete_int(param_file, "USE_LATERAL_BOUNDARY_DIFFUSION", & + hint="Use USE_HORIZONTAL_BOUNDARY_DIFFUSION instead.") call obsolete_logical(param_file, "MSTAR_FIXED", hint="Instead use MSTAR_MODE.") call obsolete_logical(param_file, "USE_VISBECK_SLOPE_BUG", .false.) diff --git a/src/framework/MOM_file_parser.F90 b/src/framework/MOM_file_parser.F90 index fd447f5193..5d658c44a4 100644 --- a/src/framework/MOM_file_parser.F90 +++ b/src/framework/MOM_file_parser.F90 @@ -124,7 +124,7 @@ module MOM_file_parser contains !> Make the contents of a parameter input file availalble in a param_file_type -subroutine open_param_file(filename, CS, checkable, component, doc_file_dir) +subroutine open_param_file(filename, CS, checkable, component, doc_file_dir, ensemble_num) character(len=*), intent(in) :: filename !< An input file name, optionally with the full path type(param_file_type), intent(inout) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters @@ -134,11 +134,13 @@ subroutine open_param_file(filename, CS, checkable, component, doc_file_dir) !! to generate parameter documentation file names; the default is"MOM" character(len=*), optional, intent(in) :: doc_file_dir !< An optional directory in which to write out !! the documentation files. The default is effectively './'. + integer, optional, intent(in) :: ensemble_num !< ensemble number to be appended to _doc filenames (optional) ! Local variables logical :: file_exists, Netcdf_file, may_check, reopened_file integer :: ios, iounit, strlen, i character(len=240) :: doc_path + character(len=5) :: ensemble_suffix type(parameter_block), pointer :: block => NULL() may_check = .true. ; if (present(checkable)) may_check = checkable @@ -211,6 +213,11 @@ subroutine open_param_file(filename, CS, checkable, component, doc_file_dir) call read_param(CS,"REPORT_UNUSED_PARAMS",CS%report_unused) call read_param(CS,"FATAL_UNUSED_PARAMS",CS%unused_params_fatal) CS%doc_file = "MOM_parameter_doc" + if (present(ensemble_num)) then + ! append instance suffix to doc_file + write(ensemble_suffix,'(A,I0.4)') '_', ensemble_num + CS%doc_file = trim(CS%doc_file)//ensemble_suffix + endif if (present(component)) CS%doc_file = trim(component)//"_parameter_doc" call read_param(CS,"DOCUMENT_FILE", CS%doc_file) if (.not.may_check) then diff --git a/src/framework/MOM_get_input.F90 b/src/framework/MOM_get_input.F90 index b6b5b89be9..4c643a5442 100644 --- a/src/framework/MOM_get_input.F90 +++ b/src/framework/MOM_get_input.F90 @@ -102,7 +102,7 @@ subroutine get_MOM_input(param_file, dirs, check_params, default_input_filename, if (len_trim(trim(parameter_filename(io))) > 0) then if (present(ensemble_num)) then call open_param_file(ensembler(parameter_filename(io),ensemble_num), param_file, & - check_params, doc_file_dir=output_dir) + check_params, doc_file_dir=output_dir, ensemble_num=ensemble_num) else call open_param_file(ensembler(parameter_filename(io)), param_file, & check_params, doc_file_dir=output_dir) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 old mode 100644 new mode 100755 index e6dd131a99..2af86488eb --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -68,6 +68,11 @@ module MOM_hor_visc logical :: use_beta_in_Leith !< If true, includes the beta term in the Leith viscosity logical :: Leith_Ah !< If true, use a biharmonic form of 2D Leith !! nonlinear eddy viscosity. AH is the background. + logical :: use_Leithy !< If true, use a biharmonic form of 2D Leith + !! nonlinear eddy viscosity with harmonic backscatter. + !! Ah is the background. Leithy = Leith+E + real :: c_K !< Fraction of energy dissipated by the biharmonic term + !! that gets backscattered in the Leith+E scheme. [nondim] logical :: use_QG_Leith_visc !< If true, use QG Leith nonlinear eddy viscosity. !! KH is the background value. logical :: bound_Coriolis !< If true & SMAGORINSKY_AH is used, the biharmonic @@ -145,10 +150,12 @@ module MOM_hor_visc n1n1_m_n2n2_q !< Factor n1**2-n2**2 in the anisotropic direction tensor at q-points [nondim] real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & - dx2h, & !< Pre-calculated dx^2 at h points [L2 ~> m2] - dy2h, & !< Pre-calculated dy^2 at h points [L2 ~> m2] - dx_dyT, & !< Pre-calculated dx/dy at h points [nondim] - dy_dxT !< Pre-calculated dy/dx at h points [nondim] + dx2h, & !< Pre-calculated dx^2 at h points [L2 ~> m2] + dy2h, & !< Pre-calculated dy^2 at h points [L2 ~> m2] + dx_dyT, & !< Pre-calculated dx/dy at h points [nondim] + dy_dxT, & !< Pre-calculated dy/dx at h points [nondim] + m_const_leithy, & !< Pre-calculated .5*sqrt(c_K)*max{dx,dy} [L ~> m] + m_leithy_max !< Pre-calculated 4./max(dx,dy)^2 at h points [L-2 ~> m-2] real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & dx2q, & !< Pre-calculated dx^2 at q points [L2 ~> m2] dy2q, & !< Pre-calculated dy^2 at q points [L2 ~> m2] @@ -257,18 +264,23 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, Del2u, & ! The u-component of the Laplacian of velocity [L-1 T-1 ~> m-1 s-1] h_u, & ! Thickness interpolated to u points [H ~> m or kg m-2]. vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1] + vort_xy_dy_smooth, & ! y-derivative of smoothed vertical vorticity [L-1 T-1 ~> m-1 s-1] div_xx_dx, & ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1] - ubtav ! zonal barotropic velocity averaged over a baroclinic time-step [L T-1 ~> m s-1] + ubtav, & ! zonal barotropic velocity averaged over a baroclinic time-step [L T-1 ~> m s-1] + u_smooth ! Zonal velocity, smoothed with a spatial low-pass filter [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G)) :: & Del2v, & ! The v-component of the Laplacian of velocity [L-1 T-1 ~> m-1 s-1] h_v, & ! Thickness interpolated to v points [H ~> m or kg m-2]. vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1] + vort_xy_dx_smooth, & ! x-derivative of smoothed vertical vorticity [L-1 T-1 ~> m-1 s-1] div_xx_dy, & ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1] - vbtav ! meridional barotropic velocity averaged over a baroclinic time-step [L T-1 ~> m s-1] + vbtav, & ! meridional barotropic velocity averaged over a baroclinic time-step [L T-1 ~> m s-1] + v_smooth ! Meridional velocity, smoothed with a spatial low-pass filter [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJ_(G)) :: & dudx_bt, dvdy_bt, & ! components in the barotropic horizontal tension [T-1 ~> s-1] div_xx, & ! Estimate of horizontal divergence at h-points [T-1 ~> s-1] sh_xx, & ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1] + sh_xx_smooth, & ! horizontal tension from smoothed velocity including metric terms [T-1 ~> s-1] sh_xx_bt, & ! barotropic horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1] str_xx,& ! str_xx is the diagonal term in the stress tensor [H L2 T-2 ~> m3 s-2 or kg s-2], but ! at some points in the code it is not yet layer integrated, so is in [L2 T-2 ~> m2 s-2]. @@ -279,23 +291,28 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] grad_div_mag_h, & ! Magnitude of divergence gradient at h-points [L-1 T-1 ~> m-1 s-1] dudx, dvdy, & ! components in the horizontal tension [T-1 ~> s-1] + dudx_smooth, dvdy_smooth, & ! components in the horizontal tension from smoothed velocity [T-1 ~> s-1] GME_effic_h, & ! The filtered efficiency of the GME terms at h points [nondim] - htot ! The total thickness of all layers [Z ~> m] + htot, & ! The total thickness of all layers [Z ~> m] + m_leithy ! Kh=m_leithy*Ah in Leith+E parameterization [L-2 ~> m-2] real :: Del2vort_h ! Laplacian of vorticity at h-points [L-2 T-1 ~> m-2 s-1] real :: grad_vel_mag_bt_h ! Magnitude of the barotropic velocity gradient tensor squared at h-points [T-2 ~> s-2] real :: boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim] real, dimension(SZIB_(G),SZJB_(G)) :: & dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1] + dvdx_smooth, dudy_smooth, & ! components in the shearing strain from smoothed velocity [T-1 ~> s-1] dDel2vdx, dDel2udy, & ! Components in the biharmonic equivalent of the shearing strain [L-2 T-1 ~> m-2 s-1] dvdx_bt, dudy_bt, & ! components in the barotropic shearing strain [T-1 ~> s-1] sh_xy, & ! horizontal shearing strain (du/dy + dv/dx) including metric terms [T-1 ~> s-1] + sh_xy_smooth, & ! horizontal shearing strain from smoothed velocity including metric terms [T-1 ~> s-1] sh_xy_bt, & ! barotropic horizontal shearing strain (du/dy + dv/dx) inc. metric terms [T-1 ~> s-1] str_xy, & ! str_xy is the cross term in the stress tensor [H L2 T-2 ~> m3 s-2 or kg s-2], but ! at some points in the code it is not yet layer integrated, so is in [L2 T-2 ~> m2 s-2]. str_xy_GME, & ! smoothed cross term in the stress tensor from GME [L2 T-2 ~> m2 s-2] bhstr_xy, & ! A copy of str_xy that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2] vort_xy, & ! Vertical vorticity (dv/dx - du/dy) including metric terms [T-1 ~> s-1] + vort_xy_smooth, & ! Vertical vorticity including metric terms, smoothed [T-1 ~> s-1] grad_vort_mag_q, & ! Magnitude of vorticity gradient at q-points [L-1 T-1 ~> m-1 s-1] grad_vort_mag_q_2d, & ! Magnitude of 2d vorticity gradient at q-points [L-1 T-1 ~> m-1 s-1] Del2vort_q, & ! Laplacian of vorticity at q-points [L-2 T-1 ~> m-2 s-1] @@ -331,6 +348,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, GME_coeff_h ! GME coefficient at h-points [L2 T-1 ~> m2 s-1] real :: AhSm ! Smagorinsky biharmonic viscosity [L4 T-1 ~> m4 s-1] real :: AhLth ! 2D Leith biharmonic viscosity [L4 T-1 ~> m4 s-1] + real :: AhLthy ! 2D Leith+E biharmonic viscosity [L4 T-1 ~> m4 s-1] real :: Shear_mag_bc ! Shear_mag value in backscatter [T-1 ~> s-1] real :: sh_xx_sq ! Square of tension (sh_xx) [T-2 ~> s-2] real :: sh_xy_sq ! Square of shearing strain (sh_xy) [T-2 ~> s-2] @@ -382,6 +400,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, Kh, & ! Laplacian viscosity (h or q) [L2 T-1 ~> m2 s-1] Shear_mag, & ! magnitude of the shear (h or q) [T-1 ~> s-1] vert_vort_mag, & ! magnitude of the vertical vorticity gradient (h or q) [L-1 T-1 ~> m-1 s-1] + vert_vort_mag_smooth, & ! magnitude of gradient of smoothed vertical vorticity (h or q) [L-1 T-1 ~> m-1 s-1] hrat_min, & ! h_min divided by the thickness at the stress point (h or q) [nondim] visc_bound_rem ! fraction of overall viscous bounds that remain to be applied (h or q) [nondim] @@ -394,6 +413,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, inv_PI2 = 1.0/((4.0*atan(1.0))**2) inv_PI6 = inv_PI3 * inv_PI3 + m_leithy(:,:) = 0. ! Initialize + if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then apply_OBC = OBC%Flather_u_BCs_exist_globally .or. OBC%Flather_v_BCs_exist_globally apply_OBC = .true. @@ -546,7 +567,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP grid_Ah, grid_Kh, d_Del2u, d_Del2v, d_str, & !$OMP Kh, Ah, AhSm, AhLth, local_strain, Sh_F_pow, & !$OMP dDel2vdx, dDel2udy, Del2vort_q, Del2vort_h, KE, & - !$OMP h2uq, h2vq, hu, hv, hq, FatH, RoScl, GME_coeff & + !$OMP h2uq, h2vq, hu, hv, hq, FatH, RoScl, GME_coeff, & + !$OMP dudx_smooth, dudy_smooth, dvdx_smooth, dvdy_smooth, & + !$OMP vort_xy_smooth, vort_xy_dx_smooth, vort_xy_dy_smooth, & + !$OMP sh_xx_smooth, sh_xy_smooth, u_smooth, v_smooth, & + !$OMP vert_vort_mag_smooth, m_leithy, AhLthy & !$OMP ) do k=1,nz @@ -575,6 +600,30 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, dudy(I,J) = CS%DX_dyBu(I,J)*(u(I,j+1,k)*G%IdxCu(I,j+1) - u(I,j,k)*G%IdxCu(I,j)) enddo ; enddo + if (CS%use_Leithy) then + ! Smooth the velocity. Right now it happens twice. In the future + ! one might make the number of smoothing cycles a user-specified parameter + u_smooth(:,:) = u(:,:,k) + v_smooth(:,:) = v(:,:,k) + call smooth_x9(CS, G, field_u=u_smooth,field_v=v_smooth) ! one call applies the filter twice + ! Calculate horizontal tension from smoothed velocity + do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 + dudx_smooth(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * u_smooth(I,j) - & + G%IdyCu(I-1,j) * u_smooth(I-1,j)) + dvdy_smooth(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * v_smooth(i,J) - & + G%IdxCv(i,J-1) * v_smooth(i,J-1)) + sh_xx_smooth(i,j) = dudx_smooth(i,j) - dvdy_smooth(i,j) + enddo ; enddo + + ! Components for the shearing strain from smoothed velocity + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 + dvdx_smooth(I,J) = CS%DY_dxBu(I,J) * & + (v_smooth(i+1,J)*G%IdyCv(i+1,J) - v_smooth(i,J)*G%IdyCv(i,J)) + dudy_smooth(I,J) = CS%DX_dyBu(I,J) * & + (u_smooth(I,j+1)*G%IdxCu(I,j+1) - u_smooth(I,j)*G%IdxCu(I,j)) + enddo ; enddo + end if ! use Leith+E + if (CS%id_normstress > 0) then do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 NoSt(i,j,k) = sh_xx(i,j) @@ -728,6 +777,20 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif + if (CS%use_Leithy) then + ! Shearing strain (including no-slip boundary conditions at the 2-D land-sea mask). + ! dudy_smooth and dvdx_smooth do not (yet) include modifications at OBCs from above. + if (CS%no_slip) then + do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + sh_xy_smooth(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_smooth(I,J) + dudy_smooth(I,J) ) + enddo ; enddo + else + do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + sh_xy_smooth(I,J) = G%mask2dBu(I,J) * ( dvdx_smooth(I,J) + dudy_smooth(I,J) ) + enddo ; enddo + endif + endif ! use Leith+E + ! Evaluate Del2u = x.Div(Grad u) and Del2v = y.Div( Grad u) if (CS%biharmonic) then do j=js-1,Jeq+1 ; do I=Isq-1,Ieq+1 @@ -765,12 +828,24 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif + if (CS%use_Leithy) then + if (CS%no_slip) then + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 + vort_xy_smooth(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_smooth(I,J) - dudy_smooth(I,J) ) + enddo ; enddo + else + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 + vort_xy_smooth(I,J) = G%mask2dBu(I,J) * ( dvdx_smooth(I,J) - dudy_smooth(I,J) ) + enddo ; enddo + endif + endif + ! Divergence do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 div_xx(i,j) = dudx(i,j) + dvdy(i,j) enddo ; enddo - if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then + if ((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then ! Vorticity gradient do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+2 @@ -783,6 +858,21 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, vort_xy_dy(I,j) = DX_dyBu * (vort_xy(I,J) * G%IdxCv(i,J) - vort_xy(I,J-1) * G%IdxCv(i,J-1)) enddo ; enddo + if (CS%use_Leithy) then + ! Gradient of smoothed vorticity + do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+2 + DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J) + vort_xy_dx_smooth(i,J) = DY_dxBu * & + (vort_xy_smooth(I,J) * G%IdyCu(I,j) - vort_xy_smooth(I-1,J) * G%IdyCu(I-1,j)) + enddo ; enddo + + do j=Jsq-1,Jeq+2 ; do I=Isq-1,Ieq+1 + DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J) + vort_xy_dy_smooth(I,j) = DX_dyBu * & + (vort_xy_smooth(I,J) * G%IdxCv(i,J) - vort_xy_smooth(I,J-1) * G%IdxCv(i,J-1)) + enddo ; enddo + endif ! If Leithy + ! Laplacian of vorticity do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J) @@ -865,6 +955,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, (0.5*(vort_xy_dy(I,j) + vort_xy_dy(I,j+1)))**2 ) enddo ; enddo + if (CS%use_Leithy) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + vert_vort_mag_smooth(i,j) = SQRT((0.5*(vort_xy_dx_smooth(i,J) + & + vort_xy_dx_smooth(i,J-1)))**2 + & + (0.5*(vort_xy_dy_smooth(I,j) + & + vort_xy_dy_smooth(I-1,j)))**2 ) + enddo ; enddo + endif ! Leithy + endif ! CS%Leith_Kh if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then @@ -890,6 +989,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif if (CS%Laplacian) then + ! Determine the Laplacian viscosity at h points, using the + ! largest value from several parameterizations. Also get + ! the Laplacian component of str_xx. + if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then if (CS%use_QG_Leith_visc) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -904,9 +1007,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif endif - ! Determine the Laplacian viscosity at h points, using the - ! largest value from several parameterizations. - ! Static (pre-computed) background viscosity do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 Kh(i,j) = CS%Kh_bg_xx(i,j) @@ -980,6 +1080,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif + ! In Leith+E parameterization Kh is computed after Ah in the biharmonic loop. + ! The harmonic component of str_xx is added in the biharmonic loop. + if (CS%use_Leithy) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Kh(i,j) = 0. + enddo ; enddo + end if + if (CS%id_Kh_h>0 .or. CS%debug) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 Kh_h(i,j,k) = Kh(i,j) @@ -1013,7 +1121,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 str_xx(i,j) = 0.0 enddo ; enddo - endif + endif ! Get Kh at h points and get Laplacian component of str_xx if (CS%anisotropic) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -1026,12 +1134,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%biharmonic) then ! Determine the biharmonic viscosity at h points, using the - ! largest value from several parameterizations. + ! largest value from several parameterizations. Also get the + ! biharmonic component of str_xx. do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 Ah(i,j) = CS%Ah_bg_xx(i,j) enddo ; enddo - if ((CS%Smagorinsky_Ah) .or. (CS%Leith_Ah)) then + if ((CS%Smagorinsky_Ah) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then if (CS%Smagorinsky_Ah) then if (CS%bound_Coriolis) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -1057,12 +1166,49 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif + if (CS%use_Leithy) then + ! Get m_leithy + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Del2vort_h = 0.25 * ((Del2vort_q(I,J) + Del2vort_q(I-1,J-1)) + & + (Del2vort_q(I-1,J) + Del2vort_q(I,J-1))) + AhLth = CS%Biharm6_const_xx(i,j) * inv_PI6 * abs(Del2vort_h) + if (AhLth .le. Ah(i,j)) then + m_leithy(i,j) = 0.0 + else + if ((CS%m_const_leithy(i,j)*vert_vort_mag(i,j)) .lt. abs(vort_xy_smooth(i,j))) then + m_leithy(i,j) = CS%c_K * (vert_vort_mag(i,j) / vort_xy_smooth(i,j))**2 + else + m_leithy(i,j) = CS%m_leithy_max(i,j) + endif + endif + enddo ; enddo + ! Smooth m_leithy + call smooth_x9(CS, G, field_h=m_leithy, zero_land=.true.) + ! Get Ah + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Del2vort_h = 0.25 * ((Del2vort_q(I,J) + Del2vort_q(I-1,J-1)) + & + (Del2vort_q(I-1,J) + Del2vort_q(I,J-1))) + AhLthy = CS%Biharm6_const_xx(i,j) * inv_PI6 * & + sqrt(max(0.,Del2vort_h**2 - m_leithy(i,j)*vert_vort_mag_smooth(i,j)**2)) + Ah(i,j) = max(Ah(i,j), AhLthy) + enddo ; enddo + ! Smooth Ah before applying upper bound + ! square, then smooth, then square root + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah_h(i,j,k) = Ah(i,j)**2 + enddo ; enddo + call smooth_x9(CS, G, field_h=Ah_h(:,:,k)) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah(i,j) = sqrt(Ah_h(i,j,k)) + enddo ; enddo + endif + if (CS%bound_Ah .and. .not. CS%better_bound_Ah) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 Ah(i,j) = min(Ah(i,j), CS%Ah_Max_xx(i,j)) enddo ; enddo endif - endif ! Smagorinsky_Ah or Leith_Ah + endif ! Smagorinsky_Ah or Leith_Ah or Leith+E if (use_MEKE_Au) then ! *Add* the MEKE contribution @@ -1096,6 +1242,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif + if (CS%use_Leithy) then + ! Compute Leith+E Kh after bounds have been applied to Ah + ! and after it has been smoothed. Kh = -m_leithy * Ah + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Kh(i,j) = -m_leithy(i,j) * Ah(i,j) + Kh_h(i,j,k) = Kh(i,j) + enddo ; enddo + endif + if (CS%id_grid_Re_Ah>0) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 KE = 0.125 * ((u(I,j,k) + u(I-1,j,k))**2 + (v(i,J,k) + v(i,J-1,k))**2) @@ -1111,10 +1266,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, str_xx(i,j) = str_xx(i,j) + d_str + if (CS%use_Leithy) str_xx(i,j) = str_xx(i,j) - Kh(i,j) * sh_xx_smooth(i,j) + ! Keep a copy of the biharmonic contribution for backscatter parameterization bhstr_xx(i,j) = d_str * (h(i,j,k) * CS%reduction_xx(i,j)) enddo ; enddo - endif + endif ! Get biharmonic coefficient at h points and biharmonic part of str_xx if (CS%biharmonic) then ! Gradient of Laplacian, for use in bi-harmonic term @@ -1203,6 +1360,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif if (CS%Laplacian) then + ! Determine the Laplacian viscosity at q points, using the + ! largest value from several parameterizations. Also get the + ! Laplacian component of str_xy. + if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then if (CS%use_QG_Leith_visc) then do J=js-1,Jeq ; do I=is-1,Ieq @@ -1217,9 +1378,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif endif - ! Determine the Laplacian viscosity at q points, using the - ! largest value from several parameterizations. - ! Static (pre-computed) background viscosity do J=js-1,Jeq ; do I=is-1,Ieq Kh(I,J) = CS%Kh_bg_xy(I,J) @@ -1286,6 +1444,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif endif + ! Leith+E doesn't recompute Kh at q points, it just interpolates it from h to q points + if (CS%use_Leithy) then + Kh(I,J) = Kh_h(i+1,j+1,k) + end if + if (CS%id_Kh_q>0 .or. CS%debug) & Kh_q(I,J,k) = Kh(I,J) @@ -1296,14 +1459,20 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, sh_xy_q(I,J,k) = sh_xy(I,J) enddo ; enddo - do J=js-1,Jeq ; do I=is-1,Ieq - str_xy(I,J) = -Kh(I,J) * sh_xy(I,J) - enddo ; enddo + if ( .not. CS%use_Leithy) then + do J=js-1,Jeq ; do I=is-1,Ieq + str_xy(I,J) = -Kh(I,J) * sh_xy(I,J) + enddo ; enddo + else + do J=js-1,Jeq ; do I=is-1,Ieq + str_xy(I,J) = -Kh(I,J) * sh_xy_smooth(I,J) + enddo ; enddo + endif else do J=js-1,Jeq ; do I=is-1,Ieq str_xy(I,J) = 0. enddo ; enddo - endif + endif ! get harmonic coefficient Kh at q points and harmonic part of str_xy if (CS%anisotropic) then do J=js-1,Jeq ; do I=is-1,Ieq @@ -1316,7 +1485,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%biharmonic) then ! Determine the biharmonic viscosity at q points, using the - ! largest value from several parameterizations. + ! largest value from several parameterizations. Also get the + ! biharmonic component of str_xy. do J=js-1,Jeq ; do I=is-1,Ieq Ah(I,J) = CS%Ah_bg_xy(I,J) enddo ; enddo @@ -1380,6 +1550,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif endif + ! Leith+E doesn't recompute Ah at q points, it just interpolates it from h to q points + if (CS%use_Leithy) then + do J=js-1,Jeq ; do I=is-1,Ieq + Ah(I,J) = Ah_h(i+1,j+1,k) + enddo ; enddo + end if + if (CS%id_Ah_q>0 .or. CS%debug) then do J=js-1,Jeq ; do I=is-1,Ieq Ah_q(I,J,k) = Ah(I,J) @@ -1395,7 +1572,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Keep a copy of the biharmonic contribution for backscatter parameterization bhstr_xy(I,J) = d_str * (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J)) enddo ; enddo - endif + endif ! Get Ah at q points and biharmonic part of str_xy if (CS%use_GME) then ! The wider halo here is to permit one pass of smoothing without a halo update. @@ -1907,6 +2084,10 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "If true, use a biharmonic Leith nonlinear eddy "//& "viscosity.", default=.false., do_not_log=.not.CS%biharmonic) if (.not.CS%biharmonic) CS%Leith_Ah = .false. + call get_param(param_file, mdl, "USE_LEITHY", CS%use_Leithy, & + "If true, use a biharmonic Leith nonlinear eddy "//& + "viscosity together with a harmonic backscatter.", & + default=.false., do_not_log=.not.(CS%biharmonic .and. CS%Laplacian)) call get_param(param_file, mdl, "BOUND_AH", CS%bound_Ah, & "If true, the biharmonic coefficient is locally limited "//& "to be stable.", default=.true., do_not_log=.not.CS%biharmonic) @@ -1965,12 +2146,11 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "Coriolis acceleration. The default is set by MAXVEL.", & units="m s-1", default=maxvel*US%L_T_to_m_s, scale=US%m_s_to_L_T, & do_not_log=.not.(CS%Smagorinsky_Ah .and. CS%bound_Coriolis)) - call get_param(param_file, mdl, "LEITH_BI_CONST", Leith_bi_const, & "The nondimensional biharmonic Leith constant, "//& "typical values are thus far undetermined.", units="nondim", default=0.0, & - fail_if_missing=CS%Leith_Ah, do_not_log=.not.CS%Leith_Ah) - + fail_if_missing=(CS%Leith_Ah .or. CS%use_Leithy), & + do_not_log=.not.(CS%Leith_Ah .or. CS%use_Leithy)) call get_param(param_file, mdl, "USE_LAND_MASK_FOR_HVISC", CS%use_land_mask, & "If true, use the land mask for the computation of thicknesses "//& "at velocity locations. This eliminates the dependence on arbitrary "//& @@ -2002,6 +2182,16 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "with the Gent and McWilliams parameterization.", default=.false.) call get_param(param_file, mdl, "SPLIT", split, & "Use the split time stepping if true.", default=.true., do_not_log=.true.) + if (CS%use_Leithy) then + if (.not.(CS%biharmonic .and. CS%Laplacian)) then + call MOM_error(FATAL, "MOM_hor_visc.F90, hor_visc_init:"//& + "LAPLACIAN and BIHARMONIC must both be True when USE_LEITHY=True.") + endif + call get_param(param_file, mdl, "LEITHY_CK", CS%c_K, & + "Fraction of biharmonic dissipation that gets backscattered, "//& + "in Leith+E.", units="nondim", default=1.0) + endif + if (CS%use_GME .and. .not.split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// & "cannot be used with SPLIT=False.") @@ -2120,9 +2310,13 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) ALLOC_(CS%Biharm_const2_xy(IsdB:IedB,JsdB:JedB)) ; CS%Biharm_const2_xy(:,:) = 0.0 endif endif - if (CS%Leith_Ah) then - ALLOC_(CS%biharm6_const_xx(isd:ied,jsd:jed)) ; CS%biharm6_const_xx(:,:) = 0.0 - ALLOC_(CS%biharm6_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%biharm6_const_xy(:,:) = 0.0 + if ((CS%Leith_Ah) .or. (CS%use_Leithy)) then + ALLOC_(CS%biharm6_const_xx(isd:ied,jsd:jed)) ; CS%biharm6_const_xx(:,:) = 0.0 + ALLOC_(CS%biharm6_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%biharm6_const_xy(:,:) = 0.0 + endif + if (CS%use_Leithy) then + ALLOC_(CS%m_const_leithy(isd:ied,jsd:jed)) ; CS%m_const_leithy(:,:) = 0.0 + ALLOC_(CS%m_leithy_max(isd:ied,jsd:jed)) ; CS%m_leithy_max(:,:) = 0.0 endif if (CS%Re_Ah > 0.0) then ALLOC_(CS%Re_Ah_const_xx(isd:ied,jsd:jed)); CS%Re_Ah_const_xx(:,:) = 0.0 @@ -2265,6 +2459,11 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) if (CS%Leith_Ah) then CS%biharm6_const_xx(i,j) = Leith_bi_const * (grid_sp_h3 * grid_sp_h3) endif + if (CS%use_Leithy) then + CS%biharm6_const_xx(i,j) = Leith_bi_const * max(G%dxT(i,j),G%dyT(i,j))**6 + CS%m_const_leithy(i,j) = 0.5 * sqrt(CS%c_K) * max(G%dxT(i,j),G%dyT(i,j)) + CS%m_leithy_max(i,j) = 4. / max(G%dxT(i,j),G%dyT(i,j))**2 + endif CS%Ah_bg_xx(i,j) = MAX(Ah, Ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2)) if (CS%Re_Ah > 0.0) CS%Re_Ah_const_xx(i,j) = grid_sp_h3 / CS%Re_Ah if (Ah_time_scale > 0.) CS%Ah_bg_xx(i,j) = & @@ -2287,7 +2486,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) (abs(G%CoriolisBu(I,J)) * BoundCorConst) endif endif - if (CS%Leith_Ah) then + if ((CS%Leith_Ah) .or. (CS%use_Leithy))then CS%biharm6_const_xy(I,J) = Leith_bi_const * (grid_sp_q3 * grid_sp_q3) endif CS%Ah_bg_xy(I,J) = MAX(Ah, Ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2)) @@ -2629,6 +2828,113 @@ subroutine smooth_GME(CS, G, GME_flux_h, GME_flux_q) enddo ! s-loop end subroutine smooth_GME +!> Apply a 9-point smoothing filter twice to reduce horizontal two-grid-point noise +!! Note that this subroutine does not conserve mass or angular momentum, so don't use it +!! in situations where you need conservation. Also can't apply it to Ah and Kh in the +!! horizontal_viscosity subroutine because they are not supposed to be halo-updated. +!! But you _can_ apply them to Kh_h and Ah_h. +subroutine smooth_x9(CS, G, field_h, field_u, field_v, field_q, zero_land) + type(hor_visc_CS), intent(in) :: CS !< Control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid + real, dimension(SZI_(G),SZJ_(G)), optional, intent(inout) :: field_h !< field to be smoothed + !! at h points + real, dimension(SZIB_(G),SZJ_(G)), optional, intent(inout) :: field_u !< field to be smoothed + !! at u points + real, dimension(SZI_(G),SZJB_(G)), optional, intent(inout) :: field_v !< field to be smoothed + !! at v points + real, dimension(SZIB_(G),SZJB_(G)), optional, intent(inout) :: field_q !< field to be smoothed + !! at q points + logical, optional, intent(in) :: zero_land !< An optional argument + !! indicating whether to set values + !! on land to zero (.true.) or + !! whether to ignore land values + !! (.false. or not present) + ! local variables. It would be good to make the _original variables allocatable. + real, dimension(SZI_(G),SZJ_(G)) :: field_h_original + real, dimension(SZIB_(G),SZJ_(G)) :: field_u_original + real, dimension(SZI_(G),SZJB_(G)) :: field_v_original + real, dimension(SZIB_(G),SZJB_(G)) :: field_q_original + real, dimension(3,3) :: weights, local_weights ! averaging weights for smoothing, nondimensional + logical :: zero_land_val ! actual value of zero_land optional argument + integer :: i, j, s + integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq + + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + + weights = reshape([1., 2., 1., 2., 4., 2., 1., 2., 1.],shape(weights))/16. + + if (present(zero_land)) then + zero_land_val = zero_land + else + zero_land_val = .false. + endif + + if (present(field_h)) then + call pass_var(field_h, G%Domain, halo=2) ! Halo size 2 ensures that you can smooth twice + do s=1,0,-1 + field_h_original(:,:) = field_h(:,:) + ! apply smoothing on field_h + do j=js-s,je+s ; do i=is-s,ie+s + ! skip land points + if (G%mask2dT(i,j)==0.) cycle + ! compute local weights + local_weights = weights*G%mask2dT(i-1:i+1,j-1:j+1) + if (zero_land_val) local_weights = local_weights/(sum(local_weights) + 1.E-16) + field_h(i,j) = sum(local_weights*field_h_original(i-1:i+1,j-1:j+1)) + enddo ; enddo + enddo + call pass_var(field_h, G%Domain) + endif + + if (present(field_u)) then + call pass_vector(field_u, field_v, G%Domain, halo=2) + do s=1,0,-1 + field_u_original(:,:) = field_u(:,:) + ! apply smoothing on field_u + do j=js-s,je+s ; do I=Isq-s,Ieq+s + ! skip land points + if (G%mask2dCu(I,j)==0.) cycle + ! compute local weights + local_weights = weights*G%mask2dCu(I-1:I+1,j-1:j+1) + if (zero_land_val) local_weights = local_weights/(sum(local_weights) + 1.E-16) + field_u(I,j) = sum(local_weights*field_u_original(I-1:I+1,j-1:j+1)) + enddo ; enddo + + field_v_original(:,:) = field_v(:,:) + ! apply smoothing on field_v + do J=Jsq-s,Jeq+s ; do i=is-s,ie+s + ! skip land points + if (G%mask2dCv(i,J)==0.) cycle + ! compute local weights + local_weights = weights*G%mask2dCv(i-1:i+1,J-1:J+1) + if (zero_land_val) local_weights = local_weights/(sum(local_weights) + 1.E-16) + field_v(i,J) = sum(local_weights*field_v_original(i-1:i+1,J-1:J+1)) + enddo ; enddo + enddo + call pass_vector(field_u, field_v, G%Domain) + endif + + if (present(field_q)) then + call pass_var(field_q, G%Domain, halo=2, position=CORNER) + do s=1,0,-1 + field_q_original(:,:) = field_q(:,:) + ! apply smoothing on field_q + do J=Jsq-s,Jeq+s ; do I=Isq-s,Ieq+s + ! skip land points + if (G%mask2dBu(I,J)==0.) cycle + ! compute local weights + local_weights = weights*G%mask2dBu(I-1:I+1,J-1:J+1) + if (zero_land_val) local_weights = local_weights/(sum(local_weights) + 1.E-16) + field_q(I,J) = sum(local_weights*field_q_original(I-1:I+1,J-1:J+1)) + enddo ; enddo + enddo + call pass_var(field_q, G%Domain, position=CORNER) + endif + +end subroutine smooth_x9 + !> Deallocates any variables allocated in hor_visc_init. subroutine hor_visc_end(CS) type(hor_visc_CS), intent(inout) :: CS !< Horizontal viscosity control structure @@ -2661,9 +2967,13 @@ subroutine hor_visc_end(CS) if (CS%Smagorinsky_Ah) then DEALLOC_(CS%Biharm_const_xx) ; DEALLOC_(CS%Biharm_const_xy) endif - if (CS%Leith_Ah) then + if ((CS%Leith_Ah) .or. (CS%use_Leithy)) then DEALLOC_(CS%Biharm6_const_xx) ; DEALLOC_(CS%Biharm6_const_xy) endif + if (CS%use_Leithy) then + DEALLOC_(CS%m_const_leithy) + DEALLOC_(CS%m_leithy_max) + endif if (CS%Re_Ah > 0.0) then DEALLOC_(CS%Re_Ah_const_xx) ; DEALLOC_(CS%Re_Ah_const_xy) endif diff --git a/src/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90 index 2a5e3f8854..4364dac0fd 100644 --- a/src/tracer/MOM_CFC_cap.F90 +++ b/src/tracer/MOM_CFC_cap.F90 @@ -5,6 +5,7 @@ module MOM_CFC_cap ! This file is part of MOM6. See LICENSE.md for the license. use MOM_coms, only : EFP_type +use MOM_debugging, only : hchksum use MOM_diag_mediator, only : diag_ctrl, register_diag_field, post_data use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_param, log_version, param_file_type @@ -18,14 +19,14 @@ module MOM_CFC_cap use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : query_initialized, set_initialized, MOM_restart_CS use MOM_spatial_means, only : global_mass_int_EFP -use MOM_time_manager, only : time_type +use MOM_time_manager, only : time_type, increment_date use time_interp_external_mod, only : init_external_field, time_interp_external use MOM_tracer_registry, only : register_tracer use MOM_tracer_types, only : tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_tracer_Z_init, only : tracer_Z_init use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : surface +use MOM_variables, only : surface, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -33,24 +34,29 @@ module MOM_CFC_cap #include public register_CFC_cap, initialize_CFC_cap, CFC_cap_unit_tests -public CFC_cap_column_physics, CFC_cap_surface_state, CFC_cap_fluxes +public CFC_cap_column_physics, CFC_cap_set_forcing public CFC_cap_stock, CFC_cap_end integer, parameter :: NTR = 2 !< the number of tracers in this module. -!> Contains the concentration array, a pointer to Tr in Tr_reg, and some metadata for a single CFC tracer +!> Contains the concentration array, surface flux, a pointer to Tr in Tr_reg, +!! and some metadata for a single CFC tracer type, private :: CFC_tracer_data - type(vardesc) :: desc !< A set of metadata for the tracer - real :: IC_val = 0.0 !< The initial value assigned to the tracer [mol kg-1]. - real :: land_val = -1.0 !< The value of the tracer used where land is masked out [mol kg-1]. - character(len=32) :: name !< Tracer variable name - integer :: id_cmor !< Diagnostic ID - real, pointer, dimension(:,:,:) :: conc !< The tracer concentration [mol kg-1]. - type(tracer_type), pointer :: tr_ptr !< pointer to tracer inside Tr_reg - end type CFC_tracer_data + type(vardesc) :: desc !< A set of metadata for the tracer + real :: IC_val = 0.0 !< The initial value assigned to the tracer [mol kg-1]. + real :: land_val = -1.0 !< The value of the tracer used where land is + !! masked out [mol kg-1]. + character(len=32) :: name !< Tracer variable name + integer :: id_cmor = -1 !< Diagnostic id + integer :: id_sfc_flux = -1 !< Surface flux id + real, pointer, dimension(:,:,:) :: conc !< The tracer concentration [mol kg-1]. + real, pointer, dimension(:,:) :: sfc_flux !< Surface flux [CU R Z T-1 ~> mol m-2 s-1] + type(tracer_type), pointer :: tr_ptr !< pointer to tracer inside Tr_reg +end type CFC_tracer_data !> The control structure for the CFC_cap tracer package type, public :: CFC_cap_CS ; private + logical :: debug !< If true, write verbose checksums for debugging purposes. character(len=200) :: IC_file !< The file in which the CFC initial values can !! be found, or an empty string for internal initilaization. logical :: Z_IC_file !< If true, the IC_file is in Z-space. The default is false. @@ -62,7 +68,12 @@ module MOM_CFC_cap !! the timing of diagnostic output. type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< Model restart control structure - type(CFC_tracer_data), dimension(2) :: CFC_data !< per-tracer parameters / metadata + type(CFC_tracer_data), dimension(NTR) :: CFC_data !< per-tracer parameters / metadata + integer :: CFC_BC_year_offset = 0 !< offset to add to model time to get time value used in CFC_BC_file + integer :: id_cfc11_atm_nh = -1 !< id number for time_interp_external. + integer :: id_cfc11_atm_sh = -1 !< id number for time_interp_external. + integer :: id_cfc12_atm_nh = -1 !< id number for time_interp_external. + integer :: id_cfc12_atm_sh = -1 !< id number for time_interp_external. end type CFC_cap_CS contains @@ -81,14 +92,17 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) ! Local variables character(len=40) :: mdl = "MOM_CFC_cap" ! This module's name. - character(len=200) :: inputdir ! The directory where NetCDF input files are. ! This include declares and sets the variable "version". # include "version_variable.h" + character(len=200) :: inputdir ! The directory where NetCDF input files are. real, dimension(:,:,:), pointer :: tr_ptr => NULL() - character(len=200) :: dummy ! Dummy variable to store params that need to be logged here. + character(len=200) :: CFC_BC_file ! filename with cfc11 and cfc12 data + character(len=30) :: CFC_BC_var_name ! varname of field in CFC_BC_file character :: m2char logical :: register_CFC_cap integer :: isd, ied, jsd, jed, nz, m + integer :: CFC_BC_data_year ! specific year in CFC BC data calendar + integer :: CFC_BC_model_year ! model year corresponding to CFC_BC_data_year isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke @@ -100,15 +114,19 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "DEBUG", CS%debug, & + "If true, write out verbose debugging data.", & + default=.false., debuggingParam=.true.) call get_param(param_file, mdl, "CFC_IC_FILE", CS%IC_file, & "The file in which the CFC initial values can be "//& "found, or an empty string for internal initialization.", & default=" ") - if ((len_trim(CS%IC_file) > 0) .and. (scan(CS%IC_file,'/') == 0)) then + if ((len_trim(CS%IC_file) > 0) .and. (scan(CS%IC_file, '/') == 0)) then ! Add the directory if CS%IC_file is not already a complete path. call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") CS%IC_file = trim(slasher(inputdir))//trim(CS%IC_file) - call log_param(param_file, mdl, "INPUTDIR/CFC_IC_FILE", CS%IC_file) + call log_param(param_file, mdl, "INPUTDIR/CFC_IC_FILE", CS%IC_file, & + "full path of CFC_IC_FILE") endif call get_param(param_file, mdl, "CFC_IC_FILE_IS_Z", CS%Z_IC_file, & "If true, CFC_IC_FILE is in depth space, not layer space", & @@ -118,7 +136,7 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) "if they are not found in the restart files. Otherwise "//& "it is a fatal error if tracers are not found in the "//& "restart files of a restarted run.", default=.false.) - do m=1,2 + do m=1,NTR write(m2char, "(I1)") m call get_param(param_file, mdl, "CFC1"//m2char//"_IC_VAL", CS%CFC_data(m)%IC_val, & "Value that CFC_1"//m2char//" is set to when it is not read from a file.", & @@ -127,28 +145,49 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) ! the following params are not used in this module. Instead, they are used in ! the cap but are logged here to keep all the CFC cap params together. - call get_param(param_file, mdl, "CFC_BC_FILE", dummy, & - "The file in which the CFC-11 and CFC-12 atm concentrations can be "//& - "found (units must be parts per trillion), or an empty string for "//& - "internal BC generation (TODO).", default=" ") - if ((len_trim(dummy) > 0) .and. (scan(dummy,'/') == 0)) then - ! Add the directory if dummy is not already a complete path. - call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") - dummy = trim(slasher(inputdir))//trim(dummy) - call log_param(param_file, mdl, "INPUTDIR/CFC_IC_FILE", dummy) + call get_param(param_file, mdl, "CFC_BC_FILE", CFC_BC_file, & + "The file in which the CFC-11 and CFC-12 atm concentrations can be "//& + "found (units must be parts per trillion).", default=" ") + if (len_trim(CFC_BC_file) == 0) then + call MOM_error(FATAL, "CFC_BC_FILE must be specified if USE_CFC_CAP=.true.") endif - if (len_trim(dummy) > 0) then - call get_param(param_file, mdl, "CFC11_VARIABLE", dummy, & - "The name of the variable representing CFC-11 in "//& - "CFC_BC_FILE.", default="CFC_11") - call get_param(param_file, mdl, "CFC12_VARIABLE", dummy, & - "The name of the variable representing CFC-12 in "//& - "CFC_BC_FILE.", default="CFC_12") + if (scan(CFC_BC_file, '/') == 0) then + ! Add the directory if CFC_BC_file is not already a complete path. + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + CFC_BC_file = trim(slasher(inputdir))//trim(CFC_BC_file) + call log_param(param_file, mdl, "INPUTDIR/CFC_BC_FILE", CFC_BC_file, & + "full path of CFC_BC_FILE") endif + call get_param(param_file, mdl, "CFC_BC_DATA_YEAR", CFC_BC_data_year, & + "Specific year in CFC_BC_FILE data calendar", default=2000) + call get_param(param_file, mdl, "CFC_BC_MODEL_YEAR", CFC_BC_model_year, & + "Model year corresponding to CFC_BC_MODEL_YEAR", default=2000) + CS%CFC_BC_year_offset = CFC_BC_data_year - CFC_BC_model_year + + call get_param(param_file, mdl, "CFC11_NH_VARIABLE", CFC_BC_var_name, & + "Variable name of NH CFC-11 atm mole fraction in CFC_BC_FILE.", & + default="cfc11_nh") + CS%id_cfc11_atm_nh = init_external_field(CFC_BC_file, CFC_BC_var_name) + + call get_param(param_file, mdl, "CFC11_SH_VARIABLE", CFC_BC_var_name, & + "Variable name of SH CFC-11 atm mole fraction in CFC_BC_FILE.", & + default="cfc11_sh") + CS%id_cfc11_atm_sh = init_external_field(CFC_BC_file, CFC_BC_var_name) + + call get_param(param_file, mdl, "CFC12_NH_VARIABLE", CFC_BC_var_name, & + "Variable name of NH CFC-12 atm mole fraction in CFC_BC_FILE.", & + default="cfc12_nh") + CS%id_cfc12_atm_nh = init_external_field(CFC_BC_file, CFC_BC_var_name) + + call get_param(param_file, mdl, "CFC12_SH_VARIABLE", CFC_BC_var_name, & + "Variable name of SH CFC-12 atm mole fraction in CFC_BC_FILE.", & + default="cfc12_sh") + CS%id_cfc12_atm_sh = init_external_field(CFC_BC_file, CFC_BC_var_name) + ! The following vardesc types contain a package of metadata about each tracer, ! including, the name; units; longname; and grid information. - do m=1,2 + do m=1,NTR write(m2char, "(I1)") m write(CS%CFC_data(m)%name, "(2A)") "CFC_1", m2char CS%CFC_data(m)%desc = var_desc(CS%CFC_data(m)%name, & @@ -157,6 +196,7 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) caller=mdl) allocate(CS%CFC_data(m)%conc(isd:ied,jsd:jed,nz), source=0.0) + allocate(CS%CFC_data(m)%sfc_flux(isd:ied,jsd:jed), source=0.0) ! This pointer assignment is needed to force the compiler not to do a copy in ! the registration calls. Curses on the designers and implementers of F90. @@ -201,7 +241,7 @@ subroutine initialize_CFC_cap(restart, day, G, GV, US, h, diag, OBC, CS) CS%Time => day CS%diag => diag - do m=1,2 + do m=1,NTR if (.not.restart .or. (CS%tracers_may_reinit .and. & .not.query_initialized(CS%CFC_data(m)%conc, CS%CFC_data(m)%name, CS%restart_CSp))) then call init_tracer_CFC(h, CS%CFC_data(m)%conc, CS%CFC_data(m)%name, CS%CFC_data(m)%land_val, & @@ -210,11 +250,23 @@ subroutine initialize_CFC_cap(restart, day, G, GV, US, h, diag, OBC, CS) endif ! cmor diagnostics + ! units for cfc11_flux and cfc12_flux are [Conc R Z T-1 ~> mol m-2 s-1] ! CFC11 cmor conventions: http://clipc-services.ceda.ac.uk/dreq/u/42625c97b8fe75124a345962c4430982.html + ! http://clipc-services.ceda.ac.uk/dreq/u/0940cbee6105037e4b7aa5579004f124.html ! CFC12 cmor conventions: http://clipc-services.ceda.ac.uk/dreq/u/3ab8e10027d7014f18f9391890369235.html + ! http://clipc-services.ceda.ac.uk/dreq/u/e9e21426e4810d0bb2d3dddb24dbf4dc.html write(m2char, "(I1)") m - CS%CFC_data(m)%id_cmor = register_diag_field('ocean_model', 'cfc1'//m2char, diag%axesTL, day, & - 'Mole Concentration of CFC1'//m2char//' in Sea Water', 'mol m-3') + CS%CFC_data(m)%id_cmor = register_diag_field('ocean_model', & + 'cfc1'//m2char, diag%axesTL, day, & + 'Mole Concentration of CFC1'//m2char//' in Sea Water', 'mol m-3') + + CS%CFC_data(m)%id_sfc_flux = register_diag_field('ocean_model', & + 'cfc1'//m2char//'_flux', diag%axesT1, day, & + 'Gas exchange flux of CFC1'//m2char//' into the ocean ', & + 'mol m-2 s-1', conversion=US%RZ_T_to_kg_m2s, & + cmor_field_name='fgcfc1'//m2char, & + cmor_long_name='Surface Downward CFC1'//m2char//' Flux', & + cmor_standard_name='surface_downward_cfc1'//m2char//'_flux') enddo @@ -308,7 +360,7 @@ subroutine CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, C ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2] real :: flux_scale - integer :: i, j, k, is, ie, js, je, nz + integer :: i, j, k, is, ie, js, je, nz, m is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -319,43 +371,43 @@ subroutine CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, C if (associated(KPP_CSp) .and. present(nonLocalTrans)) then flux_scale = GV%Z_to_H / GV%rho0 - call KPP_NonLocalTransport(KPP_CSp, G, GV, h_old, nonLocalTrans, fluxes%cfc11_flux(:,:), dt, CS%diag, & - CS%CFC_data(1)%tr_ptr, CS%CFC_data(1)%conc(:,:,:), & - flux_scale=flux_scale) - call KPP_NonLocalTransport(KPP_CSp, G, GV, h_old, nonLocalTrans, fluxes%cfc12_flux(:,:), dt, CS%diag, & - CS%CFC_data(2)%tr_ptr, CS%CFC_data(2)%conc(:,:,:), & - flux_scale=flux_scale) + do m=1,NTR + call KPP_NonLocalTransport(KPP_CSp, G, GV, h_old, nonLocalTrans, & + CS%CFC_data(m)%sfc_flux(:,:), dt, CS%diag, & + CS%CFC_data(m)%tr_ptr, CS%CFC_data(m)%conc(:,:,:), & + flux_scale=flux_scale) + enddo endif endif ! Use a tridiagonal solver to determine the concentrations after the ! surface source is applied and diapycnal advection and diffusion occurs. if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then - do k=1,nz ;do j=js,je ; do i=is,ie - h_work(i,j,k) = h_old(i,j,k) - enddo ; enddo ; enddo - call applyTracerBoundaryFluxesInOut(G, GV, CS%CFC_data(1)%conc, dt, fluxes, h_work, & - evap_CFL_limit, minimum_forcing_depth) - call tracer_vertdiff(h_work, ea, eb, dt, CS%CFC_data(1)%conc, G, GV, sfc_flux=fluxes%cfc11_flux) - - do k=1,nz ;do j=js,je ; do i=is,ie - h_work(i,j,k) = h_old(i,j,k) - enddo ; enddo ; enddo - call applyTracerBoundaryFluxesInOut(G, GV, CS%CFC_data(2)%conc, dt, fluxes, h_work, & - evap_CFL_limit, minimum_forcing_depth) - call tracer_vertdiff(h_work, ea, eb, dt, CS%CFC_data(2)%conc, G, GV, sfc_flux=fluxes%cfc12_flux) + do m=1,NTR + do k=1,nz ;do j=js,je ; do i=is,ie + h_work(i,j,k) = h_old(i,j,k) + enddo ; enddo ; enddo + call applyTracerBoundaryFluxesInOut(G, GV, CS%CFC_data(m)%conc, dt, fluxes, h_work, & + evap_CFL_limit, minimum_forcing_depth) + call tracer_vertdiff(h_work, ea, eb, dt, CS%CFC_data(m)%conc, G, GV, & + sfc_flux=CS%CFC_data(m)%sfc_flux) + enddo else - call tracer_vertdiff(h_old, ea, eb, dt, CS%CFC_data(1)%conc, G, GV, sfc_flux=fluxes%cfc11_flux) - call tracer_vertdiff(h_old, ea, eb, dt, CS%CFC_data(2)%conc, G, GV, sfc_flux=fluxes%cfc12_flux) + do m=1,NTR + call tracer_vertdiff(h_old, ea, eb, dt, CS%CFC_data(m)%conc, G, GV, & + sfc_flux=CS%CFC_data(m)%sfc_flux) + enddo endif ! If needed, write out any desired diagnostics from tracer sources & sinks here. - if (CS%CFC_data(1)%id_cmor > 0) call post_data(CS%CFC_data(1)%id_cmor, & - (GV%Rho0*US%R_to_kg_m3)*CS%CFC_data(1)%conc, & - CS%diag) - if (CS%CFC_data(2)%id_cmor > 0) call post_data(CS%CFC_data(2)%id_cmor, & - (GV%Rho0*US%R_to_kg_m3)*CS%CFC_data(2)%conc, & - CS%diag) + do m=1,NTR + if (CS%CFC_data(m)%id_cmor > 0) & + call post_data(CS%CFC_data(m)%id_cmor, & + (GV%Rho0*US%R_to_kg_m3)*CS%CFC_data(m)%conc, CS%diag) + + if (CS%CFC_data(m)%id_sfc_flux > 0) & + call post_data(CS%CFC_data(m)%id_sfc_flux, CS%CFC_data(m)%sfc_flux, CS%diag) + enddo end subroutine CFC_cap_column_physics @@ -394,95 +446,72 @@ function CFC_cap_stock(h, stocks, G, GV, CS, names, units, stock_index) return endif ; endif - do m=1,2 + do m=1,NTR call query_vardesc(CS%CFC_data(m)%desc, name=names(m), units=units(m), caller="CFC_cap_stock") units(m) = trim(units(m))//" kg" stocks(m) = global_mass_int_EFP(h, G, GV, CS%CFC_data(m)%conc, on_PE_only=.true.) enddo - CFC_cap_stock = 2 + CFC_cap_stock = NTR end function CFC_cap_stock -!> Extracts the ocean surface CFC concentrations and copies them to sfc_state. -subroutine CFC_cap_surface_state(sfc_state, G, CS) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(surface), intent(inout) :: sfc_state !< A structure containing fields that - !! describe the surface state of the ocean. - type(CFC_cap_CS), pointer :: CS!< The control structure returned by a previous - !! call to register_CFC_cap. - - ! Local variables - integer :: i, j, is, ie, js, je - - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - - if (.not.associated(CS)) return - - do j=js,je ; do i=is,ie - sfc_state%sfc_cfc11(i,j) = CS%CFC_data(1)%conc(i,j,1) - sfc_state%sfc_cfc12(i,j) = CS%CFC_data(2)%conc(i,j,1) - enddo ; enddo - -end subroutine CFC_cap_surface_state - !> Orchestrates the calculation of the CFC fluxes [mol m-2 s-1], including getting the ATM !! concentration, and calculating the solubility, Schmidt number, and gas exchange. -subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id_cfc12_atm) - type(ocean_grid_type), intent(in ) :: G !< The ocean's grid structure. - type(unit_scale_type), intent(in ) :: US !< A dimensional unit scaling type - type(surface), intent(in ) :: sfc_state !< A structure containing fields - !! that describe the surface state of the ocean. - type(forcing), intent(inout) :: fluxes !< A structure containing pointers - !! to thermodynamic and tracer forcing fields. Unused fields - !! have NULL ptrs. - real, intent(in ) :: Rho0 !< The mean ocean density [R ~> kg m-3] - type(time_type), intent(in ) :: Time !< The time of the fluxes, used for interpolating the - !! CFC's concentration in the atmosphere. - integer, optional, intent(inout):: id_cfc11_atm !< id number for time_interp_external. - integer, optional, intent(inout):: id_cfc12_atm !< id number for time_interp_external. +subroutine CFC_cap_set_forcing(sfc_state, fluxes, day_start, day_interval, G, US, Rho0, CS) + type(surface), intent(in ) :: sfc_state !< A structure containing fields + !! that describe the surface state of the ocean. + type(forcing), intent(inout) :: fluxes !< A structure containing pointers + !! to thermodynamic and tracer forcing fields. Unused fields + !! have NULL ptrs. + type(time_type), intent(in) :: day_start !< Start time of the fluxes. + type(time_type), intent(in) :: day_interval !< Length of time over which these + !! fluxes will be applied. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, intent(in) :: Rho0 !< The mean ocean density [R ~> kg m-3] + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. ! Local variables + type(time_type) :: Time_external ! time value used in CFC_BC_file real, dimension(SZI_(G),SZJ_(G)) :: & kw_wo_sc_no_term, & ! gas transfer velocity, without the Schmidt number term [Z T-1 ~> m s-1]. - kw, & ! gas transfer velocity [Z T-1 ~> m s-1]. - cair, & ! The surface gas concentration in equilibrium with the atmosphere (saturation concentration) - ! [mol kg-1]. - cfc11_atm, & !< CFC11 concentration in the atmopshere [pico mol/mol] - cfc12_atm !< CFC11 concentration in the atmopshere [pico mol/mol] - real :: ta ! Absolute sea surface temperature [hectoKelvin] - real :: sal ! Surface salinity [PSU]. - real :: alpha_11 ! The solubility of CFC 11 [mol kg-1 atm-1]. - real :: alpha_12 ! The solubility of CFC 12 [mol kg-1 atm-1]. - real :: sc_11, sc_12 ! The Schmidt numbers of CFC 11 and CFC 12 [nondim]. - real :: kw_coeff ! A coefficient used to compute the piston velocity [Z T-1 T2 L-2 = Z T L-2 ~> s / m] + kw, & ! gas transfer velocity [Z T-1 ~> m s-1]. + cair, & ! The surface gas concentration in equilibrium with the atmosphere + ! (saturation concentration) [mol kg-1]. + cfc11_atm, & ! CFC11 atm mole fraction [pico mol/mol] + cfc12_atm ! CFC12 atm mole fraction [pico mol/mol] + real :: cfc11_atm_nh ! NH value for cfc11_atm + real :: cfc11_atm_sh ! SH value for cfc11_atm + real :: cfc12_atm_nh ! NH value for cfc12_atm + real :: cfc12_atm_sh ! SH value for cfc12_atm + real :: ta ! Absolute sea surface temperature [hectoKelvin] + real :: sal ! Surface salinity [PSU]. + real :: alpha_11 ! The solubility of CFC 11 [mol kg-1 atm-1]. + real :: alpha_12 ! The solubility of CFC 12 [mol kg-1 atm-1]. + real :: sc_11, sc_12 ! The Schmidt numbers of CFC 11 and CFC 12 [nondim]. + real :: kw_coeff ! A coefficient used to compute the piston velocity [Z T-1 T2 L-2 = Z T L-2 ~> s / m] real, parameter :: pa_to_atm = 9.8692316931427e-6 ! factor for converting from Pa to atm [atm Pa-1]. - real :: press_to_atm ! converts from model pressure units to atm [atm T2 R-1 L-2 ~> atm Pa-1] - integer :: i, j, is, ie, js, je + real :: press_to_atm ! converts from model pressure units to atm [atm T2 R-1 L-2 ~> atm Pa-1] + integer :: i, j, is, ie, js, je, m is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - ! CFC11 ATM concentration - if (present(id_cfc11_atm) .and. (id_cfc11_atm /= -1)) then - call time_interp_external(id_cfc11_atm, Time, cfc11_atm) - ! convert from ppt (pico mol/mol) to mol/mol - cfc11_atm = cfc11_atm * 1.0e-12 - else - ! TODO: create cfc11_atm internally - call MOM_error(FATAL, "CFC_cap_fluxes: option to create cfc11_atm internally" //& - "has not been implemented yet.") - endif + ! Time_external = increment_date(day_start + day_interval/2, years=CS%CFC_BC_year_offset) + Time_external = increment_date(day_start, years=CS%CFC_BC_year_offset) - ! CFC12 ATM concentration - if (present(id_cfc12_atm) .and. (id_cfc12_atm /= -1)) then - call time_interp_external(id_cfc12_atm, Time, cfc12_atm) - ! convert from ppt (pico mol/mol) to mol/mol - cfc12_atm = cfc12_atm * 1.0e-12 - else - ! TODO: create cfc11_atm internally - call MOM_error(FATAL, "CFC_cap_fluxes: option to create cfc12_atm internally" //& - "has not been implemented yet.") - endif + ! CFC11 atm mole fraction, convert from ppt (pico mol/mol) to mol/mol + call time_interp_external(CS%id_cfc11_atm_nh, Time_external, cfc11_atm_nh) + cfc11_atm_nh = cfc11_atm_nh * 1.0e-12 + call time_interp_external(CS%id_cfc11_atm_sh, Time_external, cfc11_atm_sh) + cfc11_atm_sh = cfc11_atm_sh * 1.0e-12 + + ! CFC12 atm mole fraction, convert from ppt (pico mol/mol) to mol/mol + call time_interp_external(CS%id_cfc12_atm_nh, Time_external, cfc12_atm_nh) + cfc12_atm_nh = cfc12_atm_nh * 1.0e-12 + call time_interp_external(CS%id_cfc12_atm_sh, Time_external, cfc12_atm_sh) + cfc12_atm_sh = cfc12_atm_sh * 1.0e-12 !--------------------------------------------------------------------- ! Gas exchange/piston velocity parameter @@ -494,6 +523,21 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id ! set unit conversion factors press_to_atm = US%R_to_kg_m3*US%L_T_to_m_s**2 * pa_to_atm + do j=js,je ; do i=is,ie + if (G%geoLatT(i,j) < -10.0) then + cfc11_atm(i,j) = cfc11_atm_sh + cfc12_atm(i,j) = cfc12_atm_sh + elseif (G%geoLatT(i,j) <= 10.0) then + cfc11_atm(i,j) = cfc11_atm_sh + & + (0.05 * G%geoLatT(i,j) + 0.5) * (cfc11_atm_nh - cfc11_atm_sh) + cfc12_atm(i,j) = cfc12_atm_sh + & + (0.05 * G%geoLatT(i,j) + 0.5) * (cfc12_atm_nh - cfc12_atm_sh) + else + cfc11_atm(i,j) = cfc11_atm_nh + cfc12_atm(i,j) = cfc12_atm_nh + endif + enddo ; enddo + do j=js,je ; do i=is,ie ! ta in hectoKelvin ta = max(0.01, (US%C_to_degC*sfc_state%SST(i,j) + 273.15) * 0.01) @@ -512,14 +556,21 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id ! CFC flux units: CU R Z T-1 = mol kg-1 R Z T-1 ~> mol m-2 s-1 kw(i,j) = kw_wo_sc_no_term(i,j) * sqrt(660.0 / sc_11) cair(i,j) = press_to_atm * alpha_11 * cfc11_atm(i,j) * fluxes%p_surf_full(i,j) - fluxes%cfc11_flux(i,j) = kw(i,j) * (cair(i,j) - sfc_state%sfc_CFC11(i,j)) * Rho0 + CS%CFC_data(1)%sfc_flux(i,j) = kw(i,j) * (cair(i,j) - CS%CFC_data(1)%conc(i,j,1)) * Rho0 kw(i,j) = kw_wo_sc_no_term(i,j) * sqrt(660.0 / sc_12) cair(i,j) = press_to_atm * alpha_12 * cfc12_atm(i,j) * fluxes%p_surf_full(i,j) - fluxes%cfc12_flux(i,j) = kw(i,j) * (cair(i,j) - sfc_state%sfc_CFC12(i,j)) * Rho0 + CS%CFC_data(2)%sfc_flux(i,j) = kw(i,j) * (cair(i,j) - CS%CFC_data(2)%conc(i,j,1)) * Rho0 enddo ; enddo -end subroutine CFC_cap_fluxes + if (CS%debug) then + do m=1,NTR + call hchksum(CS%CFC_data(m)%sfc_flux, trim(CS%CFC_data(m)%name)//" sfc_flux", G%HI, & + scale=US%RZ_T_to_kg_m2s) + enddo + endif + +end subroutine CFC_cap_set_forcing !> Calculates the CFC's solubility function following Warner and Weiss (1985) DSR, vol 32. subroutine get_solubility(alpha_11, alpha_12, ta, sal , mask) @@ -606,8 +657,9 @@ subroutine CFC_cap_end(CS) integer :: m if (associated(CS)) then - do m=1,2 + do m=1,NTR if (associated(CS%CFC_data(m)%conc)) deallocate(CS%CFC_data(m)%conc) + if (associated(CS%CFC_data(m)%sfc_flux)) deallocate(CS%CFC_data(m)%sfc_flux) enddo deallocate(CS) @@ -684,7 +736,7 @@ logical function compare_values(verbose, test_name, calc, ans, limit) write(stdout,10) calc, ans endif -10 format("calc=",f20.16," ans",f20.16) +10 format("calc=",f22.16," ans",f22.16) end function compare_values !> \namespace mom_CFC_cap diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_hor_bnd_diffusion.F90 similarity index 80% rename from src/tracer/MOM_lateral_boundary_diffusion.F90 rename to src/tracer/MOM_hor_bnd_diffusion.F90 index f26395c119..b89552e8e4 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_hor_bnd_diffusion.F90 @@ -1,7 +1,7 @@ -!> Calculates and applies diffusive fluxes as a parameterization of lateral mixing (non-neutral) by +!> Calculates and applies diffusive fluxes as a parameterization of horizontal mixing (non-neutral) by !! mesoscale eddies near the top and bottom (to be implemented) boundary layers of the ocean. -module MOM_lateral_boundary_diffusion +module MOM_hor_bnd_diffusion ! This file is part of MOM6. See LICENSE.md for the license. @@ -28,18 +28,19 @@ module MOM_lateral_boundary_diffusion implicit none ; private -public near_boundary_unit_tests, lateral_boundary_diffusion, lateral_boundary_diffusion_init -public boundary_k_range +public near_boundary_unit_tests, hor_bnd_diffusion, hor_bnd_diffusion_init +public boundary_k_range, hor_bnd_diffusion_end ! Private parameters to avoid doing string comparisons for bottom or top boundary layer integer, public, parameter :: SURFACE = -1 !< Set a value that corresponds to the surface bopundary integer, public, parameter :: BOTTOM = 1 !< Set a value that corresponds to the bottom boundary #include -!> Sets parameters for lateral boundary mixing module. -type, public :: lbd_CS ; private +!> Sets parameters for horizontal boundary mixing module. +type, public :: hbd_CS ; private logical :: debug !< If true, write verbose checksums for debugging. integer :: deg !< Degree of polynomial reconstruction. + integer :: hbd_nk !< Maximum number of levels in the HBD grid [nondim] integer :: surface_boundary_scheme !< Which boundary layer scheme to use !! 1. ePBL; 2. KPP logical :: limiter !< Controls whether a flux limiter is applied in the @@ -51,50 +52,58 @@ module MOM_lateral_boundary_diffusion real :: H_subroundoff !< A thickness that is so small that it can be added to a thickness of !! Angstrom or larger without changing it at the bit level [H ~> m or kg m-2]. !! If Angstrom is 0 or exceedingly small, this is negligible compared to 1e-17 m. + ! HBD dynamic grids + real, allocatable, dimension(:,:,:) :: hbd_grd_u !< HBD thicknesses at t-points adjecent to + !! u-points [H ~> m or kg m-2] + real, allocatable, dimension(:,:,:) :: hbd_grd_v !< HBD thicknesses at t-points adjacent to + !! v-points (left and right) [H ~> m or kg m-2] + integer, allocatable, dimension(:,:) :: hbd_u_kmax !< Maximum vertical index in hbd_grd_u [nondim] + integer, allocatable, dimension(:,:) :: hbd_v_kmax !< Maximum vertical index in hbd_grd_v [nondim] type(remapping_CS) :: remap_CS !< Control structure to hold remapping configuration. type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to get BLD. type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure needed to get BLD. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. -end type lbd_CS +end type hbd_CS ! This include declares and sets the variable "version". #include "version_variable.h" -character(len=40) :: mdl = "MOM_lateral_boundary_diffusion" !< Name of this module -integer :: id_clock_lbd !< CPU clock for lbd +character(len=40) :: mdl = "MOM_hor_bnd_diffusion" !< Name of this module +integer :: id_clock_hbd !< CPU clock for hbd contains !> Initialization routine that reads runtime parameters and sets up pointers to other control structures that might be -!! needed for lateral boundary diffusion. -logical function lateral_boundary_diffusion_init(Time, G, GV, param_file, diag, diabatic_CSp, CS) +!! needed for horizontal boundary diffusion. +logical function hor_bnd_diffusion_init(Time, G, GV, US, param_file, diag, diabatic_CSp, CS) type(time_type), target, intent(in) :: Time !< Time structure type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< Parameter file structure type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure type(diabatic_CS), pointer :: diabatic_CSp !< KPP control structure needed to get BLD - type(lbd_CS), pointer :: CS !< Lateral boundary mixing control structure + type(hbd_CS), pointer :: CS !< Horizontal boundary mixing control structure ! local variables character(len=80) :: string ! Temporary strings - logical :: boundary_extrap ! controls if boundary extrapolation is used in the LBD code + logical :: boundary_extrap ! controls if boundary extrapolation is used in the HBD code if (ASSOCIATED(CS)) then - call MOM_error(FATAL, "lateral_boundary_diffusion_init called with associated control structure.") + call MOM_error(FATAL, "hor_bnd_diffusion_init called with associated control structure.") return endif ! Log this module and master switch for turning it on/off - call get_param(param_file, mdl, "USE_LATERAL_BOUNDARY_DIFFUSION", lateral_boundary_diffusion_init, & + call get_param(param_file, mdl, "USE_HORIZONTAL_BOUNDARY_DIFFUSION", hor_bnd_diffusion_init, & default=.false., do_not_log=.true.) call log_version(param_file, mdl, version, & - "This module implements lateral diffusion of tracers near boundaries", & - all_default=.not.lateral_boundary_diffusion_init) - call get_param(param_file, mdl, "USE_LATERAL_BOUNDARY_DIFFUSION", lateral_boundary_diffusion_init, & - "If true, enables the lateral boundary tracer's diffusion module.", & + "This module implements horizontal diffusion of tracers near boundaries", & + all_default=.not.hor_bnd_diffusion_init) + call get_param(param_file, mdl, "USE_HORIZONTAL_BOUNDARY_DIFFUSION", hor_bnd_diffusion_init, & + "If true, enables the horizonal boundary tracer's diffusion module.", & default=.false.) - if (.not. lateral_boundary_diffusion_init) return + if (.not. hor_bnd_diffusion_init) return allocate(CS) CS%diag => diag @@ -102,46 +111,55 @@ logical function lateral_boundary_diffusion_init(Time, G, GV, param_file, diag, call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp) call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp) + ! max. number of vertical layers + CS%hbd_nk = 2 + (GV%ke*2) + ! allocate the hbd grids and k_max + allocate(CS%hbd_grd_u(SZIB_(G),SZJ_(G),CS%hbd_nk), source=0.0) + allocate(CS%hbd_grd_v(SZI_(G),SZJB_(G),CS%hbd_nk), source=0.0) + allocate(CS%hbd_u_kmax(SZIB_(G),SZJ_(G)), source=0) + allocate(CS%hbd_v_kmax(SZI_(G),SZJB_(G)), source=0) + CS%surface_boundary_scheme = -1 if ( .not. ASSOCIATED(CS%energetic_PBL_CSp) .and. .not. ASSOCIATED(CS%KPP_CSp) ) then - call MOM_error(FATAL,"Lateral boundary diffusion is true, but no valid boundary layer scheme was found") + call MOM_error(FATAL,"Horizontal boundary diffusion is true, but no valid boundary layer scheme was found") endif ! Read all relevant parameters and write them to the model log. - call get_param(param_file, mdl, "LBD_LINEAR_TRANSITION", CS%linear, & + call get_param(param_file, mdl, "HBD_LINEAR_TRANSITION", CS%linear, & "If True, apply a linear transition at the base/top of the boundary. \n"//& "The flux will be fully applied at k=k_min and zero at k=k_max.", default=.false.) call get_param(param_file, mdl, "APPLY_LIMITER", CS%limiter, & "If True, apply a flux limiter in the native grid.", default=.true.) call get_param(param_file, mdl, "APPLY_LIMITER_REMAP", CS%limiter_remap, & "If True, apply a flux limiter in the remapped grid.", default=.false.) - call get_param(param_file, mdl, "LBD_BOUNDARY_EXTRAP", boundary_extrap, & - "Use boundary extrapolation in LBD code", & + call get_param(param_file, mdl, "HBD_BOUNDARY_EXTRAP", boundary_extrap, & + "Use boundary extrapolation in HBD code", & default=.false.) - call get_param(param_file, mdl, "LBD_REMAPPING_SCHEME", string, & + call get_param(param_file, mdl, "HBD_REMAPPING_SCHEME", string, & "This sets the reconstruction scheme used "//& "for vertical remapping for all variables. "//& "It can be one of the following schemes: "//& trim(remappingSchemesDoc), default=remappingDefaultScheme) - !### Revisit this hard-coded answer_date. + + ! GMM, TODO: add HBD params to control optional arguments in initialize_remapping. call initialize_remapping( CS%remap_CS, string, boundary_extrapolation = boundary_extrap ,& - check_reconstruction=.false., check_remapping=.false., answer_date=20190101) + check_reconstruction=.false., check_remapping=.false.) call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) - call get_param(param_file, mdl, "LBD_DEBUG", CS%debug, & - "If true, write out verbose debugging data in the LBD module.", & + call get_param(param_file, mdl, "HBD_DEBUG", CS%debug, & + "If true, write out verbose debugging data in the HBD module.", & default=.false.) - id_clock_lbd = cpu_clock_id('(Ocean LBD)', grain=CLOCK_MODULE) + id_clock_hbd = cpu_clock_id('(Ocean HBD)', grain=CLOCK_MODULE) -end function lateral_boundary_diffusion_init +end function hor_bnd_diffusion_init -!> Driver routine for calculating lateral diffusive fluxes near the top and bottom boundaries. +!> Driver routine for calculating horizontal diffusive fluxes near the top and bottom boundaries. !! Diffusion is applied using only information from neighboring cells, as follows: -!! 1) remap tracer to a z* grid (LBD grid) -!! 2) calculate diffusive tracer fluxes (F) in the LBD grid using a layer by layer approach +!! 1) remap tracer to a z* grid (HBD grid) +!! 2) calculate diffusive tracer fluxes (F) in the HBD grid using a layer by layer approach !! 3) remap fluxes to the native grid !! 4) update tracer by adding the divergence of F -subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) +subroutine hor_bnd_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) type(ocean_grid_type), intent(inout) :: G !< Grid type type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -152,7 +170,7 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) real, intent(in) :: dt !< Tracer time step * I_numitts !! (I_numitts in tracer_hordiff) [T ~> s] type(tracer_registry_type), pointer :: Reg !< Tracer registry - type(lbd_CS), pointer :: CS !< Control structure for this module + type(hbd_CS), pointer :: CS !< Control structure for this module ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: hbl !< Boundary layer depth [H ~> m or kg m-2] @@ -168,28 +186,32 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) type(tracer_type), pointer :: tracer => NULL() !< Pointer to the current tracer real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tracer_old !< local copy of the initial tracer concentration, !! only used to compute tendencies. - real :: tracer_int_prev !< Globally integrated tracer before LBD is applied, in mks units [conc kg] - real :: tracer_int_end !< Integrated tracer after LBD is applied, in mks units [conc kg] + real :: tracer_int_prev !< Globally integrated tracer before HBD is applied, in mks units [conc kg] + real :: tracer_int_end !< Integrated tracer after HBD is applied, in mks units [conc kg] real :: Idt !< inverse of the time step [T-1 ~> s-1] character(len=256) :: mesg !< Message for error messages. integer :: i, j, k, m !< indices to loop over - call cpu_clock_begin(id_clock_lbd) + call cpu_clock_begin(id_clock_hbd) Idt = 1./dt if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H) if (ASSOCIATED(CS%energetic_PBL_CSp)) call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, & m_to_MLD_units=GV%m_to_H) call pass_var(hbl,G%Domain) + + ! build HBD grid + call hbd_grid(SURFACE, G, GV, hbl, h, CS) + do m = 1,Reg%ntr ! current tracer tracer => Reg%tr(m) if (CS%debug) then - call hchksum(tracer%t, "before LBD "//tracer%name,G%HI) + call hchksum(tracer%t, "before HBD "//tracer%name,G%HI) endif ! for diagnostics - if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0 .or. CS%debug) then + if (tracer%id_hbdxy_conc > 0 .or. tracer%id_hbdxy_cont > 0 .or. tracer%id_hbdxy_cont_2d > 0 .or. CS%debug) then tendency(:,:,:) = 0.0 tracer_old(:,:,:) = tracer%t(:,:,:) endif @@ -198,13 +220,14 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) uFlx(:,:,:) = 0. vFlx(:,:,:) = 0. - ! LBD layer by layer + ! HBD layer by layer do j=G%jsc,G%jec do i=G%isc-1,G%iec if (G%mask2dCu(I,j)>0.) then - call fluxes_layer_method(SURFACE, G%ke, hbl(I,j), hbl(I+1,j), & + call fluxes_layer_method(SURFACE, G%ke, hbl(I,j), hbl(I+1,j), & h(I,j,:), h(I+1,j,:), tracer%t(I,j,:), tracer%t(I+1,j,:), & - Coef_x(I,j), uFlx(I,j,:), G%areaT(I,j), G%areaT(I+1,j), CS) + Coef_x(I,j), uFlx(I,j,:), G%areaT(I,j), G%areaT(I+1,j), CS%hbd_u_kmax(I,j), & + CS%hbd_grd_u(I,j,:), CS) endif enddo enddo @@ -213,7 +236,8 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) if (G%mask2dCv(i,J)>0.) then call fluxes_layer_method(SURFACE, GV%ke, hbl(i,J), hbl(i,J+1), & h(i,J,:), h(i,J+1,:), tracer%t(i,J,:), tracer%t(i,J+1,:), & - Coef_y(i,J), vFlx(i,J,:), G%areaT(i,J), G%areaT(i,J+1), CS) + Coef_y(i,J), vFlx(i,J,:), G%areaT(i,J), G%areaT(i,J+1), CS%hbd_v_kmax(i,J), & + CS%hbd_grd_v(i,J,:), CS) endif enddo enddo @@ -224,7 +248,7 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) tracer%t(i,j,k) = tracer%t(i,j,k) + (( (uFlx(I-1,j,k)-uFlx(I,j,k)) ) + ( (vFlx(i,J-1,k)-vFlx(i,J,k) ) ))* & G%IareaT(i,j) / ( h(i,j,k) + GV%H_subroundoff ) - if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0 ) then + if (tracer%id_hbdxy_conc > 0 .or. tracer%id_hbdxy_cont > 0 .or. tracer%id_hbdxy_cont_2d > 0 ) then tendency(i,j,k) = ((uFlx(I-1,j,k)-uFlx(I,j,k)) + (vFlx(i,J-1,k)-vFlx(i,J,k))) * & G%IareaT(i,j) * Idt endif @@ -239,64 +263,131 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) endif if (CS%debug) then - call hchksum(tracer%t, "after LBD "//tracer%name,G%HI) - ! tracer (native grid) integrated tracer amounts before and after LBD + call hchksum(tracer%t, "after HBD "//tracer%name,G%HI) + ! tracer (native grid) integrated tracer amounts before and after HBD tracer_int_prev = global_mass_integral(h, G, GV, tracer_old) tracer_int_end = global_mass_integral(h, G, GV, tracer%t) - write(mesg,*) 'Total '//tracer%name//' before/after LBD:', tracer_int_prev, tracer_int_end + write(mesg,*) 'Total '//tracer%name//' before/after HBD:', tracer_int_prev, tracer_int_end call MOM_mesg(mesg) endif ! Post the tracer diagnostics - if (tracer%id_lbd_dfx>0) call post_data(tracer%id_lbd_dfx, uFlx(:,:,:)*Idt, CS%diag) - if (tracer%id_lbd_dfy>0) call post_data(tracer%id_lbd_dfy, vFlx(:,:,:)*Idt, CS%diag) - if (tracer%id_lbd_dfx_2d>0) then + if (tracer%id_hbd_dfx>0) call post_data(tracer%id_hbd_dfx, uFlx(:,:,:)*Idt, CS%diag) + if (tracer%id_hbd_dfy>0) call post_data(tracer%id_hbd_dfy, vFlx(:,:,:)*Idt, CS%diag) + if (tracer%id_hbd_dfx_2d>0) then uwork_2d(:,:) = 0. do k=1,GV%ke ; do j=G%jsc,G%jec ; do I=G%isc-1,G%iec uwork_2d(I,j) = uwork_2d(I,j) + (uFlx(I,j,k) * Idt) enddo ; enddo ; enddo - call post_data(tracer%id_lbd_dfx_2d, uwork_2d, CS%diag) + call post_data(tracer%id_hbd_dfx_2d, uwork_2d, CS%diag) endif - if (tracer%id_lbd_dfy_2d>0) then + if (tracer%id_hbd_dfy_2d>0) then vwork_2d(:,:) = 0. do k=1,GV%ke ; do J=G%jsc-1,G%jec ; do i=G%isc,G%iec vwork_2d(i,J) = vwork_2d(i,J) + (vFlx(i,J,k) * Idt) enddo ; enddo ; enddo - call post_data(tracer%id_lbd_dfy_2d, vwork_2d, CS%diag) + call post_data(tracer%id_hbd_dfy_2d, vwork_2d, CS%diag) endif ! post tendency of tracer content - if (tracer%id_lbdxy_cont > 0) then - call post_data(tracer%id_lbdxy_cont, tendency, CS%diag) + if (tracer%id_hbdxy_cont > 0) then + call post_data(tracer%id_hbdxy_cont, tendency, CS%diag) endif ! post depth summed tendency for tracer content - if (tracer%id_lbdxy_cont_2d > 0) then + if (tracer%id_hbdxy_cont_2d > 0) then tendency_2d(:,:) = 0. do j=G%jsc,G%jec ; do i=G%isc,G%iec do k=1,GV%ke tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k) enddo enddo ; enddo - call post_data(tracer%id_lbdxy_cont_2d, tendency_2d, CS%diag) + call post_data(tracer%id_hbdxy_cont_2d, tendency_2d, CS%diag) endif ! post tendency of tracer concentration; this step must be ! done after posting tracer content tendency, since we alter ! the tendency array and its units. - if (tracer%id_lbdxy_conc > 0) then + if (tracer%id_hbdxy_conc > 0) then do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + CS%H_subroundoff ) enddo ; enddo ; enddo - call post_data(tracer%id_lbdxy_conc, tendency, CS%diag) + call post_data(tracer%id_hbdxy_conc, tendency, CS%diag) endif enddo - call cpu_clock_end(id_clock_lbd) + call cpu_clock_end(id_clock_hbd) + +end subroutine hor_bnd_diffusion + +!> Build the HBD grid where tracers will be rammaped to. +subroutine hbd_grid(boundary, G, GV, hbl, h, CS) + integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] + type(ocean_grid_type), intent(inout) :: G !< Grid type + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + real, dimension(SZI_(G),SZJ_(G)) :: hbl !< Boundary layer depth [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thickness in the native grid [H ~> m or kg m-2] + type(hbd_CS), pointer :: CS !< Horizontal diffusion control structure + + ! Local variables + real, allocatable :: dz_top(:) !< temporary HBD grid given by merge_interfaces [H ~> m or kg m-2] + integer :: nk, i, j, k !< number of layers in the HBD grid, and integers used in do-loops + + ! reset arrays + CS%hbd_grd_u(:,:,:) = 0.0 + CS%hbd_grd_v(:,:,:) = 0.0 + CS%hbd_u_kmax(:,:) = 0 + CS%hbd_v_kmax(:,:) = 0 + + do j=G%jsc,G%jec + do I=G%isc-1,G%iec + if (G%mask2dCu(I,j)>0.) then + call merge_interfaces(GV%ke, h(I,j,:), h(I+1,j,:), hbl(I,j), hbl(I+1,j), & + CS%H_subroundoff, dz_top) + nk = SIZE(dz_top) + if (nk > CS%hbd_nk) then + write(*,*)'nk, CS%hbd_nk', nk, CS%hbd_nk + call MOM_error(FATAL,"Houston, we've had a problem in hbd_grid, u-points (nk cannot be > CS%hbd_nk)") + endif + + CS%hbd_u_kmax(I,j) = nk + + ! set the HBD grid to dz_top + do k=1,nk + CS%hbd_grd_u(I,j,k) = dz_top(k) + enddo + deallocate(dz_top) + endif + enddo + enddo + + do J=G%jsc-1,G%jec + do i=G%isc,G%iec + if (G%mask2dCv(i,J)>0.) then + call merge_interfaces(GV%ke, h(i,J,:), h(i,J+1,:), hbl(i,J), hbl(i,J+1), & + CS%H_subroundoff, dz_top) + + nk = SIZE(dz_top) + if (nk > CS%hbd_nk) then + write(*,*)'nk, CS%hbd_nk', nk, CS%hbd_nk + call MOM_error(FATAL,"Houston, we've had a problem in hbd_grid, v-points (nk cannot be > CS%hbd_nk)") + endif + + CS%hbd_v_kmax(i,J) = nk + + ! set the HBD grid to dz_top + do k=1,nk + CS%hbd_grd_v(i,J,k) = dz_top(k) + enddo + deallocate(dz_top) + endif + enddo + enddo -end subroutine lateral_boundary_diffusion +end subroutine hbd_grid !> Calculate the harmonic mean of two quantities !! See \ref section_harmonic_mean. @@ -511,6 +602,7 @@ subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_b ! Local variables real :: htot ! Summed thickness [H ~> m or kg m-2] integer :: k + ! Surface boundary layer if ( boundary == SURFACE ) then k_top = 1 @@ -532,6 +624,7 @@ subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_b return endif enddo + ! Bottom boundary layer elseif ( boundary == BOTTOM ) then k_top = nk @@ -559,10 +652,10 @@ subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_b end subroutine boundary_k_range -!> Calculate the lateral boundary diffusive fluxes using the layer by layer method. +!> Calculate the horizontal boundary diffusive fluxes using the layer by layer method. !! See \ref section_method subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, & - khtr_u, F_layer, area_L, area_R, CS) + khtr_u, F_layer, area_L, area_R, nk, dz_top, CS) integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] integer, intent(in ) :: ke !< Number of layers in the native grid [nondim] @@ -580,10 +673,11 @@ subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_ !! in the native grid [H L2 conc ~> m3 conc] real, intent(in ) :: area_L !< Area of the horizontal grid (left) [L2 ~> m2] real, intent(in ) :: area_R !< Area of the horizontal grid (right) [L2 ~> m2] - type(lbd_CS), pointer :: CS !< Lateral diffusion control structure + integer, intent(in ) :: nk !< Number of layers in the HBD grid [nondim] + real, dimension(nk), intent(in ) :: dz_top !< The HBD z grid [H ~> m or kg m-2] + type(hbd_CS), pointer :: CS !< Horizontal diffusion control structure ! Local variables - real, allocatable :: dz_top(:) !< The LBD z grid to be created [H ~> m or kg m-2] real, allocatable :: phi_L_z(:) !< Tracer values in the ztop grid (left) [conc] real, allocatable :: phi_R_z(:) !< Tracer values in the ztop grid (right) [conc] real, allocatable :: F_layer_z(:) !< Diffusive flux at U/V-point in the ztop grid [H L2 conc ~> m3 conc] @@ -594,9 +688,6 @@ subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_ integer :: k_bot_min !< Minimum k-index for the bottom integer :: k_bot_max !< Maximum k-index for the bottom integer :: k_bot_diff !< Difference between bottom left and right k-indices - !integer :: k_top_max !< Minimum k-index for the top - !integer :: k_top_min !< Maximum k-index for the top - !integer :: k_top_diff !< Difference between top left and right k-indices integer :: k_top_L, k_bot_L !< k-indices left native grid integer :: k_top_R, k_bot_R !< k-indices right native grid real :: zeta_top_L, zeta_top_R !< distance from the top of a layer to the boundary @@ -607,17 +698,12 @@ subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_ real :: a !< coefficient to be used in the linear transition to the interior [nondim] real :: tmp1, tmp2 !< dummy variables [H ~> m or kg m-2] real :: htot_max !< depth below which no fluxes should be applied [H ~> m or kg m-2] - integer :: nk !< number of layers in the LBD grid F_layer(:) = 0.0 if (hbl_L == 0. .or. hbl_R == 0.) then return endif - ! Define vertical grid, dz_top - call merge_interfaces(ke, h_L(:), h_R(:), hbl_L, hbl_R, CS%H_subroundoff, dz_top) - nk = SIZE(dz_top) - ! allocate arrays allocate(phi_L_z(nk), source=0.0) allocate(phi_R_z(nk), source=0.0) @@ -669,27 +755,7 @@ subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_ endif endif -! TODO, boundary == BOTTOM -! if (boundary == BOTTOM) then -! ! TODO: GMM add option to apply linear decay -! k_top_max = MAX(k_top_L, k_top_R) -! ! make sure left and right k indices span same range -! if (k_top_max /= k_top_L) then -! k_top_L = k_top_max -! zeta_top_L = 1.0 -! endif -! if (k_top_max /= k_top_R) then -! k_top_R= k_top_max -! zeta_top_R = 1.0 -! endif -! -! ! tracer flux where the minimum BLD intersets layer -! F_layer(k_top_max) = (-heff * khtr_u) * (phi_R_avg - phi_L_avg) -! -! do k = k_top_max+1,nk -! F_layer_z(k) = -(heff * khtr_u) * (phi_R_z(k) - phi_L_z(k)) -! enddo -! endif + !GMM, TODO: boundary == BOTTOM ! thicknesses at velocity points do k = 1,ke @@ -723,7 +789,6 @@ subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_ enddo ! deallocated arrays - deallocate(dz_top) deallocate(phi_L_z) deallocate(phi_R_z) deallocate(F_layer_z) @@ -748,7 +813,7 @@ logical function near_boundary_unit_tests( verbose ) real :: zeta_top ! Nondimension position [nondim] integer :: k_bot ! Index of cell containing bottom of boundary real :: zeta_bot ! Nondimension position [nondim] - type(lbd_CS), pointer :: CS + type(hbd_CS), pointer :: CS allocate(CS) ! fill required fields in CS @@ -760,9 +825,11 @@ logical function near_boundary_unit_tests( verbose ) CS%debug=.false. CS%limiter=.false. CS%limiter_remap=.false. - + CS%hbd_nk = 2 + (2*2) + allocate(CS%hbd_grd_u(1,1,CS%hbd_nk), source=0.0) + allocate(CS%hbd_u_kmax(1,1), source=0) near_boundary_unit_tests = .false. - write(stdout,*) '==== MOM_lateral_boundary_diffusion =======================' + write(stdout,*) '==== MOM_hor_bnd_diffusion =======================' ! Unit tests for boundary_k_range test_name = 'Surface boundary spans the entire top cell' @@ -917,8 +984,9 @@ logical function near_boundary_unit_tests( verbose ) h_L = (/2.,2./) ; h_R = (/2.,2./) phi_L = (/0.,0./) ; phi_R = (/1.,1./) khtr_u = 1. + call hbd_grid_test(SURFACE, hbl_L, hbl_R, h_L, h_R, CS) call fluxes_layer_method(SURFACE, nk, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, & - khtr_u, F_layer, 1., 1., CS) + khtr_u, F_layer, 1., 1., CS%hbd_u_kmax(1,1), CS%hbd_grd_u(1,1,:), CS) near_boundary_unit_tests = near_boundary_unit_tests .or. & test_layer_fluxes( verbose, nk, test_name, F_layer, (/-2.0,0.0/) ) @@ -927,8 +995,9 @@ logical function near_boundary_unit_tests( verbose ) h_L = (/2.,2./) ; h_R = (/2.,2./) phi_L = (/2.,1./) ; phi_R = (/1.,1./) khtr_u = 0.5 + call hbd_grid_test(SURFACE, hbl_L, hbl_R, h_L, h_R, CS) call fluxes_layer_method(SURFACE, nk, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, & - khtr_u, F_layer, 1., 1., CS) + khtr_u, F_layer, 1., 1., CS%hbd_u_kmax(1,1), CS%hbd_grd_u(1,1,:), CS) near_boundary_unit_tests = near_boundary_unit_tests .or. & test_layer_fluxes( verbose, nk, test_name, F_layer, (/1.0,0.0/) ) @@ -937,8 +1006,9 @@ logical function near_boundary_unit_tests( verbose ) h_L = (/1.,2./) ; h_R = (/1.,2./) phi_L = (/0.,0./) ; phi_R = (/0.5,2./) khtr_u = 2. + call hbd_grid_test(SURFACE, hbl_L, hbl_R, h_L, h_R, CS) call fluxes_layer_method(SURFACE, nk, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, & - khtr_u, F_layer, 1., 1., CS) + khtr_u, F_layer, 1., 1., CS%hbd_u_kmax(1,1), CS%hbd_grd_u(1,1,:), CS) near_boundary_unit_tests = near_boundary_unit_tests .or. & test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.0,-4.0/) ) @@ -947,8 +1017,9 @@ logical function near_boundary_unit_tests( verbose ) h_L = (/6.,6./) ; h_R = (/10.,10./) phi_L = (/1.,1./) ; phi_R = (/1.,1./) khtr_u = 1. + call hbd_grid_test(SURFACE, hbl_L, hbl_R, h_L, h_R, CS) call fluxes_layer_method(SURFACE, nk, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, & - khtr_u, F_layer, 1., 1., CS) + khtr_u, F_layer, 1., 1., CS%hbd_u_kmax(1,1), CS%hbd_grd_u(1,1,:), CS) near_boundary_unit_tests = near_boundary_unit_tests .or. & test_layer_fluxes( verbose, nk, test_name, F_layer, (/0.,0./) ) @@ -958,9 +1029,9 @@ logical function near_boundary_unit_tests( verbose ) h_L = (/10.,5./) ; h_R = (/10.,0./) phi_L = (/1.,1./) ; phi_R = (/0.,0./) khtr_u = 1. + call hbd_grid_test(SURFACE, hbl_L, hbl_R, h_L, h_R, CS) call fluxes_layer_method(SURFACE, nk, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, & - khtr_u, F_layer, 1., 1., CS) - + khtr_u, F_layer, 1., 1., CS%hbd_u_kmax(1,1), CS%hbd_grd_u(1,1,:), CS) near_boundary_unit_tests = near_boundary_unit_tests .or. & test_layer_fluxes( verbose, nk, test_name, F_layer, (/10.,0.0/) ) @@ -983,7 +1054,7 @@ logical function test_layer_fluxes(verbose, nk, test_name, F_calc, F_ans) do k=1,nk if ( F_calc(k) /= F_ans(k) ) then test_layer_fluxes = .true. - write(stdout,*) "MOM_lateral_boundary_diffusion, UNIT TEST FAILED: ", test_name + write(stdout,*) "MOM_hor_bnd_diffusion, UNIT TEST FAILED: ", test_name write(stdout,10) k, F_calc(k), F_ans(k) elseif (verbose) then write(stdout,10) k, F_calc(k), F_ans(k) @@ -1026,13 +1097,55 @@ logical function test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, k_top_a end function test_boundary_k_range -!> \namespace mom_lateral_boundary_diffusion +!> Same as hbd_grid, but only used in the unit tests. +subroutine hbd_grid_test(boundary, hbl_L, hbl_R, h_L, h_R, CS) + integer, intent(in) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] + real, intent(in) :: hbl_L !< Boundary layer depth, left [H ~> m or kg m-2] + real, intent(in) :: hbl_R !< Boundary layer depth, right [H ~> m or kg m-2] + real, dimension(2), intent(in) :: h_L !< Layer thickness in the native grid, left [H ~> m or kg m-2] + real, dimension(2), intent(in) :: h_R !< Layer thickness in the native grid, right [H ~> m or kg m-2] + type(hbd_CS), pointer :: CS !< Horizontal diffusion control structure + + ! Local variables + real, allocatable :: dz_top(:) !< temporary HBD grid given by merge_interfaces [H ~> m or kg m-2] + integer :: nk, k !< number of layers in the HBD grid, and integers used in do-loops + + ! reset arrays + CS%hbd_grd_u(1,1,:) = 0.0 + CS%hbd_u_kmax(1,1) = 0 + + call merge_interfaces(2, h_L, h_R, hbl_L, hbl_R, CS%H_subroundoff, dz_top) + nk = SIZE(dz_top) + if (nk > CS%hbd_nk) then + write(*,*)'nk, CS%hbd_nk', nk, CS%hbd_nk + call MOM_error(FATAL,"Houston, we've had a problem in hbd_grid_test, (nk cannot be > CS%hbd_nk)") + endif + + CS%hbd_u_kmax(1,1) = nk + + ! set the HBD grid to dz_top + do k=1,nk + CS%hbd_grd_u(1,1,k) = dz_top(k) + enddo + deallocate(dz_top) + +end subroutine hbd_grid_test + +!> Deallocates hor_bnd_diffusion control structure +subroutine hor_bnd_diffusion_end(CS) + type(hbd_CS), pointer :: CS !< Horizontal boundary diffusion control structure + + if (associated(CS)) deallocate(CS) + +end subroutine hor_bnd_diffusion_end + +!> \namespace mom_hor_bnd_diffusion !! -!! \section section_LBD The Lateral Boundary Diffusion (LBD) framework +!! \section section_HBD The Horizontal Boundary Diffusion (HBD) framework !! -!! The LBD framework accounts for the effects of diabatic mesoscale fluxes +!! The HBD framework accounts for the effects of diabatic mesoscale fluxes !! within surface and bottom boundary layers. Unlike the equivalent adiabatic -!! fluxes, which is applied along neutral density surfaces, LBD is purely +!! fluxes, which is applied along neutral density surfaces, HBD is purely !! horizontal. To assure that diffusive fluxes are strictly horizontal !! regardless of the vertical coordinate system, this method relies on !! regridding/remapping techniques. @@ -1040,10 +1153,10 @@ end function test_boundary_k_range !! The bottom boundary layer fluxes remain to be implemented, although some !! of the steps needed to do so have already been added and tested. !! -!! Boundary lateral diffusion is applied as follows: +!! Horizontal boundary diffusion is applied as follows: !! -!! 1) remap tracer to a z* grid (LBD grid) -!! 2) calculate diffusive tracer fluxes (F) in the LBD grid using a layer by layer approach (@ref section_method) +!! 1) remap tracer to a z* grid (HBD grid) +!! 2) calculate diffusive tracer fluxes (F) in the HBD grid using a layer by layer approach (@ref section_method) !! 3) remap fluxes to the native grid !! 4) update tracer by adding the divergence of F !! @@ -1067,7 +1180,7 @@ end function test_boundary_k_range !! !! Step #3: option to linearly decay the flux from k_bot_min to k_bot_max: !! -!! If LBD_LINEAR_TRANSITION = True and k_bot_diff > 1, the diffusive flux will decay +!! If HBD_LINEAR_TRANSITION = True and k_bot_diff > 1, the diffusive flux will decay !! linearly between the top interface of the layer containing the minimum boundary !! layer depth (k_bot_min) and the lower interface of the layer containing the !! maximum layer depth (k_bot_max). @@ -1091,4 +1204,4 @@ end function test_boundary_k_range !! !! \f[ HM = \frac{2 \times h1 \times h2}{h1 + h2} \f] !! -end module MOM_lateral_boundary_diffusion +end module MOM_hor_bnd_diffusion diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index a34c2a2e58..cdabfa1277 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -27,8 +27,8 @@ module MOM_neutral_diffusion use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member -use MOM_lateral_boundary_diffusion, only : boundary_k_range, SURFACE, BOTTOM use MOM_io, only : stdout, stderr +use MOM_hor_bnd_diffusion, only : boundary_k_range, SURFACE, BOTTOM implicit none ; private diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index fc50fc4d1b..5d227defec 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -45,7 +45,7 @@ module MOM_tracer_flow_control use MOM_OCMIP2_CFC, only : OCMIP2_CFC_column_physics, OCMIP2_CFC_surface_state use MOM_OCMIP2_CFC, only : OCMIP2_CFC_stock, OCMIP2_CFC_end, OCMIP2_CFC_CS use MOM_CFC_cap, only : register_CFC_cap, initialize_CFC_cap -use MOM_CFC_cap, only : CFC_cap_column_physics, CFC_cap_surface_state +use MOM_CFC_cap, only : CFC_cap_column_physics, CFC_cap_set_forcing use MOM_CFC_cap, only : CFC_cap_stock, CFC_cap_end, CFC_cap_CS use oil_tracer, only : register_oil_tracer, initialize_oil_tracer use oil_tracer, only : oil_tracer_column_physics, oil_tracer_surface_state @@ -401,7 +401,7 @@ end subroutine get_chl_from_model !> This subroutine calls the individual tracer modules' subroutines to !! specify or read quantities related to their surface forcing. -subroutine call_tracer_set_forcing(sfc_state, fluxes, day_start, day_interval, G, CS) +subroutine call_tracer_set_forcing(sfc_state, fluxes, day_start, day_interval, G, US, Rho0, CS) type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the @@ -413,6 +413,8 @@ subroutine call_tracer_set_forcing(sfc_state, fluxes, day_start, day_interval, G type(time_type), intent(in) :: day_interval !< Length of time over which these !! fluxes will be applied. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, intent(in) :: Rho0 !< The mean ocean density [R ~> kg m-3] type(tracer_flow_control_CS), pointer :: CS !< The control structure returned by a !! previous call to call_tracer_register. @@ -421,6 +423,9 @@ subroutine call_tracer_set_forcing(sfc_state, fluxes, day_start, day_interval, G ! if (CS%use_ideal_age) & ! call ideal_age_tracer_set_forcing(sfc_state, fluxes, day_start, day_interval, & ! G, CS%ideal_age_tracer_CSp) + if (CS%use_CFC_cap) & + call CFC_cap_set_forcing(sfc_state, fluxes, day_start, day_interval, G, US, Rho0, & + CS%CFC_cap_CSp) end subroutine call_tracer_set_forcing @@ -850,8 +855,6 @@ subroutine call_tracer_surface_state(sfc_state, h, G, GV, US, CS) call advection_test_tracer_surface_state(sfc_state, h, G, GV, CS%advection_test_tracer_CSp) if (CS%use_OCMIP2_CFC) & call OCMIP2_CFC_surface_state(sfc_state, h, G, GV, US, CS%OCMIP2_CFC_CSp) - if (CS%use_CFC_cap) & - call CFC_cap_surface_state(sfc_state, G, CS%CFC_cap_CSp) if (CS%use_MOM_generic_tracer) & call MOM_generic_tracer_surface_state(sfc_state, h, G, GV, CS%MOM_generic_tracer_CSp) diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 14be088e8d..79e99f8bb7 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -3,32 +3,32 @@ module MOM_tracer_hor_diff ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end -use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE -use MOM_diag_mediator, only : post_data, diag_ctrl -use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type -use MOM_domains, only : sum_across_PEs, max_across_PEs -use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type -use MOM_domains, only : pass_vector -use MOM_debugging, only : hchksum, uvchksum -use MOM_diabatic_driver, only : diabatic_CS -use MOM_EOS, only : calculate_density, EOS_type, EOS_domain -use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe -use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery -use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint -use MOM_file_parser, only : get_param, log_version, param_file_type -use MOM_grid, only : ocean_grid_type -use MOM_lateral_mixing_coeffs, only : VarMix_CS -use MOM_MEKE_types, only : MEKE_type -use MOM_neutral_diffusion, only : neutral_diffusion_init, neutral_diffusion_end -use MOM_neutral_diffusion, only : neutral_diffusion_CS -use MOM_neutral_diffusion, only : neutral_diffusion_calc_coeffs, neutral_diffusion -use MOM_lateral_boundary_diffusion, only : lbd_CS, lateral_boundary_diffusion_init -use MOM_lateral_boundary_diffusion, only : lateral_boundary_diffusion -use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chksum -use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs -use MOM_verticalGrid, only : verticalGrid_type +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end +use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE +use MOM_diag_mediator, only : post_data, diag_ctrl +use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type +use MOM_domains, only : sum_across_PEs, max_across_PEs +use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type +use MOM_domains, only : pass_vector +use MOM_debugging, only : hchksum, uvchksum +use MOM_diabatic_driver, only : diabatic_CS +use MOM_EOS, only : calculate_density, EOS_type, EOS_domain +use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe +use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery +use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint +use MOM_file_parser, only : get_param, log_version, param_file_type +use MOM_grid, only : ocean_grid_type +use MOM_lateral_mixing_coeffs, only : VarMix_CS +use MOM_MEKE_types, only : MEKE_type +use MOM_neutral_diffusion, only : neutral_diffusion_init, neutral_diffusion_end +use MOM_neutral_diffusion, only : neutral_diffusion_CS +use MOM_neutral_diffusion, only : neutral_diffusion_calc_coeffs, neutral_diffusion +use MOM_hor_bnd_diffusion, only : hbd_CS, hor_bnd_diffusion_init +use MOM_hor_bnd_diffusion, only : hor_bnd_diffusion, hor_bnd_diffusion_end +use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chksum +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -59,13 +59,13 @@ module MOM_tracer_hor_diff !! the CFL limit is not violated. logical :: use_neutral_diffusion !< If true, use the neutral_diffusion module from within !! tracer_hor_diff. - logical :: use_lateral_boundary_diffusion !< If true, use the lateral_boundary_diffusion module from within + logical :: use_hor_bnd_diffusion !< If true, use the hor_bnd_diffusion module from within !! tracer_hor_diff. logical :: recalc_neutral_surf !< If true, recalculate the neutral surfaces if CFL has been !! exceeded type(neutral_diffusion_CS), pointer :: neutral_diffusion_CSp => NULL() !< Control structure for neutral diffusion. - type(lbd_CS), pointer :: lateral_boundary_diffusion_CSp => NULL() !< Control structure for - !! lateral boundary mixing. + type(hbd_CS), pointer :: hor_bnd_diffusion_CSp => NULL() !< Control structure for + !! horizontal boundary diffusion. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. logical :: debug !< If true, write verbose checksums for debugging purposes. @@ -387,9 +387,9 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online endif enddo - if (CS%use_lateral_boundary_diffusion) then + if (CS%use_hor_bnd_diffusion) then - if (CS%show_call_tree) call callTree_waypoint("Calling lateral boundary mixing (tracer_hordiff)") + if (CS%show_call_tree) call callTree_waypoint("Calling horizontal boundary diffusion (tracer_hordiff)") call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) @@ -403,12 +403,12 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online enddo do itt=1,num_itts - if (CS%show_call_tree) call callTree_waypoint("Calling lateral boundary diffusion (tracer_hordiff)",itt) + if (CS%show_call_tree) call callTree_waypoint("Calling horizontal boundary diffusion (tracer_hordiff)",itt) if (itt>1) then ! Update halos for subsequent iterations call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) endif - call lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, I_numitts*dt, Reg, & - CS%lateral_boundary_diffusion_CSp) + call hor_bnd_diffusion(G, GV, US, h, Coef_x, Coef_y, I_numitts*dt, Reg, & + CS%hor_bnd_diffusion_CSp) enddo ! itt endif @@ -418,7 +418,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) ! We are assuming that neutral surfaces do not evolve (much) as a result of multiple - ! lateral diffusion iterations. Otherwise the call to neutral_diffusion_calc_coeffs() + !horizontal diffusion iterations. Otherwise the call to neutral_diffusion_calc_coeffs() ! would be inside the itt-loop. -AJA if (associated(tv%p_surf)) then @@ -1541,10 +1541,10 @@ subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic diabatic_CSp, CS%neutral_diffusion_CSp ) if (CS%use_neutral_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & "USE_NEUTRAL_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") - CS%use_lateral_boundary_diffusion = lateral_boundary_diffusion_init(Time, G, GV, param_file, diag, diabatic_CSp, & - CS%lateral_boundary_diffusion_CSp) - if (CS%use_lateral_boundary_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & - "USE_LATERAL_BOUNDARY_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") + CS%use_hor_bnd_diffusion = hor_bnd_diffusion_init(Time, G, GV, US, param_file, diag, diabatic_CSp, & + CS%hor_bnd_diffusion_CSp) + if (CS%use_hor_bnd_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & + "USE_HORIZONTAL_BOUNDARY_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.) @@ -1584,6 +1584,7 @@ subroutine tracer_hor_diff_end(CS) type(tracer_hor_diff_CS), pointer :: CS !< module control structure call neutral_diffusion_end(CS%neutral_diffusion_CSp) + call hor_bnd_diffusion_end(CS%hor_bnd_diffusion_CSp) if (associated(CS)) deallocate(CS) end subroutine tracer_hor_diff_end diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index c3f5f64edf..6e1f2bee5a 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -358,14 +358,14 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u diag%axesCvL, Time, trim(flux_longname)//" diffusive meridional flux" , & trim(flux_units), v_extensive=.true., x_cell_method='sum', & conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T) - Tr%id_lbd_dfx = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffx", & - diag%axesCuL, Time, trim(flux_longname)//" diffusive zonal flux from the lateral boundary diffusion "//& - "scheme", trim(flux_units), v_extensive=.true., y_cell_method='sum', & - conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T) - Tr%id_lbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffy", & - diag%axesCvL, Time, trim(flux_longname)//" diffusive meridional flux from the lateral boundary diffusion "//& - "scheme", trim(flux_units), v_extensive=.true., x_cell_method='sum', & - conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T) + Tr%id_hbd_dfx = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffx", & + diag%axesCuL, Time, trim(flux_longname)//" diffusive zonal flux " //& + "from the horizontal boundary diffusion scheme", trim(flux_units), v_extensive=.true., & + y_cell_method='sum', conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T) + Tr%id_hbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffy", & + diag%axesCvL, Time, trim(flux_longname)//" diffusive meridional " //& + "flux from the horizontal boundary diffusion scheme", trim(flux_units), v_extensive=.true., & + x_cell_method='sum', conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T) else Tr%id_adx = register_diag_field("ocean_model", trim(shortnm)//"_adx", & diag%axesCuL, Time, "Advective (by residual mean) Zonal Flux of "//trim(flux_longname), & @@ -381,12 +381,12 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u diag%axesCvL, Time, "Diffusive Meridional Flux of "//trim(flux_longname), & flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & x_cell_method='sum') - Tr%id_lbd_dfx = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffx", & - diag%axesCuL, Time, "Lateral Boundary Diffusive Zonal Flux of "//trim(flux_longname), & + Tr%id_hbd_dfx = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffx", & + diag%axesCuL, Time, "Horizontal Boundary Diffusive Zonal Flux of "//trim(flux_longname), & flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & y_cell_method='sum') - Tr%id_lbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffy", & - diag%axesCvL, Time, "Lateral Boundary Diffusive Meridional Flux of "//trim(flux_longname), & + Tr%id_hbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffy", & + diag%axesCvL, Time, "Horizontal Boundary Diffusive Meridional Flux of "//trim(flux_longname), & flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & x_cell_method='sum') endif @@ -394,8 +394,8 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u if (Tr%id_ady > 0) call safe_alloc_ptr(Tr%ad_y,isd,ied,JsdB,JedB,nz) if (Tr%id_dfx > 0) call safe_alloc_ptr(Tr%df_x,IsdB,IedB,jsd,jed,nz) if (Tr%id_dfy > 0) call safe_alloc_ptr(Tr%df_y,isd,ied,JsdB,JedB,nz) - if (Tr%id_lbd_dfx > 0) call safe_alloc_ptr(Tr%lbd_dfx,IsdB,IedB,jsd,jed,nz) - if (Tr%id_lbd_dfy > 0) call safe_alloc_ptr(Tr%lbd_dfy,isd,ied,JsdB,JedB,nz) + if (Tr%id_hbd_dfx > 0) call safe_alloc_ptr(Tr%hbd_dfx,IsdB,IedB,jsd,jed,nz) + if (Tr%id_hbd_dfy > 0) call safe_alloc_ptr(Tr%hbd_dfy,isd,ied,JsdB,JedB,nz) Tr%id_adx_2d = register_diag_field("ocean_model", trim(shortnm)//"_adx_2d", & diag%axesCu1, Time, & @@ -415,22 +415,12 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u "Vertically Integrated Diffusive Meridional Flux of "//trim(flux_longname), & flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & x_cell_method='sum') - Tr%id_lbd_bulk_dfx = register_diag_field("ocean_model", trim(shortnm)//"_lbd_bulk_diffx", & - diag%axesCu1, Time, & - "Total Bulk Diffusive Zonal Flux of "//trim(flux_longname), & - flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & - y_cell_method='sum') - Tr%id_lbd_bulk_dfy = register_diag_field("ocean_model", trim(shortnm)//"_lbd_bulk_diffy", & - diag%axesCv1, Time, & - "Total Bulk Diffusive Meridional Flux of "//trim(flux_longname), & - flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & - x_cell_method='sum') - Tr%id_lbd_dfx_2d = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffx_2d", & - diag%axesCu1, Time, "Vertically-integrated zonal diffusive flux from the lateral boundary diffusion "//& + Tr%id_hbd_dfx_2d = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffx_2d", & + diag%axesCu1, Time, "Vertically-integrated zonal diffusive flux from the horizontal boundary diffusion "//& "scheme for "//trim(flux_longname), flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & y_cell_method='sum') - Tr%id_lbd_dfy_2d = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffy_2d", & - diag%axesCv1, Time, "Vertically-integrated meridional diffusive flux from the lateral boundary diffusion "//& + Tr%id_hbd_dfy_2d = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffy_2d", & + diag%axesCv1, Time, "Vertically-integrated meridional diffusive flux from the horizontal boundary diffusion "//& "scheme for "//trim(flux_longname), flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & x_cell_method='sum') @@ -438,10 +428,8 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u if (Tr%id_ady_2d > 0) call safe_alloc_ptr(Tr%ad2d_y,isd,ied,JsdB,JedB) if (Tr%id_dfx_2d > 0) call safe_alloc_ptr(Tr%df2d_x,IsdB,IedB,jsd,jed) if (Tr%id_dfy_2d > 0) call safe_alloc_ptr(Tr%df2d_y,isd,ied,JsdB,JedB) - if (Tr%id_lbd_bulk_dfx > 0) call safe_alloc_ptr(Tr%lbd_bulk_df_x,IsdB,IedB,jsd,jed) - if (Tr%id_lbd_bulk_dfy > 0) call safe_alloc_ptr(Tr%lbd_bulk_df_y,isd,ied,JsdB,JedB) - if (Tr%id_lbd_dfx_2d > 0) call safe_alloc_ptr(Tr%lbd_dfx_2d,IsdB,IedB,jsd,jed) - if (Tr%id_lbd_dfy_2d > 0) call safe_alloc_ptr(Tr%lbd_dfy_2d,isd,ied,JsdB,JedB) + if (Tr%id_hbd_dfx_2d > 0) call safe_alloc_ptr(Tr%hbd_dfx_2d,IsdB,IedB,jsd,jed) + if (Tr%id_hbd_dfy_2d > 0) call safe_alloc_ptr(Tr%hbd_dfy_2d,isd,ied,JsdB,JedB) Tr%id_adv_xy = register_diag_field('ocean_model', trim(shortnm)//"_advection_xy", & diag%axesTL, Time, & @@ -466,7 +454,7 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u enddo ; enddo ; enddo endif - ! Neutral/Lateral diffusion convergence tendencies + ! Neutral/Horizontal diffusion convergence tendencies if (Tr%diag_form == 1) then Tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', & diag%axesTL, Time, "Neutral diffusion tracer content tendency for "//trim(shortnm), & @@ -477,12 +465,12 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u "tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & x_cell_method='sum', y_cell_method='sum') - Tr%id_lbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency', & - diag%axesTL, Time, "Lateral diffusion tracer content tendency for "//trim(shortnm), & + Tr%id_hbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_cont_tendency', & + diag%axesTL, Time, "Horizontal boundary diffusion tracer content tendency for "//trim(shortnm), & conv_units, conversion=Tr%conv_scale*US%s_to_T, x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) - Tr%id_lbdxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency_2d', & - diag%axesT1, Time, "Depth integrated lateral diffusion tracer content "//& + Tr%id_hbdxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_cont_tendency_2d', & + diag%axesT1, Time, "Depth integrated horizontal boundary diffusion tracer content "//& "tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & x_cell_method='sum', y_cell_method='sum') else @@ -503,13 +491,13 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u cmor_long_name=trim(cmor_var_lname), cmor_standard_name=trim(cmor_long_std(cmor_var_lname)), & x_cell_method='sum', y_cell_method='sum') - Tr%id_lbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency', & - diag%axesTL, Time, "Lateral diffusion tracer content tendency for "//trim(shortnm), & + Tr%id_hbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_cont_tendency', & + diag%axesTL, Time, "Horizontal boundary diffusion tracer content tendency for "//trim(shortnm), & conv_units, conversion=Tr%conv_scale*US%s_to_T, & x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) - Tr%id_lbdxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency_2d', & - diag%axesT1, Time, "Depth integrated lateral diffusion tracer "//& + Tr%id_hbdxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_cont_tendency_2d', & + diag%axesT1, Time, "Depth integrated horizontal boundary diffusion of tracer "//& "content tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & x_cell_method='sum', y_cell_method='sum') endif @@ -517,8 +505,8 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u diag%axesTL, Time, "Neutral diffusion tracer concentration tendency for "//trim(shortnm), & trim(units)//' s-1', conversion=Tr%conc_scale*US%s_to_T) - Tr%id_lbdxy_conc = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_conc_tendency', & - diag%axesTL, Time, "Lateral diffusion tracer concentration tendency for "//trim(shortnm), & + Tr%id_hbdxy_conc = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_conc_tendency', & + diag%axesTL, Time, "Horizontal diffusion tracer concentration tendency for "//trim(shortnm), & trim(units)//' s-1', conversion=Tr%conc_scale*US%s_to_T) var_lname = "Net time tendency for "//lowercase(flux_longname) diff --git a/src/tracer/MOM_tracer_types.F90 b/src/tracer/MOM_tracer_types.F90 index 51c4508db6..bdae8bcee9 100644 --- a/src/tracer/MOM_tracer_types.F90 +++ b/src/tracer/MOM_tracer_types.F90 @@ -27,20 +27,16 @@ module MOM_tracer_types real, dimension(:,:,:), pointer :: df_x => NULL() !< diagnostic array for x-diffusive tracer flux !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), pointer :: df_y => NULL() !< diagnostic array for y-diffusive tracer flux - !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] - real, dimension(:,:,:), pointer :: lbd_dfx => NULL() !< diagnostic array for x-diffusive tracer flux - !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] - real, dimension(:,:,:), pointer :: lbd_dfy => NULL() !< diagnostic array for y-diffusive tracer flux - !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] - real, dimension(:,:), pointer :: lbd_dfx_2d => NULL() !< diagnostic array for x-diffusive tracer flux - !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] - real, dimension(:,:), pointer :: lbd_dfy_2d => NULL() !< diagnostic array for y-diffusive tracer flux - !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] - !### These two arrays may be allocated but are never used. - real, dimension(:,:), pointer :: lbd_bulk_df_x => NULL() !< diagnostic array for x-diffusive tracer flux - !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] - real, dimension(:,:), pointer :: lbd_bulk_df_y => NULL() !< diagnostic array for y-diffusive tracer flux - !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:,:), pointer :: hbd_dfx => NULL() !< diagnostic array for x-diffusive tracer flux + !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:,:), pointer :: hbd_dfy => NULL() !< diagnostic array for y-diffusive tracer flux + !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:), pointer :: hbd_dfx_2d => NULL() !< diagnostic array for x-diffusive tracer flux + !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + real, dimension(:,:), pointer :: hbd_dfy_2d => NULL() !< diagnostic array for y-diffusive tracer flux + !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: df2d_x => NULL() !< diagnostic vertical sum x-diffusive flux !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: df2d_y => NULL() !< diagnostic vertical sum y-diffusive flux @@ -106,12 +102,12 @@ module MOM_tracer_types !>@{ Diagnostic IDs integer :: id_tr = -1, id_tr_post_horzn = -1 integer :: id_adx = -1, id_ady = -1, id_dfx = -1, id_dfy = -1 - integer :: id_lbd_bulk_dfx = -1, id_lbd_bulk_dfy = -1, id_lbd_dfx = -1, id_lbd_dfy = -1 - integer :: id_lbd_dfx_2d = -1 , id_lbd_dfy_2d = -1 + integer :: id_hbd_dfx = -1, id_hbd_dfy = -1 + integer :: id_hbd_dfx_2d = -1, id_hbd_dfy_2d = -1 integer :: id_adx_2d = -1, id_ady_2d = -1, id_dfx_2d = -1, id_dfy_2d = -1 integer :: id_adv_xy = -1, id_adv_xy_2d = -1 integer :: id_dfxy_cont = -1, id_dfxy_cont_2d = -1, id_dfxy_conc = -1 - integer :: id_lbdxy_cont = -1, id_lbdxy_cont_2d = -1, id_lbdxy_conc = -1 + integer :: id_hbdxy_cont = -1, id_hbdxy_cont_2d = -1, id_hbdxy_conc = -1 integer :: id_remap_conc = -1, id_remap_cont = -1, id_remap_cont_2d = -1 integer :: id_tendency = -1, id_trxh_tendency = -1, id_trxh_tendency_2d = -1 integer :: id_tr_vardec = -1 diff --git a/src/tracer/oil_tracer.F90 b/src/tracer/oil_tracer.F90 index 20685d9711..71800284a6 100644 --- a/src/tracer/oil_tracer.F90 +++ b/src/tracer/oil_tracer.F90 @@ -113,7 +113,7 @@ function register_oil_tracer(HI, GV, US, param_file, CS, tr_Reg, restart_CS) ! Add the directory if CS%IC_file is not already a complete path. call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") CS%IC_file = trim(slasher(inputdir))//trim(CS%IC_file) - call log_param(param_file, mdl, "INPUTDIR/CFC_IC_FILE", CS%IC_file) + call log_param(param_file, mdl, "INPUTDIR/OIL_IC_FILE", CS%IC_file) endif call get_param(param_file, mdl, "OIL_IC_FILE_IS_Z", CS%Z_IC_file, & "If true, OIL_IC_FILE is in depth space, not layer space", &