Ignore:
Timestamp:
2021-05-19T17:20:25+02:00 (7 months ago)
Author:
smueller
Message:

Enabling of the extended halo option (nn_hls=2) in the implementation of the OSMOSIS boundary layer scheme (ticket #2353)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90

    r14868 r14889  
    262262      zdf_osm_alloc = zdf_osm_alloc + ierr 
    263263      ! 
    264       ALLOCATE( etmean(A2D(0),jpk), dh(jpi,jpj), r1_ft(A2D(0)), STAT=ierr ) 
     264      ALLOCATE( etmean(A2D((nn_hls-1)),jpk), dh(jpi,jpj), r1_ft(A2D((nn_hls-1))), STAT=ierr ) 
    265265      zdf_osm_alloc = zdf_osm_alloc + ierr 
    266266      ! 
    267       ALLOCATE( nbld(jpi,jpj), nmld(A2D(0)), STAT=ierr ) 
     267      ALLOCATE( nbld(jpi,jpj), nmld(A2D((nn_hls-1))), STAT=ierr ) 
    268268      zdf_osm_alloc = zdf_osm_alloc + ierr 
    269269      ! 
    270       ALLOCATE( n_ddh(A2D(0)), STAT=ierr ) 
     270      ALLOCATE( n_ddh(A2D((nn_hls-1))), STAT=ierr ) 
    271271      zdf_osm_alloc = zdf_osm_alloc + ierr 
    272272      ! 
    273       ALLOCATE( l_conv(A2D(0)), l_shear(A2D(0)), l_coup(A2D(0)), l_pyc(A2D(0)), l_flux(A2D(0)), l_mle(A2D(0)), STAT=ierr ) 
     273      ALLOCATE( l_conv(A2D((nn_hls-1))), l_shear(A2D((nn_hls-1))), l_coup(A2D((nn_hls-1))), l_pyc(A2D((nn_hls-1))), l_flux(A2D((nn_hls-1))), l_mle(A2D((nn_hls-1))), STAT=ierr ) 
    274274      zdf_osm_alloc = zdf_osm_alloc + ierr 
    275275      ! 
    276       ALLOCATE( swth0(A2D(0)),     sws0(A2D(0)),   swb0(A2D(0)),  suw0(A2D(0)),  sustar(A2D(0)), scos_wind(A2D(0)),   & 
    277          &      ssin_wind(A2D(0)), swthav(A2D(0)), swsav(A2D(0)), swbav(A2D(0)), sustke(A2D(0)), dstokes(A2D(0)),     & 
    278          &      swstrl(A2D(0)),    swstrc(A2D(0)), sla(A2D(0)),   svstr(A2D(0)), shol(A2D(0)),   STAT=ierr ) 
     276      ALLOCATE( swth0(A2D((nn_hls-1))),     sws0(A2D((nn_hls-1))),   swb0(A2D((nn_hls-1))),  suw0(A2D((nn_hls-1))),  sustar(A2D((nn_hls-1))), scos_wind(A2D((nn_hls-1))),   & 
     277         &      ssin_wind(A2D((nn_hls-1))), swthav(A2D((nn_hls-1))), swsav(A2D((nn_hls-1))), swbav(A2D((nn_hls-1))), sustke(A2D((nn_hls-1))), dstokes(A2D((nn_hls-1))),     & 
     278         &      swstrl(A2D((nn_hls-1))),    swstrc(A2D((nn_hls-1))), sla(A2D((nn_hls-1))),   svstr(A2D((nn_hls-1))), shol(A2D((nn_hls-1))),   STAT=ierr ) 
    279279      zdf_osm_alloc = zdf_osm_alloc + ierr 
    280280      ! 
    281       ALLOCATE( av_t_bl(A2D(0)), av_s_bl(A2D(0)), av_u_bl(A2D(0)), av_v_bl(A2D(0)), av_b_bl(A2D(0)), STAT=ierr) 
     281      ALLOCATE( av_t_bl(A2D((nn_hls-1))), av_s_bl(A2D((nn_hls-1))), av_u_bl(A2D((nn_hls-1))), av_v_bl(A2D((nn_hls-1))), av_b_bl(A2D((nn_hls-1))), STAT=ierr) 
    282282      zdf_osm_alloc = zdf_osm_alloc + ierr 
    283283      ! 
    284       ALLOCATE( av_dt_bl(A2D(0)), av_ds_bl(A2D(0)), av_du_bl(A2D(0)), av_dv_bl(A2D(0)), av_db_bl(A2D(0)), STAT=ierr) 
     284      ALLOCATE( av_dt_bl(A2D((nn_hls-1))), av_ds_bl(A2D((nn_hls-1))), av_du_bl(A2D((nn_hls-1))), av_dv_bl(A2D((nn_hls-1))), av_db_bl(A2D((nn_hls-1))), STAT=ierr) 
    285285      zdf_osm_alloc = zdf_osm_alloc + ierr 
    286286      ! 
    287       ALLOCATE( av_t_ml(A2D(0)), av_s_ml(A2D(0)), av_u_ml(A2D(0)), av_v_ml(A2D(0)), av_b_ml(A2D(0)), STAT=ierr) 
     287      ALLOCATE( av_t_ml(A2D((nn_hls-1))), av_s_ml(A2D((nn_hls-1))), av_u_ml(A2D((nn_hls-1))), av_v_ml(A2D((nn_hls-1))), av_b_ml(A2D((nn_hls-1))), STAT=ierr) 
    288288      zdf_osm_alloc = zdf_osm_alloc + ierr 
    289289      ! 
    290       ALLOCATE( av_dt_ml(A2D(0)), av_ds_ml(A2D(0)), av_du_ml(A2D(0)), av_dv_ml(A2D(0)), av_db_ml(A2D(0)), STAT=ierr) 
     290      ALLOCATE( av_dt_ml(A2D((nn_hls-1))), av_ds_ml(A2D((nn_hls-1))), av_du_ml(A2D((nn_hls-1))), av_dv_ml(A2D((nn_hls-1))), av_db_ml(A2D((nn_hls-1))), STAT=ierr) 
    291291      zdf_osm_alloc = zdf_osm_alloc + ierr 
    292292      ! 
    293       ALLOCATE( av_t_mle(A2D(0)), av_s_mle(A2D(0)), av_u_mle(A2D(0)), av_v_mle(A2D(0)), av_b_mle(A2D(0)), STAT=ierr) 
     293      ALLOCATE( av_t_mle(A2D((nn_hls-1))), av_s_mle(A2D((nn_hls-1))), av_u_mle(A2D((nn_hls-1))), av_v_mle(A2D((nn_hls-1))), av_b_mle(A2D((nn_hls-1))), STAT=ierr) 
    294294      zdf_osm_alloc = zdf_osm_alloc + ierr 
    295295      ! 
     
    351351      ! 
    352352      ! Scales 
    353       REAL(wp), DIMENSION(A2D(0))  ::   zrad0       ! Surface solar temperature flux (deg m/s) 
    354       REAL(wp), DIMENSION(A2D(0))  ::   zradh       ! Radiative flux at bl base (Buoyancy units) 
     353      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zrad0       ! Surface solar temperature flux (deg m/s) 
     354      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zradh       ! Radiative flux at bl base (Buoyancy units) 
    355355      REAL(wp)                     ::   zradav      ! Radiative flux, bl average (Buoyancy Units) 
    356356      REAL(wp)                     ::   zvw0        ! Surface v-momentum flux 
    357       REAL(wp), DIMENSION(A2D(0))  ::   zwb0tot     ! Total surface buoyancy flux including insolation 
    358       REAL(wp), DIMENSION(A2D(0))  ::   zwb_ent     ! Buoyancy entrainment flux 
    359       REAL(wp), DIMENSION(A2D(0))  ::   zwb_min 
    360       REAL(wp), DIMENSION(A2D(0))  ::   zwb_fk_b    ! MLE buoyancy flux averaged over OSBL 
    361       REAL(wp), DIMENSION(A2D(0))  ::   zwb_fk      ! Max MLE buoyancy flux 
    362       REAL(wp), DIMENSION(A2D(0))  ::   zdiff_mle   ! Extra MLE vertical diff 
    363       REAL(wp), DIMENSION(A2D(0))  ::   zvel_mle    ! Velocity scale for dhdt with stable ML and FK 
     357      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zwb0tot     ! Total surface buoyancy flux including insolation 
     358      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zwb_ent     ! Buoyancy entrainment flux 
     359      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zwb_min 
     360      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zwb_fk_b    ! MLE buoyancy flux averaged over OSBL 
     361      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zwb_fk      ! Max MLE buoyancy flux 
     362      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zdiff_mle   ! Extra MLE vertical diff 
     363      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zvel_mle    ! Velocity scale for dhdt with stable ML and FK 
    364364      ! 
    365365      ! mixed-layer variables 
    366       INTEGER,  DIMENSION(A2D(0)) ::   jk_ext   ! Offset for external level 
    367       ! 
    368       REAL(wp), DIMENSION(A2D(0)) ::   zhbl   ! BL depth - grid 
    369       REAL(wp), DIMENSION(A2D(0)) ::   zhml   ! ML depth - grid 
    370       ! 
    371       REAL(wp), DIMENSION(A2D(0))  ::   zhmle   ! MLE depth - grid 
     366      INTEGER,  DIMENSION(A2D((nn_hls-1))) ::   jk_ext   ! Offset for external level 
     367      ! 
     368      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zhbl   ! BL depth - grid 
     369      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zhml   ! ML depth - grid 
     370      ! 
     371      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zhmle   ! MLE depth - grid 
    372372      REAL(wp), DIMENSION(jpi,jpj) ::   zmld   ! ML depth on grid 
    373373      ! 
    374       REAL(wp), DIMENSION(A2D(0))  ::   zdh   ! Pycnocline depth - grid 
    375       REAL(wp), DIMENSION(A2D(0))  ::   zdhdt   ! BL depth tendency 
    376       REAL(wp), DIMENSION(A2D(0))  ::   zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext   ! External temperature/salinity and buoyancy gradients 
     374      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zdh   ! Pycnocline depth - grid 
     375      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zdhdt   ! BL depth tendency 
     376      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext   ! External temperature/salinity and buoyancy gradients 
    377377      REAL(wp), DIMENSION(jpi,jpj) ::   zdtdx, zdtdy, zdsdx, zdsdy   ! Horizontal gradients for Fox-Kemper parametrization 
    378378      ! 
    379       REAL(wp), DIMENSION(A2D(0)) ::   zdbds_mle   ! Magnitude of horizontal buoyancy gradient 
     379      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zdbds_mle   ! Magnitude of horizontal buoyancy gradient 
    380380      ! Flux-gradient relationship variables 
    381       REAL(wp), DIMENSION(A2D(0)) ::   zshear   ! Shear production 
    382       ! 
    383       REAL(wp), DIMENSION(A2D(0)) ::   zhbl_t   ! Holds boundary layer depth updated by full timestep 
     381      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zshear   ! Shear production 
     382      ! 
     383      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zhbl_t   ! Holds boundary layer depth updated by full timestep 
    384384      ! 
    385385      ! For calculating Ri#-dependent mixing 
     
    392392      REAL(wp)                        ::   zz0, zz1, zfac 
    393393      REAL(wp)                        ::   zus_x, zus_y   ! Temporary Stokes drift 
    394       REAL(wp), DIMENSION(A2D(0),jpk) ::   zviscos   ! Viscosity 
    395       REAL(wp), DIMENSION(A2D(0),jpk) ::   zdiffut   ! t-diffusivity 
     394      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk) ::   zviscos   ! Viscosity 
     395      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk) ::   zdiffut   ! t-diffusivity 
    396396      REAL(wp)                        ::   zabsstke 
    397397      REAL(wp)                        ::   zsqrtpi, z_two_thirds, zthickness 
     
    423423      ! 
    424424      ghamt(:,:,:)      = pp_large ; ghams(:,:,:)      = pp_large 
    425       ghamt(A2D(0),:)   = 0.0_wp   ; ghams(A2D(0),:)   = 0.0_wp 
     425      ghamt(A2D((nn_hls-1)),:)   = 0.0_wp   ; ghams(A2D((nn_hls-1)),:)   = 0.0_wp 
    426426      ghamu(:,:,:)      = pp_large ; ghamv(:,:,:)      = pp_large 
    427       ghamu(A2D(0),:)   = 0.0_wp   ; ghamv(A2D(0),:)   = 0.0_wp 
     427      ghamu(A2D((nn_hls-1)),:)   = 0.0_wp   ; ghamv(A2D((nn_hls-1)),:)   = 0.0_wp 
    428428      ! 
    429429      zdiff_mle(:,:) = 0.0_wp 
    430430      ! 
    431       ! hbl = MAX(hbl,epsln) 
     431      ! Ensure only positive hbl values are accessed when using extended halo 
     432      ! (nn_hls==2) 
     433      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     434         hbl(ji,jj) = MAX( hbl(ji,jj), epsln ) 
     435      END_2D 
     436      ! 
    432437      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    433438      ! Calculate boundary layer scales 
     
    437442      zz0 =           rn_abs   ! Assume two-band radiation model for depth of OSBL - surface equi-partition in 2-bands 
    438443      zz1 =  1.0_wp - rn_abs 
    439       DO_2D( 0, 0, 0, 0 ) 
     444      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    440445         zrad0(ji,jj)  = qsr(ji,jj) * r1_rho0_rcp   ! Surface downward irradiance (so always +ve) 
    441446         zradh(ji,jj)  = zrad0(ji,jj) *                                &   ! Downwards irradiance at base of boundary layer 
     
    448453            &                                                 zradav )                              !    over depth of OSBL 
    449454      END_2D 
    450       DO_2D( 0, 0, 0, 0 ) 
     455      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    451456         sws0(ji,jj)    = -1.0_wp * ( ( emp(ji,jj) - rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) +   &   ! Upwards surface salinity flux 
    452457            &                         sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1)                      !    for non-local term 
     
    459464            &             grav  * zbeta * swsav(ji,jj)                      ! OBSBL 
    460465      END_2D 
    461       DO_2D( 0, 0, 0, 0 ) 
     466      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    462467         suw0(ji,jj)    = -0.5_wp * (utau(ji-1,jj) + utau(ji,jj)) * r1_rho0 * tmask(ji,jj,1)   ! Surface upward velocity fluxes 
    463468         zvw0           = -0.5_wp * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1) 
     
    471476         ! Assume constant La#=0.3 
    472477      CASE(0) 
    473          DO_2D( 0, 0, 0, 0 ) 
     478         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    474479            zus_x = scos_wind(ji,jj) * sustar(ji,jj) / 0.3_wp**2 
    475480            zus_y = ssin_wind(ji,jj) * sustar(ji,jj) / 0.3_wp**2 
     
    480485         ! Assume Pierson-Moskovitz wind-wave spectrum 
    481486      CASE(1) 
    482          DO_2D( 0, 0, 0, 0 ) 
     487         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    483488            ! Use wind speed wndm included in sbc_oce module 
    484489            sustke(ji,jj)  = MAX ( 0.016_wp * wndm(ji,jj), 1e-8_wp ) 
     
    489494         zfac =  2.0_wp * rpi / 16.0_wp 
    490495         ! 
    491          DO_2D( 0, 0, 0, 0 ) 
     496         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    492497            IF ( hsw(ji,jj) > 1e-4_wp ) THEN 
    493498               ! Use  wave fields 
     
    506511      IF (ln_zdfosm_ice_shelter) THEN 
    507512         ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice 
    508          DO_2D( 0, 0, 0, 0 ) 
     513         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    509514            sustke(ji,jj)  = sustke(ji,jj)  * ( 1.0_wp - fr_i(ji,jj) ) 
    510515            dstokes(ji,jj) = dstokes(ji,jj) * ( 1.0_wp - fr_i(ji,jj) ) 
     
    519524         ! It could represent the effects of the spread of wave directions around the mean wind. The effect of this adjustment needs to be tested. 
    520525         IF(nn_osm_wave > 0) THEN 
    521             sustke(A2D(0)) = rn_zdfosm_adjust_sd * sustke(A2D(0)) 
     526            sustke(A2D((nn_hls-1))) = rn_zdfosm_adjust_sd * sustke(A2D((nn_hls-1))) 
    522527         END IF 
    523528      CASE(1) 
     
    526531         zsqrtpi = SQRT(rpi) 
    527532         z_two_thirds = 2.0_wp / 3.0_wp 
    528          DO_2D( 0, 0, 0, 0 ) 
     533         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    529534            zthickness = rn_osm_hblfrac*hbl(ji,jj) 
    530535            z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 1e-7_wp ) 
     
    540545         ! Assumes approximate depth profile of SD from Breivik (2016) 
    541546         zsqrtpi = SQRT(rpi) 
    542          DO_2D( 0, 0, 0, 0 ) 
     547         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    543548            zthickness = rn_osm_hblfrac*hbl(ji,jj) 
    544549            z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 1e-7_wp ) 
     
    562567      ! Langmuir velocity scale (swstrl), La # (sla) 
    563568      ! Mixed scale (svstr), convective velocity scale (swstrc) 
    564       DO_2D( 0, 0, 0, 0 ) 
     569      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    565570         ! Langmuir velocity scale (swstrl), at T-point 
    566571         swstrl(ji,jj) = ( sustar(ji,jj) * sustar(ji,jj) * sustke(ji,jj) )**pthird 
     
    596601      hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) ) 
    597602      nbld(:,:) = 4 
    598       DO_3D( 1, 1, 1, 1, 5, jpkm1 ) 
     603      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 5, jpkm1 ) 
    599604         IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 
    600605            nbld(ji,jj) = MIN(mbkt(ji,jj)-2, jk) 
     
    603608      ! ########################################################################## 
    604609      ! 
    605       DO_2D( 0, 0, 0, 0 ) 
     610      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    606611         zhbl(ji,jj) = gdepw(ji,jj,nbld(ji,jj),Kmm) 
    607612         nmld(ji,jj) = MAX( 3, nbld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji,jj,nbld(ji,jj)-1,Kmm) ), 1 ) ) 
     
    612617      ! Averages over well-mixed and boundary layer, note BL averages use jk_ext=2 everywhere 
    613618      jk_ext(:,:) = 1   ! ag 19/03 
    614       CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D(0)), av_t_bl, av_s_bl,      & 
     619      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D((nn_hls-1))), av_t_bl, av_s_bl,      & 
    615620         &                           av_b_bl, av_u_bl, av_v_bl, jk_ext, av_dt_bl,   & 
    616621         &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) 
    617       jk_ext(:,:) = nbld(A2D(0)) - nmld(:,:) + jk_ext(:,:) + 1   ! ag 19/03 
     622      jk_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(:,:) + jk_ext(:,:) + 1   ! ag 19/03 
    618623      CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,          & 
    619624         &                           av_b_ml, av_u_ml, av_v_ml, jk_ext, av_dt_ml,   & 
     
    632637         ! Fox-Kemper Scheme 
    633638         mld_prof = 4 
    634          DO_3D( 0, 0, 0, 0, 5, jpkm1 ) 
     639         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 ) 
    635640            IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk) 
    636641         END_3D 
    637          CALL zdf_osm_vertical_average( Kbb, Kmm, mld_prof(A2D(0)), av_t_mle, av_s_mle,   & 
     642         CALL zdf_osm_vertical_average( Kbb, Kmm, mld_prof(A2D((nn_hls-1))), av_t_mle, av_s_mle,   & 
    638643            &                           av_b_mle, av_u_mle, av_v_mle ) 
    639644         ! 
    640          DO_2D( 0, 0, 0, 0 ) 
     645         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    641646            zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 
    642647         END_2D 
     
    656661         l_flux(:,:) = .FALSE. 
    657662         l_mle(:,:)  = .FALSE. 
    658          DO_2D( 0, 0, 0, 0 ) 
     663         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    659664            IF ( l_conv(ji,jj) .AND. av_db_bl(ji,jj) < rn_osm_bl_thresh ) l_pyc(ji,jj) = .FALSE. 
    660665         END_2D 
     
    662667      ! 
    663668      !! External gradient below BL needed both with and w/o FK 
    664       CALL zdf_osm_external_gradients( Kmm, nbld(A2D(0)) + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )   ! ag 19/03 
     669      CALL zdf_osm_external_gradients( Kmm, nbld(A2D((nn_hls-1))) + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )   ! ag 19/03 
    665670      ! 
    666671      ! Test if pycnocline well resolved 
    667       !      DO_2D( 0, 0, 0, 0 )                                         Removed with ag 19/03 changes. A change in eddy diffusivity/viscosity 
     672      !      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                                         Removed with ag 19/03 changes. A change in eddy diffusivity/viscosity 
    668673      !         IF (l_conv(ji,jj) ) THEN                                  should account for this. 
    669674      !            ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,nbld(ji,jj),Kmm) 
     
    683688      ! Recalculate bl averages using jk_ext & ml averages .... note no rotation of u & v here.. 
    684689      jk_ext(:,:) = 1   ! ag 19/03 
    685       CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D(0)), av_t_bl, av_s_bl,      & 
     690      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D((nn_hls-1))), av_t_bl, av_s_bl,      & 
    686691         &                           av_b_bl, av_u_bl, av_v_bl, jk_ext, av_dt_bl,   & 
    687692         &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) 
    688       jk_ext(:,:) = nbld(A2D(0)) - nmld(:,:) + jk_ext(:,:) + 1   ! ag 19/03 
     693      jk_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(:,:) + jk_ext(:,:) + 1   ! ag 19/03 
    689694      CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,          & 
    690695         &                           av_b_ml, av_u_ml, av_v_ml, jk_ext, av_dt_ml,   & 
     
    696701      ! Test if surface boundary layer coupled to bottom 
    697702      l_coup(:,:) = .FALSE.   ! ag 19/03 
    698       DO_2D( 0, 0, 0, 0 ) 
     703      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    699704         zhbl_t(ji,jj) = hbl(ji,jj) + ( zdhdt(ji,jj) - ww(ji,jj,nbld(ji,jj)) ) * rn_Dt   ! Certainly need ww here, so subtract it 
    700705         ! Adjustment to represent limiting by ocean bottom 
     
    708713      END_2D 
    709714      ! 
    710       nmld(:,:) = nbld(A2D(0))           ! use nmld to hold previous blayer index 
     715      nmld(:,:) = nbld(A2D((nn_hls-1)))           ! use nmld to hold previous blayer index 
    711716      nbld(:,:) = 4 
    712717      ! 
    713       DO_3D( 0, 0, 0, 0, 4, jpkm1 ) 
     718      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 4, jpkm1 ) 
    714719         IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 
    715720            nbld(ji,jj) = jk 
     
    726731      ! Recalculate BL averages and differences using new BL depth 
    727732      jk_ext(:,:) = 1   ! ag 19/03 
    728       CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D(0)), av_t_bl, av_s_bl,      & 
     733      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D((nn_hls-1))), av_t_bl, av_s_bl,      & 
    729734         &                           av_b_bl, av_u_bl, av_v_bl, jk_ext, av_dt_bl,   & 
    730735         &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl ) 
     
    734739      ! 
    735740      ! Reset l_pyc before calculating terms in the flux-gradient relationship 
    736       DO_2D( 0, 0, 0, 0 ) 
     741      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    737742         IF ( av_db_bl(ji,jj) < rn_osm_bl_thresh .OR. nbld(ji,jj) >= mbkt(ji,jj) - 2 .OR.   & 
    738743            & nbld(ji,jj) - nmld(ji,jj) == 1   .OR. zdhdt(ji,jj) < 0.0_wp ) THEN   ! ag 19/03 
     
    748753      END_2D 
    749754      ! 
    750       dstokes(:,:) = MIN ( dstokes(:,:), hbl(A2D(0))/ 3.0_wp )   ! Limit delta for shallow boundary layers for calculating 
     755      dstokes(:,:) = MIN ( dstokes(:,:), hbl(A2D((nn_hls-1)))/ 3.0_wp )   ! Limit delta for shallow boundary layers for calculating 
    751756      !                                                       !    flux-gradient terms 
    752757      ! 
     
    754759      !      jk_ext = nbld - nmld + 1 
    755760      ! Recalculate ML averages and differences using new ML depth 
    756       jk_ext(:,:) = nbld(A2D(0)) - nmld(A2D(0)) + jk_ext(:,:) + 1   ! ag 19/03 
     761      jk_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(A2D((nn_hls-1))) + jk_ext(:,:) + 1   ! ag 19/03 
    757762      CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,    & 
    758763         &                           av_b_ml, av_u_ml, av_v_ml, jk_ext, av_dt_ml,   & 
    759764         &                           av_ds_ml, av_db_ml, av_du_ml, av_dv_ml ) 
    760765      ! 
    761       CALL zdf_osm_external_gradients( Kmm, nbld(A2D(0)) + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
     766      CALL zdf_osm_external_gradients( Kmm, nbld(A2D((nn_hls-1))) + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
    762767      ! Rotate mean currents and changes onto wind aligned co-ordinates 
    763768      CALL zdf_osm_velocity_rotation( av_u_ml, av_v_ml ) 
     
    786791      ! 
    787792      ! Rotate non-gradient velocity terms back to model reference frame 
    788       CALL zdf_osm_velocity_rotation( ghamu(A2D(0),:), ghamv(A2D(0),:), .FALSE., 2, nbld(A2D(0)) ) 
     793      CALL zdf_osm_velocity_rotation( ghamu(A2D((nn_hls-1)),:), ghamv(A2D((nn_hls-1)),:), .FALSE., 2, nbld(A2D((nn_hls-1))) ) 
    789794      ! 
    790795      ! KPP-style Ri# mixing 
    791796      IF ( ln_kpprimix ) THEN 
    792797         jkflt = jpk 
    793          DO_2D( 0, 0, 0, 0 ) 
     798         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    794799            IF ( nbld(ji,jj) < jkflt ) jkflt = nbld(ji,jj) 
    795800         END_2D 
    796801         DO jk = jkflt+1, jpkm1 
    797802            ! Shear production at uw- and vw-points (energy conserving form) 
    798             DO_2D( 1, 0, 1, 0 ) 
     803            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    799804               z2du(ji,jj) = 0.5_wp * ( uu(ji,jj,jk-1,Kmm) - uu(ji,jj,jk,Kmm) ) * ( uu(ji,jj,jk-1,Kbb) - uu(ji,jj,jk,Kbb) ) *   & 
    800805                  &          wumask(ji,jj,jk) / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) 
     
    802807                  &          wvmask(ji,jj,jk) / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) 
    803808            END_2D 
    804             DO_2D( 0, 0, 0, 0 ) 
     809            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    805810               IF ( jk > nbld(ji,jj) ) THEN 
    806811                  ! Shear prod. at w-point weightened by mask 
     
    821826      ! KPP-style set diffusivity large if unstable below BL 
    822827      IF ( ln_convmix) THEN 
    823          DO_2D( 0, 0, 0, 0 ) 
     828         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    824829            DO jk = nbld(ji,jj) + 1, jpkm1 
    825830               IF ( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1e-12_wp ) zdiffut(ji,jj,jk) = MAX( rn_difconv, zdiffut(ji,jj,jk) ) 
     
    829834      ! 
    830835      IF ( ln_osm_mle ) THEN   ! Set up diffusivity and non-gradient mixing 
    831          DO_2D( 0, 0, 0, 0 ) 
     836         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    832837            IF ( l_flux(ji,jj) ) THEN   ! MLE mixing extends below boundary layer 
    833838               ! Calculate MLE flux contribution from surface fluxes 
     
    862867      ! CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1.0_wp ) 
    863868      ! GN 25/8: need to change tmask --> wmask 
    864       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     869      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    865870         p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk) 
    866871         p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk) 
     
    868873      ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and 
    869874      !    v grids 
    870       CALL lbc_lnk( 'zdfosm', p_avt, 'W', 1.0_wp,   & 
    871          &                    p_avm, 'W', 1.0_wp,   & 
    872          &                    ghamu, 'W', 1.0_wp,   & 
    873          &                    ghamv, 'W', 1.0_wp ) 
     875      IF ( nn_hls == 1 ) CALL lbc_lnk( 'zdfosm', ghamu, 'W', 1.0_wp,   & 
     876         &                                       ghamv, 'W', 1.0_wp ) 
    874877      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    875878         ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) / MAX( 1.0_wp, tmask(ji,jj,jk) + tmask (ji+1,jj,jk) ) *   & 
     
    889892         CASE(0:1) 
    890893            IF ( iom_use("us_x") ) THEN                                                           ! x surface Stokes drift 
    891                osmdia2d(A2D(0)) = tmask(A2D(0),1) * sustke * scos_wind 
     894               osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustke * scos_wind 
    892895               CALL iom_put( "us_x", osmdia2d ) 
    893896            END IF 
    894897            IF ( iom_use("us_y") ) THEN                                                           ! y surface Stokes drift 
    895                osmdia2d(A2D(0)) = tmask(A2D(0),1) * sustke * ssin_wind 
     898               osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustke * ssin_wind 
    896899               CALL iom_put( "us_y", osmdia2d ) 
    897900            END IF 
    898901            IF ( iom_use("wind_wave_abs_power") ) THEN 
    899                osmdia2d(A2D(0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**2 * sustke 
     902               osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**2 * sustke 
    900903               CALL iom_put( "wind_wave_abs_power", osmdia2d ) 
    901904            END IF 
     
    913916            IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) )                      ! U_10 
    914917            IF ( iom_use("wind_wave_abs_power") ) THEN 
    915                osmdia2d(A2D(0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**2 * SQRT( ut0sd(A2D(0))**2 + vt0sd(A2D(0))**2 ) 
     918               osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**2 * SQRT( ut0sd(A2D((nn_hls-1)))**2 + vt0sd(A2D((nn_hls-1)))**2 ) 
    916919               CALL iom_put( "wind_wave_abs_power", osmdia2d ) 
    917920            END IF 
     
    922925         IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv )                      ! <vw_NL> 
    923926         IF ( iom_use("zwth0") ) THEN                                                      ! <Tw_0> 
    924             osmdia2d(A2D(0)) = tmask(A2D(0),1) * swth0;     CALL iom_put( "zwth0",     osmdia2d ) 
     927            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swth0;     CALL iom_put( "zwth0",     osmdia2d ) 
    925928         END IF 
    926929         IF ( iom_use("zws0") ) THEN                                                       ! <Sw_0> 
    927             osmdia2d(A2D(0)) = tmask(A2D(0),1) * sws0;      CALL iom_put( "zws0",      osmdia2d ) 
     930            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sws0;      CALL iom_put( "zws0",      osmdia2d ) 
    928931         END IF 
    929932         IF ( iom_use("zwb0") ) THEN                                                       ! <Sw_0> 
    930             osmdia2d(A2D(0)) = tmask(A2D(0),1) * swb0;      CALL iom_put( "zwb0",      osmdia2d ) 
     933            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swb0;      CALL iom_put( "zwb0",      osmdia2d ) 
    931934         END IF 
    932935         IF ( iom_use("zwbav") ) THEN                                                      ! Upward BL-avged turb buoyancy flux 
    933             osmdia2d(A2D(0)) = tmask(A2D(0),1) * swth0;     CALL iom_put( "zwbav",     osmdia2d ) 
     936            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swth0;     CALL iom_put( "zwbav",     osmdia2d ) 
    934937         END IF 
    935938         IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                     ! Boundary-layer depth 
    936939         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*nbld )                  ! Boundary-layer max k 
    937940         IF ( iom_use("zdt_bl") ) THEN                                                     ! dt at ml base 
    938             osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_dt_bl;  CALL iom_put( "zdt_bl", osmdia2d ) 
     941            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_dt_bl;  CALL iom_put( "zdt_bl", osmdia2d ) 
    939942         END IF 
    940943         IF ( iom_use("zds_bl") ) THEN                                                     ! ds at ml base 
    941             osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_ds_bl;  CALL iom_put( "zds_bl", osmdia2d ) 
     944            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_ds_bl;  CALL iom_put( "zds_bl", osmdia2d ) 
    942945         END IF 
    943946         IF ( iom_use("zdb_bl") ) THEN                                                     ! db at ml base 
    944             osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_db_bl;  CALL iom_put( "zdb_bl", osmdia2d ) 
     947            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_db_bl;  CALL iom_put( "zdb_bl", osmdia2d ) 
    945948         END IF 
    946949         IF ( iom_use("zdu_bl") ) THEN                                                     ! du at ml base 
    947             osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_du_bl;  CALL iom_put( "zdu_bl", osmdia2d ) 
     950            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_du_bl;  CALL iom_put( "zdu_bl", osmdia2d ) 
    948951         END IF 
    949952         IF ( iom_use("zdv_bl") ) THEN                                                     ! dv at ml base 
    950             osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_dv_bl;  CALL iom_put( "zdv_bl", osmdia2d ) 
     953            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_dv_bl;  CALL iom_put( "zdv_bl", osmdia2d ) 
    951954         END IF 
    952955         IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh )                        ! Initial boundary-layer depth 
    953956         IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml )                     ! Initial boundary-layer depth 
    954957         IF ( iom_use("zdt_ml") ) THEN                                                     ! dt at ml base 
    955             osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_dt_ml;  CALL iom_put( "zdt_ml", osmdia2d ) 
     958            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_dt_ml;  CALL iom_put( "zdt_ml", osmdia2d ) 
    956959         END IF 
    957960         IF ( iom_use("zds_ml") ) THEN                                                     ! ds at ml base 
    958             osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_ds_ml;  CALL iom_put( "zds_ml", osmdia2d ) 
     961            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_ds_ml;  CALL iom_put( "zds_ml", osmdia2d ) 
    959962         END IF 
    960963         IF ( iom_use("zdb_ml") ) THEN                                                     ! db at ml base 
    961             osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_db_ml;  CALL iom_put( "zdb_ml", osmdia2d ) 
     964            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_db_ml;  CALL iom_put( "zdb_ml", osmdia2d ) 
    962965         END IF 
    963966         IF ( iom_use("dstokes") ) THEN                                                    ! Stokes drift penetration depth 
    964             osmdia2d(A2D(0)) = tmask(A2D(0),1) * dstokes;   CALL iom_put( "dstokes",   osmdia2d ) 
     967            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * dstokes;   CALL iom_put( "dstokes",   osmdia2d ) 
    965968         END IF 
    966969         IF ( iom_use("zustke") ) THEN                                                     ! Stokes drift magnitude at T-points 
    967             osmdia2d(A2D(0)) = tmask(A2D(0),1) * sustke;    CALL iom_put( "zustke",    osmdia2d ) 
     970            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustke;    CALL iom_put( "zustke",    osmdia2d ) 
    968971         END IF 
    969972         IF ( iom_use("zwstrc") ) THEN                                                     ! Convective velocity scale 
    970             osmdia2d(A2D(0)) = tmask(A2D(0),1) * swstrc;    CALL iom_put( "zwstrc",    osmdia2d ) 
     973            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swstrc;    CALL iom_put( "zwstrc",    osmdia2d ) 
    971974         END IF 
    972975         IF ( iom_use("zwstrl") ) THEN                                                     ! Langmuir velocity scale 
    973             osmdia2d(A2D(0)) = tmask(A2D(0),1) * swstrl;    CALL iom_put( "zwstrl",    osmdia2d ) 
     976            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swstrl;    CALL iom_put( "zwstrl",    osmdia2d ) 
    974977         END IF 
    975978         IF ( iom_use("zustar") ) THEN                                                     ! Friction velocity scale 
    976             osmdia2d(A2D(0)) = tmask(A2D(0),1) * sustar;    CALL iom_put( "zustar",    osmdia2d ) 
     979            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustar;    CALL iom_put( "zustar",    osmdia2d ) 
    977980         END IF 
    978981         IF ( iom_use("zvstr") ) THEN                                                      ! Mixed velocity scale 
    979             osmdia2d(A2D(0)) = tmask(A2D(0),1) * svstr;     CALL iom_put( "zvstr",     osmdia2d ) 
     982            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * svstr;     CALL iom_put( "zvstr",     osmdia2d ) 
    980983         END IF 
    981984         IF ( iom_use("zla") ) THEN                                                        ! Langmuir # 
    982             osmdia2d(A2D(0)) = tmask(A2D(0),1) * sla;       CALL iom_put( "zla",       osmdia2d ) 
     985            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sla;       CALL iom_put( "zla",       osmdia2d ) 
    983986         END IF 
    984987         IF ( iom_use("wind_power") ) THEN                                                 ! BL depth internal to zdf_osm routine 
    985             osmdia2d(A2D(0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**3 
     988            osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**3 
    986989            CALL iom_put( "wind_power", osmdia2d ) 
    987990         END IF 
    988991         IF ( iom_use("wind_wave_power") ) THEN 
    989             osmdia2d(A2D(0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**2 * sustke 
     992            osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**2 * sustke 
    990993            CALL iom_put( "wind_wave_power", osmdia2d ) 
    991994         END IF 
    992995         IF ( iom_use("zhbl") ) THEN                                                       ! BL depth internal to zdf_osm routine 
    993             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zhbl;      CALL iom_put( "zhbl",      osmdia2d ) 
     996            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zhbl;      CALL iom_put( "zhbl",      osmdia2d ) 
    994997         END IF 
    995998         IF ( iom_use("zhml") ) THEN                                                       ! ML depth internal to zdf_osm routine 
    996             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zhml;      CALL iom_put( "zhml",      osmdia2d ) 
     999            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zhml;      CALL iom_put( "zhml",      osmdia2d ) 
    9971000         END IF 
    9981001         IF ( iom_use("imld") ) THEN                                                       ! Index for ML depth internal to zdf_osm routine 
    999             osmdia2d(A2D(0)) = tmask(A2D(0),1) * nmld;      CALL iom_put( "imld",      osmdia2d ) 
     1002            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * nmld;      CALL iom_put( "imld",      osmdia2d ) 
    10001003         END IF 
    10011004         IF ( iom_use("jp_ext") ) THEN                                                     ! =1 if pycnocline resolved internal to zdf_osm routine 
    1002             osmdia2d(A2D(0)) = tmask(A2D(0),1) * jk_ext;    CALL iom_put( "jp_ext",    osmdia2d ) 
     1005            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * jk_ext;    CALL iom_put( "jp_ext",    osmdia2d ) 
    10031006         END IF 
    10041007         IF ( iom_use("j_ddh") ) THEN                                                      ! Index forpyc thicknessh internal to zdf_osm routine 
    1005             osmdia2d(A2D(0)) = tmask(A2D(0),1) * n_ddh;     CALL iom_put( "j_ddh",     osmdia2d ) 
     1008            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * n_ddh;     CALL iom_put( "j_ddh",     osmdia2d ) 
    10061009         END IF 
    10071010         IF ( iom_use("zshear") ) THEN                                                     ! Shear production of TKE internal to zdf_osm routine 
    1008             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zshear;    CALL iom_put( "zshear",    osmdia2d ) 
     1011            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zshear;    CALL iom_put( "zshear",    osmdia2d ) 
    10091012         END IF 
    10101013         IF ( iom_use("zdh") ) THEN                                                        ! Pyc thicknessh internal to zdf_osm routine 
    1011             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zdh;       CALL iom_put( "zdh",       osmdia2d ) 
     1014            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zdh;       CALL iom_put( "zdh",       osmdia2d ) 
    10121015         END IF 
    10131016         IF ( iom_use("zhol") ) THEN                                                       ! ML depth internal to zdf_osm routine 
    1014             osmdia2d(A2D(0)) = tmask(A2D(0),1) * shol;      CALL iom_put( "zhol",      osmdia2d ) 
     1017            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * shol;      CALL iom_put( "zhol",      osmdia2d ) 
    10151018         END IF 
    10161019         IF ( iom_use("zwb_ent") ) THEN                                                    ! Upward turb buoyancy entrainment flux 
    1017             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zwb_ent;   CALL iom_put( "zwb_ent",   osmdia2d ) 
     1020            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwb_ent;   CALL iom_put( "zwb_ent",   osmdia2d ) 
    10181021         END IF 
    10191022         IF ( iom_use("zt_ml") ) THEN                                                      ! Average T in ML 
    1020             osmdia2d(A2D(0)) = tmask(A2D(0),1) * av_t_ml;   CALL iom_put( "zt_ml",     osmdia2d ) 
     1023            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_t_ml;   CALL iom_put( "zt_ml",     osmdia2d ) 
    10211024         END IF 
    10221025         IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle )                  ! FK layer depth 
    10231026         IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld )                  ! FK target layer depth 
    10241027         IF ( iom_use("zwb_fk") ) THEN                                                     ! FK b flux 
    1025             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zwb_fk;    CALL iom_put( "zwb_fk",    osmdia2d ) 
     1028            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwb_fk;    CALL iom_put( "zwb_fk",    osmdia2d ) 
    10261029         END IF 
    10271030         IF ( iom_use("zwb_fk_b") ) THEN                                                   ! FK b flux averaged over ML 
    1028             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zwb_fk_b;  CALL iom_put( "zwb_fk_b",  osmdia2d ) 
     1031            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwb_fk_b;  CALL iom_put( "zwb_fk_b",  osmdia2d ) 
    10291032         END IF 
    10301033         IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )      ! FK layer max k 
     
    10361039         IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle )      ! FK dbdy at v-pt 
    10371040         IF ( iom_use("zdiff_mle") ) THEN                                                  ! FK diff in MLE at t-pt 
    1038             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zdiff_mle; CALL iom_put( "zdiff_mle", osmdia2d ) 
     1041            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zdiff_mle; CALL iom_put( "zdiff_mle", osmdia2d ) 
    10391042         END IF 
    10401043         IF ( iom_use("zvel_mle") ) THEN                                                   ! FK diff in MLE at t-pt 
    1041             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zdiff_mle; CALL iom_put( "zvel_mle",  osmdia2d ) 
     1044            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zdiff_mle; CALL iom_put( "zvel_mle",  osmdia2d ) 
    10421045         END IF 
    10431046      END IF 
     
    10601063      !!---------------------------------------------------------------------- 
    10611064      INTEGER,                     INTENT(in   )           ::   Kbb, Kmm   ! Ocean time-level indices 
    1062       INTEGER,  DIMENSION(A2D(0)), INTENT(in   )           ::   knlev      ! Number of levels to average over. 
    1063       REAL(wp), DIMENSION(A2D(0)), INTENT(  out)           ::   pt, ps     ! Average temperature and salinity 
    1064       REAL(wp), DIMENSION(A2D(0)), INTENT(  out)           ::   pb         ! Average buoyancy 
    1065       REAL(wp), DIMENSION(A2D(0)), INTENT(  out)           ::   pu, pv     ! Average current components 
    1066       INTEGER,  DIMENSION(A2D(0)), INTENT(in   ), OPTIONAL ::   kp_ext     ! External-level offsets 
    1067       REAL(wp), DIMENSION(A2D(0)), INTENT(  out), OPTIONAL ::   pdt, pds   ! Difference between average temperature, salinity, 
    1068       REAL(wp), DIMENSION(A2D(0)), INTENT(  out), OPTIONAL ::   pdb        ! buoyancy, 
    1069       REAL(wp), DIMENSION(A2D(0)), INTENT(  out), OPTIONAL ::   pdu, pdv   ! velocity components and the OSBL 
     1065      INTEGER,  DIMENSION(A2D((nn_hls-1))), INTENT(in   )           ::   knlev      ! Number of levels to average over. 
     1066      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out)           ::   pt, ps     ! Average temperature and salinity 
     1067      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out)           ::   pb         ! Average buoyancy 
     1068      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out)           ::   pu, pv     ! Average current components 
     1069      INTEGER,  DIMENSION(A2D((nn_hls-1))), INTENT(in   ), OPTIONAL ::   kp_ext     ! External-level offsets 
     1070      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out), OPTIONAL ::   pdt, pds   ! Difference between average temperature, salinity, 
     1071      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out), OPTIONAL ::   pdb        ! buoyancy, 
     1072      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out), OPTIONAL ::   pdu, pdv   ! velocity components and the OSBL 
    10701073      ! 
    10711074      INTEGER                     ::   jk, jkflt, jkmax, ji, jj   ! Loop indices 
    10721075      INTEGER                     ::   ibld_ext                   ! External-layer index 
    1073       REAL(wp), DIMENSION(A2D(0)) ::   zthick                     ! Layer thickness 
     1076      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zthick                     ! Layer thickness 
    10741077      REAL(wp)                    ::   zthermal, zbeta            ! Thermal/haline expansion/contraction coefficients 
    10751078      !!---------------------------------------------------------------------- 
     
    10831086      jkflt = jpk 
    10841087      jkmax = 0 
    1085       DO_2D( 0, 0, 0, 0 ) 
     1088      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    10861089         IF ( knlev(ji,jj) < jkflt ) jkflt = knlev(ji,jj) 
    10871090         IF ( knlev(ji,jj) > jkmax ) jkmax = knlev(ji,jj) 
    10881091      END_2D 
    1089       DO_3D( 0, 0, 0, 0, 2, jkflt )   ! Upper, flat part of layer 
     1092      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkflt )   ! Upper, flat part of layer 
    10901093         zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm) 
    10911094         pt(ji,jj)     = pt(ji,jj)     + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) 
     
    10981101            &                               MAX( 1.0_wp , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) )          
    10991102      END_3D 
    1100       DO_3D( 0, 0, 0, 0, jkflt+1, jkmax )   ! Lower, non-flat part of layer 
     1103      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jkflt+1, jkmax )   ! Lower, non-flat part of layer 
    11011104         IF ( knlev(ji,jj) >= jk ) THEN 
    11021105            zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm) 
     
    11111114         END IF 
    11121115      END_3D 
    1113       DO_2D( 0, 0, 0, 0 ) 
     1116      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    11141117         pt(ji,jj) = pt(ji,jj) / zthick(ji,jj) 
    11151118         ps(ji,jj) = ps(ji,jj) / zthick(ji,jj) 
     
    11231126      ! Differences between vertical averages and values at an external layer 
    11241127      IF ( PRESENT( kp_ext ) ) THEN 
    1125          DO_2D( 0, 0, 0, 0 ) 
     1128         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    11261129            ibld_ext = knlev(ji,jj) + kp_ext(ji,jj) 
    11271130            IF ( ibld_ext <= mbkt(ji,jj)-1 ) THEN   ! ag 09/03 
     
    11601163      !! 
    11611164      !!----------------------------------------------------------------------       
    1162       REAL(wp),           INTENT(inout), DIMENSION(A2D(0)) ::   pu, pv           ! Components of current 
     1165      REAL(wp),           INTENT(inout), DIMENSION(A2D((nn_hls-1))) ::   pu, pv           ! Components of current 
    11631166      LOGICAL,  OPTIONAL, INTENT(in   )                    ::   fwd              ! Forward (default) or reverse rotation 
    11641167      ! 
     
    11681171      zfwd = 1.0_wp 
    11691172      IF( PRESENT(fwd) .AND. ( .NOT. fwd ) ) zfwd = -1.0_wp 
    1170       DO_2D( 0, 0, 0, 0 ) 
     1173      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    11711174         ztmp      = pu(ji,jj) 
    11721175         pu(ji,jj) = pu(ji,jj) * scos_wind(ji,jj) + zfwd * pv(ji,jj) * ssin_wind(ji,jj) 
     
    11901193      !! 
    11911194      !!----------------------------------------------------------------------       
    1192       REAL(wp),           INTENT(inout), DIMENSION(A2D(0),jpk) ::   pu, pv           ! Components of current 
     1195      REAL(wp),           INTENT(inout), DIMENSION(A2D((nn_hls-1)),jpk) ::   pu, pv           ! Components of current 
    11931196      LOGICAL,  OPTIONAL, INTENT(in   )                        ::   fwd              ! Forward (default) or reverse rotation 
    11941197      INTEGER,  OPTIONAL, INTENT(in   )                        ::   ktop             ! Minimum depth index 
    1195       INTEGER,  OPTIONAL, INTENT(in   ), DIMENSION(A2D(0))     ::   knlev            ! Array of maximum depth indices 
     1198      INTEGER,  OPTIONAL, INTENT(in   ), DIMENSION(A2D((nn_hls-1)))     ::   knlev            ! Array of maximum depth indices 
    11961199      !  
    11971200      INTEGER  ::   ji, jj, jk, jktop, jkmax   ! Loop indices 
     
    12051208      IF( PRESENT(knlev) ) THEN 
    12061209         jkmax = 0 
    1207          DO_2D( 0, 0, 0, 0 ) 
     1210         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    12081211            IF ( knlev(ji,jj) > jkmax ) jkmax = knlev(ji,jj) 
    12091212         END_2D 
     
    12131216         llkbot = .TRUE. 
    12141217      END IF 
    1215       DO_3D( 0, 0, 0, 0, jktop, jkmax ) 
     1218      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jktop, jkmax ) 
    12161219         IF ( llkbot .OR. knlev(ji,jj) >= jk ) THEN 
    12171220            ztmp         = pu(ji,jj,jk) 
     
    12371240      !!---------------------------------------------------------------------- 
    12381241      INTEGER,                     INTENT(in   ) ::   Kmm       ! Ocean time-level index 
    1239       REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pwb_ent   ! Buoyancy fluxes at base 
    1240       REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pwb_min   !    of well-mixed layer 
    1241       REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pshear    ! Production of TKE due to shear across the pycnocline 
    1242       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl      ! BL depth 
    1243       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phml      ! ML depth 
    1244       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdh       ! Pycnocline depth 
     1242      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out) ::   pwb_ent   ! Buoyancy fluxes at base 
     1243      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out) ::   pwb_min   !    of well-mixed layer 
     1244      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out) ::   pshear    ! Production of TKE due to shear across the pycnocline 
     1245      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phbl      ! BL depth 
     1246      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phml      ! ML depth 
     1247      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pdh       ! Pycnocline depth 
    12451248      ! 
    12461249      ! Local variables 
    12471250      INTEGER :: jj, ji   ! Loop indices 
    12481251      ! 
    1249       REAL(wp), DIMENSION(A2D(0)) ::   zekman 
    1250       REAL(wp), DIMENSION(A2D(0)) ::   zri_p, zri_b   ! Richardson numbers 
     1252      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zekman 
     1253      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zri_p, zri_b   ! Richardson numbers 
    12511254      REAL(wp)                    ::   zshear_u, zshear_v, zwb_shr 
    12521255      REAL(wp)                    ::   zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes 
     
    12691272      ! 
    12701273      ! Determins stability and set flag l_conv 
    1271       DO_2D( 0, 0, 0, 0 ) 
     1274      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    12721275         IF ( shol(ji,jj) < 0.0_wp ) THEN 
    12731276            l_conv(ji,jj) = .TRUE. 
     
    12781281      ! 
    12791282      pshear(:,:) = 0.0_wp 
    1280       zekman(:,:) = EXP( -1.0_wp * pp_ek * ABS( ff_t(A2D(0)) ) * phbl(:,:) / MAX( sustar(A2D(0)), 1.e-8 ) ) 
    1281       ! 
    1282       DO_2D( 0, 0, 0, 0 ) 
     1283      zekman(:,:) = EXP( -1.0_wp * pp_ek * ABS( ff_t(A2D((nn_hls-1))) ) * phbl(:,:) / MAX( sustar(A2D((nn_hls-1))), 1.e-8 ) ) 
     1284      ! 
     1285      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    12831286         IF ( l_conv(ji,jj) ) THEN 
    12841287            IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN 
     
    13321335      ! 
    13331336      ! Calculate entrainment buoyancy flux due to surface fluxes. 
    1334       DO_2D( 0, 0, 0, 0 ) 
     1337      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    13351338         IF ( l_conv(ji,jj) ) THEN 
    13361339            zwcor        = ABS( ff_t(ji,jj) ) * phbl(ji,jj) + epsln 
     
    13531356      END_2D 
    13541357      ! 
    1355       DO_2D( 0, 0, 0, 0 ) 
     1358      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    13561359         IF ( l_shear(ji,jj) ) THEN 
    13571360            IF ( l_conv(ji,jj) ) THEN 
     
    13921395      !!----------------------------------------------------------------------    
    13931396      INTEGER,                     INTENT(in   ) ::   Kmm                   ! Ocean time-level index 
    1394       INTEGER,  DIMENSION(A2D(0)), INTENT(in   ) ::   kbase                 ! OSBL base layer index 
    1395       REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pdtdz, pdsdz, pdbdz   ! External gradients of temperature, salinity and buoyancy 
     1397      INTEGER,  DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   kbase                 ! OSBL base layer index 
     1398      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out) ::   pdtdz, pdsdz, pdbdz   ! External gradients of temperature, salinity and buoyancy 
    13961399      ! 
    13971400      ! Local variables 
     
    14051408      pdbdz(:,:) = pp_large 
    14061409      ! 
    1407       DO_2D( 0, 0, 0, 0 ) 
     1410      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    14081411         IF ( kbase(ji,jj)+1 < mbkt(ji,jj) ) THEN 
    14091412            zthermal = rab_n(ji,jj,1,jp_tem)   ! Ideally use nbld not 1?? 
     
    14331436      !! 
    14341437      !!---------------------------------------------------------------------- 
    1435       REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pdhdt          ! Rate of change of hbl 
    1436       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl           ! BL depth 
    1437       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdh            ! Pycnocline depth 
    1438       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
    1439       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_min 
    1440       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
    1441       REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL 
    1442       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_fk         ! Max MLE buoyancy flux 
    1443       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pvel_mle       ! Vvelocity scale for dhdt with stable ML and FK 
     1438      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out) ::   pdhdt          ! Rate of change of hbl 
     1439      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phbl           ! BL depth 
     1440      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pdh            ! Pycnocline depth 
     1441      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
     1442      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_min 
     1443      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
     1444      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL 
     1445      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_fk         ! Max MLE buoyancy flux 
     1446      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pvel_mle       ! Vvelocity scale for dhdt with stable ML and FK 
    14441447      ! 
    14451448      ! Local variables 
     
    14551458      pwb_fk_b(:,:) = pp_large 
    14561459      ! 
    1457       DO_2D( 0, 0, 0, 0 ) 
     1460      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    14581461         ! 
    14591462         IF ( l_shear(ji,jj) ) THEN 
     
    15981601      !!---------------------------------------------------------------------- 
    15991602      INTEGER,                     INTENT(in   ) ::   Kmm        ! Ocean time-level index 
    1600       REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pdhdt      ! Rates of change of hbl 
    1601       REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   phbl       ! BL depth 
    1602       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl_t     ! BL depth 
    1603       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_ent    ! Buoyancy entrainment flux 
    1604       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_fk_b   ! MLE buoyancy flux averaged over OSBL 
     1603      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   pdhdt      ! Rates of change of hbl 
     1604      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   phbl       ! BL depth 
     1605      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phbl_t     ! BL depth 
     1606      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_ent    ! Buoyancy entrainment flux 
     1607      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_fk_b   ! MLE buoyancy flux averaged over OSBL 
    16051608      ! 
    16061609      ! Local variables 
     
    16091612      REAL(wp) ::   zthermal, zbeta 
    16101613      ! 
    1611       DO_2D( 0, 0, 0, 0 ) 
     1614      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    16121615         IF ( nbld(ji,jj) - nmld(ji,jj) > 1 ) THEN 
    16131616            ! 
     
    16981701      !!---------------------------------------------------------------------- 
    16991702      INTEGER,                     INTENT(in   ) ::   Kmm            ! Ocean time-level index 
    1700       REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pdh            ! Pycnocline thickness 
    1701       REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   phml           ! ML depth 
    1702       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdhdt          ! BL depth tendency 
    1703       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl           ! BL depth 
    1704       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
    1705       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
    1706       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL 
     1703      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   pdh            ! Pycnocline thickness 
     1704      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   phml           ! ML depth 
     1705      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pdhdt          ! BL depth tendency 
     1706      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phbl           ! BL depth 
     1707      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
     1708      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
     1709      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL 
    17071710 
    17081711      ! 
     
    17151718      REAL, PARAMETER ::   pp_ddh = 2.5_wp, pp_ddh_2 = 3.5_wp   ! Also in pycnocline_depth 
    17161719      ! 
    1717       DO_2D( 0, 0, 0, 0 ) 
     1720      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    17181721         ! 
    17191722         IF ( l_shear(ji,jj) ) THEN 
     
    18591862      !!---------------------------------------------------------------------- 
    18601863      INTEGER,                          INTENT(in   ) ::   Kmm            ! Ocean time-level index 
    1861       INTEGER,  DIMENSION(A2D(0)),      INTENT(in   ) ::   kp_ext         ! External-level offsets 
    1862       REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(  out) ::   pdbdz          ! Gradients in the pycnocline 
    1863       REAL(wp), DIMENSION(A2D(0)),      INTENT(  out) ::   palpha 
    1864       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdh            ! Pycnocline thickness 
    1865       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phbl           ! BL depth 
    1866       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
    1867       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phml           ! ML depth 
    1868       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdhdt          ! Rates of change of hbl 
     1864      INTEGER,  DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   kp_ext         ! External-level offsets 
     1865      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk),  INTENT(  out) ::   pdbdz          ! Gradients in the pycnocline 
     1866      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(  out) ::   palpha 
     1867      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdh            ! Pycnocline thickness 
     1868      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   phbl           ! BL depth 
     1869      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
     1870      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   phml           ! ML depth 
     1871      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdhdt          ! Rates of change of hbl 
    18691872      ! 
    18701873      ! Local variables 
     
    18811884      palpha(:,:)  = pp_large 
    18821885      ! 
    1883       DO_2D( 0, 0, 0, 0 ) 
     1886      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    18841887         ! 
    18851888         IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
     
    19501953      IF ( ln_dia_pyc_scl ) THEN   ! Output of pycnocline gradient profiles 
    19511954         IF ( iom_use("zdbdz_pyc") ) THEN 
    1952             osmdia3d(A2D(0),:) = wmask(A2D(0),:) * pdbdz(:,:,:); CALL iom_put( "zdbdz_pyc", osmdia3d ) 
     1955            osmdia3d(A2D((nn_hls-1)),:) = wmask(A2D((nn_hls-1)),:) * pdbdz(:,:,:); CALL iom_put( "zdbdz_pyc", osmdia3d ) 
    19531956         END IF 
    19541957      END IF 
     
    19691972      !!---------------------------------------------------------------------- 
    19701973      INTEGER,                          INTENT(in   ) ::   Kbb, Kmm       ! Ocean time-level indices 
    1971       REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(inout) ::   pdiffut        ! t-diffusivity 
    1972       REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(inout) ::   pviscos        ! Viscosity 
    1973       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phbl           ! BL depth 
    1974       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phml           ! ML depth 
    1975       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdh            ! Pycnocline depth 
    1976       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdhdt          ! BL depth tendency 
    1977       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pshear         ! Shear production 
    1978       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
    1979       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pwb_min 
     1974      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk),  INTENT(inout) ::   pdiffut        ! t-diffusivity 
     1975      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk),  INTENT(inout) ::   pviscos        ! Viscosity 
     1976      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   phbl           ! BL depth 
     1977      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   phml           ! ML depth 
     1978      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdh            ! Pycnocline depth 
     1979      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdhdt          ! BL depth tendency 
     1980      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pshear         ! Shear production 
     1981      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
     1982      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pwb_min 
    19801983      ! 
    19811984      ! Local variables 
     
    19831986      ! 
    19841987      ! Scales used to calculate eddy diffusivity and viscosity profiles 
    1985       REAL(wp), DIMENSION(A2D(0)) ::   zdifml_sc,    zvisml_sc 
    1986       REAL(wp), DIMENSION(A2D(0)) ::   zdifpyc_n_sc, zdifpyc_s_sc 
    1987       REAL(wp), DIMENSION(A2D(0)) ::   zvispyc_n_sc, zvispyc_s_sc 
    1988       REAL(wp), DIMENSION(A2D(0)) ::   zbeta_d_sc,   zbeta_v_sc 
    1989       REAL(wp), DIMENSION(A2D(0)) ::   zb_coup,      zc_coup_vis,  zc_coup_dif 
     1988      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zdifml_sc,    zvisml_sc 
     1989      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zdifpyc_n_sc, zdifpyc_s_sc 
     1990      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zvispyc_n_sc, zvispyc_s_sc 
     1991      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zbeta_d_sc,   zbeta_v_sc 
     1992      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zb_coup,      zc_coup_vis,  zc_coup_dif 
    19901993      ! 
    19911994      REAL(wp) ::   zvel_sc_pyc, zvel_sc_ml, zstab_fac, zz_b 
     
    20012004      zb_coup(:,:) = 0.0_wp 
    20022005      ! 
    2003       DO_2D( 0, 0, 0, 0 ) 
     2006      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    20042007         IF ( l_conv(ji,jj) ) THEN 
    20052008            ! 
     
    20802083      END_2D 
    20812084      ! 
    2082       DO_2D( 0, 0, 0, 0 ) 
     2085      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    20832086         IF ( l_conv(ji,jj) ) THEN 
    20842087            DO jk = 2, nmld(ji,jj)   ! Mixed layer diffusivity 
     
    21502153      END_2D 
    21512154      IF( iom_use("pb_coup") ) THEN 
    2152          osmdia2d(A2D(0)) = tmask(A2D(0),1) * zb_coup(:,:); CALL iom_put( "pb_coup", osmdia2d )   ! BBL-coupling velocity scale 
     2155         osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zb_coup(:,:); CALL iom_put( "pb_coup", osmdia2d )   ! BBL-coupling velocity scale 
    21532156      END IF 
    21542157      ! 
     
    21672170      !!---------------------------------------------------------------------- 
    21682171      INTEGER,                          INTENT(in   ) ::   Kmm            ! Time-level index 
    2169       INTEGER,  DIMENSION(A2D(0)),      INTENT(in   ) ::   kp_ext         ! Offset for external level 
    2170       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phbl           ! BL depth 
    2171       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   phml           ! ML depth 
    2172       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdh            ! Pycnocline depth 
    2173       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdhdt          ! BL depth tendency 
    2174       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pshear         ! Shear production 
    2175       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdtdz_bl_ext   ! External temperature gradients 
    2176       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdsdz_bl_ext   ! External salinity gradients 
    2177       REAL(wp), DIMENSION(A2D(0)),      INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
    2178       REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(in   ) ::   pdiffut        ! t-diffusivity 
    2179       REAL(wp), DIMENSION(A2D(0),jpk),  INTENT(in   ) ::   pviscos        ! Viscosity 
    2180       ! 
    2181       REAL(wp), DIMENSION(A2D(0))     ::   zalpha_pyc   ! 
    2182       REAL(wp), DIMENSION(A2D(0),jpk) ::   zdbdz_pyc    ! Parametrised gradient of buoyancy in the pycnocline 
     2172      INTEGER,  DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   kp_ext         ! Offset for external level 
     2173      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   phbl           ! BL depth 
     2174      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   phml           ! ML depth 
     2175      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdh            ! Pycnocline depth 
     2176      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdhdt          ! BL depth tendency 
     2177      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pshear         ! Shear production 
     2178      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdtdz_bl_ext   ! External temperature gradients 
     2179      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdsdz_bl_ext   ! External salinity gradients 
     2180      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
     2181      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk),  INTENT(in   ) ::   pdiffut        ! t-diffusivity 
     2182      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk),  INTENT(in   ) ::   pviscos        ! Viscosity 
     2183      ! 
     2184      REAL(wp), DIMENSION(A2D((nn_hls-1)))     ::   zalpha_pyc   ! 
     2185      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk) ::   zdbdz_pyc    ! Parametrised gradient of buoyancy in the pycnocline 
    21832186      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   z3ddz_pyc_1, z3ddz_pyc_2   ! Pycnocline gradient/shear profiles 
    21842187      ! 
     
    21862189      INTEGER                     ::   istat                                   ! Memory allocation status 
    21872190      REAL(wp)                    ::   zznd_d, zznd_ml, zznd_pyc, znd          ! Temporary non-dimensional depths 
    2188       REAL(wp), DIMENSION(A2D(0)) ::   zsc_wth_1,zsc_ws_1                      ! Temporary scales 
    2189       REAL(wp), DIMENSION(A2D(0)) ::   zsc_uw_1, zsc_uw_2                      ! Temporary scales 
    2190       REAL(wp), DIMENSION(A2D(0)) ::   zsc_vw_1, zsc_vw_2                      ! Temporary scales 
    2191       REAL(wp), DIMENSION(A2D(0)) ::   ztau_sc_u                               ! Dissipation timescale at base of WML 
     2191      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zsc_wth_1,zsc_ws_1                      ! Temporary scales 
     2192      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zsc_uw_1, zsc_uw_2                      ! Temporary scales 
     2193      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zsc_vw_1, zsc_vw_2                      ! Temporary scales 
     2194      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   ztau_sc_u                               ! Dissipation timescale at base of WML 
    21922195      REAL(wp)                    ::   zbuoy_pyc_sc, zdelta_pyc                ! 
    21932196      REAL(wp)                    ::   zl_c,zl_l,zl_eps                        ! Used to calculate turbulence length scale 
    2194       REAL(wp), DIMENSION(A2D(0)) ::   za_cubic, zb_cubic                      ! Coefficients in cubic polynomial specifying 
    2195       REAL(wp), DIMENSION(A2D(0)) ::   zc_cubic, zd_cubic                      ! diffusivity in pycnocline 
    2196       REAL(wp), DIMENSION(A2D(0)) ::   zwt_pyc_sc_1, zws_pyc_sc_1              ! 
    2197       REAL(wp), DIMENSION(A2D(0)) ::   zzeta_pyc                               ! 
     2197      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   za_cubic, zb_cubic                      ! Coefficients in cubic polynomial specifying 
     2198      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zc_cubic, zd_cubic                      ! diffusivity in pycnocline 
     2199      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zwt_pyc_sc_1, zws_pyc_sc_1              ! 
     2200      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zzeta_pyc                               ! 
    21982201      REAL(wp)                    ::   zomega, zvw_max                         ! 
    2199       REAL(wp), DIMENSION(A2D(0)) ::   zuw_bse,zvw_bse                         ! Momentum, heat, and salinity fluxes 
    2200       REAL(wp), DIMENSION(A2D(0)) ::   zwth_ent,zws_ent                        ! at the top of the pycnocline 
    2201       REAL(wp), DIMENSION(A2D(0)) ::   zsc_wth_pyc, zsc_ws_pyc                 ! Scales for pycnocline transport term 
     2202      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zuw_bse,zvw_bse                         ! Momentum, heat, and salinity fluxes 
     2203      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zwth_ent,zws_ent                        ! at the top of the pycnocline 
     2204      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zsc_wth_pyc, zsc_ws_pyc                 ! Scales for pycnocline transport term 
    22022205      REAL(wp)                    ::   ztmp                                    ! 
    22032206      REAL(wp)                    ::   ztgrad, zsgrad, zbgrad                  ! Variables used to calculate pycnocline gradients 
     
    22202223      jkf_mld = jpk 
    22212224      jkm_mld = 0 
    2222       DO_2D( 0, 0, 0, 0 ) 
     2225      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    22232226         IF ( nbld(ji,jj) > jkm_bld ) jkm_bld = nbld(ji,jj) 
    22242227         IF ( nmld(ji,jj) < jkf_mld ) jkf_mld = nmld(ji,jj) 
     
    22282231      ! Stokes term in scalar flux, flux-gradient relationship 
    22292232      ! ------------------------------------------------------ 
    2230       WHERE ( l_conv(A2D(0)) ) 
    2231          zsc_wth_1(:,:) = swstrl(A2D(0))**3 * swth0(A2D(0)) / ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln ) 
    2232          zsc_ws_1(:,:)  = swstrl(A2D(0))**3 * sws0(A2D(0))  / ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln ) 
     2233      WHERE ( l_conv(A2D((nn_hls-1))) ) 
     2234         zsc_wth_1(:,:) = swstrl(A2D((nn_hls-1)))**3 * swth0(A2D((nn_hls-1))) / ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln ) 
     2235         zsc_ws_1(:,:)  = swstrl(A2D((nn_hls-1)))**3 * sws0(A2D((nn_hls-1)))  / ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln ) 
    22332236      ELSEWHERE 
    2234          zsc_wth_1(:,:) = 2.0_wp * swthav(A2D(0)) 
    2235          zsc_ws_1(:,:)  = 2.0_wp * swsav(A2D(0)) 
     2237         zsc_wth_1(:,:) = 2.0_wp * swthav(A2D((nn_hls-1))) 
     2238         zsc_ws_1(:,:)  = 2.0_wp * swsav(A2D((nn_hls-1))) 
    22362239      ENDWHERE 
    2237       DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
     2240      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 
    22382241         IF ( l_conv(ji,jj) ) THEN 
    22392242            IF ( jk <= nmld(ji,jj) ) THEN 
     
    22632266      ! svstr since term needs to go to zero as swstrl goes to zero) 
    22642267      ! --------------------------------------------------------------------- 
    2265       WHERE ( l_conv(A2D(0)) ) 
    2266          zsc_uw_1(:,:) = ( swstrl(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 )**pthird * sustke(A2D(0)) /   & 
    2267             &                                  MAX( ( 1.0_wp - 1.0_wp * 6.5_wp * sla(A2D(0))**( 8.0_wp / 3.0_wp ) ), 0.2_wp ) 
    2268          zsc_uw_2(:,:) = ( swstrl(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 )**pthird * sustke(A2D(0)) /   & 
    2269             &                                  MIN( sla(A2D(0))**( 8.0_wp / 3.0_wp ) + epsln, 0.12_wp ) 
    2270          zsc_vw_1(:,:) = ff_t(A2D(0)) * phml(A2D(0)) * sustke(A2D(0))**3 * MIN( sla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) /   & 
    2271             &            ( ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 )**( 2.0_wp / 3.0_wp ) + epsln ) 
     2268      WHERE ( l_conv(A2D((nn_hls-1))) ) 
     2269         zsc_uw_1(:,:) = ( swstrl(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 )**pthird * sustke(A2D((nn_hls-1))) /   & 
     2270            &                                  MAX( ( 1.0_wp - 1.0_wp * 6.5_wp * sla(A2D((nn_hls-1)))**( 8.0_wp / 3.0_wp ) ), 0.2_wp ) 
     2271         zsc_uw_2(:,:) = ( swstrl(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 )**pthird * sustke(A2D((nn_hls-1))) /   & 
     2272            &                                  MIN( sla(A2D((nn_hls-1)))**( 8.0_wp / 3.0_wp ) + epsln, 0.12_wp ) 
     2273         zsc_vw_1(:,:) = ff_t(A2D((nn_hls-1))) * phml(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1)))**3 * MIN( sla(A2D((nn_hls-1)))**( 8.0_wp / 3.0_wp ), 0.12_wp ) /   & 
     2274            &            ( ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 )**( 2.0_wp / 3.0_wp ) + epsln ) 
    22722275      ELSEWHERE 
    2273          zsc_uw_1(:,:) = sustar(A2D(0))**2 
    2274          zsc_vw_1(:,:) = ff_t(A2D(0)) * phbl(A2D(0)) * sustke(A2D(0))**3 * MIN( sla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) /   & 
    2275             &            ( svstr(A2D(0))**2 + epsln ) 
     2276         zsc_uw_1(:,:) = sustar(A2D((nn_hls-1)))**2 
     2277         zsc_vw_1(:,:) = ff_t(A2D((nn_hls-1))) * phbl(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1)))**3 * MIN( sla(A2D((nn_hls-1)))**( 8.0_wp / 3.0_wp ), 0.12_wp ) /   & 
     2278            &            ( svstr(A2D((nn_hls-1)))**2 + epsln ) 
    22762279      ENDWHERE 
    2277       DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
     2280      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 
    22782281         IF ( l_conv(ji,jj) ) THEN 
    22792282            IF ( jk <= nmld(ji,jj) ) THEN 
     
    22972300      ! (X0.3) and pressure (X0.5)] 
    22982301      ! ---------------------------------------------------------------------- 
    2299       WHERE ( l_conv(A2D(0)) ) 
    2300          zsc_wth_1(:,:) = swbav(A2D(0)) * swth0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * shol(A2D(0)) ) ) * phml(A2D(0)) /   & 
    2301             &             ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln ) 
    2302          zsc_ws_1(:,:)  = swbav(A2D(0)) * sws0(A2D(0))  * ( 1.0_wp + EXP( 0.2_wp * shol(A2D(0)) ) ) * phml(A2D(0)) /   & 
    2303             &             ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln ) 
     2302      WHERE ( l_conv(A2D((nn_hls-1))) ) 
     2303         zsc_wth_1(:,:) = swbav(A2D((nn_hls-1))) * swth0(A2D((nn_hls-1))) * ( 1.0_wp + EXP( 0.2_wp * shol(A2D((nn_hls-1))) ) ) * phml(A2D((nn_hls-1))) /   & 
     2304            &             ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln ) 
     2305         zsc_ws_1(:,:)  = swbav(A2D((nn_hls-1))) * sws0(A2D((nn_hls-1)))  * ( 1.0_wp + EXP( 0.2_wp * shol(A2D((nn_hls-1))) ) ) * phml(A2D((nn_hls-1))) /   & 
     2306            &             ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln ) 
    23042307      ELSEWHERE 
    23052308         zsc_wth_1(:,:) = 0.0_wp 
    23062309         zsc_ws_1(:,:)  = 0.0_wp 
    23072310      ENDWHERE 
    2308       DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
     2311      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 
    23092312         IF ( l_conv(ji,jj) ) THEN 
    23102313            IF ( jk <= nmld(ji,jj) ) THEN 
     
    23272330         END IF 
    23282331      END_3D 
    2329       DO_2D( 0, 0, 0, 0 ) 
     2332      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    23302333         IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) ) THEN 
    23312334            ztau_sc_u(ji,jj)    = phml(ji,jj) / ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird *                             & 
     
    23472350         END IF 
    23482351      END_2D 
    2349       DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 
     2352      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld ) 
    23502353         IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) .AND. ( jk <= nbld(ji,jj) ) ) THEN 
    23512354            zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 
     
    23692372      IF ( ln_dia_osm ) THEN 
    23702373         IF ( iom_use("zwth_ent") ) THEN   ! Upward turb. temperature entrainment flux 
    2371             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zwth_ent(:,:); CALL iom_put( "zwth_ent", osmdia2d ) 
     2374            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwth_ent(:,:); CALL iom_put( "zwth_ent", osmdia2d ) 
    23722375         END IF 
    23732376         IF ( iom_use("zws_ent")  ) THEN   ! Upward turb. salinity entrainment flux 
    2374             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zws_ent(:,:);  CALL iom_put( "zws_ent", osmdia2d ) 
     2377            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zws_ent(:,:);  CALL iom_put( "zws_ent", osmdia2d ) 
    23752378         END IF 
    23762379      END IF 
    23772380      ! 
    23782381      zsc_vw_1(:,:) = 0.0_wp 
    2379       WHERE ( l_conv(A2D(0)) ) 
    2380          zsc_uw_1(:,:) = -1.0_wp * swb0(A2D(0)) * sustar(A2D(0))**2 * phml(A2D(0)) /   & 
    2381             &            ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln ) 
    2382          zsc_uw_2(:,:) =           swb0(A2D(0)) * sustke(A2D(0))    * phml(A2D(0)) /   & 
    2383             &            ( svstr(A2D(0))**3 + 0.5_wp * swstrc(A2D(0))**3 + epsln )**( 2.0_wp / 3.0_wp ) 
     2382      WHERE ( l_conv(A2D((nn_hls-1))) ) 
     2383         zsc_uw_1(:,:) = -1.0_wp * swb0(A2D((nn_hls-1))) * sustar(A2D((nn_hls-1)))**2 * phml(A2D((nn_hls-1))) /   & 
     2384            &            ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln ) 
     2385         zsc_uw_2(:,:) =           swb0(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1)))    * phml(A2D((nn_hls-1))) /   & 
     2386            &            ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln )**( 2.0_wp / 3.0_wp ) 
    23842387      ELSEWHERE 
    23852388         zsc_uw_1(:,:) = 0.0_wp 
    23862389      ENDWHERE 
    2387       DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
     2390      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 
    23882391         IF ( l_conv(ji,jj) ) THEN 
    23892392            IF ( jk <= nmld(ji,jj) ) THEN 
     
    24022405      END_3D 
    24032406      ! 
    2404       DO_2D( 0, 0, 0, 0 ) 
     2407      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    24052408         IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) ) THEN 
    24062409            IF ( n_ddh(ji,jj) == 0 ) THEN 
     
    24212424         END IF 
    24222425      END_2D 
    2423       DO_3D( 0, 0, 0, 0, jkf_mld, jkm_bld )   ! Need ztau_sc_u to be available. Change to array. 
     2426      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jkf_mld, jkm_bld )   ! Need ztau_sc_u to be available. Change to array. 
    24242427         IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) .AND. ( jk >= nmld(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN 
    24252428            zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj) 
     
    24362439         IF ( iom_use("ghamu_0") )    CALL iom_put( "ghamu_0", wmask * ghamu ) 
    24372440         IF ( iom_use("zsc_uw_1_0") ) THEN 
    2438             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zsc_uw_1(:,:); CALL iom_put( "zsc_uw_1_0", osmdia2d ) 
     2441            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_uw_1(:,:); CALL iom_put( "zsc_uw_1_0", osmdia2d ) 
    24392442         END IF 
    24402443      END IF 
     
    24432446      ! (X0.3) ] 
    24442447      ! ----------------------------------------------------------------------- 
    2445       WHERE ( l_conv(A2D(0)) ) 
    2446          zsc_wth_1(:,:) = swth0(A2D(0)) / ( 1.0_wp - 0.56_wp * EXP( shol(A2D(0)) ) ) 
    2447          zsc_ws_1(:,:)  = sws0(A2D(0))  / ( 1.0_wp - 0.56_wp * EXP( shol(A2D(0)) ) ) 
    2448          WHERE ( l_pyc(A2D(0)) )   ! Pycnocline scales 
    2449             zsc_wth_pyc(:,:) = -0.003_wp * swstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * av_dt_ml(A2D(0)) 
    2450             zsc_ws_pyc(:,:)  = -0.003_wp * swstrc(A2D(0)) * ( 1.0_wp - pdh(A2D(0)) / phbl(A2D(0)) ) * av_ds_ml(A2D(0)) 
     2448      WHERE ( l_conv(A2D((nn_hls-1))) ) 
     2449         zsc_wth_1(:,:) = swth0(A2D((nn_hls-1))) / ( 1.0_wp - 0.56_wp * EXP( shol(A2D((nn_hls-1))) ) ) 
     2450         zsc_ws_1(:,:)  = sws0(A2D((nn_hls-1)))  / ( 1.0_wp - 0.56_wp * EXP( shol(A2D((nn_hls-1))) ) ) 
     2451         WHERE ( l_pyc(A2D((nn_hls-1))) )   ! Pycnocline scales 
     2452            zsc_wth_pyc(:,:) = -0.003_wp * swstrc(A2D((nn_hls-1))) * ( 1.0_wp - pdh(A2D((nn_hls-1))) / phbl(A2D((nn_hls-1))) ) * av_dt_ml(A2D((nn_hls-1))) 
     2453            zsc_ws_pyc(:,:)  = -0.003_wp * swstrc(A2D((nn_hls-1))) * ( 1.0_wp - pdh(A2D((nn_hls-1))) / phbl(A2D((nn_hls-1))) ) * av_ds_ml(A2D((nn_hls-1))) 
    24512454         END WHERE 
    24522455      ELSEWHERE 
    2453          zsc_wth_1(:,:) = 2.0_wp * swthav(A2D(0)) 
    2454          zsc_ws_1(:,:)  =          sws0(A2D(0)) 
     2456         zsc_wth_1(:,:) = 2.0_wp * swthav(A2D((nn_hls-1))) 
     2457         zsc_ws_1(:,:)  =          sws0(A2D((nn_hls-1))) 
    24552458      END WHERE 
    2456       DO_3D( 0, 0, 0, 0, 1, MAX( jkm_mld, jkm_bld ) ) 
     2459      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, MAX( jkm_mld, jkm_bld ) ) 
    24572460         IF ( l_conv(ji,jj) ) THEN 
    24582461            IF ( ( jk > 1 ) .AND. ( jk <= nmld(ji,jj) ) ) THEN 
     
    24842487      END_3D 
    24852488      ! 
    2486       WHERE ( l_conv(A2D(0)) ) 
    2487          zsc_uw_1(:,:) = sustar(A2D(0))**2 
    2488          zsc_vw_1(:,:) = ff_t(A2D(0)) * sustke(A2D(0)) * phml(A2D(0)) 
     2489      WHERE ( l_conv(A2D((nn_hls-1))) ) 
     2490         zsc_uw_1(:,:) = sustar(A2D((nn_hls-1)))**2 
     2491         zsc_vw_1(:,:) = ff_t(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1))) * phml(A2D((nn_hls-1))) 
    24892492      ELSEWHERE 
    2490          zsc_uw_1(:,:) = sustar(A2D(0))**2 
     2493         zsc_uw_1(:,:) = sustar(A2D((nn_hls-1)))**2 
    24912494         zsc_uw_2(:,:) = ( 2.25_wp - 3.0_wp * ( 1.0_wp - EXP( -1.25_wp * 2.0_wp ) ) ) * ( 1.0_wp - EXP( -4.0_wp * 2.0_wp ) ) *   & 
    24922495            &            zsc_uw_1(:,:) 
    2493          zsc_vw_1(:,:) = ff_t(A2D(0)) * sustke(A2D(0)) * phbl(A2D(0)) 
     2496         zsc_vw_1(:,:) = ff_t(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1))) * phbl(A2D((nn_hls-1))) 
    24942497         zsc_vw_2(:,:) = -0.11_wp * SIN( 3.14159_wp * ( 2.0_wp + 0.4_wp ) ) * EXP( -1.0_wp * ( 1.5_wp + 2.0_wp )**2 ) *   & 
    24952498            &            zsc_vw_1(:,:) 
    24962499      ENDWHERE 
    2497       DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) ) 
     2500      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) ) 
    24982501         IF ( l_conv(ji,jj) ) THEN 
    24992502            IF ( jk <= nmld(ji,jj) ) THEN 
     
    25302533         IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask * ghamv ) 
    25312534         IF ( iom_use("zsc_uw_1_f") ) THEN 
    2532             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zsc_uw_1(:,:); CALL iom_put( "zsc_uw_1_f", osmdia2d ) 
     2535            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_uw_1(:,:); CALL iom_put( "zsc_uw_1_f", osmdia2d ) 
    25332536         END IF 
    25342537         IF ( iom_use("zsc_vw_1_f") ) THEN 
    2535             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zsc_vw_1(:,:); CALL iom_put( "zsc_vw_1_f", osmdia2d ) 
     2538            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_vw_1(:,:); CALL iom_put( "zsc_vw_1_f", osmdia2d ) 
    25362539         END IF 
    25372540         IF ( iom_use("zsc_uw_2_f") ) THEN 
    2538             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zsc_uw_2(:,:); CALL iom_put( "zsc_uw_2_f", osmdia2d ) 
     2541            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_uw_2(:,:); CALL iom_put( "zsc_uw_2_f", osmdia2d ) 
    25392542         END IF 
    25402543         IF ( iom_use("zsc_vw_2_f") ) THEN 
    2541             osmdia2d(A2D(0)) = tmask(A2D(0),1) * zsc_vw_2(:,:); CALL iom_put( "zsc_vw_2_f", osmdia2d ) 
     2544            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_vw_2(:,:); CALL iom_put( "zsc_vw_2_f", osmdia2d ) 
    25422545         END IF 
    25432546      END IF 
     
    25482551      ! Make surface forced velocity non-gradient terms go to zero at the base 
    25492552      ! of the boundary layer. 
    2550       DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 
     2553      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld ) 
    25512554         IF ( ( .NOT. l_conv(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN 
    25522555            znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / phbl(ji,jj)   ! ALMG to think about 
     
    25692572         z3ddz_pyc_2(:,:,:) = 0.0_wp 
    25702573      END IF 
    2571       DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 
     2574      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld ) 
    25722575         IF ( l_conv (ji,jj) ) THEN 
    25732576            ! Unstable conditions. Shouldn;t be needed with no pycnocline code. 
     
    26372640         IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) ) 
    26382641      END IF 
    2639       DO_3D( 0, 0, 0, 0, 2, jkm_bld ) 
     2642      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld ) 
    26402643         IF ( .NOT. l_conv (ji,jj) ) THEN 
    26412644            IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
     
    26722675      END IF 
    26732676      ! 
    2674       DO_2D( 0, 0, 0, 0 ) 
     2677      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    26752678         ghamt(ji,jj,nbld(ji,jj)) = 0.0_wp 
    26762679         ghams(ji,jj,nbld(ji,jj)) = 0.0_wp 
     
    26832686         IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask * ghamv ) 
    26842687         IF ( iom_use("zviscos") ) THEN 
    2685             osmdia3d(A2D(0),:) = wmask(A2D(0),:) * pviscos; CALL iom_put( "zviscos", osmdia3d ) 
     2688            osmdia3d(A2D((nn_hls-1)),:) = wmask(A2D((nn_hls-1)),:) * pviscos; CALL iom_put( "zviscos", osmdia3d ) 
    26862689         END IF 
    26872690      END IF 
     
    27112714      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdbdx_mle    ! MLE horiz gradients at u points 
    27122715      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdbdy_mle    ! MLE horiz gradients at v points 
    2713       REAL(wp), DIMENSION(A2D(0)),  INTENT(inout) ::   pdbds_mle    ! Magnitude of horizontal buoyancy gradient 
     2716      REAL(wp), DIMENSION(A2D((nn_hls-1))),  INTENT(inout) ::   pdbds_mle    ! Magnitude of horizontal buoyancy gradient 
    27142717      ! 
    27152718      ! Local variables 
     
    27182721      REAL(wp)                         ::   zc 
    27192722      REAL(wp)                         ::   zN2_c        ! Local buoyancy difference from 10m value 
    2720       REAL(wp), DIMENSION(A2D(1))      ::   ztm 
    2721       REAL(wp), DIMENSION(A2D(1))      ::   zsm 
    2722       REAL(wp), DIMENSION(A2D(1),jpts) ::   ztsm_midu 
    2723       REAL(wp), DIMENSION(A2D(1),jpts) ::   ztsm_midv 
    2724       REAL(wp), DIMENSION(A2D(1),jpts) ::   zabu 
    2725       REAL(wp), DIMENSION(A2D(1),jpts) ::   zabv 
    2726       REAL(wp), DIMENSION(A2D(1))      ::   zmld_midu 
    2727       REAL(wp), DIMENSION(A2D(1))      ::   zmld_midv 
     2723      REAL(wp), DIMENSION(A2D(nn_hls))      ::   ztm 
     2724      REAL(wp), DIMENSION(A2D(nn_hls))      ::   zsm 
     2725      REAL(wp), DIMENSION(A2D(nn_hls),jpts) ::   ztsm_midu 
     2726      REAL(wp), DIMENSION(A2D(nn_hls),jpts) ::   ztsm_midv 
     2727      REAL(wp), DIMENSION(A2D(nn_hls),jpts) ::   zabu 
     2728      REAL(wp), DIMENSION(A2D(nn_hls),jpts) ::   zabv 
     2729      REAL(wp), DIMENSION(A2D(nn_hls))      ::   zmld_midu 
     2730      REAL(wp), DIMENSION(A2D(nn_hls))      ::   zmld_midv 
    27282731      ! 
    27292732      ! ==  MLD used for MLE  ==! 
     
    27312734      pmld(:,:)  = 0.0_wp     ! Here hmlp used as a dummy variable, integrating vertically N^2 
    27322735      zN2_c = grav * rn_osm_mle_rho_c * r1_rho0   ! Convert density criteria into N^2 criteria 
    2733       DO_3D( 1, 1, 1, 1, nlb10, jpkm1 ) 
     2736      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, nlb10, jpkm1 ) 
    27342737         ikt = mbkt(ji,jj) 
    27352738         pmld(ji,jj) = pmld(ji,jj) + MAX( rn2b(ji,jj,jk), 0.0_wp ) * e3w(ji,jj,jk,Kmm) 
    27362739         IF( pmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level 
    27372740      END_3D 
    2738       DO_2D( 1, 1, 1, 1 ) 
     2741      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    27392742         mld_prof(ji,jj) = MAX( mld_prof(ji,jj), nbld(ji,jj) )   ! Ensure mld_prof .ge. nbld 
    27402743         pmld(ji,jj)     = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 
     
    27442747      ztm(:,:) = 0.0_wp 
    27452748      zsm(:,:) = 0.0_wp 
    2746       DO_3D( 1, 1, 1, 1, 1, ikmax ) 
     2749      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, ikmax ) 
    27472750         zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, mld_prof(ji,jj) - jk ), 1  ), KIND=wp )   ! zc being 0 outside the ML t-points 
    27482751         ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm) 
     
    27552758      zmld_midu(:,:)   =  0.0_wp 
    27562759      ztsm_midu(:,:,:) = 10.0_wp 
    2757       DO_2D( 1, 0, 0, 0 ) 
     2760      DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    27582761         pdtdx(ji,jj)            = ( ztm(ji+1,jj) - ztm(ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj) 
    27592762         pdsdx(ji,jj)            = ( zsm(ji+1,jj) - zsm(ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj) 
     
    27642767      zmld_midv(:,:)   =  0.0_wp 
    27652768      ztsm_midv(:,:,:) = 10.0_wp 
    2766       DO_2D( 0, 0, 1, 0 ) 
     2769      DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 ) 
    27672770         pdtdy(ji,jj)            = ( ztm(ji,jj+1) - ztm(ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 
    27682771         pdsdy(ji,jj)            = ( zsm(ji,jj+1) - zsm(ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj) 
     
    27732776      CALL eos_rab( ztsm_midu, zmld_midu, zabu, Kmm ) 
    27742777      CALL eos_rab( ztsm_midv, zmld_midv, zabv, Kmm ) 
    2775       DO_2D( 1, 0, 0, 0 ) 
     2778      DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    27762779         pdbdx_mle(ji,jj) = grav * ( pdtdx(ji,jj) * zabu(ji,jj,jp_tem) - pdsdx(ji,jj) * zabu(ji,jj,jp_sal) ) 
    27772780      END_2D 
    2778       DO_2D( 0, 0, 1, 0 ) 
     2781      DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 ) 
    27792782         pdbdy_mle(ji,jj) = grav * ( pdtdy(ji,jj) * zabv(ji,jj,jp_tem) - pdsdy(ji,jj) * zabv(ji,jj,jp_sal) ) 
    27802783      END_2D 
    2781       DO_2D( 0, 0, 0, 0 ) 
     2784      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    27822785         pdbds_mle(ji,jj) = SQRT( 0.5_wp * ( pdbdx_mle(ji,  jj) * pdbdx_mle(ji,  jj) + pdbdy_mle(ji,jj  ) * pdbdy_mle(ji,jj  ) +   & 
    27832786            &                                pdbdx_mle(ji-1,jj) * pdbdx_mle(ji-1,jj) + pdbdy_mle(ji,jj-1) * pdbdy_mle(ji,jj-1) ) ) 
     
    28072810      ! Outputs 
    28082811      INTEGER,                     INTENT(in   ) ::   Kmm         ! Time-level index 
    2809       REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pwb_fk 
    2810       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl        ! BL depth 
    2811       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phmle       ! MLE depth 
    2812       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb_ent     ! Buoyancy entrainment flux 
    2813       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient 
     2812      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   pwb_fk 
     2813      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phbl        ! BL depth 
     2814      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phmle       ! MLE depth 
     2815      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_ent     ! Buoyancy entrainment flux 
     2816      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient 
    28142817      ! 
    28152818      ! Local variables 
    28162819      INTEGER                     ::   ji, jj, jk        ! Dummy loop indices 
    2817       REAL(wp), DIMENSION(A2D(0)) ::   znd_param 
     2820      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   znd_param 
    28182821      REAL(wp)                    ::   zthermal, zbeta 
    28192822      REAL(wp)                    ::   zbuoy 
     
    28232826      REAL(wp)                    ::   zdbdz_mle_int 
    28242827      ! 
    2825       znd_param(A2D(0)) = 0.0_wp 
    2826       ! 
    2827       DO_2D( 0, 0, 0, 0 ) 
     2828      znd_param(A2D((nn_hls-1))) = 0.0_wp 
     2829      ! 
     2830      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    28282831         ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf 
    28292832         pwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * pdbds_mle(ji,jj) * pdbds_mle(ji,jj) 
    28302833      END_2D 
    28312834      ! 
    2832       DO_2D( 0, 0, 0, 0 ) 
     2835      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    28332836         ! 
    28342837         IF ( l_conv(ji,jj) ) THEN 
     
    28582861      ! 
    28592862      ! Diagnosis 
    2860       DO_2D( 0, 0, 0, 0 ) 
     2863      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    28612864         ! 
    28622865         IF ( l_conv(ji,jj) ) THEN 
     
    29202923      INTEGER,                     INTENT(in   ) ::   Kmm         ! Time-level index 
    29212924      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pmld        ! == Estimated FK BLD used for MLE horiz gradients == ! 
    2922       REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   phmle       ! MLE depth 
    2923       REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pvel_mle    ! Velocity scale for dhdt with stable ML and FK 
    2924       REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pdiff_mle   ! Extra MLE vertical diff 
    2925       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient 
    2926       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   phbl        ! BL depth 
    2927       REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pwb0tot     ! Total surface buoyancy flux including insolation 
     2925      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   phmle       ! MLE depth 
     2926      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   pvel_mle    ! Velocity scale for dhdt with stable ML and FK 
     2927      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   pdiff_mle   ! Extra MLE vertical diff 
     2928      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient 
     2929      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phbl        ! BL depth 
     2930      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb0tot     ! Total surface buoyancy flux including insolation 
    29282931      ! 
    29292932      ! Local variables 
     
    29392942      ! 
    29402943      ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE 
    2941       DO_2D( 0, 0, 0, 0 ) 
     2944      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    29422945         IF ( l_conv(ji,jj) ) THEN 
    29432946            ztmp =  r1_ft(ji,jj) * MIN( 111e3_wp, e1u(ji,jj) ) / rn_osm_mle_lf 
     
    29482951      END_2D 
    29492952      ! Timestep mixed layer eddy depth 
    2950       DO_2D( 0, 0, 0, 0 ) 
     2953      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    29512954         IF ( l_mle(ji,jj) ) THEN   ! MLE layer growing 
    29522955            ! Buoyancy gradient at base of MLE layer 
     
    29712974      ! 
    29722975      mld_prof(:,:) = 4 
    2973       DO_3D( 0, 0, 0, 0, 5, jpkm1 ) 
     2976      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 ) 
    29742977         IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN( mbkt(ji,jj), jk ) 
    29752978      END_3D 
    2976       DO_2D( 0, 0, 0, 0 ) 
     2979      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    29772980         phmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm) 
    29782981      END_2D 
     
    31263129         ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case) 
    31273130         z1_t2 = 2e-5_wp 
    3128          DO_2D( 0, 0, 0, 0 ) 
     3131         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    31293132            r1_ft(ji,jj) = MIN( 1.0_wp / ( ABS( ff_t(ji,jj)) + epsln ), ABS( ff_t(ji,jj) ) / z1_t2**2 ) 
    31303133         END_2D 
     
    31673170         etmean(:,:,:) = 0.0_wp 
    31683171         ! 
    3169          DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     3172         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 
    31703173            etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1.0_wp, umask(ji-1,jj,  jk) + umask(ji,jj,jk) +   & 
    31713174               &                                              vmask(ji,  jj-1,jk) + vmask(ji,jj,jk) ) 
     
    31793182         etmean(:,:,:) = 0.0_wp 
    31803183         ! 
    3181          DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     3184         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 
    31823185            etmean(ji,jj,jk) = tmask(ji, jj,jk) / MAX( 1.0_wp, 2.0_wp *   tmask(ji,jj,jk) +                               & 
    31833186               &                                               0.5_wp * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) +     & 
     
    32963299      ! 
    32973300      hbl(:,:)  = 0.0_wp   ! Here hbl used as a dummy variable, integrating vertically N^2 
    3298       DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 
     3301      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
    32993302         ikt = mbkt(ji,jj) 
    33003303         hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0.0_wp ) * e3w(ji,jj,jk,Kmm) 
     
    33023305      END_3D 
    33033306      ! 
    3304       DO_2D( 1, 1, 1, 1 ) 
     3307      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    33053308         iiki = MAX( 4, imld_rst(ji,jj) ) 
    33063309         hbl(ji,jj) = gdepw(ji,jj,iiki,Kmm  )   ! Turbocline depth 
Note: See TracChangeset for help on using the changeset viewer.