Changeset 13396

Timestamp:

2020-08-14T12:07:37+02:00 (4 years ago)

Author:

agn

Message:

Ensured no division by zero in zdf_osm_pycnocline_thickness etc.

File:

• NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90 (modified) (6 diffs)

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90

@@ -1556 +1556 @@
 
                   zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &
-                  &        ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
+                  &        MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln)
                   zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
                 ! added ajgn 23 July as temporay fix
@@ -1589 +1589 @@
                     zpert = MAX( 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )
                 ENDIF
-                zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert
+                zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln)
             ENDIF
          END DO
@@ -1722 +1722 @@
                      ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F
                      IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN
-                        IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability
+                        IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN  ! near neutral stability
                            zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
-                                (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
+                                (1.0 + zdb_bl(ji,jj)**2 / MAX( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj), epsln ) ), 0.2 )
                         ELSE                                                     ! unstable
                            zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
                                 (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
                         ENDIF
-                        ztau = 0.2 * hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird
+                        ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)
                         zddhdt = zari * hbl(ji,jj)
                      ELSE
-                        ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
+                        ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)
                         zddhdt = 0.2 * hbl(ji,jj)
                      ENDIF
                   ELSE
-                     ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
+                     ztau = 0.2 * hbl(ji,jj) /  MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)
                      zddhdt = 0.2 * hbl(ji,jj)
                   ENDIF
                ELSE ! ln_osm_mle
                   IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN
-                     IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability
-                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
-                             (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
+                     IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN  ! near neutral stability
+                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / MAX( zdbdz_ext(ji,jj) * zhbl(ji,jj), epsln ) ) / &
+                             (1.0 + zdb_bl(ji,jj)**2 /  MAX(4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj), epsln ) ), 0.2 )
                      ELSE                                                     ! unstable
-                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
-                             (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
+                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / MAX( zdbdz_ext(ji,jj) * zhbl(ji,jj), epsln ) ) / &
+                             (1.0 + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj), epsln ) ), 0.2 )
                      ENDIF
-                     ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird
+                     ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)
                      zddhdt = zari * hbl(ji,jj)
                   ELSE
-                     ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
+                     ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)
                      zddhdt = 0.2 * hbl(ji,jj)
                   ENDIF
@@ -1760 +1760 @@
                ! Alan: this hml is never defined or used
             ELSE   ! IF (lconv)
-               ztau = hbl(ji,jj) / zvstr(ji,jj)
+               ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln)
                IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here
                   ! boundary layer deepening
@@ -1766 +1766 @@
                      ! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
                      zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
-                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 )
+                          & /  MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01  , 0.2 )
                      zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj)
                   ELSE
@@ -1969 +1969 @@
                   & + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) )
              zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle
-      ! This vbelocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt.
+      ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt.
              zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 
              zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3  / rn_osm_mle_lf