Merge branch 'dev/gfdl' into dev/esmg

.testing/Makefile

@@ -310,7 +310,7 @@ $(eval $(call STAT_RULE,symmetric,symmetric,$(REPORT_COVERAGE),,,1))
 $(eval $(call STAT_RULE,asymmetric,asymmetric,,,,1))
 $(eval $(call STAT_RULE,target,target,,,,1))
 $(eval $(call STAT_RULE,repro,repro,,,,1))
-$(eval $(call STAT_RULE,openmp,openmp,,,,1))
+$(eval $(call STAT_RULE,openmp,openmp,,,GOMP_CPU_AFFINITY=0,1))
 $(eval $(call STAT_RULE,layout,symmetric,,LAYOUT=2$(,)1,,2))
 $(eval $(call STAT_RULE,rotate,symmetric,,ROTATE_INDEX=True\nINDEX_TURNS=1,,1))
 $(eval $(call STAT_RULE,nan,symmetric,,,MALLOC_PERTURB_=256,1))

src/diagnostics/MOM_diagnostics.F90

@@ -15,7 +15,7 @@ module MOM_diagnostics
 use MOM_diag_mediator,     only : diag_save_grids, diag_restore_grids, diag_copy_storage_to_diag
 use MOM_domains,           only : create_group_pass, do_group_pass, group_pass_type
 use MOM_domains,           only : To_North, To_East
-use MOM_EOS,               only : calculate_density, int_density_dz, EOS_domain
+use MOM_EOS,               only : calculate_density, calculate_density_derivs, int_density_dz, EOS_domain
 use MOM_EOS,               only : gsw_sp_from_sr, gsw_pt_from_ct
 use MOM_error_handler,     only : MOM_error, FATAL, WARNING
 use MOM_file_parser,       only : get_param, log_version, param_file_type
@@ -134,6 +134,7 @@ module MOM_diagnostics
   integer :: id_pbo            = -1
   integer :: id_thkcello       = -1, id_rhoinsitu      = -1
   integer :: id_rhopot0        = -1, id_rhopot2        = -1
+  integer :: id_drho_dT        = -1, id_drho_dS        = -1
   integer :: id_h_pre_sync     = -1
   !>@}
   !> The control structure for calculating wave speed.
@@ -619,6 +620,22 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
       enddo
       if (CS%id_rhoinsitu > 0) call post_data(CS%id_rhoinsitu, Rcv, CS%diag)
     endif
+
+    if (CS%id_drho_dT > 0 .or. CS%id_drho_dS > 0) then
+      !$OMP parallel do default(shared) private(pressure_1d)
+      do j=js,je
+        pressure_1d(:) = 0. ! Start at p=0 Pa at surface
+        do k=1,nz
+          pressure_1d(:) =  pressure_1d(:) + 0.5 * h(:,j,k) * GV%H_to_Pa ! Pressure in middle of layer k
+          ! To avoid storing more arrays, put drho_dT into Rcv, and drho_dS into work3d
+          call calculate_density_derivs(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
+                                 Rcv(:,j,k),work_3d(:,j,k),is,ie-is+1, tv%eqn_of_state)
+          pressure_1d(:) =  pressure_1d(:) + 0.5 * h(:,j,k) * GV%H_to_Pa ! Pressure at bottom of layer k
+        enddo
+      enddo
+      if (CS%id_drho_dT > 0) call post_data(CS%id_drho_dT, Rcv, CS%diag)
+      if (CS%id_drho_dS > 0) call post_data(CS%id_drho_dS, work_3d, CS%diag)
+    endif
   endif
 
   if ((CS%id_cg1>0) .or. (CS%id_Rd1>0) .or. (CS%id_cfl_cg1>0) .or. &
@@ -1600,6 +1617,10 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
       'Potential density referenced to 2000 dbar', 'kg m-3', conversion=US%R_to_kg_m3)
   CS%id_rhoinsitu = register_diag_field('ocean_model', 'rhoinsitu', diag%axesTL, Time, &
       'In situ density', 'kg m-3', conversion=US%R_to_kg_m3)
+  CS%id_drho_dT = register_diag_field('ocean_model', 'drho_dT', diag%axesTL, Time, &
+      'Partial derivative of rhoinsitu with respect to temperature (alpha)', 'kg m-3 degC-1')
+  CS%id_drho_dS = register_diag_field('ocean_model', 'drho_dS', diag%axesTL, Time, &
+      'Partial derivative of rhoinsitu with respect to salinity (beta)', 'kg^2 g-1 m-3')
 
   CS%id_du_dt = register_diag_field('ocean_model', 'dudt', diag%axesCuL, Time, &
       'Zonal Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2)

src/parameterizations/vertical/MOM_CVMix_KPP.F90

@@ -627,15 +627,13 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, &
   real :: LangEnhK     ! Langmuir enhancement for mixing coefficient
 
 
-#ifdef __DO_SAFETY_CHECKS__
   if (CS%debug) then
     call hchksum(h, "KPP in: h",G%HI,haloshift=0, scale=GV%H_to_m)
     call hchksum(uStar, "KPP in: uStar",G%HI,haloshift=0, scale=US%Z_to_m*US%s_to_T)
     call hchksum(buoyFlux, "KPP in: buoyFlux",G%HI,haloshift=0)
     call hchksum(Kt, "KPP in: Kt",G%HI,haloshift=0, scale=US%Z2_T_to_m2_s)
     call hchksum(Ks, "KPP in: Ks",G%HI,haloshift=0, scale=US%Z2_T_to_m2_s)
   endif
-#endif
 
   nonLocalTrans(:,:) = 0.0
 
@@ -862,12 +860,10 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, &
   enddo ! j
 
 
-#ifdef __DO_SAFETY_CHECKS__
   if (CS%debug) then
     call hchksum(Kt, "KPP out: Kt", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
     call hchksum(Ks, "KPP out: Ks", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
   endif
-#endif
 
   ! send diagnostics to post_data
   if (CS%id_OBLdepth > 0) call post_data(CS%id_OBLdepth, CS%OBLdepth,        CS%diag)
@@ -952,14 +948,12 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
   real :: WST
 
 
-#ifdef __DO_SAFETY_CHECKS__
   if (CS%debug) then
     call hchksum(Salt, "KPP in: S",G%HI,haloshift=0)
     call hchksum(Temp, "KPP in: T",G%HI,haloshift=0)
     call hchksum(u, "KPP in: u",G%HI,haloshift=0)
     call hchksum(v, "KPP in: v",G%HI,haloshift=0)
   endif
-#endif
 
   ! some constants
   GoRho = US%L_T_to_m_s**2*US%m_to_Z * GV%g_Earth / GV%Rho0