Changeset 5837 for branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF

Timestamp:

2015-10-26T15:59:39+01:00 (9 years ago)

Author:

timgraham

Message:

Upgraded to r5518 of trunk (NEMO 3.6 stable)

Location:

branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF

Files:

• zdf_oce.F90 (modified) (2 diffs)
• zdfbfr.F90 (modified) (15 diffs)
• zdfddm.F90 (modified) (14 diffs)
• zdfevd.F90 (modified) (2 diffs)
• zdfgls.F90 (modified) (31 diffs)
• zdfini.F90 (modified) (1 diff)
• zdfkpp.F90 (modified) (12 diffs)
• zdfmxl.F90 (modified) (6 diffs)
• zdftke.F90 (modified) (22 diffs)
• zdftmx.F90 (modified) (15 diffs)

branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90

@@ -40 +40 @@
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   avtb_2d        !: horizontal shape of background Kz profile
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   bfrua, bfrva   !: Bottom friction coefficients set in zdfbfr
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   tfrua, tfrva   !: top friction coefficients set in zdfbfr
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avmu , avmv    !: vertical viscosity coef at uw- & vw-pts       [m2/s]
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avm  , avt     !: vertical viscosity & diffusivity coef at w-pt [m2/s]
@@ -57 +58 @@
       ALLOCATE(avmb(jpk) , bfrua(jpi,jpj) ,                         &
          &     avtb(jpk) , bfrva(jpi,jpj) , avtb_2d(jpi,jpj) ,      &
+         &     tfrua(jpi, jpj), tfrva(jpi, jpj)              ,      &
          &     avmu(jpi,jpj,jpk), avm(jpi,jpj,jpk)           ,      &
          &     avmv(jpi,jpj,jpk), avt(jpi,jpj,jpk)           , STAT = zdf_oce_alloc )

branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

@@ -41 +41 @@
    REAL(wp), PUBLIC ::   rn_bfrien    ! local factor to enhance coefficient bfri (PUBLIC for TAM)
    LOGICAL , PUBLIC ::   ln_bfr2d     ! logical switch for 2D enhancement (PUBLIC for TAM)
+   REAL(wp), PUBLIC ::   rn_tfri1     ! top drag coefficient (linear case)  (PUBLIC for TAM)
+   REAL(wp), PUBLIC ::   rn_tfri2     ! top drag coefficient (non linear case) (PUBLIC for TAM)
+   REAL(wp), PUBLIC ::   rn_tfri2_max ! Maximum top drag coefficient (non linear case and ln_loglayer=T) (PUBLIC for TAM)
+   REAL(wp), PUBLIC ::   rn_tfeb2     ! background top turbulent kinetic energy  [m2/s2] (PUBLIC for TAM)
+   REAL(wp), PUBLIC ::   rn_tfrien    ! local factor to enhance coefficient tfri (PUBLIC for TAM)
+   LOGICAL , PUBLIC ::   ln_tfr2d     ! logical switch for 2D enhancement (PUBLIC for TAM)
+
    LOGICAL , PUBLIC ::   ln_loglayer  ! switch for log layer bfr coeff. (PUBLIC for TAM)
    REAL(wp), PUBLIC ::   rn_bfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM)
+   REAL(wp), PUBLIC ::   rn_tfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM)
    LOGICAL , PUBLIC ::   ln_bfrimp    ! logical switch for implicit bottom friction
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d            ! 2D bottom drag coefficient (PUBLIC for TAM)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d, tfrcoef2d   ! 2D bottom/top drag coefficient (PUBLIC for TAM)
 
    !! * Substitutions
@@ -60 +68 @@
       !!                ***  FUNCTION zdf_bfr_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( bfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc )
+      ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc )
       !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_bfr_alloc )
@@ -88 +96 @@
       INTEGER  ::   ikbt, ikbu, ikbv             ! local integers
       REAL(wp) ::   zvu, zuv, zecu, zecv, ztmp   ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:) ::  zbfrt
+      REAL(wp), POINTER, DIMENSION(:,:) ::  zbfrt, ztfrt
       !!----------------------------------------------------------------------
       !
@@ -101 +109 @@
       IF( nn_bfr == 2 ) THEN                 ! quadratic bottom friction only
          !
-         CALL wrk_alloc( jpi, jpj, zbfrt )
+         CALL wrk_alloc( jpi, jpj, zbfrt, ztfrt )
 
          IF ( ln_loglayer.AND.lk_vvl ) THEN ! "log layer" bottom friction coefficient
 
-#  if defined key_vectopt_loop
-            DO jj = 1, 1
-!CDIR NOVERRCHK
-               DO ji = 1, jpij   ! vector opt. (forced unrolling)
-#  else
-!CDIR NOVERRCHK
             DO jj = 1, jpj
-!CDIR NOVERRCHK
                DO ji = 1, jpi
-#  endif
                   ikbt = mbkt(ji,jj)
-! JC: possible WAD implementation should modify line below if layers vanish
+!! JC: possible WAD implementation should modify line below if layers vanish
                   ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp
                   zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp)
@@ -122 +122 @@
                END DO
             END DO
+! (ISF)
+            IF ( ln_isfcav ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ikbt = mikt(ji,jj)
+! JC: possible WAD implementation should modify line below if layers vanish
+                     ztmp = (1-tmask(ji,jj,1)) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp
+                     ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp)
+                     ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max)
+                  END DO
+               END DO
+            END IF
          !   
          ELSE
             zbfrt(:,:) = bfrcoef2d(:,:)
-         ENDIF
-
-# if defined key_vectopt_loop
-         DO jj = 1, 1
-!CDIR NOVERRCHK
-            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-# else
-!CDIR NOVERRCHK
+            ztfrt(:,:) = tfrcoef2d(:,:)
+         ENDIF
+
          DO jj = 2, jpjm1
-!CDIR NOVERRCHK
             DO ji = 2, jpim1
-# endif
                ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
                ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
@@ -150 +155 @@
                bfrua(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) * zecu
                bfrva(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) * zecv
+               !
+               ! in case of 2 cell water column, we assume each cell feels the top and bottom friction
+               IF ( ln_isfcav ) THEN
+                  IF ( miku(ji,jj) + 1 .GE. mbku(ji,jj) ) THEN
+                     bfrua(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) )   &
+                                  &            + ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) ) &
+                                  &          * zecu * (1._wp - umask(ji,jj,1))
+                  END IF
+                  IF ( mikv(ji,jj) + 1 .GE. mbkv(ji,jj) ) THEN
+                     bfrva(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) )   &
+                                  &            + ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) ) &
+                                  &          * zecv * (1._wp - vmask(ji,jj,1))
+                  END IF
+               END IF
             END DO
          END DO
-
-         !
          CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition
+
+         IF ( ln_isfcav ) THEN
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  ! (ISF) ========================================================================
+                  ikbu = miku(ji,jj)         ! ocean top level at u- and v-points 
+                  ikbv = mikv(ji,jj)         ! (1st wet ocean u- and v-points)
+                  !
+                  zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     &
+                     &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  )
+                  zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     &
+                     &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  )
+              !
+                  zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_tfeb2 )
+                  zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_tfeb2 )
+              !
+                  tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu * (1._wp - umask(ji,jj,1))
+                  tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) * zecv * (1._wp - vmask(ji,jj,1))
+              ! (ISF) END ====================================================================
+              ! in case of 2 cell water column, we assume each cell feels the top and bottom friction
+                  IF ( miku(ji,jj) + 1 .GE. mbku(ji,jj) ) THEN
+                     tfrua(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) )   &
+                                  &            + ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) ) &
+                                  &          * zecu * (1._wp - umask(ji,jj,1))
+                  END IF
+                  IF ( mikv(ji,jj) + 1 .GE. mbkv(ji,jj) ) THEN
+                     tfrva(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) )   &
+                                  &            + ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) ) &
+                                  &          * zecv * (1._wp - vmask(ji,jj,1))
+                  END IF
+               END DO
+            END DO
+            CALL lbc_lnk( tfrua, 'U', 1. )   ;   CALL lbc_lnk( tfrva, 'V', 1. )      ! Lateral boundary condition
+         END IF
+         !
          !
          IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        &
             &                       tab2d_2=bfrva, clinfo2=       ' v: ', mask2=vmask,ovlap=1 )
-         CALL wrk_dealloc( jpi,jpj,zbfrt )
+         CALL wrk_dealloc( jpi,jpj,zbfrt, ztfrt )
       ENDIF
       !
@@ -183 +235 @@
       INTEGER   ::   ios
       REAL(wp)  ::   zminbfr, zmaxbfr   ! temporary scalars
+      REAL(wp)  ::   zmintfr, zmaxtfr   ! temporary scalars
       REAL(wp)  ::   ztmp, zfru, zfrv   !    -         -
       !!
       NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfri2_max, rn_bfeb2, rn_bfrz0, ln_bfr2d, &
-                    &  rn_bfrien, ln_bfrimp, ln_loglayer
+                    &          rn_tfri1, rn_tfri2, rn_tfri2_max, rn_tfeb2, rn_tfrz0, ln_tfr2d, &
+                    &  rn_bfrien, rn_tfrien, ln_bfrimp, ln_loglayer
       !!----------------------------------------------------------------------
       !
@@ -215 +269 @@
          bfrua(:,:) = 0.e0
          bfrva(:,:) = 0.e0
+         tfrua(:,:) = 0.e0
+         tfrva(:,:) = 0.e0
          !
       CASE( 1 )
          IF(lwp) WRITE(numout,*) '      linear botton friction'
-         IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri1  = ', rn_bfri1
+         IF(lwp) WRITE(numout,*) '      bottom friction coef.   rn_bfri1  = ', rn_bfri1
          IF( ln_bfr2d ) THEN
             IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d
             IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien
          ENDIF
+         IF ( ln_isfcav ) THEN
+            IF(lwp) WRITE(numout,*) '      top    friction coef.   rn_bfri1  = ', rn_tfri1
+            IF( ln_tfr2d ) THEN
+               IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d
+               IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien
+            ENDIF
+         END IF
          !
          IF(ln_bfr2d) THEN
@@ -236 +299 @@
          bfrua(:,:) = - bfrcoef2d(:,:)
          bfrva(:,:) = - bfrcoef2d(:,:)
+         !
+         IF ( ln_isfcav ) THEN
+            IF(ln_tfr2d) THEN
+               ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
+               CALL iom_open('tfr_coef.nc',inum)
+               CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array
+               CALL iom_close(inum)
+               tfrcoef2d(:,:) = rn_tfri1 * ( 1 + rn_tfrien * tfrcoef2d(:,:) )
+            ELSE
+               tfrcoef2d(:,:) = rn_tfri1  ! initialize tfrcoef2d to the namelist variable
+            ENDIF
+            !
+            tfrua(:,:) = - tfrcoef2d(:,:)
+            tfrva(:,:) = - tfrcoef2d(:,:)
+         END IF
          !
       CASE( 2 )
@@ -252 +330 @@
             IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien
          ENDIF
+         IF ( ln_isfcav ) THEN
+            IF(lwp) WRITE(numout,*) '      quadratic top    friction'
+            IF(lwp) WRITE(numout,*) '      friction coef.    rn_tfri2     = ', rn_tfri2
+            IF(lwp) WRITE(numout,*) '      Max. coef. (log case)   rn_tfri2_max  = ', rn_tfri2_max
+            IF(lwp) WRITE(numout,*) '      background tke    rn_tfeb2     = ', rn_tfeb2
+            IF(lwp) WRITE(numout,*) '      log formulation   ln_tfr2d     = ', ln_loglayer
+            IF(lwp) WRITE(numout,*) '      top roughness     rn_tfrz0 [m] = ', rn_tfrz0
+            IF( rn_tfrz0<=0.e0 ) THEN
+               WRITE(ctmp1,*) '      top roughness must be strictly positive'
+               CALL ctl_stop( ctmp1 )
+            ENDIF
+            IF( ln_tfr2d ) THEN
+               IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d
+               IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien
+            ENDIF
+         END IF
          !
          IF(ln_bfr2d) THEN
@@ -263 +357 @@
             bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable
          ENDIF
+         
+         IF ( ln_isfcav ) THEN
+            IF(ln_tfr2d) THEN
+               ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
+               CALL iom_open('tfr_coef.nc',inum)
+               CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array
+               CALL iom_close(inum)
+               !
+               tfrcoef2d(:,:) = rn_tfri2 * ( 1 + rn_tfrien * tfrcoef2d(:,:) )
+            ELSE
+               tfrcoef2d(:,:) = rn_tfri2  ! initialize tfrcoef2d to the namelist variable
+            ENDIF
+         END IF
          !
          IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all
-#  if defined key_vectopt_loop
-            DO jj = 1, 1
-!CDIR NOVERRCHK
-               DO ji = 1, jpij   ! vector opt. (forced unrolling)
-#  else
-!CDIR NOVERRCHK
             DO jj = 1, jpj
-!CDIR NOVERRCHK
                DO ji = 1, jpi
-#  endif
                   ikbt = mbkt(ji,jj)
                   ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp
@@ -281 +380 @@
                END DO
             END DO
+            IF ( ln_isfcav ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ikbt = mikt(ji,jj)
+                     ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_tfrz0 ))**2._wp
+                     tfrcoef2d(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp)
+                     tfrcoef2d(ji,jj) = MIN(tfrcoef2d(ji,jj), rn_tfri2_max)
+                  END DO
+               END DO
+            END IF
          ENDIF
          !
@@ -308 +417 @@
       zminbfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient
       zmaxbfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient
-      !
-#  if defined key_vectopt_loop
-      DO jj = 1, 1
-!CDIR NOVERRCHK
-         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-#  else
-!CDIR NOVERRCHK
+      zmintfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient
+      zmaxtfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient
+      !
       DO jj = 2, jpjm1
-!CDIR NOVERRCHK
          DO ji = 2, jpim1
-#  endif
              ikbu = mbku(ji,jj)       ! deepest ocean level at u- and v-points
              ikbv = mbkv(ji,jj)
@@ -339 +442 @@
              zminbfr = MIN(  zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) )  )
              zmaxbfr = MAX(  zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) )  )
+! (ISF)
+             IF ( ln_isfcav ) THEN
+                ikbu = miku(ji,jj)       ! 1st wet ocean level at u- and v-points
+                ikbv = mikv(ji,jj)
+                zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt
+                zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt
+                IF( ABS( tfrcoef2d(ji,jj) ) > zfru ) THEN
+                   IF( ln_ctl ) THEN
+                      WRITE(numout,*) 'TFR ', narea, nimpp+ji, njmpp+jj, ikbu
+                      WRITE(numout,*) 'TFR ', ABS( tfrcoef2d(ji,jj) ), zfru
+                   ENDIF
+                   ictu = ictu + 1
+                ENDIF
+                IF( ABS( tfrcoef2d(ji,jj) ) > zfrv ) THEN
+                   IF( ln_ctl ) THEN
+                      WRITE(numout,*) 'TFR ', narea, nimpp+ji, njmpp+jj, ikbv
+                      WRITE(numout,*) 'TFR ', tfrcoef2d(ji,jj), zfrv
+                   ENDIF
+                   ictv = ictv + 1
+                ENDIF
+                zmintfr = MIN(  zmintfr, MIN( zfru, ABS( tfrcoef2d(ji,jj) ) )  )
+                zmaxtfr = MAX(  zmaxtfr, MIN( zfrv, ABS( tfrcoef2d(ji,jj) ) )  )
+             END IF
+! END ISF
          END DO
       END DO
@@ -346 +473 @@
          CALL mpp_min( zminbfr )
          CALL mpp_max( zmaxbfr )
+         IF ( ln_isfcav) CALL mpp_min( zmintfr )
+         IF ( ln_isfcav) CALL mpp_max( zmaxtfr )
       ENDIF
       IF( .NOT.ln_bfrimp) THEN
       IF( lwp .AND. ictu + ictv > 0 ) THEN
-         WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points '
-         WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points '
+         WRITE(numout,*) ' Bottom/Top friction stability check failed at ', ictu, ' U-points '
+         WRITE(numout,*) ' Bottom/Top friction stability check failed at ', ictv, ' V-points '
          WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr
-         WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary'
+         IF ( ln_isfcav ) WRITE(numout,*) ' Top friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr
+         WRITE(numout,*) ' Bottom/Top friction coefficient will be reduced where necessary'
       ENDIF
       ENDIF

branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

@@ -6 +6 @@
    !! History :  OPA  ! 2000-08  (G. Madec)  double diffusive mixing
    !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
-   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
+   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
+   !!            3.6  ! 2013-04  (G. Madec, F. Roquet) zrau compute locally using interpolation of alpha & beta
    !!----------------------------------------------------------------------
 #if defined key_zdfddm   ||   defined key_esopa
@@ -18 +19 @@
    USE dom_oce         ! ocean space and time domain variables 
    USE zdf_oce         ! ocean vertical physics variables
+   USE eosbn2         ! equation of state
+   !
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -34 +37 @@
    LOGICAL , PUBLIC, PARAMETER ::   lk_zdfddm = .TRUE.  !: double diffusive mixing flag
 
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avs    !: salinity vertical diffusivity coeff. at w-point
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   rrau   !: heat/salt buoyancy flux ratio
-
-   !                      !!* Namelist namzdf_ddm : double diffusive mixing *
-   REAL(wp) ::   rn_avts   ! maximum value of avs for salt fingering
-   REAL(wp) ::   rn_hsbfr  ! heat/salt buoyancy flux ratio
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avs   !: salinity vertical diffusivity coeff. at w-point
+
+   !                       !!* Namelist namzdf_ddm : double diffusive mixing *
+   REAL(wp) ::   rn_avts    ! maximum value of avs for salt fingering
+   REAL(wp) ::   rn_hsbfr   ! heat/salt buoyancy flux ratio
 
    !! * Substitutions
+#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -54 +57 @@
       !!                ***  ROUTINE zdf_ddm_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( avs(jpi,jpj,jpk), rrau(jpi,jpj,jpk), STAT= zdf_ddm_alloc )
-      !
+      ALLOCATE( avs(jpi,jpj,jpk) , STAT= zdf_ddm_alloc )
       IF( lk_mpp             )   CALL mpp_sum ( zdf_ddm_alloc )
       IF( zdf_ddm_alloc /= 0 )   CALL ctl_warn('zdf_ddm_alloc: failed to allocate arrays')
@@ -71 +73 @@
       !!      diffusive mixing (i.e. salt fingering and diffusive layering)
       !!      following Merryfield et al. (1999). The rate of double diffusive 
-      !!      mixing depend on the buoyancy ratio: Rrau=alpha/beta dk[T]/dk[S]
-      !!      which is computed in rn2.F
+      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]):
       !!         * salt fingering (Schmitt 1981):
-      !!      for Rrau > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (Rrau/rn_hsbfr)^6 )
-      !!      for Rrau > 1 and rn2 > 0 : zavfs = O
-      !!      otherwise                : zavft = 0.7 zavs / Rrau
+      !!      for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 )
+      !!      for R > 1 and rn2 > 0 : zavfs = O
+      !!      otherwise                : zavft = 0.7 zavs / R
       !!         * diffusive layering (Federov 1988):
-      !!      for 0< Rrau < 1 and rn2 > 0 : zavdt = 1.3635e-6  
-      !!                                 * exp( 4.6 exp(-0.54 (1/Rrau-1) ) )
+      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) )
       !!      otherwise                   : zavdt = 0 
-      !!      for .5 < Rrau < 1 and rn2 > 0 : zavds = zavdt (1.885 Rrau -0.85)
-      !!      for  0 < Rrau <.5 and rn2 > 0 : zavds = zavdt 0.15 Rrau      
+      !!      for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85)
+      !!      for  0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R      
       !!      otherwise                     : zavds = 0 
       !!         * update the eddy diffusivity:
@@ -96 +96 @@
       !
       INTEGER  ::   ji, jj , jk     ! dummy loop indices
-      REAL(wp) ::   zinr, zrr       ! temporary scalars
-      REAL(wp) ::   zavft, zavfs    !    -         -
-      REAL(wp) ::   zavdt, zavds    !    -         -
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zmsks, zmskf, zmskd1, zmskd2, zmskd3
+      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
+      REAL(wp) ::   zdt, zds
+      REAL(wp) ::   zinr, zrr       !   -      -
+      REAL(wp) ::   zavft, zavfs    !   -      -
+      REAL(wp) ::   zavdt, zavds    !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('zdf_ddm')
       !
-      CALL wrk_alloc( jpi,jpj, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
-
+      CALL wrk_alloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
+      !
       !                                                ! ===============
       DO jk = 2, jpkm1                                 ! Horizontal slab
@@ -111 +113 @@
          ! Define the mask 
          ! ---------------
-         rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) )         ! only retains positive value of rrau
+         DO jj = 1, jpj                                ! R=zrau = (alpha / beta) (dk[t] / dk[s])
+            DO ji = 1, jpi
+               zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   &
+                  &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 
+               !
+               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
+                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
+               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
+                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
+               !
+               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )
+               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 
+               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
+               zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
+            END DO
+         END DO
 
          DO jj = 1, jpj                                     ! indicators:
@@ -119 +136 @@
                ELSE                                       ;   zmsks(ji,jj) = 1._wp
                ENDIF
-               ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere            
-               IF( rrau(ji,jj,jk) <= 1.          ) THEN   ;   zmskf(ji,jj) = 0._wp
+               ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere            
+               IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
                ELSE                                       ;   zmskf(ji,jj) = 1._wp
                ENDIF
                ! diffusive layering indicators: 
-               !     ! mskdl1=1 if 0<rrau<1; 0 elsewhere
-               IF( rrau(ji,jj,jk) >= 1.          ) THEN   ;   zmskd1(ji,jj) = 0._wp
+               !     ! mskdl1=1 if 0< R <1; 0 elsewhere
+               IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
                ELSE                                       ;   zmskd1(ji,jj) = 1._wp
                ENDIF
-               !     ! mskdl2=1 if 0<rrau<0.5; 0 elsewhere
-               IF( rrau(ji,jj,jk) >= 0.5         ) THEN   ;   zmskd2(ji,jj) = 0._wp
+               !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
+               IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
                ELSE                                       ;   zmskd2(ji,jj) = 1._wp
                ENDIF
-               !   mskdl3=1 if 0.5<rrau<1; 0 elsewhere
-               IF( rrau(ji,jj,jk) <= 0.5 .OR. rrau(ji,jj,jk) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
-               ELSE                                                         ;   zmskd3(ji,jj) = 1._wp
+               !   mskdl3=1 if 0.5< R <1; 0 elsewhere
+               IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
+               ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
                ENDIF
             END DO
          END DO
          ! mask zmsk in order to have avt and avs masked
-         zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)
+         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
 
 
@@ -149 +166 @@
 !CDIR NOVERRCHK
             DO ji = 1, jpi
-               zinr = 1./rrau(ji,jj,jk)
+               zinr = 1._wp / zrau(ji,jj)
                ! salt fingering
-               zrr = rrau(ji,jj,jk)/rn_hsbfr
+               zrr = zrau(ji,jj) / rn_hsbfr
                zrr = zrr * zrr
                zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
@@ -157 +174 @@
                ! diffusive layering
                zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
-               zavds = zavdt * zmsks(ji,jj) * (  (1.85 * rrau(ji,jj,jk) - 0.85 ) * zmskd3(ji,jj)   &
-                  &                            +  0.15 * rrau(ji,jj,jk)          * zmskd2(ji,jj)  )
+               zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
+                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
                ! add to the eddy viscosity coef. previously computed
                avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
@@ -174 +191 @@
                avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    &
                   &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   &
-                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * umask(ji,jj,jk)
+                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * wumask(ji,jj,jk)
                avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    &
                   &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   &
-                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * vmask(ji,jj,jk)
+                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * wvmask(ji,jj,jk)
             END DO
          END DO
@@ -196 +213 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi,jpj, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
+      CALL wrk_dealloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
       !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_ddm')
@@ -212 +229 @@
       !!              called by zdf_ddm at the first timestep (nit000)
       !!----------------------------------------------------------------------
+      INTEGER ::   ios   ! local integer
+      !!
       NAMELIST/namzdf_ddm/ rn_avts, rn_hsbfr
-      INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!----------------------------------------------------------------------
       !
@@ -237 +255 @@
       IF( zdf_ddm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' )
       !                               ! initialization to masked Kz
-      avs(:,:,:) = rn_avt0 * tmask(:,:,:) 
+      avs(:,:,:) = rn_avt0 * wmask(:,:,:) 
       !
    END SUBROUTINE zdf_ddm_init

branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90

@@ -78 +78 @@
          !
          DO jk = 1, jpkm1 
-#if defined key_vectopt_loop
-            DO jj = 1, 1                     ! big loop forced
-               DO ji = jpi+2, jpij   
-#else
             DO jj = 2, jpj             ! no vector opt.
                DO ji = 2, jpi
-#endif
 #if defined key_zdfkpp
                   ! no evd mixing in the boundary layer with KPP
@@ -110 +105 @@
          DO jk = 1, jpkm1
 !!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL! 
-#if defined key_vectopt_loop
-            DO jj = 1, 1                     ! big loop forced
-               DO ji = 1, jpij   
-#else
             DO jj = 1, jpj             ! loop over the whole domain (no lbc_lnk call)
                DO ji = 1, jpi
-#endif
 #if defined key_zdfkpp
                   ! no evd mixing in the boundary layer with KPP

branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

@@ -20 +20 @@
    USE domvvl         ! ocean space and time domain : variable volume layer
    USE zdf_oce        ! ocean vertical physics
+   USE zdfbfr         ! bottom friction (only for rn_bfrz0)
    USE sbc_oce        ! surface boundary condition: ocean
    USE phycst         ! physical constants
@@ -52 +53 @@
 
    !                              !! ** Namelist  namzdf_gls  **
-   LOGICAL  ::   ln_crban          ! =T use Craig and Banner scheme
    LOGICAL  ::   ln_length_lim     ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988)
    LOGICAL  ::   ln_sigpsi         ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing
-   INTEGER  ::   nn_tkebc_surf     ! TKE surface boundary condition (=0/1)
-   INTEGER  ::   nn_tkebc_bot      ! TKE bottom boundary condition (=0/1)
-   INTEGER  ::   nn_psibc_surf     ! PSI surface boundary condition (=0/1)
-   INTEGER  ::   nn_psibc_bot      ! PSI bottom boundary condition (=0/1)
+   INTEGER  ::   nn_bc_surf        ! surface boundary condition (=0/1)
+   INTEGER  ::   nn_bc_bot         ! bottom boundary condition (=0/1)
+   INTEGER  ::   nn_z0_met         ! Method for surface roughness computation
    INTEGER  ::   nn_stab_func      ! stability functions G88, KC or Canuto (=0/1/2)
    INTEGER  ::   nn_clos           ! closure 0/1/2/3 MY82/k-eps/k-w/gen
@@ -66 +65 @@
    REAL(wp) ::   rn_charn          ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value)
    REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing
-
-   REAL(wp) ::   hsro          =  0.003_wp    ! Minimum surface roughness
-   REAL(wp) ::   hbro          =  0.003_wp    ! Bottom roughness (m)
+   REAL(wp) ::   rn_hsro           ! Minimum surface roughness
+   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 
+
    REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters
    REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1 
@@ -96 +95 @@
    REAL(wp) ::   rm7           =  0.0_wp
    REAL(wp) ::   rm8           =  0.318_wp
-   
+   REAL(wp) ::   rtrans        =  0.1_wp
    REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters
-   REAL(wp) ::   rc03_sqrt2_galp                                  !     -           -           -        -
-   REAL(wp) ::   rsbc_tke1, rsbc_tke2, rsbc_tke3, rfact_tke       !     -           -           -        -
-   REAL(wp) ::   rsbc_psi1, rsbc_psi2, rsbc_psi3, rfact_psi       !     -           -           -        -
-   REAL(wp) ::   rsbc_mb  , rsbc_std , rsbc_zs                    !     -           -           -        -
+   REAL(wp) ::   rsbc_tke1, rsbc_tke2, rfact_tke                  !     -           -           -        -
+   REAL(wp) ::   rsbc_psi1, rsbc_psi2, rfact_psi                  !     -           -           -        -
+   REAL(wp) ::   rsbc_zs1, rsbc_zs2                               !     -           -           -        -
    REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        -
    REAL(wp) ::   rs0, rs1, rs2, rs4, rs5, rs6                     !     -           -           -        -
@@ -147 +145 @@
       REAL(wp) ::   gh, gm, shr, dif, zsqen, zav        !   -      -
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zdep
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zkar
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves 
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhsro       ! Surface roughness (surface waves)
@@ -153 +152 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   shear       ! vertical shear
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   eps         ! dissipation rate
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi.AND.ln_crban=T)
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_a, z_elem_b, z_elem_c, psi
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi)
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   psi         ! psi at time now
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_a    ! element of the first  matrix diagonal
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_b    ! element of the second matrix diagonal
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_c    ! element of the third  matrix diagonal
       !!--------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('zdf_gls')
       !
-      CALL wrk_alloc( jpi,jpj, zdep, zflxs, zhsro )
-      CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )
-
+      CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro )
+      CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi  )
+      
       ! Preliminary computing
 
@@ -174 +176 @@
 
       ! Compute surface and bottom friction at T-points
-!CDIR NOVERRCHK
-      DO jj = 2, jpjm1
-!CDIR NOVERRCHK
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            ! 
-            ! surface friction 
+!CDIR NOVERRCHK          
+      DO jj = 2, jpjm1          
+!CDIR NOVERRCHK         
+         DO ji = fs_2, fs_jpim1   ! vector opt.         
+            !
+            ! surface friction
             ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)
-            !
-            ! bottom friction (explicit before friction)
-            ! Note that we chose here not to bound the friction as in dynbfr)
-            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &
-               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )
-            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &
-               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )
-            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)
-         END DO
-      END DO  
-
-      ! In case of breaking surface waves mixing,
-      ! Compute surface roughness length according to Charnock formula:
-      IF( ln_crban ) THEN   ;   zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro)
-      ELSE                  ;   zhsro(:,:) = hsro
-      ENDIF
+            !   
+            ! bottom friction (explicit before friction)        
+            ! Note that we chose here not to bound the friction as in dynbfr)   
+            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &         
+               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )      
+            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &         
+               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )      
+            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)         
+         END DO         
+      END DO    
+
+      ! Set surface roughness length
+      SELECT CASE ( nn_z0_met )
+      !
+      CASE ( 0 )             ! Constant roughness          
+         zhsro(:,:) = rn_hsro
+      CASE ( 1 )             ! Standard Charnock formula
+         zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro)
+      CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008)
+         zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10)
+         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11)
+      !
+      END SELECT
 
       ! Compute shear and dissipation rate
@@ -303 +312 @@
       !
       ! Set surface condition on zwall_psi (1 at the bottom)
-      IF( ln_sigpsi ) THEN
-         zcoef = rsc_psi / rsc_psi0
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zwall_psi(ji,jj,1) = zcoef
-            END DO
-         END DO
-      ENDIF
-
+      zwall_psi(:,:,1) = zwall_psi(:,:,2)
+      zwall_psi(:,:,jpk) = 1.
+      !
       ! Surface boundary condition on tke
       ! ---------------------------------
       !
-      SELECT CASE ( nn_tkebc_surf )
+      SELECT CASE ( nn_bc_surf )
       !
       CASE ( 0 )             ! Dirichlet case
-         !
-         IF (ln_crban) THEN     ! Wave induced mixing case
-            !                      ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2
-            !                      ! balance between the production and the dissipation terms including the wave effect
-            en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin )
-            z_elem_a(:,:,1) = en(:,:,1)
-            z_elem_c(:,:,1) = 0._wp
-            z_elem_b(:,:,1) = 1._wp
-            ! 
-            ! one level below
-            en(:,:,2) = MAX( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin )
-            z_elem_a(:,:,2) = 0._wp
-            z_elem_c(:,:,2) = 0._wp
-            z_elem_b(:,:,2) = 1._wp
-            !
-         ELSE                   ! No wave induced mixing case
-            !                      ! en(1) = u*^2/C0^2  &  l(1)  = K*zs
-            !                      ! balance between the production and the dissipation terms
-            en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin )
-            z_elem_a(:,:,1) = en(:,:,1) 
-            z_elem_c(:,:,1) = 0._wp
-            z_elem_b(:,:,1) = 1._wp
-            !
-            ! one level below
-            en(:,:,2) = MAX( rc02r * ustars2(:,:), rn_emin )
-            z_elem_a(:,:,2) = 0._wp
-            z_elem_c(:,:,2) = 0._wp
-            z_elem_b(:,:,2) = 1._wp
-            !
-         ENDIF
-         !
+      ! First level
+      en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
+      en(:,:,1) = MAX(en(:,:,1), rn_emin) 
+      z_elem_a(:,:,1) = en(:,:,1)
+      z_elem_c(:,:,1) = 0._wp
+      z_elem_b(:,:,1) = 1._wp
+      ! 
+      ! One level below
+      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
+      en(:,:,2) = MAX(en(:,:,2), rn_emin )
+      z_elem_a(:,:,2) = 0._wp 
+      z_elem_c(:,:,2) = 0._wp
+      z_elem_b(:,:,2) = 1._wp
+      !
+      !
       CASE ( 1 )             ! Neumann boundary condition on d(e)/dz
-         !
-         IF (ln_crban) THEN ! Shear free case: d(e)/dz= Fw
-            !
-            ! Dirichlet conditions at k=1 (Cosmetic)
-            en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin )
-            z_elem_a(:,:,1) = en(:,:,1)
-            z_elem_c(:,:,1) = 0._wp
-            z_elem_b(:,:,1) = 1._wp
-            ! at k=2, set de/dz=Fw
-            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
-            z_elem_a(:,:,2) = 0._wp        
-            zflxs(:,:) = rsbc_tke3 * ustars2(:,:)**1.5_wp * ( (zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf)
-            en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
-            !
-         ELSE                   ! No wave induced mixing case: d(e)/dz=0.
-            !
-            ! Dirichlet conditions at k=1 (Cosmetic)
-            en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin )
-            z_elem_a(:,:,1) = en(:,:,1)
-            z_elem_c(:,:,1) = 0._wp
-            z_elem_b(:,:,1) = 1._wp
-            ! at k=2 set de/dz=0.:
-            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2)  ! Remove z_elem_a from z_elem_b
-            z_elem_a(:,:,2) = 0._wp
-            !
-         ENDIF
-         !
+      !
+      ! Dirichlet conditions at k=1
+      en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
+      en(:,:,1)       = MAX(en(:,:,1), rn_emin)      
+      z_elem_a(:,:,1) = en(:,:,1)
+      z_elem_c(:,:,1) = 0._wp
+      z_elem_b(:,:,1) = 1._wp
+      !
+      ! at k=2, set de/dz=Fw
+      !cbr
+      z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
+      z_elem_a(:,:,2) = 0._wp
+      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) ))
+      zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf)
+
+      en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2)
+      !
+      !
       END SELECT
 
@@ -382 +360 @@
       ! --------------------------------
       !
-      SELECT CASE ( nn_tkebc_bot )
+      SELECT CASE ( nn_bc_bot )
       !
       CASE ( 0 )             ! Dirichlet 
@@ -457 +435 @@
       !                                            ! set the minimum value of tke 
       en(:,:,:) = MAX( en(:,:,:), rn_emin )
-      
+
       !!----------------------------------------!!
       !!   Solve prognostic equation for psi    !!
@@ -560 +538 @@
       ! ---------------------------------
       !
-      SELECT CASE ( nn_psibc_surf )
+      SELECT CASE ( nn_bc_surf )
       !
       CASE ( 0 )             ! Dirichlet boundary conditions
-         !
-         IF( ln_crban ) THEN       ! Wave induced mixing case
-            !                      ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2
-            !                      ! balance between the production and the dissipation terms including the wave effect
-            zdep(:,:) = rl_sf * zhsro(:,:)
-            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-            z_elem_a(:,:,1) = psi(:,:,1)
-            z_elem_c(:,:,1) = 0._wp
-            z_elem_b(:,:,1) = 1._wp
-            !
-            ! one level below
-            zex1 = (rmm*ra_sf+rnn)
-            zex2 = (rmm*ra_sf)
-            zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2
-            psi (:,:,2) = rsbc_psi1 * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1)
-            z_elem_a(:,:,2) = 0._wp
-            z_elem_c(:,:,2) = 0._wp
-            z_elem_b(:,:,2) = 1._wp
-            ! 
-         ELSE                   ! No wave induced mixing case
-            !                      ! en(1) = u*^2/C0^2  &  l(1)  = K*zs
-            !                      ! balance between the production and the dissipation terms
-            !
-            zdep(:,:) = vkarmn * zhsro(:,:)
-            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-            z_elem_a(:,:,1) = psi(:,:,1)
-            z_elem_c(:,:,1) = 0._wp
-            z_elem_b(:,:,1) = 1._wp
-            !
-            ! one level below
-            zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) )
-            psi (:,:,2) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-            z_elem_a(:,:,2) = 0._wp
-            z_elem_c(:,:,2) = 0._wp
-            z_elem_b(:,:,2) = 1.
-            !
-         ENDIF
-         !
+      !
+      ! Surface value
+      zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic
+      psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+      z_elem_a(:,:,1) = psi(:,:,1)
+      z_elem_c(:,:,1) = 0._wp
+      z_elem_b(:,:,1) = 1._wp
+      !
+      ! One level below
+      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) )))
+      zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:)
+      psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+      z_elem_a(:,:,2) = 0._wp
+      z_elem_c(:,:,2) = 0._wp
+      z_elem_b(:,:,2) = 1._wp
+      ! 
+      !
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
-         !
-         IF( ln_crban ) THEN     ! Wave induced mixing case
-            !
-            zdep(:,:) = rl_sf * zhsro(:,:)
-            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-            z_elem_a(:,:,1) = psi(:,:,1)
-            z_elem_c(:,:,1) = 0._wp
-            z_elem_b(:,:,1) = 1._wp
-            !
-            ! Neumann condition at k=2
-            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
-            z_elem_a(:,:,2) = 0._wp
-            !
-            ! Set psi vertical flux at the surface:
-            zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf)
-            zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 
-               &                   * en(:,:,1)**rmm * zdep         
-            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
-            !
-      ELSE                   ! No wave induced mixing
-            !
-            zdep(:,:) = vkarmn * zhsro(:,:)
-            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-            z_elem_a(:,:,1) = psi(:,:,1)
-            z_elem_c(:,:,1) = 0._wp
-            z_elem_b(:,:,1) = 1._wp
-            !
-            ! Neumann condition at k=2
-            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
-            z_elem_a(ji,jj,2) = 0._wp
-            !
-            ! Set psi vertical flux at the surface:
-            zdep(:,:)  = zhsro(:,:) + fsdept(:,:,1)
-            zflxs(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-1._wp)
-            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
-            !     
-         ENDIF
-         !
+      !
+      ! Surface value: Dirichlet
+      zdep(:,:)       = zhsro(:,:) * rl_sf
+      psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+      z_elem_a(:,:,1) = psi(:,:,1)
+      z_elem_c(:,:,1) = 0._wp
+      z_elem_b(:,:,1) = 1._wp
+      !
+      ! Neumann condition at k=2
+      z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
+      z_elem_a(:,:,2) = 0._wp
+      !
+      ! Set psi vertical flux at the surface:
+      zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
+      zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)
+      zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
+      zdep(:,:) =  rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &
+             & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.)
+      zflxs(:,:) = zdep(:,:) * zflxs(:,:)
+      psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
+
+      !   
+      !
       END SELECT
 
@@ -645 +587 @@
       ! --------------------------------
       !
-      SELECT CASE ( nn_psibc_bot )
+      SELECT CASE ( nn_bc_bot )
+      !
       !
       CASE ( 0 )             ! Dirichlet 
-         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * hbro
+         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0
          !                      ! Balance between the production and the dissipation terms
 !CDIR NOVERRCHK
@@ -656 +599 @@
                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
-               zdep(ji,jj) = vkarmn * hbro
+               zdep(ji,jj) = vkarmn * rn_bfrz0
                psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
                z_elem_a(ji,jj,ibot) = 0._wp
@@ -663 +606 @@
                !
                ! Just above last level, Dirichlet condition again (GOTM like)
-               zdep(ji,jj) = vkarmn * ( hbro + fse3t(ji,jj,ibotm1) )
+               zdep(ji,jj) = vkarmn * ( rn_bfrz0 + fse3t(ji,jj,ibotm1) )
                psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn
                z_elem_a(ji,jj,ibotm1) = 0._wp
@@ -681 +624 @@
                !
                ! Bottom level Dirichlet condition:
-               zdep(ji,jj) = vkarmn * hbro
+               zdep(ji,jj) = vkarmn * rn_bfrz0
                psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
                !
@@ -693 +636 @@
                !
                ! Set psi vertical flux at the bottom:
-               zdep(ji,jj) = hbro + 0.5_wp*fse3t(ji,jj,ibotm1)
+               zdep(ji,jj) = rn_bfrz0 + 0.5_wp*fse3t(ji,jj,ibotm1)
                zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) )   &
                   &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
@@ -736 +679 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / psi(ji,jj,jk)
+                  eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin)
                END DO
             END DO
@@ -783 +726 @@
                ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
-               mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk)  )
+               IF (ln_length_lim) mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) )
             END DO
          END DO
@@ -847 +790 @@
       ! Boundary conditions on stability functions for momentum (Neumann):
       ! Lines below are useless if GOTM style Dirichlet conditions are used
-      zcoef = rcm_sf / SQRT( 2._wp )
+
+      avmv(:,:,1) = avmv(:,:,2)
+
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            avmv(ji,jj,1) = zcoef
-         END DO
-      END DO
-      zcoef = rc0 / SQRT( 2._wp )
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            avmv(ji,jj,mbkt(ji,jj)+1) = zcoef
+            avmv(ji,jj,mbkt(ji,jj)+1) = avmv(ji,jj,mbkt(ji,jj))
          END DO
       END DO
@@ -900 +839 @@
       avmv_k(:,:,:) = avmv(:,:,:)
       !
-      CALL wrk_dealloc( jpi,jpj, zdep, zflxs, zhsro )
+      CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro )
       CALL wrk_dealloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )
       !
@@ -932 +871 @@
       !!
       NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &
-         &            rn_clim_galp, ln_crban, ln_sigpsi,     &
-         &            rn_crban, rn_charn,                    &
-         &            nn_tkebc_surf, nn_tkebc_bot,           &
-         &            nn_psibc_surf, nn_psibc_bot,           &
+         &            rn_clim_galp, ln_sigpsi, rn_hsro,      &
+         &            rn_crban, rn_charn, rn_frac_hs,        &
+         &            nn_bc_surf, nn_bc_bot, nn_z0_met,      &
          &            nn_stab_func, nn_clos
       !!----------------------------------------------------------
@@ -955 +893 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namzdf_gls : set gls mixing parameters'
-         WRITE(numout,*) '      minimum value of en                           rn_emin       = ', rn_emin
-         WRITE(numout,*) '      minimum value of eps                          rn_epsmin     = ', rn_epsmin
-         WRITE(numout,*) '      Limit dissipation rate under stable stratif.  ln_length_lim = ', ln_length_lim
-         WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp  = ', rn_clim_galp
-         WRITE(numout,*) '      TKE Surface boundary condition                nn_tkebc_surf = ', nn_tkebc_surf
-         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_tkebc_bot  = ', nn_tkebc_bot
-         WRITE(numout,*) '      PSI Surface boundary condition                nn_psibc_surf = ', nn_psibc_surf
-         WRITE(numout,*) '      PSI Bottom boundary condition                 nn_psibc_bot  = ', nn_psibc_bot
-         WRITE(numout,*) '      Craig and Banner scheme                       ln_crban      = ', ln_crban
-         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi     = ', ln_sigpsi
+         WRITE(numout,*) '      minimum value of en                           rn_emin        = ', rn_emin
+         WRITE(numout,*) '      minimum value of eps                          rn_epsmin      = ', rn_epsmin
+         WRITE(numout,*) '      Limit dissipation rate under stable stratif.  ln_length_lim  = ', ln_length_lim
+         WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp   = ', rn_clim_galp
+         WRITE(numout,*) '      TKE Surface boundary condition                nn_bc_surf     = ', nn_bc_surf
+         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot
+         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi
          WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban
          WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn
+         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met
+         WRITE(numout,*) '      Wave height frac. (used if nn_z0_met=2)       rn_frac_hs     = ', rn_frac_hs
          WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func
          WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos
-         WRITE(numout,*) '   Hard coded parameters'
-         WRITE(numout,*) '      Surface roughness (m)                         hsro          = ', hsro
-         WRITE(numout,*) '      Bottom roughness (m)                          hbro          = ', hbro
+         WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro
+         WRITE(numout,*) '      Bottom roughness (m) (nambfr namelist)        rn_bfrz0       = ', rn_bfrz0
       ENDIF
 
@@ -978 +914 @@
 
       !                                !* Check of some namelist values
-      IF( nn_tkebc_surf < 0 .OR. nn_tkebc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_surf is 0 or 1' )
-      IF( nn_psibc_surf < 0 .OR. nn_psibc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_surf is 0 or 1' )
-      IF( nn_tkebc_bot  < 0 .OR. nn_tkebc_bot  > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' )
-      IF( nn_psibc_bot  < 0 .OR. nn_psibc_bot  > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' )
+      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' )
+      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' )
+      IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' )
       IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' )
       IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' )
@@ -1001 +936 @@
          SELECT CASE ( nn_stab_func )
          CASE( 0, 1 )   ;   rpsi3m = 2.53_wp       ! G88 or KC stability functions
-         CASE( 2 )      ;   rpsi3m = 2.38_wp       ! Canuto A stability functions
+         CASE( 2 )      ;   rpsi3m = 2.62_wp       ! Canuto A stability functions
          CASE( 3 )      ;   rpsi3m = 2.38          ! Canuto B stability functions (caution : constant not identified)
          END SELECT
@@ -1012 +947 @@
          rnn     = -1._wp
          rsc_tke =  1._wp
-         rsc_psi =  1.3_wp  ! Schmidt number for psi
+         rsc_psi =  1.2_wp  ! Schmidt number for psi
          rpsi1   =  1.44_wp
          rpsi3p  =  1._wp
@@ -1140 +1075 @@
       !                                     ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
       !                                     !  or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
-      IF( ln_sigpsi .AND. ln_crban ) THEN
-         zcr = SQRT( 1.5_wp*rsc_tke ) * rcm_sf / vkarmn
-         rsc_psi0 = vkarmn*vkarmn / ( rpsi2 * rcm_sf*rcm_sf )                       & 
-        &         * ( rnn*rnn - 4._wp/3._wp * zcr*rnn*rmm - 1._wp/3._wp * zcr*rnn   &
-        &           + 2._wp/9._wp * rmm * zcr*zcr + 4._wp/9._wp * zcr*zcr * rmm*rmm )                                 
+      IF( ln_sigpsi ) THEN
+         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 
+         ! Verification: retrieve Burchard (2001) results by uncomenting the line below:
+         ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work
+         ! ra_sf = -SQRT(2./3.*rc0**3./rn_cm_sf*rn_sc_tke)/vkarmn
+         rsc_psi0 = rsc_tke/(24.*rpsi2)*(-1.+(4.*rnn + ra_sf*(1.+4.*rmm))**2./(ra_sf**2.))
       ELSE
          rsc_psi0 = rsc_psi
@@ -1151 +1087 @@
       !                                !* Shear free turbulence parameters
       !
-      ra_sf  = -4._wp * rnn * SQRT( rsc_tke ) / ( (1._wp+4._wp*rmm) * SQRT( rsc_tke )   &
-         &                                      - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
-      rl_sf  = rc0 * SQRT( rc0 / rcm_sf )                                                                   &
-         &         * SQRT(  (  (1._wp + 4._wp*rmm + 8._wp*rmm*rmm) * rsc_tke                                &
-         &                   + 12._wp * rsc_psi0 * rpsi2                                                    &
-         &                   - (1._wp + 4._wp*rmm) * SQRT( rsc_tke*(rsc_tke+ 24._wp*rsc_psi0*rpsi2) )  )    &
-         &                / ( 12._wp*rnn*rnn )                                                              )
+      ra_sf  = -4._wp*rnn*SQRT(rsc_tke) / ( (1._wp+4._wp*rmm)*SQRT(rsc_tke) &
+               &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
+
+      IF ( rn_crban==0._wp ) THEN
+         rl_sf = vkarmn
+      ELSE
+         rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp)*rsc_tke          &
+                 &                                       + 12._wp * rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm) &
+                 &                                                *SQRT(rsc_tke*(rsc_tke                 &
+                 &                                                   + 24._wp*rsc_psi0*rpsi2)) )         &
+                 &                                         /(12._wp*rnn**2.)                             &
+                 &                                       )
+      ENDIF
 
       !
@@ -1187 +1129 @@
       rc03  = rc02 * rc0
       rc04  = rc03 * rc0
-      rc03_sqrt2_galp = rc03 / SQRT(2._wp) / rn_clim_galp
-      rsbc_mb   = 0.5_wp * (15.8_wp*rn_crban)**(2._wp/3._wp)               ! Surf. bound. cond. from Mellor and Blumberg
-      rsbc_std  = 3.75_wp                                                  ! Surf. bound. cond. standard (prod=diss)
-      rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp)  ! k_eps = 53.  Dirichlet + Wave breaking 
-      rsbc_tke2 = 0.5_wp / rau0
-      rsbc_tke3 = rdt * rn_crban                                                         ! Neumann + Wave breaking
-      rsbc_zs   = rn_charn / grav                                                        ! Charnock formula
-      rsbc_psi1 = rc0**rpp * rsbc_tke1**rmm * rl_sf**rnn                           ! Dirichlet + Wave breaking
-      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi                   ! Neumann + NO Wave breaking 
-      rsbc_psi3 = -0.5_wp * rdt * rc0**rpp * rl_sf**rnn / rsc_psi  * (rnn + rmm*ra_sf) ! Neumann + Wave breaking
-      rfact_tke = -0.5_wp / rsc_tke * rdt               ! Cst used for the Diffusion term of tke
-      rfact_psi = -0.5_wp / rsc_psi * rdt               ! Cst used for the Diffusion term of tke
+      rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking
+      rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking 
+      zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )
+      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer 
+      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness
+      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness 
+      rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi
+      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 
+
+      rfact_tke = -0.5_wp / rsc_tke * rdt                                ! Cst used for the Diffusion term of tke
+      rfact_psi = -0.5_wp / rsc_psi * rdt                                ! Cst used for the Diffusion term of tke
 
       !                                !* Wall proximity function
@@ -1257 +1198 @@
                IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop'
                en  (:,:,:) = rn_emin
-               mxln(:,:,:) = 0.001        
+               mxln(:,:,:) = 0.05        
+               avt_k (:,:,:) = avt (:,:,:)
+               avm_k (:,:,:) = avm (:,:,:)
+               avmu_k(:,:,:) = avmu(:,:,:)
+               avmv_k(:,:,:) = avmv(:,:,:)
                DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO
             ENDIF
@@ -1263 +1208 @@
             IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'
             en  (:,:,:) = rn_emin
-            mxln(:,:,:) = 0.001       
+            mxln(:,:,:) = 0.05       
          ENDIF
          !
@@ -1269 +1214 @@
          !                                   ! -------------------
          IF(lwp) WRITE(numout,*) '---- gls-rst ----'
-         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
+         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     ) 
          CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
          CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
-         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
+         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 
          CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
          CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln   )

branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90

@@ -117 +117 @@
       IF( ioptio == 0 .OR. ioptio > 1 .AND. .NOT. lk_esopa )   &
          &   CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' )
+      IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. ln_isfcav )   &
+         &   CALL ctl_stop( ' only zdfcst and zdftke were tested with ice shelves cavities ' )
       !
       !                               ! ... Convection
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) '   convection :'
+      !
+#if defined key_top
+      IF( ln_zdfnpc )   CALL ctl_stop( ' zdf_init: npc scheme is not working with key_top' )
+#endif
+      !
       ioptio = 0
       IF( ln_zdfnpc ) THEN

branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90

@@ -26 +26 @@
    USE phycst         ! physical constants
    USE eosbn2         ! equation of state
-   USE zdfddm         ! double diffusion mixing
+   USE zdfddm         ! double diffusion mixing (avs array)
+   USE lib_mpp        ! MPP library
+   USE trd_oce        ! ocean trends definition
+   USE trdtra         ! tracers trends
+   !
    USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! work arrays
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE prtctl         ! Print control
-   USE trdmod_oce     ! ocean trends definition
-   USE trdtra         ! tracers trends
+   USE wrk_nemo       ! work arrays
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
 
    IMPLICIT NONE
@@ -246 +247 @@
       REAL(wp) ::   zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t
 #if defined key_zdfddm
-      REAL(wp) ::   zrrau, zds, zavdds, zavddt,zinr   ! double diffusion mixing
-      REAL(wp), POINTER, DIMENSION(:,:)   ::     zdifs
+      REAL(wp) ::   zrw, zkm1s                    ! local scalars
+      REAL(wp) ::   zrrau, zdt, zds, zavdds, zavddt, zinr   ! double diffusion mixing
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdifs
       REAL(wp), POINTER, DIMENSION(:)     ::   za2s, za3s, zkmps
-      REAL(wp) ::                            zkm1s
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zblcs
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdiffus
@@ -274 +275 @@
 #endif
 
-      zviscos(:,:,:) = 0.
-      zblcm  (:,:  ) = 0. 
-      zdiffut(:,:,:) = 0.
-      zblct  (:,:  ) = 0. 
+      zviscos(:,:,:) = 0._wp
+      zblcm  (:,:  ) = 0._wp
+      zdiffut(:,:,:) = 0._wp
+      zblct  (:,:  ) = 0._wp 
 #if defined key_zdfddm
-      zdiffus(:,:,:) = 0.
-      zblcs  (:,:  ) = 0. 
-#endif
-      ghats(:,:,:) = 0.
-     
-      zBo   (:,:) = 0.
-      zBosol(:,:) = 0.
-      zustar(:,:) = 0.
-
-
+      zdiffus(:,:,:) = 0._wp
+      zblcs  (:,:  ) = 0._wp 
+#endif
+      ghats  (:,:,:) = 0._wp
+      zBo    (:,:  ) = 0._wp
+      zBosol (:,:  ) = 0._wp
+      zustar (:,:  ) = 0._wp
+      !
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! I. Interior diffusivity and viscosity at w points ( T interfaces)
@@ -332 +331 @@
                   avt (ji,jj,jk) =  avt (ji,jj,jk) + rn_difri * zfri    
                ENDIF
+               !
 #if defined key_zdfddm 
-               avs (ji,jj,jk) =  avt (ji,jj,jk)              
+               !
                !  Double diffusion mixing ; NOT IN ROUTINE ZDFDDM.F90
-               ! ------------------------------------------------------------------
-               ! only retains positive value of rrau
-               zrrau = MAX( rrau(ji,jj,jk), epsln )
-               zds   = tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal)
-               IF( zrrau > 1. .AND. zds > 0.) THEN
-                  !
-                  ! Salt fingering case.
-                  !---------------------
-                  ! Compute interior diffusivity for double diffusive mixing of
-                  ! salinity. Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001).
-                  ! After that set interior diffusivity for double diffusive mixing
-                  ! of temperature
+               ! -------------------------
+               avs (ji,jj,jk) =  avt (ji,jj,jk)   
+
+               ! R=zrau = (alpha / beta) (dk[t] / dk[s])
+               zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   &
+                  &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 
+               !
+               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  ) * tmask(ji,jj,jk)
+               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  ) * tmask(ji,jj,jk)
+               !
+               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )
+               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 
+               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
+               zrrau = MAX(  epsln , zdt / zds  )    ! only retains positive value of zrau
+               !
+               IF( zrrau > 1. .AND. zds > 0.) THEN                        ! Salt fingering case.
+                  !                                                       !---------------------
+                  ! Compute interior diffusivity for double diffusive mixing of salinity. 
+                  ! Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001).
+                  ! After that set interior diffusivity for double diffusive mixing of temperature
                   zavdds = MIN( zrrau, Rrho0 )
                   zavdds = ( zavdds - 1.0 ) / ( Rrho0 - 1.0 )
@@ -353 +361 @@
                   zavdds = difssf * zavdds 
                   zavddt = 0.7 * zavdds
-               ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN
                   !
-                  ! Diffusive convection case.
-                  !---------------------------
-                  ! Compute interior diffusivity for double diffusive mixing of
-                  ! temperature (Marmorino and Caldwell, 1976); 
+               ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN   ! Diffusive convection case.
+                  !                                                       !---------------------------
+                  ! Compute interior diffusivity for double diffusive mixing of temperature (Marmorino and Caldwell, 1976); 
                   ! Compute interior diffusivity for double diffusive mixing of salinity 
                   zinr   = 1. / zrrau
                   zavddt = 0.909 * EXP( 4.6 * EXP( -0.54* ( zinr - 1. ) ) ) 
                   zavddt = difsdc * zavddt
-                  IF( zrrau < 0.5) THEN
-                     zavdds = zavddt * 0.15 * zrrau
-                  ELSE
-                     zavdds = zavddt * (1.85 * zrrau - 0.85 ) 
+                  IF( zrrau < 0.5) THEN   ;   zavdds = zavddt * 0.15 * zrrau
+                  ELSE                    ;   zavdds = zavddt * (1.85 * zrrau - 0.85 ) 
                   ENDIF
                ELSE
@@ -385 +389 @@
       !---------------------------------------------------------------------
       DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1     
-            IF( nn_eos < 1) THEN   
-               zt     = tsn(ji,jj,1,jp_tem)
-               zs     = tsn(ji,jj,1,jp_sal) - 35.0
-               zh     = fsdept(ji,jj,1)
-               !  potential volumic mass
-               zrhos  = rhop(ji,jj,1)
-               zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   &   ! ratio alpha/beta
-                  &                               - 0.203814e-03 ) * zt   &
-                  &                               + 0.170907e-01 ) * zt   &
-                  &   + 0.665157e-01                                      &
-                  &   +     ( - 0.678662e-05 * zs                         &
-                  &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   &
-                  &   +   ( ( - 0.302285e-13 * zh                         &
-                  &           - 0.251520e-11 * zs                         &
-                  &           + 0.512857e-12 * zt * zt           ) * zh   &
-                  &           - 0.164759e-06 * zs                         &
-                  &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   &
-                  &                               + 0.380374e-04 ) * zh
-
-               zbeta  = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      &   ! beta
-                  &                            - 0.301985e-05 ) * zt      &
-                  &   + 0.785567e-03                                      &
-                  &   + (     0.515032e-08 * zs                           &
-                  &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     &
-                  &   +(  (   0.121551e-17 * zh                           &
-                  &         - 0.602281e-15 * zs                           &
-                  &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     &
-                  &                             + 0.408195e-10   * zs     &
-                  &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     &
-                  &                             - 0.121555e-07 ) * zh
-
-               zthermal = zbeta * zalbet / ( rcp * zrhos + epsln )
-               zhalin   = zbeta * tsn(ji,jj,1,jp_sal) * rcs
-            ELSE
-               zrhos    = rhop(ji,jj,1) + rau0 * ( 1. - tmask(ji,jj,1) )
-               zthermal = rn_alpha / ( rcp * zrhos + epsln )
-               zhalin   = rn_beta * tsn(ji,jj,1,jp_sal) * rcs
-               zbeta    = rn_beta
-            ENDIF
+         DO ji = fs_2, fs_jpim1           
+            zrhos    = rau0 * ( 1._wp + rhd(ji,jj,1) ) * tmask(ji,jj,1)
+            zthermal = rab_n(ji,jj,1,jp_tem) / ( rcp * zrhos + epsln )
+            zbeta    = rab_n(ji,jj,1,jp_sal)
+            zhalin   = zbeta * tsn(ji,jj,1,jp_sal) * rcs
+            !
             ! Radiative surface buoyancy force
             zBosol(ji,jj) = grav * zthermal * qsr(ji,jj)
@@ -435 +405 @@
             ws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)                          &
                &             + sfx(ji,jj)                                     ) * rcs * tmask(ji,jj,1) 
-         ENDDO
-      ENDDO
+         END DO
+      END DO
 
       zflageos = 0.5 + SIGN( 0.5, nn_eos - 1. ) 
@@ -447 +417 @@
             ! Friction velocity (zustar), at T-point : LMD94 eq. 2
             zustar(ji,jj) = SQRT( taum(ji,jj) / ( zrhos +  epsln ) )
-         ENDDO
-      ENDDO
+         END DO
+      END DO
 
 !CDIR NOVERRCHK  
@@ -1270 +1240 @@
          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
 !!bug gm jpttdzdf ==> jpttkpp
-         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_zdf, ztrdt )
-         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_zdf, ztrds )
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )
          DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
       ENDIF
@@ -1340 +1310 @@
          IF( l_trdtrc ) THEN         ! save the non-local tracer flux trends for diagnostic
             ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:)
-            CALL trd_tra( kt, 'TRC', jn, jptra_trd_zdf, ztrtrd(:,:,:) )
+            CALL trd_tra( kt, 'TRC', jn, jptra_zdf, ztrtrd(:,:,:) )
          ENDIF
          !
@@ -1375 +1345 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk     ! dummy loop indices
+      INTEGER  ::   ios            ! local integer
 #if ! defined key_kppcustom
       INTEGER  ::   jm             ! dummy loop indices     
@@ -1382 +1353 @@
       REAL(wp) ::   zustar, zucube, zustvk, zeta, zehat   ! tempory scalars
 #endif
-      INTEGER  ::   ios            ! Local integer output status for namelist read
       REAL(wp) ::   zhbf           ! tempory scalars
       LOGICAL  ::   ll_kppcustom   ! 1st ocean level taken as surface layer

branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

@@ -6 +6 @@
    !! History :  1.0  ! 2003-08  (G. Madec)  original code
    !!            3.2  ! 2009-07  (S. Masson, G. Madec)  IOM + merge of DO-loop
+   !!            3.7  ! 2012-03  (G. Madec)  make public the density criteria for trdmxl 
+   !!             -   ! 2014-02  (F. Roquet)  mixed layer depth calculated using N2 instead of rhop 
    !!----------------------------------------------------------------------
    !!   zdf_mxl      : Compute the turbocline and mixed layer depths.
@@ -14 +16 @@
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
+   USE phycst          ! physical constants
    USE iom             ! I/O library
    USE lib_mpp         ! MPP library
@@ -25 +28 @@
    PUBLIC   zdf_mxl       ! called by step.F90
 
-   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    ! density criterion for mixed layer depth
-   REAL(wp), PUBLIC ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth
-
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by TOP)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m]
+
+   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
+   REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth
 
    !! * Substitutions
@@ -70 +73 @@
       !!      eddy diffusivity coefficient (resulting from the vertical physics
       !!      alone, not the isopycnal part, see trazdf.F) fall below a given
-      !!      value defined locally (avt_c here taken equal to 5 cm/s2)
+      !!      value defined locally (avt_c here taken equal to 5 cm/s2 by default)
       !!
       !! ** Action  :   nmln, hmld, hmlp, hmlpt
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      !!
-      INTEGER  ::   ji, jj, jk          ! dummy loop indices
-      INTEGER  ::   iikn, iiki          ! temporary integer within a do loop
-      INTEGER, POINTER, DIMENSION(:,:) ::   imld                ! temporary workspace
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   iikn, iiki, ikt, imkt   ! local integer
+      REAL(wp) ::   zN2_c        ! local scalar
+      INTEGER, POINTER, DIMENSION(:,:) ::   imld   ! 2D workspace
       !!----------------------------------------------------------------------
       !
@@ -94 +98 @@
 
       ! w-level of the mixing and mixed layers
-      nmln(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point
-      imld(:,:) = mbkt(:,:) + 1
-      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10
+      nmln(:,:)  = nlb10               ! Initialization to the number of w ocean point
+      hmlp(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2
+      zN2_c = grav * rho_c * r1_rau0   ! convert density criteria into N^2 criteria
+      DO jk = nlb10, jpkm1
+         DO jj = 1, jpj                ! Mixed layer level: w-level 
+            DO ji = 1, jpi
+               ikt = mbkt(ji,jj)
+               hmlp(ji,jj) = hmlp(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * fse3w(ji,jj,jk)
+               IF( hmlp(ji,jj) < zN2_c )   nmln(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
+            END DO
+         END DO
+      END DO
+      !
+      ! w-level of the turbocline
+      imld(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point
+      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10 
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rho_c )   nmln(ji,jj) = jk      ! Mixed layer
-               IF( avt (ji,jj,jk) < avt_c                     )   imld(ji,jj) = jk      ! Turbocline 
+               imkt = mikt(ji,jj)
+               IF( avt (ji,jj,jk) < avt_c )   imld(ji,jj) = MAX( imkt, jk )      ! Turbocline 
             END DO
          END DO
@@ -109 +126 @@
             iiki = imld(ji,jj)
             iikn = nmln(ji,jj)
-            hmld (ji,jj) = fsdepw(ji,jj,iiki  ) * tmask(ji,jj,1)    ! Turbocline depth 
-            hmlp (ji,jj) = fsdepw(ji,jj,iikn  ) * tmask(ji,jj,1)    ! Mixed layer depth
-            hmlpt(ji,jj) = fsdept(ji,jj,iikn-1)                     ! depth of the last T-point inside the mixed layer
+            imkt = mikt(ji,jj)
+            hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt )            )   * ssmask(ji,jj)    ! Turbocline depth 
+            hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj)    ! Mixed layer depth
+            hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt )            )   * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
          END DO
       END DO

branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

@@ -26 +26 @@
    !!                 !                                + cleaning of the parameters + bugs correction
    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
+   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
    !!----------------------------------------------------------------------
 #if defined key_zdftke   ||   defined key_esopa
@@ -236 +237 @@
       zfact3 = 0.5_wp       * rn_ediss
       !
+      !
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                     !  Surface boundary condition on tke
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF ( ln_isfcav ) THEN
+         DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1)
+            END DO
+         END DO
+      END IF
       DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -291 +300 @@
          END DO
          !                               ! finite LC depth
-# if defined key_vectopt_loop
-         DO jj = 1, 1
-            DO ji = 1, jpij   ! vector opt. (forced unrolling)
-# else
          DO jj = 1, jpj 
             DO ji = 1, jpi
-# endif
                zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
             END DO
@@ -313 +317 @@
                   zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
                   !                                           ! TKE Langmuir circulation source term
-                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk)
+                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
             END DO
@@ -332 +336 @@
                avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
                   &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
-                  &           / (  fse3uw_n(ji,jj,jk)         &
-                  &              * fse3uw_b(ji,jj,jk) )
+                  &                            / (  fse3uw_n(ji,jj,jk)               &
+                  &                              *  fse3uw_b(ji,jj,jk)  )
                avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
                   &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
@@ -360 +364 @@
                !                                   ! right hand side in en
                en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    &
-                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) * tmask(ji,jj,jk)
+                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) &
+                  &                                 * wmask(ji,jj,jk)
             END DO
          END DO
@@ -372 +377 @@
          END DO
       END DO
-      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
-         DO ji = fs_2, fs_jpim1    ! vector opt.
+      !
+      ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
             zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
          END DO
@@ -384 +391 @@
          END DO
       END DO
-      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
-         DO ji = fs_2, fs_jpim1    ! vector opt.
+      !
+      ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
             en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
          END DO
@@ -399 +408 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)
+               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
             END DO
          END DO
@@ -412 +421 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
-                     &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk)
+                     &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
             END DO
@@ -421 +430 @@
                jk = nmln(ji,jj)
                en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
-                  &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk)
+                  &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
             END DO
          END DO
@@ -433 +442 @@
                   ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
                   zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
-                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )    ! module of the mean stress 
+                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress 
                   zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean 
                   zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
-                     &                                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)
+                     &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
             END DO
@@ -505 +514 @@
       !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
       !
+      ! initialisation of interior minimum value (avoid a 2d loop with mikt)
+      zmxlm(:,:,:)  = rmxl_min    
+      zmxld(:,:,:)  = rmxl_min
+      !
       IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
-         zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
-         zmxlm(:,:,1) = MAX(  rn_mxl0,  zraug * taum(:,:)  )
-      ELSE                          ! surface set to the minimum value
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
+               zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) )
+            END DO
+         END DO
+      ELSE 
          zmxlm(:,:,1) = rn_mxl0
       ENDIF
-      zmxlm(:,:,jpk)  = rmxl_min     ! last level set to the interior minium value
       !
 !CDIR NOVERRCHK
@@ -520 +536 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
-               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
+               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) )
             END DO
          END DO
@@ -527 +543 @@
       !                     !* Physical limits for the mixing length
       !
-      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the zmxlm   value
+      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value 
       zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
       !
       SELECT CASE ( nn_mxl )
       !
+      ! where wmask = 0 set zmxlm == fse3w
       CASE ( 0 )           ! bounded by the distance to surface and bottom
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   &
+                  zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   &
                   &            fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) )
-                  zmxlm(ji,jj,jk) = zemxl
-                  zmxld(ji,jj,jk) = zemxl
+                  ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
+                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
+                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
                END DO
             END DO
@@ -620 +638 @@
                zsqen = SQRT( en(ji,jj,jk) )
                zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
-               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * tmask(ji,jj,jk)
-               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
+               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
+               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
                dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
             END DO
@@ -628 +646 @@
       CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
       !
-      DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points
+      DO jk = 2, jpkm1            !* vertical eddy viscosity at wu- and wv-points
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
-               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
+               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
+               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
             END DO
          END DO
@@ -653 +671 @@
 !!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!)
 !!gm              zpdlr = MAX(  0.1_wp,  ri_crit / MAX( ri_crit , zri )  )
-                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
+                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
 # if defined key_c1d
-                  e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number
-                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri
+                  e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number
+                  e_ric(ji,jj,jk) = zri   * wmask(ji,jj,jk)  ! c1d config. : save Ri
 # endif
               END DO
@@ -740 +758 @@
       !
       !                               !* Check of some namelist values
-      IF( nn_mxl  < 0  .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
-      IF( nn_pdl  < 0  .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
-      IF( nn_htau < 0  .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
-#if ! key_coupled
-      IF( nn_etau == 3 )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
-#endif
+      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
+      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
+      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
+      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
 
       IF( ln_mxl0 ) THEN
@@ -765 +781 @@
       !                               !* set vertical eddy coef. to the background value
       DO jk = 1, jpk
-         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
-         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
-         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
-         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
+         avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
+         avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
+         avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
+         avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
       END DO
       dissl(:,:,:) = 1.e-12_wp
@@ -819 +835 @@
               en (:,:,:) = rn_emin * tmask(:,:,:)
               CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
+              !
+              avt_k (:,:,:) = avt (:,:,:)
+              avm_k (:,:,:) = avm (:,:,:)
+              avmu_k(:,:,:) = avmu(:,:,:)
+              avmv_k(:,:,:) = avmv(:,:,:)
+              !
               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
            ENDIF
@@ -824 +846 @@
            en(:,:,:) = rn_emin * tmask(:,:,:)
            DO jk = 1, jpk                           ! set the Kz to the background value
-              avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
-              avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
-              avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
-              avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
+              avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
+              avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
+              avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
+              avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
            END DO
         ENDIF

branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

@@ -126 +126 @@
       zkz(:,:) = 0.e0               !* Associated potential energy consummed over the whole water column
       DO jk = 2, jpkm1
-         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk)
+         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk) * wmask(:,:,jk)
       END DO
 
@@ -135 +135 @@
       END DO
 
-      DO jk = 2, jpkm1              !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s
-         zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. )   !kz max = 300 cm2/s
+      DO jk = 2, jpkm1     !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s
+         DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
+            DO ji = 1, jpi
+               zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk)  !kz max = 300 cm2/s
+            END DO
+         END DO
       END DO
 
@@ -163 +167 @@
       !                          ! ----------------------- !
       DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing
-         avt(:,:,jk) = avt(:,:,jk) + zav_tide(:,:,jk)
-         avm(:,:,jk) = avm(:,:,jk) + zav_tide(:,:,jk)
+         DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
+            DO ji = 1, jpi
+               avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
+               avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      
+      DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1  ! vector opt.
-               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
-               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
+               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
+               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
             END DO
          END DO
@@ -237 +248 @@
       zsum2(:,:) = 0.e0
       DO jk= 2, jpk
-         zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk)
-         zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk)                
+         zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)
+         zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)               
       END DO
       DO jj = 1, jpj
@@ -274 +285 @@
       zkz(:,:) = 0.e0               ! Associated potential energy consummed over the whole water column
       DO jk = 2, jpkm1
-         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk)
+         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)
       END DO
 
@@ -284 +295 @@
 
       DO jk = 2, jpkm1              ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zavt_itf bound by 300 cm2/s
-         zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. )   ! kz max = 120 cm2/s
+         zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. ) * tmask(:,:,jk) * tmask(:,:,jk-1)   ! kz max = 120 cm2/s
       END DO
 
@@ -414 +425 @@
       ! only the energy available for mixing is taken into account,
       ! (mixing efficiency tidal dissipation efficiency)
-      en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * tmask(:,:,1)
-
+      en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * ssmask(:,:)
+
+!============
+!TG: Bug for VVL? Should this section be moved out of _init and be updated at every timestep?
       ! Vertical structure (az_tmx)
       DO jj = 1, jpj                ! part independent of the level
          DO ji = 1, jpi
-            zhdep(ji,jj) = fsdepw(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
+            zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
             zfact(ji,jj) = rau0 * rn_htmx * ( 1. - EXP( -zhdep(ji,jj) / rn_htmx ) )
             IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = en_tmx(ji,jj) / zfact(ji,jj)
@@ -427 +440 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               az_tmx(ji,jj,jk) = zfact(ji,jj) * EXP( -( zhdep(ji,jj)-fsdepw(ji,jj,jk) ) / rn_htmx ) * tmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
+               az_tmx(ji,jj,jk) = zfact(ji,jj) * EXP( -( zhdep(ji,jj)-gdepw_0(ji,jj,jk) ) / rn_htmx ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+!===========
 
       IF( nprint == 1 .AND. lwp ) THEN
@@ -446 +460 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
+                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
                END DO
             END DO
@@ -460 +474 @@
          zkz(:,:) = 0.e0
          DO jk = 2, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk)* tmask(ji,jj,jk)
-            END DO
-         END DO
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX(0.e0, rn2(ji,jj,jk)) * rau0 * zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
+               END DO
+            END DO
          END DO
          ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz
@@ -485 +499 @@
 
          DO jk = 2, jpkm1
-            zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. )   !kz max = 300 cm2/s
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk)  !kz max = 300 cm2/s
+               END DO
+            END DO
          END DO
          ztpc = 0.e0
@@ -492 +510 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
+                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
                END DO
             END DO
@@ -501 +519 @@
          DO jk = 1, jpk
             ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zav_tide(:,:,jk)     * tmask_i(:,:) )   &
-               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )
+               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) )
             ztpc = 1.E50
             DO jj = 1, jpj
@@ -522 +540 @@
             END DO
             ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zkz(:,:)     * tmask_i(:,:) )   &
-               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )
+               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) )
             WRITE(numout,*) '                jk= ', jk,'   ', ze_z * 1.e4,' cm2/s'
          END DO
@@ -528 +546 @@
             zkz(:,:) = az_tmx(:,:,jk) /rn_n2min
             ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zkz(:,:)     * tmask_i(:,:) )   &
-               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )
+               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) )
             WRITE(numout,*) 
             WRITE(numout,*) '          N2 min - jk= ', jk,'   ', ze_z * 1.e4,' cm2/s min= ',MINVAL(zkz)*1.e4,   &