Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/Venus/README
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/Venus/README	(revision 1756)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/Venus/README	(revision 1756)
@@ -0,0 +1,13 @@
+VENUS CALDENDAR : 
+
+WRF TIME COUNTER 1000-02-01_00:00:00
+                 YYYY-MM-DD_HH:mm:ss
+for Venus        0000-99-23_35:12:8
+                 100 MM = 1 Year = 1 Venus day 
+                 24 DD = 1 MM = 1/100 Venus day
+                 36 HH = 1 DD
+                 13 mm = 1 HH
+                 9 ss = 1 mm
+                 
+If GCM has 1 output each n/100 Venus day set interval (in create_readmeteo) at n*864 hours (n*24*36)
+          
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/Venus/create_readmeteo.F90
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/Venus/create_readmeteo.F90	(revision 1755)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/Venus/create_readmeteo.F90	(revision 1756)
@@ -12,8 +12,6 @@
 !------------------------------------------------------------------------!
 
-INTEGER, PARAMETER :: MONTHS_PER_YEAR = 12
-INTEGER, PARAMETER :: mday(MONTHS_PER_YEAR)   &
-           = (/61,66,66,65,60,54,50,46,47,47,51,56/)
-
+INTEGER, PARAMETER :: MONTHS_PER_YEAR = 100
+INTEGER, PARAMETER :: mday = 24 
 INTEGER :: start_day,init,i,month,day
 INTEGER :: ierr,nid,nvarid
@@ -22,5 +20,5 @@
 INTEGER :: no_please
 INTEGER :: timedim,timelen
-
+INTEGER :: inc
 REAL, DIMENSION(100) :: param
 REAL, DIMENSION(:), ALLOCATABLE :: time
@@ -32,5 +30,5 @@
 interval = 0
 
-
+!mday(:)=24
 !
 ! Open input NETCDF file
@@ -38,6 +36,5 @@
 write(*,*) "Scanning netcdf file ..."
 ierr=NF_OPEN ("input_diagfi.nc",NF_NOWRITE,nid)
-ierr=NF__OPEN ("input_diagfi.nc",NF_NOWRITE,4096,nid)
-!ierr=nc_open("input_diagfi.nc",NC_NOWRITE,nid)
+!ierr=NF90_OPEN("input_diagfi.nc",NF_NOWRITE,nid)
 IF (ierr.NE.NF_NOERR) THEN
       write(*,*)'**** Please create a symbolic link called input_diagfi.nc'
@@ -45,5 +42,5 @@
 ENDIF
 
-ierr=NF_INQ_DIMID(nid,"Time",timedim)
+ierr=NF_INQ_DIMID(nid,"time_counter",timedim)
 IF (ierr .NE. NF_NOERR) THEN
     ierr=NF_INQ_DIMID(nid,"time",timedim)
@@ -56,5 +53,5 @@
 !
 ALLOCATE(time(timelen))
-ierr = NF_INQ_VARID (nid, "Time",nvarid)
+ierr = NF_INQ_VARID (nid, "time_counter",nvarid)
 IF (ierr .NE. NF_NOERR) THEN
         ierr = NF_INQ_VARID (nid, "time",nvarid)
@@ -73,5 +70,5 @@
 IF (ierr .NE. NF_NOERR) THEN
         PRINT *, "Error: Readmeteo <ps> not found"
-        stop
+        !stop
 ENDIF
 #ifdef NC_DOUBLE
@@ -92,16 +89,19 @@
 ! beware, param(4) is the day reference of start and startfi
 ! ...have to add time(1) to get the real starting date in diagfi
-start_day=floor(param(4)+time(1))		
-start_hour=nint((param(4)-floor(param(4))+time(1))*24)	! starting hour	
-start_hour=MOD(start_hour,24)   
-IF (interval .eq. 0) interval=nint(time(1)*24)	 ! interval between each time subscript	
-
-
+!start_day=floor(param(4)+time(1))		
+!start_hour=nint((param(4)-floor(param(4))+time(1))*13)	! starting hour	
+!start_hour=MOD(start_hour,24)   
+!IF (interval .eq. 0) interval=nint(time(1)*24)	 ! interval between each time subscript	
+
+start_day=0
+start_hour=0
+interval=864 ! in hours, 864 HH = 1 MM = 1/100 Vd
 PRINT *,'*****************'
 PRINT *,'GCM data file starts at sol ',start_day,'and hour',start_hour
 PRINT *,'GCM data interval is ',interval,'hours'
 
-CALL wrf_day(start_day,month,day)
-
+month=0
+day=0
+!interval=864
 !
 ! User defined parameters
@@ -116,7 +116,7 @@
 if (no_please == 0) stop
 if (no_please == 1) then
-        my=2024
+        my=1000
         n=timelen
-        start=1
+        start=0
         interval_subs=1
 else
@@ -124,5 +124,5 @@
         write(*,*) "Starting Martian year ? ex: 24,25,26..."
         read(*,*) my
-        my=2000+my
+        my=1000+my
         write(*,*) "-- WRF data file information --"
         write(*,*) "How many files do you want to create ? at least 2, max is",timelen
@@ -177,15 +177,16 @@
 END DO
 ! WRF time reference
-hour=start_hour+(start-1)*interval
+hour=interval
 inc_hour=interval*interval_subs
-IF (hour >= 24) day=day+INT(hour/24)
-hour=MOD(hour,24)
-IF (day > mday(month)) THEN
-      day=day-mday(month)
-      month=month+1
+IF (hour >= 36) day=day+INT(hour/36)
+hour=MOD(hour,36)
+IF (day >= mday) THEN
+      inc=INT(day/mday)
+      day=day-inc*mday
+      month=month+inc
 END IF
-IF (month > MONTHS_PER_YEAR) THEN
+IF (month >= MONTHS_PER_YEAR) THEN
       my=my+1
-      month=1
+      month=0
 END IF
 DO i=1,n
@@ -207,13 +208,14 @@
         END IF
         write(6,fmt='(A1)') 'y'
-IF (hour+inc_hour >= 24) day=day+INT((hour+inc_hour)/24)
-hour=MOD(hour+inc_hour,24)
-IF (day > mday(month)) THEN
-        day=day-mday(month)
-        month=month+1
+IF (hour+inc_hour >= 36) day=day+INT((hour+inc_hour)/36)
+hour=MOD(hour+inc_hour,36)
+IF (day >= mday) THEN
+        inc=INT(day/mday)
+        day=day-inc*mday
+        month=month+inc
 END IF        
-IF (month > MONTHS_PER_YEAR) THEN
+IF (month >= MONTHS_PER_YEAR) THEN
         my=my+1
-        month=1
+        month=0
 END IF        
 END DO
@@ -221,44 +223,2 @@
 
 END
-
-
-!--------------------------------------------------
-!--------------------------------------------------
-
-SUBROUTINE wrf_day(gcm_day,wrf_month,day) 
-
-IMPLICIT NONE
-
-INTEGER, INTENT(INOUT) :: gcm_day
-INTEGER, INTENT(OUT) :: wrf_month,day
-
-INTEGER :: init,i
-INTEGER, PARAMETER :: MONTHS_PER_YEAR = 12
-INTEGER, PARAMETER :: mday(MONTHS_PER_YEAR)   &
-           = (/61,66,66,65,60,54,50,46,47,47,51,56/)
-
-!
-! Find WRF month and day
-!
-
-IF (gcm_day >= 669) THEN    !! gcm_day commence au jour 0
-        PRINT *,'out of bounds ! martian year is 669 sols !'
-        gcm_day=MOD(gcm_day,669)
-        !STOP
-ENDIF
-
-
-
-init=gcm_day+1    !!+1 sinon on decale tout
-DO i=1,MONTHS_PER_YEAR
-        wrf_month=i
-        init=init-mday(i)
-        IF (init <= 0) EXIT
-END DO
-
-PRINT *,'corresponding WRF month is ',wrf_month
-day=init+mday(wrf_month)
-PRINT *,'corresponding WRF day is ',day
-PRINT *,'*****************'
-
-END
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/Venus/readmeteo.F90
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/Venus/readmeteo.F90	(revision 1755)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/Venus/readmeteo.F90	(revision 1756)
@@ -41,5 +41,5 @@
 
 REAL :: ptop
-REAL, PARAMETER :: grav=3.72
+REAL, PARAMETER :: grav=8.89
 LOGICAL, PARAMETER :: blank=.false.
 
@@ -73,5 +73,5 @@
 
 !! Intermediate data arrays
-integer :: k,l,m,n,p
+integer :: k,l,m,n,p,i,j
 real, dimension(:), allocatable :: lat,lon,time,alt,aps,bps,levels,vertdsoil
 real, dimension(:,:), allocatable :: vide,ones,ghtsfile
@@ -81,5 +81,6 @@
 real, dimension(:,:,:), allocatable :: emissfile,co2icefile
 real, dimension(:,:,:), allocatable :: tnsfile,unsfile,vnsfile
-real, dimension(:,:,:,:), allocatable :: tfile,ufile,vfile,rfile,hfile
+real, dimension(:,:,:,:), allocatable :: tfile,ufile,vfile,rfile,hfile,ghfile
+real, dimension(:,:,:,:), allocatable :: pfile
 real, dimension(:,:,:,:), allocatable :: eta_gcm
 !real, dimension(:,:,:,:), allocatable :: tfileorig,ufileorig,vfileorig
@@ -173,5 +174,5 @@
 SELECT CASE(ident)
 CASE('LMD')
-ierr=NF_INQ_DIMID(nid,"Time",timedim)
+ierr=NF_INQ_DIMID(nid,"time_counter",timedim)
 CASE('OXF')
 ierr=NF_INQ_DIMID(nid,"time",timedim)
@@ -182,5 +183,5 @@
 SELECT CASE(ident)
 CASE('LMD')
-ierr=NF_INQ_DIMID(nid,"latitude",latdim)
+ierr=NF_INQ_DIMID(nid,"lat",latdim)
 CASE('OXF')
 ierr=NF_INQ_DIMID(nid,"lat",latdim)
@@ -191,5 +192,5 @@
 SELECT CASE(ident)
 CASE('LMD')
-ierr=NF_INQ_DIMID(nid,"longitude",londim)
+ierr=NF_INQ_DIMID(nid,"lon",londim)
 CASE('OXF')
 ierr=NF_INQ_DIMID(nid,"lon",londim)
@@ -200,5 +201,5 @@
 SELECT CASE(ident)
 CASE('LMD')
-ierr=NF_INQ_DIMID(nid,"altitude",altdim)
+ierr=NF_INQ_DIMID(nid,"presnivs",altdim)
 CASE('OXF')
 ierr=NF_INQ_DIMID(nid,"sigma",altdim)
@@ -214,4 +215,5 @@
 allocate(eta_gcm(lonlen,latlen,altlen,timelen))
 allocate(tfile(lonlen,latlen,altlen,timelen))
+allocate(pfile(lonlen,latlen,altlen,timelen))
 allocate(tsoilfile(lonlen,latlen,altlen,timelen))
 allocate(dsoilfile(lonlen,latlen,altlen,timelen))
@@ -244,4 +246,5 @@
 allocate(gwparam(lonlen,latlen,5))
 allocate(ghtsfile(lonlen,latlen))    !! no scan axis
+allocate(ghfile(lonlen,latlen,altlen,timelen))
 allocate(vide(lonlen,latlen))
 allocate(ones(lonlen,latlen))
@@ -270,4 +273,5 @@
 
 tfile(:,:,:,:)=0
+pfile(:,:,:,:)=0
 tsoilfile(:,:,:,:)=0 
 isoilfile(:,:,:,:)=0 
@@ -300,4 +304,5 @@
 gwparam(:,:,:)=0
 ghtsfile(:,:)=0
+ghfile(:,:,:,:)=0
 vide(:,:)=0
 ones(:,:)=0
@@ -313,5 +318,5 @@
 SELECT CASE(ident)
 CASE('LMD')
-   ierr = NF_INQ_VARID (nid, "Time",nvarid)
+   ierr = NF_INQ_VARID (nid, "time_counter",nvarid)
 CASE('OXF')
    ierr = NF_INQ_VARID (nid, "time",nvarid)
@@ -331,5 +336,5 @@
 SELECT CASE(ident)
 CASE('LMD')
-   ierr = NF_INQ_VARID (nid, "latitude",nvarid)
+   ierr = NF_INQ_VARID (nid, "lat",nvarid)
 CASE('OXF')
    ierr = NF_INQ_VARID (nid, "lat",nvarid)
@@ -349,5 +354,5 @@
 SELECT CASE(ident)
 CASE('LMD')
-   ierr = NF_INQ_VARID (nid, "longitude",nvarid)
+   ierr = NF_INQ_VARID (nid, "lon",nvarid)
 CASE('OXF')
    ierr = NF_INQ_VARID (nid, "lon",nvarid)
@@ -364,42 +369,42 @@
    print *,lon(1),' ... to ... ',lon(lonlen),' ... step: ',lon(lonlen)-lon(lonlen-1)
 
-SELECT CASE(ident)
-CASE('LMD')
-print *, 'Hybrid coordinates'
-   ierr = NF_INQ_VARID (nid, "aps", nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <aps> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, aps)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, aps)
-#endif
-   ierr = NF_INQ_VARID (nid, "bps", nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <bps> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bps)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, bps)
-#endif
-   print *,aps(1),' ... to ... ',aps(altlen)
-   print *,bps(1),' ... to ... ',bps(altlen)
-CASE('OXF')
-   ierr = NF_INQ_VARID (nid, "sigma", nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <sigma> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bps)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, bps)
-#endif
-   print *,bps(1),' ... to ... ',bps(altlen)
-END SELECT
+!SELECT CASE(ident)
+!CASE('LMD')
+!print *, 'Hybrid coordinates'
+!   ierr = NF_INQ_VARID (nid, "aps", nvarid)
+!   IF (ierr .NE. NF_NOERR) THEN
+!      PRINT *, "Error: Readmeteo <aps> not found"
+!      stop
+!   ENDIF
+!#ifdef NC_DOUBLE
+!   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, aps)
+!#else
+!   ierr = NF_GET_VAR_REAL(nid, nvarid, aps)
+!#endif
+!   ierr = NF_INQ_VARID (nid, "bps", nvarid)
+!   IF (ierr .NE. NF_NOERR) THEN
+!      PRINT *, "Error: Readmeteo <bps> not found"
+!      stop
+!   ENDIF
+!#ifdef NC_DOUBLE
+!   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bps)
+!#else
+!   ierr = NF_GET_VAR_REAL(nid, nvarid, bps)
+!#endif
+!   print *,aps(1),' ... to ... ',aps(altlen)
+!   print *,bps(1),' ... to ... ',bps(altlen)
+!CASE('OXF')
+!   ierr = NF_INQ_VARID (nid, "sigma", nvarid)
+!  IF (ierr .NE. NF_NOERR) THEN
+!      PRINT *, "Error: Readmeteo <sigma> not found"
+!      stop
+!   ENDIF
+!#ifdef NC_DOUBLE
+!   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bps)
+!#else
+!   ierr = NF_GET_VAR_REAL(nid, nvarid, bps)
+!#endif
+!   print *,bps(1),' ... to ... ',bps(altlen)
+!END SELECT
 
 
@@ -443,5 +448,5 @@
 
 print *,'Surface Pressure'
-   ierr = NF_INQ_VARID (nid,"ps",nvarid)
+   ierr = NF_INQ_VARID (nid,"psol",nvarid)
    IF (ierr .NE. NF_NOERR) THEN
       PRINT *, "Error: Readmeteo <ps> not found"
@@ -455,5 +460,5 @@
 
 print *,'Ground Temperature'
-   ierr = NF_INQ_VARID (nid,"tsurf",nvarid)
+   ierr = NF_INQ_VARID (nid,"tsol",nvarid)
    IF (ierr .NE. NF_NOERR) THEN
       PRINT *, "Error: Readmeteo <tsurf> not found"
@@ -518,40 +523,33 @@
 !   vnsfile=vfileorig(:,:,1,:)
 
-SELECT CASE(ident)
-CASE('LMD')
-   print *,'Geopotential height at the ground'
-   ierr = NF_INQ_VARID (nid,"phisinit",nvarid)
+print *,'>>> Read 3D meteorological fields ! - This may take some time ...'
+
+print *,'Geopotential height'
+   ierr = NF_INQ_VARID (nid,"geop",nvarid)
    IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <phisinit> not found"
+      PRINT *, "Error: Readmeteo <geop> not found"
       stop
    ENDIF
 #ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ghtsfile)
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ghfile)
 #else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, ghtsfile)
+   ierr = NF_GET_VAR_REAL(nid, nvarid, ghfile)
 #endif
-!**** from geopotential to geopotential height
-ghtsfile=ghtsfile/grav
-!****
-CASE('OXF')
-!
-! geopotential height ~ altimetry
-!
-print *,'Geopotential height at the ground from file mountain_new.nc'
-ierr=NF_OPEN("mountain_new.nc",NF_NOWRITE,nid3)
-    if (ierr.ne.NF_NOERR) then
-      write(*,*) "Error: Could not open that file either"
-      stop "Might as well stop here"
-    endif
-ierr=NF_INQ_VARID(nid3,"orography",nvarid)
-ierr=NF_GET_VAR_REAL(nid3,nvarid,ghtsfile)
-if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed reading phisinit"
-endif
-ierr=NF_CLOSE(nid3)
-END SELECT
-
-
-print *,'>>> Read 3D meteorological fields ! - This may take some time ...'
+ghtsfile=ghfile(:,:,1,1)/grav !surface geop
+
+print *,'Pressure'
+   ierr = NF_INQ_VARID (nid,"pres",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      ierr = NF_INQ_VARID (nid,"p",nvarid)
+        IF (ierr .NE. NF_NOERR) THEN
+          PRINT *, "Error: Readmeteo <p> not found"
+          stop
+        ENDIF
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, pfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, pfile)
+#endif
 
 print *,'Temperature'
@@ -572,5 +570,5 @@
 
 print *,'Zonal wind'   
-   ierr = NF_INQ_VARID (nid,"u",nvarid)
+   ierr = NF_INQ_VARID (nid,"vitu",nvarid)
    IF (ierr .NE. NF_NOERR) THEN
       PRINT *, "Error: Readmeteo <u> not found"
@@ -585,5 +583,5 @@
 
 print *,'Meridional wind'
-   ierr = NF_INQ_VARID (nid,"v",nvarid)
+   ierr = NF_INQ_VARID (nid,"vitv",nvarid)
    IF (ierr .NE. NF_NOERR) THEN
       PRINT *, "Error: Readmeteo <v> not found"
@@ -595,6 +593,5 @@
    ierr = NF_GET_VAR_REAL(nid, nvarid, vfile)
 #endif
-vnsfile=ufile(:,:,1,:)
-
+vnsfile=vfile(:,:,1,:)
 
 !!------------------------
@@ -762,4 +759,19 @@
 
 !END SELECT
+
+print*,'VENUS : rotation backward -> inversion of lat and long'
+DO i=1,latlen
+DO j=1,lonlen
+  psfile(j,i,:)=psfile(lonlen+1-j,latlen+1-i,:)
+  tsfile(i,j,:)=tsfile(lonlen+1-j,latlen+1-i,:)
+  ghtsfile(j,i)=ghtsfile(lonlen+1-j,latlen+1-i)
+  ghfile(j,i,:,:)=ghfile(lonlen+1-j,latlen+1-i,:,:)
+  pfile(j,i,:,:)=pfile(lonlen+1-j,latlen+1-i,:,:)
+  tfile(j,i,:,:)=tfile(lonlen+1-j,latlen+1-i,:,:)
+  ufile(j,i,:,:)=ufile(lonlen+1-j,latlen+1-i,:,:)
+  vfile(j,i,:,:)=vfile(lonlen+1-j,latlen+1-i,:,:)
+  tsoilfile(j,i,:,:)=tsoilfile(lonlen+1-j,latlen+1-i,:,:)
+ENDDO
+ENDDO  
 
 ierr=NF_CLOSE(nid)
@@ -876,5 +888,4 @@
 	write(1) SLAB
 !print *,'The field '//DESC//' was written to '//output
-
 !------------------------!    
 ! >>> Write a variable   !
@@ -894,5 +905,4 @@
 	write(1) SLAB
 !print *,'The field '//DESC//' was written to '//output
-
 !------------------------!    
 ! >>> Write a variable   !
@@ -903,5 +913,5 @@
 DESC='Meridional wind'
 XLVL=200100.
-SLAB=vnsfile(:,:,time_out(l))
+SLAB=-1*vnsfile(:,:,time_out(l)) !VENUS ratoting backwards : v=-1*v
 	! And now put everything in the destination file
 	! ... Header
@@ -912,5 +922,4 @@
 	write(1) SLAB
 !print *,'The field '//DESC//' was written to '//output
-
 !------------------------!    
 ! >>> Write a variable   !
@@ -1366,5 +1375,5 @@
 DO k = 1,altlen
 	XLVL=levels(k)
-	SLAB=tfile(:,:,k,time_out(l))
+        SLAB=tfile(:,:,k,time_out(l))
 		! And now put everything in the destination file
 		! ... Header
@@ -1386,5 +1395,5 @@
 DO k = 1,altlen
         XLVL=levels(k)
-	SLAB=ufile(:,:,k,time_out(l))
+        SLAB=ufile(:,:,k,time_out(l))
 		! And now put everything in the destination file
 		! ... Header
@@ -1394,4 +1403,5 @@
 		! ... Data
 		write(1) SLAB
+                !write(1) ufile(:,:,k,time_out(l))
 END DO
 !print *,'The field '//DESC//' was written to '//output
@@ -1406,5 +1416,5 @@
 DO k = 1,altlen
         XLVL=levels(k)
-	SLAB=vfile(:,:,k,time_out(l))
+        SLAB=-1*vfile(:,:,k,time_out(l)) ! VENUS rotation backward : v=-1*v
 		! And now put everything in the destination file
 		! ... Header
@@ -1430,7 +1440,7 @@
 SELECT CASE(ident)
 CASE('LMD')
-        SLAB=aps(k)+bps(k)*psfile(:,:,time_out(l))
+        SLAB=pfile(:,:,k,time_out(l))
 CASE('OXF')
-	SLAB=bps(k)*psfile(:,:,time_out(l))
+	!SLAB=bps(k)*psfile(:,:,time_out(l))
 END SELECT
 		! And now put everything in the destination file
@@ -1453,13 +1463,9 @@
 DO k = 1,altlen
         XLVL=levels(k)
-!!*******
-!! PROVISOIRE: normalement, il faudrait la hauteur geopotentielle
-!!*******
-!!however, not used by initialize_real
 SELECT CASE(ident)
 CASE('LMD')
-        SLAB=10000.*alog(610./(aps(k)+bps(k)*psfile(:,:,time_out(l))))
+        SLAB=(ghfile(:,:,k,time_out(l)))/grav
 CASE('OXF')
-	SLAB=10000.*alog(610./(bps(k)*psfile(:,:,time_out(l))))
+	!SLAB=10000.*alog(610./(bps(k)*psfile(:,:,time_out(l))))
 END SELECT
 		! And now put everything in the destination file
@@ -1563,4 +1569,5 @@
         XLVL=levels(k)
         SLAB=isoilfile(:,:,k,time_out(l))
+        
                 ! And now put everything in the destination file
                 ! ... Header
@@ -1710,4 +1717,5 @@
 deallocate(eta_gcm)
 deallocate(tfile)
+deallocate(pfile)
 deallocate(tsoilfile)
 deallocate(isoilfile)
@@ -1737,4 +1745,5 @@
 deallocate(co2icefile)
 deallocate(ghtsfile)	!! no scan axis
+deallocate(ghfile)
 deallocate(vide)
 deallocate(ones)
