Changeset 127

Timestamp:

May 24, 2011, 1:26:29 PM (14 years ago)

Author:

emillour

Message:

Ehouarn: suite de l'implementation de la discretisation verticale,
avec quelques mises a jour pour concorder avec la version terrestre.
-> Finalement, on met un flag "disvert_type" pour fixer la discretisation

disvert_type==1 (defaut si planet_type=="earth") pour cas terrestre
disvert_type==2 (defaut si planet_type!="earth") pour cas planeto (z2sig.def)

-> au passage, pour rester en phase avec modele terrestre on renomme

disvert_terre en disvert (le disvert "alternatif" demeure 'disvert_noterre')

Location:

trunk

Files:

• chantiers/commit_importants.log (modified) (2 diffs)
• chantiers/meschantiers-Ehouarn.txt (modified) (2 diffs)
• libf/dyn3d/comvert.h (modified) (1 diff)
• libf/dyn3d/disvert.F90 (moved) (moved from trunk/libf/dyn3d/disvert_terre.F90) (8 diffs)
• libf/dyn3d/disvert_noterre.F (modified) (8 diffs)
• libf/dyn3d/etat0_netcdf.F90 (modified) (5 diffs)
• libf/dyn3d/exner_hyb.F (modified) (2 diffs)
• libf/dyn3d/exner_milieu.F (modified) (2 diffs)
• libf/dyn3d/guide_mod.F90 (modified) (1 diff)
• libf/dyn3d/iniacademic.F90 (modified) (7 diffs)
• libf/dyn3d/iniconst.F (modified) (3 diffs)
• libf/dyn3d/leapfrog.F (modified) (5 diffs)
• libf/dyn3d/logic.h (modified) (3 diffs)
• libf/dyn3dpar/comvert.h (modified) (1 diff)
• libf/dyn3dpar/disvert.F90 (moved) (moved from trunk/libf/dyn3dpar/disvert_terre.F90) (8 diffs)
• libf/dyn3dpar/etat0_netcdf.F90 (modified) (5 diffs)
• libf/dyn3dpar/exner_hyb.F (modified) (2 diffs)
• libf/dyn3dpar/exner_hyb_p.F (modified) (2 diffs)
• libf/dyn3dpar/exner_milieu.F (modified) (2 diffs)
• libf/dyn3dpar/exner_milieu_p.F (modified) (2 diffs)
• libf/dyn3dpar/guide_p_mod.F90 (modified) (60 diffs)
• libf/dyn3dpar/iniacademic.F90 (modified) (7 diffs)
• libf/dyn3dpar/iniconst.F (modified) (3 diffs)
• libf/dyn3dpar/leapfrog_p.F (modified) (3 diffs)
• libf/dyn3dpar/logic.h (modified) (3 diffs)

trunk/chantiers/commit_importants.log

@@ -882 +882 @@
 
 *********************
-**** commit_v123 ****
+**** commit_v124 ****
 *********************
 
@@ -892 +892 @@
    pour atteindre la concordance entre versions séq. et //.
    
+*********************
+**** commit_v126 ****
+*********************
+
+Ehouarn: suite de l'implémentation de la discrétisation verticale, quelques
+         mises à jour pour concorder avec la version terrestre.
+-> Finalement, on met un flag "disvert_type" pour fixer la discrétisation
+   disvert_type==1 (défaut si planet_type=="earth") pour cas terrestre
+   disvert_type==2 (défaut si planet_type!="earth") pour cas planéto (z2sig.def)
+-> au passage, pour rester en phase avec modèle terrestre on renomme
+   disvert_terre en disvert (le disvert "alternatif" demeure 'disvert_noterre')

trunk/chantiers/meschantiers-Ehouarn.txt

@@ -1 +1 @@
 >>> choses à faire:
 
-- En priorité: faire qu'on puisse compiler un/des cas test (sans physique) pour pouvoir tester les modifs au fur et à mesure. OK, c'est fait, on peut compiler/tourner sans physique
+- En priorité: faire qu'on puisse compiler un/des cas test (sans physique)
+  pour pouvoir tester les modifs au fur et à mesure. OK, c'est fait, on
+  peut compiler/tourner sans physique
 
 - Uniformiser les mises à jour dyn séq. et // (commencé avec la rev. 8)
@@ -7 +9 @@
 - Peut-être revoir l'interface dynamique/physique ?
 
-- Attention à la discrétisation verticale ... 
+- Attention à la discrétisation verticale ... ( débuté avec la r109 par Seb
+  puis mes modifs en r124 et finalement r126) donc a priori
+  OK maintenant
 

trunk/libf/dyn3d/comvert.h

@@ -1 +1 @@
 !
-! $Id: comvert.h 1279 2009-12-10 09:02:56Z fairhead $
+! $Id: comvert.h 1520 2011-05-23 11:37:09Z emillour $
 !
 !-----------------------------------------------------------------------
 !   INCLUDE 'comvert.h'
 
-      COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),      &
+      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
      &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
      &               aps(llm),bps(llm),scaleheight
 
-      REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig,aps,bps
+      common/comverti/disvert_type
+
+      real ap     ! hybrid pressure contribution at interlayers
+      real bp     ! hybrid sigma contribution at interlayer
+      real presnivs ! (reference) pressure at mid-layers
+      real dpres
+      real pa     ! reference pressure (Pa) at which hybrid coordinates
+                  ! become purely pressure
+      real preff  ! reference surface pressure (Pa)
+      real nivsigs
+      real nivsig
+      real aps    ! hybrid pressure contribution at mid-layers
+      real bps    ! hybrid sigma contribution at mid-layers
       real scaleheight ! atmospheric (reference) scale height (km)
 
+      integer disvert_type ! type of vertical discretization:
+                           ! 1: Earth (default for planet_type==earth),
+                           !     automatic generation
+                           ! 2: Planets (default for planet_type!=earth),
+                           !     using 'z2sig.def' (or 'esasig.def) file
+
  !-----------------------------------------------------------------------

trunk/libf/dyn3d/disvert.F90

@@ -1 @@
-! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $
+! $Id: disvert.F90 1520 2011-05-23 11:37:09Z emillour $
 
-SUBROUTINE disvert_terre(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
+SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
 
   ! Auteur : P. Le Van
@@ -17 +17 @@
   ! ds(l) : distance entre les couches l et l-1 en coord.s
 
-  REAL pa, preff
-  REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
-  REAL presnivs(llm)
+  real,intent(in) :: pa, preff
+  real,intent(out) :: ap(llmp1), bp(llmp1)
+  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
+  real,intent(out) :: presnivs(llm)
+  real,intent(out) :: scaleheight
 
   REAL sig(llm+1), dsig(llm)
@@ -26 +28 @@
   INTEGER l
   REAL dsigmin
-  REAL alpha, beta, deltaz, h
+  REAL alpha, beta, deltaz,scaleheight
   INTEGER iostat
-  REAL pi, x
+  REAL x
+  character(len=*),parameter :: modname="disvert"
 
   !-----------------------------------------------------------------------
 
-  pi = 2 * ASIN(1.)
+  ! default scaleheight is 8km for earth
+  scaleheight=8.
 
   OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
@@ -38 +42 @@
   IF (iostat == 0) THEN
      ! cas 1 on lit les options dans sigma.def:
-     READ(99, *) h ! hauteur d'echelle 8.
+     READ(99, *) scaleheight ! hauteur d'echelle 8.
      READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
      READ(99, *) beta ! facteur d'acroissement en haut 1.3
@@ -44 +48 @@
      READ(99, *) k1 ! nombre de couches dans la transition haute
      CLOSE(99)
-     alpha=deltaz/(llm*h)
-     write(lunout, *)'disvert: h, alpha, k0, k1, beta'
+     alpha=deltaz/(llm*scaleheight)
+     write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
+                               scaleheight, alpha, k0, k1, beta
 
      alpha=deltaz/tanh(1./k0)*2.
@@ -51 +56 @@
      sig(1)=1.
      do l=1, llm
-        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
-             *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
-        zk=-h*log(sig(l+1))
+        sig(l+1)=(cosh(l/k0))**(-alpha*k0/scaleheight) &
+             *exp(-alpha/scaleheight*tanh((llm-k1)/k0) &
+                  *beta**(l-(llm-k1))/log(beta))
+        zk=-scaleheight*log(sig(l+1))
 
         dzk1=alpha*tanh(l/k0)
@@ -73 +79 @@
            dsigmin=1.
         else
-           WRITE(LUNOUT,*)'disvert: ATTENTION discretisation z a ajuster'
+           write(lunout,*) trim(modname), &
+           ' ATTENTION discretisation z a ajuster'
            dsigmin=1.
         endif
-        WRITE(LUNOUT,*) 'disvert: Discretisation verticale DSIGMIN=',dsigmin
+        write(lunout,*) trim(modname), &
+        ' Discretisation verticale DSIGMIN=',dsigmin
      endif
 
@@ -119 +127 @@
   ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
 
-  write(lunout, *)'disvert: BP '
+  write(lunout, *)  trim(modname),': BP '
   write(lunout, *) bp
-  write(lunout, *)'disvert: AP '
+  write(lunout, *)  trim(modname),': AP '
   write(lunout, *) ap
 
   write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
   write(lunout, *)'couches calcules pour une pression de surface =', preff
-  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
-  write(lunout, *)'8km'
+  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de '
+  write(lunout, *) scaleheight,' km'
   DO l = 1, llm
      dpres(l) = bp(l) - bp(l+1)
      presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
      write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
-          log(preff/presnivs(l))*8. &
-          , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
+          log(preff/presnivs(l))*scaleheight &
+          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
           max(ap(l+1)+bp(l+1)*preff, 1.e-10))
   ENDDO
 
-  write(lunout, *) 'disvert: PRESNIVS '
+  write(lunout, *) trim(modname),': PRESNIVS '
   write(lunout, *) presnivs
 
-END SUBROUTINE disvert_terre
+END SUBROUTINE disvert

trunk/libf/dyn3d/disvert_noterre.F

@@ +1 @@
+! $Id: $
       SUBROUTINE disvert_noterre
 
@@ -22 +23 @@
 c
 c=======================================================================
-c    Discretisation verticale en coordonnée hybride 
+c    Discretisation verticale en coordonnée hybride (ou sigma)
 c
 c=======================================================================
@@ -49 +50 @@
       integer iz 
       real z, ps,p
+      character(len=*),parameter :: modname="disvert_noterre"
 
 c
@@ -54 +56 @@
 c
 ! Initializations:
-      pi=2.*ASIN(1.)
+!      pi=2.*ASIN(1.) ! already done in iniconst
       
       hybrid=.true. ! default value for hybrid (ie: use hybrid coordinates)
       CALL getin('hybrid',hybrid)
-      write(lunout,*)'disvert_noterre: hybrid=',hybrid
+      write(lunout,*) trim(modname),': hybrid=',hybrid
 
 ! Ouverture possible de fichiers typiquement E.T.
@@ -156 +158 @@
 
       DO l=1,llm
-        nivsigs(l) = FLOAT(l)
+        nivsigs(l) = REAL(l)
       ENDDO
 
       DO l=1,llmp1
-        nivsig(l)= FLOAT(l)
+        nivsig(l)= REAL(l)
       ENDDO
 
@@ -199 +201 @@
       bp(llmp1) =   0.
 
-      write(lunout,*)' BP '
+      write(lunout,*) trim(modname),': BP '
       write(lunout,*)  bp
-      write(lunout,*)' AP '
+      write(lunout,*) trim(modname),': AP '
       write(lunout,*)  ap
 
@@ -225 +227 @@
       end if
 
-      write(lunout,*)' BPs '
+      write(lunout,*) trim(modname),': BPs '
       write(lunout,*)  bps
-      write(lunout,*)' APs'
+      write(lunout,*) trim(modname),': APs'
       write(lunout,*)  aps
 
@@ -235 +237 @@
       ENDDO
 
-      write(lunout,*)' PRESNIVS' 
+      write(lunout,*)trim(modname),' : PRESNIVS' 
       write(lunout,*)presnivs 
       write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ',

trunk/libf/dyn3d/etat0_netcdf.F90

@@ -1 +1 @@
 !
-! $Id: etat0_netcdf.F90 1486 2011-02-11 12:07:39Z fairhead $
+! $Id: etat0_netcdf.F90 1520 2011-05-23 11:37:09Z emillour $
 !
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE etat0_netcdf(ib, masque, letat0)
+SUBROUTINE etat0_netcdf(ib, masque, phis, letat0)
 !
 !-------------------------------------------------------------------------------
@@ -37 +37 @@
   LOGICAL,                    INTENT(IN)    :: ib     ! barycentric interpolat.
   REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
+  REAL, DIMENSION(iip1,jjp1), INTENT(OUT)   :: phis   ! geopotentiel au sol
   LOGICAL,                    INTENT(IN)    :: letat0 ! F: masque only required
 #ifndef CPP_EARTH
@@ -51 +52 @@
   REAL,    DIMENSION(klon)                 :: tsol, qsol
   REAL,    DIMENSION(klon)                 :: sn, rugmer, run_off_lic_0
-  REAL,    DIMENSION(iip1,jjp1)            :: orog, rugo, psol, phis
+  REAL,    DIMENSION(iip1,jjp1)            :: orog, rugo, psol
   REAL,    DIMENSION(iip1,jjp1,llm+1)      :: p3d
   REAL,    DIMENSION(iip1,jjp1,llm)        :: uvent, t3d, tpot, qsat, qd
@@ -138 +139 @@
                    flag_aerosol, new_aod,                               &
                    bl95_b0, bl95_b1,                                    &
-                   iflag_thermals,nsplit_thermals,tau_thermals,         &
-                   iflag_thermals_ed,iflag_thermals_optflux,            &
-                   iflag_coupl,iflag_clos,iflag_wake, read_climoz,      &
+                   read_climoz,                                         &
                    alp_offset)
 
@@ -252 +251 @@
 !*******************************************************************************
   CALL pression(ip1jmp1, ap, bp, psol, p3d)
-  CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
+  if (disvert_type.eq.1) then
+    CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
+  else ! we assume that we are in the disvert_type==2 case
+    CALL exner_milieu(ip1jmp1,psol,p3d,beta,pks,pk,y)
+  endif
   pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa)
 !  WRITE(lunout,*) 'P3D :', p3d(10,20,:)

trunk/libf/dyn3d/exner_hyb.F

@@ -51 +51 @@
       REAL SSUM
 c
+      logical,save :: firstcall=.true.
+      character(len=*),parameter :: modname="exner_hyb"
+      
+      ! Sanity check
+      if (firstcall) then
+        ! check that vertical discretization is compatible
+        ! with this routine
+        if (disvert_type.ne.1) then
+          call abort_gcm(modname,
+     &     "this routine should only be called if disvert_type==1",42)
+        endif
+        
+        ! sanity checks for Shallow Water case (1 vertical layer)
+        if (llm.eq.1) then
+          if (kappa.ne.1) then
+            call abort_gcm(modname,
+     &      "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+            call abort_gcm(modname,
+     &      "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+        endif ! of if (llm.eq.1)
+
+        firstcall=.false.
+      endif ! of if (firstcall)
 
       if (llm.eq.1) then
-        ! Specific behaviour for Shallow Water (1 vertical layer) case
-      
-        ! Sanity checks
-        if (kappa.ne.1) then
-          call abort_gcm("exner_hyb",
-     &    "kappa!=1 , but running in Shallow Water mode!!",42)
-        endif
-        if (cpp.ne.r) then
-        call abort_gcm("exner_hyb",
-     &    "cpp!=r , but running in Shallow Water mode!!",42)
-        endif
         
         ! Compute pks(:),pk(:),pkf(:)
@@ -77 +92 @@
         ! our work is done, exit routine
         return
+
       endif ! of if (llm.eq.1)
 
+!!!! General case:
      
       unpl2k    = 1.+ 2.* kappa

trunk/libf/dyn3d/exner_milieu.F

@@ -48 +48 @@
       REAL SSUM
       EXTERNAL SSUM
+      logical,save :: firstcall=.true.
+      character(len=*),parameter :: modname="exner_milieu"
 
-      if (llm.eq.1) then
-        ! Specific behaviour for Shallow Water (1 vertical layer) case
-      
-        ! Sanity checks
-        if (kappa.ne.1) then
-          call abort_gcm("exner_hyb",
-     &    "kappa!=1 , but running in Shallow Water mode!!",42)
-        endif
-        if (cpp.ne.r) then
-        call abort_gcm("exner_hyb",
-     &    "cpp!=r , but running in Shallow Water mode!!",42)
+      ! Sanity check
+      if (firstcall) then
+        ! check that vertical discretization is compatible
+        ! with this routine
+        if (disvert_type.ne.2) then
+          call abort_gcm(modname,
+     &     "this routine should only be called if disvert_type==2",42)
         endif
         
+        ! sanity checks for Shallow Water case (1 vertical layer)
+        if (llm.eq.1) then
+          if (kappa.ne.1) then
+            call abort_gcm(modname,
+     &      "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+            call abort_gcm(modname,
+     &      "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+        endif ! of if (llm.eq.1)
+
+        firstcall=.false.
+      endif ! of if (firstcall)
+
+!!!! Specific behaviour for Shallow Water (1 vertical layer) case:
+      if (llm.eq.1) then
+      
         ! Compute pks(:),pk(:),pkf(:)
         
@@ -74 +90 @@
         ! our work is done, exit routine
         return
+
       endif ! of if (llm.eq.1)
 
-      
+!!!! General case:
+
 c     -------------
 c     Calcul de pks

trunk/libf/dyn3d/guide_mod.F90

@@ -644 +644 @@
 ! -----------------------------------------------------------------
     CALL pression( ip1jmp1, ap, bp, psi, p )
-    CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
-
+    if (disvert_type==1) then
+      CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
+    else ! we assume that we are in the disvert_type==2 case
+      CALL exner_milieu(ip1jmp1,psi,p,beta,pks,pk,pkf)
+    endif
 !    ....  Calcul de pls , pression au milieu des couches ,en Pascals
     unskap=1./kappa

trunk/libf/dyn3d/iniacademic.F90

@@ -1 +1 @@
 !
-! $Id: iniacademic.F90 1474 2011-01-14 11:04:45Z lguez $
+! $Id: iniacademic.F90 1520 2011-05-23 11:37:09Z emillour $
 !
 SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
@@ -71 +71 @@
 
   REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
+  
+  character(len=*),parameter :: modname="iniacademic"
+  character(len=80) :: abort_message
 
   !-----------------------------------------------------------------------
@@ -76 +79 @@
   ! --------------------------------------
   !
-  
-  write(lunout,*) 'Iniacademic'
-  
   ! initialize planet radius, rotation rate,...
   call conf_planete
@@ -136 +136 @@
      teta0=315.     ! mean Teta (S.H. 315K)
      CALL getin('teta0',teta0)
-     write(lunout,*) 'Iniacademic - teta0 ',teta0
-     write(lunout,*) 'Iniacademic - rad ',rad
      ttp=200.       ! Tropopause temperature (S.H. 200K)
      CALL getin('ttp',ttp)
@@ -183 +181 @@
            tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
                 -delt_z*(1.-ddsin*ddsin)*log(zsig)
-           !! Aymeric -- tests particuliers
            if (planet_type=="giant") then
              tetajl(j,l)=teta0+(delt_y*                   &
@@ -207 +204 @@
         enddo
      enddo
-!     write(lunout,*) 'Iniacademic - check',tetajl(:,int(llm/2)),rlatu(:)
 
      ! 3. Initialize fields (if necessary)
@@ -217 +213 @@
 
         CALL pression ( ip1jmp1, ap, bp, ps, p       )
-        if (planet_type=="earth") then
-          CALL exner_hyb(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+        if (disvert_type.eq.1) then
+          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        elseif (disvert_type.eq.2) then
+          call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
         else
-          call exner_milieu(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+          write(abort_message,*) "Wrong value for disvert_type: ", &
+                              disvert_type
+          call abort_gcm(modname,abort_message,0)
         endif
         CALL massdair(p,masse)

trunk/libf/dyn3d/iniconst.F

@@ -1 +1 @@
 !
-! $Id: iniconst.F 1403 2010-07-01 09:02:53Z fairhead $
+! $Id: iniconst.F 1520 2011-05-23 11:37:09Z emillour $
 !
       SUBROUTINE iniconst
 
       USE control_mod
+#ifdef CPP_IOIPSL
+      use IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      use ioipsl_getincom
+#endif
 
       IMPLICIT NONE
@@ -22 +28 @@
 
 
+      character(len=*),parameter :: modname="iniconst"
+      character(len=80) :: abort_message
 c
 c
@@ -49 +57 @@
       r       = cpp * kappa
 
-      write(lunout,*)'iniconst: R  CP  Kappa ',  r , cpp,  kappa
+      write(lunout,*) trim(modname),': R  CP  Kappa ',r,cpp,kappa
 c
 c-----------------------------------------------------------------------
 
-      if (planet_type.eq."earth") then
-       CALL disvert_terre(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
+! vertical discretization: default behavior depends on planet_type flag
+      if (planet_type=="earth") then
+        disvert_type=1
       else
-       CALL disvert_noterre
+        disvert_type=2
       endif
-c
-       RETURN
-       END
+      ! but user can also specify using one or the other in run.def:
+      call getin('disvert_type',disvert_type)
+      write(lunout,*) trim(modname),': disvert_type=',disvert_type
+      
+      if (disvert_type==1) then
+       ! standard case for Earth (automatic generation of levels)
+       call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,
+     &              scaleheight)
+      else if (disvert_type==2) then
+        ! standard case for planets (levels generated using z2sig.def file)
+        call disvert_noterre
+      else
+        write(abort_message,*) "Wrong value for disvert_type: ",
+     &                        disvert_type
+        call abort_gcm(modname,abort_message,0)
+      endif
+
+      END

trunk/libf/dyn3d/leapfrog.F

@@ -178 +178 @@
 
       character*80 dynhist_file, dynhistave_file
-      character(len=20) :: modname
+      character(len=*),parameter :: modname="leapfrog"
       character*80 abort_message
 
@@ -206 +206 @@
       endif
       itaufinp1 = itaufin +1
-      modname="leapfrog"
       
 c INITIALISATIONS
@@ -238 +237 @@
       dq(:,:,:)=0.
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
-      if (planet_type.eq."earth") then
+      if (disvert_type==1) then
         CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-      else
+      else ! we assume that we are in the disvert_type==2 case
         CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
       endif
@@ -409 +408 @@
 
          CALL pression (  ip1jmp1, ap, bp, ps,  p      )
-         if (planet_type.eq."earth") then
-           CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-         else
+         if (disvert_type==1) then
+           CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
+         else ! we assume that we are in the disvert_type==2 case
            CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
          endif
@@ -531 +530 @@
 
         CALL pression ( ip1jmp1, ap, bp, ps, p                  )
-        if (planet_type.eq."earth") then
+        if (disvert_type==1) then
           CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-        else
+        else ! we assume that we are in the disvert_type==2 case
           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
         endif

trunk/libf/dyn3d/logic.h

@@ -1 +1 @@
 !
-! $Id: logic.h 1319 2010-02-23 21:29:54Z fairhead $
+! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
 !
 !
@@ -11 +11 @@
      &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
      &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
-     &  ,ok_limit,ok_etat0,hybrid
+     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf,hybrid
 
       COMMON/logici/ iflag_phys,iflag_trac
@@ -18 +18 @@
      & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
      &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
-     &  ,ok_limit,ok_etat0
-      logical hybrid ! vertcal coordinate is hybrid if true (sigma otherwise)
+     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf
+      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
+                     ! (only used if disvert_type==2)
 
-      INTEGER iflag_phys,iflag_trac
+      integer iflag_phys,iflag_trac
 !-----------------------------------------------------------------------

trunk/libf/dyn3dpar/comvert.h

@@ -1 +1 @@
 !
-! $Id: comvert.h 1279 2009-12-10 09:02:56Z fairhead $
+! $Id: comvert.h 1520 2011-05-23 11:37:09Z emillour $
 !
 !-----------------------------------------------------------------------
 !   INCLUDE 'comvert.h'
 
-      COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),      &
+      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
      &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
      &               aps(llm),bps(llm),scaleheight
 
-      REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig,aps,bps
+      common/comverti/disvert_type
+
+      real ap     ! hybrid pressure contribution at interlayers
+      real bp     ! hybrid sigma contribution at interlayer
+      real presnivs ! (reference) pressure at mid-layers
+      real dpres
+      real pa     ! reference pressure (Pa) at which hybrid coordinates
+                  ! become purely pressure
+      real preff  ! reference surface pressure (Pa)
+      real nivsigs
+      real nivsig
+      real aps    ! hybrid pressure contribution at mid-layers
+      real bps    ! hybrid sigma contribution at mid-layers
       real scaleheight ! atmospheric (reference) scale height (km)
 
+      integer disvert_type ! type of vertical discretization:
+                           ! 1: Earth (default for planet_type==earth),
+                           !     automatic generation
+                           ! 2: Planets (default for planet_type!=earth),
+                           !     using 'z2sig.def' (or 'esasig.def) file
+
  !-----------------------------------------------------------------------

trunk/libf/dyn3dpar/disvert.F90

@@ -1 @@
-! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $
+! $Id: disvert.F90 1520 2011-05-23 11:37:09Z emillour $
 
-SUBROUTINE disvert_terre(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
+SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
 
   ! Auteur : P. Le Van
@@ -17 +17 @@
   ! ds(l) : distance entre les couches l et l-1 en coord.s
 
-  REAL pa, preff
-  REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
-  REAL presnivs(llm)
+  real,intent(in) :: pa, preff
+  real,intent(out) :: ap(llmp1), bp(llmp1)
+  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
+  real,intent(out) :: presnivs(llm)
+  real,intent(out) :: scaleheight
 
   REAL sig(llm+1), dsig(llm)
@@ -26 +28 @@
   INTEGER l
   REAL dsigmin
-  REAL alpha, beta, deltaz, h
+  REAL alpha, beta, deltaz,scaleheight
   INTEGER iostat
-  REAL pi, x
+  REAL x
+  character(len=*),parameter :: modname="disvert"
 
   !-----------------------------------------------------------------------
 
-  pi = 2 * ASIN(1.)
+  ! default scaleheight is 8km for earth
+  scaleheight=8.
 
   OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
@@ -38 +42 @@
   IF (iostat == 0) THEN
      ! cas 1 on lit les options dans sigma.def:
-     READ(99, *) h ! hauteur d'echelle 8.
+     READ(99, *) scaleheight ! hauteur d'echelle 8.
      READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
      READ(99, *) beta ! facteur d'acroissement en haut 1.3
@@ -44 +48 @@
      READ(99, *) k1 ! nombre de couches dans la transition haute
      CLOSE(99)
-     alpha=deltaz/(llm*h)
-     write(lunout, *)'disvert: h, alpha, k0, k1, beta'
+     alpha=deltaz/(llm*scaleheight)
+     write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
+                               scaleheight, alpha, k0, k1, beta
 
      alpha=deltaz/tanh(1./k0)*2.
@@ -51 +56 @@
      sig(1)=1.
      do l=1, llm
-        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
-             *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
-        zk=-h*log(sig(l+1))
+        sig(l+1)=(cosh(l/k0))**(-alpha*k0/scaleheight) &
+             *exp(-alpha/scaleheight*tanh((llm-k1)/k0) &
+                  *beta**(l-(llm-k1))/log(beta))
+        zk=-scaleheight*log(sig(l+1))
 
         dzk1=alpha*tanh(l/k0)
@@ -73 +79 @@
            dsigmin=1.
         else
-           WRITE(LUNOUT,*)'disvert: ATTENTION discretisation z a ajuster'
+           write(lunout,*) trim(modname), &
+           ' ATTENTION discretisation z a ajuster'
            dsigmin=1.
         endif
-        WRITE(LUNOUT,*) 'disvert: Discretisation verticale DSIGMIN=',dsigmin
+        write(lunout,*) trim(modname), &
+        ' Discretisation verticale DSIGMIN=',dsigmin
      endif
 
@@ -119 +127 @@
   ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
 
-  write(lunout, *)'disvert: BP '
+  write(lunout, *)  trim(modname),': BP '
   write(lunout, *) bp
-  write(lunout, *)'disvert: AP '
+  write(lunout, *)  trim(modname),': AP '
   write(lunout, *) ap
 
   write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
   write(lunout, *)'couches calcules pour une pression de surface =', preff
-  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
-  write(lunout, *)'8km'
+  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de '
+  write(lunout, *) scaleheight,' km'
   DO l = 1, llm
      dpres(l) = bp(l) - bp(l+1)
      presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
      write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
-          log(preff/presnivs(l))*8. &
-          , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
+          log(preff/presnivs(l))*scaleheight &
+          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
           max(ap(l+1)+bp(l+1)*preff, 1.e-10))
   ENDDO
 
-  write(lunout, *) 'disvert: PRESNIVS '
+  write(lunout, *) trim(modname),': PRESNIVS '
   write(lunout, *) presnivs
 
-END SUBROUTINE disvert_terre
+END SUBROUTINE disvert

trunk/libf/dyn3dpar/etat0_netcdf.F90

@@ -1 +1 @@
 !
-! $Id: etat0_netcdf.F90 1486 2011-02-11 12:07:39Z fairhead $
+! $Id: etat0_netcdf.F90 1520 2011-05-23 11:37:09Z emillour $
 !
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE etat0_netcdf(ib, masque, letat0)
+SUBROUTINE etat0_netcdf(ib, masque, phis, letat0)
 !
 !-------------------------------------------------------------------------------
@@ -37 +37 @@
   LOGICAL,                    INTENT(IN)    :: ib     ! barycentric interpolat.
   REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
+  REAL, DIMENSION(iip1,jjp1), INTENT(OUT)   :: phis   ! geopotentiel au sol
   LOGICAL,                    INTENT(IN)    :: letat0 ! F: masque only required
 #ifndef CPP_EARTH
@@ -51 +52 @@
   REAL,    DIMENSION(klon)                 :: tsol, qsol
   REAL,    DIMENSION(klon)                 :: sn, rugmer, run_off_lic_0
-  REAL,    DIMENSION(iip1,jjp1)            :: orog, rugo, psol, phis
+  REAL,    DIMENSION(iip1,jjp1)            :: orog, rugo, psol
   REAL,    DIMENSION(iip1,jjp1,llm+1)      :: p3d
   REAL,    DIMENSION(iip1,jjp1,llm)        :: uvent, t3d, tpot, qsat, qd
@@ -138 +139 @@
                    flag_aerosol, new_aod,                               &
                    bl95_b0, bl95_b1,                                    &
-                   iflag_thermals,nsplit_thermals,tau_thermals,         &
-                   iflag_thermals_ed,iflag_thermals_optflux,            &
-                   iflag_coupl,iflag_clos,iflag_wake, read_climoz,      &
+                   read_climoz,                                         &
                    alp_offset)
 
@@ -252 +251 @@
 !*******************************************************************************
   CALL pression(ip1jmp1, ap, bp, psol, p3d)
-  CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
+  if (disvert_type.eq.1) then
+    CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
+  else ! we assume that we are in the disvert_type==2 case
+    CALL exner_milieu(ip1jmp1,psol,p3d,beta,pks,pk,y)
+  endif
   pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa)
 !  WRITE(lunout,*) 'P3D :', p3d(10,20,:)

trunk/libf/dyn3dpar/exner_hyb.F

@@ -51 +51 @@
       REAL SSUM
 c
+      logical,save :: firstcall=.true.
+      character(len=*),parameter :: modname="exner_hyb"
+      
+      ! Sanity check
+      if (firstcall) then
+        ! check that vertical discretization is compatible
+        ! with this routine
+        if (disvert_type.ne.1) then
+          call abort_gcm(modname,
+     &     "this routine should only be called if disvert_type==1",42)
+        endif
+        
+        ! sanity checks for Shallow Water case (1 vertical layer)
+        if (llm.eq.1) then
+          if (kappa.ne.1) then
+            call abort_gcm(modname,
+     &      "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+            call abort_gcm(modname,
+     &      "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+        endif ! of if (llm.eq.1)
+
+        firstcall=.false.
+      endif ! of if (firstcall)
 
       if (llm.eq.1) then
-        ! Specific behaviour for Shallow Water (1 vertical layer) case
-      
-        ! Sanity checks
-        if (kappa.ne.1) then
-          call abort_gcm("exner_hyb",
-     &    "kappa!=1 , but running in Shallow Water mode!!",42)
-        endif
-        if (cpp.ne.r) then
-        call abort_gcm("exner_hyb",
-     &    "cpp!=r , but running in Shallow Water mode!!",42)
-        endif
         
         ! Compute pks(:),pk(:),pkf(:)
@@ -77 +92 @@
         ! our work is done, exit routine
         return
+
       endif ! of if (llm.eq.1)
 
+!!!! General case:
      
       unpl2k    = 1.+ 2.* kappa

trunk/libf/dyn3dpar/exner_hyb_p.F

@@ -53 +53 @@
       EXTERNAL SSUM
       INTEGER ije,ijb,jje,jjb
-c
-c$OMP BARRIER
-
-      if (llm.eq.1) then
-        ! Specific behaviour for Shallow Water (1 vertical layer) case
-      
-        ! Sanity checks
-        if (kappa.ne.1) then
-          call abort_gcm("exner_hyb",
-     &    "kappa!=1 , but running in Shallow Water mode!!",42)
-        endif
-        if (cpp.ne.r) then
-        call abort_gcm("exner_hyb",
-     &    "cpp!=r , but running in Shallow Water mode!!",42)
+      logical,save :: firstcall=.true.
+!$OMP THREADPRIVATE(firstcall) 
+      character(len=*),parameter :: modname="exner_hyb_p"
+c
+
+      ! Sanity check
+      if (firstcall) then
+        ! check that vertical discretization is compatible
+        ! with this routine
+        if (disvert_type.ne.1) then
+          call abort_gcm(modname,
+     &     "this routine should only be called if disvert_type==1",42)
         endif
         
+        ! sanity checks for Shallow Water case (1 vertical layer)
+        if (llm.eq.1) then
+          if (kappa.ne.1) then
+            call abort_gcm(modname,
+     &      "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+            call abort_gcm(modname,
+     &      "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+        endif ! of if (llm.eq.1)
+
+        firstcall=.false.
+      endif ! of if (firstcall)
+
+c$OMP BARRIER
+
+! Specific behaviour for Shallow Water (1 vertical layer) case
+      if (llm.eq.1) then
+      
         ! Compute pks(:),pk(:),pkf(:)
         ijb=ij_begin
@@ -116 +134 @@
       endif ! of if (llm.eq.1)
 
+!!!! General case:
 
       unpl2k    = 1.+ 2.* kappa

trunk/libf/dyn3dpar/exner_milieu.F

@@ -48 +48 @@
       REAL SSUM
       EXTERNAL SSUM
+      logical,save :: firstcall=.true.
+      character(len=*),parameter :: modname="exner_milieu"
 
-      if (llm.eq.1) then
-        ! Specific behaviour for Shallow Water (1 vertical layer) case
-      
-        ! Sanity checks
-        if (kappa.ne.1) then
-          call abort_gcm("exner_hyb",
-     &    "kappa!=1 , but running in Shallow Water mode!!",42)
-        endif
-        if (cpp.ne.r) then
-        call abort_gcm("exner_hyb",
-     &    "cpp!=r , but running in Shallow Water mode!!",42)
+      ! Sanity check
+      if (firstcall) then
+        ! check that vertical discretization is compatible
+        ! with this routine
+        if (disvert_type.ne.2) then
+          call abort_gcm(modname,
+     &     "this routine should only be called if disvert_type==2",42)
         endif
         
+        ! sanity checks for Shallow Water case (1 vertical layer)
+        if (llm.eq.1) then
+          if (kappa.ne.1) then
+            call abort_gcm(modname,
+     &      "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+            call abort_gcm(modname,
+     &      "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+        endif ! of if (llm.eq.1)
+
+        firstcall=.false.
+      endif ! of if (firstcall)
+
+!!!! Specific behaviour for Shallow Water (1 vertical layer) case:
+      if (llm.eq.1) then
+      
         ! Compute pks(:),pk(:),pkf(:)
         
@@ -74 +90 @@
         ! our work is done, exit routine
         return
+
       endif ! of if (llm.eq.1)
 
-      
+!!!! General case:
+
 c     -------------
 c     Calcul de pks

trunk/libf/dyn3dpar/exner_milieu_p.F

@@ -50 +50 @@
       EXTERNAL SSUM
       INTEGER ije,ijb,jje,jjb
-      
-c$OMP BARRIER
-
-      if (llm.eq.1) then
-        ! Specific behaviour for Shallow Water (1 vertical layer) case
-      
-        ! Sanity checks
-        if (kappa.ne.1) then
-          call abort_gcm("exner_hyb",
-     &    "kappa!=1 , but running in Shallow Water mode!!",42)
-        endif
-        if (cpp.ne.r) then
-        call abort_gcm("exner_hyb",
-     &    "cpp!=r , but running in Shallow Water mode!!",42)
+      logical,save :: firstcall=.true.
+!$OMP THREADPRIVATE(firstcall) 
+      character(len=*),parameter :: modname="exner_milieu_p"
+
+      ! Sanity check
+      if (firstcall) then
+        ! check that vertical discretization is compatible
+        ! with this routine
+        if (disvert_type.ne.2) then
+          call abort_gcm(modname,
+     &     "this routine should only be called if disvert_type==2",42)
         endif
         
+        ! sanity checks for Shallow Water case (1 vertical layer)
+        if (llm.eq.1) then
+          if (kappa.ne.1) then
+            call abort_gcm(modname,
+     &      "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+            call abort_gcm(modname,
+     &      "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+        endif ! of if (llm.eq.1)
+
+        firstcall=.false.
+      endif ! of if (firstcall)
+      
+c$OMP BARRIER
+
+! Specific behaviour for Shallow Water (1 vertical layer) case
+      if (llm.eq.1) then
+              
         ! Compute pks(:),pk(:),pkf(:)
         ijb=ij_begin
@@ -113 +130 @@
       endif ! of if (llm.eq.1)
 
+!!!! General case:
+
 c     -------------
 c     Calcul de pks

trunk/libf/dyn3dpar/guide_p_mod.F90

@@ -19 +19 @@
 ! ---------------------------------------------
   INTEGER, PRIVATE, SAVE  :: iguide_read,iguide_int,iguide_sav
-  INTEGER, PRIVATE, SAVE  :: nlevnc
+  INTEGER, PRIVATE, SAVE  :: nlevnc, guide_plevs
   LOGICAL, PRIVATE, SAVE  :: guide_u,guide_v,guide_T,guide_Q,guide_P
   LOGICAL, PRIVATE, SAVE  :: guide_hr,guide_teta  
   LOGICAL, PRIVATE, SAVE  :: guide_BL,guide_reg,guide_add,gamma4,guide_zon 
-  LOGICAL, PRIVATE, SAVE  :: guide_modele,invert_p,invert_y,ini_anal
-  LOGICAL, PRIVATE, SAVE  :: guide_2D,guide_sav
+  LOGICAL, PRIVATE, SAVE  :: invert_p,invert_y,ini_anal
+  LOGICAL, PRIVATE, SAVE  :: guide_2D,guide_sav,guide_modele
   
   REAL, PRIVATE, SAVE     :: tau_min_u,tau_max_u
@@ -48 +48 @@
   REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: tnat1,tnat2
   REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: qnat1,qnat2
+  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: pnat1,pnat2
   REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: psnat1,psnat2
   REAL, ALLOCATABLE, DIMENSION(:),     PRIVATE, SAVE   :: apnc,bpnc
@@ -65 +66 @@
 
   SUBROUTINE guide_init
-
+    
     USE control_mod
     IMPLICIT NONE
@@ -127 +128 @@
 ! Parametres pour lecture des fichiers
     CALL getpar('iguide_read',4,iguide_read,'freq. lecture guidage')
-    CALL getpar('iguide_int',4,iguide_int,'freq. lecture guidage')
-    IF (iguide_int.GT.0) THEN
+    CALL getpar('iguide_int',4,iguide_int,'freq. interpolation vert')
+    IF (iguide_int.EQ.0) THEN
+        iguide_int=1
+    ELSEIF (iguide_int.GT.0) THEN
         iguide_int=day_step/iguide_int
     ELSE
         iguide_int=day_step*iguide_int
     ENDIF
-    CALL getpar('guide_modele',.false.,guide_modele,'guidage niveaux modele')
+    CALL getpar('guide_plevs',0,guide_plevs,'niveaux pression fichiers guidage')
+    ! Pour compatibilite avec ancienne version avec guide_modele
+    CALL getpar('guide_modele',.false.,guide_modele,'niveaux pression ap+bp*psol')
+    IF (guide_modele) THEN
+        guide_plevs=1
+    ENDIF
+    ! Fin raccord
     CALL getpar('ini_anal',.false.,ini_anal,'Etat initial = analyse')
     CALL getpar('guide_invertp',.true.,invert_p,'niveaux p inverses')
@@ -144 +153 @@
 ! ---------------------------------------------
     ncidpl=-99
-    if (guide_modele) then
+    if (guide_plevs.EQ.1) then
        if (ncidpl.eq.-99) rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
-    else
-         if (guide_u) then
-           if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
-         elseif (guide_v) then
-           if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
-         elseif (guide_T) then
-           if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
-         elseif (guide_Q) then
-           if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
-         endif
+    elseif (guide_plevs.EQ.2) then
+       if (ncidpl.EQ.-99) rcod=nf90_open('P.nc',Nf90_NOWRITe,ncidpl)
+    elseif (guide_u) then
+       if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
+    elseif (guide_v) then
+       if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
+    elseif (guide_T) then
+       if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
+    elseif (guide_Q) then
+       if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
     endif 
     error=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
@@ -240 +249 @@
     ENDIF
 
-    IF (guide_P.OR.guide_modele) THEN
+    IF (guide_plevs.EQ.2) THEN
+        ALLOCATE(pnat1(iip1,jjp1,nlevnc), stat = error)
+        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
+        ALLOCATE(pnat2(iip1,jjp1,nlevnc), stat = error)
+        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
+        pnat1=0.;pnat2=0.;
+    ENDIF
+
+    IF (guide_P.OR.guide_plevs.EQ.1) THEN
         ALLOCATE(psnat1(iip1,jjp1), stat = error)
         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
@@ -267 +284 @@
     IF (guide_T) tnat1=tnat2
     IF (guide_Q) qnat1=qnat2
-    IF (guide_P.OR.guide_modele) psnat1=psnat2
+    IF (guide_plevs.EQ.2) pnat1=pnat2
+    IF (guide_P.OR.guide_plevs.EQ.1) psnat1=psnat2
 
   END SUBROUTINE guide_init
@@ -293 +311 @@
     LOGICAL       :: f_out ! sortie guidage
     REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
-    REAL, DIMENSION (ip1jmp1,llm) :: p ! besoin si guide_P
+    ! Variables pour fonction Exner (P milieu couche)
+    REAL, DIMENSION (iip1,jjp1,llm)    :: pk, pkf
+    REAL, DIMENSION (iip1,jjp1,llm)    :: alpha, beta
+    REAL, DIMENSION (iip1,jjp1)        :: pks    
+    REAL                               :: unskap
+    REAL, DIMENSION (ip1jmp1,llmp1)    :: p ! besoin si guide_P
     ! Compteurs temps:
     INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage
@@ -300 +323 @@
     REAL, SAVE    :: factt ! pas de temps en fraction de jour
     
-    INTEGER       :: l
+    INTEGER       :: i,j,l
     
     ijb_u=ij_begin ; ije_u=ij_end ; ijn_u=ije_u-ijb_u+1  
@@ -313 +336 @@
     ENDIF
       
-    
-    
      PRINT *,'---> on rentre dans guide_main'
 !    CALL AllGather_Field(ucov,ip1jmp1,llm)
@@ -380 +401 @@
       dday_step=real(day_step)
       IF (iguide_read.LT.0) THEN
-          tau=ditau/dday_step/ REAL(iguide_read)
+          tau=ditau/dday_step/REAL(iguide_read)
       ELSE
-          tau= REAL(iguide_read)*ditau/dday_step
+          tau=REAL(iguide_read)*ditau/dday_step
       ENDIF
       reste=tau-AINT(tau)
@@ -394 +415 @@
               IF (guide_T) tnat1(:,jjb_u:jje_u,:)=tnat2(:,jjb_u:jje_u,:)
               IF (guide_Q) qnat1(:,jjb_u:jje_u,:)=qnat2(:,jjb_u:jje_u,:)
-              IF (guide_P.OR.guide_modele) psnat1(:,jjb_u:jje_u)=psnat2(:,jjb_u:jje_u)
+              IF (guide_plevs.EQ.2) pnat1(:,jjb_u:jje_u,:)=pnat2(:,jjb_u:jje_u,:)
+              IF (guide_P.OR.guide_plevs.EQ.1) psnat1(:,jjb_u:jje_u)=psnat2(:,jjb_u:jje_u)
               step_rea=step_rea+1
               itau_test=itau
@@ -430 +452 @@
 ! Sauvegarde du guidage?
     f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav)  
-    IF (f_out) CALL guide_out("S",jjp1,1,ps)
+    IF (f_out) THEN
+!       Calcul niveaux pression milieu de couches 
+        CALL pression_p( ip1jmp1, ap, bp, ps, p )
+        if (disvert_type==1) then
+          CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+        else
+          CALL exner_milieu_p(ip1jmp1,ps,p,beta,pks,pk,pkf)
+        endif
+        unskap=1./kappa
+        DO l = 1, llm
+            DO j=jjb_u,jje_u
+                DO i =1, iip1
+                    p(i+(j-1)*iip1,l) = preff * ( pk(i,j,l)/cpp) ** unskap
+                ENDDO
+            ENDDO
+        ENDDO
+        CALL guide_out("P",jjp1,llm,p(1:ip1jmp1,1:llm),1.)
+    ENDIF
     
     if (guide_u) then
@@ -441 +480 @@
         if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
-        IF (f_out) CALL guide_out("U",jjp1,llm,f_add/factt)
+        IF (f_out) CALL guide_out("U",jjp1,llm,f_add(:,:),factt)
         ucov(ijb_u:ije_u,:)=ucov(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
     endif
@@ -453 +492 @@
         if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
-        IF (f_out) CALL guide_out("T",jjp1,llm,f_add/factt)
+        IF (f_out) CALL guide_out("T",jjp1,llm,f_add(:,:),factt)
         teta(ijb_u:ije_u,:)=teta(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
     endif
@@ -465 +504 @@
         if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
         CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
-        IF (f_out) CALL guide_out("P",jjp1,1,f_add(1:ip1jmp1,1)/factt)
+        IF (f_out) CALL guide_out("SP",jjp1,1,f_add(1:ip1jmp1,1),factt)
         ps(ijb_u:ije_u)=ps(ijb_u:ije_u)+f_add(ijb_u:ije_u,1)
         CALL pression_p(ip1jmp1,ap,bp,ps,p)
@@ -479 +518 @@
         if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
-        IF (f_out) CALL guide_out("Q",jjp1,llm,f_add/factt)
+        IF (f_out) CALL guide_out("Q",jjp1,llm,f_add(:,:),factt)
         q(ijb_u:ije_u,:)=q(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
     endif
@@ -492 +531 @@
         if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
         CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
-        IF (f_out) CALL guide_out("V",jjm,llm,f_add(1:ip1jm,:)/factt)
+        IF (f_out) CALL guide_out("V",jjm,llm,f_add(1:ip1jm,:),factt)
         vcov(ijb_v:ije_v,:)=vcov(ijb_v:ije_v,:)+f_add(ijb_v:ije_v,:)
     endif
@@ -580 +619 @@
               ENDDO
           ENDDO 
-          fieldm(:,l)=fieldm(:,l)/ REAL(imax(typ)-imin(typ)+1)
+          fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
     ! Compute forcing
           DO j=jjb_v,jje_v
@@ -598 +637 @@
               ENDDO
           ENDDO 
-          fieldm(:,l)=fieldm(:,l)/ REAL(imax(typ)-imin(typ)+1)
+          fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
     ! Compute forcing
           DO j=jjb_u,jje_u
@@ -640 +679 @@
   REAL, DIMENSION (iip1,jjp1,llm)    :: alpha, beta
   REAL, DIMENSION (iip1,jjp1)        :: pks    
-  REAL                               :: prefkap,unskap
+  REAL                               :: unskap
   ! Pression de vapeur saturante
   REAL, DIMENSION (ip1jmp1,llm)      :: qsat
@@ -652 +691 @@
     print *,'Guide: conversion variables guidage'
 ! -----------------------------------------------------------------
-! Calcul des niveaux de pression champs guidage
+! Calcul des niveaux de pression champs guidage (pour T et Q)
 ! -----------------------------------------------------------------
-if (guide_modele) then
-    do i=1,iip1
-        do j=jjb_u,jje_u
-            do l=1,nlevnc
-                plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
-                plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
-            enddo
-        enddo
-    enddo
-else
-    do i=1,iip1
-        do j=jjb_u,jje_u
-            do l=1,nlevnc
-                plnc2(i,j,l)=apnc(l)
-                plnc1(i,j,l)=apnc(l)
-           enddo
-        enddo
-    enddo
-
-endif
+    IF (guide_plevs.EQ.0) THEN
+        DO l=1,nlevnc
+            DO j=jjb_u,jje_u
+                DO i=1,iip1
+                    plnc2(i,j,l)=apnc(l)
+                    plnc1(i,j,l)=apnc(l)
+               ENDDO
+            ENDDO
+        ENDDO
+    ENDIF   
+
     if (first) then
         first=.FALSE.
@@ -683 +713 @@
         enddo
         print*,'Fichiers guidage'
-        do l=1,nlevnc
-             print*,'PL(',l,')=',plnc2(1,jjb_u,l)
-        enddo
+        SELECT CASE (guide_plevs)
+        CASE (0) 
+            do l=1,nlevnc
+                 print*,'PL(',l,')=',plnc2(1,jjb_u,l)
+            enddo
+        CASE (1)
+            DO l=1,nlevnc
+                 print*,'PL(',l,')=',apnc(l)+bpnc(l)*psnat2(i,jjb_u)
+             ENDDO
+        CASE (2)
+            do l=1,nlevnc
+                 print*,'PL(',l,')=',pnat2(1,jjb_u,l)
+            enddo
+        END SELECT
         print *,'inversion de l''ordre: invert_p=',invert_p
         if (guide_u) then
@@ -702 +743 @@
 ! Calcul niveaux pression modele 
 ! -----------------------------------------------------------------
-    CALL pression_p( ip1jmp1, ap, bp, psi, p )
-    CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
 
 !    ....  Calcul de pls , pression au milieu des couches ,en Pascals
-    unskap=1./kappa
-    prefkap =  preff  ** kappa
-    DO l = 1, llm
-        DO j=jjb_u,jje_u
-            DO i =1, iip1
-                pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
-            ENDDO
+    IF (guide_plevs.EQ.1) THEN
+        DO l=1,llm
+            DO j=jjb_u,jje_u
+                DO i =1, iip1
+                    pls(i,j,l)=(ap(l)+ap(l+1))/2.+psi(i,j)*(bp(l)+bp(l+1))/2.
+                ENDDO
+            ENDDO
         ENDDO
-    ENDDO
+    ELSE
+        CALL pression_p( ip1jmp1, ap, bp, psi, p )
+        if (disvert_type==1) then
+          CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
+        else ! we assume that we are in the disvert_type==2 case
+          CALL exner_milieu_p(ip1jmp1,psi,p,beta,pks,pk,pkf)
+        endif
+        unskap=1./kappa
+        DO l = 1, llm
+            DO j=jjb_u,jje_u
+                DO i =1, iip1
+                    pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
+                ENDDO
+            ENDDO
+        ENDDO
+    ENDIF
 
 !   calcul des pressions pour les grilles u et v
@@ -747 +801 @@
 
 ! -----------------------------------------------------------------
-! Interpolation champs guidage sur niveaux modele (+inversion N/S)
+! Interpolation verticale champs guidage sur niveaux modele
 ! Conversion en variables gcm (ucov, vcov...)
 ! -----------------------------------------------------------------
@@ -762 +816 @@
     endif
 
+    IF (guide_T) THEN
+        ! Calcul des nouvelles valeurs des niveaux de pression du guidage
+        IF (guide_plevs.EQ.1) THEN
+            DO l=1,nlevnc
+                DO j=jjb_u,jje_u
+                    DO i=1,iip1
+                        plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
+                        plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
+                    ENDDO
+                ENDDO
+            ENDDO
+        ELSE IF (guide_plevs.EQ.2) THEN
+            DO l=1,nlevnc
+                DO j=jjb_u,jje_u
+                    DO i=1,iip1
+                        plnc2(i,j,l)=pnat2(i,j,l)
+                        plnc1(i,j,l)=pnat1(i,j,l)
+                    ENDDO
+                ENDDO
+            ENDDO
+        ENDIF
+
+        ! Interpolation verticale
+        CALL pres2lev(tnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm,           &
+                    plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
+        CALL pres2lev(tnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm,           &
+                    plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
+
+        ! Conversion en variables GCM
+        do l=1,llm
+            do j=jjb_u,jje_u
+                IF (guide_teta) THEN
+                    do i=1,iim
+                        ij=(j-1)*iip1+i
+                        tgui1(ij,l)=zu1(i,j,l)
+                        tgui2(ij,l)=zu2(i,j,l)
+                    enddo
+                ELSE
+                    do i=1,iim
+                        ij=(j-1)*iip1+i
+                        tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
+                        tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
+                    enddo
+                ENDIF
+                tgui1(j*iip1,l)=tgui1((j-1)*iip1+1,l)    
+                tgui2(j*iip1,l)=tgui2((j-1)*iip1+1,l)    
+            enddo
+            do i=1,iip1
+                tgui1(i,l)=tgui1(1,l)
+                tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l) 
+                tgui2(i,l)=tgui2(1,l)
+                tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l) 
+            enddo
+        enddo
+    ENDIF
+
+    IF (guide_Q) THEN
+        ! Calcul des nouvelles valeurs des niveaux de pression du guidage
+        IF (guide_plevs.EQ.1) THEN
+            DO l=1,nlevnc
+                DO j=jjb_u,jje_u
+                    DO i=1,iip1
+                        plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
+                        plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
+                    ENDDO
+                ENDDO
+            ENDDO
+        ELSE IF (guide_plevs.EQ.2) THEN
+            DO l=1,nlevnc
+                DO j=jjb_u,jje_u
+                    DO i=1,iip1
+                        plnc2(i,j,l)=pnat2(i,j,l)
+                        plnc1(i,j,l)=pnat1(i,j,l)
+                    ENDDO
+                ENDDO
+            ENDDO
+        ENDIF
+
+        ! Interpolation verticale
+        CALL pres2lev(qnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm,             &
+                      plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
+        CALL pres2lev(qnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm,             &
+                      plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
+
+        ! Conversion en variables GCM
+        ! On suppose qu'on a la bonne variable dans le fichier de guidage:
+        ! Hum.Rel si guide_hr, Hum.Spec. sinon.
+        do l=1,llm
+            do j=jjb_u,jje_u
+                do i=1,iim
+                    ij=(j-1)*iip1+i
+                    qgui1(ij,l)=zu1(i,j,l)
+                    qgui2(ij,l)=zu2(i,j,l)
+                enddo
+                qgui1(j*iip1,l)=qgui1((j-1)*iip1+1,l)    
+                qgui2(j*iip1,l)=qgui2((j-1)*iip1+1,l)    
+            enddo
+            do i=1,iip1
+                qgui1(i,l)=qgui1(1,l)
+                qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l) 
+                qgui2(i,l)=qgui2(1,l)
+                qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l) 
+            enddo
+        enddo
+        IF (guide_hr) THEN
+            CALL q_sat(iip1*jjn_u*llm,teta(:,jjb_u:jje_u,:)*pk(:,jjb_u:jje_u,:)/cpp,       &
+                       plsnc(:,jjb_u:jje_u,:),qsat(ijb_u:ije_u,:))
+            qgui1(ijb_u:ije_u,:)=qgui1(ijb_u:ije_u,:)*qsat(ijb_u:ije_u,:)*0.01 !hum. rel. en %
+            qgui2(ijb_u:ije_u,:)=qgui2(ijb_u:ije_u,:)*qsat(ijb_u:ije_u,:)*0.01 
+        ENDIF
+    ENDIF
+
     IF (guide_u) THEN
+        ! Calcul des nouvelles valeurs des niveaux de pression du guidage
+        IF (guide_plevs.EQ.1) THEN
+            DO l=1,nlevnc
+                DO j=jjb_u,jje_u
+                    DO i=1,iim
+                        plnc2(i,j,l)=apnc(l)+bpnc(l)*(psnat2(i,j)*aire(i,j)*alpha1p2(i,j) &
+                       &           +psnat2(i+1,j)*aire(i+1,j)*alpha3p4(i+1,j))/aireu(i,j)
+                        plnc1(i,j,l)=apnc(l)+bpnc(l)*(psnat1(i,j)*aire(i,j)*alpha1p2(i,j) &
+                       &           +psnat1(i+1,j)*aire(i+1,j)*alpha3p4(i+1,j))/aireu(i,j)
+                    ENDDO
+                    plnc2(iip1,j,l)=plnc2(1,j,l)
+                    plnc1(iip1,j,l)=plnc1(1,j,l)
+                ENDDO
+            ENDDO
+        ELSE IF (guide_plevs.EQ.2) THEN
+            DO l=1,nlevnc
+                DO j=jjb_u,jje_u
+                    DO i=1,iim
+                        plnc2(i,j,l)=(pnat2(i,j,l)*aire(i,j)*alpha1p2(i,j) &
+                       & +pnat2(i+1,j,l)*aire(i,j)*alpha3p4(i+1,j))/aireu(i,j)
+                        plnc1(i,j,l)=(pnat1(i,j,l)*aire(i,j)*alpha1p2(i,j) &
+                       & +pnat1(i+1,j,l)*aire(i,j)*alpha3p4(i+1,j))/aireu(i,j)
+                    ENDDO
+                    plnc2(iip1,j,l)=plnc2(1,j,l)
+                    plnc1(iip1,j,l)=plnc1(1,j,l)
+                ENDDO
+            ENDDO
+        ENDIF
+        
+        ! Interpolation verticale
         CALL pres2lev(unat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm,            &
                       plnc1(:,jjb_u:jje_u,:),plunc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
@@ -768 +964 @@
                       plnc2(:,jjb_u:jje_u,:),plunc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
 
+        ! Conversion en variables GCM
         do l=1,llm
             do j=jjb_u,jje_u
@@ -787 +984 @@
     ENDIF
     
-    IF (guide_T) THEN
-        CALL pres2lev(tnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm,           &
-                      plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
-        CALL pres2lev(tnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm,           &
-                      plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
-
-        do l=1,llm
-            do j=jjb_u,jje_u
-                IF (guide_teta) THEN
-                    do i=1,iim
-                        ij=(j-1)*iip1+i
-                        tgui1(ij,l)=zu1(i,j,l)
-                        tgui2(ij,l)=zu2(i,j,l)
-                    enddo
-                ELSE
-                    do i=1,iim
-                        ij=(j-1)*iip1+i
-                        tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
-                        tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
-                    enddo
-                ENDIF
-                tgui1(j*iip1,l)=tgui1((j-1)*iip1+1,l)    
-                tgui2(j*iip1,l)=tgui2((j-1)*iip1+1,l)    
-            enddo
-            do i=1,iip1
-                tgui1(i,l)=tgui1(1,l)
-                tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l) 
-                tgui2(i,l)=tgui2(1,l)
-                tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l) 
-            enddo
-        enddo
-    ENDIF
-
     IF (guide_v) THEN
-        
+        ! Calcul des nouvelles valeurs des niveaux de pression du guidage
+        IF (guide_plevs.EQ.1) THEN
+         CALL Register_SwapFieldHallo(psnat1,psnat1,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
+         CALL SendRequest(Req)
+         CALL WaitRequest(Req)
+         CALL Register_SwapFieldHallo(psnat2,psnat2,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
+         CALL SendRequest(Req)
+         CALL WaitRequest(Req)
+            DO l=1,nlevnc
+                DO j=jjb_v,jje_v
+                    DO i=1,iip1
+                        plnc2(i,j,l)=apnc(l)+bpnc(l)*(psnat2(i,j)*aire(i,j)*alpha2p3(i,j) &
+                       &           +psnat2(i,j+1)*aire(i,j+1)*alpha1p4(i,j+1))/airev(i,j)
+                        plnc1(i,j,l)=apnc(l)+bpnc(l)*(psnat1(i,j)*aire(i,j)*alpha2p3(i,j) &
+                       &           +psnat1(i,j+1)*aire(i,j+1)*alpha1p4(i,j+1))/airev(i,j)
+                    ENDDO
+                ENDDO
+            ENDDO
+        ELSE IF (guide_plevs.EQ.2) THEN
+         CALL Register_SwapFieldHallo(pnat1,pnat1,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
+         CALL SendRequest(Req)
+         CALL WaitRequest(Req)
+         CALL Register_SwapFieldHallo(pnat2,pnat2,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
+         CALL SendRequest(Req)
+         CALL WaitRequest(Req)
+            DO l=1,nlevnc
+                DO j=jjb_v,jje_v
+                    DO i=1,iip1
+                        plnc2(i,j,l)=(pnat2(i,j,l)*aire(i,j)*alpha2p3(i,j) &
+                       & +pnat2(i,j+1,l)*aire(i,j)*alpha1p4(i,j+1))/airev(i,j)
+                        plnc1(i,j,l)=(pnat1(i,j,l)*aire(i,j)*alpha2p3(i,j) &
+                       & +pnat1(i,j+1,l)*aire(i,j)*alpha1p4(i,j+1))/airev(i,j)
+                    ENDDO
+                ENDDO
+            ENDDO
+        ENDIF
+        ! Interpolation verticale
         CALL pres2lev(vnat1(:,jjb_v:jje_v,:),zv1(:,jjb_v:jje_v,:),nlevnc,llm,             &
                       plnc1(:,jjb_v:jje_v,:),plvnc(:,jjb_v:jje_v,:),iip1,jjn_v,invert_p)
         CALL pres2lev(vnat2(:,jjb_v:jje_v,:),zv2(:,jjb_v:jje_v,:),nlevnc,llm,             &
                       plnc2(:,jjb_v:jje_v,:),plvnc(:,jjb_v:jje_v,:),iip1,jjn_v,invert_p)
-
+        ! Conversion en variables GCM
         do l=1,llm
             do j=jjb_v,jje_v
@@ -840 +1040 @@
     ENDIF
     
-    IF (guide_Q) THEN
-        ! On suppose qu'on a la bonne variable dans le fichier de guidage:
-        ! Hum.Rel si guide_hr, Hum.Spec. sinon.
-        CALL pres2lev(qnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm,             &
-                      plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
-        CALL pres2lev(qnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm,             &
-                      plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
-
-        do l=1,llm
-            do j=jjb_u,jjb_v
-                do i=1,iim
-                    ij=(j-1)*iip1+i
-                    qgui1(ij,l)=zu1(i,j,l)
-                    qgui2(ij,l)=zu2(i,j,l)
-                enddo
-                qgui1(j*iip1,l)=qgui1((j-1)*iip1+1,l)    
-                qgui2(j*iip1,l)=qgui2((j-1)*iip1+1,l)    
-            enddo
-            do i=1,iip1
-                qgui1(i,l)=qgui1(1,l)
-                qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l) 
-                qgui2(i,l)=qgui2(1,l)
-                qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l) 
-            enddo
-        enddo
-        IF (guide_hr) THEN
-            CALL q_sat(iip1*jjn_u*llm,teta(:,jjb_u:jje_u,:)*pk(:,jjb_u:jje_u,:)/cpp,       &
-                       plsnc(:,jjb_u:jje_u,:),qsat(ijb_u:ije_u,:))
-            qgui1(ijb_u:ije_u,:)=qgui1(ijb_u:ije_u,:)*qsat(ijb_u:ije_u,:)*0.01 !hum. rel. en %
-            qgui2(ijb_u:ije_u,:)=qgui2(ijb_u:ije_u,:)*qsat(ijb_u:ije_u,:)*0.01 
-        ENDIF
-    ENDIF
 
   END SUBROUTINE guide_interp
@@ -1055 +1223 @@
     LOGICAL, SAVE         :: first=.TRUE.
 ! Identification fichiers et variables NetCDF:
-    INTEGER, SAVE         :: ncidu,varidu,ncidv,varidv,ncidQ
-    INTEGER, SAVE         :: varidQ,ncidt,varidt,ncidps,varidps
+    INTEGER, SAVE         :: ncidu,varidu,ncidv,varidv,ncidp,varidp
+    INTEGER, SAVE         :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps
     INTEGER               :: ncidpl,varidpl,varidap,varidbp
 ! Variables auxiliaires NetCDF:
@@ -1068 +1236 @@
          ncidpl=-99
          print*,'Guide: ouverture des fichiers guidage '
-! Niveaux de pression si non constants
-         if (guide_modele) then
-             print *,'Lecture du guidage sur niveaux mod�le'
+! Ap et Bp si Niveaux de pression hybrides
+         if (guide_plevs.EQ.1) then
+             print *,'Lecture du guidage sur niveaux modele'
              rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
              rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
              rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
              print*,'ncidpl,varidap',ncidpl,varidap
+         endif
+! Pression si guidage sur niveaux P variables
+         if (guide_plevs.EQ.2) then
+             rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
+             rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
+             print*,'ncidp,varidp',ncidp,varidp
+             if (ncidpl.eq.-99) ncidpl=ncidp
          endif
 ! Vent zonal
@@ -1105 +1280 @@
          endif
 ! Pression de surface
-         if ((guide_P).OR.(guide_modele)) then
+         if ((guide_P).OR.(guide_plevs.EQ.1)) then
              rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
              rcode = nf90_inq_varid(ncidps, 'SP', varidps)
@@ -1111 +1286 @@
          endif
 ! Coordonnee verticale
-         if (.not.guide_modele) then
+         if (guide_plevs.EQ.0) then
               rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
               IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
@@ -1117 +1292 @@
          endif
 ! Coefs ap, bp pour calcul de la pression aux differents niveaux
-         if (guide_modele) then
+         IF (guide_plevs.EQ.1) THEN
 #ifdef NC_DOUBLE
              status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
@@ -1125 +1300 @@
              status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
 #endif
-         else
+         ELSEIF (guide_plevs.EQ.0) THEN
 #ifdef NC_DOUBLE
              status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
@@ -1133 +1308 @@
              apnc=apnc*100.! conversion en Pascals
              bpnc(:)=0.
-         endif
+         ENDIF
          first=.FALSE.
-     endif ! (first)
+     ENDIF ! (first)
 
 ! -----------------------------------------------------------------
@@ -1152 +1327 @@
      count(4)=1
 
+! Pression 
+     if (guide_plevs.EQ.2) then
+#ifdef NC_DOUBLE
+         status=NF_GET_VARA_DOUBLE(ncidp,varidp,start,count,pnat2)
+#else
+         status=NF_GET_VARA_REAL(ncidp,varidp,start,count,pnat2)
+#endif
+         IF (invert_y) THEN
+           CALL invert_lat(iip1,jjp1,nlevnc,pnat2)
+         ENDIF
+     endif
+
 !  Vent zonal
      if (guide_u) then
@@ -1204 +1391 @@
 
 !  Pression de surface
-     if ((guide_P).OR.(guide_modele))  then
+     if ((guide_P).OR.(guide_plevs.EQ.1))  then
          start(3)=timestep
          start(4)=0
@@ -1235 +1422 @@
     LOGICAL, SAVE         :: first=.TRUE.
 ! Identification fichiers et variables NetCDF:
-    INTEGER, SAVE         :: ncidu,varidu,ncidv,varidv,ncidQ
-    INTEGER, SAVE         :: varidQ,ncidt,varidt,ncidps,varidps
+    INTEGER, SAVE         :: ncidu,varidu,ncidv,varidv,ncidp,varidp
+    INTEGER, SAVE         :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps
     INTEGER               :: ncidpl,varidpl,varidap,varidbp
 ! Variables auxiliaires NetCDF:
@@ -1252 +1439 @@
          ncidpl=-99
          print*,'Guide: ouverture des fichiers guidage '
-! Niveaux de pression si non constants
-         if (guide_modele) then
+! Ap et Bp si niveaux de pression hybrides
+         if (guide_plevs.EQ.1) then
              print *,'Lecture du guidage sur niveaux mod�le'
              rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
@@ -1259 +1446 @@
              rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
              print*,'ncidpl,varidap',ncidpl,varidap
+         endif
+! Pression
+         if (guide_plevs.EQ.2) then
+             rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
+             rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
+             print*,'ncidp,varidp',ncidp,varidp
+             if (ncidpl.eq.-99) ncidpl=ncidp
          endif
 ! Vent zonal
@@ -1289 +1483 @@
          endif
 ! Pression de surface
-         if ((guide_P).OR.(guide_modele)) then
+         if ((guide_P).OR.(guide_plevs.EQ.1)) then
              rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
              rcode = nf90_inq_varid(ncidps, 'SP', varidps)
@@ -1295 +1489 @@
          endif
 ! Coordonnee verticale
-         if (.not.guide_modele) then
+         if (guide_plevs.EQ.0) then
               rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
               IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
@@ -1301 +1495 @@
          endif
 ! Coefs ap, bp pour calcul de la pression aux differents niveaux
-         if (guide_modele) then
+         if (guide_plevs.EQ.1) then
 #ifdef NC_DOUBLE
              status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
@@ -1309 +1503 @@
              status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
 #endif
-         else
+         elseif (guide_plevs.EQ.0) THEN
 #ifdef NC_DOUBLE
              status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
@@ -1336 +1530 @@
      count(4)=1
 
+!  Pression
+     if (guide_plevs.EQ.2) then
+#ifdef NC_DOUBLE
+         status=NF_GET_VARA_DOUBLE(ncidp,varidp,start,count,zu)
+#else
+         status=NF_GET_VARA_REAL(ncidp,varidp,start,count,zu)
+#endif
+         DO i=1,iip1
+             pnat2(i,:,:)=zu(:,:)
+         ENDDO
+
+         IF (invert_y) THEN
+           CALL invert_lat(iip1,jjp1,nlevnc,pnat2)
+         ENDIF
+     endif
 !  Vent zonal
      if (guide_u) then
@@ -1350 +1559 @@
            CALL invert_lat(iip1,jjp1,nlevnc,unat2)
          ENDIF
-
      endif
 
@@ -1367 +1575 @@
            CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
          ENDIF
-
      endif
 
@@ -1384 +1591 @@
            CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
          ENDIF
-
      endif
 
@@ -1402 +1608 @@
            CALL invert_lat(iip1,jjm,nlevnc,vnat2)
          ENDIF
-
      endif
 
 !  Pression de surface
-     if ((guide_P).OR.(guide_modele))  then
+     if ((guide_P).OR.(guide_plevs.EQ.1))  then
          start(3)=timestep
          start(4)=0
@@ -1424 +1629 @@
            CALL invert_lat(iip1,jjp1,1,psnat2)
          ENDIF
-
      endif
 
@@ -1430 +1634 @@
   
 !=======================================================================
-  SUBROUTINE guide_out(varname,hsize,vsize,field)
+  SUBROUTINE guide_out(varname,hsize,vsize,field,factt)
     USE parallel
     IMPLICIT NONE
@@ -1445 +1649 @@
     INTEGER,   INTENT (IN)                         :: hsize,vsize
     REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
+    REAL, INTENT (IN)                              :: factt
 
     ! Variables locales
@@ -1507 +1712 @@
 ! --------------------------------------------------------------------
         ierr = NF_REDEF(nid)
-! Surface pressure (GCM)
-        dim3=(/id_lonv,id_latu,id_tim/)
-        ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,3,dim3,varid)
+! Pressure (GCM)
+        dim4=(/id_lonv,id_latu,id_lev,id_tim/)
+        ierr = NF_DEF_VAR(nid,"P",NF_FLOAT,4,dim4,varid)
 ! Surface pressure (guidage)
         IF (guide_P) THEN
@@ -1543 +1748 @@
 ! Enregistrement du champ
 ! --------------------------------------------------------------------
+ 
     ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
 
     SELECT CASE (varname)
-    CASE ("S")
+    CASE ("P")
         timestep=timestep+1
-        ierr = NF_INQ_VARID(nid,"SP",varid)
-        start=(/1,1,timestep,0/)
-        count=(/iip1,jjp1,1,0/)
+        ierr = NF_INQ_VARID(nid,"P",varid)
+        start=(/1,1,1,timestep/)
+        count=(/iip1,jjp1,llm,1/)
 #ifdef NC_DOUBLE
         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
@@ -1556 +1762 @@
         ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
 #endif
-    CASE ("P")
+    CASE ("SP")
         ierr = NF_INQ_VARID(nid,"ps",varid)
         start=(/1,1,timestep,0/)
         count=(/iip1,jjp1,1,0/)
 #ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
+        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
+#else
+        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
 #endif
     CASE ("U")
@@ -1570 +1776 @@
         count=(/iip1,jjp1,llm,1/)
 #ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
+        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
+#else
+        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
 #endif
     CASE ("V")
@@ -1579 +1785 @@
         count=(/iip1,jjm,llm,1/)
 #ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
+        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
+#else
+        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
 #endif
     CASE ("T")
@@ -1588 +1794 @@
         count=(/iip1,jjp1,llm,1/)
 #ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
+        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
+#else
+        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
 #endif
     CASE ("Q")
@@ -1597 +1803 @@
         count=(/iip1,jjp1,llm,1/)
 #ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
+        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
+#else
+        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
 #endif
     END SELECT

trunk/libf/dyn3dpar/iniacademic.F90

@@ -1 +1 @@
 !
-! $Id: iniacademic.F90 1474 2011-01-14 11:04:45Z lguez $
+! $Id: iniacademic.F90 1520 2011-05-23 11:37:09Z emillour $
 !
 SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
@@ -71 +71 @@
 
   REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
+  
+  character(len=*),parameter :: modname="iniacademic"
+  character(len=80) :: abort_message
 
   !-----------------------------------------------------------------------
@@ -76 +79 @@
   ! --------------------------------------
   !
-  
-  write(lunout,*) 'Iniacademic'
-  
   ! initialize planet radius, rotation rate,...
   call conf_planete
@@ -136 +136 @@
      teta0=315.     ! mean Teta (S.H. 315K)
      CALL getin('teta0',teta0)
-     write(lunout,*) 'Iniacademic - teta0 ',teta0
-     write(lunout,*) 'Iniacademic - rad ',rad
      ttp=200.       ! Tropopause temperature (S.H. 200K)
      CALL getin('ttp',ttp)
@@ -183 +181 @@
            tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
                 -delt_z*(1.-ddsin*ddsin)*log(zsig)
-           !! Aymeric -- tests particuliers
            if (planet_type=="giant") then
              tetajl(j,l)=teta0+(delt_y*                   &
@@ -207 +204 @@
         enddo
      enddo
-!     write(lunout,*) 'Iniacademic - check',tetajl(:,int(llm/2)),rlatu(:)
 
      ! 3. Initialize fields (if necessary)
@@ -217 +213 @@
 
         CALL pression ( ip1jmp1, ap, bp, ps, p       )
-        if (planet_type=="earth") then
-          CALL exner_hyb(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+        if (disvert_type.eq.1) then
+          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        elseif (disvert_type.eq.2) then
+          call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
         else
-          call exner_milieu(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+          write(abort_message,*) "Wrong value for disvert_type: ", &
+                              disvert_type
+          call abort_gcm(modname,abort_message,0)
         endif
         CALL massdair(p,masse)

trunk/libf/dyn3dpar/iniconst.F

@@ -1 +1 @@
 !
-! $Id: iniconst.F 1403 2010-07-01 09:02:53Z fairhead $
+! $Id: iniconst.F 1520 2011-05-23 11:37:09Z emillour $
 !
       SUBROUTINE iniconst
 
       USE control_mod
+#ifdef CPP_IOIPSL
+      use IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      use ioipsl_getincom
+#endif
 
       IMPLICIT NONE
@@ -22 +28 @@
 
 
+      character(len=*),parameter :: modname="iniconst"
+      character(len=80) :: abort_message
 c
 c
@@ -49 +57 @@
       r       = cpp * kappa
 
-      write(lunout,*)'iniconst: R  CP  Kappa ',  r , cpp,  kappa
+      write(lunout,*) trim(modname),': R  CP  Kappa ',r,cpp,kappa
 c
 c-----------------------------------------------------------------------
 
-      if (planet_type.eq."earth") then
-       CALL disvert_terre(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
+! vertical discretization: default behavior depends on planet_type flag
+      if (planet_type=="earth") then
+        disvert_type=1
       else
-       CALL disvert_noterre
+        disvert_type=2
       endif
-c
-       RETURN
-       END
+      ! but user can also specify using one or the other in run.def:
+      call getin('disvert_type',disvert_type)
+      write(lunout,*) trim(modname),': disvert_type=',disvert_type
+      
+      if (disvert_type==1) then
+       ! standard case for Earth (automatic generation of levels)
+       call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,
+     &              scaleheight)
+      else if (disvert_type==2) then
+        ! standard case for planets (levels generated using z2sig.def file)
+        call disvert_noterre
+      else
+        write(abort_message,*) "Wrong value for disvert_type: ",
+     &                        disvert_type
+        call abort_gcm(modname,abort_message,0)
+      endif
+
+      END

trunk/libf/dyn3dpar/leapfrog_p.F

@@ -271 +271 @@
 
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
-      if (planet_type.eq."earth") then
+      if (disvert_type==1) then
         CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-      else
+      else ! we assume that we are in the disvert_type==2 case
         CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
       endif
@@ -751 +751 @@
 
 c$OMP BARRIER
-         if (planet_type.eq."earth") then
-          CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-         else
-          CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
+         if (disvert_type==1) then
+           CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
+         else ! we assume that we are in the disvert_type==2 case
+           CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
          endif
 c$OMP BARRIER
@@ -1105 +1105 @@
         CALL pression_p ( ip1jmp1, ap, bp, ps, p                  )
 c$OMP BARRIER
-        if (planet_type.eq."earth") then
+        if (disvert_type==1) then
           CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-        else
+        else ! we assume that we are in the disvert_type==2 case
           CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
         endif

trunk/libf/dyn3dpar/logic.h

@@ -1 +1 @@
 !
-! $Id: logic.h 1319 2010-02-23 21:29:54Z fairhead $
+! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
 !
 !
@@ -11 +11 @@
      &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
      &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
-     &  ,ok_limit,ok_etat0,hybrid
+     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf,hybrid
 
       COMMON/logici/ iflag_phys,iflag_trac
@@ -18 +18 @@
      & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
      &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
-     &  ,ok_limit,ok_etat0
-      logical hybrid ! vertcal coordinate is hybrid if true (sigma otherwise)
+     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf
+      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
+                     ! (only used if disvert_type==2)
 
-      INTEGER iflag_phys,iflag_trac
+      integer iflag_phys,iflag_trac
 !$OMP THREADPRIVATE(/logicl/)
 !$OMP THREADPRIVATE(/logici/)