Changeset 2397


Ignore:
Timestamp:
2010-11-16T14:03:39+01:00 (10 years ago)
Author:
cbricaud
Message:

use stress module variable taum and simplifications for mask variables

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r2395 r2397  
    55   !!                 turbulent closure parameterization 
    66   !!====================================================================== 
    7    !! History :   3.0  ! 2009-09 (G. Reffray)  Original code 
    8    !!             3.3  ! 2010-10 (C. Bricaud)  add in the reference 
     7   !! History :   3.0  !  2009-09 (G. Reffray)  Original code 
     8   !!             3.3  !  2010-10  (C. Bricaud)  Add in the reference 
    99   !!---------------------------------------------------------------------- 
    1010#if defined key_zdfgls   ||   defined key_esopa 
     
    3434 
    3535   PUBLIC   zdf_gls        ! routine called in step module 
    36    PUBLIC   zdf_gls_init   ! routine called in step module 
     36   PUBLIC   zdf_gls_init   ! routine called in opa module 
    3737   PUBLIC   gls_rst        ! routine called in step module 
    3838 
     
    4646   !                                         !!! ** Namelist  namzdf_gls  ** 
    4747   LOGICAL  ::   ln_crban      = .FALSE.      ! =T use Craig and Banner scheme 
    48    LOGICAL  ::   ln_length_lim = .FALSE.      ! use limit on the dissipation rate under stable stratif. (Galperin et al., 1988) 
    49    LOGICAL  ::   ln_sigpsi     = .FALSE.      ! Activate Burchard (2003) modification for k-eps closure AND wave breaking mixing 
    50    REAL(wp) ::   rn_epsmin     = 1.e-12_wp    ! minimum value of dissipation (m2/s3) 
    51    REAL(wp) ::   rn_emin       = 1.e-6_wp     ! minimum value of TKE (m2/s2) 
     48   LOGICAL  ::   ln_length_lim = .FALSE.      ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988) 
     49   LOGICAL  ::   ln_sigpsi     = .FALSE.      ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing 
    5250   INTEGER  ::   nn_tkebc_surf = 0            ! TKE surface boundary condition (=0/1) 
    5351   INTEGER  ::   nn_tkebc_bot  = 0            ! TKE bottom boundary condition (=0/1) 
     
    5755   INTEGER  ::   nn_clos       = 0            ! closure 0/1/2/3 MY82/k-eps/k-w/gen 
    5856   REAL(wp) ::   rn_clim_galp  = 0.53_wp      ! Holt 2008 value for k-eps: 0.267 
    59    REAL(wp) ::   hsro          = 0.003_wp     ! Minimum surface roughness 
     57   REAL(wp) ::   rn_epsmin     = 1.e-12_wp    ! minimum value of dissipation (m2/s3) 
     58   REAL(wp) ::   rn_emin       = 1.e-6_wp     ! minimum value of TKE (m2/s2) 
    6059   REAL(wp) ::   rn_charn      = 2.e+5_wp     ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 
    6160   REAL(wp) ::   rn_crban      = 100._wp      ! Craig and Banner constant for surface breaking waves mixing 
    6261 
    63    REAL(wp) ::   rcm_sf        = 0.73_wp        ! Shear free turbulence parameters 
    64    REAL(wp) ::   ra_sf         = -2.0_wp        ! Must be negative -2 < ra_sf < -1  
    65    REAL(wp) ::   rl_sf         =  0.2_wp        ! 0 <rl_sf<vkarmn     
    66    REAL(wp) ::   rghmin        = -0.28 
    67    REAL(wp) ::   rgh0          = 0.0329 
    68    REAL(wp) ::   rghcri        = 0.03 
     62   REAL(wp) ::   hsro          =  0.003_wp    ! Minimum surface roughness 
     63   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters 
     64   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1  
     65   REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn     
     66   REAL(wp) ::   rghmin        = -0.28_wp 
     67   REAL(wp) ::   rgh0          =  0.0329_wp 
     68   REAL(wp) ::   rghcri        =  0.03_wp 
    6969   REAL(wp) ::   ra1           =  0.92_wp 
    7070   REAL(wp) ::   ra2           =  0.74_wp 
     
    8888   REAL(wp) ::   rm7           =  0.0_wp 
    8989   REAL(wp) ::   rm8           =  0.318_wp 
    90    REAL(wp) ::   rc02 
    91    REAL(wp) ::   rc02r 
    92    REAL(wp) ::   rc03 
    93    REAL(wp) ::   rc04 
    94    REAL(wp) ::   rc03_sqrt2_galp 
    95    REAL(wp) ::   rsbc_mb 
    96    REAL(wp) ::   rsbc_std 
    97    REAL(wp) ::   rsbc_tke1 
    98    REAL(wp) ::   rsbc_tke2 
    99    REAL(wp) ::   rsbc_tke3 
    100    REAL(wp) ::   rsbc_psi1 
    101    REAL(wp) ::   rsbc_psi2 
    102    REAL(wp) ::   rsbc_psi3 
    103    REAL(wp) ::   rsbc_zs 
    104    REAL(wp) ::   rfact_tke, rfact_psi 
    105    REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff 
    106    REAL(wp) ::   rs0, rs1, rs2, rs4, rs5, rs6 
    107    REAL(wp) ::   rd0, rd1, rd2, rd3, rd4, rd5 
    108    REAL(wp) ::   rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0 
    109    REAL(wp) ::   rpsi3m, rpsi3p, rpp, rmm, rnn 
     90    
     91   REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters 
     92   REAL(wp) ::   rc03_sqrt2_galp                                  !     -           -           -        - 
     93   REAL(wp) ::   rsbc_tke1, rsbc_tke2, rsbc_tke3, rfact_tke       !     -           -           -        - 
     94   REAL(wp) ::   rsbc_psi1, rsbc_psi2, rsbc_psi3, rfact_psi       !     -           -           -        - 
     95   REAL(wp) ::   rsbc_mb  , rsbc_std , rsbc_zs                    !     -           -           -        - 
     96   REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        - 
     97   REAL(wp) ::   rs0, rs1, rs2, rs4, rs5, rs6                     !     -           -           -        - 
     98   REAL(wp) ::   rd0, rd1, rd2, rd3, rd4, rd5                     !     -           -           -        - 
     99   REAL(wp) ::   rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0  !     -           -           -        - 
     100   REAL(wp) ::   rpsi3m, rpsi3p, rpp, rmm, rnn                    !     -           -           -        - 
    110101 
    111102   !! * Substitutions 
     
    124115      !! 
    125116      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity 
    126       !!      coefficients using a 2.5 turbulent closure scheme. 
     117      !!              coefficients using the GLS turbulent closure scheme. 
    127118      !!---------------------------------------------------------------------- 
    128119      USE oce,     z_elem_a  =>   ua   ! use ua as workspace 
     
    133124      INTEGER, INTENT(in) ::   kt ! ocean time step 
    134125      INTEGER  ::   ji, jj, jk, ibot, ibotm1, dir  ! dummy loop arguments 
    135       ! 
    136       REAL(wp) ::   zesh2, zsigpsi                 ! temporary scalars 
    137       REAL(wp) ::   ztx2, zty2, zup, zdown, zcof   !    -          
    138       REAL(wp) ::   zratio, zrn2, zflxb, sh        !    -         - 
    139       REAL(wp) ::   prod, buoy, diss, zdiss, sm    !    -         - 
    140       REAL(wp) ::   gh, gm, shr, dif, zsqen, zav   !    -         - 
    141       REAL(wp), DIMENSION(jpi,jpj)     ::   zdep   ! 
    142       REAL(wp), DIMENSION(jpi,jpj)     ::   zflxs  ! Turbulence fluxed induced by internal waves  
    143       REAL(wp), DIMENSION(jpi,jpj)     ::   zhsro  ! Surface roughness (surface waves) 
    144       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eb     ! tke at time before 
    145       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   mxlb   ! mixing length at time before 
    146       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   shear  ! vertical shear 
    147       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eps    ! dissipation rate 
    148       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwall_psi ! Wall function for psi schmidt number in the wb case  
    149                                                       ! (if ln_sigpsi.AND.ln_crban) 
     126      REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1, zex2   ! local scalars 
     127      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      -  
     128      REAL(wp) ::   zratio, zrn2, zflxb, sh             !   -      - 
     129      REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      - 
     130      REAL(wp) ::   gh, gm, shr, dif, zsqen, zav        !   -      - 
     131      REAL(wp), DIMENSION(jpi,jpj)     ::   zdep        ! 
     132      REAL(wp), DIMENSION(jpi,jpj)     ::   zflxs       ! Turbulence fluxed induced by internal waves  
     133      REAL(wp), DIMENSION(jpi,jpj)     ::   zhsro       ! Surface roughness (surface waves) 
     134      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eb          ! tke at time before 
     135      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   mxlb        ! mixing length at time before 
     136      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   shear       ! vertical shear 
     137      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eps         ! dissipation rate 
     138      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi.AND.ln_crban=T) 
    150139      !!-------------------------------------------------------------------- 
    151140 
    152141      ! Preliminary computing 
    153142 
    154       ustars2 = 0. ; ustarb2 = 0. ; psi  = 0. ; zwall_psi = 0. 
     143      ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp 
    155144 
    156145      ! Compute wind stress at T-points 
     
    162151          ! wind stress 
    163152          ! squared surface velocity scale 
    164           ztx2 = 2. * (utau(ji-1,jj  )*umask(ji-1,jj,1) + utau(ji,jj)*umask(ji,jj,1)) / & 
    165             &                  MAX(1., umask(ji-1,jj,1) +             umask(ji,jj,1)) 
    166           zty2 = 2. * (vtau(ji  ,jj-1)*vmask(ji,jj-1,1) + vtau(ji,jj)*vmask(ji,jj,1)) / & 
    167             &                  MAX(1., vmask(ji,jj-1,1) +             vmask(ji,jj,1)) 
    168           ustars2(ji,jj) = rsbc_tke2 * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 
     153          ustars2(ji,jj) = rau0r * taum(ji,jj) * tmask(ji,jj,1) 
    169154          ! 
    170155          ! bottom friction  
    171           ztx2 = 2. * (wbotu(ji-1,jj  )*umask(ji-1,jj,1) + wbotu(ji,jj)*umask(ji,jj,1)) / & 
    172             &                   MAX(1., umask(ji-1,jj,1) +              umask(ji,jj,1)) 
    173           zty2 = 2. * (wbotv(ji  ,jj-1)*vmask(ji,jj-1,1) + wbotv(ji,jj)*vmask(ji,jj,1)) / & 
    174             &                   MAX(1., vmask(ji,jj-1,1) +              vmask(ji,jj,1)) 
    175           ustarb2(ji,jj) = 0.5 * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 
    176         ENDDO 
    177       ENDDO   
     156          ztx2 = ( wbotu(ji-1,jj)*umask(ji-1,jj,1) + wbotu(ji,jj)*umask(ji,jj,1) )  & 
     157            &  / MAX( 1.,         umask(ji-1,jj,1) +              umask(ji,jj,1) ) 
     158          zty2 = ( wbotv(ji,jj-1)*vmask(ji,jj-1,1) + wbotv(ji,jj)*vmask(ji,jj,1) )  & 
     159            &  / MAX( 1.,         vmask(ji,jj-1,1) +              vmask(ji,jj,1) ) 
     160          ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 
     161        END DO 
     162      END DO   
    178163 
    179164      ! In case of breaking surface waves mixing, 
    180165      ! Compute surface roughness length according to Charnock formula: 
    181       IF (ln_crban) THEN 
    182         zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro) 
    183       ELSE 
    184         zhsro(:,:) = hsro 
     166      IF( ln_crban ) THEN   ;   zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro) 
     167      ELSE                  ;   zhsro(:,:) = hsro 
    185168      ENDIF 
    186169 
     
    198181                  &                            *    fse3vw_b(ji,jj,jk) ) 
    199182               eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk) 
    200             ENDDO 
    201          ENDDO 
    202       ENDDO 
     183            END DO 
     184         END DO 
     185      END DO 
    203186      ! 
    204187      ! Lateral boundary conditions (avmu,avmv) (sign unchanged) 
    205       CALL lbc_lnk( avmu, 'U', 1. )   ;    CALL lbc_lnk( avmv, 'V', 1. ) 
     188      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. ) 
    206189 
    207190      ! Save tke at before time step 
     
    209192      mxlb(:,:,:) = mxln(:,:,:) 
    210193 
    211       IF ( nn_clos .EQ. 0 ) THEN ! Mellor-Yamada 
     194      IF( nn_clos == 0 ) THEN    ! Mellor-Yamada 
    212195         DO jk = 2, jpkm1 
    213196            DO jj = 2, jpjm1  
    214197               DO ji = fs_2, fs_jpim1   ! vector opt. 
    215198                  zup   = mxln(ji,jj,jk) * fsdepw(ji,jj,mbathy(ji,jj)) 
    216                   zdown = vkarmn * fsdepw(ji,jj,jk) * (-fsdepw(ji,jj,jk)+fsdepw(ji,jj,mbathy(ji,jj))) 
    217                   zwall (ji,jj,jk) = ( 1. + re2 * ( zup / MAX( zdown, rsmall ) )**2. ) * tmask(ji,jj,jk) 
    218                ENDDO 
    219             ENDDO 
    220          ENDDO 
     199                  zdown = vkarmn * fsdepw(ji,jj,jk) * ( -fsdepw(ji,jj,jk) + fsdepw(ji,jj,mbathy(ji,jj)) ) 
     200                  zcoef = ( zup / MAX( zdown, rsmall ) ) 
     201                  zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk) 
     202               END DO 
     203            END DO 
     204         END DO 
    221205      ENDIF 
    222206 
     
    248232               diss = eps(ji,jj,jk) 
    249233               ! 
    250                dir = 0.5 + sign(0.5,shear(ji,jj,jk)+buoy)   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 
    251                ! 
    252                zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1-dir)*shear(ji,jj,jk)          ! production term 
    253                zdiss = dir*(diss/en(ji,jj,jk))   +(1-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 
     234               dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 
     235               ! 
     236               zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1._wp-dir)*shear(ji,jj,jk)          ! production term 
     237               zdiss = dir*(diss/en(ji,jj,jk))   +(1._wp-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 
    254238               ! 
    255239               ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0 
     
    257241               ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries. 
    258242               ! Otherwise, this should be rsc_psi/rsc_psi0 
    259                IF (ln_sigpsi) THEN 
    260                  zsigpsi = MIN(1., zesh2/eps(ji,jj,jk)) ! 0. <= zsigpsi <= 1. 
    261                  zwall_psi(ji,jj,jk) = rsc_psi / (zsigpsi * rsc_psi + (1.-zsigpsi) * rsc_psi0 / MAX(zwall(ji,jj,jk),1.)) 
     243               IF( ln_sigpsi ) THEN 
     244                  zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) )    ! 0. <= zsigpsi <= 1. 
     245                  zwall_psi(ji,jj,jk) = rsc_psi / (  zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp )  ) 
    262246               ELSE 
    263                  zwall_psi(ji,jj,jk) = 1. 
     247                  zwall_psi(ji,jj,jk) = 1._wp 
    264248               ENDIF 
    265249               ! 
     
    276260               ! 
    277261               ! diagonal 
    278                z_elem_b(ji,jj,jk) = 1. - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  & 
    279                   &                    + rdt * zdiss * tmask(ji,jj,jk)  
     262               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  & 
     263                  &                       + rdt * zdiss * tmask(ji,jj,jk)  
    280264               ! 
    281265               ! right hand side in en 
     
    285269      END DO 
    286270      ! 
    287       z_elem_b(:,:,jpk) = 1. 
     271      z_elem_b(:,:,jpk) = 1._wp 
    288272      ! 
    289273      ! Set surface condition on zwall_psi (1 at the bottom) 
    290       IF (ln_sigpsi) THEN 
    291         DO jj = 2, jpjm1 
    292           DO ji = fs_2, fs_jpim1   ! vector opt. 
    293             zwall_psi(ji,jj,1) = rsc_psi/rsc_psi0 
    294           END DO 
    295         END DO 
     274      IF( ln_sigpsi ) THEN 
     275         zcoef = rsc_psi / rsc_psi0 
     276         DO jj = 2, jpjm1 
     277            DO ji = fs_2, fs_jpim1   ! vector opt. 
     278               zwall_psi(ji,jj,1) = zcoef 
     279            END DO 
     280         END DO 
    296281      ENDIF 
    297282 
     
    302287      ! 
    303288      CASE ( 0 )             ! Dirichlet case 
    304       ! 
    305       IF (ln_crban) THEN     ! Wave induced mixing case 
    306       !                      ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 
    307       !                      ! balance between the production and the dissipation terms including the wave effect 
    308       ! 
    309       en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 
    310       z_elem_a(:,:,1) = en(:,:,1) 
    311       z_elem_c(:,:,1) = 0. 
    312       z_elem_b(:,:,1) = 1. 
    313       !  
    314       ! one level below 
    315       en(:,:,2) = MAX( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 
    316       z_elem_a(:,:,2) = 0. 
    317       z_elem_c(:,:,2) = 0. 
    318       z_elem_b(:,:,2) = 1. 
    319       ! 
    320       ELSE                   ! No wave induced mixing case 
    321       !                      ! en(1) = u*^2/C0^2  &  l(1)  = K*zs 
    322       !                      ! balance between the production and the dissipation terms 
    323       ! 
    324       en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 
    325       z_elem_a(:,:,1) = en(:,:,1)  
    326       z_elem_c(:,:,1) = 0. 
    327       z_elem_b(:,:,1) = 1. 
    328       ! 
    329       ! one level below 
    330       en(:,:,2) = MAX( rc02r * ustars2(:,:), rn_emin ) 
    331       z_elem_a(:,:,2) = 0. 
    332       z_elem_c(:,:,2) = 0. 
    333       z_elem_b(:,:,2) = 1. 
    334       ! 
    335       ENDIF 
    336       ! 
     289         ! 
     290         IF (ln_crban) THEN     ! Wave induced mixing case 
     291            !                      ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 
     292            !                      ! balance between the production and the dissipation terms including the wave effect 
     293            en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 
     294            z_elem_a(:,:,1) = en(:,:,1) 
     295            z_elem_c(:,:,1) = 0._wp 
     296            z_elem_b(:,:,1) = 1._wp 
     297            !  
     298            ! one level below 
     299            en(:,:,2) = MAX( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 
     300            z_elem_a(:,:,2) = 0._wp 
     301            z_elem_c(:,:,2) = 0._wp 
     302            z_elem_b(:,:,2) = 1._wp 
     303            ! 
     304         ELSE                   ! No wave induced mixing case 
     305            !                      ! en(1) = u*^2/C0^2  &  l(1)  = K*zs 
     306            !                      ! balance between the production and the dissipation terms 
     307            en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 
     308            z_elem_a(:,:,1) = en(:,:,1)  
     309            z_elem_c(:,:,1) = 0._wp 
     310            z_elem_b(:,:,1) = 1._wp 
     311            ! 
     312            ! one level below 
     313            en(:,:,2) = MAX( rc02r * ustars2(:,:), rn_emin ) 
     314            z_elem_a(:,:,2) = 0._wp 
     315            z_elem_c(:,:,2) = 0._wp 
     316            z_elem_b(:,:,2) = 1._wp 
     317            ! 
     318         ENDIF 
     319         ! 
    337320      CASE ( 1 )             ! Neumann boundary condition on d(e)/dz 
    338       ! 
    339       IF (ln_crban) THEN ! Shear free case: d(e)/dz= Fw 
    340       ! 
    341       ! Dirichlet conditions at k=1 (Cosmetic) 
    342       en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 
    343       z_elem_a(:,:,1) = en(:,:,1) 
    344       z_elem_c(:,:,1) = 0. 
    345       z_elem_b(:,:,1) = 1. 
    346       ! at k=2, set de/dz=Fw 
    347       z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    348       z_elem_a(:,:,2) = 0.            
    349       zflxs(:,:) = rsbc_tke3 * ustars2(:,:)**1.5 * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5*ra_sf) 
    350       en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 
    351       ! 
    352       ELSE                   ! No wave induced mixing case: d(e)/dz=0. 
    353       ! 
    354       ! Dirichlet conditions at k=1 (Cosmetic) 
    355       en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 
    356       z_elem_a(:,:,1) = en(:,:,1) 
    357       z_elem_c(:,:,1) = 0. 
    358       z_elem_b(:,:,1) = 1. 
    359       ! at k=2 set de/dz=0.: 
    360       z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    361       z_elem_a(:,:,2) = 0. 
    362       ! 
    363       ENDIF 
    364       ! 
     321         ! 
     322         IF (ln_crban) THEN ! Shear free case: d(e)/dz= Fw 
     323            ! 
     324            ! Dirichlet conditions at k=1 (Cosmetic) 
     325            en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 
     326            z_elem_a(:,:,1) = en(:,:,1) 
     327            z_elem_c(:,:,1) = 0._wp 
     328            z_elem_b(:,:,1) = 1._wp 
     329            ! at k=2, set de/dz=Fw 
     330            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
     331            z_elem_a(:,:,2) = 0._wp         
     332            zflxs(:,:) = rsbc_tke3 * ustars2(:,:)**1.5_wp * ( (zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 
     333            en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
     334            ! 
     335         ELSE                   ! No wave induced mixing case: d(e)/dz=0. 
     336            ! 
     337            ! Dirichlet conditions at k=1 (Cosmetic) 
     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            ! at k=2 set de/dz=0.: 
     343            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
     344            z_elem_a(:,:,2) = 0._wp 
     345            ! 
     346         ENDIF 
     347         ! 
    365348      END SELECT 
    366349 
     
    371354      ! 
    372355      CASE ( 0 )             ! Dirichlet  
    373       !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin 
    374       !                      ! Balance between the production and the dissipation terms 
    375       !                       
     356         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin 
     357         !                      ! Balance between the production and the dissipation terms 
    376358!CDIR NOVERRCHK 
    377       DO jj = 2, jpjm1 
     359         DO jj = 2, jpjm1 
    378360!CDIR NOVERRCHK 
    379          DO ji = fs_2, fs_jpim1   ! vector opt. 
    380             ibot   = mbathy(ji,jj) 
    381             ibotm1 = ibot-1 
    382             ! 
    383             ! Bottom level Dirichlet condition: 
    384             z_elem_a(ji,jj,ibot  ) = 0. 
    385             z_elem_c(ji,jj,ibot  ) = 0. 
    386             z_elem_b(ji,jj,ibot  ) = 1. 
    387             en(ji,jj,ibot  ) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
    388             ! 
    389             ! Just above last level, Dirichlet condition again 
    390             z_elem_a(ji,jj,ibotm1) = 0. 
    391             z_elem_c(ji,jj,ibotm1) = 0. 
    392             z_elem_b(ji,jj,ibotm1) = 1. 
    393             en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin )  
    394          END DO 
    395       END DO 
    396       ! 
     361            DO ji = fs_2, fs_jpim1   ! vector opt. 
     362               ibot   = mbathy(ji,jj) 
     363               ibotm1 = ibot-1 
     364               ! 
     365               ! Bottom level Dirichlet condition: 
     366               z_elem_a(ji,jj,ibot  ) = 0._wp 
     367               z_elem_c(ji,jj,ibot  ) = 0._wp 
     368               z_elem_b(ji,jj,ibot  ) = 1._wp 
     369               en(ji,jj,ibot  ) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
     370               ! 
     371               ! Just above last level, Dirichlet condition again 
     372               z_elem_a(ji,jj,ibotm1) = 0._wp 
     373               z_elem_c(ji,jj,ibotm1) = 0._wp 
     374               z_elem_b(ji,jj,ibotm1) = 1._wp 
     375               en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin )  
     376            END DO 
     377         END DO 
     378         ! 
    397379      CASE ( 1 )             ! Neumman boundary condition 
    398       !                       
     380         !                       
    399381!CDIR NOVERRCHK 
    400       DO jj = 2, jpjm1 
     382         DO jj = 2, jpjm1 
    401383!CDIR NOVERRCHK 
    402          DO ji = fs_2, fs_jpim1   ! vector opt. 
    403             ibot   = mbathy(ji,jj) 
    404             ibotm1 = ibot-1 
    405             ! 
    406             ! Bottom level Dirichlet condition: 
    407             z_elem_a(ji,jj,ibot) = 0. 
    408             z_elem_c(ji,jj,ibot) = 0. 
    409             z_elem_b(ji,jj,ibot) = 1. 
    410             en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
    411             ! 
    412             ! Just above last level: Neumann condition 
    413             z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b 
    414             z_elem_c(ji,jj,ibotm1) = 0. 
    415          END DO 
    416       END DO 
    417       ! 
     384            DO ji = fs_2, fs_jpim1   ! vector opt. 
     385               ibot   = mbathy(ji,jj) 
     386               ibotm1 = ibot-1 
     387               ! 
     388               ! Bottom level Dirichlet condition: 
     389               z_elem_a(ji,jj,ibot) = 0._wp 
     390               z_elem_c(ji,jj,ibot) = 0._wp 
     391               z_elem_b(ji,jj,ibot) = 1._wp 
     392               en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
     393               ! 
     394               ! Just above last level: Neumann condition 
     395               z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1)  ! Remove z_elem_c from z_elem_b 
     396               z_elem_c(ji,jj,ibotm1) = 0._wp 
     397            END DO 
     398         END DO 
     399         ! 
    418400      END SELECT 
    419401 
     
    442424         END DO 
    443425      END DO 
    444       ! 
    445       ! set the minimum value of tke  
     426      !                                            ! set the minimum value of tke  
    446427      en(:,:,:) = MAX( en(:,:,:), rn_emin ) 
    447428       
     
    455436      ! 
    456437      CASE( 0 )               ! k-kl  (Mellor-Yamada) 
    457       ! 
    458       DO jk = 2, jpkm1 
    459          DO jj = 2, jpjm1 
    460             DO ji = fs_2, fs_jpim1   ! vector opt. 
    461                psi(ji,jj,jk)  = en(ji,jj,jk) * mxln(ji,jj,jk) 
    462             ENDDO 
    463          ENDDO 
    464       ENDDO 
    465       ! 
     438         DO jk = 2, jpkm1 
     439            DO jj = 2, jpjm1 
     440               DO ji = fs_2, fs_jpim1   ! vector opt. 
     441                  psi(ji,jj,jk)  = en(ji,jj,jk) * mxln(ji,jj,jk) 
     442               END DO 
     443            END DO 
     444         END DO 
     445         ! 
    466446      CASE( 1 )               ! k-eps 
    467       ! 
    468       DO jk = 2, jpkm1 
    469          DO jj = 2, jpjm1 
    470             DO ji = fs_2, fs_jpim1   ! vector opt. 
    471                psi(ji,jj,jk)  = eps(ji,jj,jk) 
    472             ENDDO 
    473          ENDDO 
    474       ENDDO 
    475       ! 
     447         DO jk = 2, jpkm1 
     448            DO jj = 2, jpjm1 
     449               DO ji = fs_2, fs_jpim1   ! vector opt. 
     450                  psi(ji,jj,jk)  = eps(ji,jj,jk) 
     451               END DO 
     452            END DO 
     453         END DO 
     454         ! 
    476455      CASE( 2 )               ! k-w 
    477       ! 
    478       DO jk = 2, jpkm1 
    479          DO jj = 2, jpjm1 
    480             DO ji = fs_2, fs_jpim1   ! vector opt. 
    481                psi(ji,jj,jk)  = SQRT(en(ji,jj,jk)) / (rc0 * mxln(ji,jj,jk)) 
    482             ENDDO 
    483          ENDDO 
    484       ENDDO 
    485       ! 
    486       CASE( 3 )               ! gen 
    487       ! 
    488       DO jk = 2, jpkm1 
    489          DO jj = 2, jpjm1 
    490             DO ji = fs_2, fs_jpim1   ! vector opt. 
    491                psi(ji,jj,jk)  = rc02 * en(ji,jj,jk) * mxln(ji,jj,jk)**rnn  
    492             ENDDO 
    493          ENDDO 
    494       ENDDO 
    495       ! 
     456         DO jk = 2, jpkm1 
     457            DO jj = 2, jpjm1 
     458               DO ji = fs_2, fs_jpim1   ! vector opt. 
     459                  psi(ji,jj,jk)  = SQRT( en(ji,jj,jk) ) / ( rc0 * mxln(ji,jj,jk) ) 
     460               END DO 
     461            END DO 
     462         END DO 
     463         ! 
     464      CASE( 3 )               ! generic 
     465         DO jk = 2, jpkm1 
     466            DO jj = 2, jpjm1 
     467               DO ji = fs_2, fs_jpim1   ! vector opt. 
     468                  psi(ji,jj,jk)  = rc02 * en(ji,jj,jk) * mxln(ji,jj,jk)**rnn  
     469               END DO 
     470            END DO 
     471         END DO 
     472         ! 
    496473      END SELECT 
    497474      ! 
     
    511488               ! 
    512489               ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable) 
    513                dir = 0.5 + sign(0.5,rn2(ji,jj,jk)) 
    514                ! 
    515                rpsi3 = dir*rpsi3m+(1-dir)*rpsi3p 
     490               dir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) ) 
     491               ! 
     492               rpsi3 = dir * rpsi3m + ( 1._wp - dir ) * rpsi3p 
    516493               ! 
    517494               ! shear prod. - stratif. destruction 
     
    519496               ! 
    520497               ! stratif. destruction 
    521                buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk)) 
     498               buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk) ) 
    522499               ! 
    523500               ! shear prod. - stratif. destruction 
    524501               diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) 
    525502               ! 
    526                dir = 0.5 + sign(0.5,prod+buoy)   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 
    527                ! 
    528                zesh2 = dir*(prod+buoy)          +(1-dir)*prod                      ! production term 
    529                zdiss = dir*(diss/psi(ji,jj,jk)) +(1-dir)*(diss-buoy)/psi(ji,jj,jk) ! dissipation term 
     503               dir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 
     504               ! 
     505               zesh2 = dir * ( prod + buoy )          + (1._wp - dir ) * prod                        ! production term 
     506               zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 
    530507               !                                                         
    531508               ! building the matrix 
     
    538515                  &                      / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) ) 
    539516               ! diagonal 
    540                z_elem_b(ji,jj,jk) = 1. - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  & 
    541                   &                    + rdt * zdiss * tmask(ji,jj,jk) 
     517               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  & 
     518                  &                       + rdt * zdiss * tmask(ji,jj,jk) 
    542519               ! 
    543520               ! right hand side in psi 
     
    547524      END DO 
    548525      ! 
    549       z_elem_b(:,:,jpk) = 1. 
     526      z_elem_b(:,:,jpk) = 1._wp 
    550527 
    551528      ! Surface boundary condition on psi 
     
    555532      ! 
    556533      CASE ( 0 )             ! Dirichlet boundary conditions 
    557       ! 
    558       IF (ln_crban) THEN     ! Wave induced mixing case 
    559       !                      ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 
    560       !                      ! balance between the production and the dissipation terms including the wave effect 
    561       ! 
    562       zdep(:,:) = rl_sf * zhsro(:,:) 
    563       psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    564       z_elem_a(:,:,1) = psi(:,:,1) 
    565       z_elem_c(:,:,1) = 0. 
    566       z_elem_b(:,:,1) = 1. 
    567       ! 
    568       ! one level below 
    569       zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**(rmm*ra_sf+rnn) )  & 
    570         &                       / zhsro(:,:)**(rmm*ra_sf) 
    571       psi (:,:,2) = rsbc_psi1 * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 
    572       z_elem_a(:,:,2) = 0. 
    573       z_elem_c(:,:,2) = 0. 
    574       z_elem_b(:,:,2) = 1. 
    575       !  
    576       ELSE                   ! No wave induced mixing case 
    577       !                      ! en(1) = u*^2/C0^2  &  l(1)  = K*zs 
    578       !                      ! balance between the production and the dissipation terms 
    579       ! 
    580       zdep(:,:) = vkarmn * zhsro(:,:) 
    581       psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    582       z_elem_a(:,:,1) = psi(:,:,1) 
    583       z_elem_c(:,:,1) = 0. 
    584       z_elem_b(:,:,1) = 1. 
    585       ! 
    586       ! one level below 
    587       zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) ) 
    588       psi (:,:,2) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    589       z_elem_a(:,:,2) = 0. 
    590       z_elem_c(:,:,2) = 0. 
    591       z_elem_b(:,:,2) = 1. 
    592       ! 
    593       ENDIF 
    594       ! 
     534         ! 
     535         IF( ln_crban ) THEN       ! Wave induced mixing case 
     536            !                      ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 
     537            !                      ! balance between the production and the dissipation terms including the wave effect 
     538            zdep(:,:) = rl_sf * zhsro(:,:) 
     539            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     540            z_elem_a(:,:,1) = psi(:,:,1) 
     541            z_elem_c(:,:,1) = 0._wp 
     542            z_elem_b(:,:,1) = 1._wp 
     543            ! 
     544            ! one level below 
     545            zex1 = (rmm*ra_sf+rnn) 
     546            zex2 = (rmm*ra_sf) 
     547            zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 
     548            psi (:,:,2) = rsbc_psi1 * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 
     549            z_elem_a(:,:,2) = 0._wp 
     550            z_elem_c(:,:,2) = 0._wp 
     551            z_elem_b(:,:,2) = 1._wp 
     552            !  
     553         ELSE                   ! No wave induced mixing case 
     554            !                      ! en(1) = u*^2/C0^2  &  l(1)  = K*zs 
     555            !                      ! balance between the production and the dissipation terms 
     556            ! 
     557            zdep(:,:) = vkarmn * zhsro(:,:) 
     558            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     559            z_elem_a(:,:,1) = psi(:,:,1) 
     560            z_elem_c(:,:,1) = 0._wp 
     561            z_elem_b(:,:,1) = 1._wp 
     562            ! 
     563            ! one level below 
     564            zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) ) 
     565            psi (:,:,2) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     566            z_elem_a(:,:,2) = 0._wp 
     567            z_elem_c(:,:,2) = 0._wp 
     568            z_elem_b(:,:,2) = 1. 
     569            ! 
     570         ENDIF 
     571         ! 
    595572      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz 
    596       ! 
    597       IF (ln_crban) THEN     ! Wave induced mixing case 
    598       ! 
    599       zdep(:,:) = rl_sf * zhsro(:,:) 
    600       psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    601       z_elem_a(:,:,1) = psi(:,:,1) 
    602       z_elem_c(:,:,1) = 0. 
    603       z_elem_b(:,:,1) = 1. 
    604       ! 
    605       ! Neumann condition at k=2 
    606       z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    607       z_elem_a(:,:,2) = 0. 
    608       ! 
    609       ! Set psi vertical flux at the surface: 
    610       zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1.) / zhsro(:,:)**(rmm*ra_sf) 
    611       zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) &  
    612                  & * en(:,:,1)**rmm * zdep          
    613       psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
    614       ! 
     573         ! 
     574         IF( ln_crban ) THEN     ! Wave induced mixing case 
     575            ! 
     576            zdep(:,:) = rl_sf * zhsro(:,:) 
     577            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     578            z_elem_a(:,:,1) = psi(:,:,1) 
     579            z_elem_c(:,:,1) = 0._wp 
     580            z_elem_b(:,:,1) = 1._wp 
     581            ! 
     582            ! Neumann condition at k=2 
     583            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
     584            z_elem_a(:,:,2) = 0._wp 
     585            ! 
     586            ! Set psi vertical flux at the surface: 
     587            zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 
     588            zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) &  
     589               &                  * en(:,:,1)**rmm * zdep          
     590            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
     591            ! 
    615592      ELSE                   ! No wave induced mixing 
    616       ! 
    617       zdep(:,:) = vkarmn * zhsro(:,:) 
    618       psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    619       z_elem_a(:,:,1) = psi(:,:,1) 
    620       z_elem_c(:,:,1) = 0. 
    621       z_elem_b(:,:,1) = 1. 
    622       ! 
    623       ! Neumann condition at k=2 
    624       z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    625       z_elem_a(ji,jj,2) = 0. 
    626       ! 
    627       ! Set psi vertical flux at the surface: 
    628       zdep(:,:) = zhsro(:,:) + fsdept(:,:,1) 
    629       zflxs(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-1.) 
    630       psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
    631       !      
    632       ENDIF 
    633       ! 
     593            ! 
     594            zdep(:,:) = vkarmn * zhsro(:,:) 
     595            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     596            z_elem_a(:,:,1) = psi(:,:,1) 
     597            z_elem_c(:,:,1) = 0._wp 
     598            z_elem_b(:,:,1) = 1._wp 
     599            ! 
     600            ! Neumann condition at k=2 
     601            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
     602            z_elem_a(ji,jj,2) = 0._wp 
     603            ! 
     604            ! Set psi vertical flux at the surface: 
     605            zdep(:,:) = zhsro(:,:) + fsdept(:,:,1) 
     606            zflxs(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-1._wp) 
     607            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
     608            !      
     609         ENDIF 
     610         ! 
    634611      END SELECT 
    635612 
     
    639616      SELECT CASE ( nn_psibc_bot ) 
    640617      ! 
    641       ! 
    642618      CASE ( 0 )             ! Dirichlet  
    643       !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_hbro 
    644       !                      ! Balance between the production and the dissipation terms 
     619         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_hbro 
     620         !                      ! Balance between the production and the dissipation terms 
    645621!CDIR NOVERRCHK 
    646       DO jj = 2, jpjm1 
     622         DO jj = 2, jpjm1 
    647623!CDIR NOVERRCHK 
    648          DO ji = fs_2, fs_jpim1   ! vector opt. 
    649             ibot = mbathy(ji,jj) 
    650             ibotm1 = ibot-1 
    651             zdep(ji,jj) = vkarmn * rn_hbro 
    652             psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 
    653             z_elem_a(ji,jj,ibot) = 0. 
    654             z_elem_c(ji,jj,ibot) = 0. 
    655             z_elem_b(ji,jj,ibot) = 1. 
    656             ! 
    657             ! Just above last level, Dirichlet condition again (GOTM like) 
    658             zdep(ji,jj) = vkarmn * (rn_hbro + fse3t(ji,jj,ibotm1)) 
    659             psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn 
    660             z_elem_a(ji,jj,ibotm1) = 0. 
    661             z_elem_c(ji,jj,ibotm1) = 0. 
    662             z_elem_b(ji,jj,ibotm1) = 1. 
    663          END DO 
    664       END DO 
    665       ! 
    666       ! 
     624            DO ji = fs_2, fs_jpim1   ! vector opt. 
     625               ibot = mbathy(ji,jj) 
     626               ibotm1 = ibot-1 
     627               zdep(ji,jj) = vkarmn * rn_hbro 
     628               psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 
     629               z_elem_a(ji,jj,ibot) = 0._wp 
     630               z_elem_c(ji,jj,ibot) = 0._wp 
     631               z_elem_b(ji,jj,ibot) = 1._wp 
     632               ! 
     633               ! Just above last level, Dirichlet condition again (GOTM like) 
     634               zdep(ji,jj) = vkarmn * (rn_hbro + fse3t(ji,jj,ibotm1)) 
     635               psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn 
     636               z_elem_a(ji,jj,ibotm1) = 0._wp 
     637               z_elem_c(ji,jj,ibotm1) = 0._wp 
     638               z_elem_b(ji,jj,ibotm1) = 1._wp 
     639            END DO 
     640         END DO 
     641         ! 
    667642      CASE ( 1 )             ! Neumman boundary condition 
    668       !                       
     643         !                       
    669644!CDIR NOVERRCHK 
    670       DO jj = 2, jpjm1 
     645         DO jj = 2, jpjm1 
    671646!CDIR NOVERRCHK 
    672          DO ji = fs_2, fs_jpim1   ! vector opt. 
    673             ibot   = mbathy(ji,jj) 
    674             ibotm1 = ibot-1 
    675             ! 
    676             ! Bottom level Dirichlet condition: 
    677             zdep(ji,jj) = vkarmn * rn_hbro 
    678             psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 
    679             ! 
    680             z_elem_a(ji,jj,ibot) = 0. 
    681             z_elem_c(ji,jj,ibot) = 0. 
    682             z_elem_b(ji,jj,ibot) = 1. 
    683             ! 
    684             ! Just above last level: Neumann condition with flux injection 
    685             z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b 
    686             z_elem_c(ji,jj,ibotm1) = 0. 
    687             ! 
    688             ! Set psi vertical flux at the bottom: 
    689             zdep(ji,jj) = rn_hbro + 0.5*fse3t(ji,jj,ibotm1) 
    690             zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) * & 
    691               & (0.5*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1.) 
    692             psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / fse3w(ji,jj,ibotm1) 
    693          END DO 
    694       END DO 
    695       ! 
     647            DO ji = fs_2, fs_jpim1   ! vector opt. 
     648               ibot   = mbathy(ji,jj) 
     649               ibotm1 = ibot-1 
     650               ! 
     651               ! Bottom level Dirichlet condition: 
     652               zdep(ji,jj) = vkarmn * rn_hbro 
     653               psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 
     654               ! 
     655               z_elem_a(ji,jj,ibot) = 0._wp 
     656               z_elem_c(ji,jj,ibot) = 0._wp 
     657               z_elem_b(ji,jj,ibot) = 1._wp 
     658               ! 
     659               ! Just above last level: Neumann condition with flux injection 
     660               z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b 
     661               z_elem_c(ji,jj,ibotm1) = 0. 
     662               ! 
     663               ! Set psi vertical flux at the bottom: 
     664               zdep(ji,jj) = rn_hbro + 0.5_wp*fse3t(ji,jj,ibotm1) 
     665               zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) )  & 
     666                  &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 
     667               psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / fse3w(ji,jj,ibotm1) 
     668            END DO 
     669         END DO 
     670         ! 
    696671      END SELECT 
    697672 
     
    727702      ! 
    728703      CASE( 0 )               ! k-kl  (Mellor-Yamada) 
    729       ! 
    730       DO jk = 1, jpkm1 
    731          DO jj = 2, jpjm1 
    732             DO ji = fs_2, fs_jpim1   ! vector opt. 
    733               eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / psi(ji,jj,jk) 
    734             ENDDO 
    735          ENDDO 
    736       ENDDO 
    737       ! 
     704         DO jk = 1, jpkm1 
     705            DO jj = 2, jpjm1 
     706               DO ji = fs_2, fs_jpim1   ! vector opt. 
     707                  eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / psi(ji,jj,jk) 
     708               END DO 
     709            END DO 
     710         END DO 
     711         ! 
    738712      CASE( 1 )               ! k-eps 
    739       ! 
    740       DO jk = 1, jpkm1 
    741          DO jj = 2, jpjm1 
    742             DO ji = fs_2, fs_jpim1   ! vector opt. 
    743                eps(ji,jj,jk)  = psi(ji,jj,jk) 
    744             ENDDO 
    745          ENDDO 
    746       ENDDO 
    747       ! 
     713         DO jk = 1, jpkm1 
     714            DO jj = 2, jpjm1 
     715               DO ji = fs_2, fs_jpim1   ! vector opt. 
     716                  eps(ji,jj,jk) = psi(ji,jj,jk) 
     717               END DO 
     718            END DO 
     719         END DO 
     720         ! 
    748721      CASE( 2 )               ! k-w 
    749       ! 
    750       DO jk = 1, jpkm1 
    751          DO jj = 2, jpjm1 
    752             DO ji = fs_2, fs_jpim1   ! vector opt. 
    753               eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)  
    754             ENDDO 
    755          ENDDO 
    756       ENDDO 
    757       ! 
    758       CASE( 3 )               ! gen 
    759       ! 
    760       DO jk = 1, jpkm1 
    761          DO jj = 2, jpjm1 
    762             DO ji = fs_2, fs_jpim1   ! vector opt. 
    763               eps(ji,jj,jk) = rc0**(3.+rpp/rnn) * en(ji,jj,jk)**(1.5+rmm/rnn) * psi(ji,jj,jk)**(-1./rnn) 
    764             ENDDO 
    765          ENDDO 
    766       ENDDO 
    767       ! 
     722         DO jk = 1, jpkm1 
     723            DO jj = 2, jpjm1 
     724               DO ji = fs_2, fs_jpim1   ! vector opt. 
     725                  eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)  
     726               END DO 
     727            END DO 
     728         END DO 
     729         ! 
     730      CASE( 3 )               ! generic 
     731         zcoef = rc0**( 3._wp  + rpp/rnn ) 
     732         zex1  =      ( 1.5_wp + rmm/rnn ) 
     733         zex2  = -1._wp / rnn 
     734         DO jk = 1, jpkm1 
     735            DO jj = 2, jpjm1 
     736               DO ji = fs_2, fs_jpim1   ! vector opt. 
     737                  eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 
     738               END DO 
     739            END DO 
     740         END DO 
     741         ! 
    768742      END SELECT 
    769743 
     
    775749               ! limitation 
    776750               eps(ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin ) 
    777                mxln(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / eps(ji,jj,jk) 
     751               mxln(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 
    778752               ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)  
    779753               zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    780                mxln(ji,jj,jk) = MIN( rn_clim_galp*SQRT( 2.*en(ji,jj,jk)/zrn2 ), mxln(ji,jj,jk) ) 
     754               mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) ) 
    781755            END DO 
    782756         END DO 
     
    790764      ! 
    791765      CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions 
    792       ! 
    793       DO jk = 2, jpkm1 
    794          DO jj = 2, jpjm1 
    795             DO ji = fs_2, fs_jpim1   ! vector opt. 
    796                ! zcof =  l²/q² 
    797                zcof = mxlb(ji,jj,jk)**2. / ( 2.*eb(ji,jj,jk) ) 
    798                ! Gh = -N²l²/q² 
    799                gh = - rn2(ji,jj,jk) * zcof 
    800                gh = MIN( gh, rgh0   ) 
    801                gh = MAX( gh, rghmin ) 
    802                ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin) 
    803                sh = ra2*( 1.-6.*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1.-rc3 ) ) ) 
    804                sm = ( rb1**(-1./3.) + ( 18*ra1*ra1 + 9.*ra1*ra2*(1.-rc2) )*sh*gh ) / (1.-9.*ra1*ra2*gh) 
    805                ! 
    806                ! Store stability function in avmu and avmv 
    807                avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 
    808                avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 
    809             END DO 
    810          END DO 
    811       END DO 
    812       ! 
     766         DO jk = 2, jpkm1 
     767            DO jj = 2, jpjm1 
     768               DO ji = fs_2, fs_jpim1   ! vector opt. 
     769                  ! zcof =  l²/q² 
     770                  zcof = mxlb(ji,jj,jk) * mxlb(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) 
     771                  ! Gh = -N²l²/q² 
     772                  gh = - rn2(ji,jj,jk) * zcof 
     773                  gh = MIN( gh, rgh0   ) 
     774                  gh = MAX( gh, rghmin ) 
     775                  ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin) 
     776                  sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) ) 
     777                  sm = ( rb1**(-1._wp/3._wp) + ( 18._wp*ra1*ra1 + 9._wp*ra1*ra2*(1._wp-rc2) )*sh*gh ) / (1._wp-9._wp*ra1*ra2*gh) 
     778                  ! 
     779                  ! Store stability function in avmu and avmv 
     780                  avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 
     781                  avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 
     782               END DO 
     783            END DO 
     784         END DO 
     785         ! 
    813786      CASE ( 2, 3 )               ! Canuto stability functions 
    814       ! 
    815       DO jk = 2, jpkm1 
    816          DO jj = 2, jpjm1 
    817             DO ji = fs_2, fs_jpim1   ! vector opt. 
    818                ! zcof =  l²/q² 
    819                zcof = mxlb(ji,jj,jk)**2. / ( 2.*eb(ji,jj,jk) ) 
    820                ! Gh = -N²l²/q² 
    821                gh = - rn2(ji,jj,jk) * zcof 
    822                gh = MIN( gh, rgh0   ) 
    823                gh = MAX( gh, rghmin ) 
    824                gh = gh*rf6 
    825                ! Gm =  M²l²/q² Shear number 
    826                shr = shear(ji,jj,jk)/MAX(avm(ji,jj,jk), rsmall) 
    827                gm = shr * zcof 
    828                gm = MAX(gm, 1.e-10) 
    829                gm = gm*rf6 
    830                gm = MIN ( (rd0 - rd1*gh + rd3*gh**2.)/(rd2-rd4*gh) , gm ) 
    831                ! Stability functions from Canuto 
    832                rcff = rd0 - rd1*gh +rd2*gm + rd3*gh**2. - rd4*gh*gm + rd5*gm**2. 
    833                sm = (rs0 - rs1*gh + rs2*gm) / rcff 
    834                sh = (rs4 - rs5*gh + rs6*gm) / rcff 
    835                ! 
    836                ! Store stability function in avmu and avmv 
    837                avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 
    838                avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 
    839             END DO 
    840          END DO 
    841       END DO 
    842       ! 
     787         DO jk = 2, jpkm1 
     788            DO jj = 2, jpjm1 
     789               DO ji = fs_2, fs_jpim1   ! vector opt. 
     790                  ! zcof =  l²/q² 
     791                  zcof = mxlb(ji,jj,jk)*mxlb(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) 
     792                  ! Gh = -N²l²/q² 
     793                  gh = - rn2(ji,jj,jk) * zcof 
     794                  gh = MIN( gh, rgh0   ) 
     795                  gh = MAX( gh, rghmin ) 
     796                  gh = gh * rf6 
     797                  ! Gm =  M²l²/q² Shear number 
     798                  shr = shear(ji,jj,jk) / MAX( avm(ji,jj,jk), rsmall ) 
     799                  gm = MAX( shr * zcof , 1.e-10 ) 
     800                  gm = gm * rf6 
     801                  gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm ) 
     802                  ! Stability functions from Canuto 
     803                  rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm 
     804                  sm = (rs0 - rs1*gh + rs2*gm) / rcff 
     805                  sh = (rs4 - rs5*gh + rs6*gm) / rcff 
     806                  ! 
     807                  ! Store stability function in avmu and avmv 
     808                  avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 
     809                  avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 
     810               END DO 
     811            END DO 
     812         END DO 
     813         ! 
    843814      END SELECT 
    844815 
    845816      ! Boundary conditions on stability functions for momentum (Neumann): 
    846817      ! Lines below are useless if GOTM style Dirichlet conditions are used 
     818      zcoef = rcm_sf / SQRT( 2._wp ) 
    847819      DO jj = 2, jpjm1 
    848820         DO ji = fs_2, fs_jpim1   ! vector opt. 
    849             avmv(ji,jj,1) = rcm_sf / SQRT(2.) 
    850          END DO 
    851       END DO 
     821            avmv(ji,jj,1) = zcoef 
     822         END DO 
     823      END DO 
     824      zcoef = rc0 / SQRT( 2._wp ) 
    852825      DO jj = 2, jpjm1 
    853826         DO ji = fs_2, fs_jpim1   ! vector opt. 
    854             ibot=mbathy(ji,jj) 
    855             ibotm1=ibot-1 
    856             avmv(ji,jj,ibot) = rc0 / SQRT(2.) 
     827            ibot   = mbathy(ji,jj) 
     828            avmv(ji,jj,ibot) = zcoef 
    857829         END DO 
    858830      END DO 
     
    863835         DO jj = 2, jpjm1 
    864836            DO ji = fs_2, fs_jpim1   ! vector opt. 
    865                zsqen           = SQRT(2.*en(ji,jj,jk)) * mxln(ji,jj,jk) 
    866                zav             = zsqen * avmu(ji,jj,jk) 
    867                avt(ji,jj,jk)   = MAX(zav,avtb(jk))*tmask(ji,jj,jk) ! apply mask for zdfmxl routine 
    868                zav             = zsqen * avmv(ji,jj,jk) 
    869                avm(ji,jj,jk)   = MAX(zav,avmb(jk)) ! Note that avm is not masked at the surface and the bottom 
    870             END DO 
    871          END DO 
    872       END DO 
    873  
     837               zsqen         = SQRT( 2._wp * en(ji,jj,jk) ) * mxln(ji,jj,jk) 
     838               zav           = zsqen * avmu(ji,jj,jk) 
     839               avt(ji,jj,jk) = MAX( zav, avtb(jk) )*tmask(ji,jj,jk) ! apply mask for zdfmxl routine 
     840               zav           = zsqen * avmv(ji,jj,jk) 
     841               avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 
     842            END DO 
     843         END DO 
     844      END DO 
    874845      ! 
    875846      ! Lateral boundary conditions (sign unchanged) 
    876       ! 
    877       avt(:,:,1)  = 0. 
     847      avt(:,:,1)  = 0._wp 
    878848      CALL lbc_lnk( avm, 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. ) 
    879849 
     
    881851         DO jj = 2, jpjm1 
    882852            DO ji = fs_2, fs_jpim1   ! vector opt. 
    883                avmu(ji,jj,jk) = ( avm  (ji,jj,jk)*tmask(ji,jj,jk) + avm (ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) * umask(ji,jj,jk)   & 
    884                &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) ) 
    885                avmv(ji,jj,jk) = ( avm  (ji,jj,jk)*tmask(ji,jj,jk) + avm (ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)   & 
    886                &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) ) 
    887             END DO 
    888          END DO 
    889       END DO 
    890  
    891       avmu(:,:,1) = 0. 
    892       avmv(:,:,1) = 0. 
    893  
    894       CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
     853               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
     854               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
     855            END DO 
     856         END DO 
     857      END DO 
     858      avmu(:,:,1) = 0._wp             ;   avmv(:,:,1) = 0._wp                 ! set surface to zero 
     859      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )       ! Lateral boundary conditions 
    895860 
    896861      IF(ln_ctl) THEN 
     
    918883      !! 
    919884      !!---------------------------------------------------------------------- 
     885      USE dynzdf_exp 
     886      USE trazdf_exp 
     887      ! 
    920888      INTEGER ::   jk    ! dummy loop indices 
    921889      REAL(wp)::   zcr   ! local scalar 
     
    929897      !!---------------------------------------------------------- 
    930898 
    931       ! Read Namelist namzdf_gls 
    932       ! ------------------------ 
    933       REWIND ( numnam ) 
     899      REWIND ( numnam )                !* Read Namelist namzdf_gls 
    934900      READ   ( numnam, namzdf_gls ) 
    935901 
    936       ! Parameter control and print 
    937       ! --------------------------- 
    938       ! Control print 
    939       IF(lwp) THEN 
     902      IF(lwp) THEN                     !* Control print 
    940903         WRITE(numout,*) 
    941904         WRITE(numout,*) 'zdf_gls_init : gls turbulent closure scheme' 
    942905         WRITE(numout,*) '~~~~~~~~~~~~' 
    943          WRITE(numout,*) '          Namelist namzdf_gls : set gls mixing parameters' 
    944          WRITE(numout,*) '             minimum value of en                       rn_emin  = ', rn_emin 
    945          WRITE(numout,*) '             minimum value of eps                    rn_epsmin  = ', rn_epsmin 
    946          WRITE(numout,*) '             Surface roughness (m)                         hsro = ', hsro 
    947          WRITE(numout,*) '             Bottom roughness (m)                        rn_hbro = ', rn_hbro 
    948          WRITE(numout,*) '             Limit dissipation rate under stable stratification ln_length_lim = ',ln_length_lim 
    949          WRITE(numout,*) '             Galperin length scale limitation coef (Standard: 0.53, Holt: 0.26) rn_clim_galp = ', rn_clim_galp 
    950          WRITE(numout,*) '             TKE Surface boundary condition       nn_tkebc_surf = ', nn_tkebc_surf 
    951          WRITE(numout,*) '             TKE Bottom boundary condition        nn_tkebc_bot  = ', nn_tkebc_bot 
    952          WRITE(numout,*) '             PSI Surface boundary condition       nn_psibc_surf = ', nn_psibc_surf 
    953          WRITE(numout,*) '             PSI Bottom boundary condition        nn_psibc_bot  = ', nn_psibc_bot 
    954          WRITE(numout,*) '             Craig and Banner scheme                  ln_crban  = ', ln_crban 
    955          WRITE(numout,*) '             Modify psi Schmidt number (wb case)     ln_sigpsi  = ', ln_sigpsi 
    956          WRITE(numout,*) '             Craig and Banner coef.                   rn_crban  = ', rn_crban 
    957          WRITE(numout,*) '             Charnock coef.                           rn_charn  = ', rn_charn 
    958          WRITE(numout,*) '             Stability functions                   nn_stab_func = ',nn_stab_func 
    959          WRITE(numout,*) '             Closure                                 nn_clos    = ',nn_clos 
     906         WRITE(numout,*) '   Namelist namzdf_gls : set gls mixing parameters' 
     907         WRITE(numout,*) '      minimum value of en                           rn_emin       = ', rn_emin 
     908         WRITE(numout,*) '      minimum value of eps                          rn_epsmin     = ', rn_epsmin 
     909         WRITE(numout,*) '      Surface roughness (m)                         hsro          = ', hsro 
     910         WRITE(numout,*) '      Bottom roughness (m)                          rn_hbro      = ', rn_hbro 
     911         WRITE(numout,*) '      Limit dissipation rate under stable stratif.  ln_length_lim = ', ln_length_lim 
     912         WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp = ', rn_clim_galp 
     913         WRITE(numout,*) '      TKE Surface boundary condition                nn_tkebc_surf = ', nn_tkebc_surf 
     914         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_tkebc_bot  = ', nn_tkebc_bot 
     915         WRITE(numout,*) '      PSI Surface boundary condition                nn_psibc_surf = ', nn_psibc_surf 
     916         WRITE(numout,*) '      PSI Bottom boundary condition                 nn_psibc_bot  = ', nn_psibc_bot 
     917         WRITE(numout,*) '      Craig and Banner scheme                       ln_crban      = ', ln_crban 
     918         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi     = ', ln_sigpsi 
     919         WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban 
     920         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn 
     921         WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func 
     922         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos 
    960923         WRITE(numout,*) 
    961924      ENDIF 
    962925 
    963       ! Check of some namelist values 
     926      !                                !* Check of some namelist values 
    964927      IF( nn_tkebc_surf < 0 .OR. nn_tkebc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_surf is 0 or 1' ) 
    965928      IF( nn_psibc_surf < 0 .OR. nn_psibc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_surf is 0 or 1' ) 
    966       IF( nn_tkebc_bot < 0 .OR. nn_tkebc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' ) 
    967       IF( nn_psibc_bot < 0 .OR. nn_psibc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' ) 
    968       IF( nn_stab_func < 0 .OR. nn_stab_func > 3) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) 
    969       IF( nn_clos < 0 .OR. nn_clos > 3) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) 
     929      IF( nn_tkebc_bot  < 0 .OR. nn_tkebc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' ) 
     930      IF( nn_psibc_bot  < 0 .OR. nn_psibc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' ) 
     931      IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) 
     932      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) 
    970933 
    971934      ! Initialisation of the parameters for the choosen closure 
     
    975938      ! 
    976939      CASE( 0 )               ! k-kl  (Mellor-Yamada) 
    977       ! 
    978       IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada' 
    979       rpp       = 0. 
    980       rmm       = 1. 
    981       rnn       = 1. 
    982       rsc_tke = 1.96 
    983       rsc_psi = 1.96 
    984       rpsi1 = 0.9 
    985       rpsi3p = 1. 
    986       rpsi2  = 0.5 
    987       ! 
     940         ! 
     941         IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada' 
     942         rpp     = 0._wp 
     943         rmm     = 1._wp 
     944         rnn     = 1._wp 
     945         rsc_tke = 1.96_wp 
     946         rsc_psi = 1.96_wp 
     947         rpsi1   = 0.9_wp 
     948         rpsi3p  = 1._wp 
     949         rpsi2   = 0.5_wp 
     950         ! 
    988951         SELECT CASE ( nn_stab_func ) 
    989          ! 
    990          CASE( 0, 1 )               ! G88 or KC stability functions 
    991             rpsi3m = 2.53 
    992          CASE( 2 )                  ! Canuto A stability functions 
    993             rpsi3m = 2.38 
    994          CASE( 3 )                  ! Canuto B stability functions 
    995             rpsi3m = 2.38         ! caution : constant not identified  
     952         CASE( 0, 1 )   ;   rpsi3m = 2.53_wp       ! G88 or KC stability functions 
     953         CASE( 2 )      ;   rpsi3m = 2.38_wp       ! Canuto A stability functions 
     954         CASE( 3 )      ;   rpsi3m = 2.38          ! Canuto B stability functions (caution : constant not identified) 
    996955         END SELECT 
    997       ! 
     956         ! 
    998957      CASE( 1 )               ! k-eps 
    999       ! 
    1000       IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps' 
    1001       rpp       = 3. 
    1002       rmm       = 1.5 
    1003       rnn       = -1. 
    1004       rsc_tke = 1. 
    1005       rsc_psi = 1.3  ! Schmidt number for psi 
    1006       rpsi1 = 1.44 
    1007       rpsi3p = 1. 
    1008       rpsi2  = 1.92 
    1009       ! 
    1010         SELECT CASE ( nn_stab_func ) 
    1011         !  
    1012         CASE( 0, 1 )               ! G88 or KC stability functions 
    1013            rpsi3m = -0.52 
    1014         CASE( 2 )                  ! Canuto A stability functions 
    1015            rpsi3m = -0.629  
    1016         CASE( 3 )                  ! Canuto B stability functions 
    1017            rpsi3m = -0.566 
    1018         END SELECT 
    1019       ! 
     958         ! 
     959         IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps' 
     960         rpp     =  3._wp 
     961         rmm     =  1.5_wp 
     962         rnn     = -1._wp 
     963         rsc_tke =  1._wp 
     964         rsc_psi =  1.3_wp  ! Schmidt number for psi 
     965         rpsi1   =  1.44_wp 
     966         rpsi3p  =  1._wp 
     967         rpsi2   =  1.92_wp 
     968         ! 
     969         SELECT CASE ( nn_stab_func ) 
     970         CASE( 0, 1 )   ;   rpsi3m = -0.52_wp      ! G88 or KC stability functions 
     971         CASE( 2 )      ;   rpsi3m = -0.629_wp     ! Canuto A stability functions 
     972         CASE( 3 )      ;   rpsi3m = -0.566        ! Canuto B stability functions 
     973         END SELECT 
     974         ! 
    1020975      CASE( 2 )               ! k-omega 
    1021       ! 
    1022       IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega' 
    1023       rpp       = -1. 
    1024       rmm       = 0.5 
    1025       rnn       = -1. 
    1026       rsc_tke = 2. 
    1027       rsc_psi = 2. 
    1028       rpsi1 = 0.555 
    1029       rpsi3p = 1. 
    1030       rpsi2  = 0.833 
    1031       ! 
    1032         SELECT CASE ( nn_stab_func ) 
    1033         ! 
    1034         CASE( 0, 1 )               ! G88 or KC stability functions 
    1035            rpsi3m = -0.58 
    1036         CASE( 2 )                  ! Canuto A stability functions 
    1037            rpsi3m = -0.64 
    1038         CASE( 3 )                  ! Canuto B stability functions 
    1039            rpsi3m = -0.64        ! caution : constant not identified 
    1040         END SELECT 
    1041       ! 
    1042       CASE( 3 )               ! gen 
    1043       ! 
    1044       IF(lwp) WRITE(numout,*) 'The choosen closure is generic' 
    1045       rpp       = 2. 
    1046       rmm       = 1. 
    1047       rnn       = -0.67 
    1048       rsc_tke = 0.8 
    1049       rsc_psi = 1.07 
    1050       rpsi1 = 1. 
    1051       rpsi3p = 1. 
    1052       rpsi2  = 1.22 
    1053       ! 
    1054         SELECT CASE ( nn_stab_func ) 
    1055         ! 
    1056         CASE( 0, 1 )               ! G88 or KC stability functions 
    1057            rpsi3m = 0.1 
    1058         CASE( 2 )                  ! Canuto A stability functions 
    1059            rpsi3m = 0.05 
    1060         CASE( 3 )                  ! Canuto B stability functions 
    1061            rpsi3m = 0.05         ! caution : constant not identified 
    1062         END SELECT 
    1063       ! 
     976         ! 
     977         IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega' 
     978         rpp     = -1._wp 
     979         rmm     =  0.5_wp 
     980         rnn     = -1._wp 
     981         rsc_tke =  2._wp 
     982         rsc_psi =  2._wp 
     983         rpsi1   =  0.555_wp 
     984         rpsi3p  =  1._wp 
     985         rpsi2   =  0.833_wp 
     986         ! 
     987         SELECT CASE ( nn_stab_func ) 
     988         CASE( 0, 1 )   ;   rpsi3m = -0.58_wp       ! G88 or KC stability functions 
     989         CASE( 2 )      ;   rpsi3m = -0.64_wp       ! Canuto A stability functions 
     990         CASE( 3 )      ;   rpsi3m = -0.64_wp       ! Canuto B stability functions caution : constant not identified) 
     991         END SELECT 
     992         ! 
     993      CASE( 3 )               ! generic 
     994         ! 
     995         IF(lwp) WRITE(numout,*) 'The choosen closure is generic' 
     996         rpp     = 2._wp 
     997         rmm     = 1._wp 
     998         rnn     = -0.67_wp 
     999         rsc_tke = 0.8_wp 
     1000         rsc_psi = 1.07_wp 
     1001         rpsi1   = 1._wp 
     1002         rpsi3p  = 1._wp 
     1003         rpsi2   = 1.22_wp 
     1004         ! 
     1005         SELECT CASE ( nn_stab_func ) 
     1006         CASE( 0, 1 )   ;   rpsi3m = 0.1_wp         ! G88 or KC stability functions 
     1007         CASE( 2 )      ;   rpsi3m = 0.05_wp        ! Canuto A stability functions 
     1008         CASE( 3 )      ;   rpsi3m = 0.05_wp        ! Canuto B stability functions caution : constant not identified) 
     1009         END SELECT 
     1010         ! 
    10641011      END SELECT 
    10651012 
     
    10701017      ! 
    10711018      CASE ( 0 )             ! Galperin stability functions 
    1072       ! 
    1073       IF(lwp) WRITE(numout,*) 'Stability functions from Galperin' 
    1074       rc2 = 0. 
    1075       rc3 = 0. 
    1076       rc_diff = 1. 
    1077       rc0     = 0.5544 
    1078       rcm_sf = 0.9884 
    1079       rghmin = -0.28 
    1080       rgh0   = 0.0233 
    1081       rghcri = 0.02 
    1082       ! 
     1019         ! 
     1020         IF(lwp) WRITE(numout,*) 'Stability functions from Galperin' 
     1021         rc2     =  0._wp 
     1022         rc3     =  0._wp 
     1023         rc_diff =  1._wp 
     1024         rc0     =  0.5544_wp 
     1025         rcm_sf  =  0.9884_wp 
     1026         rghmin  = -0.28_wp 
     1027         rgh0    =  0.0233_wp 
     1028         rghcri  =  0.02_wp 
     1029         ! 
    10831030      CASE ( 1 )             ! Kantha-Clayson stability functions 
    1084       ! 
    1085       IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson' 
    1086       rc2 = 0.7 
    1087       rc3 = 0.2 
    1088       rc_diff = 1. 
    1089       rc0     = 0.5544 
    1090       rcm_sf = 0.9884 
    1091       rghmin = -0.28 
    1092       rgh0   = 0.0233 
    1093       rghcri = 0.02 
    1094       ! 
     1031         ! 
     1032         IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson' 
     1033         rc2     =  0.7_wp 
     1034         rc3     =  0.2_wp 
     1035         rc_diff =  1._wp 
     1036         rc0     =  0.5544_wp 
     1037         rcm_sf  =  0.9884_wp 
     1038         rghmin  = -0.28_wp 
     1039         rgh0    =  0.0233_wp 
     1040         rghcri  =  0.02_wp 
     1041         ! 
    10951042      CASE ( 2 )             ! Canuto A stability functions 
    1096       ! 
    1097       IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A' 
    1098       rs0 = 1.5*rl1*rl5**2. 
    1099       rs1 = -rl4*(rl6+rl7) + 2.*rl4*rl5*(rl1-(1./3.)*rl2-rl3)+1.5*rl1*rl5*rl8 
    1100       rs2 = -(3./8.)*rl1*(rl6**2.-rl7**2.) 
    1101       rs4 = 2.*rl5 
    1102       rs5 = 2.*rl4 
    1103       rs6 = (2./3.)*rl5*(3.*rl3**2.-rl2**2.)-0.5*rl5*rl1*(3.*rl3-rl2)+0.75*rl1*(rl6-rl7) 
    1104       rd0 = 3*rl5**2. 
    1105       rd1 = rl5*(7.*rl4+3.*rl8) 
    1106       rd2 = rl5**2.*(3.*rl3**2.-rl2**2.)-0.75*(rl6**2.-rl7**2.) 
    1107       rd3 = rl4*(4.*rl4+3.*rl8) 
    1108       rd4 = rl4*(rl2*rl6-3.*rl3*rl7-rl5*(rl2**2.-rl3**2.))+rl5*rl8*(3.*rl3**2.-rl2**2.) 
    1109       rd5 = 0.25*(rl2**2.-3.*rl3**2.)*(rl6**2.-rl7**2.) 
    1110       rc0 = 0.5268 
    1111       rf6 = 8. / (rc0**6.) 
    1112       rc_diff = SQRT(2.)/(rc0**3.) 
    1113       rcm_sf = 0.7310 
    1114       rghmin = -0.28 
    1115       rgh0   = 0.0329 
    1116       rghcri = 0.03 
    1117       ! 
     1043         ! 
     1044         IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A' 
     1045         rs0 = 1.5_wp * rl1 * rl5*rl5 
     1046         rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8 
     1047         rs2 = -(3._wp/8._wp) * rl1*(rl6*rl6-rl7*rl7) 
     1048         rs4 = 2._wp * rl5 
     1049         rs5 = 2._wp * rl4 
     1050         rs6 = (2._wp/3._wp) * rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.5_wp * rl5*rl1 * (3._wp*rl3-rl2)   & 
     1051            &                                                    + 0.75_wp * rl1 * ( rl6 - rl7 ) 
     1052         rd0 = 3._wp * rl5*rl5 
     1053         rd1 = rl5 * ( 7._wp*rl4 + 3._wp*rl8 ) 
     1054         rd2 = rl5*rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.75_wp*(rl6*rl6 - rl7*rl7 ) 
     1055         rd3 = rl4 * ( 4._wp*rl4 + 3._wp*rl8) 
     1056         rd4 = rl4 * ( rl2 * rl6 - 3._wp*rl3*rl7 - rl5*(rl2*rl2 - rl3*rl3 ) ) + rl5*rl8 * ( 3._wp*rl3*rl3 - rl2*rl2 ) 
     1057         rd5 = 0.25_wp * ( rl2*rl2 - 3._wp *rl3*rl3 ) * ( rl6*rl6 - rl7*rl7 ) 
     1058         rc0 = 0.5268_wp 
     1059         rf6 = 8._wp / (rc0**6._wp) 
     1060         rc_diff = SQRT(2._wp) / (rc0**3._wp) 
     1061         rcm_sf  =  0.7310_wp 
     1062         rghmin  = -0.28_wp 
     1063         rgh0    =  0.0329_wp 
     1064         rghcri  =  0.03_wp 
     1065         ! 
    11181066      CASE ( 3 )             ! Canuto B stability functions 
    1119       ! 
    1120       IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B' 
    1121       rs0 = 1.5*rm1*rm5**2. 
    1122       rs1 = -rm4*(rm6+rm7) + 2.*rm4*rm5*(rm1-(1./3.)*rm2-rm3)+1.5*rm1*rm5*rm8 
    1123       rs2 = -(3./8.)*rm1*(rm6**2.-rm7**2.) 
    1124       rs4 = 2.*rm5 
    1125       rs5 = 2.*rm4 
    1126       rs6 = (2./3.)*rm5*(3.*rm3**2.-rm2**2.)-0.5*rm5*rm1*(3.*rm3-rm2)+0.75*rm1*(rm6-rm7) 
    1127       rd0 = 3*rm5**2. 
    1128       rd1 = rm5*(7.*rm4+3.*rm8) 
    1129       rd2 = rm5**2.*(3.*rm3**2.-rm2**2.)-0.75*(rm6**2.-rm7**2.) 
    1130       rd3 = rm4*(4.*rm4+3.*rm8) 
    1131       rd4 = rm4*(rm2*rm6-3.*rm3*rm7-rm5*(rm2**2.-rm3**2.))+rm5*rm8*(3.*rm3**2.-rm2**2.) 
    1132       rd5 = 0.25*(rm2**2.-3.*rm3**2.)*(rm6**2.-rm7**2.) 
    1133       rc0 = 0.5268  !!       rc0 = 0.5540 (Warner ...) to verify ! 
    1134       rf6 = 8. / (rc0**6.) 
    1135       rc_diff = SQRT(2.)/(rc0**3.) 
    1136       rcm_sf = 0.7470 
    1137       rghmin = -0.28 
    1138       rgh0   = 0.0444 
    1139       rghcri = 0.0414 
    1140       ! 
     1067         ! 
     1068         IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B' 
     1069         rs0 = 1.5_wp * rm1 * rm5*rm5 
     1070         rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8 
     1071         rs2 = -(3._wp/8._wp) * rm1 * (rm6*rm6-rm7*rm7 ) 
     1072         rs4 = 2._wp * rm5 
     1073         rs5 = 2._wp * rm4 
     1074         rs6 = (2._wp/3._wp) * rm5 * (3._wp*rm3*rm3-rm2*rm2) - 0.5_wp * rm5*rm1*(3._wp*rm3-rm2) + 0.75_wp * rm1*(rm6-rm7) 
     1075         rd0 = 3._wp * rm5*rm5 
     1076         rd1 = rm5 * (7._wp*rm4 + 3._wp*rm8) 
     1077         rd2 = rm5*rm5 * (3._wp*rm3*rm3 - rm2*rm2) - 0.75_wp * (rm6*rm6 - rm7*rm7) 
     1078         rd3 = rm4 * ( 4._wp*rm4 + 3._wp*rm8 ) 
     1079         rd4 = rm4 * ( rm2*rm6 -3._wp*rm3*rm7 - rm5*(rm2*rm2 - rm3*rm3) ) + rm5 * rm8 * ( 3._wp*rm3*rm3 - rm2*rm2 ) 
     1080         rd5 = 0.25_wp * ( rm2*rm2 - 3._wp*rm3*rm3 ) * ( rm6*rm6 - rm7*rm7 ) 
     1081         rc0 = 0.5268_wp            !!       rc0 = 0.5540_wp (Warner ...) to verify ! 
     1082         rf6 = 8._wp / ( rc0**6._wp ) 
     1083         rc_diff = SQRT(2._wp)/(rc0**3.) 
     1084         rcm_sf  =  0.7470_wp 
     1085         rghmin  = -0.28_wp 
     1086         rgh0    =  0.0444_wp 
     1087         rghcri  =  0.0414_wp 
     1088         ! 
    11411089      END SELECT 
    11421090     
    1143       ! Set Schmidt number for psi diffusion 
    1144       ! In the wave breaking case 
     1091      ! Set Schmidt number for psi diffusion in the wave breaking case 
    11451092      ! See equation 13 of Carniel et al, Ocean modelling, 30, 225-239, 2009 
    11461093      ! or equation (17) of Burchard, JPO, 31, 3133-3145, 2001 
    1147       IF ((ln_sigpsi).AND.(ln_crban)) THEN 
    1148          zcr = SQRT(1.5*rsc_tke) * rcm_sf /vkarmn 
    1149          rsc_psi0 = vkarmn**2/(rpsi2*rcm_sf**2) *              &  
    1150         &             (  rnn**2 - 4./3.*zcr*rnn*rmm - 1./3.*zcr*rnn  & 
    1151         &              + 2./9.*rmm*zcr**2 + 4./9.*zcr**2*rmm**2)                                  
     1094      IF( ln_sigpsi .AND. ln_crban ) THEN 
     1095         zcr = SQRT( 1.5_wp*rsc_tke ) * rcm_sf / vkarmn 
     1096         rsc_psi0 = vkarmn*vkarmn / ( rpsi2 * rcm_sf*rcm_sf )                       &  
     1097        &         * ( rnn*rnn - 4._wp/3._wp * zcr*rnn*rmm - 1._wp/3._wp * zcr*rnn   & 
     1098        &           + 2._wp/9._wp * rmm * zcr*zcr + 4._wp/9._wp * zcr*zcr * rmm*rmm )                                  
    11521099      ELSE 
    11531100         rsc_psi0 = rsc_psi 
     
    11561103      ! Shear free turbulence parameters: 
    11571104      ! 
    1158       ra_sf  = -4.*rnn*SQRT(rsc_tke) / ( (1.+4.*rmm)*SQRT(rsc_tke) & 
    1159                &                              - SQRT(rsc_tke + 24.*rsc_psi0*rpsi2 ) ) 
    1160       rl_sf  = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1. + 4.*rmm + 8.*rmm**2)*rsc_tke       & 
    1161                &                                          + 12. * rsc_psi0*rpsi2 - (1. + 4.*rmm) & 
    1162                &                                          *SQRT(rsc_tke*(rsc_tke                & 
    1163                &                                                + 24.*rsc_psi0*rpsi2)) )         & 
    1164                &                                          /(12.*rnn**2.) & 
    1165                &                                       ) 
    1166  
    1167       ! Control print 
    1168       ! 
    1169       IF(lwp) THEN 
     1105      ra_sf  = -4._wp * rnn * SQRT( rsc_tke ) / ( (1._wp+4._wp*rmm) * SQRT( rsc_tke )   & 
     1106         &                                      - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 
     1107      rl_sf  = rc0 * SQRT( rc0 / rcm_sf )                                                                   & 
     1108         &         * SQRT(  (  (1._wp + 4._wp*rmm + 8._wp*rmm*rmm) * rsc_tke                                & 
     1109         &                   + 12._wp * rsc_psi0 * rpsi2                                                    & 
     1110         &                   - (1._wp + 4._wp*rmm) * SQRT( rsc_tke*(rsc_tke+ 24._wp*rsc_psi0*rpsi2) )  )    & 
     1111         &                / ( 12._wp*rnn*rnn )                                                              ) 
     1112 
     1113      ! 
     1114      IF(lwp) THEN      ! Control print 
    11701115         WRITE(numout,*) 
    11711116         WRITE(numout,*) 'Limit values' 
     
    11911136 
    11921137      ! Constants initialization 
    1193       rc02r = 1. / rc0**2. 
    1194       rc02  = rc0**2._wp 
    1195       rc03  = rc0**3._wp 
    1196       rc04  = rc0**4._wp 
     1138      rc02  = rc0  * rc0   ;   rc02r = 1. / rc02 
     1139      rc03  = rc02 * rc0 
     1140      rc04  = rc03 * rc0 
    11971141      rc03_sqrt2_galp = rc03 / SQRT(2._wp) / rn_clim_galp 
    1198       rsbc_mb   = 0.5 * (15.8*rn_crban)**(2./3.) ! Surf. bound. cond. from Mellor and Blumberg 
    1199       rsbc_std  = 3.75                           ! Surf. bound. cond. standard (prod=diss) 
    1200       rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2./3.) ! k_eps = 53.  Dirichlet + Wave breaking  
    1201       rsbc_tke2 = 0.5 / rau0 
     1142      rsbc_mb   = 0.5_wp * (15.8_wp*rn_crban)**(2._wp/3._wp)              ! Surf. bound. cond. from Mellor and Blumberg 
     1143      rsbc_std  = 3.75_wp                                                  ! Surf. bound. cond. standard (prod=diss) 
     1144      rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp) ! k_eps = 53.  Dirichlet + Wave breaking  
     1145      rsbc_tke2 = 0.5_wp / rau0 
    12021146      rsbc_tke3 = rdt * rn_crban                                                         ! Neumann + Wave breaking 
    1203       rsbc_zs   = rn_charn/grav                                                          ! Charnock formula 
     1147      rsbc_zs   = rn_charn / grav                                                        ! Charnock formula 
    12041148      rsbc_psi1 = rc0**rpp * rsbc_tke1**rmm * rl_sf**rnn                           ! Dirichlet + Wave breaking 
    1205       rsbc_psi2 = -0.5 * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi                   ! Neumann + NO Wave breaking  
    1206       rsbc_psi3 = -0.5 * rdt * rc0**rpp * rl_sf**rnn / rsc_psi  * (rnn + rmm*ra_sf) ! Neumann + Wave breaking 
    1207       rfact_tke = -0.5 / rsc_tke * rdt               ! Cst used for the Diffusion term of tke 
    1208       rfact_psi = -0.5 / rsc_psi * rdt               ! Cst used for the Diffusion term of tke 
    1209  
    1210       ! Wall proximity function 
     1149      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi                   ! Neumann + NO Wave breaking  
     1150      rsbc_psi3 = -0.5_wp * rdt * rc0**rpp * rl_sf**rnn / rsc_psi  * (rnn + rmm*ra_sf) ! Neumann + Wave breaking 
     1151      rfact_tke = -0.5_wp / rsc_tke * rdt               ! Cst used for the Diffusion term of tke 
     1152      rfact_psi = -0.5_wp / rsc_psi * rdt               ! Cst used for the Diffusion term of tke 
     1153 
     1154      !                                !* Wall proximity function 
    12111155      zwall (:,:,:) = 1._wp * tmask(:,:,:) 
    12121156 
    1213       !                               !* set vertical eddy coef. to the background value 
     1157      !                                !* set vertical eddy coef. to the background value 
    12141158      DO jk = 1, jpk 
    12151159         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 
     
    12181162         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 
    12191163      END DO 
    1220       !                               !* read or initialize all required files  
     1164      !                                !* read or initialize all required files  
    12211165      CALL gls_rst( nit000, 'READ' ) 
    12221166      ! 
     
    12361180     INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
    12371181     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    1238      !! 
     1182     ! 
    12391183     INTEGER ::   jit, jk   ! dummy loop indices 
    12401184     INTEGER ::   id1, id2, id3, id4, id5, id6, id7, id8 
     
    13271271      WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt 
    13281272   END SUBROUTINE zdf_gls 
    1329    SUBROUTINE zdf_gls_init          ! Empty routine 
    1330    END SUBROUTINE zdf_gls_init 
    1331    SUBROUTINE gls_rst( kt,cdrw )          ! Empty routine 
     1273   SUBROUTINE gls_rst( kt, cdrw )          ! Empty routine 
    13321274      INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
    13331275      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    1334       WRITE(*,*) 'gls_Rst: You should not have seen this print! error?', kt, cdrw 
     1276      WRITE(*,*) 'gls_rst: You should not have seen this print! error?', kt, cdrw 
    13351277   END SUBROUTINE gls_rst 
    13361278#endif 
     
    13381280   !!====================================================================== 
    13391281END MODULE zdfgls 
     1282 
Note: See TracChangeset for help on using the changeset viewer.