Changeset 14441

Timestamp:

2021-02-11T17:34:32+01:00 (3 years ago)

Author:

agn

Message:

Included Alan's fix for shear code

File:

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

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

@@ -2103 +2103 @@
                    zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
                 ENDIF
-                zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
+                zvel_max = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
                 IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN
                    ! OSBL is deepening, entrainment > restratification
-                   IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN
-                      zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj)
-                      zpsi = -zalpha_pyc(ji,jj) * ( zwb0tot(ji,jj) - ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) ) * zdh(ji,jj) / zhbl(ji,jj)
-                      zpsi = zpsi - 4.0 * ( 1.0 + zdh(ji,jj) /zhbl(ji,jj) ) * zgamma_b_nd * ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) )
+                   IF ( zdb_bl(ji,jj) > 1.0e-15 ) THEN
+                      zgamma_b_nd = MAX( zdbdz_bl_ext(ji,jj), 0._wp ) * zdh(ji,jj) / zdb_ml(ji,jj)
+                      zpsi = ( 1.0 - 0.5 * zdh(ji,jj) / zhbl(ji,jj) ) * ( zwb0(ji,jj) - MIN( ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ), 0._wp ) ) * zdh(ji,jj) / zhbl(ji,jj)
+                      zpsi = zpsi + 1.75 * ( 1.0 - 0.5 * zdh(ji,jj) / zhbl(ji,jj) )*( zdh(ji,jj) / zhbl(ji,jj) + zgamma_b_nd ) * MIN( ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ), 0._wp )
                       zpsi = zalpha_b * MAX ( zpsi, 0._wp )
-                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) + zpsi / MAX( zdb_bl(ji,jj), 1.e-15 )
+                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) + zpsi / ( zvel_max + MAX( zdb_bl(ji,jj), 1.e-15 ) )
                       IF ( j_ddh(ji,jj) == 1 ) THEN
                         IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
@@ -2128 +2128 @@
                         zddhdt = 0._wp
                       ENDIF ! j_ddh
-                      zdhdt(ji,jj) = zdhdt(ji,jj) + zalpha_b * ( zalpha_pyc(ji,jj) * zdb_ml(ji,jj) + 2.0 * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) ) * zddhdt / zdb_bl(ji,jj)
+                      zdhdt(ji,jj) = zdhdt(ji,jj) + zalpha_b * ( 1.0 -0.5 * zdh(ji,jj) / zhbl(ji,jj) ) * zddhdt / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) )
                    ELSE    ! zdb_bl >0
                        zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15)
@@ -2334 +2334 @@
       INTEGER :: jj, ji
       INTEGER :: inhml
-      REAL(wp) :: zari, ztau, zdh_ref, zddhdt
+      REAL(wp) :: zari, ztau, zdh_ref, zddhdt, zvel_max
       REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth
 
@@ -2342 +2342 @@
          IF ( lshear(ji,jj) ) THEN
             IF ( lconv(ji,jj) ) THEN
-             IF ( zdb_bl(ji,jj) > 0._wp ) THEN
+             IF ( zdb_bl(ji,jj) > 1.0e-15) THEN
               IF ( j_ddh(ji,jj) == 0 ) THEN
+                zvel_max = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
 ! ddhdt for pycnocline determined in osm_calculate_dhdt
-!                zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)
-!                zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt
+                zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) )
+                zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt
 ! maximum limit for how thick the shear layer can grow relative to the thickness of the boundary kayer
-!                dh(ji,jj) = MIN( dh(ji,jj) + zddhdt * rn_rdt, 0.625 * hbl(ji,jj) )
-                ztau = MAX( zdb_bl(ji,jj) * ( 0.625 * hbl(ji,jj) ) / ( a_ddh * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_rdt )
-                dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + 0.625 * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) )
-                IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = 0.8 * zhbl(ji,jj)
+                dh(ji,jj) = MIN( dh(ji,jj) + zddhdt * rn_rdt, 0.625 * hbl(ji,jj) )
               ELSE
-! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt 
+! Need to recalculate because hbl has been updated.
                 IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
                   zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )