NCAR to main candidate branch (2023-3-02)

config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90

@@ -53,7 +53,7 @@ module MOM_ocean_model_nuopc
 use MOM_coupler_types,       only : coupler_type_set_diags, coupler_type_send_data
 use mpp_domains_mod,         only : domain2d, mpp_get_layout, mpp_get_global_domain
 use mpp_domains_mod,         only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain
-use fms_mod,                 only : stdout
+use MOM_io,                  only : stdout
 use MOM_EOS,                 only : gsw_sp_from_sr, gsw_pt_from_ct
 use MOM_wave_interface,      only : wave_parameters_CS, MOM_wave_interface_init
 use MOM_wave_interface,      only : Update_Surface_Waves, query_wave_properties
@@ -446,7 +446,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
   call diag_mediator_close_registration(OS%diag)
 
   if (is_root_pe()) &
-    write(*,'(/12x,a/)') '======== COMPLETED MOM INITIALIZATION ========'
+    write(stdout,'(/12x,a/)') '======== COMPLETED MOM INITIALIZATION ========'
 
   call callTree_leave("ocean_model_init(")
 end subroutine ocean_model_init
@@ -1123,20 +1123,18 @@ subroutine ocean_public_type_chksum(id, timestep, ocn)
   ! Local variables
   integer(kind=int64) :: chks ! A checksum for the field
   logical :: root    ! True only on the root PE
-  integer :: outunit ! The output unit to write to
 
-  outunit = stdout()
   root = is_root_pe()
 
-  if (root) write(outunit,*) "BEGIN CHECKSUM(ocean_type):: ", id, timestep
-  chks = field_chksum(ocn%t_surf ) ; if (root) write(outunit,100) 'ocean%t_surf   ', chks
-  chks = field_chksum(ocn%s_surf ) ; if (root) write(outunit,100) 'ocean%s_surf   ', chks
-  chks = field_chksum(ocn%u_surf ) ; if (root) write(outunit,100) 'ocean%u_surf   ', chks
-  chks = field_chksum(ocn%v_surf ) ; if (root) write(outunit,100) 'ocean%v_surf   ', chks
-  chks = field_chksum(ocn%sea_lev) ; if (root) write(outunit,100) 'ocean%sea_lev  ', chks
-  chks = field_chksum(ocn%frazil ) ; if (root) write(outunit,100) 'ocean%frazil   ', chks
-  chks = field_chksum(ocn%melt_potential) ; if (root) write(outunit,100) 'ocean%melt_potential   ', chks
-  call coupler_type_write_chksums(ocn%fields, outunit, 'ocean%')
+  if (root) write(stdout,*) "BEGIN CHECKSUM(ocean_type):: ", id, timestep
+  chks = field_chksum(ocn%t_surf ) ; if (root) write(stdout,100) 'ocean%t_surf   ', chks
+  chks = field_chksum(ocn%s_surf ) ; if (root) write(stdout,100) 'ocean%s_surf   ', chks
+  chks = field_chksum(ocn%u_surf ) ; if (root) write(stdout,100) 'ocean%u_surf   ', chks
+  chks = field_chksum(ocn%v_surf ) ; if (root) write(stdout,100) 'ocean%v_surf   ', chks
+  chks = field_chksum(ocn%sea_lev) ; if (root) write(stdout,100) 'ocean%sea_lev  ', chks
+  chks = field_chksum(ocn%frazil ) ; if (root) write(stdout,100) 'ocean%frazil   ', chks
+  chks = field_chksum(ocn%melt_potential) ; if (root) write(stdout,100) 'ocean%melt_potential   ', chks
+  call coupler_type_write_chksums(ocn%fields, stdout, 'ocean%')
 100 FORMAT("   CHECKSUM::",A20," = ",Z20)
 
 end subroutine ocean_public_type_chksum

src/diagnostics/MOM_obsolete_params.F90

@@ -97,6 +97,8 @@ subroutine find_obsolete_params(param_file)
   ! This parameter is on the to-do list to be obsoleted.
   ! call obsolete_logical(param_file, "NEW_SPONGES", hint="Use INTERPOLATE_SPONGE_TIME_SPACE instead.")
 
+  call obsolete_logical(param_file, "SMOOTH_RI", hint="Instead use N_SMOOTH_RI.")
+
   ! Write the file version number to the model log.
   call log_version(param_file, mdl, version)
 

src/framework/_Diagnostics.dox

@@ -90,7 +90,7 @@ An arbitrary number of lines, one per diagnostic field:
   "average" or "mean" performs a time-average.
   "min" or "max" diagnose the minium or maxium over each time period.
 
-- `regional_section` : "none" means global output. A string of six space separated numbers, "lat_min, lat_max, lon_min, lon_max, vert_min, vert_max", limits the diagnostic to a region.
+- `regional_section` : "none" means global output. A string of six space separated numbers, "lon_min lon_max lat_min lat_max vert_min vert_max", limits the diagnostic to a region.
 
 - `packing` : Data representation in the file. 1 means "real*8", 2 means "real*4", 4 mean 16-bit integers, 8 means 1-byte.
 

src/framework/posix.F90

@@ -188,7 +188,7 @@ end subroutine longjmp_posix
 
   !> C interface to POSIX siglongjmp()
   !! Users should use the Fortran-defined siglongjmp() function.
-  subroutine siglongjmp_posix(env, val) bind(c, name="longjmp")
+  subroutine siglongjmp_posix(env, val) bind(c, name="siglongjmp")
     ! #include <setjmp.h>
     ! int siglongjmp(jmp_buf env, int val);
     import :: sigjmp_buf, c_int
@@ -360,6 +360,9 @@ function sigsetjmp_missing(env, savesigs) result(rc) bind(c)
   print '(a)', 'ERROR: sigsetjmp() is not implemented in this build.'
   print '(a)', 'Recompile with autoconf or -DSIGSETJMP_NAME=\"<symbol name>\".'
   error stop
+
+  ! NOTE: Compilers may expect a return value, even if it is unreachable
+  rc = -1
 end function sigsetjmp_missing
 
 end module posix

src/parameterizations/vertical/MOM_CVMix_shear.F90

@@ -31,22 +31,22 @@ module MOM_CVMix_shear
 type, public :: CVMix_shear_cs ! TODO: private
   logical :: use_LMD94                      !< Flags to use the LMD94 scheme
   logical :: use_PP81                       !< Flags to use Pacanowski and Philander (JPO 1981)
-  logical :: smooth_ri                      !< If true, smooth Ri using a 1-2-1 filter
+  integer :: n_smooth_ri                    !< Number of times to smooth Ri using a 1-2-1 filter
   real    :: Ri_zero                        !< LMD94 critical Richardson number
   real    :: Nu_zero                        !< LMD94 maximum interior diffusivity
   real    :: KPP_exp                        !< Exponent of unitless factor of diff.
                                             !! for KPP internal shear mixing scheme.
   real, allocatable, dimension(:,:,:) :: N2 !< Squared Brunt-Vaisala frequency [T-2 ~> s-2]
   real, allocatable, dimension(:,:,:) :: S2 !< Squared shear frequency [T-2 ~> s-2]
   real, allocatable, dimension(:,:,:) :: ri_grad !< Gradient Richardson number
-  real, allocatable, dimension(:,:,:) :: ri_grad_smooth !< Gradient Richardson number
-                                                        !! after smoothing
+  real, allocatable, dimension(:,:,:) :: ri_grad_orig !< Gradient Richardson number
+                                                        !! before smoothing
   character(10) :: Mix_Scheme               !< Mixing scheme name (string)
 
   type(diag_ctrl), pointer :: diag => NULL() !< Pointer to the diagnostics control structure
   !>@{ Diagnostic handles
   integer :: id_N2 = -1, id_S2 = -1, id_ri_grad = -1, id_kv = -1, id_kd = -1
-  integer :: id_ri_grad_smooth = -1
+  integer :: id_ri_grad_orig = -1
   !>@}
 
 end type CVMix_shear_cs
@@ -72,7 +72,7 @@ subroutine calculate_CVMix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS )
   type(CVMix_shear_cs),                       pointer     :: CS  !< The control structure returned by a previous
                                                                  !! call to CVMix_shear_init.
   ! Local variables
-  integer :: i, j, k, kk, km1
+  integer :: i, j, k, kk, km1, s
   real :: GoRho  ! Gravitational acceleration divided by density [Z T-2 R-1 ~> m4 s-2 kg-1]
   real :: pref   ! Interface pressures [R L2 T-2 ~> Pa]
   real :: DU, DV ! Velocity differences [L T-1 ~> m s-1]
@@ -85,7 +85,8 @@ subroutine calculate_CVMix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS )
   real, dimension(2*(GV%ke)) :: temp_1d ! A column of temperatures [C ~> degC]
   real, dimension(2*(GV%ke)) :: salt_1d ! A column of salinities [S ~> ppt]
   real, dimension(2*(GV%ke)) :: rho_1d  ! A column of densities at interface pressures [R ~> kg m-3]
-  real, dimension(GV%ke+1) :: Ri_Grad !< Gradient Richardson number [nondim]
+  real, dimension(GV%ke+1) :: Ri_Grad   !< Gradient Richardson number [nondim]
+  real, dimension(GV%ke+1) :: Ri_Grad_prev !< Gradient Richardson number before s.th smoothing iteration [nondim]
   real, dimension(GV%ke+1) :: Kvisc   !< Vertical viscosity at interfaces [m2 s-1]
   real, dimension(GV%ke+1) :: Kdiff   !< Diapycnal diffusivity at interfaces [m2 s-1]
   real :: epsln  !< Threshold to identify vanished layers [H ~> m or kg m-2]
@@ -145,27 +146,35 @@ subroutine calculate_CVMix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS )
 
       Ri_grad(GV%ke+1) = Ri_grad(GV%ke)
 
-      if (CS%id_ri_grad > 0) CS%ri_grad(i,j,:) = Ri_Grad(:)
+      if (CS%n_smooth_ri > 0) then
+
+        if (CS%id_ri_grad_orig > 0) CS%ri_grad_orig(i,j,:) = Ri_Grad(:)
 
-      if (CS%smooth_ri) then
         ! 1) fill Ri_grad in vanished layers with adjacent value
         do k = 2, GV%ke
           if (h(i,j,k) <= epsln) Ri_grad(k) = Ri_grad(k-1)
         enddo
 
         Ri_grad(GV%ke+1) = Ri_grad(GV%ke)
 
-        ! 2) vertically smooth Ri with 1-2-1 filter
-        dummy =  0.25 * Ri_grad(2)
-        Ri_grad(GV%ke+1) = Ri_grad(GV%ke)
-        do k = 3, GV%ke
-          Ri_Grad(k) = dummy + 0.5 * Ri_Grad(k) + 0.25 * Ri_grad(k+1)
-          dummy = 0.25 * Ri_grad(k)
+        do s=1,CS%n_smooth_ri
+
+          Ri_Grad_prev(:) = Ri_Grad(:)
+
+          ! 2) vertically smooth Ri with 1-2-1 filter
+          dummy =  0.25 * Ri_grad_prev(2)
+          do k = 3, GV%ke
+            Ri_Grad(k) = dummy + 0.5 * Ri_Grad_prev(k) + 0.25 * Ri_grad_prev(k+1)
+            dummy = 0.25 * Ri_grad(k)
+          enddo
         enddo
 
-        if (CS%id_ri_grad_smooth > 0) CS%ri_grad_smooth(i,j,:) = Ri_Grad(:)
+        Ri_grad(GV%ke+1) = Ri_grad(GV%ke)
+
       endif
 
+      if (CS%id_ri_grad > 0) CS%ri_grad(i,j,:) = Ri_Grad(:)
+
       do K=1,GV%ke+1
         Kvisc(K) = US%Z2_T_to_m2_s * kv(i,j,K)
         Kdiff(K) = US%Z2_T_to_m2_s * kd(i,j,K)
@@ -190,7 +199,7 @@ subroutine calculate_CVMix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS )
   if (CS%id_N2 > 0) call post_data(CS%id_N2, CS%N2, CS%diag)
   if (CS%id_S2 > 0) call post_data(CS%id_S2, CS%S2, CS%diag)
   if (CS%id_ri_grad > 0) call post_data(CS%id_ri_grad, CS%ri_grad, CS%diag)
-  if (CS%id_ri_grad_smooth > 0) call post_data(CS%id_ri_grad_smooth ,CS%ri_grad_smooth, CS%diag)
+  if (CS%id_ri_grad_orig > 0) call post_data(CS%id_ri_grad_orig ,CS%ri_grad_orig, CS%diag)
 
 end subroutine calculate_CVMix_shear
 
@@ -274,10 +283,10 @@ logical function CVMix_shear_init(Time, G, GV, US, param_file, diag, CS)
                  "Exponent of unitless factor of diffusivities, "// &
                  "for KPP internal shear mixing scheme." &
                  ,units="nondim", default=3.0)
-  call get_param(param_file, mdl, "SMOOTH_RI", CS%smooth_ri, &
-                 "If true, vertically smooth the Richardson "// &
-                 "number by applying a 1-2-1 filter once.", &
-                 default = .false.)
+  call get_param(param_file, mdl, "N_SMOOTH_RI", CS%n_smooth_ri, &
+                 "If > 0, vertically smooth the Richardson "// &
+                 "number by applying a 1-2-1 filter N_SMOOTH_RI times.", &
+                 default = 0)
   call cvmix_init_shear(mix_scheme=CS%Mix_Scheme, &
                         KPP_nu_zero=CS%Nu_Zero,   &
                         KPP_Ri_zero=CS%Ri_zero,   &
@@ -304,11 +313,14 @@ logical function CVMix_shear_init(Time, G, GV, US, param_file, diag, CS)
     allocate( CS%ri_grad( SZI_(G), SZJ_(G), SZK_(GV)+1 ), source=1.e8 )
   endif
 
-  CS%id_ri_grad_smooth = register_diag_field('ocean_model', 'ri_grad_shear_smooth', &
-       diag%axesTi, Time, &
-      'Smoothed gradient Richarson number used by MOM_CVMix_shear module','nondim')
-  if (CS%id_ri_grad_smooth > 0) then !Initialize w/ large Richardson value
-    allocate( CS%ri_grad_smooth( SZI_(G), SZJ_(G), SZK_(GV)+1 ), source=1.e8 )
+  if (CS%n_smooth_ri > 0) then
+    CS%id_ri_grad_orig = register_diag_field('ocean_model', 'ri_grad_shear_orig', &
+         diag%axesTi, Time, &
+        'Original gradient Richarson number, before smoothing was applied. This is '//&
+        'part of the MOM_CVMix_shear module and only available when N_SMOOTH_RI > 0','nondim')
+  endif
+  if (CS%id_ri_grad_orig > 0 .or. CS%n_smooth_ri > 0) then !Initialize w/ large Richardson value
+    allocate( CS%ri_grad_orig( SZI_(G), SZJ_(G), SZK_(GV)+1 ), source=1.e8 )
   endif
 
   CS%id_kd = register_diag_field('ocean_model', 'kd_shear_CVMix', diag%axesTi, Time, &

src/tracer/MOM_tracer_flow_control.F90

@@ -468,7 +468,8 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV,
       call ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
                                            G, GV, US, CS%ideal_age_tracer_CSp, &
                                            evap_CFL_limit=evap_CFL_limit, &
-                                           minimum_forcing_depth=minimum_forcing_depth)
+                                           minimum_forcing_depth=minimum_forcing_depth, &
+                                           Hbl=Hml)
     if (CS%use_regional_dyes) &
       call dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
                                      G, GV, US, CS%dye_tracer_CSp, &
@@ -544,7 +545,7 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV,
                                       G, GV, US, CS%RGC_tracer_CSp)
     if (CS%use_ideal_age) &
       call ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
-                                           G, GV, US, CS%ideal_age_tracer_CSp)
+                                           G, GV, US, CS%ideal_age_tracer_CSp, Hbl=Hml)
     if (CS%use_regional_dyes) &
       call dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
                                            G, GV, US, CS%dye_tracer_CSp)