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 12377 for NEMO/trunk/src/OCE/ZDF/zdfosm.F90 – NEMO

Ignore:
Timestamp:
2020-02-12T15:39:06+01:00 (4 years ago)
Author:
acc
Message:

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk

    • Property svn:externals

      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 

  • NEMO/trunk/src/OCE/ZDF/zdfosm.F90

    r11536 r12377  
    4242   !!---------------------------------------------------------------------- 
    4343   USE oce            ! ocean dynamics and active tracers 
    44                       ! uses wn from previous time step (which is now wb) to calculate hbl 
     44                      ! uses ww from previous time step (which is now wb) to calculate hbl 
    4545   USE dom_oce        ! ocean space and time domain 
    4646   USE zdf_oce        ! ocean vertical physics 
    … …  
    103103   INTEGER :: idebug = 236 
    104104   INTEGER :: jdebug = 228 
     105   !! * Substitutions 
     106#  include "do_loop_substitute.h90" 
    105107   !!---------------------------------------------------------------------- 
    106108   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    122124 
    123125 
    124    SUBROUTINE zdf_osm( kt, p_avm, p_avt ) 
     126   SUBROUTINE zdf_osm( kt, Kbb, Kmm, Krhs, p_avm, p_avt ) 
    125127      !!---------------------------------------------------------------------- 
    126128      !!                   ***  ROUTINE zdf_osm  *** 
    … …  
    157159      !!         the equation number. (LMD94, here after) 
    158160      !!---------------------------------------------------------------------- 
    159       INTEGER                   , INTENT(in   ) ::   kt            ! ocean time step 
     161      INTEGER                   , INTENT(in   ) ::  kt             ! ocean time step 
     162      INTEGER                   , INTENT(in   ) ::  Kbb, Kmm, Krhs ! ocean time level indices 
    160163      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::  p_avm, p_avt   ! momentum and tracer Kz (w-points) 
    161164      !! 
    … …  
    295298     zz0 =       rn_abs       ! surface equi-partition in 2-bands 
    296299     zz1 =  1. - rn_abs 
    297      DO jj = 2, jpjm1 
    298         DO ji = 2, jpim1 
    299            ! Surface downward irradiance (so always +ve) 
    300            zrad0(ji,jj) = qsr(ji,jj) * r1_rau0_rcp 
    301            ! Downwards irradiance at base of boundary layer 
    302            zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) ) 
    303            ! Downwards irradiance averaged over depth of the OSBL 
    304            zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 & 
    305                  &                         + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj) 
    306         END DO 
    307      END DO 
     300     DO_2D_00_00 
     301        ! Surface downward irradiance (so always +ve) 
     302        zrad0(ji,jj) = qsr(ji,jj) * r1_rau0_rcp 
     303        ! Downwards irradiance at base of boundary layer 
     304        zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) ) 
     305        ! Downwards irradiance averaged over depth of the OSBL 
     306        zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 & 
     307              &                         + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj) 
     308     END_2D 
    308309     ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL 
    309      DO jj = 2, jpjm1 
    310         DO ji = 2, jpim1 
    311            zthermal = rab_n(ji,jj,1,jp_tem) 
    312            zbeta    = rab_n(ji,jj,1,jp_sal) 
    313            ! Upwards surface Temperature flux for non-local term 
    314            zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1) 
    315            ! Upwards surface salinity flux for non-local term 
    316            zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)  + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1) 
    317            ! Non radiative upwards surface buoyancy flux 
    318            zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) -  grav * zbeta * zws0(ji,jj) 
    319            ! turbulent heat flux averaged over depth of OSBL 
    320            zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) ) 
    321            ! turbulent salinity flux averaged over depth of the OBSL 
    322            zwsav(ji,jj) = 0.5 * zws0(ji,jj) 
    323            ! turbulent buoyancy flux averaged over the depth of the OBSBL 
    324            zwbav(ji,jj) = grav  * zthermal * zwthav(ji,jj) - grav  * zbeta * zwsav(ji,jj) 
    325            ! Surface upward velocity fluxes 
    326            zuw0(ji,jj) = -utau(ji,jj) * r1_rau0 * tmask(ji,jj,1) 
    327            zvw0(ji,jj) = -vtau(ji,jj) * r1_rau0 * tmask(ji,jj,1) 
    328            ! Friction velocity (zustar), at T-point : LMD94 eq. 2 
    329            zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) 
    330            zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 
    331            zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 
    332         END DO 
    333      END DO 
     310     DO_2D_00_00 
     311        zthermal = rab_n(ji,jj,1,jp_tem) 
     312        zbeta    = rab_n(ji,jj,1,jp_sal) 
     313        ! Upwards surface Temperature flux for non-local term 
     314        zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1) 
     315        ! Upwards surface salinity flux for non-local term 
     316        zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm)  + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1) 
     317        ! Non radiative upwards surface buoyancy flux 
     318        zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) -  grav * zbeta * zws0(ji,jj) 
     319        ! turbulent heat flux averaged over depth of OSBL 
     320        zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) ) 
     321        ! turbulent salinity flux averaged over depth of the OBSL 
     322        zwsav(ji,jj) = 0.5 * zws0(ji,jj) 
     323        ! turbulent buoyancy flux averaged over the depth of the OBSBL 
     324        zwbav(ji,jj) = grav  * zthermal * zwthav(ji,jj) - grav  * zbeta * zwsav(ji,jj) 
     325        ! Surface upward velocity fluxes 
     326        zuw0(ji,jj) = -utau(ji,jj) * r1_rau0 * tmask(ji,jj,1) 
     327        zvw0(ji,jj) = -vtau(ji,jj) * r1_rau0 * tmask(ji,jj,1) 
     328        ! Friction velocity (zustar), at T-point : LMD94 eq. 2 
     329        zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) 
     330        zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 
     331        zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) ) 
     332     END_2D 
    334333     ! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes) 
    335334     SELECT CASE (nn_osm_wave) 
    336335     ! Assume constant La#=0.3 
    337336     CASE(0) 
    338         DO jj = 2, jpjm1 
    339            DO ji = 2, jpim1 
    340               zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2 
    341               zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 
    342               zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 ) 
    343               ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init 
    344            END DO 
    345         END DO 
     337        DO_2D_00_00 
     338           zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2 
     339           zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 
     340           zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 ) 
     341           ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init 
     342        END_2D 
    346343     ! Assume Pierson-Moskovitz wind-wave spectrum 
    347344     CASE(1) 
    348         DO jj = 2, jpjm1 
    349            DO ji = 2, jpim1 
    350               ! Use wind speed wndm included in sbc_oce module 
    351               zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 
    352               dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav 
    353            END DO 
    354         END DO 
     345        DO_2D_00_00 
     346           ! Use wind speed wndm included in sbc_oce module 
     347           zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 
     348           dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav 
     349        END_2D 
    355350     ! Use ECMWF wave fields as output from SBCWAVE 
    356351     CASE(2) 
    357352        zfac =  2.0_wp * rpi / 16.0_wp 
    358         DO jj = 2, jpjm1 
    359            DO ji = 2, jpim1 
    360               ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas. 
    361               !    The coefficient 0.8 gives La=0.3  in this situation. 
    362               ! It could represent the effects of the spread of wave directions 
    363               ! around the mean wind. The effect of this adjustment needs to be tested. 
    364               zustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), & 
    365                    &                zustar(ji,jj) / ( 0.45 * 0.45 )                                                  ) 
    366               dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 
    367            END DO 
    368         END DO 
     353        DO_2D_00_00 
     354           ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas. 
     355           !    The coefficient 0.8 gives La=0.3  in this situation. 
     356           ! It could represent the effects of the spread of wave directions 
     357           ! around the mean wind. The effect of this adjustment needs to be tested. 
     358           zustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), & 
     359                &                zustar(ji,jj) / ( 0.45 * 0.45 )                                                  ) 
     360           dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes ! 
     361        END_2D 
    369362     END SELECT 
    370363 
    371364     ! Langmuir velocity scale (zwstrl), La # (zla) 
    372365     ! mixed scale (zvstr), convective velocity scale (zwstrc) 
    373      DO jj = 2, jpjm1 
    374         DO ji = 2, jpim1 
    375            ! Langmuir velocity scale (zwstrl), at T-point 
    376            zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird 
    377            ! Modify zwstrl to allow for small and large values of dstokes/hbl. 
    378            ! Intended as a possible test. Doesn't affect LES results for entrainment, 
    379            !  but hasn't been shown to be correct as dstokes/h becomes large or small. 
    380            zwstrl(ji,jj) = zwstrl(ji,jj) *  & 
    381                 & (1.12 * ( 1.0 - ( 1.0 - EXP( -hbl(ji,jj) / dstokes(ji,jj) ) ) * dstokes(ji,jj) / hbl(ji,jj) ))**pthird * & 
    382                 & ( 1.0 - EXP( -15.0 * dstokes(ji,jj) / hbl(ji,jj) )) 
    383            ! define La this way so effects of Stokes penetration depth on velocity scale are included 
    384            zla(ji,jj) = SQRT ( zustar(ji,jj) / zwstrl(ji,jj) )**3 
    385            ! Velocity scale that tends to zustar for large Langmuir numbers 
    386            zvstr(ji,jj) = ( zwstrl(ji,jj)**3  + & 
    387                 & ( 1.0 - EXP( -0.5 * zla(ji,jj)**2 ) ) * zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird 
    388  
    389            ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 
    390            ! Note zustke and zwstrl are not amended. 
    391            IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45 
    392            ! 
    393            ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv 
    394            IF ( zwbav(ji,jj) > 0.0) THEN 
    395               zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird 
    396               zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln ) 
    397               lconv(ji,jj) = .TRUE. 
    398            ELSE 
    399               zhol(ji,jj) = -hbl(ji,jj) *  2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3  + epsln ) 
    400               lconv(ji,jj) = .FALSE. 
    401            ENDIF 
    402         END DO 
    403      END DO 
     366     DO_2D_00_00 
     367        ! Langmuir velocity scale (zwstrl), at T-point 
     368        zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird 
     369        ! Modify zwstrl to allow for small and large values of dstokes/hbl. 
     370        ! Intended as a possible test. Doesn't affect LES results for entrainment, 
     371        !  but hasn't been shown to be correct as dstokes/h becomes large or small. 
     372        zwstrl(ji,jj) = zwstrl(ji,jj) *  & 
     373             & (1.12 * ( 1.0 - ( 1.0 - EXP( -hbl(ji,jj) / dstokes(ji,jj) ) ) * dstokes(ji,jj) / hbl(ji,jj) ))**pthird * & 
     374             & ( 1.0 - EXP( -15.0 * dstokes(ji,jj) / hbl(ji,jj) )) 
     375        ! define La this way so effects of Stokes penetration depth on velocity scale are included 
     376        zla(ji,jj) = SQRT ( zustar(ji,jj) / zwstrl(ji,jj) )**3 
     377        ! Velocity scale that tends to zustar for large Langmuir numbers 
     378        zvstr(ji,jj) = ( zwstrl(ji,jj)**3  + & 
     379             & ( 1.0 - EXP( -0.5 * zla(ji,jj)**2 ) ) * zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird 
     380 
     381        ! limit maximum value of Langmuir number as approximate treatment for shear turbulence. 
     382        ! Note zustke and zwstrl are not amended. 
     383        IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45 
     384        ! 
     385        ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv 
     386        IF ( zwbav(ji,jj) > 0.0) THEN 
     387           zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird 
     388           zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln ) 
     389           lconv(ji,jj) = .TRUE. 
     390        ELSE 
     391           zhol(ji,jj) = -hbl(ji,jj) *  2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3  + epsln ) 
     392           lconv(ji,jj) = .FALSE. 
     393        ENDIF 
     394     END_2D 
    404395 
    405396     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    407398     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    408399     ! BL must be always 2 levels deep. 
    409       hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,3) ) 
     400      hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,3,Kmm) ) 
    410401      ibld(:,:) = 3 
    411       DO jk = 4, jpkm1 
    412          DO jj = 2, jpjm1 
    413             DO ji = 2, jpim1 
    414                IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 
    415                   ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 
    416                ENDIF 
     402      DO_3D_00_00( 4, jpkm1 ) 
     403         IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 
     404            ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 
     405         ENDIF 
     406      END_3D 
     407 
     408      DO_2D_00_00 
     409            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
     410            zbeta    = rab_n(ji,jj,1,jp_sal) 
     411            zt   = 0._wp 
     412            zs   = 0._wp 
     413            zu   = 0._wp 
     414            zv   = 0._wp 
     415            ! average over depth of boundary layer 
     416            zthick=0._wp 
     417            DO jm = 2, ibld(ji,jj) 
     418               zthick=zthick+e3t(ji,jj,jm,Kmm) 
     419               zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 
     420               zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 
     421               zu   = zu  + e3t(ji,jj,jm,Kmm) & 
     422                  &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 
     423                  &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
     424               zv   = zv  + e3t(ji,jj,jm,Kmm) & 
     425                  &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 
     426                  &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    417427            END DO 
    418          END DO 
    419       END DO 
    420  
    421       DO jj = 2, jpjm1                                 !  Vertical slab 
    422          DO ji = 2, jpim1 
    423                zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    424                zbeta    = rab_n(ji,jj,1,jp_sal) 
    425                zt   = 0._wp 
    426                zs   = 0._wp 
    427                zu   = 0._wp 
    428                zv   = 0._wp 
    429                ! average over depth of boundary layer 
    430                zthick=0._wp 
    431                DO jm = 2, ibld(ji,jj) 
    432                   zthick=zthick+e3t_n(ji,jj,jm) 
    433                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    434                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    435                   zu   = zu  + e3t_n(ji,jj,jm) & 
    436                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    437                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    438                   zv   = zv  + e3t_n(ji,jj,jm) & 
    439                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    440                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    441                END DO 
    442                zt_bl(ji,jj) = zt / zthick 
    443                zs_bl(ji,jj) = zs / zthick 
    444                zu_bl(ji,jj) = zu / zthick 
    445                zv_bl(ji,jj) = zv / zthick 
    446                zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    447                zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    448                zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    449                      &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    450                zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    451                      &   / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    452                zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 
    453                IF ( lconv(ji,jj) ) THEN    ! Convective 
    454                       zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 
    455                            &            + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 
    456  
    457                       zvel_max =  - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 
    458                            &   * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     428            zt_bl(ji,jj) = zt / zthick 
     429            zs_bl(ji,jj) = zs / zthick 
     430            zu_bl(ji,jj) = zu / zthick 
     431            zv_bl(ji,jj) = zv / zthick 
     432            zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 
     433            zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 
     434            zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 
     435                  &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
     436            zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 
     437                  &   / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
     438            zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 
     439            IF ( lconv(ji,jj) ) THEN    ! Convective 
     440                   zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 
     441                        &            + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 
     442 
     443                   zvel_max =  - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 
     444                        &   * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    459445! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment. 
    460446!                      zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) & 
    … …  
    463449!                      zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 
    464450!                           &       ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    465                       zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) 
    466                ELSE                        ! Stable 
    467                       zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) & 
    468                            &   + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) & 
    469                            & - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) & 
    470                            &  * zwstrl(ji,jj)**3 / hbli(ji,jj) 
    471                       zzdhdt = zzdhdt + zwbav(ji,jj) 
    472                       IF ( zzdhdt < 0._wp ) THEN 
    473                       ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
    474                          zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 
    475                       ELSE 
    476                          zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 
    477                               &  + MAX( zdb_bl(ji,jj), 0.0 ) 
    478                       ENDIF 
    479                       zzdhdt = 2.0 * zzdhdt / zpert 
    480                ENDIF 
    481                zdhdt(ji,jj) = zzdhdt 
    482            END DO 
    483       END DO 
     451                   zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) 
     452            ELSE                        ! Stable 
     453                   zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) & 
     454                        &   + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) & 
     455                        & - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) & 
     456                        &  * zwstrl(ji,jj)**3 / hbli(ji,jj) 
     457                   zzdhdt = zzdhdt + zwbav(ji,jj) 
     458                   IF ( zzdhdt < 0._wp ) THEN 
     459                   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
     460                      zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 
     461                   ELSE 
     462                      zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 
     463                           &  + MAX( zdb_bl(ji,jj), 0.0 ) 
     464                   ENDIF 
     465                   zzdhdt = 2.0 * zzdhdt / zpert 
     466            ENDIF 
     467            zdhdt(ji,jj) = zzdhdt 
     468      END_2D 
    484469 
    485470      ! Calculate averages over depth of boundary layer 
    … …  
    487472      ibld(:,:) = 3 
    488473 
    489       zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 
    490       zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_n(:,:)) 
    491       zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + wn(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 
    492  
    493       DO jk = 4, jpkm1 
    494          DO jj = 2, jpjm1 
    495             DO ji = 2, jpim1 
    496                IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 
    497                   ibld(ji,jj) =  MIN(mbkt(ji,jj), jk) 
    498                ENDIF 
    499             END DO 
    500          END DO 
    501       END DO 
     474      zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 
     475      zhbl_t(:,:) = MIN(zhbl_t(:,:), ht(:,:)) 
     476      zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 
     477 
     478      DO_3D_00_00( 4, jpkm1 ) 
     479         IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 
     480            ibld(ji,jj) =  MIN(mbkt(ji,jj), jk) 
     481         ENDIF 
     482      END_3D 
    502483 
    503484! 
    504485! Step through model levels taking account of buoyancy change to determine the effect on dhdt 
    505486! 
    506       DO jj = 2, jpjm1 
    507          DO ji = 2, jpim1 
    508             IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 
     487      DO_2D_00_00 
     488         IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 
    509489! 
    510490! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 
    511491! 
    512                zhbl_s = hbl(ji,jj) 
    513                jm = imld(ji,jj) 
    514                zthermal = rab_n(ji,jj,1,jp_tem) 
    515                zbeta = rab_n(ji,jj,1,jp_sal) 
    516                IF ( lconv(ji,jj) ) THEN 
     492            zhbl_s = hbl(ji,jj) 
     493            jm = imld(ji,jj) 
     494            zthermal = rab_n(ji,jj,1,jp_tem) 
     495            zbeta = rab_n(ji,jj,1,jp_sal) 
     496            IF ( lconv(ji,jj) ) THEN 
    517497!unstable 
    518                   zvel_max =  - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 
    519                        &   * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    520  
    521                   DO jk = imld(ji,jj), ibld(ji,jj) 
    522                      zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 
    523                           & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max 
    524  
    525                      zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) ) 
    526                      zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 
    527  
    528                      IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 
    529                   END DO 
    530                   hbl(ji,jj) = zhbl_s 
    531                   ibld(ji,jj) = jm 
    532                   hbli(ji,jj) = hbl(ji,jj) 
    533                ELSE 
     498               zvel_max =  - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 
     499                    &   * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     500 
     501               DO jk = imld(ji,jj), ibld(ji,jj) 
     502                  zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) & 
     503                       & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) + zvel_max 
     504 
     505                  zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w(ji,jj,jk,Kmm) ) 
     506                  zhbl_s = MIN(zhbl_s, ht(ji,jj)) 
     507 
     508                  IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1 
     509               END DO 
     510               hbl(ji,jj) = zhbl_s 
     511               ibld(ji,jj) = jm 
     512               hbli(ji,jj) = hbl(ji,jj) 
     513            ELSE 
    534514! stable 
    535                   DO jk = imld(ji,jj), ibld(ji,jj) 
    536                      zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )          & 
    537                           &               - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) & 
    538                           & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s 
    539  
    540                      zhbl_s = zhbl_s +  (                                                                                & 
    541                           &                     0.32         *                         ( hbli(ji,jj) / zhbl_s -1.0 )     & 
    542                           &               * zwstrl(ji,jj)**3 / hbli(ji,jj)                                               & 
    543                           &               + ( ( 0.32 / 3.0 )           * EXP( -  2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) )   & 
    544                           &               -   ( 0.32 / 3.0  - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s      ) ) ) & 
    545                           &          * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w_n(ji,jj,jk) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include wn here 
    546  
    547                      zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 
    548                      IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 
    549                   END DO 
    550                   hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,3) ) 
    551                   ibld(ji,jj) = MAX(jm, 3 ) 
    552                   IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 
    553                ENDIF   ! IF ( lconv ) 
     515               DO jk = imld(ji,jj), ibld(ji,jj) 
     516                  zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) )          & 
     517                       &               - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) & 
     518                       & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s 
     519 
     520                  zhbl_s = zhbl_s +  (                                                                                & 
     521                       &                     0.32         *                         ( hbli(ji,jj) / zhbl_s -1.0 )     & 
     522                       &               * zwstrl(ji,jj)**3 / hbli(ji,jj)                                               & 
     523                       &               + ( ( 0.32 / 3.0 )           * EXP( -  2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) )   & 
     524                       &               -   ( 0.32 / 3.0  - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s      ) ) ) & 
     525                       &          * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w(ji,jj,jk,Kmm) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include ww here 
     526 
     527                  zhbl_s = MIN(zhbl_s, ht(ji,jj)) 
     528                  IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1 
     529               END DO 
     530               hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,3,Kmm) ) 
     531               ibld(ji,jj) = MAX(jm, 3 ) 
     532               IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 
     533            ENDIF   ! IF ( lconv ) 
     534         ELSE 
     535! change zero or one model level. 
     536            hbl(ji,jj) = zhbl_t(ji,jj) 
     537            IF ( lconv(ji,jj) ) THEN 
     538               hbli(ji,jj) = hbl(ji,jj) 
    554539            ELSE 
    555 ! change zero or one model level. 
    556                hbl(ji,jj) = zhbl_t(ji,jj) 
    557                IF ( lconv(ji,jj) ) THEN 
    558                   hbli(ji,jj) = hbl(ji,jj) 
    559                ELSE 
    560                   hbl(ji,jj) = MAX(hbl(ji,jj), gdepw_n(ji,jj,3) ) 
    561                   IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 
    562                ENDIF 
     540               hbl(ji,jj) = MAX(hbl(ji,jj), gdepw(ji,jj,3,Kmm) ) 
     541               IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 
    563542            ENDIF 
    564             zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
    565          END DO 
    566       END DO 
     543         ENDIF 
     544         zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 
     545      END_2D 
    567546      dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. )  !  Limit delta for shallow boundary layers for calculating flux-gradient terms. 
    568547 
    … …  
    570549     ! Consider later  combining this into the loop above and looking for columns 
    571550     ! where the index for base of the boundary layer have changed 
    572       DO jj = 2, jpjm1                                 !  Vertical slab 
    573          DO ji = 2, jpim1 
    574                zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    575                zbeta    = rab_n(ji,jj,1,jp_sal) 
    576                zt   = 0._wp 
    577                zs   = 0._wp 
    578                zu   = 0._wp 
    579                zv   = 0._wp 
    580                ! average over depth of boundary layer 
    581                zthick=0._wp 
    582                DO jm = 2, ibld(ji,jj) 
    583                   zthick=zthick+e3t_n(ji,jj,jm) 
    584                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    585                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    586                   zu   = zu  + e3t_n(ji,jj,jm) & 
    587                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    588                      &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    589                   zv   = zv  + e3t_n(ji,jj,jm) & 
    590                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    591                      &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    592                END DO 
    593                zt_bl(ji,jj) = zt / zthick 
    594                zs_bl(ji,jj) = zs / zthick 
    595                zu_bl(ji,jj) = zu / zthick 
    596                zv_bl(ji,jj) = zv / zthick 
    597                zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    598                zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    599                zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    600                       &   / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    601                zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    602                       &  / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    603                zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 
    604                zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
    605                IF ( lconv(ji,jj) ) THEN 
    606                   IF ( zdb_bl(ji,jj) > 0._wp )THEN 
    607                      IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability 
    608                            zari = 4.5 * ( zvstr(ji,jj)**2 ) & 
    609                              & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 
    610                      ELSE                                                     ! unstable 
    611                            zari = 4.5 * ( zwstrc(ji,jj)**2 ) & 
    612                              & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 
    613                      ENDIF 
    614                      IF ( zari > 0.2 ) THEN                                                ! This test checks for weakly stratified pycnocline 
    615                         zari = 0.2 
    616                         zwb_ent(ji,jj) = 0._wp 
    617                      ENDIF 
    618                      inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
     551      DO_2D_00_00 
     552            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
     553            zbeta    = rab_n(ji,jj,1,jp_sal) 
     554            zt   = 0._wp 
     555            zs   = 0._wp 
     556            zu   = 0._wp 
     557            zv   = 0._wp 
     558            ! average over depth of boundary layer 
     559            zthick=0._wp 
     560            DO jm = 2, ibld(ji,jj) 
     561               zthick=zthick+e3t(ji,jj,jm,Kmm) 
     562               zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 
     563               zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 
     564               zu   = zu  + e3t(ji,jj,jm,Kmm) & 
     565                  &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 
     566                  &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
     567               zv   = zv  + e3t(ji,jj,jm,Kmm) & 
     568                  &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 
     569                  &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
     570            END DO 
     571            zt_bl(ji,jj) = zt / zthick 
     572            zs_bl(ji,jj) = zs / zthick 
     573            zu_bl(ji,jj) = zu / zthick 
     574            zv_bl(ji,jj) = zv / zthick 
     575            zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 
     576            zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 
     577            zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 
     578                   &   / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
     579            zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 
     580                   &  / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
     581            zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 
     582            zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 
     583            IF ( lconv(ji,jj) ) THEN 
     584               IF ( zdb_bl(ji,jj) > 0._wp )THEN 
     585                  IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability 
     586                        zari = 4.5 * ( zvstr(ji,jj)**2 ) & 
     587                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 
     588                  ELSE                                                     ! unstable 
     589                        zari = 4.5 * ( zwstrc(ji,jj)**2 ) & 
     590                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 
     591                  ENDIF 
     592                  IF ( zari > 0.2 ) THEN                                                ! This test checks for weakly stratified pycnocline 
     593                     zari = 0.2 
     594                     zwb_ent(ji,jj) = 0._wp 
     595                  ENDIF 
     596                  inhml = MAX( INT( zari * zhbl(ji,jj) / e3t(ji,jj,ibld(ji,jj),Kmm) ) , 1 ) 
     597                  imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 
     598                  zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
     599                  zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
     600               ELSE  ! IF (zdb_bl) 
     601                  imld(ji,jj) = ibld(ji,jj) - 1 
     602                  zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
     603                  zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
     604               ENDIF 
     605            ELSE   ! IF (lconv) 
     606               IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here 
     607               ! boundary layer deepening 
     608                  IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     609               ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
     610                     zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
     611                       & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 ) 
     612                     inhml = MAX( INT( zari * zhbl(ji,jj) / e3t(ji,jj,ibld(ji,jj),Kmm) ) , 1 ) 
    619613                     imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 
    620                      zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     614                     zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
    621615                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    622                   ELSE  ! IF (zdb_bl) 
     616                  ELSE 
    623617                     imld(ji,jj) = ibld(ji,jj) - 1 
    624                      zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     618                     zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
    625619                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    626                   ENDIF 
    627                ELSE   ! IF (lconv) 
    628                   IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here 
    629                   ! boundary layer deepening 
    630                      IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    631                   ! pycnocline thickness set by stratification - use same relationship as for neutral conditions. 
    632                         zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
    633                           & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 ) 
    634                         inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
    635                         imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 
    636                         zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    637                         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    638                      ELSE 
    639                         imld(ji,jj) = ibld(ji,jj) - 1 
    640                         zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
    641                         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    642                      ENDIF ! IF (zdb_bl > 0.0) 
    643                   ELSE     ! IF(dhdt >= 0) 
    644                   ! boundary layer collapsing. 
    645                      imld(ji,jj) = ibld(ji,jj) 
    646                      zhml(ji,jj) = zhbl(ji,jj) 
    647                      zdh(ji,jj) = 0._wp 
    648                   ENDIF    ! IF (dhdt >= 0) 
    649                ENDIF       ! IF (lconv) 
    650          END DO 
    651       END DO 
     620                  ENDIF ! IF (zdb_bl > 0.0) 
     621               ELSE     ! IF(dhdt >= 0) 
     622               ! boundary layer collapsing. 
     623                  imld(ji,jj) = ibld(ji,jj) 
     624                  zhml(ji,jj) = zhbl(ji,jj) 
     625                  zdh(ji,jj) = 0._wp 
     626               ENDIF    ! IF (dhdt >= 0) 
     627            ENDIF       ! IF (lconv) 
     628      END_2D 
    652629 
    653630      ! Average over the depth of the mixed layer in the convective boundary layer 
    654631      ! Also calculate entrainment fluxes for temperature and salinity 
    655       DO jj = 2, jpjm1                                 !  Vertical slab 
    656          DO ji = 2, jpim1 
    657             zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    658             zbeta    = rab_n(ji,jj,1,jp_sal) 
    659             IF ( lconv(ji,jj) ) THEN 
     632      DO_2D_00_00 
     633         zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
     634         zbeta    = rab_n(ji,jj,1,jp_sal) 
     635         IF ( lconv(ji,jj) ) THEN 
     636            zt   = 0._wp 
     637            zs   = 0._wp 
     638            zu   = 0._wp 
     639            zv   = 0._wp 
     640            ! average over depth of boundary layer 
     641            zthick=0._wp 
     642            DO jm = 2, imld(ji,jj) 
     643               zthick=zthick+e3t(ji,jj,jm,Kmm) 
     644               zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 
     645               zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 
     646               zu   = zu  + e3t(ji,jj,jm,Kmm) & 
     647                  &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 
     648                  &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
     649               zv   = zv  + e3t(ji,jj,jm,Kmm) & 
     650                  &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 
     651                  &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
     652            END DO 
     653            zt_ml(ji,jj) = zt / zthick 
     654            zs_ml(ji,jj) = zs / zthick 
     655            zu_ml(ji,jj) = zu / zthick 
     656            zv_ml(ji,jj) = zv / zthick 
     657            zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 
     658            zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 
     659            zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 
     660                  &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
     661            zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 
     662                  &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
     663            zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 
     664         ELSE 
     665         ! stable, if entraining calulate average below interface layer. 
     666            IF ( zdhdt(ji,jj) >= 0._wp ) THEN 
    660667               zt   = 0._wp 
    661668               zs   = 0._wp 
    … …  
    665672               zthick=0._wp 
    666673               DO jm = 2, imld(ji,jj) 
    667                   zthick=zthick+e3t_n(ji,jj,jm) 
    668                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    669                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    670                   zu   = zu  + e3t_n(ji,jj,jm) & 
    671                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
     674                  zthick=zthick+e3t(ji,jj,jm,Kmm) 
     675                  zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 
     676                  zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 
     677                  zu   = zu  + e3t(ji,jj,jm,Kmm) & 
     678                     &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 
    672679                     &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    673                   zv   = zv  + e3t_n(ji,jj,jm) & 
    674                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
     680                  zv   = zv  + e3t(ji,jj,jm,Kmm) & 
     681                     &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 
    675682                     &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    676683               END DO 
    … …  
    679686               zu_ml(ji,jj) = zu / zthick 
    680687               zv_ml(ji,jj) = zv / zthick 
    681                zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    682                zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    683                zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
     688               zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 
     689               zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 
     690               zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 
    684691                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    685                zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
     692               zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 
    686693                     &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    687694               zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 
    688             ELSE 
    689             ! stable, if entraining calulate average below interface layer. 
    690                IF ( zdhdt(ji,jj) >= 0._wp ) THEN 
    691                   zt   = 0._wp 
    692                   zs   = 0._wp 
    693                   zu   = 0._wp 
    694                   zv   = 0._wp 
    695                   ! average over depth of boundary layer 
    696                   zthick=0._wp 
    697                   DO jm = 2, imld(ji,jj) 
    698                      zthick=zthick+e3t_n(ji,jj,jm) 
    699                      zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    700                      zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    701                      zu   = zu  + e3t_n(ji,jj,jm) & 
    702                         &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
    703                         &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    704                      zv   = zv  + e3t_n(ji,jj,jm) & 
    705                         &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
    706                         &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    707                   END DO 
    708                   zt_ml(ji,jj) = zt / zthick 
    709                   zs_ml(ji,jj) = zs / zthick 
    710                   zu_ml(ji,jj) = zu / zthick 
    711                   zv_ml(ji,jj) = zv / zthick 
    712                   zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    713                   zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    714                   zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
    715                         &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    716                   zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
    717                         &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    718                   zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 
    719                ENDIF 
    720695            ENDIF 
    721          END DO 
    722       END DO 
     696         ENDIF 
     697      END_2D 
    723698    ! 
    724699    ! rotate mean currents and changes onto wind align co-ordinates 
    725700    ! 
    726701 
    727       DO jj = 2, jpjm1 
    728          DO ji = 2, jpim1 
    729             ztemp = zu_ml(ji,jj) 
    730             zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj) 
    731             zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    732             ztemp = zdu_ml(ji,jj) 
    733             zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj) 
    734             zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    735     ! 
    736             ztemp = zu_bl(ji,jj) 
    737             zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj) 
    738             zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    739             ztemp = zdu_bl(ji,jj) 
    740             zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj) 
    741             zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
    742          END DO 
    743       END DO 
     702      DO_2D_00_00 
     703         ztemp = zu_ml(ji,jj) 
     704         zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj) 
     705         zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
     706         ztemp = zdu_ml(ji,jj) 
     707         zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj) 
     708         zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
     709 ! 
     710         ztemp = zu_bl(ji,jj) 
     711         zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj) 
     712         zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
     713         ztemp = zdu_bl(ji,jj) 
     714         zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj) 
     715         zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj) 
     716      END_2D 
    744717 
    745718     zuw_bse = 0._wp 
    746719     zvw_bse = 0._wp 
    747      DO jj = 2, jpjm1 
    748         DO ji = 2, jpim1 
    749  
    750            IF ( lconv(ji,jj) ) THEN 
    751               IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    752                  zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
    753                  zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
    754               ENDIF 
    755            ELSE 
    756               zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 
    757               zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 
     720     DO_2D_00_00 
     721 
     722        IF ( lconv(ji,jj) ) THEN 
     723           IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     724              zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
     725              zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln) 
    758726           ENDIF 
    759         END DO 
    760      END DO 
     727        ELSE 
     728           zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) ) 
     729           zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) ) 
     730        ENDIF 
     731     END_2D 
    761732 
    762733      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    764735      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    765736 
    766        DO jj = 2, jpjm1 
    767           DO ji = 2, jpim1 
    768           ! 
    769              IF ( lconv (ji,jj) ) THEN 
    770              ! Unstable conditions 
    771                 IF( zdb_bl(ji,jj) > 0._wp ) THEN 
    772                 ! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero 
    773                    ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) ) 
    774                    zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) ) 
    775                    zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) ) 
    776                    DO jk = 2 , ibld(ji,jj) 
    777                       znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    778                       zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    779                       zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    780                       zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    781                    END DO 
    782                 ENDIF 
    783              ELSE 
    784              ! stable conditions 
    785              ! if pycnocline profile only defined when depth steady of increasing. 
    786                 IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady. 
    787                    IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
    788                      IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline 
    789                          ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 
    790                          zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 
    791                          zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 
    792                          DO jk = 2, ibld(ji,jj) 
    793                             znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    794                             zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    795                             zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    796                             zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    797                          END DO 
    798                      ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
    799                          ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 
    800                          zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 
    801                          zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 
    802                          DO jk = 2, ibld(ji,jj) 
    803                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    804                             zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    805                             zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    806                             zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    807                          END DO 
    808                       ENDIF ! IF (zhol >=0.5) 
    809                    ENDIF    ! IF (zdb_bl> 0.) 
    810                 ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero 
    811              ENDIF          ! IF (lconv) 
    812             ! 
    813           END DO 
    814        END DO 
    815 ! 
    816        DO jj = 2, jpjm1 
    817           DO ji = 2, jpim1 
    818           ! 
    819              IF ( lconv (ji,jj) ) THEN 
    820              ! Unstable conditions 
    821                  zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / & 
    822                & (( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0) 
    823                 zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / & 
    824               & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    825                 DO jk = 2 , ibld(ji,jj)-1 
    826                    znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    827                    zdudz_pyc(ji,jj,jk) =  zugrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    828                    zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    829                 END DO 
    830              ELSE 
    831              ! stable conditions 
    832                 zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 
    833                 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 
    834                 DO jk = 2, ibld(ji,jj) 
    835                    znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    836                    IF ( znd < 1.0 ) THEN 
    837                       zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 
    838                    ELSE 
    839                       zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 
    840                    ENDIF 
    841                    zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 
     737       DO_2D_00_00 
     738       ! 
     739          IF ( lconv (ji,jj) ) THEN 
     740          ! Unstable conditions 
     741             IF( zdb_bl(ji,jj) > 0._wp ) THEN 
     742             ! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero 
     743                ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) ) 
     744                zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) ) 
     745                zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) ) 
     746                DO jk = 2 , ibld(ji,jj) 
     747                   znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
     748                   zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     749                   zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     750                   zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    842751                END DO 
    843752             ENDIF 
    844             ! 
    845           END DO 
    846        END DO 
     753          ELSE 
     754          ! stable conditions 
     755          ! if pycnocline profile only defined when depth steady of increasing. 
     756             IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady. 
     757                IF ( zdb_bl(ji,jj) > 0._wp ) THEN 
     758                  IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline 
     759                      ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj) 
     760                      zsgrad = zds_bl(ji,jj) / zhbl(ji,jj) 
     761                      zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 
     762                      DO jk = 2, ibld(ji,jj) 
     763                         znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
     764                         zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     765                         zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     766                         zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
     767                      END DO 
     768                  ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form. 
     769                      ztgrad = zdt_bl(ji,jj) / zdh(ji,jj) 
     770                      zsgrad = zds_bl(ji,jj) / zdh(ji,jj) 
     771                      zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 
     772                      DO jk = 2, ibld(ji,jj) 
     773                         znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
     774                         zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     775                         zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     776                         zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     777                      END DO 
     778                   ENDIF ! IF (zhol >=0.5) 
     779                ENDIF    ! IF (zdb_bl> 0.) 
     780             ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero 
     781          ENDIF          ! IF (lconv) 
     782         ! 
     783       END_2D 
     784! 
     785       DO_2D_00_00 
     786       ! 
     787          IF ( lconv (ji,jj) ) THEN 
     788          ! Unstable conditions 
     789              zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / & 
     790            & (( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0) 
     791             zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / & 
     792           & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     793             DO jk = 2 , ibld(ji,jj)-1 
     794                znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
     795                zdudz_pyc(ji,jj,jk) =  zugrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     796                zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
     797             END DO 
     798          ELSE 
     799          ! stable conditions 
     800             zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj) 
     801             zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 
     802             DO jk = 2, ibld(ji,jj) 
     803                znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
     804                IF ( znd < 1.0 ) THEN 
     805                   zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 
     806                ELSE 
     807                   zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 ) 
     808                ENDIF 
     809                zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 ) 
     810             END DO 
     811          ENDIF 
     812         ! 
     813       END_2D 
    847814       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    848815       ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship 
    … …  
    860827      !     zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 
    861828      !  ENDWHERE 
    862        DO jj = 2, jpjm1 
    863           DO ji = 2, jpim1 
    864              IF ( lconv(ji,jj) ) THEN 
    865                zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    866                zvisml_sc(ji,jj) = zdifml_sc(ji,jj) 
    867                zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 
    868                zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 
    869                zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 
    870                zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj) 
    871              ELSE 
    872                zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 
    873                zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 
    874             END IF 
    875         END DO 
    876     END DO 
     829       DO_2D_00_00 
     830          IF ( lconv(ji,jj) ) THEN 
     831            zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
     832            zvisml_sc(ji,jj) = zdifml_sc(ji,jj) 
     833            zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 
     834            zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj) 
     835            zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third 
     836            zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj) 
     837          ELSE 
     838            zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 
     839            zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ) 
     840         END IF 
     841       END_2D 
    877842! 
    878        DO jj = 2, jpjm1 
    879           DO ji = 2, jpim1 
    880              IF ( lconv(ji,jj) ) THEN 
    881                 DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity 
    882                     zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 
     843       DO_2D_00_00 
     844          IF ( lconv(ji,jj) ) THEN 
     845             DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity 
     846                 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 
     847                 ! 
     848                 zdiffut(ji,jj,jk) = 0.8   * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml    )**1.5 
     849                 ! 
     850                 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml    ) & 
     851                      &            *                                      ( 1.0 -               0.5 * zznd_ml**2 ) 
     852             END DO 
     853             ! pycnocline - if present linear profile 
     854             IF ( zdh(ji,jj) > 0._wp ) THEN 
     855                DO jk = imld(ji,jj)+1 , ibld(ji,jj) 
     856                    zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
    883857                    ! 
    884                     zdiffut(ji,jj,jk) = 0.8   * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml    )**1.5 
     858                    zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
    885859                    ! 
    886                     zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml    ) & 
    887                          &            *                                      ( 1.0 -               0.5 * zznd_ml**2 ) 
     860                    zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
    888861                END DO 
    889                 ! pycnocline - if present linear profile 
    890                 IF ( zdh(ji,jj) > 0._wp ) THEN 
    891                    DO jk = imld(ji,jj)+1 , ibld(ji,jj) 
    892                        zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    893                        ! 
    894                        zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
    895                        ! 
    896                        zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
    897                    END DO 
    898                 ENDIF 
    899                 ! Temporay fix to ensure zdiffut is +ve; won't be necessary with wn taken out 
    900                 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t_n(ji,jj,ibld(ji,jj)) 
    901                 ! could be taken out, take account of entrainment represents as a diffusivity 
    902                 ! should remove w from here, represents entrainment 
    903              ELSE 
    904              ! stable conditions 
    905                 DO jk = 2, ibld(ji,jj) 
    906                    zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    907                    zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 
    908                    zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 
    909                 END DO 
    910              ENDIF   ! end if ( lconv ) 
     862             ENDIF 
     863             ! Temporay fix to ensure zdiffut is +ve; won't be necessary with ww taken out 
     864             zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t(ji,jj,ibld(ji,jj),Kmm) 
     865             ! could be taken out, take account of entrainment represents as a diffusivity 
     866             ! should remove w from here, represents entrainment 
     867          ELSE 
     868          ! stable conditions 
     869             DO jk = 2, ibld(ji,jj) 
     870                zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
     871                zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 
     872                zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 
     873             END DO 
     874          ENDIF   ! end if ( lconv ) 
    911875! 
    912           END DO  ! end of ji loop 
    913        END DO     ! end of jj loop 
     876       END_2D 
    914877 
    915878       ! 
    … …  
    928891 
    929892 
    930        DO jj = 2, jpjm1 
    931           DO ji = 2, jpim1 
    932             IF ( lconv(ji,jj) ) THEN 
    933               DO jk = 2, imld(ji,jj) 
    934                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
    935                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 
    936                  ! 
    937                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) *  zsc_ws_1(ji,jj) 
    938               END DO ! end jk loop 
    939             ELSE     ! else for if (lconv) 
     893       DO_2D_00_00 
     894         IF ( lconv(ji,jj) ) THEN 
     895           DO jk = 2, imld(ji,jj) 
     896              zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     897              ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 
     898              ! 
     899              ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) *  zsc_ws_1(ji,jj) 
     900           END DO ! end jk loop 
     901         ELSE     ! else for if (lconv) 
    940902 ! Stable conditions 
    941                DO jk = 2, ibld(ji,jj) 
    942                   zznd_d=gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
    943                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 
    944                        &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 
    945                   ! 
    946                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 
    947                        &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) *  zsc_ws_1(ji,jj) 
    948                END DO 
    949             ENDIF               ! endif for check on lconv 
    950  
    951           END DO  ! end of ji loop 
    952        END DO     ! end of jj loop 
     903            DO jk = 2, ibld(ji,jj) 
     904               zznd_d=gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     905               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 
     906                    &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 
     907               ! 
     908               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) & 
     909                    &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) *  zsc_ws_1(ji,jj) 
     910            END DO 
     911         ENDIF               ! endif for check on lconv 
     912 
     913       END_2D 
    953914 
    954915 
    … …  
    963924       ENDWHERE 
    964925 
    965        DO jj = 2, jpjm1 
    966           DO ji = 2, jpim1 
    967              IF ( lconv(ji,jj) ) THEN 
    968                 DO jk = 2, imld(ji,jj) 
    969                    zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
    970                    ghamu(ji,jj,jk) = ghamu(ji,jj,jk) +      ( -0.05 * EXP ( -0.4 * zznd_d )   * zsc_uw_1(ji,jj)   & 
    971                         &          +                        0.00125 * EXP (      - zznd_d )   * zsc_uw_2(ji,jj) ) & 
    972                         &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) ) 
     926       DO_2D_00_00 
     927          IF ( lconv(ji,jj) ) THEN 
     928             DO jk = 2, imld(ji,jj) 
     929                zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     930                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) +      ( -0.05 * EXP ( -0.4 * zznd_d )   * zsc_uw_1(ji,jj)   & 
     931                     &          +                        0.00125 * EXP (      - zznd_d )   * zsc_uw_2(ji,jj) ) & 
     932                     &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) ) 
    973933! 
    974                    ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 *  0.15 * EXP (      - zznd_d )                       & 
    975                         &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj) 
    976                 END DO   ! end jk loop 
    977              ELSE 
     934                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 *  0.15 * EXP (      - zznd_d )                       & 
     935                     &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj) 
     936             END DO   ! end jk loop 
     937          ELSE 
    978938! Stable conditions 
    979                 DO jk = 2, ibld(ji,jj) ! corrected to ibld 
    980                    zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
    981                    ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 *   1.3 * EXP ( -0.5 * zznd_d )                       & 
    982                         &                                   * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj) 
    983                    ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp 
    984                 END DO   ! end jk loop 
    985              ENDIF 
    986           END DO  ! ji loop 
    987        END DO     ! jj loo 
     939             DO jk = 2, ibld(ji,jj) ! corrected to ibld 
     940                zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     941                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 *   1.3 * EXP ( -0.5 * zznd_d )                       & 
     942                     &                                   * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj) 
     943                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp 
     944             END DO   ! end jk loop 
     945          ENDIF 
     946       END_2D 
    988947 
    989948! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)] 
    … …  
    997956       ENDWHERE 
    998957 
    999        DO jj = 2, jpjm1 
    1000           DO ji = 2, jpim1 
    1001              IF (lconv(ji,jj) ) THEN 
    1002                 DO jk = 2, imld(ji,jj) 
    1003                    zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 
    1004                    ! calculate turbulent length scale 
    1005                    zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           & 
    1006                         &     * ( 1.0 - EXP ( -15.0 * (     1.1 - zznd_ml          ) ) ) 
    1007                    zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           & 
    1008                         &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 
    1009                    zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 
    1010                    ! non-gradient buoyancy terms 
    1011                    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 
    1012                    ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 *  zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 
    1013                 END DO 
    1014              ELSE 
    1015                 DO jk = 2, ibld(ji,jj) 
    1016                    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) 
    1017                    ghams(ji,jj,jk) = ghams(ji,jj,jk) +  zsc_ws_1(ji,jj) 
    1018                 END DO 
    1019              ENDIF 
    1020           END DO   ! ji loop 
    1021        END DO      ! jj loop 
     958       DO_2D_00_00 
     959          IF (lconv(ji,jj) ) THEN 
     960             DO jk = 2, imld(ji,jj) 
     961                zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 
     962                ! calculate turbulent length scale 
     963                zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           & 
     964                     &     * ( 1.0 - EXP ( -15.0 * (     1.1 - zznd_ml          ) ) ) 
     965                zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           & 
     966                     &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) ) 
     967                zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0) 
     968                ! non-gradient buoyancy terms 
     969                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 
     970                ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 *  zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml ) 
     971             END DO 
     972          ELSE 
     973             DO jk = 2, ibld(ji,jj) 
     974                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj) 
     975                ghams(ji,jj,jk) = ghams(ji,jj,jk) +  zsc_ws_1(ji,jj) 
     976             END DO 
     977          ENDIF 
     978       END_2D 
    1022979 
    1023980 
    … …  
    1031988       ENDWHERE 
    1032989 
    1033        DO jj = 2, jpjm1 
    1034           DO ji = 2, jpim1 
    1035              IF ( lconv(ji,jj) ) THEN 
    1036                 DO jk = 2 , imld(ji,jj) 
    1037                    zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
    1038                    ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) +   0.125 * EXP( -0.5 * zznd_d )     & 
    1039                         &                                                            * (   1.0 - EXP( -0.5 * zznd_d ) )   & 
    1040                         &                                          * zsc_uw_2(ji,jj)                                    ) 
    1041                    ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 
    1042                 END DO  ! jk loop 
    1043              ELSE 
    1044              ! stable conditions 
    1045                 DO jk = 2, ibld(ji,jj) 
    1046                    ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) 
    1047                    ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 
    1048                 END DO 
    1049              ENDIF 
    1050           END DO        ! ji loop 
    1051        END DO           ! jj loop 
     990       DO_2D_00_00 
     991          IF ( lconv(ji,jj) ) THEN 
     992             DO jk = 2 , imld(ji,jj) 
     993                zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     994                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) +   0.125 * EXP( -0.5 * zznd_d )     & 
     995                     &                                                            * (   1.0 - EXP( -0.5 * zznd_d ) )   & 
     996                     &                                          * zsc_uw_2(ji,jj)                                    ) 
     997                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 
     998             END DO  ! jk loop 
     999          ELSE 
     1000          ! stable conditions 
     1001             DO jk = 2, ibld(ji,jj) 
     1002                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) 
     1003                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj) 
     1004             END DO 
     1005          ENDIF 
     1006       END_2D 
    10521007 
    10531008! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ] 
    … …  
    10611016       ENDWHERE 
    10621017 
    1063        DO jj = 2, jpjm1 
    1064           DO ji = 2, jpim1 
    1065             IF ( lconv(ji,jj) ) THEN 
    1066                DO jk = 2, imld(ji,jj) 
    1067                   zznd_ml=gdepw_n(ji,jj,jk) / zhml(ji,jj) 
    1068                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj)                & 
    1069                        &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      & 
    1070                        &                               - EXP(     - 6.0 * zznd_ml    ) ) )  & 
    1071                        &          * ( 1.0 - EXP( - 15.0 * (         1.0 - zznd_ml    ) ) ) 
    1072                   ! 
    1073                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj)  & 
    1074                        &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      & 
    1075                        &                               - EXP(     - 6.0 * zznd_ml    ) ) )  & 
    1076                        &          * ( 1.0 - EXP ( -15.0 * (         1.0 - zznd_ml    ) ) ) 
    1077                END DO 
    1078             ELSE 
    1079                DO jk = 2, ibld(ji,jj) 
    1080                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
    1081                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    1082                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
    1083                &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 
    1084                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
    1085                &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 
    1086                END DO 
    1087             ENDIF 
    1088           ENDDO    ! ji loop 
    1089        END DO      ! jj loop 
     1018       DO_2D_00_00 
     1019         IF ( lconv(ji,jj) ) THEN 
     1020            DO jk = 2, imld(ji,jj) 
     1021               zznd_ml=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 
     1022               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj)                & 
     1023                    &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      & 
     1024                    &                               - EXP(     - 6.0 * zznd_ml    ) ) )  & 
     1025                    &          * ( 1.0 - EXP( - 15.0 * (         1.0 - zznd_ml    ) ) ) 
     1026               ! 
     1027               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj)  & 
     1028                    &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      & 
     1029                    &                               - EXP(     - 6.0 * zznd_ml    ) ) )  & 
     1030                    &          * ( 1.0 - EXP ( -15.0 * (         1.0 - zznd_ml    ) ) ) 
     1031            END DO 
     1032         ELSE 
     1033            DO jk = 2, ibld(ji,jj) 
     1034               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     1035               znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
     1036               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
     1037            &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 
     1038               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
     1039            &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj) 
     1040            END DO 
     1041         ENDIF 
     1042       END_2D 
    10901043 
    10911044 
    … …  
    11001053       ENDWHERE 
    11011054 
    1102        DO jj = 2, jpjm1 
    1103           DO ji = 2, jpim1 
    1104              IF ( lconv(ji,jj) ) THEN 
    1105                DO jk = 2, imld(ji,jj) 
    1106                   zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 
    1107                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
     1055       DO_2D_00_00 
     1056          IF ( lconv(ji,jj) ) THEN 
     1057            DO jk = 2, imld(ji,jj) 
     1058               zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 
     1059               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     1060               ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 
     1061                    & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj) 
     1062               ! 
     1063               ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 
     1064                    & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj) 
     1065            END DO 
     1066          ELSE 
     1067            DO jk = 2, ibld(ji,jj) 
     1068               znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
     1069               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     1070               IF ( zznd_d <= 2.0 ) THEN 
     1071                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 & 
     1072                       &*  ( 2.25 - 3.0  * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj) 
     1073                  ! 
     1074               ELSE 
    11081075                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 
    1109                        & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj) 
     1076                       & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj) 
    11101077                  ! 
    1111                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 
    1112                        & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj) 
    1113                END DO 
    1114              ELSE 
    1115                DO jk = 2, ibld(ji,jj) 
    1116                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    1117                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
    1118                   IF ( zznd_d <= 2.0 ) THEN 
    1119                      ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 & 
    1120                           &*  ( 2.25 - 3.0  * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj) 
    1121                      ! 
    1122                   ELSE 
    1123                      ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 
    1124                           & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj) 
    1125                      ! 
    1126                   ENDIF 
    1127  
    1128                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 
    1129                        & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj) 
    1130                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 
    1131                        & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj) 
    1132                END DO 
    1133              ENDIF 
    1134           END DO 
    1135        END DO 
     1078               ENDIF 
     1079 
     1080               ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 
     1081                    & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj) 
     1082               ghamv(ji,jj,jk) = ghamv(ji,jj,jk)& 
     1083                    & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj) 
     1084            END DO 
     1085          ENDIF 
     1086       END_2D 
    11361087! 
    11371088! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 
    11381089 
    1139       DO jj = 2, jpjm1 
    1140          DO ji = 2, jpim1 
    1141             IF ( lconv(ji,jj) ) THEN 
    1142                DO jk = 2, ibld(ji,jj) 
    1143                   znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 
    1144                   IF ( znd >= 0.0 ) THEN 
    1145                      ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 
    1146                      ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 
    1147                   ELSE 
    1148                      ghamu(ji,jj,jk) = 0._wp 
    1149                      ghamv(ji,jj,jk) = 0._wp 
    1150                   ENDIF 
    1151                END DO 
    1152             ELSE 
    1153                DO jk = 2, ibld(ji,jj) 
    1154                   znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 
    1155                   IF ( znd >= 0.0 ) THEN 
    1156                      ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 
    1157                      ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 
    1158                   ELSE 
    1159                      ghamu(ji,jj,jk) = 0._wp 
    1160                      ghamv(ji,jj,jk) = 0._wp 
    1161                   ENDIF 
    1162                END DO 
    1163             ENDIF 
    1164          END DO 
    1165       END DO 
     1090      DO_2D_00_00 
     1091         IF ( lconv(ji,jj) ) THEN 
     1092            DO jk = 2, ibld(ji,jj) 
     1093               znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 
     1094               IF ( znd >= 0.0 ) THEN 
     1095                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 
     1096                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 
     1097               ELSE 
     1098                  ghamu(ji,jj,jk) = 0._wp 
     1099                  ghamv(ji,jj,jk) = 0._wp 
     1100               ENDIF 
     1101            END DO 
     1102         ELSE 
     1103            DO jk = 2, ibld(ji,jj) 
     1104               znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 
     1105               IF ( znd >= 0.0 ) THEN 
     1106                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 
     1107                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 
     1108               ELSE 
     1109                  ghamu(ji,jj,jk) = 0._wp 
     1110                  ghamv(ji,jj,jk) = 0._wp 
     1111               ENDIF 
     1112            END DO 
     1113         ENDIF 
     1114      END_2D 
    11661115 
    11671116      ! pynocline contributions 
    11681117       ! Temporary fix to avoid instabilities when zdb_bl becomes very very small 
    11691118       zsc_uw_1 = 0._wp ! 50.0 * zla**(8.0/3.0) * zustar**2 * zhbl / ( zdb_bl + epsln ) 
    1170        DO jj = 2, jpjm1 
    1171           DO ji = 2, jpim1 
    1172              DO jk= 2, ibld(ji,jj) 
    1173                 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    1174                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 
    1175                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
    1176                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 
    1177                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 
    1178                 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 
     1119       DO_2D_00_00 
     1120          DO jk= 2, ibld(ji,jj) 
     1121             znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
     1122             ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 
     1123             ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
     1124             ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk) 
     1125             ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk) 
     1126             ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk) 
     1127          END DO 
     1128       END_2D 
     1129 
     1130! Entrainment contribution. 
     1131 
     1132       DO_2D_00_00 
     1133          IF ( lconv(ji,jj) ) THEN 
     1134            DO jk = 1, imld(ji,jj) - 1 
     1135               znd=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 
     1136               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 
     1137               ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 
     1138               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 
     1139               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd 
     1140            END DO 
     1141            DO jk = imld(ji,jj), ibld(ji,jj) 
     1142               znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
     1143               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 
     1144               ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 
     1145               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 
     1146               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 
    11791147             END DO 
    1180            END DO 
    1181        END DO 
    1182  
    1183 ! Entrainment contribution. 
    1184  
    1185        DO jj=2, jpjm1 
    1186           DO ji = 2, jpim1 
    1187              IF ( lconv(ji,jj) ) THEN 
    1188                DO jk = 1, imld(ji,jj) - 1 
    1189                   znd=gdepw_n(ji,jj,jk) / zhml(ji,jj) 
    1190                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd 
    1191                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd 
    1192                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd 
    1193                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd 
    1194                END DO 
    1195                DO jk = imld(ji,jj), ibld(ji,jj) 
    1196                   znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
    1197                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 
    1198                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 
    1199                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd ) 
    1200                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd ) 
    1201                 END DO 
    1202              ENDIF 
    1203              ghamt(ji,jj,ibld(ji,jj)) = 0._wp 
    1204              ghams(ji,jj,ibld(ji,jj)) = 0._wp 
    1205              ghamu(ji,jj,ibld(ji,jj)) = 0._wp 
    1206              ghamv(ji,jj,ibld(ji,jj)) = 0._wp 
    1207           END DO       ! ji loop 
    1208        END DO          ! jj loop 
     1148          ENDIF 
     1149          ghamt(ji,jj,ibld(ji,jj)) = 0._wp 
     1150          ghams(ji,jj,ibld(ji,jj)) = 0._wp 
     1151          ghamu(ji,jj,ibld(ji,jj)) = 0._wp 
     1152          ghamv(ji,jj,ibld(ji,jj)) = 0._wp 
     1153       END_2D 
    12091154 
    12101155 
    … …  
    12201165       ! rotate non-gradient velocity terms back to model reference frame 
    12211166 
    1222        DO jj = 2, jpjm1 
    1223           DO ji = 2, jpim1 
    1224              DO jk = 2, ibld(ji,jj) 
    1225                 ztemp = ghamu(ji,jj,jk) 
    1226                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj) 
    1227                 ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj) 
    1228              END DO 
     1167       DO_2D_00_00 
     1168          DO jk = 2, ibld(ji,jj) 
     1169             ztemp = ghamu(ji,jj,jk) 
     1170             ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj) 
     1171             ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj) 
    12291172          END DO 
    1230        END DO 
     1173       END_2D 
    12311174 
    12321175       IF(ln_dia_osm) THEN 
    … …  
    12361179! KPP-style Ri# mixing 
    12371180       IF( ln_kpprimix) THEN 
    1238           DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form) 
    1239              DO jj = 1, jpjm1 
    1240                 DO ji = 1, jpim1   ! vector opt. 
    1241                    z3du(ji,jj,jk) = 0.5 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   & 
    1242                         &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) * wumask(ji,jj,jk) & 
    1243                         &                 / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 
    1244                    z3dv(ji,jj,jk) = 0.5 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   & 
    1245                         &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) * wvmask(ji,jj,jk) & 
    1246                         &                 / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 
    1247                 END DO 
     1181          DO_3D_10_10( 2, jpkm1 ) 
     1182             z3du(ji,jj,jk) = 0.5 * (  uu(ji,jj,jk-1,Kmm) -  uu(ji  ,jj,jk,Kmm) )   & 
     1183                  &                 * (  uu(ji,jj,jk-1,Kbb) -  uu(ji  ,jj,jk,Kbb) ) * wumask(ji,jj,jk) & 
     1184                  &                 / (  e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) 
     1185             z3dv(ji,jj,jk) = 0.5 * (  vv(ji,jj,jk-1,Kmm) -  vv(ji,jj  ,jk,Kmm) )   & 
     1186                  &                 * (  vv(ji,jj,jk-1,Kbb) -  vv(ji,jj  ,jk,Kbb) ) * wvmask(ji,jj,jk) & 
     1187                  &                 / (  e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) 
     1188          END_3D 
     1189      ! 
     1190         DO_3D_00_00( 2, jpkm1 ) 
     1191            !                                          ! shear prod. at w-point weightened by mask 
     1192            zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
     1193               &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 
     1194            !                                          ! local Richardson number 
     1195            zri   = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln) 
     1196            zfri =  MIN( zri / rn_riinfty , 1.0_wp ) 
     1197            zfri  = ( 1.0_wp - zfri * zfri ) 
     1198            zrimix(ji,jj,jk)  =  zfri * zfri  * zfri * wmask(ji, jj, jk) 
     1199         END_3D 
     1200 
     1201          DO_2D_00_00 
     1202             DO jk = ibld(ji,jj) + 1, jpkm1 
     1203                zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 
     1204                zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 
    12481205             END DO 
    1249           END DO 
    1250       ! 
    1251          DO jk = 2, jpkm1 
    1252             DO jj = 2, jpjm1 
    1253                DO ji = 2, jpim1   ! vector opt. 
    1254                   !                                          ! shear prod. at w-point weightened by mask 
    1255                   zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
    1256                      &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 
    1257                   !                                          ! local Richardson number 
    1258                   zri   = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln) 
    1259                   zfri =  MIN( zri / rn_riinfty , 1.0_wp ) 
    1260                   zfri  = ( 1.0_wp - zfri * zfri ) 
    1261                   zrimix(ji,jj,jk)  =  zfri * zfri  * zfri * wmask(ji, jj, jk) 
    1262                 END DO 
    1263              END DO 
    1264           END DO 
    1265  
    1266           DO jj = 2, jpjm1 
    1267              DO ji = 2, jpim1 
    1268                 DO jk = ibld(ji,jj) + 1, jpkm1 
    1269                    zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 
    1270                    zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri 
    1271                 END DO 
    1272              END DO 
    1273           END DO 
     1206          END_2D 
    12741207 
    12751208       END IF ! ln_kpprimix = .true. 
    … …  
    12771210! KPP-style set diffusivity large if unstable below BL 
    12781211       IF( ln_convmix) THEN 
    1279           DO jj = 2, jpjm1 
    1280              DO ji = 2, jpim1 
    1281                 DO jk = ibld(ji,jj) + 1, jpkm1 
    1282                   IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv 
    1283                 END DO 
     1212          DO_2D_00_00 
     1213             DO jk = ibld(ji,jj) + 1, jpkm1 
     1214               IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv 
    12841215             END DO 
    1285           END DO 
     1216          END_2D 
    12861217       END IF ! ln_convmix = .true. 
    12871218 
    … …  
    12911222       ! GN 25/8: need to change tmask --> wmask 
    12921223 
    1293      DO jk = 2, jpkm1 
    1294          DO jj = 2, jpjm1 
    1295              DO ji = 2, jpim1 
    1296                 p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk) 
    1297                 p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk) 
    1298              END DO 
    1299          END DO 
    1300      END DO 
     1224     DO_3D_00_00( 2, jpkm1 ) 
     1225          p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk) 
     1226          p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk) 
     1227     END_3D 
    13011228      ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid  (sign unchanged), needed to caclulate gham[uv] on u and v grids 
    13021229     CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1. , p_avm, 'W', 1.,   & 
    13031230      &                  ghamu, 'W', 1. , ghamv, 'W', 1. ) 
    1304        DO jk = 2, jpkm1 
    1305            DO jj = 2, jpjm1 
    1306                DO ji = 2, jpim1 
    1307                   ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) & 
    1308                      &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk) 
    1309  
    1310                   ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) & 
    1311                       &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk) 
    1312  
    1313                   ghamt(ji,jj,jk) =  ghamt(ji,jj,jk) * tmask(ji,jj,jk) 
    1314                   ghams(ji,jj,jk) =  ghams(ji,jj,jk) * tmask(ji,jj,jk) 
    1315                END DO 
    1316            END DO 
    1317         END DO 
     1231       DO_3D_00_00( 2, jpkm1 ) 
     1232            ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) & 
     1233               &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk) 
     1234 
     1235            ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) & 
     1236                &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk) 
     1237 
     1238            ghamt(ji,jj,jk) =  ghamt(ji,jj,jk) * tmask(ji,jj,jk) 
     1239            ghams(ji,jj,jk) =  ghams(ji,jj,jk) * tmask(ji,jj,jk) 
     1240       END_3D 
    13181241        ! Lateral boundary conditions on final outputs for gham[ts],  on W-grid  (sign unchanged) 
    13191242        ! Lateral boundary conditions on final outputs for gham[uv],  on [UV]-grid  (sign unchanged) 
    … …  
    13641287 
    13651288 
    1366    SUBROUTINE zdf_osm_init 
     1289   SUBROUTINE zdf_osm_init( Kmm )  
    13671290     !!---------------------------------------------------------------------- 
    13681291     !!                  ***  ROUTINE zdf_osm_init  *** 
    … …  
    13761299     !! ** input   :   Namlist namosm 
    13771300     !!---------------------------------------------------------------------- 
     1301     INTEGER, INTENT(in)    :: Kmm ! time level index (middle) 
     1302     ! 
    13781303     INTEGER  ::   ios            ! local integer 
    13791304     INTEGER  ::   ji, jj, jk     ! dummy loop indices 
    … …  
    13841309     !!---------------------------------------------------------------------- 
    13851310     ! 
    1386      REWIND( numnam_ref )              ! Namelist namzdf_osm in reference namelist : Osmosis ML model 
    13871311     READ  ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901) 
    13881312901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' ) 
    13891313 
    1390      REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy 
    13911314     READ  ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 ) 
    13921315902  IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' ) 
    … …  
    14231346     IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 
    14241347 
    1425      call osm_rst( nit000, 'READ' ) !* read or initialize hbl 
     1348     call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl 
    14261349 
    14271350     IF( ln_zdfddm) THEN 
    … …  
    14591382        etmean(:,:,:) = 0.e0 
    14601383 
    1461         DO jk = 1, jpkm1 
    1462            DO jj = 2, jpjm1 
    1463               DO ji = 2, jpim1   ! vector opt. 
    1464                  etmean(ji,jj,jk) = tmask(ji,jj,jk)                     & 
    1465                       &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   & 
    1466                       &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  ) 
    1467               END DO 
    1468            END DO 
    1469         END DO 
     1384        DO_3D_00_00( 1, jpkm1 ) 
     1385           etmean(ji,jj,jk) = tmask(ji,jj,jk)                     & 
     1386                &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   & 
     1387                &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  ) 
     1388        END_3D 
    14701389 
    14711390     CASE ( 1 )                ! horizontal average 
    … …  
    14771396        etmean(:,:,:) = 0.e0 
    14781397 
    1479         DO jk = 1, jpkm1 
    1480            DO jj = 2, jpjm1 
    1481               DO ji = 2, jpim1   ! vector opt. 
    1482                  etmean(ji,jj,jk) = tmask(ji, jj,jk)                           & 
    1483                       & / MAX( 1., 2.* tmask(ji,jj,jk)                           & 
    1484                       &      +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk)   & 
    1485                       &             +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) & 
    1486                       &      +1. * ( tmask(ji-1,jj  ,jk) + tmask(ji  ,jj+1,jk)   & 
    1487                       &             +tmask(ji  ,jj-1,jk) + tmask(ji+1,jj  ,jk) ) ) 
    1488               END DO 
    1489            END DO 
    1490         END DO 
     1398        DO_3D_00_00( 1, jpkm1 ) 
     1399           etmean(ji,jj,jk) = tmask(ji, jj,jk)                           & 
     1400                & / MAX( 1., 2.* tmask(ji,jj,jk)                           & 
     1401                &      +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk)   & 
     1402                &             +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) & 
     1403                &      +1. * ( tmask(ji-1,jj  ,jk) + tmask(ji  ,jj+1,jk)   & 
     1404                &             +tmask(ji  ,jj-1,jk) + tmask(ji+1,jj  ,jk) ) ) 
     1405        END_3D 
    14911406 
    14921407     CASE DEFAULT 
    … …  
    15171432 
    15181433 
    1519    SUBROUTINE osm_rst( kt, cdrw ) 
     1434   SUBROUTINE osm_rst( kt, Kmm, cdrw ) 
    15201435     !!--------------------------------------------------------------------- 
    15211436     !!                   ***  ROUTINE osm_rst  *** 
    … …  
    15271442     !!---------------------------------------------------------------------- 
    15281443 
    1529      INTEGER, INTENT(in) :: kt 
     1444     INTEGER         , INTENT(in) ::   kt     ! ocean time step index 
     1445     INTEGER         , INTENT(in) ::   Kmm    ! ocean time level index (middle) 
    15301446     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    15311447 
    … …  
    15451461        id1 = iom_varid( numror, 'wn'   , ldstop = .FALSE. ) 
    15461462        IF( id1 > 0 ) THEN                       ! 'wn' exists; read 
    1547            CALL iom_get( numror, jpdom_autoglo, 'wn', wn, ldxios = lrxios ) 
    1548            WRITE(numout,*) ' ===>>>> :  wn read from restart file' 
     1463           CALL iom_get( numror, jpdom_autoglo, 'wn', ww, ldxios = lrxios ) 
     1464           WRITE(numout,*) ' ===>>>> :  ww read from restart file' 
    15491465        ELSE 
    1550            wn(:,:,:) = 0._wp 
    1551            WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially' 
     1466           ww(:,:,:) = 0._wp 
     1467           WRITE(numout,*) ' ===>>>> :  ww not in restart file, set to zero initially' 
    15521468        END IF 
    15531469        id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. ) 
    … …  
    15681484     IF( TRIM(cdrw) == 'WRITE') THEN     !* Write hbli into the restart file, then return 
    15691485        IF(lwp) WRITE(numout,*) '---- osm-rst ----' 
    1570          CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn  , ldxios = lwxios ) 
     1486         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , ww  , ldxios = lwxios ) 
    15711487         CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl , ldxios = lwxios ) 
    15721488         CALL iom_rstput( kt, nitrst, numrow, 'hbli'   , hbli, ldxios = lwxios ) 
    … …  
    15801496     ALLOCATE( imld_rst(jpi,jpj) ) 
    15811497     ! w-level of the mixing and mixed layers 
    1582      CALL eos_rab( tsn, rab_n ) 
    1583      CALL bn2(tsn, rab_n, rn2) 
     1498     CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm ) 
     1499     CALL bn2(ts(:,:,:,:,Kmm), rab_n, rn2, Kmm) 
    15841500     imld_rst(:,:)  = nlb10         ! Initialization to the number of w ocean point 
    15851501     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2 
    … …  
    15871503     ! 
    15881504     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2 
    1589      DO jk = 1, jpkm1 
    1590         DO jj = 1, jpj              ! Mixed layer level: w-level 
    1591            DO ji = 1, jpi 
    1592               ikt = mbkt(ji,jj) 
    1593               hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk) 
    1594               IF( hbl(ji,jj) < zN2_c )   imld_rst(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level 
    1595            END DO 
    1596         END DO 
    1597      END DO 
     1505     DO_3D_11_11( 1, jpkm1 ) 
     1506        ikt = mbkt(ji,jj) 
     1507        hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) 
     1508        IF( hbl(ji,jj) < zN2_c )   imld_rst(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level 
     1509     END_3D 
    15981510     ! 
    1599      DO jj = 1, jpj 
    1600         DO ji = 1, jpi 
    1601            iiki = imld_rst(ji,jj) 
    1602            hbl (ji,jj) = gdepw_n(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth 
    1603         END DO 
    1604      END DO 
     1511     DO_2D_11_11 
     1512        iiki = imld_rst(ji,jj) 
     1513        hbl (ji,jj) = gdepw(ji,jj,iiki  ,Kmm) * ssmask(ji,jj)    ! Turbocline depth 
     1514     END_2D 
    16051515     hbl = MAX(hbl,epsln) 
    16061516     hbli(:,:) = hbl(:,:) 
    … …  
    16101520 
    16111521 
    1612    SUBROUTINE tra_osm( kt ) 
     1522   SUBROUTINE tra_osm( kt, Kmm, pts, Krhs ) 
    16131523      !!---------------------------------------------------------------------- 
    16141524      !!                  ***  ROUTINE tra_osm  *** 
    … …  
    16201530      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace 
    16211531      !!---------------------------------------------------------------------- 
    1622       INTEGER, INTENT(in) :: kt 
     1532      INTEGER                                  , INTENT(in)    :: kt        ! time step index 
     1533      INTEGER                                  , INTENT(in)    :: Kmm, Krhs ! time level indices 
     1534      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts       ! active tracers and RHS of tracer equation 
     1535      ! 
    16231536      INTEGER :: ji, jj, jk 
    16241537      ! 
    … …  
    16301543 
    16311544      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
    1632          ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    1633          ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
     1545         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 
     1546         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 
    16341547      ENDIF 
    16351548 
    16361549      ! add non-local temperature and salinity flux 
    1637       DO jk = 1, jpkm1 
    1638          DO jj = 2, jpjm1 
    1639             DO ji = 2, jpim1 
    1640                tsa(ji,jj,jk,jp_tem) =  tsa(ji,jj,jk,jp_tem)                      & 
    1641                   &                 - (  ghamt(ji,jj,jk  )  & 
    1642                   &                    - ghamt(ji,jj,jk+1) ) /e3t_n(ji,jj,jk) 
    1643                tsa(ji,jj,jk,jp_sal) =  tsa(ji,jj,jk,jp_sal)                      & 
    1644                   &                 - (  ghams(ji,jj,jk  )  & 
    1645                   &                    - ghams(ji,jj,jk+1) ) / e3t_n(ji,jj,jk) 
    1646             END DO 
    1647          END DO 
    1648       END DO 
     1550      DO_3D_00_00( 1, jpkm1 ) 
     1551         pts(ji,jj,jk,jp_tem,Krhs) =  pts(ji,jj,jk,jp_tem,Krhs)                      & 
     1552            &                 - (  ghamt(ji,jj,jk  )  & 
     1553            &                    - ghamt(ji,jj,jk+1) ) /e3t(ji,jj,jk,Kmm) 
     1554         pts(ji,jj,jk,jp_sal,Krhs) =  pts(ji,jj,jk,jp_sal,Krhs)                      & 
     1555            &                 - (  ghams(ji,jj,jk  )  & 
     1556            &                    - ghams(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm) 
     1557      END_3D 
    16491558 
    16501559 
    16511560      ! save the non-local tracer flux trends for diagnostic 
    16521561      IF( l_trdtra )   THEN 
    1653          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    1654          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
     1562         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 
     1563         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 
    16551564!!bug gm jpttdzdf ==> jpttosm 
    1656          CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
    1657          CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 
     1565         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
     1566         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdf, ztrds ) 
    16581567         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds ) 
    16591568      ENDIF 
    16601569 
    1661       IF(ln_ctl) THEN 
    1662          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' osm  - Ta: ', mask1=tmask,   & 
    1663          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     1570      IF(sn_cfctl%l_prtctl) THEN 
     1571         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' osm  - Ta: ', mask1=tmask,   & 
     1572         &             tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    16641573      ENDIF 
    16651574      ! 
    … …  
    16841593 
    16851594 
    1686    SUBROUTINE dyn_osm( kt ) 
     1595   SUBROUTINE dyn_osm( kt, Kmm, puu, pvv, Krhs ) 
    16871596      !!---------------------------------------------------------------------- 
    16881597      !!                  ***  ROUTINE dyn_osm  *** 
    … …  
    16931602      !! ** Method  :   ??? 
    16941603      !!---------------------------------------------------------------------- 
    1695       INTEGER, INTENT(in) ::   kt   ! 
     1604      INTEGER                             , INTENT( in )  ::  kt          ! ocean time step index 
     1605      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     1606      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    16961607      ! 
    16971608      INTEGER :: ji, jj, jk   ! dummy loop indices 
    … …  
    17051616      !code saving tracer trends removed, replace with trdmxl_oce 
    17061617 
    1707       DO jk = 1, jpkm1           ! add non-local u and v fluxes 
    1708          DO jj = 2, jpjm1 
    1709             DO ji = 2, jpim1 
    1710                ua(ji,jj,jk) =  ua(ji,jj,jk)                      & 
    1711                   &                 - (  ghamu(ji,jj,jk  )  & 
    1712                   &                    - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk) 
    1713                va(ji,jj,jk) =  va(ji,jj,jk)                      & 
    1714                   &                 - (  ghamv(ji,jj,jk  )  & 
    1715                   &                    - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk) 
    1716             END DO 
    1717          END DO 
    1718       END DO 
     1618      DO_3D_00_00( 1, jpkm1 ) 
     1619         puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs)                      & 
     1620            &                 - (  ghamu(ji,jj,jk  )  & 
     1621            &                    - ghamu(ji,jj,jk+1) ) / e3u(ji,jj,jk,Kmm) 
     1622         pvv(ji,jj,jk,Krhs) =  pvv(ji,jj,jk,Krhs)                      & 
     1623            &                 - (  ghamv(ji,jj,jk  )  & 
     1624            &                    - ghamv(ji,jj,jk+1) ) / e3v(ji,jj,jk,Kmm) 
     1625      END_3D 
    17191626      ! 
    17201627      ! code for saving tracer trends removed 
Note: See TracChangeset for help on using the changeset viewer.