Fixes for a 2-d simulation

kshedstrom · Nov 18, 2019 · 92785e5 · 92785e5
1 parent b98a6aa
commit 92785e5
Show file tree

Hide file tree

Showing 7 changed files with 99 additions and 52 deletions.
diff --git a/ROMS/Modules/mod_forces.F b/ROMS/Modules/mod_forces.F
@@ -871,6 +871,42 @@ SUBROUTINE allocate_forces (ng, LBi, UBi, LBj, UBj)
       Dmem(ng)=Dmem(ng)+size2d
 #endif
 
+#if defined BULK_FLUXES || defined BULK_FLUXES2D || defined ECOSIM \
+    || defined CCSM_FLUXES2D
+      allocate ( FORCES(ng) % Uwind(LBi:UBi,LBj:UBj) )
+      Dmem(ng)=Dmem(ng)+size2d
+
+      allocate ( FORCES(ng) % Vwind(LBi:UBi,LBj:UBj) )
+      Dmem(ng)=Dmem(ng)+size2d
+
+# ifndef ANA_WINDS
+      allocate ( FORCES(ng) % UwindG(LBi:UBi,LBj:UBj,2) )
+      Dmem(ng)=Dmem(ng)+2.0_r8*size2d
+
+      allocate ( FORCES(ng) % VwindG(LBi:UBi,LBj:UBj,2) )
+      Dmem(ng)=Dmem(ng)+2.0_r8*size2d
+# endif
+#endif
+
+#if defined BULK_FLUXES || defined ECOSIM || defined CCSM_FLUXES2D ||  \
+    (defined SHORTWAVE && defined ANA_SRFLUX)
+      allocate ( FORCES(ng) % Hair(LBi:UBi,LBj:UBj) )
+      Dmem(ng)=Dmem(ng)+size2d
+
+# ifndef ANA_HUMIDITY
+      allocate ( FORCES(ng) % HairG(LBi:UBi,LBj:UBj,2) )
+      Dmem(ng)=Dmem(ng)+2.0_r8*size2d
+# endif
+
+      allocate ( FORCES(ng) % Tair(LBi:UBi,LBj:UBj) )
+      Dmem(ng)=Dmem(ng)+size2d
+
+# ifndef ANA_TAIR
+      allocate ( FORCES(ng) % TairG(LBi:UBi,LBj:UBj,2) )
+      Dmem(ng)=Dmem(ng)+2.0_r8*size2d
+# endif
+#endif
+
 #ifdef SOLVE3D
 
 # ifdef SHORTWAVE
@@ -966,41 +1002,6 @@ SUBROUTINE allocate_forces (ng, LBi, UBi, LBj, UBj)
       Dmem(ng)=Dmem(ng)+size2d
 # endif
 
-#if defined BULK_FLUXES || defined ECOSIM || defined CCSM_FLUXES2D ||  \
-    (defined SHORTWAVE && defined ANA_SRFLUX)
-      allocate ( FORCES(ng) % Hair(LBi:UBi,LBj:UBj) )
-      Dmem(ng)=Dmem(ng)+size2d
-
-#  ifndef ANA_HUMIDITY
-      allocate ( FORCES(ng) % HairG(LBi:UBi,LBj:UBj,2) )
-      Dmem(ng)=Dmem(ng)+2.0_r8*size2d
-#  endif
-
-      allocate ( FORCES(ng) % Tair(LBi:UBi,LBj:UBj) )
-      Dmem(ng)=Dmem(ng)+size2d
-
-#  ifndef ANA_TAIR
-      allocate ( FORCES(ng) % TairG(LBi:UBi,LBj:UBj,2) )
-      Dmem(ng)=Dmem(ng)+2.0_r8*size2d
-#  endif
-# endif
-
-# if defined BULK_FLUXES || defined BULK_FLUXES2D || defined ECOSIM
-      allocate ( FORCES(ng) % Uwind(LBi:UBi,LBj:UBj) )
-      Dmem(ng)=Dmem(ng)+size2d
-
-      allocate ( FORCES(ng) % Vwind(LBi:UBi,LBj:UBj) )
-      Dmem(ng)=Dmem(ng)+size2d
-
-#  ifndef ANA_WINDS
-      allocate ( FORCES(ng) % UwindG(LBi:UBi,LBj:UBj,2) )
-      Dmem(ng)=Dmem(ng)+2.0_r8*size2d
-
-      allocate ( FORCES(ng) % VwindG(LBi:UBi,LBj:UBj,2) )
-      Dmem(ng)=Dmem(ng)+2.0_r8*size2d
-#  endif
-# endif
-
 # ifdef BULK_FLUXES
       allocate ( FORCES(ng) % rain(LBi:UBi,LBj:UBj) )
       Dmem(ng)=Dmem(ng)+size2d
@@ -1523,7 +1524,8 @@ SUBROUTINE initialize_forces (ng, tile, model)
 # endif
 #endif
 
-#if defined BULK_FLUXES || defined BULK_FLUXES2D || defined ECOSIM
+#if defined BULK_FLUXES || defined BULK_FLUXES2D || defined ECOSIM \
+    || defined CCSM_FLUXES2D
             FORCES(ng) % Uwind(i,j) = IniVal
             FORCES(ng) % Vwind(i,j) = IniVal
 # ifndef ANA_WINDS

diff --git a/ROMS/Modules/mod_grid.F b/ROMS/Modules/mod_grid.F
@@ -244,7 +244,7 @@ MODULE mod_grid
           real(r8), pointer :: rmask(:,:)
           real(r8), pointer :: umask(:,:)
           real(r8), pointer :: vmask(:,:)
-# ifdef BERING_STRAIT
+# ifdef MASK_HACK
           real(r8), pointer :: mask2(:,:)
 # endif
 # ifdef ALBEDO_HACK
@@ -589,7 +589,7 @@ SUBROUTINE allocate_grid (ng, LBi, UBi, LBj, UBj, LBij, UBij)
       allocate ( GRID(ng) % vmask(LBi:UBi,LBj:UBj) )
       Dmem(ng)=Dmem(ng)+size2d
 
-#ifdef BERING_STRAIT
+#ifdef MASK_HACK
       allocate ( GRID(ng) % mask2(LBi:UBi,LBj:UBj) )
       Dmem(ng)=Dmem(ng)+size2d
 #endif
@@ -920,7 +920,7 @@ SUBROUTINE initialize_grid (ng, tile, model)
             GRID(ng) % rmask(i,j) = IniMetricVal
             GRID(ng) % umask(i,j) = IniMetricVal
             GRID(ng) % vmask(i,j) = IniMetricVal
-#ifdef BERING_STRAIT
+#ifdef MASK_HACK
             GRID(ng) % mask2(i,j) = IniVal
 #endif
 #ifdef ALBEDO_HACK

diff --git a/ROMS/Modules/mod_scalars.F b/ROMS/Modules/mod_scalars.F
@@ -3189,7 +3189,7 @@ SUBROUTINE allocate_scalars
         Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
       END IF
 
-#ifdef BULK_FLUXES
+#if defined BULK_FLUXES || defined BULK_FLUXES2D
       IF (.not.allocated(blk_ZQ)) THEN
         allocate ( blk_ZQ(Ngrids) )
         Dmem(1)=Dmem(1)+REAL(Ngrids,r8)

diff --git a/ROMS/Nonlinear/bulk_flux.F b/ROMS/Nonlinear/bulk_flux.F
@@ -1964,7 +1964,7 @@ SUBROUTINE bulk_flux (ng, tile)
      &                     GRID(ng) % rmask,                            &
      &                     GRID(ng) % umask,                            &
      &                     GRID(ng) % vmask,                            &
-     &                     GRID(ng) % mask2,                            &
+!    &                     GRID(ng) % mask2,                            &
 # endif
      &                     FORCES(ng) % Uwind,                          &
      &                     FORCES(ng) % Vwind,                          &
@@ -1987,7 +1987,7 @@ SUBROUTINE bulk_flux_tile (ng, tile,                              &
 # endif
 # ifdef MASKING
      &                           rmask, umask, vmask,                   &
-     &                           mask2,                                 &
+!    &                           mask2,                                 &
 # endif
      &                           Uwind, Vwind, sustr, svstr)
 !***********************************************************************
@@ -2016,21 +2016,21 @@ SUBROUTINE bulk_flux_tile (ng, tile,                              &
       real(r8), intent(in) :: rmask(LBi:,LBj:)
       real(r8), intent(in) :: umask(LBi:,LBj:)
       real(r8), intent(in) :: vmask(LBi:,LBj:)
-      real(r8), intent(in) :: mask2(LBi:,LBj:)
+!     real(r8), intent(in) :: mask2(LBi:,LBj:)
 #  endif
-      real(r8), intent(out) :: Uwind(LBi:,LBj:)
-      real(r8), intent(out) :: Vwind(LBi:,LBj:)
+      real(r8), intent(in) :: Uwind(LBi:,LBj:)
+      real(r8), intent(in) :: Vwind(LBi:,LBj:)
       real(r8), intent(out) :: sustr(LBi:,LBj:)
       real(r8), intent(out) :: svstr(LBi:,LBj:)
 # else
 #  ifdef MASKING
       real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
       real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
       real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
-      real(r8), intent(in) :: mask2(LBi:UBi,LBj:UBj)
+!     real(r8), intent(in) :: mask2(LBi:UBi,LBj:UBj)
 #  endif
-      real(r8), intent(out) :: Uwind(LBi:,LBj:)
-      real(r8), intent(out) :: Vwind(LBi:,LBj:)
+      real(r8), intent(in) :: Uwind(LBi:,LBj:)
+      real(r8), intent(in) :: Vwind(LBi:,LBj:)
       real(r8), intent(out) :: sustr(LBi:UBi,LBj:UBj)
       real(r8), intent(out) :: svstr(LBi:UBi,LBj:UBj)
 # endif

diff --git a/ROMS/Nonlinear/ccsm_flux.F b/ROMS/Nonlinear/ccsm_flux.F
@@ -1833,6 +1833,8 @@ SUBROUTINE ccsm_flux_tile (ng, tile,                              &
           ssq    = 0.98_r8 * qsat(TairK(i)) / rhoAir(i)  ! sea surf hum (kg/kg)
 !          ssq    = 0.98_r8 * qsat(TseaK) / rhoAir(i)    ! sea surf hum (kg/kg)
           delq   = Hair(i,j) - ssq                     ! spec hum dif (kg/kg)
+          ! Note: this whole option is untested and unimplemented due to this.
+          delt = ???
           alz    = log(blk_ZW(ng)/zref)
           cpair     = shr_const_cpdair*(1.0 + shr_const_cpvir*ssq)
 

diff --git a/ROMS/Nonlinear/set_vbc.F b/ROMS/Nonlinear/set_vbc.F
@@ -1148,13 +1148,15 @@ SUBROUTINE set_vbc (ng, tile)
      &                   LBi, UBi, LBj, UBj,                            &
      &                   IminS, ImaxS, JminS, JmaxS,                    &
      &                   krhs(ng), kstp(ng), knew(ng),                  &
+     &                   GRID(ng) % h,                                  &
 #  if defined UV_LDRAG
      &                   GRID(ng) % rdrag,                              &
 #  elif defined UV_QDRAG
      &                   GRID(ng) % rdrag2,                             &
 #  endif
      &                   OCEAN(ng) % ubar,                              &
      &                   OCEAN(ng) % vbar,                              &
+     &                   OCEAN(ng) % zeta,                              &
      &                   FORCES(ng) % bustr,                            &
      &                   FORCES(ng) % bvstr)
 #  ifdef PROFILE
@@ -1169,12 +1171,13 @@ SUBROUTINE set_vbc_tile (ng, tile,                                &
      &                         LBi, UBi, LBj, UBj,                      &
      &                         IminS, ImaxS, JminS, JmaxS,              &
      &                         krhs, kstp, knew,                        &
+     &                         h,                                       &
 #  if defined UV_LDRAG
      &                         rdrag,                                   &
 #  elif defined UV_QDRAG
      &                         rdrag2,                                  &
 #  endif
-     &                         ubar, vbar, bustr, bvstr)
+     &                         ubar, vbar, zeta, bustr, bvstr)
 !***********************************************************************
 !
       USE mod_param
@@ -1193,6 +1196,7 @@ SUBROUTINE set_vbc_tile (ng, tile,                                &
       integer, intent(in) :: krhs, kstp, knew
 !
 #  ifdef ASSUMED_SHAPE
+      real(r8), intent(in) :: h(LBi:,LBj:)
 #   ifdef UV_LDRAG
       real(r8), intent(in) :: rdrag(LBi:,LBj:)
 #   endif
@@ -1201,9 +1205,11 @@ SUBROUTINE set_vbc_tile (ng, tile,                                &
 #   endif
       real(r8), intent(in) :: ubar(LBi:,LBj:,:)
       real(r8), intent(in) :: vbar(LBi:,LBj:,:)
+      real(r8), intent(in) :: zeta(LBi:,LBj:,:)
       real(r8), intent(inout) :: bustr(LBi:,LBj:)
       real(r8), intent(inout) :: bvstr(LBi:,LBj:)
 #  else
+      real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
 #   ifdef UV_LDRAG
       real(r8), intent(in) :: rdrag(LBi:UBi,LBj:UBj)
 #   endif
@@ -1212,6 +1218,7 @@ SUBROUTINE set_vbc_tile (ng, tile,                                &
 #   endif
       real(r8), intent(in) :: ubar(LBi:UBi,LBj:UBj,3)
       real(r8), intent(in) :: vbar(LBi:UBi,LBj:UBj,3)
+      real(r8), intent(in) :: zeta(LBi:UBi,LBj:UBj,3)
       real(r8), intent(inout) :: bustr(LBi:UBi,LBj:UBj)
       real(r8), intent(inout) :: bvstr(LBi:UBi,LBj:UBj)
 #  endif
@@ -1220,14 +1227,22 @@ SUBROUTINE set_vbc_tile (ng, tile,                                &
 !
       integer :: i, j
 
-      real(r8) :: cff1, cff2
+      real(r8) :: cff, cff1, cff2, cff3
 
 #  include "set_bounds.h"
 !
 !-----------------------------------------------------------------------
 !  Set kinematic barotropic bottom momentum stress (m2/s2).
 !-----------------------------------------------------------------------
 
+#  ifdef LIMIT_BSTRESS
+!
+!  Set limiting factor for bottom stress. The bottom stress is adjusted
+!  to not change the direction of momentum.  It only should slow down
+!  to zero.  The value of 0.75 is arbitrary limitation assigment.
+!
+      cff=0.3_r8/dt(ng)
+#  endif
 #  if defined UV_LDRAG
 !
 !  Set linear bottom stress.
@@ -1236,12 +1251,26 @@ SUBROUTINE set_vbc_tile (ng, tile,                                &
         DO i=IstrU,Iend
           bustr(i,j)=0.5_r8*(rdrag(i-1,j)+rdrag(i,j))*                  &
      &               ubar(i,j,krhs)
+#    ifdef LIMIT_BSTRESS
+          cff3=cff*0.5_r8*(h(i-1,j)+h(i,j)+                             &
+     &             zeta(i-1,j,krhs)+zeta(i,j,krhs))
+          bustr(i,j)=SIGN(1.0_r8, bustr(i,j))*                          &
+     &               MIN(ABS(bustr(i,j)),                               &
+     &                   ABS(ubar(i,j,krhs))*cff3)
+#    endif
         END DO
       END DO
       DO j=JstrV,Jend
         DO i=Istr,Iend
           bvstr(i,j)=0.5_r8*(rdrag(i,j-1)+rdrag(i,j))*                  &
      &               vbar(i,j,krhs)
+#    ifdef LIMIT_BSTRESS
+          cff3=cff*0.5_r8*(h(i,j-1)+h(i,j)+                             &
+     &             zeta(i,j-1,krhs)+zeta(i,j,krhs))
+          bvstr(i,j)=SIGN(1.0_r8, bvstr(i,j))*                          &
+     &               MIN(ABS(bvstr(i,j)),                               &
+     &                   ABS(vbar(i,j,krhs))*cff3)
+#    endif
         END DO
       END DO
 #  elif defined UV_QDRAG
@@ -1257,6 +1286,13 @@ SUBROUTINE set_vbc_tile (ng, tile,                                &
           cff2=SQRT(ubar(i,j,krhs)*ubar(i,j,krhs)+cff1*cff1)
           bustr(i,j)=0.5_r8*(rdrag2(i-1,j)+rdrag2(i,j))*                &
      &               ubar(i,j,krhs)*cff2
+#    ifdef LIMIT_BSTRESS
+          cff3=cff*0.5_r8*(h(i-1,j)+h(i,j)+                             &
+     &             zeta(i-1,j,krhs)+zeta(i,j,krhs))
+          bustr(i,j)=SIGN(1.0_r8, bustr(i,j))*                          &
+     &               MIN(ABS(bustr(i,j)),                               &
+     &                   ABS(ubar(i,j,krhs))*cff3)
+#    endif
         END DO
       END DO
       DO j=JstrV,Jend
@@ -1268,6 +1304,13 @@ SUBROUTINE set_vbc_tile (ng, tile,                                &
           cff2=SQRT(cff1*cff1+vbar(i,j,krhs)*vbar(i,j,krhs))
           bvstr(i,j)=0.5_r8*(rdrag2(i,j-1)+rdrag2(i,j))*                &
      &               vbar(i,j,krhs)*cff2
+#    ifdef LIMIT_BSTRESS
+          cff3=cff*0.5_r8*(h(i,j-1)+h(i,j)+                             &
+     &             zeta(i,j-1,krhs)+zeta(i,j,krhs))
+          bvstr(i,j)=SIGN(1.0_r8, bvstr(i,j))*                          &
+     &               MIN(ABS(bvstr(i,j)),                               &
+     &                   ABS(vbar(i,j,krhs))*cff3)
+#    endif
         END DO
       END DO
 #  endif

diff --git a/ROMS/Utility/read_phypar.F b/ROMS/Utility/read_phypar.F
@@ -867,7 +867,7 @@ SUBROUTINE read_PhyPar (model, inp, out, Lwrite)
               Npts=load_r(Nval, Rval, Ngrids, Zob)
             CASE ('Zos')
               Npts=load_r(Nval, Rval, Ngrids, Zos)
-#ifdef BULK_FLUXES
+#if defined BULK_FLUXES || defined BULK_FLUXES2D
             CASE ('BLK_ZQ')
               Npts=load_r(Nval, Rval, Ngrids, blk_ZQ)
             CASE ('BLK_ZT')
@@ -4404,7 +4404,7 @@ SUBROUTINE read_PhyPar (model, inp, out, Lwrite)
           WRITE (out,200) Zos(ng), 'Zos',                               &
      &          'Surface roughness (m).'
 # endif
-# ifdef BULK_FLUXES
+# if defined BULK_FLUXES || defined BULK_FLUXES2D
           WRITE (out,200) blk_ZQ(ng), 'blk_ZQ',                         &
      &          'Height (m) of surface air humidity measurement.'
           IF (blk_ZQ(ng).le.0.0_r8) THEN