Changeset 13191 for branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ZDF

Timestamp:

2020-07-01T15:01:22+02:00 (4 years ago)

Author:

jwhile

Message:

Updates for 1d runnig

Location:

branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ZDF

Files:

• zdfbfr.F90 (modified) (4 diffs)
• zdfevd.F90 (modified) (7 diffs)
• zdfgls.F90 (modified) (39 diffs)
• zdftke.F90 (modified) (37 diffs)

branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

@@ -107 +107 @@
          WRITE(numout,*) '~~~~~~~~'
       ENDIF
-      ! 
+      !
+#if defined key_traldf_c2d || key_traldf_c3d
       IF( ln_stopack .AND. ( nn_spp_bfr > 0 ) ) THEN
          bfrcoef2d = bfrcoef2d0
          CALL spp_gen(kt, bfrcoef2d, nn_spp_bfr, rn_bfr_sd, jk_spp_bfr)
       ENDIF
+#else
+      IF ( ln_stopack .AND. ( nn_spp_bfr > 0 ) ) &
+         & CALL ctl_stop( 'zdf_bfr: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
       !
       IF( nn_bfr == 2 ) THEN                 ! quadratic bottom friction only
@@ -140 +146 @@
                END DO
             END IF
-         !   
+         !
          ELSE
             zbfrt(:,:) = bfrcoef2d(:,:)
@@ -183 +189 @@
                DO ji = 2, jpim1
                   ! (ISF) ========================================================================
-                  ikbu = miku(ji,jj)         ! ocean top level at u- and v-points 
+                  ikbu = miku(ji,jj)         ! ocean top level at u- and v-points
                   ikbv = mikv(ji,jj)         ! (1st wet ocean u- and v-points)
                   !
@@ -363 +369 @@
             bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable
          ENDIF
-         
+
          IF ( ln_isfcav ) THEN
             IF(ln_tfr2d) THEN

branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90

@@ -20 +20 @@
    USE zdfkpp          ! KPP vertical mixing
    USE trd_oce         ! trends: ocean variables
-   USE trdtra          ! trends manager: tracers 
+   USE trdtra          ! trends manager: tracers
    USE in_out_manager  ! I/O manager
    USE iom             ! for iom_put
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE timing          ! Timing
-   USE stopack         
+   USE stopack
 
    IMPLICIT NONE
@@ -45 +45 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_evd  ***
-      !!                   
+      !!
       !! ** Purpose :   Local increased the vertical eddy viscosity and diffu-
       !!      sivity coefficients when a static instability is encountered.
       !!
       !! ** Method  :   avt, avm, and the 4 neighbouring avmu, avmv coefficients
-      !!      are set to avevd (namelist parameter) if the water column is 
+      !!      are set to avevd (namelist parameter) if the water column is
       !!      statically unstable (i.e. if rn2 < -1.e-12 )
       !!
@@ -77 +77 @@
       zavt_evd(:,:,:) = avt(:,:,:)           ! set avt prior to evd application
 
+#if defined key_traldf_c2d || key_traldf_c3d
       IF( ln_stopack .AND. ( nn_spp_aevd > 0 ) ) THEN
          rn_avevd0(:,:) = rn_avevd
          CALL spp_gen(kt, rn_avevd0, nn_spp_aevd, rn_aevd_sd, jk_spp_aevd)
       ENDIF
+#else
+      IF ( ln_stopack .AND. ( nn_spp_aevd > 0 ) ) &
+         & CALL ctl_stop( 'zdf_evd: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
 
       SELECT CASE ( nn_evdm )
@@ -88 +94 @@
          zavm_evd(:,:,:) = avm(:,:,:)           ! set avm prior to evd application
          !
-         DO jk = 1, jpkm1 
+         DO jk = 1, jpkm1
             DO jj = 2, jpj             ! no vector opt.
                DO ji = 2, jpi
@@ -106 +112 @@
                END DO
             END DO
-         END DO 
+         END DO
          CALL lbc_lnk( avt , 'W', 1. )   ;   CALL lbc_lnk( avm , 'W', 1. )   ! Lateral boundary conditions
          CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
@@ -113 +119 @@
          CALL iom_put( "avm_evd", zavm_evd )              ! output this change
          !
-      CASE DEFAULT         ! enhance vertical eddy diffusivity only (if rn2<-1.e-12) 
+      CASE DEFAULT         ! enhance vertical eddy diffusivity only (if rn2<-1.e-12)
          DO jk = 1, jpkm1
-!!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL! 
+!!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL!
             DO jj = 1, jpj             ! loop over the whole domain (no lbc_lnk call)
                DO ji = 1, jpi
 #if defined key_zdfkpp
                   ! no evd mixing in the boundary layer with KPP
-                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12  .AND.  fsdepw(ji,jj,jk) > hkpp(ji,jj)  )   &          
+                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12  .AND.  fsdepw(ji,jj,jk) > hkpp(ji,jj)  )   &
 #else
                   IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 )   &
@@ -129 +135 @@
          END DO
          !
-      END SELECT 
+      END SELECT
 
       zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd

branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

@@ -2 +2 @@
    !!======================================================================
    !!                       ***  MODULE  zdfgls  ***
-   !! Ocean physics:  vertical mixing coefficient computed from the gls 
+   !! Ocean physics:  vertical mixing coefficient computed from the gls
    !!                 turbulent closure parameterization
    !!======================================================================
@@ -16 +16 @@
    !!   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 +31 @@
    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)
    USE stopack
 
@@ -62 +62 @@
    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
@@ -73 +73 @@
    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
@@ -134 +134 @@
       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         !   -      -
@@ -140 +140 @@
       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
@@ -157 +157 @@
       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
 
@@ -166 +166 @@
          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         
-         DO ji = fs_2, fs_jpim1   ! vector opt.         
+!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
@@ -227 +227 @@
       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)
@@ -276 +276 @@
                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
@@ -295 +295 @@
                ! 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
@@ -317 +317 @@
       ! 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
@@ -335 +335 @@
       ! 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
@@ -358 +358 @@
       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
@@ -378 +378 @@
                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
@@ -429 +429 @@
          END DO
       END DO
-      !                                            ! set the minimum value of tke 
+      !                                            ! set the minimum value of tke
       en(:,:,:) = MAX( en(:,:,:), rn_emin )
 
@@ -471 +471 @@
             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
@@ -490 +490 @@
                !
                ! 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)
@@ -510 +510 @@
                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)
@@ -552 +552 @@
       z_elem_c(:,:,2) = 0._wp
       z_elem_b(:,:,2) = 1._wp
-      ! 
+      !
       !
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
@@ -576 +576 @@
       psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
 
-      !   
+      !
       !
       END SELECT
@@ -586 +586 @@
       !
       !
-      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
@@ -611 +611 @@
          !
       CASE ( 1 )             ! Neumman boundary condition
-         !                      
+         !
 !CDIR NOVERRCHK
          DO jj = 2, jpjm1
@@ -693 +693 @@
             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
@@ -720 +720 @@
                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
 
       !
@@ -807 +807 @@
             END DO
          END DO
+#if defined key_traldf_c2d || key_traldf_c3d
          IF( ln_stopack) THEN
             IF( nn_spp_avt > 0 ) CALL spp_gen(kt, avt(:,:,jk), nn_spp_avt, rn_avt_sd,jk)
             IF( nn_spp_avm > 0 ) CALL spp_gen(kt, avm(:,:,jk), nn_spp_avm, rn_avm_sd,jk)
-         ENDIF 
+         ENDIF
+#else
+         IF ( ln_stopack ) &
+            & CALL ctl_stop( 'zdf_gls: parameter perturbation will only work with '// &
+                             'key_traldf_c2d or key_traldf_c3d')
+#endif
       END DO
       !
@@ -851 +857 @@
       !!----------------------------------------------------------------------
       !!                  ***  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
       !!
@@ -1071 +1077 @@
          !
       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
@@ -1084 +1090 @@
          rsc_psi0 = rsc_psi
       ENDIF
- 
+
       !                                !* Shear free turbulence parameters
       !
@@ -1130 +1136 @@
       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
@@ -1151 +1157 @@
          avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
       END DO
-      !                              
+      !
       CALL gls_rst( nit000, 'READ' )   !* read or initialize all required files
       !
@@ -1162 +1168 @@
       !!---------------------------------------------------------------------
       !!                   ***  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)
       !!----------------------------------------------------------------------
@@ -1178 +1184 @@
       !!----------------------------------------------------------------------
       !
-      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
          !                                   ! ---------------
          IF( ln_rstart ) THEN                   !* Read the restart file
@@ -1197 +1203 @@
                CALL iom_get( numror, jpdom_autoglo, 'mxln'  , mxln   )
                IF(nn_timing == 2)  CALL timing_stop('iom_rstget')
-            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 (:,:,:)
@@ -1210 +1216 @@
             IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'
             en  (:,:,:) = rn_emin
-            mxln(:,:,:) = 0.05       
+            mxln(:,:,:) = 0.05
          ENDIF
          !
@@ -1217 +1223 @@
          IF(lwp) WRITE(numout,*) '---- gls-rst ----'
          IF(nn_timing == 2)  CALL timing_start('iom_rstput')
-         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   )
@@ -1239 +1245 @@
       WRITE(*,*) 'zdf_gls_init: You should not have seen this print! error?'
    END SUBROUTINE zdf_gls_init
-   
+
    SUBROUTINE zdf_gls( kt )          ! Empty routine
    IMPLICIT NONE
-      INTEGER, INTENT(in) ::   kt ! ocean time step   
+      INTEGER, INTENT(in) ::   kt ! ocean time step
       WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt
    END SUBROUTINE zdf_gls
-   
+
    SUBROUTINE gls_rst( kt, cdrw )          ! Empty routine
    IMPLICIT NONE
@@ -1256 +1262 @@
    !!======================================================================
 END MODULE zdfgls
-

branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

@@ -2 +2 @@
    !!======================================================================
    !!                       ***  MODULE  zdftke  ***
-   !! Ocean physics:  vertical mixing coefficient computed from the tke 
+   !! Ocean physics:  vertical mixing coefficient computed from the tke
    !!                 turbulent closure parameterization
    !!=====================================================================
@@ -22 +22 @@
    !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
    !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
-   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term 
-   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl 
+   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term
+   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl
    !!                 !                                + cleaning of the parameters + bugs correction
    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
@@ -52 +52 @@
    USE wrk_nemo       ! work arrays
    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)
 #if defined key_agrif
    USE agrif_opa_interp
@@ -74 +74 @@
    INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
    REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
-   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation 
+   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation
    REAL(wp) ::   rn_ebb    ! coefficient of the surface input of tke
    REAL(wp) ::   rn_emin   ! minimum value of tke           [m2/s2]
@@ -102 +102 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4
-   
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -118 +118 @@
       !!----------------------------------------------------------------------
       ALLOCATE(                                                                    &
-         &      efr  (jpi,jpj)     , e_niw(jpi,jpj,jpk) ,                         &      
+         &      efr  (jpi,jpj)     , e_niw(jpi,jpj,jpk) ,                         &
 #if defined key_c1d
          &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          &
@@ -147 +147 @@
       !!         surface: en = max( rn_emin0, rn_ebb * taum )
       !!         bottom : en = rn_emin
-      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 
-      !!
-      !!        The now Turbulent kinetic energy is computed using the following 
+      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
+      !!
+      !!        The now Turbulent kinetic energy is computed using the following
       !!      time stepping: implicit for vertical diffusion term, linearized semi
-      !!      implicit for kolmogoroff dissipation term, and explicit forward for 
-      !!      both buoyancy and shear production terms. Therefore a tridiagonal 
+      !!      implicit for kolmogoroff dissipation term, and explicit forward for
+      !!      both buoyancy and shear production terms. Therefore a tridiagonal
       !!      linear system is solved. Note that buoyancy and shear terms are
       !!      discretized in a energy conserving form (Bruchard 2002).
@@ -160 +160 @@
       !!
       !!        The now vertical eddy vicosity and diffusivity coefficients are
-      !!      given by: 
+      !!      given by:
       !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
-      !!              avt = max( avmb, pdl * avm                 )  
+      !!              avt = max( avmb, pdl * avm                 )
       !!              eav = max( avmb, avm )
       !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
-      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1 
+      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1
       !!
       !! ** Action  :   compute en (now turbulent kinetic energy)
@@ -180 +180 @@
       !
       IF( kt /= nit000 ) THEN   ! restore before value to compute tke
-         avt (:,:,:) = avt_k (:,:,:) 
-         avm (:,:,:) = avm_k (:,:,:) 
-         avmu(:,:,:) = avmu_k(:,:,:) 
-         avmv(:,:,:) = avmv_k(:,:,:) 
-      ENDIF
-      !
+         avt (:,:,:) = avt_k (:,:,:)
+         avm (:,:,:) = avm_k (:,:,:)
+         avmu(:,:,:) = avmu_k(:,:,:)
+         avmv(:,:,:) = avmv_k(:,:,:)
+      ENDIF
+      !
+#if defined key_traldf_c2d || key_traldf_c3d
       IF( ln_stopack) THEN
          IF( nn_spp_tkelc > 0 ) THEN
@@ -208 +209 @@
          ENDIF
       ENDIF
+#else
+      IF ( ln_stopack ) &
+         & CALL ctl_stop( 'zdf_tke: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
       !
       CALL tke_tke      ! now tke (en)
@@ -213 +219 @@
       CALL tke_avn      ! now avt, avm, avmu, avmv
       !
-      avt_k (:,:,:) = avt (:,:,:) 
-      avm_k (:,:,:) = avm (:,:,:) 
-      avmu_k(:,:,:) = avmu(:,:,:) 
-      avmv_k(:,:,:) = avmv(:,:,:) 
+      avt_k (:,:,:) = avt (:,:,:)
+      avm_k (:,:,:) = avm (:,:,:)
+      avmu_k(:,:,:) = avmu(:,:,:)
+      avmv_k(:,:,:) = avmv(:,:,:)
       !
 #if defined key_agrif
-      ! Update child grid f => parent grid 
+      ! Update child grid f => parent grid
       IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only
-#endif      
+#endif
       IF (  kt == nitend ) THEN
         DEALLOCATE ( rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 )
@@ -271 +277 @@
       CALL wrk_alloc( jpi,jpj, zfact2 )
       CALL wrk_alloc( jpi,jpj, zfact3 )
-      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
+      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )
       !
       zbbrau = rn_ebb0 / rau0       ! Local constant initialisation
-      zfact1 = -.5_wp * rdt 
+      zfact1 = -.5_wp * rdt
       zfact2 = 1.5_wp * rdt * rn_ediss0
       zfact3 = 0.5_wp       * rn_ediss0
@@ -294 +300 @@
          END DO
       END DO
-      
+
 !!bfr   - start commented area
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -333 +339 @@
          imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
          DO jk = jpkm1, 2, -1
-            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us 
+            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us
                DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
                   zus  = zcof * taum(ji,jj)
@@ -341 +347 @@
          END DO
          !                               ! finite LC depth
-         DO jj = 1, jpj 
+         DO jj = 1, jpj
             DO ji = 1, jpi
                zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
@@ -378 +384 @@
             DO ji = 1, jpi
                avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
-                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
+                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   &
                   &                            / (  fse3uw_n(ji,jj,jk)               &
                   &                              *  fse3uw_b(ji,jj,jk)  )
@@ -407 +413 @@
                   !                                                           ! shear prod. at w-point weightened by mask
                zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
-                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
+                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
                   !
                zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
@@ -464 +470 @@
       END DO
 
-      !                                 ! Save TKE prior to nn_etau addition  
-      e_niw(:,:,:) = en(:,:,:)  
-      !  
+      !                                 ! Save TKE prior to nn_etau addition
+      e_niw(:,:,:) = en(:,:,:)
+      !
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  TKE due to surface and internal wave breaking
@@ -508 +514 @@
                   ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
                   zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
-                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress 
-                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean 
+                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress
+                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean
                   zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau(ji,jj) * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
@@ -517 +523 @@
          END DO
       ELSEIF( nn_etau == 4 ) THEN       !* column integral independant of htau (rn_efr must be scaled up)
-         IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau 
+         IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -535 +541 @@
       CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
       !
-      DO jk = 2, jpkm1                             ! TKE budget: near-inertial waves term  
-         DO jj = 2, jpjm1  
-            DO ji = fs_2, fs_jpim1   ! vector opt.  
-               e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk)  
-            END DO  
-         END DO  
-      END DO  
-      !  
-      CALL lbc_lnk( e_niw, 'W', 1. )  
+      DO jk = 2, jpkm1                             ! TKE budget: near-inertial waves term
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( e_niw, 'W', 1. )
       !
       CALL wrk_dealloc( jpi,jpj, imlc )    ! integer
-      CALL wrk_dealloc( jpi,jpj, zhlc ) 
-      CALL wrk_dealloc( jpi,jpj, zbbrau ) 
-      CALL wrk_dealloc( jpi,jpj, zfact2 ) 
-      CALL wrk_dealloc( jpi,jpj, zfact3 ) 
-      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
+      CALL wrk_dealloc( jpi,jpj, zhlc )
+      CALL wrk_dealloc( jpi,jpj, zbbrau )
+      CALL wrk_dealloc( jpi,jpj, zfact2 )
+      CALL wrk_dealloc( jpi,jpj, zfact3 )
+      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tke_tke')
@@ -563 +569 @@
       !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
       !!
-      !! ** Method  :   At this stage, en, the now TKE, is known (computed in 
-      !!              the tke_tke routine). First, the now mixing lenth is 
+      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
+      !!              the tke_tke routine). First, the now mixing lenth is
       !!      computed from en and the strafification (N^2), then the mixings
       !!      coefficients are computed.
       !!              - Mixing length : a first evaluation of the mixing lengh
       !!      scales is:
-      !!                      mxl = sqrt(2*en) / N  
+      !!                      mxl = sqrt(2*en) / N
       !!      where N is the brunt-vaisala frequency, with a minimum value set
       !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
-      !!        The mixing and dissipative length scale are bound as follow : 
+      !!        The mixing and dissipative length scale are bound as follow :
       !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
       !!                        zmxld = zmxlm = mxl
       !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
-      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 
+      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
       !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
       !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
-      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to 
-      !!                    the surface to obtain ldown. the resulting length 
+      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
+      !!                    the surface to obtain ldown. the resulting length
       !!                    scales are:
-      !!                        zmxld = sqrt( lup * ldown ) 
+      !!                        zmxld = sqrt( lup * ldown )
       !!                        zmxlm = min ( lup , ldown )
       !!              - Vertical eddy viscosity and diffusivity:
       !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
-      !!                      avt = max( avmb, pdlr * avm )  
+      !!                      avt = max( avmb, pdlr * avm )
       !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
       !!
@@ -601 +607 @@
       IF( nn_timing == 1 )  CALL timing_start('tke_avn')
 
-      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
+      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )
 
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -610 +616 @@
       !
       ! initialisation of interior minimum value (avoid a 2d loop with mikt)
-      zmxlm(:,:,:)  = rmxl_min    
+      zmxlm(:,:,:)  = rmxl_min
       zmxld(:,:,:)  = rmxl_min
       !
@@ -620 +626 @@
             END DO
          END DO
-      ELSE 
+      ELSE
          zmxlm(:,:,1) = rn_mxl0
       ENDIF
@@ -638 +644 @@
       !                     !* Physical limits for the mixing length
       !
-      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value 
+      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value
       zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
       !
@@ -738 +744 @@
             END DO
          END DO
+#if defined key_traldf_c2d || key_traldf_c3d
          IF( ln_stopack) THEN
             IF(nn_spp_avt > 0 ) CALL spp_gen( 1, avt(:,:,jk), nn_spp_avt, rn_avt_sd, jk_spp_avt, jk)
             IF(nn_spp_avm > 0 ) CALL spp_gen( 1, avm(:,:,jk), nn_spp_avm, rn_avm_sd, jk_spp_avm, jk)
          ENDIF
+#else
+         IF ( ln_stopack  ) &
+            & CALL ctl_stop( 'tke_avn: parameter perturbation will only work with '// &
+                             'key_traldf_c2d or key_traldf_c3d')
+#endif
       END DO
       CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
@@ -787 +799 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
+      CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tke_avn')
@@ -797 +809 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_tke_init  ***
-      !!                     
-      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
+      !!
+      !! ** Purpose :   Initialization of the vertical eddy diffivity and
       !!              viscosity when using a tke turbulent closure scheme
       !!
@@ -814 +826 @@
          &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   &
          &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   &
-         &                 nn_etau , nn_htau  , rn_efr , rn_c   
+         &                 nn_etau , nn_htau  , rn_efr , rn_c
       !!----------------------------------------------------------------------
 
@@ -884 +896 @@
       ALLOCATE( rn_ebb0  (jpi,jpj) ) ; rn_ebb0   = rn_ebb
       ALLOCATE( rn_efr0  (jpi,jpj) ) ; rn_efr0   = rn_efr
-      
+
       IF( nn_etau == 2  ) THEN
           ierr = zdf_mxl_alloc()
@@ -896 +908 @@
 
       !                               !* depth of penetration of surface tke
-      IF( nn_etau /= 0 ) THEN      
+      IF( nn_etau /= 0 ) THEN
          htau(:,:) = 0._wp
          SELECT CASE( nn_htau )             ! Choice of the depth of penetration
@@ -902 +914 @@
             htau(:,:) = 10._wp
          CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
-            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )            
+            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )
          CASE( 2 )                                 ! fraction of depth-integrated TKE within mixed-layer
             rhtau = -1._wp / LOG( 1._wp - rn_c )
@@ -945 +957 @@
       END DO
       dissl(:,:,:) = 1.e-12_wp
-      !                              
+      !
       CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files
       !
@@ -954 +966 @@
      !!---------------------------------------------------------------------
      !!                   ***  ROUTINE tke_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 
-     !!                set to rn_emin or recomputed 
+     !!                if the restart does not contain TKE, en is either
+     !!                set to rn_emin or recomputed
      !!----------------------------------------------------------------------
      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
@@ -968 +980 @@
      !!----------------------------------------------------------------------
      !
-     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
         !                                   ! ---------------
         IF( ln_rstart ) THEN                   !* Read the restart file
@@ -1006 +1018 @@
         avmu_k(:,:,:) = avmu(:,:,:)
         avmv_k(:,:,:) = avmv(:,:,:)
-        
+
      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
         !                                   ! -------------------