Calculate a photoacclimation mixed layer depth in MOM6 (NOAA-GFDL#16)

* Calculate a photoacclimation MLD in MOM6

Using the diagnose MLD routines calculate a MLD that can be used
in COBALTv3.

* Update interfaces for diagnoseMLD

Some previously optional arguments were made nonoptinal in the
final version merged into dev/gfdl. So these areguents have been
added to the call in MOM_generic_tracers.

* Add unset variable in diagnoseMLDbyenergy

* Remove timescale variable

This parameter seems unneeded at this time.

* Change paramter name and though up PR

Change "MLD_PHA_FIXED" to "MLD_PHA_CALC" which is default false.
Touch up a few formatting issues.
Add initialization and comment for mld_pha.

* Add fatal errors

Add fatal errors if PHA_MLD_CALC is true, but a method is not specified.

* Add additional argument to generic tracer stub

Add the photo_acc_mld argument to the generic tracer stub in MOM6.

---------

Co-authored-by: Theresa Morrison <Theresa.Morrison@gaea57.ncrc.gov>
Co-authored-by: Theresa Morrison <Theresa.Morrison@gaea56.ncrc.gov>

Author: 3 people

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

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

src/tracer/MOM_generic_tracer.F90

@@ -31,6 +31,7 @@ module MOM_generic_tracer
 
   use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS
   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
@@ -83,6 +84,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
@@ -423,6 +432,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
@@ -507,6 +552,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
@@ -572,6 +619,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
     !
@@ -582,7 +640,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.
@@ -593,14 +652,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