Changeset 14357 for NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/LDF/ldfslp.F90

Timestamp:

2021-01-28T11:45:54+01:00 (3 years ago)

Author:

cbricaud

Message:

add in dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING branch ISF top limit condition for triads (ticket #2601)

File:

• NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/LDF/ldfslp.F90 (modified) (7 diffs)

NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/LDF/ldfslp.F90

@@ -148 +148 @@
          END_2D
       ENDIF
-      IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the bottom ocean level
+      IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the top ocean level
          DO_2D( 1, 0, 1, 0 )
             IF( miku(ji,jj) > 1 )   zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 
@@ -356 +356 @@
             END_2D
          ENDIF
+         IF( ln_zps .AND. ln_isfcav .AND. l_grad_zps) THEN           ! partial steps correction at the top ocean level
+            DO_2D( 1, 0, 1, 0 )
+               iku  = miku(ji,jj)          ;   ikv  = mikv(ji,jj)        ! first wet ocean level (u- & v-points)
+               zdit = gtui(ji,jj,jp_tem)   ;   zdjt = gtvi(ji,jj,jp_tem) ! i- & j-gradient of Temperature
+               zdis = gtui(ji,jj,jp_sal)   ;   zdjs = gtvi(ji,jj,jp_sal) ! i- & j-gradient of Salinity
+               zdxrho_raw = ( - rab_b(ji+ip,jj   ,iku,jp_tem) * zdit + rab_b(ji+ip,jj   ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj)
+               zdyrho_raw = ( - rab_b(ji   ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji   ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj)
+               zdxrho(ji+ip,jj   ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign
+               zdyrho(ji   ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw )
+            END_2D
+         ENDIF
          !
       END DO
@@ -368 +379 @@
                zdks = 0._wp
             ENDIF
+            IF ( ln_isfcav ) THEN
+               IF( mikt(ji,jj) > 1 .AND.  jk+kp <= mikt(ji,jj) )THEN
+                  zdkt = 0._wp                                          ! 1st wet level gradient set to zero
+                  zdks = 0._wp
+               ENDIF
+            ENDIF
+
             zdzrho_raw = ( - rab_b(ji,jj,jk   ,jp_tem) * zdkt & 
                        &   + rab_b(ji,jj,jk   ,jp_sal) * zdks &
@@ -382 +400 @@
       !                                       !==  intialisations to zero  ==!
       !
-      wslp2  (:,:,:)     = 0._wp              ! wslp2 will be cumulated 3D field set to zero
-      triadi_g(:,:,1,:,:) = 0._wp   ;   triadi_g(:,:,jpk,:,:) = 0._wp   ! set surface and bottom slope to zero
-      triadj_g(:,:,1,:,:) = 0._wp   ;   triadj_g(:,:,jpk,:,:) = 0._wp
-      !!gm _iso set to zero missing
-      triadi  (:,:,1,:,:) = 0._wp   ;   triadj  (:,:,jpk,:,:) = 0._wp   ! set surface and bottom slope to zero
-      triadj  (:,:,1,:,:) = 0._wp   ;   triadj  (:,:,jpk,:,:) = 0._wp
+      wslp2   (:,:,:)     = 0._wp              ! wslp2 will be cumulated 3D field
+      triadi_g(:,:,:,:,:) = 0._wp ; triadj_g(:,:,:,:,:) = 0._wp
+      triadi  (:,:,:,:,:) = 0._wp ; triadj  (:,:,:,:,:) = 0._wp
 
       !-------------------------------------!
@@ -445 +460 @@
                   !
                   ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface
-                  zti_coord = znot_thru_surface * ( gdept(ji+1,jj  ,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj)
-                  ztj_coord = znot_thru_surface * ( gdept(ji  ,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj)     ! unmasked
+                  zti_coord =  wmask(ji,jj,jk+kp) * znot_thru_surface * ( gdept(ji+1,jj  ,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj)
+                  ztj_coord =  wmask(ji,jj,jk+kp) * znot_thru_surface * ( gdept(ji  ,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj)     ! unmasked
+
                   zti_g_raw = zti_raw - zti_coord      ! ref to geopot surfaces
                   ztj_g_raw = ztj_raw - ztj_coord
@@ -464 +480 @@
                   zti_g_lim =          ( zfacti   * zti_g_lim                       &
                      &      + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp)   &
-                     &                           * gdepw(ji+ip,jj,jk+kp,Kmm) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp)
+                     &                           * (gdepw(ji+ip,jj,jk+kp,Kmm)-risfdep(ji+ip,jj)) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp)
                   ztj_g_lim =          ( zfactj   * ztj_g_lim                       &
                      &      + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp)   &
-                     &                           * gdepw(ji,jj+jp,jk+kp,Kmm) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp)
+                     &                           * (gdepw(ji,jj+jp,jk+kp,Kmm)-risfdep(ji,jj+jp)) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp)
                   !
                   triadi_g(ji+ip,jj   ,jk,1-ip,kp) = zti_g_lim
@@ -478 +494 @@
                   ! To do this, retain limited sx**2  in vertical flux, but divide by real slope for 13/31 terms
                   ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2
+                  !
+                  zti_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS(zti_lim ) ), zti_lim )
+                  ztj_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS(ztj_lim ) ), ztj_lim )
                   !
                   zti_lim  = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp)    ! remove coordinate slope => relative to coordinate surfaces