Index: /LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90	(revision 4803)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90	(revision 4804)
@@ -23,9 +23,38 @@
 !=======================================================================
 
-
-
 USE lmdz_atke_turbulence_ini, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual, ri0, ri1
 USE lmdz_atke_turbulence_ini, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes, rg, rd, rv, atke_ok_vdiff
 USE lmdz_atke_turbulence_ini, ONLY : viscom, viscoh, clmix, clmixshear, iflag_atke_lmix, lmin, smmin, cn
+
+!!-------------------------------------------------------------------------------------------------------------
+! integer :: iflag_atke        ! flag that controls options in atke_compute_km_kh
+! integer :: iflag_atke_lmix   ! flag that controls the calculation of mixing length in atke
+! integer :: iflag_num_atke    ! flag that controls the numerical treatment of diffusion coeffiient calculation
+
+! logical :: atke_ok_vdiff    ! activate vertical diffusion of TKE or not
+! logical :: atke_ok_virtual  ! account for vapor for flottability
+
+! real :: kappa = 0.4         ! Von Karman constant
+! real :: cn                  ! Sm value at Ri=0
+! real :: cinf                ! Sm value at Ri=-Inf
+! real :: ri0                 ! Richardson number near zero to guarantee continuity in slope of Sm (stability function) at Ri=0
+! real :: ri1                 ! Richardson number near zero to guarantee continuity in slope of Pr (Prandlt's number) at Ri=0
+! real :: lmin                ! minimum mixing length corresponding to the Kolmogov dissipation scale  in planetary atmospheres (Chen et al 2016, JGR Atmos)
+! real :: ctkes               ! coefficient for surface TKE
+! real :: clmixshear          ! coefficient for mixing length depending on local wind shear
+
+! Tunable parameters for the ATKE scheme and their range of values
+!!-------------------------------------------------------------------------------------------------------------
+! real :: cepsilon            ! controls the value of the dissipation length scale, range [1.2 - 10]
+! real :: cke                 ! controls the value of the diffusion coefficient of TKE, range [1 - 5]
+! real :: l0                  ! asymptotic mixing length far from the ground [m] (Sun et al 2011, JAMC), range [15 - 75]
+! real :: clmix               ! controls the value of the mixing length in stratified conditions, range [0.1 - 2]
+! real :: ric                 ! critical Richardson number controlling the slope of Sm in stable conditions, range [0.19 - 0.25]
+! real :: smmin               ! minimum value of Sm in very stable conditions (defined here as minsqrt(Ez/Ek)) at large Ri, range [0.025 - 0.1]
+! real :: pr_neut             ! neutral value of the Prandtl number (Ri=0), range [0.7 - 1]
+! real :: pr_slope            ! linear slope of Pr with Ri in the very stable regime, range [3 - 5]
+! real :: cinffac             ! ratio between cinf and cn controlling the convective limit of Sm, range [1.2 - 5.0]
+! real :: pr_asym             ! value of Prandlt in the convective limit(Ri=-Inf), range [0.3 - 0.5]
+!!-------------------------------------------------------------------------------------------------------------
 
 implicit none
@@ -96,17 +125,17 @@
 ! results should not depend on the choice of preff
 DO ilay=1,nlay
-     DO igrid = 1, ngrid
+    DO igrid = 1, ngrid
         theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd)
-     END DO
+    END DO
 END DO
 
 ! account for water vapor mass for buoyancy calculation
 IF (atke_ok_virtual) THEN
-  DO ilay=1,nlay
-     DO igrid = 1, ngrid
-        rvap=max(0.,qvap(igrid,ilay)/(1.-qvap(igrid,ilay)))
-        theta(igrid,ilay)=theta(igrid,ilay)*(1.+rvap/(RD/RV))/(1.+rvap)
-     END DO
-  END DO
+    DO ilay=1,nlay
+        DO igrid = 1, ngrid
+            rvap=max(0.,qvap(igrid,ilay)/(1.-qvap(igrid,ilay)))
+            theta(igrid,ilay)=theta(igrid,ilay)*(1.+rvap/(RD/RV))/(1.+rvap)
+        END DO
+    END DO
 ENDIF
 
@@ -123,7 +152,7 @@
 
 DO ilay=1,nlay
-   DO igrid=1,ngrid 
-      z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay))
-   ENDDO
+    DO igrid=1,ngrid 
+        z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay))
+    ENDDO
 ENDDO
 
@@ -150,6 +179,6 @@
                                 + Ri(igrid,ilay) * pr_slope
             IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN
-               call abort_physic("atke_compute_km_kh", &
-               'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1)
+                call abort_physic("atke_compute_km_kh", &
+                'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1)
             ENDIF
         END IF
@@ -166,59 +195,55 @@
 
 IF (iflag_atke_lmix .EQ. 1 ) THEN
-! Blackadar formulation (~kappa l) + buoyancy length scale (Deardoff 1980) for very stable conditions
-   DO ilay=2,nlay
-      DO igrid=1,ngrid
-          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
-          IF (N2(igrid,ilay) .GT. 0.) THEN
-             lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay))
-             lstrat=max(lstrat,lmin)
-             !Inverse interpolation, Van de Wiel et al. 2010
-             l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
-          ENDIF 
-      ENDDO
-   ENDDO
-
-ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN
-! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms
-! implies 2 tuning coefficients clmix and clmixshear
-DO ilay=2,nlay
-      DO igrid=1,ngrid
-          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
-          IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN
-             lstrat=min(clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)), &
+    ! Blackadar formulation (~kappa l) + buoyancy length scale (Deardoff 1980) for very stable conditions
+        DO ilay=2,nlay
+            DO igrid=1,ngrid
+                l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
+                IF (N2(igrid,ilay) .GT. 0.) THEN
+                    lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay))
+                    lstrat=max(lstrat,lmin)
+                    !Inverse interpolation, Van de Wiel et al. 2010
+                    l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
+                ENDIF 
+            ENDDO
+        ENDDO
+
+    ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN
+    ! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms
+    ! implies 2 tuning coefficients clmix and clmixshear
+    DO ilay=2,nlay
+        DO igrid=1,ngrid
+            l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
+            IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN
+                lstrat=min(clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)), &
                     clmixshear*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay)))
-             lstrat=max(lstrat,lmin)
-             !Inverse interpolation, Van de Wiel et al. 2010   
-             l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
-          ENDIF
-      ENDDO
-   ENDDO
-
-ELSE IF (iflag_atke_lmix .EQ. 3 ) THEN
-! add effect of wind shear on lstrat following grisogono 2010, qjrms
-! keeping a single tuning coefficient clmix
-DO ilay=2,nlay
-      DO igrid=1,ngrid
-          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
-          IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN
-             lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))*(1.0+Ri(igrid,ilay)/(2.*Prandtl(igrid,ilay)))
-             lstrat=max(lstrat,lmin)
-             !Inverse interpolation, Van de Wiel et al. 2010   
-             l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
-          ENDIF
-      ENDDO
-   ENDDO
-
-
+                lstrat=max(lstrat,lmin)
+                !Inverse interpolation, Van de Wiel et al. 2010   
+                l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
+            ENDIF
+        ENDDO
+    ENDDO
+
+    ELSE IF (iflag_atke_lmix .EQ. 3 ) THEN
+    ! add effect of wind shear on lstrat following grisogono 2010, qjrms
+    ! keeping a single tuning coefficient clmix
+    DO ilay=2,nlay
+        DO igrid=1,ngrid
+            l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
+            IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN
+                lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))*(1.0+Ri(igrid,ilay)/(2.*Prandtl(igrid,ilay)))
+                lstrat=max(lstrat,lmin)
+                !Inverse interpolation, Van de Wiel et al. 2010   
+                l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
+            ENDIF
+        ENDDO
+    ENDDO
 
 ELSE
-! default Blackadar formulation: neglect effect of local stratification and shear
-
-   DO ilay=2,nlay+1
-      DO igrid=1,ngrid
-          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
-      ENDDO
-
-   ENDDO
+    ! default Blackadar formulation: neglect effect of local stratification and shear
+    DO ilay=2,nlay+1
+        DO igrid=1,ngrid
+            l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
+        ENDDO
+    ENDDO
 ENDIF
 
@@ -228,84 +253,84 @@
 IF (iflag_atke == 0) THEN
 
-! stationary solution (dtke/dt=0)
-
-   DO ilay=2,nlay
-        DO igrid=1,ngrid
-            tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
-            shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
-        ENDDO
-    ENDDO
+    ! stationary solution (dtke/dt=0)
+
+    DO ilay=2,nlay
+            DO igrid=1,ngrid
+                tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
+                shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
+            ENDDO
+        ENDDO
 
 ELSE IF (iflag_atke == 1) THEN
 
-! full implicit scheme resolved with a second order polynomial equation
-! default solution which shows fair convergence properties
+    ! full implicit scheme resolved with a second order polynomial equation
+    ! default solution which shows fair convergence properties
+        DO ilay=2,nlay
+            DO igrid=1,ngrid
+            qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
+            delta=(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime)**2. &
+                    +4.*(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime*qq + &
+                    2.*l_exchange(igrid,ilay)*l_exchange(igrid,ilay)*cepsilon*Sm(igrid,ilay) &
+                    *shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)))
+            qq=(-2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime + sqrt(delta))/2.
+            qq=max(0.,qq)
+            tke(igrid,ilay)=0.5*(qq**2) 
+            ENDDO
+        ENDDO
+
+
+ELSE IF (iflag_atke == 2) THEN
+
+    ! semi implicit scheme when l does not depend on tke
+    ! positive-guaranteed if pr slope in stable condition >1
+
     DO ilay=2,nlay
-        DO igrid=1,ngrid
-           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
-           delta=(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime)**2. &
-                 +4.*(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime*qq + &
-                 2.*l_exchange(igrid,ilay)*l_exchange(igrid,ilay)*cepsilon*Sm(igrid,ilay) &
-                 *shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)))
-           qq=(-2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime + sqrt(delta))/2.
-           qq=max(0.,qq)
-           tke(igrid,ilay)=0.5*(qq**2) 
-        ENDDO
-    ENDDO
-
-
-ELSE IF (iflag_atke == 2) THEN
-
-! semi implicit scheme when l does not depend on tke
-! positive-guaranteed if pr slope in stable condition >1
-
-   DO ilay=2,nlay
-        DO igrid=1,ngrid
-           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
-           qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.)      & 
-               *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
-               /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) 
-           tke(igrid,ilay)=0.5*(qq**2) 
-        ENDDO
-    ENDDO
+            DO igrid=1,ngrid
+            qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
+            qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.)      & 
+                *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
+                /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) 
+            tke(igrid,ilay)=0.5*(qq**2) 
+            ENDDO
+        ENDDO
 
 
 ELSE IF (iflag_atke == 3) THEN
-! numerical resolution adapted from that in MAR (Deleersnijder 1992)
-! positively defined by construction
-
-    DO ilay=2,nlay
-        DO igrid=1,ngrid
-           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
-           IF (Ri(igrid,ilay) .LT. 0.) THEN
-              netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))
-              netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))
-           ELSE
-              netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ &
-                      l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay)
-              netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)
-           ENDIF
-           qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss)
-           tke(igrid,ilay)=0.5*(qq**2)
-        ENDDO
-    ENDDO
+    ! numerical resolution adapted from that in MAR (Deleersnijder 1992)
+    ! positively defined by construction
+
+        DO ilay=2,nlay
+            DO igrid=1,ngrid
+            qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
+            IF (Ri(igrid,ilay) .LT. 0.) THEN
+                netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))
+                netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))
+            ELSE
+                netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ &
+                        l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay)
+                netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)
+            ENDIF
+            qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss)
+            tke(igrid,ilay)=0.5*(qq**2)
+            ENDDO
+        ENDDO
 
 ELSE IF (iflag_atke == 4) THEN
-! semi implicit scheme from Arpege (V. Masson methodology with 
-! Taylor expansion of the dissipation term)
-    DO ilay=2,nlay
-        DO igrid=1,ngrid
-           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
-           qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
-             +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) &
-             /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
-           qq=max(0.,qq)
-           tke(igrid,ilay)=0.5*(qq**2)
-        ENDDO
-    ENDDO
+    ! semi implicit scheme from Arpege (V. Masson methodology with 
+    ! Taylor expansion of the dissipation term)
+        DO ilay=2,nlay
+            DO igrid=1,ngrid
+            qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
+            qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
+                +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) &
+                /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
+            qq=max(0.,qq)
+            tke(igrid,ilay)=0.5*(qq**2)
+            ENDDO
+        ENDDO
 
 
 ELSE
-   call abort_physic("atke_compute_km_kh", &
+    call abort_physic("atke_compute_km_kh", &
         'numerical treatment of TKE not possible yet', 1)
 
@@ -316,5 +341,5 @@
 
 DO igrid=1,ngrid
- tke(igrid,nlay+1)=0.
+    tke(igrid,nlay+1)=0.
 END DO
 
@@ -324,6 +349,6 @@
 ! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note)
 DO igrid=1,ngrid
- ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2))
- tke(igrid,1)=ctkes*(ustar**2)
+    ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2))
+    tke(igrid,1)=ctkes*(ustar**2)
 END DO
 
@@ -332,5 +357,5 @@
 !==========================
 IF (atke_ok_vdiff) THEN
-   CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
+    CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
 ENDIF
 
@@ -395,5 +420,5 @@
 DO ilay=2,nlay
     DO igrid=1,ngrid
-       Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke 
+        Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke 
     ENDDO
 ENDDO
@@ -423,6 +448,6 @@
 
 DO igrid=1,ngrid
-  CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay)
-  DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay)
+    CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay)
+    DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay)
 ENDDO
 
@@ -431,5 +456,5 @@
     DO igrid=1,ngrid
         CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) &
-                       / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
+                        / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
         DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
     ENDDO
Index: /LMDZ6/trunk/libf/phylmd/lmdz_atke_turbulence_ini.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_atke_turbulence_ini.F90	(revision 4803)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_atke_turbulence_ini.F90	(revision 4804)
@@ -4,176 +4,175 @@
 
 ! declaration of constants and parameters
-save
 
-integer :: iflag_atke, iflag_num_atke, iflag_atke_lmix
-  !$OMP THREADPRIVATE(iflag_atke, iflag_num_atke, iflag_atke_lmix)
-  real :: kappa = 0.4 ! Von Karman constant
-  !$OMP THREADPRIVATE(kappa)
-  real :: cinffac
-  !$OMP THREADPRIVATE(cinffac)
-  real :: l0,ric,ri0,ri1,cinf,cn,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke
-  !$OMP THREADPRIVATE(l0,ri0,ri1,ric,cinf,cn,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke)
-  integer :: lunout,prt_level
-  !$OMP THREADPRIVATE(lunout,prt_level)
-  real :: rg, rd, rpi, rcpd, rv
-  !$OMP THREADPRIVATE(rg, rd, rpi, rcpd, rv)
-  real :: viscom, viscoh
-  !$OMP THREADPRIVATE(viscom,viscoh)
-  real :: lmin=0.01              ! minimum mixing length corresponding to the Kolmogov dissipation scale 
-                                 ! in planetary atmospheres (Chen et al 2016, JGR Atmos)
-  !$OMP THREADPRIVATE(lmin)
-  logical :: atke_ok_vdiff, atke_ok_virtual
-  !$OMP THREADPRIVATE(atke_ok_vdiff,atke_ok_virtual)
+real, save, protected  :: rg, rd, rpi, rcpd, rv, viscom, viscoh        
+!$OMP THREADPRIVATE(rg, rd, rpi, rcpd, rv, viscom, viscoh)     
+
+integer, save, protected :: iflag_atke        ! flag that controls options in atke_compute_km_kh
+integer, save, protected :: iflag_atke_lmix   ! flag that controls the calculation of mixing length in atke
+integer, save, protected :: iflag_num_atke    ! flag that controls the numerical treatment of diffusion coeffiient calculation
+!$OMP THREADPRIVATE(iflag_atke, iflag_atke_lmix, iflag_num_atke)
+
+logical, save, protected :: atke_ok_vdiff    ! activate vertical diffusion of TKE or not
+logical, save, protected :: atke_ok_virtual  ! account for vapor for flottability
+!$OMP THREADPRIVATE(atke_ok_vdiff, atke_ok_virtual)
+
+real, save, protected :: kappa = 0.4         ! Von Karman constant
+real, save, protected :: cn                  ! Sm value at Ri=0
+real, save, protected :: cinf                ! Sm value at Ri=-Inf
+real, save, protected :: ri0                 ! Richardson number near zero to guarantee continuity in slope of Sm (stability function) at Ri=0
+real, save, protected :: ri1                 ! Richardson number near zero to guarantee continuity in slope of Pr (Prandlt's number) at Ri=0
+real, save, protected :: lmin = 0.01         ! minimum mixing length corresponding to the Kolmogov dissipation scale  in planetary atmospheres (Chen et al 2016, JGR Atmos)
+real, save, protected :: ctkes               ! coefficient for surface TKE
+real, save, protected :: clmixshear          ! coefficient for mixing length depending on local wind shear
+!$OMP THREADPRIVATE(kappa, cn, cinf, ri0, ri1, lmin, ctkes, clmixshear)
+
+
+! Tunable parameters for the ATKE scheme and their range of values
+!!-------------------------------------------------------------------------------------------------------------
+real, save, protected :: cepsilon            ! controls the value of the dissipation length scale, range [1.2 - 10]
+real, save, protected :: cke                 ! controls the value of the diffusion coefficient of TKE, range [1 - 5]
+real, save, protected :: l0                  ! asymptotic mixing length far from the ground [m] (Sun et al 2011, JAMC), range [15 - 75]
+real, save, protected :: clmix               ! controls the value of the mixing length in stratified conditions, range [0.1 - 2]
+real, save, protected :: ric                 ! critical Richardson number controlling the slope of Sm in stable conditions, range [0.19 - 0.25]
+real, save, protected :: smmin               ! minimum value of Sm in very stable conditions (defined here as minsqrt(Ez/Ek)) at large Ri, range [0.025 - 0.1]
+real, save, protected :: pr_neut             ! neutral value of the Prandtl number (Ri=0), range [0.7 - 1]
+real, save, protected :: pr_slope            ! linear slope of Pr with Ri in the very stable regime, range [3 - 5]
+real, save, protected :: cinffac             ! ratio between cinf and cn controlling the convective limit of Sm, range [1.2 - 5.0]
+real, save, protected :: pr_asym             ! value of Prandlt in the convective limit(Ri=-Inf), range [0.3 - 0.5]
+!$OMP THREADPRIVATE(cepsilon, cke, l0, clmix, ric, smmin, pr_neut, pr_slope, cinffac, pr_asym)
+!!-------------------------------------------------------------------------------------------------------------
 
 CONTAINS
 
 SUBROUTINE atke_ini(rg_in, rd_in, rpi_in, rcpd_in, rv_in, viscom_in, viscoh_in)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE  atke_ini  ***
+      !!
+      !! ** Purpose :   Initialization of the atke module and choice of some constants
+      !!               
+      !!----------------------------------------------------------------------
 
-   USE ioipsl_getin_p_mod, ONLY : getin_p
+        USE ioipsl_getin_p_mod, ONLY : getin_p                                                                                                      
 
-  real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in, rv_in, viscom_in, viscoh_in
+      ! input arguments (universal constants for planet)
+      !-------------------------------------------------
+      real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in, rv_in, viscom_in, viscoh_in
+      !!----------------------------------------------------------------------
+      
+      rg=rg_in          ! gravity acceleration
+      rd=rd_in          ! dry gas constant (R/M, R=perfect gas constant and M is the molar mass of the fluid)
+      rpi=rpi_in        ! Pi number
+      rcpd=rcpd_in      ! cp per unit mass of the gas
+      rv=rv_in          ! water vapor constant (for simulations in Earth's atmosphere)
+      viscom=viscom_in  ! kinematic molecular viscosity for momentum
+      viscoh=viscoh_in  ! kinematic molecular viscosity for heat
 
 
-  ! input arguments (universal constants for planet)
-  !-------------------------------------------------
- 
-  ! gravity acceleration
-  rg=rg_in
-  ! dry gas constant (R/M, R=perfect gas constant and M is the molar mass of the fluid)
-  rd=rd_in
-  ! Pi number
-  rpi=rpi_in
-  ! cp per unit mass of the gas
-  rcpd=rcpd_in
-  ! water vapor constant (for simulations in Earth's atmosphere)
-  rv=rv_in
-  ! kinematic molecular viscosity for momentum
-  viscom=viscom_in
-  ! kinematic molecular viscosity for heat
-  viscoh=viscoh_in
+      ! Read flag values in .def files
+      !-------------------------------
+
+      ! flag that controls options in atke_compute_km_kh
+      iflag_atke=0
+      CALL getin_p('iflag_atke',iflag_atke)
+
+      ! flag that controls the calculation of mixing length in atke
+      iflag_atke_lmix=0
+      CALL getin_p('iflag_atke_lmix',iflag_atke_lmix)
+
+      if (iflag_atke .eq. 0 .and. iflag_atke_lmix>0) then
+            call abort_physic("atke_turbulence_ini", &
+            'stationary scheme must use mixing length formulation not depending on tke', 1)
+      endif
+
+      ! activate vertical diffusion of TKE or not
+      atke_ok_vdiff=.false.
+      CALL getin_p('atke_ok_vdiff',atke_ok_vdiff)
+
+      ! account for vapor for flottability
+      atke_ok_virtual=.true.
+      CALL getin_p('atke_ok_virtual',atke_ok_virtual)
 
 
-  !viscom=1.46E-5
-  !viscoh=2.06E-5
+      ! flag that controls the numerical treatment of diffusion coeffiient calculation
+      iflag_num_atke=0
+      CALL getin_p('iflag_num_atke',iflag_num_atke)
 
+      ! asymptotic mixing length in neutral conditions [m]
+      ! Sun et al 2011, JAMC 
+      ! between 10 and 40
+      l0=15.0 
+      CALL getin_p('atke_l0',l0)
 
-  ! Read flag values in .def files
-  !-------------------------------
+      ! critical Richardson number 
+      ric=0.25
+      CALL getin_p('atke_ric',ric)
 
+      ! constant for tke dissipation calculation
+      cepsilon=5.87 ! default value as in yamada4
+      CALL getin_p('atke_cepsilon',cepsilon)
 
-  ! flag that controls options in atke_compute_km_kh
-  iflag_atke=0
-  CALL getin_p('iflag_atke',iflag_atke)
+      ! calculation of cn = Sm value at Ri=0
+      ! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0
+      cn=(1./sqrt(cepsilon))**(2/3)
 
-  ! flag that controls the calculation of mixing length in atke
-  iflag_atke_lmix=0
-  CALL getin_p('iflag_atke_lmix',iflag_atke_lmix)
+      ! asymptotic value of Sm for Ri=-Inf
+      cinffac=2.0
+      CALL getin_p('atke_cinffac',cinffac)
+      cinf=cinffac*cn
+      if (cinf .le. cn) then
+            call abort_physic("atke_turbulence_ini", &
+            'cinf cannot be lower than cn', 1)
+      endif
 
-  if (iflag_atke .eq. 0 .and. iflag_atke_lmix>0) then
-        call abort_physic("atke_turbulence_ini", &
-        'stationary scheme must use mixing length formulation not depending on tke', 1)
-  endif
+      ! coefficient for surface TKE 
+      ! following Lenderink & Holtslag 2004, ctkes=(cepsilon)**(2/3)
+      ! (provided by limit condition in neutral conditions)
+      ctkes=(cepsilon)**(2./3.)
 
-  ! activate vertical diffusion of TKE or not
-  atke_ok_vdiff=.false.
-  CALL getin_p('atke_ok_vdiff',atke_ok_vdiff)
+      ! slope of Pr=f(Ri) for stable conditions
+      pr_slope=5.0 ! default value from Zilitinkevich et al. 2005
+      CALL getin_p('atke_pr_slope',pr_slope)
+      if (pr_slope .le. 1) then
+            call abort_physic("atke_turbulence_ini", &
+            'pr_slope has to be greater than 1 for consistency of the tke scheme', 1)
+      endif
 
+      ! value of turbulent prandtl number in neutral conditions (Ri=0)
+      pr_neut=0.8
+      CALL getin_p('atke_pr_neut',pr_neut)
 
-  ! account for vapor for flottability
-  atke_ok_virtual=.true.
-  CALL getin_p('atke_ok_virtual',atke_ok_virtual)
+      ! asymptotic turbulent prandt number value for Ri=-Inf
+      pr_asym=0.4
+      CALL getin_p('atke_pr_asym',pr_asym)
+      if (pr_asym .ge. pr_neut) then
+            call abort_physic("atke_turbulence_ini", &
+            'pr_asym must be be lower than pr_neut', 1)
+      endif
 
+      ! coefficient for mixing length depending on local stratification
+      clmix=0.5
+      CALL getin_p('atke_clmix',clmix)
 
-  ! flag that controls the numerical treatment of diffusion coeffiient calculation
-  iflag_num_atke=0
-  CALL getin_p('iflag_num_atke',iflag_num_atke)
+      ! coefficient for mixing length depending on local wind shear
+      clmixshear=0.5
+      CALL getin_p('atke_clmixshear',clmixshear)
 
-  ! asymptotic mixing length in neutral conditions [m]
-  ! Sun et al 2011, JAMC 
-  ! between 10 and 40
-  l0=15.0 
-  CALL getin_p('atke_l0',l0)
+      ! minimum anisotropy coefficient (defined here as minsqrt(Ez/Ek)) at large Ri. 
+      ! From Zilitinkevich et al. 2013, it equals sqrt(0.03)~0.17  
+      smmin=0.17
+      CALL getin_p('atke_smmin',smmin)
 
-  ! critical Richardson number 
-  ric=0.25
-  CALL getin_p('atke_ric',ric)
+      ! ratio between the eddy diffusivity coeff for tke wrt that for momentum
+      ! default value from Lenderink et al. 2004
+      cke=2.
+      CALL getin_p('atke_cke',cke)
 
+      ! calculation of Ri0 such that continuity in slope of Sm at Ri=0
+      ri0=2./rpi*(cinf - cn)*ric/cn
 
-  ! constant for tke dissipation calculation
-  cepsilon=5.87 ! default value as in yamada4
-  CALL getin_p('atke_cepsilon',cepsilon)
+      ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0
+      ri1 = -2./rpi * (pr_asym - pr_neut)
 
-
-  ! calculation of cn = Sm value at Ri=0
-  ! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0
-  cn=(1./sqrt(cepsilon))**(2/3)
-
-  ! asymptotic value of Sm for Ri=-Inf
-  cinffac=2.0
-  CALL getin_p('atke_cinffac',cinffac)
-  cinf=cinffac*cn
-  if (cinf .le. cn) then
-        call abort_physic("atke_turbulence_ini", &
-        'cinf cannot be lower than cn', 1)
-  endif
-
-
- ! coefficient for surface TKE 
- ! following Lenderink & Holtslag 2004, ctkes=(cepsilon)**(2/3)
- ! (provided by limit condition in neutral conditions)
-  ctkes=(cepsilon)**(2./3.)
-
-  ! slope of Pr=f(Ri) for stable conditions
-  pr_slope=5.0 ! default value from Zilitinkevich et al. 2005
-  CALL getin_p('atke_pr_slope',pr_slope)
-  if (pr_slope .le. 1) then
-        call abort_physic("atke_turbulence_ini", &
-        'pr_slope has to be greater than 1 for consistency of the tke scheme', 1)
-  endif
-
-  ! value of turbulent prandtl number in neutral conditions (Ri=0)
-  pr_neut=0.8
-  CALL getin_p('atke_pr_neut',pr_neut)
-
-
- ! asymptotic turbulent prandt number value for Ri=-Inf
-  pr_asym=0.4
-  CALL getin_p('atke_pr_asym',pr_asym)
-  if (pr_asym .ge. pr_neut) then
-        call abort_physic("atke_turbulence_ini", &
-        'pr_asym must be be lower than pr_neut', 1)
-  endif
-
-
-
-  ! coefficient for mixing length depending on local stratification
-  clmix=0.5
-  CALL getin_p('atke_clmix',clmix)
- 
-  ! coefficient for mixing length depending on local wind shear
-  clmixshear=0.5
-  CALL getin_p('atke_clmixshear',clmixshear)
-
-
-  ! minimum anisotropy coefficient (defined here as minsqrt(Ez/Ek)) at large Ri. 
-  ! From Zilitinkevich et al. 2013, it equals sqrt(0.03)~0.17  
-  smmin=0.17
-  CALL getin_p('atke_smmin',smmin)
-
-
-  ! ratio between the eddy diffusivity coeff for tke wrt that for momentum
-  ! default value from Lenderink et al. 2004
-  cke=2.
-  CALL getin_p('atke_cke',cke)
-
-
-  ! calculation of Ri0 such that continuity in slope of Sm at Ri=0
-  ri0=2./rpi*(cinf - cn)*ric/cn
-
-  ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0
-  ri1 = -2./rpi * (pr_asym - pr_neut)
-
-
- RETURN
+RETURN
 
 END SUBROUTINE atke_ini
