Changeset 14358

Timestamp:

2021-01-28T12:05:04+01:00 (3 years ago)

Author:

cbricaud

Message:

add ISF top limit condition for GLS in dev_r14318_MOI-02_MOI-GLS-ICE (ticket #2604 )

Location:

NEMO/branches/2021/dev_r14318_MOI-02_MOI-GLS-ICE/src/OCE/ZDF

Files:

• zdfgls.F90 (modified) (5 diffs)
• zdfphy.F90 (modified) (1 diff)

NEMO/branches/2021/dev_r14318_MOI-02_MOI-GLS-ICE/src/OCE/ZDF/zdfgls.F90

@@ -340 +340 @@
       zdiag(:,:,2) = 1._wp
       !
+      IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
+         DO_2D( 0, 0, 0, 0 )
+            IF( mikt(ji,jj) > 1 )THEN
+               itop   = mikt(ji,jj)       ! k   top w-point
+               itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
+               !
+               z_en = MAX( rc02r *  ustar2_top(ji,jj), rn_emin )
+               !
+               ! Dirichlet condition applied at: 
+               !     top level (itop)         &      Just below it (itopp1)   
+               zd_lw(ji,jj,itop) = 0._wp   ;   zd_lw(ji,jj,itopp1) = 0._wp
+               zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
+               zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = 1._wp
+               en   (ji,jj,itop) = z_en    ;   en   (ji,jj,itopp1) = z_en
+            ENDIF
+         END_2D
+      ENDIF
       !
       CASE ( 1 )             ! Neumann boundary condition (set d(e)/dz)
@@ -361 +378 @@
       en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm)
       !
+      IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
+         DO_2D( 0, 0, 0, 0 )
+            IF( mikt(ji,jj) > 1 )THEN
+               itop   = mikt(ji,jj)       ! k   top w-point
+               itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
+               !
+               zd_lw(ji,jj,itop) = en(ji,jj,itop)
+               zd_up(ji,jj,itop) = 0._wp
+               zdiag(ji,jj,itop) = 1._wp
+               !
+               zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_lw(ji,jj,itopp1) !Remove zd_lw from zdiag
+               zd_lw(ji,jj,itopp1) = 0._wp
+
+               zkar (ji,jj)   = (rl_sf +(vkarmn-rl_sf)*(1.-EXP(-rtrans*(gdept(ji,jj,itop,Kmm)-gdepw(ji,jj,itop,Kmm))/r_z0_top)))
+               zflxs(ji,jj)   = rsbc_tke2 * ustar2_top(ji,jj)**1.5_wp * zkar(ji,jj) &
+                            &                    * (  (r_z0_top+(gdept(ji,jj,itop,Kmm)-gdepw(ji,jj,itop,Kmm)) ) / r_z0_top)**(1.5_wp*ra_sf)
+               en(ji,jj,itopp1) = en(ji,jj,itopp1) + zflxs(ji,jj)/e3w(ji,jj,itopp1,Kmm)
+            ENDIF
+         END_2D
+      ENDIF
       !
       END SELECT
@@ -422 +459 @@
             en   (ji,jj,ibot) = z_en
          END_2D
+         !
          IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
             DO_2D( 0, 0, 0, 0 )
@@ -555 +593 @@
          zdiag(:,:,2) = 1._wp
          !
+         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
+            DO_2D( 0, 0, 0, 0 )
+               IF( mikt(ji,jj) > 1 )THEN
+                  itop   = mikt(ji,jj)       ! k   top w-point
+                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
+                  !
+                  zdep(ji,jj)         = vkarmn * r_z0_top
+                  psi  (ji,jj,itop)   = rc0**rpp * en(ji,jj,itop)**rmm*zdep(ji,jj)**rnn
+                  zd_lw(ji,jj,itop)   = 0._wp
+                  zd_up(ji,jj,itop)   = 0._wp
+                  zdiag(ji,jj,itop)   = 1._wp
+                  !
+                  zdep(ji,jj)         = vkarmn * ( r_z0_top + e3t(ji,jj,itopp1,Kmm) )
+                  psi  (ji,jj,itopp1) = rc0**rpp * en(ji,jj,itopp1)**rmm *zdep(ji,jj)**rnn
+                  zd_lw(ji,jj,itopp1) = 0._wp
+                  zd_up(ji,jj,itopp1) = 0._wp
+                  zdiag(ji,jj,itopp1) = 1._wp
+               ENDIF
+            END_2D
+         ENDIF
+         !
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
          !
@@ -579 +638 @@
          zflxs(:,:)   = zdep(:,:) * zflxs(:,:)
          psi  (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm)
+         !
+         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
+            DO_2D( 0, 0, 0, 0 )
+               IF( mikt(ji,jj) > 1 )THEN
+                  itop   = mikt(ji,jj)       ! k   top w-point
+                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
+                  !
+                  zdep(ji,jj) = vkarmn * r_z0_top
+                  psi  (ji,jj,itop) = rc0**rpp * en(ji,jj,itop)**rmm * zdep(ji,jj)**rnn
+                  zd_lw(ji,jj,itop) = psi(ji,jj,itop)
+                  zd_up(ji,jj,itop) = 0._wp
+                  zdiag(ji,jj,itop) = 1._wp
+                  !
+                  zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_lw(ji,jj,itopp1)
+                  zd_lw(ji,jj,itopp1) = 0._wp
+                  zdep(ji,jj) = vkarmn * ( r_z0_top + e3t(ji,jj,itopp1,Kmm) )
+                  zflxb = rsbc_psi2 * ( p_avm(ji,jj,itop) + p_avm(ji,jj,itopp1)) &
+                        &  * (0.5_wp*(en(ji,jj,itop)+en(ji,jj,itopp1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
+                  psi(ji,jj,itopp1) = psi(ji,jj,itopp1) + zflxb / e3w(ji,jj,itopp1,Kmm)
+               ENDIF
+            END_2D
+         ENDIF
+
          !
       END SELECT

NEMO/branches/2021/dev_r14318_MOI-02_MOI-GLS-ICE/src/OCE/ZDF/zdfphy.F90

@@ -207 +207 @@
       IF( ioptio /= 1 )    CALL ctl_stop( 'zdf_phy_init: one and only one vertical diffusion option has to be defined ' )
       IF( ln_isfcav ) THEN
-      IF( ln_zdfric .OR. ln_zdfgls )    CALL ctl_stop( 'zdf_phy_init: zdfric and zdfgls never tested with ice shelves cavities ' )
+      IF( ln_zdfric )    CALL ctl_stop( 'zdf_phy_init: zdfric never tested with ice shelves cavities ' )
       ENDIF
       !                                ! shear production term flag