New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 14986 for NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/ZDF/zdfgls.F90 – NEMO

Ignore:
Timestamp:
2021-06-14T13:34:08+02:00 (3 years ago)
Author:
sparonuz
Message:

Merge trunk -r14984:HEAD

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/ZDF/zdfgls.F90

    r14200 r14986  
    2626   USE zdfmxl         ! mixed layer 
    2727   USE sbcwave , ONLY : hsw   ! significant wave height 
     28#if defined key_si3 
     29   USE ice, ONLY: hm_i, h_i 
     30#endif 
     31#if defined key_cice 
     32   USE sbc_ice, ONLY: h_i 
     33#endif 
    2834   ! 
    2935   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    … …  
    5157   LOGICAL  ::   ln_length_lim     ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988) 
    5258   LOGICAL  ::   ln_sigpsi         ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing 
     59   INTEGER  ::   nn_mxlice         ! type of scaling under sea-ice (=0/1/2/3) 
    5360   INTEGER  ::   nn_bc_surf        ! surface boundary condition (=0/1) 
    5461   INTEGER  ::   nn_bc_bot         ! bottom boundary condition (=0/1) 
    … …  
    137144      USE zdf_oce , ONLY : en, avtb, avmb   ! ocean vertical physics 
    138145      !! 
    139       INTEGER                   , INTENT(in   ) ::   kt             ! ocean time step 
    140       INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices 
    141       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term 
    142       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points) 
     146      INTEGER                             , INTENT(in   ) ::   kt             ! ocean time step 
     147      INTEGER                             , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices 
     148      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   p_sh2          ! shear production term 
     149      REAL(wp), DIMENSION(:,:,:)          , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points) 
    143150      ! 
    144151      INTEGER  ::   ji, jj, jk    ! dummy loop arguments 
    … …  
    151158      REAL(wp) ::   gh, gm, shr, dif, zsqen, zavt, zavm !   -      - 
    152159      REAL(wp) ::   zmsku, zmskv                        !   -      - 
    153       REAL(wp), DIMENSION(jpi,jpj)     ::   zdep 
    154       REAL(wp), DIMENSION(jpi,jpj)     ::   zkar 
    155       REAL(wp), DIMENSION(jpi,jpj)     ::   zflxs       ! Turbulence fluxed induced by internal waves 
    156       REAL(wp), DIMENSION(jpi,jpj)     ::   zhsro       ! Surface roughness (surface waves) 
    157       REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra    ! Tapering of wave breaking under sea ice 
    158       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eb          ! tke at time before 
    159       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   hmxl_b      ! mixing length at time before 
    160       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eps         ! dissipation rate 
    161       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi) 
    162       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   psi         ! psi at time now 
    163       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zd_lw, zd_up, zdiag   ! lower, upper  and diagonal of the matrix 
    164       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zstt, zstm  ! stability function on tracer and momentum 
     160      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zdep 
     161      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zkar 
     162      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zflxs                 ! Turbulence fluxed induced by internal waves 
     163      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zhsro                 ! Surface roughness (surface waves) 
     164      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zice_fra              ! Tapering of wave breaking under sea ice 
     165      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   eb                    ! tke at time before 
     166      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   hmxl_b                ! mixing length at time before 
     167      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   eps                   ! dissipation rate 
     168      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zwall_psi             ! Wall function use in the wb case (ln_sigpsi) 
     169      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   psi                   ! psi at time now 
     170      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zd_lw, zd_up, zdiag   ! lower, upper  and diagonal of the matrix 
     171      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zstt, zstm            ! stability function on tracer and momentum 
    165172      !!-------------------------------------------------------------------- 
    166173      ! 
    167174      ! Preliminary computing 
    168  
    169       ustar2_surf(:,:) = 0._wp   ;         psi(:,:,:) = 0._wp 
    170       ustar2_top (:,:) = 0._wp   ;   zwall_psi(:,:,:) = 0._wp 
    171       ustar2_bot (:,:) = 0._wp 
     175      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     176         ustar2_surf(ji,jj) = 0._wp   ;   ustar2_top(ji,jj) = 0._wp   ;   ustar2_bot(ji,jj) = 0._wp 
     177      END_2D 
     178      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 
     179         psi(ji,jj,jk) = 0._wp   ;   zwall_psi(ji,jj,jk) = 0._wp 
     180      END_3D 
    172181 
    173182      SELECT CASE ( nn_z0_ice ) 
    174183      CASE( 0 )   ;   zice_fra(:,:) = 0._wp 
    175       CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp ) 
    176       CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:) 
    177       CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp ) 
     184      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(A2D(nn_hls)) * 10._wp ) 
     185      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(A2D(nn_hls)) 
     186      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp ) 
    178187      END SELECT 
    179188 
    180189      ! Compute surface, top and bottom friction at T-points 
    181       DO_2D( 0, 0, 0, 0 )          !==  surface ocean friction 
     190      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )          !==  surface ocean friction 
    182191         ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1)   ! surface friction 
    183192      END_2D 
    … …  
    186195      ! 
    187196      IF( .NOT.ln_drg_OFF ) THEN     !== top/bottom friction   (explicit before friction) 
    188          DO_2D( 0, 0, 0, 0 )         ! bottom friction (explicit before friction) 
     197         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )          ! bottom friction (explicit before friction) 
    189198            zmsku = 0.5_wp * ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 
    190199            zmskv = 0.5_wp * ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0) 
    … …  
    193202         END_2D 
    194203         IF( ln_isfcav ) THEN 
    195             DO_2D( 0, 0, 0, 0 )      ! top friction 
     204            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )      ! top friction 
    196205               zmsku = 0.5_wp * ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 
    197206               zmskv = 0.5_wp * ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0) 
    … …  
    206215         zhsro(:,:) = rn_hsro 
    207216      CASE ( 1 )             ! Standard Charnock formula 
    208          zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(:,:) , rn_hsro ) 
     217         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     218            zhsro(ji,jj) = MAX( rsbc_zs1 * ustar2_surf(ji,jj) , rn_hsro ) 
     219         END_2D 
    209220      CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 
    210221!!gm faster coding : the 2 comment lines should be used 
    211222!!gm         zcof = 2._wp * 0.6_wp / 28._wp 
    212223!!gm         zdep(:,:)  = 30._wp * TANH(  zcof/ SQRT( MAX(ustar2_surf(:,:),rsmall) )  )       ! Wave age (eq. 10) 
    213          zdep (:,:) = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(:,:),rsmall))) )         ! Wave age (eq. 10) 
    214          zhsro(:,:) = MAX(rsbc_zs2 * ustar2_surf(:,:) * zdep(:,:)**1.5, rn_hsro)          ! zhsro = rn_frac_hs * Hsw (eq. 11) 
     224         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     225            zcof = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(ji,jj),rsmall))) )          ! Wave age (eq. 10) 
     226            zhsro(ji,jj) = MAX(rsbc_zs2 * ustar2_surf(ji,jj) * zcof**1.5, rn_hsro)        ! zhsro = rn_frac_hs * Hsw (eq. 11) 
     227         END_2D 
    215228      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file) 
    216          zhsro(:,:) = MAX(rn_frac_hs * hsw(:,:), rn_hsro)   ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 ) 
     229         zhsro(:,:) = MAX(rn_frac_hs * hsw(A2D(nn_hls)), rn_hsro)   ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 ) 
    217230      END SELECT 
    218231      ! 
    219232      ! adapt roughness where there is sea ice 
    220       zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1)  + (1._wp - tmask(:,:,1))*rn_hsro 
    221       ! 
    222       DO_3D( 0, 0, 0, 0, 2, jpkm1 )  !==  Compute dissipation rate  ==! 
     233      SELECT CASE( nn_mxlice )       ! Type of scaling under sea-ice 
     234      ! 
     235      CASE( 1 )                      ! scaling with constant sea-ice roughness 
     236         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     237            zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * rn_hsri )*tmask(ji,jj,1)  + (1._wp - tmask(ji,jj,1))*rn_hsro 
     238         END_2D 
     239         ! 
     240      CASE( 2 )                      ! scaling with mean sea-ice thickness 
     241#if defined key_si3 
     242         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     243            zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * hm_i(ji,jj) )*tmask(ji,jj,1)  + (1._wp - tmask(ji,jj,1))*rn_hsro 
     244         END_2D 
     245#endif 
     246         ! 
     247      CASE( 3 )                      ! scaling with max sea-ice thickness 
     248#if defined key_si3 || defined key_cice 
     249         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     250            zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * MAXVAL(h_i(ji,jj,:)) )*tmask(ji,jj,1)  + (1._wp - tmask(ji,jj,1))*rn_hsro 
     251         END_2D 
     252#endif 
     253         ! 
     254      END SELECT 
     255      ! 
     256      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )  !==  Compute dissipation rate  ==! 
    223257         eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 
    224258      END_3D 
    225259 
    226260      ! Save tke at before time step 
    227       eb    (:,:,:) = en    (:,:,:) 
    228       hmxl_b(:,:,:) = hmxl_n(:,:,:) 
     261      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 
     262         eb    (ji,jj,jk) = en    (ji,jj,jk) 
     263         hmxl_b(ji,jj,jk) = hmxl_n(ji,jj,jk) 
     264      END_3D 
    229265 
    230266      IF( nn_clos == 0 ) THEN    ! Mellor-Yamada 
    231          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     267         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    232268            zup   = hmxl_n(ji,jj,jk) * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 
    233269            zdown = vkarmn * gdepw(ji,jj,jk,Kmm) * ( -gdepw(ji,jj,jk,Kmm) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ) 
    … …  
    250286      ! Warning : after this step, en : right hand side of the matrix 
    251287 
    252       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     288      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    253289         ! 
    254290         buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)     ! stratif. destruction 
    … …  
    303339      ! 
    304340      CASE ( 0 )             ! Dirichlet boundary condition (set e at k=1 & 2) 
    305       ! First level 
    306       en   (:,:,1) = MAX(  rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3  ) 
    307       zd_lw(:,:,1) = en(:,:,1) 
    308       zd_up(:,:,1) = 0._wp 
    309       zdiag(:,:,1) = 1._wp 
    310       ! 
    311       ! One level below 
    312       en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm)) & 
    313          &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp) , rn_emin   ) 
    314       zd_lw(:,:,2) = 0._wp 
    315       zd_up(:,:,2) = 0._wp 
    316       zdiag(:,:,2) = 1._wp 
    317       ! 
    318       ! 
     341         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     342            ! First level 
     343            en   (ji,jj,1) = MAX(  rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3  ) 
     344            zd_lw(ji,jj,1) = en(ji,jj,1) 
     345            zd_up(ji,jj,1) = 0._wp 
     346            zdiag(ji,jj,1) = 1._wp 
     347            ! 
     348            ! One level below 
     349            en   (ji,jj,2) =  MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1          & 
     350               &                             * ((zhsro(ji,jj)+gdepw(ji,jj,2,Kmm)) / zhsro(ji,jj) )**(1.5_wp*ra_sf)  )**r2_3 ) 
     351            zd_lw(ji,jj,2) = 0._wp 
     352            zd_up(ji,jj,2) = 0._wp 
     353            zdiag(ji,jj,2) = 1._wp 
     354         END_2D 
     355         ! 
     356         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity) 
     357            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     358               IF( mikt(ji,jj) > 1 )THEN 
     359                  itop   = mikt(ji,jj)       ! k   top w-point 
     360                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one 
     361                  !                                                ! mask at the 
     362                  !                                                ocean surface 
     363                  !                                                points 
     364                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 
     365                  ! 
     366                  ! Dirichlet condition applied at: 
     367                  !     top level (itop)         &      Just below it (itopp1) 
     368                  zd_lw(ji,jj,itop) = 0._wp   ;   zd_lw(ji,jj,itopp1) = 0._wp 
     369                  zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp 
     370                  zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = 1._wp 
     371                  en   (ji,jj,itop) = z_en    ;   en   (ji,jj,itopp1) = z_en 
     372               ENDIF 
     373            END_2D 
     374         ENDIF 
     375         ! 
    319376      CASE ( 1 )             ! Neumann boundary condition (set d(e)/dz) 
    320       ! 
    321       ! Dirichlet conditions at k=1 
    322       en   (:,:,1) = MAX(  rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 , rn_emin  ) 
    323       zd_lw(:,:,1) = en(:,:,1) 
    324       zd_up(:,:,1) = 0._wp 
    325       zdiag(:,:,1) = 1._wp 
    326       ! 
    327       ! at k=2, set de/dz=Fw 
    328       !cbr 
    329       DO_2D( 0, 0, 0, 0 )   ! zdiag zd_lw not defined/used on the halo 
    330          zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 
    331          zd_lw(ji,jj,2) = 0._wp 
    332       END_2D 
    333       zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) )) 
    334       zflxs(:,:)   = rsbc_tke2 * (1._wp-zice_fra(:,:)) * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 
    335           &                    * (  ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:)  )**(1.5_wp*ra_sf) 
     377         ! 
     378         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     379            ! Dirichlet conditions at k=1 
     380            en   (ji,jj,1) = MAX(  rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3  ) 
     381            zd_lw(ji,jj,1) = en(ji,jj,1) 
     382            zd_up(ji,jj,1) = 0._wp 
     383            zdiag(ji,jj,1) = 1._wp 
     384            ! 
     385            ! at k=2, set de/dz=Fw 
     386            !cbr 
     387            ! zdiag zd_lw not defined/used on the halo 
     388            zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 
     389            zd_lw(ji,jj,2) = 0._wp 
     390            ! 
     391            zkar (ji,jj)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj)) )) 
     392            zflxs(ji,jj)   = rsbc_tke2 * (1._wp-zice_fra(ji,jj)) * ustar2_surf(ji,jj)**1.5_wp * zkar(ji,jj) & 
     393                &                    * (  ( zhsro(ji,jj)+gdept(ji,jj,1,Kmm) ) / zhsro(ji,jj)  )**(1.5_wp*ra_sf) 
    336394!!gm why not   :                        * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf) 
    337       en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm) 
    338       ! 
    339       ! 
     395            en(ji,jj,2) = en(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm) 
     396         END_2D 
     397         ! 
     398         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity) 
     399            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     400               IF( mikt(ji,jj) > 1 )THEN 
     401                  itop   = mikt(ji,jj)       ! k   top w-point 
     402                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one 
     403                  !                                                ! mask at the 
     404                  !                                                ocean surface 
     405                  !                                                points 
     406                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 
     407                  ! 
     408                  ! Bottom level Dirichlet condition: 
     409                  !     Bottom level (ibot)      &      Just above it (ibotm1) 
     410                  !         Dirichlet            !         Neumann 
     411                  zd_lw(ji,jj,itop) = 0._wp   !   ! Remove zd_up from zdiag 
     412                  zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) 
     413                  zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp 
     414                  en   (ji,jj,itop) = z_en 
     415               ENDIF 
     416            END_2D 
     417         ENDIF 
     418         ! 
    340419      END SELECT 
    341420 
    … …  
    348427         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin 
    349428         !                      ! Balance between the production and the dissipation terms 
    350          DO_2D( 0, 0, 0, 0 ) 
     429         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    351430!!gm This means that bottom and ocean w-level above have a specified "en" value.   Sure ???? 
    352431!!   With thick deep ocean level thickness, this may be quite large, no ??? 
    … …  
    365444         END_2D 
    366445         ! 
     446         ! NOTE: ctl_stop with ln_isfcav when using GLS 
    367447         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity) 
    368             DO_2D( 0, 0, 0, 0 ) 
     448            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    369449               itop   = mikt(ji,jj)       ! k   top w-point 
    370450               itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one 
    … …  
    384464      CASE ( 1 )             ! Neumman boundary condition 
    385465         ! 
    386          DO_2D( 0, 0, 0, 0 ) 
     466         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    387467            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
    388468            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
    … …  
    398478            en   (ji,jj,ibot) = z_en 
    399479         END_2D 
     480         ! NOTE: ctl_stop with ln_isfcav when using GLS 
    400481         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity) 
    401             DO_2D( 0, 0, 0, 0 ) 
     482            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    402483               itop   = mikt(ji,jj)       ! k   top w-point 
    403484               itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one 
    … …  
    420501      ! ---------------------------------------------------------- 
    421502      ! 
    422       DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
     503      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    423504         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    424505      END_3D 
    425       DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
     506      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    426507         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) 
    427508      END_3D 
    428       DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )           ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
     509      DO_3DS_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    429510         en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    430511      END_3D 
    431512      !                                            ! set the minimum value of tke 
    432       en(:,:,:) = MAX( en(:,:,:), rn_emin ) 
     513      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 
     514         en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) 
     515      END_3D 
    433516 
    434517      !!----------------------------------------!! 
    … …  
    441524      ! 
    442525      CASE( 0 )               ! k-kl  (Mellor-Yamada) 
    443          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     526         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    444527            psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 
    445528         END_3D 
    446529         ! 
    447530      CASE( 1 )               ! k-eps 
    448          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     531         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    449532            psi(ji,jj,jk)  = eps(ji,jj,jk) 
    450533         END_3D 
    451534         ! 
    452535      CASE( 2 )               ! k-w 
    453          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     536         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    454537            psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 
    455538         END_3D 
    456539         ! 
    457540      CASE( 3 )               ! generic 
    458          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     541         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    459542            psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 
    460543         END_3D 
    … …  
    469552      ! Warning : after this step, en : right hand side of the matrix 
    470553 
    471       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     554      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    472555         ! 
    473556         ! psi / k 
    … …  
    516599      CASE ( 0 )             ! Dirichlet boundary conditions 
    517600         ! 
    518          ! Surface value 
    519          zdep    (:,:)   = zhsro(:,:) * rl_sf ! Cosmetic 
    520          psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    521          zd_lw(:,:,1) = psi(:,:,1) 
    522          zd_up(:,:,1) = 0._wp 
    523          zdiag(:,:,1) = 1._wp 
    524          ! 
    525          ! One level below 
    526          zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(:,:,2,Kmm)/zhsro(:,:) ))) 
    527          zdep    (:,:)   = (zhsro(:,:) + gdepw(:,:,2,Kmm)) * zkar(:,:) 
    528          psi     (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    529          zd_lw(:,:,2) = 0._wp 
    530          zd_up(:,:,2) = 0._wp 
    531          zdiag(:,:,2) = 1._wp 
     601         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     602            ! Surface value 
     603            zdep    (ji,jj)   = zhsro(ji,jj) * rl_sf ! Cosmetic 
     604            psi     (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1) 
     605            zd_lw(ji,jj,1) = psi(ji,jj,1) 
     606            zd_up(ji,jj,1) = 0._wp 
     607            zdiag(ji,jj,1) = 1._wp 
     608            ! 
     609            ! One level below 
     610            zkar    (ji,jj)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(ji,jj,2,Kmm)/zhsro(ji,jj) ))) 
     611            zdep    (ji,jj)   = (zhsro(ji,jj) + gdepw(ji,jj,2,Kmm)) * zkar(ji,jj) 
     612            psi     (ji,jj,2) = rc0**rpp * en(ji,jj,2)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1) 
     613            zd_lw(ji,jj,2) = 0._wp 
     614            zd_up(ji,jj,2) = 0._wp 
     615            zdiag(ji,jj,2) = 1._wp 
     616         END_2D 
    532617         ! 
    533618      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz 
    534619         ! 
    535          ! Surface value: Dirichlet 
    536          zdep    (:,:)   = zhsro(:,:) * rl_sf 
    537          psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    538          zd_lw(:,:,1) = psi(:,:,1) 
    539          zd_up(:,:,1) = 0._wp 
    540          zdiag(:,:,1) = 1._wp 
    541          ! 
    542          ! Neumann condition at k=2 
    543          DO_2D( 0, 0, 0, 0 )   ! zdiag zd_lw not defined/used on the halo 
     620         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     621            ! Surface value: Dirichlet 
     622            zdep    (ji,jj)   = zhsro(ji,jj) * rl_sf 
     623            psi     (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1) 
     624            zd_lw(ji,jj,1) = psi(ji,jj,1) 
     625            zd_up(ji,jj,1) = 0._wp 
     626            zdiag(ji,jj,1) = 1._wp 
     627            ! 
     628            ! Neumann condition at k=2, zdiag zd_lw not defined/used on the halo 
    544629            zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 
    545630            zd_lw(ji,jj,2) = 0._wp 
    546          END_2D 
    547          ! 
    548          ! Set psi vertical flux at the surface: 
    549          zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:) )) ! Lengh scale slope 
    550          zdep (:,:)   = ((zhsro(:,:) + gdept(:,:,1,Kmm)) / zhsro(:,:))**(rmm*ra_sf) 
    551          zflxs(:,:)   = (rnn + (1._wp-zice_fra(:,:))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:)) & 
    552             &           *(1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
    553          zdep (:,:)   = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * & 
    554             &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept(:,:,1,Kmm))**(rnn-1.) 
    555          zflxs(:,:)   = zdep(:,:) * zflxs(:,:) 
    556          psi  (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm) 
     631            ! 
     632            ! Set psi vertical flux at the surface: 
     633            zkar (ji,jj)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj) )) ! Lengh scale slope 
     634            zdep (ji,jj)   = ((zhsro(ji,jj) + gdept(ji,jj,1,Kmm)) / zhsro(ji,jj))**(rmm*ra_sf) 
     635            zflxs(ji,jj)   = (rnn + (1._wp-zice_fra(ji,jj))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(ji,jj)) & 
     636               &           *(1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1*zdep(ji,jj))**(2._wp*rmm/3._wp-1_wp) 
     637            zdep (ji,jj)   = rsbc_psi1 * (zwall_psi(ji,jj,1)*p_avm(ji,jj,1)+zwall_psi(ji,jj,2)*p_avm(ji,jj,2)) * & 
     638               &           ustar2_surf(ji,jj)**rmm * zkar(ji,jj)**rnn * (zhsro(ji,jj) + gdept(ji,jj,1,Kmm))**(rnn-1.) 
     639            zflxs(ji,jj)   = zdep(ji,jj) * zflxs(ji,jj) 
     640            psi  (ji,jj,2) = psi(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm) 
     641         END_2D 
    557642         ! 
    558643      END SELECT 
    … …  
    569654         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot 
    570655         !                      ! Balance between the production and the dissipation terms 
    571          DO_2D( 0, 0, 0, 0 ) 
     656         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    572657            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
    573658            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
    … …  
    586671         END_2D 
    587672         ! 
     673         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity) 
     674            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     675               IF ( mikt(ji,jj) > 1 ) THEN 
     676                  itop   = mikt(ji,jj)       ! k   top w-point 
     677                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one 
     678                  ! 
     679                  zdep(ji,jj) = vkarmn * r_z0_top 
     680                  psi (ji,jj,itop) = rc0**rpp * en(ji,jj,itop)**rmm *zdep(ji,jj)**rnn 
     681                  zd_lw(ji,jj,itop) = 0._wp 
     682                  zd_up(ji,jj,itop) = 0._wp 
     683                  zdiag(ji,jj,itop) = 1._wp 
     684                  ! 
     685                  ! Just above last level, Dirichlet condition again (GOTM like) 
     686                  zdep(ji,jj) = vkarmn * ( r_z0_top + e3t(ji,jj,itopp1,Kmm) ) 
     687                  psi (ji,jj,itopp1) = rc0**rpp * en(ji,jj,itop  )**rmm *zdep(ji,jj)**rnn 
     688                  zd_lw(ji,jj,itopp1) = 0._wp 
     689                  zd_up(ji,jj,itopp1) = 0._wp 
     690                  zdiag(ji,jj,itopp1) = 1._wp 
     691               END IF 
     692            END_2D 
     693         END IF 
     694         ! 
    588695      CASE ( 1 )             ! Neumman boundary condition 
    589696         ! 
    590          DO_2D( 0, 0, 0, 0 ) 
     697         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    591698            ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
    592699            ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
    … …  
    611718         END_2D 
    612719         ! 
     720         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity) 
     721            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     722               IF ( mikt(ji,jj) > 1 ) THEN 
     723                  itop   = mikt(ji,jj)       ! k   top w-point 
     724                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one 
     725                  ! 
     726                  ! Bottom level Dirichlet condition: 
     727                  zdep(ji,jj) = vkarmn * r_z0_top 
     728                  psi (ji,jj,itop) = rc0**rpp * en(ji,jj,itop)**rmm *zdep(ji,jj)**rnn 
     729                  ! 
     730                  zd_lw(ji,jj,itop) = 0._wp 
     731                  zd_up(ji,jj,itop) = 0._wp 
     732                  zdiag(ji,jj,itop) = 1._wp 
     733                  ! 
     734                  ! Just below cavity level: Neumann condition with flux 
     735                  ! injection 
     736                  zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) ! Remove zd_up from zdiag 
     737                  zd_up(ji,jj,itopp1) = 0._wp 
     738                  ! 
     739                  ! Set psi vertical flux below cavity: 
     740                  zdep(ji,jj) = r_z0_top + 0.5_wp*e3t(ji,jj,itopp1,Kmm) 
     741                  zflxb = rsbc_psi2 * ( p_avm(ji,jj,itop) + p_avm(ji,jj,itopp1))   & 
     742                     &  * (0.5_wp*(en(ji,jj,itop)+en(ji,jj,itopp1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 
     743                  psi(ji,jj,itopp1) = psi(ji,jj,itopp1) + zflxb / e3w(ji,jj,itopp1,Kmm) 
     744               END IF 
     745            END_2D 
     746         END IF 
     747 
     748         ! 
    613749      END SELECT 
    614750 
    … …  
    616752      ! ---------------- 
    617753      ! 
    618       DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
     754      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    619755         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    620756      END_3D 
    621       DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
     757      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    622758         zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 
    623759      END_3D 
    624       DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
     760      DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    625761         psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    626762      END_3D 
    … …  
    632768      ! 
    633769      CASE( 0 )               ! k-kl  (Mellor-Yamada) 
    634          DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     770         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 
    635771            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) 
    636772         END_3D 
    637773         ! 
    638774      CASE( 1 )               ! k-eps 
    639          DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     775         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 
    640776            eps(ji,jj,jk) = psi(ji,jj,jk) 
    641777         END_3D 
    642778         ! 
    643779      CASE( 2 )               ! k-w 
    644          DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     780         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 
    645781            eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 
    646782         END_3D 
    … …  
    650786         zex1  =      ( 1.5_wp + rmm/rnn ) 
    651787         zex2  = -1._wp / rnn 
    652          DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     788         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 
    653789            eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 
    654790         END_3D 
    … …  
    658794      ! Limit dissipation rate under stable stratification 
    659795      ! -------------------------------------------------- 
    660       DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! Note that this set boundary conditions on hmxl_n at the same time 
     796      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )   ! Note that this set boundary conditions on hmxl_n at the same time 
    661797         ! limitation 
    662798         eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin ) 
    663799         hmxl_n(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 
    664          ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 
    665          zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    666          IF( ln_length_lim )   hmxl_n(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 
    667       END_3D 
     800      END_3D 
     801      IF( ln_length_lim ) THEN        ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 
     802         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 
     803            zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
     804            hmxl_n(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 
     805         END_3D 
     806      ENDIF 
    668807 
    669808      ! 
    … …  
    674813      ! 
    675814      CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions 
    676          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     815         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    677816            ! zcof =  l²/q² 
    678817            zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) 
    … …  
    691830         ! 
    692831      CASE ( 2, 3 )               ! Canuto stability functions 
    693          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     832         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    694833            ! zcof =  l²/q² 
    695834            zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) 
    … …  
    723862      ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0) 
    724863      zstm(:,:,jpk) = 0. 
    725       DO_2D( 0, 0, 0, 0 )             ! update bottom with good values 
     864      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )             ! update bottom with good values 
    726865         zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 
    727866      END_2D 
    728867 
    729       zstt(:,:,  1) = wmask(:,:,  1)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0) 
    730       zstt(:,:,jpk) = wmask(:,:,jpk)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0) 
     868      zstt(:,:,  1) = wmask(A2D(nn_hls),  1)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0) 
     869      zstt(:,:,jpk) = wmask(A2D(nn_hls),jpk)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0) 
    731870 
    732871!!gm should be done for ISF (top boundary cond.) 
    … …  
    738877      !     later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk) 
    739878      !     for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0) 
    740       DO_3D( 0, 0, 0, 0, 1, jpk ) 
     879      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 
    741880         zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 
    742881         zavt  = zsqen * zstt(ji,jj,jk) 
    … …  
    745884         p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom 
    746885      END_3D 
    747       p_avt(:,:,1) = 0._wp 
     886      p_avt(A2D(nn_hls),1) = 0._wp 
    748887      ! 
    749888      IF(sn_cfctl%l_prtctl) THEN 
    … …  
    775914      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim,       & 
    776915         &            rn_clim_galp, ln_sigpsi, rn_hsro, rn_hsri,   & 
    777          &            rn_crban, rn_charn, rn_frac_hs,              & 
     916         &            nn_mxlice, rn_crban, rn_charn, rn_frac_hs,   & 
    778917         &            nn_bc_surf, nn_bc_bot, nn_z0_met, nn_z0_ice, & 
    779918         &            nn_stab_func, nn_clos 
    … …  
    815954         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos 
    816955         WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro 
    817          WRITE(numout,*) '      Ice-ocean roughness (used if nn_z0_ice/=0)    rn_hsri        = ', rn_hsri 
     956         WRITE(numout,*) '      type of scaling under sea-ice                 nn_mxlice      = ', nn_mxlice 
     957         IF( nn_mxlice == 1 ) & 
     958            WRITE(numout,*) '      Ice-ocean roughness (used if nn_z0_ice/=0) rn_hsri        = ', rn_hsri 
     959         SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice 
     960            CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   No scaling under sea-ice' 
     961            CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   scaling with constant sea-ice thickness' 
     962            CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   scaling with mean     sea-ice thickness' 
     963            CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   scaling with max      sea-ice thickness' 
     964            CASE DEFAULT 
     965               CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 ') 
     966         END SELECT 
    818967         WRITE(numout,*) 
    819968      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.