diff --git a/config_src/external/GFDL_ocean_BGC/generic_tracer.F90 b/config_src/external/GFDL_ocean_BGC/generic_tracer.F90 index 75cd57a08d..1e97d91eb4 100644 --- a/config_src/external/GFDL_ocean_BGC/generic_tracer.F90 +++ b/config_src/external/GFDL_ocean_BGC/generic_tracer.F90 @@ -70,7 +70,7 @@ end subroutine generic_tracer_coupler_accumulate !> Calls the corresponding generic_X_update_from_source routine for each package X subroutine generic_tracer_source(Temp,Salt,rho_dzt,dzt,hblt_depth,ilb,jlb,tau,dtts,& grid_dat,model_time,nbands,max_wavelength_band,sw_pen_band,opacity_band,internal_heat,& - frunoff,grid_ht, current_wave_stress, sosga, geolat, eqn_of_state) + frunoff,grid_ht, current_wave_stress, sosga, geolat, eqn_of_state, photo_acc_dpth) integer, intent(in) :: ilb !< Lower bounds of x extent of input arrays on data domain integer, intent(in) :: jlb !< Lower bounds of y extent of input arrays on data domain real, dimension(ilb:,jlb:,:), intent(in) :: Temp !< Potential temperature [deg C] @@ -98,6 +98,7 @@ subroutine generic_tracer_source(Temp,Salt,rho_dzt,dzt,hblt_depth,ilb,jlb,tau,dt real, optional , intent(in) :: sosga !< Global average sea surface salinity [ppt] real, dimension(ilb:,jlb:),optional, intent(in) :: geolat !< Latitude type(EOS_type), optional, intent(in) :: eqn_of_state !< A pointer to the equation of state + real, dimension(ilb:,jlb:), optional, intent(in) :: photo_acc_dpth end subroutine generic_tracer_source !> Update the tracers from bottom fluxes diff --git a/src/diagnostics/MOM_diagnose_MLD.F90 b/src/diagnostics/MOM_diagnose_MLD.F90 index 29b66ef6ac..b42c974a95 100644 --- a/src/diagnostics/MOM_diagnose_MLD.F90 +++ b/src/diagnostics/MOM_diagnose_MLD.F90 @@ -30,7 +30,7 @@ module MOM_diagnose_mld !> Diagnose a mixed layer depth (MLD) determined by a given density difference with the surface. !> This routine is appropriate in MOM_diabatic_aux due to its position within the time stepping. subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, diagPtr, & - ref_h_mld, id_ref_z, id_ref_rho, id_N2subML, id_MLDsq, dz_subML) + ref_h_mld, id_ref_z, id_ref_rho, id_N2subML, id_MLDsq, dz_subML, MLD_out) type(ocean_grid_type), intent(in) :: 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 @@ -44,6 +44,8 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, real, intent(in) :: ref_h_mld !< Depth of the calculated "surface" densisty [Z ~> m] integer, intent(in) :: id_ref_z !< Handle (ID) of reference depth diagnostic integer, intent(in) :: id_ref_rho !< Handle (ID) of reference density diagnostic + real, dimension(SZI_(G),SZJ_(G)), & + optional, intent(inout) :: MLD_out !< Send MLD to other routines [Z ~> m] integer, optional, intent(in) :: id_N2subML !< Optional handle (ID) of subML stratification integer, optional, intent(in) :: id_MLDsq !< Optional handle (ID) of squared MLD real, optional, intent(in) :: dz_subML !< The distance over which to calculate N2subML @@ -234,11 +236,13 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, if ((id_ref_z > 0) .and. (pRef_MLD(is)/=0.)) call post_data(id_ref_z, z_ref_diag , diagPtr) if (id_ref_rho > 0) call post_data(id_ref_rho, rhoSurf_2d , diagPtr) + if (present(MLD_out)) MLD_out(:,:) = MLD(:,:) + end subroutine diagnoseMLDbyDensityDifference !> Diagnose a mixed layer depth (MLD) determined by the depth a given energy value would mix. !> This routine is appropriate in MOM_diabatic_aux due to its position within the time stepping. -subroutine diagnoseMLDbyEnergy(id_MLD, h, tv, G, GV, US, Mixing_Energy, diagPtr) +subroutine diagnoseMLDbyEnergy(id_MLD, h, tv, G, GV, US, Mixing_Energy, diagPtr, MLD_out) ! Author: Brandon Reichl ! Date: October 2, 2020 ! // @@ -270,6 +274,8 @@ subroutine diagnoseMLDbyEnergy(id_MLD, h, tv, G, GV, US, Mixing_Energy, diagPtr) type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any !! available thermodynamic fields. type(diag_ctrl), pointer :: diagPtr !< Diagnostics structure + real, dimension(SZI_(G),SZJ_(G)), & + optional, intent(inout) :: MLD_out !< Send MLD to other routines [Z ~> m] ! Local variables real, dimension(SZI_(G),SZJ_(G),3) :: MLD ! Diagnosed mixed layer depth [Z ~> m]. @@ -467,6 +473,8 @@ subroutine diagnoseMLDbyEnergy(id_MLD, h, tv, G, GV, US, Mixing_Energy, diagPtr) if (id_MLD(2) > 0) call post_data(id_MLD(2), MLD(:,:,2), diagPtr) if (id_MLD(3) > 0) call post_data(id_MLD(3), MLD(:,:,3), diagPtr) + if (present(MLD_out)) MLD_out(:,:) = MLD(:,:,1) + end subroutine diagnoseMLDbyEnergy !> \namespace mom_diagnose_mld diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 943e1f9f48..d943886c8a 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -32,6 +32,7 @@ module MOM_generic_tracer use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS use MOM_ALE_sponge, only : ALE_sponge_CS, initialize_ALE_sponge use MOM_coms, only : EFP_type, max_across_PEs, min_across_PEs, PE_here + use MOM_diagnose_mld, only : diagnoseMLDbyDensityDifference, diagnoseMLDbyEnergy use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_ctrl, get_diag_time_end use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, is_root_pe @@ -84,6 +85,14 @@ module MOM_generic_tracer logical :: tracers_may_reinit !< If true, tracers may go through the !! initialization code if they are not found in the restart files. + logical :: mld_pha_calc = .False. !< If true, use a fixed value for photoacclimation MLD + real :: mld_pha_val = 0.0 !< The value of fixed photoacclimation MLD + logical :: mld_pha_use_delta_rho = .False. !< If true, use a density diference to find the MLD + real :: mld_pha_href = 0.0 !< The reference depth for density difference based MLD + real :: mld_pha_drho = 0.03 !< The density thershold for a density difference based MLD + logical :: mld_pha_use_delta_eng = .False. !< If true, use an energy diference to find the MLD + real :: mld_pha_deng = 25.0 !< The energy threshold for an energy d ifference based MLD + type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< Restart control structure @@ -425,6 +434,42 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, tv, param_f enddo !! end section to re-initialize generic tracers + call get_param(param_file, "MOM", "PHA_MLD_CALC", CS%mld_pha_calc, & + "If false, use a fixed value for the photoacclimation mixed layer depth within the "//& + "generic tracer update. This MLD is only used for photoacclimation. This variable should "//& + "be set to true if using COBALTv3 for the BGC.", default=.false.) + if (CS%mld_pha_calc) then + call get_param(param_file, "MOM", "PHA_MLD_USE_DELTA_RHO", CS%mld_pha_use_delta_rho, & + "If true, use a density difference to find the photoacclimation mixed layer depth "//& + "within the generic tracer update. This MLD is only used for photoacclimation.", default=.false.) + call get_param(param_file, "MOM", "PHA_MLD_USE_DELTA_ENG", CS%mld_pha_use_delta_eng, & + "If true, use an energy difference to find the photoacclimation mixed layer depth "//& + "with the generic tracer update. This MLD is only used for photoacclimation.", default=.false.) + if (CS%mld_pha_use_delta_rho .and. CS%mld_pha_use_delta_eng) then + call MOM_error(FATAL, "PHA_MLD_CALC is set to true and PHA_MLD_USE_DELTA_RHO and PHA_MLD_USE_DELTA_ENG "// & + "are both true. Choose only one option for the calculated photoacclimation MLD!") + elseif ((.not.CS%mld_pha_use_delta_rho) .and. (.not.CS%mld_pha_use_delta_eng)) then + call MOM_error(FATAL, "PHA_MLD_CALC is set to true but PHA_MLD_USE_DELTA_RHO and PHA_MLD_USE_DELTA_ENG "// & + "are both false. Choose an option for the calculated photoacclimation MLD!") + endif + if (CS%mld_pha_use_delta_rho) then + call get_param(param_file, "MOM", "PHA_MLD_HREF", CS%mld_pha_href, & + "The reference depth for a density difference based photoacclimation MLD [m].", & + units='m', default=0.0, scale=US%m_to_Z, do_not_log=.not.CS%mld_pha_use_delta_rho) + call get_param(param_file, "MOM", "PHA_MLD_DRHO", CS%mld_pha_drho, & + "The density difference for a density difference based photoacclimation MLD [kg m-3].", & + units='kg/m3', default=0.03, scale=US%kg_m3_to_R, do_not_log=.not.CS%mld_pha_use_delta_rho) + elseif (CS%mld_pha_use_delta_eng) then + call get_param(param_file, "MOM", "PHA_MLD_DENG", CS%mld_pha_deng, & + "The energy for an energy difference based photoacclimation MLD.", default=25.0, & + units='J/m2',scale=US%W_m2_to_RZ3_T3*US%s_to_T, do_not_log=.not.CS%mld_pha_use_delta_eng) + endif + else + call get_param(param_file, "MOM", "PHA_MLD_VAL", CS%mld_pha_val, & + "The depth of photoacclimation if fixed depth is used [m].", & + units='m', default=0.0, scale=US%m_to_Z) + endif + !Now we can reset the grid mask, axes and time to their true values !Note that grid_tmask must be set correctly on the data domain boundary @@ -654,6 +699,8 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: rho_dzt ! Layer mass per unit area [kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dzt ! Layer vertical extents [m] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! A work array of thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)) :: mld_pha ! The mixed layer depth calculated for photoacclimation + ! that is used in COBALTv3 integer :: i, j, k, isc, iec, jsc, jec, nk isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke @@ -719,6 +766,17 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, enddo ; enddo sosga = global_area_mean(surface_field, G, unscale=US%S_to_ppt) + mld_pha(:,:) = 0.0 + if (.not.CS%mld_pha_calc) then + mld_pha(:,:) = CS%mld_pha_val + else + if (CS%mld_pha_use_delta_rho) then + call diagnoseMLDbyDensityDifference(-1, h_old, tv, CS%mld_pha_drho, G, GV, US, CS%diag, CS%mld_pha_href, id_ref_z=-1, id_ref_rho=-1, MLD_out=mld_pha) + elseif (CS%mld_pha_use_delta_eng) then + call diagnoseMLDbyEnergy((/-1, -1, -1/), h_old, tv, G, GV, US, (/CS%mld_pha_deng, CS%mld_pha_deng, CS%mld_pha_deng/), CS%diag, MLD_out=mld_pha) + endif + endif + ! !Calculate tendencies (i.e., field changes at dt) from the sources / sinks ! @@ -729,7 +787,8 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & G%areaT, get_diag_time_end(CS%diag), & optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band, & - internal_heat=tv%internal_heat, frunoff=fluxes%frunoff, sosga=sosga, geolat=G%geolatT, eqn_of_state=tv%eqn_of_state) + internal_heat=tv%internal_heat, frunoff=fluxes%frunoff, sosga=sosga, geolat=G%geolatT, eqn_of_state=tv%eqn_of_state, & + photo_acc_dpth=mld_pha) else ! tv%internal_heat is a null pointer unless DO_GEOTHERMAL = True, ! so we have to check and only do the scaling if it is associated. @@ -740,14 +799,16 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, sw_pen_band=G%US%QRZ_T_to_W_m2*optics%sw_pen_band(:,:,:), & opacity_band=G%US%m_to_Z*optics%opacity_band(:,:,:,:), & internal_heat=G%US%RZ_to_kg_m2*US%C_to_degC*tv%internal_heat(:,:), & - frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga, geolat=G%geolatT, eqn_of_state=tv%eqn_of_state) + frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga, geolat=G%geolatT, eqn_of_state=tv%eqn_of_state, & + photo_acc_dpth=mld_pha*US%Z_to_m) else call generic_tracer_source(US%C_to_degC*tv%T, US%S_to_ppt*tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), & optics%nbands, optics%max_wavelength_band, & sw_pen_band=G%US%QRZ_T_to_W_m2*optics%sw_pen_band(:,:,:), & opacity_band=G%US%m_to_Z*optics%opacity_band(:,:,:,:), & - frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga, geolat=G%geolatT, eqn_of_state=tv%eqn_of_state) + frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga, geolat=G%geolatT, eqn_of_state=tv%eqn_of_state, & + photo_acc_dpth=mld_pha*US%Z_to_m) endif endif