Index: /codes/icosagcm/trunk/src/disvert.f90
===================================================================
--- /codes/icosagcm/trunk/src/disvert.f90	(revision 466)
+++ /codes/icosagcm/trunk/src/disvert.f90	(revision 467)
@@ -23,4 +23,5 @@
   SUBROUTINE init_disvert
   USE disvert_std_mod, ONLY: ap_std=>ap, bp_std=>bp, presnivs_std=>presnivs, init_disvert_std=>init_disvert
+  USE disvert_strato_mod, ONLY: ap_strato=>ap, bp_strato=>bp, presnivs_strato=>presnivs, init_disvert_strato=>init_disvert
   USE disvert_apbp_mod, ONLY: ap_apbp=>ap, bp_apbp=>bp, presnivs_apbp=>presnivs, init_disvert_apbp=>init_disvert
   USE disvert_ncar_mod, ONLY: ap_ncar=>ap, bp_ncar=>bp, presnivs_ncar=>presnivs, init_disvert_ncar=>init_disvert
@@ -65,4 +66,11 @@
         presnivs=>presnivs_std
       
+      CASE('strato')
+        
+        CALL init_disvert_strato
+        ap=>ap_strato
+        bp=>bp_strato
+        presnivs=>presnivs_strato
+
       CASE('read_apbp')
     
Index: /codes/icosagcm/trunk/src/disvert_strato.f90
===================================================================
--- /codes/icosagcm/trunk/src/disvert_strato.f90	(revision 467)
+++ /codes/icosagcm/trunk/src/disvert_strato.f90	(revision 467)
@@ -0,0 +1,84 @@
+MODULE disvert_strato_mod
+  USE icosa
+  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: ap(:)
+!$OMP THREADPRIVATE(ap)
+  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: bp(:)
+!$OMP THREADPRIVATE(bp)
+  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:)
+!$OMP THREADPRIVATE(presnivs)
+
+CONTAINS
+
+  SUBROUTINE init_disvert
+  USE icosa
+  IMPLICIT NONE
+  
+    ALLOCATE(ap(llm+1))
+    ALLOCATE(bp(llm+1))
+    ALLOCATE(presnivs(llm))
+
+    CALL disvert(ap,bp,presnivs)    
+
+  END SUBROUTINE init_disvert  
+
+
+  SUBROUTINE disvert(ap,bp,presnivs)
+  USE icosa
+  USE mpipara
+  USE earth_const
+  IMPLICIT NONE
+  REAL(rstd),INTENT(OUT) :: ap(:)
+  REAL(rstd),INTENT(OUT) :: bp(:)
+  REAL(rstd),INTENT(OUT) :: presnivs(:)
+  
+  REAL(rstd) :: dsig(llm)
+  REAL(rstd) :: sig(llm+1)
+  REAL(rstd) :: snorm
+  REAL(rstd) :: x
+  INTEGER :: l
+  
+    snorm  = 0.
+    DO l = 1, llm
+      x = Pi*(l-0.5)/(llm+1)
+      dsig(l) = (1.0 + 7.0 * SIN(x)**2) * (0.5* (1.-tanh( (x-Pi/2)/(Pi/2) ) ))**2
+      snorm = snorm + dsig(l)
+    ENDDO    
+    
+    DO l = 1, llm
+      dsig(l) = dsig(l)/snorm
+    ENDDO
+
+    sig(llm+1) = 0.
+    DO l = llm, 1, -1
+      sig(l) = sig(l+1) + dsig(l)
+    ENDDO
+
+    bp(llm+1) =   0.
+    DO l = 1, llm
+      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
+      ap(l) = pa * ( sig(l) - bp(l) )
+    ENDDO
+    bp(1)=1.
+    ap(1)=0.
+    ap(llm+1) = pa * ( sig(llm+1) - bp(llm+1) )
+    DO l = 1, llm
+      presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
+    ENDDO
+    
+    ! tell the world about it
+    IF (is_mpi_root) THEN
+!$OMP MASTER
+      WRITE(*,*) "ap()=",ap
+      WRITE(*,*) "bp()=",bp
+      WRITE(*,*) "Approximative mid-layer pressure, assuming a surface pressure preff=",preff," Pa"
+      WRITE(*,*) "and approximative mid-layer height, assuming an atmospheric scale height of ",scale_height/1000," (km)"
+      DO l=1,llm
+        WRITE(*,*) 'PRESNIVS(',l,')=',presnivs(l),'  Z ~ ',log(preff/presnivs(l))*scale_height/1000,       &
+                   ' DZ ~ ',scale_height/1000*log((ap(l)+bp(l)*preff)/ max(ap(l+1)+bp(l+1)*preff,1.e-10))
+      ENDDO
+!$OMP END MASTER
+    ENDIF
+  
+  END SUBROUTINE disvert
+  
+END MODULE disvert_strato_mod
Index: /codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- /codes/icosagcm/trunk/src/etat0.f90	(revision 466)
+++ /codes/icosagcm/trunk/src/etat0.f90	(revision 467)
@@ -9,8 +9,38 @@
     REAL(rstd) :: etat0_temp
 
-    PUBLIC :: etat0, etat0_type
+    PUBLIC :: etat0, init_etat0, etat0_type
 
 CONTAINS
   
+  SUBROUTINE init_etat0
+  USE etat0_database_mod
+  IMPLICIT NONE
+
+    CALL getin("etat0",etat0_type)
+
+    SELECT CASE (TRIM(etat0_type))
+      CASE ('isothermal')
+      CASE ('temperature_profile')
+      CASE ('jablonowsky06')
+      CASE ('dcmip5')
+      CASE ('williamson91.6')
+      CASE ('start_file')
+      CASE ('database')
+        CALL init_etat0_database
+      CASE ('academic')
+      CASE ('held_suarez')
+      CASE ('venus')
+      CASE ('dcmip1')
+      CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
+      CASE ('dcmip3')
+      CASE ('dcmip4')
+      CASE DEFAULT
+         PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
+            '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> '
+         STOP
+    END SELECT
+
+  END SUBROUTINE init_etat0
+
   SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_w, f_q)
     USE mpipara, ONLY : is_mpi_root
@@ -28,7 +58,8 @@
     USE etat0_dcmip2016_supercell_mod, ONLY : getin_etat0_dcmip2016_supercell=>getin_etat0
     ! Ad hoc interfaces
-    USE etat0_academic_mod, ONLY : etat0_academic=>etat0  
-    USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0  
-    USE etat0_venus_mod,  ONLY : etat0_venus=>etat0  
+    USE etat0_academic_mod, ONLY : etat0_academic=>etat0
+    USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0
+    USE etat0_venus_mod,  ONLY : etat0_venus=>etat0
+    USE etat0_database_mod, ONLY : etat0_database=>etat0
     USE etat0_start_file_mod, ONLY : etat0_start_file=>etat0  
 
@@ -93,4 +124,6 @@
     !------------------- Ad hoc interfaces --------------------
     SELECT CASE (TRIM(etat0_type))
+     CASE ('database')
+        CALL etat0_database(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
     CASE ('start_file')
        CALL etat0_start_file(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
@@ -109,5 +142,5 @@
          PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
               '> options are <isothermal>, <temperature_profile>, <held_suarez>, &
-              &<bubble>, <venus>, <jablonowsky06>, <academic>, <dcmip[1-4]> '
+              &<bubble>, <venus>, <jablonowsky06>, <start_file>, <database>, <academic>, <dcmip[1-4]> '
          STOP
       END IF
Index: /codes/icosagcm/trunk/src/etat0_database.f90
===================================================================
--- /codes/icosagcm/trunk/src/etat0_database.f90	(revision 467)
+++ /codes/icosagcm/trunk/src/etat0_database.f90	(revision 467)
@@ -0,0 +1,152 @@
+MODULE etat0_database_mod
+
+
+CONTAINS
+
+  SUBROUTINE init_etat0_database
+  USE xios
+  IMPLICIT NONE
+  
+    CALL xios_set_fieldgroup_attr("read_fields",enabled=.TRUE.)
+    CALL xios_set_filegroup_attr("read_files",enabled=.TRUE.)
+
+  END SUBROUTINE init_etat0_database
+
+  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+  USE icosa
+  USE restart_mod
+  USE wind_mod
+  USE write_field_mod
+  USE time_mod
+  USE transfert_mod
+  USE xios_mod
+  USE write_field_mod
+  USE vertical_remap_mod
+  USE theta2theta_rhodz_mod
+  USE qsat_mod
+  USE pression_mod
+  USE omp_para
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_phis(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_q(:)
+
+    TYPE(t_field),POINTER,SAVE :: f_ulon_reg(:)
+    TYPE(t_field),POINTER,SAVE :: f_ulat_reg(:)
+    TYPE(t_field),POINTER,SAVE :: f_temp_reg(:)
+    TYPE(t_field),POINTER,SAVE :: f_q_reg(:)
+
+    TYPE(t_field),POINTER,SAVE :: f_ts(:)
+    TYPE(t_field),POINTER,SAVE :: f_z(:)
+    TYPE(t_field),POINTER,SAVE :: f_ulon(:)
+    TYPE(t_field),POINTER,SAVE :: f_ulat(:)
+    TYPE(t_field),POINTER,SAVE :: f_temp(:)
+    TYPE(t_field),POINTER,SAVE :: f_q1(:)
+    TYPE(t_field),POINTER,SAVE :: f_qsat(:)
+    TYPE(t_field),POINTER,SAVE :: f_p(:)
+    INTEGER :: nb_level
+    REAL,ALLOCATABLE:: levels(:)
+    INTEGER :: ind
+
+    CALL xios_read_field("relief",f_phis)
+    
+    CALL writeField("relief_out",f_phis,once=.TRUE.)
+
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      f_phis(ind)%rval2d(:)=f_phis(ind)%rval2d(:)*g      
+    ENDDO
+
+
+    IF (is_omp_master) CALL xios_get_axis_attr("lev_ecdyn",n_glo=nb_level)
+    CALL bcast_omp(nb_level)
+    ALLOCATE(levels(nb_level))
+
+    IF (is_omp_master) CALL xios_get_axis_attr("lev_ecdyn",value=levels)
+    CALL bcast_omp(levels)
+    
+    levels=levels*100  ! hectoPascal -> Pascal
+ 
+    CALL allocate_field(f_ts, field_t, type_real, name="ts")
+    CALL allocate_field(f_z, field_t, type_real, name="z")
+    CALL allocate_field(f_ulon_reg, field_t, type_real,nb_level)
+    CALL allocate_field(f_ulat_reg, field_t, type_real,nb_level)
+    CALL allocate_field(f_temp_reg, field_t, type_real,nb_level)
+    CALL allocate_field(f_q_reg,    field_t, type_real,nb_level)
+    
+    CALL allocate_field(f_q1,    field_t, type_real,llm)
+    CALL allocate_field(f_qsat,  field_t, type_real,llm)
+    CALL allocate_field(f_p,  field_t, type_real,llm+1)
+    CALL allocate_field(f_temp,  field_t, type_real,llm)
+    CALL allocate_field(f_ulon,  field_t, type_real,llm)
+    CALL allocate_field(f_ulat,  field_t, type_real,llm)
+
+    CALL xios_read_field("z",f_z)
+    CALL xios_read_field("ps",f_ps)
+    CALL xios_read_field("ts",f_ts)
+    CALL writeField("ps_out",f_ps)
+
+!$OMP BARRIER
+    
+!    CALL writeField("phis_out",f_phis,once=.TRUE.)
+!    CALL writeField("ts_out",f_ts,once=.TRUE.)
+
+! make correction to ps due to relief at higher resolution
+! difference with LMDZ : tsol is taken from ECDYN.NC and not from ECPHY
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      f_ps(ind)%rval2d(:)=f_ps(ind)%rval2d(:)*(1.+(f_z(ind)%rval2d(:)-f_phis(ind)%rval2d(:))/287.0/f_ts(ind)%rval2d(:))
+    ENDDO
+    CALL transfert_request(f_ps,req_i0) 
+    CALL writeField("ps_out",f_ps)
+    
+    
+
+    CALL xios_read_field("temp",f_temp_reg)
+    CALL vertical_remap(levels,f_temp_reg,f_ps,f_temp)
+    CALL transfert_request(f_temp,req_i0) 
+    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
+
+    CALL xios_read_field("u",f_ulon_reg)
+    CALL vertical_remap(levels,f_ulon_reg,f_ps,f_ulon)
+    CALL xios_read_field("v",f_ulat_reg)
+    CALL vertical_remap(levels,f_ulat_reg,f_ps,f_ulat)
+    CALL transfert_request(f_ulat,req_i0) 
+    CALL transfert_request(f_ulon,req_i0) 
+    CALL ulonlat2un(f_ulon, f_ulat,f_u)
+
+    CALL xios_read_field("q",f_q_reg)
+    CALL vertical_remap(levels,f_q_reg,f_ps,f_q1)
+
+    CALL pression(f_ps,f_p)
+! difference with LMDZ : for qsat, pressure at mid layer is computed as a mean value pmid=(p(l)+p(l+1))/2    
+    CALL qsat(f_temp,f_p,f_qsat)
+    CALL transfert_request(f_qsat,req_i0) 
+
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      f_q(ind)%rval4d(:,:,:)=0
+      f_q(ind)%rval4d(:,:,1)=f_q1(ind)%rval3d(:,:)*f_qsat(ind)%rval3d(:,:)*0.01
+      WHERE(f_q(ind)%rval4d(:,:,1)<0) f_q(ind)%rval4d(:,:,1)=0
+    ENDDO
+
+
+    CALL deallocate_field(f_ts)
+    CALL deallocate_field(f_z)
+    CALL deallocate_field(f_ulon_reg)
+    CALL deallocate_field(f_ulat_reg)
+    CALL deallocate_field(f_temp_reg)
+    CALL deallocate_field(f_q_reg)
+    
+    CALL deallocate_field(f_q1)
+    CALL deallocate_field(f_qsat)
+    CALL deallocate_field(f_p)
+    CALL deallocate_field(f_temp)
+    CALL deallocate_field(f_ulon)
+    CALL deallocate_field(f_ulat)   
+ 
+  END SUBROUTINE etat0
+
+END MODULE etat0_database_mod
Index: /codes/icosagcm/trunk/src/icosagcm.f90
===================================================================
--- /codes/icosagcm/trunk/src/icosagcm.f90	(revision 466)
+++ /codes/icosagcm/trunk/src/icosagcm.f90	(revision 467)
@@ -31,4 +31,5 @@
     CALL compute_domain
     CALL init_transfert
+    CALL init_etat0
     CALL init_writefield
     CALL init_trace
Index: /codes/icosagcm/trunk/src/q_sat.f90
===================================================================
--- /codes/icosagcm/trunk/src/q_sat.f90	(revision 467)
+++ /codes/icosagcm/trunk/src/q_sat.f90	(revision 467)
@@ -0,0 +1,98 @@
+MODULE qsat_mod
+
+
+
+
+CONTAINS
+
+  SUBROUTINE qsat(f_temp,f_p,f_qs)
+  USE icosa
+  IMPLICIT NONE
+    TYPE(t_field), POINTER :: f_temp(:)  ! IN    : temperature
+    TYPE(t_field), POINTER :: f_p(:)     ! IN    : pressure at mid-levels
+    TYPE(t_field), POINTER :: f_qs(:)    ! OUT   : vapeur d'eau saturante en kg/kg
+    
+    REAL(rstd),POINTER :: temp(:,:),  p(:,:), qs(:,:)
+    INTEGER :: ind
+
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       temp=f_temp(ind)
+       p=f_p(ind)
+       qs=f_qs(ind)
+       CALL compute_qsat(temp,p,qs)
+    END DO
+  
+  END SUBROUTINE qsat
+
+  
+  SUBROUTINE compute_qsat(temp,p,qsat)
+  USE icosa
+  USE omp_para
+  IMPLICIT NONE
+  
+!======================================================================
+! 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
+!======================================================================
+!
+      REAL,INTENT(IN)  :: temp(iim*jjm,llm)
+      REAL,INTENT(IN)  :: p   (iim*jjm,llm+1)
+      REAL,INTENT(OUT) :: qsat(iim*jjm,llm)
+
+      REAL, PARAMETER  ::  r2es=611.14 *18.0153/28.9644
+
+      REAL :: r3es
+      REAL, PARAMETER ::  r3les=17.269
+      REAL, PARAMETER ::  r3ies=21.875
+
+      REAL :: r4es
+      REAL, PARAMETER :: r4les=35.86
+      REAL, PARAMETER :: r4ies=7.66
+
+      REAL, PARAMETER :: rtt=273.16
+      REAL, PARAMETER :: retv=28.9644/18.0153 - 1.0
+
+      REAL :: zqsat, pmid
+      INTEGER :: l,i,j,ij
+!
+!     ------------------------------------------------------------------
+!
+!
+
+    DO l=ll_begin,ll_end
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+
+          IF (temp(ij,l) .LE. rtt) THEN
+            r3es = r3ies
+            r4es = r4ies
+          ELSE
+            r3es = r3les
+            r4es = r4les
+          ENDIF
+          pmid=0.5*(p(ij,l)+p(ij,l+1))
+          zqsat=r2es/pmid*EXP(r3es*(temp(ij,l)-rtt)/(temp(ij,l)-r4es))
+          zqsat=MIN(0.5,zqsat)
+          zqsat=zqsat/(1.-retv *zqsat)
+
+          qsat(ij,l)= zqsat
+
+        ENDDO
+      ENDDO
+    ENDDO
+  
+  
+  END SUBROUTINE compute_qsat
+
+END MODULE qsat_mod
Index: /codes/icosagcm/trunk/src/vertical_remap.f90
===================================================================
--- /codes/icosagcm/trunk/src/vertical_remap.f90	(revision 467)
+++ /codes/icosagcm/trunk/src/vertical_remap.f90	(revision 467)
@@ -0,0 +1,110 @@
+MODULE vertical_remap_mod
+  USE icosa
+  
+
+CONTAINS
+
+   
+  SUBROUTINE vertical_remap(pressure_level,field_in,f_ps,field_out)
+  USE icosa
+  USE pression_mod
+  USE omp_para
+  IMPLICIT NONE
+    REAL(rstd), INTENT(IN) :: pressure_level(:)
+    TYPE(t_field),POINTER :: field_in(:)
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: field_out(:)
+
+    TYPE(t_field),POINTER,SAVE :: f_p(:)
+    REAL(rstd),POINTER :: in(:,:)
+    REAL(rstd),POINTER :: out(:,:)
+    REAL(rstd),POINTER :: p(:,:)
+    
+    INTEGER :: ind
+
+    CALL allocate_field(f_p,field_t,type_real,llm+1)
+    CALL pression(f_ps,f_p)
+ 
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      p=f_p(ind)
+      in=field_in(ind)
+      out=field_out(ind)
+      CALL compute_vertical_remap(pressure_level,in,p,out)
+    ENDDO
+    
+  END SUBROUTINE  vertical_remap
+
+  SUBROUTINE compute_vertical_remap(pressure_level,in,p,out)
+  USE omp_para
+  IMPLICIT NONE
+    REAL(rstd),INTENT(IN)  :: pressure_level(:)
+    REAL(rstd),INTENT(IN)  :: in(:,:)
+    REAL(rstd),INTENT(IN)  :: p(iim*jjm,llm+1)
+    REAL(rstd),INTENT(OUT) :: out(iim*jjm,llm)
+    REAL(rstd) :: coeff, pmid
+    INTEGER :: i,j,ij,l,n,nb_level
+    INTEGER :: a
+    INTEGER :: b
+    LOGICAL :: positive
+    
+    nb_level=size(pressure_level)
+    IF (pressure_level(1)>=pressure_level(nb_level)) THEN 
+      positive=.FALSE.
+    ELSE 
+      positive=.TRUE.
+    ENDIF
+      
+ !$OMP BARRIER    
+    IF (is_omp_level_master) THEN
+    
+    DO l=1,llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          pmid=0.5*(p(ij,l)+p(ij,l+1))
+          IF (positive) THEN
+            a=0
+            DO n=1,nb_level-1
+              IF ( (pmid>=pressure_level(n) .AND. pmid<pressure_level(n+1))) THEN
+               a=n ; b=n+1 ; EXIT
+              ENDIF
+            ENDDO
+            IF (a==0) THEN
+              IF (pmid<=pressure_level(1)) THEN
+                a=1 ; b=2
+              ELSE
+                a=nb_level-1 ; b=nb_level
+              ENDIF
+            ENDIF
+          ELSE
+            a=0
+            DO n=1,nb_level-1
+              IF ( (pmid<=pressure_level(n) .AND. pmid>pressure_level(n+1))) THEN
+               a=n ; b=n+1 ; EXIT
+              ENDIF
+            ENDDO
+            
+            IF (a==0) THEN
+              IF (pmid >= pressure_level(1)) THEN
+                a=1 ; b=2
+              ELSE
+                a=nb_level-1 ; b=nb_level
+              ENDIF
+            ENDIF
+          ENDIF
+                  
+          coeff=(pmid-pressure_level(a))/(pressure_level(a)-pressure_level(b))
+          out(ij,l)=in(ij,a)+coeff*(in(ij,a)-in(ij,b))
+        ENDDO
+      ENDDO
+    ENDDO
+ 
+    ENDIF
+ !$OMP BARRIER    
+        
+  END SUBROUTINE compute_vertical_remap
+
+END MODULE vertical_remap_mod
