Merge pull request NOAA-EMC#20 from NOAA-GFDL/dev/gfdl

merge in latest dev/gfdl changes

src/core/MOM.F90

@@ -2748,10 +2748,10 @@ subroutine extract_surface_state(CS, sfc_state)
       sfc_state%SST(i,j) = CS%tv%T(i,j,1)
       sfc_state%SSS(i,j) = CS%tv%S(i,j,1)
     enddo ; enddo ; endif
-    do j=js,je ; do I=IscB,IecB
+    do j=js,je ; do I=is-1,ie
       sfc_state%u(I,j) = u(I,j,1)
     enddo ; enddo
-    do J=JscB,JecB ; do i=is,ie
+    do J=js-1,je ; do i=is,ie
       sfc_state%v(i,J) = v(i,J,1)
     enddo ; enddo
 
@@ -2803,12 +2803,15 @@ subroutine extract_surface_state(CS, sfc_state)
     enddo ! end of j loop
 
 !   Determine the mean velocities in the uppermost depth_ml fluid.
+    ! NOTE: Velocity loops start on `[ij]s-1` in order to update halo values
+    !       required by the speed diagnostic on the non-symmetric grid.
+    !       This assumes that u and v halos have already been updated.
     if (CS%Hmix_UV>0.) then
       !### This calculation should work in thickness (H) units instead of Z, but that
       !### would change answers at roundoff in non-Boussinesq cases.
       depth_ml = CS%Hmix_UV
       !$OMP parallel do default(shared) private(depth,dh,hv)
-      do J=jscB,jecB
+      do J=js-1,ie
         do i=is,ie
           depth(i) = 0.0
           sfc_state%v(i,J) = 0.0
@@ -2835,11 +2838,11 @@ subroutine extract_surface_state(CS, sfc_state)
 
       !$OMP parallel do default(shared) private(depth,dh,hu)
       do j=js,je
-        do I=iscB,iecB
+        do I=is-1,ie
           depth(I) = 0.0
           sfc_state%u(I,j) = 0.0
         enddo
-        do k=1,nz ; do I=iscB,iecB
+        do k=1,nz ; do I=is-1,ie
           hu = 0.5 * (h(i,j,k) + h(i+1,j,k)) * GV%H_to_Z
           if (depth(i) + hu < depth_ml) then
             dh = hu
@@ -2852,17 +2855,17 @@ subroutine extract_surface_state(CS, sfc_state)
           depth(I) = depth(I) + dh
         enddo ; enddo
         ! Calculate the average properties of the mixed layer depth.
-        do I=iscB,iecB
+        do I=is-1,ie
           if (depth(I) < GV%H_subroundoff*GV%H_to_Z) &
               depth(I) = GV%H_subroundoff*GV%H_to_Z
           sfc_state%u(I,j) = sfc_state%u(I,j) / depth(I)
         enddo
       enddo ! end of j loop
     else ! Hmix_UV<=0.
-      do j=js,je ; do I=IscB,IecB
+      do j=js,je ; do I=is-1,ie
         sfc_state%u(I,j) = u(I,j,1)
       enddo ; enddo
-      do J=JscB,JecB ; do i=is,ie
+      do J=js-1,je ; do i=is,ie
         sfc_state%v(i,J) = v(i,J,1)
       enddo ; enddo
     endif

src/diagnostics/MOM_diagnostics.F90

@@ -1076,16 +1076,21 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, CS)
 end subroutine calculate_energy_diagnostics
 
 !> This subroutine registers fields to calculate a diagnostic time derivative.
-subroutine register_time_deriv(f_ptr, deriv_ptr, CS)
-  real, dimension(:,:,:), target :: f_ptr     !< Field whose derivative is taken.
-  real, dimension(:,:,:), target :: deriv_ptr !< Field in which the calculated time derivatives
-                                              !! will be placed.
+subroutine register_time_deriv(lb, f_ptr, deriv_ptr, CS)
+  integer, intent(in), dimension(3) :: lb     !< Lower index bound of f_ptr
+  real, dimension(lb(1):,lb(2):,:), target :: f_ptr
+                                              !< Time derivative operand
+  real, dimension(lb(1):,lb(2):,:), target :: deriv_ptr
+                                              !< Time derivative of f_ptr
   type(diagnostics_CS),  pointer :: CS        !< Control structure returned by previous call to
                                               !! diagnostics_init.
 
   ! This subroutine registers fields to calculate a diagnostic time derivative.
+  ! NOTE: Lower bound is required for grid indexing in calculate_derivs().
+  !       We assume that the vertical axis is 1-indexed.
 
-  integer :: m
+  integer :: m      !< New index of deriv_ptr in CS%deriv
+  integer :: ub(3)  !< Upper index bound of f_ptr, based on shape.
 
   if (.not.associated(CS)) call MOM_error(FATAL, &
          "register_time_deriv: Module must be initialized before it is used.")
@@ -1098,9 +1103,11 @@ subroutine register_time_deriv(f_ptr, deriv_ptr, CS)
 
   m = CS%num_time_deriv+1 ; CS%num_time_deriv = m
 
-  CS%nlay(m) = size(f_ptr(:,:,:),3)
+  ub(:) = lb(:) + shape(f_ptr) - 1
+
+  CS%nlay(m) = size(f_ptr, 3)
   CS%deriv(m)%p => deriv_ptr
-  allocate(CS%prev_val(m)%p(size(f_ptr(:,:,:),1), size(f_ptr(:,:,:),2), CS%nlay(m)) )
+  allocate(CS%prev_val(m)%p(lb(1):ub(1), lb(2):ub(2), CS%nlay(m)))
 
   CS%var_ptr(m)%p => f_ptr
   CS%prev_val(m)%p(:,:,:) = f_ptr(:,:,:)
@@ -1551,21 +1558,21 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
       'Zonal Acceleration', 'm s-2')
   if ((CS%id_du_dt>0) .and. .not.associated(CS%du_dt)) then
     call safe_alloc_ptr(CS%du_dt,IsdB,IedB,jsd,jed,nz)
-    call register_time_deriv(MIS%u, CS%du_dt, CS)
+    call register_time_deriv(lbound(MIS%u), MIS%u, CS%du_dt, CS)
   endif
 
   CS%id_dv_dt = register_diag_field('ocean_model', 'dvdt', diag%axesCvL, Time, &
       'Meridional Acceleration', 'm s-2')
   if ((CS%id_dv_dt>0) .and. .not.associated(CS%dv_dt)) then
     call safe_alloc_ptr(CS%dv_dt,isd,ied,JsdB,JedB,nz)
-    call register_time_deriv(MIS%v, CS%dv_dt, CS)
+    call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS)
   endif
 
   CS%id_dh_dt = register_diag_field('ocean_model', 'dhdt', diag%axesTL, Time, &
       'Thickness tendency', trim(thickness_units)//" s-1", v_extensive = .true.)
   if ((CS%id_dh_dt>0) .and. .not.associated(CS%dh_dt)) then
     call safe_alloc_ptr(CS%dh_dt,isd,ied,jsd,jed,nz)
-    call register_time_deriv(MIS%h, CS%dh_dt, CS)
+    call register_time_deriv(lbound(MIS%h), MIS%h, CS%dh_dt, CS)
   endif
 
   ! layer thickness variables
@@ -2009,15 +2016,15 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, CS)
   if (associated(CS%dKE_dt)) then
     if (.not.associated(CS%du_dt)) then
       call safe_alloc_ptr(CS%du_dt,IsdB,IedB,jsd,jed,nz)
-      call register_time_deriv(MIS%u, CS%du_dt, CS)
+      call register_time_deriv(lbound(MIS%u), MIS%u, CS%du_dt, CS)
     endif
     if (.not.associated(CS%dv_dt)) then
       call safe_alloc_ptr(CS%dv_dt,isd,ied,JsdB,JedB,nz)
-      call register_time_deriv(MIS%v, CS%dv_dt, CS)
+      call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS)
     endif
     if (.not.associated(CS%dh_dt)) then
       call safe_alloc_ptr(CS%dh_dt,isd,ied,jsd,jed,nz)
-      call register_time_deriv(MIS%h, CS%dh_dt, CS)
+      call register_time_deriv(lbound(MIS%h), MIS%h, CS%dh_dt, CS)
     endif
   endif
 

src/parameterizations/vertical/MOM_CVMix_KPP.F90

@@ -147,6 +147,7 @@ module MOM_CVMix_KPP
   real, allocatable, dimension(:,:)   :: OBLdepth_original  !< Depth (positive) of OBL [m] without smoothing
   real, allocatable, dimension(:,:)   :: kOBL      !< Level (+fraction) of OBL extent
   real, allocatable, dimension(:,:)   :: OBLdepthprev !< previous Depth (positive) of OBL [m]
+  real, allocatable, dimension(:,:)   :: La_SL     !< Langmuir number used in KPP
   real, allocatable, dimension(:,:,:) :: dRho      !< Bulk difference in density [kg m-3]
   real, allocatable, dimension(:,:,:) :: Uz2       !< Square of bulk difference in resolved velocity [m2 s-2]
   real, allocatable, dimension(:,:,:) :: BulkRi    !< Bulk Richardson number for each layer (dimensionless)
@@ -536,6 +537,8 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves)
   CS%OBLdepth(:,:) = 0.
   allocate( CS%kOBL( SZI_(G), SZJ_(G) ) )
   CS%kOBL(:,:) = 0.
+  allocate( CS%La_SL( SZI_(G), SZJ_(G) ) )
+  CS%La_SL(:,:) = 0.
   allocate( CS%Vt2( SZI_(G), SZJ_(G), SZK_(G) ) )
   CS%Vt2(:,:,:) = 0.
   if (CS%id_OBLdepth_original > 0) allocate( CS%OBLdepth_original( SZI_(G), SZJ_(G) ) )
@@ -578,7 +581,7 @@ end function KPP_init
 !> KPP vertical diffusivity/viscosity and non-local tracer transport
 subroutine KPP_calculate(CS, G, GV, US, h, uStar, &
                          buoyFlux, Kt, Ks, Kv, nonLocalTransHeat,&
-                         nonLocalTransScalar, Waves)
+                         nonLocalTransScalar, waves)
 
   ! Arguments
   type(KPP_CS),                               pointer       :: CS    !< Control structure
@@ -718,11 +721,11 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, &
            LangEnhK = CS%KPP_K_ENH_FAC
         elseif (CS%LT_K_METHOD==LT_K_MODE_VR12) then
            ! Added minimum value for La_SL, so removed maximum value for LangEnhK.
-           LangEnhK = sqrt(1.+(1.5*WAVES%La_SL(i,j))**(-2) + &
-                (5.4*WAVES%La_SL(i,j))**(-4))
+           LangEnhK = sqrt(1.+(1.5*CS%La_SL(i,j))**(-2) + &
+                (5.4*CS%La_SL(i,j))**(-4))
         elseif (CS%LT_K_METHOD==LT_K_MODE_RW16) then
           !This maximum value is proposed in Reichl et al., 2016 JPO formula
-          LangEnhK = min(2.25, 1. + 1./WAVES%La_SL(i,j))
+          LangEnhK = min(2.25, 1. + 1./CS%La_SL(i,j))
         else
            !This shouldn't be reached.
            !call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in LT_K_ENHANCEMENT")
@@ -1069,15 +1072,10 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyF
       enddo ! k-loop finishes
 
       if (CS%LT_K_ENHANCEMENT .or. CS%LT_VT2_ENHANCEMENT) then
-        if (.not.(present(WAVES).and.associated(WAVES))) then
-          call MOM_error(FATAL,"Trying to use input WAVES information in KPP\n"//&
-               "without activating USEWAVES")
-        endif
-        !For now get Langmuir number based on prev. MLD (otherwise must compute 3d LA)
         MLD_GUESS = max( 1.*US%m_to_Z, abs(US%m_to_Z*CS%OBLdepthprev(i,j) ) )
-        call get_Langmuir_Number( LA, G, GV, US, MLD_guess, uStar(i,j), i, j, &
+        call get_Langmuir_Number( LA, G, GV, US, MLD_guess, surfFricVel, i, j, &
              H=H(i,j,:), U_H=U_H, V_H=V_H, WAVES=WAVES)
-        WAVES%La_SL(i,j)=LA
+        CS%La_SL(i,j)=LA
       endif
 
 
@@ -1125,14 +1123,14 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyF
           enddo
         elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_VR12) then
           !Introduced minimum value for La_SL, so maximum value for enhvt2 is removed.
-          enhvt2 = sqrt(1.+(1.5*WAVES%La_SL(i,j))**(-2) + &
-                   (5.4*WAVES%La_SL(i,j))**(-4))
+          enhvt2 = sqrt(1.+(1.5*CS%La_SL(i,j))**(-2) + &
+                   (5.4*CS%La_SL(i,j))**(-4))
           do k=1,G%ke
              LangEnhVT2(k) = enhvt2
           enddo
         elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_RW16) then
           !Introduced minimum value for La_SL, so maximum value for enhvt2 is removed.
-          enhvt2 = 1. + 2.3*WAVES%La_SL(i,j)**(-0.5)
+          enhvt2 = 1. + 2.3*CS%La_SL(i,j)**(-0.5)
           do k=1,G%ke
             LangEnhVT2(k) = enhvt2
           enddo
@@ -1141,7 +1139,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyF
           do k=1,G%ke
             WST = (max(0.,-buoyflux(i,j,1))*(-cellHeight(k)))**(1./3.)
             LangEnhVT2(k) = sqrt((0.15*WST**3. + 0.17*surfFricVel**3.* &
-                 (1.+0.49*WAVES%La_SL(i,j)**(-2.)))  / &
+                 (1.+0.49*CS%La_SL(i,j)**(-2.)))  / &
                  (0.2*ws_1d(k)**3/(CS%cs*CS%surf_layer_ext*CS%vonKarman**4.)))
           enddo
         else
@@ -1314,11 +1312,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyF
   if (CS%id_BulkUz2  > 0) call post_data(CS%id_BulkUz2,  CS%Uz2,             CS%diag)
   if (CS%id_EnhK     > 0) call post_data(CS%id_EnhK,     CS%EnhK,            CS%diag)
   if (CS%id_EnhVt2   > 0) call post_data(CS%id_EnhVt2,   CS%EnhVt2,          CS%diag)
-  if (present(WAVES)) then
-    if ((CS%id_La_SL>0) .and. associated(WAVES)) then
-      call post_data(CS%id_La_SL,WAVES%La_SL,CS%diag)
-    endif
-  endif
+  if (CS%id_La_SL    > 0) call post_data(CS%id_La_SL,    CS%La_SL,           CS%diag)
 
   ! BLD smoothing:
   if (CS%n_smooth > 0) call KPP_smooth_BLD(CS,G,GV,h)