Changeset 12384

Timestamp:

2020-02-13T17:21:21+01:00 (3 years ago)

Author:

agn

Message:

corrections to tramle & zdfosm to avoid division by 0

Location:

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE

Files:

• TRA/tramle.F90 (modified) (1 diff)
• ZDF/zdfosm.F90 (modified) (4 diffs)

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/tramle.F90

@@ -254 +254 @@
     DO jj = 1, jpjm1
        DO ji = 1, fs_jpim1   ! vector opt.
-          zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu
-          zhv(ji,jj) = 1._wp / zhv(ji,jj)
+          zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall)                   ! hu --> 1/hu
+          zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall) 
        END DO
     END DO

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdfosm.F90

@@ -1358 +1358 @@
      REAL(wp) :: ztgrad, zsgrad, zbgrad
      REAL(wp) :: zgamma_b_nd, zgamma_c, znd
-     REAL(wp) :: zzeta_s=0.3
+     REAL(wp) :: zzeta_s=0.3, ztmp
 
      DO jj = 2, jpjm1
@@ -1367 +1367 @@
                          & (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh &
                          &  .OR.  zdb_bl(ji,jj) > 0._wp)) THEN  ! zdhdt could be <0 due to FK, hence check zdhdt>0
-                       ztgrad = 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_ext(ji,jj)
-                       zsgrad = 0.5 * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_ext(ji,jj)
-                       zbgrad = 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_ext(ji,jj)
-                       zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj)
+                       ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
+                       ztgrad = 0.5 * zdt_ml(ji,jj) * ztmp + zdtdz_ext(ji,jj)
+                       zsgrad = 0.5 * zds_ml(ji,jj) * ztmp + zdsdz_ext(ji,jj)
+                       zbgrad = 0.5 * zdb_ml(ji,jj) * ztmp + zdbdz_ext(ji,jj)
+                       zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln)
                        zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /&
                             & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check
                        DO jk = 2, ibld(ji,jj)+ibld_ext
-                          znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
+                          znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp
                           IF ( znd <= zzeta_s ) THEN
-                             zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) * &
+                             zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) * ztmp * &
                                   & EXP( -6.0 * ( znd -zzeta_s )**2 )
-                             zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) / zdh(ji,jj) * &
+                             zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) * ztmp * &
                                   & EXP( -6.0 * ( znd -zzeta_s )**2 )
-                             zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) * &
+                             zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) * ztmp * &
                                   & EXP( -6.0 * ( znd -zzeta_s )**2 )
                           ELSE
@@ -1395 +1396 @@
                     IF ( zdb_bl(ji,jj) > 0._wp ) THEN
                        IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline
-                          ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj)
-                          zsgrad = zds_bl(ji,jj) / zhbl(ji,jj)
-                          zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj)
+                          ztmp = 1._wp/MAX(zhbl(ji,jj), epsln)
+                          ztgrad = zdt_bl(ji,jj) * ztmp
+                          zsgrad = zds_bl(ji,jj) * ztmp
+                          zbgrad = zdb_bl(ji,jj) * ztmp
                           DO jk = 2, ibld(ji,jj)
-                             znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
+                             znd = gdepw_n(ji,jj,jk) * ztmp
                              zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
                              zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
@@ -1405 +1407 @@
                           END DO
                        ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
-                          ztgrad = zdt_bl(ji,jj) / zdh(ji,jj)
-                          zsgrad = zds_bl(ji,jj) / zdh(ji,jj)
-                          zbgrad = zdb_bl(ji,jj) / zdh(ji,jj)
+                          ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
+                          ztgrad = zdt_bl(ji,jj) * ztmp
+                          zsgrad = zds_bl(ji,jj) * ztmp
+                          zbgrad = zdb_bl(ji,jj) * ztmp
                           DO jk = 2, ibld(ji,jj)
-                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
+                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) * ztmp
                              zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
                              zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )