Merge pull request mom-ocean#1026 from marshallward/nan_merge

NaN initialization fixes

src/core/MOM_dynamics_split_RK2.F90

@@ -1105,6 +1105,9 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
 !  Accel_diag%u_accel_bt => CS%u_accel_bt ; Accel_diag%v_accel_bt => CS%v_accel_bt
 !  Accel_diag%u_av => CS%u_av ; Accel_diag%v_av => CS%v_av
 
+  id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', &
+                                     grain=CLOCK_ROUTINE)
+
   call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
   call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv_CSp)
   if (use_tides) call tidal_forcing_init(Time, G, param_file, CS%tides_CSp)
@@ -1244,7 +1247,6 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   id_clock_horvisc    = cpu_clock_id('(Ocean horizontal viscosity)',     grain=CLOCK_MODULE)
   id_clock_mom_update = cpu_clock_id('(Ocean momentum increments)',      grain=CLOCK_MODULE)
   id_clock_pass       = cpu_clock_id('(Ocean message passing)',          grain=CLOCK_MODULE)
-  id_clock_pass_init  = cpu_clock_id('(Ocean init message passing)',     grain=CLOCK_ROUTINE)
   id_clock_btcalc     = cpu_clock_id('(Ocean barotropic mode calc)',     grain=CLOCK_MODULE)
   id_clock_btstep     = cpu_clock_id('(Ocean barotropic mode stepping)', grain=CLOCK_MODULE)
   id_clock_btforce    = cpu_clock_id('(Ocean barotropic forcing calc)',  grain=CLOCK_MODULE)

src/core/MOM_forcing_type.F90

@@ -150,12 +150,12 @@ module MOM_forcing_type
                                  !! or freezing (negative) [m year-1]
 
   ! Scalars set by surface forcing modules
-  real :: vPrecGlobalAdj     !< adjustment to restoring vprec to zero out global net [kg m-2 s-1]
-  real :: saltFluxGlobalAdj  !< adjustment to restoring salt flux to zero out global net [kgSalt m-2 s-1]
-  real :: netFWGlobalAdj     !< adjustment to net fresh water to zero out global net [kg m-2 s-1]
-  real :: vPrecGlobalScl     !< scaling of restoring vprec to zero out global net ( -1..1 ) [nondim]
-  real :: saltFluxGlobalScl  !< scaling of restoring salt flux to zero out global net ( -1..1 ) [nondim]
-  real :: netFWGlobalScl     !< scaling of net fresh water to zero out global net ( -1..1 ) [nondim]
+  real :: vPrecGlobalAdj = 0.     !< adjustment to restoring vprec to zero out global net [kg m-2 s-1]
+  real :: saltFluxGlobalAdj = 0.  !< adjustment to restoring salt flux to zero out global net [kgSalt m-2 s-1]
+  real :: netFWGlobalAdj = 0.     !< adjustment to net fresh water to zero out global net [kg m-2 s-1]
+  real :: vPrecGlobalScl = 0.     !< scaling of restoring vprec to zero out global net ( -1..1 ) [nondim]
+  real :: saltFluxGlobalScl = 0.  !< scaling of restoring salt flux to zero out global net ( -1..1 ) [nondim]
+  real :: netFWGlobalScl = 0.     !< scaling of net fresh water to zero out global net ( -1..1 ) [nondim]
 
   logical :: fluxes_used = .true. !< If true, all of the heat, salt, and mass
                                   !! fluxes have been applied to the ocean.

src/diagnostics/MOM_diagnostics.F90

@@ -829,15 +829,18 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
     call post_data(CS%id_col_ht, z_bot, CS%diag)
   endif
 
+  ! NOTE: int_density_z expects z_top and z_btm values from [ij]sq to [ij]eq+1
   if (CS%id_col_mass > 0 .or. CS%id_pbo > 0) then
     do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo
     if (GV%Boussinesq) then
       if (associated(tv%eqn_of_state)) then
         IG_Earth = 1.0 / GV%mks_g_Earth
 !       do j=js,je ; do i=is,ie ; z_bot(i,j) = -P_SURF(i,j)/GV%H_to_Pa ; enddo ; enddo
-        do j=js,je ; do i=is,ie ; z_bot(i,j) = 0.0 ; enddo ; enddo
+        do j=G%jscB,G%jecB+1 ; do i=G%iscB,G%iecB+1
+          z_bot(i,j) = 0.0
+        enddo ; enddo
         do k=1,nz
-          do j=js,je ; do i=is,ie
+          do j=G%jscB,G%jecB+1 ; do i=G%iscB,G%iecB+1
             z_top(i,j) = z_bot(i,j)
             z_bot(i,j) = z_top(i,j) - GV%H_to_Z*h(i,j,k)
           enddo ; enddo

src/diagnostics/MOM_wave_speed.F90

@@ -75,12 +75,14 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, &
     pres, T_int, S_int, &
     gprime        ! The reduced gravity across each interface [m2 Z-1 s-2 ~> m s-2].
   real, dimension(SZK_(G)) :: &
-    Igl, Igu      ! The inverse of the reduced gravity across an interface times
-                  ! the thickness of the layer below (Igl) or above (Igu) it [s2 m-2].
+    Igl, Igu, Igd ! The inverse of the reduced gravity across an interface times
+                  ! the thickness of the layer below (Igl) or above (Igu) it.
+                  ! Igd is provided for the tridiagonal solver.  [s2 m-2]
   real, dimension(SZK_(G),SZI_(G)) :: &
     Hf, Tf, Sf, Rf
   real, dimension(SZK_(G)) :: &
     Hc, Tc, Sc, Rc
+  real, dimension(SZK_(G)) :: Hc_H        ! Hc(:) rescaled from Z to thickness
   real det, ddet, detKm1, detKm2, ddetKm1, ddetKm2
   real :: lam, dlam, lam0
   real :: min_h_frac
@@ -145,8 +147,8 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, &
 !$OMP                                  Z_to_Pa,cg1,g_Rho0,rescale,I_rescale,L2_to_Z2)  &
 !$OMP                          private(htot,hmin,kf,H_here,HxT_here,HxS_here,HxR_here, &
 !$OMP                                  Hf,Tf,Sf,Rf,pres,T_int,S_int,drho_dT,           &
-!$OMP                                  drho_dS,drxh_sum,kc,Hc,Tc,Sc,I_Hnew,gprime,     &
-!$OMP                                  Rc,speed2_tot,Igl,Igu,lam0,lam,lam_it,dlam,     &
+!$OMP                                  drho_dS,drxh_sum,kc,Hc,Hc_H,Tc,Sc,I_Hnew,gprime,&
+!$OMP                                  Rc,speed2_tot,Igl,Igu,Igd,lam0,lam,lam_it,dlam, &
 !$OMP                                  mode_struct,sum_hc,N2min,gp,hw,                 &
 !$OMP                                  ms_min,ms_max,ms_sq,                            &
 !$OMP                                  det,ddet,detKm1,ddetKm1,detKm2,ddetKm2,det_it,ddet_it)
@@ -424,7 +426,10 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, &
             endif
 
             if (calc_modal_structure) then
-              call tdma6(kc, -igu, igu+igl, -igl, lam, mode_struct)
+              do k = 1,kc
+                Igd(k) = Igu(k) + Igl(k)
+              enddo
+              call tdma6(kc, -Igu, Igd, -Igl, lam, mode_struct)
               ms_min = mode_struct(1)
               ms_max = mode_struct(1)
               ms_sq = mode_struct(1)**2
@@ -456,8 +461,12 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, &
             endif
             ! Note that remapping_core_h requires that the same units be used
             ! for both the source and target grid thicknesses, here [H ~> m or kg m-2].
-            call remapping_core_h(CS%remapping_CS, kc, GV%Z_to_H*Hc(:), mode_struct, &
-                                  nz, h(i,j,:), modal_structure(i,j,:), 1.0e-30*GV%m_to_H, 1.0e-10*GV%m_to_H)
+            do k = 1,kc
+              Hc_H(k) = GV%Z_to_H * Hc(k)
+            enddo
+            call remapping_core_h(CS%remapping_CS, kc, Hc_H(:), mode_struct, &
+                                  nz, h(i,j,:), modal_structure(i,j,:), &
+                                  1.0e-30*GV%m_to_H, 1.0e-10*GV%m_to_H)
           endif
         else
           cg1(i,j) = 0.0

src/parameterizations/lateral/MOM_mixed_layer_restrat.F90

@@ -214,7 +214,9 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt_in_T, MLD_in
         ! Mixed-layer depth, using sigma-0 (surface reference pressure)
         deltaRhoAtKm1(:) = deltaRhoAtK(:) ! Store value from previous iteration of K
         call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, deltaRhoAtK, is-1, ie-is+3, tv%eqn_of_state)
-        deltaRhoAtK(:) = deltaRhoAtK(:) - rhoSurf(:) ! Density difference between layer K and surface
+        do i = is-1,ie+1
+          deltaRhoAtK(i) = deltaRhoAtK(i) - rhoSurf(i) ! Density difference between layer K and surface
+        enddo
         do i = is-1, ie+1
           ddRho = deltaRhoAtK(i) - deltaRhoAtKm1(i)
           if ((MLD_fast(i,j)==0.) .and. (ddRho>0.) .and. &

src/parameterizations/vertical/MOM_entrain_diffusive.F90

@@ -196,6 +196,7 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, &
   real :: Idt        ! The inverse of the time step [T-1 ~> s-1].
 
   logical :: do_any
+  logical :: do_entrain_eakb    ! True if buffer layer is entrained
   logical :: do_i(SZI_(G)), did_i(SZI_(G)), reiterate, correct_density
   integer :: it, i, j, k, is, ie, js, je, nz, K2, kmb
   integer :: kb(SZI_(G))  ! The value of kb in row j.
@@ -254,7 +255,7 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, &
   !$OMP                          private(dtKd,dtKd_int,do_i,Ent_bl,dtKd_kb,h_bl,     &
   !$OMP                                  I2p2dsp1_ds,grats,htot,max_eakb,I_dSkbp1,   &
   !$OMP                                  zeros,maxF_kb,maxF,ea_kbp1,eakb,Sref,       &
-  !$OMP                                  maxF_correct,do_any,                        &
+  !$OMP                                  maxF_correct,do_any,do_entrain_eakb,        &
   !$OMP                                  err_min_eakb0,err_max_eakb0,eakb_maxF,      &
   !$OMP                                  min_eakb,err_eakb0,F,minF,hm,fk,F_kb_maxent,&
   !$OMP                                  F_kb,is1,ie1,kb_min_act,dFdfm_kb,b1,dFdfm,  &
@@ -355,10 +356,16 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, &
                         kmb, is, ie, G, GV, CS, F_kb_maxEnt, do_i_in = do_i)
 
       do i=is,ie
-        if ((.not.do_i(i)) .or. (err_max_eakb0(i) >= 0.0)) then
-          eakb(i) = 0.0 ; min_eakb(i) = 0.0
-        else ! If error_max_eakb0 < 0 the buffer layers are always all entrained.
+        do_entrain_eakb = .false.
+        ! If error_max_eakb0 < 0, then buffer layers are always all entrained
+        if (do_i(i)) then ; if (err_max_eakb0(i) < 0.0) then
+          do_entrain_eakb = .true.
+        endif ; endif
+
+        if (do_entrain_eakb) then
           eakb(i) = max_eakb(i) ; min_eakb(i) = max_eakb(i)
+        else
+          eakb(i) = 0.0 ; min_eakb(i) = 0.0
         endif
       enddo
 
@@ -413,11 +420,13 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, &
       endif
     endif
     do k=nz-1,kb_min,-1 ; do i=is,ie ; if (do_i(i)) then
-      if (k>=kb(i)) then
+      if (k >= kb(i)) then
         maxF(i,k) = MIN(maxF(i,k),dsp1_ds(i,k+1)*maxF(i,k+1) + htot(i))
         htot(i) = htot(i) + (h(i,j,k) - Angstrom)
-        if ( (k == kb(i)) .and. ((maxF(i,k) < F_kb(i)) .or. &
-            (maxF(i,k) < maxF_kb(i)) .and. (eakb_maxF(i) <= max_eakb(i))) ) then
+      endif
+      if (k == kb(i)) then
+        if ((maxF(i,k) < F_kb(i)) .or. (maxF(i,k) < maxF_kb(i)) &
+            .and. (eakb_maxF(i) <= max_eakb(i))) then
           !   In this case, too much was being entrained by the topmost interior
           ! layer, even with the minimum initial estimate.  The buffer layer
           ! will always entrain the maximum amount.

src/parameterizations/vertical/MOM_opacity.F90

@@ -841,7 +841,12 @@ subroutine sumSWoverBands(G, GV, US, h, nsw, optics, j, dt, &
 
   pen_SW_bnd(:,:) = iPen_SW_bnd(:,:)
   do i=is,ie ; h_heat(i) = 0.0 ; enddo
-  netPen(:,1) = sum( pen_SW_bnd(:,:), dim=1 ) ! Surface interface
+  do i=is,ie
+    netPen(i,1) = 0.
+    do n=1,max(nsw,1)
+      netPen(i,1) = netPen(i,1) + pen_SW_bnd(n,i)   ! Surface interface
+    enddo
+  enddo
 
   ! Apply penetrating SW radiation to remaining parts of layers.
   ! Excessively thin layers are not heated to avoid runaway temps.

src/parameterizations/vertical/MOM_set_viscosity.F90

@@ -260,6 +260,9 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize)
   real :: C2pi_3           ! An irrational constant, 2/3 pi.
   real :: tmp              ! A temporary variable.
   real :: tmp_val_m1_to_p1
+  real :: curv_tol         ! Numerator of curvature cubed, used to estimate
+                           ! accuracy of a single L(:) Newton iteration
+  logical :: use_L0, do_one_L_iter    ! Control flags for L(:) Newton iteration
   logical :: use_BBL_EOS, do_i(SZIB_(G))
   integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m, n, K2, nkmb, nkml
   integer :: itt, maxitt=20
@@ -773,19 +776,29 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize)
               dV_dL2 = 0.5*(slope+a) - a*L0 ; dVol = (vol-Vol_0)
            !  dV_dL2 = 0.5*(slope+a) - a*L0 ; dVol = max(vol-Vol_0, 0.0)
 
-             ! The following code is more robust when GV%Angstrom_H=0, but it changes answers.
-             if (.not.CS%answers_2018) then
-               Vol_tol = max(0.5*GV%Angstrom_H + GV%H_subroundoff, 1e-14*vol)
-               Vol_quit = max(0.9*GV%Angstrom_H + GV%H_subroundoff, 1e-14*vol)
-             endif
+              use_L0 = .false.
+              do_one_L_iter = .false.
+              if (CS%answers_2018) then
+                curv_tol = GV%Angstrom_H*dV_dL2**2 &
+                           * (0.25 * dV_dL2 * GV%Angstrom_H - a * L0 * dVol)
+                do_one_L_iter = (a * a * dVol**3) < curv_tol
+              else
+                ! The following code is more robust when GV%Angstrom_H=0, but
+                ! it changes answers.
+                use_L0 = (dVol <= 0.)
+
+                Vol_tol = max(0.5 * GV%Angstrom_H + GV%H_subroundoff, 1e-14 * vol)
+                Vol_quit = max(0.9 * GV%Angstrom_H + GV%H_subroundoff, 1e-14 * vol)
+
+                curv_tol = Vol_tol * dV_dL2**2 &
+                           * (dV_dL2 * Vol_tol - 2.0 * a * L0 * dVol)
+                do_one_L_iter = (a * a * dVol**3) < curv_tol
+              endif
 
-              if ((.not.CS%answers_2018) .and. (dVol <= 0.0)) then
+              if (use_L0) then
                 L(K) = L0
                 Vol_err = 0.5*(L(K)*L(K))*(slope + a_3*(3.0-4.0*L(K))) - vol
-              elseif ( ((.not.CS%answers_2018) .and. &
-                  (a*a*dVol**3 < Vol_tol*dV_dL2**2 *(dV_dL2*Vol_tol - 2.0*a*L0*dVol))) .or. &
-                  (CS%answers_2018 .and. (a*a*dVol**3 < GV%Angstrom_H*dV_dL2**2 * &
-                                          (0.25*dV_dL2*GV%Angstrom_H - a*L0*dVol) )) ) then
+              elseif (do_one_L_iter) then
                 ! One iteration of Newton's method should give an estimate
                 ! that is accurate to within Vol_tol.
                 L(K) = sqrt(L0*L0 + dVol / dV_dL2)

src/tracer/MOM_neutral_diffusion.F90

@@ -2401,6 +2401,9 @@ logical function ndiff_unit_tests_discontinuous(verbose)
   ! Tests for linearized version of searching the layer for neutral surface position
   ! EOS linear in T, uniform alpha
   CS%max_iter = 10
+  ! Unit tests require explicit initialization of tolerance
+  CS%Drho_tol = 0.
+  CS%x_tol = 0.
   ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,        &
              find_neutral_pos_linear(CS, 0., 10., 35., 0., -0.2, 0.,                      &
                                      0., -0.2, 0., 10., -0.2, 0.,                     &

src/tracer/MOM_tracer_advect.F90

@@ -652,10 +652,14 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
     do m=1,ntr
 
       ! update tracer
-      do i=is,ie ; if ((do_i(i)) .and. (Ihnew(i) > 0.0)) then
-        Tr(m)%t(i,j,k) = (Tr(m)%t(i,j,k) * hlst(i) - &
-                          (flux_x(I,m) - flux_x(I-1,m))) * Ihnew(i)
-      endif ; enddo
+      do i=is,ie
+        if (do_i(i)) then
+          if (Ihnew(i) > 0.0) then
+            Tr(m)%t(i,j,k) = (Tr(m)%t(i,j,k) * hlst(i) - &
+                              (flux_x(I,m) - flux_x(I-1,m))) * Ihnew(i)
+          endif
+        endif
+      enddo
 
       ! diagnostics
       if (associated(Tr(m)%ad_x)) then ; do i=is,ie ; if (do_i(i)) then