Index: /dynamico_lmdz/simple_physics/bash/beautify.sh
===================================================================
--- /dynamico_lmdz/simple_physics/bash/beautify.sh	(revision 4229)
+++ /dynamico_lmdz/simple_physics/bash/beautify.sh	(revision 4229)
@@ -0,0 +1,8 @@
+#!/bin/bash
+for x in $* ; do
+    # indent source file using emacs in batch mode
+    # following standard settings of emacs f90-mode
+    emacs -batch $x --eval '(f90-mode)' --eval '(indent-region (point-min) (point-max) nil)' --eval '(delete-trailing-whitespace)' -f save-buffer
+    # remove blank lines at end of file
+    sed -i '' -e :a -e '/^\n*$/{$d;N;};/\n$/ba' $x 
+done
Index: /dynamico_lmdz/simple_physics/phyparam/Makefile
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/Makefile	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/Makefile	(revision 4229)
@@ -22,5 +22,5 @@
 clean :
 	rm -f obj/* include/* lib/* xml/*
-
+	../bash/beautify.sh physics/*.F90
 %.so : $(OBJECTS)
 	$(F90) -shared $^ -o $@ 
Index: /dynamico_lmdz/simple_physics/phyparam/physics/astronomy.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/astronomy.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/astronomy.F90	(revision 4229)
@@ -1,4 +1,4 @@
 MODULE astronomy
-  
+
 #include "use_logging.h"
 
@@ -12,5 +12,5 @@
 
 CONTAINS
- 
+
   SUBROUTINE solarlong(pday,psollong)
     REAL, INTENT(IN) :: pday               ! jour de l'annee (le jour 0 correspondant a l'equinoxe)
@@ -18,14 +18,14 @@
     LOGICAL, PARAMETER ::  lwrite=.TRUE.
 
-! Local:
-! ------
+    ! Local:
+    ! ------
     REAL zanom,xref,zx0,zdx,zteta,zz
     INTEGER iter
 
-!--------------------------------------------------------
-! calcul de l'angle polaire et de la distance au soleil :
-! -------------------------------------------------------
+    !--------------------------------------------------------
+    ! calcul de l'angle polaire et de la distance au soleil :
+    ! -------------------------------------------------------
 
-!  calcul de l'zanomalie moyenne
+    !  calcul de l'zanomalie moyenne
 
     zz=(pday-peri_day)/year_day
@@ -33,7 +33,7 @@
     xref=abs(zanom)
 
-!  resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
-!  methode de Newton
-    
+    !  resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
+    !  methode de Newton
+
     zx0=xref+e_elips*sin(xref)
     DO iter=1,10
@@ -43,11 +43,11 @@
     zx0=zx0+zdx
     if(zanom.lt.0.) zx0=-zx0
-    
+
     ! zteta est la longitude solaire
-    
+
     zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
-    
+
     psollong=zteta-timeperi
-    
+
     IF(psollong.LT.0.) psollong=psollong+2.*pi
     IF(psollong.GT.2.*pi) psollong=psollong-2.*pi
@@ -55,5 +55,5 @@
     !   sorties eventuelles:
     !   ---------------------
-    
+
     IF (lwrite) THEN
        WRITELOG(*,*) 'day of the year  :',pday
@@ -61,7 +61,7 @@
        LOG_DBG('solarlong')
     ENDIF
-    
+
   END SUBROUTINE solarlong
-  
+
   SUBROUTINE iniorbit
     !=======================================================================
@@ -92,21 +92,21 @@
     !
     !=======================================================================
-    
+
     !-----------------------------------------------------------------------
-    
+
     !   Local:
     !   ------
     REAL zxref,zanom,zz,zx0,zdx
     INTEGER iter
-    
+
     !-----------------------------------------------------------------------
-    
+
     WRITELOG(*,*) 'Perihelie en Mkm  ',periheli
-    WRITELOG(*,*) 'Aphelise  en Mkm  ',aphelie 
+    WRITELOG(*,*) 'Aphelise  en Mkm  ',aphelie
     WRITELOG(*,*) 'obliquite en degres  :',obliquit
-    
+
     e_elips=(aphelie-periheli)/(periheli+aphelie)
     p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr
-    
+
     WRITELOG(*,*) 'e_elips',e_elips
     WRITELOG(*,*) 'p_elips',p_elips
@@ -115,33 +115,32 @@
     ! calcul de l'angle polaire et de la distance au soleil :
     ! -------------------------------------------------------
-    
+
     !  calcul de l'zanomalie moyenne
-    
+
     zz=(year_day-peri_day)/year_day
     zanom=2.*pi*(zz-nint(zz))
     zxref=abs(zanom)
     WRITELOG(*,*) 'zanom  ',zanom
-    
+
     !  resolution de l'equation horaire  zx0 - e * sin (zx0) = zxref
     !  methode de Newton
-    
+
     zx0=zxref+e_elips*sin(zxref)
     DO  iter=1,100
        zdx=-(zx0-e_elips*sin(zx0)-zxref)/(1.-e_elips*cos(zx0))
-!       if(abs(zdx).le.(1.e-12)) goto 120
        zx0=zx0+zdx
     END DO
-    
+
     zx0=zx0+zdx
     if(zanom.lt.0.) zx0=-zx0
     WRITELOG(*,*) 'zx0   ',zx0
-    
+
     ! zteta est la longitude solaire
-    
+
     timeperi=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
     WRITELOG(*,*) 'longitude solaire du perihelie timeperi = ',timeperi
-    
+
     LOG_INFO('iniorbit')
-    
+
   END SUBROUTINE iniorbit
 
@@ -171,5 +170,5 @@
     !   Declarations:
     !   -------------
-    
+
     ! arguments:
     ! ----------
@@ -177,17 +176,16 @@
     REAL, INTENT(IN) :: pls
     REAL, INTENT(OUT) :: pdist_sol,pdecli
-    
+
     !-----------------------------------------------------------------------
-    
+
     ! Distance Sun-Planet
-    
+
     pdist_sol=p_elips/(1.+e_elips*cos(pls+timeperi))
-    
+
     ! Solar declination
-    
+
     pdecli= asin (sin(pls)*sin(obliquit*pi/180.))
-    
+
   END SUBROUTINE orbite
 
 END MODULE astronomy
-
Index: /dynamico_lmdz/simple_physics/phyparam/physics/callkeys.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/callkeys.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/callkeys.F90	(revision 4229)
@@ -2,5 +2,5 @@
   IMPLICIT NONE
   SAVE
-  
+
   LOGICAL callrad,calldifv,calladj,callcond,callsoil,season,diurnal,lverbose
   INTEGER iradia
Index: /dynamico_lmdz/simple_physics/phyparam/physics/comgeomfi.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/comgeomfi.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/comgeomfi.F90	(revision 4229)
@@ -5,9 +5,9 @@
   REAL, ALLOCATABLE :: long(:), lati(:), sinlon(:), coslon(:), sinlat(:), coslat(:)
   INTEGER :: ngridmax, nlayermx, nsoilmx
-!$OMP THREADPRIVATE(long,lati,sinlon,coslon,sinlat,coslat,totarea)
-!$OMP THREADPRIVATE(ngridmax,nlayermx,nsoilmx)
+  !$OMP THREADPRIVATE(long,lati,sinlon,coslon,sinlat,coslat,totarea)
+  !$OMP THREADPRIVATE(ngridmax,nlayermx,nsoilmx)
 
 CONTAINS
-    
+
   SUBROUTINE init_comgeomfi(klon, klev, longitude, latitude)
     INTEGER, INTENT(IN) :: klon, klev
@@ -29,4 +29,4 @@
     coslon(:)=cos(long(:))
   END SUBROUTINE init_comgeomfi
-  
+
 END MODULE comgeomfi
Index: /dynamico_lmdz/simple_physics/phyparam/physics/convection.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/convection.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/convection.F90	(revision 4229)
@@ -76,5 +76,5 @@
           ! * momentum u,v is mixed with weight zalpha, i.e. u:=zalpha*mean(u)+(1-zalpha)*u
           ! where zalpha = sum( abs(theta-mean(theta))*dsig) / mean(theta)*sum(dsig)
-          ! for large deviations of theta from its mean, this formula could in principe yield zalpha>1. 
+          ! for large deviations of theta from its mean, this formula could in principe yield zalpha>1.
           zalpha=0.
           zum=0.
@@ -84,5 +84,5 @@
              zh2(i, l) = zhm
              ! we must use compute zum, zvm from zu2, zv2 to conserve momentum (see below)
-             zum=zum+dsig(l)*zu2(i,l) 
+             zum=zum+dsig(l)*zu2(i,l)
              zvm=zvm+dsig(l)*zv2(i,l)
           END DO
@@ -101,9 +101,9 @@
              zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
              zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
-          ENDDO          
+          ENDDO
        END IF
-       
+
     END DO
-    
+
   END SUBROUTINE adjust_onecolumn
 
@@ -113,5 +113,5 @@
        &                  pdufi,pdvfi,pdhfi,     &
        &                  pduadj,pdvadj,pdhadj)
-    
+
     !=======================================================================
     !
@@ -124,5 +124,5 @@
     !
     !=======================================================================
-    
+
     INTEGER, INTENT(IN) :: ngrid,nlay
     REAL, INTENT(IN) ::  ptimestep
@@ -132,12 +132,12 @@
     REAL, INTENT(OUT) :: pdhadj(ngrid,nlay), pduadj(ngrid,nlay), pdvadj(ngrid,nlay)
 
-    !   local:    
+    !   local:
     REAL sig(nlay+1),sdsig(nlay),dsig(nlay)
     REAL zu(ngrid,nlay), zv(ngrid,nlay), zh(ngrid,nlay)
-    REAL zu2(ngrid,nlay), zv2(ngrid,nlay), zh2(ngrid,nlay)    
+    REAL zu2(ngrid,nlay), zv2(ngrid,nlay), zh2(ngrid,nlay)
     LOGICAL vtest(ngrid)
     INTEGER jadrs(ngrid)
     INTEGER jcnt, ig, l, jj
-    
+
     ! Add physics tendencies to u,v,h
     DO l=1,nlay
@@ -151,5 +151,5 @@
     zv2(:,:)=zv(:,:)
     zh2(:,:)=zh(:,:)
-    
+
     !-----------------------------------------------------------------------
     !   detection des profils a modifier:
@@ -160,5 +160,5 @@
     !   sinon, il reste a sa valeur initiale (.FALSE.)
     !   cette operation est vectorisable
-    
+
     vtest(:)=.FALSE.
     DO l=2,nlay
@@ -186,5 +186,5 @@
        END DO
     END DO
-    
+
   END SUBROUTINE convadj
 
Index: /dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90	(revision 4229)
@@ -9,5 +9,5 @@
 
 CONTAINS
-  
+
   SUBROUTINE iniphyparam(ngrid,nlayer, &
        &           punjours,              &
@@ -18,5 +18,5 @@
     USE planet, ONLY : coefir, coefvis
     USE astronomy
-    USE turbulence, ONLY : lmixmin, emin_turb 
+    USE turbulence, ONLY : lmixmin, emin_turb
     USE surface
     USE read_param_mod
@@ -27,5 +27,5 @@
          pdayref        ! Day of reference for the simulation
     REAL, INTENT(IN)  :: ptimestep, prad, pg, pr, pcpp, punjours
-    
+
     CALL read_param('unjours', 86400., unjours,'unjours')
     CALL read_param('planet_rad',prad,planet_rad,'planet_rad')
@@ -53,5 +53,5 @@
     CALL read_param('coefvis',.99     ,coefvis,'coefvis')
     CALL read_param('coefir',.08      ,coefir,'coefir')
-    
+
     CALL read_param('callrad',.true.,callrad,'appel rayonnemen')
     CALL read_param('calldifv',.true.,calldifv,'appel difv')
@@ -63,5 +63,5 @@
     CALL read_param('lverbose',.true.,lverbose,'appel soil')
     CALL read_param('period_sort',1.,period_sort,'period sorties en jour')
-    
+
     CALL check_mismatch('unjours', punjours, unjours)
     CALL check_mismatch('rad', prad, planet_rad)
@@ -69,5 +69,5 @@
     CALL check_mismatch('R', pr, r)
     CALL check_mismatch('cpp', pcpp, cpp)
-    
+
     WRITELOG(*,*) 'Activation de la physique:'
     WRITELOG(*,*) ' R=',r
Index: /dynamico_lmdz/simple_physics/phyparam/physics/logging.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/logging.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/logging.F90	(revision 4229)
@@ -16,5 +16,5 @@
      END SUBROUTINE plugin
 
-     ! Plugin that writes into string 'line' information about the gridpoint of index 'index' 
+     ! Plugin that writes into string 'line' information about the gridpoint of index 'index'
      SUBROUTINE plugin_log_gridpoint(index, line)
        INTEGER, INTENT(IN) :: index ! index of gridpoint
@@ -25,5 +25,5 @@
 
 
-  ! This module provides a default implementations of plugins but the top-level driver is welcome to override them.  
+  ! This module provides a default implementations of plugins but the top-level driver is welcome to override them.
   ! Note F2003/F2008: pgfortran (F2003) accepts to initialize pointers only to NULL()
   ! => plugins are initialzed to NULL() and set to default values in flush_log and log_gridpoint
@@ -39,5 +39,5 @@
   INTEGER :: logging_lineno=0
 
-  ! messages with a log level > max_log_level are not printed 
+  ! messages with a log level > max_log_level are not printed
   INTEGER, PARAMETER, PUBLIC :: log_level_fatal=1, log_level_error=2, log_level_warn=3, log_level_info=4, log_level_dbg=5
   INTEGER, PUBLIC :: max_log_level = log_level_info
@@ -81,5 +81,5 @@
     logging_lineno = logging_lineno+1
 #ifndef XCODEML
-    IF(.NOT.ASSOCIATED(log_gridpoint_plugin)) log_gridpoint_plugin => default_log_gridpoint 
+    IF(.NOT.ASSOCIATED(log_gridpoint_plugin)) log_gridpoint_plugin => default_log_gridpoint
     CALL log_gridpoint_plugin(index, logging_buf(logging_lineno))
 #endif
@@ -87,5 +87,5 @@
 
   SUBROUTINE default_log_gridpoint(index, line)
-    INTEGER, INTENT(IN) :: index ! index of gridpoint                                                                                                                                    
+    INTEGER, INTENT(IN) :: index ! index of gridpoint
     CHARACTER(*), INTENT(OUT) :: line
     line=''
Index: /dynamico_lmdz/simple_physics/phyparam/physics/mellor_yamada.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/mellor_yamada.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/mellor_yamada.F90	(revision 4229)
@@ -5,27 +5,27 @@
 CONTAINS
 
-  PURE SUBROUTINE my_25(imax,kmax,dt,gp,zi,z,u,v,teta,cd,q2,long,km,kh) 
-    
-    !********************************************************************** 
-    !****** FERMETURE MELLOR-YAMADA DE NIVEAU 2.5 (QUASI-EQUILIBRE) ******* 
-    !* q2 au interfaces entre mailles.                                      
-    !**********************************************************************
-    
-    INTEGER, INTENT(IN) :: imax,kmax 
+  PURE SUBROUTINE my_25(imax,kmax,dt,gp,zi,z,u,v,teta,cd,q2,long,km,kh)
+
+    !**********************************************************************
+    !****** FERMETURE MELLOR-YAMADA DE NIVEAU 2.5 (QUASI-EQUILIBRE) *******
+    !* q2 au interfaces entre mailles.
+    !**********************************************************************
+
+    INTEGER, INTENT(IN) :: imax,kmax
     REAL, INTENT(IN) :: dt, gp
     REAL, INTENT(IN) :: z(imax,kmax), &
-         u(imax,kmax), v(imax,kmax), teta(imax,kmax), cd(imax)         
+         u(imax,kmax), v(imax,kmax), teta(imax,kmax), cd(imax)
     REAL, INTENT(INOUT) :: zi(imax,kmax+1), q2(imax,kmax+1)
     REAL, INTENT(OUT) :: long(imax,kmax+1), km(imax,kmax+1), kh(imax,kmax+1)
 
-    !************* DECLARATIONS ******************************************* 
-   
-    INTEGER, PARAMETER :: impmax=5 
+    !************* DECLARATIONS *******************************************
+
+    INTEGER, PARAMETER :: impmax=5
     REAL, PARAMETER :: kappa=0.4, &
          a1=0.92, a2=0.74, b1=16.6, b2=10.1, c1=0.08,           &
-         e1=1.8, e2=1.33, &                                       
+         e1=1.8, e2=1.33, &
          khmin=1.0e-5, kmmin=1.0e-5, kqmin=1.0e-5,              &
-         q2min=0.001, q2lmin=0.001, &                             
-         ghmax=0.023, ghmin=-0.28 
+         q2min=0.001, q2lmin=0.001, &
+         ghmax=0.023, ghmin=-0.28
 
     REAL longc, au, ateta, az, adz, akq, acd, &
@@ -33,221 +33,221 @@
          am, azi, bet, beta, dest, du, dv, gh, &
          prod, q2c, sh, sm, sq, us, vt, vt1, vt2
-         
-    REAL unsdz(imax,kmax),unsdzi(imax,kmax+1) 
+
+    REAL unsdz(imax,kmax),unsdzi(imax,kmax+1)
     REAL kq(imax,kmax),                                               &
-         &     m2(imax,kmax+1),n2(imax,kmax+1),ri(imax,kmax+1)              
+         &     m2(imax,kmax+1),n2(imax,kmax+1),ri(imax,kmax+1)
     REAL a(imax,kmax),b(imax,kmax),c(imax,kmax),f(imax,kmax),         &
-         &     alph(imax,kmax)                                              
-    REAL ksdz2inf,ksdz2sup 
+         &     alph(imax,kmax)
+    REAL ksdz2inf,ksdz2sup
 
     INTEGER :: i,k
-    
-    !************* INCREMENTS VERTICAUX *********************************** 
-    
-    DO i=1,imax 
-       zi(i,kmax+1)=zi(i,kmax)+2.0*(z(i,kmax)-zi(i,kmax)) 
-    END DO
-    DO k=1,kmax 
-       DO  i=1,imax 
-          unsdz(i,k)=1.0/(zi(i,k+1)-zi(i,k)) 
-       END DO
-    END DO
-
-    DO k=2,kmax 
-       DO i=1,imax 
-          unsdzi(i,k)=1.0/(z(i,k)-z(i,k-1)) 
+
+    !************* INCREMENTS VERTICAUX ***********************************
+
+    DO i=1,imax
+       zi(i,kmax+1)=zi(i,kmax)+2.0*(z(i,kmax)-zi(i,kmax))
+    END DO
+    DO k=1,kmax
+       DO  i=1,imax
+          unsdz(i,k)=1.0/(zi(i,k+1)-zi(i,k))
+       END DO
+    END DO
+
+    DO k=2,kmax
+       DO i=1,imax
+          unsdzi(i,k)=1.0/(z(i,k)-z(i,k-1))
        END do
     END DO
 
-    DO i=1,imax 
-       unsdzi(i,1)=0.5/(z(i,1)-zi(i,1)) 
-       unsdzi(i,kmax+1)=0.5/(zi(i,kmax+1)-z(i,kmax)) 
+    DO i=1,imax
+       unsdzi(i,1)=0.5/(z(i,1)-zi(i,1))
+       unsdzi(i,kmax+1)=0.5/(zi(i,kmax+1)-z(i,kmax))
     END do
-    
-    !********************************************************************** 
-                      
-    
-    !************* DIFFUSIVITES KH, KM et KQ ****************************** 
-    ! Ci-dessous, une premiere estimation des diffusivites turbulentes km * 
-    ! et kh est effectuee pour utilisation dans les taux de production    * 
-    ! et de destruction de q2 et q2l. On calcule aussi kq.                * 
-    
-    DO k=2,kmax 
-       DO i=1,imax 
-          beta=2.0/(teta(i,k)+teta(i,k-1)) 
-          n2(i,k)=beta*gp*unsdzi(i,k)*(teta(i,k)-teta(i,k-1)) 
-          n2(i,k)=amax1(0.0,n2(i,k)) 
-          du=unsdzi(i,k)*(u(i,k)-u(i,k-1)) 
-          dv=unsdzi(i,k)*(v(i,k)-v(i,k-1)) 
-          m2(i,k)=du*du+dv*dv 
-          ri(i,k)=n2(i,k)/(m2(i,k)+1.0e-10) 
-          ri(i,k)=amax1(-0.1,min(4.0,ri(i,k))) 
-       END DO
-    END DO
-
-    DO k=2,kmax 
-       DO i=1,imax 
-          vt=kappa*(zi(i,k)-zi(i,1)) 
-          long(i,k)=vt/(1.0+vt/160.0) 
-          IF(n2(i,k).gt.0.0) THEN 
-             long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k)) 
+
+    !**********************************************************************
+
+
+    !************* DIFFUSIVITES KH, KM et KQ ******************************
+    ! Ci-dessous, une premiere estimation des diffusivites turbulentes km *
+    ! et kh est effectuee pour utilisation dans les taux de production    *
+    ! et de destruction de q2 et q2l. On calcule aussi kq.                *
+
+    DO k=2,kmax
+       DO i=1,imax
+          beta=2.0/(teta(i,k)+teta(i,k-1))
+          n2(i,k)=beta*gp*unsdzi(i,k)*(teta(i,k)-teta(i,k-1))
+          n2(i,k)=amax1(0.0,n2(i,k))
+          du=unsdzi(i,k)*(u(i,k)-u(i,k-1))
+          dv=unsdzi(i,k)*(v(i,k)-v(i,k-1))
+          m2(i,k)=du*du+dv*dv
+          ri(i,k)=n2(i,k)/(m2(i,k)+1.0e-10)
+          ri(i,k)=amax1(-0.1,min(4.0,ri(i,k)))
+       END DO
+    END DO
+
+    DO k=2,kmax
+       DO i=1,imax
+          vt=kappa*(zi(i,k)-zi(i,1))
+          long(i,k)=vt/(1.0+vt/160.0)
+          IF(n2(i,k).gt.0.0) THEN
+             long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k))
           END IF
           gh=amax1(ghmin,                                                 &
-               &     min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))             
+               &     min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))
           sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh*                          &
                &           ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/     &
-               &        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))          
-          km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm) 
-          sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2)) 
-          kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh) 
-       END DO
-    END DO
-
-    DO i=1,imax 
-       us=sqrt(cd(i)*(u(i,1)*u(i,1)+v(i,1)*v(i,1))) 
-       vt1=(b1**0.666667)*us*us 
+               &        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))
+          km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm)
+          sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2))
+          kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh)
+       END DO
+    END DO
+
+    DO i=1,imax
+       us=sqrt(cd(i)*(u(i,1)*u(i,1)+v(i,1)*v(i,1)))
+       vt1=(b1**0.666667)*us*us
        vt2=(b1**0.6666667)*kappa*kappa*                                 &
-            &     m2(i,2)*(zi(i,2)-zi(i,1))*(zi(i,2)-zi(i,1))                  
-       !       q2(i,1)=amax1(q2min,vt1,vt2)                                    
-       q2(i,1)=amax1(q2min,vt1) 
-       long(i,1)=0.0 
-       long(i,kmax+1)=long(i,kmax) 
-       sq=0.2 
-       kq(i,1)=amax1(kqmin,kappa*(z(i,1)-zi(i,1))*us*sq) 
-    END DO
-    
-    DO k=2,kmax 
-       DO i=1,imax 
-          longc=0.5*(long(i,k)+long(i,k+1)) 
-          q2c=0.5*(q2(i,k)+q2(i,k+1)) 
-          sq=0.2 
-          kq(i,k)=amax1(kqmin,longc*sqrt(q2c)*sq) 
-       END DO
-    END DO
-
-    !********************************************************************** 
-                                                                        
-    
-    !************* CALCUL DE Q2 ******************************************* 
-    
-    DO k=2,kmax 
-       DO i=1,imax 
-          prod=2.0*(km(i,k)*m2(i,k)+amax1(0.0,-kh(i,k)*n2(i,k))) 
+            &     m2(i,2)*(zi(i,2)-zi(i,1))*(zi(i,2)-zi(i,1))
+       !       q2(i,1)=amax1(q2min,vt1,vt2)
+       q2(i,1)=amax1(q2min,vt1)
+       long(i,1)=0.0
+       long(i,kmax+1)=long(i,kmax)
+       sq=0.2
+       kq(i,1)=amax1(kqmin,kappa*(z(i,1)-zi(i,1))*us*sq)
+    END DO
+
+    DO k=2,kmax
+       DO i=1,imax
+          longc=0.5*(long(i,k)+long(i,k+1))
+          q2c=0.5*(q2(i,k)+q2(i,k+1))
+          sq=0.2
+          kq(i,k)=amax1(kqmin,longc*sqrt(q2c)*sq)
+       END DO
+    END DO
+
+    !**********************************************************************
+
+
+    !************* CALCUL DE Q2 *******************************************
+
+    DO k=2,kmax
+       DO i=1,imax
+          prod=2.0*(km(i,k)*m2(i,k)+amax1(0.0,-kh(i,k)*n2(i,k)))
           dest=2.0*(amax1(0.0,kh(i,k)*n2(i,k))                            &
-               &           +q2(i,k)*sqrt(q2(i,k))/(b1*long(i,k)))                 
-          IF(k.lt.kmax) THEN 
-             ksdz2sup=unsdzi(i,k)*unsdz(i,k)*kq(i,k) 
-          ELSE 
-             ksdz2sup=0.0 
+               &           +q2(i,k)*sqrt(q2(i,k))/(b1*long(i,k)))
+          IF(k.lt.kmax) THEN
+             ksdz2sup=unsdzi(i,k)*unsdz(i,k)*kq(i,k)
+          ELSE
+             ksdz2sup=0.0
           END IF
-          ksdz2inf=unsdzi(i,k)*unsdz(i,k-1)*kq(i,k-1) 
-          b(i,k)=-ksdz2inf*dt 
-          a(i,k)=1.0+dt*(dest/q2(i,k)+ksdz2inf+ksdz2sup) 
-          c(i,k)=-ksdz2sup*dt 
-          f(i,k)=q2(i,k)+dt*prod 
-       END DO
-    END DO
-    DO i=1,imax 
+          ksdz2inf=unsdzi(i,k)*unsdz(i,k-1)*kq(i,k-1)
+          b(i,k)=-ksdz2inf*dt
+          a(i,k)=1.0+dt*(dest/q2(i,k)+ksdz2inf+ksdz2sup)
+          c(i,k)=-ksdz2sup*dt
+          f(i,k)=q2(i,k)+dt*prod
+       END DO
+    END DO
+    DO i=1,imax
        f(i,2)=f(i,2)                                                    &
-            &       +dt*unsdzi(i,2)*unsdz(i,1)*kq(i,1)*q2(i,1)                 
-    END DO
-    
-    DO i=1,imax 
-       alph(i,2)=a(i,2) 
-    END DO
-    
-    DO k=3,kmax 
-       DO i=1,imax 
-          bet=b(i,k)/alph(i,k-1) 
-          alph(i,k)=a(i,k)-bet*c(i,k-1) 
-          f(i,k)=f(i,k)-bet*f(i,k-1) 
-       END DO
-    END DO
-    
-    DO i=1,imax 
-       q2(i,kmax)=amax1(q2min,f(i,kmax)/alph(i,kmax)) 
-       q2(i,kmax+1)=q2(i,kmax) 
-    END DO
-    
-    DO k=kmax-1,2,-1 
-       DO i=1,imax 
-          q2(i,k)=amax1(q2min,(f(i,k)-c(i,k)*q2(i,k+1))/alph(i,k)) 
-       END DO
-    END DO
-    
-    DO i=1,imax 
-       q2(i,2)=amax1(q2(i,2),q2(i,1)) 
-    END DO
-    
-!********************************************************************** 
-                                                                        
-                                                                        
-!************* EVALUATION FINALE DE KH ET KM ************************** 
-                                                                        
-    DO k=2,kmax 
-       DO i=1,imax 
-          IF(n2(i,k).gt.0.0) THEN 
-             long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k))	 
+            &       +dt*unsdzi(i,2)*unsdz(i,1)*kq(i,1)*q2(i,1)
+    END DO
+
+    DO i=1,imax
+       alph(i,2)=a(i,2)
+    END DO
+
+    DO k=3,kmax
+       DO i=1,imax
+          bet=b(i,k)/alph(i,k-1)
+          alph(i,k)=a(i,k)-bet*c(i,k-1)
+          f(i,k)=f(i,k)-bet*f(i,k-1)
+       END DO
+    END DO
+
+    DO i=1,imax
+       q2(i,kmax)=amax1(q2min,f(i,kmax)/alph(i,kmax))
+       q2(i,kmax+1)=q2(i,kmax)
+    END DO
+
+    DO k=kmax-1,2,-1
+       DO i=1,imax
+          q2(i,k)=amax1(q2min,(f(i,k)-c(i,k)*q2(i,k+1))/alph(i,k))
+       END DO
+    END DO
+
+    DO i=1,imax
+       q2(i,2)=amax1(q2(i,2),q2(i,1))
+    END DO
+
+    !**********************************************************************
+
+
+    !************* EVALUATION FINALE DE KH ET KM **************************
+
+    DO k=2,kmax
+       DO i=1,imax
+          IF(n2(i,k).gt.0.0) THEN
+             long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k))
           END IF
           gh=amax1(ghmin,                                                 &
-               &         min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))         
+               &         min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))
           sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh*                          &
                &           ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/     &
-               &        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))          
-          km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm) 
-          sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2)) 
-          kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh) 
-       END DO
-    END DO
-    
-    DO i=1,imax 
-       km(i,1)=kmmin 
-       km(i,kmax+1)=km(i,kmax) 
-       kh(i,1)=khmin 
-       kh(i,kmax+1)=kh(i,kmax) 
-    END DO
-    
-    !********************************************************************** 
-    
-    am=1.0/float(imax) 
-    DO k=kmax,1,-1 
-       au=0.0 
-       ateta=0.0 
-       az=0.0 
-       adz=0.0 
-       akq=0.0 
-       acd=0.0 
-       DO i=1,imax 
-          au=au+am*sqrt(u(i,k)*u(i,k)+v(i,k)*v(i,k)) 
-          ateta=ateta+am*teta(i,k) 
-          az=az+am*z(i,k) 
-          adz=adz+am*(1.0/unsdz(i,k)) 
-          akq=akq+am*kq(i,k) 
-          acd=acd+am*cd(i) 
-       END DO
-    END DO
-    
-    DO k=kmax+1,1,-1 
-       azi=0.0 
-       adzi=0.0 
-       aq2=0.0 
-       al=0.0 
-       akm=0.0 
-       akh=0.0 
-       am2=0.0 
-       al0=0.0 
-       DO i=1,imax 
-          azi=azi+am*zi(i,k) 
-          adzi=adzi+am*(1.0/unsdzi(i,k)) 
-          aq2=aq2+am*q2(i,k) 
-          al=al+am*long(i,k) 
-          akm=akm+am*km(i,k) 
-          akh=akh+am*kh(i,k) 
-          am2=am2+am*m2(i,k) 
-          !         al0=al0+am*long0d(i)                                          
-       END DO
-    END DO
-    
+               &        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))
+          km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm)
+          sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2))
+          kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh)
+       END DO
+    END DO
+
+    DO i=1,imax
+       km(i,1)=kmmin
+       km(i,kmax+1)=km(i,kmax)
+       kh(i,1)=khmin
+       kh(i,kmax+1)=kh(i,kmax)
+    END DO
+
+    !**********************************************************************
+
+    am=1.0/float(imax)
+    DO k=kmax,1,-1
+       au=0.0
+       ateta=0.0
+       az=0.0
+       adz=0.0
+       akq=0.0
+       acd=0.0
+       DO i=1,imax
+          au=au+am*sqrt(u(i,k)*u(i,k)+v(i,k)*v(i,k))
+          ateta=ateta+am*teta(i,k)
+          az=az+am*z(i,k)
+          adz=adz+am*(1.0/unsdz(i,k))
+          akq=akq+am*kq(i,k)
+          acd=acd+am*cd(i)
+       END DO
+    END DO
+
+    DO k=kmax+1,1,-1
+       azi=0.0
+       adzi=0.0
+       aq2=0.0
+       al=0.0
+       akm=0.0
+       akh=0.0
+       am2=0.0
+       al0=0.0
+       DO i=1,imax
+          azi=azi+am*zi(i,k)
+          adzi=adzi+am*(1.0/unsdzi(i,k))
+          aq2=aq2+am*q2(i,k)
+          al=al+am*long(i,k)
+          akm=akm+am*km(i,k)
+          akh=akh+am*kh(i,k)
+          am2=am2+am*m2(i,k)
+          !         al0=al0+am*long0d(i)
+       END DO
+    END DO
+
   END SUBROUTINE my_25
-  
+
 END MODULE mellor_yamada
Index: /dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90	(revision 4229)
@@ -12,5 +12,5 @@
        ps_rad=1.e5, height_scale=10000., ref_temp=285., &
        capcal_nosoil=1e5, tsoil_init=300.
-    
+
   ! internal state, written to / read from disk at checkpoint / restart
   REAL, ALLOCATABLE :: tsurf(:), tsoil(:,:), capcal(:), fluxgrd(:)
@@ -19,5 +19,5 @@
   ! precomputed arrays
   REAL, ALLOCATABLE :: rnatur(:), albedo(:),emissiv(:), z0(:), inertie(:)
-  !$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie)          
+  !$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie)
 
   INTEGER :: icount
@@ -42,5 +42,5 @@
 
     !=======================================================================
-    !   Top routine of the physical parametrisations of the LMD 
+    !   Top routine of the physical parametrisations of the LMD
     !   20 parameters GCM for planetary atmospheres.
     !   It includes:
@@ -56,5 +56,5 @@
          nlayer                   ! Number of vertical layers.
     LOGICAL, INTENT(IN) ::      &
-         firstcall,             & ! True at the first call 
+         firstcall,             & ! True at the first call
          lastcall                 ! True at the last call
     REAL, INTENT(IN)    ::      &
@@ -73,5 +73,5 @@
          pdt(ngrid,nlayer),     &
          pdpsrf(ngrid)
-    
+
     !    Local variables :
     REAL, DIMENSION(ngrid) :: mu0
@@ -90,5 +90,5 @@
     REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer)
     REAL zfluxsw(ngrid),zfluxlw(ngrid), fluxrad(ngrid)
-    
+
     WRITELOG(*,*) 'latitude0', ngrid, lati(1:2), lati(ngrid-1:ngrid)
     WRITELOG(*,*) 'nlayer',nlayer
@@ -104,5 +104,5 @@
 
     nlevel=nlayer+1
-    igout=ngrid/2+1    
+    igout=ngrid/2+1
     zday=rjourvrai+gmtime
 
@@ -110,5 +110,5 @@
     !    0. Allocate and initialize at first call
     !    --------------------
-    
+
     IF(firstcall) THEN
        !         zday_last=rjourvrai
@@ -122,7 +122,7 @@
     !    1. Initialisations :
     !    --------------------
-    
+
     icount=icount+1
-    
+
     pdv(:,:)  = 0.
     pdu(:,:)  = 0.
@@ -131,9 +131,9 @@
     zflubid(:)= 0.
     zdtsrf(:) = 0.
-    
+
     !-----------------------------------------------------------------------
     !   calcul du geopotentiel aux niveaux intercouches
     !   ponderation des altitudes au niveau des couches en dp/p
-    
+
     DO l=1,nlayer
        DO ig=1,ngrid
@@ -151,5 +151,5 @@
        ENDDO
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !   Transformation de la temperature en temperature potentielle
@@ -160,12 +160,12 @@
        ENDDO
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !    2. Compute radiative tendencies :
     !    ---------------------------------------
-    
+
     IF(callrad) CALL radiative_tendencies(ngrid, igout, nlayer, &
-            gmtime, ptimestep*float(iradia), zday, pplev, pplay, pt, &
-            pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, fluxrad, mu0)
+         gmtime, ptimestep*float(iradia), zday, pplev, pplay, pt, &
+         pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, fluxrad, mu0)
 
     !-----------------------------------------------------------------------
@@ -178,8 +178,8 @@
           zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
        ENDDO
-       
+
        zdum1(:,:)=0.
        zdum2(:,:)=0.
-       
+
        do l=1,nlayer
           do ig=1,ngrid
@@ -187,12 +187,12 @@
           enddo
        enddo
-       
+
        CALL vdif(ngrid,nlayer,zday,        &
-       &        ptimestep,capcal,z0,       &
-       &        pplay,pplev,zzlay,zzlev,   &
-       &        pu,pv,zh,tsurf,emissiv,    &
-       &        zdum1,zdum2,zdum3,zflubid, &
-       &        zdufr,zdvfr,zdhfr,zdtsrfr, &
-       &        lverbose)
+            &        ptimestep,capcal,z0,       &
+            &        pplay,pplev,zzlay,zzlev,   &
+            &        pu,pv,zh,tsurf,emissiv,    &
+            &        zdum1,zdum2,zdum3,zflubid, &
+            &        zdufr,zdvfr,zdhfr,zdtsrfr, &
+            &        lverbose)
 
        DO l=1,nlayer
@@ -203,13 +203,13 @@
           ENDDO
        ENDDO
-       
+
        DO ig=1,ngrid
           zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig)
        ENDDO
-       
+
     ELSE
        DO ig=1,ngrid
           zdtsrf(ig)=zdtsrf(ig)+ &
-          &      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
+               &      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
        ENDDO
     ENDIF
@@ -218,7 +218,7 @@
     !   4. Dry convective adjustment:
     !   -----------------------------
-    
+
     IF(calladj) THEN
-       
+
        DO l=1,nlayer
           DO ig=1,ngrid
@@ -226,15 +226,15 @@
           ENDDO
        ENDDO
-       
+
        zdufr(:,:)=0.
        zdvfr(:,:)=0.
        zdhfr(:,:)=0.
-       
+
        CALL convadj(ngrid,nlayer,ptimestep, &
-       &                pplay,pplev,zpopsk, &
-       &                pu,pv,zh,           &
-       &                pdu,pdv,zdum1,      &
-       &                zdufr,zdvfr,zdhfr)
-       
+            &                pplay,pplev,zpopsk, &
+            &                pu,pv,zh,           &
+            &                pdu,pdv,zdum1,      &
+            &                zdufr,zdvfr,zdhfr)
+
        DO l=1,nlayer
           DO ig=1,ngrid
@@ -244,5 +244,5 @@
           ENDDO
        ENDDO
-       
+
     ENDIF
 
@@ -250,18 +250,18 @@
     !   On ajoute les tendances physiques a la temperature du sol:
     !   ---------------------------------------------------------------
-    
+
     DO ig=1,ngrid
-       tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig) 
-    ENDDO
-    
+       tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)
+    ENDDO
+
     WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1)
-    
+
     !-----------------------------------------------------------------------
     !   soil temperatures:
     !   --------------------
-    
+
     IF (callsoil) THEN
        CALL soil(ngrid,nsoilmx,.false.,inertie, &
-       &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+            &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
        IF(lverbose) THEN
           WRITELOG(*,*) 'Surface Heat capacity,conduction Flux, Ts, dTs, dt'
@@ -271,5 +271,5 @@
        ENDIF
     ENDIF
-    
+
     !-----------------------------------------------------------------------
     !   sorties:
@@ -285,9 +285,9 @@
        zday_last=zday
        !  Ecriture/extension de la coordonnee temps
-       
+
        do ig=1,ngridmax
           zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(r*ref_temp))
        enddo
-       
+
        call writefield('u','Vent zonal moy','m/s',pu)
        call writefield('v','Vent meridien moy','m/s',pv)
@@ -295,5 +295,5 @@
        call writefield('geop','Geopotential','m2/s2',pphi)
        call writefield('plev','plev','Pa',pplev(:,1:nlayer))
-       
+
        call writefield('du','du',' ',pdu)
        call writefield('dv','du',' ',pdv)
@@ -301,5 +301,5 @@
        call writefield('dtsw','dtsw',' ',zdtsw)
        call writefield('dtlw','dtlw',' ',zdtlw)
-       
+
        call writefield('ts','Surface temper','K',tsurf)
        call writefield('coslon','coslon',' ',coslon)
@@ -313,7 +313,7 @@
        call writefield('swsurf','SW surf','Pa',zfluxsw)
        call writefield('lwsurf','LW surf','Pa',zfluxlw)
-       
+
     endif
-    
+
   END SUBROUTINE phyparam
 
@@ -338,5 +338,5 @@
          &                 zfluxsw(ngrid),      & ! short-wave flux at surface
          &                 fluxrad(ngrid),      & ! surface flux
-                           mu0(ngrid)             ! ??
+         mu0(ngrid)             ! ??
     ! local variables
     REAL :: fract(ngrid), dtrad(ngrid, nlayer)
@@ -347,8 +347,8 @@
     !    2.1 Insolation
     !    --------------------------------------------------
-    
+
     CALL solarlong(zday,zls)
     CALL orbite(zls,dist_sol,declin)
-    
+
     IF(diurnal) THEN
        IF ( .TRUE. ) then
@@ -362,5 +362,5 @@
           CALL zenang(ngrid,zls,gmtime,zdtime,lati,long,mu0,fract)
        ENDIF
-       
+
        IF(lverbose) THEN
           WRITELOG(*,*) 'day, declin, sinlon,coslon,sinlat,coslat'
@@ -375,10 +375,10 @@
        CALL mucorr(ngrid,declin,lati,mu0,fract,height_scale,planet_rad)
     ENDIF
-    
+
     zinsol=solarcst/(dist_sol*dist_sol)
-    
+
     !    2.2 Radiative tendencies and fluxes:
     !    --------------------------------------------------
-    
+
     CALL sw(ngrid,nlayer,diurnal,coefvis,albedo, &
          &              pplev,ps_rad,                 &
@@ -386,13 +386,13 @@
          &              zfluxsw,zdtsw,                &
          &              lverbose)
-    
+
     CALL lw(ngrid,nlayer,coefir,emissiv, &
          &             pplev,ps_rad,tsurf,pt, &
          &             zfluxlw,zdtlw,         &
          &             lverbose)
-    
+
     !    2.4 surface fluxes
     !    ------------------------------
-    
+
     DO ig=1,ngrid
        fluxrad(ig)=emissiv(ig)*zfluxlw(ig) &
@@ -402,8 +402,8 @@
        fluxrad(ig)=fluxrad(ig)-zplanck(ig)
     ENDDO
-    
+
     !    2.5 Temperature tendencies
     !    --------------------------
-    
+
     DO l=1,nlayer
        DO ig=1,ngrid
@@ -412,5 +412,5 @@
        ENDDO
     ENDDO
-    
+
     IF(lverbose) THEN
        WRITELOG(*,*) 'Diagnostics for radiation'
@@ -428,5 +428,5 @@
        LOG_DBG('radiative_tendencies')
     ENDIF
-    
+
   END SUBROUTINE radiative_tendencies
 
Index: /dynamico_lmdz/simple_physics/phyparam/physics/phys_const.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/phys_const.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/phys_const.F90	(revision 4229)
@@ -1,4 +1,4 @@
 MODULE phys_const
-  
+
   REAL planet_rad,g,r,cpp,rcp,dtphys,unjours,mugaz,omeg
 
Index: /dynamico_lmdz/simple_physics/phyparam/physics/planet.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/planet.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/planet.F90	(revision 4229)
@@ -4,4 +4,4 @@
 
   REAL :: coefvis, coefir
-  
+
 END MODULE planet
Index: /dynamico_lmdz/simple_physics/phyparam/physics/radiative_lw.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/radiative_lw.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/radiative_lw.F90	(revision 4229)
@@ -5,14 +5,14 @@
   IMPLICIT NONE
   SAVE
-  
+
   PRIVATE
-  
+
   PUBLIC :: lw
-  
+
   LOGICAL, PARAMETER :: lstrong=.TRUE.
   REAL, PARAMETER :: stephan=5.67e-08
-  
+
 CONTAINS
-  
+
   SUBROUTINE lw(ngrid,nlayer,coefir,emissiv, &
        pp,ps_rad,ptsurf,pt,              &
@@ -45,11 +45,11 @@
     !
     !=======================================================================
-    
+
     !   declarations:
     !   -------------
-    
+
     !   arguments:
     !   ----------
-    
+
     INTEGER, INTENT(IN)  ::  ngrid,nlayer
     REAL,    INTENT(IN)  :: coefir,emissiv(ngrid),ps_rad
@@ -57,8 +57,8 @@
     REAL,    INTENT(OUT) ::  pdtlw(ngrid,nlayer),pfluxir(ngrid)
     LOGICAL, INTENT(IN)  :: lwrite
-    
+
     !   variables locales:
     !   ------------------
-    
+
     INTEGER nlevel,ilev,ig,i,il
     REAL zplanck(ngrid,nlayer+1),zcoef
@@ -72,11 +72,11 @@
     !   initialisations:
     !   ----------------
-    
+
     nlevel=nlayer+1
-    
+
     !-----------------------------------------------------------------------
     !   2. calcul des quantites d'absorbants:
     !   -------------------------------------
-    
+
     !   absorption forte
     IF(lstrong) THEN
@@ -93,5 +93,5 @@
        ENDIF
        zcoef=-log(coefir)/sqrt(ps_rad*ps_rad/(2.*g))
-       
+
        !   absorption faible
     ELSE
@@ -103,10 +103,10 @@
        zcoef=-log(coefir)/ps_rad
     ENDIF
-    
-    
+
+
     !-----------------------------------------------------------------------
     !   2. calcul de la fonction de corps noir:
     !   ---------------------------------------
-    
+
     DO ilev=1,nlayer
        DO ig=1,ngrid
@@ -116,9 +116,9 @@
        ENDDO
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !   4. flux descendants:
     !   --------------------
-    
+
     DO ilev=1,nlayer
        DO ig=1,ngrid
@@ -129,5 +129,5 @@
        ENDDO
        CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
-       
+
        DO il=nlayer,ilev,-1
           zlwtr2(:)=zlwtr1(:)
@@ -142,10 +142,10 @@
        ENDDO
     ENDDO
-    
+
     DO ig=1,ngrid
        zfluxdn(ig,nlevel)=0.
        pfluxir(ig)=emissiv(ig)*zfluxdn(ig,1)
     ENDDO
-    
+
     DO ig=1,ngrid
        zfluxup(ig,1)=ptsurf(ig)*ptsurf(ig)
@@ -153,9 +153,9 @@
             +(1.-emissiv(ig))*zfluxdn(ig,1)
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !   3. flux montants:
     !   ------------------
-    
+
     DO ilev=1,nlayer
        DO ig=1,ngrid
@@ -177,5 +177,5 @@
           ENDDO
        ENDDO
-       
+
     ENDDO
 
@@ -183,5 +183,5 @@
     !   5. calcul des flux nets:
     !   ------------------------
-    
+
     DO ilev=1,nlevel
        DO ig=1,ngrid
@@ -189,9 +189,9 @@
        ENDDO
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !   6. Calcul des taux de refroidissement:
     !   --------------------------------------
-    
+
     DO ilev=1,nlayer
        DO ig=1,ngrid
@@ -200,13 +200,13 @@
        ENDDO
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !   10. sorties eventuelles:
     !   ------------------------
-    
-    IF (lwrite) THEN       
+
+    IF (lwrite) THEN
        WRITELOG(*,*) 'Diagnostique rayonnement thermique'
        WRITELOG(*,*) 'temperature     ', &
-          'flux montant    flux desc.     taux de refroid.'
+            'flux montant    flux desc.     taux de refroid.'
        i=ngrid/2+1
        WRITELOG(6,'(4e18.4)') ptsurf(i)
@@ -215,10 +215,10 @@
                zfluxup(i,ilev),zfluxdn(i,ilev),pdtlw(i,ilev)
        ENDDO
-       WRITELOG(6,'(4e18.4)') zfluxup(i,nlevel),zfluxdn(i,nlevel)       
+       WRITELOG(6,'(4e18.4)') zfluxup(i,nlevel),zfluxdn(i,nlevel)
        LOG_DBG(tag)
     ENDIF
 
-!-----------------------------------------------------------------------
-    
+    !-----------------------------------------------------------------------
+
   END SUBROUTINE lw
 
@@ -239,6 +239,6 @@
        ENDDO
     ENDIF
-    
+
   END SUBROUTINE lwtr
-  
+
 END MODULE radiative_lw
Index: /dynamico_lmdz/simple_physics/phyparam/physics/radiative_sw.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/radiative_sw.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/radiative_sw.F90	(revision 4229)
@@ -5,5 +5,5 @@
   IMPLICIT NONE
   SAVE
-  
+
   PRIVATE
 
@@ -11,6 +11,6 @@
 
 CONTAINS
-  
-  PURE SUBROUTINE monGATHER(n,a,b,index)      
+
+  PURE SUBROUTINE monGATHER(n,a,b,index)
     INTEGER, INTENT(IN) ::  n,index(n)
     REAL, INTENT(IN) :: b(n)
@@ -32,5 +32,5 @@
     END DO
   end subroutine monscatter
-  
+
   SUBROUTINE sw(ngrid,nlayer,ldiurn, &
        coefvis,albedo, &
@@ -41,5 +41,5 @@
     !=======================================================================
     !
-    !   Rayonnement solaire en atmosphere non diffusante avec un 
+    !   Rayonnement solaire en atmosphere non diffusante avec un
     !   coefficient d'absoprption gris.
     !
@@ -60,9 +60,9 @@
     REAL, INTENT(IN)    :: psolarf0
     REAL, INTENT(OUT)   :: fsrfvis(ngrid),dtsw(ngrid,nlayer)
-    
+
     REAL zalb(ngrid),zmu(ngrid),zfract(ngrid)
     REAL zplev(ngrid,nlayer+1)
     REAL zflux(ngrid),zdtsw(ngrid,nlayer)
-    
+
     INTEGER ig,l,nlevel,index(ngrid),ncount,igout
     REAL ztrdir(ngrid,nlayer+1),ztrref(ngrid,nlayer+1)
@@ -71,9 +71,9 @@
     REAL zu(ngrid,nlayer+1)
     REAL tau0
-    
+
     !-----------------------------------------------------------------------
     !   1. initialisations:
     !   -------------------
-    
+
     nlevel=nlayer+1
 
@@ -81,5 +81,5 @@
     !   Definitions des tableaux locaux pour les points ensoleilles:
     !   ------------------------------------------------------------
-    
+
     IF (ldiurn) THEN
        ncount=0
@@ -106,11 +106,11 @@
        zplev(:,:)=plevel(:,:)
     ENDIF
-    
+
     !-----------------------------------------------------------------------
     !   calcul des profondeurs optiques integres depuis p=0:
     !   ----------------------------------------------------
-    
+
     tau0=-.5*log(coefvis)
-    
+
     ! calcul de la partie homogene de l'opacite
     tau0=tau0/ps_rad
@@ -120,9 +120,9 @@
        ENDDO
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !   2. calcul de la transmission depuis le sommet de l'atmosphere:
     !   -----------------------------------------------------------
-    
+
     DO l=1,nlevel
        DO ig=1,ncount
@@ -130,8 +130,8 @@
        ENDDO
     ENDDO
-    
+
     IF (lwrite) THEN
        igout=ncount/2+1
-       WRITELOG(*,*) 
+       WRITELOG(*,*)
        WRITELOG(*,*) 'Diagnostique des transmission dans le spectre solaire'
        WRITELOG(*,*) 'zfract, zmu, zalb'
@@ -142,9 +142,9 @@
        ENDDO
     ENDIF
-    
+
     !-----------------------------------------------------------------------
     !   3. taux de chauffage, ray. solaire direct:
     !   ------------------------------------------
-    
+
     DO l=1,nlayer
        DO ig=1,ncount
@@ -155,5 +155,5 @@
     ENDDO
     IF (lwrite) THEN
-       WRITELOG(*,*) 
+       WRITELOG(*,*)
        WRITELOG(*,*) 'Diagnostique des taux de chauffage solaires:'
        WRITELOG(*,*) ' 1 taux de chauffage lie au ray. solaire  direct'
@@ -162,10 +162,10 @@
        ENDDO
     ENDIF
-    
-    
+
+
     !-----------------------------------------------------------------------
     !   4. calcul du flux solaire arrivant sur le sol:
     !   ----------------------------------------------
-    
+
     DO ig=1,ncount
        z1(ig)=zfract(ig)*zmu(ig)*psolarf0*ztrdir(ig,1)
@@ -174,15 +174,15 @@
     ENDDO
     IF (lwrite) THEN
-       WRITELOG(*,*) 
+       WRITELOG(*,*)
        WRITELOG(*,*) 'Diagnostique des taux de chauffage solaires:'
        WRITELOG(*,*) ' 2 flux solaire net incident sur le sol'
        WRITELOG(*,*) zflux(igout)
     ENDIF
-    
-    
+
+
     !-----------------------------------------------------------------------
     !   5.calcul des traansmissions depuis le sol, cas diffus:
     !   ------------------------------------------------------
-    
+
     DO l=1,nlevel
        DO ig=1,ncount
@@ -190,7 +190,7 @@
        ENDDO
     ENDDO
-    
-    IF (lwrite) THEN
-       WRITELOG(*,*) 
+
+    IF (lwrite) THEN
+       WRITELOG(*,*)
        WRITELOG(*,*) 'Diagnostique des taux de chauffage solaires'
        WRITELOG(*,*) ' 3 transmission avec les sol'
@@ -200,9 +200,9 @@
        ENDDO
     ENDIF
-    
-    !-----------------------------------------------------------------------
-    !   6.ajout a l'echauffement de la contribution du ray. sol. reflechit: 
+
+    !-----------------------------------------------------------------------
+    !   6.ajout a l'echauffement de la contribution du ray. sol. reflechit:
     !   -------------------------------------------------------------------
-    
+
     DO l=1,nlayer
        DO ig=1,ncount
@@ -212,11 +212,11 @@
        ENDDO
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !   10. sorites eventuelles:
     !   ------------------------
-    
-    IF (lwrite) THEN
-       WRITELOG(*,*) 
+
+    IF (lwrite) THEN
+       WRITELOG(*,*)
        WRITELOG(*,*) 'Diagnostique des taux de chauffage solaires:'
        WRITELOG(*,*) ' 3 taux de chauffage total'
@@ -225,5 +225,5 @@
        ENDDO
     ENDIF
-    
+
     IF (ldiurn) THEN
        fsrfvis(:)=0.
@@ -246,7 +246,6 @@
     !        call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir     ')
     LOG_INFO('rad_sw')
-    
+
   END SUBROUTINE sw
-  
+
 END MODULE radiative_sw
-
Index: /dynamico_lmdz/simple_physics/phyparam/physics/read_param_mod.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/read_param_mod.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/read_param_mod.F90	(revision 4229)
@@ -5,5 +5,5 @@
   SAVE
 
-  INTERFACE 
+  INTERFACE
 
      ! each of these plugins reads a parameter of a certain type and stores the value in 'val'
@@ -43,5 +43,5 @@
   INTERFACE read_param
      PROCEDURE read_paramr, read_parami, read_paramb
-  END INTERFACE read_param
+  END INTERFACE
 
   PUBLIC :: read_param
@@ -68,5 +68,5 @@
 #endif
   END SUBROUTINE read_parami
-  
+
   SUBROUTINE read_paramb(name, defval, val, comment)
     CHARACTER(*), INTENT(IN) :: name, comment
@@ -78,4 +78,4 @@
 #endif
   END SUBROUTINE read_paramb
-  
+
 END MODULE read_param_mod
Index: /dynamico_lmdz/simple_physics/phyparam/physics/solar.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/solar.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/solar.F90	(revision 4229)
@@ -6,5 +6,5 @@
   PRIVATE
 
-  REAL, PARAMETER :: pi_local = 4.0 * ATAN(1.0), deux_pi_local = 2.0 * pi_local 
+  REAL, PARAMETER :: pi_local = 4.0 * ATAN(1.0), deux_pi_local = 2.0 * pi_local
 
   PUBLIC :: solang, zenang, mucorr
@@ -14,5 +14,5 @@
   PURE SUBROUTINE solang ( kgrid,psilon,pcolon,psilat,pcolat, &
        ptim1,ptim2,ptim3,pmu0,pfract )
-    
+
     !-----------------------------------------------------------------------
     !          CALCULATES THE SOLAR ANGLE FOR ALL THE POINTS OF THE GRID
@@ -93,176 +93,176 @@
        ENDIF
     ENDDO
-    
+
   END SUBROUTINE solang
 
   SUBROUTINE zenang(klon,longi,gmtime,pdtrad,lat,long,              &
-       &                  pmu0,frac)                                      
-    USE astronomy 
-    
-    !=============================================================          
-    ! Auteur : O. Boucher (LMD/CNRS)                                        
-    !          d'apres les routines zenith et angle de Z.X. Li              
-    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal     
-    !          et l'ensoleillement moyen entre gmtime1 et gmtime2           
-    !          connaissant la declinaison, la latitude et la longitude.     
-    ! Rque   : Different de la routine angle en ce sens que zenang          
-    !          fournit des moyennes de pmu0 et non des valeurs              
-    !          instantanees, du coup frac prend toutes les valeurs          
-    !          entre 0 et 1.                                                
-    ! Date   : premiere version le 13 decembre 1994                         
-    !          revu pour  GCM  le 30 septembre 1996                         
-    !===============================================================        
-    ! longi----INPUT : la longitude vraie de la terre dans son plan         
-    !                  solaire a partir de l'equinoxe de printemps (degre)  
-    ! gmtime---INPUT : temps universel en fraction de jour                  
-    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)               
-    ! lat------INPUT : latitude en degres                                   
-    ! long-----INPUT : longitude en degres                                  
-    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad   
-    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad   
-    !================================================================       
-    integer klon 
-    !================================================================       
-    real longi, gmtime, pdtrad 
-    real lat(klon), long(klon), pmu0(klon), frac(klon) 
-    !================================================================       
-    integer i 
-    real gmtime1, gmtime2 
-    real incl 
-    real omega1, omega2, omega 
-    ! omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi.        
-    ! omega : heure en radian du coucher de soleil                          
-    ! -omega est donc l'heure en radian de lever du soleil                  
-    real omegadeb, omegafin 
-    real zfrac1, zfrac2, z1_mu, z2_mu 
-    ! declinaison en radian                     
-    real lat_sun 
-    ! longitude solaire en radian               
-    real lon_sun 
-    ! latitude du pt de grille en radian        
-    real latr 
-    !================================================================       
-    !                                                                       
-    !     incl=R_incl * pi_local / 180.                                     
+       &                  pmu0,frac)
+    USE astronomy
+
+    !=============================================================
+    ! Auteur : O. Boucher (LMD/CNRS)
+    !          d'apres les routines zenith et angle de Z.X. Li
+    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
+    !          et l'ensoleillement moyen entre gmtime1 et gmtime2
+    !          connaissant la declinaison, la latitude et la longitude.
+    ! Rque   : Different de la routine angle en ce sens que zenang
+    !          fournit des moyennes de pmu0 et non des valeurs
+    !          instantanees, du coup frac prend toutes les valeurs
+    !          entre 0 et 1.
+    ! Date   : premiere version le 13 decembre 1994
+    !          revu pour  GCM  le 30 septembre 1996
+    !===============================================================
+    ! longi----INPUT : la longitude vraie de la terre dans son plan
+    !                  solaire a partir de l'equinoxe de printemps (degre)
+    ! gmtime---INPUT : temps universel en fraction de jour
+    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
+    ! lat------INPUT : latitude en degres
+    ! long-----INPUT : longitude en degres
+    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
+    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
+    !================================================================
+    integer klon
+    !================================================================
+    real longi, gmtime, pdtrad
+    real lat(klon), long(klon), pmu0(klon), frac(klon)
+    !================================================================
+    integer i
+    real gmtime1, gmtime2
+    real incl
+    real omega1, omega2, omega
+    ! omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi.
+    ! omega : heure en radian du coucher de soleil
+    ! -omega est donc l'heure en radian de lever du soleil
+    real omegadeb, omegafin
+    real zfrac1, zfrac2, z1_mu, z2_mu
+    ! declinaison en radian
+    real lat_sun
+    ! longitude solaire en radian
+    real lon_sun
+    ! latitude du pt de grille en radian
+    real latr
+    !================================================================
+    !
+    !     incl=R_incl * pi_local / 180.
     WRITELOG(*,*) 'Obliquite =' ,obliquit
     LOG_INFO('solar')
 
-    incl=obliquit * pi_local / 180. 
-    !                                                                       
-    !     lon_sun = longi * pi_local / 180.0                                
-    lon_sun = longi 
-    lat_sun = ASIN (SIN(lon_sun)*SIN(incl) ) 
-    !                                                                       
-    gmtime1=gmtime*86400. 
-    gmtime2=gmtime*86400.+pdtrad 
-    !                                                                       
-    DO i = 1, klon 
-       !                                                                       
-       !     latr = lat(i) * pi_local / 180.                                   
-       latr = lat(i) 
-       !                                                                       
-       !--pose probleme quand lat=+/-90 degres                                 
-       !                                                                       
-       !      omega = -TAN(latr)*TAN(lat_sun)                                  
-       !      omega = ACOS(omega)                                              
-       !      IF (latr.GE.(pi_local/2.+lat_sun)                                
-       !     .    .OR. latr.LE.(-pi_local/2.+lat_sun)) THEN                    
-       !         omega = 0.0       ! nuit polaire                              
-       !      ENDIF                                                            
-       !      IF (latr.GE.(pi_local/2.-lat_sun)                                
-       !     .          .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN              
-       !         omega = pi_local  ! journee polaire                           
-       !      ENDIF                                                            
-       !                                                                       
-       !--remplace par cela (le cas par defaut est different)                  
-       !                                                                       
-       !--nuit polaire                                        
-       omega=0.0 
+    incl=obliquit * pi_local / 180.
+    !
+    !     lon_sun = longi * pi_local / 180.0
+    lon_sun = longi
+    lat_sun = ASIN (SIN(lon_sun)*SIN(incl) )
+    !
+    gmtime1=gmtime*86400.
+    gmtime2=gmtime*86400.+pdtrad
+    !
+    DO i = 1, klon
+       !
+       !     latr = lat(i) * pi_local / 180.
+       latr = lat(i)
+       !
+       !--pose probleme quand lat=+/-90 degres
+       !
+       !      omega = -TAN(latr)*TAN(lat_sun)
+       !      omega = ACOS(omega)
+       !      IF (latr.GE.(pi_local/2.+lat_sun)
+       !     .    .OR. latr.LE.(-pi_local/2.+lat_sun)) THEN
+       !         omega = 0.0       ! nuit polaire
+       !      ENDIF
+       !      IF (latr.GE.(pi_local/2.-lat_sun)
+       !     .          .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN
+       !         omega = pi_local  ! journee polaire
+       !      ENDIF
+       !
+       !--remplace par cela (le cas par defaut est different)
+       !
+       !--nuit polaire
+       omega=0.0
        IF (latr.GE.(pi_local/2.-lat_sun)                                 &
-            &          .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN               
-          ! journee polaire                            
-          omega = pi_local 
+            &          .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN
+          ! journee polaire
+          omega = pi_local
        ENDIF
        IF (latr.LT.(pi_local/2.+lat_sun).AND.                            &
             &    latr.GT.(-pi_local/2.+lat_sun).AND.                           &
             &    latr.LT.(pi_local/2.-lat_sun).AND.                            &
-            &    latr.GT.(-pi_local/2.-lat_sun)) THEN                          
-          omega = -TAN(latr)*TAN(lat_sun) 
-          omega = ACOS(omega) 
-       ENDIF
-       !                                                                       
-       omega1 = gmtime1 + long(i)*86400.0/360.0 
-       omega1 = omega1 / 86400.0*deux_pi_local 
-       omega1 = MOD (omega1+deux_pi_local, deux_pi_local) 
-       omega1 = omega1 - pi_local 
-       !                                                                       
-       omega2 = gmtime2 + long(i)*86400.0/360.0 
-       omega2 = omega2 / 86400.0*deux_pi_local 
-       omega2 = MOD (omega2+deux_pi_local, deux_pi_local) 
-       omega2 = omega2 - pi_local 
-       !                                                                       
-       !--on est dans la meme journee locale 
-       IF (omega1.LE.omega2) THEN 
-          !                                                                       
+            &    latr.GT.(-pi_local/2.-lat_sun)) THEN
+          omega = -TAN(latr)*TAN(lat_sun)
+          omega = ACOS(omega)
+       ENDIF
+       !
+       omega1 = gmtime1 + long(i)*86400.0/360.0
+       omega1 = omega1 / 86400.0*deux_pi_local
+       omega1 = MOD (omega1+deux_pi_local, deux_pi_local)
+       omega1 = omega1 - pi_local
+       !
+       omega2 = gmtime2 + long(i)*86400.0/360.0
+       omega2 = omega2 / 86400.0*deux_pi_local
+       omega2 = MOD (omega2+deux_pi_local, deux_pi_local)
+       omega2 = omega2 - pi_local
+       !
+       !--on est dans la meme journee locale
+       IF (omega1.LE.omega2) THEN
+          !
           IF (omega2.LE.-omega .OR. omega1.GE.omega .OR. omega.LT.1e-5)  &
-               THEN                                                           
-             !--nuit           
-             frac(i)=0.0 
-             pmu0(i)=0.0 
+               THEN
+             !--nuit
+             frac(i)=0.0
+             pmu0(i)=0.0
              !--jour+nuit/jou
-          ELSE 
-             omegadeb=MAX(-omega,omega1) 
-             omegafin=MIN(omega,omega2) 
-             frac(i)=(omegafin-omegadeb)/(omega2-omega1) 
+          ELSE
+             omegadeb=MAX(-omega,omega1)
+             omegafin=MIN(omega,omega2)
+             frac(i)=(omegafin-omegadeb)/(omega2-omega1)
              pmu0(i)=SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)*    &
-                  (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)          
+                  (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)
           ENDIF
-          !                                                                       
-          !---omega1 GT omega2 -- a cheval sur deux journees          
-       ELSE 
-          !                                                                       
-          !-------------------entre omega1 et pi                                  
-          !--nuit                                
-          IF (omega1.GE.omega) THEN 
-             zfrac1=0.0 
-             z1_mu =0.0 
-             !--jour+nuit                           
-          ELSE 
-             omegadeb=MAX(-omega,omega1) 
-             omegafin=omega 
-             zfrac1=omegafin-omegadeb 
+          !
+          !---omega1 GT omega2 -- a cheval sur deux journees
+       ELSE
+          !
+          !-------------------entre omega1 et pi
+          !--nuit
+          IF (omega1.GE.omega) THEN
+             zfrac1=0.0
+             z1_mu =0.0
+             !--jour+nuit
+          ELSE
+             omegadeb=MAX(-omega,omega1)
+             omegafin=omega
+             zfrac1=omegafin-omegadeb
              z1_mu =SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)*     &
-                  (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)          
+                  (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)
           ENDIF
-          !---------------------entre -pi et omega2                               
-          !--nuit                              
-          IF (omega2.LE.-omega) THEN 
-             zfrac2=0.0 
-             z2_mu =0.0 
-             !--jour+nuit                         
-          ELSE 
-             omegadeb=-omega 
-             omegafin=MIN(omega,omega2) 
-             zfrac2=omegafin-omegadeb 
+          !---------------------entre -pi et omega2
+          !--nuit
+          IF (omega2.LE.-omega) THEN
+             zfrac2=0.0
+             z2_mu =0.0
+             !--jour+nuit
+          ELSE
+             omegadeb=-omega
+             omegafin=MIN(omega,omega2)
+             zfrac2=omegafin-omegadeb
              z2_mu =SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)*     &
-                  (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)          
-             !                                                                       
+                  (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)
+             !
           ENDIF
-          !-----------------------moyenne                                         
-          frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi_local-omega1) 
-          pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2,1.E-10) 
-          !                                                                       
-          !---comparaison omega1 et omega2                          
-       ENDIF
-       !                                                                       
+          !-----------------------moyenne
+          frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi_local-omega1)
+          pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2,1.E-10)
+          !
+          !---comparaison omega1 et omega2
+       ENDIF
+       !
     ENDDO
-    !                                                                       
+    !
   END SUBROUTINE zenang
-  
+
   PURE SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad)
-    
+
     !=======================================================================
     !
-    !   Calcul of equivalent solar angle and and fraction of day whithout 
+    !   Calcul of equivalent solar angle and and fraction of day whithout
     !   diurnal cycle.
     !
@@ -284,10 +284,10 @@
     !
     !=======================================================================
-    
+
     !-----------------------------------------------------------------------
     !
     !    0. Declarations :
     !    -----------------
-    
+
     !     Arguments :
     !     -----------
@@ -306,11 +306,11 @@
 
     !-----------------------------------------------------------------------
-    
+
     z = pdeclin
     cz = cos (z)
     sz = sin (z)
-    
+
     DO j = 1, npts
-       
+
        phi = plat(j)
        cphi = cos(phi)
@@ -322,5 +322,5 @@
        a = 1.0 - t*t
        ap = a
-       
+
        IF(t.eq.0.) then
           t=0.5*pi
@@ -334,5 +334,5 @@
           ENDIF
        ENDIF
-       
+
        pmu(j) = (sphi*sz*t) / pi + b*sin(t)/pi
        pfract(j) = t / pi
@@ -344,11 +344,11 @@
        pmu(j) = pmu(j) / pfract(j)
        IF (pmu(j).eq.0.) pfract(j) = 0.
-       
+
     END DO
-    
+
     !-----------------------------------------------------------------------
     !   correction de rotondite:
     !   ------------------------
-    
+
     alph=phaut/prad
     DO j=1,npts
@@ -357,5 +357,5 @@
        !    $          (sqrt(alph*alph*pmu(j)*pmu(j)+2.*alph+1.)-alph*pmu(j))
     END DO
-    
+
   END SUBROUTINE mucorr
 
Index: /dynamico_lmdz/simple_physics/phyparam/physics/surface.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/surface.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/surface.F90	(revision 4229)
@@ -14,10 +14,10 @@
 
   ! precomputed variables
-  REAL :: lambda 
+  REAL :: lambda
   REAL,ALLOCATABLE :: dz1(:),dz2(:)
   !$OMP THREADPRIVATE(dz1,dz2,zc,lambda)
 
   ! internal state variables needed for checkpoint/restart
-  REAL, ALLOCATABLE :: zc(:,:),zd(:,:) 
+  REAL, ALLOCATABLE :: zc(:,:),zd(:,:)
   !$OMP THREADPRIVATE(zc,zd)
 
@@ -28,167 +28,167 @@
   SUBROUTINE init_soil(ngrid,nsoil)
     INTEGER, INTENT(IN) :: ngrid, nsoil
-    REAL min_period,dalph_soil 
-    REAL fz,rk,fz1,rk1,rk2 
+    REAL min_period,dalph_soil
+    REAL fz,rk,fz1,rk1,rk2
     INTEGER :: jk
 
     ! this is a function definition
-    fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.) 
-    
-    !-----------------------------------------------------------------------
-    !   ground levels                                                       
-    !   grnd=z/l where l is the skin depth of the diurnal cycle:            
-    !   --------------------------------------------------------            
-    
+    fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
+
+    !-----------------------------------------------------------------------
+    !   ground levels
+    !   grnd=z/l where l is the skin depth of the diurnal cycle:
+    !   --------------------------------------------------------
+
     WRITELOG(*,*) 'nsoil,ngrid,firstcall=',nsoil,ngrid, .TRUE.
 
-    ALLOCATE(dz1(nsoil),dz2(nsoil)) 
-    
-    min_period=20000. 
-    dalph_soil=2. 
-    
-    !   la premiere couche represente un dixieme de cycle diurne            
-    fz1=sqrt(min_period/pi) 
-    
-    DO jk=1,nsoil 
-       rk1=jk 
-       rk2=jk-1 
-       dz2(jk)=fz(rk1)-fz(rk2) 
-    ENDDO
-    DO jk=1,nsoil-1 
-       rk1=jk+.5 
-       rk2=jk-.5 
-       dz1(jk)=1./(fz(rk1)-fz(rk2)) 
-    ENDDO
-    lambda=fz(.5)*dz1(1) 
-
-    WRITELOG(*,*) 'full layers, intermediate layers (secoonds)' 
-    DO jk=1,nsoil 
-       rk=jk 
-       rk1=jk+.5 
-       rk2=jk-.5 
+    ALLOCATE(dz1(nsoil),dz2(nsoil))
+
+    min_period=20000.
+    dalph_soil=2.
+
+    !   la premiere couche represente un dixieme de cycle diurne
+    fz1=sqrt(min_period/pi)
+
+    DO jk=1,nsoil
+       rk1=jk
+       rk2=jk-1
+       dz2(jk)=fz(rk1)-fz(rk2)
+    ENDDO
+    DO jk=1,nsoil-1
+       rk1=jk+.5
+       rk2=jk-.5
+       dz1(jk)=1./(fz(rk1)-fz(rk2))
+    ENDDO
+    lambda=fz(.5)*dz1(1)
+
+    WRITELOG(*,*) 'full layers, intermediate layers (secoonds)'
+    DO jk=1,nsoil
+       rk=jk
+       rk1=jk+.5
+       rk2=jk-.5
        WRITELOG(*,*) fz(rk1)*fz(rk2)*pi,        &
-            &        fz(rk)*fz(rk)*pi                                          
+            &        fz(rk)*fz(rk)*pi
     ENDDO
     LOG_INFO('init_soil')
 
   END SUBROUTINE init_soil
-  
+
   SUBROUTINE soil(ngrid,nsoil,firstcall,ptherm_i,          &
        &          ptimestep,ptsrf,ptsoil,                  &
-       &          pcapcal,pfluxgrd)                                       
-    
+       &          pcapcal,pfluxgrd)
+
     !=======================================================================
-    !                                                                       
-    !   Auteur:  Frederic Hourdin     30/01/92                              
-    !   -------                                                             
-    !                                                                       
-    !   objet:  computation of : the soil temperature evolution             
-    !   ------                   the surfacic heat capacity "Capcal"        
-    !                            the surface conduction flux pcapcal        
-    !                                                                       
-    !                                                                       
-    !   Method: implicit time integration                                   
-    !   -------                                                             
-    !   Consecutive ground temperatures are related by:                     
-    !           T(k+1) = C(k) + D(k)*T(k)  (1)                              
-    !   the coefficients C and D are computed at the t-dt time-step.        
-    !   Routine structure:                                                  
-    !   1)new temperatures are computed  using (1)                          
-    !   2)C and D coefficients are computed from the new temperature        
-    !     profile for the t+dt time-step                                    
-    !   3)the coefficients A and B are computed where the diffusive         
-    !     fluxes at the t+dt time-step is given by                          
-    !            Fdiff = A + B Ts(t+dt)                                     
-    !     or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt                    
-    !            with F0 = A + B (Ts(t))                                    
-    !                 Capcal = B*dt                                         
-    !                                                                       
-    !   Interface:                                                          
-    !   ----------                                                          
-    !                                                                       
-    !   Arguments:                                                          
-    !   ----------                                                          
-    !   ngrid               number of grid-points                           
-    !   ptimestep              physical timestep (s)                        
-    !   pto(ngrid,nsoil)     temperature at time-step t (K)                 
-    !   ptn(ngrid,nsoil)     temperature at time step t+dt (K)              
-    !   pcapcal(ngrid)      specific heat (W*m-2*s*K-1)                     
-    !   pfluxgrd(ngrid)      surface diffusive flux from ground (Wm-2)      
-    !                                                                       
+    !
+    !   Auteur:  Frederic Hourdin     30/01/92
+    !   -------
+    !
+    !   objet:  computation of : the soil temperature evolution
+    !   ------                   the surfacic heat capacity "Capcal"
+    !                            the surface conduction flux pcapcal
+    !
+    !
+    !   Method: implicit time integration
+    !   -------
+    !   Consecutive ground temperatures are related by:
+    !           T(k+1) = C(k) + D(k)*T(k)  (1)
+    !   the coefficients C and D are computed at the t-dt time-step.
+    !   Routine structure:
+    !   1)new temperatures are computed  using (1)
+    !   2)C and D coefficients are computed from the new temperature
+    !     profile for the t+dt time-step
+    !   3)the coefficients A and B are computed where the diffusive
+    !     fluxes at the t+dt time-step is given by
+    !            Fdiff = A + B Ts(t+dt)
+    !     or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
+    !            with F0 = A + B (Ts(t))
+    !                 Capcal = B*dt
+    !
+    !   Interface:
+    !   ----------
+    !
+    !   Arguments:
+    !   ----------
+    !   ngrid               number of grid-points
+    !   ptimestep              physical timestep (s)
+    !   pto(ngrid,nsoil)     temperature at time-step t (K)
+    !   ptn(ngrid,nsoil)     temperature at time step t+dt (K)
+    !   pcapcal(ngrid)      specific heat (W*m-2*s*K-1)
+    !   pfluxgrd(ngrid)      surface diffusive flux from ground (Wm-2)
+    !
     !=======================================================================
-    !   declarations:                                                       
-    !   -------------                                                       
-    
-    
-    !-----------------------------------------------------------------------
-    !  arguments                                                            
-    !  ---------                                                            
-    
-    INTEGER ngrid,nsoil 
-    REAL ptimestep 
-    REAL ptsrf(ngrid),ptsoil(ngrid,nsoil),ptherm_i(ngrid) 
-    REAL pcapcal(ngrid),pfluxgrd(ngrid) 
-    LOGICAL firstcall 
-    
-    
-    !-----------------------------------------------------------------------
-    !  local arrays                                                         
-    !  ------------                                                         
-    
-    INTEGER ig,jk 
-    REAL zdz2(nsoil),z1(ngrid) 
-        
-    IF (firstcall) THEN 
+    !   declarations:
+    !   -------------
+
+
+    !-----------------------------------------------------------------------
+    !  arguments
+    !  ---------
+
+    INTEGER ngrid,nsoil
+    REAL ptimestep
+    REAL ptsrf(ngrid),ptsoil(ngrid,nsoil),ptherm_i(ngrid)
+    REAL pcapcal(ngrid),pfluxgrd(ngrid)
+    LOGICAL firstcall
+
+
+    !-----------------------------------------------------------------------
+    !  local arrays
+    !  ------------
+
+    INTEGER ig,jk
+    REAL zdz2(nsoil),z1(ngrid)
+
+    IF (firstcall) THEN
        CALL init_soil(ngrid, nsoil)
     ELSE
        !-----------------------------------------------------------------------
-       !   Computation of the soil temperatures using the Cgrd and Dgrd        
-       !  coefficient computed at the previous time-step:                      
-       !  -----------------------------------------------                      
-       
-       !    surface temperature                                                
-       DO ig=1,ngrid 
+       !   Computation of the soil temperatures using the Cgrd and Dgrd
+       !  coefficient computed at the previous time-step:
+       !  -----------------------------------------------
+
+       !    surface temperature
+       DO ig=1,ngrid
           ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/                   &
-               &      (lambda*(1.-zd(ig,1))+1.)                                   
+               &      (lambda*(1.-zd(ig,1))+1.)
        ENDDO
-       
-       !   other temperatures                                                  
-       DO jk=1,nsoil-1 
-          DO ig=1,ngrid 
+
+       !   other temperatures
+       DO jk=1,nsoil-1
+          DO ig=1,ngrid
              ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
           ENDDO
        ENDDO
-       
+
     ENDIF
 
     !-----------------------------------------------------------------------
-    !   Computation of the Cgrd and Dgrd coefficient for the next step:     
-    !   ---------------------------------------------------------------     
-    
-    DO jk=1,nsoil 
-       zdz2(jk)=dz2(jk)/ptimestep 
-    ENDDO
-    
-    DO ig=1,ngrid 
-       z1(ig)=zdz2(nsoil)+dz1(nsoil-1) 
-       zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig) 
-       zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig) 
-    ENDDO
-    
-    DO jk=nsoil-1,2,-1 
-       DO ig=1,ngrid 
-          z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk))) 
+    !   Computation of the Cgrd and Dgrd coefficient for the next step:
+    !   ---------------------------------------------------------------
+
+    DO jk=1,nsoil
+       zdz2(jk)=dz2(jk)/ptimestep
+    ENDDO
+
+    DO ig=1,ngrid
+       z1(ig)=zdz2(nsoil)+dz1(nsoil-1)
+       zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig)
+       zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig)
+    ENDDO
+
+    DO jk=nsoil-1,2,-1
+       DO ig=1,ngrid
+          z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
           zc(ig,jk-1)=                                                &
-               &      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig)           
-          zd(ig,jk-1)=dz1(jk-1)*z1(ig) 
+               &      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig)
+          zd(ig,jk-1)=dz1(jk-1)*z1(ig)
        ENDDO
     ENDDO
-    
-    !-----------------------------------------------------------------------
-    !   computation of the surface diffusive flux from ground and           
-    !   calorific capacity of the ground:                                   
-    !   ---------------------------------                                   
-    
-    DO ig=1,ngrid 
+
+    !-----------------------------------------------------------------------
+    !   computation of the surface diffusive flux from ground and
+    !   calorific capacity of the ground:
+    !   ---------------------------------
+
+    DO ig=1,ngrid
        pfluxgrd(ig)=ptherm_i(ig)*dz1(1)*                              &
             &   (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
@@ -198,8 +198,8 @@
        pfluxgrd(ig)=pfluxgrd(ig)                                      &
             &   +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig))   &
-            &   /ptimestep                                                     
-    ENDDO
-    
+            &   /ptimestep
+    ENDDO
+
   END SUBROUTINE soil
-  
+
 END MODULE surface
Index: /dynamico_lmdz/simple_physics/phyparam/physics/turbulence.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/turbulence.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/turbulence.F90	(revision 4229)
@@ -1,4 +1,4 @@
 MODULE turbulence
-  
+
 #include "use_logging.h"
 
@@ -8,10 +8,10 @@
 
   REAL, PARAMETER :: karman=0.4
-  REAL :: lmixmin=100., emin_turb=1e-8  
+  REAL :: lmixmin=100., emin_turb=1e-8
 
   PUBLIC :: vdif, lmixmin, emin_turb
-  
+
 CONTAINS
-  
+
   PURE SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
     !=======================================================================
@@ -47,25 +47,25 @@
     !   Declarations:
     !   -------------
-    
+
     !   Arguments:
     !   ----------
-    
+
     INTEGER, INTENT(IN) ::  ngrid
     REAL, INTENT(IN)  ::  pg, pz0(ngrid), pz(ngrid), &
          pu(ngrid),pv(ngrid), pts(ngrid),ph(ngrid)
     REAL, INTENT(OUT) :: pcdv(ngrid),pcdh(ngrid)
-    
+
     !   Local:
     !   ------
-    
+
     REAL, PARAMETER :: b=5., c=5., d=5., umin2=1e-12, &
          c2b=2.*b, c3bc=3.*b*c, c3b=3.*b
     INTEGER ig
-    REAL zu2,z1,zri,zcd0,zz    
-    
+    REAL zu2,z1,zri,zcd0,zz
+
     !-----------------------------------------------------------------------
     !   couche de surface:
     !   ------------------
-    
+
     !     DO ig=1,ngrid
     !        zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
@@ -74,6 +74,6 @@
     !     ENDDO
     !     RETURN
-    
-    
+
+
 !!!! WARNING, verifier la formule originale de Louis!
     DO ig=1,ngrid
@@ -93,9 +93,9 @@
        ENDIF
     ENDDO
-    
-    !-----------------------------------------------------------------------
-    
+
+    !-----------------------------------------------------------------------
+
   END SUBROUTINE vdif_cd
-  
+
   PURE SUBROUTINE vdif_k(ngrid,nlay,   &
        ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pkv,pkh)
@@ -106,8 +106,8 @@
          pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
     REAL, INTENT(OUT) :: pkv(ngrid,nlay+1),pkh(ngrid,nlay+1)
-    
+
     INTEGER ig,il
     REAL zdu,zdv,zri,zdvodz2,zdz,z1,lmix
-    
+
     DO ig=1,ngrid
        pkv(ig,1)=0.
@@ -116,5 +116,5 @@
        pkh(ig,nlay+1)=0.
     ENDDO
-    
+
     DO il=2,nlay
        DO ig=1,ngrid
@@ -138,12 +138,12 @@
        ENDDO
     ENDDO
-    
+
   END SUBROUTINE vdif_k
 
-  
+
   SUBROUTINE multipl(n,x1,x2,y)
-!=======================================================================
-!   multiplication de deux vecteurs
-!=======================================================================
+    !=======================================================================
+    !   multiplication de deux vecteurs
+    !=======================================================================
     INTEGER n,i
     REAL x1(n),x2(n),y(n)
@@ -162,25 +162,25 @@
        lwrite)
     USE phys_const
-    
-!=======================================================================
-!
-!   Diffusion verticale
-!   Shema implicite
-!   On commence par rajouter au variables x la tendance physique
-!   et on resoult en fait:
-!      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
-!
-!   arguments:
-!   ----------
-!
-!   entree:
-!   -------
-!
-!
-!=======================================================================
-
-!
-!   arguments:
-!   ----------
+
+    !=======================================================================
+    !
+    !   Diffusion verticale
+    !   Shema implicite
+    !   On commence par rajouter au variables x la tendance physique
+    !   et on resoult en fait:
+    !      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
+    !
+    !   arguments:
+    !   ----------
+    !
+    !   entree:
+    !   -------
+    !
+    !
+    !=======================================================================
+
+    !
+    !   arguments:
+    !   ----------
 
     INTEGER ngrid,nlay
@@ -198,5 +198,5 @@
     !   local:
     !   ------
-    
+
     INTEGER ilev,ig,ilay,nlev
     INTEGER unit,ierr,it1,it2
@@ -213,16 +213,16 @@
     REAL zc(ngrid,nlay),zd(ngrid,nlay)
     REAL zcst1
-    
+
     !
     !-----------------------------------------------------------------------
     !   initialisations:
     !   ----------------
-    
+
     nlev=nlay+1
-    
+
     !   computation of rho*dz and dt*rho/dz=dt*rho**2 g/dp:
     !   with rho=p/RT=p/ (R Theta) (p/ps)**kappa
     !   ---------------------------------
-    
+
     DO ilay=1,nlay
        DO ig=1,ngrid
@@ -230,5 +230,5 @@
        ENDDO
     ENDDO
-    
+
     zcst1=4.*g*ptimestep/(r*r)
     DO ilev=2,nlev-1
@@ -258,9 +258,9 @@
        LOG_DBG('vdif')
     ENDIF
-    
+
     !-----------------------------------------------------------------------
     !   2. ajout des tendances physiques:
     !   ------------------------------
-    
+
     DO ilev=1,nlay
        DO ig=1,ngrid
@@ -270,5 +270,5 @@
        ENDDO
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !   3. calcul de  cd :
@@ -277,5 +277,5 @@
     CALL vdif_cd( ngrid,pz0,g,pzlay,pu,pv,ptsrf,ph,zcdv,zcdh)
     CALL vdif_k(ngrid,nlay,ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zkv,zkh)
-    
+
     DO ig=1,ngrid
        zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
@@ -283,6 +283,6 @@
        zcdh(ig)=zcdh(ig)*sqrt(zu2)
     ENDDO
-    
-    IF(lwrite) THEN       
+
+    IF(lwrite) THEN
        WRITELOG(*,*) 'Diagnostique diffusion verticale'
        WRITELOG(*,*) 'LMIXMIN',lmixmin
@@ -295,9 +295,9 @@
        LOG_DBG('vdif')
     ENDIF
-    
+
     !-----------------------------------------------------------------------
     !   integration verticale pour u:
     !   -----------------------------
-    
+
     CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
     CALL multipl(ngrid,zcdv,zb0,zb)
@@ -307,5 +307,5 @@
        zd(ig,nlay)=zb(ig,nlay)*z1(ig)
     ENDDO
-    
+
     DO ilay=nlay-1,1,-1
        DO ig=1,ngrid
@@ -317,5 +317,5 @@
        ENDDO
     ENDDO
-    
+
     DO ig=1,ngrid
        zu(ig,1)=zc(ig,1)
@@ -326,5 +326,5 @@
        ENDDO
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !   integration verticale pour v:
@@ -336,5 +336,5 @@
        zd(ig,nlay)=zb(ig,nlay)*z1(ig)
     ENDDO
-    
+
     DO ilay=nlay-1,1,-1
        DO ig=1,ngrid
@@ -346,5 +346,5 @@
        ENDDO
     ENDDO
-    
+
     DO ig=1,ngrid
        zv(ig,1)=zc(ig,1)
@@ -355,5 +355,5 @@
        ENDDO
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !   integration verticale pour h:
@@ -367,5 +367,5 @@
        zd(ig,nlay)=zb(ig,nlay)*z1(ig)
     ENDDO
-    
+
     DO ilay=nlay-1,1,-1
        DO ig=1,ngrid
@@ -377,18 +377,18 @@
        ENDDO
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !   rajout eventuel de planck dans le shema implicite:
     !   --------------------------------------------------
-    
+
     z4st=4.*5.67e-8*ptimestep
     DO ig=1,ngrid
        zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !   calcul le l'evolution de la temperature du sol:
     !   -----------------------------------------------
-    
+
     DO ig=1,ngrid
        z1(ig) = pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1) &
@@ -399,9 +399,9 @@
        pdtsrf(ig) = (ztsrf2(ig)-ptsrf(ig))/ptimestep
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !   integration verticale finale:
     !   -----------------------------
-    
+
     DO ilay=2,nlay
        DO ig=1,ngrid
@@ -409,9 +409,9 @@
        ENDDO
     ENDDO
-    
+
     !-----------------------------------------------------------------------
     !   calcul final des tendances de la diffusion verticale:
     !   -----------------------------------------------------
-    
+
     DO ilev = 1, nlay
        DO ig=1,ngrid
@@ -424,7 +424,7 @@
        ENDDO
     ENDDO
-    
+
     IF(lwrite) THEN
-       WRITELOG(*,*) 
+       WRITELOG(*,*)
        WRITELOG(*,*) 'Diagnostique de la diffusion verticale'
        WRITELOG(*,*) 'h avant et apres diffusion verticale'
@@ -437,4 +437,4 @@
 
   END SUBROUTINE vdif
-  
+
 END MODULE turbulence
Index: /dynamico_lmdz/simple_physics/phyparam/physics/use_logging.h
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/use_logging.h	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/use_logging.h	(revision 4229)
@@ -1,3 +1,3 @@
-   USE logging
+  USE logging
 #define WRITELOG(junk, fmt) logging_lineno = MIN(logging_bufsize, logging_lineno+1) ; WRITE(logging_buf(logging_lineno), fmt)
 #define LOG_WARN(tag) CALL flush_log(log_level_warn, tag)
Index: /dynamico_lmdz/simple_physics/phyparam/physics/writefield_mod.F90
===================================================================
--- /dynamico_lmdz/simple_physics/phyparam/physics/writefield_mod.F90	(revision 4228)
+++ /dynamico_lmdz/simple_physics/phyparam/physics/writefield_mod.F90	(revision 4229)
@@ -24,5 +24,5 @@
   INTERFACE writefield
      PROCEDURE writefield1, writefield2
-  END INTERFACE writefield
+  END INTERFACE
 
   PUBLIC :: writefield
@@ -53,4 +53,4 @@
 #endif
   END SUBROUTINE writefield1
-  
+
 END MODULE writefield_mod
