Index: /branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- /branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 15796)
+++ /branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 15797)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                       ***  MODULE  zdfgls  ***
-   !! Ocean physics:  vertical mixing coefficient computed from the gls 
+   !! Ocean physics:  vertical mixing coefficient computed from the gls
    !!                 turbulent closure parameterization
    !!======================================================================
@@ -16,5 +16,5 @@
    !!   gls_rst       : read/write gls restart in ocean restart file
    !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and active tracers 
+   USE oce            ! ocean dynamics and active tracers
    USE dom_oce        ! ocean space and time domain
    USE domvvl         ! ocean space and time domain : variable volume layer
@@ -31,5 +31,5 @@
    USE iom            ! I/O manager library
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
 
    IMPLICIT NONE
@@ -61,9 +61,9 @@
    REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing
    REAL(wp) ::   rn_hsro           ! Minimum surface roughness
-   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 
+   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1)
 
    REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters
-   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1 
-   REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn    
+   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1
+   REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn
    REAL(wp) ::   rghmin        = -0.28_wp
    REAL(wp) ::   rgh0          =  0.0329_wp
@@ -72,6 +72,6 @@
    REAL(wp) ::   ra2           =  0.74_wp
    REAL(wp) ::   rb1           = 16.60_wp
-   REAL(wp) ::   rb2           = 10.10_wp         
-   REAL(wp) ::   re2           =  1.33_wp         
+   REAL(wp) ::   rb2           = 10.10_wp
+   REAL(wp) ::   re2           =  1.33_wp
    REAL(wp) ::   rl1           =  0.107_wp
    REAL(wp) ::   rl2           =  0.0032_wp
@@ -133,5 +133,5 @@
       INTEGER  ::   ji, jj, jk, ibot, ibotm1, dir  ! dummy loop arguments
       REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1, zex2   ! local scalars
-      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      - 
+      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      -
       REAL(wp) ::   zratio, zrn2, zflxb, sh             !   -      -
       REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      -
@@ -139,5 +139,5 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zdep
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zkar
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves 
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhsro       ! Surface roughness (surface waves)
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   eb          ! tke at time before
@@ -156,5 +156,5 @@
       CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro )
       CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi  )
-      
+
       ! Preliminary computing
 
@@ -165,30 +165,30 @@
          avm (:,:,:) = avm_k (:,:,:)
          avmu(:,:,:) = avmu_k(:,:,:)
-         avmv(:,:,:) = avmv_k(:,:,:) 
+         avmv(:,:,:) = avmv_k(:,:,:)
       ENDIF
 
       ! Compute surface and bottom friction at T-points
-!CDIR NOVERRCHK          
-      DO jj = 2, jpjm1          
-!CDIR NOVERRCHK         
+!CDIR NOVERRCHK
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
          DO ji = fs_2, fs_jpim1   ! vector opt.         
             !
             ! surface friction
             ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)
-            !   
-            ! bottom friction (explicit before friction)        
-            ! Note that we chose here not to bound the friction as in dynbfr)   
-            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &         
-               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )      
-            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &         
-               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )      
-            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)         
-         END DO         
-      END DO    
+            !
+            ! bottom friction (explicit before friction)
+            ! Note that we chose here not to bound the friction as in dynbfr)
+            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &
+               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )
+            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &
+               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )
+            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)
+         END DO
+      END DO
 
       ! Set surface roughness length
       SELECT CASE ( nn_z0_met )
       !
-      CASE ( 0 )             ! Constant roughness          
+      CASE ( 0 )             ! Constant roughness
          zhsro(:,:) = rn_hsro
       CASE ( 1 )             ! Standard Charnock formula
@@ -226,5 +226,5 @@
       IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1 
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zup   = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1)
@@ -275,5 +275,5 @@
                IF( ln_sigpsi ) THEN
                   zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) )     ! 0. <= zsigpsi <= 1.
-                  zwall_psi(ji,jj,jk) = rsc_psi /   & 
+                  zwall_psi(ji,jj,jk) = rsc_psi /   &
                      &     (  zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp )  )
                ELSE
@@ -294,5 +294,5 @@
                ! diagonal
                z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  &
-                  &                       + rdt * zdiss * tmask(ji,jj,jk) 
+                  &                       + rdt * zdiss * tmask(ji,jj,jk)
                !
                ! right hand side in en
@@ -316,14 +316,14 @@
       ! First level
       en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
-      en(:,:,1) = MAX(en(:,:,1), rn_emin) 
+      en(:,:,1) = MAX(en(:,:,1), rn_emin)
       z_elem_a(:,:,1) = en(:,:,1)
       z_elem_c(:,:,1) = 0._wp
       z_elem_b(:,:,1) = 1._wp
-      ! 
+      !
       ! One level below
       en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) &
           &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
       en(:,:,2) = MAX(en(:,:,2), rn_emin )
-      z_elem_a(:,:,2) = 0._wp 
+      z_elem_a(:,:,2) = 0._wp
       z_elem_c(:,:,2) = 0._wp
       z_elem_b(:,:,2) = 1._wp
@@ -334,5 +334,5 @@
       ! Dirichlet conditions at k=1
       en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
-      en(:,:,1)       = MAX(en(:,:,1), rn_emin)      
+      en(:,:,1)       = MAX(en(:,:,1), rn_emin)
       z_elem_a(:,:,1) = en(:,:,1)
       z_elem_c(:,:,1) = 0._wp
@@ -357,5 +357,5 @@
       SELECT CASE ( nn_bc_bot )
       !
-      CASE ( 0 )             ! Dirichlet 
+      CASE ( 0 )             ! Dirichlet
          !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin
          !                      ! Balance between the production and the dissipation terms
@@ -377,10 +377,10 @@
                z_elem_c(ji,jj,ibotm1) = 0._wp
                z_elem_b(ji,jj,ibotm1) = 1._wp
-               en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
+               en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
             END DO
          END DO
          !
       CASE ( 1 )             ! Neumman boundary condition
-         !                      
+         !
 !CDIR NOVERRCHK
          DO jj = 2, jpjm1
@@ -428,5 +428,5 @@
          END DO
       END DO
-      !                                            ! set the minimum value of tke 
+      !                                            ! set the minimum value of tke
       en(:,:,:) = MAX( en(:,:,:), rn_emin )
 
@@ -470,5 +470,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn 
+                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn
                END DO
             END DO
@@ -489,5 +489,5 @@
                !
                ! psi / k
-               zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 
+               zratio = psi(ji,jj,jk) / eb(ji,jj,jk)
                !
                ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable)
@@ -509,5 +509,5 @@
                zesh2 = dir * ( prod + buoy )          + (1._wp - dir ) * prod                        ! production term
                zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
-               !                                                        
+               !
                ! building the matrix
                zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
@@ -551,5 +551,5 @@
       z_elem_c(:,:,2) = 0._wp
       z_elem_b(:,:,2) = 1._wp
-      ! 
+      !
       !
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
@@ -575,5 +575,5 @@
       psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
 
-      !   
+      !
       !
       END SELECT
@@ -585,5 +585,5 @@
       !
       !
-      CASE ( 0 )             ! Dirichlet 
+      CASE ( 0 )             ! Dirichlet
          !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0
          !                      ! Balance between the production and the dissipation terms
@@ -610,5 +610,5 @@
          !
       CASE ( 1 )             ! Neumman boundary condition
-         !                      
+         !
 !CDIR NOVERRCHK
          DO jj = 2, jpjm1
@@ -692,5 +692,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 
+                  eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)
                END DO
             END DO
@@ -719,10 +719,10 @@
                eps(ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
                mxln(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
-               ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 
+               ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
                IF (ln_length_lim) mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) )
             END DO
          END DO
-      END DO 
+      END DO
 
       !
@@ -846,6 +846,6 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_gls_init  ***
-      !!                     
-      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
+      !!
+      !! ** Purpose :   Initialization of the vertical eddy diffivity and
       !!      viscosity when using a gls turbulent closure scheme
       !!
@@ -1066,10 +1066,10 @@
          !
       END SELECT
-    
+
       !                                !* Set Schmidt number for psi diffusion in the wave breaking case
       !                                     ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
       !                                     !  or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
       IF( ln_sigpsi ) THEN
-         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 
+         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf
          ! Verification: retrieve Burchard (2001) results by uncomenting the line below:
          ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work
@@ -1079,5 +1079,5 @@
          rsc_psi0 = rsc_psi
       ENDIF
- 
+
       !                                !* Shear free turbulence parameters
       !
@@ -1125,11 +1125,11 @@
       rc04  = rc03 * rc0
       rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking
-      rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking 
+      rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking
       zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )
-      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer 
+      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer
       rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness
-      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness 
+      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness
       rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi
-      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 
+      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking
 
       rfact_tke = -0.5_wp / rsc_tke * rdt                                ! Cst used for the Diffusion term of tke
@@ -1146,5 +1146,5 @@
          avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
       END DO
-      !                              
+      !
       CALL gls_rst( nit000, 'READ' )   !* read or initialize all required files
       !
@@ -1157,9 +1157,9 @@
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE ts_rst  ***
-      !!                     
+      !!
       !! ** Purpose :   Read or write TKE file (en) in restart file
       !!
       !! ** Method  :   use of IOM library
-      !!                if the restart does not contain TKE, en is either 
+      !!                if the restart does not contain TKE, en is either
       !!                set to rn_emin or recomputed (nn_igls/=0)
       !!----------------------------------------------------------------------
@@ -1173,5 +1173,5 @@
       !!----------------------------------------------------------------------
       !
-      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
          !                                   ! ---------------
          IF( ln_rstart ) THEN                   !* Read the restart file
@@ -1190,8 +1190,8 @@
                CALL iom_get( numror, jpdom_autoglo, 'avmv'  , avmv   )
                CALL iom_get( numror, jpdom_autoglo, 'mxln'  , mxln   )
-            ELSE                        
+            ELSE
                IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop'
                en  (:,:,:) = rn_emin
-               mxln(:,:,:) = 0.05        
+               mxln(:,:,:) = 0.05
                avt_k (:,:,:) = avt (:,:,:)
                avm_k (:,:,:) = avm (:,:,:)
@@ -1203,5 +1203,5 @@
             IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'
             en  (:,:,:) = rn_emin
-            mxln(:,:,:) = 0.05       
+            mxln(:,:,:) = 0.05
          ENDIF
          !
@@ -1209,8 +1209,8 @@
          !                                   ! -------------------
          IF(lwp) WRITE(numout,*) '---- gls-rst ----'
-         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     ) 
+         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
          CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
          CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
-         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 
+         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
          CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
          CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln   )
@@ -1241,3 +1241,2 @@
    !!======================================================================
 END MODULE zdfgls
-
