Index: NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90
===================================================================
--- NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90	(revision 14858)
+++ NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90	(revision 14859)
@@ -364,5 +364,5 @@
       !
       ! mixed-layer variables
-      INTEGER,  DIMENSION(A2D(0)) ::   jp_ext   ! Offset for external level
+      INTEGER,  DIMENSION(A2D(0)) ::   jk_ext   ! Offset for external level
       !
       REAL(wp), DIMENSION(A2D(0)) ::   zhbl   ! BL depth - grid
@@ -610,12 +610,12 @@
       END_2D
       !
-      ! Averages over well-mixed and boundary layer, note BL averages use jp_ext=2 everywhere
-      jp_ext(:,:) = 1   ! ag 19/03
+      ! Averages over well-mixed and boundary layer, note BL averages use jk_ext=2 everywhere
+      jk_ext(:,:) = 1   ! ag 19/03
       CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D(0)), av_t_bl, av_s_bl,      &
-         &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, av_dt_bl,   &
+         &                           av_b_bl, av_u_bl, av_v_bl, jk_ext, av_dt_bl,   &
          &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl )
-      jp_ext(:,:) = nbld(A2D(0)) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03
+      jk_ext(:,:) = nbld(A2D(0)) - nmld(:,:) + jk_ext(:,:) + 1   ! ag 19/03
       CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,          &
-         &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, av_dt_ml,   &
+         &                           av_b_ml, av_u_ml, av_v_ml, jk_ext, av_dt_ml,   &
          &                           av_ds_ml, av_db_ml, av_du_ml, av_dv_ml )
       ! Velocity components in frame aligned with surface stress
@@ -649,6 +649,6 @@
             &                        zdbds_mle )
          ! Recalculate hmle, zmle, zvel_mle, zdiff_mle & redefine mld_proc to be index for new hmle
-         CALL zdf_osm_mle_parameters( Kmm, mld_prof, zmld, zhmle, zvel_mle,   &
-            &                         zdiff_mle, zdbds_mle, zhbl, zwb0tot )
+         CALL zdf_osm_mle_parameters( Kmm, zmld, zhmle, zvel_mle, zdiff_mle,   &
+            &                         zdbds_mle, zhbl, zwb0tot )
       ELSE    ! ln_osm_mle
          ! FK not selected, Boundary Layer only.
@@ -670,23 +670,23 @@
       !            IF ( ztmp > 6 ) THEN
       !               ! pycnocline well resolved
-      !               jp_ext(ji,jj) = 1
+      !               jk_ext(ji,jj) = 1
       !            ELSE
       !               ! pycnocline poorly resolved
-      !               jp_ext(ji,jj) = 0
+      !               jk_ext(ji,jj) = 0
       !            ENDIF
       !         ELSE
       !            ! Stable conditions
-      !            jp_ext(ji,jj) = 0
+      !            jk_ext(ji,jj) = 0
       !         ENDIF
       !      END_2D
       !
-      ! Recalculate bl averages using jp_ext & ml averages .... note no rotation of u & v here..
-      jp_ext(:,:) = 1   ! ag 19/03
+      ! Recalculate bl averages using jk_ext & ml averages .... note no rotation of u & v here..
+      jk_ext(:,:) = 1   ! ag 19/03
       CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D(0)), av_t_bl, av_s_bl,      &
-         &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, av_dt_bl,   &
+         &                           av_b_bl, av_u_bl, av_v_bl, jk_ext, av_dt_bl,   &
          &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl )
-      jp_ext(:,:) = nbld(A2D(0)) - nmld(:,:) + jp_ext(:,:) + 1   ! ag 19/03
+      jk_ext(:,:) = nbld(A2D(0)) - nmld(:,:) + jk_ext(:,:) + 1   ! ag 19/03
       CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,          &
-         &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, av_dt_ml,   &
+         &                           av_b_ml, av_u_ml, av_v_ml, jk_ext, av_dt_ml,   &
          &                           av_ds_ml, av_db_ml, av_du_ml, av_dv_ml )   ! ag 19/03
       !
@@ -725,7 +725,7 @@
       !
       ! Recalculate BL averages and differences using new BL depth
-      jp_ext(:,:) = 1   ! ag 19/03
+      jk_ext(:,:) = 1   ! ag 19/03
       CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D(0)), av_t_bl, av_s_bl,      &
-         &                           av_b_bl, av_u_bl, av_v_bl, jp_ext, av_dt_bl,   &
+         &                           av_b_bl, av_u_bl, av_v_bl, jk_ext, av_dt_bl,   &
          &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl )
       !
@@ -752,9 +752,9 @@
       !
       ! Average over the depth of the mixed layer in the convective boundary layer
-      !      jp_ext = nbld - nmld + 1
+      !      jk_ext = nbld - nmld + 1
       ! Recalculate ML averages and differences using new ML depth
-      jp_ext(:,:) = nbld(A2D(0)) - nmld(A2D(0)) + jp_ext(:,:) + 1   ! ag 19/03
+      jk_ext(:,:) = nbld(A2D(0)) - nmld(A2D(0)) + jk_ext(:,:) + 1   ! ag 19/03
       CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, av_t_ml, av_s_ml,    &
-         &                           av_b_ml, av_u_ml, av_v_ml, jp_ext, av_dt_ml,   &
+         &                           av_b_ml, av_u_ml, av_v_ml, jk_ext, av_dt_ml,   &
          &                           av_ds_ml, av_db_ml, av_du_ml, av_dv_ml )
       !
@@ -774,6 +774,6 @@
       ! Calculate non-gradient components of the flux-gradient relationships
       ! --------------------------------------------------------------------
-      jp_ext(:,:) = 1   ! ag 19/03
-      CALL zdf_osm_fgr_terms( Kmm, jp_ext, zhbl, zhml, zdh,                              &
+      jk_ext(:,:) = 1   ! ag 19/03
+      CALL zdf_osm_fgr_terms( Kmm, jk_ext, zhbl, zhml, zdh,                              &
          &                    zdhdt, zshear, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext,   &
          &                    zdiffut, zviscos )
@@ -920,5 +920,5 @@
             IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) )                      ! U_10
             IF ( iom_use("wind_wave_abs_power") ) THEN
-               osmdia2d(A2D(0)) = 1000.0_wp * rho0 * tmask(:,:,1) * sustar**2 * SQRT( ut0sd**2 + vt0sd**2 )
+               osmdia2d(A2D(0)) = 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar**2 * SQRT( ut0sd(A2D(0))**2 + vt0sd(A2D(0))**2 )
                CALL iom_put( "wind_wave_abs_power", osmdia2d )
             END IF
@@ -1007,5 +1007,5 @@
          END IF
          IF ( iom_use("jp_ext") ) THEN                                                     ! =1 if pycnocline resolved internal to zdf_osm routine
-            osmdia2d(A2D(0)) = tmask(A2D(0),1) * jp_ext;    CALL iom_put( "jp_ext",    osmdia2d )
+            osmdia2d(A2D(0)) = tmask(A2D(0),1) * jk_ext;    CALL iom_put( "jp_ext",    osmdia2d )
          END IF
          IF ( iom_use("j_ddh") ) THEN                                                      ! Index forpyc thicknessh internal to zdf_osm routine
@@ -1277,5 +1277,5 @@
       ! Determins stability and set flag l_conv
       DO_2D( 0, 0, 0, 0 )
-         IF ( shol(ji,jj) < 0._wp ) THEN
+         IF ( shol(ji,jj) < 0.0_wp ) THEN
             l_conv(ji,jj) = .TRUE.
          ELSE
@@ -1659,5 +1659,5 @@
                   zdb = MAX(  grav * ( zthermal * ( av_t_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) -               &
                      &                 zbeta    * ( av_s_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0_wp ) +   &
-                     &  2.0 * svstr(ji,jj)**2 / zhbl_s
+                     &  2.0_wp * svstr(ji,jj)**2 / zhbl_s
                   !
                   ! Alan is thuis right? I have simply changed hbli to hbl
@@ -1956,5 +1956,7 @@
       !
       IF ( ln_dia_pyc_scl ) THEN   ! Output of pycnocline gradient profiles
-         IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask(:,:,:) * pdbdz(:,:,:) )
+         IF ( iom_use("zdbdz_pyc") ) THEN
+            osmdia3d(A2D(0),:) = wmask(A2D(0),:) * pdbdz(:,:,:); CALL iom_put( "zdbdz_pyc", osmdia3d )
+         END IF
       END IF
       !
@@ -2154,5 +2156,7 @@
          !
       END_2D
-      IF( iom_use("pb_coup") ) CALL iom_put( "pb_coup", tmask(:,:,1) * zb_coup(:,:) )   ! BBL-coupling velocity scale
+      IF( iom_use("pb_coup") ) THEN
+         osmdia2d(A2D(0)) = tmask(A2D(0),1) * zb_coup(:,:); CALL iom_put( "pb_coup", osmdia2d )   ! BBL-coupling velocity scale
+      END IF
       !
    END SUBROUTINE zdf_osm_diffusivity_viscosity
@@ -2259,6 +2263,6 @@
       !
       IF ( ln_dia_osm ) THEN
-         IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu )
-         IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv )
+         IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask * ghamu )
+         IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask * ghamv )
       END IF
       !
@@ -2371,6 +2375,10 @@
       !
       IF ( ln_dia_osm ) THEN
-         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! Upward turb. temperature entrainment flux
-         IF ( iom_use("zws_ent")  ) CALL iom_put( "zws_ent",  tmask(:,:,1)*zws_ent  )   ! Upward turb. salinity    entrainment flux
+         IF ( iom_use("zwth_ent") ) THEN   ! Upward turb. temperature entrainment flux
+            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zwth_ent(:,:); CALL iom_put( "zwth_ent", osmdia2d )
+         END IF
+         IF ( iom_use("zws_ent")  ) THEN   ! Upward turb. salinity entrainment flux
+            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zws_ent(:,:);  CALL iom_put( "zws_ent", osmdia2d )
+         END IF
       END IF
       !
@@ -2413,5 +2421,5 @@
                zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * av_dv_ml(ji,jj)
             ENDIF
-            zb_cubic(ji,jj) = pdh(ji,jj) / phbl(ji,jj) * suw0(ji,jj) - ( 2.0 + pdh(ji,jj) / phml(ji,jj) ) * zuw_bse(ji,jj)
+            zb_cubic(ji,jj) = pdh(ji,jj) / phbl(ji,jj) * suw0(ji,jj) - ( 2.0_wp + pdh(ji,jj) / phml(ji,jj) ) * zuw_bse(ji,jj)
             za_cubic(ji,jj) = zuw_bse(ji,jj) - zb_cubic(ji,jj)
             zvw_max = 0.7_wp * ff_t(ji,jj) * ( sustke(ji,jj) * dstokes(ji,jj) + 0.7_wp * sustar(ji,jj) * phml(ji,jj) )
@@ -2433,6 +2441,8 @@
       !
       IF ( ln_dia_osm ) THEN
-         IF ( iom_use("ghamu_0") )    CALL iom_put( "ghamu_0",    wmask*ghamu           )
-         IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 )
+         IF ( iom_use("ghamu_0") )    CALL iom_put( "ghamu_0", wmask * ghamu )
+         IF ( iom_use("zsc_uw_1_0") ) THEN
+            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zsc_uw_1(:,:); CALL iom_put( "zsc_uw_1_0", osmdia2d )
+         END IF
       END IF
       !
@@ -2448,6 +2458,6 @@
          END WHERE
       ELSEWHERE
-         zsc_wth_1(:,:) = 2.0 * swthav(A2D(0))
-         zsc_ws_1(:,:)  =       sws0(A2D(0))
+         zsc_wth_1(:,:) = 2.0_wp * swthav(A2D(0))
+         zsc_ws_1(:,:)  =          sws0(A2D(0))
       END WHERE
       DO_3D( 0, 0, 0, 0, 1, MAX( jkm_mld, jkm_bld ) )
@@ -2508,5 +2518,5 @@
                znd    = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj)
                zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
-               IF ( zznd_d <= 2.0 ) THEN
+               IF ( zznd_d <= 2.0_wp ) THEN
                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5_wp * 0.3_wp *                                              &
                      &                                ( 2.25_wp - 3.0_wp * ( 1.0_wp - EXP( -1.25_wp * zznd_d ) ) *   &
@@ -2524,10 +2534,18 @@
       !
       IF ( ln_dia_osm ) THEN
-         IF ( iom_use("ghamu_f") )    CALL iom_put( "ghamu_f",    wmask       *ghamu    )
-         IF ( iom_use("ghamv_f") )    CALL iom_put( "ghamv_f",    wmask       *ghamv    )
-         IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 )
-         IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 )
-         IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 )
-         IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 )
+         IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask * ghamu )
+         IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask * ghamv )
+         IF ( iom_use("zsc_uw_1_f") ) THEN
+            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zsc_uw_1(:,:); CALL iom_put( "zsc_uw_1_f", osmdia2d )
+         END IF
+         IF ( iom_use("zsc_vw_1_f") ) THEN
+            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zsc_vw_1(:,:); CALL iom_put( "zsc_vw_1_f", osmdia2d )
+         END IF
+         IF ( iom_use("zsc_uw_2_f") ) THEN
+            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zsc_uw_2(:,:); CALL iom_put( "zsc_uw_2_f", osmdia2d )
+         END IF
+         IF ( iom_use("zsc_vw_2_f") ) THEN
+            osmdia2d(A2D(0)) = tmask(A2D(0),1) * zsc_vw_2(:,:); CALL iom_put( "zsc_vw_2_f", osmdia2d )
+         END IF
       END IF
       !
@@ -2654,6 +2672,6 @@
       END IF
       IF ( ln_dia_osm ) THEN
-         IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu )
-         IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv )
+         IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask * ghamu )
+         IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask * ghamv )
       END IF
       IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN   ! Deallocate arrays used for output of pycnocline gradient/shear profiles
@@ -2669,6 +2687,6 @@
       !
       IF ( ln_dia_osm ) THEN
-         IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu )
-         IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv )
+         IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask * ghamu )
+         IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask * ghamv )
          IF ( iom_use("zviscos") ) THEN
             osmdia3d(A2D(0),:) = wmask(A2D(0),:) * pviscos; CALL iom_put( "zviscos", osmdia3d )
@@ -2892,6 +2910,6 @@
    END SUBROUTINE zdf_osm_osbl_state_fk
 
-   SUBROUTINE zdf_osm_mle_parameters( Kmm, kmld_prof, pmld, phmle, pvel_mle,   &
-      &                               pdiff_mle, pdbds_mle, phbl, pwb0tot )
+   SUBROUTINE zdf_osm_mle_parameters( Kmm, pmld, phmle, pvel_mle, pdiff_mle,   &
+      &                               pdbds_mle, phbl, pwb0tot )
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE zdf_osm_mle_parameters  ***
@@ -2908,5 +2926,4 @@
       !!----------------------------------------------------------------------
       INTEGER,                     INTENT(in   ) ::   Kmm         ! Time-level index
-      INTEGER,  DIMENSION(:,:),    INTENT(inout) ::   kmld_prof
       REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pmld        ! == Estimated FK BLD used for MLE horiz gradients == !
       REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   phmle       ! MLE depth
@@ -2943,6 +2960,6 @@
             zthermal = rab_n(ji,jj,1,jp_tem)
             zbeta    = rab_n(ji,jj,1,jp_sal)
-            zbuoy = grav * ( zthermal * ts(ji,jj,kmld_prof(ji,jj)+2,jp_tem,Kmm) -   &
-               &             zbeta    * ts(ji,jj,kmld_prof(ji,jj)+2,jp_sal,Kmm) )
+            zbuoy = grav * ( zthermal * ts(ji,jj,mld_prof(ji,jj)+2,jp_tem,Kmm) -   &
+               &             zbeta    * ts(ji,jj,mld_prof(ji,jj)+2,jp_sal,Kmm) )
             zdb_mle = av_b_bl(ji,jj) - zbuoy
             ! Timestep hmle
@@ -2960,10 +2977,10 @@
       END_2D
       !
-      kmld_prof(:,:) = 4
+      mld_prof(:,:) = 4
       DO_3D( 0, 0, 0, 0, 5, jpkm1 )
-         IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) kmld_prof(ji,jj) = MIN( mbkt(ji,jj), jk )
+         IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN( mbkt(ji,jj), jk )
       END_3D
       DO_2D( 0, 0, 0, 0 )
-         phmle(ji,jj) = gdepw(ji,jj,kmld_prof(ji,jj),Kmm)
+         phmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm)
       END_2D
       !
@@ -3194,5 +3211,8 @@
       ghamv(:,:,:) = 0.0_wp
       !
-      IF ( ln_dia_osm ) osmdia2d(:,:) = 0.0_wp   ! Initialise auxiliary array for diagnostic output
+      IF ( ln_dia_osm ) THEN   ! Initialise auxiliary arrays for diagnostic output
+         osmdia2d(:,:)   = 0.0_wp
+         osmdia3d(:,:,:) = 0.0_wp
+      END IF
       !
    END SUBROUTINE zdf_osm_init
@@ -3303,5 +3323,5 @@
       END IF
       !
-      ww(:,:,:) = 0._wp
+      ww(:,:,:) = 0.0_wp
       WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
       !
