Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 1282)
+++ /trunk/LMDZ.GENERIC/README	(revision 1283)
@@ -1024,2 +1024,10 @@
   over the whole planet. This should be further imroved with computation of
   means (possibly area weighed), etc.
+
+== 05/06/2014 == EM
+Bug fixes:
+- hice() in physiq.F90 must be a saved array.
+- bad use of min/max on arrays in h2o_cloudrad (radii_mod.F90) which actually
+  ended up setting all array elements to the same value.
+And some cosmetic cleanup in rain.F90, vdif_kc.F and turbdiff.F90
+
Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 1282)
+++ /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 1283)
@@ -348,5 +348,5 @@
 
 !     included by BC for hydrology
-      real hice(ngrid)
+      real,allocatable,save :: hice(:)
 
 !     included by RW to test water conservation (by routine)
@@ -454,4 +454,5 @@
         ALLOCATE(cloudfrac(ngrid,nlayermx))
         ALLOCATE(totcloudfrac(ngrid))
+        ALLOCATE(hice(ngrid))
         ALLOCATE(qsurf_hist(ngrid,nq))
         ALLOCATE(reffrad(ngrid,nlayermx,naerkind))
@@ -976,5 +977,5 @@
          zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
          if (tracer) then 
-           pdq(1:ngrid,1:nlayermx,1:nq)=pdq(1:ngrid,1:nlayermx,1:nq)+ zdqdif(1:ngrid,1:nlayermx,1:nq) 
+           pdq(1:ngrid,1:nlayermx,1:nq)=pdq(1:ngrid,1:nlayermx,1:nq)+ zdqdif(1:ngrid,1:nlayermx,1:nq)
            dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
          end if ! of if (tracer)
@@ -1064,5 +1065,5 @@
          pdt(1:ngrid,1:nlayermx)    = pdt(1:ngrid,1:nlayermx) + zdhadj(1:ngrid,1:nlayermx)*zpopsk(1:ngrid,1:nlayermx)
          zdtadj(1:ngrid,1:nlayermx) = zdhadj(1:ngrid,1:nlayermx)*zpopsk(1:ngrid,1:nlayermx) ! for diagnostic only
- 
+
          if(tracer) then 
             pdq(1:ngrid,1:nlayermx,1:nq) = pdq(1:ngrid,1:nlayermx,1:nq) + zdqadj(1:ngrid,1:nlayermx,1:nq) 
@@ -1112,5 +1113,4 @@
               zdtc,zdtsurfc,pdpsrf,zduc,zdvc,               &
               zdqc)
-
 
          pdt(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)+zdtc(1:ngrid,1:nlayermx)
@@ -1246,5 +1246,4 @@
                dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
                rainout(1:ngrid)             = zdqsrain(1:ngrid)+zdqssnow(1:ngrid) ! diagnostic   
-
 
                !-------------------------
Index: /trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90	(revision 1282)
+++ /trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90	(revision 1283)
@@ -38,7 +38,7 @@
       Implicit none
 
-#include "callkeys.h"
-#include "dimensions.h"
-#include "dimphys.h"
+      include "callkeys.h"
+      include "dimensions.h"
+      include "dimphys.h"
 
       integer,intent(in) :: ngrid
@@ -140,8 +140,8 @@
       Implicit none
 
-#include "callkeys.h"
-#include "dimensions.h"
-#include "dimphys.h"
-#include "comcstfi.h"
+      include "callkeys.h"
+      include "dimensions.h"
+      include "dimphys.h"
+      include "comcstfi.h"
 
       integer,intent(in) :: ngrid
@@ -200,16 +200,16 @@
       Implicit none
 
-#include "callkeys.h"
-#include "dimensions.h"
-#include "dimphys.h"
-#include "comcstfi.h"
+      include "callkeys.h"
+      include "dimensions.h"
+      include "dimphys.h"
+      include "comcstfi.h"
 
       integer,intent(in) :: ngrid
 
       real, intent(in) :: pql(ngrid,nlayermx) !condensed water mixing ratios (kg/kg)
-      real, intent(out) :: reffliq(ngrid,nlayermx),reffice(ngrid,nlayermx)     !liquid and ice water particle radii (K)
+      real, intent(out) :: reffliq(ngrid,nlayermx),reffice(ngrid,nlayermx)     !liquid and ice water particle radii (m)
 
       real,external :: CBRT            
-      
+      integer :: i,k
 
       if (radfixed) then
@@ -217,9 +217,14 @@
 	 reffice(1:ngrid,1:nlayermx)= rad_h2o_ice
       else
-	 reffliq(1:ngrid,1:nlayermx)  = CBRT( 3*pql(1:ngrid,1:nlayermx)/(4*Nmix_h2o*pi*rhowater) )
-	 reffliq(1:ngrid,1:nlayermx)  = min(max(reffliq(1:ngrid,1:nlayermx),1.e-6),1000.e-6)
-	 reffice(1:ngrid,1:nlayermx)  = CBRT( 3*pql(1:ngrid,1:nlayermx)/(4*Nmix_h2o_ice*pi*rhowaterice) )
-	 reffice(1:ngrid,1:nlayermx)  = min(max(reffice(1:ngrid,1:nlayermx),1.e-6),1000.e-6)
-      end if
+	 do k=1,nlayermx
+	   do i=1,ngrid
+	     reffliq(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o*pi*rhowater))
+	     reffliq(i,k) = min(max(reffliq(i,k),1.e-6),1000.e-6)
+	   
+	     reffice(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o_ice*pi*rhowaterice))
+	     reffice(i,k) = min(max(reffice(i,k),1.e-6),1000.e-6)
+	   enddo
+	 enddo
+      endif
 
    end subroutine h2o_cloudrad
@@ -243,8 +248,8 @@
       Implicit none
 
-#include "callkeys.h"
-#include "dimensions.h"
-#include "dimphys.h"
-#include "comcstfi.h"
+      include "callkeys.h"
+      include "dimensions.h"
+      include "dimphys.h"
+      include "comcstfi.h"
 
       integer,intent(in) :: ngrid,nq
@@ -289,7 +294,7 @@
       Implicit none
 
-#include "callkeys.h"
-#include "dimensions.h"
-#include "dimphys.h"
+      include "callkeys.h"
+      include "dimensions.h"
+      include "dimphys.h"
 
       integer,intent(in) :: ngrid
@@ -317,7 +322,7 @@
       Implicit none
 
-#include "callkeys.h"
-#include "dimensions.h"
-#include "dimphys.h"
+      include "callkeys.h"
+      include "dimensions.h"
+      include "dimphys.h"
 
       integer,intent(in) :: ngrid
@@ -346,7 +351,7 @@
      Implicit none
 
-#include "callkeys.h"
-#include "dimensions.h"
-#include "dimphys.h"
+     include "callkeys.h"
+     include "dimensions.h"
+     include "dimphys.h"
 
       integer,intent(in) :: ngrid
Index: /trunk/LMDZ.GENERIC/libf/phystd/rain.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/rain.F90	(revision 1282)
+++ /trunk/LMDZ.GENERIC/libf/phystd/rain.F90	(revision 1283)
@@ -2,9 +2,8 @@
 
 
-! to use  'getin'
-  use ioipsl_getincom 
+  use ioipsl_getincom, only: getin
   use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2,Psat_water,Tsat_water,rhowater
   use radii_mod, only: h2o_cloudrad
-  USE tracer_h
+  USE tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
   implicit none
 
@@ -23,8 +22,8 @@
 !==================================================================
 
-#include "dimensions.h"
-#include "dimphys.h"
-#include "comcstfi.h"
-#include "callkeys.h"
+  include "dimensions.h"
+  include "dimphys.h"
+  include "comcstfi.h"
+  include "callkeys.h"
 
 !     Arguments
@@ -149,5 +148,5 @@
 
          firstcall = .false.
-      ENDIF
+      ENDIF ! of IF (firstcall)
 
 !     GCM -----> subroutine variables
@@ -173,15 +172,9 @@
 
 !     Initialise the outputs
-      DO k = 1, nlayermx
-      DO i = 1, ngrid
-         d_t(i,k) = 0.0
-         d_q(i,k) = 0.0
-         d_ql(i,k) = 0.0
-      ENDDO
-      ENDDO
-      DO i = 1, ngrid
-         zrfl(i)  = 0.0
-         zrfln(i) = 0.0
-      ENDDO
+      d_t(1:ngrid,1:nlayermx) = 0.0
+      d_q(1:ngrid,1:nlayermx) = 0.0
+      d_ql(1:ngrid,1:nlayermx) = 0.0
+      zrfl(1:ngrid) = 0.0
+      zrfln(1:ngrid) = 0.0
 
       ! calculate saturation mixing ratio
@@ -203,9 +196,8 @@
       ENDDO
 
-
       ! Vertical loop (from top to bottom)
       ! We carry the rain with us and calculate that added by warm/cold precipitation
       ! processes and that subtracted by evaporation at each level.
-      DO 9999 k = nlayermx, 1, -1
+      DO k = nlayermx, 1, -1
 
          IF (evap_prec) THEN ! note no rneb dependence!
@@ -236,11 +228,9 @@
 		     
 
-               ENDIF
+               ENDIF ! of IF (zrfl(i) .GT.0.)
             ENDDO
-         ENDIF
-
-         DO i = 1, ngrid
-            zoliq(i) = 0.0
-         ENDDO
+         ENDIF ! of IF (evap_prec)
+
+         zoliq(1:ngrid) = 0.0
 
 
@@ -367,5 +357,5 @@
          endif ! if precip_scheme=1
 
- 9999 continue
+      ENDDO ! of DO k = nlayermx, 1, -1
 
 !     Rain or snow on the ground
@@ -385,19 +375,12 @@
 
 !     now subroutine -----> GCM variables
-      DO k = 1, nlayermx
-         DO i = 1, ngrid
-            
-            if(evap_prec)then
-               dqrain(i,k,i_vap) = d_q(i,k)/ptimestep
-               d_t(i,k)          = d_t(i,k)/ptimestep
-            else
-               dqrain(i,k,i_vap) = 0.0
-               d_t(i,k)          = 0.0
-            endif
-            dqrain(i,k,i_ice) = d_ql(i,k)/ptimestep
-         
-         ENDDO
-      ENDDO
-
-      RETURN
+      if (evap_prec) then
+        dqrain(1:ngrid,1:nlayermx,i_vap)=dqrain(1:ngrid,1:nlayermx,i_vap)/ptimestep
+        d_t(1:ngrid,1:nlayermx)=d_t(1:ngrid,1:nlayermx)/ptimestep
+      else
+        dqrain(1:ngrid,1:nlayermx,i_vap)=0.0
+        d_t(1:ngrid,1:nlayermx)=0.0
+      endif
+      dqrain(1:ngrid,1:nlayermx,i_ice) = d_ql(1:ngrid,1:nlayermx)/ptimestep
+
     end subroutine rain
Index: /trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90	(revision 1282)
+++ /trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90	(revision 1283)
@@ -9,7 +9,6 @@
       use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol,Psat_water
       use radcommon_h, only : sigma, glat
-      USE surfdat_h
-      USE comgeomfi_h
-      USE tracer_h
+      use surfdat_h, only: dryness
+      use tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
 
       implicit none
@@ -42,39 +41,43 @@
 !     ------------
 
-#include "dimensions.h"
-#include "dimphys.h"
-#include "comcstfi.h"
-#include "callkeys.h"
+      include "dimensions.h"
+      include "dimphys.h"
+      include "comcstfi.h"
+      include "callkeys.h"
 
 !     arguments
 !     ---------
-      INTEGER ngrid,nlay
-      REAL ptimestep
-      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
-      REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
-      REAL pu(ngrid,nlay),pv(ngrid,nlay)
-      REAL pt(ngrid,nlay),ppopsk(ngrid,nlay)
-      REAL ptsrf(ngrid),pemis(ngrid)
-      REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay)
-      REAL pdtfi(ngrid,nlay)
-      REAL pfluxsrf(ngrid)
-      REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
-      REAL pdtdif(ngrid,nlay)
-      REAL pdtsrf(ngrid),sensibFlux(ngrid),pcapcal(ngrid)
-      REAL pq2(ngrid,nlay+1)
+      INTEGER,INTENT(IN) :: ngrid,nlay
+      REAL,INTENT(IN) :: ptimestep
+      REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
+      REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
+      REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
+      REAL,INTENT(IN) :: pt(ngrid,nlay),ppopsk(ngrid,nlay)
+      REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
+      REAL,INTENT(IN) :: pemis(ngrid)
+      REAL,INTENT(IN) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay)
+      REAL,INTENT(IN) :: pdtfi(ngrid,nlay)
+      REAL,INTENT(IN) :: pfluxsrf(ngrid)
+      REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
+      REAL,INTENT(OUT) :: pdtdif(ngrid,nlay)
+      REAL,INTENT(OUT) :: pdtsrf(ngrid) ! tendency (K/s) on surface temperature
+      REAL,INTENT(OUT) :: sensibFlux(ngrid)
+      REAL,INTENT(IN) :: pcapcal(ngrid)
+      REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
       
-      integer rnat(ngrid)      
+      integer,intent(in) :: rnat(ngrid)      
+      LOGICAL,INTENT(IN) :: lastcall ! not used
 
 !     Arguments added for condensation
-      logical lecrit
-      REAL pz0
+      logical,intent(in) :: lecrit ! not used.
+      REAL,INTENT(IN) :: pz0
 
 !     Tracers
 !     --------
-      integer nq 
-      real pqsurf(ngrid,nq)
-      real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) 
-      real pdqdif(ngrid,nlay,nq) 
-      real pdqsdif(ngrid,nq) 
+      integer,intent(in) :: nq 
+      real,intent(in) :: pqsurf(ngrid,nq)
+      real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) 
+      real,intent(out) :: pdqdif(ngrid,nlay,nq) 
+      real,intent(out) :: pdqsdif(ngrid,nq) 
       
 !     local
@@ -102,8 +105,6 @@
       REAL zx_alf1(ngrid),zx_alf2(ngrid)
 
-      LOGICAL firstcall
-      SAVE firstcall
+      LOGICAL,SAVE :: firstcall=.true.
       
-      LOGICAL lastcall
 !     Tracers
 !     -------
@@ -115,5 +116,4 @@
       REAL rho(ngrid)         ! near-surface air density
       REAL qsat(ngrid),psat(ngrid)
-      DATA firstcall/.true./
       REAL kmixmin
 
@@ -242,5 +242,5 @@
 
       call vdif_kc(ngrid,ptimestep,g,pzlev,pzlay,pu,pv,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature 
-
+      
 !     Adding eddy mixing to mimic 3D general circulation in 1D
 !     R. Wordsworth & F. Forget (2010)
@@ -731,5 +731,3 @@
 !      endif
 
-
-      return
       end
Index: /trunk/LMDZ.GENERIC/libf/phystd/vdif_kc.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/vdif_kc.F	(revision 1282)
+++ /trunk/LMDZ.GENERIC/libf/phystd/vdif_kc.F	(revision 1283)
@@ -27,15 +27,15 @@
 c 
 c.......................................................................
-      INTEGER ngrid
-      REAL dt,g
-      REAL zlev(ngrid,nlayermx+1)
-      REAL zlay(ngrid,nlayermx)
-      REAL u(ngrid,nlayermx)
-      REAL v(ngrid,nlayermx)
-      REAL teta(ngrid,nlayermx)
-      REAL cd(ngrid)
-      REAL q2(ngrid,nlayermx+1)
-      REAL km(ngrid,nlayermx+1)
-      REAL kn(ngrid,nlayermx+1)
+      INTEGER,INTENT(IN) :: ngrid
+      REAL,INTENT(IN) :: dt,g
+      REAL,INTENT(IN) :: zlev(ngrid,nlayermx+1)
+      REAL,INTENT(IN) :: zlay(ngrid,nlayermx)
+      REAL,INTENT(IN) :: u(ngrid,nlayermx)
+      REAL,INTENT(IN) :: v(ngrid,nlayermx)
+      REAL,INTENT(IN) :: teta(ngrid,nlayermx)
+      REAL,INTENT(IN) :: cd(ngrid)
+      REAL,INTENT(INOUT) :: q2(ngrid,nlayermx+1)
+      REAL,INTENT(OUT) :: km(ngrid,nlayermx+1)
+      REAL,INTENT(OUT) :: kn(ngrid,nlayermx+1)
 c.......................................................................
 c
@@ -50,5 +50,6 @@
 c
 c.......................................................................
-      INTEGER nlay,nlev
+      INTEGER,PARAMETER :: nlay=nlayermx
+      INTEGER,PARAMETER :: nlev=nlayermx+1
       REAL unsdz(ngrid,nlayermx)
       REAL unsdzdec(ngrid,nlayermx+1)
@@ -114,6 +115,6 @@
 c.......................................................................
       REAL gn
-      REAL gnmin
-      REAL gnmax
+      REAL,PARAMETER :: gnmin=-10.E+0
+      REAL,PARAMETER :: gnmax=0.0233E+0
       LOGICAL gninf
       LOGICAL gnsup
@@ -134,9 +135,15 @@
 c
 c.......................................................................
-      REAL kappa
-      REAL long0
-      REAL a1,a2,b1,b2,c1
-      REAL cn1,cn2
-      REAL cm1,cm2,cm3,cm4
+      REAL,PARAMETER :: kappa=0.4E+0
+      REAL,PARAMETER :: long0=160.E+0
+      REAL,PARAMETER :: a1=0.92E+0,a2=0.74E+0
+      REAL,PARAMETER :: b1=16.6E+0,b2=10.1E+0,c1=0.08E+0
+      REAL,PARAMETER :: cn1=a2*(1.E+0 -6.E+0 *a1/b1)
+      REAL,PARAMETER :: cn2=-3.E+0 *a2*(6.E+0 *a1+b2)
+      REAL,PARAMETER :: cm1=a1*(1.E+0 -3.E+0 *c1-6.E+0 *a1/b1)
+      REAL,PARAMETER :: cm2=a1*(-3.E+0 *a2*((b2-3.E+0 *a2)*
+     &          (1.E+0 -6.E+0 *a1/b1)-3.E+0 *c1*(b2+6.E+0 *a1)))
+      REAL,PARAMETER :: cm3=-3.E+0 *a2*(6.E+0 *a1+b2)
+      REAL,PARAMETER :: cm4=-9.E+0 *a1*a2
 c.......................................................................
 c
@@ -157,51 +164,26 @@
 c
 c.......................................................................
-      REAL q2min
-      REAL q2max
+      REAL,PARAMETER :: q2min=1.E-3
+      REAL,PARAMETER :: q2max=1.E+2
 c.......................................................................
 c knmin : borne inferieure de kn
 c kmmin : borne inferieure de km
 c.......................................................................
-      REAL knmin
-      REAL kmmin
+      REAL,PARAMETER :: knmin=1.E-5
+      REAL,PARAMETER :: kmmin=1.E-5
 c.......................................................................
       INTEGER ilay,ilev,igrid
       REAL tmp1,tmp2
 c.......................................................................
-      PARAMETER (kappa=0.4E+0)
-      PARAMETER (long0=160.E+0)
-      PARAMETER (gnmin=-10.E+0)
-      PARAMETER (gnmax=0.0233E+0)
-      PARAMETER (a1=0.92E+0)
-      PARAMETER (a2=0.74E+0)
-      PARAMETER (b1=16.6E+0)
-      PARAMETER (b2=10.1E+0)
-      PARAMETER (c1=0.08E+0)
-      PARAMETER (knmin=1.E-5)
-      PARAMETER (kmmin=1.E-5)
-      PARAMETER (q2min=1.E-3)
-      PARAMETER (q2max=1.E+2)
-      PARAMETER (nlay=nlayermx)
-      PARAMETER (nlev=nlayermx+1)
-c
-      PARAMETER (
-     &  cn1=a2*(1.E+0 -6.E+0 *a1/b1)
-     &          )
-      PARAMETER (
-     &  cn2=-3.E+0 *a2*(6.E+0 *a1+b2)
-     &          )
-      PARAMETER (
-     &  cm1=a1*(1.E+0 -3.E+0 *c1-6.E+0 *a1/b1)
-     &          )
-      PARAMETER (
-     &  cm2=a1*(-3.E+0 *a2*((b2-3.E+0 *a2)*(1.E+0 -6.E+0 *a1/b1)
-     &          -3.E+0 *c1*(b2+6.E+0 *a1)))
-     &          )
-      PARAMETER (
-     &  cm3=-3.E+0 *a2*(6.E+0 *a1+b2)
-     &          )
-      PARAMETER (
-     &  cm4=-9.E+0 *a1*a2
-     &          )
+c
+
+! initialization of local variables:
+      long(:,:)=0.
+      n2(:,:)=0.
+      sn(:,:)=0.
+      snq2(:,:)=0.
+      sm(:,:)=0.
+      smq2(:,:)=0.
+
 c.......................................................................
 c  traitment des valeur de q2 en entree
@@ -209,16 +191,16 @@
 c
       DO ilev=1,nlev
-                                                 DO igrid=1,ngrid 
+       DO igrid=1,ngrid 
         q2(igrid,ilev)=amax1(q2(igrid,ilev),q2min)
         q(igrid,ilev)=sqrt(q2(igrid,ilev))
-                                                 ENDDO
-      ENDDO
-c
-                                                 DO igrid=1,ngrid
-      tmp1=cd(igrid)*(u(igrid,1)**2+v(igrid,1)**2)
-      q2(igrid,1)=b1**(2.E+0/3.E+0)*tmp1
-      q2(igrid,1)=amax1(q2(igrid,1),q2min)
-      q(igrid,1)=sqrt(q2(igrid,1))
-                                                 ENDDO
+       ENDDO
+      ENDDO
+c
+      DO igrid=1,ngrid
+        tmp1=cd(igrid)*(u(igrid,1)**2+v(igrid,1)**2)
+        q2(igrid,1)=b1**(2.E+0/3.E+0)*tmp1
+        q2(igrid,1)=amax1(q2(igrid,1),q2min)
+        q(igrid,1)=sqrt(q2(igrid,1))
+      ENDDO
 c
 c.......................................................................
@@ -237,19 +219,24 @@
 c
       DO ilay=1,nlay
-                                                 DO igrid=1,ngrid
+       DO igrid=1,ngrid
         unsdz(igrid,ilay)=1.E+0/(zlev(igrid,ilay+1)-zlev(igrid,ilay))
-                                                 ENDDO
-      ENDDO
-                                                 DO igrid=1,ngrid
-      unsdzdec(igrid,1)=1.E+0/(zlay(igrid,1)-zlev(igrid,1))
-                                                 ENDDO
+       ENDDO
+      ENDDO
+
+      DO igrid=1,ngrid
+        unsdzdec(igrid,1)=1.E+0/(zlay(igrid,1)-zlev(igrid,1))
+      ENDDO
+
       DO ilay=2,nlay
-                                                 DO igrid=1,ngrid
-        unsdzdec(igrid,ilay)=1.E+0/(zlay(igrid,ilay)-zlay(igrid,ilay-1))
-                                                 ENDDO
-      ENDDO
-                                                 DO igrid=1,ngrid
-      unsdzdec(igrid,nlay+1)=1.E+0/(zlev(igrid,nlay+1)-zlay(igrid,nlay))
-                                                 ENDDO
+        DO igrid=1,ngrid
+          unsdzdec(igrid,ilay)=1.E+0/
+     &                          (zlay(igrid,ilay)-zlay(igrid,ilay-1))
+        ENDDO
+      ENDDO
+      
+      DO igrid=1,ngrid
+        unsdzdec(igrid,nlay+1)=1.E+0/
+     &                          (zlev(igrid,nlay+1)-zlay(igrid,nlay))
+      ENDDO
 c
 c.......................................................................
@@ -257,16 +244,16 @@
 c.......................................................................
 c
-                                                 DO igrid=1,ngrid
-      m2(igrid,1)=(unsdzdec(igrid,1)
+      DO igrid=1,ngrid
+        m2(igrid,1)=(unsdzdec(igrid,1)
      &                   *u(igrid,1))**2
      &                 +(unsdzdec(igrid,1)
      &                   *v(igrid,1))**2
-      m(igrid,1)=sqrt(m2(igrid,1))
-      mpre(igrid,1)=m(igrid,1)
-                                                 ENDDO
+        m(igrid,1)=sqrt(m2(igrid,1))
+        mpre(igrid,1)=m(igrid,1)
+      ENDDO
 c
 c-----------------------------------------------------------------------
       DO ilev=2,nlev-1
-                                                 DO igrid=1,ngrid
+       DO igrid=1,ngrid
 c-----------------------------------------------------------------------
 c
@@ -295,13 +282,14 @@
 c
 c-----------------------------------------------------------------------
-                                                 ENDDO
-      ENDDO
-c-----------------------------------------------------------------------
-c
-                                                 DO igrid=1,ngrid
-      m2(igrid,nlev)=m2(igrid,nlev-1)
-      m(igrid,nlev)=m(igrid,nlev-1)
-      mpre(igrid,nlev)=m(igrid,nlev)
-                                                 ENDDO
+       ENDDO
+      ENDDO
+c-----------------------------------------------------------------------
+c
+      DO igrid=1,ngrid
+        m2(igrid,nlev)=m2(igrid,nlev-1)
+        m(igrid,nlev)=m(igrid,nlev-1)
+        mpre(igrid,nlev)=m(igrid,nlev)
+      ENDDO
+      
 c
 c.......................................................................
@@ -311,5 +299,5 @@
 c-----------------------------------------------------------------------
       DO ilev=2,nlev-1
-                                                 DO igrid=1,ngrid
+       DO igrid=1,ngrid
 c-----------------------------------------------------------------------
 c
@@ -375,6 +363,6 @@
 c
 c-----------------------------------------------------------------------
-                                                  ENDDO
-      ENDDO
+       ENDDO ! of DO igrid=1,ngrid
+      ENDDO ! of DO ilev=2,nlev-1
 c-----------------------------------------------------------------------
 c
@@ -383,15 +371,13 @@
 c.......................................................................
 c
-                                                 DO igrid=1,ngrid
-      kn(igrid,1)=knmin
-      km(igrid,1)=kmmin
-      kmpre(igrid,1)=km(igrid,1)
-                                                 ENDDO
+      DO igrid=1,ngrid
+        kn(igrid,1)=knmin
+        km(igrid,1)=kmmin
+        kmpre(igrid,1)=km(igrid,1)
+      ENDDO
 c
 c-----------------------------------------------------------------------
       DO ilev=2,nlev-1
-                                                 DO igrid=1,ngrid
-c-----------------------------------------------------------------------
-c
+       DO igrid=1,ngrid
         kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
      &                                         *sn(igrid,ilev)
@@ -399,15 +385,13 @@
      &                                         *sm(igrid,ilev)
         kmpre(igrid,ilev)=km(igrid,ilev)
-c
-c-----------------------------------------------------------------------
-                                                 ENDDO
-      ENDDO
-c-----------------------------------------------------------------------
-c
-                                                 DO igrid=1,ngrid
-      kn(igrid,nlev)=kn(igrid,nlev-1)
-      km(igrid,nlev)=km(igrid,nlev-1)
-      kmpre(igrid,nlev)=km(igrid,nlev)
-                                                 ENDDO
+       ENDDO
+      ENDDO
+c-----------------------------------------------------------------------
+c
+      DO igrid=1,ngrid
+        kn(igrid,nlev)=kn(igrid,nlev-1)
+        km(igrid,nlev)=km(igrid,nlev-1)
+        kmpre(igrid,nlev)=km(igrid,nlev)
+      ENDDO
 c
 c.......................................................................
@@ -539,16 +523,12 @@
 c
 c
-                                                 DO igrid=1,ngrid
-      kn(igrid,1)=knmin
-      km(igrid,1)=kmmin
-      q2(igrid,nlev)=q2(igrid,nlev-1)
-      q(igrid,nlev)=q(igrid,nlev-1)
-      kn(igrid,nlev)=kn(igrid,nlev-1)
-      km(igrid,nlev)=km(igrid,nlev-1)
-                                                 ENDDO
+      DO igrid=1,ngrid
+        kn(igrid,1)=knmin
+        km(igrid,1)=kmmin
+        q2(igrid,nlev)=q2(igrid,nlev-1)
+        q(igrid,nlev)=q(igrid,nlev-1)
+        kn(igrid,nlev)=kn(igrid,nlev-1)
+        km(igrid,nlev)=km(igrid,nlev-1)
+      ENDDO
 
-c
-c.......................................................................
-c
-      RETURN
       END
