Changeset 13851

Timestamp:

2020-11-23T15:24:10+01:00 (4 years ago)

Author:

emanuelaclementi

Message:

corrections for zdftke Neumann BC plus loop alignement - ticket #2155 #2339

File:

• NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/ZDF/zdftke.F90 (modified) (4 diffs)

NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/ZDF/zdftke.F90

@@ -239 +239 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
-      DO_2D( 0, 0, 0, 0 )         ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
-!! clem: this should be the right formulation but it makes the model unstable unless drags are calculated implicitly
-!!       one way around would be to increase zbbirau 
-!!          en(ji,jj,1) = MAX( rn_emin0, ( ( 1._wp - fr_i(ji,jj) ) * zbbrau + &
-!!             &                                     fr_i(ji,jj)   * zbbirau ) * taum(ji,jj) ) * tmask(ji,jj,1)
-         en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
-      END_2D
+      zd_up(ji,jj,1) = 0._wp
+      zd_lw(ji,jj,1) = 1._wp
+
+      IF ( cpl_phioc .and. ln_phioc )  THEN
+         SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves
+
+         CASE ( 0 ) ! Dirichlet BC
+
+            DO_2D( 0, 0, 0, 0 )    ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
+               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
+               en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) )  * tmask(ji,jj,1)
+               zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)  ! choose to keep coherence with former estimation of
+                                                   ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)
+            END_2D
+
+         CASE ( 1 ) ! Neumann BC
+            DO_2D( 0, 0, 0, 0 )
+               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
+               en(ji,jj,2)    = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 )/e3w(ji,jj,2,Kmm)
+               en(ji,jj,1)    = en(ji,jj,2) + (2 * e3t(ji,jj,1,Kmm) * phioc(ji,jj)/rho0) / ( p_avm(ji,jj,1) + p_avm(ji,jj,2) )
+               zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2)
+               zd_lw(ji,jj,2) = 0._wp
+               zdiag(ji,jj,1) = 1._wp
+            END_2D
+         
+         END SELECT
+
+      ELSE  ! TKE Dirichlet boundary condition (without wave coupling)
+
+         DO_2D( 0, 0, 0, 0 )       ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
+            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
+            zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) ! choose to keep coherence with former estimation of
+                                               ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)
+         END_2D
+
+      ENDIF
       !
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -400 +429 @@
       !                            !  Surface boundary condition on tke
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! [EC] Should we keep this??
+      ![EC] should it be removed?
       IF ( ln_isfcav ) THEN
          DO_2D( 1, 1, 1, 1 )       ! en(mikt(ji,jj))   = rn_emin
@@ -406 +435 @@
          END_2D
       END IF
-
-      IF ( cpl_phioc .and. ln_phioc )  THEN
-         SELECT CASE (nn_bc_surf)  ! Boundary Condition using surface TKE flux from waves
-
-         CASE ( 0 ) ! Dirichlet BC
-
-            DO_2D( 0, 0, 0, 0 )    ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
-               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
-               en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) )  * tmask(ji,jj,1)
-               zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)  ! choose to keep coherence with former estimation of
-               zd_lw(ji,jj,1) = 1._wp              ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)
-               zd_up(ji,jj,1) = 0._wp
-            END_2D
-
-         CASE ( 1 ) ! Neumann BC
-            DO_2D( 0, 0, 0, 0 )
-               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
-               en(ji,jj,2)    = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm)
-               en(ji,jj,1)    = en(ji,jj,2) + (2 * e3t(ji,jj,1,Kmm) * phioc(ji,jj)) / ( p_avm(ji,jj,1) + p_avm(ji,jj,2) )
-               zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2)
-               zd_lw(ji,jj,2) = 0._wp
-               zd_up(ji,jj,1) = 0._wp
-            END_2D
-!
-         END SELECT
-
-      ELSE  ! TKE Dirichlet boundary condition (without wave coupling)
-
-         DO_2D( 0, 0, 0, 0 )       ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
-            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
-            zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) ! choose to keep coherence with former estimation of
-            zd_lw(ji,jj,1) = 1._wp             ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)
-            zd_up(ji,jj,1) = 0._wp
-         END_2D
-
-      ENDIF
       !
       !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
@@ -771 +764 @@
          IF ( cpl_phioc .and. ln_phioc )  THEN
             SELECT CASE( nn_bc_surf)             ! Type of scaling under sea-ice
-            CASE( 0 )   ;   WRITE(numout,*) '  nn_bc_surf=0  ==>>>   DIRICHLET Surface Boundary Condition using surface TKE flux from waves'
-            CASE( 1 )   ;   WRITE(numout,*) '  nn_bc_surf=1  ==>>>   NEUMANN Surface Boundary Condition using surface TKE flux from waves'
+            CASE( 0 )   ;   WRITE(numout,*) '  nn_bc_surf=0 ==>>> DIRICHLET SBC using surface TKE flux from waves'
+            CASE( 1 )   ;   WRITE(numout,*) '  nn_bc_surf=1 ==>>> NEUMANN SBC using surface TKE flux from waves'
             END SELECT
-         ENDIF 
+         ENDIF
          WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
          WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau