Index: NEMO/branches/2021/dev_r15388_updated_zdfiwm/src/OCE/ZDF/zdfiwm.F90
===================================================================
--- NEMO/branches/2021/dev_r15388_updated_zdfiwm/src/OCE/ZDF/zdfiwm.F90	(revision 15514)
+++ NEMO/branches/2021/dev_r15388_updated_zdfiwm/src/OCE/ZDF/zdfiwm.F90	(revision 15517)
@@ -42,6 +42,6 @@
    LOGICAL ::  ln_tsdiff   ! account for differential T/S wave-driven mixing (=T) or not (=F)
 
-   REAL(wp)::  r1_6 = 1._wp / 6._wp
-   REAL(wp)::  rnu  = 1.4e-6_wp   ! molecular kinematic viscosity
+   REAL(wp)::  z1_6 = 1._wp / 6._wp
+   REAL(wp)::  znu  = 1.4e-6_wp   ! molecular kinematic viscosity
 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ebot_iwm   ! bottom-intensified dissipation above abyssal hills (W/m2)
@@ -135,5 +135,5 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp), SAVE :: zztmp
-      REAL(wp)       :: ztmp1, ztmp2        ! scalar workspace
+      !
       REAL(wp), DIMENSION(A2D(nn_hls))     ::   zfact       ! Used for vertical structure
       REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zReb        ! Turbulence intensity parameter
@@ -145,14 +145,8 @@
       !!----------------------------------------------------------------------
       !
-      !                       !* Set to zero the 1st and last vertical levels of appropriate variables
-      IF( iom_use("emix_iwm") ) THEN
-         zemx_iwm(:,:,:) = 0._wp
-      ENDIF
-      IF( iom_use("av_ratio") ) THEN
-         zav_ratio(:,:,:) = 0._wp
-      ENDIF
-      IF( iom_use("av_wave") .OR. sn_cfctl%l_prtctl ) THEN
-         zav_wave(:,:,:) = 0._wp
-      ENDIF
+      !                       !* Initialize appropriately certain variables
+      zav_ratio(:,:,:) = 1._wp * wmask(:,:,:)  ! important to set it to 1 here 
+      IF( iom_use("emix_iwm") )                         zemx_iwm (:,:,:) = 0._wp
+      IF( iom_use("av_wave") .OR. sn_cfctl%l_prtctl )   zav_wave (:,:,:) = 0._wp
       !
       !                       ! ----------------------------- !
@@ -163,32 +157,28 @@
       !                       !* ocean depth using an exponential decay from the seafloor.
       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                ! part independent of the level
-         zfact(ji,jj) = rho0 * (  1._wp - EXP( -ht(ji,jj) * hcri_iwm(ji,jj) )  )
-         IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj)
+         IF( ht(ji,jj) /= 0._wp ) THEN ; zfact(ji,jj) = ecri_iwm(ji,jj) * r1_rho0 / ( 1._wp - EXP( -ht(ji,jj) * hcri_iwm(ji,jj) ) )
+         ELSE                          ; zfact(ji,jj) = 0._wp
+         ENDIF
       END_2D
 
       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )      ! complete with the level-dependent part
-         IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization
-            zemx_iwm(ji,jj,jk) = 0._wp
-         ELSE
-            zemx_iwm(ji,jj,jk) = zfact(ji,jj) * (  EXP( ( gdept(ji,jj,jk  ,Kmm) - ht(ji,jj) ) * hcri_iwm(ji,jj) )     &
-                 &                               - EXP( ( gdept(ji,jj,jk-1,Kmm) - ht(ji,jj) ) * hcri_iwm(ji,jj) ) )   &
-                 &                            / e3w(ji,jj,jk,Kmm)
-         ENDIF
+         zemx_iwm(ji,jj,jk) = zfact(ji,jj) * (  EXP( ( gdept(ji,jj,jk  ,Kmm) - ht(ji,jj) ) * hcri_iwm(ji,jj) )   &
+            &                                 - EXP( ( gdept(ji,jj,jk-1,Kmm) - ht(ji,jj) ) * hcri_iwm(ji,jj) )   &
+            &                                ) * wmask(ji,jj,jk) / e3w(ji,jj,jk,Kmm)
       END_3D
 
                                !* 'bot' component: distribute energy over the time-varying
                                !* ocean depth using an algebraic decay above the seafloor.
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
-         zfact(ji,jj) = 0._wp
-      END_2D
       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )               ! part independent of the level
-         IF( ht(ji,jj) /= 0._wp )  zfact(ji,jj) = ebot_iwm(ji,jj) * (  1._wp +  hbot_iwm(ji,jj) / ht(ji,jj)  ) / rho0
+         IF( ht(ji,jj) /= 0._wp ) THEN ; zfact(ji,jj) = ebot_iwm(ji,jj) * ( 1._wp +  hbot_iwm(ji,jj) / ht(ji,jj)  ) * r1_rho0
+         ELSE                          ; zfact(ji,jj) = 0._wp
+         ENDIF
       END_2D
 
       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )     ! complete with the level-dependent part
-         zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) +                                                                                          & 
-            &         zfact(ji,jj) * (   1._wp / ( 1._wp + ( ht(ji,jj) - gdept(ji,jj,jk  ,Kmm) ) / hbot_iwm(ji,jj) )                        &
-            &                         -  1._wp / ( 1._wp + ( ht(ji,jj) - gdept(ji,jj,jk-1,Kmm) ) / hbot_iwm(ji,jj) )   ) * wmask(ji,jj,jk)  &
-            &                      / e3w(ji,jj,jk,Kmm)
+         zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) +                                                                           & 
+            &                 zfact(ji,jj) * (  1._wp / ( 1._wp + ( ht(ji,jj) - gdept(ji,jj,jk  ,Kmm) ) / hbot_iwm(ji,jj) )  &
+            &                                 - 1._wp / ( 1._wp + ( ht(ji,jj) - gdept(ji,jj,jk-1,Kmm) ) / hbot_iwm(ji,jj) )  &
+            &                                ) * wmask(ji,jj,jk) / e3w(ji,jj,jk,Kmm)
       END_3D
 
@@ -203,5 +193,5 @@
       !
       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
-         IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ensq_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
+         IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ensq_iwm(ji,jj) * r1_rho0 / zfact(ji,jj)
       END_2D
       !
@@ -216,23 +206,23 @@
       END_2D
       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )     ! part independent of the level
-         zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  )
+         zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) )
       END_3D
       !
       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
-         IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = esho_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
+         IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = esho_iwm(ji,jj) * r1_rho0 / zfact(ji,jj)
       END_2D
       !
       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )     ! complete with the level-dependent part
-         zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  )
+         zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) )
       END_3D
 
       ! Calculate turbulence intensity parameter Reb
       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
-         zReb(ji,jj,jk) = zemx_iwm(ji,jj,jk) / MAX( 1.e-20_wp, rnu * rn2(ji,jj,jk) )
+         zReb(ji,jj,jk) = zemx_iwm(ji,jj,jk) / MAX( 1.e-20_wp, znu * rn2(ji,jj,jk) )
       END_3D
       !
       ! Define internal wave-induced diffusivity
       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
-         zav_wave(ji,jj,jk) = zReb(ji,jj,jk) * r1_6 * rnu  ! This corresponds to a constant mixing efficiency of 1/6
+         zav_wave(ji,jj,jk) = zReb(ji,jj,jk) * z1_6 * znu  ! This corresponds to a constant mixing efficiency of 1/6
       END_3D
       !
@@ -240,7 +230,7 @@
          DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224) regimes
             IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
-               zav_wave(ji,jj,jk) = 3.6515_wp * rnu * SQRT( zReb(ji,jj,jk) )
+               zav_wave(ji,jj,jk) = 3.6515_wp * znu * SQRT( zReb(ji,jj,jk) )
             ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN
-               zav_wave(ji,jj,jk) = 0.052125_wp * rnu * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
+               zav_wave(ji,jj,jk) = 0.052125_wp * znu * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
             ENDIF
          END_3D
@@ -248,6 +238,42 @@
       !
       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )    ! Bound diffusivity by molecular value and 100 cm2/s
-         zav_wave(ji,jj,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(ji,jj,jk) ), 1.e-2_wp  ) * wmask(ji,jj,jk)
-      END_3D
+         zav_wave(ji,jj,jk) = MIN( MAX( 1.4e-7_wp, zav_wave(ji,jj,jk) ), 1.e-2_wp ) * wmask(ji,jj,jk)
+      END_3D
+      !
+      !                          ! ----------------------- !
+      !                          !   Update  mixing coefs  !                          
+      !                          ! ----------------------- !
+      !
+      IF( ln_tsdiff ) THEN                !* Option for differential mixing of salinity and temperature
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! Calculate S/T diffusivity ratio as a function of Reb (else it is set to 1)
+            zav_ratio(ji,jj,jk) = ( 0.505_wp + &
+               &                    0.495_wp * TANH( 0.92_wp * ( LOG10( MAX( 1.e-20, zReb(ji,jj,jk) * 5._wp * z1_6 ) ) - 0.60_wp ) ) &
+               &                  ) * wmask(ji,jj,jk)
+         END_3D
+      ENDIF
+      CALL iom_put( "av_ratio", zav_ratio )
+      !
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* update momentum & tracer diffusivity with wave-driven mixing
+         p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) * zav_ratio(ji,jj,jk)
+         p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk)
+         p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk)
+      END_3D
+      !                             !* output internal wave-driven mixing coefficient
+      CALL iom_put( "av_wave", zav_wave )
+                                    !* output useful diagnostics: Kz*N^2 , 
+                                    !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm)
+      IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN
+         ALLOCATE( z2d(A2D(nn_hls)) , z3d(A2D(nn_hls),jpk) )
+         z2d(:,:) = 0._wp ; z3d(:,:,:) = 0._wp   ! Initialisation for iom_put
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+            z3d(ji,jj,jk) = MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk)
+            z2d(ji,jj) = z2d(ji,jj) + rho0 * e3w(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * wmask(ji,jj,jk)
+         END_3D
+         CALL iom_put(  "bflx_iwm", z3d )
+         CALL iom_put( "pcmap_iwm", z2d )
+         DEALLOCATE( z2d , z3d )
+      ENDIF
+      CALL iom_put( "emix_iwm", zemx_iwm )
+
       !
       IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave
@@ -260,5 +286,5 @@
          IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
             CALL mpp_sum( 'zdfiwm', zztmp )
-            zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing
+            zztmp = rho0 * zztmp ! Global integral of rho0 * Kz * N^2 = power contributing to mixing
             !
             IF(lwp) THEN
@@ -272,54 +298,4 @@
       ENDIF
 
-      !                          ! ----------------------- !
-      !                          !   Update  mixing coefs  !                          
-      !                          ! ----------------------- !
-      !
-      IF( ln_tsdiff ) THEN                !* Option for differential mixing of salinity and temperature
-         ztmp1 = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10( 1.e-20_wp ) - 0.60_wp ) )
-         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! Calculate S/T diffusivity ratio as a function of Reb
-            ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6
-            IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN
-               zav_ratio(ji,jj,jk) = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10(ztmp2) - 0.60_wp ) )
-            ELSE
-               zav_ratio(ji,jj,jk) = ztmp1 * wmask(ji,jj,jk)
-            ENDIF
-         END_3D
-         CALL iom_put( "av_ratio", zav_ratio )
-         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* update momentum & tracer diffusivity with wave-driven mixing
-            p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) * zav_ratio(ji,jj,jk)
-            p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk)
-            p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk)
-         END_3D
-         !
-      ELSE                                !* update momentum & tracer diffusivity with wave-driven mixing
-         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
-            p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk)
-            p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk)
-            p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk)
-         END_3D
-      ENDIF
-      !                             !* output internal wave-driven mixing coefficient
-      CALL iom_put( "av_wave", zav_wave )
-                                    !* output useful diagnostics: Kz*N^2 , 
-                                    !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm)
-      IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN
-         ALLOCATE( z2d(A2D(nn_hls)) , z3d(A2D(nn_hls),jpk) )
-         ! Initialisation for iom_put
-         z2d(:,:) = 0._wp ; z3d(:,:,:) = 0._wp
-
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-            z3d(ji,jj,jk) = MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk)
-            z2d(ji,jj) = z2d(ji,jj) + e3w(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * wmask(ji,jj,jk)
-         END_3D
-         DO_2D( 0, 0, 0, 0 )
-            z2d(ji,jj) = rho0 * z2d(ji,jj)
-         END_2D
-         CALL iom_put(  "bflx_iwm", z3d )
-         CALL iom_put( "pcmap_iwm", z2d )
-         DEALLOCATE( z2d , z3d )
-      ENDIF
-      CALL iom_put( "emix_iwm", zemx_iwm )
-
       IF(sn_cfctl%l_prtctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', kdim=jpk)
       !
@@ -356,5 +332,4 @@
       INTEGER  ::   inum               ! local integer
       INTEGER  ::   ios
-      REAL(wp) ::   zbot, zcri, znsq, zsho   ! local scalars
       !
       CHARACTER(len=256)            ::   cn_dir                 ! Root directory for location of ssr files
@@ -372,4 +347,7 @@
       TYPE(FLD  ), DIMENSION(jpiwm) ::   sf_iwm                         ! structure of input fields (file informations, fields read)
       !
+      REAL(wp), DIMENSION(jpi,jpj,4) ::   ztmp
+      REAL(wp), DIMENSION(4)         ::   zdia
+      !
       NAMELIST/namzdf_iwm/ ln_mevar, ln_tsdiff, &
           &                cn_dir, sn_mpb, sn_mpc, sn_mpn, sn_mps, sn_dsb, sn_dsc
@@ -395,7 +373,7 @@
       ! ensures that avt remains larger than its molecular value (=1.4e-7). Therefore, avtb should 
       ! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6).
-      avmb(:) = rnu              ! molecular value
+      avmb(:) = znu              ! molecular value
       avtb(:) = 1.e-10_wp        ! very small diffusive minimum (background avt is specified in zdf_iwm)    
-      avtb_2d(:,:) = 1._wp     ! uniform 
+      avtb_2d(:,:) = 1._wp       ! uniform 
       IF(lwp) THEN                  ! Control print
          WRITE(numout,*)
@@ -422,5 +400,5 @@
       sf_iwm(jp_mpb)%fnow(:,:,1) = 1.e-10_wp
       sf_iwm(jp_mpc)%fnow(:,:,1) = 1.e-10_wp
-      sf_iwm(jp_mpn)%fnow(:,:,1) = 1.e-6_wp
+      sf_iwm(jp_mpn)%fnow(:,:,1) = 1.e-5_wp
       sf_iwm(jp_mps)%fnow(:,:,1) = 1.e-10_wp
       sf_iwm(jp_dsb)%fnow(:,:,1) = 100._wp
@@ -439,14 +417,17 @@
       hcri_iwm(:,:) = 1._wp / hcri_iwm(:,:) ! only the inverse height is used, hence calculated here once for all
 
-      zbot = glob_sum( 'zdfiwm', e1e2t(:,:) * ebot_iwm(:,:) )
-      zcri = glob_sum( 'zdfiwm', e1e2t(:,:) * ecri_iwm(:,:) )
-      znsq = glob_sum( 'zdfiwm', e1e2t(:,:) * ensq_iwm(:,:) )
-      zsho = glob_sum( 'zdfiwm', e1e2t(:,:) * esho_iwm(:,:) )
+      ! diags
+      ztmp(:,:,1) = e1e2t(:,:) * ebot_iwm(:,:)
+      ztmp(:,:,2) = e1e2t(:,:) * ecri_iwm(:,:)
+      ztmp(:,:,3) = e1e2t(:,:) * ensq_iwm(:,:)
+      ztmp(:,:,4) = e1e2t(:,:) * esho_iwm(:,:)
+
+      zdia(1:4) = glob_sum_vec( 'zdfiwm', ztmp(:,:,1:4) )
 
       IF(lwp) THEN
-         WRITE(numout,*) '      Dissipation above abyssal hills:        ', zbot * 1.e-12_wp, 'TW'
-         WRITE(numout,*) '      Dissipation along topographic slopes:   ', zcri * 1.e-12_wp, 'TW'
-         WRITE(numout,*) '      Dissipation scaling with N^2:           ', znsq * 1.e-12_wp, 'TW'
-         WRITE(numout,*) '      Dissipation due to shoaling:            ', zsho * 1.e-12_wp, 'TW'
+         WRITE(numout,*) '      Dissipation above abyssal hills:        ', zdia(1) * 1.e-12_wp, 'TW'
+         WRITE(numout,*) '      Dissipation along topographic slopes:   ', zdia(2) * 1.e-12_wp, 'TW'
+         WRITE(numout,*) '      Dissipation scaling with N^2:           ', zdia(3) * 1.e-12_wp, 'TW'
+         WRITE(numout,*) '      Dissipation due to shoaling:            ', zdia(4) * 1.e-12_wp, 'TW'
       ENDIF
       !
