Merge pull request #63 from NOAA-GFDL/dev/gfdl

Sync with NOAA-GFDL dev/gfdl

.testing/Makefile

@@ -92,11 +92,34 @@ TARGET_SOURCE = $(call SOURCE,build/target_codebase/src) \
   $(wildcard build/target_codebase/config_src/ext*/*.F90)
 FMS_SOURCE = $(call SOURCE,$(DEPS)/fms/src)
 
+#---
+# Python preprocessing environment configuration
+
+HAS_NUMPY = $(shell python -c "import numpy" 2> /dev/null && echo "yes")
+HAS_NETCDF4 = $(shell python -c "import netCDF4" 2> /dev/null && echo "yes")
+
+USE_VENV =
+ifneq ($(HAS_NUMPY), yes)
+  USE_VENV = yes
+endif
+ifneq ($(HAS_NETCDF4), yes)
+  USE_VENV = yes
+endif
+
+# When disabled, activation is a null operation (`true`)
+VENV_PATH =
+VENV_ACTIVATE = true
+ifeq ($(USE_VENV), yes)
+  VENV_PATH = work/local-env
+  VENV_ACTIVATE = . $(VENV_PATH)/bin/activate
+endif
+
+
 #---
 # Rules
 
 .PHONY: all build.regressions
-all: $(foreach b,$(BUILDS),build/$(b)/MOM6)
+all: $(foreach b,$(BUILDS),build/$(b)/MOM6) $(VENV_PATH)
 build.regressions: $(foreach b,symmetric target,build/$(b)/MOM6)
 
 # Executable
@@ -184,6 +207,18 @@ $(LIST_PATHS) $(MKMF):
 	cd $(DEPS)/mkmf; git checkout $(MKMF_COMMIT)
 
 
+#---
+# Python preprocessing
+# NOTE: Some less mature environments (e.g. Arm64 Ubuntu) require explicit
+#   installation of numpy before netCDF4, as well as wheel and cython support.
+work/local-env:
+	python3 -m venv $@
+	. $@/bin/activate \
+	  && pip3 install wheel \
+	  && pip3 install cython \
+	  && pip3 install numpy \
+	  && pip3 install netCDF4
+
 #----
 # Testing
 
@@ -264,7 +299,6 @@ $(eval $(call CMP_RULE,regression,symmetric target))
 
 # TODO: chksum_diag parsing of restart files
 
-
 #---
 # Test run output files
 
@@ -281,7 +315,13 @@ work/%/$(1)/ocean.stats work/%/$(1)/chksum_diag: build/$(2)/MOM6
 	if [ $(3) ]; then find build/$(2) -name *.gcda -exec rm -f '{}' \; ; fi
 	mkdir -p $$(@D)
 	cp -rL $$*/* $$(@D)
-	cd $$(@D) && if [ -f Makefile ]; then $(MAKE); fi
+	if [ -f $$(@D)/Makefile ]; then \
+	  $$(VENV_ACTIVATE) \
+	    && cd $$(@D) \
+	    && $(MAKE); \
+	else \
+	  cd $$(@D); \
+	fi
 	mkdir -p $$(@D)/RESTART
 	echo -e "$(4)" > $$(@D)/MOM_override
 	cd $$(@D) \
@@ -327,7 +367,13 @@ work/%/restart/ocean.stats: build/symmetric/MOM6
 	rm -rf $(@D)
 	mkdir -p $(@D)
 	cp -rL $*/* $(@D)
-	cd work/$*/restart && if [ -f Makefile ]; then $(MAKE); fi
+	if [ -f $(@D)/Makefile ]; then \
+	  $(VENV_ACTIVATE) \
+	    && cd work/$*/restart \
+	    && $(MAKE); \
+	else \
+	  cd work/$*/restart; \
+	fi
 	mkdir -p $(@D)/RESTART
 	# Generate the half-period input namelist
 	# TODO: Assumes that runtime set by DAYMAX, will fail if set by input.nml

.testing/tc4/.gitignore

@@ -0,0 +1,4 @@
+ocean_hgrid.nc
+sponge.nc
+temp_salt_ic.nc
+topog.nc

.testing/tc4/Makefile

@@ -1,3 +1,8 @@
-ocean_hgrid.nc topog.nc temp_salt_ic.nc sponge.nc:
+OUT=ocean_hgrid.nc sponge.nc temp_salt_ic.nc topog.nc
+
+$(OUT):
 	python build_grid.py
 	python build_data.py
+
+clean:
+	rm -rf $(OUT)

.travis.yml

@@ -14,6 +14,7 @@ addons:
     - mpich libmpich-dev
     - doxygen graphviz flex bison cmake
     - python-numpy python-netcdf4
+    - python3 python3-dev python3-venv python3-pip
     - bc
 
 jobs:

src/core/MOM_CoriolisAdv.F90

@@ -256,7 +256,7 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
   enddo ; enddo
 
   !$OMP parallel do default(private) shared(u,v,h,uh,vh,CAu,CAv,G,CS,AD,Area_h,Area_q,&
-  !$OMP                        RV,PV,is,ie,js,je,Isq,Ieq,Jsq,Jeq,nz,h_neglect,h_tiny,OBC)
+  !$OMP                        RV,PV,is,ie,js,je,Isq,Ieq,Jsq,Jeq,nz,h_neglect,h_tiny,OBC,eps_vel)
   do k=1,nz
 
     ! Here the second order accurate layer potential vorticities, q,

src/core/MOM_PressureForce_Montgomery.F90

@@ -507,7 +507,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
     ! This no longer includes any pressure dependency, since this routine
     ! will come down with a fatal error if there is any compressibility.
     !$OMP parallel do default(shared)
-    do k=1,nz+1 ; do j=Jsq,Jeq+1
+    do k=1,nz ; do j=Jsq,Jeq+1
       call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), &
                              tv%eqn_of_state, EOSdom)
       do i=Isq,Ieq+1 ; rho_star(i,j,k) = G_Rho0*rho_star(i,j,k) ; enddo

src/core/MOM_open_boundary.F90

@@ -1590,7 +1590,7 @@ subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CSp)
   ! Local variables
   real :: vel2_rescale ! A rescaling factor for squared velocities from the representation in
                        ! a restart file to the internal representation in this run.
-  integer :: i, j, k, isd, ied, jsd, jed, nz
+  integer :: i, j, k, isd, ied, jsd, jed, nz, m
   integer :: IsdB, IedB, JsdB, JedB
   isd  = G%isd  ; ied  = G%ied  ; jsd  = G%jsd  ; jed  = G%jed ; nz = GV%ke
   IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
@@ -1603,6 +1603,11 @@ subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CSp)
   if (OBC%oblique_BCs_exist_globally) call pass_vector(OBC%rx_oblique, OBC%ry_oblique, G%Domain, &
                      To_All+Scalar_Pair)
   if (associated(OBC%cff_normal)) call pass_var(OBC%cff_normal, G%Domain, position=CORNER)
+  if (associated(OBC%tres_x) .or. associated(OBC%tres_y)) then
+    do m=1,OBC%ntr
+      call pass_vector(OBC%tres_x(:,:,:,m), OBC%tres_y(:,:,:,m), G%Domain, To_All+Scalar_Pair)
+    enddo
+  endif
 
   ! The rx_normal and ry_normal arrays used with radiation OBCs are currently in units of grid
   ! points per timestep, but if this were to be corrected to [L T-1 ~> m s-1] or [T-1 ~> s-1] to
@@ -4711,7 +4716,8 @@ subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart
   endif
 
   ! Still painfully inefficient, now in four dimensions.
-  if (any(OBC%tracer_x_reservoirs_used)) then
+  ! Allocating both for now so that the pass_vector works.
+  if (any(OBC%tracer_x_reservoirs_used) .or. any(OBC%tracer_y_reservoirs_used)) then
     allocate(OBC%tres_x(HI%isdB:HI%iedB,HI%jsd:HI%jed,GV%ke,OBC%ntr))
     OBC%tres_x(:,:,:,:) = 0.0
     do m=1,OBC%ntr
@@ -4727,8 +4733,8 @@ subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart
         endif
       endif
     enddo
-  endif
-  if (any(OBC%tracer_y_reservoirs_used)) then
+! endif
+! if (any(OBC%tracer_y_reservoirs_used)) then
     allocate(OBC%tres_y(HI%isd:HI%ied,HI%jsdB:HI%jedB,GV%ke,OBC%ntr))
     OBC%tres_y(:,:,:,:) = 0.0
     do m=1,OBC%ntr

src/diagnostics/MOM_wave_speed.F90

@@ -806,7 +806,8 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos, better_spee
 
   min_h_frac = tol_Hfrac / real(nz)
   !$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,GV,US,min_h_frac,use_EOS,T,S, &
-  !$OMP                                     Z_to_pres,tv,cn,g_Rho0,nmodes)
+  !$OMP                                     Z_to_pres,tv,cn,g_Rho0,nmodes,cg1_min2,better_est, &
+  !$OMP                                     c1_thresh,tol_solve,tol_merge)
   do j=js,je
     !   First merge very thin layers with the one above (or below if they are
     ! at the top).  This also transposes the row order so that columns can

src/ice_shelf/MOM_ice_shelf_dynamics.F90

@@ -1405,6 +1405,7 @@ subroutine ice_shelf_advect_thickness_y(CS, G, LB, time_step, hmask, h0, h_after
   real :: h_face     ! Thickness at a face for transport [Z ~> m]
   real :: slope_lim  ! The value of the slope limiter, in the range of 0 to 2 [nondim]
 
+  ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh
 
   ! hmask coded values: 1) fully covered; 2) partly covered - no export; 3) Specified boundary condition
   ! relevant u_face_mask coded values: 1) Normal interior point; 4) Specified flux BC

src/parameterizations/vertical/MOM_kappa_shear.F90

@@ -396,8 +396,8 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_
   real :: dz_in_lay     !   The running sum of the thickness in a layer [Z ~> m].
   real :: k0dt          ! The background diffusivity times the timestep [Z2 ~> m2].
   real :: dz_massless   ! A layer thickness that is considered massless [Z ~> m].
-  real :: I_hwt ! The inverse of the masked thickness weights [H-1 ~> m-1 or m2 kg-1].
-  real :: I_Prandtl
+  real :: I_hwt         ! The inverse of the masked thickness weights [H-1 ~> m-1 or m2 kg-1].
+  real :: I_Prandtl     ! The inverse of the turbulent Prandtl number [nondim].
   logical :: use_temperature  !  If true, temperature and salinity have been
                         ! allocated and are being used as state variables.
   logical :: new_kappa = .true. ! If true, ignore the value of kappa from the
@@ -422,7 +422,7 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_
   I_Prandtl = 0.0 ; if (CS%Prandtl_turb > 0.0) I_Prandtl = 1.0 / CS%Prandtl_turb
 
   !$OMP parallel do default(private) shared(jsB,jeB,isB,ieB,nz,h,u_in,v_in,use_temperature,new_kappa, &
-  !$OMP                                tv,G,GV,US,CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io)
+  !$OMP                                tv,G,GV,US,CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io,I_Prandtl)
   do J=JsB,JeB
     J2 = mod(J,2)+1 ; J2m1 = 3-J2 ! = mod(J-1,2)+1
 

src/parameterizations/vertical/MOM_vert_friction.F90

@@ -629,6 +629,8 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
   real :: z2      ! The distance from the bottom, normalized by Hbbl, nondim.
   real :: z2_wt   ! A nondimensional (0-1) weight used when calculating z2.
   real :: z_clear ! The clearance of an interface above the surrounding topography [H ~> m or kg m-2].
+  real :: a_cpl_max  ! The maximum drag doefficient across interfaces, set so that it will be
+                     ! representable as a 32-bit float in MKS units  [Z T-1 ~> m s-1]
   real :: h_neglect  ! A thickness that is so small it is usually lost
                      ! in roundoff and can be neglected [H ~> m or kg m-2].
 
@@ -648,6 +650,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
          "Module must be initialized before it is used.")
 
   h_neglect = GV%H_subroundoff
+  a_cpl_max = 1.0e37 * US%m_to_Z * US%T_to_s
   I_Hbbl(:) = 1.0 / (CS%Hbbl + h_neglect)
   I_valBL = 0.0 ; if (CS%harm_BL_val > 0.0) I_valBL = 1.0 / CS%harm_BL_val
 
@@ -676,7 +679,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
   endif
 
   !$OMP parallel do default(private) shared(G,GV,US,CS,visc,Isq,Ieq,nz,u,h,forces,hML_u, &
-  !$OMP                                     OBC,h_neglect,dt,I_valBL,Kv_u) &
+  !$OMP                                     OBC,h_neglect,dt,I_valBL,Kv_u,a_cpl_max) &
   !$OMP                     firstprivate(i_hbbl)
   do j=G%Jsc,G%Jec
     do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0) ; enddo
@@ -811,13 +814,13 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
 
     if (do_any_shelf) then
       do K=1,nz+1 ; do I=Isq,Ieq ; if (do_i_shelf(I)) then
-        CS%a_u(I,j,K) = forces%frac_shelf_u(I,j)  * a_shelf(I,K) + &
-                   (1.0-forces%frac_shelf_u(I,j)) * a_cpl(I,K)
+        CS%a_u(I,j,K) = min(a_cpl_max, forces%frac_shelf_u(I,j)  * a_shelf(I,K) + &
+                                       (1.0-forces%frac_shelf_u(I,j)) * a_cpl(I,K))
 ! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH
-!        CS%a_u(I,j,K) = forces%frac_shelf_u(I,j)  * max(a_shelf(I,K), a_cpl(I,K)) + &
-!                   (1.0-forces%frac_shelf_u(I,j)) * a_cpl(I,K)
+!        CS%a_u(I,j,K) = min(a_cpl_max, forces%frac_shelf_u(I,j)  * max(a_shelf(I,K), a_cpl(I,K)) + &
+!                                       (1.0-forces%frac_shelf_u(I,j)) * a_cpl(I,K))
       elseif (do_i(I)) then
-        CS%a_u(I,j,K) = a_cpl(I,K)
+        CS%a_u(I,j,K) = min(a_cpl_max, a_cpl(I,K))
       endif ; enddo ; enddo
       do k=1,nz ; do I=Isq,Ieq ; if (do_i_shelf(I)) then
         ! Should we instead take the inverse of the average of the inverses?
@@ -827,7 +830,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
         CS%h_u(I,j,k) = hvel(I,k) + h_neglect
       endif ; enddo ; enddo
     else
-      do K=1,nz+1 ; do I=Isq,Ieq ; if (do_i(I)) CS%a_u(I,j,K) = a_cpl(I,K) ; enddo ; enddo
+      do K=1,nz+1 ; do I=Isq,Ieq ; if (do_i(I)) CS%a_u(I,j,K) = min(a_cpl_max, a_cpl(I,K)) ; enddo ; enddo
       do k=1,nz ; do I=Isq,Ieq ; if (do_i(I)) CS%h_u(I,j,k) = hvel(I,k) + h_neglect ; enddo ; enddo
     endif
 
@@ -843,7 +846,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
 
   ! Now work on v-points.
   !$OMP parallel do default(private) shared(G,GV,CS,US,visc,is,ie,Jsq,Jeq,nz,v,h,forces,hML_v, &
-  !$OMP                                  OBC,h_neglect,dt,I_valBL,Kv_v) &
+  !$OMP                                  OBC,h_neglect,dt,I_valBL,Kv_v,a_cpl_max) &
   !$OMP                     firstprivate(i_hbbl)
   do J=Jsq,Jeq
     do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0) ; enddo
@@ -979,13 +982,13 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
 
     if (do_any_shelf) then
       do K=1,nz+1 ; do i=is,ie ; if (do_i_shelf(i)) then
-        CS%a_v(i,J,K) = forces%frac_shelf_v(i,J)  * a_shelf(i,k) + &
-                   (1.0-forces%frac_shelf_v(i,J)) * a_cpl(i,K)
+        CS%a_v(i,J,K) = min(a_cpl_max, forces%frac_shelf_v(i,J)  * a_shelf(i,k) + &
+                                       (1.0-forces%frac_shelf_v(i,J)) * a_cpl(i,K))
 ! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH
-!        CS%a_v(i,J,K) = forces%frac_shelf_v(i,J)  * max(a_shelf(i,K), a_cpl(i,K)) + &
-!                   (1.0-forces%frac_shelf_v(i,J)) * a_cpl(i,K)
+!        CS%a_v(i,J,K) = min(a_cpl_max, forces%frac_shelf_v(i,J)  * max(a_shelf(i,K), a_cpl(i,K)) + &
+                    !                   (1.0-forces%frac_shelf_v(i,J)) * a_cpl(i,K))
       elseif (do_i(i)) then
-        CS%a_v(i,J,K) = a_cpl(i,K)
+        CS%a_v(i,J,K) = min(a_cpl_max, a_cpl(i,K))
       endif ; enddo ; enddo
       do k=1,nz ; do i=is,ie ; if (do_i_shelf(i)) then
         ! Should we instead take the inverse of the average of the inverses?
@@ -995,7 +998,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
         CS%h_v(i,J,k) = hvel(i,k) + h_neglect
       endif ; enddo ; enddo
     else
-      do K=1,nz+1 ; do i=is,ie ; if (do_i(i)) CS%a_v(i,J,K) = a_cpl(i,K) ; enddo ; enddo
+      do K=1,nz+1 ; do i=is,ie ; if (do_i(i)) CS%a_v(i,J,K) = min(a_cpl_max, a_cpl(i,K)) ; enddo ; enddo
       do k=1,nz ; do i=is,ie ; if (do_i(i)) CS%h_v(i,J,k) = hvel(i,k) + h_neglect ; enddo ; enddo
     endif
 
@@ -1109,10 +1112,10 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i,
   nz = G%ke
   h_neglect = GV%H_subroundoff
 
-  !   The maximum coupling coefficent was originally introduced to avoid
-  ! truncation error problems in the tridiagonal solver. Effectively, the 1e-10
-  ! sets the maximum coupling coefficient increment to 1e10 m per timestep.
   if (CS%answers_2018) then
+    !   The maximum coupling coefficent was originally introduced to avoid
+    ! truncation error problems in the tridiagonal solver. Effectively, the 1e-10
+    ! sets the maximum coupling coefficient increment to 1e10 m per timestep.
     I_amax = (1.0e-10*US%Z_to_m) * dt
   else
     I_amax = 0.0