Changeset 2019

Timestamp:

Apr 16, 2014, 7:16:58 AM (10 years ago)

Author:

fhourdin

Message:

Passage au format libre pour inclure les ficheirs du 1D dans lmdz1d.F90

Location:

LMDZ5/trunk/libf

Files:

• dyn3d/conf_gcm.F (modified) (32 diffs)
• dyn3d_common/q_sat.F (modified) (3 diffs)
• phylmd/1DUTILS.h (modified) (88 diffs)
• phylmd/1D_decl_cases.h (modified) (1 diff)
• phylmd/1D_interp_cases.h (modified) (17 diffs)
• phylmd/1D_nudge_sandu_astex.h (modified) (1 diff)
• phylmd/1D_read_forc_cases.h (modified) (22 diffs)
• phylmd/1Dconv.h (modified) (39 diffs)
• phylmd/compar1d.h (modified) (1 diff)
• phylmd/fcg_gcssold.h (modified) (1 diff)
• phylmd/lmdz1d.F90 (modified) (1 diff)

LMDZ5/trunk/libf/dyn3d/conf_gcm.F

@@ -2 +2 @@
 ! $Id$
 !
-c
-c
+!
+!
       SUBROUTINE conf_gcm( tapedef, etatinit, clesphy0 )
-c
+!
       USE control_mod
 #ifdef CPP_IOIPSL
@@ -17 +17 @@
 
       IMPLICIT NONE
-c-----------------------------------------------------------------------
-c     Auteurs :   L. Fairhead , P. Le Van  .
-c
-c     Arguments :
-c
-c     tapedef   :
-c     etatinit  :     = TRUE   , on ne  compare pas les valeurs des para- 
-c     -metres  du zoom  avec  celles lues sur le fichier start .
-c      clesphy0 :  sortie  .
-c
+!-----------------------------------------------------------------------
+!     Auteurs :   L. Fairhead , P. Le Van  .
+!
+!     Arguments :
+!
+!     tapedef   :
+!     etatinit  :     = TRUE   , on ne  compare pas les valeurs des para- 
+!     -metres  du zoom  avec  celles lues sur le fichier start .
+!      clesphy0 :  sortie  .
+!
        LOGICAL etatinit
        INTEGER tapedef
@@ -33 +33 @@
        PARAMETER(     longcles = 20 )
        REAL clesphy0( longcles )
-c
-c   Declarations :
-c   --------------
+!
+!   Declarations :
+!   --------------
 #include "dimensions.h"
 #include "paramet.h"
@@ -47 +47 @@
 ! #include "clesphys.h"
 #include "iniprint.h"
-c
-c
-c   local:
-c   ------
+!
+!
+!   local:
+!   ------
 
       CHARACTER ch1*72,ch2*72,ch3*72,ch4*12
@@ -58 +58 @@
       INTEGER i
       LOGICAL use_filtre_fft
-c
-c  -------------------------------------------------------------------
-c
-c       .........     Version  du 29/04/97       ..........
-c
-c   Nouveaux parametres nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,
-c      tetatemp   ajoutes  pour la dissipation   .
-c
-c   Autre parametre ajoute en fin de liste de tapedef : ** fxyhypb ** 
-c
-c  Si fxyhypb = .TRUE. , choix de la fonction a derivee tangente hyperb.
-c    Sinon , choix de fxynew  , a derivee sinusoidale  ..
-c
-c   ......  etatinit = . TRUE. si defrun  est appele dans ETAT0_LMD  ou
-c         LIMIT_LMD  pour l'initialisation de start.dat (dic) et
-c                de limit.dat ( dic)                        ...........
-c           Sinon  etatinit = . FALSE .
-c
-c   Donc etatinit = .F.  si on veut comparer les valeurs de  grossismx ,
-c    grossismy,clon,clat, fxyhypb  lues sur  le fichier  start  avec
-c   celles passees  par run.def ,  au debut du gcm, apres l'appel a 
-c    lectba .  
-c   Ces parmetres definissant entre autres la grille et doivent etre
-c   pareils et coherents , sinon il y aura  divergence du gcm .
-c
-c-----------------------------------------------------------------------
-c   initialisations:
-c   ----------------
+!
+!  -------------------------------------------------------------------
+!
+!       .........     Version  du 29/04/97       ..........
+!
+!   Nouveaux parametres nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,
+!      tetatemp   ajoutes  pour la dissipation   .
+!
+!   Autre parametre ajoute en fin de liste de tapedef : ** fxyhypb ** 
+!
+!  Si fxyhypb = .TRUE. , choix de la fonction a derivee tangente hyperb.
+!    Sinon , choix de fxynew  , a derivee sinusoidale  ..
+!
+!   ......  etatinit = . TRUE. si defrun  est appele dans ETAT0_LMD  ou
+!         LIMIT_LMD  pour l'initialisation de start.dat (dic) et
+!                de limit.dat ( dic)                        ...........
+!           Sinon  etatinit = . FALSE .
+!
+!   Donc etatinit = .F.  si on veut comparer les valeurs de  grossismx ,
+!    grossismy,clon,clat, fxyhypb  lues sur  le fichier  start  avec
+!   celles passees  par run.def ,  au debut du gcm, apres l'appel a 
+!    lectba .  
+!   Ces parmetres definissant entre autres la grille et doivent etre
+!   pareils et coherents , sinon il y aura  divergence du gcm .
+!
+!-----------------------------------------------------------------------
+!   initialisations:
+!   ----------------
 
 !Config  Key  = lunout
@@ -95 +95 @@
       CALL getin('lunout', lunout)
       IF (lunout /= 5 .and. lunout /= 6) THEN
-        OPEN(UNIT=lunout,FILE='lmdz.out',ACTION='write', 
+        OPEN(UNIT=lunout,FILE='lmdz.out',ACTION='write',                     &
      &          STATUS='unknown',FORM='formatted')
       ENDIF
@@ -107 +107 @@
       CALL getin('prt_level',prt_level)
 
-c-----------------------------------------------------------------------
-c  Parametres de controle du run:
-c-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
+!  Parametres de controle du run:
+!-----------------------------------------------------------------------
 !Config  Key  = planet_type
 !Config  Desc = planet type ("earth", "mars", "venus", ...)
@@ -232 +232 @@
        CALL getin('dissip_period',dissip_period)
 
-ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
-ccc
+!cc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
+!cc
 
 !Config  Key  = lstardis
@@ -348 +348 @@
        CALL getin('ok_guide',ok_guide)
 
-c    ...............................................................
+!    ...............................................................
 
 !Config  Key  =  read_start
@@ -390 +390 @@
       ENDDO
 
-ccc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
-c     .........   (  modif  le 17/04/96 )   .........
-c
+!cc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
+!     .........   (  modif  le 17/04/96 )   .........
+!
       IF( etatinit ) GO TO 100
 
@@ -411 +411 @@
        CALL getin('clat',clatt)
 
-c
-c
+!
+!
       IF( ABS(clat - clatt).GE. 0.001 )  THEN
-        write(lunout,*)'conf_gcm: La valeur de clat passee par run.def',
+        write(lunout,*)'conf_gcm: La valeur de clat passee par run.def',     &
      &    ' est differente de celle lue sur le fichier  start '
         STOP
@@ -429 +429 @@
 
       IF( ABS(grossismx - grossismxx).GE. 0.001 )  THEN
-        write(lunout,*)'conf_gcm: La valeur de grossismx passee par ',
+        write(lunout,*)'conf_gcm: La valeur de grossismx passee par ',       &
      &  'run.def est differente de celle lue sur le fichier  start '
         STOP
@@ -443 +443 @@
 
       IF( ABS(grossismy - grossismyy).GE. 0.001 )  THEN
-        write(lunout,*)'conf_gcm: La valeur de grossismy passee par ',
+        write(lunout,*)'conf_gcm: La valeur de grossismy passee par ',        &
      & 'run.def est differente de celle lue sur le fichier  start '
         STOP
@@ -449 +449 @@
       
       IF( grossismx.LT.1. )  THEN
-        write(lunout,*)
+        write(lunout,*)                                                        &
      &       'conf_gcm: ***  ATTENTION !! grossismx < 1 .   *** '
          STOP
@@ -458 +458 @@
 
       IF( grossismy.LT.1. )  THEN
-        write(lunout,*)
+        write(lunout,*)                                                        &
      &       'conf_gcm: ***  ATTENTION !! grossismy < 1 .   *** '
          STOP
@@ -466 +466 @@
 
       write(lunout,*)'conf_gcm: alphax alphay',alphax,alphay
-c
-c    alphax et alphay sont les anciennes formulat. des grossissements
-c
-c
+!
+!    alphax et alphay sont les anciennes formulat. des grossissements
+!
+!
 
 !Config  Key  = fxyhypb
@@ -482 +482 @@
          IF( fxyhypbb )     THEN
             write(lunout,*)' ********  PBS DANS  CONF_GCM  ******** '
-            write(lunout,*)' *** fxyhypb lu sur le fichier start est ',
-     *       'F alors  qu il est  T  sur  run.def  ***'
+            write(lunout,*)' *** fxyhypb lu sur le fichier start est ',     &
+     &       'F alors  qu il est  T  sur  run.def  ***'
               STOP
          ENDIF
@@ -489 +489 @@
          IF( .NOT.fxyhypbb )   THEN
             write(lunout,*)' ********  PBS DANS  CONF_GCM  ******** '
-            write(lunout,*)' ***  fxyhypb lu sur le fichier start est ',
-     *        'T alors  qu il est  F  sur  run.def  ****  '
+            write(lunout,*)' ***  fxyhypb lu sur le fichier start est ',     &
+     &        'T alors  qu il est  F  sur  run.def  ****  '
               STOP
          ENDIF
       ENDIF
-c
+!
 !Config  Key  = dzoomx
 !Config  Desc = extension en longitude
@@ -505 +505 @@
       IF( fxyhypb )  THEN
        IF( ABS(dzoomx - dzoomxx).GE. 0.001 )  THEN
-        write(lunout,*)'conf_gcm: La valeur de dzoomx passee par ',
-     *  'run.def est differente de celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de dzoomx passee par ',         &
+     &  'run.def est differente de celle lue sur le fichier  start '
         STOP
        ENDIF
@@ -521 +521 @@
       IF( fxyhypb )  THEN
        IF( ABS(dzoomy - dzoomyy).GE. 0.001 )  THEN
-        write(lunout,*)'conf_gcm: La valeur de dzoomy passee par ',
-     * 'run.def est differente de celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de dzoomy passee par ',          &
+     & 'run.def est differente de celle lue sur le fichier  start '
         STOP
        ENDIF
@@ -536 +536 @@
       IF( fxyhypb )  THEN
        IF( ABS(taux - tauxx).GE. 0.001 )  THEN
-        write(lunout,*)'conf_gcm: La valeur de taux passee par ',
-     * 'run.def est differente de celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de taux passee par ',           &
+     & 'run.def est differente de celle lue sur le fichier  start '
         STOP
        ENDIF
@@ -551 +551 @@
       IF( fxyhypb )  THEN
        IF( ABS(tauy - tauyy).GE. 0.001 )  THEN
-        write(lunout,*)'conf_gcm: La valeur de tauy passee par ',
-     * 'run.def est differente de celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de tauy passee par ',           &
+     & 'run.def est differente de celle lue sur le fichier  start '
         STOP
        ENDIF
       ENDIF
 
-cc
+!c
       IF( .NOT.fxyhypb  )  THEN
 
@@ -572 +572 @@
           IF( ysinuss )     THEN
             write(lunout,*)' ********  PBS DANS  CONF_GCM  ******** '
-            write(lunout,*)' *** ysinus lu sur le fichier start est F',
-     *       ' alors  qu il est  T  sur  run.def  ***'
+            write(lunout,*)' *** ysinus lu sur le fichier start est F',     &
+     &       ' alors  qu il est  T  sur  run.def  ***'
             STOP
           ENDIF
@@ -579 +579 @@
           IF( .NOT.ysinuss )   THEN
             write(lunout,*)' ********  PBS DANS  CONF_GCM  ******** '
-            write(lunout,*)' *** ysinus lu sur le fichier start est T',
-     *        ' alors  qu il est  F  sur  run.def  ****  '
+            write(lunout,*)' *** ysinus lu sur le fichier start est T',     &
+     &        ' alors  qu il est  F  sur  run.def  ****  '
               STOP
           ENDIF
         ENDIF
       ENDIF ! of IF( .NOT.fxyhypb  )
-c
+!
 !Config  Key  = offline
 !Config  Desc = Nouvelle eau liquide
@@ -682 +682 @@
 
       RETURN
-c   ...............................................
-c
+!   ...............................................
+!
 100   CONTINUE
 !Config  Key  = clon
@@ -718 +718 @@
 
       IF( grossismx.LT.1. )  THEN
-        write(lunout,*)
-     &   'conf_gcm: ***  ATTENTION !! grossismx < 1 .   *** '
+        write(lunout,*)'conf_gcm: ***ATTENTION !! grossismx < 1 . *** '
          STOP
       ELSE
@@ -727 +726 @@
 
       IF( grossismy.LT.1. )  THEN
-        write(lunout,*)
-     &  'conf_gcm: ***  ATTENTION !! grossismy < 1 .   *** '
+        write(lunout,*) 'conf_gcm: ***ATTENTION !! grossismy < 1 . *** '
          STOP
       ELSE
@@ -735 +733 @@
 
       write(lunout,*)'conf_gcm: alphax alphay ',alphax,alphay
-c
-c    alphax et alphay sont les anciennes formulat. des grossissements
-c
-c
+!
+!    alphax et alphay sont les anciennes formulat. des grossissements
+!
+!
 
 !Config  Key  = fxyhypb
@@ -786 +784 @@
        ysinus = .TRUE.
        CALL getin('ysinus',ysinus)
-c
+!
 !Config  Key  = offline
 !Config  Desc = Nouvelle eau liquide
@@ -864 +862 @@
       vert_prof_dissip = merge(1, 0, ok_strato .and. llm==39)
       CALL getin('vert_prof_dissip', vert_prof_dissip)
-      call assert(vert_prof_dissip == 0 .or. vert_prof_dissip ==  1,
-     $     "bad value for vert_prof_dissip")
+      call assert(vert_prof_dissip == 0 .or. vert_prof_dissip ==  1,        &
+     &     "bad value for vert_prof_dissip")
 
 !Config  Key  = ok_gradsfile
@@ -892 +890 @@
 
       write(lunout,*)' #########################################'
-      write(lunout,*)' Configuration des parametres de cel0'
+      write(lunout,*)' Configuration des parametres de cel0'                &
      &             //'_limit: '
       write(lunout,*)' planet_type = ', planet_type
@@ -937 +935 @@
       write(lunout,*)' ok_limit = ', ok_limit
       write(lunout,*)' ok_etat0 = ', ok_etat0
-c
+!
       RETURN
       END

LMDZ5/trunk/libf/dyn3d_common/q_sat.F

@@ -2 +2 @@
 ! $Header$
 !
-c
-c
+!
+!
 
       subroutine q_sat(np,temp,pres,qsat)
-c
+!
       IMPLICIT none
-c======================================================================
-c Autheur(s): Z.X. Li (LMD/CNRS)
-c  reecriture vectorisee par F. Hourdin.
-c Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
-c======================================================================
-c Arguments:
-c kelvin---input-R: temperature en Kelvin
-c millibar--input-R: pression en mb
-c
-c q_sat----output-R: vapeur d'eau saturante en kg/kg
-c======================================================================
-c
+!======================================================================
+! Autheur(s): Z.X. Li (LMD/CNRS)
+!  reecriture vectorisee par F. Hourdin.
+! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
+!======================================================================
+! Arguments:
+! kelvin---input-R: temperature en Kelvin
+! millibar--input-R: pression en mb
+!
+! q_sat----output-R: vapeur d'eau saturante en kg/kg
+!======================================================================
+!
       integer np
       REAL temp(np),pres(np),qsat(np)
-c
+!
       REAL r2es
       PARAMETER (r2es=611.14 *18.0153/28.9644)
-c
+!
       REAL r3les, r3ies, r3es
       PARAMETER (R3LES=17.269)
       PARAMETER (R3IES=21.875)
-c
+!
       REAL r4les, r4ies, r4es
       PARAMETER (R4LES=35.86)
       PARAMETER (R4IES=7.66)
-c
+!
       REAL rtt
       PARAMETER (rtt=273.16)
-c
+!
       REAL retv
       PARAMETER (retv=28.9644/18.0153 - 1.0)
@@ -42 +42 @@
       real zqsat
       integer ip
-c
-C     ------------------------------------------------------------------
-c
-c
+!
+!     ------------------------------------------------------------------
+!
+!
 
       do ip=1,np
 
-c      write(*,*)'kelvin,millibar=',kelvin,millibar
-c       write(*,*)'temp,pres=',temp(ip),pres(ip)
-c
+!      write(*,*)'kelvin,millibar=',kelvin,millibar
+!       write(*,*)'temp,pres=',temp(ip),pres(ip)
+!
          IF (temp(ip) .LE. rtt) THEN
             r3es = r3ies
@@ -59 +59 @@
             r4es = r4les
          ENDIF
-c
+!
          zqsat=r2es/pres(ip)*EXP(r3es*(temp(ip)-rtt)/(temp(ip)-r4es))
          zqsat=MIN(0.5,ZQSAT)
          zqsat=zqsat/(1.-retv *zqsat)
-c
+!
          qsat(ip)= zqsat
-c      write(*,*)'qsat=',qsat(ip)
+!      write(*,*)'qsat=',qsat(ip)
 
       enddo
-c
+!
       RETURN
       END

LMDZ5/trunk/libf/phylmd/1DUTILS.h

@@ -5 +5 @@
 ! $Id: conf_unicol.F 1279 2010-08-04 17:20:56Z lahellec $
 !
-c
-c
+!
+!
       SUBROUTINE conf_unicol
-c
+!
 #ifdef CPP_IOIPSL
       use IOIPSL
@@ -16 +16 @@
 #endif
       IMPLICIT NONE
-c-----------------------------------------------------------------------
-c     Auteurs :   A. Lahellec  .
-c
-c   Declarations :
-c   --------------
+!-----------------------------------------------------------------------
+!     Auteurs :   A. Lahellec  .
+!
+!   Declarations :
+!   --------------
 
 #include "compar1d.h"
@@ -28 +28 @@
 #include "fcg_racmo.h"
 #include "iniprint.h"
-c
-c
-c   local:
-c   ------
-
-c      CHARACTER ch1*72,ch2*72,ch3*72,ch4*12
+!
+!
+!   local:
+!   ------
+
+!      CHARACTER ch1*72,ch2*72,ch3*72,ch4*12
       
-c
-c  -------------------------------------------------------------------
-c
-c      .........    Initilisation parametres du lmdz1D      ..........
-c
-c---------------------------------------------------------------------
-c   initialisations:
-c   ----------------
+!
+!  -------------------------------------------------------------------
+!
+!      .........    Initilisation parametres du lmdz1D      ..........
+!
+!---------------------------------------------------------------------
+!   initialisations:
+!   ----------------
 
 !Config  Key  = lunout
@@ -60 +60 @@
 !Config  Help = Niveau d'impression pour le débogage
 !Config         (0 = minimum d'impression)
-c      prt_level = 0
-c      CALL getin('prt_level',prt_level)
-
-c-----------------------------------------------------------------------
-c  Parametres de controle du run:
-c-----------------------------------------------------------------------
+!      prt_level = 0
+!      CALL getin('prt_level',prt_level)
+
+!-----------------------------------------------------------------------
+!  Parametres de controle du run:
+!-----------------------------------------------------------------------
 
 !Config  Key  = restart
@@ -71 +71 @@
 !Config  Def  = false
 !Config  Help = les fichiers restart doivent etre renomme en start
-       restart = .FALSE.
+       restart =.false.
        CALL getin('restart',restart)
 
@@ -137 +137 @@
 !Config  Def  = false
 !Config  Help = forcage ou non par les flux de surface
-       ok_flux_surf = .FALSE.
+       ok_flux_surf =.false.
        CALL getin('ok_flux_surf',ok_flux_surf)
 
@@ -319 +319 @@
       write(lunout,*)' +++++++++++++++++++++++++++++++++++++++'
       write(lunout,*)
-c
+!
       RETURN
       END
@@ -325 +325 @@
 ! $Id: dyn1deta0.F 1279 2010/07/30 A Lahellec$
 !
-c
-      SUBROUTINE dyn1deta0(fichnom,plev,play,phi,phis,presnivs,
+!
+      SUBROUTINE dyn1deta0(fichnom,plev,play,phi,phis,presnivs,                 &
      &                          ucov,vcov,temp,q,omega2)
       USE dimphy
@@ -339 +339 @@
 
       IMPLICIT NONE
-c=======================================================
-c Ecriture du fichier de redemarrage sous format NetCDF
-c=======================================================
-c   Declarations:
-c   -------------
+!=======================================================
+! Ecriture du fichier de redemarrage sous format NetCDF
+!=======================================================
+!   Declarations:
+!   -------------
 #include "dimensions.h"
 #include "comconst.h"
@@ -351 +351 @@
 #include "netcdf.inc"
 
-c   Arguments:
-c   ----------
+!   Arguments:
+!   ----------
       CHARACTER*(*) fichnom
-cAl1 plev tronque pour .nc mais plev(klev+1):=0
+!Al1 plev tronque pour .nc mais plev(klev+1):=0
       real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
       real :: presnivs(klon,klev)
       real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
       real :: q(klon,klev,nqtot),omega2(klon,klev)
-c      real :: ug(klev),vg(klev),fcoriolis
+!      real :: ug(klev),vg(klev),fcoriolis
       real :: phis(klon) 
 
-c   Variables locales pour NetCDF:
-c   ------------------------------
+!   Variables locales pour NetCDF:
+!   ------------------------------
       INTEGER iq
       INTEGER length
@@ -383 +383 @@
       print*,'after open startphy ',fichnom,nmq
 
-c
-c Lecture des parametres de controle:
-c
+!
+! Lecture des parametres de controle:
+!
       CALL get_var("controle",tab_cntrl)
        
@@ -394 +394 @@
       day_ref    = tab_cntrl(4)
       annee_ref  = tab_cntrl(5)
-c      rad        = tab_cntrl(6)
-c      omeg       = tab_cntrl(7)
-c      g          = tab_cntrl(8)
-c      cpp        = tab_cntrl(9)
-c      kappa      = tab_cntrl(10)
-c      daysec     = tab_cntrl(11)
-c      dtvr       = tab_cntrl(12)
-c      etot0      = tab_cntrl(13)
-c      ptot0      = tab_cntrl(14)
-c      ztot0      = tab_cntrl(15)
-c      stot0      = tab_cntrl(16)
-c      ang0       = tab_cntrl(17)
-c      pa         = tab_cntrl(18)
-c      preff      = tab_cntrl(19)
-c
-c      clon       = tab_cntrl(20)
-c      clat       = tab_cntrl(21)
-c      grossismx  = tab_cntrl(22)
-c      grossismy  = tab_cntrl(23)
-c
+!      rad        = tab_cntrl(6)
+!      omeg       = tab_cntrl(7)
+!      g          = tab_cntrl(8)
+!      cpp        = tab_cntrl(9)
+!      kappa      = tab_cntrl(10)
+!      daysec     = tab_cntrl(11)
+!      dtvr       = tab_cntrl(12)
+!      etot0      = tab_cntrl(13)
+!      ptot0      = tab_cntrl(14)
+!      ztot0      = tab_cntrl(15)
+!      stot0      = tab_cntrl(16)
+!      ang0       = tab_cntrl(17)
+!      pa         = tab_cntrl(18)
+!      preff      = tab_cntrl(19)
+!
+!      clon       = tab_cntrl(20)
+!      clat       = tab_cntrl(21)
+!      grossismx  = tab_cntrl(22)
+!      grossismy  = tab_cntrl(23)
+!
       IF ( tab_cntrl(24).EQ.1. )  THEN
-        fxyhypb  = . TRUE .
-c        dzoomx   = tab_cntrl(25)
-c        dzoomy   = tab_cntrl(26)
-c        taux     = tab_cntrl(28)
-c        tauy     = tab_cntrl(29)
+        fxyhypb  =.true.
+!        dzoomx   = tab_cntrl(25)
+!        dzoomy   = tab_cntrl(26)
+!        taux     = tab_cntrl(28)
+!        tauy     = tab_cntrl(29)
       ELSE
-        fxyhypb = . FALSE .
-        ysinus  = . FALSE .
-        IF( tab_cntrl(27).EQ.1. ) ysinus = . TRUE. 
+        fxyhypb = .false.
+        ysinus  = .false.
+        IF( tab_cntrl(27).EQ.1. ) ysinus =.true. 
       ENDIF
 
       day_ini = tab_cntrl(30)
       itau_dyn = tab_cntrl(31)
-c   .................................................................
-c
-c
-c      PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
-cAl1
-       Print*,'day_ref,annee_ref,day_ini,itau_dyn',
+!   .................................................................
+!
+!
+!      PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
+!Al1
+       Print*,'day_ref,annee_ref,day_ini,itau_dyn',                         &
      &              day_ref,annee_ref,day_ini,itau_dyn
 
-c  Lecture des champs
-c
+!  Lecture des champs
+!
       plev(1,klev+1)=0.
       CALL get_field("plev",plev,found)
@@ -460 +460 @@
       Do iq=1,nqtot
         CALL get_field("q"//nmq(iq),q(:,:,iq),found)
-        IF (.NOT. found)
-     &  PRINT*, modname//'Le champ <q'//nmq//'> est absent'
+        IF (.NOT.found)PRINT*, modname//'Le champ <q'//nmq//'> est absent'
       EndDo
 
       CALL close_startphy
-      print*,' close startphy'
-     .      ,fichnom,play(1,1),play(1,klev),temp(1,klev)
-c
+      print*,' close startphy',fichnom,play(1,1),play(1,klev),temp(1,klev)
+!
       RETURN
       END
@@ -473 +471 @@
 ! $Id: dyn1dredem.F 1279 2010/07/29 A Lahellec$
 !
-c
-      SUBROUTINE dyn1dredem(fichnom,plev,play,phi,phis,presnivs,
+!
+      SUBROUTINE dyn1dredem(fichnom,plev,play,phi,phis,presnivs,           &
      &                          ucov,vcov,temp,q,omega2)
       USE dimphy
@@ -485 +483 @@
 
       IMPLICIT NONE
-c=======================================================
-c Ecriture du fichier de redemarrage sous format NetCDF
-c=======================================================
-c   Declarations:
-c   -------------
+!=======================================================
+! Ecriture du fichier de redemarrage sous format NetCDF
+!=======================================================
+!   Declarations:
+!   -------------
 #include "dimensions.h"
 #include "comconst.h"
@@ -497 +495 @@
 #include "netcdf.inc"
 
-c   Arguments:
-c   ----------
+!   Arguments:
+!   ----------
       CHARACTER*(*) fichnom
-cAl1 plev tronque pour .nc mais plev(klev+1):=0
+!Al1 plev tronque pour .nc mais plev(klev+1):=0
       real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
       real :: presnivs(klon,klev)
@@ -506 +504 @@
       real :: q(klon,klev,nqtot)
       real :: omega2(klon,klev),rho(klon,klev+1)
-c      real :: ug(klev),vg(klev),fcoriolis
+!      real :: ug(klev),vg(klev),fcoriolis
       real :: phis(klon) 
 
-c   Variables locales pour NetCDF:
-c   ------------------------------
+!   Variables locales pour NetCDF:
+!   ------------------------------
       INTEGER nid
       INTEGER ierr
@@ -520 +518 @@
       character*20 modname
       character*80 abort_message
-c
+!
       INTEGER nb
       SAVE nb
@@ -554 +552 @@
        tab_cntrl(11) = daysec
        tab_cntrl(12) = dtvr
-c       tab_cntrl(13) = etot0
-c       tab_cntrl(14) = ptot0
-c       tab_cntrl(15) = ztot0
-c       tab_cntrl(16) = stot0
-c       tab_cntrl(17) = ang0
-c       tab_cntrl(18) = pa
-c       tab_cntrl(19) = preff
-c
-c    .....    parametres  pour le zoom      ......   
-
-c       tab_cntrl(20)  = clon
-c       tab_cntrl(21)  = clat
-c       tab_cntrl(22)  = grossismx
-c       tab_cntrl(23)  = grossismy
-c
+!       tab_cntrl(13) = etot0
+!       tab_cntrl(14) = ptot0
+!       tab_cntrl(15) = ztot0
+!       tab_cntrl(16) = stot0
+!       tab_cntrl(17) = ang0
+!       tab_cntrl(18) = pa
+!       tab_cntrl(19) = preff
+!
+!    .....    parametres  pour le zoom      ......   
+
+!       tab_cntrl(20)  = clon
+!       tab_cntrl(21)  = clat
+!       tab_cntrl(22)  = grossismx
+!       tab_cntrl(23)  = grossismy
+!
       IF ( fxyhypb )   THEN
        tab_cntrl(24) = 1.
-c       tab_cntrl(25) = dzoomx
-c       tab_cntrl(26) = dzoomy
+!       tab_cntrl(25) = dzoomx
+!       tab_cntrl(26) = dzoomy
        tab_cntrl(27) = 0.
-c       tab_cntrl(28) = taux
-c       tab_cntrl(29) = tauy
+!       tab_cntrl(28) = taux
+!       tab_cntrl(29) = tauy
       ELSE
        tab_cntrl(24) = 0.
-c       tab_cntrl(25) = dzoomx
-c       tab_cntrl(26) = dzoomy
+!       tab_cntrl(25) = dzoomx
+!       tab_cntrl(26) = dzoomy
        tab_cntrl(27) = 0.
        tab_cntrl(28) = 0.
@@ -585 +583 @@
        IF( ysinus )  tab_cntrl(27) = 1.
       ENDIF
-CAl1 iday_end -> day_end
+!Al1 iday_end -> day_end
        tab_cntrl(30) = FLOAT(day_end)
        tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
-c
+!
       CALL put_var("controle","Param. de controle Dyn1D",tab_cntrl)
-c
-
-c  Ecriture/extension de la coordonnee temps
+!
+
+!  Ecriture/extension de la coordonnee temps
 
       nb = nb + 1
 
-c  Ecriture des champs
-c
+!  Ecriture des champs
+!
       CALL put_field("plev","p interfaces sauf la nulle",plev)
       CALL put_field("play","",play)
@@ -609 +607 @@
 
       Do iq=1,nqtot
-        CALL put_field("q"//nmq(iq),"eau vap ou condens et traceurs",
-     .                                                      q(:,:,iq))
+        CALL put_field("q"//nmq(iq),"eau vap ou condens et traceurs",           &
+     &                                                      q(:,:,iq))
       EndDo
       CALL close_restartphy
 
-c
+!
       RETURN
       END
@@ -770 +768 @@
 !   -------
  
-      IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1)
-     s    STOP 'probleme de dim'
+      IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1)                          &
+     &    STOP 'probleme de dim'
 !   traitement des poles
       CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
@@ -829 +827 @@
       pi=2.*ASIN(1.)
  
-      OPEN(99,file='sigma.def',status='old',form='formatted',
-     s   iostat=ierr)
+      OPEN(99,file='sigma.def',status='old',form='formatted',                   &
+     &   iostat=ierr)
  
 !-----------------------------------------------------------------------
@@ -850 +848 @@
  
        DO 1  l = 1, llm
-       dsig(l) = (alpha+(1.-alpha)*exp(-beta*(llm-l)))*
-     $          ( (tanh(gama*l)/tanh(gama*llm))**np +
-     $            (1.-l/FLOAT(llm))*delta )
+       dsig(l) = (alpha+(1.-alpha)*exp(-beta*(llm-l)))*                    &
+     &          ( (tanh(gama*l)/tanh(gama*llm))**np +                      &
+     &            (1.-l/FLOAT(llm))*delta )
    1   CONTINUE
  
@@ -1005 +1003 @@
 
 
-       SUBROUTINE advect_va(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,
-     s                q,temp,u,v,play)
+       SUBROUTINE advect_va(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,              &
+     &                q,temp,u,v,play)
 !itlmd 
 !----------------------------------------------------------------------
@@ -1031 +1029 @@
 !si omgup pour la couche 1, alors tendance nulle
         omgdown=max(omega(2),0.0)
-        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*
-     &           (1.+q(l,1)))
-        d_t_va(l)= alpha*(omgdown)- 
-     &              omgdown*(temp(l)-temp(l+1))
-     &              /(play(l)-play(l+1))              
-
-        d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))
-     &              /(play(l)-play(l+1))              
-
-        d_u_va(l)= -omgdown*(u(l)-u(l+1))
-     &              /(play(l)-play(l+1))              
-        d_v_va(l)= -omgdown*(v(l)-v(l+1))
-     &              /(play(l)-play(l+1))              
+        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*(1.+q(l,1)))
+        d_t_va(l)= alpha*(omgdown)-omgdown*(temp(l)-temp(l+1))             &
+     &       /(play(l)-play(l+1))
+
+        d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))/(play(l)-play(l+1))              
+
+        d_u_va(l)= -omgdown*(u(l)-u(l+1))/(play(l)-play(l+1))              
+        d_v_va(l)= -omgdown*(v(l)-v(l+1))/(play(l)-play(l+1))              
 
         
        elseif(l.eq.llm) then
         omgup=min(omega(l),0.0)
-        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*
-     &           (1.+q(l,1)))
-        d_t_va(l)= alpha*(omgup)- 
+        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*(1.+q(l,1)))
+        d_t_va(l)= alpha*(omgup)-                                          &
+
 !bug?     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
      &              omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
@@ -1060 +1053 @@
         omgup=min(omega(l),0.0)
         omgdown=max(omega(l+1),0.0)
-        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*
-     &           (1.+q(l,1)))
-        d_t_va(l)= alpha*(omgup+omgdown)- 
-     &              omgdown*(temp(l)-temp(l+1))
-     &              /(play(l)-play(l+1))-              
+        alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)*(1.+q(l,1)))
+        d_t_va(l)= alpha*(omgup+omgdown)-omgdown*(temp(l)-temp(l+1))       &
+     &              /(play(l)-play(l+1))-                                  &
 !bug?     &              omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l))
      &              omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l))
 !      print*, '  ??? '
 
-        d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))
-     &              /(play(l)-play(l+1))-              
-     &              omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
-        d_u_va(l)= -omgdown*(u(l)-u(l+1))
-     &              /(play(l)-play(l+1))-              
-     &              omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
-        d_v_va(l)= -omgdown*(v(l)-v(l+1))
-     &              /(play(l)-play(l+1))-              
+        d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))                            &
+     &              /(play(l)-play(l+1))-                                  &
+     &              omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l)) 
+        d_u_va(l)= -omgdown*(u(l)-u(l+1))                                  &
+     &              /(play(l)-play(l+1))-                                  &
+     &              omgup*(u(l-1)-u(l))/(play(l-1)-play(l)) 
+        d_v_va(l)= -omgdown*(v(l)-v(l+1))                                  &
+     &              /(play(l)-play(l+1))-                                  &
      &              omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
        
@@ -1086 +1077 @@
         end
 !       SUBROUTINE lstendH(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,
-       SUBROUTINE lstendH(llm,nqtot,omega,d_t_va,d_q_va,
-     !                q,temp,u,v,play)
+       SUBROUTINE lstendH(llm,nqtot,omega,d_t_va,d_q_va,                        &
+     &                q,temp,u,v,play)
 !itlmd 
 !----------------------------------------------------------------------
@@ -1143 +1134 @@
 !      d_v_va(k)  = -omdn*dvdp(k)-omup*dvdp(k-1)
       d_q_va(k,1)= -omdn*dqdp(k)-omup*dqdp(k-1)
-      d_t_va(k)  = -omdn*dtdp(k)-omup*dtdp(k-1)
-     :              +(omup+omdn)*cor(k)
+      d_t_va(k)  = -omdn*dtdp(k)-omup*dtdp(k-1)+(omup+omdn)*cor(k)
       enddo
 
@@ -1170 +1160 @@
 
 !======================================================================
-      SUBROUTINE read_togacoare(fich_toga,nlev_toga,nt_toga
-     :             ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga
-     :             ,ht_toga,vt_toga,hq_toga,vq_toga)
+      SUBROUTINE read_togacoare(fich_toga,nlev_toga,nt_toga                     &
+     &             ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga        &
+     &             ,ht_toga,vt_toga,hq_toga,vq_toga)
       implicit none
 
-c-------------------------------------------------------------------------
-c Read TOGA-COARE forcing data 
-c-------------------------------------------------------------------------
+!-------------------------------------------------------------------------
+! Read TOGA-COARE forcing data 
+!-------------------------------------------------------------------------
 
       integer nlev_toga,nt_toga
@@ -1207 +1197 @@
 
        do k = 1, nlev_toga
-         read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip)
-     :       ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip)
-     :       ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip)
+         read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip)          &
+     &       ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip)                     &
+     &       ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip)
 
 ! conversion in SI units:
@@ -1233 +1223 @@
           end
 
-c-------------------------------------------------------------------------
+!-------------------------------------------------------------------------
       SUBROUTINE read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
       implicit none
 
-c-------------------------------------------------------------------------
-c Read I.SANDU case forcing data
-c-------------------------------------------------------------------------
+!-------------------------------------------------------------------------
+! Read I.SANDU case forcing data
+!-------------------------------------------------------------------------
 
       integer nlev_sandu,nt_sandu
@@ -1269 +1259 @@
 
 !=====================================================================
-c-------------------------------------------------------------------------
-      SUBROUTINE read_astex(fich_astex,nlev_astex,nt_astex,div_astex,
-     : ts_astex,ug_astex,vg_astex,ufa_astex,vfa_astex)
+!-------------------------------------------------------------------------
+      SUBROUTINE read_astex(fich_astex,nlev_astex,nt_astex,div_astex,      &
+     & ts_astex,ug_astex,vg_astex,ufa_astex,vfa_astex)
       implicit none
 
-c-------------------------------------------------------------------------
-c Read Astex case forcing data
-c-------------------------------------------------------------------------
+!-------------------------------------------------------------------------
+! Read Astex case forcing data
+!-------------------------------------------------------------------------
 
       integer nlev_astex,nt_astex
@@ -1297 +1287 @@
       read(21,'(a)')
       read(21,'(a)')
-      read(21,223) iy, im, id, ih, div_astex(ip),ts_astex(ip),
-     :ug_astex(ip),vg_astex(ip),ufa_astex(ip),vfa_astex(ip)
+      read(21,223) iy, im, id, ih, div_astex(ip),ts_astex(ip),             &
+     &ug_astex(ip),vg_astex(ip),ufa_astex(ip),vfa_astex(ip)
       ts_astex(ip)=ts_astex(ip)+273.15
-      print *,'ts=',iy,im,id,ih,ip,div_astex(ip),ts_astex(ip),
-     :ug_astex(ip),vg_astex(ip),ufa_astex(ip),vg_astex(ip)
+      print *,'ts=',iy,im,id,ih,ip,div_astex(ip),ts_astex(ip),             &
+     &ug_astex(ip),vg_astex(ip),ufa_astex(ip),vg_astex(ip)
       enddo
       close(21)
@@ -1310 +1300 @@
           end
 !=====================================================================
-      subroutine read_twpice(fich_twpice,nlevel,ntime
-     :     ,T_srf,plev,T,q,u,v,omega
-     :     ,T_adv_h,T_adv_v,q_adv_h,q_adv_v)
+      subroutine read_twpice(fich_twpice,nlevel,ntime                       &
+     &     ,T_srf,plev,T,q,u,v,omega                                       &
+     &     ,T_adv_h,T_adv_v,q_adv_h,q_adv_v)
 
 !program reading forcings of the TWP-ICE experiment
@@ -1792 +1782 @@
 !=====================================================================
 
-       SUBROUTINE interp_sandu_vertical(play,nlev_sandu,plev_prof
-     :         ,t_prof,thl_prof,q_prof,u_prof,v_prof,w_prof
-     :         ,omega_prof,o3mmr_prof
-     :         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod
-     :         ,omega_mod,o3mmr_mod,mxcalc)
+       SUBROUTINE interp_sandu_vertical(play,nlev_sandu,plev_prof          &
+     &         ,t_prof,thl_prof,q_prof,u_prof,v_prof,w_prof                &
+     &         ,omega_prof,o3mmr_prof                                      &
+     &         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod                      &
+     &         ,omega_mod,o3mmr_mod,mxcalc)
 
        implicit none
@@ -1802 +1792 @@
 #include "dimensions.h"
 
-c-------------------------------------------------------------------------
-c Vertical interpolation of SANDUREF forcing data onto model levels
-c-------------------------------------------------------------------------
+!-------------------------------------------------------------------------
+! Vertical interpolation of SANDUREF forcing data onto model levels
+!-------------------------------------------------------------------------
 
        integer nlevmax
@@ -1838 +1828 @@
 
          do k = 1, nlev_sandu-1
-          if (play(l).le.plev_prof(k)
-     :       .and. play(l).gt.plev_prof(k+1)) then
+          if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
             k1=k
             k2=k+1
@@ -1860 +1849 @@
          v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
          w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
-         omega_mod(l)=omega_prof(k2)-
-     :      frac*(omega_prof(k2)-omega_prof(k1))
-         o3mmr_mod(l)=o3mmr_prof(k2)-
-     :      frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
+         omega_mod(l)=omega_prof(k2)-frac*(omega_prof(k2)-omega_prof(k1))
+         o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
 
          else !play>plev_prof(1)
@@ -1884 +1871 @@
         else ! above max altitude of forcing file
 
-cjyg
+!jyg
          fact=20.*(plev_prof(nlev_sandu)-play(l))/plev_prof(nlev_sandu) !jyg
          fact = max(fact,0.)                                           !jyg
@@ -1909 +1896 @@
           end
 !=====================================================================
-       SUBROUTINE interp_astex_vertical(play,nlev_astex,plev_prof
-     :         ,t_prof,thl_prof,qv_prof,ql_prof,qt_prof,u_prof,v_prof
-     :         ,w_prof,tke_prof,o3mmr_prof
-     :         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod
-     :         ,tke_mod,o3mmr_mod,mxcalc)
+       SUBROUTINE interp_astex_vertical(play,nlev_astex,plev_prof          &
+     &         ,t_prof,thl_prof,qv_prof,ql_prof,qt_prof,u_prof,v_prof      &
+     &         ,w_prof,tke_prof,o3mmr_prof                                 &
+     &         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod       &
+     &         ,tke_mod,o3mmr_mod,mxcalc)
 
        implicit none
@@ -1919 +1906 @@
 #include "dimensions.h"
 
-c-------------------------------------------------------------------------
-c Vertical interpolation of Astex forcing data onto model levels
-c-------------------------------------------------------------------------
+!-------------------------------------------------------------------------
+! Vertical interpolation of Astex forcing data onto model levels
+!-------------------------------------------------------------------------
 
        integer nlevmax
@@ -1956 +1943 @@
 
          do k = 1, nlev_astex-1
-          if (play(l).le.plev_prof(k)
-     :       .and. play(l).gt.plev_prof(k+1)) then
+          if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
             k1=k
             k2=k+1
@@ -1981 +1967 @@
          w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
          tke_mod(l)= tke_prof(k2) - frac*(tke_prof(k2)-tke_prof(k1))
-         o3mmr_mod(l)=o3mmr_prof(k2)-
-     :      frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
+         o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
 
          else !play>plev_prof(1)
@@ -2005 +1990 @@
         else ! above max altitude of forcing file
 
-cjyg
+!jyg
          fact=20.*(plev_prof(nlev_astex)-play(l))/plev_prof(nlev_astex) !jyg
          fact = max(fact,0.)                                           !jyg
@@ -2033 +2018 @@
 
 !======================================================================
-      SUBROUTINE read_rico(fich_rico,nlev_rico,ps_rico,play
-     :             ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico
-     :             ,dth_dyn,dqh_dyn)
+      SUBROUTINE read_rico(fich_rico,nlev_rico,ps_rico,play                &
+     &             ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico             &
+     &             ,dth_dyn,dqh_dyn)
       implicit none
 
-c-------------------------------------------------------------------------
-c Read RICO forcing data 
-c-------------------------------------------------------------------------
+!-------------------------------------------------------------------------
+! Read RICO forcing data 
+!-------------------------------------------------------------------------
 #include "dimensions.h"
 
@@ -2081 +2066 @@
           if(prico(l)>play(k)) then
               if(play(k)>prico(l+1)) then
-                zlay(k)=zrico(l)+(play(k)-prico(l)) * 
+                zlay(k)=zrico(l)+(play(k)-prico(l)) *                      &
      &              (zrico(l+1)-zrico(l))/(prico(l+1)-prico(l))
-              else
-                zlay(k)=zrico(l)+(play(k)-prico(80))* 
+              else 
+                zlay(k)=zrico(l)+(play(k)-prico(80))*                      &
      &              (zrico(81)-zrico(80))/(prico(81)-prico(80))
               endif
@@ -2094 +2079 @@
           u_rico(k)=-9.9 + (-1.9 + 9.9)*zlay(k)/4000
         elseif(4000 < zlay(k) .and. zlay(k) < 12000) then
-       u_rico(k)=  -1.9 + (30.0 + 1.9) / 
-     :          (12000 - 4000) * (zlay(k) - 4000)
+       u_rico(k)=  -1.9 + (30.0 + 1.9) /                                   &
+     &          (12000 - 4000) * (zlay(k) - 4000)
         elseif(12000 < zlay(k) .and. zlay(k) < 13000) then
           u_rico(k)=30.0
         elseif(13000 < zlay(k) .and. zlay(k) < 20000) then
-          u_rico(k)=30.0 - (30.0) / 
-     : (20000 - 13000) * (zlay(k) - 13000)
+          u_rico(k)=30.0 - (30.0) /                                        &
+     & (20000 - 13000) * (zlay(k) - 13000)
         else
           u_rico(k)=0.0
@@ -2109 +2094 @@
           q_rico(k)=16.0 + (13.8 - 16.0) / (740) * zlay(k)
         elseif(740 < zlay(k) .and. zlay(k) < 3260) then
-          q_rico(k)=13.8 + (2.4 - 13.8) / 
-     :          (3260 - 740) * (zlay(k) - 740)
+          q_rico(k)=13.8 + (2.4 - 13.8) /                                   &
+     &          (3260 - 740) * (zlay(k) - 740)
         elseif(3260 < zlay(k) .and. zlay(k) < 4000) then
-          q_rico(k)=2.4 + (1.8 - 2.4) / 
-     :               (4000 - 3260) * (zlay(k) - 3260)
+          q_rico(k)=2.4 + (1.8 - 2.4) /                                    &
+     &               (4000 - 3260) * (zlay(k) - 3260)
         elseif(4000 < zlay(k) .and. zlay(k) < 9000) then
-          q_rico(k)=1.8 + (0 - 1.8) / 
-     :             (10000 - 4000) * (zlay(k) - 4000)
+          q_rico(k)=1.8 + (0 - 1.8) /                                      &
+     &             (10000 - 4000) * (zlay(k) - 4000)
         else
           q_rico(k)=0.0
@@ -2125 +2110 @@
           t_rico(k)=299.2 + (292.0 - 299.2) / (740) * zlay(k)
         elseif(740 < zlay(k) .and. zlay(k) < 4000) then
-          t_rico(k)=292.0 + (278.0 - 292.0) /
-     :       (4000 - 740) * (zlay(k) - 740)
+          t_rico(k)=292.0 + (278.0 - 292.0) /                              &                       
+     &       (4000 - 740) * (zlay(k) - 740)
         elseif(4000 < zlay(k) .and. zlay(k) < 15000) then
-          t_rico(k)=278.0 + (203.0 - 278.0) /
-     :       (15000 - 4000) * (zlay(k) - 4000)
+          t_rico(k)=278.0 + (203.0 - 278.0) /                              &
+     &       (15000 - 4000) * (zlay(k) - 4000)
         elseif(15000 < zlay(k) .and. zlay(k) < 17500) then
-          t_rico(k)=203.0 + (194.0 - 203.0) / 
-     :       (17500 - 15000)* (zlay(k) - 15000)
+          t_rico(k)=203.0 + (194.0 - 203.0) /                              & 
+     &       (17500 - 15000)* (zlay(k) - 15000)
         elseif(17500 < zlay(k) .and. zlay(k) < 20000) then
-          t_rico(k)=194.0 + (206.0 - 194.0) /
-     :       (20000 - 17500)* (zlay(k) - 17500)
+          t_rico(k)=194.0 + (206.0 - 194.0) /                              &
+     &       (20000 - 17500)* (zlay(k) - 17500)
         elseif(20000 < zlay(k) .and. zlay(k) < 60000) then
-          t_rico(k)=206.0 + (270.0 - 206.0) / 
-     :        (60000 - 20000)* (zlay(k) - 20000)
+          t_rico(k)=206.0 + (270.0 - 206.0) /                              & 
+     &        (60000 - 20000)* (zlay(k) - 20000)
         endif
 
@@ -2154 +2139 @@
 ! dThrz+dTsw0+dTlw0
         if(0 < zlay(k) .and. zlay(k) < 4000) then
-          dth_dyn(k)=- 2.51 / 86400 + (-2.18 + 2.51 )/ 
-     :               (86400*4000) * zlay(k)
+          dth_dyn(k)=- 2.51 / 86400 + (-2.18 + 2.51 )/                     &
+     &               (86400*4000) * zlay(k)
         elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
-          dth_dyn(k)=- 2.18 / 86400 + ( 2.18 ) /
-     :           (86400*(5000 - 4000)) * (zlay(k) - 4000)
+          dth_dyn(k)=- 2.18 / 86400 + ( 2.18 ) /                           &
+     &           (86400*(5000 - 4000)) * (zlay(k) - 4000)
         else
           dth_dyn(k)=0.0
@@ -2164 +2149 @@
 ! dQhrz
         if(0 < zlay(k) .and. zlay(k) < 3000) then
-          dqh_dyn(k)=-1.0 / 86400 + (0.345 + 1.0)/
-     :                    (86400*3000) * (zlay(k))
+          dqh_dyn(k)=-1.0 / 86400 + (0.345 + 1.0)/                         &
+     &                    (86400*3000) * (zlay(k))
         elseif(3000 < zlay(k) .and. zlay(k) < 4000) then
           dqh_dyn(k)=0.345 / 86400
         elseif(4000 < zlay(k) .and. zlay(k) < 5000) then
-          dqh_dyn(k)=0.345 / 86400 + 
-     +   (-0.345)/(86400 * (5000 - 4000)) * (zlay(k)-4000)
+          dqh_dyn(k)=0.345 / 86400 +                                       &
+     &   (-0.345)/(86400 * (5000 - 4000)) * (zlay(k)-4000)
         else
           dqh_dyn(k)=0.0
@@ -2196 +2181 @@
 
 !======================================================================
-        SUBROUTINE interp_sandu_time(day,day1,annee_ref
-     i             ,year_ini_sandu,day_ini_sandu,nt_sandu,dt_sandu
-     i             ,nlev_sandu,ts_sandu,ts_prof)
+        SUBROUTINE interp_sandu_time(day,day1,annee_ref                    &
+     &             ,year_ini_sandu,day_ini_sandu,nt_sandu,dt_sandu         &
+     &             ,nlev_sandu,ts_sandu,ts_prof)
         implicit none
 
@@ -2247 +2232 @@
        time_sandu2=(it_sandu2-1)*dt_sandu
        print *,'timeit day day_ini_sandu',timeit,day,day_ini_sandu
-       print *,'it_sandu1,it_sandu2,time_sandu1,time_sandu2',
-     .          it_sandu1,it_sandu2,time_sandu1,time_sandu2
+       print *,'it_sandu1,it_sandu2,time_sandu1,time_sandu2',              &
+     &          it_sandu1,it_sandu2,time_sandu1,time_sandu2
 
        if (it_sandu1 .ge. nt_sandu) then
-        write(*,*) 'PB-stop: day, it_sandu1, it_sandu2, timeit: '
-     :        ,day,it_sandu1,it_sandu2,timeit/86400.
+        write(*,*) 'PB-stop: day, it_sandu1, it_sandu2, timeit: '          &
+     &        ,day,it_sandu1,it_sandu2,timeit/86400.
         stop
        endif
@@ -2260 +2245 @@
        frac=max(frac,0.0)
 
-       ts_prof = ts_sandu(it_sandu2)
-     :          -frac*(ts_sandu(it_sandu2)-ts_sandu(it_sandu1))
-
-         print*,
-     :'day,annee_ref,day_ini_sandu,timeit,it_sandu1,it_sandu2,SST:',
-     :day,annee_ref,day_ini_sandu,timeit/86400.,it_sandu1,
-     :it_sandu2,ts_prof
+       ts_prof = ts_sandu(it_sandu2)                                       &
+     &          -frac*(ts_sandu(it_sandu2)-ts_sandu(it_sandu1))
+
+         print*,                                                           &
+     &'day,annee_ref,day_ini_sandu,timeit,it_sandu1,it_sandu2,SST:',       &
+     &day,annee_ref,day_ini_sandu,timeit/86400.,it_sandu1,                  &
+     &it_sandu2,ts_prof
 
         return
         END
 !=====================================================================
-c-------------------------------------------------------------------------
-      SUBROUTINE read_armcu(fich_armcu,nlev_armcu,nt_armcu,
-     : sens,flat,adv_theta,rad_theta,adv_qt)
+!-------------------------------------------------------------------------
+      SUBROUTINE read_armcu(fich_armcu,nlev_armcu,nt_armcu,                &
+     & sens,flat,adv_theta,rad_theta,adv_qt)
       implicit none
 
-c-------------------------------------------------------------------------
-c Read ARM_CU case forcing data
-c-------------------------------------------------------------------------
+!-------------------------------------------------------------------------
+! Read ARM_CU case forcing data
+!-------------------------------------------------------------------------
 
       integer nlev_armcu,nt_armcu
@@ -2296 +2281 @@
       read(21,'(a)')
       read(21,'(a)')
-      read(21,223) iy, im, id, ih, in, sens(ip),flat(ip),
-     :             adv_theta(ip),rad_theta(ip),adv_qt(ip)
-      print *,'forcages=',iy,im,id,ih,in, sens(ip),flat(ip),
-     :             adv_theta(ip),rad_theta(ip),adv_qt(ip)
+      read(21,223) iy, im, id, ih, in, sens(ip),flat(ip),                  &
+     &             adv_theta(ip),rad_theta(ip),adv_qt(ip)
+      print *,'forcages=',iy,im,id,ih,in, sens(ip),flat(ip),               &
+     &             adv_theta(ip),rad_theta(ip),adv_qt(ip)
       enddo
       close(21)
@@ -2309 +2294 @@
 
 !=====================================================================
-       SUBROUTINE interp_toga_vertical(play,nlev_toga,plev_prof
-     :         ,t_prof,q_prof,u_prof,v_prof,w_prof
-     :         ,ht_prof,vt_prof,hq_prof,vq_prof
-     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
-     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
+       SUBROUTINE interp_toga_vertical(play,nlev_toga,plev_prof            &
+     &         ,t_prof,q_prof,u_prof,v_prof,w_prof                         &
+     &         ,ht_prof,vt_prof,hq_prof,vq_prof                            &
+     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                              &
+     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
  
        implicit none
@@ -2319 +2304 @@
 #include "dimensions.h"
 
-c-------------------------------------------------------------------------
-c Vertical interpolation of TOGA-COARE forcing data onto model levels
-c-------------------------------------------------------------------------
+!-------------------------------------------------------------------------
+! Vertical interpolation of TOGA-COARE forcing data onto model levels
+!-------------------------------------------------------------------------
  
        integer nlevmax
@@ -2357 +2342 @@
 
          do k = 1, nlev_toga-1
-          if (play(l).le.plev_prof(k) 
-     :       .and. play(l).gt.plev_prof(k+1)) then
+          if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then
             k1=k
             k2=k+1
@@ -2403 +2387 @@
         else ! above max altitude of forcing file
  
-cjyg
+!jyg
          fact=20.*(plev_prof(nlev_toga)-play(l))/plev_prof(nlev_toga) !jyg
          fact = max(fact,0.)                                           !jyg
@@ -2430 +2414 @@
  
 !======================================================================
-        SUBROUTINE interp_astex_time(day,day1,annee_ref
-     i             ,year_ini_astex,day_ini_astex,nt_astex,dt_astex
-     i             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex
-     i             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof
-     i             ,ufa_prof,vfa_prof)
+        SUBROUTINE interp_astex_time(day,day1,annee_ref                    &
+     &             ,year_ini_astex,day_ini_astex,nt_astex,dt_astex         &
+     &             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex        &
+     &             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof   &
+     &             ,ufa_prof,vfa_prof)
         implicit none
 
@@ -2486 +2470 @@
        time_astex2=(it_astex2-1)*dt_astex
        print *,'timeit day day_ini_astex',timeit,day,day_ini_astex
-       print *,'it_astex1,it_astex2,time_astex1,time_astex2',
-     .          it_astex1,it_astex2,time_astex1,time_astex2
+       print *,'it_astex1,it_astex2,time_astex1,time_astex2',              &
+     &          it_astex1,it_astex2,time_astex1,time_astex2
 
        if (it_astex1 .ge. nt_astex) then
-        write(*,*) 'PB-stop: day, it_astex1, it_astex2, timeit: '
-     :        ,day,it_astex1,it_astex2,timeit/86400.
+        write(*,*) 'PB-stop: day, it_astex1, it_astex2, timeit: '          &
+     &        ,day,it_astex1,it_astex2,timeit/86400.
         stop
        endif
@@ -2499 +2483 @@
        frac=max(frac,0.0)
 
-       div_prof = div_astex(it_astex2)
-     :          -frac*(div_astex(it_astex2)-div_astex(it_astex1))
-       ts_prof = ts_astex(it_astex2)
-     :          -frac*(ts_astex(it_astex2)-ts_astex(it_astex1))
-       ug_prof = ug_astex(it_astex2)
-     :          -frac*(ug_astex(it_astex2)-ug_astex(it_astex1))
-       vg_prof = vg_astex(it_astex2)
-     :          -frac*(vg_astex(it_astex2)-vg_astex(it_astex1))
-       ufa_prof = ufa_astex(it_astex2)
-     :          -frac*(ufa_astex(it_astex2)-ufa_astex(it_astex1))
-       vfa_prof = vfa_astex(it_astex2)
-     :          -frac*(vfa_astex(it_astex2)-vfa_astex(it_astex1))
-
-         print*,
-     :'day,annee_ref,day_ini_astex,timeit,it_astex1,it_astex2,SST:',
-     :day,annee_ref,day_ini_astex,timeit/86400.,it_astex1,
-     :it_astex2,div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
+       div_prof = div_astex(it_astex2)                                     &
+     &          -frac*(div_astex(it_astex2)-div_astex(it_astex1))
+       ts_prof = ts_astex(it_astex2)                                        &
+     &          -frac*(ts_astex(it_astex2)-ts_astex(it_astex1))
+       ug_prof = ug_astex(it_astex2)                                       &
+     &          -frac*(ug_astex(it_astex2)-ug_astex(it_astex1))
+       vg_prof = vg_astex(it_astex2)                                       &
+     &          -frac*(vg_astex(it_astex2)-vg_astex(it_astex1))
+       ufa_prof = ufa_astex(it_astex2)                                     &
+     &          -frac*(ufa_astex(it_astex2)-ufa_astex(it_astex1))
+       vfa_prof = vfa_astex(it_astex2)                                     &
+     &          -frac*(vfa_astex(it_astex2)-vfa_astex(it_astex1))
+
+         print*,                                                           &
+     &'day,annee_ref,day_ini_astex,timeit,it_astex1,it_astex2,SST:',       &
+     &day,annee_ref,day_ini_astex,timeit/86400.,it_astex1,                 &
+     &it_astex2,div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof 
 
         return
@@ -2521 +2505 @@
 
 !======================================================================
-        SUBROUTINE interp_toga_time(day,day1,annee_ref
-     i             ,year_ini_toga,day_ini_toga,nt_toga,dt_toga,nlev_toga
-     i             ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga
-     i             ,ht_toga,vt_toga,hq_toga,vq_toga
-     o             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof
-     o             ,ht_prof,vt_prof,hq_prof,vq_prof)
+        SUBROUTINE interp_toga_time(day,day1,annee_ref                     &
+     &             ,year_ini_toga,day_ini_toga,nt_toga,dt_toga,nlev_toga   &
+     &             ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga   &
+     &             ,ht_toga,vt_toga,hq_toga,vq_toga                        &
+     &             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof   &
+     &             ,ht_prof,vt_prof,hq_prof,vq_prof)
         implicit none
 
@@ -2618 +2602 @@
 
        if (it_toga1 .ge. nt_toga) then
-        write(*,*) 'PB-stop: day, it_toga1, it_toga2, timeit: '
-     :        ,day,it_toga1,it_toga2,timeit/86400.
+        write(*,*) 'PB-stop: day, it_toga1, it_toga2, timeit: '            &
+     &        ,day,it_toga1,it_toga2,timeit/86400.
         stop
        endif
@@ -2627 +2611 @@
        frac=max(frac,0.0)
 
-       ts_prof = ts_toga(it_toga2)
-     :          -frac*(ts_toga(it_toga2)-ts_toga(it_toga1))
+       ts_prof = ts_toga(it_toga2)                                         &
+     &          -frac*(ts_toga(it_toga2)-ts_toga(it_toga1))
 
 !        print*,
@@ -2635 +2619 @@
 
        do k=1,nlev_toga
-        plev_prof(k) = 100.*(plev_toga(k,it_toga2)
-     :          -frac*(plev_toga(k,it_toga2)-plev_toga(k,it_toga1)))
-        t_prof(k) = t_toga(k,it_toga2)
-     :          -frac*(t_toga(k,it_toga2)-t_toga(k,it_toga1))
-        q_prof(k) = q_toga(k,it_toga2)
-     :          -frac*(q_toga(k,it_toga2)-q_toga(k,it_toga1))
-        u_prof(k) = u_toga(k,it_toga2)
-     :          -frac*(u_toga(k,it_toga2)-u_toga(k,it_toga1))
-        v_prof(k) = v_toga(k,it_toga2)
-     :          -frac*(v_toga(k,it_toga2)-v_toga(k,it_toga1))
-        w_prof(k) = w_toga(k,it_toga2)
-     :          -frac*(w_toga(k,it_toga2)-w_toga(k,it_toga1))
-        ht_prof(k) = ht_toga(k,it_toga2)
-     :          -frac*(ht_toga(k,it_toga2)-ht_toga(k,it_toga1))
-        vt_prof(k) = vt_toga(k,it_toga2)
-     :          -frac*(vt_toga(k,it_toga2)-vt_toga(k,it_toga1))
-        hq_prof(k) = hq_toga(k,it_toga2)
-     :          -frac*(hq_toga(k,it_toga2)-hq_toga(k,it_toga1))
-        vq_prof(k) = vq_toga(k,it_toga2)
-     :          -frac*(vq_toga(k,it_toga2)-vq_toga(k,it_toga1))
+        plev_prof(k) = 100.*(plev_toga(k,it_toga2)                         &
+     &          -frac*(plev_toga(k,it_toga2)-plev_toga(k,it_toga1)))
+        t_prof(k) = t_toga(k,it_toga2)                                     &
+     &          -frac*(t_toga(k,it_toga2)-t_toga(k,it_toga1))
+        q_prof(k) = q_toga(k,it_toga2)                                     &
+     &          -frac*(q_toga(k,it_toga2)-q_toga(k,it_toga1))
+        u_prof(k) = u_toga(k,it_toga2)                                     &
+     &          -frac*(u_toga(k,it_toga2)-u_toga(k,it_toga1))
+        v_prof(k) = v_toga(k,it_toga2)                                     &
+     &          -frac*(v_toga(k,it_toga2)-v_toga(k,it_toga1))
+        w_prof(k) = w_toga(k,it_toga2)                                     &
+     &          -frac*(w_toga(k,it_toga2)-w_toga(k,it_toga1))
+        ht_prof(k) = ht_toga(k,it_toga2)                                   &
+     &          -frac*(ht_toga(k,it_toga2)-ht_toga(k,it_toga1))
+        vt_prof(k) = vt_toga(k,it_toga2)                                   &
+     &          -frac*(vt_toga(k,it_toga2)-vt_toga(k,it_toga1))
+        hq_prof(k) = hq_toga(k,it_toga2)                                   &
+     &          -frac*(hq_toga(k,it_toga2)-hq_toga(k,it_toga1))
+        vq_prof(k) = vq_toga(k,it_toga2)                                   &
+     &          -frac*(vq_toga(k,it_toga2)-vq_toga(k,it_toga1))
         enddo
 
@@ -2661 +2645 @@
 
 !======================================================================
-        SUBROUTINE interp_armcu_time(day,day1,annee_ref
-     i             ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu
-     i             ,nlev_armcu,fs_armcu,fl_armcu,at_armcu,rt_armcu
-     i             ,aqt_armcu,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof)
+        SUBROUTINE interp_armcu_time(day,day1,annee_ref                    &
+     &             ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu         &
+     &             ,nlev_armcu,fs_armcu,fl_armcu,at_armcu,rt_armcu         &
+     &             ,aqt_armcu,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof)
         implicit none
 
@@ -2707 +2691 @@
        time_armcu2=(it_armcu2-1)*dt_armcu
        print *,'timeit day day_ini_armcu',timeit,day,day_ini_armcu
-       print *,'it_armcu1,it_armcu2,time_armcu1,time_armcu2',
-     .          it_armcu1,it_armcu2,time_armcu1,time_armcu2
+       print *,'it_armcu1,it_armcu2,time_armcu1,time_armcu2',              &
+     &          it_armcu1,it_armcu2,time_armcu1,time_armcu2
 
        if (it_armcu1 .ge. nt_armcu) then
-        write(*,*) 'PB-stop: day, it_armcu1, it_armcu2, timeit: '
-     :        ,day,it_armcu1,it_armcu2,timeit/86400.
+        write(*,*) 'PB-stop: day, it_armcu1, it_armcu2, timeit: '          &
+     &        ,day,it_armcu1,it_armcu2,timeit/86400.
         stop
        endif
@@ -2720 +2704 @@
        frac=max(frac,0.0)
 
-       fs_prof = fs_armcu(it_armcu2)
-     :          -frac*(fs_armcu(it_armcu2)-fs_armcu(it_armcu1))
-       fl_prof = fl_armcu(it_armcu2)
-     :          -frac*(fl_armcu(it_armcu2)-fl_armcu(it_armcu1))
-       at_prof = at_armcu(it_armcu2)
-     :          -frac*(at_armcu(it_armcu2)-at_armcu(it_armcu1))
-       rt_prof = rt_armcu(it_armcu2)
-     :          -frac*(rt_armcu(it_armcu2)-rt_armcu(it_armcu1))
-       aqt_prof = aqt_armcu(it_armcu2)
-     :          -frac*(aqt_armcu(it_armcu2)-aqt_armcu(it_armcu1))
-
-         print*,
-     :'day,annee_ref,day_ini_armcu,timeit,it_armcu1,it_armcu2,SST:',
-     :day,annee_ref,day_ini_armcu,timeit/86400.,it_armcu1,
-     :it_armcu2,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
+       fs_prof = fs_armcu(it_armcu2)                                       &
+     &          -frac*(fs_armcu(it_armcu2)-fs_armcu(it_armcu1))
+       fl_prof = fl_armcu(it_armcu2)                                       &
+     &          -frac*(fl_armcu(it_armcu2)-fl_armcu(it_armcu1))
+       at_prof = at_armcu(it_armcu2)                                       &
+     &          -frac*(at_armcu(it_armcu2)-at_armcu(it_armcu1))
+       rt_prof = rt_armcu(it_armcu2)                                       &
+     &          -frac*(rt_armcu(it_armcu2)-rt_armcu(it_armcu1))
+       aqt_prof = aqt_armcu(it_armcu2)                                       &
+     &          -frac*(aqt_armcu(it_armcu2)-aqt_armcu(it_armcu1))
+
+         print*,                                                           &
+     &'day,annee_ref,day_ini_armcu,timeit,it_armcu1,it_armcu2,SST:',       &
+     &day,annee_ref,day_ini_armcu,timeit/86400.,it_armcu1,                 &
+     &it_armcu2,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof
 
         return
@@ -2740 +2724 @@
 
 !=====================================================================
-      subroutine readprofiles(nlev_max,kmax,ntrac,height,
-     .           thlprof,qtprof,uprof,
-     .           vprof,e12prof,ugprof,vgprof,
-     .           wfls,dqtdxls,dqtdyls,dqtdtls,
-     .           thlpcar,tracer,nt1,nt2)
+      subroutine readprofiles(nlev_max,kmax,ntrac,height,                  &
+     &           thlprof,qtprof,uprof,                                     &
+     &           vprof,e12prof,ugprof,vgprof,                              &
+     &           wfls,dqtdxls,dqtdyls,dqtdtls,                             &
+     &           thlpcar,tracer,nt1,nt2)
       implicit none
 
@@ -2750 +2734 @@
         logical :: llesread = .true.
 
-        real height(nlev_max),thlprof(nlev_max),qtprof(nlev_max),
-     .       uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max),
-     .       ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max),
-     .       dqtdxls(nlev_max),dqtdyls(nlev_max),dqtdtls(nlev_max),
-     .           thlpcar(nlev_max),tracer(nlev_max,ntrac)
+        real height(nlev_max),thlprof(nlev_max),qtprof(nlev_max),          &
+     &       uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max),            &
+     &       ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max),             &
+     &       dqtdxls(nlev_max),dqtdyls(nlev_max),dqtdtls(nlev_max),        &
+     &           thlpcar(nlev_max),tracer(nlev_max,ntrac)
 
         integer, parameter :: ilesfile=1
@@ -2765 +2749 @@
         read (ilesfile,*) kmax
         do k=1,kmax
-          read (ilesfile,*) height(k),thlprof(k),qtprof (k),
-     .                      uprof (k),vprof  (k),e12prof(k)
+          read (ilesfile,*) height(k),thlprof(k),qtprof (k),               &
+     &                      uprof (k),vprof  (k),e12prof(k)
         enddo
         close(ilesfile)
@@ -2779 +2763 @@
         endif
         do k=1,kmax
-          read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k),
-     .                      dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)
+          read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k),         &
+     &                      dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)
         end do
         close(ilesfile)
@@ -2806 +2790 @@
         end
 !======================================================================
-      subroutine readprofile_sandu(nlev_max,kmax,height,pprof,tprof,
-     .       thlprof,qprof,uprof,vprof,wprof,omega,o3mmr)
+      subroutine readprofile_sandu(nlev_max,kmax,height,pprof,tprof,       &
+     &       thlprof,qprof,uprof,vprof,wprof,omega,o3mmr)
 !======================================================================
       implicit none
@@ -2814 +2798 @@
         logical :: llesread = .true.
 
-        real height(nlev_max),pprof(nlev_max),tprof(nlev_max),
-     .  thlprof(nlev_max),
-     .  qprof(nlev_max),uprof(nlev_max),vprof(nlev_max),
-     .  wprof(nlev_max),omega(nlev_max),o3mmr(nlev_max)
+        real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
+        real thlprof(nlev_max)
+        real qprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
+        real wprof(nlev_max),omega(nlev_max),o3mmr(nlev_max)
 
         integer, parameter :: ilesfile=1
@@ -2828 +2812 @@
         read (ilesfile,*) kmax
         do k=1,kmax
-          read (ilesfile,*) height(k),pprof(k),  tprof(k),thlprof(k),
-     .                      qprof (k),uprof(k),  vprof(k),  wprof(k),
-     .                      omega (k),o3mmr(k)
+          read (ilesfile,*) height(k),pprof(k),  tprof(k),thlprof(k),      &
+     &                      qprof (k),uprof(k),  vprof(k),  wprof(k),      &
+     &                      omega (k),o3mmr(k)
         enddo
         close(ilesfile)
@@ -2838 +2822 @@
 
 !======================================================================
-      subroutine readprofile_astex(nlev_max,kmax,height,pprof,tprof,
-     .    thlprof,qvprof,qlprof,qtprof,uprof,vprof,wprof,tkeprof,o3mmr)
+      subroutine readprofile_astex(nlev_max,kmax,height,pprof,tprof,       &
+     &    thlprof,qvprof,qlprof,qtprof,uprof,vprof,wprof,tkeprof,o3mmr)
 !======================================================================
       implicit none
@@ -2846 +2830 @@
         logical :: llesread = .true.
 
-        real height(nlev_max),pprof(nlev_max),tprof(nlev_max),
-     .  thlprof(nlev_max),qlprof(nlev_max),qtprof(nlev_max),
-     .  qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max),
-     .  wprof(nlev_max),tkeprof(nlev_max),o3mmr(nlev_max)
+        real height(nlev_max),pprof(nlev_max),tprof(nlev_max),             &
+     &  thlprof(nlev_max),qlprof(nlev_max),qtprof(nlev_max),               &
+     &  qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max),                  &
+     &  wprof(nlev_max),tkeprof(nlev_max),o3mmr(nlev_max)
 
         integer, parameter :: ilesfile=1
@@ -2860 +2844 @@
         read (ilesfile,*) kmax
         do k=1,kmax
-          read (ilesfile,*) height(k),pprof(k),  tprof(k),thlprof(k),
-     .                qvprof (k),qlprof (k),qtprof (k),
-     .                uprof(k),  vprof(k),  wprof(k),tkeprof(k),o3mmr(k)
+          read (ilesfile,*) height(k),pprof(k),  tprof(k),thlprof(k),      &
+     &                qvprof (k),qlprof (k),qtprof (k),                    &
+     &                uprof(k),  vprof(k),  wprof(k),tkeprof(k),o3mmr(k)
         enddo
         close(ilesfile)
@@ -2872 +2856 @@
 
 !======================================================================
-      subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof,
-     .       vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)
+      subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof,       &
+     &       vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)
 !======================================================================
       implicit none
@@ -2880 +2864 @@
         logical :: llesread = .true.
 
-        real height(nlev_max),pprof(nlev_max),tprof(nlev_max),
-     .  thetaprof(nlev_max),rvprof(nlev_max),
-     .  qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max),
-     .  aprof(nlev_max+1),bprof(nlev_max+1)
+        real height(nlev_max),pprof(nlev_max),tprof(nlev_max)
+        real thetaprof(nlev_max),rvprof(nlev_max)
+        real qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max)
+        real aprof(nlev_max+1),bprof(nlev_max+1)
 
         integer, parameter :: ilesfile=1
@@ -2902 +2886 @@
         read (ilesfile,*) kmax
         do k=1,kmax
-          read (ilesfile,*) height(k)    ,pprof(k),  uprof(k), vprof(k),
-     .                      thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)
+          read (ilesfile,*) height(k)    ,pprof(k),  uprof(k), vprof(k),   &
+     &                      thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)
         enddo
         close(ilesfile)
@@ -2926 +2910 @@
         end
 !=====================================================================
-      subroutine read_amma(fich_amma,nlevel,ntime
-     :     ,zz,pp,temp,qv,u,v,dw
-     :     ,dt,dq,sens,flat)
+      subroutine read_amma(fich_amma,nlevel,ntime                          &
+     &     ,zz,pp,temp,qv,u,v,dw                                           &
+     &     ,dt,dq,sens,flat)
 
 !program reading forcings of the AMMA case study
@@ -3156 +3140 @@
          end subroutine read_amma
 !======================================================================
-        SUBROUTINE interp_amma_time(day,day1,annee_ref
-     i         ,year_ini_amma,day_ini_amma,nt_amma,dt_amma,nlev_amma
-     i         ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma
-     o         ,vitw_prof,ht_prof,hq_prof,lat_prof,sens_prof)
+        SUBROUTINE interp_amma_time(day,day1,annee_ref                     &
+     &         ,year_ini_amma,day_ini_amma,nt_amma,dt_amma,nlev_amma       &
+     &         ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma               &
+     &         ,vitw_prof,ht_prof,hq_prof,lat_prof,sens_prof)
         implicit none
 
@@ -3237 +3221 @@
 
        if (it_amma1 .gt. nt_amma) then
-        write(*,*) 'PB-stop: day, it_amma1, it_amma2, timeit: '
-     :        ,day,day_ini_amma,it_amma1,it_amma2,timeit/86400.
+        write(*,*) 'PB-stop: day, it_amma1, it_amma2, timeit: '            &
+     &        ,day,day_ini_amma,it_amma1,it_amma2,timeit/86400.
         stop
        endif
@@ -3246 +3230 @@
        frac=max(frac,0.0)
 
-       lat_prof = lat_amma(it_amma2)
-     :          -frac*(lat_amma(it_amma2)-lat_amma(it_amma1))
-       sens_prof = sens_amma(it_amma2)
-     :          -frac*(sens_amma(it_amma2)-sens_amma(it_amma1))
+       lat_prof = lat_amma(it_amma2)                                       &
+     &          -frac*(lat_amma(it_amma2)-lat_amma(it_amma1)) 
+       sens_prof = sens_amma(it_amma2)                                     &
+     &          -frac*(sens_amma(it_amma2)-sens_amma(it_amma1))
 
        do k=1,nlev_amma
-        vitw_prof(k) = vitw_amma(k,it_amma2)
-     :          -frac*(vitw_amma(k,it_amma2)-vitw_amma(k,it_amma1))
-        ht_prof(k) = ht_amma(k,it_amma2)
-     :          -frac*(ht_amma(k,it_amma2)-ht_amma(k,it_amma1))
-        hq_prof(k) = hq_amma(k,it_amma2)
-     :          -frac*(hq_amma(k,it_amma2)-hq_amma(k,it_amma1))
+        vitw_prof(k) = vitw_amma(k,it_amma2)                               &
+     &          -frac*(vitw_amma(k,it_amma2)-vitw_amma(k,it_amma1))
+        ht_prof(k) = ht_amma(k,it_amma2)                                   &
+     &          -frac*(ht_amma(k,it_amma2)-ht_amma(k,it_amma1))
+        hq_prof(k) = hq_amma(k,it_amma2)                                   &
+     &          -frac*(hq_amma(k,it_amma2)-hq_amma(k,it_amma1))
         enddo
 
@@ -3264 +3248 @@
 
 !=====================================================================
-      subroutine read_fire(fich_fire,nlevel,ntime
-     :     ,zz,thl,qt,u,v,tke
-     :     ,ug,vg,wls,dqtdx,dqtdy,dqtdt,thl_rad)
+      subroutine read_fire(fich_fire,nlevel,ntime                          &
+     &     ,zz,thl,qt,u,v,tke                                              &
+     &     ,ug,vg,wls,dqtdx,dqtdy,dqtdt,thl_rad)
 
 !program reading forcings of the FIRE case study

LMDZ5/trunk/libf/phylmd/1D_decl_cases.h

@@ -157 +157 @@
         logical :: Turb_fcg_gcssold
 
-        common /turb_forcing/
-     s  dtime_frcg,hthturb_gcssold, hqturb_gcssold,Turb_fcg_gcssold
+        common /turb_forcing/                                                   &
+     &  dtime_frcg,hthturb_gcssold, hqturb_gcssold,Turb_fcg_gcssold
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! Declarations specifiques au cas Arm_cu

LMDZ5/trunk/libf/phylmd/1D_interp_cases.h

@@ -4 +4 @@
       if (forcing_GCSSold) then
 
-       call get_uvd(it,timestep,fich_gcssold_ctl,fich_gcssold_dat,
-     :               ht_gcssold,hq_gcssold,hw_gcssold,
-     :               hu_gcssold,hv_gcssold,
-     :               hthturb_gcssold,hqturb_gcssold,Ts_gcssold,
-     :               imp_fcg_gcssold,ts_fcg_gcssold,
-     :               Tp_fcg_gcssold,Turb_fcg_gcssold)
+       call get_uvd(it,timestep,fich_gcssold_ctl,fich_gcssold_dat,              &
+     &               ht_gcssold,hq_gcssold,hw_gcssold,                          &
+     &               hu_gcssold,hv_gcssold,                                     &
+     &               hthturb_gcssold,hqturb_gcssold,Ts_gcssold,                 &
+     &               imp_fcg_gcssold,ts_fcg_gcssold,                            &
+     &               Tp_fcg_gcssold,Turb_fcg_gcssold)
        if (prt_level.ge.1) then
          print *,' get_uvd -> hqturb_gcssold ',it,hqturb_gcssold
@@ -38 +38 @@
 
        if (prt_level.ge.1) then
-        print*,
-     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_toga=',
-     :    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_toga
+        print*,                                                             &
+     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_toga=',     &
+     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_toga
        endif
 
 ! time interpolation:
-        CALL interp_toga_time(daytime,day1,annee_ref
-     i             ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga
-     i             ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga
-     i             ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga
-     o             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof
-     o             ,ht_prof,vt_prof,hq_prof,vq_prof)
+        CALL interp_toga_time(daytime,day1,annee_ref                        &
+     &             ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga           &
+     &             ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga        &
+     &             ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga           &
+     &             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof    &
+     &             ,ht_prof,vt_prof,hq_prof,vq_prof)
 
         if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
 
 ! vertical interpolation:
-      CALL interp_toga_vertical(play,nlev_toga,plev_prof
-     :         ,t_prof,q_prof,u_prof,v_prof,w_prof
-     :         ,ht_prof,vt_prof,hq_prof,vq_prof
-     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
-     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
+      CALL interp_toga_vertical(play,nlev_toga,plev_prof                    &
+     &         ,t_prof,q_prof,u_prof,v_prof,w_prof                          &
+     &         ,ht_prof,vt_prof,hq_prof,vq_prof                             &
+     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                               &
+     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
 
 ! large-scale forcing :
@@ -84 +84 @@
       if (forcing_twpice) then
 
-        print*,
-     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_twpi=',
-     :    daytime,day1,(daytime-day1)*86400.,
-     :    (daytime-day1)*86400/dt_twpi
+        print*,                                                             &
+     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_twpi=',     &
+     &    daytime,day1,(daytime-day1)*86400.,                               &
+     &    (daytime-day1)*86400/dt_twpi
 
 ! time interpolation:
-        CALL interp_toga_time(daytime,day1,annee_ref
-     i       ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi
-     i       ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi
-     i       ,ht_twpi,vt_twpi,hq_twpi,vq_twpi
-     o       ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp,u_proftwp
-     o       ,v_proftwp,w_proftwp
-     o       ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp)
+        CALL interp_toga_time(daytime,day1,annee_ref                        &
+     &       ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi       &
+     &       ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi          &
+     &       ,ht_twpi,vt_twpi,hq_twpi,vq_twpi                               &
+     &       ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp,u_proftwp         &
+     &       ,v_proftwp,w_proftwp                                           &
+     &       ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp)
 
 ! vertical interpolation:
-      CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp
-     :         ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp
-     :         ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp
-     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
-     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
+      CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp                 &
+     &         ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp           &
+     &         ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp                 &
+     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                               &
+     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
 
 
 !calcul de l'advection verticale a partir du omega
-cCalcul des gradients verticaux
-cinitialisation
+!Calcul des gradients verticaux
+!initialisation
       d_t_z(:)=0.
       d_q_z(:)=0.
@@ -114 +114 @@
       d_q_dyn_z(:)=0.
       DO l=2,llm-1
-       d_t_z(l)=(temp(l+1)-temp(l-1))
-     &          /(play(l+1)-play(l-1))
-       d_q_z(l)=(q(l+1,1)-q(l-1,1))
-     &          /(play(l+1)-play(l-1))
+       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
+       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
       ENDDO
       d_t_z(1)=d_t_z(2)
@@ -124 +122 @@
       d_q_z(llm)=d_q_z(llm-1)
 
-cCalcul de l advection verticale
+!Calcul de l advection verticale
       d_t_dyn_z(:)=w_mod(:)*d_t_z(:)
       d_q_dyn_z(:)=w_mod(:)*d_q_z(:)
@@ -133 +131 @@
 !           if (phi(l).gt.5000.) then
         if (phi(l).gt.0.) then
-        u(l)=u(l)
-     .      +timestep*(u_mod(l)-u(l))/(2.*3600.)
-        v(l)=v(l) 
-     .      +timestep*(v_mod(l)-v(l))/(2.*3600.)
+        u(l)=u(l)+timestep*(u_mod(l)-u(l))/(2.*3600.)
+        v(l)=v(l)+timestep*(v_mod(l)-v(l))/(2.*3600.)
            endif    
         else
@@ -150 +146 @@
            if ((zz(l).le.16000.).and.(zz(l).gt.15000.)) then
              zfact=(zz(l)-15000.)/1000.
-        q(l,1)=q(l,1)
-     .      +timestep*(q_mod(l)-q(l,1))/(6.*3600.)*zfact
-        temp(l)=temp(l) 
-     .      +timestep*(t_mod(l)-temp(l))/(6.*3600.)*zfact
+        q(l,1)=q(l,1)+timestep*(q_mod(l)-q(l,1))/(6.*3600.)*zfact
+        temp(l)=temp(l)+timestep*(t_mod(l)-temp(l))/(6.*3600.)*zfact
            else if (zz(l).gt.16000.) then 
-        q(l,1)=q(l,1)
-     .      +timestep*(q_mod(l)-q(l,1))/(6.*3600.)
-        temp(l)=temp(l) 
-     .      +timestep*(t_mod(l)-temp(l))/(6.*3600.)
+        q(l,1)=q(l,1)+timestep*(q_mod(l)-q(l,1))/(6.*3600.)
+        temp(l)=temp(l)+timestep*(t_mod(l)-temp(l))/(6.*3600.)
            endif
         enddo    
@@ -189 +181 @@
        if (forcing_amma) then
 
-        print*,
-     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_amma=',
-     :    daytime,day1,(daytime-day1)*86400.,
-     :    (daytime-day1)*86400/dt_amma
+        print*,                                                             &
+     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_amma=',     &
+     &    daytime,day1,(daytime-day1)*86400.,                               &
+     &    (daytime-day1)*86400/dt_amma
 
 ! time interpolation using TOGA interpolation routine
-        CALL interp_amma_time(daytime,day1,annee_ref
-     i       ,year_ini_amma,day_ju_ini_amma,nt_amma,dt_amma,nlev_amma
-     i       ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma
-     o       ,vitw_profamma,ht_profamma,hq_profamma,lat_profamma
-     :       ,sens_profamma)
+        CALL interp_amma_time(daytime,day1,annee_ref                        &
+     &       ,year_ini_amma,day_ju_ini_amma,nt_amma,dt_amma,nlev_amma       &
+     &       ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma                  &
+     &       ,vitw_profamma,ht_profamma,hq_profamma,lat_profamma            &
+     &       ,sens_profamma)
 
       print*,'apres interpolation temporelle AMMA'
@@ -213 +205 @@
 ! vertical interpolation using TOGA interpolation routine:
 !      write(*,*)'avant interp vert', t_proftwp
-      CALL interp_toga_vertical(play,nlev_amma,plev_amma
-     :         ,th_profamma,q_profamma,u_profamma,v_profamma
-     :         ,vitw_profamma
-     :         ,ht_profamma,vt_profamma,hq_profamma,vq_profamma
-     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
-     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
+      CALL interp_toga_vertical(play,nlev_amma,plev_amma                      &
+     &         ,th_profamma,q_profamma,u_profamma,v_profamma                 &
+     &         ,vitw_profamma                                               &
+     &         ,ht_profamma,vt_profamma,hq_profamma,vq_profamma             &
+     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                               &
+     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
        write(*,*) 'Profil initial forcing AMMA interpole'
 
 
 !calcul de l'advection verticale a partir du omega
-cCalcul des gradients verticaux
-cinitialisation
+!Calcul des gradients verticaux
+!initialisation
       do l=1,llm
       d_t_z(l)=0.
@@ -231 +223 @@
 
       DO l=2,llm-1
-       d_t_z(l)=(temp(l+1)-temp(l-1))
-     &          /(play(l+1)-play(l-1))
-       d_q_z(l)=(q(l+1,1)-q(l-1,1))
-     &          /(play(l+1)-play(l-1))
+       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
+       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
       ENDDO
       d_t_z(1)=d_t_z(2)
@@ -250 +240 @@
 !        d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-omega(l)*d_t_z(l)
 !attention: on impose dth
-        d_th_adv(l) = alpha*omega(l)/rcpd+
+        d_th_adv(l) = alpha*omega(l)/rcpd+                                  &
      &         ht_mod(l)*(play(l)/pzero)**rkappa-omega(l)*d_t_z(l)
 !        d_th_adv(l) = 0.
@@ -275 +265 @@
 !       call lstendH(llm,omega,dt_dyn,dq_dyn,du_dyn, dv_dyn,
 !     :  q,temp,u,v,play)
-       call lstendH(llm,nqtot,omega,dt_dyn,dq_dyn,
-     :  q,temp,u,v,play)
+       call lstendH(llm,nqtot,omega,dt_dyn,dq_dyn,q,temp,u,v,play)
 
         do l=1,llm
@@ -290 +279 @@
       if (forcing_armcu) then
 
-        print*,
-     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_armcu=',
-     :    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_armcu
+        print*,                                                             &
+     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_armcu=',    &
+     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_armcu
 
 ! time interpolation:
 ! ATTENTION, cet appel ne convient pas pour TOGA !!
 ! revoir 1DUTILS.h et les arguments
-      CALL interp_armcu_time(daytime,day1,annee_ref
-     i            ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu
-     i            ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu
-     i            ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof
-     i            ,adv_theta_prof,rad_theta_prof,adv_qt_prof)
+      CALL interp_armcu_time(daytime,day1,annee_ref                         &
+     &            ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu        &
+     &            ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu          &
+     &            ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof         &
+     &            ,adv_theta_prof,rad_theta_prof,adv_qt_prof)
 
 ! vertical interpolation:
@@ -346 +335 @@
       if (forcing_sandu) then
 
-        print*,
-     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_sandu=',
-     :    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_sandu
+        print*,                                                             &
+     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_sandu=',    &
+     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_sandu 
 
 ! time interpolation:
 ! ATTENTION, cet appel ne convient pas pour TOGA !!
 ! revoir 1DUTILS.h et les arguments
-      CALL interp_sandu_time(daytime,day1,annee_ref
-     i             ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu
-     i             ,nlev_sandu
-     i             ,ts_sandu,ts_prof)
+      CALL interp_sandu_time(daytime,day1,annee_ref                         &
+     &             ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu       &
+     &             ,nlev_sandu                                              &
+     &             ,ts_sandu,ts_prof)
 
         if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
 
 ! vertical interpolation:
-      CALL interp_sandu_vertical(play,nlev_sandu,plev_profs
-     :         ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs
-     :         ,omega_profs,o3mmr_profs
-     :         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod
-     :         ,omega_mod,o3mmr_mod,mxcalc)
-!calcul de l'advection verticale
-cCalcul des gradients verticaux
-cinitialisation
-      d_t_z(:)=0.
-      d_q_z(:)=0.
-      d_t_dyn_z(:)=0.
-      d_q_dyn_z(:)=0.
-! schema centre
-!     DO l=2,llm-1
-!      d_t_z(l)=(temp(l+1)-temp(l-1))
-!    &          /(play(l+1)-play(l-1))
-!      d_q_z(l)=(q(l+1,1)-q(l-1,1))
-!    &          /(play(l+1)-play(l-1))
-! schema amont
-      DO l=2,llm-1
-       d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l))
-       d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l))
-!     print *,'l temp2 temp0 play2 play0 omega_mod',
-!    & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l)
-      ENDDO
-      d_t_z(1)=d_t_z(2)
-      d_q_z(1)=d_q_z(2)
-      d_t_z(llm)=d_t_z(llm-1)
-      d_q_z(llm)=d_q_z(llm-1)
-
-!  calcul de l advection verticale
-! Confusion w (m/s) et omega (Pa/s) !!
-      d_t_dyn_z(:)=omega_mod(:)*d_t_z(:)
-      d_q_dyn_z(:)=omega_mod(:)*d_q_z(:)
-!     do l=1,llm
-!      print *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z',
-!    :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l)
-!     enddo
-
-
-! large-scale forcing : pour le cas Sandu ces forcages sont la SST
-! et une divergence constante -> profil de omega
-      tsurf = ts_prof
-      write(*,*) 'SST suivante: ',tsurf
-      do l = 1, llm
-       omega(l) = omega_mod(l)
-       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
-
-       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
-!
-!      d_th_adv(l) = 0.0
-!      d_q_adv(l,1) = 0.0
-!CR:test advection=0
-!calcul de l'advection verticale
-        d_th_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
-!        print*,'temp adv',l,-d_t_dyn_z(l)
-        d_q_adv(l,1) = -d_q_dyn_z(l)
-!        print*,'q adv',l,-d_q_dyn_z(l)
-       dt_cooling(l) = 0.0
-      enddo
-      endif ! forcing_sandu
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!---------------------------------------------------------------------
-! Interpolation forcing in time and onto model levels
-!---------------------------------------------------------------------
-      if (forcing_astex) then
-
-        print*,
-     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_astex=',
-     :    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_astex
-
-! time interpolation:
-! ATTENTION, cet appel ne convient pas pour TOGA !!
-! revoir 1DUTILS.h et les arguments
-      CALL interp_astex_time(daytime,day1,annee_ref
-     i             ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex
-     i             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex
-     i             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof
-     i             ,ufa_prof,vfa_prof)
-
-        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
-
-! vertical interpolation:
-      CALL interp_astex_vertical(play,nlev_astex,plev_profa
-     :         ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa
-     :         ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa
-     :         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod
-     :         ,tke_mod,o3mmr_mod,mxcalc)
+      CALL interp_sandu_vertical(play,nlev_sandu,plev_profs                 &
+     &         ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs           &
+     &         ,omega_profs,o3mmr_profs                                     &
+     &         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod                       &
+     &         ,omega_mod,o3mmr_mod,mxcalc)
 !calcul de l'advection verticale
 !Calcul des gradients verticaux
@@ -476 +382 @@
 !  calcul de l advection verticale
 ! Confusion w (m/s) et omega (Pa/s) !!
-      d_t_dyn_z(:)=w_mod(:)*d_t_z(:)
-      d_q_dyn_z(:)=w_mod(:)*d_q_z(:)
+      d_t_dyn_z(:)=omega_mod(:)*d_t_z(:)
+      d_q_dyn_z(:)=omega_mod(:)*d_q_z(:)
 !     do l=1,llm
 !      print *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z',
@@ -484 +390 @@
 
 
-! large-scale forcing : pour le cas Astex ces forcages sont la SST
-! la divergence,ug,vg,ufa,vfa
+! large-scale forcing : pour le cas Sandu ces forcages sont la SST
+! et une divergence constante -> profil de omega
       tsurf = ts_prof
       write(*,*) 'SST suivante: ',tsurf
       do l = 1, llm
-       omega(l) = w_mod(l)
+       omega(l) = omega_mod(l)
        omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
 
@@ -504 +410 @@
        dt_cooling(l) = 0.0
       enddo
+      endif ! forcing_sandu
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!---------------------------------------------------------------------
+! Interpolation forcing in time and onto model levels
+!---------------------------------------------------------------------
+      if (forcing_astex) then
+
+        print*,                                                             &
+     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_astex=',    &
+     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_astex
+
+! time interpolation:
+! ATTENTION, cet appel ne convient pas pour TOGA !!
+! revoir 1DUTILS.h et les arguments
+      CALL interp_astex_time(daytime,day1,annee_ref                         &
+     &             ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex       &
+     &             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex         &
+     &             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof    &
+     &             ,ufa_prof,vfa_prof)
+
+        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
+
+! vertical interpolation:
+      CALL interp_astex_vertical(play,nlev_astex,plev_profa                 &
+     &         ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa                &
+     &         ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa               &
+     &         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod        &
+     &         ,tke_mod,o3mmr_mod,mxcalc)
+!calcul de l'advection verticale
+!Calcul des gradients verticaux
+!initialisation
+      d_t_z(:)=0.
+      d_q_z(:)=0.
+      d_t_dyn_z(:)=0.
+      d_q_dyn_z(:)=0.
+! schema centre
+!     DO l=2,llm-1
+!      d_t_z(l)=(temp(l+1)-temp(l-1))
+!    &          /(play(l+1)-play(l-1))
+!      d_q_z(l)=(q(l+1,1)-q(l-1,1))
+!    &          /(play(l+1)-play(l-1))
+! schema amont
+      DO l=2,llm-1
+       d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l))
+       d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l))
+!     print *,'l temp2 temp0 play2 play0 omega_mod',
+!    & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l)
+      ENDDO
+      d_t_z(1)=d_t_z(2)
+      d_q_z(1)=d_q_z(2)
+      d_t_z(llm)=d_t_z(llm-1)
+      d_q_z(llm)=d_q_z(llm-1)
+
+!  calcul de l advection verticale
+! Confusion w (m/s) et omega (Pa/s) !!
+      d_t_dyn_z(:)=w_mod(:)*d_t_z(:)
+      d_q_dyn_z(:)=w_mod(:)*d_q_z(:)
+!     do l=1,llm
+!      print *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z',
+!    :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l)
+!     enddo
+
+
+! large-scale forcing : pour le cas Astex ces forcages sont la SST
+! la divergence,ug,vg,ufa,vfa
+      tsurf = ts_prof
+      write(*,*) 'SST suivante: ',tsurf
+      do l = 1, llm
+       omega(l) = w_mod(l)
+       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
+
+       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
+!
+!      d_th_adv(l) = 0.0
+!      d_q_adv(l,1) = 0.0
+!CR:test advection=0
+!calcul de l'advection verticale
+        d_th_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
+!        print*,'temp adv',l,-d_t_dyn_z(l)
+        d_q_adv(l,1) = -d_q_dyn_z(l)
+!        print*,'q adv',l,-d_q_dyn_z(l)
+       dt_cooling(l) = 0.0
+      enddo
       endif ! forcing_astex
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

LMDZ5/trunk/libf/phylmd/1D_nudge_sandu_astex.h

@@ -20 +20 @@
 !     print *,'l dt  relax dtadv',l,dt_phys(l),relax_thl(l),d_th_adv(l)
       enddo
-        u(1:mxcalc)=u(1:mxcalc) + timestep*(
-     .              du_phys(1:mxcalc) - relax_u(1:mxcalc))
-        v(1:mxcalc)=v(1:mxcalc) + timestep*(
-     .              dv_phys(1:mxcalc) - relax_v(1:mxcalc))
-c       q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(
-c    .              dq(1:mxcalc,:) - relax_q(1:mxcalc,:)+
-c    .              d_q_adv(1:mxcalc,:))
-        q(1:mxcalc,1)=q(1:mxcalc,1)+timestep*(
-     .      dq(1:mxcalc,1) - relax_q(1:mxcalc,1)+d_q_adv(1:mxcalc,1))
-        q(1:mxcalc,2)=q(1:mxcalc,2)+timestep*(
-     .      dq(1:mxcalc,2) - relax_q(1:mxcalc,2)+d_q_adv(1:mxcalc,2))
-        temp(1:mxcalc)=temp(1:mxcalc)+timestep*(
-     .      dt_phys(1:mxcalc)-relax_thl(1:mxcalc)+d_th_adv(1:mxcalc))
+        u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &          
+     &              du_phys(1:mxcalc) - relax_u(1:mxcalc))
+        v(1:mxcalc)=v(1:mxcalc) + timestep*(                                &
+     &               dv_phys(1:mxcalc) - relax_v(1:mxcalc))
+!       q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(
+!    .              dq(1:mxcalc,:) - relax_q(1:mxcalc,:)+
+!    .              d_q_adv(1:mxcalc,:))
+        q(1:mxcalc,1)=q(1:mxcalc,1)+timestep*(                              &
+     &      dq(1:mxcalc,1) - relax_q(1:mxcalc,1)+d_q_adv(1:mxcalc,1))
+        q(1:mxcalc,2)=q(1:mxcalc,2)+timestep*(                              &
+     &      dq(1:mxcalc,2) - relax_q(1:mxcalc,2)+d_q_adv(1:mxcalc,2))
+        temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
+     &      dt_phys(1:mxcalc)-relax_thl(1:mxcalc)+d_th_adv(1:mxcalc))

LMDZ5/trunk/libf/phylmd/1D_read_forc_cases.h

@@ -8 +8 @@
       nq2=0
 
-      if (forcing_les .or. forcing_radconv
-     :    .or. forcing_GCSSold .or. forcing_fire) then
+      if (forcing_les .or. forcing_radconv                                      &
+     &    .or. forcing_GCSSold .or. forcing_fire) then
 
       if (forcing_fire) then
@@ -16 +16 @@
 !----------------------------------------------------------------------
       fich_fire='fire.nc'
-      call read_fire(fich_fire,nlev_fire,nt_fire
-     :     ,height,tttprof,qtprof,uprof,vprof,e12prof
-     :     ,ugprof,vgprof,wfls,dqtdxls
-     :     ,dqtdyls,dqtdtls,thlpcar)
+      call read_fire(fich_fire,nlev_fire,nt_fire                                &
+     &     ,height,tttprof,qtprof,uprof,vprof,e12prof                           &
+     &     ,ugprof,vgprof,wfls,dqtdxls                                          &
+     &     ,dqtdyls,dqtdtls,thlpcar)
       write(*,*) 'Forcing FIRE lu'
       kmax=120            ! nombre de niveaux dans les profils et forcages
@@ -28 +28 @@
 !----------------------------------------------------------------------
 
-      call readprofiles(nlev_max,kmax,nqtot,height,
-     .           tttprof,qtprof,uprof,vprof,
-     .           e12prof,ugprof,vgprof,
-     .           wfls,dqtdxls,dqtdyls,dqtdtls,
-     .           thlpcar,qprof,nq1,nq2)
+      call readprofiles(nlev_max,kmax,nqtot,height,                             &
+     &           tttprof,qtprof,uprof,vprof,                                    &
+     &           e12prof,ugprof,vgprof,                                         &
+     &           wfls,dqtdxls,dqtdyls,dqtdtls,                                  &
+     &           thlpcar,qprof,nq1,nq2)
       endif
 
@@ -66 +66 @@
         ug(l)   = ugprof(kmax)-frac*( ugprof(kmax)- ugprof(kmax-1))
         vg(l)   = vgprof(kmax)-frac*( vgprof(kmax)- vgprof(kmax-1))
-        IF (nq2>0) q(l,nq1:nq2)=qprof(kmax,nq1:nq2)
-     s               -frac*(qprof(kmax,nq1:nq2)-qprof(kmax-1,nq1:nq2))
+        IF (nq2>0) q(l,nq1:nq2)=qprof(kmax,nq1:nq2)                         &
+     &               -frac*(qprof(kmax,nq1:nq2)-qprof(kmax-1,nq1:nq2))
         omega(l)=   wfls(kmax)-frac*(   wfls(kmax)-   wfls(kmax-1))
 
         dq_dyn(l,1) = dqtdtls(kmax)-frac*(dqtdtls(kmax)-dqtdtls(kmax-1))
-        dt_cooling(l)
-     .          =thlpcar(kmax)-frac*(thlpcar(kmax)-thlpcar(kmax-1))
+        dt_cooling(l)=thlpcar(kmax)-frac*(thlpcar(kmax)-thlpcar(kmax-1))
         do k=2,kmax
           frac = (height(k)-zlay(l))/(height(k)-height(k-1))
@@ -91 +90 @@
             ug(l)   = ugprof(k)-frac*( ugprof(k)- ugprof(k-1))
             vg(l)   = vgprof(k)-frac*( vgprof(k)- vgprof(k-1))
-            IF (nq2>0) q(l,nq1:nq2)=qprof(k,nq1:nq2)
-     s                   -frac*(qprof(k,nq1:nq2)-qprof(k-1,nq1:nq2))
+            IF (nq2>0) q(l,nq1:nq2)=qprof(k,nq1:nq2)                        &
+     &                   -frac*(qprof(k,nq1:nq2)-qprof(k-1,nq1:nq2))
             omega(l)=   wfls(k)-frac*(   wfls(k)-   wfls(k-1))
             dq_dyn(l,1)=dqtdtls(k)-frac*(dqtdtls(k)-dqtdtls(k-1))
-            dt_cooling(l)
-     .              =thlpcar(k)-frac*(thlpcar(k)-thlpcar(k-1))
+            dt_cooling(l)=thlpcar(k)-frac*(thlpcar(k)-thlpcar(k-1))
           elseif(zlay(l)<height(1)) then ! profils uniformes pour z<height(1)
             ttt =tttprof(1)
@@ -144 +142 @@
        fich_gcssold_dat = './forcing8.dat'
        call copie(llm,play,psurf,fich_gcssold_ctl)
-       call get_uvd2(it,timestep,fich_gcssold_ctl,fich_gcssold_dat,
-     :               ht_gcssold,hq_gcssold,hw_gcssold,
-     :               hu_gcssold,hv_gcssold,
-     :               hthturb_gcssold,hqturb_gcssold,Ts_gcssold,
-     :               imp_fcg_gcssold,ts_fcg_gcssold,
-     :               Tp_fcg_gcssold,Turb_fcg_gcssold)
+       call get_uvd2(it,timestep,fich_gcssold_ctl,fich_gcssold_dat,         &
+     &               ht_gcssold,hq_gcssold,hw_gcssold,                      &
+     &               hu_gcssold,hv_gcssold,                                 &
+     &               hthturb_gcssold,hqturb_gcssold,Ts_gcssold,             &
+     &               imp_fcg_gcssold,ts_fcg_gcssold,                        &
+     &               Tp_fcg_gcssold,Turb_fcg_gcssold)
        print *,' get_uvd2 -> hqturb_gcssold ',hqturb_gcssold
       endif ! forcing_GCSSold
@@ -160 +158 @@
 !       call writefield_phy('omega', omega,llm+1)
       fich_rico = 'rico.txt'
-       call read_rico(fich_rico,nlev_rico,ps_rico,play
-     :             ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico
-     :             ,dth_rico,dqh_rico)
+       call read_rico(fich_rico,nlev_rico,ps_rico,play                      &
+     &             ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico              &
+     &             ,dth_rico,dqh_rico)
         print*, ' on a lu et prepare RICO'
 
@@ -189 +187 @@
 ! read TOGA-COARE forcing (native vertical grid, nt_toga timesteps):
       fich_toga = './d_toga/ifa_toga_coare_v21_dime.txt'
-      CALL read_togacoare(fich_toga,nlev_toga,nt_toga
-     :         ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga
-     :         ,ht_toga,vt_toga,hq_toga,vq_toga)
+      CALL read_togacoare(fich_toga,nlev_toga,nt_toga                       &
+     &         ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga        &
+     &         ,ht_toga,vt_toga,hq_toga,vq_toga)
 
        write(*,*) 'Forcing TOGA lu'
@@ -197 +195 @@
 ! time interpolation for initial conditions:
       write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
-      CALL interp_toga_time(daytime,day1,annee_ref
-     i             ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga
-     i             ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga
-     i             ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga
-     o             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof
-     o             ,ht_prof,vt_prof,hq_prof,vq_prof)
+      CALL interp_toga_time(daytime,day1,annee_ref                          &
+     &             ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga           &
+     &             ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga        &
+     &             ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga           &
+     &             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof    &
+     &             ,ht_prof,vt_prof,hq_prof,vq_prof)
 
 ! vertical interpolation:
-      CALL interp_toga_vertical(play,nlev_toga,plev_prof
-     :         ,t_prof,q_prof,u_prof,v_prof,w_prof
-     :         ,ht_prof,vt_prof,hq_prof,vq_prof
-     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
-     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
+      CALL interp_toga_vertical(play,nlev_toga,plev_prof                    &
+     &         ,t_prof,q_prof,u_prof,v_prof,w_prof                          &
+     &         ,ht_prof,vt_prof,hq_prof,vq_prof                             &
+     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                               &
+     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
        write(*,*) 'Profil initial forcing TOGA interpole'
 
@@ -239 +237 @@
       if (forcing_twpice) then
 !read TWP-ICE forcings
-      fich_twpice=
-     : 'd_twpi/twp180iopsndgvarana_v2.1_C3.c1.20060117.000000.cdf'
-      call read_twpice(fich_twpice,nlev_twpi,nt_twpi
-     :     ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi
-     :     ,ht_twpi,vt_twpi,hq_twpi,vq_twpi)
+     fich_twpice='d_twpi/twp180iopsndgvarana_v2.1_C3.c1.20060117.000000.cdf'
+      call read_twpice(fich_twpice,nlev_twpi,nt_twpi                        &
+     &     ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi            &
+     &     ,ht_twpi,vt_twpi,hq_twpi,vq_twpi)
 
       write(*,*) 'Forcing TWP-ICE lu'
 !Time interpolation for initial conditions using TOGA interpolation routine
          write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1   
-      CALL interp_toga_time(daytime,day1,annee_ref
-     i          ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi
-     i             ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi
-     i             ,ht_twpi,vt_twpi,hq_twpi,vq_twpi
-     o             ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp
-     o             ,u_proftwp,v_proftwp,w_proftwp
-     o             ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp)
+      CALL interp_toga_time(daytime,day1,annee_ref                          &
+     &          ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi    &
+     &             ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi    &
+     &             ,ht_twpi,vt_twpi,hq_twpi,vq_twpi                         &
+     &             ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp             &
+     &             ,u_proftwp,v_proftwp,w_proftwp                           &
+     &             ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp)
 
 ! vertical interpolation using TOGA interpolation routine:
 !      write(*,*)'avant interp vert', t_proftwp
-      CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp
-     :         ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp
-     :         ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp
-     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
-     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
+      CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp                 &
+     &         ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp           &
+     &         ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp                 &
+     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                               &
+     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
 !       write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod
 
@@ -295 +292 @@
 !read AMMA forcings
       fich_amma='amma.nc'
-      call read_amma(fich_amma,nlev_amma,nt_amma
-     :     ,z_amma,plev_amma,th_amma,q_amma,u_amma,v_amma,vitw_amma
-     :     ,ht_amma,hq_amma,sens_amma,lat_amma)
+      call read_amma(fich_amma,nlev_amma,nt_amma                            &
+     &     ,z_amma,plev_amma,th_amma,q_amma,u_amma,v_amma,vitw_amma         &
+     &     ,ht_amma,hq_amma,sens_amma,lat_amma)
 
       write(*,*) 'Forcing AMMA lu'
@@ -318 +315 @@
 ! vertical interpolation using TOGA interpolation routine:
 !      write(*,*)'avant interp vert', t_proftwp
-      CALL interp_toga_vertical(play,nlev_amma,plev_amma
-     :         ,th_ammai,q_ammai,u_ammai,v_ammai,vitw_ammai
-     :         ,ht_ammai,vt_ammai,hq_ammai,vq_ammai
-     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
-     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
+      CALL interp_toga_vertical(play,nlev_amma,plev_amma                    &
+     &         ,th_ammai,q_ammai,u_ammai,v_ammai,vitw_ammai                 &
+     &         ,ht_ammai,vt_ammai,hq_ammai,vq_ammai                         &
+     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                               &
+     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
 !       write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod
 
@@ -351 +348 @@
        dt_cooling(l)=0.
       enddo     
-       write(*,*) 'Profil initial forcing AMMA interpole temp39',
-     &           temp(39)
+       write(*,*) 'Prof initeforcing AMMA interpole temp39',temp(39)
       
 
@@ -374 +370 @@
        write(*,*) 'Avant lecture Forcing Arm_Cu'
       fich_armcu = './ifa_armcu.txt'
-      CALL read_armcu(fich_armcu,nlev_armcu,nt_armcu,
-     : sens_armcu,flat_armcu,adv_theta_armcu,
-     : rad_theta_armcu,adv_qt_armcu)
+      CALL read_armcu(fich_armcu,nlev_armcu,nt_armcu,                       &
+     & sens_armcu,flat_armcu,adv_theta_armcu,                               &
+     & rad_theta_armcu,adv_qt_armcu)
        write(*,*) 'Forcing Arm_Cu lu'
 
@@ -391 +387 @@
 !----------------------------------------------------------------------
 
-      call readprofile_armcu(nlev_max,kmax,height,play_mod,u_mod,
-     .           v_mod,theta_mod,t_mod,qv_mod,rv_mod,ap,bp)
+      call readprofile_armcu(nlev_max,kmax,height,play_mod,u_mod,           &
+     &           v_mod,theta_mod,t_mod,qv_mod,rv_mod,ap,bp)
 
 ! time interpolation for initial conditions:
@@ -406 +402 @@
       print *,'dt_armcu=',dt_armcu
       print *,'nlev_armcu=',nlev_armcu
-      CALL interp_armcu_time(daytime,day1,annee_ref
-     i            ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu
-     i            ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu
-     i            ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof
-     i            ,adv_theta_prof,rad_theta_prof,adv_qt_prof)
+      CALL interp_armcu_time(daytime,day1,annee_ref                         &
+     &            ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu        &
+     &            ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu         &
+     &            ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof         &
+     &            ,adv_theta_prof,rad_theta_prof,adv_qt_prof)
        write(*,*) 'Forcages interpoles dans temps'
 
@@ -456 +452 @@
       do l = 1, llm
       plev(l+1) = ap(llm-l+1)+bp(llm-l+1)*psurf
-      print *,'Read_forc: l height play plev zlay temp',
-     :   l,height(l),play(l),plev(l),zlay(l),temp(l)
+      print *,'Read_forc: l height play plev zlay temp',                    &
+     &   l,height(l),play(l),plev(l),zlay(l),temp(l)
       enddo
 ! For this case, fluxes are imposed
@@ -482 +478 @@
 !----------------------------------------------------------------------
 
-      call readprofile_sandu(nlev_max,kmax,height,plev_profs,t_profs,
-     .           thl_profs,q_profs,u_profs,v_profs,
-     .           w_profs,omega_profs,o3mmr_profs)
+      call readprofile_sandu(nlev_max,kmax,height,plev_profs,t_profs,       &
+     &           thl_profs,q_profs,u_profs,v_profs,                         &
+     &           w_profs,omega_profs,o3mmr_profs)
 
 ! time interpolation for initial conditions:
@@ -500 +496 @@
       print *,'dt_sandu=',dt_sandu
       print *,'nlev_sandu=',nlev_sandu
-      CALL interp_sandu_time(daytime,day1,annee_ref
-     i             ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu
-     i             ,nlev_sandu
-     i             ,ts_sandu,ts_prof)
+      CALL interp_sandu_time(daytime,day1,annee_ref                         &
+     &             ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu       &
+     &             ,nlev_sandu                                              &
+     &             ,ts_sandu,ts_prof)
 
 ! vertical interpolation:
       print *,'Avant interp_vertical: nlev_sandu=',nlev_sandu
-      CALL interp_sandu_vertical(play,nlev_sandu,plev_profs
-     :         ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs
-     :         ,omega_profs,o3mmr_profs
-     :         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod
-     :         ,omega_mod,o3mmr_mod,mxcalc)
+      CALL interp_sandu_vertical(play,nlev_sandu,plev_profs                 &
+     &         ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs           &
+     &         ,omega_profs,o3mmr_profs                                     &
+     &         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod                       &
+     &         ,omega_mod,o3mmr_mod,mxcalc)
        write(*,*) 'Profil initial forcing SANDU interpole'
 
@@ -548 +544 @@
 ! read astex forcing :
       fich_astex = './ifa_astex.txt'
-      CALL read_astex(fich_astex,nlev_astex,nt_astex,div_astex,ts_astex,
-     :  ug_astex,vg_astex,ufa_astex,vfa_astex)
+      CALL read_astex(fich_astex,nlev_astex,nt_astex,div_astex,ts_astex,    &
+     &  ug_astex,vg_astex,ufa_astex,vfa_astex)
 
        write(*,*) 'Forcing Astex lu'
@@ -557 +553 @@
 !----------------------------------------------------------------------
 
-      call readprofile_astex(nlev_max,kmax,height,plev_profa,t_profa,
-     .           thl_profa,qv_profa,ql_profa,qt_profa,u_profa,v_profa,
-     .           w_profa,tke_profa,o3mmr_profa)
+      call readprofile_astex(nlev_max,kmax,height,plev_profa,t_profa,       &
+     &           thl_profa,qv_profa,ql_profa,qt_profa,u_profa,v_profa,      &
+     &           w_profa,tke_profa,o3mmr_profa)
 
 ! time interpolation for initial conditions:
@@ -575 +571 @@
       print *,'dt_astex=',dt_astex
       print *,'nlev_astex=',nlev_astex
-      CALL interp_astex_time(daytime,day1,annee_ref
-     i             ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex
-     i             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex
-     i             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof
-     i             ,ufa_prof,vfa_prof)
+      CALL interp_astex_time(daytime,day1,annee_ref                         &
+     &             ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex       &
+     &             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex         &
+     &             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof    &
+     &             ,ufa_prof,vfa_prof)
 
 ! vertical interpolation:
       print *,'Avant interp_vertical: nlev_astex=',nlev_astex
-      CALL interp_astex_vertical(play,nlev_astex,plev_profa
-     :         ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa
-     :         ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa
-     :         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod
-     :         ,tke_mod,o3mmr_mod,mxcalc)
+      CALL interp_astex_vertical(play,nlev_astex,plev_profa                 &
+     &         ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa                &
+     &         ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa               &
+     &         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod        &
+     &         ,tke_mod,o3mmr_mod,mxcalc)
        write(*,*) 'Profil initial forcing Astex interpole'
 

LMDZ5/trunk/libf/phylmd/1Dconv.h

@@ -1 @@
-        subroutine get_uvd(itap,dtime,file_forctl,file_fordat,
-     :       ht,hq,hw,hu,hv,hthturb,hqturb,
-     :       Ts,imp_fcg,ts_fcg,Tp_fcg,Turb_fcg)
-c
+        subroutine get_uvd(itap,dtime,file_forctl,file_fordat,                  &
+     &       ht,hq,hw,hu,hv,hthturb,hqturb,                                     &
+     &       Ts,imp_fcg,ts_fcg,Tp_fcg,Turb_fcg)                                 
+!
         implicit none
  
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c cette routine permet d obtenir u_convg,v_convg,ht,hq et ainsi de
-c pouvoir calculer la convergence et le cisaillement dans la physiq
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! cette routine permet d obtenir u_convg,v_convg,ht,hq et ainsi de
+! pouvoir calculer la convergence et le cisaillement dans la physiq
+!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
 #include "YOMCST.h"
@@ -29 +29 @@
       COMMON/com2_phys_gcss/playm,hplaym,nblvlm
 
-c======================================================================
-c methode: on va chercher les donnees du mesoNH de meteo france, on y
-c          a acces a tout pas detemps grace a la routine rdgrads qui
-c          est une boucle lisant dans ces fichiers.
-c          Puis on interpole ces donnes sur les 11 niveaux du gcm et
-c          et sur les pas de temps de ce meme gcm
-c----------------------------------------------------------------------
-c input:
-c       pasmax     :nombre de pas de temps maximum du mesoNH
-c       dt         :pas de temps du meso_NH (en secondes)
-c----------------------------------------------------------------------
+!======================================================================
+! methode: on va chercher les donnees du mesoNH de meteo france, on y
+!          a acces a tout pas detemps grace a la routine rdgrads qui
+!          est une boucle lisant dans ces fichiers.
+!          Puis on interpole ces donnes sur les 11 niveaux du gcm et
+!          et sur les pas de temps de ce meme gcm
+!----------------------------------------------------------------------
+! input:
+!       pasmax     :nombre de pas de temps maximum du mesoNH
+!       dt         :pas de temps du meso_NH (en secondes)
+!----------------------------------------------------------------------
       integer pasmax,dt
       save pasmax,dt
-c----------------------------------------------------------------------
-c arguments:
-c           itap   :compteur de la physique(le nombre de ces pas est
-c                   fixe dans la subroutine calcul_ini_gcm de interpo
-c                   -lation
-c           dtime  :pas detemps du gcm (en secondes)
-c           ht     :convergence horizontale de temperature(K/s)
-c           hq     :    "         "       d humidite (kg/kg/s)
-c           hw     :vitesse verticale moyenne (m/s**2)
-c           hu     :convergence horizontale d impulsion le long de x
-c                  (kg/(m^2 s^2)
-c           hv     : idem le long de y.
-c           Ts     : Temperature de surface (K)
-c           imp_fcg: var. logical .eq. T si forcage en impulsion
-c           ts_fcg: var. logical .eq. T si forcage en Ts present dans fichier
-c           Tp_fcg: var. logical .eq. T si forcage donne en Temp potentielle
-c           Turb_fcg: var. logical .eq. T si forcage turbulent present dans fichier
-c----------------------------------------------------------------------
+!----------------------------------------------------------------------
+! arguments:
+!           itap   :compteur de la physique(le nombre de ces pas est
+!                   fixe dans la subroutine calcul_ini_gcm de interpo
+!                   -lation
+!           dtime  :pas detemps du gcm (en secondes)
+!           ht     :convergence horizontale de temperature(K/s)
+!           hq     :    "         "       d humidite (kg/kg/s)
+!           hw     :vitesse verticale moyenne (m/s**2)
+!           hu     :convergence horizontale d impulsion le long de x
+!                  (kg/(m^2 s^2)
+!           hv     : idem le long de y.
+!           Ts     : Temperature de surface (K)
+!           imp_fcg: var. logical .eq. T si forcage en impulsion
+!           ts_fcg: var. logical .eq. T si forcage en Ts present dans fichier
+!           Tp_fcg: var. logical .eq. T si forcage donne en Temp potentielle
+!           Turb_fcg: var. logical .eq. T si forcage turbulent present dans fichier
+!----------------------------------------------------------------------
         integer itap
         real dtime
@@ -74 +74 @@
         logical Tp_fcg
         logical Turb_fcg
-c----------------------------------------------------------------------
-c Variables internes de get_uvd (note : l interpolation temporelle
-c est faite entre les pas de temps before et after, sur les variables
-c definies sur la grille du SCM; on atteint exactement les valeurs Meso
-c aux milieux des pas de temps Meso)
-c     time0     :date initiale en secondes
-c     time      :temps associe a chaque pas du SCM
-c     pas       :numero du pas du meso_NH (on lit en pas : le premier pas
-c                 des donnees est duplique)
-c     pasprev   :numero du pas de lecture precedent
-c     htaft     :advection horizontale de temp. au pas de temps after
-c     hqaft     :    "         "      d humidite        "
-c     hwaft     :vitesse verticalle moyenne  au pas de temps after
-c     huaft,hvaft :advection horizontale d impulsion au pas de temps after
-c     tsaft     : surface temperature 'after time step'
-c     htbef     :idem htaft, mais pour le pas de temps before
-c     hqbef     :voir hqaft
-c     hwbef     :voir hwaft
-c     hubef,hvbef : idem huaft,hvaft, mais pour before
-c     tsbef     : surface temperature 'before time step'
-c----------------------------------------------------------------------
+!----------------------------------------------------------------------
+! Variables internes de get_uvd (note : l interpolation temporelle
+! est faite entre les pas de temps before et after, sur les variables
+! definies sur la grille du SCM; on atteint exactement les valeurs Meso
+! aux milieux des pas de temps Meso)
+!     time0     :date initiale en secondes
+!     time      :temps associe a chaque pas du SCM
+!     pas       :numero du pas du meso_NH (on lit en pas : le premier pas
+!                 des donnees est duplique)
+!     pasprev   :numero du pas de lecture precedent
+!     htaft     :advection horizontale de temp. au pas de temps after
+!     hqaft     :    "         "      d humidite        "
+!     hwaft     :vitesse verticalle moyenne  au pas de temps after
+!     huaft,hvaft :advection horizontale d impulsion au pas de temps after
+!     tsaft     : surface temperature 'after time step'
+!     htbef     :idem htaft, mais pour le pas de temps before
+!     hqbef     :voir hqaft
+!     hwbef     :voir hwaft
+!     hubef,hvbef : idem huaft,hvaft, mais pour before
+!     tsbef     : surface temperature 'before time step'
+!----------------------------------------------------------------------
         integer time0,pas,pasprev
         save time0,pas,pasprev
@@ -106 +106 @@
         real Tsbef
         save htbef,hqbef,hwbef,hubef,hvbef,hthturbbef,hqturbbef
-c
+!
         real timeaft,timebef
         save timeaft,timebef
         integer temps
         character*4 string
-c----------------------------------------------------------------------
-c variables arguments de la subroutine rdgrads
-c---------------------------------------------------------------------
+!----------------------------------------------------------------------
+! variables arguments de la subroutine rdgrads
+!---------------------------------------------------------------------
         integer icompt,icomp1 !compteurs de rdgrads
         real z(100)         ! altitude (grille Meso)
@@ -128 +128 @@
         real hqturb_mes(100) !tendance horizontale d humidite, due aux
                               !flux turbulents
-c
-c---------------------------------------------------------------------
-c variable argument de la subroutine copie
-c---------------------------------------------------------------------
-c SB        real pplay(100)    !pression en milieu de couche du gcm
-c SB                            !argument de la physique
-c---------------------------------------------------------------------
-c variables destinees a la lecture du pas de temps du fichier de donnees
-c---------------------------------------------------------------------
+!
+!---------------------------------------------------------------------
+! variable argument de la subroutine copie
+!---------------------------------------------------------------------
+! SB        real pplay(100)    !pression en milieu de couche du gcm
+! SB                            !argument de la physique
+!---------------------------------------------------------------------
+! variables destinees a la lecture du pas de temps du fichier de donnees
+!---------------------------------------------------------------------
        character*80 aaa,atemps,spaces,apasmax
        integer nch,imn,ipa
-c---------------------------------------------------------------------
-c  procedures appelees
+!---------------------------------------------------------------------
+!  procedures appelees
         external rdgrads    !lire en iterant dans forcing.dat
-c---------------------------------------------------------------------
+!---------------------------------------------------------------------
                print*,'le pas itap est:',itap
-c*** on determine le pas du meso_NH correspondant au nouvel itap ***
-c*** pour aller chercher les champs dans rdgrads                 ***
-c
+!*** on determine le pas du meso_NH correspondant au nouvel itap ***
+!*** pour aller chercher les champs dans rdgrads                 ***
+!
         time=time0+itap*dtime
-cc        temps=int(time/dt+1)
-cc        pas=min(temps,pasmax)
+!c        temps=int(time/dt+1)
+!c        pas=min(temps,pasmax)
         temps = 1 + int((dt + 2*time)/(2*dt))
         pas=min(temps,pasmax-1)
              print*,'le pas Meso est:',pas
-c
-c
-c===================================================================
-c
-c*** on remplit les champs before avec les champs after du pas   ***
-c*** precedent en format gcm                                     ***
+!
+!
+!===================================================================
+!
+!*** on remplit les champs before avec les champs after du pas   ***
+!*** precedent en format gcm                                     ***
         if(pas.gt.pasprev)then
           do i=1,klev
@@ -180 +180 @@
          print *, imp_fcg,ts_fcg,Turb_fcg,pas,nblvlm,icompt
                        print*,'le pas pas est:',pas
-c*** on va chercher les nouveaux champs after dans toga.dat     ***
-c*** champs en format meso_NH                                   ***
-          open(99,FILE=file_fordat,FORM='UNFORMATTED',
-     .             ACCESS='DIRECT',RECL=8)
-          call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes
-     .                  ,hu_mes,hv_mes,hthturb_mes,hqturb_mes
-     .                  ,ts_fcg,ts_subr,imp_fcg,Turb_fcg)
-c
+!*** on va chercher les nouveaux champs after dans toga.dat     ***
+!*** champs en format meso_NH                                   ***
+          open(99,FILE=file_fordat,FORM='UNFORMATTED',                        &
+     &             ACCESS='DIRECT',RECL=8)
+          call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes                &
+     &                  ,hu_mes,hv_mes,hthturb_mes,hqturb_mes                 &
+     &                  ,ts_fcg,ts_subr,imp_fcg,Turb_fcg)
+!
 
                if(Tp_fcg) then
-c     (le forcage est donne en temperature potentielle)
+!     (le forcage est donne en temperature potentielle)
          do i = 1,nblvlm
            ht_mes(i) = ht_mes(i)*(hplaym(i)/1000.)**rkappa
@@ -200 +200 @@
          enddo
         endif  ! Turb_fcg
-c
+!
                print*,'ht_mes ',(ht_mes(i),i=1,nblvlm)
                print*,'hq_mes ',(hq_mes(i),i=1,nblvlm)
@@ -213 +213 @@
                   endif
           IF (ts_fcg) print*,'ts_subr', ts_subr
-c*** on interpole les champs meso_NH sur les niveaux de pression***
-c*** gcm . on obtient le nouveau champ after                    ***
+!*** on interpole les champs meso_NH sur les niveaux de pression***
+!*** gcm . on obtient le nouveau champ after                    ***
             do k=1,klev
              if (JM(k) .eq. 0) then
@@ -237 +237 @@
                endif ! imp_fcg
                if(Turb_fcg) then
-           hThTurbaft(k)=coef1(k)*hThTurb_mes(jm(k))
-     $               +coef2(k)*hThTurb_mes(jm(k)+1)
-           hqTurbaft(k) =coef1(k)*hqTurb_mes(jm(k))
-     $               +coef2(k)*hqTurb_mes(jm(k)+1)
+           hThTurbaft(k)=coef1(k)*hThTurb_mes(jm(k))                            &
+     &               +coef2(k)*hThTurb_mes(jm(k)+1)
+           hqTurbaft(k) =coef1(k)*hqTurb_mes(jm(k))                             &
+     &               +coef2(k)*hqTurb_mes(jm(k)+1)
                endif ! Turb_fcg
              endif ! JM(k) .eq. 0
@@ -250 +250 @@
          endif ! pas.gt.pasprev    fin du bloc relatif au passage au pas
                                   !de temps (meso) suivant
-c*** si on atteint le pas max des donnees experimentales ,on     ***
-c*** on conserve les derniers champs calcules                    ***
+!*** si on atteint le pas max des donnees experimentales ,on     ***
+!*** on conserve les derniers champs calcules                    ***
       if(temps.ge.pasmax)then
           do ll=1,klev
@@ -264 +264 @@
           ts_subr = tsaft
       else ! temps.ge.pasmax
-c*** on interpole sur les pas de temps de 10mn du gcm a partir   ***
-c** des pas de temps de 1h du meso_NH                            ***
+!*** on interpole sur les pas de temps de 10mn du gcm a partir   ***
+!** des pas de temps de 1h du meso_NH                            ***
          do j=1,klev
          ht(j)=((timeaft-time)*htbef(j)+(time-timebef)*htaft(j))/dt
@@ -275 +275 @@
              endif ! imp_fcg
              if(Turb_fcg) then
-         hThTurb(j)=((timeaft-time)*hThTurbbef(j)
-     $           +(time-timebef)*hThTurbaft(j))/dt
-         hqTurb(j)= ((timeaft-time)*hqTurbbef(j)
-     $           +(time-timebef)*hqTurbaft(j))/dt
+         hThTurb(j)=((timeaft-time)*hThTurbbef(j)                           &
+     &           +(time-timebef)*hThTurbaft(j))/dt
+         hqTurb(j)= ((timeaft-time)*hqTurbbef(j)                            &
+     &           +(time-timebef)*hqTurbaft(j))/dt
              endif ! Turb_fcg
          enddo
          ts_subr = ((timeaft-time)*tsbef + (time-timebef)*tsaft)/dt
        endif ! temps.ge.pasmax
-c
+!
         print *,' time,timebef,timeaft',time,timebef,timeaft
         print *,' ht,htbef,htaft,hthturb,hthturbbef,hthturbaft'
         do j= 1,klev
-           print *, j,ht(j),htbef(j),htaft(j),
-     $             hthturb(j),hthturbbef(j),hthturbaft(j)
+           print *, j,ht(j),htbef(j),htaft(j),                              &
+     &             hthturb(j),hthturbbef(j),hthturbaft(j)
         enddo
         print *,' hq,hqbef,hqaft,hqturb,hqturbbef,hqturbaft'
         do j= 1,klev
-           print *, j,hq(j),hqbef(j),hqaft(j),
-     $             hqturb(j),hqturbbef(j),hqturbaft(j)
+           print *, j,hq(j),hqbef(j),hqaft(j),                              &
+     &             hqturb(j),hqturbbef(j),hqturbaft(j)
         enddo
-c
-c-------------------------------------------------------------------
-c
+!
+!-------------------------------------------------------------------
+!
          IF (Ts_fcg) Ts = Ts_subr
          return
-c
-c-----------------------------------------------------------------------
-c on sort les champs de "convergence" pour l instant initial 'in'
-c ceci se passe au pas temps itap=0 de la physique
-c-----------------------------------------------------------------------
-        entry get_uvd2(itap,dtime,file_forctl,file_fordat,
-     .           ht,hq,hw,hu,hv,hThTurb,hqTurb,ts,
-     .           imp_fcg,ts_fcg,Tp_fcg,Turb_fcg)
+!
+!-----------------------------------------------------------------------
+! on sort les champs de "convergence" pour l instant initial 'in'
+! ceci se passe au pas temps itap=0 de la physique
+!-----------------------------------------------------------------------
+        entry get_uvd2(itap,dtime,file_forctl,file_fordat,                  &
+     &           ht,hq,hw,hu,hv,hThTurb,hqTurb,ts,                          &
+     &           imp_fcg,ts_fcg,Tp_fcg,Turb_fcg)
              print*,'le pas itap est:',itap
-c
-c===================================================================
-c
+!
+!===================================================================
+!
        write(*,'(a)') 'OPEN '//file_forctl
        open(97,FILE=file_forctl,FORM='FORMATTED')
-c
-c------------------
+!
+!------------------
       do i=1,1000
       read(97,1000,end=999) string
@@ -322 +322 @@
       enddo
  50   backspace(97)
-c-------------------------------------------------------------------
-c   *** on lit le pas de temps dans le fichier de donnees ***
-c   *** "forcing.ctl" et pasmax                           ***
-c-------------------------------------------------------------------
+!-------------------------------------------------------------------
+!   *** on lit le pas de temps dans le fichier de donnees ***
+!   *** "forcing.ctl" et pasmax                           ***
+!-------------------------------------------------------------------
       read(97,2000) aaa
  2000  format (a80)
@@ -344 +344 @@
          print*,'pasmax est',pasmax
       CLOSE(97)
-c------------------------------------------------------------------
-c *** on lit le pas de temps initial de la simulation ***
-c------------------------------------------------------------------
+!------------------------------------------------------------------
+! *** on lit le pas de temps initial de la simulation ***
+!------------------------------------------------------------------
                   in=itap
-cc                  pasprev=in
-cc                  time0=dt*(pasprev-1)
+!c                  pasprev=in
+!c                  time0=dt*(pasprev-1)
                   pasprev=in-1
                   time0=dt*pasprev
-C
+!
           close(98)
-c
+!
       write(*,'(a)') 'OPEN '//file_fordat
-      open(99,FILE=file_fordat,FORM='UNFORMATTED',
-     .          ACCESS='DIRECT',RECL=8)
+      open(99,FILE=file_fordat,FORM='UNFORMATTED',                          &
+     &          ACCESS='DIRECT',RECL=8)
           icomp1 = nblvlm*4
           IF (ts_fcg) icomp1 = icomp1 + 1
@@ -363 +363 @@
           IF (Turb_fcg) icomp1 = icomp1 + nblvlm*2
           icompt = icomp1*(in-1)
-          call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes
-     .                   ,hu_mes,hv_mes,hthturb_mes,hqturb_mes
-     .                   ,ts_fcg,ts_subr,imp_fcg,Turb_fcg)
+          call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes              &
+     &                   ,hu_mes,hv_mes,hthturb_mes,hqturb_mes              &
+     &                   ,ts_fcg,ts_subr,imp_fcg,Turb_fcg)
           print *, 'get_uvd : rdgrads ->'
           print *, tp_fcg
-c
-c following commented out because we have temperature already in ARM case
-c   (otherwise this is the potential temperature )
-c------------------------------------------------------------------------
+!
+! following commented out because we have temperature already in ARM case
+!   (otherwise this is the potential temperature )
+!------------------------------------------------------------------------
                if(Tp_fcg) then
           do i = 1,nblvlm
@@ -389 +389 @@
                  print*,'hqTurb ',     (hqTurb_mes(i),i=1,nblvlm)
               endif ! Turb_fcg
-c----------------------------------------------------------------------
-c on a obtenu des champs initiaux sur les niveaux du meso_NH
-c on interpole sur les niveaux du gcm(niveau pression bien sur!)
-c-----------------------------------------------------------------------
+!----------------------------------------------------------------------
+! on a obtenu des champs initiaux sur les niveaux du meso_NH
+! on interpole sur les niveaux du gcm(niveau pression bien sur!)
+!-----------------------------------------------------------------------
             do k=1,klev
              if (JM(k) .eq. 0) then
-cFKC bug? ne faut il pas convertir tsol en tendance ????
-cRT bug taken care of by removing the stuff
+!FKC bug? ne faut il pas convertir tsol en tendance ????
+!RT bug taken care of by removing the stuff
            htaft(k)=              ht_mes(jm(k)+1)
            hqaft(k)=              hq_mes(jm(k)+1)
@@ -417 +417 @@
                endif ! imp_fcg
                if(Turb_fcg) then
-           hThTurbaft(k)=coef1(k)*hThTurb_mes(jm(k))
-     $                  +coef2(k)*hThTurb_mes(jm(k)+1)
-           hqTurbaft(k) =coef1(k)*hqTurb_mes(jm(k))
-     $                  +coef2(k)*hqTurb_mes(jm(k)+1)
+           hThTurbaft(k)=coef1(k)*hThTurb_mes(jm(k))                        &
+     &                  +coef2(k)*hThTurb_mes(jm(k)+1)
+           hqTurbaft(k) =coef1(k)*hqTurb_mes(jm(k))                         &
+     &                  +coef2(k)*hqTurb_mes(jm(k)+1)
                endif ! Turb_fcg
              endif ! JM(k) .eq. 0
             enddo
             tsaft = ts_subr
-c valeurs initiales des champs de convergence
+! valeurs initiales des champs de convergence
           do k=1,klev
              ht(k)=htaft(k)
@@ -442 +442 @@
           close(99)
           close(98)
-c
-c-------------------------------------------------------------------
-c
-c
+!
+!-------------------------------------------------------------------
+!
+!
  100      IF (Ts_fcg) Ts = Ts_subr
         return
-c
+!
 999     continue
         stop 'erreur lecture, file forcing.ctl'
         end
 
-      SUBROUTINE advect_tvl(dtime,zt,zq,vu_f,vv_f,t_f,q_f
-     :                     ,d_t_adv,d_q_adv)
+      SUBROUTINE advect_tvl(dtime,zt,zq,vu_f,vv_f,t_f,q_f                   &
+     &                     ,d_t_adv,d_q_adv)
       use dimphy
       implicit none
 
 #include "dimensions.h"
-ccccc#include "dimphy.h"
+!cccc#include "dimphy.h"
 
       integer k
       real dtime, fact, du, dv, cx, cy, alx, aly
       real zt(klev), zq(klev,3)
-     :   , vu_f(klev), vv_f(klev), t_f(klev), q_f(klev,3)
+      real vu_f(klev), vv_f(klev), t_f(klev), q_f(klev,3)
 
       real d_t_adv(klev), d_q_adv(klev,3)
 
-c Velocity of moving cell
+! Velocity of moving cell
       data cx,cy /12., -2./
 
-c Dimensions of moving cell
-      data alx,aly /100 000.,150 000./
+! Dimensions of moving cell
+      data alx,aly /100000.,150000./
 
       do k = 1, klev
@@ -490 +490 @@
       implicit none
 
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c cette routine remplit les COMMON com1_phys_gcss et com2_phys_gcss.h
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! cette routine remplit les COMMON com1_phys_gcss et com2_phys_gcss.h
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
       INTEGER klev !nombre de niveau de pression du GCM
@@ -514 +514 @@
       klev = klevgcm
 
-c---------------------------------------------------------------------
-c pression au milieu des couches du gcm dans la physiq
-c (SB: remplace le call conv_lipress_gcm(playgcm) )
-c---------------------------------------------------------------------
+!---------------------------------------------------------------------
+! pression au milieu des couches du gcm dans la physiq
+! (SB: remplace le call conv_lipress_gcm(playgcm) )
+!---------------------------------------------------------------------
 
        do k = 1, klev
@@ -524 +524 @@
        enddo
 
-c----------------------------------------------------------------------
-c lecture du descripteur des donnees Meso-NH (forcing.ctl):
-c  -> nb niveaux du meso.NH (nblvlm) + pressions meso.NH
-c (on remplit le COMMON com2_phys_gcss)
-c----------------------------------------------------------------------
+!----------------------------------------------------------------------
+! lecture du descripteur des donnees Meso-NH (forcing.ctl):
+!  -> nb niveaux du meso.NH (nblvlm) + pressions meso.NH
+! (on remplit le COMMON com2_phys_gcss)
+!----------------------------------------------------------------------
 
       call mesolupbis(file_forctl)
@@ -534 +534 @@
       print*,'la valeur de nblvlm est:',nblvlm
 
-c----------------------------------------------------------------------
-c etude de la correspondance entre les niveaux meso.NH et GCM;
-c calcul des coefficients d interpolation coef1 et coef2
-c (on remplit le COMMON com1_phys_gcss)
-c----------------------------------------------------------------------
+!----------------------------------------------------------------------
+! etude de la correspondance entre les niveaux meso.NH et GCM;
+! calcul des coefficients d interpolation coef1 et coef2
+! (on remplit le COMMON com1_phys_gcss)
+!----------------------------------------------------------------------
 
       call corresbis(psolgcm)
 
-c---------------------------------------------------------
-c TEST sur le remplissage de com1_phys_gcss et com2_phys_gcss:
-c---------------------------------------------------------
+!---------------------------------------------------------
+! TEST sur le remplissage de com1_phys_gcss et com2_phys_gcss:
+!---------------------------------------------------------
  
       write(*,*) ' '
@@ -562 +562 @@
       SUBROUTINE mesolupbis(file_forctl)
       implicit none 
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-c Lecture descripteur des donnees MESO-NH (forcing.ctl):
-c -------------------------------------------------------
-c
-c     Cette subroutine lit dans le fichier de controle "essai.ctl"
-c     et affiche le nombre de niveaux du Meso-NH ainsi que les valeurs
-c     des pressions en milieu de couche du Meso-NH (en Pa puis en hPa).
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
+!
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!
+! Lecture descripteur des donnees MESO-NH (forcing.ctl):
+! -------------------------------------------------------
+!
+!     Cette subroutine lit dans le fichier de controle "essai.ctl"
+!     et affiche le nombre de niveaux du Meso-NH ainsi que les valeurs
+!     des pressions en milieu de couche du Meso-NH (en Pa puis en hPa).
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!
       INTEGER nblvlm !nombre de niveau de pression du mesoNH
       REAL playm(100)  !pression en Pa milieu de chaque couche Meso-NH
@@ -588 +588 @@
       lu=9
       open(lu,file=file_forctl,form='formatted')
-c
+!
       do i=1,1000
       read(lu,1000,end=999) a
       if (a .eq. 'ZDEF') go to 100
       enddo
-c
+!
  100  backspace(lu)
       print*,'  DESCRIPTION DES 2 MODELES : '
       print*,' '
-c
+!
       read(lu,2000) aaa
  2000  format (a80)
@@ -604 +604 @@
          read(anblvl,*) nblvlm
 
-c
+!
       print*,'nbre de niveaux de pression Meso-NH :',nblvlm
       print*,' '
       print*,'pression en Pa de chaque couche du meso-NH :'
-c
+!
       read(lu,*) (playm(mlz),mlz=1,nblvlm)
-c      Si la pression est en HPa, la multiplier par 100
+!      Si la pression est en HPa, la multiplier par 100
       if (playm(1) .lt. 10000.) then
         do mlz = 1,nblvlm
@@ -617 +617 @@
       endif
       print*,(playm(mlz),mlz=1,nblvlm)
-c
+!
  1000 format (a4)
  1001 format(5x,i2)
-c
+!
       print*,' '
       do mlzh=1,nblvlm
       hplaym(mlzh)=playm(mlzh)/100.
       enddo
-c
+!
       print*,'pression en hPa de chaque couche du meso-NH: '
       print*,(hplaym(mlzh),mlzh=1,nblvlm)
-c
+!
       close (lu)
       return
-c
+!
  999  stop 'erreur lecture des niveaux pression des donnees'
       end
 
-      SUBROUTINE rdgrads(itape,icount,nl,z,ht,hq,hw,hu,hv,hthtur,hqtur,
-     :  ts_fcg,ts,imp_fcg,Turb_fcg)
+      SUBROUTINE rdgrads(itape,icount,nl,z,ht,hq,hw,hu,hv,hthtur,hqtur,     &
+     &  ts_fcg,ts,imp_fcg,Turb_fcg)
       IMPLICIT none
       INTEGER itape,icount,icomp, nl
@@ -642 +642 @@
       real hthtur(nl),hqtur(nl)
       real ts
-c
+!
       INTEGER k
-c
+!
       LOGICAL imp_fcg,ts_fcg,Turb_fcg
-c
+!
       icomp = icount
-c
-c
+!
+!
          do k=1,nl
             icomp=icomp+1
@@ -664 +664 @@
             read(itape,rec=icomp)hQ(k)
          enddo
-c
+!
              if(turb_fcg) then
          do k=1,nl
@@ -676 +676 @@
              endif
          print *,' apres lecture hthtur, hqtur'
-c
+!
           if(imp_fcg) then
 
@@ -689 +689 @@
 
           endif
-c
+!
          do k=1,nl
             icomp=icomp+1
             read(itape,rec=icomp)hw(k)
          enddo
-c
+!
               if(ts_fcg) then
          icomp=icomp+1
          read(itape,rec=icomp)ts
               endif
-c
+!
       print *,' rdgrads ->'
 
@@ -708 +708 @@
       implicit none
 
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c Cette subroutine calcule et affiche les valeurs des coefficients
-c d interpolation qui serviront dans la formule d interpolation elle-
-c meme.
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! Cette subroutine calcule et affiche les valeurs des coefficients
+! d interpolation qui serviront dans la formule d interpolation elle-
+! meme.
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
       INTEGER klev    !nombre de niveau de pression du GCM
@@ -737 +737 @@
           mlz = 0
           JM(1) = mlz
-          coef1(1)=(playm(mlz+1)-val)
-     *             /(playm(mlz+1)-psol)
-          coef2(1)=(val-psol)
-     *             /(playm(mlz+1)-psol)
+          coef1(1)=(playm(mlz+1)-val)/(playm(mlz+1)-psol)
+          coef2(1)=(val-psol)/(playm(mlz+1)-psol)
        else if (val .gt. playm(nblvlm)) then
          do mlz=1,nblvlm
-          if (     val .le. playm(mlz)
-     *       .and. val .gt. playm(mlz+1))then
+          if (     val .le. playm(mlz).and. val .gt. playm(mlz+1))then
            JM(k)=mlz
-           coef1(k)=(playm(mlz+1)-val)
-     *              /(playm(mlz+1)-playm(mlz))
-           coef2(k)=(val-playm(mlz))
-     *              /(playm(mlz+1)-playm(mlz))
+           coef1(k)=(playm(mlz+1)-val)/(playm(mlz+1)-playm(mlz))
+           coef2(k)=(val-playm(mlz))/(playm(mlz+1)-playm(mlz))
           endif
          enddo
@@ -758 +753 @@
        endif
       enddo
-c
-cc      if (play(klev) .le. playm(nblvlm)) then
-cc         mlz=nblvlm-1
-cc         JM(klev)=mlz
-cc         coef1(klev)=(playm(mlz+1)-val)
-cc     *            /(playm(mlz+1)-playm(mlz))
-cc         coef2(klev)=(val-playm(mlz))
-cc     *            /(playm(mlz+1)-playm(mlz))
-cc      endif
-c
+!
+!c      if (play(klev) .le. playm(nblvlm)) then
+!c         mlz=nblvlm-1
+!c         JM(klev)=mlz
+!c         coef1(klev)=(playm(mlz+1)-val)
+!c     *            /(playm(mlz+1)-playm(mlz))
+!c         coef2(klev)=(val-playm(mlz))
+!c     *            /(playm(mlz+1)-playm(mlz))
+!c      endif
+!
       print*,' '
       print*,'         INTERPOLATION  : '
@@ -776 +771 @@
       print*,(JM(k),k=1,klev)
       print*,' '
-      print*,'valeurs du premier coef d"interpolation pour les 9 niveaux
-     *: '
+      print*,'vals du premier coef d"interpolation pour les 9 niveaux: '
       print*,(coef1(k),k=1,klev)
       print*,' '
-      print*,'valeurs du deuxieme coef d"interpolation pour les 9 niveau
-     *x: '
+      print*,'valeurs du deuxieme coef d"interpolation pour les 9 niveaux:'
       print*,(coef2(k),k=1,klev)
-c
+!
       return
       end
       SUBROUTINE GETSCH(STR,DEL,TRM,NTH,SST,NCH)
-C***************************************************************
-C*                                                             *
-C*                                                             *
-C* GETSCH                                                      *
-C*                                                             *
-C*                                                             *
-C* modified by :                                               *
-C***************************************************************
-C*   Return in SST the character string found between the NTH-1 and NTH
-C*   occurence of the delimiter 'DEL' but before the terminator 'TRM' in
-C*   the input string 'STR'. If TRM=DEL then STR is considered unlimited.
-C*   NCH=Length of the string returned in SST or =-1 if NTH is <1 or if
-C*   NTH is greater than the number of delimiters in STR.
+!***************************************************************
+!*                                                             *
+!*                                                             *
+!* GETSCH                                                      *
+!*                                                             *
+!*                                                             *
+!* modified by :                                               *
+!***************************************************************
+!*   Return in SST the character string found between the NTH-1 and NTH
+!*   occurence of the delimiter 'DEL' but before the terminator 'TRM' in
+!*   the input string 'STR'. If TRM=DEL then STR is considered unlimited.
+!*   NCH=Length of the string returned in SST or =-1 if NTH is <1 or if
+!*   NTH is greater than the number of delimiters in STR.
       IMPLICIT INTEGER (A-Z)
       CHARACTER STR*(*),DEL*1,TRM*1,SST*(*)
@@ -811 +804 @@
           IF(LENGTH.LT.0) LENGTH=LEN(STR)
         ENDIF
-C*     Find beginning and end of the NTH DEL-limited substring in STR
+!*     Find beginning and end of the NTH DEL-limited substring in STR
         END=-1
         DO 1,N=1,NTH
@@ -818 +811 @@
         END=BEG+INDEX(STR(BEG:LENGTH),DEL)-2
         IF(END.EQ.BEG-2) END=LENGTH
-C*        PRINT *,'NTH,LENGTH,N,BEG,END=',NTH,LENGTH,N,BEG,END
+!*        PRINT *,'NTH,LENGTH,N,BEG,END=',NTH,LENGTH,N,BEG,END
     1   CONTINUE
         NCH=END-BEG+1
@@ -825 +818 @@
       END
       CHARACTER*(*) FUNCTION SPACES(STR,NSPACE)
-C
-C CERN PROGLIB# M433    SPACES          .VERSION KERNFOR  4.14  860211
-C ORIG.  6/05/86 M.GOOSSENS/DD
-C
-C-    The function value SPACES returns the character string STR with
-C-    leading blanks removed and each occurence of one or more blanks
-C-    replaced by NSPACE blanks inside the string STR
-C
+!
+! CERN PROGLIB# M433    SPACES          .VERSION KERNFOR  4.14  860211
+! ORIG.  6/05/86 M.GOOSSENS/DD
+!
+!-    The function value SPACES returns the character string STR with
+!-    leading blanks removed and each occurence of one or more blanks
+!-    replaced by NSPACE blanks inside the string STR
+!
       CHARACTER*(*) STR
-C
+!
       LENSPA = LEN(SPACES)
       SPACES = ' '
@@ -857 +850 @@
   999 END
       FUNCTION INDEXC(STR,SSTR)
-C
-C CERN PROGLIB# M433    INDEXC          .VERSION KERNFOR  4.14  860211
-C ORIG. 26/03/86 M.GOOSSENS/DD
-C
-C-    Find the leftmost position where substring SSTR does not match
-C-    string STR scanning forward
-C
+!
+! CERN PROGLIB# M433    INDEXC          .VERSION KERNFOR  4.14  860211
+! ORIG. 26/03/86 M.GOOSSENS/DD
+!
+!-    Find the leftmost position where substring SSTR does not match
+!-    string STR scanning forward
+!
       CHARACTER*(*) STR,SSTR
-C
+!
       LENS   = LEN(STR)
       LENSS  = LEN(SSTR)
-C
+!
       DO 10 I=1,LENS-LENSS+1
           IF (STR(I:I+LENSS-1).NE.SSTR) THEN
@@ -876 +869 @@
    10 CONTINUE
       INDEXC = 0
-C
+!
   999 END

LMDZ5/trunk/libf/phylmd/compar1d.h

@@ -27 +27 @@
       logical :: ok_old_disvert
 
-      common/com_par1d/
+      common/com_par1d/                                                 &
      & nat_surf,tsurf,rugos,                                            &
      & qsol,qsurf,psurf,zsurf,albedo,time,time_ini,xlat,xlon,airefi,    &
      & wtsurf,wqsurf,restart_runoff,xagesno,qsolinp,zpicinp,            &
-     & forcing_type,
+     & forcing_type,                                                    &
      & restart,ok_old_disvert
 

LMDZ5/trunk/libf/phylmd/fcg_gcssold.h

@@ -2 +2 @@
 ! $Id: fcg_gcssold.h 2010-08-10 17:02:56Z lahellec $
 !
-      logical :: imp_fcg_gcssold,ts_fcg_gcssold,Tp_fcg_gcssold,
-     $ Tp_ini_gcssold,
-     $ xTurb_fcg_gcssold
+      logical :: imp_fcg_gcssold,ts_fcg_gcssold,Tp_fcg_gcssold
+      logical :: Tp_ini_gcssold
+      logical :: xTurb_fcg_gcssold
 
-      common /fcg_gcssold/imp_fcg_gcssold,ts_fcg_gcssold,Tp_fcg_gcssold,
-     $ Tp_ini_gcssold,
-     $ xTurb_fcg_gcssold
+      common /fcg_gcssold/imp_fcg_gcssold,ts_fcg_gcssold,Tp_fcg_gcssold,        & 
+     & Tp_ini_gcssold,                                                          &
+     & xTurb_fcg_gcssold
 
 !$OMP THREADPRIVATE(/fcg_gcssold/)

LMDZ5/trunk/libf/phylmd/lmdz1d.F90

@@ -4 +4 @@
 #include "../dyn3d_common/disvert.F90"
 
-      CALL lmdz1d
-
-      END
+
+      PROGRAM lmdz1d
+
+      USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar
+      use phys_state_var_mod
+      use comgeomphy
+      use dimphy
+      use surface_data, only : type_ocean,ok_veget
+      use pbl_surface_mod, only : ftsoil, pbl_surface_init,                     &
+     &                            pbl_surface_final
+      use fonte_neige_mod, only : fonte_neige_init, fonte_neige_final 
+
+      use infotrac ! new
+      use control_mod
+      USE indice_sol_mod
+      USE phyaqua_mod
+
+      implicit none
+#include "dimensions.h"
+#include "YOMCST.h"
+#include "temps.h"
+!!#include "control.h"
+#include "iniprint.h"
+#include "clesphys.h"
+#include "dimsoil.h"
+!#include "indicesol.h"
+
+#include "comvert.h"
+#include "compar1d.h"
+#include "flux_arp.h"
+#include "tsoilnudge.h"
+#include "fcg_gcssold.h"
+!!!#include "fbforcing.h"
+
+!=====================================================================
+! DECLARATIONS
+!=====================================================================
+
+!---------------------------------------------------------------------
+!  Externals
+!---------------------------------------------------------------------
+      external fq_sat
+      real fq_sat
+
+!---------------------------------------------------------------------
+!  Arguments d' initialisations de la physique (USER DEFINE)
+!---------------------------------------------------------------------
+
+      integer, parameter :: ngrid=1
+      real :: zcufi    = 1.
+      real :: zcvfi    = 1.
+
+!-      real :: nat_surf
+!-      logical :: ok_flux_surf
+!-      real :: fsens
+!-      real :: flat
+!-      real :: tsurf
+!-      real :: rugos
+!-      real :: qsol(1:2)
+!-      real :: qsurf
+!-      real :: psurf
+!-      real :: zsurf
+!-      real :: albedo
+!-
+!-      real :: time     = 0.
+!-      real :: time_ini
+!-      real :: xlat 
+!-      real :: xlon 
+!-      real :: wtsurf 
+!-      real :: wqsurf 
+!-      real :: restart_runoff 
+!-      real :: xagesno 
+!-      real :: qsolinp 
+!-      real :: zpicinp 
+!-
+      real :: fnday 
+      real :: day, daytime 
+      real :: day1
+      real :: heure
+      integer :: jour
+      integer :: mois
+      integer :: an
+ 
+!---------------------------------------------------------------------
+!  Declarations related to forcing and initial profiles 
+!---------------------------------------------------------------------
+
+        integer :: kmax = llm
+        integer llm700,nq1,nq2
+        INTEGER, PARAMETER :: nlev_max=1000, nqmx=1000
+        real timestep, frac
+        real height(nlev_max),tttprof(nlev_max),qtprof(nlev_max)
+        real  uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max)
+        real  ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max)
+        real  dqtdxls(nlev_max),dqtdyls(nlev_max)
+        real  dqtdtls(nlev_max),thlpcar(nlev_max)
+        real  qprof(nlev_max,nqmx)
+
+!        integer :: forcing_type
+        logical :: forcing_les     = .false.
+        logical :: forcing_armcu   = .false.
+        logical :: forcing_rico    = .false.
+        logical :: forcing_radconv = .false.
+        logical :: forcing_toga    = .false.
+        logical :: forcing_twpice  = .false.
+        logical :: forcing_amma    = .false.
+        logical :: forcing_GCM2SCM = .false.
+        logical :: forcing_GCSSold = .false.
+        logical :: forcing_sandu   = .false.
+        logical :: forcing_astex   = .false.
+        logical :: forcing_fire    = .false.
+        integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
+!                                                            (cf read_tsurf1d.F)
+
+!vertical advection computation
+!       real d_t_z(llm), d_q_z(llm)
+!       real d_t_dyn_z(llm), d_q_dyn_z(llm)
+!       real zz(llm)
+!       real zfact
+
+!flag forcings
+        logical :: nudge_wind=.true.
+        logical :: nudge_thermo=.false.
+        logical :: cptadvw=.true.
+!=====================================================================
+! DECLARATIONS FOR EACH CASE
+!=====================================================================
+!
+#include "1D_decl_cases.h"
+!
+!---------------------------------------------------------------------
+!  Declarations related to vertical discretization:
+!---------------------------------------------------------------------
+      real :: pzero=1.e5
+      real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1)
+      real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1),poub
+
+!---------------------------------------------------------------------
+!  Declarations related to variables
+!---------------------------------------------------------------------
+
+      real :: phi(llm)
+      real :: teta(llm),tetal(llm),temp(llm),u(llm),v(llm),w(llm)
+      real :: rlat_rad(1),rlon_rad(1)
+      real :: omega(llm+1),omega2(llm),rho(llm+1)
+      real :: ug(llm),vg(llm),fcoriolis
+      real :: sfdt, cfdt
+      real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
+      real :: dt_dyn(llm)
+      real :: dt_cooling(llm),d_th_adv(llm)
+      real :: alpha
+      real :: ttt
+
+      REAL, ALLOCATABLE, DIMENSION(:,:):: q
+      REAL, ALLOCATABLE, DIMENSION(:,:):: dq
+      REAL, ALLOCATABLE, DIMENSION(:,:):: dq_dyn
+      REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
+!     REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv
+
+!---------------------------------------------------------------------
+!  Initialization of surface variables
+!---------------------------------------------------------------------
+      real :: run_off_lic_0(1)
+      real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf)
+      real :: evap(1,nbsrf),frugs(1,nbsrf)
+      real :: tsoil(1,nsoilmx,nbsrf)
+      real :: agesno(1,nbsrf)
+
+!---------------------------------------------------------------------
+!  Call to phyredem
+!---------------------------------------------------------------------
+      logical :: ok_writedem =.true.
+      
+!---------------------------------------------------------------------
+!  Call to physiq
+!---------------------------------------------------------------------
+      integer, parameter :: longcles=20
+      logical :: firstcall=.true.
+      logical :: lastcall=.false.
+      real :: phis    = 0.0
+      real :: clesphy0(longcles) = 0.0
+      real :: dpsrf 
+
+!---------------------------------------------------------------------
+!  Initializations of boundary conditions 
+!---------------------------------------------------------------------
+      integer, parameter :: yd = 360
+      real :: phy_nat (yd) = 0.0 ! 0=ocean libre,1=land,2=glacier,3=banquise
+      real :: phy_alb (yd)  ! Albedo land only (old value condsurf_jyg=0.3)
+      real :: phy_sst (yd)  ! SST (will not be used; cf read_tsurf1d.F)
+      real :: phy_bil (yd) = 1.0 ! Ne sert que pour les slab_ocean
+      real :: phy_rug (yd) ! Longueur rugosite utilisee sur land only
+      real :: phy_ice (yd) = 0.0 ! Fraction de glace
+      real :: phy_fter(yd) = 0.0 ! Fraction de terre
+      real :: phy_foce(yd) = 0.0 ! Fraction de ocean
+      real :: phy_fsic(yd) = 0.0 ! Fraction de glace
+      real :: phy_flic(yd) = 0.0 ! Fraction de glace
+
+!---------------------------------------------------------------------
+!  Fichiers et d'autres variables
+!---------------------------------------------------------------------
+      integer :: k,l,i,it=1,mxcalc
+      integer jjmp1
+      parameter (jjmp1=jjm+1-1/jjm)
+      INTEGER nbteta
+      PARAMETER(nbteta=3)
+      REAL dudyn(iim+1,jjmp1,llm)
+      REAL PVteta(1,nbteta)
+      INTEGER read_climoz 
+!Al1
+      integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
+      data ecrit_slab_oc/-1/
+
+!=====================================================================
+! INITIALIZATIONS 
+!=====================================================================
+! Initialization of Common turb_forcing
+       dtime_frcg = 0.
+       Turb_fcg_gcssold=.false.
+       hthturb_gcssold = 0.
+       hqturb_gcssold = 0.
+
+!---------------------------------------------------------------------
+! OPTIONS OF THE 1D SIMULATION (lmdz1d.def => unicol.def)
+!---------------------------------------------------------------------
+!Al1
+        call conf_unicol
+!Al1 moves this gcssold var from common fcg_gcssold to 
+        Turb_fcg_gcssold = xTurb_fcg_gcssold
+! --------------------------------------------------------------------
+        close(1)
+!Al1
+        write(*,*) 'lmdz1d.def lu => unicol.def'
+
+! forcing_type defines the way the SCM is forced:
+!forcing_type = 0 ==> forcing_les = .true.
+!             initial profiles from file prof.inp.001
+!             no forcing by LS convergence ; 
+!             surface temperature imposed ;
+!             radiative cooling may be imposed (iflag_radia=0 in physiq.def)
+!forcing_type = 1 ==> forcing_radconv = .true.
+!             idem forcing_type = 0, but the imposed radiative cooling 
+!             is set to 0 (hence, if iflag_radia=0 in physiq.def, 
+!             then there is no radiative cooling at all)
+!forcing_type = 2 ==> forcing_toga = .true.
+!             initial profiles from TOGA-COARE IFA files 
+!             LS convergence and SST imposed from TOGA-COARE IFA files 
+!forcing_type = 3 ==> forcing_GCM2SCM = .true.
+!             initial profiles from the GCM output
+!             LS convergence imposed from the GCM output
+!forcing_type = 4 ==> forcing_twpice = .true.
+!             initial profiles from TWP-ICE cdf file 
+!             LS convergence, omega and SST imposed from TWP-ICE files 
+!forcing_type = 5 ==> forcing_rico = .true.
+!             initial profiles from RICO files 
+!             LS convergence imposed from RICO files 
+!forcing_type = 6 ==> forcing_amma = .true.
+!             initial profiles from AMMA nc file 
+!             LS convergence, omega and surface fluxes imposed from AMMA file  
+!forcing_type = 40 ==> forcing_GCSSold = .true.
+!             initial profile from GCSS file
+!             LS convergence imposed from GCSS file
+!forcing_type = 50 ==> forcing_fire = .true.
+!             forcing from fire.nc
+!forcing_type = 59 ==> forcing_sandu = .true.
+!             initial profiles from sanduref file: see prof.inp.001
+!             SST varying with time and divergence constante: see ifa_sanduref.txt file
+!             Radiation has to be computed interactively
+!forcing_type = 60 ==> forcing_astex = .true.
+!             initial profiles from file: see prof.inp.001 
+!             SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
+!             Radiation has to be computed interactively
+!forcing_type = 61 ==> forcing_armcu = .true.
+!             initial profiles from file: see prof.inp.001 
+!             sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
+!             large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
+!             use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s 
+!             Radiation to be switched off
+!
+      if (forcing_type .eq.0) THEN
+       forcing_les = .true.
+      elseif (forcing_type .eq.1) THEN
+       forcing_radconv = .true.
+      elseif (forcing_type .eq.2) THEN
+       forcing_toga    = .true.
+      elseif (forcing_type .eq.3) THEN
+       forcing_GCM2SCM = .true.
+      elseif (forcing_type .eq.4) THEN
+       forcing_twpice = .true.
+      elseif (forcing_type .eq.5) THEN
+       forcing_rico = .true.
+      elseif (forcing_type .eq.6) THEN
+       forcing_amma = .true.
+      elseif (forcing_type .eq.40) THEN
+       forcing_GCSSold = .true.
+      elseif (forcing_type .eq.50) THEN
+       forcing_fire = .true.
+      elseif (forcing_type .eq.59) THEN
+       forcing_sandu   = .true.
+      elseif (forcing_type .eq.60) THEN
+       forcing_astex   = .true.
+      elseif (forcing_type .eq.61) THEN
+       forcing_armcu = .true.
+       IF(llm.NE.19.AND.llm.NE.40) stop 'Erreur nombre de niveaux !!'
+      else
+       write (*,*) 'ERROR : unknown forcing_type ', forcing_type
+       stop 'Forcing_type should be 0,1,2,3,4,5,6 or 40,59,60,61'
+      ENDIF
+      print*,"forcing type=",forcing_type
+
+! if type_ts_forcing=0, the surface temp of 1D simulation is constant in time
+! (specified by tsurf in lmdz1d.def); if type_ts_forcing=1, the surface temperature
+! varies in time according to a forcing (e.g. forcing_toga) and is passed to read_tsurf1d.F
+! through the common sst_forcing.
+
+        type_ts_forcing = 0
+        if (forcing_toga.or.forcing_sandu.or.forcing_astex)                 &
+     &    type_ts_forcing = 1
+
+!---------------------------------------------------------------------
+!  Definition of the run
+!---------------------------------------------------------------------
+
+      call conf_gcm( 99, .TRUE. , clesphy0 )
+!-----------------------------------------------------------------------
+!   Choix du calendrier
+!   -------------------
+
+!      calend = 'earth_365d'
+      if (calend == 'earth_360d') then
+        call ioconf_calendar('360d')
+        write(*,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
+      else if (calend == 'earth_365d') then
+        call ioconf_calendar('noleap')
+        write(*,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
+      else if (calend == 'earth_366d') then
+        call ioconf_calendar('all_leap')
+        write(*,*)'CALENDRIER CHOISI: Terrestre bissextile'
+      else if (calend == 'gregorian') then
+        call ioconf_calendar('gregorian') ! not to be used by normal users
+        write(*,*)'CALENDRIER CHOISI: Gregorien'
+      else
+        write (*,*) 'ERROR : unknown calendar ', calend
+        stop 'calend should be 360d,earth_365d,earth_366d,gregorian'
+      endif
+!-----------------------------------------------------------------------
+!
+!c Date :
+!      La date est supposee donnee sous la forme [annee, numero du jour dans 
+!      l annee] ; l heure est donnee dans time_ini, lu dans lmdz1d.def.
+!      On appelle ymds2ju pour convertir [annee, jour] en [jour Julien].
+!      Le numero du jour est dans "day". L heure est traitee separement.
+!      La date complete est dans "daytime" (l'unite est le jour).
+      if (nday>0) then
+         fnday=nday
+      else
+         fnday=-nday/float(day_step)
+      endif
+
+! Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026)
+      IF(forcing_type .EQ. 61) fnday=53100./86400.
+! Special case for amma which lasts less than one day : 64800s !! (MPL 20120216)
+      IF(forcing_type .EQ. 6) fnday=64800./86400.
+      annee_ref = anneeref
+      mois = 1
+      day_ref = dayref
+      heure = 0.
+      itau_dyn = 0
+      itau_phy = 0
+      call ymds2ju(annee_ref,mois,day_ref,heure,day)
+      day_ini = int(day)
+      day_end = day_ini + fnday
+
+      IF (forcing_type .eq.2) THEN
+! Convert the initial date of Toga-Coare to Julian day
+      call ymds2ju                                                          &
+     & (year_ini_toga,mth_ini_toga,day_ini_toga,heure,day_ju_ini_toga)
+
+      ELSEIF (forcing_type .eq.4) THEN
+! Convert the initial date of TWPICE to Julian day
+      call ymds2ju                                                          &
+     & (year_ini_twpi,mth_ini_twpi,day_ini_twpi,heure_ini_twpi              &
+     & ,day_ju_ini_twpi)
+      ELSEIF (forcing_type .eq.6) THEN
+! Convert the initial date of AMMA to Julian day
+      call ymds2ju                                                          &
+     & (year_ini_amma,mth_ini_amma,day_ini_amma,heure_ini_amma              &
+     & ,day_ju_ini_amma)
+
+      ELSEIF (forcing_type .eq.59) THEN
+! Convert the initial date of Sandu case to Julian day
+      call ymds2ju                                                          &
+     &   (year_ini_sandu,mth_ini_sandu,day_ini_sandu,                       &
+     &    time_ini*3600.,day_ju_ini_sandu)
+
+      ELSEIF (forcing_type .eq.60) THEN
+! Convert the initial date of Astex case to Julian day
+      call ymds2ju                                                          &
+     &   (year_ini_astex,mth_ini_astex,day_ini_astex,                        &
+     &    time_ini*3600.,day_ju_ini_astex)
+
+      ELSEIF (forcing_type .eq.61) THEN
+
+! Convert the initial date of Arm_cu case to Julian day
+      call ymds2ju                                                          &
+     & (year_ini_armcu,mth_ini_armcu,day_ini_armcu,heure_ini_armcu          &
+     & ,day_ju_ini_armcu)
+      ENDIF
+
+      daytime = day + time_ini/24. ! 1st day and initial time of the simulation 
+! Print out the actual date of the beginning of the simulation :
+      call ju2ymds(daytime,year_print, month_print,day_print,sec_print)
+      print *,' Time of beginning : ',                                      &
+     &        year_print, month_print, day_print, sec_print
+
+!---------------------------------------------------------------------
+! Initialization of dimensions, geometry and initial state
+!---------------------------------------------------------------------
+      call init_phys_lmdz(1,1,llm,1,(/1/))
+      call suphel
+      call initcomgeomphy
+      call infotrac_init 
+
+      if (nqtot>nqmx) STOP'Augmenter nqmx dans lmdz1d.F'
+      allocate(q(llm,nqtot)) ; q(:,:)=0.
+      allocate(dq(llm,nqtot)) 
+      allocate(dq_dyn(llm,nqtot)) 
+      allocate(d_q_adv(llm,nqtot)) 
+!     allocate(d_th_adv(llm)) 
+
+!
+!   No ozone climatology need be read in this pre-initialization
+!          (phys_state_var_init is called again in physiq)
+      read_climoz = 0
+!
+      call phys_state_var_init(read_climoz)
+
+      if (ngrid.ne.klon) then
+         print*,'stop in inifis'
+         print*,'Probleme de dimensions :'
+         print*,'ngrid = ',ngrid
+         print*,'klon  = ',klon
+         stop
+      endif
+!!!=====================================================================
+!!! Feedback forcing values for Gateaux differentiation (al1)
+!!!=====================================================================
+!!! Surface Planck forcing bracketing call radiation
+!!      surf_Planck = 0.
+!!      surf_Conv   = 0.
+!!      write(*,*) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv
+!!! a mettre dans le lmdz1d.def ou autre
+!!
+!!
+      qsol = qsolinp
+      qsurf = fq_sat(tsurf,psurf/100.)
+      rlat=xlat
+      rlon=xlon
+      day1= day_ini
+      time=daytime-day
+      ts_toga(1)=tsurf ! needed by read_tsurf1d.F
+      rho(1)=psurf/(rd*tsurf*(1.+(rv/rd-1.)*qsurf)) 
+
+!
+!! mpl et jyg le 22/08/2012 : 
+!!  pour que les cas a flux de surface imposes marchent
+      IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN
+       fsens=-wtsurf*rcpd*rho(1)
+       flat=-wqsurf*rlvtt*rho(1)
+       print *,'Flux: ok_flux wtsurf wqsurf',ok_flux_surf,wtsurf,wqsurf
+      ENDIF
+      print*,'Flux sol ',fsens,flat
+!!      ok_flux_surf=.false.
+!!      fsens=-wtsurf*rcpd*rho(1)
+!!      flat=-wqsurf*rlvtt*rho(1)
+!!!!
+
+! Vertical discretization and pressure levels at half and mid levels:
+
+      pa   = 5e4
+!!      preff= 1.01325e5
+      preff = psurf
+      IF (ok_old_disvert) THEN
+        call disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
+        print *,'On utilise disvert0'
+      ELSE
+        call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,          &
+     &                 scaleheight)
+        print *,'On utilise disvert'
+!       Nouvelle version disvert permettant d imposer ap,bp (modif L.Guez) MPL 18092012
+!       Dans ce cas, on lit ap,bp dans le fichier hybrid.txt
+      ENDIF
+      sig_s=presnivs/preff
+      plev =ap+bp*psurf
+      play = 0.5*(plev(1:llm)+plev(2:llm+1))
+!cc      zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles
+
+      IF (forcing_type .eq. 59) THEN
+! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m
+      write(*,*) '***********************'
+      do l = 1, llm
+       write(*,*) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l)
+       if (trouve_700 .and. play(l).le.70000) then
+         llm700=l
+         print *,'llm700,play=',llm700,play(l)/100.
+         trouve_700= .false.
+       endif
+      enddo
+      write(*,*) '***********************'
+      ENDIF
+
+!
+!=====================================================================
+! EVENTUALLY, READ FORCING DATA :
+!=====================================================================
+
+#include "1D_read_forc_cases.h"
+
+      if (forcing_GCM2SCM) then
+        write (*,*) 'forcing_GCM2SCM not yet implemented'
+        stop 'in initialization'
+      endif ! forcing_GCM2SCM
+
+      print*,'mxcalc=',mxcalc
+      print*,'zlay=',zlay(mxcalc)
+      print*,'play=',play(mxcalc)
+
+!Al1 pour SST forced, appellé depuis ocean_forced_noice
+      ts_cur = tsurf ! SST used in read_tsurf1d
+!=====================================================================
+! Initialisation de la physique : 
+!=====================================================================
+
+!  Rq: conf_phys.F90 lit tous les flags de physiq.def; conf_phys appele depuis physiq.F
+!
+! day_step, iphysiq lus dans gcm.def ci-dessus
+! timestep: calcule ci-dessous from rday et day_step
+! ngrid=1
+! llm: defini dans .../modipsl/modeles/LMDZ4/libf/grid/dimension
+! rday: defini dans suphel.F (86400.)
+! day_ini: lu dans run.def (dayref)
+! rlat_rad,rlon-rad: transformes en radian de rlat,rlon lus dans lmdz1d.def (en degres)
+! airefi,zcufi,zcvfi initialises au debut de ce programme
+! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
+      day_step = float(nsplit_phys)*day_step/float(iphysiq)
+      write (*,*) 'Time step divided by nsplit_phys (=',nsplit_phys,')'
+      timestep =rday/day_step
+      dtime_frcg = timestep
+!
+      zcufi=airefi
+      zcvfi=airefi
+!
+      rlat_rad(:)=rlat(:)*rpi/180.
+      rlon_rad(:)=rlon(:)*rpi/180.
+
+      call iniphysiq(ngrid,llm,rday,day_ini,timestep,                        &
+     &     rlat_rad,rlon_rad,airefi,zcufi,zcvfi,ra,rg,rd,rcpd,(/1/))
+      print*,'apres iniphysiq'
+
+! 2 PARAMETRES QUI DEVRAIENT ETRE LUS DANS run.def MAIS NE LE SONT PAS ICI:
+      co2_ppm= 330.0
+      solaire=1370.0
+
+! Ecriture du startphy avant le premier appel a la physique.
+! On le met juste avant pour avoir acces a tous les champs
+! NB: les clesphy0 seront remplies dans phyredem d'apres les flags lus dans gcm.def
+
+      if (ok_writedem) then
+
+!--------------------------------------------------------------------------
+! pbl_surface_init (called here) and pbl_surface_final (called by phyredem)
+! need : qsol fder snow qsurf evap rugos agesno ftsoil
+!--------------------------------------------------------------------------
+
+        type_ocean = "force"
+        run_off_lic_0(1) = restart_runoff 
+        call fonte_neige_init(run_off_lic_0)
+
+        fder=0.
+        snsrf(1,:)=0.        ! couverture de neige des sous surface
+        qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
+        evap=0.
+        frugs(1,:)=rugos     ! couverture de neige des sous surface
+        agesno  = xagesno
+        tsoil(:,:,:)=tsurf
+!------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012)
+!       tsoil(1,1,1)=299.18
+!       tsoil(1,2,1)=300.08
+!       tsoil(1,3,1)=301.88
+!       tsoil(1,4,1)=305.48
+!       tsoil(1,5,1)=308.00
+!       tsoil(1,6,1)=308.00
+!       tsoil(1,7,1)=308.00
+!       tsoil(1,8,1)=308.00
+!       tsoil(1,9,1)=308.00
+!       tsoil(1,10,1)=308.00
+!       tsoil(1,11,1)=308.00
+!-----------------------------------------------------------------------
+        call pbl_surface_init(qsol, fder, snsrf, qsurfsrf,                  &
+     &                                    evap, frugs, agesno, tsoil)
+
+!------------------ prepare limit conditions for limit.nc -----------------
+!--   Ocean force
+
+        print*,'avant phyredem'
+        pctsrf(1,:)=0.
+        if (nat_surf.eq.0.) then
+          pctsrf(1,is_oce)=1.
+          pctsrf(1,is_ter)=0.
+        else
+          pctsrf(1,is_oce)=0.
+          pctsrf(1,is_ter)=1.
+        end if
+
+        print*,'nat_surf,pctsrf(1,is_oce),pctsrf(1,is_ter)',nat_surf         &
+     &        ,pctsrf(1,is_oce),pctsrf(1,is_ter)
+
+        zmasq=pctsrf(1,is_ter)+pctsrf(1,is_lic)
+        zpic = zpicinp
+        ftsol=tsurf
+        falb1 = albedo                            
+        falb2 = albedo                            
+        rugoro=rugos 
+        t_ancien(1,:)=temp(:)
+        q_ancien(1,:)=q(:,1)
+        pbl_tke=1.e-8
+
+        rain_fall=0.
+        snow_fall=0.
+        solsw=0.
+        sollw=0.
+        radsol=0.
+        rnebcon=0.
+        ratqs=0.
+        clwcon=0.
+        zmea=0.
+        zstd=0.
+        zsig=0.
+        zgam=0.
+        zval=0. 
+        zthe=0. 
+        sig1=0.
+        w01=0.
+        u_ancien(1,:)=u(:)
+        v_ancien(1,:)=v(:)
+ 
+!------------------------------------------------------------------------
+! Make file containing restart for the physics (startphy.nc)
+!
+! NB: List of the variables to be written by phyredem (via put_field):
+! rlon,rlat,zmasq,pctsrf(:,is_ter),pctsrf(:,is_lic),pctsrf(:,is_oce)
+! pctsrf(:,is_sic),ftsol(:,nsrf),tsoil(:,isoil,nsrf),qsurf(:,nsrf)
+! qsol,falb1(:,nsrf),falb2(:,nsrf),evap(:,nsrf),snow(:,nsrf)
+! radsol,solsw,sollw,fder,rain_fall,snow_fall,frugs(:,nsrf)
+! agesno(:,nsrf),zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro
+! t_ancien,q_ancien,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
+! run_off_lic_0,pbl_tke(:,1:klev,nsrf),zmax0,f0,sig1,w01
+! wake_deltat,wake_deltaq,wake_s,wake_cstar,wake_fip
+!------------------------------------------------------------------------
+!Al1 =============== restart option ==========================
+        if (.not.restart) then
+          call phyredem ("startphy.nc")
+        else
+! (desallocations)
+        print*,'callin surf final'
+          call pbl_surface_final(qsol, fder, snsrf, qsurfsrf,               &
+     &                                    evap, frugs, agesno, tsoil)
+        print*,'after surf final'
+          CALL fonte_neige_final(run_off_lic_0)
+        endif
+
+        ok_writedem=.false.
+        print*,'apres phyredem'
+
+      endif ! ok_writedem
+      
+!------------------------------------------------------------------------
+! Make file containing boundary conditions (limit.nc) **Al1->restartdyn***
+! --------------------------------------------------
+! NB: List of the variables to be written in limit.nc 
+!     (by writelim.F, subroutine of 1DUTILS.h):
+!        phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
+!        phy_fter,phy_foce,phy_flic,phy_fsic)
+!------------------------------------------------------------------------
+      do i=1,yd
+        phy_nat(i)  = nat_surf
+        phy_alb(i)  = albedo
+        phy_sst(i)  = tsurf ! read_tsurf1d will be used instead
+        phy_rug(i)  = rugos
+        phy_fter(i) = pctsrf(1,is_ter)
+        phy_foce(i) = pctsrf(1,is_oce)
+        phy_fsic(i) = pctsrf(1,is_sic)
+        phy_flic(i) = pctsrf(1,is_lic)
+      enddo
+
+! fabrication de limit.nc
+      call writelim (1,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,             &
+     &               phy_ice,phy_fter,phy_foce,phy_flic,phy_fsic)
+
+
+      call phys_state_var_end
+!Al1
+      if (restart) then
+        print*,'call to restart dyn 1d'
+        Call dyn1deta0("start1dyn.nc",plev,play,phi,phis,presnivs,          &
+     &              u,v,temp,q,omega2)
+
+       print*,'fnday,annee_ref,day_ref,day_ini',                            &
+     &     fnday,annee_ref,day_ref,day_ini
+!**      call ymds2ju(annee_ref,mois,day_ini,heure,day)
+       day = day_ini
+       day_end = day_ini + nday
+       daytime = day + time_ini/24. ! 1st day and initial time of the simulation 
+
+! Print out the actual date of the beginning of the simulation :
+       call ju2ymds(daytime, an, mois, jour, heure)
+       print *,' Time of beginning : y m d h',an, mois,jour,heure/3600.
+
+       day = int(daytime)
+       time=daytime-day
+  
+       print*,'****** intialised fields from restart1dyn *******'
+       print*,'plev,play,phi,phis,presnivs,u,v,temp,q,omega2'
+       print*,'temp(1),q(1,1),u(1),v(1),plev(1),phis :'
+       print*,temp(1),q(1,1),u(1),v(1),plev(1),phis
+! raz for safety
+       do l=1,llm
+         dq_dyn(l,1) = 0.
+       enddo
+      endif
+!Al1 ================  end restart =================================
+      IF (ecrit_slab_oc.eq.1) then
+         open(97,file='div_slab.dat',STATUS='UNKNOWN')
+       elseif (ecrit_slab_oc.eq.0) then
+         open(97,file='div_slab.dat',STATUS='OLD')
+       endif
+!=====================================================================
+! START OF THE TEMPORAL LOOP :
+!=====================================================================
+           
+      do while(it.le.nint(fnday*day_step))
+
+       if (prt_level.ge.1) then
+         print*,'XXXXXXXXXXXXXXXXXXX ITAP,day,time=',                       &
+     &                                it,day,time,nint(fnday*day_step)
+         print*,'PAS DE TEMPS ',timestep
+       endif
+!Al1 demande de restartphy.nc
+       if (it.eq.nint(fnday*day_step)) lastcall=.True.
+
+!---------------------------------------------------------------------
+! Interpolation of forcings in time and onto model levels
+!---------------------------------------------------------------------
+
+#include "1D_interp_cases.h"
+
+      if (forcing_GCM2SCM) then
+        write (*,*) 'forcing_GCM2SCM not yet implemented'
+        stop 'in time loop'
+      endif ! forcing_GCM2SCM
+
+!---------------------------------------------------------------------
+!  Geopotential :
+!---------------------------------------------------------------------
+
+        phi(1)=RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
+        do l = 1, llm-1
+          phi(l+1)=phi(l)+RD*(temp(l)+temp(l+1))*                           &
+     &    (play(l)-play(l+1))/(play(l)+play(l+1))
+        enddo
+
+!---------------------------------------------------------------------
+! Listing output for debug prt_level>=1
+!---------------------------------------------------------------------
+       if (prt_level>=1) then
+         print *,' avant physiq : -------- day time ',day,time
+         write(*,*) 'firstcall,lastcall,phis',                               &
+     &               firstcall,lastcall,phis
+         write(*,'(a10,2a4,4a13)') 'BEFOR1 IT=','it','l',                   &
+     &        'presniv','plev','play','phi'
+         write(*,'(a10,2i4,4f13.2)') ('BEFOR1 IT= ',it,l,                   &
+     &         presnivs(l),plev(l),play(l),phi(l),l=1,llm)
+         write(*,'(a11,2a4,a11,6a8)') 'BEFOR2','it','l',                    &
+     &         'presniv','u','v','temp','q1','q2','omega2'
+         write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ',it,l,         &
+     &   presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
+       endif
+
+!---------------------------------------------------------------------
+!   Call physiq :
+!---------------------------------------------------------------------
+
+        call physiq(ngrid,llm,                                              &
+     &              firstcall,lastcall,                                     &
+     &              day,time,timestep,                                      &
+     &              plev,play,phi,phis,presnivs,clesphy0,                   &
+     &              u,v,temp,q,omega2,                                      &
+     &              du_phys,dv_phys,dt_phys,dq,dpsrf,                        &
+     &              dudyn,PVteta)
+        firstcall=.false.
+
+!---------------------------------------------------------------------
+! Listing output for debug prt_level>=1
+!---------------------------------------------------------------------
+        if (prt_level>=1) then
+          write(*,'(a11,2a4,4a13)') 'AFTER1 IT=','it','l',                  &
+     &        'presniv','plev','play','phi'
+          write(*,'(a11,2i4,4f13.2)') ('AFTER1 it= ',it,l,                  &
+     &    presnivs(l),plev(l),play(l),phi(l),l=1,llm)
+          write(*,'(a11,2a4,a11,6a8)') 'AFTER2','it','l',                   &
+     &         'presniv','u','v','temp','q1','q2','omega2'
+          write(*,'(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ',it,l,       &
+     &    presnivs(l),u(l),v(l),temp(l),q(l,1),q(l,2),omega2(l),l=1,llm)
+          write(*,'(a11,2a4,a11,5a8)') 'AFTER3','it','l',                   &
+     &         'presniv','du_phys','dv_phys','dt_phys','dq1','dq2'    
+           write(*,'(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ',it,l,            &
+     &      presnivs(l),86400*du_phys(l),86400*dv_phys(l),                   &
+     &       86400*dt_phys(l),86400*dq(l,1),dq(l,2),l=1,llm)
+          write(*,*) 'dpsrf',dpsrf
+        endif
+!---------------------------------------------------------------------
+!   Add physical tendencies :
+!---------------------------------------------------------------------
+
+       fcoriolis=2.*sin(rpi*xlat/180.)*romega
+       if (forcing_radconv .or. forcing_fire) then
+         fcoriolis=0.0
+         dt_cooling=0.0
+         d_th_adv=0.0
+         d_q_adv=0.0
+       endif
+
+       if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice            &
+     &    .or.forcing_amma) then
+         fcoriolis=0.0 ; ug=0. ; vg=0.
+       endif
+         if(forcing_rico) then
+          dt_cooling=0.
+        endif
+
+      print*, 'fcoriolis ', fcoriolis, xlat
+
+        du_age(1:mxcalc)=fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
+       dv_age(1:mxcalc)=-fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
+
+!!!!!!!!!!!!!!!!!!!!!!!!
+! Geostrophic wind
+!!!!!!!!!!!!!!!!!!!!!!!!
+       sfdt = sin(0.5*fcoriolis*timestep)
+       cfdt = cos(0.5*fcoriolis*timestep)
+!
+        du_age(1:mxcalc)= -2.*sfdt/timestep*                                &
+     &          (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -                          &
+     &           cfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
+!!     : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
+!
+       dv_age(1:mxcalc)= -2.*sfdt/timestep*                                 &
+     &          (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) +                           &
+     &           sfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
+!!     : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!         call  writefield_phy('dv_age' ,dv_age,llm)
+!         call  writefield_phy('du_age' ,du_age,llm)
+!         call  writefield_phy('du_phys' ,du_phys,llm)
+!         call  writefield_phy('u_tend' ,u,llm)
+!         call  writefield_phy('u_g' ,ug,llm)
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! Increment state variables
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! pour les cas sandu et astex, on reclacule u,v,q,temp et teta dans 1D_nudge_sandu_astex.h
+! au dessus de 700hpa, on relaxe vers les profils initiaux
+      if (forcing_sandu .OR. forcing_astex) then
+#include "1D_nudge_sandu_astex.h"
+      else
+        u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &
+     &              du_phys(1:mxcalc)                                       &
+     &             +du_age(1:mxcalc) )            
+        v(1:mxcalc)=v(1:mxcalc) + timestep*(                                 &
+     &              dv_phys(1:mxcalc)                                       &
+     &             +dv_age(1:mxcalc) )
+        q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(                              &
+     &                dq(1:mxcalc,:)                                        &
+     &               +d_q_adv(1:mxcalc,:) )
+
+        if (prt_level.ge.1) then
+          print *,                                                          &
+     &    'physiq-> temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1) ',         &
+     &              temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1)
+           print*,du_phys
+           print*, v
+           print*, vg
+        endif
+
+        temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
+     &              dt_phys(1:mxcalc)                                       &
+     &             +d_th_adv(1:mxcalc)                                      &
+     &             +dt_cooling(1:mxcalc))  ! Taux de chauffage ou refroid.
+
+      endif  ! forcing_sandu or forcing_astex
+
+        teta=temp*(pzero/play)**rkappa
+!
+!---------------------------------------------------------------------
+!   Nudge soil temperature if requested
+!---------------------------------------------------------------------
+
+      IF (nudge_tsoil) THEN
+       ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:)                     &
+     &  -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
+      ENDIF
+
+!---------------------------------------------------------------------
+!   Add large-scale tendencies (advection, etc) :
+!---------------------------------------------------------------------
+
+!cc nrlmd
+!cc        tmpvar=teta
+!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
+!cc
+!cc        teta(1:mxcalc)=tmpvar(1:mxcalc)
+!cc        tmpvar(:)=q(:,1)
+!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
+!cc        q(1:mxcalc,1)=tmpvar(1:mxcalc)
+!cc        tmpvar(:)=q(:,2)
+!cc        call advect_vert(llm,omega,timestep,tmpvar,plev)
+!cc        q(1:mxcalc,2)=tmpvar(1:mxcalc)
+
+!---------------------------------------------------------------------
+!   Air temperature :
+!---------------------------------------------------------------------        
+        if (lastcall) then
+          print*,'Pas de temps final ',it
+          call ju2ymds(daytime, an, mois, jour, heure)
+          print*,'a la date : a m j h',an, mois, jour ,heure/3600.
+        endif
+
+!  incremente day time
+!        print*,'daytime bef',daytime,1./day_step
+        daytime = daytime+1./day_step
+!Al1dbg
+        day = int(daytime+0.1/day_step)
+!        time = max(daytime-day,0.0)
+!Al1&jyg: correction de bug
+!cc        time = real(mod(it,day_step))/day_step
+        time = time_ini/24.+real(mod(it,day_step))/day_step
+!        print*,'daytime nxt time',daytime,time
+        it=it+1
+
+      enddo
+
+!Al1
+      if (ecrit_slab_oc.ne.-1) close(97)
+
+!Al1 Call to 1D equivalent of dynredem (an,mois,jour,heure ?)
+! -------------------------------------
+       call dyn1dredem("restart1dyn.nc",                                    &
+     &              plev,play,phi,phis,presnivs,                            &
+     &              u,v,temp,q,omega2)
+
+        CALL abort_gcm ('lmdz1d   ','The End  ',0)
+
+      end
+
+#include "1DUTILS.h"
+#include "1Dconv.h"
+
+