Changeset 14566

Timestamp:

2021-03-02T11:42:20+01:00 (3 years ago)

Author:

smueller

Message:

Synchronisation of the OSMOSIS boundary layer scheme with the version developed in branch ^{/NEMO/branches/NERC/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0: transfer of changeset [14516] (ticket #2353)}

File:

• NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90 (modified) (12 diffs)

NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90

@@ -550 +550 @@
       DO_2D( 0, 0, 0, 0 )
          zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm)
-         imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji, jj, ibld(ji,jj), Kmm )) , 1 ))
+         imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji, jj, ibld(ji,jj) - 1, Kmm )) , 1 ))
          zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
          zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
@@ -582 +582 @@
          END_2D
 
-!! External gradient
-         CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )
+!! Calculate fairly-well-mixed depth zmld & its index mld_prof + lateral zmld-averaged gradients
          CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle )
+!! Calculate vertical gradients immediately below zmld
          CALL zdf_osm_external_gradients( mld_prof, zdtdz_mle_ext, zdsdz_mle_ext, zdbdz_mle_ext )
+!! Calculate max vertical FK flux zwb_fk & set logical descriptors
          CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk )
+!! recalculate hmle, zmle, zvel_mle, zdiff_mle & redefine mld_proc to be index for new hmle
          CALL zdf_osm_mle_parameters( zmld, mld_prof, hmle, zhmle, zvel_mle, zdiff_mle )
       ELSE    ! ln_osm_mle
@@ -597 +599 @@
          END_2D
       ENDIF   ! ln_osm_mle
+
+!! External gradient below BL needed both with and w/o FK
+      CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )
 
 ! Test if pycnocline well resolved
@@ -615 +620 @@
       END_2D
 
+      ! Recalculate bl averages using jp_ext & ml averages .... note no rotation of u & v here..
       CALL zdf_osm_vertical_average( Kbb, Kmm,                                          &
          &                           ibld, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl,           &
@@ -648 +654 @@
 ! is external level in bounds?
 
+!   Recalculate BL averages and differences using new BL depth
       CALL zdf_osm_vertical_average( Kbb, Kmm,                                          &
          &                           ibld, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl,           &
@@ -663 +670 @@
       dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. )  !  Limit delta for shallow boundary layers for calculating flux-gradient terms.
 !
-    ! Average over the depth of the mixed layer in the convective boundary layer
+      ! Average over the depth of the mixed layer in the convective boundary layer
 !      jp_ext = ibld - imld +1
+!     Recalculate ML averages and differences using new ML depth
       CALL zdf_osm_vertical_average( Kbb, Kmm,                                               &
          &                           imld-1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml,              &
          &                           ibld-imld+1, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml )
-    ! rotate mean currents and changes onto wind align co-ordinates
-    !
-     CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml )
-     CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl )
+      ! rotate mean currents and changes onto wind align co-ordinates
+      !
+      CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml )
+      CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl )
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       !  Pycnocline gradients for scalars and velocity
@@ -993 +1001 @@
                 DO jk = imld(ji,jj) , ibld(ji,jj)
                    zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6)
-                    zviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 )
-                    zviscos(ji,jj,jk) = zviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 -0.2 * zznd_pyc**3 )
+                   zviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 )
+                   zviscos(ji,jj,jk) = zviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 -0.2 * zznd_pyc**3 )
                 END DO
                 IF ( zdhdt(ji,jj) > 0._wp ) THEN
@@ -1338 +1346 @@
            jkb1 = MIN(jkb + 1, mbkt(ji,jj))
            zdtdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_tem,Kmm) - ts(ji,jj,jkb,jp_tem,Kmm ) ) &
-                &   / e3t(ji,jj,jkb,Kmm)
+                &   / e3w(ji,jj,jkb1,Kmm)
            zdsdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_sal,Kmm) - ts(ji,jj,jkb,jp_sal,Kmm ) ) &
-                &   / e3t(ji,jj,jkb,Kmm)
+                &   / e3w(ji,jj,jkb1,Kmm)
            zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj)
         ELSE
@@ -1634 +1642 @@
                  IF ( ln_osm_mle ) THEN
                     zhbl_s = zhbl_s + MIN( &
-                     & rn_Dt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
-                     & e3w(ji,jj,jm,Kmm) )
+                       & rn_Dt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
+                       & e3w(ji,jj,jm,Kmm) )
                  ELSE
-                   zhbl_s = zhbl_s + MIN( &
-                     & rn_Dt * (  -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
-                     & e3w(ji,jj,jm,Kmm) )
+                    zhbl_s = zhbl_s + MIN( &
+                       & rn_Dt * (  -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
+                       & e3w(ji,jj,jm,Kmm) )
                  ENDIF
 
-!                    zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol)
+                 !                    zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol)
                  IF ( zhbl_s >= gdepw(ji,jj,mbkt(ji,jj) + 1,Kmm) ) THEN
-                   zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol)
-                   lpyc(ji,jj) = .FALSE.
+                    zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol)
+                    lpyc(ji,jj) = .FALSE.
                  ENDIF
                  IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1
@@ -1680 +1688 @@
 ! change zero or one model level.
           hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw(ji,jj,4,Kmm) )
-        ENDIF
-        zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm)
+       ENDIF
+       zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm)
      END_2D
      IF( ln_timing ) CALL timing_stop('zdf_osm_th')
@@ -1754 +1762 @@
                zdh_ref = 0.2 * hbl(ji,jj)
             ENDIF    ! IF (dhdt >= 0)
-            dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) )
+            dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) )
             IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref       ! can be a problem with dh>hbl for rapid collapse
-            ! Alan: this hml is never defined or used -- do we need it?
          ENDIF
 
@@ -1824 +1831 @@
 
       hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)
-      inhml = MAX( INT( dh(ji,jj) / MAX(e3t(ji,jj,ibld(ji,jj),Kmm), 1.e-3) ) , 1 )
+      inhml = MAX( INT( dh(ji,jj) / MAX(e3t(ji,jj,ibld(ji,jj) - 1,Kmm), 1.e-3) ) , 1 )
       imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3)
       zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)