Merge branch 'IS_converge_crit' into IS_coulomb

NOAA-GFDL · Aug 24, 2023 · e6f30eb · e6f30eb
2 parents ae2acb7 + 9874b2d
commit e6f30eb
Show file tree

Hide file tree

Showing 2 changed files with 71 additions and 21 deletions.
diff --git a/src/parameterizations/vertical/MOM_sponge.F90 b/src/parameterizations/vertical/MOM_sponge.F90
@@ -347,10 +347,14 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml)
                    ! give 0 at the surface [nondim].
 
   real :: e(SZK_(GV)+1)  ! The interface heights [Z ~> m], usually negative.
+  real :: dz_to_h(SZK_(GV)+1)  ! Factors used to convert interface height movement
+                   ! to thickness fluxes [H Z-1 ~> nondim or kg m-3]
   real :: e0       ! The height of the free surface [Z ~> m].
   real :: e_str    ! A nondimensional amount by which the reference
                    ! profile must be stretched for the free surfaces
                    ! heights in the two profiles to agree [nondim].
+  real :: w_mean   ! The vertical displacement of water moving upward through an
+                   ! interface within 1 timestep [Z ~> m].
   real :: w        ! The thickness of water moving upward through an
                    ! interface within 1 timestep [H ~> m or kg m-2].
   real :: wm       ! wm is w if w is negative and 0 otherwise [H ~> m or kg m-2].
@@ -384,9 +388,15 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml)
                   "work properly with i-mean sponges and a bulk mixed layer.")
 
     do j=js,je ; do i=is,ie ; e_D(i,j,nz+1) = -G%bathyT(i,j) ; enddo ; enddo
-    do k=nz,1,-1 ; do j=js,je ; do i=is,ie
-      e_D(i,j,K) = e_D(i,j,K+1) + h(i,j,k)*GV%H_to_Z
-    enddo ; enddo ; enddo
+    if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then
+      do k=nz,1,-1 ; do j=js,je ; do i=is,ie
+        e_D(i,j,K) = e_D(i,j,K+1) + GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k)
+      enddo ; enddo ; enddo
+    else
+      do k=nz,1,-1 ; do j=js,je ; do i=is,ie
+        e_D(i,j,K) = e_D(i,j,K+1) + h(i,j,k)*GV%H_to_Z
+      enddo ; enddo ; enddo
+    endif
     do j=js,je
       do i=is,ie
         dilate(i) = (G%bathyT(i,j) + G%Z_ref) / (e_D(i,j,1) + G%bathyT(i,j))
@@ -424,20 +434,39 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml)
       do K=2,nz+1 ; do i=is,ie
         h_above(i,K) = h_above(i,K-1) + max(h(i,j,k-1)-GV%Angstrom_H, 0.0)
       enddo ; enddo
-      do K=2,nz
-        ! w is positive for an upward (lightward) flux of mass, resulting
-        ! in the downward movement of an interface.
-        w = damp_1pdamp * eta_mean_anom(j,K) * GV%Z_to_H
-        do i=is,ie
+
+      ! In both blocks below, w is positive for an upward (lightward) flux of mass,
+      ! resulting in the downward movement of an interface.
+      if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then
+        do K=2,nz
+          w_mean = damp_1pdamp * eta_mean_anom(j,K)
+          do i=is,ie
+            w = w_mean * 2.0*GV%RZ_to_H / (tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j,k))
+            if (w > 0.0) then
+              w_int(i,j,K) = min(w, h_below(i,K))
+              eb(i,j,k-1) = eb(i,j,k-1) + w_int(i,j,K)
+            else
+              w_int(i,j,K) = max(w, -h_above(i,K))
+              ea(i,j,k) = ea(i,j,k) - w_int(i,j,K)
+            endif
+          enddo
+        enddo
+      else
+        do K=2,nz
+          w = damp_1pdamp * eta_mean_anom(j,K) * GV%Z_to_H
           if (w > 0.0) then
-            w_int(i,j,K) = min(w, h_below(i,K))
-            eb(i,j,k-1) = eb(i,j,k-1) + w_int(i,j,K)
+            do i=is,ie
+              w_int(i,j,K) = min(w, h_below(i,K))
+              eb(i,j,k-1) = eb(i,j,k-1) + w_int(i,j,K)
+            enddo
           else
-            w_int(i,j,K) = max(w, -h_above(i,K))
-            ea(i,j,k) = ea(i,j,k) - w_int(i,j,K)
+            do i=is,ie
+              w_int(i,j,K) = max(w, -h_above(i,K))
+              ea(i,j,k) = ea(i,j,k) - w_int(i,j,K)
+            enddo
           endif
         enddo
-      enddo
+      endif
       do k=1,nz ; do i=is,ie
         ea_k = max(0.0, -w_int(i,j,K))
         eb_k = max(0.0, w_int(i,j,K+1))
@@ -462,9 +491,20 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml)
     damp = dt * CS%Iresttime_col(c)
 
     e(1) = 0.0 ; e0 = 0.0
-    do K=1,nz
-      e(K+1) = e(K) - h(i,j,k)*GV%H_to_Z
-    enddo
+    if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then
+      do K=1,nz
+        e(K+1) = e(K) - GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k)
+      enddo
+      dz_to_h(1) = GV%RZ_to_H / tv%SpV_avg(i,j,1)
+      do K=2,nz
+        dz_to_h(K) = 2.0*GV%RZ_to_H / (tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j,k))
+      enddo
+    else
+      do K=1,nz
+        e(K+1) = e(K) - h(i,j,k)*GV%H_to_Z
+        dz_to_h(K) = GV%Z_to_H
+      enddo
+    endif
     e_str = e(nz+1) / CS%Ref_eta(nz+1,c)
 
     if ( CS%bulkmixedlayer ) then
@@ -481,7 +521,7 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml)
       wpb = 0.0; wb = 0.0
       do k=nz,nkmb+1,-1
         if (GV%Rlay(k) > Rcv_ml(i,j)) then
-          w = MIN((((e(K)-e0) - e_str*CS%Ref_eta(K,c)) * damp)*GV%Z_to_H, &
+          w = MIN((((e(K)-e0) - e_str*CS%Ref_eta(K,c)) * damp)*dz_to_h(K), &
                     ((wb + h(i,j,k)) - GV%Angstrom_H))
           wm = 0.5*(w-ABS(w))
           do m=1,CS%fldno
@@ -537,7 +577,7 @@ subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml)
       wpb = 0.0
       wb = 0.0
       do k=nz,1,-1
-        w = MIN((((e(K)-e0) - e_str*CS%Ref_eta(K,c)) * damp)*GV%Z_to_H, &
+        w = MIN((((e(K)-e0) - e_str*CS%Ref_eta(K,c)) * damp)*dz_to_h(K), &
                   ((wb + h(i,j,k)) - GV%Angstrom_H))
         wm = 0.5*(w - ABS(w))
         do m=1,CS%fldno

diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90
@@ -322,6 +322,7 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg)
   real :: D_edge            ! The thickness [Z ~> m] of the dense fluid at the
                             ! inner edge of the inflow
   real :: RLay_range        ! The range of densities [R ~> kg m-3].
+  real :: Rlay_Ref          ! The surface layer's target density [R ~> kg m-3].
   real :: f_0               ! The reference value of the Coriolis parameter [T-1 ~> s-1]
   real :: f_inflow          ! The value of the Coriolis parameter used to determine DOME inflow
                             ! properties [T-1 ~> s-1]
@@ -351,6 +352,9 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg)
   call get_param(PF, mdl, "DENSITY_RANGE", Rlay_range, &
                  "The range of reference potential densities in the layers.", &
                  units="kg m-3", default=2.0, scale=US%kg_m3_to_R)
+  call get_param(PF, mdl, "LIGHTEST_DENSITY", Rlay_Ref, &
+                 "The reference potential density used for layer 1.", &
+                 units="kg m-3", default=US%R_to_kg_m3*GV%Rho0, scale=US%kg_m3_to_R)
   call get_param(PF, mdl, "F_0", f_0, &
                  "The reference value of the Coriolis parameter with the betaplane option.", &
                  units="s-1", default=0.0, scale=US%T_to_s)
@@ -369,9 +373,15 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg)
 
   if (.not.associated(OBC)) return
 
-  g_prime_tot = (GV%g_Earth / GV%Rho0) * Rlay_range
-  Def_Rad = sqrt(D_edge*g_prime_tot) / abs(f_inflow)
-  tr_0 = (-D_edge*sqrt(D_edge*g_prime_tot)*0.5*Def_Rad) * GV%Z_to_H
+  if (GV%Boussinesq) then
+    g_prime_tot = (GV%g_Earth / GV%Rho0) * Rlay_range
+    Def_Rad = sqrt(D_edge*g_prime_tot) / abs(f_inflow)
+    tr_0 = (-D_edge*sqrt(D_edge*g_prime_tot)*0.5*Def_Rad) * GV%Z_to_H
+  else
+    g_prime_tot = (GV%g_Earth / (Rlay_Ref + 0.5*Rlay_range)) * Rlay_range
+    Def_Rad = sqrt(D_edge*g_prime_tot) / abs(f_inflow)
+    tr_0 = (-D_edge*sqrt(D_edge*g_prime_tot)*0.5*Def_Rad) * (Rlay_Ref + 0.5*Rlay_range) * GV%RZ_to_H
+  endif
 
   I_Def_Rad = 1.0 / (1.0e-3*US%L_to_m*Def_Rad)