Ignore:
Timestamp:
2015-07-09T12:14:37+02:00 (5 years ago)
Author:
davestorkey
Message:

Update UKMO/dev_r5107_hadgem3_cplseq branch to trunk revision 5518
(= branching point of NEMO 3.6_stable).

Location:
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/ZDF
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r5477 r5572  
    120120                  zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp) 
    121121                  zbfrt(ji,jj) = MIN(zbfrt(ji,jj), rn_bfri2_max) 
    122 ! (ISF) 
    123                   ikbt = mikt(ji,jj) 
    124 ! JC: possible WAD implementation should modify line below if layers vanish 
    125                   ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 
    126                   ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 
    127                   ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max) 
    128  
    129122               END DO 
    130123            END DO 
     124! (ISF) 
     125            IF ( ln_isfcav ) THEN 
     126               DO jj = 1, jpj 
     127                  DO ji = 1, jpi 
     128                     ikbt = mikt(ji,jj) 
     129! JC: possible WAD implementation should modify line below if layers vanish 
     130                     ztmp = (1-tmask(ji,jj,1)) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 
     131                     ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 
     132                     ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max) 
     133                  END DO 
     134               END DO 
     135            END IF 
    131136         !    
    132137         ELSE 
     
    152157               ! 
    153158               ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 
    154                IF ( miku(ji,jj) + 2 .GE. mbku(ji,jj) ) THEN 
    155                   bfrua(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) )   & 
    156                                &            + ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) ) & 
    157                                &          * zecu * (1._wp - umask(ji,jj,1)) 
    158                END IF 
    159                IF ( mikv(ji,jj) + 2 .GE. mbkv(ji,jj) ) THEN 
    160                   bfrva(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) )   & 
    161                                &            + ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) ) & 
    162                                &          * zecv * (1._wp - vmask(ji,jj,1)) 
    163                END IF 
    164                ! (ISF) ======================================================================== 
    165                ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points  
    166                ikbv = mikv(ji,jj)         ! (deepest ocean u- and v-points) 
    167                ! 
    168                zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     & 
    169                   &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  ) 
    170                zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     & 
    171                   &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  ) 
    172                ! 
    173                zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 ) 
    174                zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 ) 
    175                ! 
    176                tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu * (1._wp - umask(ji,jj,1)) 
    177                tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) * zecv * (1._wp - vmask(ji,jj,1)) 
    178                ! (ISF) END ==================================================================== 
    179                ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 
    180                IF ( miku(ji,jj) + 2 .GE. mbku(ji,jj) ) THEN 
    181                   tfrua(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) )   & 
    182                                &            + ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) ) & 
    183                                &          * zecu * (1._wp - umask(ji,jj,1)) 
    184                END IF 
    185                IF ( mikv(ji,jj) + 2 .GE. mbkv(ji,jj) ) THEN 
    186                   tfrva(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) )   & 
    187                                &            + ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) ) & 
    188                                &          * zecv * (1._wp - vmask(ji,jj,1)) 
     159               IF ( ln_isfcav ) THEN 
     160                  IF ( miku(ji,jj) + 1 .GE. mbku(ji,jj) ) THEN 
     161                     bfrua(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) )   & 
     162                                  &            + ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) ) & 
     163                                  &          * zecu * (1._wp - umask(ji,jj,1)) 
     164                  END IF 
     165                  IF ( mikv(ji,jj) + 1 .GE. mbkv(ji,jj) ) THEN 
     166                     bfrva(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) )   & 
     167                                  &            + ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) ) & 
     168                                  &          * zecv * (1._wp - vmask(ji,jj,1)) 
     169                  END IF 
    189170               END IF 
    190171            END DO 
    191172         END DO 
    192          ! 
    193173         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition 
     174 
     175         IF ( ln_isfcav ) THEN 
     176            DO jj = 2, jpjm1 
     177               DO ji = 2, jpim1 
     178                  ! (ISF) ======================================================================== 
     179                  ikbu = miku(ji,jj)         ! ocean top level at u- and v-points  
     180                  ikbv = mikv(ji,jj)         ! (1st wet ocean u- and v-points) 
     181                  ! 
     182                  zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     & 
     183                     &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  ) 
     184                  zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     & 
     185                     &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  ) 
     186              ! 
     187                  zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_tfeb2 ) 
     188                  zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_tfeb2 ) 
     189              ! 
     190                  tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu * (1._wp - umask(ji,jj,1)) 
     191                  tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) * zecv * (1._wp - vmask(ji,jj,1)) 
     192              ! (ISF) END ==================================================================== 
     193              ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 
     194                  IF ( miku(ji,jj) + 1 .GE. mbku(ji,jj) ) THEN 
     195                     tfrua(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) )   & 
     196                                  &            + ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) ) & 
     197                                  &          * zecu * (1._wp - umask(ji,jj,1)) 
     198                  END IF 
     199                  IF ( mikv(ji,jj) + 1 .GE. mbkv(ji,jj) ) THEN 
     200                     tfrva(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) )   & 
     201                                  &            + ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) ) & 
     202                                  &          * zecv * (1._wp - vmask(ji,jj,1)) 
     203                  END IF 
     204               END DO 
     205            END DO 
     206            CALL lbc_lnk( tfrua, 'U', 1. )   ;   CALL lbc_lnk( tfrva, 'V', 1. )      ! Lateral boundary condition 
     207         END IF 
     208         ! 
    194209         ! 
    195210         IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        & 
     
    264279            IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien 
    265280         ENDIF 
    266          IF(lwp) WRITE(numout,*) '      top    friction coef.   rn_bfri1  = ', rn_bfri1 
    267          IF( ln_tfr2d ) THEN 
    268             IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d 
    269             IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien 
    270          ENDIF 
     281         IF ( ln_isfcav ) THEN 
     282            IF(lwp) WRITE(numout,*) '      top    friction coef.   rn_bfri1  = ', rn_tfri1 
     283            IF( ln_tfr2d ) THEN 
     284               IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d 
     285               IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien 
     286            ENDIF 
     287         END IF 
    271288         ! 
    272289         IF(ln_bfr2d) THEN 
     
    282299         bfrua(:,:) = - bfrcoef2d(:,:) 
    283300         bfrva(:,:) = - bfrcoef2d(:,:) 
     301         ! 
     302         IF ( ln_isfcav ) THEN 
     303            IF(ln_tfr2d) THEN 
     304               ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 
     305               CALL iom_open('tfr_coef.nc',inum) 
     306               CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array 
     307               CALL iom_close(inum) 
     308               tfrcoef2d(:,:) = rn_tfri1 * ( 1 + rn_tfrien * tfrcoef2d(:,:) ) 
     309            ELSE 
     310               tfrcoef2d(:,:) = rn_tfri1  ! initialize tfrcoef2d to the namelist variable 
     311            ENDIF 
     312            ! 
     313            tfrua(:,:) = - tfrcoef2d(:,:) 
     314            tfrva(:,:) = - tfrcoef2d(:,:) 
     315         END IF 
    284316         ! 
    285317      CASE( 2 ) 
     
    298330            IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien 
    299331         ENDIF 
    300          IF(lwp) WRITE(numout,*) '      quadratic top    friction' 
    301          IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri2  = ', rn_tfri2 
    302          IF(lwp) WRITE(numout,*) '      Max. coef. (log case)   rn_tfri2_max  = ', rn_tfri2_max 
    303          IF(lwp) WRITE(numout,*) '      background tke   rn_tfeb2  = ', rn_tfeb2 
    304          IF(lwp) WRITE(numout,*) '      log formulation   ln_tfr2d = ', ln_loglayer 
    305          IF(lwp) WRITE(numout,*) '      bottom roughness  rn_tfrz0 [m] = ', rn_tfrz0 
    306          IF( rn_tfrz0<=0.e0 ) THEN 
    307             WRITE(ctmp1,*) '      bottom roughness must be strictly positive' 
    308             CALL ctl_stop( ctmp1 ) 
    309          ENDIF 
    310          IF( ln_tfr2d ) THEN 
    311             IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d 
    312             IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien 
    313          ENDIF 
     332         IF ( ln_isfcav ) THEN 
     333            IF(lwp) WRITE(numout,*) '      quadratic top    friction' 
     334            IF(lwp) WRITE(numout,*) '      friction coef.    rn_tfri2     = ', rn_tfri2 
     335            IF(lwp) WRITE(numout,*) '      Max. coef. (log case)   rn_tfri2_max  = ', rn_tfri2_max 
     336            IF(lwp) WRITE(numout,*) '      background tke    rn_tfeb2     = ', rn_tfeb2 
     337            IF(lwp) WRITE(numout,*) '      log formulation   ln_tfr2d     = ', ln_loglayer 
     338            IF(lwp) WRITE(numout,*) '      top roughness     rn_tfrz0 [m] = ', rn_tfrz0 
     339            IF( rn_tfrz0<=0.e0 ) THEN 
     340               WRITE(ctmp1,*) '      top roughness must be strictly positive' 
     341               CALL ctl_stop( ctmp1 ) 
     342            ENDIF 
     343            IF( ln_tfr2d ) THEN 
     344               IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d 
     345               IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien 
     346            ENDIF 
     347         END IF 
    314348         ! 
    315349         IF(ln_bfr2d) THEN 
     
    323357            bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable 
    324358         ENDIF 
     359          
     360         IF ( ln_isfcav ) THEN 
     361            IF(ln_tfr2d) THEN 
     362               ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 
     363               CALL iom_open('tfr_coef.nc',inum) 
     364               CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array 
     365               CALL iom_close(inum) 
     366               ! 
     367               tfrcoef2d(:,:) = rn_tfri2 * ( 1 + rn_tfrien * tfrcoef2d(:,:) ) 
     368            ELSE 
     369               tfrcoef2d(:,:) = rn_tfri2  ! initialize tfrcoef2d to the namelist variable 
     370            ENDIF 
     371         END IF 
    325372         ! 
    326373         IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all 
     
    333380               END DO 
    334381            END DO 
     382            IF ( ln_isfcav ) THEN 
     383               DO jj = 1, jpj 
     384                  DO ji = 1, jpi 
     385                     ikbt = mikt(ji,jj) 
     386                     ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_tfrz0 ))**2._wp 
     387                     tfrcoef2d(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 
     388                     tfrcoef2d(ji,jj) = MIN(tfrcoef2d(ji,jj), rn_tfri2_max) 
     389                  END DO 
     390               END DO 
     391            END IF 
    335392         ENDIF 
    336393         ! 
     
    385442             zminbfr = MIN(  zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) )  ) 
    386443             zmaxbfr = MAX(  zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) )  ) 
     444! (ISF) 
     445             IF ( ln_isfcav ) THEN 
     446                ikbu = miku(ji,jj)       ! 1st wet ocean level at u- and v-points 
     447                ikbv = mikv(ji,jj) 
     448                zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt 
     449                zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt 
     450                IF( ABS( tfrcoef2d(ji,jj) ) > zfru ) THEN 
     451                   IF( ln_ctl ) THEN 
     452                      WRITE(numout,*) 'TFR ', narea, nimpp+ji, njmpp+jj, ikbu 
     453                      WRITE(numout,*) 'TFR ', ABS( tfrcoef2d(ji,jj) ), zfru 
     454                   ENDIF 
     455                   ictu = ictu + 1 
     456                ENDIF 
     457                IF( ABS( tfrcoef2d(ji,jj) ) > zfrv ) THEN 
     458                   IF( ln_ctl ) THEN 
     459                      WRITE(numout,*) 'TFR ', narea, nimpp+ji, njmpp+jj, ikbv 
     460                      WRITE(numout,*) 'TFR ', tfrcoef2d(ji,jj), zfrv 
     461                   ENDIF 
     462                   ictv = ictv + 1 
     463                ENDIF 
     464                zmintfr = MIN(  zmintfr, MIN( zfru, ABS( tfrcoef2d(ji,jj) ) )  ) 
     465                zmaxtfr = MAX(  zmaxtfr, MIN( zfrv, ABS( tfrcoef2d(ji,jj) ) )  ) 
     466             END IF 
     467! END ISF 
    387468         END DO 
    388469      END DO 
     
    392473         CALL mpp_min( zminbfr ) 
    393474         CALL mpp_max( zmaxbfr ) 
     475         IF ( ln_isfcav) CALL mpp_min( zmintfr ) 
     476         IF ( ln_isfcav) CALL mpp_max( zmaxtfr ) 
    394477      ENDIF 
    395478      IF( .NOT.ln_bfrimp) THEN 
    396479      IF( lwp .AND. ictu + ictv > 0 ) THEN 
    397          WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points ' 
    398          WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points ' 
     480         WRITE(numout,*) ' Bottom/Top friction stability check failed at ', ictu, ' U-points ' 
     481         WRITE(numout,*) ' Bottom/Top friction stability check failed at ', ictv, ' V-points ' 
    399482         WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr 
    400          WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr 
    401          WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary' 
     483         IF ( ln_isfcav ) WRITE(numout,*) ' Top friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr 
     484         WRITE(numout,*) ' Bottom/Top friction coefficient will be reduced where necessary' 
    402485      ENDIF 
    403486      ENDIF 
  • branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

    r5477 r5572  
    156156         END DO 
    157157         ! mask zmsk in order to have avt and avs masked 
    158          zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk) 
     158         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk) 
    159159 
    160160 
     
    191191               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    & 
    192192                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   & 
    193                   &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * umask(ji,jj,jk) 
     193                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * wumask(ji,jj,jk) 
    194194               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    & 
    195195                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   & 
    196                   &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * vmask(ji,jj,jk) 
     196                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * wvmask(ji,jj,jk) 
    197197            END DO 
    198198         END DO 
     
    255255      IF( zdf_ddm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' ) 
    256256      !                               ! initialization to masked Kz 
    257       avs(:,:,:) = rn_avt0 * tmask(:,:,:)  
     257      avs(:,:,:) = rn_avt0 * wmask(:,:,:)  
    258258      ! 
    259259   END SUBROUTINE zdf_ddm_init 
  • branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r5477 r5572  
    2020   USE domvvl         ! ocean space and time domain : variable volume layer 
    2121   USE zdf_oce        ! ocean vertical physics 
     22   USE zdfbfr         ! bottom friction (only for rn_bfrz0) 
    2223   USE sbc_oce        ! surface boundary condition: ocean 
    2324   USE phycst         ! physical constants 
     
    5253 
    5354   !                              !! ** Namelist  namzdf_gls  ** 
    54    LOGICAL  ::   ln_crban          ! =T use Craig and Banner scheme 
    5555   LOGICAL  ::   ln_length_lim     ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988) 
    5656   LOGICAL  ::   ln_sigpsi         ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing 
    57    INTEGER  ::   nn_tkebc_surf     ! TKE surface boundary condition (=0/1) 
    58    INTEGER  ::   nn_tkebc_bot      ! TKE bottom boundary condition (=0/1) 
    59    INTEGER  ::   nn_psibc_surf     ! PSI surface boundary condition (=0/1) 
    60    INTEGER  ::   nn_psibc_bot      ! PSI bottom boundary condition (=0/1) 
     57   INTEGER  ::   nn_bc_surf        ! surface boundary condition (=0/1) 
     58   INTEGER  ::   nn_bc_bot         ! bottom boundary condition (=0/1) 
     59   INTEGER  ::   nn_z0_met         ! Method for surface roughness computation 
    6160   INTEGER  ::   nn_stab_func      ! stability functions G88, KC or Canuto (=0/1/2) 
    6261   INTEGER  ::   nn_clos           ! closure 0/1/2/3 MY82/k-eps/k-w/gen 
     
    6665   REAL(wp) ::   rn_charn          ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 
    6766   REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing 
    68  
    69    REAL(wp) ::   hsro          =  0.003_wp    ! Minimum surface roughness 
    70    REAL(wp) ::   hbro          =  0.003_wp    ! Bottom roughness (m) 
     67   REAL(wp) ::   rn_hsro           ! Minimum surface roughness 
     68   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1)  
     69 
    7170   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters 
    7271   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1  
     
    9695   REAL(wp) ::   rm7           =  0.0_wp 
    9796   REAL(wp) ::   rm8           =  0.318_wp 
    98     
     97   REAL(wp) ::   rtrans        =  0.1_wp 
    9998   REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters 
    100    REAL(wp) ::   rc03_sqrt2_galp                                  !     -           -           -        - 
    101    REAL(wp) ::   rsbc_tke1, rsbc_tke2, rsbc_tke3, rfact_tke       !     -           -           -        - 
    102    REAL(wp) ::   rsbc_psi1, rsbc_psi2, rsbc_psi3, rfact_psi       !     -           -           -        - 
    103    REAL(wp) ::   rsbc_mb  , rsbc_std , rsbc_zs                    !     -           -           -        - 
     99   REAL(wp) ::   rsbc_tke1, rsbc_tke2, rfact_tke                  !     -           -           -        - 
     100   REAL(wp) ::   rsbc_psi1, rsbc_psi2, rfact_psi                  !     -           -           -        - 
     101   REAL(wp) ::   rsbc_zs1, rsbc_zs2                               !     -           -           -        - 
    104102   REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        - 
    105103   REAL(wp) ::   rs0, rs1, rs2, rs4, rs5, rs6                     !     -           -           -        - 
     
    147145      REAL(wp) ::   gh, gm, shr, dif, zsqen, zav        !   -      - 
    148146      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zdep 
     147      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zkar 
    149148      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves  
    150149      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhsro       ! Surface roughness (surface waves) 
     
    153152      REAL(wp), POINTER, DIMENSION(:,:,:) ::   shear       ! vertical shear 
    154153      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eps         ! dissipation rate 
    155       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi.AND.ln_crban=T) 
    156       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_a, z_elem_b, z_elem_c, psi 
     154      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi) 
     155      REAL(wp), POINTER, DIMENSION(:,:,:) ::   psi         ! psi at time now 
     156      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_a    ! element of the first  matrix diagonal 
     157      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_b    ! element of the second matrix diagonal 
     158      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_c    ! element of the third  matrix diagonal 
    157159      !!-------------------------------------------------------------------- 
    158160      ! 
    159161      IF( nn_timing == 1 )  CALL timing_start('zdf_gls') 
    160162      ! 
    161       CALL wrk_alloc( jpi,jpj, zdep, zflxs, zhsro ) 
    162       CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 
    163  
     163      CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 
     164      CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi  ) 
     165       
    164166      ! Preliminary computing 
    165167 
     
    174176 
    175177      ! Compute surface and bottom friction at T-points 
    176 !CDIR NOVERRCHK 
    177       DO jj = 2, jpjm1 
    178 !CDIR NOVERRCHK 
    179          DO ji = fs_2, fs_jpim1   ! vector opt. 
    180             !  
    181             ! surface friction  
     178!CDIR NOVERRCHK           
     179      DO jj = 2, jpjm1           
     180!CDIR NOVERRCHK          
     181         DO ji = fs_2, fs_jpim1   ! vector opt.          
     182            ! 
     183            ! surface friction 
    182184            ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 
    183             ! 
    184             ! bottom friction (explicit before friction) 
    185             ! Note that we chose here not to bound the friction as in dynbfr) 
    186             ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   & 
    187                & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  ) 
    188             zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   & 
    189                & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  ) 
    190             ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 
    191          END DO 
    192       END DO   
    193  
    194       ! In case of breaking surface waves mixing, 
    195       ! Compute surface roughness length according to Charnock formula: 
    196       IF( ln_crban ) THEN   ;   zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro) 
    197       ELSE                  ;   zhsro(:,:) = hsro 
    198       ENDIF 
     185            !    
     186            ! bottom friction (explicit before friction)         
     187            ! Note that we chose here not to bound the friction as in dynbfr)    
     188            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &          
     189               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )       
     190            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &          
     191               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )       
     192            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)          
     193         END DO          
     194      END DO     
     195 
     196      ! Set surface roughness length 
     197      SELECT CASE ( nn_z0_met ) 
     198      ! 
     199      CASE ( 0 )             ! Constant roughness           
     200         zhsro(:,:) = rn_hsro 
     201      CASE ( 1 )             ! Standard Charnock formula 
     202         zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro) 
     203      CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 
     204         zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10) 
     205         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 
     206      ! 
     207      END SELECT 
    199208 
    200209      ! Compute shear and dissipation rate 
     
    303312      ! 
    304313      ! Set surface condition on zwall_psi (1 at the bottom) 
    305       IF( ln_sigpsi ) THEN 
    306          zcoef = rsc_psi / rsc_psi0 
    307          DO jj = 2, jpjm1 
    308             DO ji = fs_2, fs_jpim1   ! vector opt. 
    309                zwall_psi(ji,jj,1) = zcoef 
    310             END DO 
    311          END DO 
    312       ENDIF 
    313  
     314      zwall_psi(:,:,1) = zwall_psi(:,:,2) 
     315      zwall_psi(:,:,jpk) = 1. 
     316      ! 
    314317      ! Surface boundary condition on tke 
    315318      ! --------------------------------- 
    316319      ! 
    317       SELECT CASE ( nn_tkebc_surf ) 
     320      SELECT CASE ( nn_bc_surf ) 
    318321      ! 
    319322      CASE ( 0 )             ! Dirichlet case 
    320          ! 
    321          IF (ln_crban) THEN     ! Wave induced mixing case 
    322             !                      ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 
    323             !                      ! balance between the production and the dissipation terms including the wave effect 
    324             en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 
    325             z_elem_a(:,:,1) = en(:,:,1) 
    326             z_elem_c(:,:,1) = 0._wp 
    327             z_elem_b(:,:,1) = 1._wp 
    328             !  
    329             ! one level below 
    330             en(:,:,2) = MAX( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 
    331             z_elem_a(:,:,2) = 0._wp 
    332             z_elem_c(:,:,2) = 0._wp 
    333             z_elem_b(:,:,2) = 1._wp 
    334             ! 
    335          ELSE                   ! No wave induced mixing case 
    336             !                      ! en(1) = u*^2/C0^2  &  l(1)  = K*zs 
    337             !                      ! balance between the production and the dissipation terms 
    338             en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 
    339             z_elem_a(:,:,1) = en(:,:,1)  
    340             z_elem_c(:,:,1) = 0._wp 
    341             z_elem_b(:,:,1) = 1._wp 
    342             ! 
    343             ! one level below 
    344             en(:,:,2) = MAX( rc02r * ustars2(:,:), rn_emin ) 
    345             z_elem_a(:,:,2) = 0._wp 
    346             z_elem_c(:,:,2) = 0._wp 
    347             z_elem_b(:,:,2) = 1._wp 
    348             ! 
    349          ENDIF 
    350          ! 
     323      ! First level 
     324      en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
     325      en(:,:,1) = MAX(en(:,:,1), rn_emin)  
     326      z_elem_a(:,:,1) = en(:,:,1) 
     327      z_elem_c(:,:,1) = 0._wp 
     328      z_elem_b(:,:,1) = 1._wp 
     329      !  
     330      ! One level below 
     331      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
     332      en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
     333      z_elem_a(:,:,2) = 0._wp  
     334      z_elem_c(:,:,2) = 0._wp 
     335      z_elem_b(:,:,2) = 1._wp 
     336      ! 
     337      ! 
    351338      CASE ( 1 )             ! Neumann boundary condition on d(e)/dz 
    352          ! 
    353          IF (ln_crban) THEN ! Shear free case: d(e)/dz= Fw 
    354             ! 
    355             ! Dirichlet conditions at k=1 (Cosmetic) 
    356             en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 
    357             z_elem_a(:,:,1) = en(:,:,1) 
    358             z_elem_c(:,:,1) = 0._wp 
    359             z_elem_b(:,:,1) = 1._wp 
    360             ! at k=2, set de/dz=Fw 
    361             z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    362             z_elem_a(:,:,2) = 0._wp         
    363             zflxs(:,:) = rsbc_tke3 * ustars2(:,:)**1.5_wp * ( (zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 
    364             en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
    365             ! 
    366          ELSE                   ! No wave induced mixing case: d(e)/dz=0. 
    367             ! 
    368             ! Dirichlet conditions at k=1 (Cosmetic) 
    369             en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 
    370             z_elem_a(:,:,1) = en(:,:,1) 
    371             z_elem_c(:,:,1) = 0._wp 
    372             z_elem_b(:,:,1) = 1._wp 
    373             ! at k=2 set de/dz=0.: 
    374             z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2)  ! Remove z_elem_a from z_elem_b 
    375             z_elem_a(:,:,2) = 0._wp 
    376             ! 
    377          ENDIF 
    378          ! 
     339      ! 
     340      ! Dirichlet conditions at k=1 
     341      en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
     342      en(:,:,1)       = MAX(en(:,:,1), rn_emin)       
     343      z_elem_a(:,:,1) = en(:,:,1) 
     344      z_elem_c(:,:,1) = 0._wp 
     345      z_elem_b(:,:,1) = 1._wp 
     346      ! 
     347      ! at k=2, set de/dz=Fw 
     348      !cbr 
     349      z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
     350      z_elem_a(:,:,2) = 0._wp 
     351      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 
     352      zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 
     353 
     354      en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 
     355      ! 
     356      ! 
    379357      END SELECT 
    380358 
     
    382360      ! -------------------------------- 
    383361      ! 
    384       SELECT CASE ( nn_tkebc_bot ) 
     362      SELECT CASE ( nn_bc_bot ) 
    385363      ! 
    386364      CASE ( 0 )             ! Dirichlet  
     
    457435      !                                            ! set the minimum value of tke  
    458436      en(:,:,:) = MAX( en(:,:,:), rn_emin ) 
    459        
     437 
    460438      !!----------------------------------------!! 
    461439      !!   Solve prognostic equation for psi    !! 
     
    560538      ! --------------------------------- 
    561539      ! 
    562       SELECT CASE ( nn_psibc_surf ) 
     540      SELECT CASE ( nn_bc_surf ) 
    563541      ! 
    564542      CASE ( 0 )             ! Dirichlet boundary conditions 
    565          ! 
    566          IF( ln_crban ) THEN       ! Wave induced mixing case 
    567             !                      ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 
    568             !                      ! balance between the production and the dissipation terms including the wave effect 
    569             zdep(:,:) = rl_sf * zhsro(:,:) 
    570             psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    571             z_elem_a(:,:,1) = psi(:,:,1) 
    572             z_elem_c(:,:,1) = 0._wp 
    573             z_elem_b(:,:,1) = 1._wp 
    574             ! 
    575             ! one level below 
    576             zex1 = (rmm*ra_sf+rnn) 
    577             zex2 = (rmm*ra_sf) 
    578             zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 
    579             psi (:,:,2) = rsbc_psi1 * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 
    580             z_elem_a(:,:,2) = 0._wp 
    581             z_elem_c(:,:,2) = 0._wp 
    582             z_elem_b(:,:,2) = 1._wp 
    583             !  
    584          ELSE                   ! No wave induced mixing case 
    585             !                      ! en(1) = u*^2/C0^2  &  l(1)  = K*zs 
    586             !                      ! balance between the production and the dissipation terms 
    587             ! 
    588             zdep(:,:) = vkarmn * zhsro(:,:) 
    589             psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    590             z_elem_a(:,:,1) = psi(:,:,1) 
    591             z_elem_c(:,:,1) = 0._wp 
    592             z_elem_b(:,:,1) = 1._wp 
    593             ! 
    594             ! one level below 
    595             zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) ) 
    596             psi (:,:,2) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    597             z_elem_a(:,:,2) = 0._wp 
    598             z_elem_c(:,:,2) = 0._wp 
    599             z_elem_b(:,:,2) = 1. 
    600             ! 
    601          ENDIF 
    602          ! 
     543      ! 
     544      ! Surface value 
     545      zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
     546      psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     547      z_elem_a(:,:,1) = psi(:,:,1) 
     548      z_elem_c(:,:,1) = 0._wp 
     549      z_elem_b(:,:,1) = 1._wp 
     550      ! 
     551      ! One level below 
     552      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 
     553      zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 
     554      psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     555      z_elem_a(:,:,2) = 0._wp 
     556      z_elem_c(:,:,2) = 0._wp 
     557      z_elem_b(:,:,2) = 1._wp 
     558      !  
     559      ! 
    603560      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz 
    604          ! 
    605          IF( ln_crban ) THEN     ! Wave induced mixing case 
    606             ! 
    607             zdep(:,:) = rl_sf * zhsro(:,:) 
    608             psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    609             z_elem_a(:,:,1) = psi(:,:,1) 
    610             z_elem_c(:,:,1) = 0._wp 
    611             z_elem_b(:,:,1) = 1._wp 
    612             ! 
    613             ! Neumann condition at k=2 
    614             z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    615             z_elem_a(:,:,2) = 0._wp 
    616             ! 
    617             ! Set psi vertical flux at the surface: 
    618             zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 
    619             zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) &  
    620                &                   * en(:,:,1)**rmm * zdep          
    621             psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
    622             ! 
    623       ELSE                   ! No wave induced mixing 
    624             ! 
    625             zdep(:,:) = vkarmn * zhsro(:,:) 
    626             psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    627             z_elem_a(:,:,1) = psi(:,:,1) 
    628             z_elem_c(:,:,1) = 0._wp 
    629             z_elem_b(:,:,1) = 1._wp 
    630             ! 
    631             ! Neumann condition at k=2 
    632             z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    633             z_elem_a(ji,jj,2) = 0._wp 
    634             ! 
    635             ! Set psi vertical flux at the surface: 
    636             zdep(:,:)  = zhsro(:,:) + fsdept(:,:,1) 
    637             zflxs(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-1._wp) 
    638             psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
    639             !      
    640          ENDIF 
    641          ! 
     561      ! 
     562      ! Surface value: Dirichlet 
     563      zdep(:,:)       = zhsro(:,:) * rl_sf 
     564      psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     565      z_elem_a(:,:,1) = psi(:,:,1) 
     566      z_elem_c(:,:,1) = 0._wp 
     567      z_elem_b(:,:,1) = 1._wp 
     568      ! 
     569      ! Neumann condition at k=2 
     570      z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
     571      z_elem_a(:,:,2) = 0._wp 
     572      ! 
     573      ! Set psi vertical flux at the surface: 
     574      zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
     575      zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
     576      zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
     577      zdep(:,:) =  rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 
     578             & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 
     579      zflxs(:,:) = zdep(:,:) * zflxs(:,:) 
     580      psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
     581 
     582      !    
     583      ! 
    642584      END SELECT 
    643585 
     
    645587      ! -------------------------------- 
    646588      ! 
    647       SELECT CASE ( nn_psibc_bot ) 
     589      SELECT CASE ( nn_bc_bot ) 
     590      ! 
    648591      ! 
    649592      CASE ( 0 )             ! Dirichlet  
    650          !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * hbro 
     593         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0 
    651594         !                      ! Balance between the production and the dissipation terms 
    652595!CDIR NOVERRCHK 
     
    656599               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
    657600               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
    658                zdep(ji,jj) = vkarmn * hbro 
     601               zdep(ji,jj) = vkarmn * rn_bfrz0 
    659602               psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 
    660603               z_elem_a(ji,jj,ibot) = 0._wp 
     
    663606               ! 
    664607               ! Just above last level, Dirichlet condition again (GOTM like) 
    665                zdep(ji,jj) = vkarmn * ( hbro + fse3t(ji,jj,ibotm1) ) 
     608               zdep(ji,jj) = vkarmn * ( rn_bfrz0 + fse3t(ji,jj,ibotm1) ) 
    666609               psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn 
    667610               z_elem_a(ji,jj,ibotm1) = 0._wp 
     
    681624               ! 
    682625               ! Bottom level Dirichlet condition: 
    683                zdep(ji,jj) = vkarmn * hbro 
     626               zdep(ji,jj) = vkarmn * rn_bfrz0 
    684627               psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 
    685628               ! 
     
    693636               ! 
    694637               ! Set psi vertical flux at the bottom: 
    695                zdep(ji,jj) = hbro + 0.5_wp*fse3t(ji,jj,ibotm1) 
     638               zdep(ji,jj) = rn_bfrz0 + 0.5_wp*fse3t(ji,jj,ibotm1) 
    696639               zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) )   & 
    697640                  &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 
     
    736679            DO jj = 2, jpjm1 
    737680               DO ji = fs_2, fs_jpim1   ! vector opt. 
    738                   eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / psi(ji,jj,jk) 
     681                  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) 
    739682               END DO 
    740683            END DO 
     
    783726               ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)  
    784727               zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    785                mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) ) 
     728               IF (ln_length_lim) mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) ) 
    786729            END DO 
    787730         END DO 
     
    847790      ! Boundary conditions on stability functions for momentum (Neumann): 
    848791      ! Lines below are useless if GOTM style Dirichlet conditions are used 
    849       zcoef = rcm_sf / SQRT( 2._wp ) 
     792 
     793      avmv(:,:,1) = avmv(:,:,2) 
     794 
    850795      DO jj = 2, jpjm1 
    851796         DO ji = fs_2, fs_jpim1   ! vector opt. 
    852             avmv(ji,jj,1) = zcoef 
    853          END DO 
    854       END DO 
    855       zcoef = rc0 / SQRT( 2._wp ) 
    856       DO jj = 2, jpjm1 
    857          DO ji = fs_2, fs_jpim1   ! vector opt. 
    858             avmv(ji,jj,mbkt(ji,jj)+1) = zcoef 
     797            avmv(ji,jj,mbkt(ji,jj)+1) = avmv(ji,jj,mbkt(ji,jj)) 
    859798         END DO 
    860799      END DO 
     
    900839      avmv_k(:,:,:) = avmv(:,:,:) 
    901840      ! 
    902       CALL wrk_dealloc( jpi,jpj, zdep, zflxs, zhsro ) 
     841      CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 
    903842      CALL wrk_dealloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 
    904843      ! 
     
    932871      !! 
    933872      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 
    934          &            rn_clim_galp, ln_crban, ln_sigpsi,     & 
    935          &            rn_crban, rn_charn,                    & 
    936          &            nn_tkebc_surf, nn_tkebc_bot,           & 
    937          &            nn_psibc_surf, nn_psibc_bot,           & 
     873         &            rn_clim_galp, ln_sigpsi, rn_hsro,      & 
     874         &            rn_crban, rn_charn, rn_frac_hs,        & 
     875         &            nn_bc_surf, nn_bc_bot, nn_z0_met,      & 
    938876         &            nn_stab_func, nn_clos 
    939877      !!---------------------------------------------------------- 
     
    955893         WRITE(numout,*) '~~~~~~~~~~~~' 
    956894         WRITE(numout,*) '   Namelist namzdf_gls : set gls mixing parameters' 
    957          WRITE(numout,*) '      minimum value of en                           rn_emin       = ', rn_emin 
    958          WRITE(numout,*) '      minimum value of eps                          rn_epsmin     = ', rn_epsmin 
    959          WRITE(numout,*) '      Limit dissipation rate under stable stratif.  ln_length_lim = ', ln_length_lim 
    960          WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp  = ', rn_clim_galp 
    961          WRITE(numout,*) '      TKE Surface boundary condition                nn_tkebc_surf = ', nn_tkebc_surf 
    962          WRITE(numout,*) '      TKE Bottom boundary condition                 nn_tkebc_bot  = ', nn_tkebc_bot 
    963          WRITE(numout,*) '      PSI Surface boundary condition                nn_psibc_surf = ', nn_psibc_surf 
    964          WRITE(numout,*) '      PSI Bottom boundary condition                 nn_psibc_bot  = ', nn_psibc_bot 
    965          WRITE(numout,*) '      Craig and Banner scheme                       ln_crban      = ', ln_crban 
    966          WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi     = ', ln_sigpsi 
     895         WRITE(numout,*) '      minimum value of en                           rn_emin        = ', rn_emin 
     896         WRITE(numout,*) '      minimum value of eps                          rn_epsmin      = ', rn_epsmin 
     897         WRITE(numout,*) '      Limit dissipation rate under stable stratif.  ln_length_lim  = ', ln_length_lim 
     898         WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp   = ', rn_clim_galp 
     899         WRITE(numout,*) '      TKE Surface boundary condition                nn_bc_surf     = ', nn_bc_surf 
     900         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot 
     901         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi 
    967902         WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban 
    968903         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn 
     904         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met 
     905         WRITE(numout,*) '      Wave height frac. (used if nn_z0_met=2)       rn_frac_hs     = ', rn_frac_hs 
    969906         WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func 
    970907         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos 
    971          WRITE(numout,*) '   Hard coded parameters' 
    972          WRITE(numout,*) '      Surface roughness (m)                         hsro          = ', hsro 
    973          WRITE(numout,*) '      Bottom roughness (m)                          hbro          = ', hbro 
     908         WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro 
     909         WRITE(numout,*) '      Bottom roughness (m) (nambfr namelist)        rn_bfrz0       = ', rn_bfrz0 
    974910      ENDIF 
    975911 
     
    978914 
    979915      !                                !* Check of some namelist values 
    980       IF( nn_tkebc_surf < 0 .OR. nn_tkebc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_surf is 0 or 1' ) 
    981       IF( nn_psibc_surf < 0 .OR. nn_psibc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_surf is 0 or 1' ) 
    982       IF( nn_tkebc_bot  < 0 .OR. nn_tkebc_bot  > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' ) 
    983       IF( nn_psibc_bot  < 0 .OR. nn_psibc_bot  > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' ) 
     916      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 
     917      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 
     918      IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' ) 
    984919      IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) 
    985920      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) 
     
    1001936         SELECT CASE ( nn_stab_func ) 
    1002937         CASE( 0, 1 )   ;   rpsi3m = 2.53_wp       ! G88 or KC stability functions 
    1003          CASE( 2 )      ;   rpsi3m = 2.38_wp       ! Canuto A stability functions 
     938         CASE( 2 )      ;   rpsi3m = 2.62_wp       ! Canuto A stability functions 
    1004939         CASE( 3 )      ;   rpsi3m = 2.38          ! Canuto B stability functions (caution : constant not identified) 
    1005940         END SELECT 
     
    1012947         rnn     = -1._wp 
    1013948         rsc_tke =  1._wp 
    1014          rsc_psi =  1.3_wp  ! Schmidt number for psi 
     949         rsc_psi =  1.2_wp  ! Schmidt number for psi 
    1015950         rpsi1   =  1.44_wp 
    1016951         rpsi3p  =  1._wp 
     
    11401075      !                                     ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009 
    11411076      !                                     !  or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001 
    1142       IF( ln_sigpsi .AND. ln_crban ) THEN 
    1143          zcr = SQRT( 1.5_wp*rsc_tke ) * rcm_sf / vkarmn 
    1144          rsc_psi0 = vkarmn*vkarmn / ( rpsi2 * rcm_sf*rcm_sf )                       &  
    1145         &         * ( rnn*rnn - 4._wp/3._wp * zcr*rnn*rmm - 1._wp/3._wp * zcr*rnn   & 
    1146         &           + 2._wp/9._wp * rmm * zcr*zcr + 4._wp/9._wp * zcr*zcr * rmm*rmm )                                  
     1077      IF( ln_sigpsi ) THEN 
     1078         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf  
     1079         ! Verification: retrieve Burchard (2001) results by uncomenting the line below: 
     1080         ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work 
     1081         ! ra_sf = -SQRT(2./3.*rc0**3./rn_cm_sf*rn_sc_tke)/vkarmn 
     1082         rsc_psi0 = rsc_tke/(24.*rpsi2)*(-1.+(4.*rnn + ra_sf*(1.+4.*rmm))**2./(ra_sf**2.)) 
    11471083      ELSE 
    11481084         rsc_psi0 = rsc_psi 
     
    11511087      !                                !* Shear free turbulence parameters 
    11521088      ! 
    1153       ra_sf  = -4._wp * rnn * SQRT( rsc_tke ) / ( (1._wp+4._wp*rmm) * SQRT( rsc_tke )   & 
    1154          &                                      - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 
    1155       rl_sf  = rc0 * SQRT( rc0 / rcm_sf )                                                                   & 
    1156          &         * SQRT(  (  (1._wp + 4._wp*rmm + 8._wp*rmm*rmm) * rsc_tke                                & 
    1157          &                   + 12._wp * rsc_psi0 * rpsi2                                                    & 
    1158          &                   - (1._wp + 4._wp*rmm) * SQRT( rsc_tke*(rsc_tke+ 24._wp*rsc_psi0*rpsi2) )  )    & 
    1159          &                / ( 12._wp*rnn*rnn )                                                              ) 
     1089      ra_sf  = -4._wp*rnn*SQRT(rsc_tke) / ( (1._wp+4._wp*rmm)*SQRT(rsc_tke) & 
     1090               &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 
     1091 
     1092      IF ( rn_crban==0._wp ) THEN 
     1093         rl_sf = vkarmn 
     1094      ELSE 
     1095         rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp)*rsc_tke          & 
     1096                 &                                       + 12._wp * rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm) & 
     1097                 &                                                *SQRT(rsc_tke*(rsc_tke                 & 
     1098                 &                                                   + 24._wp*rsc_psi0*rpsi2)) )         & 
     1099                 &                                         /(12._wp*rnn**2.)                             & 
     1100                 &                                       ) 
     1101      ENDIF 
    11601102 
    11611103      ! 
     
    11871129      rc03  = rc02 * rc0 
    11881130      rc04  = rc03 * rc0 
    1189       rc03_sqrt2_galp = rc03 / SQRT(2._wp) / rn_clim_galp 
    1190       rsbc_mb   = 0.5_wp * (15.8_wp*rn_crban)**(2._wp/3._wp)               ! Surf. bound. cond. from Mellor and Blumberg 
    1191       rsbc_std  = 3.75_wp                                                  ! Surf. bound. cond. standard (prod=diss) 
    1192       rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp)  ! k_eps = 53.  Dirichlet + Wave breaking  
    1193       rsbc_tke2 = 0.5_wp / rau0 
    1194       rsbc_tke3 = rdt * rn_crban                                                         ! Neumann + Wave breaking 
    1195       rsbc_zs   = rn_charn / grav                                                        ! Charnock formula 
    1196       rsbc_psi1 = rc0**rpp * rsbc_tke1**rmm * rl_sf**rnn                           ! Dirichlet + Wave breaking 
    1197       rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi                   ! Neumann + NO Wave breaking  
    1198       rsbc_psi3 = -0.5_wp * rdt * rc0**rpp * rl_sf**rnn / rsc_psi  * (rnn + rmm*ra_sf) ! Neumann + Wave breaking 
    1199       rfact_tke = -0.5_wp / rsc_tke * rdt               ! Cst used for the Diffusion term of tke 
    1200       rfact_psi = -0.5_wp / rsc_psi * rdt               ! Cst used for the Diffusion term of tke 
     1131      rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking 
     1132      rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking  
     1133      zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
     1134      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer  
     1135      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness 
     1136      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness  
     1137      rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 
     1138      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking  
     1139 
     1140      rfact_tke = -0.5_wp / rsc_tke * rdt                                ! Cst used for the Diffusion term of tke 
     1141      rfact_psi = -0.5_wp / rsc_psi * rdt                                ! Cst used for the Diffusion term of tke 
    12011142 
    12021143      !                                !* Wall proximity function 
     
    12571198               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 
    12581199               en  (:,:,:) = rn_emin 
    1259                mxln(:,:,:) = 0.001         
     1200               mxln(:,:,:) = 0.05         
    12601201               avt_k (:,:,:) = avt (:,:,:) 
    12611202               avm_k (:,:,:) = avm (:,:,:) 
     
    12671208            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values' 
    12681209            en  (:,:,:) = rn_emin 
    1269             mxln(:,:,:) = 0.001        
     1210            mxln(:,:,:) = 0.05        
    12701211         ENDIF 
    12711212         ! 
     
    12731214         !                                   ! ------------------- 
    12741215         IF(lwp) WRITE(numout,*) '---- gls-rst ----' 
    1275          CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     ) 
     1216         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )  
    12761217         CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  ) 
    12771218         CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  ) 
    1278          CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 
     1219         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )  
    12791220         CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 
    12801221         CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln   ) 
  • branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90

    r5477 r5572  
    1414   !!---------------------------------------------------------------------- 
    1515   USE par_oce         ! mesh and scale factors 
    16    USE sbc_oce         ! surface module (only for nn_isf in the option compatibility test) 
    1716   USE ldftra_oce      ! ocean active tracers: lateral physics 
    1817   USE ldfdyn_oce      ! ocean dynamics lateral physics 
     
    118117      IF( ioptio == 0 .OR. ioptio > 1 .AND. .NOT. lk_esopa )   & 
    119118         &   CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' ) 
    120       IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. nn_isf .NE. 0 )   & 
     119      IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. ln_isfcav )   & 
    121120         &   CALL ctl_stop( ' only zdfcst and zdftke were tested with ice shelves cavities ' ) 
    122121      ! 
     
    125124      IF(lwp) WRITE(numout,*) '   convection :' 
    126125      ! 
    127       IF( ln_zdfnpc )   CALL ctl_stop( ' zdf_init: non penetrative convective scheme is not working',   & 
    128          &                                       ' set ln_zdfnpc to FALSE' ) 
     126#if defined key_top 
     127      IF( ln_zdfnpc )   CALL ctl_stop( ' zdf_init: npc scheme is not working with key_top' ) 
     128#endif 
    129129      ! 
    130130      ioptio = 0 
  • branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r5477 r5572  
    2626   !!                 !                                + cleaning of the parameters + bugs correction 
    2727   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     28   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
    2829   !!---------------------------------------------------------------------- 
    2930#if defined key_zdftke   ||   defined key_esopa 
     
    236237      zfact3 = 0.5_wp       * rn_ediss 
    237238      ! 
     239      ! 
    238240      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    239241      !                     !  Surface boundary condition on tke 
    240242      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     243      IF ( ln_isfcav ) THEN 
     244         DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin 
     245            DO ji = fs_2, fs_jpim1   ! vector opt. 
     246               en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1) 
     247            END DO 
     248         END DO 
     249      END IF 
    241250      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    242251         DO ji = fs_2, fs_jpim1   ! vector opt. 
    243             IF (mikt(ji,jj) .GT. 1) THEN 
    244                en(ji,jj,mikt(ji,jj))=rn_emin 
    245             ELSE 
    246                en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
    247             END IF 
     252            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
    248253         END DO 
    249254      END DO 
     
    301306         END DO 
    302307         zcof = 0.016 / SQRT( zrhoa * zcdrag ) 
     308!CDIR NOVERRCHK 
    303309         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en 
    304             DO jj = 2, jpjm1 
     310!CDIR NOVERRCHK 
     311            DO jj = 2, jpjm1 
     312!CDIR NOVERRCHK 
    305313               DO ji = fs_2, fs_jpim1   ! vector opt. 
    306314                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
     
    309317                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 
    310318                  !                                           ! TKE Langmuir circulation source term 
    311                   en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk) 
     319                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    312320               END DO 
    313321            END DO 
     
    328336               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   & 
    329337                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   &  
    330                   &           / (  fse3uw_n(ji,jj,jk)         & 
    331                   &              * fse3uw_b(ji,jj,jk) ) 
     338                  &                            / (  fse3uw_n(ji,jj,jk)               & 
     339                  &                              *  fse3uw_b(ji,jj,jk) ) 
    332340               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   & 
    333341                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   & 
     
    338346      END DO 
    339347      ! 
    340       DO jj = 2, jpjm1 
    341          DO ji = fs_2, fs_jpim1   ! vector opt. 
    342             DO jk = mikt(ji,jj)+1, jpkm1           !* Matrix and right hand side in en 
     348      DO jk = 2, jpkm1           !* Matrix and right hand side in en 
     349         DO jj = 2, jpjm1 
     350            DO ji = fs_2, fs_jpim1   ! vector opt. 
    343351               zcof   = zfact1 * tmask(ji,jj,jk) 
    344352               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal 
     
    357365               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    & 
    358366                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
    359                   &                                 * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
    360             END DO 
    361             !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
    362             DO jk = mikt(ji,jj)+2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
     367                  &                                 * wmask(ji,jj,jk) 
     368            END DO 
     369         END DO 
     370      END DO 
     371      !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
     372      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
     373         DO jj = 2, jpjm1 
     374            DO ji = fs_2, fs_jpim1    ! vector opt. 
    363375               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    364376            END DO 
    365             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    366             zd_lw(ji,jj,mikt(ji,jj)+1) = en(ji,jj,mikt(ji,jj)+1) - zd_lw(ji,jj,mikt(ji,jj)+1) * en(ji,jj,mikt(ji,jj))    ! Surface boudary conditions on tke 
    367             ! 
    368             DO jk = mikt(ji,jj)+2, jpkm1 
     377         END DO 
     378      END DO 
     379      ! 
     380      ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
     381      DO jj = 2, jpjm1 
     382         DO ji = fs_2, fs_jpim1   ! vector opt. 
     383            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
     384         END DO 
     385      END DO 
     386      DO jk = 3, jpkm1 
     387         DO jj = 2, jpjm1 
     388            DO ji = fs_2, fs_jpim1    ! vector opt. 
    369389               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 
    370390            END DO 
    371             ! 
    372             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
     391         END DO 
     392      END DO 
     393      ! 
     394      ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
     395      DO jj = 2, jpjm1 
     396         DO ji = fs_2, fs_jpim1   ! vector opt. 
    373397            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 
    374             ! 
    375             DO jk = jpk-2, mikt(ji,jj)+1, -1 
     398         END DO 
     399      END DO 
     400      DO jk = jpk-2, 2, -1 
     401         DO jj = 2, jpjm1 
     402            DO ji = fs_2, fs_jpim1    ! vector opt. 
    376403               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    377404            END DO 
    378             ! 
    379             DO jk = mikt(ji,jj), jpkm1                             ! set the minimum value of tke 
    380                en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 
     405         END DO 
     406      END DO 
     407      DO jk = 2, jpkm1                             ! set the minimum value of tke 
     408         DO jj = 2, jpjm1 
     409            DO ji = fs_2, fs_jpim1   ! vector opt. 
     410               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 
    381411            END DO 
    382412         END DO 
     
    391421               DO ji = fs_2, fs_jpim1   ! vector opt. 
    392422                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    393                      &                                 * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) * tmask(ji,jj,1) 
     423                     &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    394424               END DO 
    395425            END DO 
     
    400430               jk = nmln(ji,jj) 
    401431               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    402                   &                                 * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) * tmask(ji,jj,1) 
     432                  &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    403433            END DO 
    404434         END DO 
     
    416446                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    417447                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    418                      &                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * tmask(ji,jj,1) 
     448                     &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    419449               END DO 
    420450            END DO 
     
    484514      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2) 
    485515      ! 
     516      ! initialisation of interior minimum value (avoid a 2d loop with mikt) 
     517      zmxlm(:,:,:)  = rmxl_min     
     518      zmxld(:,:,:)  = rmxl_min 
     519      ! 
    486520      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
    487521         DO jj = 2, jpjm1 
    488522            DO ji = fs_2, fs_jpim1 
    489                IF (mikt(ji,jj) .GT. 1) THEN 
    490                   zmxlm(ji,jj,mikt(ji,jj)) = rmxl_min 
    491                ELSE 
    492                   zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
    493                   zmxlm(ji,jj,mikt(ji,jj)) = MAX( rn_mxl0, zraug * taum(ji,jj) ) 
    494                END IF 
     523               zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
     524               zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 
    495525            END DO 
    496526         END DO 
    497527      ELSE  
    498          DO jj = 2, jpjm1 
    499             DO ji = fs_2, fs_jpim1                         ! surface set to the minimum value 
    500                zmxlm(ji,jj,mikt(ji,jj)) = MAX( tmask(ji,jj,1) * rn_mxl0, rmxl_min) 
    501             END DO 
    502          END DO 
     528         zmxlm(:,:,1) = rn_mxl0 
    503529      ENDIF 
    504       zmxlm(:,:,jpk)  = rmxl_min     ! last level set to the interior minium value 
    505       ! 
    506 !CDIR NOVERRCHK 
    507       DO jj = 2, jpjm1 
    508 !CDIR NOVERRCHK 
    509          DO ji = fs_2, fs_jpim1   ! vector opt. 
    510             !CDIR NOVERRCHK 
    511             DO jk = mikt(ji,jj)+1, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
     530      ! 
     531!CDIR NOVERRCHK 
     532      DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
     533!CDIR NOVERRCHK 
     534         DO jj = 2, jpjm1 
     535!CDIR NOVERRCHK 
     536            DO ji = fs_2, fs_jpim1   ! vector opt. 
    512537               zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    513                zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  ) 
    514             END DO 
    515             zmxld(ji,jj,mikt(ji,jj)) = zmxlm(ji,jj,mikt(ji,jj))   ! surface set to the minimum value  
     538               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 
     539            END DO 
    516540         END DO 
    517541      END DO 
     
    519543      !                     !* Physical limits for the mixing length 
    520544      ! 
    521       zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the zmxlm   value 
     545      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value  
    522546      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value 
    523547      ! 
    524548      SELECT CASE ( nn_mxl ) 
    525549      ! 
     550      ! where wmask = 0 set zmxlm == fse3w 
    526551      CASE ( 0 )           ! bounded by the distance to surface and bottom 
    527          DO jj = 2, jpjm1 
    528             DO ji = fs_2, fs_jpim1   ! vector opt. 
    529                DO jk = mikt(ji,jj)+1, jpkm1 
     552         DO jk = 2, jpkm1 
     553            DO jj = 2, jpjm1 
     554               DO ji = fs_2, fs_jpim1   ! vector opt. 
    530555                  zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   & 
    531556                  &            fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 
    532                   zmxlm(ji,jj,jk) = zemxl 
    533                   zmxld(ji,jj,jk) = zemxl 
     557                  ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 
     558                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 
     559                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 
    534560               END DO 
    535561            END DO 
     
    537563         ! 
    538564      CASE ( 1 )           ! bounded by the vertical scale factor 
    539          DO jj = 2, jpjm1 
    540             DO ji = fs_2, fs_jpim1   ! vector opt. 
    541                DO jk = mikt(ji,jj)+1, jpkm1 
     565         DO jk = 2, jpkm1 
     566            DO jj = 2, jpjm1 
     567               DO ji = fs_2, fs_jpim1   ! vector opt. 
    542568                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 
    543569                  zmxlm(ji,jj,jk) = zemxl 
     
    548574         ! 
    549575      CASE ( 2 )           ! |dk[xml]| bounded by e3t : 
    550          DO jj = 2, jpjm1 
    551             DO ji = fs_2, fs_jpim1   ! vector opt. 
    552                DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom : 
     576         DO jk = 2, jpkm1         ! from the surface to the bottom : 
     577            DO jj = 2, jpjm1 
     578               DO ji = fs_2, fs_jpim1   ! vector opt. 
    553579                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    554580               END DO 
    555                DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface : 
     581            END DO 
     582         END DO 
     583         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : 
     584            DO jj = 2, jpjm1 
     585               DO ji = fs_2, fs_jpim1   ! vector opt. 
    556586                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    557587                  zmxlm(ji,jj,jk) = zemxl 
     
    562592         ! 
    563593      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t : 
    564          DO jj = 2, jpjm1 
    565             DO ji = fs_2, fs_jpim1   ! vector opt. 
    566                DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom : lup 
     594         DO jk = 2, jpkm1         ! from the surface to the bottom : lup 
     595            DO jj = 2, jpjm1 
     596               DO ji = fs_2, fs_jpim1   ! vector opt. 
    567597                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    568598               END DO 
    569                DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface : ldown 
     599            END DO 
     600         END DO 
     601         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown 
     602            DO jj = 2, jpjm1 
     603               DO ji = fs_2, fs_jpim1   ! vector opt. 
    570604                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    571605               END DO 
     
    604638               zsqen = SQRT( en(ji,jj,jk) ) 
    605639               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen 
    606                avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * tmask(ji,jj,jk) 
    607                avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
     640               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk) 
     641               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
    608642               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 
    609643            END DO 
     
    612646      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
    613647      ! 
    614       DO jj = 2, jpjm1 
    615          DO ji = fs_2, fs_jpim1   ! vector opt. 
    616             DO jk = miku(ji,jj)+1, jpkm1            !* vertical eddy viscosity at u- and v-points 
    617                avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
    618             END DO 
    619             DO jk = mikv(ji,jj)+1, jpkm1 
    620                avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
     648      DO jk = 2, jpkm1            !* vertical eddy viscosity at wu- and wv-points 
     649         DO jj = 2, jpjm1 
     650            DO ji = fs_2, fs_jpim1   ! vector opt. 
     651               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk) 
     652               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk) 
    621653            END DO 
    622654         END DO 
     
    625657      ! 
    626658      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt 
    627          DO jj = 2, jpjm1 
    628             DO ji = fs_2, fs_jpim1   ! vector opt. 
    629                DO jk = mikt(ji,jj)+1, jpkm1 
     659         DO jk = 2, jpkm1 
     660            DO jj = 2, jpjm1 
     661               DO ji = fs_2, fs_jpim1   ! vector opt. 
    630662                  zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 
    631663                  !                                          ! shear 
     
    639671!!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!) 
    640672!!gm              zpdlr = MAX(  0.1_wp,  ri_crit / MAX( ri_crit , zri )  ) 
    641                   avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
     673                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
    642674# if defined key_c1d 
    643                   e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number 
    644                   e_ric(ji,jj,jk) = zri   * tmask(ji,jj,jk)  ! c1d config. : save Ri 
     675                  e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number 
     676                  e_ric(ji,jj,jk) = zri   * wmask(ji,jj,jk)  ! c1d config. : save Ri 
    645677# endif 
    646678              END DO 
     
    729761      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' ) 
    730762      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 
    731       IF( nn_etau == 3 .AND. .NOT. lk_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 
     763      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 
    732764 
    733765      IF( ln_mxl0 ) THEN 
     
    749781      !                               !* set vertical eddy coef. to the background value 
    750782      DO jk = 1, jpk 
    751          avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 
    752          avm (:,:,jk) = avmb(jk) * tmask(:,:,jk) 
    753          avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) 
    754          avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 
     783         avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 
     784         avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 
     785         avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 
     786         avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 
    755787      END DO 
    756788      dissl(:,:,:) = 1.e-12_wp 
     
    803835              en (:,:,:) = rn_emin * tmask(:,:,:) 
    804836              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation) 
     837              ! 
     838              avt_k (:,:,:) = avt (:,:,:) 
     839              avm_k (:,:,:) = avm (:,:,:) 
     840              avmu_k(:,:,:) = avmu(:,:,:) 
     841              avmv_k(:,:,:) = avmv(:,:,:) 
     842              ! 
    805843              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO 
    806844           ENDIF 
     
    808846           en(:,:,:) = rn_emin * tmask(:,:,:) 
    809847           DO jk = 1, jpk                           ! set the Kz to the background value 
    810               avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 
    811               avm (:,:,jk) = avmb(jk) * tmask(:,:,jk) 
    812               avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) 
    813               avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 
     848              avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 
     849              avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 
     850              avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 
     851              avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 
    814852           END DO 
    815853        ENDIF 
  • branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

    r5477 r5572  
    126126      zkz(:,:) = 0.e0               !* Associated potential energy consummed over the whole water column 
    127127      DO jk = 2, jpkm1 
    128          zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk) * tmask(:,:,jk-1) 
     128         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk) * wmask(:,:,jk) 
    129129      END DO 
    130130 
     
    135135      END DO 
    136136 
    137       DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
    138          DO ji = 1, jpi 
    139             DO jk = mikt(ji,jj)+1, jpkm1     !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 
    140                zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. )   !kz max = 300 cm2/s 
     137      DO jk = 2, jpkm1     !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 
     138         DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
     139            DO ji = 1, jpi 
     140               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 
    141141            END DO 
    142142         END DO 
     
    166166      !                          !   Update  mixing coefs  !                           
    167167      !                          ! ----------------------- ! 
    168       DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
    169          DO ji = 1, jpi 
    170             DO jk = mikt(ji,jj)+1, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
    171                avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) 
    172                avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) 
     168      DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
     169         DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
     170            DO ji = 1, jpi 
     171               avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk) 
     172               avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk) 
    173173            END DO 
    174174         END DO 
    175175      END DO 
    176176       
    177       DO jj = 2, jpjm1 
    178          DO ji = fs_2, fs_jpim1  ! vector opt. 
    179             DO jk = mikt(ji,jj)+1, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
    180                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) 
    181                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) 
     177      DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
     178         DO jj = 2, jpjm1 
     179            DO ji = fs_2, fs_jpim1  ! vector opt. 
     180               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) 
     181               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) 
    182182            END DO 
    183183         END DO 
     
    457457         ztpc = 0.e0 
    458458         zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:) 
    459          DO jj = 1, jpj 
    460             DO ji = 1, jpi 
    461                DO jk= mikt(ji,jj)+1, jpkm1 
    462                   ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 
     459         DO jk= 2, jpkm1 
     460            DO jj = 1, jpj 
     461               DO ji = 1, jpi 
     462                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
    463463               END DO 
    464464            END DO 
     
    473473         zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 )    
    474474         zkz(:,:) = 0.e0 
    475          DO jj = 1, jpj 
    476             DO ji = 1, jpi 
    477                DO jk = mikt(ji,jj)+1, jpkm1 
    478                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) 
     475         DO jk = 2, jpkm1 
     476            DO jj = 1, jpj 
     477               DO ji = 1, jpi 
     478                  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) 
    479479               END DO 
    480480            END DO 
     
    498498         WRITE(numout,*) '          Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) ) 
    499499 
    500          DO jj = 1, jpj 
    501             DO ji = 1, jpi 
    502                DO jk = mikt(ji,jj)+1, jpkm1 
    503                   zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. )   !kz max = 300 cm2/s 
     500         DO jk = 2, jpkm1 
     501            DO jj = 1, jpj 
     502               DO ji = 1, jpi 
     503                  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 
    504504               END DO 
    505505            END DO 
     
    510510            DO jj = 1, jpj 
    511511               DO ji = 1, jpi 
    512                   ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 
     512                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
    513513               END DO 
    514514            END DO 
     
    519519         DO jk = 1, jpk 
    520520            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zav_tide(:,:,jk)     * tmask_i(:,:) )   & 
    521                &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) ) 
     521               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 
    522522            ztpc = 1.E50 
    523523            DO jj = 1, jpj 
     
    540540            END DO 
    541541            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zkz(:,:)     * tmask_i(:,:) )   & 
    542                &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) ) 
     542               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 
    543543            WRITE(numout,*) '                jk= ', jk,'   ', ze_z * 1.e4,' cm2/s' 
    544544         END DO 
     
    546546            zkz(:,:) = az_tmx(:,:,jk) /rn_n2min 
    547547            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zkz(:,:)     * tmask_i(:,:) )   & 
    548                &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) ) 
     548               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 
    549549            WRITE(numout,*)  
    550550            WRITE(numout,*) '          N2 min - jk= ', jk,'   ', ze_z * 1.e4,' cm2/s min= ',MINVAL(zkz)*1.e4,   & 
Note: See TracChangeset for help on using the changeset viewer.