Merge pull request #80 from gustavo-marques/pass_OBLD_caps

Makes OBLD available to MCT (and other caps)

config_src/mct_driver/MOM_ocean_model.F90

@@ -140,7 +140,8 @@ module MOM_ocean_model
     frazil =>NULL(), &   !< Accumulated heating (in Joules/m^2) from frazil
                          !! formation in the ocean.
     melt_potential => NULL(), & !< Accumulated heat used to melt sea ice (in W/m^2)
-    area => NULL()              !< cell area of the ocean surface, in m2.
+    area => NULL(), &    !< cell area of the ocean surface, in m2.
+    OBLD => NULL()       !< Ocean boundary layer depth, in m.
   type(coupler_2d_bc_type) :: fields    !< A structure that may contain an
                                         !! array of named tracer-related fields.
   integer                  :: avg_kount !< Used for accumulating averages of this type.
@@ -732,6 +733,7 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap,
             Ocean_sfc%v_surf (isc:iec,jsc:jec), &
             Ocean_sfc%sea_lev(isc:iec,jsc:jec), &
             Ocean_sfc%area   (isc:iec,jsc:jec), &
+            Ocean_sfc%OBLD   (isc:iec,jsc:jec), &
             Ocean_sfc%melt_potential(isc:iec,jsc:jec), &
             Ocean_sfc%frazil (isc:iec,jsc:jec))
 
@@ -742,6 +744,7 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap,
   Ocean_sfc%sea_lev = 0.0  ! time averaged thickness of top model grid cell (m) plus patm/rho0/grav
   Ocean_sfc%frazil  = 0.0  ! time accumulated frazil (J/m^2) passed to ice model
   Ocean_sfc%melt_potential  = 0.0  ! time accumulated melt potential (J/m^2) passed to ice model
+  Ocean_sfc%OBLD    = 0.0  ! ocean boundary layer depth, in m
   Ocean_sfc%area    = 0.0
   Ocean_sfc%axes    = diag%axesT1%handles !diag axes to be used by coupler tracer flux diagnostics
 
@@ -812,6 +815,8 @@ subroutine convert_state_to_ocean_type(state, Ocean_sfc, G, patm, press_to_z)
       Ocean_sfc%frazil(i,j) = state%frazil(i+i0,j+j0)
     if (allocated(state%melt_potential)) &
       Ocean_sfc%melt_potential(i,j) = state%melt_potential(i+i0,j+j0)
+    if (allocated(state%Hml)) &
+      Ocean_sfc%OBLD(i,j) = state%Hml(i+i0,j+j0)
   enddo ; enddo
 
   if (Ocean_sfc%stagger == AGRID) then
@@ -1036,6 +1041,7 @@ subroutine ocean_public_type_chksum(id, timestep, ocn)
   write(outunit,100) 'ocean%v_surf   ',mpp_chksum(ocn%v_surf )
   write(outunit,100) 'ocean%sea_lev  ',mpp_chksum(ocn%sea_lev)
   write(outunit,100) 'ocean%frazil   ',mpp_chksum(ocn%frazil )
+  write(outunit,100) 'ocean%OBLD     ',mpp_chksum(ocn%OBLD   )
   write(outunit,100) 'ocean%melt_potential  ',mpp_chksum(ocn%melt_potential)
 
   call coupler_type_write_chksums(ocn%fields, outunit, 'ocean%')

config_src/mct_driver/ocn_cap_methods.F90

@@ -176,6 +176,7 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
         o2x(ind%o2x_So_s, n) = ocn_public%s_surf(ig,jg) * grid%mask2dT(i,j)
         o2x(ind%o2x_So_u, n) = ocn_public%u_surf(ig,jg) * grid%mask2dT(i,j)
         o2x(ind%o2x_So_v, n) = ocn_public%v_surf(ig,jg) * grid%mask2dT(i,j)
+        o2x(ind%o2x_So_bldepth, n) = ocn_public%OBLD(ig,jg) * grid%mask2dT(i,j)
         ! ocean melt and freeze potential (o2x_Fioo_q), W m-2
         if (ocn_public%frazil(ig,jg) > 0.0) then
           ! Frazil: change from J/m^2 to W/m^2

src/core/MOM.F90

@@ -548,7 +548,6 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, &
 
   if (therm_reset) then
     CS%time_in_thermo_cycle = 0.0
-    if (allocated(sfc_state%melt_potential)) sfc_state%melt_potential(:,:)  = 0.0
     if (associated(CS%tv%frazil))        CS%tv%frazil(:,:)        = 0.0
     if (associated(CS%tv%salt_deficit))  CS%tv%salt_deficit(:,:)  = 0.0
     if (associated(CS%tv%TempxPmE))      CS%tv%TempxPmE(:,:)      = 0.0
@@ -2717,6 +2716,13 @@ subroutine extract_surface_state(CS, sfc_state)
     sfc_state%sea_lev(i,j) = CS%ave_ssh_ibc(i,j)
   enddo ; enddo
 
+  ! copy Hml into sfc_state, so that caps can access it
+  if (associated(CS%Hml)) then
+    do j=js,je ; do i=is,ie
+      sfc_state%Hml(i,j) = CS%Hml(i,j)
+    enddo ; enddo
+  endif
+
   if (CS%Hmix < 0.0) then  ! A bulk mixed layer is in use, so layer 1 has the properties
     if (use_temperature) then ; do j=js,je ; do i=is,ie
       sfc_state%SST(i,j) = CS%tv%T(i,j,1)
@@ -2729,9 +2735,6 @@ subroutine extract_surface_state(CS, sfc_state)
       sfc_state%v(i,J) = v(i,J,1)
     enddo ; enddo
 
-    if (associated(CS%Hml)) then ; do j=js,je ; do i=is,ie
-      sfc_state%Hml(i,j) = CS%Hml(i,j)
-    enddo ; enddo ; endif
   else  ! (CS%Hmix >= 0.0)
 
     depth_ml = CS%Hmix
@@ -2773,7 +2776,6 @@ subroutine extract_surface_state(CS, sfc_state)
         else
           sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) / depth(i)
         endif
-        sfc_state%Hml(i,j) = depth(i)
       enddo
     enddo ! end of j loop
 
@@ -2867,12 +2869,13 @@ subroutine extract_surface_state(CS, sfc_state)
       enddo ; enddo
 
       do i=is,ie
+       ! set melt_potential to zero to avoid passing previous values
+       sfc_state%melt_potential(i,j) = 0.0
+
        if (G%mask2dT(i,j)>0.) then
-         ! time accumulated melt_potential, in J/m^2
-         sfc_state%melt_potential(i,j) = sfc_state%melt_potential(i,j) +  (CS%tv%C_p * CS%GV%Rho0 * delT(i))
-       else
-         sfc_state%melt_potential(i,j) = 0.0
-       endif! G%mask2dT
+         ! instantaneous melt_potential, in J/m^2
+         sfc_state%melt_potential(i,j) = CS%tv%C_p * CS%GV%Rho0 * delT(i)
+       endif
       enddo
     enddo ! end of j loop
   endif   ! melt_potential

src/core/MOM_variables.F90

@@ -40,8 +40,8 @@ module MOM_variables
     v, &        !< The mixed layer meridional velocity in m s-1.
     sea_lev, &  !< The sea level in m.  If a reduced surface gravity is
                 !! used, that is compensated for in sea_lev.
-    melt_potential, & !< Amount of heat that can be used to melt sea ice, in J m-2.
-                      !! This is computed w.r.t. surface freezing temperature.
+    melt_potential, & !< instantaneous amount of heat that can be used to melt sea ice,
+                      !! in J m-2. This is computed w.r.t. surface freezing temperature.
     ocean_mass, &  !< The total mass of the ocean in kg m-2.
     ocean_heat, &  !< The total heat content of the ocean in C kg m-2.
     ocean_salt, &  !< The total salt content of the ocean in kgSalt m-2.

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -92,6 +92,8 @@ module MOM_diabatic_driver
   logical :: use_energetic_PBL       !< If true, use the implicit energetics planetary
                                      !! boundary layer scheme to determine the diffusivity
                                      !! in the surface boundary layer.
+  logical :: use_KPP                 !< If true, use CVMix/KPP boundary layer scheme to determine the
+                                     !! OBLD and the diffusivities within this layer.
   logical :: use_kappa_shear         !< If true, use the kappa_shear module to find the
                                      !! shear-driven diapycnal diffusivity.
   logical :: use_CVMix_shear         !< If true, use the CVMix module to find the
@@ -266,7 +268,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h         !< thickness (m for Bouss / kg/m2 for non-Bouss)
   type(thermo_var_ptrs),                     intent(inout) :: tv        !< points to thermodynamic fields
                                                                         !! unused have NULL ptrs
-  real, dimension(:,:),                      pointer       :: Hml       !< active mixed layer depth
+  real, dimension(:,:),                      pointer       :: Hml       !< mixed layer depth, m
   type(forcing),                             intent(inout) :: fluxes    !< points to forcing fields
                                                                         !! unused fields have NULL ptrs
   type(vertvisc_type),                       intent(inout) :: visc      !< vertical viscosities, BBL properies, and
@@ -736,11 +738,11 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
     call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, &
          CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves)
 
-    ! If visc%MLD exists, copy the ePBL's MLD into it
-    if (associated(visc%MLD)) then
-      call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G)
-      call pass_var(visc%MLD, G%domain, halo=1)
-      Hml(:,:) = visc%MLD(:,:)
+    if (associated(Hml)) then
+      call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G)
+      call pass_var(Hml, G%domain, halo=1)
+      ! If visc%MLD exists, copy KPP's BLD into it
+      if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:)
     endif
 
     ! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL.
@@ -2815,7 +2817,10 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag,
   call get_param(param_file, mod, "DOUBLE_DIFFUSION", differentialDiffusion, &
                  "If true, apply parameterization of double-diffusion.", &
                  default=.false. )
-
+  call get_param(param_file, mod, "USE_KPP", CS%use_KPP, &
+                 "If true, turns on the [CVMix] KPP scheme of Large et al., 1994,\n"// &
+                 "to calculate diffusivities and non-local transport in the OBL.",     &
+                 default=.false., do_not_log=.true.)
   CS%use_CVMix_ddiff = CVMix_ddiff_is_used(param_file)
 
   if (CS%use_CVMix_ddiff .and. differentialDiffusion) then