Changeset 1302 for trunk/LMDZ.COMMON/libf

Timestamp:

Jun 26, 2014, 6:07:05 PM (11 years ago)

Author:

emillour

Message:

Common dynamics:
Some updates to keep up with LMDZ5 Earth model evolution
(up to LMDZ5 rev 2070). See file "DOC/chantiers/commit_importants.log"
for detailed list of changes.
Note that the updates of exner* routines change (as expected) results
at numerical roundoff level.
EM

Location:

trunk/LMDZ.COMMON/libf

Files:

• bibio/pres2lev_mod.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/pres2lev_mod.F90)
• bibio/wxios.F90 (modified) (13 diffs)
• dyn3d/calfis.F (modified) (3 diffs)
• dyn3d/ce0l.F90 (modified) (4 diffs)
• dyn3d/conf_gcm.F (modified) (29 diffs)
• dyn3d/gcm.F (modified) (3 diffs)
• dyn3d/guide_mod.F90 (modified) (17 diffs)
• dyn3d/leapfrog.F (modified) (6 diffs)
• dyn3d/logic.h (modified) (2 diffs)
• dyn3d_common/PVtheta.F (deleted)
• dyn3d_common/disvert.F90 (modified) (8 diffs)
• dyn3d_common/exner_hyb_m.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/exner_hyb.F) (1 diff)
• dyn3d_common/exner_milieu_m.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/exner_milieu.F) (1 diff)
• dyn3d_common/grilles_gcm_netcdf.F (deleted)
• dyn3d_common/iniacademic.F90 (modified) (3 diffs)
• dyn3d_common/iniconst.F90 (modified) (1 diff)
• dyn3d_common/q_sat.F (modified) (3 diffs)
• dyn3d_common/tourabs.F (deleted)
• dyn3dpar/calfis_p.F (modified) (4 diffs)
• dyn3dpar/ce0l.F90 (modified) (3 diffs)
• dyn3dpar/conf_gcm.F (modified) (16 diffs)
• dyn3dpar/exner_hyb_p_m.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p.F) (1 diff)
• dyn3dpar/exner_milieu_p_m.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p.F) (1 diff)
• dyn3dpar/gcm.F (modified) (1 diff)
• dyn3dpar/guide_p_mod.F90 (modified) (25 diffs)
• dyn3dpar/leapfrog_p.F (modified) (8 diffs)
• dyn3dpar/logic.h (modified) (2 diffs)
• dyn3dpar/mod_const_mpi.F90 (modified) (3 diffs)
• dyn3dpar/parallel_lmdz.F90 (modified) (3 diffs)

trunk/LMDZ.COMMON/libf/bibio/wxios.F90

@@ -26 +26 @@
     
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    !   str + i   =>   str_i   !!!!!!!!!!!!!!!!!!!!
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    
-    SUBROUTINE concat(str, str2, str_str2)
-        CHARACTER(len=*), INTENT(IN) :: str, str2
-        CHARACTER(len=20), INTENT(OUT) :: str_str2
-        
-        
-        str_str2 = TRIM(ADJUSTL(str//"_"//TRIM(ADJUSTL(str2))))
-        !IF (prt_level >= 10) WRITE(lunout,*) "Xios. ",str,"+",str2,"=",str_str2
-    END SUBROUTINE concat
-    
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !   36day => 36d etc     !!!!!!!!!!!!!!!!!!!!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -109 +96 @@
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-    SUBROUTINE wxios_init(xios_ctx_name, locom, outcom)
+    SUBROUTINE wxios_init(xios_ctx_name, locom, outcom, type_ocean)
         IMPLICIT NONE
         INCLUDE 'iniprint.h'
@@ -116 +103 @@
       INTEGER, INTENT(IN), OPTIONAL :: locom
       INTEGER, INTENT(OUT), OPTIONAL :: outcom
+      CHARACTER(len=6), INTENT(IN), OPTIONAL :: type_ocean
 
     
@@ -142 +130 @@
         g_ctx_name = xios_ctx_name
         
-        CALL wxios_context_init()
-        
+        ! Si couple alors init fait dans cpl_init
+        IF (.not. PRESENT(type_ocean)) THEN
+            CALL wxios_context_init()
+        ENDIF
+
     END SUBROUTINE wxios_init
 
@@ -158 +149 @@
         g_ctx = xios_ctx
 
-        IF (prt_level >= 10) WRITE(lunout,*) "wxios_context_init: Current context is ",trim(g_ctx_name)
-
+        IF (prt_level >= 10) THEN
+          WRITE(lunout,*) "wxios_context_init: Current context is ",trim(g_ctx_name)
+          WRITE(lunout,*) "     now call xios_solve_inheritance()"
+        ENDIF
         !Une première analyse des héritages:
         CALL xios_solve_inheritance()
@@ -303 +296 @@
     ! Pour déclarer un axe vertical !!!!!!!!!!!!!!!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    SUBROUTINE wxios_add_vaxis(axisgroup_id, axis_file, axis_size, axis_value)
-        IMPLICIT NONE
-        INCLUDE 'iniprint.h'
-
-        CHARACTER (len=*), INTENT(IN) :: axisgroup_id, axis_file
+    SUBROUTINE wxios_add_vaxis(axis_id, axis_size, axis_value)
+        IMPLICIT NONE
+        INCLUDE 'iniprint.h'
+
+        CHARACTER (len=*), INTENT(IN) :: axis_id
         INTEGER, INTENT(IN) :: axis_size
         REAL, DIMENSION(axis_size), INTENT(IN) :: axis_value
         
-        TYPE(xios_axisgroup) :: axgroup
-        TYPE(xios_axis) :: ax
-        CHARACTER(len=20) :: axis_id
-        
-        
-        !Préparation du nom de l'axe:
-        CALL concat(axisgroup_id, axis_file, axis_id)
+!        TYPE(xios_axisgroup) :: axgroup
+!        TYPE(xios_axis) :: ax
+!        CHARACTER(len=50) :: axis_id 
+        
+!        IF (len_trim(axisgroup_id).gt.len(axis_id)) THEN
+!          WRITE(lunout,*) "wxios_add_vaxis: error, size of axis_id too small!!"
+!          WRITE(lunout,*) "     increase it to at least ",len_trim(axisgroup_id)
+!          CALL abort_gcm("wxios_add_vaxis","len(axis_id) too small",1)
+!        ENDIF
+!        axis_id=trim(axisgroup_id)
         
         !On récupère le groupe d'axes qui va bien:
-        CALL xios_get_axisgroup_handle(axisgroup_id, axgroup)
+        !CALL xios_get_axisgroup_handle(axisgroup_id, axgroup)
         
         !On ajoute l'axe correspondant à ce fichier:
-        CALL xios_add_axis(axgroup, ax, TRIM(ADJUSTL(axis_id)))
+        !CALL xios_add_axis(axgroup, ax, TRIM(ADJUSTL(axis_id)))
         
         !Et on le parametrise:
-        CALL xios_set_axis_attr_hdl(ax, size=axis_size, value=axis_value)
+        !CALL xios_set_axis_attr_hdl(ax, size=axis_size, value=axis_value)
+        
+        ! Ehouarn: New way to declare axis, without axis_group:
+        CALL xios_set_axis_attr(trim(axis_id),size=axis_size,value=axis_value)
         
         !Vérification:
@@ -332 +331 @@
             IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_vaxis: Axis created: ", TRIM(ADJUSTL(axis_id))
         ELSE
-            WRITE(*,*) "wxios_add_vaxis: Invalid axis: ", TRIM(ADJUSTL(axis_id))
+            WRITE(lunout,*) "wxios_add_vaxis: Invalid axis: ", TRIM(ADJUSTL(axis_id))
         END IF
 
@@ -367 +366 @@
         
             IF (xios_is_valid_file("X"//fname)) THEN
-                IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_file: New file: ", "X"//fname
-                IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_file: output_freq=",TRIM(ADJUSTL(nffreq)),"; output_lvl=",flvl
+                IF (prt_level >= 10) THEN
+                  WRITE(lunout,*) "wxios_add_file: New file: ", "X"//fname
+                  WRITE(lunout,*) "wxios_add_file: output_freq=",TRIM(ADJUSTL(nffreq)),"; output_lvl=",flvl
+                ENDIF
             ELSE
-                WRITE(*,*) "wxios_add_file: Error, invalid file: ", "X"//trim(fname)
-                WRITE(*,*) "wxios_add_file: output_freq=",TRIM(ADJUSTL(nffreq)),"; output_lvl=",flvl
+                WRITE(lunout,*) "wxios_add_file: Error, invalid file: ", "X"//trim(fname)
+                WRITE(lunout,*) "wxios_add_file: output_freq=",TRIM(ADJUSTL(nffreq)),"; output_lvl=",flvl
             END IF
         ELSE
-            IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_file: File ",trim(fname), " défined using XML."
-                CALL xios_set_file_attr(fname, enabled=.TRUE.)
+            IF (prt_level >= 10) THEN
+              WRITE(lunout,*) "wxios_add_file: File ",trim(fname), " défined using XML."
+            ENDIF
+            ! Ehouarn: add an enable=.true. on top of xml definitions... why???
+            CALL xios_set_file_attr(fname, enabled=.TRUE.)
         END IF
     END SUBROUTINE wxios_add_file
     
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    ! Pour créer un champ      !!!!!!!!!!!!!!!
+    ! Pour créer un champ      !!!!!!!!!!!!!!!!!!!!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     SUBROUTINE wxios_add_field(fieldname, fieldgroup, fieldlongname, fieldunit)
@@ -418 +422 @@
     
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    ! Pour déclarer un champ      !!!!!!!!!!!!!!!
+    ! Pour déclarer un champ      !!!!!!!!!!!!!!!!!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     SUBROUTINE wxios_add_field_to_file(fieldname, fdim, fid, fname, fieldlongname, fieldunit, field_level, op)
@@ -432 +436 @@
         CHARACTER(len=*), INTENT(IN) :: op
         
-        CHARACTER(len=20) :: axis_id
+        CHARACTER(len=20) :: axis_id ! Ehouarn: dangerous...
         CHARACTER(len=100) :: operation
         TYPE(xios_file) :: f
@@ -441 +445 @@
         
         
-        !Préparation du nom de l'axe:
-        CALL concat("presnivs", fname, axis_id)
+        ! Ajout Abd pour NMC:
+        IF (fid.LE.6) THEN
+          axis_id="presnivs"
+        ELSE
+          axis_id="plev"
+        ENDIF
         
         !on prépare le nom de l'opération:
@@ -448 +456 @@
         
         
-        
         !On selectionne le bon groupe de champs:
         IF (fdim.EQ.2) THEN
-            CALL xios_get_fieldgroup_handle("fields_2D", fieldgroup)
+          CALL xios_get_fieldgroup_handle("fields_2D", fieldgroup)
         ELSE
           CALL xios_get_fieldgroup_handle("fields_3D", fieldgroup)
@@ -515 +522 @@
             !Sinon on se contente de l'activer:
             CALL xios_set_field_attr(fieldname, enabled=.TRUE.)
+            !NB: This will override an enable=.false. set by a user in the xml file;
+            !   then the only way to not output the field is by changing its
+            !   output level
         ENDIF       
         
     END SUBROUTINE wxios_add_field_to_file
     
-    SUBROUTINE wxios_update_calendar(ito)
-        INTEGER, INTENT(IN) :: ito
-        CALL xios_update_calendar(ito)
-    END SUBROUTINE wxios_update_calendar
-    
-    SUBROUTINE wxios_write_2D(fieldname, fdata)
-        CHARACTER(len=*), INTENT(IN) :: fieldname
-        REAL, DIMENSION(:,:), INTENT(IN) :: fdata
-        
-        CALL xios_send_field(fieldname, fdata)
-    END SUBROUTINE wxios_write_2D
-    
-    SUBROUTINE wxios_write_3D(fieldname, fdata)
-        CHARACTER(len=*), INTENT(IN) :: fieldname
-        REAL, DIMENSION(:,:,:), INTENT(IN) :: fdata
-        
-        CALL xios_send_field(fieldname, fdata)
-    END SUBROUTINE wxios_write_3D
+!    SUBROUTINE wxios_update_calendar(ito)
+!        INTEGER, INTENT(IN) :: ito
+!        CALL xios_update_calendar(ito)
+!    END SUBROUTINE wxios_update_calendar
+!    
+!    SUBROUTINE wxios_write_2D(fieldname, fdata)
+!        CHARACTER(len=*), INTENT(IN) :: fieldname
+!        REAL, DIMENSION(:,:), INTENT(IN) :: fdata
+!
+!        CALL xios_send_field(fieldname, fdata)
+!    END SUBROUTINE wxios_write_2D
+    
+!    SUBROUTINE wxios_write_3D(fieldname, fdata)
+!        CHARACTER(len=*), INTENT(IN) :: fieldname
+!        REAL, DIMENSION(:,:,:), INTENT(IN) :: fdata
+!        
+!        CALL xios_send_field(fieldname, fdata)
+!    END SUBROUTINE wxios_write_3D
     
     SUBROUTINE wxios_closedef()

trunk/LMDZ.COMMON/libf/dyn3d/calfis.F

@@ -172 +172 @@
       REAL unskap, pksurcp
       save unskap
-
-cIM diagnostique PVteta, Amip2
-      INTEGER,PARAMETER :: ntetaSTD=3
-      REAL,SAVE :: rtetaSTD(ntetaSTD)=(/350.,380.,405./) ! Earth-specific, beware !!
-      REAL PVteta(ngridmx,ntetaSTD)
 
       REAL flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
@@ -651 +646 @@
       ENDDO
 c
-      if (planet_type=="earth") then
-#ifdef CPP_EARTH
-! PVtheta calls tetalevel, which is in the (Earth) physics
-cIM calcul PV a teta=350, 380, 405K
-      CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
-     $           ztfi,zplay,zplev,
-     $           ntetaSTD,rtetaSTD,PVteta)
-#endif
-      endif
-c
 c On change de grille, dynamique vers physiq, pour le flux de masse verticale
       CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
@@ -713 +698 @@
      .             zdqfi,
      .             zdpsrf,
-cIM diagnostique PVteta, Amip2          
-     .             pducov,
-     .             PVteta)
+     .             pducov)
 
       else if ( planet_type=="generic" ) then

trunk/LMDZ.COMMON/libf/dyn3d/ce0l.F90

@@ -1 +1 @@
 !
-! $Id: ce0l.F90 1511 2011-04-28 15:21:47Z jghattas $
+! $Id: ce0l.F90 1984 2014-02-18 09:59:29Z emillour $
 !
 !-------------------------------------------------------------------------------
@@ -30 +30 @@
 #ifndef CPP_EARTH
 #include "iniprint.h"
-  WRITE(lunout,*)'limit_netcdf: Earth-specific program, needs Earth physics'
+  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
 #else
 !-------------------------------------------------------------------------------
@@ -94 +94 @@
   END IF
 
-  IF (grilles_gcm_netcdf) THEN
-     WRITE(lunout,'(//)')
-     WRITE(lunout,*) '  ***************************  '
-     WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
-     WRITE(lunout,*) '  ***************************  '
-     WRITE(lunout,'(//)')
-     CALL grilles_gcm_netcdf_sub(masque,phis)
-  END IF
+  
+  WRITE(lunout,'(//)')
+  WRITE(lunout,*) '  ***************************  '
+  WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
+  WRITE(lunout,*) '  ***************************  '
+  WRITE(lunout,'(//)')
+  CALL grilles_gcm_netcdf_sub(masque,phis)
+
 #endif
 ! of #ifndef CPP_EARTH #else
@@ -108 +108 @@
 !
 !-------------------------------------------------------------------------------
+

trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F

@@ -2 +2 @@
 ! $Id: conf_gcm.F 1418 2010-07-19 15:11:24Z jghattas $
 !
-c
-c
+!
+!
       SUBROUTINE conf_gcm( tapedef, etatinit )
-c
+!
       USE control_mod
 #ifdef CPP_IOIPSL
@@ -18 +18 @@
 
       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
+!-----------------------------------------------------------------------
+!     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 .
+!
        LOGICAL etatinit
        INTEGER tapedef
 
-c   Declarations :
-c   --------------
+!   Declarations :
+!   --------------
 #include "dimensions.h"
 #include "paramet.h"
@@ -43 +43 @@
 ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
 ! #include "clesphys.h"
-c
-c
-c   local:
-c   ------
+!
+!
+!   local:
+!   ------
 
       CHARACTER ch1*72,ch2*72,ch3*72,ch4*12
@@ -54 +54 @@
       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
@@ -91 +91 @@
       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
@@ -103 +103 @@
       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", ...)
@@ -264 +264 @@
        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
@@ -430 +430 @@
        CALL getin('ok_guide',ok_guide)
 
-c    ...............................................................
+!    ...............................................................
 
 !Config  Key  =  read_start
@@ -587 +587 @@
       CALL getin('ok_etat0',ok_etat0)
 
-!Config  Key  = grilles_gcm_netcdf
-!Config  Desc = creation de fichier grilles_gcm.nc dans create_etat0_limit
-!Config  Def  = n
-      grilles_gcm_netcdf = .FALSE.
-      CALL getin('grilles_gcm_netcdf',grilles_gcm_netcdf)
-
-c----------------------------------------
-c Parameters for zonal averages in the case of Titan
+!----------------------------------------
+! Parameters for zonal averages in the case of Titan
       moyzon_mu = .false.
       moyzon_ch = .false.
@@ -601 +595 @@
        CALL getin('moyzon_ch', moyzon_ch)
       endif
-c----------------------------------------
-
-c----------------------------------------
-ccc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
-c     .........   (  modif  le 17/04/96 )   .........
-c
-C ZOOM PARAMETERS ... the ones read in start.nc prevail anyway ! (SL, 2012)
-c
-c----------------------------------------
+!----------------------------------------
+
+!----------------------------------------
+!cc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
+!     .........   (  modif  le 17/04/96 )   .........
+!
+! ZOOM PARAMETERS ... the ones read in start.nc prevail anyway ! (SL, 2012)
+!
+!----------------------------------------
       IF( etatinit ) then
 
@@ -645 +639 @@
 
       IF( grossismx.LT.1. )  THEN
-        write(lunout,*)
-     &   'conf_gcm: ***  ATTENTION !! grossismx < 1 .   *** '
+        write(lunout,*)                                                        &
+     &       'conf_gcm: ***  ATTENTION !! grossismx < 1 .   *** '
          STOP
       ELSE
@@ -654 +648 @@
 
       IF( grossismy.LT.1. )  THEN
-        write(lunout,*)
-     &  'conf_gcm: ***  ATTENTION !! grossismy < 1 .   *** '
+        write(lunout,*)                                                        &
+     &       'conf_gcm: ***  ATTENTION !! grossismy < 1 .   *** '
          STOP
       ELSE
@@ -662 +656 @@
 
       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
@@ -737 +731 @@
 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
@@ -752 +746 @@
 
       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
@@ -766 +760 @@
 
       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
@@ -772 +766 @@
       
       IF( grossismx.LT.1. )  THEN
-        write(lunout,*)
+        write(lunout,*)                                                        &
      &       'conf_gcm: ***  ATTENTION !! grossismx < 1 .   *** '
          STOP
@@ -781 +775 @@
 
       IF( grossismy.LT.1. )  THEN
-        write(lunout,*)
+        write(lunout,*)                                                        &
      &       'conf_gcm: ***  ATTENTION !! grossismy < 1 .   *** '
          STOP
@@ -789 +783 @@
 
       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
@@ -805 +799 @@
          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
@@ -812 +806 @@
          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
@@ -828 +822 @@
       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
@@ -844 +838 @@
       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
@@ -859 +853 @@
       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
@@ -874 +868 @@
       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
@@ -895 +889 @@
           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
@@ -902 +896 @@
           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
@@ -910 +904 @@
 
       endif ! etatinit
-c----------------------------------------
+!----------------------------------------
 
 
@@ -962 +956 @@
       write(lunout,*)' ok_limit = ', ok_limit
       write(lunout,*)' ok_etat0 = ', ok_etat0
-      write(lunout,*)' grilles_gcm_netcdf = ', grilles_gcm_netcdf
       if (planet_type=="titan") then
        write(lunout,*)' moyzon_mu = ', moyzon_mu

trunk/LMDZ.COMMON/libf/dyn3d/gcm.F

@@ -107 +107 @@
       REAL ps(ip1jmp1)                       ! pression  au sol
       REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
-      REAL pks(ip1jmp1)                      ! exner au  sol
-      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
       REAL masse(ip1jmp1,llm)                ! masse d'air
       REAL phis(ip1jmp1)                     ! geopotentiel au sol
@@ -133 +130 @@
       data call_iniphys/.true./
 
-      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
 c+jld variables test conservation energie
 c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
@@ -500 +496 @@
 c   ------------------------
 
-      day_end=day_ini+nday
+      if (nday>=0) then ! standard case
+        day_end=day_ini+nday
+      else ! special case when nday <0, run -nday dynamical steps
+        day_end=day_ini-nday/day_step
+      endif
       if (less1day) then
         day_end=day_ini+floor(time_0+fractday)

trunk/LMDZ.COMMON/libf/dyn3d/guide_mod.F90

@@ -437 +437 @@
 ! 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) CALL guide_out("SP",jjp1,1,ps)
     
     if (guide_u) then
@@ -447 +447 @@
         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("ua",jjp1,llm,(1.-tau)*ugui1+tau*ugui2)
+        IF (f_out) CALL guide_out("u",jjp1,llm,ucov)
+        IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add/factt)
         ucov=ucov+f_add
     endif
@@ -459 +461 @@
         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("teta",jjp1,llm,f_add/factt)
         teta=teta+f_add
     endif
@@ -471 +473 @@
         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("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt)
         ps=ps+f_add(1:ip1jmp1,1)
         CALL pression(ip1jmp1,ap,bp,ps,p)
@@ -485 +487 @@
         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=q+f_add
     endif
@@ -497 +499 @@
         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,vcov)
+        IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)*vgui1+tau*vgui2)
+        IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt)
         vcov=vcov+f_add(1:ip1jm,:)
     endif
@@ -589 +593 @@
   SUBROUTINE guide_interp(psi,teta)
   
+  use exner_hyb_m, only: exner_hyb
+  use exner_milieu_m, only: exner_milieu
   IMPLICIT NONE
 
@@ -610 +616 @@
   REAL, DIMENSION (iip1,jjm,llm)     :: pbary 
   ! 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,llm)    :: pk
   REAL, DIMENSION (iip1,jjp1)        :: pks    
   REAL                               :: prefkap,unskap
@@ -676 +681 @@
     CALL pression( ip1jmp1, ap, bp, psi, p )
     if (pressure_exner) then
-      CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
+      CALL exner_hyb(ip1jmp1,psi,p,pks,pk)
     else
-      CALL exner_milieu(ip1jmp1,psi,p,beta,pks,pk,pkf)
+      CALL exner_milieu(ip1jmp1,psi,p,pks,pk)
     endif
 !    ....  Calcul de pls , pression au milieu des couches ,en Pascals
@@ -1507 +1512 @@
     
     ! Variables entree
-    CHARACTER, INTENT(IN)                          :: varname
+    CHARACTER*(*), INTENT(IN)                          :: varname
     INTEGER,   INTENT (IN)                         :: hsize,vsize
     REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
@@ -1516 +1521 @@
     INTEGER       :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
     INTEGER       :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
+    INTEGER       :: vid_au,vid_av
     INTEGER, DIMENSION (3) :: dim3
     INTEGER, DIMENSION (4) :: dim4,count,start
-    INTEGER                :: ierr, varid
+    INTEGER                :: ierr, varid,l
+    REAL, DIMENSION (iip1,hsize,vsize) :: field2
 
     print *,'Guide: output timestep',timestep,'var ',varname
@@ -1542 +1549 @@
         ierr=NF_DEF_VAR(nid,"LEVEL",NF_FLOAT,1,id_lev,vid_lev)
         ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
+        ierr=NF_DEF_VAR(nid,"au",NF_FLOAT,2,(/id_lonu,id_latu/),vid_au)
         ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
+        ierr=NF_DEF_VAR(nid,"av",NF_FLOAT,2,(/id_lonv,id_latv/),vid_av)
         
         ierr=NF_ENDDEF(nid)
@@ -1555 +1564 @@
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
+        ierr = NF_PUT_VAR_DOUBLE(nid,vid_au,alpha_u)
+        ierr = NF_PUT_VAR_DOUBLE(nid,vid_av,alpha_v)
 #else
         ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
@@ -1563 +1574 @@
         ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
         ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
+        ierr = NF_PUT_VAR_REAL(nid,vid_au,alpha_u)
+        ierr = NF_PUT_VAR_REAL(nid,vid_av,alpha_v)
 #endif
 ! --------------------------------------------------------------------
@@ -1579 +1592 @@
         IF (guide_u) THEN
             dim4=(/id_lonu,id_latu,id_lev,id_tim/)
+            ierr = NF_DEF_VAR(nid,"u",NF_FLOAT,4,dim4,varid)
+            ierr = NF_DEF_VAR(nid,"ua",NF_FLOAT,4,dim4,varid)
             ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
         ENDIF
@@ -1584 +1599 @@
         IF (guide_v) THEN
             dim4=(/id_lonv,id_latv,id_lev,id_tim/)
+            ierr = NF_DEF_VAR(nid,"v",NF_FLOAT,4,dim4,varid)
+            ierr = NF_DEF_VAR(nid,"va",NF_FLOAT,4,dim4,varid)
             ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
         ENDIF
@@ -1606 +1623 @@
     ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
 
+    IF (varname=="SP") timestep=timestep+1
+
+    ierr = NF_INQ_VARID(nid,varname,varid)
     SELECT CASE (varname)
-    CASE ("S")
-        timestep=timestep+1
-        ierr = NF_INQ_VARID(nid,"SP",varid)
+    CASE ("SP","ps")
         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)
-#endif
-    CASE ("P")
-        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)
-#endif
-    CASE ("U")
-        ierr = NF_INQ_VARID(nid,"ucov",varid)
+    CASE ("v","va","vcov")
+        start=(/1,1,1,timestep/)
+        count=(/iip1,jjm,llm,1/)
+    CASE DEFAULT
         start=(/1,1,1,timestep/)
         count=(/iip1,jjp1,llm,1/)
+    END SELECT
+
+    SELECT CASE (varname)
+    CASE("u","ua")
+        DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO
+        field2(:,1,:)=0. ; field2(:,jjp1,:)=0.
+    CASE("v","va")
+        DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO
+    CASE DEFAULT
+        field2=field
+    END SELECT
+
+
 #ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
+    ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field2)
 #else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
+    ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field2)
 #endif
-    CASE ("V")
-        ierr = NF_INQ_VARID(nid,"vcov",varid)
-        start=(/1,1,1,timestep/)
-        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)
-#endif
-    CASE ("T")
-        ierr = NF_INQ_VARID(nid,"teta",varid)
-        start=(/1,1,1,timestep/)
-        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)
-#endif
-    CASE ("Q")
-        ierr = NF_INQ_VARID(nid,"q",varid)
-        start=(/1,1,1,timestep/)
-        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)
-#endif
-    END SELECT
- 
+
     ierr = NF_CLOSE(nid)
 

trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F

@@ -20 +20 @@
      &                       ok_dynzon,periodav,ok_dyn_ave,iecri,
      &                       ok_dyn_ins,output_grads_dyn
+      use exner_hyb_m, only: exner_hyb
+      use exner_milieu_m, only: exner_milieu
       use cpdet_mod, only: cpdet,tpot2t,t2tpot
       use sponge_mod, only: callsponge,mode_sponge,sponge
@@ -217 +219 @@
       endif
 
-      itaufin   = nday*day_step
+      if (nday>=0) then
+         itaufin   = nday*day_step
+      else
+         ! to run a given (-nday) number of dynamical steps
+         itaufin   = -nday
+      endif
       if (less1day) then
 c MODIF VENUS: to run less than one day:
@@ -262 +269 @@
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
       if (pressure_exner) then
-        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
       else
-        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+        CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
       endif
 
@@ -476 +483 @@
          CALL pression (  ip1jmp1, ap, bp, ps,  p      )
          if (pressure_exner) then
-           CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
+           CALL exner_hyb(  ip1jmp1, ps, p,pks, pk, pkf )
          else
-           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+           CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
          endif
+
+! Compute geopotential (physics might need it)
+         CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
 
            jD_cur = jD_ref + day_ini - day_ref +                        &
@@ -551 +561 @@
           CALL massdair(p,masse)
           if (pressure_exner) then
-            CALL exner_hyb(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+            CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
           else
-            CALL exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
+            CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
           endif
           
@@ -604 +614 @@
         CALL pression ( ip1jmp1, ap, bp, ps, p                  )
         if (pressure_exner) then
-          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+          CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
         else
-          CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+          CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
         endif
         CALL massdair(p,masse)

trunk/LMDZ.COMMON/libf/dyn3d/logic.h

@@ -11 +11 @@
      &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
      &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
-     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf,hybrid                    &
+     &  ,ok_limit,ok_etat0,hybrid                                       &
      &  ,moyzon_mu,moyzon_ch
 
@@ -19 +19 @@
      & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
      &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
-     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf
+     &  ,ok_limit,ok_etat0
+
       logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
                      ! (only used if disvert_type==2)

trunk/LMDZ.COMMON/libf/dyn3d_common/disvert.F90

@@ -1 +1 @@
 ! $Id: disvert.F90 1645 2012-07-30 16:01:50Z lguez $
 
-SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
+SUBROUTINE disvert()
 
   ! Auteur : P. Le Van
 
+#ifdef CPP_IOIPSL
+  use ioipsl, only: getin
+#else
+  USE ioipsl_getincom, only: getin
+#endif
   use new_unit_m, only: new_unit
-  use ioipsl, only: getin
   use assert_m, only: assert
 
@@ -13 +17 @@
   include "dimensions.h"
   include "paramet.h"
+  include "comvert.h"
+  include "comconst.h"
   include "iniprint.h"
   include "logic.h"
 
-  ! s = sigma ** kappa : coordonnee verticale
-  ! dsig(l) : epaisseur de la couche l ds la coord. s
-  ! sig(l) : sigma a l'interface des couches l et l-1
-  ! ds(l) : distance entre les couches l et l-1 en coord.s
-
-  real,intent(in) :: pa, preff
-  real,intent(out) :: ap(llmp1) ! in Pa
-  real, intent(out):: bp(llmp1)
-  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
-  real,intent(out) :: presnivs(llm)
-  real,intent(out) :: scaleheight
-
+!-------------------------------------------------------------------------------
+! Purpose: Vertical distribution functions for LMDZ.
+!          Triggered by the levels number llm.
+!-------------------------------------------------------------------------------
+! Read    in "comvert.h":
+! pa                         !--- PURE PRESSURE COORDINATE FOR P<pa (in Pascals)
+! preff                      !--- REFERENCE PRESSURE                 (101325 Pa)
+! Written in "comvert.h":
+! ap(llm+1), bp(llm+1)       !--- Ap, Bp HYBRID COEFFICIENTS AT INTERFACES
+! aps(llm),  bps(llm)        !--- Ap, Bp HYBRID COEFFICIENTS AT MID-LAYERS
+! dpres(llm)                 !--- PRESSURE DIFFERENCE FOR EACH LAYER
+! presnivs(llm)              !--- PRESSURE AT EACH MID-LAYER
+! scaleheight                !--- VERTICAL SCALE HEIGHT            (Earth: 8kms)
+! nivsig(llm+1)              !--- SIGMA INDEX OF EACH LAYER INTERFACE
+! nivsigs(llm)               !--- SIGMA INDEX OF EACH MID-LAYER
+!-------------------------------------------------------------------------------
+! Local variables:
   REAL sig(llm+1), dsig(llm)
-  real zk, zkm1, dzk1, dzk2, k0, k1
+  REAL sig0(llm+1), zz(llm+1)
+  REAL zk, zkm1, dzk1, dzk2, z, k0, k1
 
   INTEGER l, unit
   REAL dsigmin
+  REAL vert_scale,vert_dzmin,vert_dzlow,vert_z0low,vert_dzmid,vert_z0mid,vert_h_mid,vert_dzhig,vert_z0hig,vert_h_hig
+
   REAL alpha, beta, deltaz
   REAL x
   character(len=*),parameter :: modname="disvert"
 
-  character(len=6):: vert_sampling
+  character(len=24):: vert_sampling
   ! (allowed values are "param", "tropo", "strato" and "read")
 
   !-----------------------------------------------------------------------
 
-  print *, "Call sequence information: disvert"
+  WRITE(lunout,*) TRIM(modname)//" starts"
 
   ! default scaleheight is 8km for earth
@@ -49 +63 @@
   vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
   call getin('vert_sampling', vert_sampling)
-  print *, 'vert_sampling = ' // vert_sampling
+  WRITE(lunout,*) TRIM(modname)//' vert_sampling = ' // vert_sampling
+  if (llm==39 .and. vert_sampling=="strato") then
+     dsigmin=0.3 ! Vieille option par défaut pour CMIP5
+  else
+     dsigmin=1.
+  endif
+  call getin('dsigmin', dsigmin)
+  WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
+
 
   select case (vert_sampling)
@@ -86 +108 @@
 
      ap = pa * (sig - bp)
-  case("tropo")
+  case("sigma")
      DO l = 1, llm
         x = 2*asin(1.) * (l - 0.5) / (llm + 1)
-        dsig(l) = 1.0 + 7.0 * SIN(x)**2
+        dsig(l) = dsigmin + 7.0 * SIN(x)**2
      ENDDO
      dsig = dsig / sum(dsig)
@@ -98 +120 @@
 
      bp(1)=1.
+     bp(2: llm) = sig(2:llm)
+     bp(llmp1) = 0.
+     ap(:)=0.
+  case("tropo")
+     DO l = 1, llm
+        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
+        dsig(l) = dsigmin + 7.0 * SIN(x)**2
+     ENDDO
+     dsig = dsig / sum(dsig)
+     sig(llm+1) = 0.
+     DO l = llm, 1, -1
+        sig(l) = sig(l+1) + dsig(l)
+     ENDDO
+
+     bp(1)=1.
      bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
      bp(llmp1) = 0.
@@ -104 +141 @@
      ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
   case("strato")
-     if (llm==39) then
-        dsigmin=0.3
-     else if (llm==50) then
-        dsigmin=1.
-     else
-        write(lunout,*) trim(modname), ' ATTENTION discretisation z a ajuster'
-        dsigmin=1.
-     endif
-     WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
-
      DO l = 1, llm
         x = 2*asin(1.) * (l - 0.5) / (llm + 1)
@@ -131 +158 @@
      ap(1)=0.
      ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
+  case("strato_correct")
+!==================================================================
+! Fredho 2014/05/18, Saint-Louis du Senegal
+! Cette version de la discretisation strato est corrige au niveau
+! du passage des sig aux ap, bp
+! la version precedente donne un coude dans l'epaisseur des couches
+! vers la tropopause
+!==================================================================
+
+
+     DO l = 1, llm
+        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
+        dsig(l) =(dsigmin + 7. * SIN(x)**2) &
+             *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
+     ENDDO
+     dsig = dsig / sum(dsig)
+     sig0(llm+1) = 0.
+     DO l = llm, 1, -1
+        sig0(l) = sig0(l+1) + dsig(l)
+     ENDDO
+     sig=racinesig(sig0)
+
+     bp(1)=1.
+     bp(2:llm)=EXP(1.-1./sig(2: llm)**2)
+     bp(llmp1)=0.
+
+     ap(1)=0.
+     ap(2:llm)=pa*(sig(2:llm)-bp(2:llm))
+     ap(llm+1)=0.
+
+  CASE("strato_custom0") 
+!=======================================================
+! Version Transitoire
+    ! custumize strato distribution with specific alpha & beta values and function
+    ! depending on llm (experimental and temporary)!
+    SELECT CASE (llm)
+      CASE(55)
+        alpha=0.45
+        beta=4.0
+      CASE(63)
+        alpha=0.45
+        beta=5.0
+      CASE(71)
+        alpha=3.05
+        beta=65.
+      CASE(79)
+        alpha=3.20
+        ! alpha=2.05 ! FLOTT 79 (PLANTE)
+        beta=70.
+    END SELECT
+    ! Or used values provided by user in def file:
+    CALL getin("strato_alpha",alpha)
+    CALL getin("strato_beta",beta)
+    
+    ! Build geometrical distribution
+    scaleheight=7.
+    zz(1)=0.
+    IF (llm==55.OR.llm==63) THEN
+      DO l=1,llm
+        z=zz(l)/scaleheight
+        zz(l+1)=zz(l)+0.03+z*1.5*(1.-TANH(z-0.5))+alpha*(1.+TANH(z-1.5))     &
+                            +5.0*EXP((l-llm)/beta)
+      ENDDO
+    ELSEIF (llm==71) THEN !.OR.llm==79) THEN      ! FLOTT 79 (PLANTE)
+      DO l=1,llm
+        z=zz(l)
+        zz(l+1)=zz(l)+0.02+0.88*TANH(z/2.5)+alpha*(1.+TANH((z-beta)/15.))
+      ENDDO
+    ELSEIF (llm==79) THEN
+      DO l=1,llm
+        z=zz(l)
+        zz(l+1)=zz(l)+0.02+0.80*TANH(z/3.8)+alpha*(1+TANH((z-beta)/17.))     &
+                            +0.03*TANH(z/.25)
+      ENDDO
+    ENDIF ! of IF (llm==55.OR.llm==63) ...
+    
+
+    ! Build sigma distribution
+    sig0=EXP(-zz(:)/scaleheight)
+    sig0(llm+1)=0.
+!    sig=ridders(sig0)
+    sig=racinesig(sig0)
+    
+    ! Compute ap() and bp()
+    bp(1)=1.
+    bp(2:llm)=EXP(1.-1./sig(2:llm)**2)
+    bp(llm+1)=0.
+    ap=pa*(sig-bp)
+
+  CASE("strato_custom") 
+!===================================================================
+! David Cugnet, François Lott, Lionel Guez, Ehouoarn Millour, Fredho
+! 2014/05
+! custumize strato distribution
+! Al the parameter are given in km assuming a given scalehigh
+    vert_scale=7.     ! scale hight
+    vert_dzmin=0.02   ! width of first layer
+    vert_dzlow=1.     ! dz in the low atmosphere
+    vert_z0low=8.     ! height at which resolution recches dzlow
+    vert_dzmid=3.     ! dz in the mid atmsophere
+    vert_z0mid=70.    ! height at which resolution recches dzmid
+    vert_h_mid=20.    ! width of the transition
+    vert_dzhig=11.    ! dz in the high atmsophere
+    vert_z0hig=80.    ! height at which resolution recches dz
+    vert_h_hig=20.    ! width of the transition
+!===================================================================
+
+    call getin('vert_scale',vert_scale)
+    call getin('vert_dzmin',vert_dzmin)
+    call getin('vert_dzlow',vert_dzlow)
+    call getin('vert_z0low',vert_z0low)
+    CALL getin('vert_dzmid',vert_dzmid)
+    CALL getin('vert_z0mid',vert_z0mid)
+    call getin('vert_h_mid',vert_h_mid)
+    call getin('vert_dzhig',vert_dzhig)
+    call getin('vert_z0hig',vert_z0hig)
+    call getin('vert_h_hig',vert_h_hig)
+
+    scaleheight=vert_scale ! for consistency with further computations
+    ! Build geometrical distribution
+    zz(1)=0.
+    DO l=1,llm
+       z=zz(l)
+       zz(l+1)=zz(l)+vert_dzmin+vert_dzlow*TANH(z/vert_z0low)+                &
+&      (vert_dzmid-vert_dzlow)* &
+&           (TANH((z-vert_z0mid)/vert_h_mid)-TANH((-vert_z0mid)/vert_h_mid)) &
+&     +(vert_dzhig-vert_dzmid-vert_dzlow)*                                  &
+&           (TANH((z-vert_z0hig)/vert_h_hig)-TANH((-vert_z0hig)/vert_h_hig))
+    ENDDO
+
+
+!===================================================================
+! Comment added Fredho 2014/05/18, Saint-Louis, Senegal
+! From approximate z to ap, bp, so that p=ap+bp*p0 and p/p0=exp(-z/H)
+! sig0 is p/p0
+! sig is an intermediate distribution introduce to estimate bp
+! 1.   sig0=exp(-z/H)
+! 2.   inversion of sig0=(1-pa/p0)*sig+(1-pa/p0)*exp(1-1/sig**2)
+! 3.   bp=exp(1-1/sig**2)
+! 4.   ap deduced from  the combination of 2 and 3 so that sig0=ap/p0+bp
+!===================================================================
+
+    sig0=EXP(-zz(:)/vert_scale)
+    sig0(llm+1)=0.
+    sig=racinesig(sig0)
+    bp(1)=1.
+    bp(2:llm)=EXP(1.-1./sig(2:llm)**2)
+    bp(llm+1)=0.
+    ap=pa*(sig-bp)
+
   case("read")
      ! Read "ap" and "bp". First line is skipped (title line). "ap"
@@ -178 +355 @@
   write(lunout, *) presnivs
 
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+FUNCTION ridders(sig) RESULT(sg)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Purpose: Search for s solving (Pa/Preff)*s+(1-Pa/Preff)*EXP(1-1./s**2)=sg
+! Notes:   Uses Ridders' method, quite robust. Initial bracketing: 0<=sg<=1.
+! Reference: Ridders, C. F. J. "A New Algorithm for Computing a Single Root of a
+!       Real Continuous Function" IEEE Trans. Circuits Systems 26, 979-980, 1979
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)  :: sig(:)
+  REAL              :: sg(SIZE(sig))
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: it, ns, maxit
+  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
+!-------------------------------------------------------------------------------
+  ns=SIZE(sig); maxit=9999
+  c1=Pa/Preff; c2=1.-c1
+  DO l=1,ns
+    xx=HUGE(1.)
+    x1=0.0; f1=distrib(x1,c1,c2,sig(l))
+    x2=1.0; f2=distrib(x2,c1,c2,sig(l))
+    DO it=1,maxit
+      x3=0.5*(x1+x2); f3=distrib(x3,c1,c2,sig(l))
+      s=SQRT(f3**2-f1*f2);                 IF(s==0.) EXIT
+      x4=x3+(x3-x1)*(SIGN(1.,f1-f2)*f3/s); IF(ABS(10.*LOG(x4-xx))<=1E-5) EXIT
+      xx=x4; f4=distrib(x4,c1,c2,sig(l));  IF(f4==0.) EXIT
+      IF(SIGN(f3,f4)/=f3) THEN;      x1=x3; f1=f3; x2=xx; f2=f4
+      ELSE IF(SIGN(f1,f4)/=f1) THEN; x2=xx; f2=f4
+      ELSE IF(SIGN(f2,f4)/=f2) THEN; x1=xx; f1=f4
+      ELSE; CALL abort_gcm("ridders",'Algorithm failed (which is odd...')
+      END IF
+      IF(ABS(10.*LOG(ABS(x2-x1)))<=1E-5) EXIT       !--- ERROR ON SIG <= 0.01m            
+    END DO
+    IF(it==maxit+1) WRITE(lunout,'(a,i3)')'WARNING in ridder: failed to converg&
+     &e for level ',l
+    sg(l)=xx
+  END DO
+  sg(1)=1.; sg(ns)=0.
+
+END FUNCTION ridders
+
+FUNCTION racinesig(sig) RESULT(sg)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Fredho 2014/05/18
+! Purpose: Search for s solving (Pa/Preff)*sg+(1-Pa/Preff)*EXP(1-1./sg**2)=s
+! Notes:   Uses Newton Raphson search
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)  :: sig(:)
+  REAL              :: sg(SIZE(sig))
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: it, ns, maxit
+  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
+!-------------------------------------------------------------------------------
+  ns=SIZE(sig); maxit=100
+  c1=Pa/Preff; c2=1.-c1
+  DO l=2,ns-1
+    sg(l)=sig(l)
+    DO it=1,maxit
+       f1=exp(1-1./sg(l)**2)*(1.-c1)
+       sg(l)=sg(l)-(c1*sg(l)+f1-sig(l))/(c1+2*f1*sg(l)**(-3))
+    ENDDO
+!   print*,'SSSSIG ',sig(l),sg(l),c1*sg(l)+exp(1-1./sg(l)**2)*(1.-c1)
+  ENDDO
+  sg(1)=1.; sg(ns)=0.
+
+END FUNCTION racinesig
+
+
+
+
 END SUBROUTINE disvert
+
+!-------------------------------------------------------------------------------
+
+FUNCTION distrib(x,c1,c2,x0) RESULT(res)
+!
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN) :: x, c1, c2, x0
+  REAL             :: res
+!-------------------------------------------------------------------------------
+  res=c1*x+c2*EXP(1-1/(x**2))-x0
+
+END FUNCTION distrib
+
+
+

trunk/LMDZ.COMMON/libf/dyn3d_common/exner_hyb_m.F90

@@ -1 @@
-!
-! $Id $
-!
-      SUBROUTINE  exner_hyb ( ngrid, ps, p,alpha,beta, pks, pk, pkf )
-c
-c     Auteurs :  P.Le Van  , Fr. Hourdin  .
-c    ..........
-c
-c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
-c    .... alpha,beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
-c
-c   ************************************************************************
-c    Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des 
-c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
-c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
-c   ************************************************************************
-c  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
-c    la pression et la fonction d'Exner  au  sol  .
-c
-c                                 -------- z                                   
-c    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
-c                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
-c    ( voir note de Fr.Hourdin )  ,
-c
-c    on determine successivement , du haut vers le bas des couches, les 
-c    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
-c    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
-c     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
-c
-c
-      IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-#include "comvert.h"
-#include "serre.h"
+module exner_hyb_m
 
-      INTEGER  ngrid
-      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
-      REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
+  IMPLICIT NONE
 
-c    .... variables locales   ...
+contains
 
-      INTEGER l, ij
-      REAL unpl2k,dellta
+  SUBROUTINE  exner_hyb ( ngrid, ps, p, pks, pk, pkf )
 
-      REAL ppn(iim),pps(iim)
-      REAL xpn, xps
-      REAL SSUM
-c
-      logical,save :: firstcall=.true.
-      character(len=*),parameter :: modname="exner_hyb"
-      
-      ! Sanity check
-      if (firstcall) then
-        ! sanity checks for Shallow Water case (1 vertical layer)
-        if (llm.eq.1) then
+    !     Auteurs :  P.Le Van  , Fr. Hourdin  .
+    !    ..........
+    !
+    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
+    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
+    !
+    !   ************************************************************************
+    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
+    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
+    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
+    !   ************************************************************************
+    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
+    !    la pression et la fonction d'Exner  au  sol  .
+    !
+    !                                 -------- z
+    !    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
+    !                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
+    !    ( voir note de Fr.Hourdin )  ,
+    !
+    !    on determine successivement , du haut vers le bas des couches, les 
+    !    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
+    !    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
+    !     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
+    !
+    !
+    !
+    include "dimensions.h"
+    include "paramet.h"
+    include "comconst.h"
+    include "comgeom.h"
+    include "comvert.h"
+    include "serre.h"
+
+    INTEGER  ngrid
+    REAL p(ngrid,llmp1),pk(ngrid,llm)
+    real, optional:: pkf(ngrid,llm)
+    REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
+
+    !    .... variables locales   ...
+
+    INTEGER l, ij
+    REAL unpl2k,dellta
+
+    logical,save :: firstcall=.true.
+    character(len=*),parameter :: modname="exner_hyb"
+
+    ! Sanity check
+    if (firstcall) then
+       ! 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)
+             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)
+             call abort_gcm(modname, &
+                  "cpp!=r , but running in Shallow Water mode!!",42)
           endif
-        endif ! of if (llm.eq.1)
+       endif ! of if (llm.eq.1)
 
-        firstcall=.false.
-      endif ! of if (firstcall)
+       firstcall=.false.
+    endif ! of if (firstcall)
 
-      if (llm.eq.1) then
-        
-        ! Compute pks(:),pk(:),pkf(:)
-        
-        DO   ij  = 1, ngrid
-          pks(ij) = (cpp/preff) * ps(ij) 
+    ! Specific behaviour for Shallow Water (1 vertical layer) case:
+    if (llm.eq.1) then
+
+       ! Compute pks(:),pk(:),pkf(:)
+
+       DO   ij  = 1, ngrid
+          pks(ij) = (cpp/preff) * ps(ij)
           pk(ij,1) = .5*pks(ij)
-        ENDDO
-        
-        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
-        
-        ! our work is done, exit routine
-        return
+       ENDDO
 
-      endif ! of if (llm.eq.1)
+       if (present(pkf)) then
+          pkf = pk
+          CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
+       end if
 
-!!!! General case:
-     
-      unpl2k    = 1.+ 2.* kappa
-c
-      DO   ij  = 1, ngrid
-        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-      ENDDO
+       ! our work is done, exit routine
+       return
+    endif ! of if (llm.eq.1)
 
-      DO  ij   = 1, iim
-        ppn(ij) = aire(   ij   ) * pks(  ij     )
-        pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-      ENDDO
-      xpn      = SSUM(iim,ppn,1) /apoln
-      xps      = SSUM(iim,pps,1) /apols
+    ! General case:
 
-      DO ij   = 1, iip1
-        pks(   ij     )  =  xpn
-        pks( ij+ip1jm )  =  xps
-      ENDDO
-c
-c
-c    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
-c
-      DO     ij      = 1, ngrid
+    unpl2k    = 1.+ 2.* kappa
+
+    !     -------------
+    !     Calcul de pks
+    !     -------------
+
+    DO   ij  = 1, ngrid
+       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+    ENDDO
+
+    !    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
+    !
+    DO     ij      = 1, ngrid
        alpha(ij,llm) = 0.
        beta (ij,llm) = 1./ unpl2k
-      ENDDO
-c
-c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
-c
-      DO l = llm -1 , 2 , -1
-c
-        DO ij = 1, ngrid
-        dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
-        alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
-        beta (ij,l)  =   p(ij,l  ) / dellta   
-        ENDDO
-c
-      ENDDO
-c
-c  ***********************************************************************
-c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
-c
-      DO   ij   = 1, ngrid
-       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  /
-     *    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
-      ENDDO
-c
-c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
-c
-      DO l = 2, llm
-        DO   ij   = 1, ngrid
-         pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
-        ENDDO
-      ENDDO
-c
-c
-      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-      CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
-      
+    ENDDO
+    !
+    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
+    !
+    DO l = llm -1 , 2 , -1
+       !
+       DO ij = 1, ngrid
+          dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
+          alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
+          beta (ij,l)  =   p(ij,l  ) / dellta   
+       ENDDO
+    ENDDO
 
-      RETURN
-      END
+    !  ***********************************************************************
+    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....
+    !
+    DO   ij   = 1, ngrid
+       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  / &
+            (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
+    ENDDO
+    !
+    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
+    !
+    DO l = 2, llm
+       DO   ij   = 1, ngrid
+          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
+       ENDDO
+    ENDDO
+
+    if (present(pkf)) then
+       !    calcul de pkf
+       pkf = pk
+       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
+    end if
+
+  END SUBROUTINE exner_hyb
+
+end module exner_hyb_m
+

trunk/LMDZ.COMMON/libf/dyn3d_common/exner_milieu_m.F90

@@ -1 @@
-!
-! $Id $
-!
-      SUBROUTINE  exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf )
-c
-c     Auteurs :  F. Forget , Y. Wanherdrick
-c P.Le Van  , Fr. Hourdin  .
-c    ..........
-c
-c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
-c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
-c
-c   ************************************************************************
-c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
-c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
-c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
-c   ************************************************************************
-c    .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
-c    la pression et la fonction d'Exner  au  sol  .
-c
-c     WARNING : CECI est une version speciale de exner_hyb originale
-c               Utilise dans la version martienne pour pouvoir 
-c               tourner avec des coordonnees verticales complexe
-c              => Il ne verifie PAS la condition la proportionalite en 
-c              energie totale/ interne / potentielle (F.Forget 2001)
-c    ( voir note de Fr.Hourdin )  ,
-c
-      IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-#include "comvert.h"
-#include "serre.h"
+module exner_milieu_m
 
-      INTEGER  ngrid
-      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
-      REAL ps(ngrid),pks(ngrid), beta(ngrid,llm)
+  IMPLICIT NONE
 
-c    .... variables locales   ...
+contains
 
-      INTEGER l, ij
-      REAL dum1
+  SUBROUTINE  exner_milieu ( ngrid, ps, p, pks, pk, pkf )
+    !
+    !     Auteurs :  F. Forget , Y. Wanherdrick
+    ! P.Le Van  , Fr. Hourdin  .
+    !    ..........
+    !
+    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
+    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
+    !
+    !   ************************************************************************
+    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
+    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
+    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
+    !   ************************************************************************
+    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
+    !    la pression et la fonction d'Exner  au  sol  .
+    !
+    !     WARNING : CECI est une version speciale de exner_hyb originale
+    !               Utilise dans la version martienne pour pouvoir 
+    !               tourner avec des coordonnees verticales complexe
+    !              => Il ne verifie PAS la condition la proportionalite en 
+    !              energie totale/ interne / potentielle (F.Forget 2001)
+    !    ( voir note de Fr.Hourdin )  ,
+    !
+    !
+    include "dimensions.h"
+    include "paramet.h"
+    include "comconst.h"
+    include "comgeom.h"
+    include "comvert.h"
+    include "serre.h"
 
-      REAL ppn(iim),pps(iim)
-      REAL xpn, xps
-      REAL SSUM
-      EXTERNAL SSUM
-      logical,save :: firstcall=.true.
-      character(len=*),parameter :: modname="exner_milieu"
+    INTEGER  ngrid
+    REAL p(ngrid,llmp1),pk(ngrid,llm)
+    real, optional:: pkf(ngrid,llm)
+    REAL ps(ngrid),pks(ngrid)
 
-      ! Sanity check
-      if (firstcall) then
-        ! sanity checks for Shallow Water case (1 vertical layer)
-        if (llm.eq.1) then
+    !    .... variables locales   ...
+
+    INTEGER l, ij
+    REAL dum1
+
+    logical,save :: firstcall=.true.
+    character(len=*),parameter :: modname="exner_milieu"
+
+    ! Sanity check
+    if (firstcall) then
+       ! 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)
+             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)
+             call abort_gcm(modname, &
+                  "cpp!=r , but running in Shallow Water mode!!",42)
           endif
-        endif ! of if (llm.eq.1)
+       endif ! of if (llm.eq.1)
 
-        firstcall=.false.
-      endif ! of if (firstcall)
+       firstcall=.false.
+    endif ! of if (firstcall)
 
-!!!! Specific behaviour for Shallow Water (1 vertical layer) case:
-      if (llm.eq.1) then
-      
-        ! Compute pks(:),pk(:),pkf(:)
-        
-        DO   ij  = 1, ngrid
-          pks(ij) = (cpp/preff) * ps(ij) 
+    ! Specific behaviour for Shallow Water (1 vertical layer) case:
+    if (llm.eq.1) then
+
+       ! Compute pks(:),pk(:),pkf(:)
+
+       DO   ij  = 1, ngrid
+          pks(ij) = (cpp/preff) * ps(ij)
           pk(ij,1) = .5*pks(ij)
-        ENDDO
-        
-        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
-        
-        ! our work is done, exit routine
-        return
+       ENDDO
 
-      endif ! of if (llm.eq.1)
+       if (present(pkf)) then
+          pkf = pk
+          CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
+       end if
 
-!!!! General case:
+       ! our work is done, exit routine
+       return
+    endif ! of if (llm.eq.1)
 
-c     -------------
-c     Calcul de pks
-c     -------------
-   
-      DO   ij  = 1, ngrid
-        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-      ENDDO
+    ! General case:
 
-      DO  ij   = 1, iim
-        ppn(ij) = aire(   ij   ) * pks(  ij     )
-        pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-      ENDDO
-      xpn      = SSUM(iim,ppn,1) /apoln
-      xps      = SSUM(iim,pps,1) /apols
+    !     -------------
+    !     Calcul de pks
+    !     -------------
 
-      DO ij   = 1, iip1
-        pks(   ij     )  =  xpn
-        pks( ij+ip1jm )  =  xps
-      ENDDO
-c
-c
-c    .... Calcul de pk  pour la couche l 
-c    --------------------------------------------
-c
-      dum1 = cpp * (2*preff)**(-kappa) 
-      DO l = 1, llm-1
-        DO   ij   = 1, ngrid
-         pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
-        ENDDO
-      ENDDO
+    DO   ij  = 1, ngrid
+       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+    ENDDO
 
-c    .... Calcul de pk  pour la couche l = llm ..
-c    (on met la meme distance (en log pression)  entre Pk(llm)
-c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
+    !    .... Calcul de pk  pour la couche l 
+    !    --------------------------------------------
+    !
+    dum1 = cpp * (2*preff)**(-kappa) 
+    DO l = 1, llm-1
+       DO   ij   = 1, ngrid
+          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
+       ENDDO
+    ENDDO
 
-      DO   ij   = 1, ngrid
-         pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
-      ENDDO
+    !    .... Calcul de pk  pour la couche l = llm ..
+    !    (on met la meme distance (en log pression)  entre Pk(llm)
+    !    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
 
+    DO   ij   = 1, ngrid
+       pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
+    ENDDO
 
-c    calcul de pkf
-c    -------------
-      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-      CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
-      
-c    EST-CE UTILE ?? : calcul de beta
-c    --------------------------------
-      DO l = 2, llm
-        DO   ij   = 1, ngrid
-          beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
-        ENDDO
-      ENDDO
+    if (present(pkf)) then
+       !    calcul de pkf
+       pkf = pk
+       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
+    end if
 
-      RETURN
-      END
+  END SUBROUTINE exner_milieu
+
+end module exner_milieu_m
+

trunk/LMDZ.COMMON/libf/dyn3d_common/iniacademic.F90

@@ -14 +14 @@
 #endif
   USE Write_Field
+  use exner_hyb_m, only: exner_hyb
+  use exner_milieu_m, only: exner_milieu
 
   !   Author:    Frederic Hourdin      original: 15/01/93
@@ -54 +56 @@
   REAL pks(ip1jmp1)                      ! exner au  sol
   REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-  REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
   REAL phi(ip1jmp1,llm)                  ! geopotentiel
   REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
@@ -223 +224 @@
         CALL pression ( ip1jmp1, ap, bp, ps, p       )
         if (pressure_exner) then
-          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-        else
-          call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
+          CALL exner_hyb( ip1jmp1, ps, p, pks, pk)
+        else
+          call exner_milieu(ip1jmp1,ps,p,pks,pk)
         endif
         CALL massdair(p,masse)

trunk/LMDZ.COMMON/libf/dyn3d_common/iniconst.F90

@@ -73 +73 @@
   if (disvert_type==1) then
      ! standard case for Earth (automatic generation of levels)
-     call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig, scaleheight)
+     call disvert()
   else if (disvert_type==2) then
      ! standard case for planets (levels generated using z2sig.def file)

trunk/LMDZ.COMMON/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
+

trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F

@@ -225 +225 @@
       save unskap
 
-cIM diagnostique PVteta, Amip2
-      INTEGER,PARAMETER :: ntetaSTD=3
-      REAL,SAVE :: rtetaSTD(ntetaSTD)=(/350.,380.,405./) ! Earth-specific, beware !!
-      REAL PVteta(klon,ntetaSTD)
-      
-      
       REAL SSUM
 
@@ -267 +261 @@
       klon=klon_mpi
       
-      PVteta(:,:)=0.
-            
       IF ( firstcal )  THEN
         debut = .TRUE.
@@ -684 +676 @@
       endif
 
-
-      IF (is_sequential.and.(planet_type=="earth")) THEN
-#ifdef CPP_EARTH
-! PVtheta calls tetalevel, which is in the physics
-cIM calcul PV a teta=350, 380, 405K
-        CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
-     $           ztfi,zplay,zplev,
-     $           ntetaSTD,rtetaSTD,PVteta)
-c
-#endif
-      ENDIF
-
 c On change de grille, dynamique vers physiq, pour le flux de masse verticale
       CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
@@ -923 +903 @@
      .             zdqfi_omp,
      .             zdpsrf_omp,
-cIM diagnostique PVteta, Amip2          
-     .             pducov,
-     .             PVteta)
+     .             pducov)
 
       else if ( planet_type=="generic" ) then

trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90

@@ -1 +1 @@
 !
-! $Id: ce0l.F90 1615 2012-02-10 15:42:26Z emillour $
+! $Id: ce0l.F90 1984 2014-02-18 09:59:29Z emillour $
 !
 !-------------------------------------------------------------------------------
@@ -115 +115 @@
   END IF
 
-  IF (grilles_gcm_netcdf) THEN
-     WRITE(lunout,'(//)')
-     WRITE(lunout,*) '  ***************************  '
-     WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
-     WRITE(lunout,*) '  ***************************  '
-     WRITE(lunout,'(//)')
-     CALL grilles_gcm_netcdf_sub(masque,phis)
-  END IF
+  WRITE(lunout,'(//)')
+  WRITE(lunout,*) '  ***************************  '
+  WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
+  WRITE(lunout,*) '  ***************************  '
+  WRITE(lunout,'(//)')
+  CALL grilles_gcm_netcdf_sub(masque,phis)
   
 #ifdef CPP_MPI
@@ -137 +135 @@
 !
 !-------------------------------------------------------------------------------
+

trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F

@@ -2 +2 @@
 ! $Id: conf_gcm.F 1438 2010-10-08 10:19:34Z jghattas $
 !
-c
-c
+!
+!
       SUBROUTINE conf_gcm( tapedef, etatinit )
-c
+!
 #ifdef CPP_IOIPSL
       use IOIPSL
@@ -20 +20 @@
       use sponge_mod_p, only: callsponge,mode_sponge,nsponge,tetasponge
       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
+!-----------------------------------------------------------------------
+!     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 .
+!
        LOGICAL etatinit
        INTEGER tapedef
 
-c   Declarations :
-c   --------------
+!   Declarations :
+!   --------------
 #include "dimensions.h"
 #include "paramet.h"
@@ -45 +45 @@
 ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
 ! #include "clesphys.h"
-c
-c
-c   local:
-c   ------
+!
+!
+!   local:
+!   ------
 
       CHARACTER ch1*72,ch2*72,ch3*72,ch4*12
@@ -60 +60 @@
       integer,external :: OMP_GET_NUM_THREADS
 #endif      
-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
@@ -290 +290 @@
        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
@@ -456 +456 @@
        CALL getin('ok_guide',ok_guide)
 
-c    ...............................................................
+!     ...............................................................
 
 !Config  Key  =  read_start
@@ -632 +632 @@
       CALL getin('ok_etat0',ok_etat0)
 
-!Config  Key  = grilles_gcm_netcdf
-!Config  Desc = creation de fichier grilles_gcm.nc dans create_etat0_limit
-!Config  Def  = n
-      grilles_gcm_netcdf = .FALSE.
-      CALL getin('grilles_gcm_netcdf',grilles_gcm_netcdf)
-
-c----------------------------------------
-c Parameters for zonal averages in the case of Titan
+!----------------------------------------
+! Parameters for zonal averages in the case of Titan
       moyzon_mu = .false.
       moyzon_ch = .false.
@@ -646 +640 @@
        CALL getin('moyzon_ch', moyzon_ch)
       endif
-c----------------------------------------
-
-c----------------------------------------
-ccc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
-c     .........   (  modif  le 17/04/96 )   .........
-c
-C ZOOM PARAMETERS ... the ones read in start.nc prevail anyway ! (SL, 2012)
-c
-c----------------------------------------
+!----------------------------------------
+
+!----------------------------------------
+!cc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
+!     .........   (  modif  le 17/04/96 )   .........
+!
+! ZOOM PARAMETERS ... the ones read in start.nc prevail anyway ! (SL, 2012)
+!
+!----------------------------------------
       IF( etatinit ) then
 
@@ -707 +701 @@
 
       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
@@ -758 +752 @@
        ysinus = .TRUE.
        CALL getin('ysinus',ysinus)
-c
-c----------------------------------------
+!
+!----------------------------------------
        else ! etatinit=false
-c----------------------------------------
+!----------------------------------------
 
 !Config  Key  = clon
@@ -779 +773 @@
        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',
@@ -834 +828 @@
 
       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
@@ -862 +856 @@
          ENDIF
       ENDIF
-c
+!
 !Config  Key  = dzoomx
 !Config  Desc = extension en longitude
@@ -925 +919 @@
       ENDIF
 
-cc
+!c
       IF( .NOT.fxyhypb  )  THEN
 
@@ -955 +949 @@
 
       endif ! etatinit
-c----------------------------------------
+!----------------------------------------
 
 
@@ -1009 +1003 @@
       write(lunout,*)' ok_limit = ', ok_limit
       write(lunout,*)' ok_etat0 = ', ok_etat0
-      write(lunout,*)' grilles_gcm_netcdf = ', grilles_gcm_netcdf
       if (planet_type=="titan") then
        write(lunout,*)' moyzon_mu = ', moyzon_mu

trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p_m.F90

@@ -1 @@
-!
-! $Id $
-!
-      SUBROUTINE  exner_hyb_p ( ngrid, ps, p,alpha,beta, pks, pk, pkf )
-c
-c     Auteurs :  P.Le Van  , Fr. Hourdin  .
-c    ..........
-c
-c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
-c    .... alpha,beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
-c
-c   ************************************************************************
-c    Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des 
-c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
-c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
-c   ************************************************************************
-c  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
-c    la pression et la fonction d'Exner  au  sol  .
-c
-c                                 -------- z                                   
-c    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
-c                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
-c    ( voir note de Fr.Hourdin )  ,
-c
-c    on determine successivement , du haut vers le bas des couches, les 
-c    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
-c    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
-c     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
-c
-c
-      USE parallel_lmdz
-      IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-#include "comvert.h"
-#include "serre.h"
+module exner_hyb_p_m
 
-      INTEGER  ngrid
-      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
-      REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
+  IMPLICIT NONE
 
-c    .... variables locales   ...
+contains
 
-      INTEGER l, ij
-      REAL unpl2k,dellta
+  SUBROUTINE  exner_hyb_p ( ngrid, ps, p, pks, pk, pkf )
 
-      REAL ppn(iim),pps(iim)
-      REAL xpn, xps
-      REAL SSUM
-      EXTERNAL SSUM
-      INTEGER ije,ijb,jje,jjb
-      logical,save :: firstcall=.true.
-!$OMP THREADPRIVATE(firstcall) 
-      character(len=*),parameter :: modname="exner_hyb_p"
-c
+    !     Auteurs :  P.Le Van  , Fr. Hourdin  .
+    !    ..........
+    !
+    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
+    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
+    !
+    !   ************************************************************************
+    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
+    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
+    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
+    !   ************************************************************************
+    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
+    !    la pression et la fonction d'Exner  au  sol  .
+    !
+    !                                 -------- z
+    !    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
+    !                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
+    !    ( voir note de Fr.Hourdin )  ,
+    !
+    !    on determine successivement , du haut vers le bas des couches, les 
+    !    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
+    !    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
+    !     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
+    !
+    !
+    USE parallel_lmdz
+    !
+    include "dimensions.h"
+    include "paramet.h"
+    include "comconst.h"
+    include "comgeom.h"
+    include "comvert.h"
+    include "serre.h"
 
-      ! Sanity check
-      if (firstcall) then
-        ! sanity checks for Shallow Water case (1 vertical layer)
-        if (llm.eq.1) then
+    INTEGER  ngrid
+    REAL p(ngrid,llmp1),pk(ngrid,llm)
+    REAL, optional:: pkf(ngrid,llm)
+    REAL ps(ngrid),pks(ngrid)
+    REAL alpha(ngrid,llm),beta(ngrid,llm)
+
+    !    .... variables locales   ...
+
+    INTEGER l, ij
+    REAL unpl2k,dellta
+
+    INTEGER ije,ijb,jje,jjb
+    logical,save :: firstcall=.true.
+    !$OMP THREADPRIVATE(firstcall) 
+    character(len=*),parameter :: modname="exner_hyb_p"
+
+    ! Sanity check
+    if (firstcall) then
+       ! 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)
+             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)
+             call abort_gcm(modname, &
+                  "cpp!=r , but running in Shallow Water mode!!",42)
           endif
-        endif ! of if (llm.eq.1)
+       endif ! of if (llm.eq.1)
 
-        firstcall=.false.
-      endif ! of if (firstcall)
+       firstcall=.false.
+    endif ! of if (firstcall)
 
-c$OMP BARRIER
+    !$OMP BARRIER
 
-! Specific behaviour for Shallow Water (1 vertical layer) case
-      if (llm.eq.1) then
-      
-        ! Compute pks(:),pk(:),pkf(:)
-        ijb=ij_begin
-        ije=ij_end
-!$OMP DO SCHEDULE(STATIC)
-        DO ij=ijb, ije
-          pks(ij)=(cpp/preff)*ps(ij)
+    ! Specific behaviour for Shallow Water (1 vertical layer) case:
+    if (llm.eq.1) then
+
+       ! Compute pks(:),pk(:),pkf(:)
+       ijb=ij_begin
+       ije=ij_end
+       !$OMP DO SCHEDULE(STATIC)
+       DO ij=ijb, ije
+          pks(ij) = (cpp/preff) * ps(ij)
           pk(ij,1) = .5*pks(ij)
-          pkf(ij,1)=pk(ij,1)
-        ENDDO
-!$OMP ENDDO
+          if (present(pkf)) pkf(ij,1)=pk(ij,1)
+       ENDDO
+       !$OMP ENDDO
 
-!$OMP MASTER
-      if (pole_nord) then
-        DO  ij   = 1, iim
-          ppn(ij) = aire(   ij   ) * pks(  ij     )
-        ENDDO
-        xpn      = SSUM(iim,ppn,1) /apoln
-  
-        DO ij   = 1, iip1
-          pks(   ij     )  =  xpn
-          pk(ij,1) = .5*pks(ij)
-          pkf(ij,1)=pk(ij,1)
-        ENDDO
-      endif
-      
-      if (pole_sud) then
-        DO  ij   = 1, iim
-          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-        ENDDO
-        xps      = SSUM(iim,pps,1) /apols 
-  
-        DO ij   = 1, iip1
-          pks( ij+ip1jm )  =  xps
-          pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
-          pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
-        ENDDO
-      endif
-!$OMP END MASTER
-!$OMP BARRIER
-        jjb=jj_begin
-        jje=jj_end
-        CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+       !$OMP BARRIER
+       if (present(pkf)) then
+          jjb=jj_begin
+          jje=jj_end
+          CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+       end if
 
-        ! our work is done, exit routine
-        return
-      endif ! of if (llm.eq.1)
+       ! our work is done, exit routine
+       return
+    endif ! of if (llm.eq.1)
 
-!!!! General case:
+    ! General case:
 
-      unpl2k    = 1.+ 2.* kappa
-c
-      ijb=ij_begin
-      ije=ij_end
+    unpl2k    = 1.+ 2.* kappa
 
-c$OMP DO SCHEDULE(STATIC)
-      DO   ij  = ijb, ije
-        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-      ENDDO
-c$OMP ENDDO
-c Synchro OPENMP ici
+    !     -------------
+    !     Calcul de pks
+    !     -------------
 
-c$OMP MASTER
-      if (pole_nord) then
-        DO  ij   = 1, iim
-          ppn(ij) = aire(   ij   ) * pks(  ij     )
-        ENDDO
-        xpn      = SSUM(iim,ppn,1) /apoln
-  
-        DO ij   = 1, iip1
-          pks(   ij     )  =  xpn
-        ENDDO
-      endif
-      
-      if (pole_sud) then
-        DO  ij   = 1, iim
-          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-        ENDDO
-        xps      = SSUM(iim,pps,1) /apols 
-  
-        DO ij   = 1, iip1
-          pks( ij+ip1jm )  =  xps
-        ENDDO
-      endif
-c$OMP END MASTER
-c$OMP BARRIER
-c
-c
-c    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
-c
-c$OMP DO SCHEDULE(STATIC)
-      DO     ij      = ijb,ije
+    ijb=ij_begin
+    ije=ij_end
+
+    !$OMP DO SCHEDULE(STATIC)
+    DO   ij  = ijb, ije
+       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+    ENDDO
+    !$OMP ENDDO
+    ! Synchro OPENMP ici
+
+    !$OMP BARRIER
+    !
+    !
+    !    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
+    !
+    !$OMP DO SCHEDULE(STATIC)
+    DO     ij      = ijb,ije
        alpha(ij,llm) = 0.
        beta (ij,llm) = 1./ unpl2k
-      ENDDO
-c$OMP ENDDO NOWAIT
-c
-c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
-c
-      DO l = llm -1 , 2 , -1
-c
-c$OMP DO SCHEDULE(STATIC)
-        DO ij = ijb, ije
-        dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
-        alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
-        beta (ij,l)  =   p(ij,l  ) / dellta   
-        ENDDO
-c$OMP ENDDO NOWAIT
-c
-      ENDDO
+    ENDDO
+    !$OMP ENDDO NOWAIT
+    !
+    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
+    !
+    DO l = llm -1 , 2 , -1
+       !
+       !$OMP DO SCHEDULE(STATIC)
+       DO ij = ijb, ije
+          dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
+          alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
+          beta (ij,l)  =   p(ij,l  ) / dellta   
+       ENDDO
+       !$OMP ENDDO NOWAIT
+    ENDDO
 
-c
-c  ***********************************************************************
-c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
-c
-c$OMP DO SCHEDULE(STATIC)
-      DO   ij   = ijb, ije
-       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  /
-     *    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
-      ENDDO
-c$OMP ENDDO NOWAIT
-c
-c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
-c
-      DO l = 2, llm
-c$OMP DO SCHEDULE(STATIC)
-        DO   ij   = ijb, ije
-         pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
-        ENDDO
-c$OMP ENDDO NOWAIT        
-      ENDDO
-c
-c
-c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-      DO l = 1, llm
-c$OMP DO SCHEDULE(STATIC)
-         DO   ij   = ijb, ije
-           pkf(ij,l)=pk(ij,l)
-         ENDDO
-c$OMP ENDDO NOWAIT             
-      ENDDO
+    !  ***********************************************************************
+    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....
+    !
+    !$OMP DO SCHEDULE(STATIC)
+    DO   ij   = ijb, ije
+       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  / &
+            (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
+    ENDDO
+    !$OMP ENDDO NOWAIT
+    !
+    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
+    !
+    DO l = 2, llm
+       !$OMP DO SCHEDULE(STATIC)
+       DO   ij   = ijb, ije
+          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
+       ENDDO
+       !$OMP ENDDO NOWAIT        
+    ENDDO
 
-c$OMP BARRIER
-      
-      jjb=jj_begin
-      jje=jj_end
-      CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
-      
+    if (present(pkf)) then
+       !    calcul de pkf
 
-      RETURN
-      END
+       DO l = 1, llm
+          !$OMP DO SCHEDULE(STATIC)
+          DO   ij   = ijb, ije
+             pkf(ij,l)=pk(ij,l)
+          ENDDO
+          !$OMP ENDDO NOWAIT             
+       ENDDO
+
+       !$OMP BARRIER
+
+       jjb=jj_begin
+       jje=jj_end
+       CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+    end if
+
+  END SUBROUTINE exner_hyb_p
+
+end module exner_hyb_p_m
+

trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p_m.F90

@@ -1 @@
-!
-! $Id $
-!
-      SUBROUTINE  exner_milieu_p ( ngrid, ps, p,beta, pks, pk, pkf )
-c
-c     Auteurs :  F. Forget , Y. Wanherdrick
-c P.Le Van  , Fr. Hourdin  .
-c    ..........
-c
-c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
-c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
-c
-c   ************************************************************************
-c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
-c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
-c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
-c   ************************************************************************
-c    .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
-c    la pression et la fonction d'Exner  au  sol  .
-c
-c     WARNING : CECI est une version speciale de exner_hyb originale
-c               Utilise dans la version martienne pour pouvoir 
-c               tourner avec des coordonnees verticales complexe
-c              => Il ne verifie PAS la condition la proportionalite en 
-c              energie totale/ interne / potentielle (F.Forget 2001)
-c    ( voir note de Fr.Hourdin )  ,
-c
-      USE parallel_lmdz
-      IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-#include "comvert.h"
-#include "serre.h"
+module exner_milieu_p_m
 
-      INTEGER  ngrid
-      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
-      REAL ps(ngrid),pks(ngrid), beta(ngrid,llm)
+  IMPLICIT NONE
 
-c    .... variables locales   ...
+contains
 
-      INTEGER l, ij
-      REAL dum1
+  SUBROUTINE  exner_milieu_p ( ngrid, ps, p, pks, pk, pkf )
+    !
+    !     Auteurs :  F. Forget , Y. Wanherdrick
+    ! P.Le Van  , Fr. Hourdin  .
+    !    ..........
+    !
+    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
+    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
+    !
+    !   ************************************************************************
+    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
+    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
+    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
+    !   ************************************************************************
+    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
+    !    la pression et la fonction d'Exner  au  sol  .
+    !
+    !     WARNING : CECI est une version speciale de exner_hyb originale
+    !               Utilise dans la version martienne pour pouvoir 
+    !               tourner avec des coordonnees verticales complexe
+    !              => Il ne verifie PAS la condition la proportionalite en 
+    !              energie totale/ interne / potentielle (F.Forget 2001)
+    !    ( voir note de Fr.Hourdin )  ,
+    !
+    USE parallel_lmdz
+    !
+    include "dimensions.h"
+    include "paramet.h"
+    include "comconst.h"
+    include "comgeom.h"
+    include "comvert.h"
+    include "serre.h"
 
-      REAL ppn(iim),pps(iim)
-      REAL xpn, xps
-      REAL SSUM
-      EXTERNAL SSUM
-      INTEGER ije,ijb,jje,jjb
-      logical,save :: firstcall=.true.
-!$OMP THREADPRIVATE(firstcall) 
-      character(len=*),parameter :: modname="exner_milieu_p"
+    INTEGER  ngrid
+    REAL p(ngrid,llmp1),pk(ngrid,llm)
+    REAL, optional:: pkf(ngrid,llm)
+    REAL ps(ngrid),pks(ngrid)
 
-      ! Sanity check
-      if (firstcall) then
-        ! sanity checks for Shallow Water case (1 vertical layer)
-        if (llm.eq.1) then
+    !    .... variables locales   ...
+
+    INTEGER l, ij,ijb,ije,jjb,jje
+    REAL dum1
+
+    logical,save :: firstcall=.true.
+    !$OMP THREADPRIVATE(firstcall) 
+    character(len=*),parameter :: modname="exner_milieu_p"
+
+    ! Sanity check
+    if (firstcall) then
+       ! 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)
+             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)
+             call abort_gcm(modname, &
+                  "cpp!=r , but running in Shallow Water mode!!",42)
           endif
-        endif ! of if (llm.eq.1)
+       endif ! of if (llm.eq.1)
 
-        firstcall=.false.
-      endif ! of if (firstcall)
-      
-c$OMP BARRIER
+       firstcall=.false.
+    endif ! of if (firstcall)
 
-! Specific behaviour for Shallow Water (1 vertical layer) case
-      if (llm.eq.1) then
-              
-        ! Compute pks(:),pk(:),pkf(:)
-        ijb=ij_begin
-        ije=ij_end
-!$OMP DO SCHEDULE(STATIC)
-        DO ij=ijb, ije
-          pks(ij)=(cpp/preff)*ps(ij)
+    !$OMP BARRIER
+
+    ! Specific behaviour for Shallow Water (1 vertical layer) case:
+    if (llm.eq.1) then
+
+       ! Compute pks(:),pk(:),pkf(:)
+       ijb=ij_begin
+       ije=ij_end
+       !$OMP DO SCHEDULE(STATIC)
+       DO ij=ijb, ije
+          pks(ij) = (cpp/preff) * ps(ij)
           pk(ij,1) = .5*pks(ij)
-          pkf(ij,1)=pk(ij,1)
-        ENDDO
-!$OMP ENDDO
+          if (present(pkf)) pkf(ij,1)=pk(ij,1)
+       ENDDO
+       !$OMP ENDDO
 
-!$OMP MASTER
-      if (pole_nord) then
-        DO  ij   = 1, iim
-          ppn(ij) = aire(   ij   ) * pks(  ij     )
-        ENDDO
-        xpn      = SSUM(iim,ppn,1) /apoln
-  
-        DO ij   = 1, iip1
-          pks(   ij     )  =  xpn
-          pk(ij,1) = .5*pks(ij)
-          pkf(ij,1)=pk(ij,1)
-        ENDDO
-      endif
-      
-      if (pole_sud) then
-        DO  ij   = 1, iim
-          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-        ENDDO
-        xps      = SSUM(iim,pps,1) /apols 
-  
-        DO ij   = 1, iip1
-          pks( ij+ip1jm )  =  xps
-          pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
-          pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
-        ENDDO
-      endif
-!$OMP END MASTER
-!$OMP BARRIER
-        jjb=jj_begin
-        jje=jj_end
-        CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+       !$OMP BARRIER
+       if (present(pkf)) then
+          jjb=jj_begin
+          jje=jj_end
+          CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+       end if
 
-        ! our work is done, exit routine
-        return
-      endif ! of if (llm.eq.1)
+       ! our work is done, exit routine
+       return
+    endif ! of if (llm.eq.1)
 
-!!!! General case:
+    ! General case:
 
-c     -------------
-c     Calcul de pks
-c     -------------
-   
-      ijb=ij_begin
-      ije=ij_end
+    !     -------------
+    !     Calcul de pks
+    !     -------------
 
-c$OMP DO SCHEDULE(STATIC)
-      DO   ij  = ijb, ije
-        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-      ENDDO
-c$OMP ENDDO
-c Synchro OPENMP ici
+    ijb=ij_begin
+    ije=ij_end
 
-c$OMP MASTER
-      if (pole_nord) then
-        DO  ij   = 1, iim
-          ppn(ij) = aire(   ij   ) * pks(  ij     )
-        ENDDO
-        xpn      = SSUM(iim,ppn,1) /apoln
-  
-        DO ij   = 1, iip1
-          pks(   ij     )  =  xpn
-        ENDDO
-      endif
-      
-      if (pole_sud) then
-        DO  ij   = 1, iim
-          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-        ENDDO
-        xps      = SSUM(iim,pps,1) /apols 
-  
-        DO ij   = 1, iip1
-          pks( ij+ip1jm )  =  xps
-        ENDDO
-      endif
-c$OMP END MASTER
-c$OMP BARRIER
-c
-c
-c    .... Calcul de pk  pour la couche l 
-c    --------------------------------------------
-c
-      dum1 = cpp * (2*preff)**(-kappa) 
-      DO l = 1, llm-1
-c$OMP DO SCHEDULE(STATIC)
-        DO   ij   = ijb, ije
-         pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
-        ENDDO
-c$OMP ENDDO NOWAIT        
-      ENDDO
+    !$OMP DO SCHEDULE(STATIC)
+    DO   ij  = ijb, ije
+       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+    ENDDO
+    !$OMP ENDDO
+    ! Synchro OPENMP ici
 
-c    .... Calcul de pk  pour la couche l = llm ..
-c    (on met la meme distance (en log pression)  entre Pk(llm)
-c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
+    !$OMP BARRIER
+    !
+    !
+    !    .... Calcul de pk  pour la couche l 
+    !    --------------------------------------------
+    !
+    dum1 = cpp * (2*preff)**(-kappa) 
+    DO l = 1, llm-1
+       !$OMP DO SCHEDULE(STATIC)
+       DO   ij   = ijb, ije
+          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
+       ENDDO
+       !$OMP ENDDO NOWAIT
+    ENDDO
 
-c$OMP DO SCHEDULE(STATIC)
-      DO   ij   = ijb, ije
-         pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
-      ENDDO
-c$OMP ENDDO NOWAIT        
+    !    .... Calcul de pk  pour la couche l = llm ..
+    !    (on met la meme distance (en log pression)  entre Pk(llm)
+    !    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
 
+    !$OMP DO SCHEDULE(STATIC)
+    DO   ij   = ijb, ije
+       pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
+    ENDDO
+    !$OMP ENDDO NOWAIT        
 
-c    calcul de pkf
-c    -------------
-c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-      DO l = 1, llm
-c$OMP DO SCHEDULE(STATIC)
-         DO   ij   = ijb, ije
-           pkf(ij,l)=pk(ij,l)
-         ENDDO
-c$OMP ENDDO NOWAIT             
-      ENDDO
+    if (present(pkf)) then
+       !    calcul de pkf
 
-c$OMP BARRIER
-      
-      jjb=jj_begin
-      jje=jj_end
-      CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
-      
-c    EST-CE UTILE ?? : calcul de beta
-c    --------------------------------
-      DO l = 2, llm
-c$OMP DO SCHEDULE(STATIC)
-        DO   ij   = ijb, ije
-          beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
-        ENDDO
-c$OMP ENDDO NOWAIT             
-      ENDDO
+       DO l = 1, llm
+          !$OMP DO SCHEDULE(STATIC)
+          DO   ij   = ijb, ije
+             pkf(ij,l)=pk(ij,l)
+          ENDDO
+          !$OMP ENDDO NOWAIT
+       ENDDO
 
-      RETURN
-      END
+       !$OMP BARRIER
+
+       jjb=jj_begin
+       jje=jj_end
+       CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+    end if
+
+  END SUBROUTINE exner_milieu_p
+
+end module exner_milieu_p_m
+

trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F

@@ -516 +516 @@
 
 
-      day_end = day_ini + nday
+      if (nday>=0) then ! standard case
+        day_end=day_ini+nday
+      else ! special case when nday <0, run -nday dynamical steps
+        day_end=day_ini-nday/day_step
+      endif
       if (less1day) then
         day_end=day_ini+floor(time_0+fractday)

trunk/LMDZ.COMMON/libf/dyn3dpar/guide_p_mod.F90

@@ -328 +328 @@
 !=======================================================================
   SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
+    use exner_hyb_p_m, only: exner_hyb_p
+    use exner_milieu_p_m, only: exner_milieu_p
     USE parallel_lmdz
     USE control_mod
@@ -349 +351 @@
     REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
     ! 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,llm)    :: pk
     REAL, DIMENSION (iip1,jjp1)        :: pks    
     REAL                               :: unskap
@@ -491 +492 @@
     IF (f_out) THEN
 !       Calcul niveaux pression milieu de couches 
-        CALL pression_p( ip1jmp1, ap, bp, ps, p )
-        if (pressure_exner) 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)
+        CALL pression_p( ip1jmp1, ap, bp, ps, p )
+        if (pressure_exner) then
+          CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk)
+        else
+          CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk)
         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.)
+        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("SP",jjp1,llm,p(1:ip1jmp1,1:llm),1.)
     ENDIF
     
@@ -517 +518 @@
         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,ucov,factt)
+        IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)*ugui1(:,:)+tau*ugui2(:,:),factt)
+        IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add(:,:)/factt,factt)
         ucov(ijb_u:ije_u,:)=ucov(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
     endif
@@ -529 +532 @@
         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("teta",jjp1,llm,f_add(:,:)/factt,factt)
         teta(ijb_u:ije_u,:)=teta(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
     endif
@@ -541 +544 @@
         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("SP",jjp1,1,f_add(1:ip1jmp1,1),factt)
+        IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt,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)
@@ -555 +558 @@
         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,factt)
         q(ijb_u:ije_u,:)=q(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
     endif
@@ -568 +571 @@
         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,vcov(1:ip1jm,:),factt)
+        IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)*vgui1(:,:)+tau*vgui2(:,:),factt)
+        IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt,factt)
         vcov(ijb_v:ije_v,:)=vcov(ijb_v:ije_v,:)+f_add(ijb_v:ije_v,:)
     endif
@@ -689 +694 @@
 !=======================================================================
   SUBROUTINE guide_interp(psi,teta)
+    use exner_hyb_p_m, only: exner_hyb_p
+    use exner_milieu_p_m, only: exner_milieu_p
   USE parallel_lmdz
   USE mod_hallo
@@ -713 +720 @@
   REAL, DIMENSION (iip1,jjm,llm)     :: pbary 
   ! 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,llm)    :: pk
   REAL, DIMENSION (iip1,jjp1)        :: pks    
   REAL                               :: unskap
@@ -784 +790 @@
     IF (guide_plevs.EQ.1) THEN
         DO l=1,llm
-            DO j=jjb_u,jje_u
-                DO i =1, iip1
+            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
         ENDDO
     ELSE
-        CALL pression_p( ip1jmp1, ap, bp, psi, p )
-        if (pressure_exner) then
-          CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
+        CALL pression_p( ip1jmp1, ap, bp, psi, p )
+        if (pressure_exner) then
+          CALL exner_hyb_p(ip1jmp1,psi,p,pks,pk)
         else
-          CALL exner_milieu_p(ip1jmp1,psi,p,beta,pks,pk,pkf)
+          CALL exner_milieu_p(ip1jmp1,psi,p,pks,pk)
         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
+        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
 
@@ -1024 +1030 @@
         ! 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)
+         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
@@ -1041 +1047 @@
             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)
+         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
@@ -1807 +1813 @@
     
     ! Variables entree
-    CHARACTER, INTENT(IN)                          :: varname
+    CHARACTER*(*), INTENT(IN)                          :: varname
     INTEGER,   INTENT (IN)                         :: hsize,vsize
     REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
@@ -1817 +1823 @@
     INTEGER       :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
     INTEGER       :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
+    INTEGER       :: vid_au,vid_av
+    INTEGER       :: l
     INTEGER, DIMENSION (3) :: dim3
     INTEGER, DIMENSION (4) :: dim4,count,start
     INTEGER                :: ierr, varid
+    REAL, DIMENSION (iip1,hsize,vsize) :: field2
     
     CALL gather_field(field,iip1*hsize,vsize,0)
@@ -1834 +1843 @@
 ! Definition des dimensions
         ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu) 
+        print*,'id_lonu 1 ',id_lonu
         ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv) 
         ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu) 
@@ -1842 +1852 @@
 ! Creation des variables dimensions
         ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)
+        print*,'id_lonu 2 ',id_lonu
         ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)
         ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)
@@ -1848 +1859 @@
         ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
         ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
+        ierr=NF_DEF_VAR(nid,"au",NF_FLOAT,2,(/id_lonu,id_latu/),vid_au)
+        ierr=NF_DEF_VAR(nid,"av",NF_FLOAT,2,(/id_lonv,id_latv/),vid_av)
         
         ierr=NF_ENDDEF(nid)
@@ -1853 +1866 @@
 ! Enregistrement des variables dimensions
 #ifdef NC_DOUBLE
+        print*,'id_lonu DOUBLE ',id_lonu,rlonu*180./pi
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi)
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi)
@@ -1860 +1874 @@
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
+        ierr = NF_PUT_VAR_DOUBLE(nid,vid_au,alpha_u)
+        ierr = NF_PUT_VAR_DOUBLE(nid,vid_av,alpha_v)
 #else
+        print*,'id_lonu 3 ',id_lonu,rlonu*180./pi
         ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
         ierr = NF_PUT_VAR_REAL(nid,vid_lonv,rlonv*180./pi)
@@ -1868 +1885 @@
         ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
         ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
+        ierr = NF_PUT_VAR_REAL(nid,vid_au,alpha_u)
+        ierr = NF_PUT_VAR_REAL(nid,vid_av,alpha_v)
 #endif
 ! --------------------------------------------------------------------
@@ -1875 +1894 @@
 ! Pressure (GCM)
         dim4=(/id_lonv,id_latu,id_lev,id_tim/)
-        ierr = NF_DEF_VAR(nid,"P",NF_FLOAT,4,dim4,varid)
+        ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,4,dim4,varid)
 ! Surface pressure (guidage)
         IF (guide_P) THEN
@@ -1883 +1902 @@
 ! Zonal wind
         IF (guide_u) THEN
+        print*,'id_lonu 4 ',id_lonu,varname
             dim4=(/id_lonu,id_latu,id_lev,id_tim/)
+            ierr = NF_DEF_VAR(nid,"u",NF_FLOAT,4,dim4,varid)
+            ierr = NF_DEF_VAR(nid,"ua",NF_FLOAT,4,dim4,varid)
             ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
         ENDIF
@@ -1889 +1911 @@
         IF (guide_v) THEN
             dim4=(/id_lonv,id_latv,id_lev,id_tim/)
+            ierr = NF_DEF_VAR(nid,"v",NF_FLOAT,4,dim4,varid)
+            ierr = NF_DEF_VAR(nid,"va",NF_FLOAT,4,dim4,varid)
             ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
         ENDIF
@@ -1912 +1936 @@
     ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
 
+    IF (varname=="SP") timestep=timestep+1
+
+    IF (varname=="SP") THEN
+      print*,'varname=SP=',varname
+    ELSE
+      print*,'varname diff SP=',varname
+    ENDIF
+
+
+    ierr = NF_INQ_VARID(nid,varname,varid)
     SELECT CASE (varname)
-    CASE ("P")
-        timestep=timestep+1
-        ierr = NF_INQ_VARID(nid,"P",varid)
+    CASE ("SP","ps")
         start=(/1,1,1,timestep/)
         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)
-#endif
-    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/factt)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
-#endif
-    CASE ("U")
-        ierr = NF_INQ_VARID(nid,"ucov",varid)
+    CASE ("v","va","vcov")
+        start=(/1,1,1,timestep/)
+        count=(/iip1,jjm,llm,1/)
+    CASE DEFAULT
         start=(/1,1,1,timestep/)
         count=(/iip1,jjp1,llm,1/)
+    END SELECT
+
+    SELECT CASE (varname)
+    CASE("u","ua")
+        DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO
+        field2(:,1,:)=0. ; field2(:,jjp1,:)=0.
+    CASE("v","va")
+        DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO
+    CASE DEFAULT
+        field2=field
+    END SELECT
+
 #ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
+    ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field2)
 #else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
+    ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field2)
 #endif
-    CASE ("V")
-        ierr = NF_INQ_VARID(nid,"vcov",varid)
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjm,llm,1/)
-#ifdef NC_DOUBLE
-        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")
-        ierr = NF_INQ_VARID(nid,"teta",varid)
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjp1,llm,1/)
-#ifdef NC_DOUBLE
-        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")
-        ierr = NF_INQ_VARID(nid,"q",varid)
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjp1,llm,1/)
-#ifdef NC_DOUBLE
-        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
  
     ierr = NF_CLOSE(nid)

trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F

@@ -8 +8 @@
      &                    time_0)
 
+      use exner_hyb_m, only: exner_hyb
+      use exner_milieu_m, only: exner_milieu
+      use exner_hyb_p_m, only: exner_hyb_p
+      use exner_milieu_p_m, only: exner_milieu_p
        USE misc_mod
        USE parallel_lmdz
@@ -17 +21 @@
        USE vampir
        USE timer_filtre, ONLY : print_filtre_timer
-       USE infotrac
+       USE infotrac, ONLY: nqtot
        USE guide_p_mod, ONLY : guide_main
        USE getparam
@@ -165 +169 @@
       character*10 string10
 
-      REAL,SAVE :: alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
       REAL,SAVE :: flxw(ip1jmp1,llm) ! flux de masse verticale
 
@@ -236 +239 @@
       lafin=.false.
       
-      itaufin   = nday*day_step
+      if (nday>=0) then
+         itaufin   = nday*day_step
+      else
+         ! to run a given (-nday) number of dynamical steps
+         itaufin   = -nday
+      endif
       if (less1day) then
 c MODIF VENUS: to run less than one day:
@@ -292 +300 @@
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
       if (pressure_exner) then
-        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
       else
-        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+        CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
       endif
 c$OMP END MASTER
@@ -830 +838 @@
 c$OMP BARRIER
          if (pressure_exner) then
-           CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
+           CALL exner_hyb_p(  ip1jmp1, ps, p,pks, pk, pkf )
          else
-           CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
+           CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
          endif
+! Compute geopotential (physics might need it)
+         CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
 c$OMP BARRIER
            jD_cur = jD_ref + day_ini - day_ref
@@ -1040 +1050 @@
 c$OMP BARRIER
           if (pressure_exner) then
-            CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+            CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk,pkf)
           else
-            CALL exner_milieu_p(ip1jmp1,ps,p,beta,pks,pk,pkf)
+            CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk,pkf)
           endif
 c$OMP BARRIER
@@ -1202 +1212 @@
 c$OMP BARRIER
         if (pressure_exner) then
-          CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+          CALL exner_hyb_p( ip1jmp1, ps, p, pks, pk, pkf )
         else
-          CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
+          CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
         endif
 c$OMP BARRIER

trunk/LMDZ.COMMON/libf/dyn3dpar/logic.h

@@ -11 +11 @@
      &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
      &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
-     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf,hybrid                    &
+     &  ,ok_limit,ok_etat0,hybrid                                       &
      &  ,moyzon_mu,moyzon_ch
 
@@ -19 +19 @@
      & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
      &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
-     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf
+     &  ,ok_limit,ok_etat0
       logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
                      ! (only used if disvert_type==2)

trunk/LMDZ.COMMON/libf/dyn3dpar/mod_const_mpi.F90

@@ -1 +1 @@
 ! 
-! $Id: mod_const_mpi.F90 1700 2012-12-20 14:43:19Z lguez $
+! $Id: mod_const_mpi.F90 2055 2014-06-04 12:33:27Z acaubel $
 !
 MODULE mod_const_mpi
@@ -17 +17 @@
     USE ioipsl_getincom, only: getin
 #endif
-
+! Use of Oasis-MCT coupler 
+#ifdef CPP_OMCT
+    USE mod_prism
+#endif
+#ifdef CPP_XIOS
+    USE wxios, only: wxios_init
+#endif
     IMPLICIT NONE
 #ifdef CPP_MPI
@@ -38 +44 @@
 #ifdef CPP_COUPLE
 !$OMP MASTER
-       CALL prism_init_comp_proto (comp_id, 'lmdz.x', ierr)
+#ifdef CPP_XIOS
+      CALL wxios_init("LMDZ", outcom=COMM_LMDZ, type_ocean=type_ocean)
+#else
+       CALL prism_init_comp_proto (comp_id, 'LMDZ', ierr)
        CALL prism_get_localcomm_proto(COMM_LMDZ,ierr)
+#endif
 !$OMP END MASTER
 #endif

trunk/LMDZ.COMMON/libf/dyn3dpar/parallel_lmdz.F90

@@ -225 +225 @@
 #endif
 #ifdef CPP_COUPLE
+! Use of Oasis-MCT coupler 
+#if defined CPP_OMCT
+    use mod_prism
+#else
     use mod_prism_proto
 #endif
-#ifdef CPP_EARTH
 ! Ehouarn: surface_data module is in 'phylmd' ...
       use surface_data, only : type_ocean
@@ -252 +255 @@
 
       if (type_ocean == 'couple') then
+#ifdef CPP_XIOS
+    !Fermeture propre de XIOS
+      CALL wxios_close()
+#else
 #ifdef CPP_COUPLE
          call prism_terminate_proto(ierr)
@@ -258 +265 @@
          endif
 #endif 
+#endif
       else
 #ifdef CPP_XIOS