Changeset 1862 for LMDZ5/trunk/libf/phylmd/physiq.F90

Timestamp:

Sep 10, 2013, 2:18:14 PM (11 years ago)

Author:

lguez

Message:

Converted physiq.F to free source form (on mandate of the United Poihl).

File:

• LMDZ5/trunk/libf/phylmd/physiq.F90 (moved) (moved from LMDZ5/trunk/libf/phylmd/physiq.F) (149 diffs)

LMDZ5/trunk/libf/phylmd/physiq.F90

@@ -1 +1 @@
 ! $Id$
-c#define IO_DEBUG
-
-      SUBROUTINE physiq (nlon,nlev,
-     .            debut,lafin,jD_cur, jH_cur,pdtphys,
-     .            paprs,pplay,pphi,pphis,presnivs,clesphy0,
-     .            u,v,t,qx,
-     .            flxmass_w,
-     .            d_u, d_v, d_t, d_qx, d_ps
-     .            , dudyn
-     .            , PVteta)
-
-      USE ioipsl, only: histbeg, histvert, histdef, histend, histsync,
-     $     histwrite, ju2ymds, ymds2ju, ioget_year_len
+!#define IO_DEBUG
+
+      SUBROUTINE physiq (nlon,nlev, &
+                  debut,lafin,jD_cur, jH_cur,pdtphys, &
+                  paprs,pplay,pphi,pphis,presnivs,clesphy0, &
+                  u,v,t,qx, &
+                  flxmass_w, &
+                  d_u, d_v, d_t, d_qx, d_ps &
+                  , dudyn &
+                  , PVteta)
+
+      USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
+           histwrite, ju2ymds, ymds2ju, ioget_year_len
       USE comgeomphy
       USE phys_cal_mod
@@ -37 +37 @@
       use regr_pr_av_m, only: regr_pr_av
       use netcdf95, only: nf95_close
-cIM for NMC files
-c     use netcdf, only: nf90_fill_real
+!IM for NMC files
+!     use netcdf, only: nf90_fill_real
       use netcdf
       use mod_phys_lmdz_mpi_data, only: is_mpi_root
@@ -69 +69 @@
 !!   =====================
 #define histNMC
-c#define histISCCP
+!#define histISCCP
 !!======================================================================
 !!    modif   ( P. Le Van ,  12/10/98 )
@@ -114 +114 @@
 #include "iniprint.h"
 #include "thermcell.h"
-c======================================================================
+!======================================================================
       LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE
       PARAMETER (ok_cvl=.TRUE.)
@@ -121 +121 @@
       integer iflag_radia     ! active ou non le rayonnement (MPL)
       save iflag_radia
-c$OMP THREADPRIVATE(iflag_radia)
-c======================================================================
+!$OMP THREADPRIVATE(iflag_radia)
+!======================================================================
       LOGICAL check ! Verifier la conservation du modele en eau
       PARAMETER (check=.FALSE.)
       LOGICAL ok_stratus ! Ajouter artificiellement les stratus
       PARAMETER (ok_stratus=.FALSE.)
-c======================================================================
+!======================================================================
       REAL amn, amx
       INTEGER igout
-c======================================================================
-c Clef controlant l'activation du cycle diurne:
-ccc      LOGICAL cycle_diurne
-ccc      PARAMETER (cycle_diurne=.FALSE.)
-c======================================================================
-c Modele thermique du sol, a activer pour le cycle diurne:
-ccc      LOGICAL soil_model
-ccc      PARAMETER (soil_model=.FALSE.)
-c======================================================================
-c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
-c le calcul du rayonnement est celle apres la precipitation des nuages.
-c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
-c la condensation et la precipitation. Cette cle augmente les impacts
-c radiatifs des nuages.
-ccc      LOGICAL new_oliq
-ccc      PARAMETER (new_oliq=.FALSE.)
-c======================================================================
-c Clefs controlant deux parametrisations de l'orographie:
-cc      LOGICAL ok_orodr
-ccc      PARAMETER (ok_orodr=.FALSE.)
-ccc      LOGICAL ok_orolf
-ccc      PARAMETER (ok_orolf=.FALSE.)
-c======================================================================
+!======================================================================
+! Clef controlant l'activation du cycle diurne:
+!cc      LOGICAL cycle_diurne
+!cc      PARAMETER (cycle_diurne=.FALSE.)
+!======================================================================
+! Modele thermique du sol, a activer pour le cycle diurne:
+!cc      LOGICAL soil_model
+!cc      PARAMETER (soil_model=.FALSE.)
+!======================================================================
+! Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
+! le calcul du rayonnement est celle apres la precipitation des nuages.
+! Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
+! la condensation et la precipitation. Cette cle augmente les impacts
+! radiatifs des nuages.
+!cc      LOGICAL new_oliq
+!cc      PARAMETER (new_oliq=.FALSE.)
+!======================================================================
+! Clefs controlant deux parametrisations de l'orographie:
+!c      LOGICAL ok_orodr
+!cc      PARAMETER (ok_orodr=.FALSE.)
+!cc      LOGICAL ok_orolf
+!cc      PARAMETER (ok_orolf=.FALSE.)
+!======================================================================
       LOGICAL ok_journe ! sortir le fichier journalier
       save ok_journe
-c$OMP THREADPRIVATE(ok_journe)
-c
+!$OMP THREADPRIVATE(ok_journe)
+!
       LOGICAL ok_mensuel ! sortir le fichier mensuel
       save ok_mensuel
-c$OMP THREADPRIVATE(ok_mensuel)
-c
+!$OMP THREADPRIVATE(ok_mensuel)
+!
       LOGICAL ok_instan ! sortir le fichier instantane
       save ok_instan
-c$OMP THREADPRIVATE(ok_instan)
-c
+!$OMP THREADPRIVATE(ok_instan)
+!
       LOGICAL ok_LES ! sortir le fichier LES 
       save ok_LES                            
-c$OMP THREADPRIVATE(ok_LES)                  
-c
+!$OMP THREADPRIVATE(ok_LES)                  
+!
       LOGICAL callstats ! sortir le fichier stats 
       save callstats                            
-c$OMP THREADPRIVATE(callstats)                  
-c
+!$OMP THREADPRIVATE(callstats)                  
+!
       LOGICAL ok_region ! sortir le fichier regional
       PARAMETER (ok_region=.FALSE.)
-c======================================================================
+!======================================================================
       real seuil_inversion
       save seuil_inversion
-c$OMP THREADPRIVATE(seuil_inversion)
+!$OMP THREADPRIVATE(seuil_inversion)
       integer iflag_ratqs
       save iflag_ratqs
-c$OMP THREADPRIVATE(iflag_ratqs)
+!$OMP THREADPRIVATE(iflag_ratqs)
       real facteur
 
@@ -194 +194 @@
       real zqsatth(klon,klev)
 
-c======================================================================
-c
+!======================================================================
+!
       INTEGER ivap          ! indice de traceurs pour vapeur d'eau
       PARAMETER (ivap=1)
@@ -201 +201 @@
       PARAMETER (iliq=2)
 
-c
-c
-c Variables argument:
-c
+!
+!
+! Variables argument:
+!
       INTEGER nlon
       INTEGER nlev
@@ -222 +222 @@
       REAL v(klon,klev)
       REAL t(klon,klev),thetal(klon,klev)
-c thetal: ligne suivante a decommenter si vous avez les fichiers     MPL 20130625
-c fth_fonctions.F90 et parkind1.F90
-c sinon thetal=theta
-c     REAL fth_thetae,fth_thetav,fth_thetal
+! thetal: ligne suivante a decommenter si vous avez les fichiers     MPL 20130625
+! fth_fonctions.F90 et parkind1.F90
+! sinon thetal=theta
+!     REAL fth_thetae,fth_thetav,fth_thetal
       REAL qx(klon,klev,nqtot)
       REAL flxmass_w(klon,klev)
@@ -245 +245 @@
 !IM definition dynamique o_trac dans phys_output_open
 !      type(ctrl_out) :: o_trac(nqtot)
-c
-cIM Amip2 PV a theta constante 
-c
+!
+!IM Amip2 PV a theta constante 
+!
       INTEGER nbteta
       PARAMETER(nbteta=3)
@@ -253 +253 @@
       DATA ctetaSTD/'350','380','405'/
       SAVE ctetaSTD
-c$OMP THREADPRIVATE(ctetaSTD)
+!$OMP THREADPRIVATE(ctetaSTD)
       REAL rtetaSTD(nbteta)
       DATA rtetaSTD/350., 380., 405./
       SAVE rtetaSTD
-c$OMP THREADPRIVATE(rtetaSTD)     
-c
+!$OMP THREADPRIVATE(rtetaSTD)     
+!
       REAL PVteta(klon,nbteta)
       REAL zx_tmp_3dte(iim,jjmp1,nbteta)
-c
-cMI Amip2 PV a theta constante
-
-cym      INTEGER klevp1, klevm1
-cym      PARAMETER(klevp1=klev+1,klevm1=klev-1)
-cym#include "raddim.h"
-c
-c
-cIM Amip2
-c variables a une pression donnee
-c
+!
+!MI Amip2 PV a theta constante
+
+!ym      INTEGER klevp1, klevm1
+!ym      PARAMETER(klevp1=klev+1,klevm1=klev-1)
+!ym#include "raddim.h"
+!
+!
+!IM Amip2
+! variables a une pression donnee
+!
 #include "declare_STDlev.h"
-c
+!
       CHARACTER*4 bb2
       CHARACTER*2 bb3
-c
+!
 #include "radopt.h"
-c
-c
+!
+!
 
       REAL convliq(klon,klev)  ! eau liquide nuageuse convective
@@ -290 +290 @@
 
       INTEGER linv, kp1
-c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
-c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
+! flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
+! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
       REAL flwp_c(klon), fiwp_c(klon)
       REAL flwc_c(klon,klev), fiwc_c(klon,klev)
@@ -298 +298 @@
 
 
-cIM ISCCP simulator v3.4
-c dans clesphys.h top_height, overlap
-cv3.4
+!IM ISCCP simulator v3.4
+! dans clesphys.h top_height, overlap
+!v3.4
       INTEGER debug, debugcol
-cym      INTEGER npoints
-cym      PARAMETER(npoints=klon) 
-c
+!ym      INTEGER npoints
+!ym      PARAMETER(npoints=klon) 
+!
       INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
       INTEGER nregISCtot
       PARAMETER(nregISCtot=1) 
-c
-c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
-c y compris pour 1 point
-c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
-c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
+!
+! imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
+! y compris pour 1 point
+! imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
+! jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
       INTEGER imin_debut, nbpti
       INTEGER jmin_debut, nbptj 
-cIM parametres ISCCP BEG
+!IM parametres ISCCP BEG
       INTEGER nbapp_isccp
 !     INTEGER nbapp_isccp,isccppas
@@ -324 +324 @@
       DATA ifreq_isccp/3/
       SAVE ifreq_isccp
-c$OMP THREADPRIVATE(ifreq_isccp)
+!$OMP THREADPRIVATE(ifreq_isccp)
       CHARACTER*5 typinout(napisccp)
       DATA typinout/'i3od'/
       SAVE typinout
-c$OMP THREADPRIVATE(typinout)
-cIM verif boxptop BEG
+!$OMP THREADPRIVATE(typinout)
+!IM verif boxptop BEG
       CHARACTER*1 verticaxe(napisccp)
       DATA verticaxe/'1'/ 
       SAVE verticaxe
-c$OMP THREADPRIVATE(verticaxe)
-cIM verif boxptop END
+!$OMP THREADPRIVATE(verticaxe)
+!IM verif boxptop END
       INTEGER nvlev(napisccp)
-c     INTEGER nvlev
+!     INTEGER nvlev
       REAL t1, aa
       REAL seed_re(klon,napisccp)
-cym !!!! A voir plus tard 
-cym      INTEGER iphy(iim,jjmp1)
-cIM parametres ISCCP END
-c
-c ncol = nb. de sous-colonnes pour chaque maille du GCM 
-c ncolmx = No. max. de sous-colonnes pour chaque maille du GCM 
-c      INTEGER ncol(napisccp), ncolmx, seed(klon,napisccp)
+!ym !!!! A voir plus tard 
+!ym      INTEGER iphy(iim,jjmp1)
+!IM parametres ISCCP END
+!
+! ncol = nb. de sous-colonnes pour chaque maille du GCM 
+! ncolmx = No. max. de sous-colonnes pour chaque maille du GCM 
+!      INTEGER ncol(napisccp), ncolmx, seed(klon,napisccp)
       INTEGER,SAVE :: ncol(napisccp)
-c$OMP THREADPRIVATE(ncol)
+!$OMP THREADPRIVATE(ncol)
       INTEGER ncolmx, seed(klon,napisccp)
       REAL nbsunlit(nregISCtot,klon,napisccp)  !nbsunlit : moyenne de sunlit
-c     PARAMETER(ncolmx=1500)
+!     PARAMETER(ncolmx=1500)
       PARAMETER(ncolmx=300)
-c
-cIM verif boxptop BEG
+!
+!IM verif boxptop BEG
       REAL vertlev(ncolmx,napisccp)
-cIM verif boxptop END
-c
+!IM verif boxptop END
+!
       REAL,SAVE :: tautab_omp(0:255),tautab(0:255)
       INTEGER,SAVE :: invtau_omp(-20:45000),invtau(-20:45000)
-c$OMP THREADPRIVATE(tautab,invtau)
+!$OMP THREADPRIVATE(tautab,invtau)
       REAL emsfc_lw
       PARAMETER(emsfc_lw=0.99)
-c     REAL    ran0                      ! type for random number fuction
-c
+!     REAL    ran0                      ! type for random number fuction
+!
       REAL cldtot(klon,klev)
-c variables de haut en bas pour le simulateur ISCCP
+! variables de haut en bas pour le simulateur ISCCP
       REAL dtau_s(klon,klev) !tau nuages startiformes
       REAL dtau_c(klon,klev) !tau nuages convectifs
       REAL dem_s(klon,klev)  !emissivite nuages startiformes 
       REAL dem_c(klon,klev)  !emissivite nuages convectifs
-c
-c variables de haut en bas pour le simulateur ISCCP
+!
+! variables de haut en bas pour le simulateur ISCCP
       REAL pfull(klon,klev)
       REAL phalf(klon,klev+1)
@@ -382 +382 @@
       REAL dem_sH2B(klon,klev)
       REAL dem_cH2B(klon,klev)
-c
+!
       INTEGER kmax, lmax, lmax3
       PARAMETER(kmax=8, lmax=8, lmax3=3)
@@ -388 +388 @@
       PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
       INTEGER iimx7, jjmx7, jjmp1x7
-      PARAMETER(iimx7=iim*kmaxm1, jjmx7=jjm*lmaxm1, 
-     .jjmp1x7=jjmp1*lmaxm1)
-c
-c output from ISCCP simulator
+      PARAMETER(iimx7=iim*kmaxm1, jjmx7=jjm*lmaxm1,  &
+      jjmp1x7=jjmp1*lmaxm1)
+!
+! output from ISCCP simulator
       REAL fq_isccp(klon,kmaxm1,lmaxm1,napisccp)
       REAL fq_is_true(klon,kmaxm1,lmaxm1,napisccp)
@@ -401 +401 @@
       REAL zx_tmp_fi3d_bx(klon,ncolmx)
       REAL zx_tmp_3d_bx(iim,jjmp1,ncolmx)
-c
+!
       REAL cld_fi3d(klon,lmax3)
       REAL cld_3d(iim,jjmp1,lmax3)
-c
+!
       INTEGER iw, iwmax
       REAL wmin, pas_w
-c     PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
-cIM 051005     PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
+!     PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
+!IM 051005     PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
       PARAMETER(wmin=-100.,pas_w=10.,iwmax=20)
       REAL o500(klon)
-c
-
-c sorties ISCCP
+!
+
+! sorties ISCCP
 
       integer nid_isccp
       save nid_isccp        
-c$OMP THREADPRIVATE(nid_isccp)
+!$OMP THREADPRIVATE(nid_isccp)
 
       REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
@@ -424 +424 @@
       DATA zx_pc/180., 310., 440., 560., 680., 800., 1000./
       SAVE zx_pc
-c$OMP THREADPRIVATE(zx_tau,zx_pc)
-c cldtopres pression au sommet des nuages
+!$OMP THREADPRIVATE(zx_tau,zx_pc)
+! cldtopres pression au sommet des nuages
       REAL cldtopres(lmaxm1), cldtopres3(lmax3)
       DATA cldtopres/180., 310., 440., 560., 680., 800., 1000./
       DATA cldtopres3/440., 680., 1000./
       SAVE cldtopres,cldtopres3
-c$OMP THREADPRIVATE(cldtopres,cldtopres3)
-cIM 051005 BEG
+!$OMP THREADPRIVATE(cldtopres,cldtopres3)
+!IM 051005 BEG
       INTEGER komega, nhoriRD 
 
-c taulev: numero du niveau de tau dans les sorties ISCCP
+! taulev: numero du niveau de tau dans les sorties ISCCP
       CHARACTER *4 taulev(kmaxm1)
-c     DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
+!     DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
       DATA taulev/'tau0','tau1','tau2','tau3','tau4','tau5','tau6'/
       CHARACTER *3 pclev(lmaxm1)
       DATA pclev/'pc1','pc2','pc3','pc4','pc5','pc6','pc7'/
       SAVE taulev,pclev
-c$OMP THREADPRIVATE(taulev,pclev)
-c
-c cnameisccp
+!$OMP THREADPRIVATE(taulev,pclev)
+!
+! cnameisccp
       CHARACTER *29 cnameisccp(lmaxm1,kmaxm1)
-cIM bad 151205     DATA cnameisccp/'pc< 50hPa, tau< 0.3', 
-      DATA cnameisccp/'pc= 50-180hPa, tau< 0.3',
-     .                'pc= 180-310hPa, tau< 0.3',
-     .                'pc= 310-440hPa, tau< 0.3',
-     .                'pc= 440-560hPa, tau< 0.3',
-     .                'pc= 560-680hPa, tau< 0.3',
-     .                'pc= 680-800hPa, tau< 0.3',
-     .                'pc= 800-1000hPa, tau< 0.3',
-     .                'pc= 50-180hPa, tau= 0.3-1.3',
-     .                'pc= 180-310hPa, tau= 0.3-1.3',
-     .                'pc= 310-440hPa, tau= 0.3-1.3',
-     .                'pc= 440-560hPa, tau= 0.3-1.3',
-     .                'pc= 560-680hPa, tau= 0.3-1.3',
-     .                'pc= 680-800hPa, tau= 0.3-1.3',
-     .                'pc= 800-1000hPa, tau= 0.3-1.3',
-     .                'pc= 50-180hPa, tau= 1.3-3.6',
-     .                'pc= 180-310hPa, tau= 1.3-3.6',
-     .                'pc= 310-440hPa, tau= 1.3-3.6',
-     .                'pc= 440-560hPa, tau= 1.3-3.6',
-     .                'pc= 560-680hPa, tau= 1.3-3.6',
-     .                'pc= 680-800hPa, tau= 1.3-3.6',
-     .                'pc= 800-1000hPa, tau= 1.3-3.6',
-     .                'pc= 50-180hPa, tau= 3.6-9.4',
-     .                'pc= 180-310hPa, tau= 3.6-9.4',
-     .                'pc= 310-440hPa, tau= 3.6-9.4',
-     .                'pc= 440-560hPa, tau= 3.6-9.4',
-     .                'pc= 560-680hPa, tau= 3.6-9.4',
-     .                'pc= 680-800hPa, tau= 3.6-9.4',
-     .                'pc= 800-1000hPa, tau= 3.6-9.4',
-     .                'pc= 50-180hPa, tau= 9.4-23',
-     .                'pc= 180-310hPa, tau= 9.4-23',
-     .                'pc= 310-440hPa, tau= 9.4-23',
-     .                'pc= 440-560hPa, tau= 9.4-23',
-     .                'pc= 560-680hPa, tau= 9.4-23',
-     .                'pc= 680-800hPa, tau= 9.4-23',
-     .                'pc= 800-1000hPa, tau= 9.4-23',
-     .                'pc= 50-180hPa, tau= 23-60',
-     .                'pc= 180-310hPa, tau= 23-60',
-     .                'pc= 310-440hPa, tau= 23-60',
-     .                'pc= 440-560hPa, tau= 23-60',
-     .                'pc= 560-680hPa, tau= 23-60',
-     .                'pc= 680-800hPa, tau= 23-60',
-     .                'pc= 800-1000hPa, tau= 23-60',
-     .                'pc= 50-180hPa, tau> 60.',
-     .                'pc= 180-310hPa, tau> 60.',
-     .                'pc= 310-440hPa, tau> 60.',
-     .                'pc= 440-560hPa, tau> 60.',
-     .                'pc= 560-680hPa, tau> 60.',
-     .                'pc= 680-800hPa, tau> 60.',
-     .                'pc= 800-1000hPa, tau> 60.'/
+!IM bad 151205     DATA cnameisccp/'pc< 50hPa, tau< 0.3', 
+      DATA cnameisccp/'pc= 50-180hPa, tau< 0.3', &
+                      'pc= 180-310hPa, tau< 0.3', &
+                      'pc= 310-440hPa, tau< 0.3', &
+                      'pc= 440-560hPa, tau< 0.3', &
+                      'pc= 560-680hPa, tau< 0.3', &
+                      'pc= 680-800hPa, tau< 0.3', &
+                      'pc= 800-1000hPa, tau< 0.3', &
+                      'pc= 50-180hPa, tau= 0.3-1.3', &
+                      'pc= 180-310hPa, tau= 0.3-1.3', &
+                      'pc= 310-440hPa, tau= 0.3-1.3', &
+                      'pc= 440-560hPa, tau= 0.3-1.3', &
+                      'pc= 560-680hPa, tau= 0.3-1.3', &
+                      'pc= 680-800hPa, tau= 0.3-1.3', &
+                      'pc= 800-1000hPa, tau= 0.3-1.3', &
+                      'pc= 50-180hPa, tau= 1.3-3.6', &
+                      'pc= 180-310hPa, tau= 1.3-3.6', &
+                      'pc= 310-440hPa, tau= 1.3-3.6', &
+                      'pc= 440-560hPa, tau= 1.3-3.6', &
+                      'pc= 560-680hPa, tau= 1.3-3.6', &
+                      'pc= 680-800hPa, tau= 1.3-3.6', &
+                      'pc= 800-1000hPa, tau= 1.3-3.6', &
+                      'pc= 50-180hPa, tau= 3.6-9.4', &
+                      'pc= 180-310hPa, tau= 3.6-9.4', &
+                      'pc= 310-440hPa, tau= 3.6-9.4', &
+                      'pc= 440-560hPa, tau= 3.6-9.4', &
+                      'pc= 560-680hPa, tau= 3.6-9.4', &
+                      'pc= 680-800hPa, tau= 3.6-9.4', &
+                      'pc= 800-1000hPa, tau= 3.6-9.4', &
+                      'pc= 50-180hPa, tau= 9.4-23', &
+                      'pc= 180-310hPa, tau= 9.4-23', &
+                      'pc= 310-440hPa, tau= 9.4-23', &
+                      'pc= 440-560hPa, tau= 9.4-23', &
+                      'pc= 560-680hPa, tau= 9.4-23', &
+                      'pc= 680-800hPa, tau= 9.4-23', &
+                      'pc= 800-1000hPa, tau= 9.4-23', &
+                      'pc= 50-180hPa, tau= 23-60', &
+                      'pc= 180-310hPa, tau= 23-60', &
+                      'pc= 310-440hPa, tau= 23-60', &
+                      'pc= 440-560hPa, tau= 23-60', &
+                      'pc= 560-680hPa, tau= 23-60', &
+                      'pc= 680-800hPa, tau= 23-60', &
+                      'pc= 800-1000hPa, tau= 23-60', &
+                      'pc= 50-180hPa, tau> 60.', &
+                      'pc= 180-310hPa, tau> 60.', &
+                      'pc= 310-440hPa, tau> 60.', &
+                      'pc= 440-560hPa, tau> 60.', &
+                      'pc= 560-680hPa, tau> 60.', &
+                      'pc= 680-800hPa, tau> 60.', &
+                      'pc= 800-1000hPa, tau> 60.'/
        SAVE cnameisccp
-c$OMP THREADPRIVATE(cnameisccp)
-c
-c     REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
-c     INTEGER nhorix7
-cIM: region='3d' <==> sorties en global
+!$OMP THREADPRIVATE(cnameisccp)
+!
+!     REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
+!     INTEGER nhorix7
+!IM: region='3d' <==> sorties en global
       CHARACTER*3 region
       PARAMETER(region='3d')
-c
-cIM ISCCP simulator v3.4
-c
+!
+!IM ISCCP simulator v3.4
+!
       logical ok_hf
-c
+!
       integer nid_hf, nid_hf3d
       save ok_hf, nid_hf, nid_hf3d
-c$OMP THREADPRIVATE(ok_hf, nid_hf, nid_hf3d)
-c  QUESTION : noms de variables ?
+!$OMP THREADPRIVATE(ok_hf, nid_hf, nid_hf3d)
+!  QUESTION : noms de variables ?
 
       INTEGER        longcles
       PARAMETER    ( longcles = 20 )
       REAL clesphy0( longcles      )
-c
-c Variables propres a la physique
+!
+! Variables propres a la physique
       INTEGER itap
       SAVE itap                   ! compteur pour la physique
-c$OMP THREADPRIVATE(itap)
-c
+!$OMP THREADPRIVATE(itap)
+!
       REAL,save ::  solarlong0
-c$OMP THREADPRIVATE(solarlong0)
-
-c
-c  Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
-c
-cIM 141004     REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
+!$OMP THREADPRIVATE(solarlong0)
+
+!
+!  Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
+!
+!IM 141004     REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
       REAL zulow(klon),zvlow(klon)
-c
+!
       INTEGER igwd,idx(klon),itest(klon)
-c
-c      REAL,allocatable,save :: run_off_lic_0(:)
-cc$OMP THREADPRIVATE(run_off_lic_0)
-cym      SAVE run_off_lic_0
-cKE43
-c Variables liees a la convection de K. Emanuel (sb):
-c
+!
+!      REAL,allocatable,save :: run_off_lic_0(:)
+!!$OMP THREADPRIVATE(run_off_lic_0)
+!ym      SAVE run_off_lic_0
+!KE43
+! Variables liees a la convection de K. Emanuel (sb):
+!
       REAL bas, top             ! cloud base and top levels
       SAVE bas
       SAVE top
-c$OMP THREADPRIVATE(bas, top)
+!$OMP THREADPRIVATE(bas, top)
 
       REAL wdn(klon), tdn(klon), qdn(klon)
-c
-c=================================================================================================
-cCR04.12.07: on ajoute les nouvelles variables du nouveau schema de convection avec poches froides
-c Variables li\'ees \`a la poche froide (jyg)
+!
+!=================================================================================================
+!CR04.12.07: on ajoute les nouvelles variables du nouveau schema de convection avec poches froides
+! Variables li\'ees \`a la poche froide (jyg)
 
       REAL mip(klon,klev)  ! mass flux shed by the adiab ascent at each level
-c
+!
       REAL wape_prescr, fip_prescr
       INTEGER it_wape_prescr
       SAVE wape_prescr, fip_prescr, it_wape_prescr
-c$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
-c
-c variables supplementaires de concvl
+!$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
+!
+! variables supplementaires de concvl
       REAL Tconv(klon,klev)
       REAL ment(klon,klev,klev),sij(klon,klev,klev)
@@ -570 +570 @@
       real, save :: wake_s_min_lsp=0.1
 
-c$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
-c$OMP THREADPRIVATE(ale_max,alp_max)
-c$OMP THREADPRIVATE(wake_s_min_lsp)
+!$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
+!$OMP THREADPRIVATE(ale_max,alp_max)
+!$OMP THREADPRIVATE(wake_s_min_lsp)
 
 
       real ok_wk_lsp(klon)
 
-cRC
-c Variables li\'ees \`a la poche froide (jyg et rr)
-c Version diagnostique pour l'instant : pas de r\'etroaction sur la convection
+!RC
+! Variables li\'ees \`a la poche froide (jyg et rr)
+! Version diagnostique pour l'instant : pas de r\'etroaction sur la convection
 
       REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
@@ -595 +595 @@
       REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
       REAL wake_spread(klon,klev)     ! spreading term in wake_delt
-c
-cpourquoi y'a pas de save??
-c
+!
+!pourquoi y'a pas de save??
+!
       INTEGER wake_k(klon)            ! Wake sommet
-c
+!
       REAL t_undi(klon,klev)               ! temperature moyenne dans la zone non perturbee
       REAL q_undi(klon,klev)               ! humidite moyenne dans la zone non perturbee
-c
-cjyg
-ccc      REAL wake_pe(klon)              ! Wake potential energy - WAPE 
+!
+!jyg
+!cc      REAL wake_pe(klon)              ! Wake potential energy - WAPE 
 
       REAL wake_gfl(klon)             ! Gust Front Length
       REAL wake_dens(klon)
-c
-c
+!
+!
       REAL dt_dwn(klon,klev)
       REAL dq_dwn(klon,klev)
@@ -621 +621 @@
       REAL dq_a(klon,klev)
       REAL, SAVE :: alp_offset
-c$OMP THREADPRIVATE(alp_offset)
-
-c
-cRR:fin declarations poches froides
-c=======================================================================================================
+!$OMP THREADPRIVATE(alp_offset)
+
+!
+!RR:fin declarations poches froides
+!=======================================================================================================
         
       REAL ztv(klon,klev),ztva(klon,klev)
@@ -632 +632 @@
       REAL zthl(klon,klev)
 
-ccc nrlmd le 10/04/2012
-
-c--------Stochastic Boundary Layer Triggering: ALE_BL--------
-c---Propri\'et\'es du thermiques au LCL 
+!cc nrlmd le 10/04/2012
+
+!--------Stochastic Boundary Layer Triggering: ALE_BL--------
+!---Propri\'et\'es du thermiques au LCL 
       real zlcl_th(klon)                                     ! Altitude du LCL calcul\'e continument (pcon dans thermcell_main.F90)
       real fraca0(klon)                                      ! Fraction des thermiques au LCL
@@ -644 +644 @@
       real env_tke_max0(klon)                                ! TKE dans l'environnement au LCL 
 
-c---D\'eclenchement stochastique
+!---D\'eclenchement stochastique
       integer :: tau_trig(klon)
 
-c--------Statistical Boundary Layer Closure: ALP_BL--------
-c---Profils de TKE dans et hors du thermique
+!--------Statistical Boundary Layer Closure: ALP_BL--------
+!---Profils de TKE dans et hors du thermique
       real pbl_tke_input(klon,klev+1,nbsrf)
       real therm_tke_max(klon,klev)                          ! Profil de TKE dans les thermiques
@@ -654 +654 @@
 
 
-ccc fin nrlmd le 10/04/2012
-
-c Variables locales pour la couche limite (al1):
-c
-cAl1      REAL pblh(klon)           ! Hauteur de couche limite
-cAl1      SAVE pblh
-c34EK
-c
-c Variables locales:
-c
-cAA
-cAA  Pour phytrac 
+!cc fin nrlmd le 10/04/2012
+
+! Variables locales pour la couche limite (al1):
+!
+!Al1      REAL pblh(klon)           ! Hauteur de couche limite
+!Al1      SAVE pblh
+!34EK
+!
+! Variables locales:
+!
+!AA
+!AA  Pour phytrac 
       REAL u1(klon)             ! vents dans la premiere couche U
       REAL v1(klon)             ! vents dans la premiere couche V
 
-c@$$      LOGICAL offline           ! Controle du stockage ds "physique"
-c@$$      PARAMETER (offline=.false.)
-c@$$      INTEGER physid
+!@$$      LOGICAL offline           ! Controle du stockage ds "physique"
+!@$$      PARAMETER (offline=.false.)
+!@$$      INTEGER physid
       REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
       REAL frac_nucl(klon,klev) ! idem (nucleation)
@@ -680 +680 @@
       REAL          :: calday
 
-cIM cf FH pour Tiedtke 080604
+!IM cf FH pour Tiedtke 080604
       REAL rain_tiedtke(klon),snow_tiedtke(klon)
-c
-cIM 050204 END
+!
+!IM 050204 END
       REAL devap(klon) ! evaporation et sa derivee
       REAL dsens(klon) ! chaleur sensible et sa derivee
 
-c
-c Conditions aux limites
-c
+!
+! Conditions aux limites
+!
 !
       REAL :: day_since_equinox
@@ -698 +698 @@
       LOGICAL, parameter :: new_orbit = .true.
 
-c
+!
       INTEGER lmt_pas
       SAVE lmt_pas                ! frequence de mise a jour
-c$OMP THREADPRIVATE(lmt_pas) 
+!$OMP THREADPRIVATE(lmt_pas) 
       real zmasse(klon, llm),exner(klon, llm) 
-C     (column-density of mass of air in a cell, in kg m-2)
+!     (column-density of mass of air in a cell, in kg m-2)
       real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
 
-cIM sorties
+!IM sorties
       REAL un_jour
       PARAMETER(un_jour=86400.)
-c======================================================================
-c
-c Declaration des procedures appelees
-c
+!======================================================================
+!
+! Declaration des procedures appelees
+!
       EXTERNAL angle     ! calculer angle zenithal du soleil
       EXTERNAL alboc     ! calculer l'albedo sur ocean
       EXTERNAL ajsec     ! ajustement sec
       EXTERNAL conlmd    ! convection (schema LMD)
-cKE43
+!KE43
       EXTERNAL conema3  ! convect4.3
       EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)
-cAA 
+!AA 
       EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
-c                          ! stockage des coefficients necessaires au
-c                          ! lessivage OFF-LINE et ON-LINE
+!                          ! stockage des coefficients necessaires au
+!                          ! lessivage OFF-LINE et ON-LINE
       EXTERNAL hgardfou  ! verifier les temperatures
       EXTERNAL nuage     ! calculer les proprietes radiatives
-CC      EXTERNAL o3cm      ! initialiser l'ozone
+!C      EXTERNAL o3cm      ! initialiser l'ozone
       EXTERNAL orbite    ! calculer l'orbite terrestre
       EXTERNAL phyetat0  ! lire l'etat initial de la physique
@@ -736 +736 @@
       EXTERNAL ecrirega  ! ecrire le fichier binaire regional
       EXTERNAL ecriregs  ! ecrire le fichier binaire regional
-cIM
+!IM
       EXTERNAL haut2bas  !variables de haut en bas
       EXTERNAL ini_undefSTD  !initialise a 0 une variable a 1 niveau de pression
       EXTERNAL undefSTD      !somme les valeurs definies d'1 var a 1 niveau de pression
-c     EXTERNAL moy_undefSTD  !moyenne d'1 var a 1 niveau de pression
-c     EXTERNAL moyglo_aire   !moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
-c                            !par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
-c
-c Variables locales
-c
+!     EXTERNAL moy_undefSTD  !moyenne d'1 var a 1 niveau de pression
+!     EXTERNAL moyglo_aire   !moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
+!                            !par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
+!
+! Variables locales
+!
       REAL rhcl(klon,klev)    ! humiditi relative ciel clair
       REAL dialiq(klon,klev)  ! eau liquide nuageuse
       REAL diafra(klon,klev)  ! fraction nuageuse
       REAL cldliq(klon,klev)  ! eau liquide nuageuse
-c
-CXXX PB 
+!
+!XXX PB 
       REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
-c
+!
       REAL zxfluxt(klon, klev)
       REAL zxfluxq(klon, klev)
@@ -759 +759 @@
       REAL zxfluxv(klon, klev)
 
-c Le rayonnement n'est pas calcule tous les pas, il faut donc
-c                      sauvegarder les sorties du rayonnement
-cym      SAVE  heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
-cym      SAVE  sollwdownclr, toplwdown, toplwdownclr
-cym      SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
-c
+! Le rayonnement n'est pas calcule tous les pas, il faut donc
+!                      sauvegarder les sorties du rayonnement
+!ym      SAVE  heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
+!ym      SAVE  sollwdownclr, toplwdown, toplwdownclr
+!ym      SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
+!
       INTEGER itaprad
       SAVE itaprad
-c$OMP THREADPRIVATE(itaprad)
-c
+!$OMP THREADPRIVATE(itaprad)
+!
       REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
       REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
 
-c
+!
       REAL zxsnow(klon)
       REAL zxsnow_dummy(klon)
-c
+!
       REAL dist, rmu0(klon), fract(klon)
       REAL zdtime, zlongi
-c
+!
       CHARACTER*2 str2
       CHARACTER*2 iqn
-c
+!
       REAL qcheck
       REAL z_avant(klon), z_apres(klon), z_factor(klon)
       LOGICAL zx_ajustq
-c
+!
       REAL za, zb
       REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
@@ -793 +793 @@
       PARAMETER (t_coup=234.0)
 
-cym A voir plus tard !!
-cym      REAL zx_relief(iim,jjmp1)
-cym      REAL zx_aire(iim,jjmp1)
-c
-c Grandeurs de sorties
+!ym A voir plus tard !!
+!ym      REAL zx_relief(iim,jjmp1)
+!ym      REAL zx_aire(iim,jjmp1)
+!
+! Grandeurs de sorties
       REAL s_capCL(klon)
       REAL s_oliqCL(klon), s_cteiCL(klon)
       REAL s_trmb1(klon), s_trmb2(klon)
       REAL s_trmb3(klon)
-cKE43
-c Variables locales pour la convection de K. Emanuel (sb):
+!KE43
+! Variables locales pour la convection de K. Emanuel (sb):
 
       REAL tvp(klon,klev)       ! virtual temp of lifted parcel
@@ -811 +811 @@
       INTEGER iflagctrl(klon)          ! flag fonctionnement de convect
 
-c -- convect43:
+! -- convect43:
       INTEGER ntra              ! nb traceurs pour convect4.3
       REAL pori_con(klon)    ! pressure at the origin level of lifted parcel
@@ -817 +817 @@
       REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
       REAL dplcldt(klon), dplcldr(klon)
-c?     .     condm_con(klon,klev),conda_con(klon,klev),
-c?     .     mr_con(klon,klev),ep_con(klon,klev)
-c?     .    ,sadiab(klon,klev),wadiab(klon,klev)
-c --
-c34EK
-c
-c Variables du changement
-c
-c con: convection
-c lsc: condensation a grande echelle (Large-Scale-Condensation)
-c ajs: ajustement sec
-c eva: evaporation de l'eau liquide nuageuse
-c vdf: couche limite (Vertical DiFfusion)
+!?     .     condm_con(klon,klev),conda_con(klon,klev),
+!?     .     mr_con(klon,klev),ep_con(klon,klev)
+!?     .    ,sadiab(klon,klev),wadiab(klon,klev)
+! --
+!34EK
+!
+! Variables du changement
+!
+! con: convection
+! lsc: condensation a grande echelle (Large-Scale-Condensation)
+! ajs: ajustement sec
+! eva: evaporation de l'eau liquide nuageuse
+! vdf: couche limite (Vertical DiFfusion)
 
 ! tendance nulles
       REAL du0(klon,klev),dv0(klon,klev),dq0(klon,klev),dql0(klon,klev)
 
-c
-*********************************************************
-*     declarations
+!
+!********************************************************
+!     declarations
       
-*********************************************************
-cIM 081204 END
-c
+!********************************************************
+!IM 081204 END
+!
       REAL pen_u(klon,klev), pen_d(klon,klev)
       REAL pde_u(klon,klev), pde_d(klon,klev)
       INTEGER kcbot(klon), kctop(klon), kdtop(klon)
-c
+!
       REAL ratqsc(klon,klev)
       real ratqsbas,ratqshaut,tau_ratqs
       save ratqsbas,ratqshaut,tau_ratqs
-c$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
+!$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
       real zpt_conv(klon,klev)
 
-c Parametres lies au nouveau schema de nuages (SB, PDF)
+! Parametres lies au nouveau schema de nuages (SB, PDF)
       real fact_cldcon
       real facttemps
       logical ok_newmicro
       save ok_newmicro
-c$OMP THREADPRIVATE(ok_newmicro)
+!$OMP THREADPRIVATE(ok_newmicro)
       save fact_cldcon,facttemps
-c$OMP THREADPRIVATE(fact_cldcon,facttemps)
+!$OMP THREADPRIVATE(fact_cldcon,facttemps)
 
       integer iflag_cldcon
       save iflag_cldcon
-c$OMP THREADPRIVATE(iflag_cldcon)
+!$OMP THREADPRIVATE(iflag_cldcon)
       logical ptconv(klon,klev)
-cIM cf. AM 081204 BEG
+!IM cf. AM 081204 BEG
       logical ptconvth(klon,klev)
-cIM cf. AM 081204 END
-c
-c Variables liees a l'ecriture de la bande histoire physique
-c
-c======================================================================
-c
-cIM cf. AM 081204 BEG
-c   declarations pour sortir sur une sous-region
+!IM cf. AM 081204 END
+!
+! Variables liees a l'ecriture de la bande histoire physique
+!
+!======================================================================
+!
+!IM cf. AM 081204 BEG
+!   declarations pour sortir sur une sous-region
       integer imin_ins,imax_ins,jmin_ins,jmax_ins
       save imin_ins,imax_ins,jmin_ins,jmax_ins
-c$OMP THREADPRIVATE(imin_ins,imax_ins,jmin_ins,jmax_ins)
-c      real lonmin_ins,lonmax_ins,latmin_ins
-c     s  ,latmax_ins
-c     data lonmin_ins,lonmax_ins,latmin_ins
-c    s  ,latmax_ins/
-c valeurs initiales     s   -5.,20.,41.,55./   
-c    s   100.,130.,-20.,20./
-c    s   -180.,180.,-90.,90./
-c======================================================================
-cIM cf. AM 081204 END
-
-c
+!$OMP THREADPRIVATE(imin_ins,imax_ins,jmin_ins,jmax_ins)
+!      real lonmin_ins,lonmax_ins,latmin_ins
+!     s  ,latmax_ins
+!     data lonmin_ins,lonmax_ins,latmin_ins
+!    s  ,latmax_ins/
+! valeurs initiales     s   -5.,20.,41.,55./   
+!    s   100.,130.,-20.,20./
+!    s   -180.,180.,-90.,90./
+!======================================================================
+!IM cf. AM 081204 END
+
+!
       integer itau_w   ! pas de temps ecriture = itap + itau_phy
-c
-c
-c Variables locales pour effectuer les appels en serie
-c
-cIM RH a 2m (la surface)
+!
+!
+! Variables locales pour effectuer les appels en serie
+!
+!IM RH a 2m (la surface)
       REAL Lheat
 
@@ -899 +899 @@
       PARAMETER    ( length = 100 )
       REAL tabcntr0( length       )
-c
+!
       INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev)
-cIM
+!IM
       INTEGER ndex2d1(iwmax)
-c
-cIM AMIP2 BEG
+!
+!IM AMIP2 BEG
       REAL moyglo, mountor
-cIM 141004 BEG
+!IM 141004 BEG
       REAL zustrdr(klon), zvstrdr(klon)
       REAL zustrli(klon), zvstrli(klon)
@@ -912 +912 @@
       REAL zustrhi(klon), zvstrhi(klon)
       REAL aam, torsfc
-cIM 141004 END
-cIM 190504 BEG
+!IM 141004 END
+!IM 190504 BEG
       INTEGER ij, imp1jmp1
       PARAMETER(imp1jmp1=(iim+1)*jjmp1)
-cym A voir plus tard
+!ym A voir plus tard
       REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
       REAL padyn(iim+1,jjmp1,klev+1)
       REAL dudyn(iim+1,jjmp1,klev)
       REAL rlatdyn(iim+1,jjmp1)
-cIM 190504 END
+!IM 190504 END
       LOGICAL ok_msk
       REAL msk(klon)
-cIM 
+!IM 
       REAL airetot, pi
-cym A voir plus tard
-cym      REAL zm_wo(jjmp1, klev)
-cIM AMIP2 END
-c
+!ym A voir plus tard
+!ym      REAL zm_wo(jjmp1, klev)
+!IM AMIP2 END
+!
       REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
       REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D 
@@ -936 +936 @@
       REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
       REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
-c
+!
       INTEGER nid_day, nid_mth, nid_ins, nid_mthnmc, nid_daynmc
       INTEGER nid_hfnmc, nid_day_seri, nid_ctesGCM
       SAVE nid_day, nid_mth, nid_ins, nid_mthnmc, nid_daynmc
       SAVE nid_hfnmc, nid_day_seri, nid_ctesGCM
-c$OMP THREADPRIVATE(nid_day, nid_mth, nid_ins)
-c$OMP THREADPRIVATE(nid_mthnmc, nid_daynmc, nid_hfnmc)
-c$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
-c
-cIM 280405 BEG
+!$OMP THREADPRIVATE(nid_day, nid_mth, nid_ins)
+!$OMP THREADPRIVATE(nid_mthnmc, nid_daynmc, nid_hfnmc)
+!$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
+!
+!IM 280405 BEG
       INTEGER nid_bilKPins, nid_bilKPave
       SAVE nid_bilKPins, nid_bilKPave
-c$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
-c
+!$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
+!
       REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
       REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
       REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
       REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
-c
+!
       INTEGER nhori, nvert, nvert1, nvert3
       REAL zsto, zsto1, zsto2
@@ -961 +961 @@
       REAL zout_isccp(napisccp)
       SAVE zcals, zcalh, zoutj, zout_isccp
-c$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
+!$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
 
       real zjulian
       save zjulian
-c$OMP THREADPRIVATE(zjulian)
+!$OMP THREADPRIVATE(zjulian)
 
       character*20 modname
@@ -973 +973 @@
       integer idayref
 
-C essai writephys
+! essai writephys
       integer fid_day, fid_mth, fid_ins
       parameter (fid_ins = 1, fid_day = 2, fid_mth = 3) 
       integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
-      parameter (prof2d_on = 1, prof3d_on = 2,
-     .           prof2d_av = 3, prof3d_av = 4)
+      parameter (prof2d_on = 1, prof3d_on = 2, &
+                 prof2d_av = 3, prof3d_av = 4)
       character*30 nom_fichier
       character*40 varname
       character*40 vartitle
       character*20 varunits
-C     Variables liees au bilan d'energie et d'enthalpi
+!     Variables liees au bilan d'energie et d'enthalpi
       REAL ztsol(klon)
-      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
-     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
-      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
-     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
-c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot,
-c$OMP+              h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
+      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot &
+              , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
+      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot &
+              , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
+!$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot)
+!$OMP THREADPRIVATE(h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
       REAL      d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
       REAL      d_h_vcol_phy
       REAL      fs_bound, fq_bound
       SAVE      d_h_vcol_phy
-c$OMP THREADPRIVATE(d_h_vcol_phy)
+!$OMP THREADPRIVATE(d_h_vcol_phy)
       REAL      zero_v(klon)
       CHARACTER*15 ztit
@@ -1001 +1001 @@
       SAVE      ip_ebil
       DATA      ip_ebil/0/
-c$OMP THREADPRIVATE(ip_ebil)
+!$OMP THREADPRIVATE(ip_ebil)
       INTEGER   if_ebil ! level for energy conserv. dignostics
       SAVE      if_ebil
-c$OMP THREADPRIVATE(if_ebil)
-c+jld ec_conser
+!$OMP THREADPRIVATE(if_ebil)
+!+jld ec_conser
       REAL ZRCPD
-c-jld ec_conser
+!-jld ec_conser
       REAL q2m(klon,nbsrf)  ! humidite a 2m
 
-cIM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
+!IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
       CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
       CHARACTER*40 tinst, tave, typeval
@@ -1030 +1030 @@
       REAL bl95_b0, bl95_b1   ! Parameter in Boucher and Lohmann (1995)
       SAVE ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1
-c$OMP THREADPRIVATE(ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1)
+!$OMP THREADPRIVATE(ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1)
       LOGICAL, SAVE :: aerosol_couple ! true  : calcul des aerosols dans INCA
                                       ! false : lecture des aerosol dans un fichier
-c$OMP THREADPRIVATE(aerosol_couple)    
+!$OMP THREADPRIVATE(aerosol_couple)    
       INTEGER, SAVE :: flag_aerosol 
-c$OMP THREADPRIVATE(flag_aerosol) 
+!$OMP THREADPRIVATE(flag_aerosol) 
       LOGICAL, SAVE :: new_aod
-c$OMP THREADPRIVATE(new_aod) 
-c
-c--STRAT AEROSOL
+!$OMP THREADPRIVATE(new_aod) 
+!
+!--STRAT AEROSOL
       LOGICAL, SAVE :: flag_aerosol_strat
-c$OMP THREADPRIVATE(flag_aerosol_strat)
-cc-fin STRAT AEROSOL
-c
-c Declaration des constantes et des fonctions thermodynamiques
-c
+!$OMP THREADPRIVATE(flag_aerosol_strat)
+!c-fin STRAT AEROSOL
+!
+! Declaration des constantes et des fonctions thermodynamiques
+!
       LOGICAL,SAVE :: first=.true.
-c$OMP THREADPRIVATE(first)
+!$OMP THREADPRIVATE(first)
 
       integer iunit
 
       integer, save::  read_climoz ! read ozone climatology
-C     (let it keep the default OpenMP shared attribute)
-C     Allowed values are 0, 1 and 2
-C     0: do not read an ozone climatology
-C     1: read a single ozone climatology that will be used day and night
-C     2: read two ozone climatologies, the average day and night
-C     climatology and the daylight climatology
+!     (let it keep the default OpenMP shared attribute)
+!     Allowed values are 0, 1 and 2
+!     0: do not read an ozone climatology
+!     1: read a single ozone climatology that will be used day and night
+!     2: read two ozone climatologies, the average day and night
+!     climatology and the daylight climatology
 
       integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
-C     (let it keep the default OpenMP shared attribute)
+!     (let it keep the default OpenMP shared attribute)
 
       real, pointer, save:: press_climoz(:)
-C     (let it keep the default OpenMP shared attribute)
+!     (let it keep the default OpenMP shared attribute)
 !     edges of pressure intervals for ozone climatologies, in Pa, in strictly
 !     ascending order
@@ -1069 +1069 @@
       integer, save:: co3i = 0
 !     time index in NetCDF file of current ozone fields
-c$OMP THREADPRIVATE(co3i) 
+!$OMP THREADPRIVATE(co3i) 
 
       integer ro3i
@@ -1079 +1079 @@
 #include "YOETHF.h"
 #include "FCTTRE.h"
-cIM 100106 BEG : pouvoir sortir les ctes de la physique
+!IM 100106 BEG : pouvoir sortir les ctes de la physique
 #include "conema3.h"
 #include "fisrtilp.h"
 #include "nuage.h"
 #include "compbl.h"
-cIM 100106 END : pouvoir sortir les ctes de la physique
-c
+!IM 100106 END : pouvoir sortir les ctes de la physique
+!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-c Declarations pour Simulateur COSP
-c============================================================
+! Declarations pour Simulateur COSP
+!============================================================
       real :: mr_ozone(klon,klev)
 
-cIM sorties fichier 1D paramLMDZ_phy.nc
+!IM sorties fichier 1D paramLMDZ_phy.nc
       REAL :: zx_tmp_0d(1,1)
       INTEGER, PARAMETER :: np=1
@@ -1099 +1099 @@
       REAL grain(1), gtsol(1), gt2m(1), gprw(1)
 
-cIM stations CFMIP
+!IM stations CFMIP
       INTEGER, SAVE :: nCFMIP
-c$OMP THREADPRIVATE(nCFMIP)
+!$OMP THREADPRIVATE(nCFMIP)
       INTEGER, PARAMETER :: npCFMIP=120
       INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
       REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
-c$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP) 
+!$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP) 
       INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
       REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
-c$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
+!$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
       INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
-c$OMP THREADPRIVATE(iGCM, jGCM)
+!$OMP THREADPRIVATE(iGCM, jGCM)
       logical, dimension(nfiles)            :: phys_out_filestations
       logical, parameter :: lNMC=.FALSE.
 
-cIM betaCRF
+!IM betaCRF
       REAL, SAVE :: pfree, beta_pbl, beta_free
-c$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
+!$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
       REAL, SAVE :: lon1_beta,  lon2_beta, lat1_beta, lat2_beta
-c$OMP THREADPRIVATE(lon1_beta,  lon2_beta, lat1_beta, lat2_beta)
+!$OMP THREADPRIVATE(lon1_beta,  lon2_beta, lat1_beta, lat2_beta)
       LOGICAL, SAVE :: mskocean_beta
-c$OMP THREADPRIVATE(mskocean_beta)
+!$OMP THREADPRIVATE(mskocean_beta)
       REAL, dimension(klon, klev) :: beta         ! facteur sur cldtaurad et cldemirad pour evaluer les retros liees aux CRF
       REAL, dimension(klon, klev) :: cldtaurad    ! epaisseur optique pour radlwsw pour tester "CRF off"
@@ -1131 +1131 @@
       integer iostat
 
-c======================================================================
+!======================================================================
 ! Gestion calendrier : mise a jour du module phys_cal_mod
 !
       CALL phys_cal_update(jD_cur,jH_cur)
 
-c======================================================================
+!======================================================================
 ! Ecriture eventuelle d'un profil verticale en entree de la physique.
 ! Utilise notamment en 1D mais peut etre active egalement en 3D
 ! en imposant la valeur de igout.
-c======================================================================d
+!======================================================================d
       if (prt_level.ge.1) then
           igout=klon/2+1/klon
          write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
-         write(lunout,*)
-     s 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
-         write(lunout,*)
-     s  nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys 
+         write(lunout,*) &
+       'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
+         write(lunout,*) &
+        nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys 
 
          write(lunout,*) 'paprs, play, phi, u, v, t'
          do k=1,klev
-            write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
-     s   u(igout,k),v(igout,k),t(igout,k)
+            write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
+         u(igout,k),v(igout,k),t(igout,k)
          enddo
          write(lunout,*) 'ovap (g/kg),  oliq (g/kg)'
@@ -1160 +1160 @@
       endif
 
-c======================================================================
+!======================================================================
 
       if (first) then 
       
-cCR:nvelles variables convection/poches froides
+!CR:nvelles variables convection/poches froides
       
       print*, '================================================='
       print*, 'Allocation des variables locales et sauvegardees'
       call phys_local_var_init
-c
+!
       pasphys=pdtphys
-c     appel a la lecture du run.def physique
-      call conf_phys(ok_journe, ok_mensuel,
-     .     ok_instan, ok_hf,
-     .     ok_LES,
-     .     callstats,
-     .     solarlong0,seuil_inversion,
-     .     fact_cldcon, facttemps,ok_newmicro,iflag_radia,
-     .     iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,
-     .     ok_ade, ok_aie, ok_cdnc, aerosol_couple, 
-     .     flag_aerosol, flag_aerosol_strat, new_aod,
-     .     bl95_b0, bl95_b1,
-c     nv flags pour la convection et les poches froides
-     .     read_climoz,
-     &     alp_offset)
+!     appel a la lecture du run.def physique
+      call conf_phys(ok_journe, ok_mensuel, &
+           ok_instan, ok_hf, &
+           ok_LES, &
+           callstats, &
+           solarlong0,seuil_inversion, &
+           fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
+           iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
+           ok_ade, ok_aie, ok_cdnc, aerosol_couple,  &
+           flag_aerosol, flag_aerosol_strat, new_aod, &
+           bl95_b0, bl95_b1, &
+!     nv flags pour la convection et les poches froides
+           read_climoz, &
+           alp_offset)
       call phys_state_var_init(read_climoz)
       call phys_output_var_init
       print*, '================================================='
-c
+!
           dnwd0=0.0
           ftd=0.0
           fqd=0.0
           cin=0.
-cym Attention pbase pas initialise dans concvl !!!!
+!ym Attention pbase pas initialise dans concvl !!!!
           pbase=0
-cIM 180608
+!IM 180608
 
         itau_con=0
@@ -1202 +1202 @@
       endif  ! first
 
-cym => necessaire pour iflag_con != 2    
+!ym => necessaire pour iflag_con != 2    
       pmfd(:,:) = 0.
       pen_u(:,:) = 0.
@@ -1216 +1216 @@
 
        modname = 'physiq'
-cIM
+!IM
       IF (ip_ebil_phy.ge.1) THEN
         DO i=1,klon
@@ -1231 +1231 @@
 
 
-c======================================================================
+!======================================================================
 ! Gestion calendrier : mise a jour du module phys_cal_mod
 !
-c     CALL phys_cal_update(jD_cur,jH_cur)
-
-c
-c Si c'est le debut, il faut initialiser plusieurs choses
-c          ********
-c
+!     CALL phys_cal_update(jD_cur,jH_cur)
+
+!
+! Si c'est le debut, il faut initialiser plusieurs choses
+!          ********
+!
        IF (debut) THEN
 !rv
-cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation 
-cde la convection a partir des caracteristiques du thermique
+!CRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation 
+!de la convection a partir des caracteristiques du thermique
          wght_th(:,:)=1.
          lalim_conv(:)=1 
-cRC
+!RC
          ustar(:,:)=0.
          u10m(:,:)=0.
@@ -1273 +1273 @@
          clwcon(:,:) = 0.0
 
-cIM      
+!IM      
          IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
-c
-      print*,'iflag_coupl,iflag_clos,iflag_wake',
-     .   iflag_coupl,iflag_clos,iflag_wake
+!
+      print*,'iflag_coupl,iflag_clos,iflag_wake', &
+         iflag_coupl,iflag_clos,iflag_wake
       print*,'CYCLE_DIURNE', cycle_diurne
-c
+!
       IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
          abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
          CALL abort_gcm (modname,abort_message,1)
       ENDIF
-c
+!
       IF(ok_isccp.AND.iflag_con.LE.2) THEN
-         abort_message = 'ISCCP-like outputs may be available for KE
-     .(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
+         abort_message = 'ISCCP-like outputs may be available for KE' // &
+     '(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
          CALL abort_gcm (modname,abort_message,1)
       ENDIF
-c
-c Initialiser les compteurs:
-c
+!
+! Initialiser les compteurs:
+!
          itap    = 0
          itaprad = 0
@@ -1320 +1320 @@
          PRINT*,'FH WARNING : lignes a supprimer'
          ENDIF
-cIM begin
-          print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
-     $,ratqs(1,1)
-cIM end
+!IM begin
+          print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
+      ,ratqs(1,1)
+!IM end
 
 
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-c
-C on remet le calendrier a zero
-c
+!
+! on remet le calendrier a zero
+!
          IF (raz_date .eq. 1) THEN
            itau_phy = 0
          ENDIF
 
-cIM cf. AM 081204 BEG
+!IM cf. AM 081204 BEG
          PRINT*,'cycle_diurne3 =',cycle_diurne
-cIM cf. AM 081204 END
-c
-         CALL printflag( tabcntr0,radpas,ok_journe,
-     ,                    ok_instan, ok_region )
-c
+!IM cf. AM 081204 END
+!
+         CALL printflag( tabcntr0,radpas,ok_journe, &
+                          ok_instan, ok_region )
+!
          IF (ABS(dtime-pdtphys).GT.0.001) THEN
-            WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
-     .                        pdtphys
+            WRITE(lunout,*) 'Pas physique n est pas correct',dtime, &
+                              pdtphys
             abort_message='Pas physique n est pas correct '
 !           call abort_gcm(modname,abort_message,1)
@@ -1350 +1350 @@
          ENDIF
          IF (nlon .NE. klon) THEN
-            WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, 
-     .                      klon
+            WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,  &
+                            klon
             abort_message='nlon et klon ne sont pas coherents'
             call abort_gcm(modname,abort_message,1)
          ENDIF
          IF (nlev .NE. klev) THEN
-            WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
-     .                       klev
+            WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
+                             klev
             abort_message='nlev et klev ne sont pas coherents'
             call abort_gcm(modname,abort_message,1)
          ENDIF
-c
+!
          IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN 
            WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
@@ -1369 +1369 @@
          ENDIF
          WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
-         WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
-     .                   ok_cvl
-c
-cKE43
-c Initialisation pour la convection de K.E. (sb):
+         WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
+                         ok_cvl
+!
+!KE43
+! Initialisation pour la convection de K.E. (sb):
          IF (iflag_con.GE.3) THEN
 
          WRITE(lunout,*)"*** Convection de Kerry Emanuel 4.3  "
-         WRITE(lunout,*)
-     .      "On va utiliser le melange convectif des traceurs qui"
+         WRITE(lunout,*) &
+            "On va utiliser le melange convectif des traceurs qui"
          WRITE(lunout,*)"est calcule dans convect4.3"
          WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
@@ -1386 +1386 @@
            ema_pcb(i)  = 0.
            ema_pct(i)  = 0.
-c          ema_workcbmf(i) = 0.
+!          ema_workcbmf(i) = 0.
           ENDDO
-cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
+!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
           DO i = 1, klon
            ibas_con(i) = 1
            itop_con(i) = 1
           ENDDO
-cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
-c===============================================================================
-cCR:04.12.07: initialisations poches froides
-c Controle de ALE et ALP pour la fermeture convective (jyg)
+!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
+!===============================================================================
+!CR:04.12.07: initialisations poches froides
+! Controle de ALE et ALP pour la fermeture convective (jyg)
           if (iflag_wake>=1) then
-            CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
-     s                  ,alp_bl_prescr, ale_bl_prescr)
-c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
-c        print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
+            CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
+                        ,alp_bl_prescr, ale_bl_prescr)
+! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
+!        print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
           endif
 
@@ -1409 +1409 @@
         enddo
 
-c================================================================================
-cIM stations CFMIP
+!================================================================================
+!IM stations CFMIP
       nCFMIP=npCFMIP
-      OPEN(98,file='npCFMIP_param.data',status='old',
-     $          form='formatted',iostat=iostat)
+      OPEN(98,file='npCFMIP_param.data',status='old', &
+                form='formatted',iostat=iostat)
             if (iostat == 0) then
       READ(98,*,end=998) nCFMIP
@@ -1426 +1426 @@
       ENDIF
 
-c
+!
       ALLOCATE(tabCFMIP(nCFMIP))
       ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
@@ -1432 +1432 @@
       ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
       ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
-c
-c lecture des nCFMIP stations CFMIP, de leur numero 
-c et des coordonnees geographiques lonCFMIP, latCFMIP
-c
-         CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, 
-     $lonCFMIP, latCFMIP)
-c
-c identification des
-c 1) coordonnees lonGCM, latGCM des points CFMIP dans la grille de LMDZ
-c 2) indices points tabijGCM de la grille physique 1d sur klon points
-c 3) indices iGCM, jGCM de la grille physique 2d
-c
-         CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP,
-     $tabijGCM, lonGCM, latGCM, iGCM, jGCM)
-c
+!
+! lecture des nCFMIP stations CFMIP, de leur numero 
+! et des coordonnees geographiques lonCFMIP, latCFMIP
+!
+         CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP,  &
+      lonCFMIP, latCFMIP)
+!
+! identification des
+! 1) coordonnees lonGCM, latGCM des points CFMIP dans la grille de LMDZ
+! 2) indices points tabijGCM de la grille physique 1d sur klon points
+! 3) indices iGCM, jGCM de la grille physique 2d
+!
+         CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
+      tabijGCM, lonGCM, latGCM, iGCM, jGCM)
+!
             else
                ALLOCATE(tabijGCM(0))
@@ -1462 +1462 @@
            ENDDO
 
-c34EK
+!34EK
          IF (ok_orodr) THEN
 
@@ -1489 +1489 @@
            ENDDO
          ENDIF
-c
-c
+!
+!
          lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
-         WRITE(lunout,*)'La frequence de lecture surface est de ', 
-     .                   lmt_pas
-c
+         WRITE(lunout,*)'La frequence de lecture surface est de ',  &
+                         lmt_pas
+!
       capemaxcels = 't_max(X)'
       t2mincels = 't_min(X)'
@@ -1500 +1500 @@
       tinst = 'inst(X)'
       tave = 'ave(X)'
-cIM cf. AM 081204 BEG
+!IM cf. AM 081204 BEG
       write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
-cIM cf. AM 081204 END
-c
-c=============================================================
-c   Initialisation des sorties
-c=============================================================
+!IM cf. AM 081204 END
+!
+!=============================================================
+!   Initialisation des sorties
+!=============================================================
 
 #ifdef CPP_IOIPSL
 
-c$OMP MASTER
-       call phys_output_open(rlon,rlat,nCFMIP,tabijGCM,
-     &                       iGCM,jGCM,lonGCM,latGCM,
-     &                       jjmp1,nlevSTD,clevSTD,rlevSTD,
-     &                       nbteta, ctetaSTD, dtime,ok_veget,
-     &                       type_ocean,iflag_pbl,ok_mensuel,ok_journe,
-     &                       ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, 
-     &                       read_climoz, phys_out_filestations,
-     &                       new_aod, aerosol_couple,
-     &                       flag_aerosol_strat, pdtphys, paprs, pphis, 
-     &                       pplay, lmax_th, ptconv, ptconvth, ivap, 
-     &                       d_t, qx, d_qx, zmasse, ok_sync)
-c$OMP END MASTER
-c$OMP BARRIER
+!$OMP MASTER
+       call phys_output_open(rlon,rlat,nCFMIP,tabijGCM, &
+                             iGCM,jGCM,lonGCM,latGCM, &
+                             jjmp1,nlevSTD,clevSTD,rlevSTD, &
+                             nbteta, ctetaSTD, dtime,ok_veget, &
+                             type_ocean,iflag_pbl,ok_mensuel,ok_journe, &
+                             ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,  &
+                             read_climoz, phys_out_filestations, &
+                             new_aod, aerosol_couple, &
+                             flag_aerosol_strat, pdtphys, paprs, pphis,  &
+                             pplay, lmax_th, ptconv, ptconvth, ivap,  &
+                             d_t, qx, d_qx, zmasse, ok_sync)
+!$OMP END MASTER
+!$OMP BARRIER
 
 #undef histISCCP
@@ -1538 +1538 @@
          ecrit_tra = ecrit_tra * un_jour
        
-cXXXPB Positionner date0 pour initialisation de ORCHIDEE
+!XXXPB Positionner date0 pour initialisation de ORCHIDEE
       date0 = jD_ref 
       WRITE(*,*) 'physiq date0 : ',date0
-c
-c
-c
-c Prescrire l'ozone dans l'atmosphere
-c
-c
-cc         DO i = 1, klon
-cc         DO k = 1, klev
-cc            CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
-cc         ENDDO
-cc         ENDDO
-c
+!
+!
+!
+! Prescrire l'ozone dans l'atmosphere
+!
+!
+!c         DO i = 1, klon
+!c         DO k = 1, klev
+!c            CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
+!c         ENDDO
+!c         ENDDO
+!
       IF (type_trac == 'inca') THEN
 #ifdef INCA
@@ -1562 +1562 @@
          WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
 
-         CALL chemini( 
-     $                   rg,
-     $                   ra,
-     $                   airephy,
-     $                   rlat,
-     $                   rlon,
-     $                   presnivs,
-     $                   calday,
-     $                   klon,
-     $                   nqtot,
-     $                   pdtphys,
-     $                   annee_ref,
-     $                   day_ref, 
-     $                   itau_phy)
+         CALL chemini(  &
+                         rg, &
+                         ra, &
+                         airephy, &
+                         rlat, &
+                         rlon, &
+                         presnivs, &
+                         calday, &
+                         klon, &
+                         nqtot, &
+                         pdtphys, &
+                         annee_ref, &
+                         day_ref,  &
+                         itau_phy)
 
          CALL VTe(VTinca)
@@ -1581 +1581 @@
 #endif
       END IF
-c
-c
+!
+!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! Nouvelle initialisation pour le rayonnement RRTM
@@ -1589 +1589 @@
       call iniradia(klon,klev,paprs(1,1:klev+1))
 
-C$omp single
+!$omp single
       if (read_climoz >= 1) then
          call open_climoz(ncid_climoz, press_climoz)
       END IF
-C$omp end single
-c
-cIM betaCRF
+!$omp end single
+!
+!IM betaCRF
       pfree=70000. !Pa
       beta_pbl=1.
@@ -1605 +1605 @@
       mskocean_beta=.FALSE.
 
-      OPEN(99,file='beta_crf.data',status='old',
-     $          form='formatted',err=9999)
+      OPEN(99,file='beta_crf.data',status='old', &
+                form='formatted',err=9999)
       READ(99,*,end=9998) pfree
       READ(99,*,end=9998) beta_pbl
@@ -1634 +1634 @@
 !
       itap   = itap + 1
-c
+!
 !
 ! Update fraction of the sub-surfaces (pctsrf) and 
@@ -1640 +1640 @@
 ! on the surface fraction.
 !
-      CALL change_srf_frac(itap, dtime, days_elapsed+1, 
-     *     pctsrf, falb1, falb2, ftsol, ustar, u10m, v10m, pbl_tke)
+      CALL change_srf_frac(itap, dtime, days_elapsed+1,  &
+           pctsrf, falb1, falb2, ftsol, ustar, u10m, v10m, pbl_tke)
 
 
@@ -1660 +1660 @@
       dq0(:,:)=0.
       dql0(:,:)=0.
-c
-c Mettre a zero des variables de sortie (pour securite)
-c
+!
+! Mettre a zero des variables de sortie (pour securite)
+!
       DO i = 1, klon
          d_ps(i) = 0.0
@@ -1695 +1695 @@
 ! RomP <<<
 
-c
-c Ne pas affecter les valeurs entrees de u, v, h, et q
-c
+!
+! Ne pas affecter les valeurs entrees de u, v, h, et q
+!
       DO k = 1, klev
       DO i = 1, klon
@@ -1724 +1724 @@
       ENDDO
       ENDIF
-C
+!
       DO i = 1, klon
         ztsol(i) = 0.
@@ -1733 +1733 @@
         ENDDO
       ENDDO
-cIM
+!IM
       IF (ip_ebil_phy.ge.1) THEN 
         ztit='after dynamic'
-        CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
-     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
-     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-C     Comme les tendances de la physique sont ajoute dans la dynamique,
-C     on devrait avoir que la variation d'entalpie par la dynamique
-C     est egale a la variation de la physique au pas de temps precedent.
-C     Donc la somme de ces 2 variations devrait etre nulle.
-        call diagphy(airephy,ztit,ip_ebil_phy
-     e      , zero_v, zero_v, zero_v, zero_v, zero_v
-     e      , zero_v, zero_v, zero_v, ztsol
-     e      , d_h_vcol+d_h_vcol_phy, d_qt, 0.
-     s      , fs_bound, fq_bound )
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
+            , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
+            , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+!     Comme les tendances de la physique sont ajoute dans la dynamique,
+!     on devrait avoir que la variation d'entalpie par la dynamique
+!     est egale a la variation de la physique au pas de temps precedent.
+!     Donc la somme de ces 2 variations devrait etre nulle.
+        call diagphy(airephy,ztit,ip_ebil_phy &
+            , zero_v, zero_v, zero_v, zero_v, zero_v &
+            , zero_v, zero_v, zero_v, ztsol &
+            , d_h_vcol+d_h_vcol_phy, d_qt, 0. &
+            , fs_bound, fq_bound )
       END IF 
 
-c Diagnostiquer la tendance dynamique
-c
+! Diagnostiquer la tendance dynamique
+!
       IF (ancien_ok) THEN
          DO k = 1, klev
@@ -1766 +1766 @@
           DO k = 1, klev
           DO i = 1, klon
-            d_tr_dyn(i,k,iq-2)=
-     $       (tr_seri(i,k,iq-2)-tr_ancien(i,k,iq-2))/dtime
+            d_tr_dyn(i,k,iq-2)= &
+             (tr_seri(i,k,iq-2)-tr_ancien(i,k,iq-2))/dtime
 !         iiq=niadv(iq)
 !         print*,i,k," d_tr_dyn",d_tr_dyn(i,k,iq-2),"tra:",iq,tname(iiq)
@@ -1797 +1797 @@
          ancien_ok = .TRUE.
       ENDIF
-c
-c Ajouter le geopotentiel du sol:
-c
+!
+! Ajouter le geopotentiel du sol:
+!
       DO k = 1, klev
       DO i = 1, klon
@@ -1805 +1805 @@
       ENDDO
       ENDDO
-c
-c Verifier les temperatures
-c
-cIM BEG
+!
+! Verifier les temperatures
+!
+!IM BEG
       IF (check) THEN
        amn=MIN(ftsol(1,is_ter),1000.)
@@ -1816 +1816 @@
         amx=MAX(ftsol(i,is_ter),amx)
        ENDDO
-c
+!
        PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
       ENDIF !(check) THEN
-cIM END
-c
+!IM END
+!
       CALL hgardfou(t_seri,ftsol,'debutphy')
-c
-cIM BEG
+!
+!IM BEG
       IF (check) THEN
        amn=MIN(ftsol(1,is_ter),1000.)
@@ -1831 +1831 @@
         amx=MAX(ftsol(i,is_ter),amx)
        ENDDO
-c
+!
        PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
       ENDIF !(check) THEN
-cIM END
-c
-c Mettre en action les conditions aux limites (albedo, sst, etc.).
-c Prescrire l'ozone et calculer l'albedo sur l'ocean.
-c
+!IM END
+!
+! Mettre en action les conditions aux limites (albedo, sst, etc.).
+! Prescrire l'ozone et calculer l'albedo sur l'ocean.
+!
       if (read_climoz >= 1) then
-C        Ozone from a file
+!        Ozone from a file
 !        Update required ozone index:
-         ro3i = int((days_elapsed + jh_cur - jh_1jan)
-     $        / ioget_year_len(year_cur) * 360.) + 1
+         ro3i = int((days_elapsed + jh_cur - jh_1jan) &
+              / ioget_year_len(year_cur) * 360.) + 1
          if (ro3i == 361) ro3i = 360
-C        (This should never occur, except perhaps because of roundup
-C        error. See documentation.)
+!        (This should never occur, except perhaps because of roundup
+!        error. See documentation.)
          if (ro3i /= co3i) then
-C           Update ozone field:
+!           Update ozone field:
             if (read_climoz == 1) then
-               call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,
-     $              press_in_edg=press_climoz, paprs=paprs, v3=wo)
+               call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
+                    press_in_edg=press_climoz, paprs=paprs, v3=wo)
             else
-C              read_climoz == 2
-               call regr_pr_av(ncid_climoz,
-     $              (/"tro3         ", "tro3_daylight"/),
-     $              julien=ro3i, press_in_edg=press_climoz, paprs=paprs,
-     $              v3=wo)
+!              read_climoz == 2
+               call regr_pr_av(ncid_climoz, &
+                    (/"tro3         ", "tro3_daylight"/), &
+                    julien=ro3i, press_in_edg=press_climoz, paprs=paprs, &
+                    v3=wo)
             end if
 !           Convert from mole fraction of ozone to column density of ozone in a
 !           cell, in kDU:
-            forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)
-     $           * rmo3 / rmd * zmasse / dobson_u / 1e3
-C           (By regridding ozone values for LMDZ only once every 360th of
-C           year, we have already neglected the variation of pressure in one
-C           360th of year. So do not recompute "wo" at each time step even if
-C           "zmasse" changes a little.)
+            forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) &
+                 * rmo3 / rmd * zmasse / dobson_u / 1e3
+!           (By regridding ozone values for LMDZ only once every 360th of
+!           year, we have already neglected the variation of pressure in one
+!           360th of year. So do not recompute "wo" at each time step even if
+!           "zmasse" changes a little.)
             co3i = ro3i
          end if
       elseif (MOD(itap-1,lmt_pas) == 0) THEN
-C        Once per day, update ozone from Royer:
+!        Once per day, update ozone from Royer:
          wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
       ENDIF
-c
-c Re-evaporer l'eau liquide nuageuse
-c
+!
+! Re-evaporer l'eau liquide nuageuse
+!
       DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
       DO i = 1, klon
          zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
-cjyg<
-c      Attention : Arnaud a propose des formules completement differentes
-c                  A verifier !!!
+!jyg<
+!      Attention : Arnaud a propose des formules completement differentes
+!                  A verifier !!!
          zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
          IF (iflag_ice_thermo .EQ. 0) THEN
            zlsdcp=zlvdcp
          ENDIF 
-c>jyg
+!>jyg
 
          zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
          zb = MAX(0.0,ql_seri(i,k))
-         za = - MAX(0.0,ql_seri(i,k))
-     .                  * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
+         za = - MAX(0.0,ql_seri(i,k)) &
+                        * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
          t_seri(i,k) = t_seri(i,k) + za
          q_seri(i,k) = q_seri(i,k) + zb
@@ -1899 +1899 @@
       ENDDO
       ENDDO
-cIM
+!IM
       IF (ip_ebil_phy.ge.2) THEN 
         ztit='after reevap'
-        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
-     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
-     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-         call diagphy(airephy,ztit,ip_ebil_phy
-     e      , zero_v, zero_v, zero_v, zero_v, zero_v
-     e      , zero_v, zero_v, zero_v, ztsol
-     e      , d_h_vcol, d_qt, d_ec
-     s      , fs_bound, fq_bound )
-C
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime &
+            , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
+            , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+         call diagphy(airephy,ztit,ip_ebil_phy &
+            , zero_v, zero_v, zero_v, zero_v, zero_v &
+            , zero_v, zero_v, zero_v, ztsol &
+            , d_h_vcol, d_qt, d_ec &
+            , fs_bound, fq_bound )
+!
       END IF 
 
-c
-c=========================================================================
+!
+!=========================================================================
 ! Calculs de l'orbite.
 ! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
@@ -1938 +1938 @@
       endif
       if(prt_level.ge.1)                                                &
-     &    write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
+         write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
 
 
@@ -1967 +1967 @@
       endif
 
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c Appel au pbl_surface : Planetary Boudary Layer et Surface
-c Cela implique tous les interactions des sous-surfaces et la partie diffusion 
-c turbulent du couche limit. 
-c 
-c Certains varibales de sorties de pbl_surface sont utiliser que pour 
-c ecriture des fihiers hist_XXXX.nc, ces sont :
-c   qsol,      zq2m,      s_pblh,  s_lcl,
-c   s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
-c   s_therm,   s_trmb1,   s_trmb2, s_trmb3,
-c   zxrugs,    zu10m,     zv10m,   fder,
-c   zxqsurf,   rh2m,      zxfluxu, zxfluxv,
-c   frugs,     agesno,    fsollw,  fsolsw,
-c   d_ts,      fevap,     fluxlat, t2m,
-c   wfbils,    wfbilo,    fluxt,   fluxu, fluxv,
-c
-c Certains ne sont pas utiliser du tout : 
-c   dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
-c
-
-c Calcul de l'humidite de saturation au niveau du sol
+!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! Appel au pbl_surface : Planetary Boudary Layer et Surface
+! Cela implique tous les interactions des sous-surfaces et la partie diffusion 
+! turbulent du couche limit. 
+! 
+! Certains varibales de sorties de pbl_surface sont utiliser que pour 
+! ecriture des fihiers hist_XXXX.nc, ces sont :
+!   qsol,      zq2m,      s_pblh,  s_lcl,
+!   s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
+!   s_therm,   s_trmb1,   s_trmb2, s_trmb3,
+!   zxrugs,    zu10m,     zv10m,   fder,
+!   zxqsurf,   rh2m,      zxfluxu, zxfluxv,
+!   frugs,     agesno,    fsollw,  fsolsw,
+!   d_ts,      fevap,     fluxlat, t2m,
+!   wfbils,    wfbilo,    fluxt,   fluxu, fluxv,
+!
+! Certains ne sont pas utiliser du tout : 
+!   dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
+!
+
+! Calcul de l'humidite de saturation au niveau du sol
 
 
@@ -1993 +1993 @@
       if (iflag_pbl/=0) then
 
-      CALL pbl_surface( 
-     e     dtime,     date0,     itap,    days_elapsed+1,
-     e     debut,     lafin,
-     e     rlon,      rlat,      rugoro,  rmu0,     
-     e     rain_fall, snow_fall, solsw,   sollw,    
-     e     t_seri,    q_seri,    u_seri,  v_seri,   
-     e     pplay,     paprs,     pctsrf,            
-     +     ftsol,falb1,falb2,ustar,u10m,v10m,wstar,
-     s     sollwdown, cdragh,    cdragm,  u1,    v1,
-     s     albsol1,   albsol2,   sens,    evap,  
-     s     zxtsol,    zxfluxlat, zt2m,    qsat2m, 
-     s     d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf, d_t_diss,
-     s     coefh,     coefm,     slab_wfbils,                
-     d     qsol,      zq2m,      s_pblh,  s_lcl,
-     d     s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
-     d     s_therm,   s_trmb1,   s_trmb2, s_trmb3,
-     d     zxrugs,    zustar, zu10m,     zv10m,   fder,
-     d     zxqsurf,   rh2m,      zxfluxu, zxfluxv,
-     d     frugs,     agesno,    fsollw,  fsolsw,
-     d     d_ts,      fevap,     fluxlat, t2m,
-     d     wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
-     -     dsens,     devap,     zxsnow,
-     -     zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
+      CALL pbl_surface(  &
+           dtime,     date0,     itap,    days_elapsed+1, &
+           debut,     lafin, &
+           rlon,      rlat,      rugoro,  rmu0,      &
+           rain_fall, snow_fall, solsw,   sollw,     &
+           t_seri,    q_seri,    u_seri,  v_seri,    &
+           pplay,     paprs,     pctsrf,             &
+           ftsol,falb1,falb2,ustar,u10m,v10m,wstar, &
+           sollwdown, cdragh,    cdragm,  u1,    v1, &
+           albsol1,   albsol2,   sens,    evap,   &
+           zxtsol,    zxfluxlat, zt2m,    qsat2m,  &
+           d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf, d_t_diss, &
+           coefh,     coefm,     slab_wfbils,                 &
+           qsol,      zq2m,      s_pblh,  s_lcl, &
+           s_capCL,   s_oliqCL,  s_cteiCL,s_pblT, &
+           s_therm,   s_trmb1,   s_trmb2, s_trmb3, &
+           zxrugs,    zustar, zu10m,     zv10m,   fder, &
+           zxqsurf,   rh2m,      zxfluxu, zxfluxv, &
+           frugs,     agesno,    fsollw,  fsolsw, &
+           d_ts,      fevap,     fluxlat, t2m, &
+           wfbils,    wfbilo,    fluxt,   fluxu,  fluxv, &
+           dsens,     devap,     zxsnow, &
+           zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
 
 
 !-----------------------------------------------------------------------------------------
 ! ajout des tendances de la diffusion turbulente
-      CALL add_phys_tend
-     s     (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,'vdf')
+      CALL add_phys_tend &
+           (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,'vdf')
 !-----------------------------------------------------------------------------------------
 
@@ -2031 +2031 @@
       endif
 
-         CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh,
-     e      t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
+         CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
+            t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
 
 
       IF (ip_ebil_phy.ge.2) THEN 
         ztit='after surface_main'
-        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
-     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
-     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-         call diagphy(airephy,ztit,ip_ebil_phy
-     e      , zero_v, zero_v, zero_v, zero_v, sens
-     e      , evap  , zero_v, zero_v, ztsol
-     e      , d_h_vcol, d_qt, d_ec
-     s      , fs_bound, fq_bound )
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+            , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
+            , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+         call diagphy(airephy,ztit,ip_ebil_phy &
+            , zero_v, zero_v, zero_v, zero_v, sens &
+            , evap  , zero_v, zero_v, ztsol &
+            , d_h_vcol, d_qt, d_ec &
+            , fs_bound, fq_bound )
       END IF 
 
       ENDIF
-c =================================================================== c
-c   Calcul de Qsat
+! =================================================================== c
+!   Calcul de Qsat
 
       DO k = 1, klev
@@ -2075 +2075 @@
       write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
       endif
-c
-c Appeler la convection (au choix)
-c
+!
+! Appeler la convection (au choix)
+!
       DO k = 1, klev
       DO i = 1, klon
-         conv_q(i,k) = d_q_dyn(i,k) 
-     .               + d_q_vdf(i,k)/dtime
-         conv_t(i,k) = d_t_dyn(i,k) 
-     .               + d_t_vdf(i,k)/dtime
+         conv_q(i,k) = d_q_dyn(i,k)  &
+                     + d_q_vdf(i,k)/dtime
+         conv_t(i,k) = d_t_dyn(i,k)  &
+                     + d_t_vdf(i,k)/dtime
       ENDDO
       ENDDO
@@ -2098 +2098 @@
          DO k = 1, klev
          DO i = 1, klon
-            z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
-     .                        *(paprs(i,k)-paprs(i,k+1))/RG
+            z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
+                              *(paprs(i,k)-paprs(i,k+1))/RG
          ENDDO
          ENDDO
       ENDIF
 
-c Calcule de vitesse verticale a partir de flux de masse verticale
+! Calcule de vitesse verticale a partir de flux de masse verticale
       DO k = 1, klev
          DO i = 1, klon
@@ -2110 +2110 @@
          END DO
       END DO
-      if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ',
-     $     omega(igout, :)
+      if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
+           omega(igout, :)
 
       IF (iflag_con.EQ.1) THEN
         abort_message ='reactiver le call conlmd dans physiq.F'
         CALL abort_gcm (modname,abort_message,1)
-c     CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
-c    .             d_t_con, d_q_con,
-c    .             rain_con, snow_con, ibas_con, itop_con)
+!     CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
+!    .             d_t_con, d_q_con,
+!    .             rain_con, snow_con, ibas_con, itop_con)
       ELSE IF (iflag_con.EQ.2) THEN
-      CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
-     e            conv_t, conv_q, -evap, omega,
-     s            d_t_con, d_q_con, rain_con, snow_con,
-     s            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
-     s            kcbot, kctop, kdtop, pmflxr, pmflxs)
+      CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
+                  conv_t, conv_q, -evap, omega, &
+                  d_t_con, d_q_con, rain_con, snow_con, &
+                  pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
+                  kcbot, kctop, kdtop, pmflxr, pmflxs)
       d_u_con = 0.
       d_v_con = 0.
@@ -2135 +2135 @@
       ENDDO
       ELSE IF (iflag_con.GE.3) THEN
-c nb of tracers for the KE convection:
-c MAF la partie traceurs est faite dans phytrac
-c on met ntra=1 pour limiter les appels mais on peut
-c supprimer les calculs / ftra.
+! nb of tracers for the KE convection:
+! MAF la partie traceurs est faite dans phytrac
+! on met ntra=1 pour limiter les appels mais on peut
+! supprimer les calculs / ftra.
               ntra = 1
 
-c=====================================================================================
-cajout pour la parametrisation des poches froides: 
-ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri 
+!=====================================================================================
+!ajout pour la parametrisation des poches froides: 
+!calcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri 
       do k=1,klev
             do i=1,klon
              if (iflag_wake>=1) then
-             t_wake(i,k) = t_seri(i,k)
-     .           +(1-wake_s(i))*wake_deltat(i,k)
-             q_wake(i,k) = q_seri(i,k)
-     .           +(1-wake_s(i))*wake_deltaq(i,k)
-             t_undi(i,k) = t_seri(i,k)
-     .           -wake_s(i)*wake_deltat(i,k)
-             q_undi(i,k) = q_seri(i,k)
-     .           -wake_s(i)*wake_deltaq(i,k)
+             t_wake(i,k) = t_seri(i,k) &
+                 +(1-wake_s(i))*wake_deltat(i,k)
+             q_wake(i,k) = q_seri(i,k) &
+                 +(1-wake_s(i))*wake_deltaq(i,k)
+             t_undi(i,k) = t_seri(i,k) &
+                 -wake_s(i)*wake_deltat(i,k)
+             q_undi(i,k) = q_seri(i,k) &
+                 -wake_s(i)*wake_deltaq(i,k)
              else
              t_wake(i,k) = t_seri(i,k)
@@ -2164 +2164 @@
          enddo
       
-cc--   Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
-cc--    pour le soulevement des particules dans le modele convectif
-c
+!c--   Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
+!c--    pour le soulevement des particules dans le modele convectif
+!
       do i = 1,klon
         ALE(i) = 0.
         ALP(i) = 0.
       enddo
-c
-ccalcul de ale_wake et alp_wake
+!
+!calcul de ale_wake et alp_wake
        if (iflag_wake>=1) then
          if (itap .le. it_wape_prescr) then
@@ -2181 +2181 @@
          else
           do i = 1,klon
-cjyg  ALE=WAPE au lieu de ALE = 1/2 Cstar**2
-ccc           ale_wake(i) = 0.5*wake_cstar(i)**2
+!jyg  ALE=WAPE au lieu de ALE = 1/2 Cstar**2
+!cc           ale_wake(i) = 0.5*wake_cstar(i)**2
            ale_wake(i) = wake_pe(i)
            alp_wake(i) = wake_fip(i)
@@ -2193 +2193 @@
          enddo
        endif
-ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
-cdans le thermique sinon
+!combinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
+!dans le thermique sinon
        if (iflag_coupl.eq.0) then
-          if (debut.and.prt_level.gt.9)
-     $                     WRITE(lunout,*)'ALE et ALP imposes'
+          if (debut.and.prt_level.gt.9) &
+                           WRITE(lunout,*)'ALE et ALP imposes'
           do i = 1,klon
-con ne couple que ale
-c           ALE(i) = max(ale_wake(i),Ale_bl(i))
+!on ne couple que ale
+!           ALE(i) = max(ale_wake(i),Ale_bl(i))
             ALE(i) = max(ale_wake(i),ale_bl_prescr)
-con ne couple que alp
-c           ALP(i) = alp_wake(i) + Alp_bl(i)
+!on ne couple que alp
+!           ALP(i) = alp_wake(i) + Alp_bl(i)
             ALP(i) = alp_wake(i) + alp_bl_prescr
           enddo
@@ -2223 +2223 @@
        do i = 1,klon
           ALE(i) = max(ale_wake(i),Ale_bl(i))
-ccc nrlmd le 10/04/2012----------Stochastic triggering--------------
+!cc nrlmd le 10/04/2012----------Stochastic triggering--------------
           if (iflag_trig_bl.ge.1) then
              ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
           endif
-ccc fin nrlmd le 10/04/2012
+!cc fin nrlmd le 10/04/2012
           if (alp_offset>=0.) then
             ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
@@ -2233 +2233 @@
             ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
             if (alp(i)<0.) then
-               print*,'ALP ',alp(i),alp_wake(i)
-     s         ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
+               print*,'ALP ',alp(i),alp_wake(i) &
+               ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
             endif
           endif
@@ -2244 +2244 @@
           if (alp(i)>alp_max) then
              IF(prt_level>9)WRITE(lunout,*)                             &
-     &       'WARNING SUPER ALP (seuil=',alp_max,
-     ,       '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
+            'WARNING SUPER ALP (seuil=',alp_max, &
+             '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
              alp(i)=alp_max
           endif
           if (ale(i)>ale_max) then
              IF(prt_level>9)WRITE(lunout,*)                             &
-     &       'WARNING SUPER ALE (seuil=',ale_max,
-     ,       '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
+            'WARNING SUPER ALE (seuil=',ale_max, &
+             '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
              ale(i)=ale_max
           endif
        enddo
 
-cfin calcul ale et alp
-c=================================================================================================
-
-
-c sb, oct02:
-c Schema de convection modularise et vectorise:
-c (driver commun aux versions 3 et 4)
-c
+!fin calcul ale et alp
+!=================================================================================================
+
+
+! sb, oct02:
+! Schema de convection modularise et vectorise:
+! (driver commun aux versions 3 et 4)
+!
           IF (ok_cvl) THEN ! new driver for convectL
 
@@ -2271 +2271 @@
                 nbtr_tmp=nbtr
              END IF
-cjyg   iflag_con est dans clesphys
-cc          CALL concvl (iflag_con,iflag_clos,
-          CALL concvl (iflag_clos,
-     .        dtime,paprs,pplay,t_undi,q_undi,
-     .        t_wake,q_wake,wake_s,
-     .        u_seri,v_seri,tr_seri,nbtr_tmp,
-     .        ALE,ALP,
-     .        sig1,w01,
-     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
-     .        rain_con, snow_con, ibas_con, itop_con, sigd,
-     .        ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0,
-     .        Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
-     .        pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
+!jyg   iflag_con est dans clesphys
+!c          CALL concvl (iflag_con,iflag_clos,
+          CALL concvl (iflag_clos, &
+              dtime,paprs,pplay,t_undi,q_undi, &
+              t_wake,q_wake,wake_s, &
+              u_seri,v_seri,tr_seri,nbtr_tmp, &
+              ALE,ALP, &
+              sig1,w01, &
+              d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
+              rain_con, snow_con, ibas_con, itop_con, sigd, &
+              ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0, &
+              Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
+              pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
 ! RomP >>>
 !!     .        pmflxr,pmflxs,da,phi,mp,
 !!     .        ftd,fqd,lalim_conv,wght_th)
-     .        pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij,
-     .        ftd,fqd,lalim_conv,wght_th,
-     .        ev, ep,epmlmMm,eplaMm,
-     .        wdtrainA,wdtrainM)
+              pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
+              ftd,fqd,lalim_conv,wght_th, &
+              ev, ep,epmlmMm,eplaMm, &
+              wdtrainA,wdtrainM)
 ! RomP <<<
 
-cIM begin
-c       print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
-c    .dnwd0(1,1),ftd(1,1),fqd(1,1)
-cIM end
-cIM cf. FH
+!IM begin
+!       print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
+!    .dnwd0(1,1),ftd(1,1),fqd(1,1)
+!IM end
+!IM cf. FH
               clwcon0=qcondc
               pmfu(:,:)=upwd(:,:)+dnwd(:,:)
@@ -2307 +2307 @@
           ELSE ! ok_cvl
 
-c MAF conema3 ne contient pas les traceurs
-          CALL conema3 (dtime,
-     .        paprs,pplay,t_seri,q_seri,
-     .        u_seri,v_seri,tr_seri,ntra,
-     .        sig1,w01,
-     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
-     .        rain_con, snow_con, ibas_con, itop_con,
-     .        upwd,dnwd,dnwd0,bas,top,
-     .        Ma,cape,tvp,rflag,
-     .        pbase
-     .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
-     .        ,clwcon0)
+! MAF conema3 ne contient pas les traceurs
+          CALL conema3 (dtime, &
+              paprs,pplay,t_seri,q_seri, &
+              u_seri,v_seri,tr_seri,ntra, &
+              sig1,w01, &
+              d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
+              rain_con, snow_con, ibas_con, itop_con, &
+              upwd,dnwd,dnwd0,bas,top, &
+              Ma,cape,tvp,rflag, &
+              pbase &
+              ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
+              ,clwcon0)
 
           ENDIF ! ok_cvl
 
-c
-c Correction precip
+!
+! Correction precip
           rain_con = rain_con * cvl_corr
           snow_con = snow_con * cvl_corr
-c
+!
 
            IF (.NOT. ok_gust) THEN
@@ -2334 +2334 @@
            ENDIF
 
-c =================================================================== c
-c Calcul des proprietes des nuages convectifs
-c
-
-c   calcul des proprietes des nuages convectifs
+! =================================================================== c
+! Calcul des proprietes des nuages convectifs
+!
+
+!   calcul des proprietes des nuages convectifs
              clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
-             call clouds_gno
-     s       (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
-
-c =================================================================== c
+             call clouds_gno &
+             (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
+
+! =================================================================== c
 
           DO i = 1, klon
@@ -2382 +2382 @@
       ENDIF
 
-c     CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
-c    .              d_u_con, d_v_con)
+!     CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
+!    .              d_u_con, d_v_con)
 
 !-----------------------------------------------------------------------------------------
@@ -2397 +2397 @@
       endif
 
-cIM
+!IM
       IF (ip_ebil_phy.ge.2) THEN 
         ztit='after convect'
-        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
-     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
-     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-         call diagphy(airephy,ztit,ip_ebil_phy
-     e      , zero_v, zero_v, zero_v, zero_v, zero_v
-     e      , zero_v, rain_con, snow_con, ztsol
-     e      , d_h_vcol, d_qt, d_ec
-     s      , fs_bound, fq_bound )
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+            , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
+            , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+         call diagphy(airephy,ztit,ip_ebil_phy &
+            , zero_v, zero_v, zero_v, zero_v, zero_v &
+            , zero_v, rain_con, snow_con, ztsol &
+            , d_h_vcol, d_qt, d_ec &
+            , fs_bound, fq_bound )
       END IF 
-C
+!
       IF (check) THEN
           za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
@@ -2417 +2417 @@
           DO i = 1, klon
             za = za + airephy(i)/REAL(klon)
-            zx_t = zx_t + (rain_con(i)+
-     .                   snow_con(i))*airephy(i)/REAL(klon)
+            zx_t = zx_t + (rain_con(i)+ &
+                         snow_con(i))*airephy(i)/REAL(klon)
           ENDDO
           zx_t = zx_t/za*dtime
@@ -2429 +2429 @@
           DO k = 1, klev
             DO i = 1, klon
-              z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
-     .            *(paprs(i,k)-paprs(i,k+1))/RG
+              z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
+                  *(paprs(i,k)-paprs(i,k+1))/RG
             ENDDO
           ENDDO
           DO i = 1, klon
-            z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
-     .          /z_apres(i)
+            z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
+                /z_apres(i)
           ENDDO
           DO k = 1, klev
             DO i = 1, klon
-              IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
-     .            z_factor(i).LT.(1.0-1.0E-08)) THEN
+              IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
+                  z_factor(i).LT.(1.0-1.0E-08)) THEN
                   q_seri(i,k) = q_seri(i,k) * z_factor(i)
               ENDIF
@@ -2448 +2448 @@
       zx_ajustq=.FALSE.
 
-c
-c=============================================================================
-cRR:Evolution de la poche froide: on ne fait pas de separation wake/env 
-cpour la couche limite diffuse pour l instant
-c
+!
+!=============================================================================
+!RR:Evolution de la poche froide: on ne fait pas de separation wake/env 
+!pour la couche limite diffuse pour l instant
+!
       if (iflag_wake>=1) then
       DO k=1,klev
@@ -2472 +2472 @@
         ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
         DO k = 1,klev
-         dt_dwn(:,k)= dt_dwn(:,k)+
-     :            ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
-         dq_dwn(:,k)= dq_dwn(:,k)+
-     :            ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
+         dt_dwn(:,k)= dt_dwn(:,k)+ &
+                  ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
+         dq_dwn(:,k)= dq_dwn(:,k)+ &
+                  ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
         ENDDO
       endif
-c
-ccalcul caracteristiques de la poche froide
-      call calWAKE (paprs,pplay,dtime
-     :               ,t_seri,q_seri,omega
-     :               ,dt_dwn,dq_dwn,M_dwn,M_up
-     :               ,dt_a,dq_a,sigd
-     :               ,wdt_PBL,wdq_PBL
-     :               ,udt_PBL,udq_PBL
-     o               ,wake_deltat,wake_deltaq,wake_dth
-     o               ,wake_h,wake_s,wake_dens
-     o               ,wake_pe,wake_fip,wake_gfl
-     o               ,dt_wake,dq_wake
-     o               ,wake_k, t_undi,q_undi
-     o               ,wake_omgbdth,wake_dp_omgb
-     o               ,wake_dtKE,wake_dqKE
-     o               ,wake_dtPBL,wake_dqPBL
-     o               ,wake_omg,wake_dp_deltomg
-     o               ,wake_spread,wake_Cstar,wake_d_deltat_gw
-     o               ,wake_ddeltat,wake_ddeltaq)
-c
+!
+!calcul caracteristiques de la poche froide
+      call calWAKE (paprs,pplay,dtime &
+                     ,t_seri,q_seri,omega &
+                     ,dt_dwn,dq_dwn,M_dwn,M_up &
+                     ,dt_a,dq_a,sigd &
+                     ,wdt_PBL,wdq_PBL &
+                     ,udt_PBL,udq_PBL &
+                     ,wake_deltat,wake_deltaq,wake_dth &
+                     ,wake_h,wake_s,wake_dens &
+                     ,wake_pe,wake_fip,wake_gfl &
+                     ,dt_wake,dq_wake &
+                     ,wake_k, t_undi,q_undi &
+                     ,wake_omgbdth,wake_dp_omgb &
+                     ,wake_dtKE,wake_dqKE &
+                     ,wake_dtPBL,wake_dqPBL &
+                     ,wake_omg,wake_dp_deltomg &
+                     ,wake_spread,wake_Cstar,wake_d_deltat_gw &
+                     ,wake_ddeltat,wake_ddeltaq)
+!
 !-----------------------------------------------------------------------------------------
 ! ajout des tendances des poches froides
@@ -2508 +2508 @@
 
       endif
-c
-c===================================================================
-cJYG
+!
+!===================================================================
+!JYG
       IF (ip_ebil_phy.ge.2) THEN 
         ztit='after wake'
-        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
-     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
-     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-        call diagphy(airephy,ztit,ip_ebil_phy
-     e      , zero_v, zero_v, zero_v, zero_v, zero_v
-     e      , zero_v, zero_v, zero_v, ztsol
-     e      , d_h_vcol, d_qt, d_ec
-     s      , fs_bound, fq_bound )
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+            , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
+            , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+        call diagphy(airephy,ztit,ip_ebil_phy &
+            , zero_v, zero_v, zero_v, zero_v, zero_v &
+            , zero_v, zero_v, zero_v, ztsol &
+            , d_h_vcol, d_qt, d_ec &
+            , fs_bound, fq_bound )
       END IF 
 
-c      print*,'apres callwake iflag_cldcon=', iflag_cldcon
-c
-c===================================================================
-c Convection seche (thermiques ou ajustement)
-c===================================================================
-c
-       call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
-     s ,seuil_inversion,weak_inversion,dthmin)
+!      print*,'apres callwake iflag_cldcon=', iflag_cldcon
+!
+!===================================================================
+! Convection seche (thermiques ou ajustement)
+!===================================================================
+!
+       call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
+       ,seuil_inversion,weak_inversion,dthmin)
 
 
@@ -2541 +2541 @@
       d_q_ajs(:,:)=0.
       clwcon0th(:,:)=0.
-c
-c      fm_therm(:,:)=0.
-c      entr_therm(:,:)=0.
-c      detr_therm(:,:)=0.
-c
-      IF(prt_level>9)WRITE(lunout,*)
-     .    'AVANT LA CONVECTION SECHE , iflag_thermals='
-     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
+!
+!      fm_therm(:,:)=0.
+!      entr_therm(:,:)=0.
+!      detr_therm(:,:)=0.
+!
+      IF(prt_level>9)WRITE(lunout,*) &
+          'AVANT LA CONVECTION SECHE , iflag_thermals=' &
+         ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
       if(iflag_thermals.lt.0) then
-c  Rien
-c  ====
+!  Rien
+!  ====
          IF(prt_level>9)WRITE(lunout,*)'pas de convection'
 
@@ -2557 +2557 @@
       else
 
-c  Thermiques
-c  ==========
-         IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
-     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
-
-
-ccc nrlmd le 10/04/2012
+!  Thermiques
+!  ==========
+         IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
+         ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
+
+
+!cc nrlmd le 10/04/2012
          DO k=1,klev+1
            DO i=1,klon
@@ -2572 +2572 @@
            ENDDO
          ENDDO
-ccc fin nrlmd le 10/04/2012
+!cc fin nrlmd le 10/04/2012
 
          if (iflag_thermals>=1) then
-         call calltherm(pdtphys
-     s      ,pplay,paprs,pphi,weak_inversion
-     s      ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
-     s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
-     s      ,fm_therm,entr_therm,detr_therm
-     s      ,zqasc,clwcon0th,lmax_th,ratqscth
-     s      ,ratqsdiff,zqsatth
-con rajoute ale et alp, et les caracteristiques de la couche alim
-     s      ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca
-     s      ,ztv,zpspsk,ztla,zthl
-ccc nrlmd le 10/04/2012
-     e      ,pbl_tke_input,pctsrf,omega,airephy
-     s      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0
-     s      ,n2,s2,ale_bl_stat
-     s      ,therm_tke_max,env_tke_max
-     s      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke
-     s      ,alp_bl_conv,alp_bl_stat
-ccc fin nrlmd le 10/04/2012
-     s      ,zqla,ztva )
-
-ccc nrlmd le 10/04/2012
-c-----------Stochastic triggering-----------
+         call calltherm(pdtphys &
+            ,pplay,paprs,pphi,weak_inversion &
+            ,u_seri,v_seri,t_seri,q_seri,zqsat,debut &
+            ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
+            ,fm_therm,entr_therm,detr_therm &
+            ,zqasc,clwcon0th,lmax_th,ratqscth &
+            ,ratqsdiff,zqsatth &
+!on rajoute ale et alp, et les caracteristiques de la couche alim
+            ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
+            ,ztv,zpspsk,ztla,zthl &
+!cc nrlmd le 10/04/2012
+            ,pbl_tke_input,pctsrf,omega,airephy &
+            ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
+            ,n2,s2,ale_bl_stat &
+            ,therm_tke_max,env_tke_max &
+            ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
+            ,alp_bl_conv,alp_bl_stat &
+!cc fin nrlmd le 10/04/2012
+            ,zqla,ztva )
+
+!cc nrlmd le 10/04/2012
+!-----------Stochastic triggering-----------
       if (iflag_trig_bl.ge.1) then
-c
+!
         IF (prt_level .GE. 10) THEN
-         print *,'cin, ale_bl_stat, alp_bl_stat ',
-     $            cin, ale_bl_stat, alp_bl_stat
+         print *,'cin, ale_bl_stat, alp_bl_stat ', &
+                  cin, ale_bl_stat, alp_bl_stat
         ENDIF
 
-c----Initialisations
+!----Initialisations
       do i=1,klon
       proba_notrig(i)=1.
@@ -2614 +2614 @@
         endif
       enddo
-c
+!
         IF (prt_level .GE. 10) THEN
-         print *,'random_notrig, tau_trig ',
-     $            random_notrig, tau_trig
-          print *,'s_trig,s2,n2 ',
-     $             s_trig,s2,n2
+         print *,'random_notrig, tau_trig ', &
+                  random_notrig, tau_trig
+          print *,'s_trig,s2,n2 ', &
+                   s_trig,s2,n2
         ENDIF
 
-c----Tirage al\'eatoire et calcul de ale_bl_trig
+!----Tirage al\'eatoire et calcul de ale_bl_trig
       do i=1,klon
         if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) )  then
-        proba_notrig(i)=(1.-exp(-s_trig/s2(i)))**
-     $                  (n2(i)*dtime/tau_trig(i))
-c        print *, 'proba_notrig(i) ',proba_notrig(i)
+        proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
+                        (n2(i)*dtime/tau_trig(i))
+!        print *, 'proba_notrig(i) ',proba_notrig(i)
           if (random_notrig(i) .ge. proba_notrig(i)) then 
           ale_bl_trig(i)=ale_bl_stat(i)
@@ -2639 +2639 @@
         endif
       enddo
-c
+!
         IF (prt_level .GE. 10) THEN
-         print *,'proba_notrig, ale_bl_trig ',
-     $            proba_notrig, ale_bl_trig
+         print *,'proba_notrig, ale_bl_trig ', &
+                  proba_notrig, ale_bl_trig
         ENDIF
 
       endif !(iflag_trig_bl)
 
-c-----------Statistical closure-----------
+!-----------Statistical closure-----------
       if (iflag_clos_bl.ge.1) then 
 
@@ -2664 +2664 @@
         ENDIF
 
-ccc fin nrlmd le 10/04/2012
+!cc fin nrlmd le 10/04/2012
 
 ! ----------------------------------------------------------------------
@@ -2692 +2692 @@
 
 
-c  Ajustement sec
-c  ==============
+!  Ajustement sec
+!  ==============
 
 ! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
@@ -2715 +2715 @@
 
          if (iflag_thermals.eq.0) then
-            CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
-     s      , d_t_ajsb, d_q_ajsb)
+            CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
+            , d_t_ajsb, d_q_ajsb)
          else
-            CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
-     s      , d_t_ajsb, d_q_ajsb)
+            CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
+            , d_t_ajsb, d_q_ajsb)
          endif
 
@@ -2733 +2733 @@
 
       endif
-c
-c===================================================================
-cIM
+!
+!===================================================================
+!IM
       IF (ip_ebil_phy.ge.2) THEN 
         ztit='after dry_adjust'
-        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
-     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
-     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-        call diagphy(airephy,ztit,ip_ebil_phy
-     e      , zero_v, zero_v, zero_v, zero_v, zero_v
-     e      , zero_v, zero_v, zero_v, ztsol
-     e      , d_h_vcol, d_qt, d_ec
-     s      , fs_bound, fq_bound )
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+            , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
+            , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+        call diagphy(airephy,ztit,ip_ebil_phy &
+            , zero_v, zero_v, zero_v, zero_v, zero_v &
+            , zero_v, zero_v, zero_v, ztsol &
+            , d_h_vcol, d_qt, d_ec &
+            , fs_bound, fq_bound )
       END IF 
 
 
-c-------------------------------------------------------------------------
+!-------------------------------------------------------------------------
 ! Computation of ratqs, the width (normalized) of the subrid scale 
 ! water distribution
-      CALL  calcratqs(klon,klev,prt_level,lunout,       
-     s     iflag_ratqs,iflag_con,iflag_cldcon,pdtphys, 
-     s     ratqsbas,ratqshaut,tau_ratqs,fact_cldcon,  
-     s     ptconv,ptconvth,clwcon0th, rnebcon0th,    
-     s     paprs,pplay,q_seri,zqsat,fm_therm,
-     s     ratqs,ratqsc)
-
-
-c
-c Appeler le processus de condensation a grande echelle
-c et le processus de precipitation
-c-------------------------------------------------------------------------
+      CALL  calcratqs(klon,klev,prt_level,lunout,        &
+           iflag_ratqs,iflag_con,iflag_cldcon,pdtphys,  &
+           ratqsbas,ratqshaut,tau_ratqs,fact_cldcon,   &
+           ptconv,ptconvth,clwcon0th, rnebcon0th,     &
+           paprs,pplay,q_seri,zqsat,fm_therm, &
+           ratqs,ratqsc)
+
+
+!
+! Appeler le processus de condensation a grande echelle
+! et le processus de precipitation
+!-------------------------------------------------------------------------
       IF (prt_level .GE.10) THEN
        print *,' ->fisrtilp '
       ENDIF
-c-------------------------------------------------------------------------
-      CALL fisrtilp(dtime,paprs,pplay,
-     .           t_seri, q_seri,ptconv,ratqs,
-     .           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
-     .           rain_lsc, snow_lsc,
-     .           pfrac_impa, pfrac_nucl, pfrac_1nucl,
-     .           frac_impa, frac_nucl, beta_prec_fisrt,
-     .           prfl, psfl, rhcl, 
-     .           zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon,
-     .           iflag_ice_thermo)
+!-------------------------------------------------------------------------
+      CALL fisrtilp(dtime,paprs,pplay, &
+                 t_seri, q_seri,ptconv,ratqs, &
+                 d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, &
+                 rain_lsc, snow_lsc, &
+                 pfrac_impa, pfrac_nucl, pfrac_1nucl, &
+                 frac_impa, frac_nucl, beta_prec_fisrt, &
+                 prfl, psfl, rhcl,  &
+                 zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon, &
+                 iflag_ice_thermo)
 
       WHERE (rain_lsc < 0) rain_lsc = 0.
@@ -2797 +2797 @@
          DO i = 1, klon
             za = za + airephy(i)/REAL(klon)
-            zx_t = zx_t + (rain_lsc(i)
-     .                  + snow_lsc(i))*airephy(i)/REAL(klon)
+            zx_t = zx_t + (rain_lsc(i) &
+                        + snow_lsc(i))*airephy(i)/REAL(klon)
         ENDDO
          zx_t = zx_t/za*dtime
          WRITE(lunout,*)"Precip=", zx_t
       ENDIF
-cIM
+!IM
       IF (ip_ebil_phy.ge.2) THEN 
         ztit='after fisrt'
-        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
-     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
-     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-        call diagphy(airephy,ztit,ip_ebil_phy
-     e      , zero_v, zero_v, zero_v, zero_v, zero_v
-     e      , zero_v, rain_lsc, snow_lsc, ztsol
-     e      , d_h_vcol, d_qt, d_ec
-     s      , fs_bound, fq_bound )
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+            , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
+            , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+        call diagphy(airephy,ztit,ip_ebil_phy &
+            , zero_v, zero_v, zero_v, zero_v, zero_v &
+            , zero_v, rain_lsc, snow_lsc, ztsol &
+            , d_h_vcol, d_qt, d_ec &
+            , fs_bound, fq_bound )
       END IF 
 
@@ -2823 +2823 @@
       endif
 
-c
-c-------------------------------------------------------------------
-c  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
-c-------------------------------------------------------------------
-
-c 1. NUAGES CONVECTIFS
-c
-cIM cf FH
-c     IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
+!
+!-------------------------------------------------------------------
+!  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
+!-------------------------------------------------------------------
+
+! 1. NUAGES CONVECTIFS
+!
+!IM cf FH
+!     IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
       IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
        snow_tiedtke=0.
-c     print*,'avant calcul de la pseudo precip '
-c     print*,'iflag_cldcon',iflag_cldcon
+!     print*,'avant calcul de la pseudo precip '
+!     print*,'iflag_cldcon',iflag_cldcon
        if (iflag_cldcon.eq.-1) then
           rain_tiedtke=rain_con
        else
-c       print*,'calcul de la pseudo precip '
+!       print*,'calcul de la pseudo precip '
           rain_tiedtke=0.
-c         print*,'calcul de la pseudo precip 0'
+!         print*,'calcul de la pseudo precip 0'
           do k=1,klev
           do i=1,klon
              if (d_q_con(i,k).lt.0.) then
-                rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
-     s         *(paprs(i,k)-paprs(i,k+1))/rg
+                rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
+               *(paprs(i,k)-paprs(i,k+1))/rg
              endif
           enddo
           enddo
        endif
-c
-c     call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
-c
-
-c Nuages diagnostiques pour Tiedtke
-      CALL diagcld1(paprs,pplay,
-cIM cf FH  .             rain_con,snow_con,ibas_con,itop_con,
-     .             rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
-     .             diafra,dialiq)
+!
+!     call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
+!
+
+! Nuages diagnostiques pour Tiedtke
+      CALL diagcld1(paprs,pplay, &
+!IM cf FH  .             rain_con,snow_con,ibas_con,itop_con,
+                   rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
+                   diafra,dialiq)
       DO k = 1, klev
       DO i = 1, klon
@@ -2870 +2870 @@
 
       ELSE IF (iflag_cldcon.ge.3) THEN
-c  On prend pour les nuages convectifs le max du calcul de la
-c  convection et du calcul du pas de temps precedent diminue d'un facteur
-c  facttemps
+!  On prend pour les nuages convectifs le max du calcul de la
+!  convection et du calcul du pas de temps precedent diminue d'un facteur
+!  facttemps
       facteur = pdtphys *facttemps
       do k=1,klev
          do i=1,klon
             rnebcon(i,k)=rnebcon(i,k)*facteur
-            if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k))
-     s      then
+            if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k)) &
+            then
                 rnebcon(i,k)=rnebcon0(i,k)
                 clwcon(i,k)=clwcon0(i,k)
@@ -2885 +2885 @@
       enddo
 
-c
-cjq - introduce the aerosol direct and first indirect radiative forcings
-cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
+!
+!jq - introduce the aerosol direct and first indirect radiative forcings
+!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
       IF (flag_aerosol .gt. 0) THEN
-         IF (.NOT. aerosol_couple)
-     &        CALL readaerosol_optic(
-     &        debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref,
-     &        pdtphys, pplay, paprs, t_seri, rhcl, presnivs, 
-     &        mass_solu_aero, mass_solu_aero_pi, 
-     &        tau_aero, piz_aero, cg_aero, 
-     &        tausum_aero, tau3d_aero)
+         IF (.NOT. aerosol_couple) &
+              CALL readaerosol_optic( &
+              debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
+              pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
+              mass_solu_aero, mass_solu_aero_pi,  &
+              tau_aero, piz_aero, cg_aero,  &
+              tausum_aero, tau3d_aero)
       ELSE
          tausum_aero(:,:,:) = 0.
@@ -2902 +2902 @@
          cg_aero(:,:,:,:)  = 0.
       ENDIF
-c
-c--STRAT AEROSOL
-c--updates tausum_aero,tau_aero,piz_aero,cg_aero
+!
+!--STRAT AEROSOL
+!--updates tausum_aero,tau_aero,piz_aero,cg_aero
       IF (flag_aerosol_strat) THEN
          PRINT *,'appel a readaerosolstrat', mth_cur
          CALL readaerosolstrato(debut)
       ENDIF
-c--fin STRAT AEROSOL
-
-cIM calcul nuages par le simulateur ISCCP
-c
+!--fin STRAT AEROSOL
+
+!IM calcul nuages par le simulateur ISCCP
+!
 #ifdef histISCCP
       IF (ok_isccp) THEN
-c
-cIM lecture invtau, tautab des fichiers formattes
-c
+!
+!IM lecture invtau, tautab des fichiers formattes
+!
       IF (debut) THEN
-c$OMP MASTER
-c
+!$OMP MASTER
+!
       open(99,file='tautab.formatted', FORM='FORMATTED')
       read(99,'(f30.20)') tautab_omp
       close(99)
-c
+!
       open(99,file='invtau.formatted',form='FORMATTED')
       read(99,'(i10)') invtau_omp
 
-c     print*,'calcul_simulISCCP invtau_omp',invtau_omp
-c     write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
+!     print*,'calcul_simulISCCP invtau_omp',invtau_omp
+!     write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
 
       close(99)
-c$OMP END MASTER
-c$OMP BARRIER 
+!$OMP END MASTER
+!$OMP BARRIER 
       tautab=tautab_omp
       invtau=invtau_omp
-c
+!
       ENDIF !debut
-c
-cIM appel simulateur toutes les  NINT(freq_ISCCP/dtime) heures
+!
+!IM appel simulateur toutes les  NINT(freq_ISCCP/dtime) heures
        IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
 #include "calcul_simulISCCP.h"
@@ -2946 +2946 @@
 #endif
 
-c   On prend la somme des fractions nuageuses et des contenus en eau
+!   On prend la somme des fractions nuageuses et des contenus en eau
 
       if (iflag_cldcon>=5) then
@@ -3015 +3015 @@
 !        enddo prfl, psfl,
 !     enddo
-c
-c 2. NUAGES STARTIFORMES
-c
+!
+! 2. NUAGES STARTIFORMES
+!
       IF (ok_stratus) THEN
       CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
@@ -3029 +3029 @@
       ENDDO
       ENDIF
-c
-c Precipitation totale
-c
+!
+! Precipitation totale
+!
       DO i = 1, klon
          rain_fall(i) = rain_con(i) + rain_lsc(i)
          snow_fall(i) = snow_con(i) + snow_lsc(i)
       ENDDO
-cIM
+!IM
       IF (ip_ebil_phy.ge.2) THEN 
         ztit="after diagcld"
-        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
-     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
-     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-        call diagphy(airephy,ztit,ip_ebil_phy
-     e      , zero_v, zero_v, zero_v, zero_v, zero_v
-     e      , zero_v, zero_v, zero_v, ztsol
-     e      , d_h_vcol, d_qt, d_ec
-     s      , fs_bound, fq_bound )
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+            , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
+            , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+        call diagphy(airephy,ztit,ip_ebil_phy &
+            , zero_v, zero_v, zero_v, zero_v, zero_v &
+            , zero_v, zero_v, zero_v, ztsol &
+            , d_h_vcol, d_qt, d_ec &
+            , fs_bound, fq_bound )
       END IF 
-c
-c Calculer l'humidite relative pour diagnostique
-c
+!
+! Calculer l'humidite relative pour diagnostique
+!
       DO k = 1, klev
       DO i = 1, klon
@@ -3072 +3072 @@
       ENDDO
 
-cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle 
-c   equivalente a 2m (tpote) pour diagnostique
-c
+!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle 
+!   equivalente a 2m (tpote) pour diagnostique
+!
       DO i = 1, klon
        tpot(i)=zt2m(i)*(100000./paprs(i,1))**RKAPPA
@@ -3090 +3090 @@
         ENDIF
        ENDIF
-       tpote(i) = tpot(i)*     
-     . EXP((Lheat *qsat2m(i))/(RCPD*zt2m(i)))
+       tpote(i) = tpot(i)*      &
+       EXP((Lheat *qsat2m(i))/(RCPD*zt2m(i)))
       ENDDO
 
@@ -3102 +3102 @@
          call chemtime(itap+itau_phy-1, date0, dtime)
          IF (config_inca == 'aero') THEN
-            CALL AEROSOL_METEO_CALC(
-     $           calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
-     $           prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
+            CALL AEROSOL_METEO_CALC( &
+                 calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
+                 prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
          END IF
 
          zxsnow_dummy(:) = 0.0
 
-         CALL chemhook_begin (calday,
-     $                          days_elapsed+1,
-     $                          jH_cur,
-     $                          pctsrf(1,1),
-     $                          rlat,
-     $                          rlon,
-     $                          airephy,
-     $                          paprs,
-     $                          pplay,
-     $                          coefh(:,:,is_ave),
-     $                          pphi,
-     $                          t_seri,
-     $                          u,
-     $                          v,
-     $                          wo(:, :, 1),
-     $                          q_seri,
-     $                          zxtsol,
-     $                          zxsnow_dummy,
-     $                          solsw,
-     $                          albsol1,
-     $                          rain_fall,
-     $                          snow_fall,
-     $                          itop_con,
-     $                          ibas_con,
-     $                          cldfra,
-     $                          iim,
-     $                          jjm,
-     $                          tr_seri,
-     $                          ftsol,
-     $                          paprs,
-     $                          cdragh,
-     $                          cdragm,
-     $                          pctsrf,
-     $                          pdtphys,
-     $                            itap)
+         CALL chemhook_begin (calday, &
+                                days_elapsed+1, &
+                                jH_cur, &
+                                pctsrf(1,1), &
+                                rlat, &
+                                rlon, &
+                                airephy, &
+                                paprs, &
+                                pplay, &
+                                coefh(:,:,is_ave), &
+                                pphi, &
+                                t_seri, &
+                                u, &
+                                v, &
+                                wo(:, :, 1), &
+                                q_seri, &
+                                zxtsol, &
+                                zxsnow_dummy, &
+                                solsw, &
+                                albsol1, &
+                                rain_fall, &
+                                snow_fall, &
+                                itop_con, &
+                                ibas_con, &
+                                cldfra, &
+                                iim, &
+                                jjm, &
+                                tr_seri, &
+                                ftsol, &
+                                paprs, &
+                                cdragh, &
+                                cdragm, &
+                                pctsrf, &
+                                pdtphys, &
+                                  itap)
 
          CALL VTe(VTinca)
@@ -3149 +3149 @@
 #endif 
       END IF !type_trac = inca
-c     
-c Calculer les parametres optiques des nuages et quelques
-c parametres pour diagnostiques:
-c
+!     
+! Calculer les parametres optiques des nuages et quelques
+! parametres pour diagnostiques:
+!
 
       IF (aerosol_couple) THEN 
@@ -3160 +3160 @@
 
       if (ok_newmicro) then
-      CALL newmicro (ok_cdnc, bl95_b0, bl95_b1,
-     .              paprs, pplay, t_seri, cldliq, cldfra,
-     .              cldtau, cldemi, cldh, cldl, cldm, cldt, cldq,
-     e              flwp, fiwp, flwc, fiwc,
-     e              mass_solu_aero, mass_solu_aero_pi,
-     s              cldtaupi, re, fl, ref_liq, ref_ice) 
+      CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
+                    paprs, pplay, t_seri, cldliq, cldfra, &
+                    cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
+                    flwp, fiwp, flwc, fiwc, &
+                    mass_solu_aero, mass_solu_aero_pi, &
+                    cldtaupi, re, fl, ref_liq, ref_ice) 
       else
-      CALL nuage (paprs, pplay,
-     .            t_seri, cldliq, cldfra, cldtau, cldemi,
-     .            cldh, cldl, cldm, cldt, cldq,
-     e            ok_aie,
-     e            mass_solu_aero, mass_solu_aero_pi,
-     e            bl95_b0, bl95_b1,
-     s            cldtaupi, re, fl)
+      CALL nuage (paprs, pplay, &
+                  t_seri, cldliq, cldfra, cldtau, cldemi, &
+                  cldh, cldl, cldm, cldt, cldq, &
+                  ok_aie, &
+                  mass_solu_aero, mass_solu_aero_pi, &
+                  bl95_b0, bl95_b1, &
+                  cldtaupi, re, fl)
       endif
-c
-cIM betaCRF
-c
+!
+!IM betaCRF
+!
       cldtaurad   = cldtau
       cldtaupirad = cldtaupi
       cldemirad   = cldemi
      
-c
-      if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND.
-     $lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
-c
-c global
-c
+!
+      if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
+      lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
+!
+! global
+!
        DO k=1, klev
        DO i=1, klon
@@ -3204 +3204 @@
        ENDDO
        ENDDO 
-c
+!
       else
-c
-c regional
-c
+!
+! regional
+!
        DO k=1, klev
        DO i=1,klon
-c
-        if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND.
-     $      rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
+!
+        if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND. &
+            rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
          if (pplay(i,k).GE.pfree) THEN
           beta(i,k) = beta_pbl
@@ -3227 +3227 @@
         cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
         endif
-c
+!
        ENDDO
        ENDDO
-c
+!
       endif
-c
-c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
-c
+!
+! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
+!
       IF (MOD(itaprad,radpas).EQ.0) THEN
 
       DO i = 1, klon
-         albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
-     .             + falb1(i,is_lic) * pctsrf(i,is_lic)
-     .             + falb1(i,is_ter) * pctsrf(i,is_ter)
-     .             + falb1(i,is_sic) * pctsrf(i,is_sic)
-         albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
-     .               + falb2(i,is_lic) * pctsrf(i,is_lic)
-     .               + falb2(i,is_ter) * pctsrf(i,is_ter)
-     .               + falb2(i,is_sic) * pctsrf(i,is_sic)
+         albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce) &
+                   + falb1(i,is_lic) * pctsrf(i,is_lic) &
+                   + falb1(i,is_ter) * pctsrf(i,is_ter) &
+                   + falb1(i,is_sic) * pctsrf(i,is_sic)
+         albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce) &
+                     + falb2(i,is_lic) * pctsrf(i,is_lic) &
+                     + falb2(i,is_ter) * pctsrf(i,is_ter) &
+                     + falb2(i,is_sic) * pctsrf(i,is_sic)
       ENDDO
 
@@ -3257 +3257 @@
       IF (aerosol_couple) THEN 
 #ifdef INCA
-         CALL radlwsw_inca 
-     e        (kdlon,kflev,dist, rmu0, fract, solaire,
-     e        paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
-     e        wo(:, :, 1),
-     e        cldfrarad, cldemirad, cldtaurad,
-     s        heat,heat0,cool,cool0,radsol,albpla,
-     s        topsw,toplw,solsw,sollw,
-     s        sollwdown,
-     s        topsw0,toplw0,solsw0,sollw0,
-     s        lwdn0, lwdn, lwup0, lwup, 
-     s        swdn0, swdn, swup0, swup,
-     e        ok_ade, ok_aie,
-     e        tau_aero, piz_aero, cg_aero,
-     s        topswad_aero, solswad_aero,
-     s        topswad0_aero, solswad0_aero,
-     s        topsw_aero, topsw0_aero,
-     s        solsw_aero, solsw0_aero,
-     e        cldtaupirad,
-     s        topswai_aero, solswai_aero)
+         CALL radlwsw_inca  &
+              (kdlon,kflev,dist, rmu0, fract, solaire, &
+              paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
+              wo(:, :, 1), &
+              cldfrarad, cldemirad, cldtaurad, &
+              heat,heat0,cool,cool0,radsol,albpla, &
+              topsw,toplw,solsw,sollw, &
+              sollwdown, &
+              topsw0,toplw0,solsw0,sollw0, &
+              lwdn0, lwdn, lwup0, lwup,  &
+              swdn0, swdn, swup0, swup, &
+              ok_ade, ok_aie, &
+              tau_aero, piz_aero, cg_aero, &
+              topswad_aero, solswad_aero, &
+              topswad0_aero, solswad0_aero, &
+              topsw_aero, topsw0_aero, &
+              solsw_aero, solsw0_aero, &
+              cldtaupirad, &
+              topswai_aero, solswai_aero)
             
 #endif
       ELSE
-c
-cIM calcul radiatif pour le cas actuel
-c
+!
+!IM calcul radiatif pour le cas actuel
+!
        RCO2 = RCO2_act
        RCH4 = RCH4_act
@@ -3287 +3287 @@
        RCFC11 = RCFC11_act
        RCFC12 = RCFC12_act
-c
+!
       IF (prt_level .GE.10) THEN
        print *,' ->radlwsw, number 1 '
       ENDIF
-c
-         CALL radlwsw
-     e        (dist, rmu0, fract, 
-     e        paprs, pplay,zxtsol,albsol1, albsol2, 
-     e        t_seri,q_seri,wo,
-     e        cldfrarad, cldemirad, cldtaurad,
-     e        ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol,
-     e        flag_aerosol_strat,
-     e        tau_aero, piz_aero, cg_aero,
-     e        cldtaupirad,new_aod,
-     e        zqsat, flwc, fiwc,
-     s        heat,heat0,cool,cool0,radsol,albpla,
-     s        topsw,toplw,solsw,sollw,
-     s        sollwdown,
-     s        topsw0,toplw0,solsw0,sollw0,
-     s        lwdn0, lwdn, lwup0, lwup, 
-     s        swdn0, swdn, swup0, swup,
-     s        topswad_aero, solswad_aero,
-     s        topswai_aero, solswai_aero,
-     o        topswad0_aero, solswad0_aero,
-     o        topsw_aero, topsw0_aero,
-     o        solsw_aero, solsw0_aero,
-     o        topswcf_aero, solswcf_aero)
+!
+         CALL radlwsw &
+              (dist, rmu0, fract,  &
+              paprs, pplay,zxtsol,albsol1, albsol2,  &
+              t_seri,q_seri,wo, &
+              cldfrarad, cldemirad, cldtaurad, &
+              ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
+              flag_aerosol_strat, &
+              tau_aero, piz_aero, cg_aero, &
+              cldtaupirad,new_aod, &
+              zqsat, flwc, fiwc, &
+              heat,heat0,cool,cool0,radsol,albpla, &
+              topsw,toplw,solsw,sollw, &
+              sollwdown, &
+              topsw0,toplw0,solsw0,sollw0, &
+              lwdn0, lwdn, lwup0, lwup,  &
+              swdn0, swdn, swup0, swup, &
+              topswad_aero, solswad_aero, &
+              topswai_aero, solswai_aero, &
+              topswad0_aero, solswad0_aero, &
+              topsw_aero, topsw0_aero, &
+              solsw_aero, solsw0_aero, &
+              topswcf_aero, solswcf_aero)
          
-c
-cIM 2eme calcul radiatif pour le cas perturbe ou au moins un
-cIM des taux doit etre different du taux actuel
-cIM Par defaut on a les taux perturbes egaux aux taux actuels
-c
+!
+!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
+!IM des taux doit etre different du taux actuel
+!IM Par defaut on a les taux perturbes egaux aux taux actuels
+!
       if (ok_4xCO2atm) then
-       if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR.
-     $RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR.
-     $RCFC12_per.NE.RCFC12_act) THEN
-c
+       if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
+      RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
+      RCFC12_per.NE.RCFC12_act) THEN
+!
        RCO2 = RCO2_per
        RCH4 = RCH4_per
@@ -3330 +3330 @@
        RCFC11 = RCFC11_per
        RCFC12 = RCFC12_per
-c
+!
       IF (prt_level .GE.10) THEN
        print *,' ->radlwsw, number 2 '
       ENDIF
-c
-         CALL radlwsw
-     e        (dist, rmu0, fract, 
-     e        paprs, pplay,zxtsol,albsol1, albsol2, 
-     e        t_seri,q_seri,wo,
-     e        cldfra, cldemi, cldtau,
-     e        ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol,
-     e        flag_aerosol_strat,
-     e        tau_aero, piz_aero, cg_aero,
-     e        cldtaupi,new_aod,
-     e        zqsat, flwc, fiwc,
-     s        heatp,heat0p,coolp,cool0p,radsolp,albplap,
-     s        topswp,toplwp,solswp,sollwp,
-     s        sollwdownp,
-     s        topsw0p,toplw0p,solsw0p,sollw0p,
-     s        lwdn0p, lwdnp, lwup0p, lwupp, 
-     s        swdn0p, swdnp, swup0p, swupp,
-     s        topswad_aerop, solswad_aerop,
-     s        topswai_aerop, solswai_aerop,
-     o        topswad0_aerop, solswad0_aerop,
-     o        topsw_aerop, topsw0_aerop,
-     o        solsw_aerop, solsw0_aerop,
-     o        topswcf_aerop, solswcf_aerop)
+!
+         CALL radlwsw &
+              (dist, rmu0, fract,  &
+              paprs, pplay,zxtsol,albsol1, albsol2,  &
+              t_seri,q_seri,wo, &
+              cldfra, cldemi, cldtau, &
+              ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
+              flag_aerosol_strat, &
+              tau_aero, piz_aero, cg_aero, &
+              cldtaupi,new_aod, &
+              zqsat, flwc, fiwc, &
+              heatp,heat0p,coolp,cool0p,radsolp,albplap, &
+              topswp,toplwp,solswp,sollwp, &
+              sollwdownp, &
+              topsw0p,toplw0p,solsw0p,sollw0p, &
+              lwdn0p, lwdnp, lwup0p, lwupp,  &
+              swdn0p, swdnp, swup0p, swupp, &
+              topswad_aerop, solswad_aerop, &
+              topswai_aerop, solswai_aerop, &
+              topswad0_aerop, solswad0_aerop, &
+              topsw_aerop, topsw0_aerop, &
+              solsw_aerop, solsw0_aerop, &
+              topswcf_aerop, solswcf_aerop)
        endif
       endif
-c
+!
       ENDIF ! aerosol_couple
       itaprad = 0
@@ -3387 +3387 @@
       END IF
 
-c
-c Ajouter la tendance des rayonnements (tous les pas)
-c
+!
+! Ajouter la tendance des rayonnements (tous les pas)
+!
       DO k = 1, klev
       DO i = 1, klon
-         t_seri(i,k) = t_seri(i,k)
-     .               + (heat(i,k)-cool(i,k)) * dtime/RDAY
-      ENDDO
-      ENDDO
-c
+         t_seri(i,k) = t_seri(i,k) &
+                     + (heat(i,k)-cool(i,k)) * dtime/RDAY
+      ENDDO
+      ENDDO
+!
       if (mydebug) then
         call writefield_phy('u_seri',u_seri,llm)
@@ -3404 +3404 @@
       endif
  
-cIM
+!IM
       IF (ip_ebil_phy.ge.2) THEN 
         ztit='after rad'
-        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
-     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
-     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-        call diagphy(airephy,ztit,ip_ebil_phy
-     e      , topsw, toplw, solsw, sollw, zero_v
-     e      , zero_v, zero_v, zero_v, ztsol
-     e      , d_h_vcol, d_qt, d_ec
-     s      , fs_bound, fq_bound )
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+            , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
+            , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+        call diagphy(airephy,ztit,ip_ebil_phy &
+            , topsw, toplw, solsw, sollw, zero_v &
+            , zero_v, zero_v, zero_v, ztsol &
+            , d_h_vcol, d_qt, d_ec &
+            , fs_bound, fq_bound )
       END IF 
-c
-c
-c Calculer l'hydrologie de la surface
-c
-c      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
-c     .            agesno, ftsol,fqsurf,fsnow, ruis)
-c
-
-c
-c Calculer le bilan du sol et la derive de temperature (couplage)
-c
+!
+!
+! Calculer l'hydrologie de la surface
+!
+!      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
+!     .            agesno, ftsol,fqsurf,fsnow, ruis)
+!
+
+!
+! Calculer le bilan du sol et la derive de temperature (couplage)
+!
       DO i = 1, klon
-c         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
-c a la demande de JLD
+!         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
+! a la demande de JLD
          bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
       ENDDO
-c
-cmoddeblott(jan95)
-c Appeler le programme de parametrisation de l'orographie
-c a l'echelle sous-maille:
-c
+!
+!moddeblott(jan95)
+! Appeler le programme de parametrisation de l'orographie
+! a l'echelle sous-maille:
+!
       IF (prt_level .GE.10) THEN
        print *,' call orography ? ', ok_orodr
       ENDIF
-c
+!
       IF (ok_orodr) THEN
-c
-c  selection des points pour lesquels le shema est actif:
+!
+!  selection des points pour lesquels le shema est actif:
         igwd=0
         DO i=1,klon
         itest(i)=0
-c        IF ((zstd(i).gt.10.0)) THEN
+!        IF ((zstd(i).gt.10.0)) THEN
         IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
           itest(i)=1
@@ -3454 +3454 @@
         ENDIF
         ENDDO
-c        igwdim=MAX(1,igwd)
-c
+!        igwdim=MAX(1,igwd)
+!
         IF (ok_strato) THEN
         
-          CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,
-     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
-     e                   igwd,idx,itest,
-     e                   t_seri, u_seri, v_seri,
-     s                   zulow, zvlow, zustrdr, zvstrdr,
-     s                   d_t_oro, d_u_oro, d_v_oro)
+          CALL drag_noro_strato(klon,klev,dtime,paprs,pplay, &
+                         zmea,zstd, zsig, zgam, zthe,zpic,zval, &
+                         igwd,idx,itest, &
+                         t_seri, u_seri, v_seri, &
+                         zulow, zvlow, zustrdr, zvstrdr, &
+                         d_t_oro, d_u_oro, d_v_oro)
 
        ELSE
-        CALL drag_noro(klon,klev,dtime,paprs,pplay,
-     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
-     e                   igwd,idx,itest,
-     e                   t_seri, u_seri, v_seri,
-     s                   zulow, zvlow, zustrdr, zvstrdr,
-     s                   d_t_oro, d_u_oro, d_v_oro)
+        CALL drag_noro(klon,klev,dtime,paprs,pplay, &
+                         zmea,zstd, zsig, zgam, zthe,zpic,zval, &
+                         igwd,idx,itest, &
+                         t_seri, u_seri, v_seri, &
+                         zulow, zvlow, zustrdr, zvstrdr, &
+                         d_t_oro, d_u_oro, d_v_oro)
        ENDIF
-c
-c  ajout des tendances
+!
+!  ajout des tendances
 !-----------------------------------------------------------------------------------------
 ! ajout des tendances de la trainee de l'orographie
       CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
 !-----------------------------------------------------------------------------------------
-c
+!
       ENDIF ! fin de test sur ok_orodr
-c
+!
       if (mydebug) then
         call writefield_phy('u_seri',u_seri,llm)
@@ -3490 +3490 @@
       
       IF (ok_orolf) THEN
-c
-c  selection des points pour lesquels le shema est actif:
+!
+!  selection des points pour lesquels le shema est actif:
         igwd=0
         DO i=1,klon
@@ -3501 +3501 @@
         ENDIF
         ENDDO
-c        igwdim=MAX(1,igwd)
-c
+!        igwdim=MAX(1,igwd)
+!
         IF (ok_strato) THEN
 
-          CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,
-     e                   rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,
-     e                   igwd,idx,itest,
-     e                   t_seri, u_seri, v_seri,
-     s                   zulow, zvlow, zustrli, zvstrli,
-     s                   d_t_lif, d_u_lif, d_v_lif               )
+          CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
+                         rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
+                         igwd,idx,itest, &
+                         t_seri, u_seri, v_seri, &
+                         zulow, zvlow, zustrli, zvstrli, &
+                         d_t_lif, d_u_lif, d_v_lif               )
         
         ELSE
-          CALL lift_noro(klon,klev,dtime,paprs,pplay,
-     e                   rlat,zmea,zstd,zpic,
-     e                   itest,
-     e                   t_seri, u_seri, v_seri,
-     s                   zulow, zvlow, zustrli, zvstrli,
-     s                   d_t_lif, d_u_lif, d_v_lif)
+          CALL lift_noro(klon,klev,dtime,paprs,pplay, &
+                         rlat,zmea,zstd,zpic, &
+                         itest, &
+                         t_seri, u_seri, v_seri, &
+                         zulow, zvlow, zustrli, zvstrli, &
+                         d_t_lif, d_u_lif, d_v_lif)
        ENDIF
-c   
+!   
 !-----------------------------------------------------------------------------------------
 ! ajout des tendances de la portance de l'orographie
       CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
 !-----------------------------------------------------------------------------------------
-c
+!
       ENDIF ! fin de test sur ok_orolf
-C  HINES GWD PARAMETRIZATION
+!  HINES GWD PARAMETRIZATION
 
        IF (ok_hines) then
 
-         CALL hines_gwd(klon,klev,dtime,paprs,pplay,
-     i                  rlat,t_seri,u_seri,v_seri,
-     o                  zustrhi,zvstrhi,
-     o                  d_t_hin, d_u_hin, d_v_hin)
-c
-c  ajout des tendances
+         CALL hines_gwd(klon,klev,dtime,paprs,pplay, &
+                        rlat,t_seri,u_seri,v_seri, &
+                        zustrhi,zvstrhi, &
+                        d_t_hin, d_u_hin, d_v_hin)
+!
+!  ajout des tendances
         CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'hin')
 
       ENDIF
-c
-
-c
-cIM cf. FLott BEG
-C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
+!
+
+!
+!IM cf. FLott BEG
+! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
 
       if (mydebug) then
@@ -3559 +3559 @@
       DO k = 1, klev
       DO i = 1, klon
-       zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
-     c            (paprs(i,k)-paprs(i,k+1))/rg
-       zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
-     c            (paprs(i,k)-paprs(i,k+1))/rg
-      ENDDO
-      ENDDO
-c
-cIM calcul composantes axiales du moment angulaire et couple des montagnes
-c
+       zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
+                  (paprs(i,k)-paprs(i,k+1))/rg
+       zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
+                  (paprs(i,k)-paprs(i,k+1))/rg
+      ENDDO
+      ENDDO
+!
+!IM calcul composantes axiales du moment angulaire et couple des montagnes
+!
       IF (is_sequential .and. ok_orodr) THEN 
-        CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
-     C                 ra,rg,romega,
-     C                 rlat,rlon,pphis,
-     C                 zustrdr,zustrli,zustrph,
-     C                 zvstrdr,zvstrli,zvstrph,
-     C                 paprs,u,v,
-     C                 aam, torsfc)
+        CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
+                       ra,rg,romega, &
+                       rlat,rlon,pphis, &
+                       zustrdr,zustrli,zustrph, &
+                       zvstrdr,zvstrli,zvstrph, &
+                       paprs,u,v, &
+                       aam, torsfc)
        ENDIF
-cIM cf. FLott END
-cIM
+!IM cf. FLott END
+!IM
       IF (ip_ebil_phy.ge.2) THEN 
         ztit='after orography'
-        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
-     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
-     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-         call diagphy(airephy,ztit,ip_ebil_phy
-     e      , zero_v, zero_v, zero_v, zero_v, zero_v
-     e      , zero_v, zero_v, zero_v, ztsol
-     e      , d_h_vcol, d_qt, d_ec
-     s      , fs_bound, fq_bound )
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
+            , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
+            , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+         call diagphy(airephy,ztit,ip_ebil_phy &
+            , zero_v, zero_v, zero_v, zero_v, zero_v &
+            , zero_v, zero_v, zero_v, ztsol &
+            , d_h_vcol, d_qt, d_ec &
+            , fs_bound, fq_bound )
       END IF 
-c
-c
+!
+!
 !====================================================================
 ! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
@@ -3606 +3606 @@
 !       print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
 !     s        ref_liq,ref_ice
-          call phys_cosp(itap,dtime,freq_cosp,
-     $                   ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP,
-     $                   ecrit_mth,ecrit_day,ecrit_hf,
-     $                   klon,klev,rlon,rlat,presnivs,overlap,
-     $                   ref_liq,ref_ice,
-     $                   pctsrf(:,is_ter)+pctsrf(:,is_lic),
-     $                   zu10m,zv10m,pphis,
-     $                   zphi,paprs(:,1:klev),pplay,zxtsol,t_seri,
-     $                   qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc,
-     $                   prfl(:,1:klev),psfl(:,1:klev),
-     $                   pmflxr(:,1:klev),pmflxs(:,1:klev),
-     $                   mr_ozone,cldtau, cldemi)
+          call phys_cosp(itap,dtime,freq_cosp, &
+                         ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
+                         ecrit_mth,ecrit_day,ecrit_hf, &
+                         klon,klev,rlon,rlat,presnivs,overlap, &
+                         ref_liq,ref_ice, &
+                         pctsrf(:,is_ter)+pctsrf(:,is_lic), &
+                         zu10m,zv10m,pphis, &
+                         zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
+                         qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
+                         prfl(:,1:klev),psfl(:,1:klev), &
+                         pmflxr(:,1:klev),pmflxs(:,1:klev), &
+                         mr_ozone,cldtau, cldemi)
 
 !     L          calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
@@ -3630 +3630 @@
        ENDIF  !ok_cosp
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-cAA
-cAA Installation de l'interface online-offline pour traceurs
-cAA
-c====================================================================
-c   Calcul  des tendances traceurs
-c====================================================================
-C
+!AA
+!AA Installation de l'interface online-offline pour traceurs
+!AA
+!====================================================================
+!   Calcul  des tendances traceurs
+!====================================================================
+!
 
        IF (type_trac=='repr') THEN
@@ -3644 +3644 @@
        END IF
 
-      call phytrac (
-     I     itap,     days_elapsed+1,    jH_cur,   debut,
-     I     lafin,    dtime,     u, v,     t,
-     I     paprs,    pplay,     pmfu,     pmfd,
-     I     pen_u,    pde_u,     pen_d,    pde_d,
-     I     cdragh,   coefh(:,:,is_ave),   fm_therm, entr_therm,
-     I     u1,       v1,        ftsol,    pctsrf,
-     I     zustar,   zu10m,     zv10m,
-     I     wstar(:,is_ave),    ale_bl,         ale_wake,
-     I     rlat,     rlon,
-     I     frac_impa,frac_nucl, beta_prec_fisrt,beta_prec,
-     I     presnivs, pphis,     pphi,     albsol1,
-     I     sh_in,    rhcl,      cldfra,   rneb,
-     I     diafra,   cldliq,    itop_con, ibas_con,
-     I     pmflxr,   pmflxs,    prfl,     psfl,
-     I     da,       phi,       mp,       upwd,
-     I     phi2,     d1a,       dam,      sij,         !<<RomP
-     I     wdtrainA, wdtrainM,  sigd,     clw,elij,    !<<RomP
-     I     ev,       ep,        epmlmMm,  eplaMm,      !<<RomP
-     I     dnwd,     aerosol_couple,      flxmass_w,
-     I     tau_aero, piz_aero,  cg_aero,  ccm,
-     I     rfname,
-     I     d_tr_dyn,                                   !<<RomP
-     O     tr_seri)
+      call phytrac ( &
+           itap,     days_elapsed+1,    jH_cur,   debut, &
+           lafin,    dtime,     u, v,     t, &
+           paprs,    pplay,     pmfu,     pmfd, &
+           pen_u,    pde_u,     pen_d,    pde_d, &
+           cdragh,   coefh(:,:,is_ave),   fm_therm, entr_therm, &
+           u1,       v1,        ftsol,    pctsrf, &
+           zustar,   zu10m,     zv10m, &
+           wstar(:,is_ave),    ale_bl,         ale_wake, &
+           rlat,     rlon, &
+           frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
+           presnivs, pphis,     pphi,     albsol1, &
+           sh_in,    rhcl,      cldfra,   rneb, &
+           diafra,   cldliq,    itop_con, ibas_con, &
+           pmflxr,   pmflxs,    prfl,     psfl, &
+           da,       phi,       mp,       upwd, &
+           phi2,     d1a,       dam,      sij, &        !<<RomP
+           wdtrainA, wdtrainM,  sigd,     clw,elij, &   !<<RomP
+           ev,       ep,        epmlmMm,  eplaMm, &     !<<RomP
+           dnwd,     aerosol_couple,      flxmass_w, &
+           tau_aero, piz_aero,  cg_aero,  ccm, &
+           rfname, &
+           d_tr_dyn, &                                 !<<RomP
+           tr_seri)
 
       IF (offline) THEN
 
-       IF (prt_level.ge.9)
-     $    print*,'Attention on met a 0 les thermiques pour phystoke'
-         call phystokenc (
-     I                   nlon,klev,pdtphys,rlon,rlat,
-     I                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
-     I                   fm_therm,entr_therm,
-     I                   cdragh,coefh(:,:,is_ave),u1,v1,ftsol,pctsrf,
-     I                   frac_impa, frac_nucl,
-     I                   pphis,airephy,dtime,itap,
-     I                   qx(:,:,ivap),da,phi,mp,upwd,dnwd)
+       IF (prt_level.ge.9) &
+          print*,'Attention on met a 0 les thermiques pour phystoke'
+         call phystokenc ( &
+                         nlon,klev,pdtphys,rlon,rlat, &
+                         t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
+                         fm_therm,entr_therm, &
+                         cdragh,coefh(:,:,is_ave),u1,v1,ftsol,pctsrf, &
+                         frac_impa, frac_nucl, &
+                         pphis,airephy,dtime,itap, &
+                         qx(:,:,ivap),da,phi,mp,upwd,dnwd)
 
 
       ENDIF
 
-c
-c Calculer le transport de l'eau et de l'energie (diagnostique)
-c
-      CALL transp (paprs,zxtsol,
-     e                   t_seri, q_seri, u_seri, v_seri, zphi,
-     s                   ve, vq, ue, uq)
-c
-cIM global posePB BEG
+!
+! Calculer le transport de l'eau et de l'energie (diagnostique)
+!
+      CALL transp (paprs,zxtsol, &
+                         t_seri, q_seri, u_seri, v_seri, zphi, &
+                         ve, vq, ue, uq)
+!
+!IM global posePB BEG
       IF(1.EQ.0) THEN
-c
-      CALL transp_lay (paprs,zxtsol,
-     e                   t_seri, q_seri, u_seri, v_seri, zphi,
-     s                   ve_lay, vq_lay, ue_lay, uq_lay)
-c
+!
+      CALL transp_lay (paprs,zxtsol, &
+                         t_seri, q_seri, u_seri, v_seri, zphi, &
+                         ve_lay, vq_lay, ue_lay, uq_lay)
+!
       ENDIF !(1.EQ.0) THEN
-cIM global posePB END
-c Accumuler les variables a stocker dans les fichiers histoire:
-c
+!IM global posePB END
+! Accumuler les variables a stocker dans les fichiers histoire:
+!
 
 !================================================================
@@ -3711 +3711 @@
       d_t_ec(:,:)=0.
       forall (k=1: llm) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
-      CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),
-     s        u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:),
-     s        zmasse,exner,d_t_ec)
+      CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap), &
+              u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
+              zmasse,exner,d_t_ec)
       t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
 
-cIM
+!IM
       IF (ip_ebil_phy.ge.1) THEN 
         ztit='after physic'
-        CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
-     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
-     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-C     Comme les tendances de la physique sont ajoute dans la dynamique,
-C     on devrait avoir que la variation d'entalpie par la dynamique
-C     est egale a la variation de la physique au pas de temps precedent.
-C     Donc la somme de ces 2 variations devrait etre nulle.
-
-        call diagphy(airephy,ztit,ip_ebil_phy
-     e      , topsw, toplw, solsw, sollw, sens
-     e      , evap, rain_fall, snow_fall, ztsol
-     e      , d_h_vcol, d_qt, d_ec
-     s      , fs_bound, fq_bound )
-C
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
+            , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
+            , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+!     Comme les tendances de la physique sont ajoute dans la dynamique,
+!     on devrait avoir que la variation d'entalpie par la dynamique
+!     est egale a la variation de la physique au pas de temps precedent.
+!     Donc la somme de ces 2 variations devrait etre nulle.
+
+        call diagphy(airephy,ztit,ip_ebil_phy &
+            , topsw, toplw, solsw, sollw, sens &
+            , evap, rain_fall, snow_fall, ztsol &
+            , d_h_vcol, d_qt, d_ec &
+            , fs_bound, fq_bound )
+!
       d_h_vcol_phy=d_h_vcol
-C
+!
       END IF 
-C
-c=======================================================================
-c   SORTIES
-c=======================================================================
-
-cIM Interpolation sur les niveaux de pression du NMC
-c   -------------------------------------------------
-c
+!
+!=======================================================================
+!   SORTIES
+!=======================================================================
+
+!IM Interpolation sur les niveaux de pression du NMC
+!   -------------------------------------------------
+!
 #include "calcul_STDlev.h"
-c
-c slp sea level pressure
+!
+! slp sea level pressure
       slp(:) = paprs(:,1)*exp(pphis(:)/(RD*t_seri(:,1)))
-c
-ccc prw = eau precipitable
+!
+!cc prw = eau precipitable
       DO i = 1, klon
        prw(i) = 0.
        DO k = 1, klev
-        prw(i) = prw(i) +
-     .           q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
+        prw(i) = prw(i) + &
+                 q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
        ENDDO
       ENDDO
-c
-cIM initialisation + calculs divers diag AMIP2
-c
+!
+!IM initialisation + calculs divers diag AMIP2
+!
 #include "calcul_divers.h"
-c
+!
       IF (type_trac == 'inca') THEN
 #ifdef INCA
@@ -3767 +3767 @@
          CALL VTb(VTinca)
 
-         CALL chemhook_end (
-     $                        dtime,
-     $                        pplay,
-     $                        t_seri,
-     $                        tr_seri,
-     $                        nbtr,
-     $                        paprs,
-     $                        q_seri,
-     $                        airephy,
-     $                        pphi,
-     $                        pphis,
-     $                        zx_rh)
+         CALL chemhook_end ( &
+                              dtime, &
+                              pplay, &
+                              t_seri, &
+                              tr_seri, &
+                              nbtr, &
+                              paprs, &
+                              q_seri, &
+                              airephy, &
+                              pphi, &
+                              pphis, &
+                              zx_rh)
 
          CALL VTe(VTinca)
@@ -3786 +3786 @@
 
 
-c
-c Convertir les incrementations en tendances
-c
+!
+! Convertir les incrementations en tendances
+!
       IF (prt_level .GE.10) THEN
         print *,'Convertir les incrementations en tendances '
       ENDIF
-c
+!
       if (mydebug) then
         call writefield_phy('u_seri',u_seri,llm)
@@ -3809 +3809 @@
       ENDDO
       ENDDO
-c
+!
       IF (nqtot.GE.3) THEN
       DO iq = 3, nqtot
@@ -3819 +3819 @@
       ENDDO
       ENDIF
-c
-cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
-cIM global posePB#include "write_bilKP_ins.h"
-cIM global posePB#include "write_bilKP_ave.h"
-c
-
-c Sauvegarder les valeurs de t et q a la fin de la physique:
-c
+!
+!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
+!IM global posePB#include "write_bilKP_ins.h"
+!IM global posePB#include "write_bilKP_ave.h"
+!
+
+! Sauvegarder les valeurs de t et q a la fin de la physique:
+!
       DO k = 1, klev
       DO i = 1, klon
@@ -3856 +3856 @@
       if (prt_level.ge.1) then
       write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
-      write(lunout,*)
-     s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
-      write(lunout,*)
-     s  nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
-     s  pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
-     s  pctsrf(igout,is_sic)
+      write(lunout,*) &
+       'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
+      write(lunout,*) &
+        nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
+        pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
+        pctsrf(igout,is_sic)
       write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
       do k=1,klev
-         write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
-     s   d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
-     s   d_t_eva(igout,k)
+         write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
+         d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
+         d_t_eva(igout,k)
       enddo
       write(lunout,*) 'cool,heat'
@@ -3875 +3875 @@
       write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
       do k=1,klev
-         write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
-     s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
+         write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
+       d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
       enddo
 
@@ -3882 +3882 @@
       write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
       do k=1,klev
-         write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
-     s  d_qx(igout,k,1),d_qx(igout,k,2)
+         write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
+        d_qx(igout,k,1),d_qx(igout,k,2)
       enddo
       endif
@@ -3889 +3889 @@
 !==========================================================================
 
-c============================================================
-c   Calcul de la temperature potentielle
-c============================================================
+!============================================================
+!   Calcul de la temperature potentielle
+!============================================================
       DO k = 1, klev
       DO i = 1, klon
-cJYG/IM theta en debut du pas de temps
-cJYG/IM       theta(i,k)=t(i,k)*(100000./pplay(i,k))**(RD/RCPD)
-cJYG/IM theta en fin de pas de temps de physique
+!JYG/IM theta en debut du pas de temps
+!JYG/IM       theta(i,k)=t(i,k)*(100000./pplay(i,k))**(RD/RCPD)
+!JYG/IM theta en fin de pas de temps de physique
         theta(i,k)=t_seri(i,k)*(100000./pplay(i,k))**(RD/RCPD)
-c thetal: 2 lignes suivantes a decommenter si vous avez les fichiers     MPL 20130625
-c fth_fonctions.F90 et parkind1.F90
-c sinon thetal=theta
-c       thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
-c    :         ql_seri(i,k))
+! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers     MPL 20130625
+! fth_fonctions.F90 et parkind1.F90
+! sinon thetal=theta
+!       thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
+!    :         ql_seri(i,k))
         thetal(i,k)=theta(i,k)
       ENDDO
       ENDDO
-c
-
-c 22.03.04 BEG
-c=============================================================
-c   Ecriture des sorties
-c=============================================================
+!
+
+! 22.03.04 BEG
+!=============================================================
+!   Ecriture des sorties
+!=============================================================
 #ifdef CPP_IOIPSL
  
-c Recupere des varibles calcule dans differents modules
-c pour ecriture dans histxxx.nc 
+! Recupere des varibles calcule dans differents modules
+! pour ecriture dans histxxx.nc 
 
       ! Get some variables from module fonte_neige_mod
-      CALL fonte_neige_get_vars(pctsrf, 
-     .     zxfqcalving, zxfqfonte, zxffonte)
-
-
-
-
-c=============================================================
+      CALL fonte_neige_get_vars(pctsrf,  &
+           zxfqcalving, zxfqfonte, zxffonte)
+
+
+
+
+!=============================================================
 ! Separation entre thermiques et non thermiques dans les sorties
 ! de fisrtilp
-c=============================================================
+!=============================================================
 
       if (iflag_thermals>=1) then
@@ -3956 +3956 @@
 
       CALL phys_output_write(itap, pdtphys, paprs, pphis,               &
-     &                  pplay, lmax_th, aerosol_couple,                 &
-     &                  ok_ade, ok_aie, ivap, new_aod, ok_sync,         &
-     &                  ptconv, read_climoz, clevSTD, freq_moyNMC,      &
-     &                  ptconvth, d_t, qx, d_qx, zmasse,                &
-     &                  flag_aerosol_strat)
+                        pplay, lmax_th, aerosol_couple,                 &
+                        ok_ade, ok_aie, ivap, new_aod, ok_sync,         &
+                        ptconv, read_climoz, clevSTD, freq_moyNMC,      &
+                        ptconvth, d_t, qx, d_qx, zmasse,                &
+                        flag_aerosol_strat)
 
 
@@ -3975 +3975 @@
 #endif
 
-c 22.03.04 END
-c
-c====================================================================
-c Si c'est la fin, il faut conserver l'etat de redemarrage
-c====================================================================
-c
-
-c        -----------------------------------------------------------------
-c        WSTATS: Saving statistics
-c        -----------------------------------------------------------------
-c        ("stats" stores and accumulates 8 key variables in file "stats.nc"
-c        which can later be used to make the statistic files of the run:
-c        "stats")          only possible in 3D runs !
+! 22.03.04 END
+!
+!====================================================================
+! Si c'est la fin, il faut conserver l'etat de redemarrage
+!====================================================================
+!
+
+!        -----------------------------------------------------------------
+!        WSTATS: Saving statistics
+!        -----------------------------------------------------------------
+!        ("stats" stores and accumulates 8 key variables in file "stats.nc"
+!        which can later be used to make the statistic files of the run:
+!        "stats")          only possible in 3D runs !
 
          
          IF (callstats) THEN
 
-           call wstats(klon,o_psol%name,"Surface pressure","Pa"
-     &                 ,2,paprs(:,1))
-           call wstats(klon,o_tsol%name,"Surface temperature","K",
-     &                 2,zxtsol)
+           call wstats(klon,o_psol%name,"Surface pressure","Pa" &
+                       ,2,paprs(:,1))
+           call wstats(klon,o_tsol%name,"Surface temperature","K", &
+                       2,zxtsol)
            zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
-           call wstats(klon,o_precip%name,"Precip Totale liq+sol",
-     &                 "kg/(s*m2)",2,zx_tmp_fi2d)
+           call wstats(klon,o_precip%name,"Precip Totale liq+sol", &
+                       "kg/(s*m2)",2,zx_tmp_fi2d)
            zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
-           call wstats(klon,o_plul%name,"Large-scale Precip",
-     &                 "kg/(s*m2)",2,zx_tmp_fi2d)
+           call wstats(klon,o_plul%name,"Large-scale Precip", &
+                       "kg/(s*m2)",2,zx_tmp_fi2d)
            zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
-           call wstats(klon,o_pluc%name,"Convective Precip",
-     &                 "kg/(s*m2)",2,zx_tmp_fi2d)
-           call wstats(klon,o_sols%name,"Solar rad. at surf.",
-     &                 "W/m2",2,solsw)
-           call wstats(klon,o_soll%name,"IR rad. at surf.",
-     &                 "W/m2",2,sollw)
+           call wstats(klon,o_pluc%name,"Convective Precip", &
+                       "kg/(s*m2)",2,zx_tmp_fi2d)
+           call wstats(klon,o_sols%name,"Solar rad. at surf.", &
+                       "W/m2",2,solsw)
+           call wstats(klon,o_soll%name,"IR rad. at surf.", &
+                       "W/m2",2,sollw)
           zx_tmp_fi2d(:) = topsw(:)-toplw(:)
-          call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA",
-     &                 "W/m2",2,zx_tmp_fi2d)
-
-
-
-           call wstats(klon,o_temp%name,"Air temperature","K",
-     &                 3,t_seri)
-           call wstats(klon,o_vitu%name,"Zonal wind","m.s-1",
-     &                 3,u_seri)
-           call wstats(klon,o_vitv%name,"Meridional wind",
-     &                "m.s-1",3,v_seri)
-           call wstats(klon,o_vitw%name,"Vertical wind",
-     &                "m.s-1",3,omega)
-           call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg",
-     &                 3,q_seri)
+          call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA", &
+                       "W/m2",2,zx_tmp_fi2d)
+
+
+
+           call wstats(klon,o_temp%name,"Air temperature","K", &
+                       3,t_seri)
+           call wstats(klon,o_vitu%name,"Zonal wind","m.s-1", &
+                       3,u_seri)
+           call wstats(klon,o_vitv%name,"Meridional wind", &
+                      "m.s-1",3,v_seri)
+           call wstats(klon,o_vitw%name,"Vertical wind", &
+                      "m.s-1",3,omega)
+           call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg", &
+                       3,q_seri)
  
 
@@ -4041 +4041 @@
 !         write(97) u_seri,v_seri,t_seri,q_seri
 !         close(97)
-C$OMP MASTER
+!$OMP MASTER
          if (read_climoz >= 1) then
             if (is_mpi_root) then
@@ -4048 +4048 @@
             deallocate(press_climoz) ! pointer
          end if
-C$OMP END MASTER
+!$OMP END MASTER
       ENDIF
       
@@ -4057 +4057 @@
       FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
       IMPLICIT none
-c
-c Calculer et imprimer l'eau totale. A utiliser pour verifier
-c la conservation de l'eau
-c
+!
+! Calculer et imprimer l'eau totale. A utiliser pour verifier
+! la conservation de l'eau
+!
 #include "YOMCST.h"
       INTEGER klon,klev
@@ -4067 +4067 @@
       REAL qtotal, zx, qcheck
       INTEGER i, k
-c
+!
       zx = 0.0
       DO i = 1, klon
@@ -4075 +4075 @@
       DO k = 1, klev
       DO i = 1, klon
-         qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
-     .                     *(paprs(i,k)-paprs(i,k+1))/RG
-      ENDDO
-      ENDDO
-c
+         qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i) &
+                           *(paprs(i,k)-paprs(i,k+1))/RG
+      ENDDO
+      ENDDO
+!
       qcheck = qtotal/zx
-c
+!
       RETURN
       END
       SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
       IMPLICIT none
-c
-c Tranformer une variable de la grille physique a
-c la grille d'ecriture
-c
+!
+! Tranformer une variable de la grille physique a
+! la grille d'ecriture
+!
       INTEGER nfield,nlon,iim,jjmp1, jjm
       REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
-c
+!
       INTEGER i, n, ig
-c
+!
       jjm = jjmp1 - 1
       DO n = 1, nfield