Changeset 5044 for branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

Timestamp:

2015-01-26T15:48:14+01:00 (9 years ago)

Author:

hliu

Message:

correct three potential bugs in vertical scale calculation when wetting/drying is on

File:

• branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (9 diffs)

branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

@@ -1564 +1564 @@
          DO jj = 1, jpj
 
-            IF( hbatt(ji,jj) > 0._wp) THEN
+            IF( bathy(ji,jj) > 0._wp) THEN
                DO jk = 1, mbathy(ji,jj)
                  ! check coordinate is monotonically increasing
@@ -1628 +1628 @@
       INTEGER  ::   ji, jj, jk           ! dummy loop argument
       REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
-      REAL(wp) ::   ztmpu, ztmpv, ztmpf
+      REAL(wp) ::   ztmpu,  ztmpv,  ztmpf
+      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
       !
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
@@ -1686 +1687 @@
             ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
             ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
+            ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
+            ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
+            ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
+                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
             DO jk = 1, jpk
-               IF(ABS(ztmpu) < rn_wdmin1.AND.ln_wd) THEN
+               IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN
                  z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
                  z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
@@ -1697 +1702 @@
                END IF
 
-               IF(ABS(ztmpv) < rn_wdmin1.AND.ln_wd) THEN
+               IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN
                  z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
                  z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
@@ -1707 +1712 @@
                END IF
 
-               IF(ABS(ztmpf) < rn_wdmin1.AND.ln_wd) THEN
+               IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN
                  z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)  &
                     &                            + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) )
@@ -1768 +1773 @@
       REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
       REAL(wp) ::   ztmpu, ztmpv, ztmpf
+      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
       !
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
@@ -1850 +1856 @@
            ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
            ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
+           ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
+           ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
+           ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
+                  & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
            DO jk = 1, jpk
-              IF(ABS(ztmpu) < 1.e-5.AND.ln_wd) THEN
+              IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN
                 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
                 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
@@ -1861 +1871 @@
               END IF
 
-              IF(ABS(ztmpv) < 1.e-5.AND.ln_wd) THEN
+              IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN
                 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
                 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
@@ -1871 +1881 @@
               END IF
 
-              IF(ABS(ztmpf) < 1.e-5.AND.ln_wd) THEN
+              IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN
                 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)  &
                    &                            + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) )