Index: trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/obsolete/readmeteo.F90_backup
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/obsolete/readmeteo.F90_backup	(revision 72)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/obsolete/readmeteo.F90_backup	(revision 72)
@@ -0,0 +1,974 @@
+program readmeteo
+
+implicit none
+include "netcdf.inc"
+
+!------------------------------------------------------------------------!
+! Readmeteo prepares an initial state for the WPS pre-processor of WRF   !
+!                                                                        !
+! Input is a diagfi.nc NETCDF file from a LMD/GCM simulation             !
+!                                                                        !
+! Outputs are WRFSI intermediate format files ready for metgrid.exe      !
+!       http://www.mmm.ucar.edu/wrf/OnLineTutorial/WPS/IM_si.htm         !
+!                                                                        !
+! A. Spiga - 16/04/2007                                                  !
+!	     22/05/2007 : need to run zrecast as a preliminary           !
+!            07/06/2007 : external parameter file 'readmeteo.def         !
+!	     30/07/2007 : manage additional variables (tsoil, emiss,...) !	
+!            19/10/2007 : no need to run zrecast at all (eta levels)     !   
+!		12/2007 : include co2ice and emissivity			 ! 
+!									 !
+! spiga@lmd.jussieu.fr							 !	
+!------------------------------------------------------------------------!
+
+
+
+!***************************************************************************
+!***************************************************************************
+REAL, PARAMETER :: emiss_prescribed=0.95
+REAL, PARAMETER :: co2ice_prescribed=0.
+
+
+
+
+!***************************************************************************
+!***************************************************************************
+
+REAL :: ptop
+REAL, PARAMETER :: grav=3.72
+LOGICAL, PARAMETER :: blank=.false.
+CHARACTER (LEN=3), PARAMETER :: ident='LMD'
+
+
+!! Variables to be written in the output file
+!! *** NB: conformity with metgrid/src/read_met_module.F90
+CHARACTER (LEN=33) :: output
+INTEGER :: IFV
+CHARACTER (LEN=24) :: HDATE
+REAL :: XFCST
+CHARACTER (LEN=32) :: SOURCE
+CHARACTER (LEN=9) :: FIELD
+CHARACTER (LEN=25) :: UNITS
+CHARACTER (LEN=46) :: DESC
+REAL :: XLVL
+INTEGER :: NX,NY,IPROJ
+CHARACTER (LEN=8) :: STARTLOC
+REAL :: STARTLAT,STARTLON,DELTALAT,DELTALON
+REAL, POINTER, DIMENSION(:,:) :: SLAB
+
+!! Variables related to NETCDF format
+integer :: nid,nid2,nvarid,ierr,ierr2
+integer :: timedim,londim,latdim,altdim,altdim2
+integer :: rlatvdim,rlonudim
+integer :: timelen,lonlen,latlen,altlen,altlen2
+integer :: rlatvlen,rlonulen
+integer :: soillen,soildim
+integer, dimension(4) :: corner,edges
+
+!! Intermediate data arrays
+integer :: k,l
+real, dimension(:), allocatable :: lat,lon,time,alt,aps,bps,levels
+real, dimension(:,:), allocatable :: vide,ones,ghtsfile
+real, dimension(:,:,:), allocatable :: psfile,tsfile
+real, dimension(:,:,:), allocatable :: emissfile,co2icefile
+real, dimension(:,:,:), allocatable :: tnsfile,unsfile,vnsfile
+real, dimension(:,:,:,:), allocatable :: tfile,ufile,vfile,rfile,hfile
+real, dimension(:,:,:,:), allocatable :: eta_gcm
+real, dimension(:,:,:,:), allocatable :: tfileorig,ufileorig,vfileorig
+real, dimension(:,:,:,:), allocatable :: tsoilfile
+
+!! Reading the parameter file 
+integer :: tmp,increment,FILES
+integer :: tmp2,tmp3,tmp4
+character*1 :: answer
+character*4 :: str
+character*2 :: str2, str3, str4
+integer, dimension(:), allocatable :: time_out
+character*13, dimension(:), allocatable :: date_out
+character*19, dimension(:), allocatable :: date_out2
+!***************************************************************************
+!***************************************************************************
+
+
+!!---------------------
+!! Open the datafile
+!!---------------------
+ierr=NF_OPEN ("input_diagfi.nc",NF_NOWRITE,nid)
+IF (ierr.NE.NF_NOERR) THEN
+   write(*,*)'**** Please create a symbolic link called input_diagfi.nc'
+   CALL ABORT
+ENDIF
+
+
+!!--------------------------
+!! Ask for data references
+!!--------------------------
+write(*,*) "Create n files. What is n ?"
+read(*,*) FILES
+allocate(time_out(FILES))
+allocate(date_out(FILES))
+allocate(date_out2(FILES))
+
+write(*,*) ""
+write(*,*) "INPUT: Diagfi time"
+write(*,*) "list of n subscripts, separated by <Return>s"
+increment=1
+do while (increment.ne.FILES+1)
+    read(*,*) tmp
+    time_out(increment)=tmp
+    increment=increment+1
+enddo
+
+write(*,*) ""
+write(*,*) "OUTPUT: WRF time"
+write(*,*) "fill 4 lines indicating"
+write(*,*) "-year (4 digits)"
+write(*,*) "-month (2 digits)" 
+write(*,*) "-day (2 digits)"
+write(*,*) "-hour (2 digits)"
+increment=1
+do while (increment.ne.FILES+1)
+    read(*,*) str
+    read(*,*) str2
+    read(*,*) str3
+    read(*,*) str4
+    date_out(increment)=str//'-'//str2//'-'//str3//'_'//str4
+    date_out2(increment)=str//'-'//str2//'-'//str3//'_'//str4//':00:00'
+    print *,date_out(increment)
+    write(*,*) "ok? (y/n)"
+    read(*,*) answer
+    if (answer.eq.'n') then
+        write(*,*) "please write again"
+    else 
+        write(*,*) "next one, please"    
+        increment=increment+1
+    endif
+enddo
+
+
+!!-------------------
+!! Get array sizes
+!!-------------------
+ierr=NF_INQ_DIMID(nid,"Time",timedim)
+ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
+  write(*,*) "timelen: ",timelen
+ierr=NF_INQ_DIMID(nid,"latitude",latdim)
+ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
+  write(*,*) "latlen: ",latlen
+ierr=NF_INQ_DIMID(nid,"longitude",londim)
+ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
+  write(*,*) "lonlen: ",lonlen
+ierr=NF_INQ_DIMID(nid,"altitude",altdim)
+ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
+  write(*,*) "altlen: ",altlen
+
+
+
+!!-------------------------
+!! Allocate local arrays
+!!-------------------------
+allocate(eta_gcm(lonlen,latlen,altlen,timelen))
+allocate(tfile(lonlen,latlen,altlen,timelen))
+allocate(tfileorig(lonlen,latlen,altlen,timelen))
+allocate(ufile(lonlen,latlen,altlen,timelen))
+allocate(ufileorig(lonlen,latlen,altlen,timelen))
+allocate(vfile(lonlen,latlen,altlen,timelen))
+allocate(vfileorig(lonlen,latlen,altlen,timelen))
+allocate(rfile(lonlen,latlen,altlen,timelen))
+allocate(hfile(lonlen,latlen,altlen,timelen))  
+allocate(psfile(lonlen,latlen,timelen)) 
+allocate(tsfile(lonlen,latlen,timelen))
+allocate(tnsfile(lonlen,latlen,timelen))
+allocate(unsfile(lonlen,latlen,timelen))
+allocate(vnsfile(lonlen,latlen,timelen))
+allocate(emissfile(lonlen,latlen,timelen))
+allocate(co2icefile(lonlen,latlen,timelen))
+allocate(ghtsfile(lonlen,latlen))	!! no scan axis
+allocate(vide(lonlen,latlen))
+allocate(ones(lonlen,latlen))
+allocate(lat(latlen), lon(lonlen), alt(altlen), time(timelen))
+allocate(aps(altlen),bps(altlen),levels(altlen))
+
+tfile(:,:,:,:)=0 
+tfileorig(:,:,:,:)=0
+ufileorig(:,:,:,:)=0
+vfileorig(:,:,:,:)=0
+ufile(:,:,:,:)=0 
+vfile(:,:,:,:)=0 
+rfile(:,:,:,:)=0 
+hfile(:,:,:,:)=0
+psfile(:,:,:)=0 
+tsfile(:,:,:)=0
+emissfile(:,:,:)=0
+co2icefile(:,:,:)=0
+tnsfile(:,:,:)=0
+unsfile(:,:,:)=0
+vnsfile(:,:,:)=0
+ghtsfile(:,:)=0
+vide(:,:)=0
+ones(:,:)=0
+
+
+!!-------------------
+!! Read dimensions
+!!-------------------
+
+print *,'>>> Read dimensions !'
+
+print *,'Time'
+   ierr = NF_INQ_VARID (nid, "Time",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <Time> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, time)
+#endif
+   print *,time(1),' ... to ... ',time(timelen)
+
+print *,'Latitude'
+   ierr = NF_INQ_VARID (nid, "latitude",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <latitude> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lat)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, lat)
+#endif
+   print *,lat(1),' ... to ... ',lat(latlen),' ... step: ',lat(latlen)-lat(latlen-1)
+
+print *,'Longitude'
+   ierr = NF_INQ_VARID (nid, "longitude",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <longitude> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lon)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, lon)
+#endif
+   print *,lon(1),' ... to ... ',lon(lonlen),' ... step: ',lon(lonlen)-lon(lonlen-1)
+
+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)
+
+
+
+
+!!-------------------------------------------
+!! Reading 2D and 3D meteorological fields
+!!-------------------------------------------
+
+IF (blank .EQV. .false.) THEN
+
+print *,'>>> Read 2D optional fields !'
+
+print *,'Emissivity'
+   ierr = NF_INQ_VARID (nid,"emiss",nvarid)
+IF (ierr .NE. NF_NOERR) THEN
+	PRINT *, '...warning: not found in diagfi !'
+        PRINT *, '...will be filled with a prescribed value', emiss_prescribed
+	emissfile(:,:,:)=emiss_prescribed	
+ELSE
+#ifdef NC_DOUBLE
+	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, emissfile)
+#else
+	ierr = NF_GET_VAR_REAL(nid, nvarid, emissfile)
+#endif
+ENDIF   
+
+print *,'CO2 ice'
+   ierr = NF_INQ_VARID (nid,"co2ice",nvarid)
+IF (ierr .NE. NF_NOERR) THEN
+	PRINT *, '...warning: not found in diagfi !'
+	PRINT *, '...will be filled with a prescribed value', co2ice_prescribed
+	co2icefile(:,:,:)=co2ice_prescribed	
+ELSE
+#ifdef NC_DOUBLE
+	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, co2icefile)
+#else
+	ierr = NF_GET_VAR_REAL(nid, nvarid, co2icefile)
+#endif
+ENDIF
+
+print *,'>>> Read 2D surface fields !'
+
+print *,'Surface Pressure'
+   ierr = NF_INQ_VARID (nid,"ps",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <ps> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, psfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, psfile)
+#endif
+
+print *,'Ground Temperature'
+   ierr = NF_INQ_VARID (nid,"tsurf",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <tsurf> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, tsfile)
+#endif
+
+
+!!"atmospheric" surface temperature is taken
+!!from original diagfi.nc first level
+!!... level is ~3-5 meters
+print *,'Near-Surface Temperature'
+   ierr = NF_INQ_VARID (nid,"temp",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <temp> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfileorig)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, tfileorig)
+#endif
+   tnsfile=tfileorig(:,:,1,:)
+
+!!"atmospheric" surface u is taken
+!!from original diagfi.nc first level
+!!... level is ~3-5 meters
+print *,'Near-Surface Zonal Wind'
+   ierr = NF_INQ_VARID (nid,"u",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+     PRINT *, "Error: Readmeteo <u> not found"
+     stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufileorig)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, ufileorig)
+#endif
+   unsfile=ufileorig(:,:,1,:)
+
+!!"atmospheric" surface v is taken
+!!from original diagfi.nc first level
+!!... level is ~3-5 meters
+print *,'Near-Surface Meridional Wind'
+   ierr = NF_INQ_VARID (nid,"v",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+     PRINT *, "Error: Readmeteo <v> not found"
+     stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfileorig)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, vfileorig)
+#endif
+   vnsfile=vfileorig(:,:,1,:)
+
+
+   print *,'Geopotential height at the ground'
+   ierr = NF_INQ_VARID (nid,"phisinit",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <phisinit> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ghtsfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, ghtsfile)
+#endif
+!**** from geopotential to geopotential height
+ghtsfile=ghtsfile/grav
+!****
+
+print *,'>>> Read 3D meteorological fields ! - This may take some time ...'
+
+print *,'Temperature'
+   ierr = NF_INQ_VARID (nid,"temp",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <t> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, tfile)
+#endif
+
+print *,'Zonal wind'   
+   ierr = NF_INQ_VARID (nid,"u",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <u> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, ufile)
+#endif
+
+print *,'Meridional wind'
+   ierr = NF_INQ_VARID (nid,"v",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <v> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, vfile)
+#endif
+
+ENDIF
+
+
+!!-----------------------------
+!! Loop on the written files
+!! >>> what follows is pretty repetitive ... 
+!!        gotta do something (one more loop?)
+!!-----------------------------
+
+!!! Equivalent eta levels for WRF interpolation in initialize_real
+!print *,'Computing equivalent eta levels'
+!	
+!	!ptop will be passed through RH(surface)
+!	ptop=MINVAL(aps(altlen)+bps(altlen)*psfile(:,:,:))
+!	print *, 'ptop', ptop
+!
+!print *, 'sample: equivalent eta levels at i=1,j=1'
+!DO k = 1,altlen
+!        levels(k)=-k
+!        eta_gcm(:,:,k,:)=(aps(k)+bps(k)*psfile(:,:,:)-ptop)/(psfile(:,:,:)-ptop)
+!        print *, eta_gcm(1,1,k,1) 
+!END DO
+
+!!---better to pass pressure values 
+!!---the eta treatment is done in module_initialize_real
+DO k = 1,altlen
+        levels(k)=-k
+        !!dummy decreasing levels
+END DO
+
+
+
+
+!! Dummy surface level is XLVL=200100.
+
+
+DO l=1,FILES
+
+
+!!---------------------------------------------
+!! Write file in intermediate format for WPS
+!! 1. Surface data
+!!---------------------------------------------
+
+!
+! These variables remain the same in the "loop"
+!
+HDATE=date_out(l)
+IFV=4
+XFCST=0.
+SOURCE='LMD GCM'
+NX=lonlen
+NY=latlen
+ALLOCATE(SLAB(NX,NY))
+IPROJ=0
+STARTLOC='SWCORNER'
+STARTLAT=lat(1) 
+STARTLON=lon(1)
+DELTALAT=lat(latlen)-lat(latlen-1)
+DELTALON=lon(lonlen)-lon(lonlen-1)
+!
+! Open the data destination file
+!
+output=ident//"/"//ident//":"//date_out2(l)       
+open(UNIT=1,			&
+	FILE=output,		&
+	STATUS='REPLACE',	&
+	FORM='UNFORMATTED',	&
+	ACCESS='SEQUENTIAL')
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='T'
+UNITS='K'
+DESC='Atmospheric temperature'
+XLVL=200100.
+!SLAB=tsfile(:,:,time_out(l))
+!SLAB=tfileorig(:,:,1,time_out(l))
+SLAB=tnsfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='U'
+UNITS='m s{-1}'
+DESC='Zonal wind'
+XLVL=200100.
+!SLAB=ufile(:,:,1,time_out(l))
+!SLAB=ufileorig(:,:,1,time_out(l))
+SLAB=unsfile(:,:,time_out(l))
+        ! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='V'
+UNITS='m s{-1}'
+DESC='Meridional wind'
+XLVL=200100.
+!SLAB=vfile(:,:,1,time_out(l))
+!SLAB=vfileorig(:,:,1,time_out(l))
+SLAB=vnsfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='RH'
+UNITS=''
+DESC='Customized 2D static field'
+XLVL=200100.
+!SLAB=co2icefile(:,:,time_out(l))
+SLAB=vide(:,:)
+!SLAB=vide(:,:)+ptop
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SOILHGT'
+UNITS='m' 
+DESC='Terrain field of source analysis'  !Ground geopotential height
+XLVL=200100.
+SLAB=ghtsfile(:,:)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='PSFC'
+UNITS='Pa'
+DESC='Surface Pressure'
+XLVL=200100.
+SLAB=psfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='ST000010'
+UNITS=''
+DESC='Emissivity'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=emissfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='ST010040'
+UNITS=''
+DESC='CO2 ice'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=co2icefile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='ST040100'
+UNITS='K'
+DESC='T of 40-100 cm ground layer'
+XLVL=200100.
+SLAB=vide(:,:)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='ST100200'
+UNITS='K'
+DESC='T of 100-200 cm ground layer'
+XLVL=200100.
+SLAB=vide(:,:)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='LANDSEA'
+UNITS='0/1 Flag'
+DESC='Land/Sea flag'
+XLVL=200100.
+SLAB=ones(:,:)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SKINTEMP'
+UNITS='K'
+DESC='Ground temperature'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=tsfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SM000010'
+UNITS='fraction'
+DESC='Relative humidity'
+XLVL=200100.
+SLAB=vide(:,:)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SM010040'
+UNITS='fraction'
+DESC='Relative humidity'
+XLVL=200100.
+SLAB=vide(:,:)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SM040100'
+UNITS='fraction'
+DESC='Relative humidity'
+XLVL=200100.
+SLAB=vide(:,:)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SM100200'
+UNITS='fraction'
+DESC='Relative humidity'
+XLVL=200100.
+SLAB=vide(:,:)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
+	! ... Data
+	write(1) SLAB
+print *,'The field '//DESC//' was written to '//output
+
+
+
+!!----------------------------------------------------
+!! Write file in intermediate format for WPS
+!! 2. 3D meteorological data
+!! >>> same stuff, but handle vertical profile
+!! >>> !! XLVL must be decreasing (pressure levels) 
+!!----------------------------------------------------
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='T'
+UNITS='K'
+DESC='Atmospheric temperature'
+DO k = 1,altlen
+	XLVL=levels(k)
+	SLAB=tfile(:,:,k,time_out(l))
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='U'
+UNITS='m s{-1}'
+DESC='Zonal wind'
+DO k = 1,altlen
+        XLVL=levels(k)
+	SLAB=ufile(:,:,k,time_out(l))
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='V'
+UNITS='m s{-1}'
+DESC='Meridional wind'
+DO k = 1,altlen
+        XLVL=levels(k)
+	SLAB=vfile(:,:,k,time_out(l))
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='RH'
+UNITS=''
+DESC='Customized 2D static field'
+DESC='Eta-levels from the GCM'
+DESC='Pressure values from the GCM'
+DO k = 1,altlen
+        XLVL=levels(k)
+	!SLAB=rfile(:,:,k,time_out(l))
+        !SLAB=eta_gcm(:,:,k,time_out(l))
+        SLAB=aps(k)+bps(k)*psfile(:,:,time_out(l))
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='HGT'
+UNITS='m'
+DESC='Height'
+DO k = 1,altlen
+        XLVL=levels(k)
+!!*******
+!! PROVISOIRE: normalement, il faudrait la hauteur geopotentielle
+!!*******
+!!SLAB=hfile(:,:,k,time_out(l))
+!!SLAB=vide(:,:)
+        SLAB=10000.*alog(610./(aps(k)+bps(k)*psfile(:,:,time_out(l))))
+!!however, not used by initialize_real
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+print *,'The field '//DESC//' was written to '//output
+
+print *,'****Close file '//output
+close(1)
+
+DEALLOCATE(SLAB)
+
+END DO !! end of the files loop
+
+
+!!-------------------------
+!! DeAllocate local arrays
+!!-------------------------
+deallocate(eta_gcm)
+deallocate(tfile)
+deallocate(tfileorig)
+deallocate(ufile)
+deallocate(ufileorig)
+deallocate(vfile)
+deallocate(vfileorig)
+deallocate(rfile)
+deallocate(hfile)  
+deallocate(psfile) 
+deallocate(tsfile)
+deallocate(tnsfile)
+deallocate(unsfile)
+deallocate(vnsfile)
+deallocate(emissfile)
+deallocate(co2icefile)
+deallocate(ghtsfile)	!! no scan axis
+deallocate(vide)
+deallocate(ones)
+deallocate(lat, lon, alt, time)
+deallocate(aps,bps,levels)
+
+
+print *, 'End of pre-WPS process !'
+
+
+end
Index: trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/obsolete/readmeteo_newphys.F90
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/obsolete/readmeteo_newphys.F90	(revision 72)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/obsolete/readmeteo_newphys.F90	(revision 72)
@@ -0,0 +1,1491 @@
+program readmeteo
+
+implicit none
+include "netcdf.inc"
+
+!------------------------------------------------------------------------!
+! Readmeteo prepares an initial state for the WPS pre-processor of WRF   !
+!                                                                        !
+! Input is a diagfi.nc NETCDF file from a LMD/GCM simulation             !
+!                                                                        !
+! Outputs are WRFSI intermediate format files ready for metgrid.exe      !
+!       http://www.mmm.ucar.edu/wrf/OnLineTutorial/WPS/IM_si.htm         !
+!                                                                        !
+! **** Please create a WPSFEED folder (or a symbolic link) ****          !       
+!                                                                        !
+! A. Spiga - 16/04/2007                                                  !
+!	     22/05/2007 : need to run zrecast as a preliminary           !
+!            07/06/2007 : external parameter file 'readmeteo.def         !
+!	     30/07/2007 : manage additional variables (tsoil, emiss,...) !	
+!            19/10/2007 : no need to run zrecast at all (eta levels)     !   
+!		12/2007 : include co2ice and emissivity			 ! 
+!               02/2008 : include water vapor and ice                    !
+!               01/2010 : possible use of diagsoil for new physics       ! 
+!               02/2011 : add dust tracers                               !
+!									 !
+! spiga@lmd.jussieu.fr							 !	
+!------------------------------------------------------------------------!
+
+
+
+!***************************************************************************
+!***************************************************************************
+REAL, PARAMETER :: emiss_prescribed=0.95
+REAL, PARAMETER :: co2ice_prescribed=0.
+CHARACTER (LEN=3), PARAMETER :: ident='LMD'
+
+
+
+!***************************************************************************
+!***************************************************************************
+
+REAL :: ptop
+REAL, PARAMETER :: grav=3.72
+LOGICAL, PARAMETER :: blank=.false.
+
+
+
+!! Variables to be written in the output file
+!! *** NB: conformity with metgrid/src/read_met_module.F90
+CHARACTER (LEN=33) :: output
+INTEGER :: IFV
+CHARACTER (LEN=24) :: HDATE
+REAL :: XFCST
+CHARACTER (LEN=32) :: SOURCE
+CHARACTER (LEN=9) :: FIELD
+CHARACTER (LEN=25) :: UNITS
+CHARACTER (LEN=46) :: DESC
+REAL :: XLVL
+INTEGER :: NX,NY,IPROJ
+CHARACTER (LEN=8) :: STARTLOC
+REAL :: STARTLAT,STARTLON,DELTALAT,DELTALON
+REAL, POINTER, DIMENSION(:,:) :: SLAB
+
+!! Variables related to NETCDF format
+integer :: nid,nid2,nid3,nvarid,ierr,ierr2
+integer :: timedim,londim,latdim,altdim,altdim2
+integer :: rlatvdim,rlonudim
+integer :: timelen,lonlen,latlen,altlen,altlen2
+integer :: rlatvlen,rlonulen
+integer :: soillen,soildim
+integer :: nphys
+integer, dimension(4) :: corner,edges
+
+!! Intermediate data arrays
+integer :: k,l,m,n,p
+real, dimension(:), allocatable :: lat,lon,time,alt,aps,bps,levels
+real, dimension(:,:), allocatable :: vide,ones,ghtsfile
+real, dimension(:,:), allocatable :: interm
+real, dimension(:,:,:), allocatable :: gwparam
+real, dimension(:,:,:), allocatable :: psfile,tsfile
+real, dimension(:,:,:), allocatable :: emissfile,co2icefile
+real, dimension(:,:,:), allocatable :: tnsfile,unsfile,vnsfile
+real, dimension(:,:,:,:), allocatable :: tfile,ufile,vfile,rfile,hfile
+real, dimension(:,:,:,:), allocatable :: eta_gcm
+real, dimension(:,:,:,:), allocatable :: tfileorig,ufileorig,vfileorig
+real, dimension(:,:,:,:), allocatable :: tsoilfile, dsoilfile, isoilfile
+real, dimension(:,:,:,:), allocatable :: waterfile, watericefile
+real, dimension(:,:,:,:), allocatable :: dustfile,dustnfile
+
+!! Reading the parameter file 
+integer :: tmp,increment,FILES
+integer :: tmp2,tmp3,tmp4
+character*1 :: answer
+character*4 :: str
+character*2 :: str2, str3, str4
+integer, dimension(:), allocatable :: time_out
+character*13, dimension(:), allocatable :: date_out
+character*19, dimension(:), allocatable :: date_out2
+!***************************************************************************
+!***************************************************************************
+
+
+!!---------------------
+!! Open the datafile
+!!---------------------
+ierr=NF_OPEN ("input_diagfi.nc",NF_NOWRITE,nid)
+IF (ierr.NE.NF_NOERR) THEN
+   write(*,*)'**** Please create a symbolic link called input_diagfi.nc'
+   CALL ABORT
+ENDIF
+
+
+!!--------------------------
+!! Ask for data references
+!!--------------------------
+!write(*,*) "Create n files. What is n ?"
+read(*,*) FILES
+allocate(time_out(FILES))
+allocate(date_out(FILES))
+allocate(date_out2(FILES))
+
+!write(*,*) ""
+!write(*,*) "INPUT: Diagfi time"
+!write(*,*) "list of n subscripts, separated by <Return>s"
+increment=1
+do while (increment.ne.FILES+1)
+    read(*,*) tmp
+    time_out(increment)=tmp
+    increment=increment+1
+enddo
+
+!write(*,*) ""
+!write(*,*) "OUTPUT: WRF time"
+!write(*,*) "fill 4 lines indicating"
+!write(*,*) "-year (4 digits)"
+!write(*,*) "-month (2 digits)" 
+!write(*,*) "-day (2 digits)"
+!write(*,*) "-hour (2 digits)"
+increment=1
+do while (increment.ne.FILES+1)
+    read(*,*) str
+    read(*,*) str2
+    read(*,*) str3
+    read(*,*) str4
+    date_out(increment)=str//'-'//str2//'-'//str3//'_'//str4
+    date_out2(increment)=str//'-'//str2//'-'//str3//'_'//str4//':00:00'
+    !print *,date_out(increment)
+    !write(*,*) "ok? (y/n)"
+    read(*,*) answer
+    if (answer.eq.'n') then
+        !write(*,*) "please write again"
+    else 
+        !write(*,*) "next one, please"    
+        increment=increment+1
+    endif
+enddo
+
+
+!!-------------------
+!! Get array sizes
+!!-------------------
+SELECT CASE(ident)
+CASE('LMD')
+ierr=NF_INQ_DIMID(nid,"Time",timedim)
+CASE('OXF')
+ierr=NF_INQ_DIMID(nid,"time",timedim)
+END SELECT
+ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
+  write(*,*) "timelen: ",timelen
+
+SELECT CASE(ident)
+CASE('LMD')
+ierr=NF_INQ_DIMID(nid,"latitude",latdim)
+CASE('OXF')
+ierr=NF_INQ_DIMID(nid,"lat",latdim)
+END SELECT
+ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
+  write(*,*) "latlen: ",latlen
+
+SELECT CASE(ident)
+CASE('LMD')
+ierr=NF_INQ_DIMID(nid,"longitude",londim)
+CASE('OXF')
+ierr=NF_INQ_DIMID(nid,"lon",londim)
+END SELECT
+ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
+  write(*,*) "lonlen: ",lonlen
+
+SELECT CASE(ident)
+CASE('LMD')
+ierr=NF_INQ_DIMID(nid,"altitude",altdim)
+CASE('OXF')
+ierr=NF_INQ_DIMID(nid,"sigma",altdim)
+END SELECT
+ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
+  write(*,*) "altlen: ",altlen
+
+
+
+!!-------------------------
+!! Allocate local arrays
+!!-------------------------
+allocate(eta_gcm(lonlen,latlen,altlen,timelen))
+allocate(tfile(lonlen,latlen,altlen,timelen))
+allocate(tsoilfile(lonlen,latlen,altlen,timelen))
+allocate(dsoilfile(lonlen,latlen,altlen,timelen))
+allocate(isoilfile(lonlen,latlen,altlen,timelen))
+allocate(tfileorig(lonlen,latlen,altlen,timelen))
+allocate(ufile(lonlen,latlen,altlen,timelen))
+allocate(ufileorig(lonlen,latlen,altlen,timelen))
+allocate(vfile(lonlen,latlen,altlen,timelen))
+allocate(vfileorig(lonlen,latlen,altlen,timelen))
+allocate(rfile(lonlen,latlen,altlen,timelen))
+allocate(hfile(lonlen,latlen,altlen,timelen))  
+allocate(waterfile(lonlen,latlen,altlen,timelen))
+allocate(watericefile(lonlen,latlen,altlen,timelen))
+allocate(dustfile(lonlen,latlen,altlen,timelen))
+allocate(dustnfile(lonlen,latlen,altlen,timelen))
+allocate(psfile(lonlen,latlen,timelen)) 
+allocate(tsfile(lonlen,latlen,timelen))
+allocate(tnsfile(lonlen,latlen,timelen))
+allocate(unsfile(lonlen,latlen,timelen))
+allocate(vnsfile(lonlen,latlen,timelen))
+allocate(emissfile(lonlen,latlen,timelen))
+allocate(co2icefile(lonlen,latlen,timelen))
+allocate(interm(lonlen,latlen))
+allocate(gwparam(lonlen,latlen,5))
+allocate(ghtsfile(lonlen,latlen))  !! no scan axis
+allocate(vide(lonlen,latlen))
+allocate(ones(lonlen,latlen))
+allocate(lat(latlen), lon(lonlen), alt(altlen), time(timelen))
+allocate(aps(altlen),bps(altlen),levels(altlen))
+
+tfile(:,:,:,:)=0
+tsoilfile(:,:,:,:)=0 
+isoilfile(:,:,:,:)=0 
+dsoilfile(:,:,:,:)=0 
+tfileorig(:,:,:,:)=0
+ufileorig(:,:,:,:)=0
+vfileorig(:,:,:,:)=0
+ufile(:,:,:,:)=0 
+vfile(:,:,:,:)=0 
+rfile(:,:,:,:)=0 
+hfile(:,:,:,:)=0
+waterfile(:,:,:,:)=0
+watericefile(:,:,:,:)=0
+dustfile(:,:,:,:)=0
+dustnfile(:,:,:,:)=0
+psfile(:,:,:)=0 
+tsfile(:,:,:)=0
+emissfile(:,:,:)=0
+co2icefile(:,:,:)=0
+tnsfile(:,:,:)=0
+unsfile(:,:,:)=0
+vnsfile(:,:,:)=0
+interm(:,:)=0
+gwparam(:,:,:)=0
+ghtsfile(:,:)=0
+vide(:,:)=0
+ones(:,:)=0
+
+
+!!-------------------
+!! Read dimensions
+!!-------------------
+
+print *,'>>> Read dimensions !'
+
+print *,'Time'
+SELECT CASE(ident)
+CASE('LMD')
+   ierr = NF_INQ_VARID (nid, "Time",nvarid)
+CASE('OXF')
+   ierr = NF_INQ_VARID (nid, "time",nvarid)
+END SELECT
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <Time> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, time)
+#endif
+   print *,time(1),' ... to ... ',time(timelen)
+
+print *,'Latitude'
+SELECT CASE(ident)
+CASE('LMD')
+   ierr = NF_INQ_VARID (nid, "latitude",nvarid)
+CASE('OXF')
+   ierr = NF_INQ_VARID (nid, "lat",nvarid)
+END SELECT
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <latitude> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lat)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, lat)
+#endif
+   print *,lat(1),' ... to ... ',lat(latlen),' ... step: ',lat(latlen)-lat(latlen-1)
+
+print *,'Longitude'
+SELECT CASE(ident)
+CASE('LMD')
+   ierr = NF_INQ_VARID (nid, "longitude",nvarid)
+CASE('OXF')
+   ierr = NF_INQ_VARID (nid, "lon",nvarid)
+END SELECT
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <longitude> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lon)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, lon)
+#endif
+   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
+
+
+!!-------------------------------------------
+!! Reading 2D and 3D meteorological fields
+!!-------------------------------------------
+
+IF (blank .EQV. .false.) THEN
+
+print *,'>>> Read 2D optional fields !'
+
+print *,'Emissivity'
+   ierr = NF_INQ_VARID (nid,"emis",nvarid)
+IF (ierr .NE. NF_NOERR) THEN
+        PRINT *, '...warning: not found in diagfi !'
+        PRINT *, '...will be filled with a prescribed value', emiss_prescribed
+        emissfile(:,:,:)=emiss_prescribed
+ELSE
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, emissfile)
+#else
+        ierr = NF_GET_VAR_REAL(nid, nvarid, emissfile)
+#endif
+ENDIF   
+
+print *,'CO2 ice'
+   ierr = NF_INQ_VARID (nid,"co2ice",nvarid)
+IF (ierr .NE. NF_NOERR) THEN
+        PRINT *, '...warning: not found in diagfi !'
+        PRINT *, '...will be filled with a prescribed value', co2ice_prescribed
+        co2icefile(:,:,:)=co2ice_prescribed
+ELSE
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, co2icefile)
+#else
+        ierr = NF_GET_VAR_REAL(nid, nvarid, co2icefile)
+#endif
+ENDIF
+
+print *,'>>> Read 2D surface fields !'
+
+print *,'Surface Pressure'
+   ierr = NF_INQ_VARID (nid,"ps",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <ps> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, psfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, psfile)
+#endif
+
+print *,'Ground Temperature'
+   ierr = NF_INQ_VARID (nid,"tsurf",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <tsurf> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, tsfile)
+#endif
+
+
+!!"atmospheric" surface temperature is taken
+!!from original diagfi.nc first level
+!!... level is ~3-5 meters
+print *,'Near-Surface Temperature'
+   ierr = NF_INQ_VARID (nid,"temp",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      ierr = NF_INQ_VARID (nid,"t",nvarid)
+        IF (ierr .NE. NF_NOERR) THEN
+           PRINT *, "Error: Readmeteo <temp> not found"
+           stop
+        ENDIF
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfileorig)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, tfileorig)
+#endif
+   tnsfile=tfileorig(:,:,1,:)
+
+!!"atmospheric" surface u is taken
+!!from original diagfi.nc first level
+!!... level is ~3-5 meters
+print *,'Near-Surface Zonal Wind'
+   ierr = NF_INQ_VARID (nid,"u",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+     PRINT *, "Error: Readmeteo <u> not found"
+     stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufileorig)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, ufileorig)
+#endif
+   unsfile=ufileorig(:,:,1,:)
+
+!!"atmospheric" surface v is taken
+!!from original diagfi.nc first level
+!!... level is ~3-5 meters
+print *,'Near-Surface Meridional Wind'
+   ierr = NF_INQ_VARID (nid,"v",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+     PRINT *, "Error: Readmeteo <v> not found"
+     stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfileorig)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, vfileorig)
+#endif
+   vnsfile=vfileorig(:,:,1,:)
+
+SELECT CASE(ident)
+CASE('LMD')
+   print *,'Geopotential height at the ground'
+   ierr = NF_INQ_VARID (nid,"phisinit",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <phisinit> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ghtsfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, ghtsfile)
+#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 ...'
+
+print *,'Temperature'
+   ierr = NF_INQ_VARID (nid,"temp",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      ierr = NF_INQ_VARID (nid,"t",nvarid)
+        IF (ierr .NE. NF_NOERR) THEN
+          PRINT *, "Error: Readmeteo <t> not found"
+          stop
+        ENDIF
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, tfile)
+#endif
+
+print *,'Zonal wind'   
+   ierr = NF_INQ_VARID (nid,"u",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <u> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, ufile)
+#endif
+
+print *,'Meridional wind'
+   ierr = NF_INQ_VARID (nid,"v",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <v> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, vfile)
+#endif
+
+
+!!------------------------
+!! special water stuff
+!!------------------------
+    print *,'Water vapor'
+    ierr=NF_INQ_VARID(nid,"q02",nvarid)
+    if (ierr.ne.NF_NOERR) then
+      write(*,*) "...No q02 - Water vapor set to 0"
+      waterfile(:,:,:,:)=0.  
+    endif
+    ierr=NF_GET_VAR_REAL(nid,nvarid,waterfile)
+
+    print *,'Water ice'
+    ierr=NF_INQ_VARID(nid,"q01",nvarid)
+    if (ierr.ne.NF_NOERR) then
+      write(*,*) "...No q01 - Water ice set to 0"  
+      watericefile(:,:,:,:)=0.
+    endif
+    ierr=NF_GET_VAR_REAL(nid,nvarid,watericefile)
+!!------------------------
+!! special water stuff
+!!------------------------
+
+!!------------------------
+!! special dust stuff
+!!------------------------
+    print *,'Dust mass'
+    ierr=NF_INQ_VARID(nid,"dustq",nvarid)
+    if (ierr.ne.NF_NOERR) then
+      write(*,*) "...No dustq - Dust mass set to 0"
+      dustfile(:,:,:,:)=0.
+    endif
+    ierr=NF_GET_VAR_REAL(nid,nvarid,dustfile)
+
+    print *,'Dust number'
+    ierr=NF_INQ_VARID(nid,"dustN",nvarid)
+    if (ierr.ne.NF_NOERR) then
+      write(*,*) "...No dustN - Dust number set to 0"
+      dustnfile(:,:,:,:)=0.
+    endif
+    ierr=NF_GET_VAR_REAL(nid,nvarid,dustnfile)
+!!------------------------
+!! special dust stuff
+!!------------------------
+
+
+!SELECT CASE(ident)
+!CASE('LMD')
+
+    print *,'Soil temperatures'
+    ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
+    if (ierr.ne.NF_NOERR) then
+        write(*,*) "...No tsoil - Soil temperatures set isothermal with surface temperature"
+        DO l=1,altlen
+                tsoilfile(:,:,l,:)=tsfile(:,:,:)
+        ENDDO
+    endif
+    ierr=NF_GET_VAR_REAL(nid,nvarid,tsoilfile)
+  
+!!!!!!!!
+!!!!!!!! new physics (but still compatible with old physics)
+
+    print *,'Soil depths'
+    ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
+    if (ierr.ne.NF_NOERR) then
+        write(*,*) "...No soildepth - Set to default"  !!! see soil_settings in LMD physics
+        DO l=1,altlen
+                dsoilfile(1,1,l,:)=-999.
+        ENDDO
+    endif
+    ierr=NF_GET_VAR_REAL(nid,nvarid,dsoilfile(1,1,:,1))
+    print *, dsoilfile(1,1,:,1)
+    DO m=1,lonlen
+    DO n=1,latlen
+    DO p=1,timelen
+      dsoilfile(m,n,:,p) = dsoilfile(1,1,:,1) 
+    ENDDO
+    ENDDO
+    ENDDO
+
+    print *,'Soil thermal inertia'
+    ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
+    if (ierr.ne.NF_NOERR) then
+        write(*,*) "...No soil therm. inert. - Set to -999"
+        DO l=1,altlen
+                isoilfile(:,:,l,:)=-999.
+        ENDDO
+    endif
+    ierr=NF_GET_VAR_REAL(nid,nvarid,isoilfile)
+
+!!!!!!!!
+!!!!!!!! new physics
+
+
+!END SELECT
+
+ierr=NF_CLOSE(nid)
+
+!!!!lott stuff
+!print *,'get lott'
+!ierr=NF_OPEN("init_lott.nc",NF_NOWRITE,nid3)
+!ierr=NF_INQ_VARID(nid3,"ZMEA",nvarid)
+!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
+!gwparam(:,:,1)=interm(:,:)
+!ierr=NF_INQ_VARID(nid3,"ZSTD",nvarid)
+!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
+!gwparam(:,:,2)=interm(:,:)
+!ierr=NF_INQ_VARID(nid3,"ZSIG",nvarid)
+!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
+!gwparam(:,:,3)=interm(:,:)
+!ierr=NF_INQ_VARID(nid3,"ZGAM",nvarid)
+!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
+!gwparam(:,:,4)=interm(:,:)
+!ierr=NF_INQ_VARID(nid3,"ZTHE",nvarid)
+!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
+!gwparam(:,:,5)=interm(:,:)
+!ierr=NF_CLOSE(nid3)
+
+ENDIF
+
+
+!!-----------------------------
+!! Loop on the written files
+!! >>> what follows is pretty repetitive ... 
+!!        gotta do something (one more loop?)
+!!-----------------------------
+
+!!! Equivalent eta levels for WRF interpolation in initialize_real
+!print *,'Computing equivalent eta levels'
+!	
+!	!ptop will be passed through RH(surface)
+!	ptop=MINVAL(aps(altlen)+bps(altlen)*psfile(:,:,:))
+!	print *, 'ptop', ptop
+!
+!print *, 'sample: equivalent eta levels at i=1,j=1'
+!DO k = 1,altlen
+!        levels(k)=-k
+!        eta_gcm(:,:,k,:)=(aps(k)+bps(k)*psfile(:,:,:)-ptop)/(psfile(:,:,:)-ptop)
+!        print *, eta_gcm(1,1,k,1) 
+!END DO
+
+!!---better to pass pressure values 
+!!---the eta treatment is done in module_initialize_real
+DO k = 1,altlen
+        levels(k)=-k
+        !!dummy decreasing levels
+END DO
+
+
+
+
+!! Dummy surface level is XLVL=200100.
+
+
+DO l=1,FILES
+
+
+!!---------------------------------------------
+!! Write file in intermediate format for WPS
+!! 1. Surface data
+!!---------------------------------------------
+!!
+!! a mettre pour tous sinon
+!!     WRF_DEBUG: Warning DIM             4 , NAME num_metgrid_levels REDIFINED  by
+!!            var DUSTN            26           25  in wrf_io.F90 line         2349
+!!
+
+!
+! These variables remain the same in the "loop"
+!
+HDATE=date_out(l)
+IFV=4
+XFCST=0.
+SOURCE=ident
+NX=lonlen
+NY=latlen
+ALLOCATE(SLAB(NX,NY))
+IPROJ=0
+STARTLOC='SWCORNER'
+STARTLAT=lat(1) 
+STARTLON=lon(1)
+DELTALAT=lat(latlen)-lat(latlen-1)
+DELTALON=lon(lonlen)-lon(lonlen-1)
+!
+! Open the data destination file
+!
+output="./WPSFEED/"//ident//":"//date_out2(l)       
+open(UNIT=1,			&
+	FILE=output,		&
+	STATUS='REPLACE',	&
+	FORM='UNFORMATTED',	&
+	ACCESS='SEQUENTIAL')
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='T'
+UNITS='K'
+DESC='Atmospheric temperature'
+XLVL=200100.
+!SLAB=tsfile(:,:,time_out(l))
+!SLAB=tfileorig(:,:,1,time_out(l))
+SLAB=tnsfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='U'
+UNITS='m s{-1}'
+DESC='Zonal wind'
+XLVL=200100.
+!SLAB=ufile(:,:,1,time_out(l))
+!SLAB=ufileorig(:,:,1,time_out(l))
+SLAB=unsfile(:,:,time_out(l))
+        ! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='V'
+UNITS='m s{-1}'
+DESC='Meridional wind'
+XLVL=200100.
+!SLAB=vfile(:,:,1,time_out(l))
+!SLAB=vfileorig(:,:,1,time_out(l))
+SLAB=vnsfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='RH'
+UNITS=''
+DESC='Customized 2D static field'
+XLVL=200100.
+!SLAB=co2icefile(:,:,time_out(l))
+SLAB=vide(:,:)
+!SLAB=vide(:,:)+ptop
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SOILHGT'
+UNITS='m' 
+DESC='Terrain field of source analysis'  !Ground geopotential height
+XLVL=200100.
+SLAB=ghtsfile(:,:)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='PSFC'
+UNITS='Pa'
+DESC='Surface Pressure'
+XLVL=200100.
+SLAB=psfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='ST000010'
+UNITS=''
+DESC='Emissivity'
+XLVL=200100.
+SLAB=emissfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='ST010040'
+UNITS=''
+DESC='CO2 ice'
+XLVL=200100.
+SLAB=co2icefile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='ST040100'
+UNITS=''
+DESC='ZMEA'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=gwparam(:,:,1)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='ST100200'
+UNITS=''
+DESC='ZSTD'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=gwparam(:,:,2)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='LANDSEA'
+UNITS='0/1 Flag'
+DESC='Land/Sea flag'
+XLVL=200100.
+SLAB=ones(:,:)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SKINTEMP'
+UNITS='K'
+DESC='Ground temperature'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=tsfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SM000010'
+UNITS=''
+DESC='ZSIG'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=gwparam(:,:,3)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SM010040'
+UNITS=''
+DESC='ZGAM'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=gwparam(:,:,4)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SM040100'
+UNITS=''
+DESC='ZTHE'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=gwparam(:,:,5)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SM100200'
+UNITS='fraction'
+DESC='Relative humidity'
+XLVL=200100.
+SLAB=vide(:,:)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='HV'
+UNITS='kg/kg'
+DESC='Water vapor'
+XLVL=200100.
+SLAB=waterfile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='HI'
+UNITS='kg/kg'
+DESC='Water ice'
+XLVL=200100.
+SLAB=watericefile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='TSOIL'
+UNITS='K'
+DESC='Soil temperature'
+XLVL=200100.
+SLAB=tsoilfile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='DSOIL'
+UNITS='m'
+DESC='Soil depths'
+XLVL=200100.
+SLAB=dsoilfile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='ISOIL'
+UNITS='tiu'
+DESC='Soil thermal inertia'
+XLVL=200100.
+SLAB=isoilfile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='DUSTQ'
+UNITS='kg/kg'
+DESC='Dust mixing ratio'
+XLVL=200100.
+SLAB=dustfile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='DUSTN'
+UNITS='part/kg' 
+DESC='Dust number density'
+XLVL=200100.
+SLAB=dustnfile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+
+!!----------------------------------------------------
+!! Write file in intermediate format for WPS
+!! 2. 3D meteorological data
+!! >>> same stuff, but handle vertical profile
+!! >>> !! XLVL must be decreasing (pressure levels) 
+!!----------------------------------------------------
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='T'
+UNITS='K'
+DESC='Atmospheric temperature'
+DO k = 1,altlen
+	XLVL=levels(k)
+	SLAB=tfile(:,:,k,time_out(l))
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='U'
+UNITS='m s{-1}'
+DESC='Zonal wind'
+DO k = 1,altlen
+        XLVL=levels(k)
+	SLAB=ufile(:,:,k,time_out(l))
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='V'
+UNITS='m s{-1}'
+DESC='Meridional wind'
+DO k = 1,altlen
+        XLVL=levels(k)
+	SLAB=vfile(:,:,k,time_out(l))
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='RH'
+UNITS=''
+DESC='Customized 2D static field'
+DESC='Eta-levels from the GCM'
+DESC='Pressure values from the GCM'
+DO k = 1,altlen
+        XLVL=levels(k)
+SELECT CASE(ident)
+CASE('LMD')
+        SLAB=aps(k)+bps(k)*psfile(:,:,time_out(l))
+CASE('OXF')
+	SLAB=bps(k)*psfile(:,:,time_out(l))
+END SELECT
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='HGT'
+UNITS='m'
+DESC='Height'
+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))))
+CASE('OXF')
+	SLAB=10000.*alog(610./(bps(k)*psfile(:,:,time_out(l))))
+END SELECT
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='HV'
+UNITS='kg/kg'
+DESC='Water vapor'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=waterfile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='HI'
+UNITS='kg/kg'
+DESC='Water ice'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=watericefile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='TSOIL'
+UNITS='K'
+DESC='Soil temperature'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=tsoilfile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='DSOIL'
+UNITS='m'
+DESC='Soil depths'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=dsoilfile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='ISOIL'
+UNITS='tiu'
+DESC='Soil thermal inertia'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=isoilfile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='DUSTQ'
+UNITS='kg/kg'
+DESC='Dust mixing ratio'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=dustfile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='DUSTN'
+UNITS='part/kg'
+DESC='Dust number density'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=dustnfile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+print *,'****done file '//output, int(100.*float(l)/float(FILES)), ' % '
+close(1)
+
+DEALLOCATE(SLAB)
+
+END DO !! end of the files loop
+
+
+!!-------------------------
+!! DeAllocate local arrays
+!!-------------------------
+deallocate(eta_gcm)
+deallocate(tfile)
+deallocate(tsoilfile)
+deallocate(isoilfile)
+deallocate(dsoilfile)
+deallocate(tfileorig)
+deallocate(ufile)
+deallocate(ufileorig)
+deallocate(vfile)
+deallocate(vfileorig)
+deallocate(rfile)
+deallocate(hfile)  
+deallocate(waterfile)
+deallocate(watericefile)
+deallocate(dustfile)
+deallocate(dustnfile)
+deallocate(psfile) 
+deallocate(tsfile)
+deallocate(tnsfile)
+deallocate(unsfile)
+deallocate(vnsfile)
+deallocate(emissfile)
+deallocate(co2icefile)
+deallocate(ghtsfile)	!! no scan axis
+deallocate(vide)
+deallocate(ones)
+deallocate(lat, lon, alt, time)
+deallocate(aps,bps,levels)
+
+
+print *, '------------------------'
+print *, 'End of pre-WPS process !'
+print *, '------------------------'
+print *, 'Here is what you should set in namelist.wps:' 
+print *, " start_date = '"//date_out2(1)//"'"
+print *, " end_date = '"//date_out2(FILES)//"'"
+print *, '------------------------'
+print *, 'Any date between those is ok'
+print *, 'for example, second available date is ... '//date_out2(2)
+print *, '------------------------'
+
+end
Index: trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/readmeteo.F90
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/readmeteo.F90	(revision 70)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/readmeteo.F90	(revision 72)
@@ -21,4 +21,6 @@
 !		12/2007 : include co2ice and emissivity			 ! 
 !               02/2008 : include water vapor and ice                    !
+!               01/2010 : possible use of diagsoil for new physics       ! 
+!               02/2011 : add dust tracers, correct surface              !
 !									 !
 ! spiga@lmd.jussieu.fr							 !	
@@ -71,6 +73,6 @@
 
 !! Intermediate data arrays
-integer :: k,l
-real, dimension(:), allocatable :: lat,lon,time,alt,aps,bps,levels
+integer :: k,l,m,n,p
+real, dimension(:), allocatable :: lat,lon,time,alt,aps,bps,levels,dsoilvert
 real, dimension(:,:), allocatable :: vide,ones,ghtsfile
 real, dimension(:,:), allocatable :: interm
@@ -81,8 +83,9 @@
 real, dimension(:,:,:,:), allocatable :: tfile,ufile,vfile,rfile,hfile
 real, dimension(:,:,:,:), allocatable :: eta_gcm
-real, dimension(:,:,:,:), allocatable :: tfileorig,ufileorig,vfileorig
-real, dimension(:,:,:,:), allocatable :: tsoilfile
+!real, dimension(:,:,:,:), allocatable :: tfileorig,ufileorig,vfileorig
+real, dimension(:,:,:,:), allocatable :: tsoilfile, dsoilfile, isoilfile
 real, dimension(:,:,:,:), allocatable :: waterfile, watericefile
-real, dimension(:,:,:), allocatable :: swaterfile, swatericefile
+real, dimension(:,:,:), allocatable :: swatericefile!, swaterfile
+real, dimension(:,:,:,:), allocatable :: dustfile,dustnfile
 
 !! Reading the parameter file 
@@ -202,15 +205,20 @@
 allocate(tfile(lonlen,latlen,altlen,timelen))
 allocate(tsoilfile(lonlen,latlen,altlen,timelen))
-allocate(tfileorig(lonlen,latlen,altlen,timelen))
+allocate(dsoilfile(lonlen,latlen,altlen,timelen))
+allocate(isoilfile(lonlen,latlen,altlen,timelen))
+allocate(dsoilvert(altlen))
+!allocate(tfileorig(lonlen,latlen,altlen,timelen))
 allocate(ufile(lonlen,latlen,altlen,timelen))
-allocate(ufileorig(lonlen,latlen,altlen,timelen))
+!allocate(ufileorig(lonlen,latlen,altlen,timelen))
 allocate(vfile(lonlen,latlen,altlen,timelen))
-allocate(vfileorig(lonlen,latlen,altlen,timelen))
+!allocate(vfileorig(lonlen,latlen,altlen,timelen))
 allocate(rfile(lonlen,latlen,altlen,timelen))
 allocate(hfile(lonlen,latlen,altlen,timelen))  
 allocate(waterfile(lonlen,latlen,altlen,timelen))
 allocate(watericefile(lonlen,latlen,altlen,timelen))
-allocate(swaterfile(lonlen,latlen,timelen))
+!allocate(swaterfile(lonlen,latlen,timelen))
 allocate(swatericefile(lonlen,latlen,timelen))
+allocate(dustfile(lonlen,latlen,altlen,timelen))
+allocate(dustnfile(lonlen,latlen,altlen,timelen))
 allocate(psfile(lonlen,latlen,timelen)) 
 allocate(tsfile(lonlen,latlen,timelen))
@@ -230,7 +238,10 @@
 tfile(:,:,:,:)=0
 tsoilfile(:,:,:,:)=0 
-tfileorig(:,:,:,:)=0
-ufileorig(:,:,:,:)=0
-vfileorig(:,:,:,:)=0
+isoilfile(:,:,:,:)=0 
+dsoilfile(:,:,:,:)=0 
+dsoilvert(:)=0.
+!tfileorig(:,:,:,:)=0
+!ufileorig(:,:,:,:)=0
+!vfileorig(:,:,:,:)=0
 ufile(:,:,:,:)=0 
 vfile(:,:,:,:)=0 
@@ -239,6 +250,8 @@
 waterfile(:,:,:,:)=0
 watericefile(:,:,:,:)=0
-swaterfile(:,:,:)=0
+!swaterfile(:,:,:)=0
 swatericefile(:,:,:)=0
+dustfile(:,:,:,:)=0
+dustnfile(:,:,:,:)=0
 psfile(:,:,:)=0 
 tsfile(:,:,:)=0
@@ -366,12 +379,12 @@
    ierr = NF_INQ_VARID (nid,"emis",nvarid)
 IF (ierr .NE. NF_NOERR) THEN
-	PRINT *, '...warning: not found in diagfi !'
+        PRINT *, '...warning: not found in diagfi !'
         PRINT *, '...will be filled with a prescribed value', emiss_prescribed
-	emissfile(:,:,:)=emiss_prescribed	
+  emissfile(:,:,:)=emiss_prescribed
 ELSE
 #ifdef NC_DOUBLE
-	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, emissfile)
+  ierr = NF_GET_VAR_DOUBLE(nid, nvarid, emissfile)
 #else
-	ierr = NF_GET_VAR_REAL(nid, nvarid, emissfile)
+  ierr = NF_GET_VAR_REAL(nid, nvarid, emissfile)
 #endif
 ENDIF   
@@ -380,12 +393,12 @@
    ierr = NF_INQ_VARID (nid,"co2ice",nvarid)
 IF (ierr .NE. NF_NOERR) THEN
-	PRINT *, '...warning: not found in diagfi !'
-	PRINT *, '...will be filled with a prescribed value', co2ice_prescribed
-	co2icefile(:,:,:)=co2ice_prescribed	
+  PRINT *, '...warning: not found in diagfi !'
+  PRINT *, '...will be filled with a prescribed value', co2ice_prescribed
+  co2icefile(:,:,:)=co2ice_prescribed
 ELSE
 #ifdef NC_DOUBLE
-	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, co2icefile)
+  ierr = NF_GET_VAR_DOUBLE(nid, nvarid, co2icefile)
 #else
-	ierr = NF_GET_VAR_REAL(nid, nvarid, co2icefile)
+  ierr = NF_GET_VAR_REAL(nid, nvarid, co2icefile)
 #endif
 ENDIF
@@ -418,54 +431,54 @@
 
 
-!!"atmospheric" surface temperature is taken
-!!from original diagfi.nc first level
-!!... level is ~3-5 meters
-print *,'Near-Surface Temperature'
-   ierr = NF_INQ_VARID (nid,"temp",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      ierr = NF_INQ_VARID (nid,"t",nvarid)
-	IF (ierr .NE. NF_NOERR) THEN
-           PRINT *, "Error: Readmeteo <temp> not found"
-           stop
-	ENDIF
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfileorig)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, tfileorig)
-#endif
-   tnsfile=tfileorig(:,:,1,:)
-
-!!"atmospheric" surface u is taken
-!!from original diagfi.nc first level
-!!... level is ~3-5 meters
-print *,'Near-Surface Zonal Wind'
-   ierr = NF_INQ_VARID (nid,"u",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-     PRINT *, "Error: Readmeteo <u> not found"
-     stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufileorig)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, ufileorig)
-#endif
-   unsfile=ufileorig(:,:,1,:)
-
-!!"atmospheric" surface v is taken
-!!from original diagfi.nc first level
-!!... level is ~3-5 meters
-print *,'Near-Surface Meridional Wind'
-   ierr = NF_INQ_VARID (nid,"v",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-     PRINT *, "Error: Readmeteo <v> not found"
-     stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfileorig)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, vfileorig)
-#endif
-   vnsfile=vfileorig(:,:,1,:)
+!!!"atmospheric" surface temperature is taken
+!!!from original diagfi.nc first level
+!!!... level is ~3-5 meters
+!print *,'Near-Surface Temperature'
+!   ierr = NF_INQ_VARID (nid,"temp",nvarid)
+!   IF (ierr .NE. NF_NOERR) THEN
+!      ierr = NF_INQ_VARID (nid,"t",nvarid)
+!	IF (ierr .NE. NF_NOERR) THEN
+!           PRINT *, "Error: Readmeteo <temp> not found"
+!           stop
+!	ENDIF
+!   ENDIF
+!#ifdef NC_DOUBLE
+!   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfileorig)
+!#else
+!   ierr = NF_GET_VAR_REAL(nid, nvarid, tfileorig)
+!#endif
+!   tnsfile=tfileorig(:,:,1,:)
+
+!!!"atmospheric" surface u is taken
+!!!from original diagfi.nc first level
+!!!... level is ~3-5 meters
+!print *,'Near-Surface Zonal Wind'
+!   ierr = NF_INQ_VARID (nid,"u",nvarid)
+!   IF (ierr .NE. NF_NOERR) THEN
+!     PRINT *, "Error: Readmeteo <u> not found"
+!     stop
+!   ENDIF
+!#ifdef NC_DOUBLE
+!   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufileorig)
+!#else
+!   ierr = NF_GET_VAR_REAL(nid, nvarid, ufileorig)
+!#endif
+!   unsfile=ufileorig(:,:,1,:)
+
+!!!"atmospheric" surface v is taken
+!!!from original diagfi.nc first level
+!!!... level is ~3-5 meters
+!print *,'Near-Surface Meridional Wind'
+!   ierr = NF_INQ_VARID (nid,"v",nvarid)
+!   IF (ierr .NE. NF_NOERR) THEN
+!     PRINT *, "Error: Readmeteo <v> not found"
+!     stop
+!   ENDIF
+!#ifdef NC_DOUBLE
+!   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfileorig)
+!#else
+!   ierr = NF_GET_VAR_REAL(nid, nvarid, vfileorig)
+!#endif
+!   vnsfile=vfileorig(:,:,1,:)
 
 SELECT CASE(ident)
@@ -520,4 +533,5 @@
    ierr = NF_GET_VAR_REAL(nid, nvarid, tfile)
 #endif
+tnsfile=tfile(:,:,1,:)
 
 print *,'Zonal wind'   
@@ -532,4 +546,5 @@
    ierr = NF_GET_VAR_REAL(nid, nvarid, ufile)
 #endif
+unsfile=ufile(:,:,1,:)
 
 print *,'Meridional wind'
@@ -544,4 +559,5 @@
    ierr = NF_GET_VAR_REAL(nid, nvarid, vfile)
 #endif
+vnsfile=ufile(:,:,1,:)
 
 
@@ -554,6 +570,7 @@
       write(*,*) "...No q02 - Water vapor set to 0"
       waterfile(:,:,:,:)=0.  
-    endif
-    ierr=NF_GET_VAR_REAL(nid,nvarid,waterfile)
+    else
+      ierr=NF_GET_VAR_REAL(nid,nvarid,waterfile)
+    endif   
 
     print *,'Water ice'
@@ -562,14 +579,14 @@
       write(*,*) "...No q01 - Water ice set to 0"  
       watericefile(:,:,:,:)=0.
+    else
+      ierr=NF_GET_VAR_REAL(nid,nvarid,watericefile)
     endif
-    ierr=NF_GET_VAR_REAL(nid,nvarid,watericefile)
-
-    print *,'Surface Water vapor'
-    ierr=NF_INQ_VARID(nid,"qsurf02",nvarid)
-    if (ierr.ne.NF_NOERR) then
-      write(*,*) "...No qsurf02 - surface Water vapor set to 0"
-      swaterfile(:,:,:)=0.
-    endif
-    ierr=NF_GET_VAR_REAL(nid,nvarid,swaterfile)
+!    print *,'Surface Water vapor'
+!    ierr=NF_INQ_VARID(nid,"qsurf02",nvarid)
+!    if (ierr.ne.NF_NOERR) then
+!      write(*,*) "...No qsurf02 - surface Water vapor set to 0"
+!      swaterfile(:,:,:)=0.
+!    endif
+!    ierr=NF_GET_VAR_REAL(nid,nvarid,swaterfile)
 
     print *,'Surface Water ice'
@@ -578,8 +595,33 @@
       write(*,*) "...No qsurf01 - surface Water ice set to 0"
       swatericefile(:,:,:)=0.
+    else
+      ierr=NF_GET_VAR_REAL(nid,nvarid,swatericefile)
     endif
-    ierr=NF_GET_VAR_REAL(nid,nvarid,swatericefile)
 !!------------------------
 !! special water stuff
+!!------------------------
+
+!!------------------------
+!! special dust stuff
+!!------------------------
+    print *,'Dust mass'
+    ierr=NF_INQ_VARID(nid,"dustq",nvarid)
+    if (ierr.ne.NF_NOERR) then
+      write(*,*) "...No dustq - Dust mass set to 0"
+      dustfile(:,:,:,:)=0.
+    else
+      ierr=NF_GET_VAR_REAL(nid,nvarid,dustfile)
+    endif
+
+    print *,'Dust number'
+    ierr=NF_INQ_VARID(nid,"dustN",nvarid)
+    if (ierr.ne.NF_NOERR) then
+      write(*,*) "...No dustN - Dust number set to 0"
+      dustnfile(:,:,:,:)=0.
+    else
+      ierr=NF_GET_VAR_REAL(nid,nvarid,dustnfile)
+    endif
+!!------------------------
+!! special dust stuff
 !!------------------------
 
@@ -595,6 +637,42 @@
                 tsoilfile(:,:,l,:)=tsfile(:,:,:)
         ENDDO
+    else
+        ierr=NF_GET_VAR_REAL(nid,nvarid,tsoilfile)
     endif
-    ierr=NF_GET_VAR_REAL(nid,nvarid,tsoilfile)
+
+!!!!!!!!
+!!!!!!!! new physics (but still compatible with old physics)
+    print *,'Soil depths'
+    ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
+    if (ierr.ne.NF_NOERR) then
+        write(*,*) "...No soildepth - Set to -999"  !!! see soil_settings in LMD physics
+        DO l=1,altlen
+                dsoilvert(l)=-999.
+        ENDDO
+    else
+        ierr=NF_GET_VAR_REAL(nid,nvarid,dsoilvert)
+    endif
+    print *, 'wait a minute' !! AS: I know this could be better
+    DO m=1,lonlen
+     DO n=1,latlen
+      DO p=1,timelen
+       dsoilfile(m,n,:,p) = dsoilvert(:) 
+      ENDDO
+     ENDDO
+    ENDDO 
+    DEALLOCATE(vertdsoil)
+
+    print *,'Soil thermal inertia'
+    ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
+    if (ierr.ne.NF_NOERR) then
+        write(*,*) "...No soil therm. inert. - Set to -999"
+        DO l=1,altlen
+                isoilfile(:,:,l,:)=-999.
+        ENDDO
+    else
+        ierr=NF_GET_VAR_REAL(nid,nvarid,isoilfile)
+    endif
+!!!!!!!!
+!!!!!!!! new physics
 
 
@@ -625,5 +703,4 @@
 ENDIF
 
-
 !!-----------------------------
 !! Loop on the written files
@@ -666,4 +743,9 @@
 !! 1. Surface data
 !!---------------------------------------------
+!!
+!! a mettre pour tous sinon
+!!     WRF_DEBUG: Warning DIM             4 , NAME num_metgrid_levels REDIFINED  by
+!!            var DUSTN            26           25  in wrf_io.F90 line         2349
+!!
 
 !
@@ -701,6 +783,4 @@
 DESC='Atmospheric temperature'
 XLVL=200100.
-!SLAB=tsfile(:,:,time_out(l))
-!SLAB=tfileorig(:,:,1,time_out(l))
 SLAB=tnsfile(:,:,time_out(l))
 	! And now put everything in the destination file
@@ -721,6 +801,4 @@
 DESC='Zonal wind'
 XLVL=200100.
-!SLAB=ufile(:,:,1,time_out(l))
-!SLAB=ufileorig(:,:,1,time_out(l))
 SLAB=unsfile(:,:,time_out(l))
         ! And now put everything in the destination file
@@ -741,6 +819,4 @@
 DESC='Meridional wind'
 XLVL=200100.
-!SLAB=vfile(:,:,1,time_out(l))
-!SLAB=vfileorig(:,:,1,time_out(l))
 SLAB=vnsfile(:,:,time_out(l))
 	! And now put everything in the destination file
@@ -761,7 +837,5 @@
 DESC='Customized 2D static field'
 XLVL=200100.
-!SLAB=co2icefile(:,:,time_out(l))
 SLAB=vide(:,:)
-!SLAB=vide(:,:)+ptop
 	! And now put everything in the destination file
 	! ... Header
@@ -853,5 +927,4 @@
 DESC='ZMEA'
 XLVL=200100.
-!SLAB=vide(:,:)
 SLAB=gwparam(:,:,1)
 	! And now put everything in the destination file
@@ -872,5 +945,4 @@
 DESC='ZSTD'
 XLVL=200100.
-!SLAB=vide(:,:)
 SLAB=gwparam(:,:,2)
 	! And now put everything in the destination file
@@ -909,5 +981,4 @@
 DESC='Ground temperature'
 XLVL=200100.
-!SLAB=vide(:,:)
 SLAB=tsfile(:,:,time_out(l))
 	! And now put everything in the destination file
@@ -928,5 +999,4 @@
 DESC='ZSIG'
 XLVL=200100.
-!SLAB=vide(:,:)
 SLAB=gwparam(:,:,3)
 	! And now put everything in the destination file
@@ -947,5 +1017,4 @@
 DESC='ZGAM'
 XLVL=200100.
-!SLAB=vide(:,:)
 SLAB=gwparam(:,:,4)
 	! And now put everything in the destination file
@@ -966,5 +1035,4 @@
 DESC='ZTHE'
 XLVL=200100.
-!SLAB=vide(:,:)
 SLAB=gwparam(:,:,5)
 	! And now put everything in the destination file
@@ -982,8 +1050,8 @@
 !------------------------! 
 FIELD='SM100200'
-UNITS='fraction'
-DESC='Relative humidity'
-XLVL=200100.
-SLAB=vide(:,:)
+UNITS='kg/kg'
+DESC='Surf water ice'
+XLVL=200100.
+SLAB=swatericefile(:,:,time_out(l))
 	! And now put everything in the destination file
 	! ... Header
@@ -1004,6 +1072,5 @@
 DESC='Water vapor'
 XLVL=200100.
-!SLAB=waterfile(:,:,1,time_out(l))
-SLAB=swaterfile(:,:,time_out(l))
+SLAB=waterfile(:,:,1,time_out(l))
         ! And now put everything in the destination file
         ! ... Header
@@ -1023,6 +1090,5 @@
 DESC='Water ice'
 XLVL=200100.
-!SLAB=watericefile(:,:,1,time_out(l))
-SLAB=swatericefile(:,:,time_out(l))
+SLAB=watericefile(:,:,1,time_out(l))
         ! And now put everything in the destination file
         ! ... Header
@@ -1052,4 +1118,75 @@
 !print *,'The field '//DESC//' was written to '//output
 
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='DSOIL'
+UNITS='m'
+DESC='Soil depths'
+XLVL=200100.
+SLAB=dsoilfile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='ISOIL'
+UNITS='tiu'
+DESC='Soil thermal inertia'
+XLVL=200100.
+SLAB=isoilfile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='DUSTQ'
+UNITS='kg/kg'
+DESC='Dust mixing ratio'
+XLVL=200100.
+SLAB=dustfile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='DUSTN'
+UNITS='part/kg' 
+DESC='Dust number density'
+XLVL=200100.
+SLAB=dustnfile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
 
 
@@ -1237,4 +1374,83 @@
 !print *,'The field '//DESC//' was written to '//output
 
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='DSOIL'
+UNITS='m'
+DESC='Soil depths'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=dsoilfile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='ISOIL'
+UNITS='tiu'
+DESC='Soil thermal inertia'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=isoilfile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='DUSTQ'
+UNITS='kg/kg'
+DESC='Dust mixing ratio'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=dustfile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='DUSTN'
+UNITS='part/kg'
+DESC='Dust number density'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=dustnfile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
 
 print *,'****done file '//output, int(100.*float(l)/float(FILES)), ' % '
@@ -1252,15 +1468,20 @@
 deallocate(tfile)
 deallocate(tsoilfile)
-deallocate(tfileorig)
+deallocate(isoilfile)
+deallocate(dsoilfile)
+!deallocate(dsoilvert)
+!deallocate(tfileorig)
 deallocate(ufile)
-deallocate(ufileorig)
+!deallocate(ufileorig)
 deallocate(vfile)
-deallocate(vfileorig)
+!deallocate(vfileorig)
 deallocate(rfile)
 deallocate(hfile)  
 deallocate(waterfile)
 deallocate(watericefile)
-deallocate(swaterfile)
+!deallocate(swaterfile)
 deallocate(swatericefile)
+deallocate(dustfile)
+deallocate(dustnfile)
 deallocate(psfile) 
 deallocate(tsfile)
Index: trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/readmeteo.F90_backup
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/readmeteo.F90_backup	(revision 70)
+++ 	(revision )
@@ -1,974 +1,0 @@
-program readmeteo
-
-implicit none
-include "netcdf.inc"
-
-!------------------------------------------------------------------------!
-! Readmeteo prepares an initial state for the WPS pre-processor of WRF   !
-!                                                                        !
-! Input is a diagfi.nc NETCDF file from a LMD/GCM simulation             !
-!                                                                        !
-! Outputs are WRFSI intermediate format files ready for metgrid.exe      !
-!       http://www.mmm.ucar.edu/wrf/OnLineTutorial/WPS/IM_si.htm         !
-!                                                                        !
-! A. Spiga - 16/04/2007                                                  !
-!	     22/05/2007 : need to run zrecast as a preliminary           !
-!            07/06/2007 : external parameter file 'readmeteo.def         !
-!	     30/07/2007 : manage additional variables (tsoil, emiss,...) !	
-!            19/10/2007 : no need to run zrecast at all (eta levels)     !   
-!		12/2007 : include co2ice and emissivity			 ! 
-!									 !
-! spiga@lmd.jussieu.fr							 !	
-!------------------------------------------------------------------------!
-
-
-
-!***************************************************************************
-!***************************************************************************
-REAL, PARAMETER :: emiss_prescribed=0.95
-REAL, PARAMETER :: co2ice_prescribed=0.
-
-
-
-
-!***************************************************************************
-!***************************************************************************
-
-REAL :: ptop
-REAL, PARAMETER :: grav=3.72
-LOGICAL, PARAMETER :: blank=.false.
-CHARACTER (LEN=3), PARAMETER :: ident='LMD'
-
-
-!! Variables to be written in the output file
-!! *** NB: conformity with metgrid/src/read_met_module.F90
-CHARACTER (LEN=33) :: output
-INTEGER :: IFV
-CHARACTER (LEN=24) :: HDATE
-REAL :: XFCST
-CHARACTER (LEN=32) :: SOURCE
-CHARACTER (LEN=9) :: FIELD
-CHARACTER (LEN=25) :: UNITS
-CHARACTER (LEN=46) :: DESC
-REAL :: XLVL
-INTEGER :: NX,NY,IPROJ
-CHARACTER (LEN=8) :: STARTLOC
-REAL :: STARTLAT,STARTLON,DELTALAT,DELTALON
-REAL, POINTER, DIMENSION(:,:) :: SLAB
-
-!! Variables related to NETCDF format
-integer :: nid,nid2,nvarid,ierr,ierr2
-integer :: timedim,londim,latdim,altdim,altdim2
-integer :: rlatvdim,rlonudim
-integer :: timelen,lonlen,latlen,altlen,altlen2
-integer :: rlatvlen,rlonulen
-integer :: soillen,soildim
-integer, dimension(4) :: corner,edges
-
-!! Intermediate data arrays
-integer :: k,l
-real, dimension(:), allocatable :: lat,lon,time,alt,aps,bps,levels
-real, dimension(:,:), allocatable :: vide,ones,ghtsfile
-real, dimension(:,:,:), allocatable :: psfile,tsfile
-real, dimension(:,:,:), allocatable :: emissfile,co2icefile
-real, dimension(:,:,:), allocatable :: tnsfile,unsfile,vnsfile
-real, dimension(:,:,:,:), allocatable :: tfile,ufile,vfile,rfile,hfile
-real, dimension(:,:,:,:), allocatable :: eta_gcm
-real, dimension(:,:,:,:), allocatable :: tfileorig,ufileorig,vfileorig
-real, dimension(:,:,:,:), allocatable :: tsoilfile
-
-!! Reading the parameter file 
-integer :: tmp,increment,FILES
-integer :: tmp2,tmp3,tmp4
-character*1 :: answer
-character*4 :: str
-character*2 :: str2, str3, str4
-integer, dimension(:), allocatable :: time_out
-character*13, dimension(:), allocatable :: date_out
-character*19, dimension(:), allocatable :: date_out2
-!***************************************************************************
-!***************************************************************************
-
-
-!!---------------------
-!! Open the datafile
-!!---------------------
-ierr=NF_OPEN ("input_diagfi.nc",NF_NOWRITE,nid)
-IF (ierr.NE.NF_NOERR) THEN
-   write(*,*)'**** Please create a symbolic link called input_diagfi.nc'
-   CALL ABORT
-ENDIF
-
-
-!!--------------------------
-!! Ask for data references
-!!--------------------------
-write(*,*) "Create n files. What is n ?"
-read(*,*) FILES
-allocate(time_out(FILES))
-allocate(date_out(FILES))
-allocate(date_out2(FILES))
-
-write(*,*) ""
-write(*,*) "INPUT: Diagfi time"
-write(*,*) "list of n subscripts, separated by <Return>s"
-increment=1
-do while (increment.ne.FILES+1)
-    read(*,*) tmp
-    time_out(increment)=tmp
-    increment=increment+1
-enddo
-
-write(*,*) ""
-write(*,*) "OUTPUT: WRF time"
-write(*,*) "fill 4 lines indicating"
-write(*,*) "-year (4 digits)"
-write(*,*) "-month (2 digits)" 
-write(*,*) "-day (2 digits)"
-write(*,*) "-hour (2 digits)"
-increment=1
-do while (increment.ne.FILES+1)
-    read(*,*) str
-    read(*,*) str2
-    read(*,*) str3
-    read(*,*) str4
-    date_out(increment)=str//'-'//str2//'-'//str3//'_'//str4
-    date_out2(increment)=str//'-'//str2//'-'//str3//'_'//str4//':00:00'
-    print *,date_out(increment)
-    write(*,*) "ok? (y/n)"
-    read(*,*) answer
-    if (answer.eq.'n') then
-        write(*,*) "please write again"
-    else 
-        write(*,*) "next one, please"    
-        increment=increment+1
-    endif
-enddo
-
-
-!!-------------------
-!! Get array sizes
-!!-------------------
-ierr=NF_INQ_DIMID(nid,"Time",timedim)
-ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
-  write(*,*) "timelen: ",timelen
-ierr=NF_INQ_DIMID(nid,"latitude",latdim)
-ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
-  write(*,*) "latlen: ",latlen
-ierr=NF_INQ_DIMID(nid,"longitude",londim)
-ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
-  write(*,*) "lonlen: ",lonlen
-ierr=NF_INQ_DIMID(nid,"altitude",altdim)
-ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
-  write(*,*) "altlen: ",altlen
-
-
-
-!!-------------------------
-!! Allocate local arrays
-!!-------------------------
-allocate(eta_gcm(lonlen,latlen,altlen,timelen))
-allocate(tfile(lonlen,latlen,altlen,timelen))
-allocate(tfileorig(lonlen,latlen,altlen,timelen))
-allocate(ufile(lonlen,latlen,altlen,timelen))
-allocate(ufileorig(lonlen,latlen,altlen,timelen))
-allocate(vfile(lonlen,latlen,altlen,timelen))
-allocate(vfileorig(lonlen,latlen,altlen,timelen))
-allocate(rfile(lonlen,latlen,altlen,timelen))
-allocate(hfile(lonlen,latlen,altlen,timelen))  
-allocate(psfile(lonlen,latlen,timelen)) 
-allocate(tsfile(lonlen,latlen,timelen))
-allocate(tnsfile(lonlen,latlen,timelen))
-allocate(unsfile(lonlen,latlen,timelen))
-allocate(vnsfile(lonlen,latlen,timelen))
-allocate(emissfile(lonlen,latlen,timelen))
-allocate(co2icefile(lonlen,latlen,timelen))
-allocate(ghtsfile(lonlen,latlen))	!! no scan axis
-allocate(vide(lonlen,latlen))
-allocate(ones(lonlen,latlen))
-allocate(lat(latlen), lon(lonlen), alt(altlen), time(timelen))
-allocate(aps(altlen),bps(altlen),levels(altlen))
-
-tfile(:,:,:,:)=0 
-tfileorig(:,:,:,:)=0
-ufileorig(:,:,:,:)=0
-vfileorig(:,:,:,:)=0
-ufile(:,:,:,:)=0 
-vfile(:,:,:,:)=0 
-rfile(:,:,:,:)=0 
-hfile(:,:,:,:)=0
-psfile(:,:,:)=0 
-tsfile(:,:,:)=0
-emissfile(:,:,:)=0
-co2icefile(:,:,:)=0
-tnsfile(:,:,:)=0
-unsfile(:,:,:)=0
-vnsfile(:,:,:)=0
-ghtsfile(:,:)=0
-vide(:,:)=0
-ones(:,:)=0
-
-
-!!-------------------
-!! Read dimensions
-!!-------------------
-
-print *,'>>> Read dimensions !'
-
-print *,'Time'
-   ierr = NF_INQ_VARID (nid, "Time",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <Time> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, time)
-#endif
-   print *,time(1),' ... to ... ',time(timelen)
-
-print *,'Latitude'
-   ierr = NF_INQ_VARID (nid, "latitude",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <latitude> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lat)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, lat)
-#endif
-   print *,lat(1),' ... to ... ',lat(latlen),' ... step: ',lat(latlen)-lat(latlen-1)
-
-print *,'Longitude'
-   ierr = NF_INQ_VARID (nid, "longitude",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <longitude> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lon)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, lon)
-#endif
-   print *,lon(1),' ... to ... ',lon(lonlen),' ... step: ',lon(lonlen)-lon(lonlen-1)
-
-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)
-
-
-
-
-!!-------------------------------------------
-!! Reading 2D and 3D meteorological fields
-!!-------------------------------------------
-
-IF (blank .EQV. .false.) THEN
-
-print *,'>>> Read 2D optional fields !'
-
-print *,'Emissivity'
-   ierr = NF_INQ_VARID (nid,"emiss",nvarid)
-IF (ierr .NE. NF_NOERR) THEN
-	PRINT *, '...warning: not found in diagfi !'
-        PRINT *, '...will be filled with a prescribed value', emiss_prescribed
-	emissfile(:,:,:)=emiss_prescribed	
-ELSE
-#ifdef NC_DOUBLE
-	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, emissfile)
-#else
-	ierr = NF_GET_VAR_REAL(nid, nvarid, emissfile)
-#endif
-ENDIF   
-
-print *,'CO2 ice'
-   ierr = NF_INQ_VARID (nid,"co2ice",nvarid)
-IF (ierr .NE. NF_NOERR) THEN
-	PRINT *, '...warning: not found in diagfi !'
-	PRINT *, '...will be filled with a prescribed value', co2ice_prescribed
-	co2icefile(:,:,:)=co2ice_prescribed	
-ELSE
-#ifdef NC_DOUBLE
-	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, co2icefile)
-#else
-	ierr = NF_GET_VAR_REAL(nid, nvarid, co2icefile)
-#endif
-ENDIF
-
-print *,'>>> Read 2D surface fields !'
-
-print *,'Surface Pressure'
-   ierr = NF_INQ_VARID (nid,"ps",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <ps> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, psfile)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, psfile)
-#endif
-
-print *,'Ground Temperature'
-   ierr = NF_INQ_VARID (nid,"tsurf",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <tsurf> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsfile)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, tsfile)
-#endif
-
-
-!!"atmospheric" surface temperature is taken
-!!from original diagfi.nc first level
-!!... level is ~3-5 meters
-print *,'Near-Surface Temperature'
-   ierr = NF_INQ_VARID (nid,"temp",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <temp> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfileorig)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, tfileorig)
-#endif
-   tnsfile=tfileorig(:,:,1,:)
-
-!!"atmospheric" surface u is taken
-!!from original diagfi.nc first level
-!!... level is ~3-5 meters
-print *,'Near-Surface Zonal Wind'
-   ierr = NF_INQ_VARID (nid,"u",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-     PRINT *, "Error: Readmeteo <u> not found"
-     stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufileorig)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, ufileorig)
-#endif
-   unsfile=ufileorig(:,:,1,:)
-
-!!"atmospheric" surface v is taken
-!!from original diagfi.nc first level
-!!... level is ~3-5 meters
-print *,'Near-Surface Meridional Wind'
-   ierr = NF_INQ_VARID (nid,"v",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-     PRINT *, "Error: Readmeteo <v> not found"
-     stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfileorig)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, vfileorig)
-#endif
-   vnsfile=vfileorig(:,:,1,:)
-
-
-   print *,'Geopotential height at the ground'
-   ierr = NF_INQ_VARID (nid,"phisinit",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <phisinit> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ghtsfile)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, ghtsfile)
-#endif
-!**** from geopotential to geopotential height
-ghtsfile=ghtsfile/grav
-!****
-
-print *,'>>> Read 3D meteorological fields ! - This may take some time ...'
-
-print *,'Temperature'
-   ierr = NF_INQ_VARID (nid,"temp",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <t> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfile)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, tfile)
-#endif
-
-print *,'Zonal wind'   
-   ierr = NF_INQ_VARID (nid,"u",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <u> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufile)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, ufile)
-#endif
-
-print *,'Meridional wind'
-   ierr = NF_INQ_VARID (nid,"v",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <v> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfile)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, vfile)
-#endif
-
-ENDIF
-
-
-!!-----------------------------
-!! Loop on the written files
-!! >>> what follows is pretty repetitive ... 
-!!        gotta do something (one more loop?)
-!!-----------------------------
-
-!!! Equivalent eta levels for WRF interpolation in initialize_real
-!print *,'Computing equivalent eta levels'
-!	
-!	!ptop will be passed through RH(surface)
-!	ptop=MINVAL(aps(altlen)+bps(altlen)*psfile(:,:,:))
-!	print *, 'ptop', ptop
-!
-!print *, 'sample: equivalent eta levels at i=1,j=1'
-!DO k = 1,altlen
-!        levels(k)=-k
-!        eta_gcm(:,:,k,:)=(aps(k)+bps(k)*psfile(:,:,:)-ptop)/(psfile(:,:,:)-ptop)
-!        print *, eta_gcm(1,1,k,1) 
-!END DO
-
-!!---better to pass pressure values 
-!!---the eta treatment is done in module_initialize_real
-DO k = 1,altlen
-        levels(k)=-k
-        !!dummy decreasing levels
-END DO
-
-
-
-
-!! Dummy surface level is XLVL=200100.
-
-
-DO l=1,FILES
-
-
-!!---------------------------------------------
-!! Write file in intermediate format for WPS
-!! 1. Surface data
-!!---------------------------------------------
-
-!
-! These variables remain the same in the "loop"
-!
-HDATE=date_out(l)
-IFV=4
-XFCST=0.
-SOURCE='LMD GCM'
-NX=lonlen
-NY=latlen
-ALLOCATE(SLAB(NX,NY))
-IPROJ=0
-STARTLOC='SWCORNER'
-STARTLAT=lat(1) 
-STARTLON=lon(1)
-DELTALAT=lat(latlen)-lat(latlen-1)
-DELTALON=lon(lonlen)-lon(lonlen-1)
-!
-! Open the data destination file
-!
-output=ident//"/"//ident//":"//date_out2(l)       
-open(UNIT=1,			&
-	FILE=output,		&
-	STATUS='REPLACE',	&
-	FORM='UNFORMATTED',	&
-	ACCESS='SEQUENTIAL')
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='T'
-UNITS='K'
-DESC='Atmospheric temperature'
-XLVL=200100.
-!SLAB=tsfile(:,:,time_out(l))
-!SLAB=tfileorig(:,:,1,time_out(l))
-SLAB=tnsfile(:,:,time_out(l))
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='U'
-UNITS='m s{-1}'
-DESC='Zonal wind'
-XLVL=200100.
-!SLAB=ufile(:,:,1,time_out(l))
-!SLAB=ufileorig(:,:,1,time_out(l))
-SLAB=unsfile(:,:,time_out(l))
-        ! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='V'
-UNITS='m s{-1}'
-DESC='Meridional wind'
-XLVL=200100.
-!SLAB=vfile(:,:,1,time_out(l))
-!SLAB=vfileorig(:,:,1,time_out(l))
-SLAB=vnsfile(:,:,time_out(l))
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='RH'
-UNITS=''
-DESC='Customized 2D static field'
-XLVL=200100.
-!SLAB=co2icefile(:,:,time_out(l))
-SLAB=vide(:,:)
-!SLAB=vide(:,:)+ptop
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='SOILHGT'
-UNITS='m' 
-DESC='Terrain field of source analysis'  !Ground geopotential height
-XLVL=200100.
-SLAB=ghtsfile(:,:)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='PSFC'
-UNITS='Pa'
-DESC='Surface Pressure'
-XLVL=200100.
-SLAB=psfile(:,:,time_out(l))
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='ST000010'
-UNITS=''
-DESC='Emissivity'
-XLVL=200100.
-!SLAB=vide(:,:)
-SLAB=emissfile(:,:,time_out(l))
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='ST010040'
-UNITS=''
-DESC='CO2 ice'
-XLVL=200100.
-!SLAB=vide(:,:)
-SLAB=co2icefile(:,:,time_out(l))
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='ST040100'
-UNITS='K'
-DESC='T of 40-100 cm ground layer'
-XLVL=200100.
-SLAB=vide(:,:)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='ST100200'
-UNITS='K'
-DESC='T of 100-200 cm ground layer'
-XLVL=200100.
-SLAB=vide(:,:)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='LANDSEA'
-UNITS='0/1 Flag'
-DESC='Land/Sea flag'
-XLVL=200100.
-SLAB=ones(:,:)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='SKINTEMP'
-UNITS='K'
-DESC='Ground temperature'
-XLVL=200100.
-!SLAB=vide(:,:)
-SLAB=tsfile(:,:,time_out(l))
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='SM000010'
-UNITS='fraction'
-DESC='Relative humidity'
-XLVL=200100.
-SLAB=vide(:,:)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='SM010040'
-UNITS='fraction'
-DESC='Relative humidity'
-XLVL=200100.
-SLAB=vide(:,:)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='SM040100'
-UNITS='fraction'
-DESC='Relative humidity'
-XLVL=200100.
-SLAB=vide(:,:)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='SM100200'
-UNITS='fraction'
-DESC='Relative humidity'
-XLVL=200100.
-SLAB=vide(:,:)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
-	! ... Data
-	write(1) SLAB
-print *,'The field '//DESC//' was written to '//output
-
-
-
-!!----------------------------------------------------
-!! Write file in intermediate format for WPS
-!! 2. 3D meteorological data
-!! >>> same stuff, but handle vertical profile
-!! >>> !! XLVL must be decreasing (pressure levels) 
-!!----------------------------------------------------
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='T'
-UNITS='K'
-DESC='Atmospheric temperature'
-DO k = 1,altlen
-	XLVL=levels(k)
-	SLAB=tfile(:,:,k,time_out(l))
-		! And now put everything in the destination file
-		! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-		! ... Data
-		write(1) SLAB
-END DO
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='U'
-UNITS='m s{-1}'
-DESC='Zonal wind'
-DO k = 1,altlen
-        XLVL=levels(k)
-	SLAB=ufile(:,:,k,time_out(l))
-		! And now put everything in the destination file
-		! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-		! ... Data
-		write(1) SLAB
-END DO
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='V'
-UNITS='m s{-1}'
-DESC='Meridional wind'
-DO k = 1,altlen
-        XLVL=levels(k)
-	SLAB=vfile(:,:,k,time_out(l))
-		! And now put everything in the destination file
-		! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-		! ... Data
-		write(1) SLAB
-END DO
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='RH'
-UNITS=''
-DESC='Customized 2D static field'
-DESC='Eta-levels from the GCM'
-DESC='Pressure values from the GCM'
-DO k = 1,altlen
-        XLVL=levels(k)
-	!SLAB=rfile(:,:,k,time_out(l))
-        !SLAB=eta_gcm(:,:,k,time_out(l))
-        SLAB=aps(k)+bps(k)*psfile(:,:,time_out(l))
-		! And now put everything in the destination file
-		! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-		! ... Data
-		write(1) SLAB
-END DO
-print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='HGT'
-UNITS='m'
-DESC='Height'
-DO k = 1,altlen
-        XLVL=levels(k)
-!!*******
-!! PROVISOIRE: normalement, il faudrait la hauteur geopotentielle
-!!*******
-!!SLAB=hfile(:,:,k,time_out(l))
-!!SLAB=vide(:,:)
-        SLAB=10000.*alog(610./(aps(k)+bps(k)*psfile(:,:,time_out(l))))
-!!however, not used by initialize_real
-		! And now put everything in the destination file
-		! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-		! ... Data
-		write(1) SLAB
-END DO
-print *,'The field '//DESC//' was written to '//output
-
-print *,'****Close file '//output
-close(1)
-
-DEALLOCATE(SLAB)
-
-END DO !! end of the files loop
-
-
-!!-------------------------
-!! DeAllocate local arrays
-!!-------------------------
-deallocate(eta_gcm)
-deallocate(tfile)
-deallocate(tfileorig)
-deallocate(ufile)
-deallocate(ufileorig)
-deallocate(vfile)
-deallocate(vfileorig)
-deallocate(rfile)
-deallocate(hfile)  
-deallocate(psfile) 
-deallocate(tsfile)
-deallocate(tnsfile)
-deallocate(unsfile)
-deallocate(vnsfile)
-deallocate(emissfile)
-deallocate(co2icefile)
-deallocate(ghtsfile)	!! no scan axis
-deallocate(vide)
-deallocate(ones)
-deallocate(lat, lon, alt, time)
-deallocate(aps,bps,levels)
-
-
-print *, 'End of pre-WPS process !'
-
-
-end
Index: trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/readmeteo_newphys.F90
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/readmeteo_newphys.F90	(revision 70)
+++ 	(revision )
@@ -1,1491 +1,0 @@
-program readmeteo
-
-implicit none
-include "netcdf.inc"
-
-!------------------------------------------------------------------------!
-! Readmeteo prepares an initial state for the WPS pre-processor of WRF   !
-!                                                                        !
-! Input is a diagfi.nc NETCDF file from a LMD/GCM simulation             !
-!                                                                        !
-! Outputs are WRFSI intermediate format files ready for metgrid.exe      !
-!       http://www.mmm.ucar.edu/wrf/OnLineTutorial/WPS/IM_si.htm         !
-!                                                                        !
-! **** Please create a WPSFEED folder (or a symbolic link) ****          !       
-!                                                                        !
-! A. Spiga - 16/04/2007                                                  !
-!	     22/05/2007 : need to run zrecast as a preliminary           !
-!            07/06/2007 : external parameter file 'readmeteo.def         !
-!	     30/07/2007 : manage additional variables (tsoil, emiss,...) !	
-!            19/10/2007 : no need to run zrecast at all (eta levels)     !   
-!		12/2007 : include co2ice and emissivity			 ! 
-!               02/2008 : include water vapor and ice                    !
-!               01/2010 : possible use of diagsoil for new physics       ! 
-!               02/2011 : add dust tracers                               !
-!									 !
-! spiga@lmd.jussieu.fr							 !	
-!------------------------------------------------------------------------!
-
-
-
-!***************************************************************************
-!***************************************************************************
-REAL, PARAMETER :: emiss_prescribed=0.95
-REAL, PARAMETER :: co2ice_prescribed=0.
-CHARACTER (LEN=3), PARAMETER :: ident='LMD'
-
-
-
-!***************************************************************************
-!***************************************************************************
-
-REAL :: ptop
-REAL, PARAMETER :: grav=3.72
-LOGICAL, PARAMETER :: blank=.false.
-
-
-
-!! Variables to be written in the output file
-!! *** NB: conformity with metgrid/src/read_met_module.F90
-CHARACTER (LEN=33) :: output
-INTEGER :: IFV
-CHARACTER (LEN=24) :: HDATE
-REAL :: XFCST
-CHARACTER (LEN=32) :: SOURCE
-CHARACTER (LEN=9) :: FIELD
-CHARACTER (LEN=25) :: UNITS
-CHARACTER (LEN=46) :: DESC
-REAL :: XLVL
-INTEGER :: NX,NY,IPROJ
-CHARACTER (LEN=8) :: STARTLOC
-REAL :: STARTLAT,STARTLON,DELTALAT,DELTALON
-REAL, POINTER, DIMENSION(:,:) :: SLAB
-
-!! Variables related to NETCDF format
-integer :: nid,nid2,nid3,nvarid,ierr,ierr2
-integer :: timedim,londim,latdim,altdim,altdim2
-integer :: rlatvdim,rlonudim
-integer :: timelen,lonlen,latlen,altlen,altlen2
-integer :: rlatvlen,rlonulen
-integer :: soillen,soildim
-integer :: nphys
-integer, dimension(4) :: corner,edges
-
-!! Intermediate data arrays
-integer :: k,l,m,n,p
-real, dimension(:), allocatable :: lat,lon,time,alt,aps,bps,levels
-real, dimension(:,:), allocatable :: vide,ones,ghtsfile
-real, dimension(:,:), allocatable :: interm
-real, dimension(:,:,:), allocatable :: gwparam
-real, dimension(:,:,:), allocatable :: psfile,tsfile
-real, dimension(:,:,:), allocatable :: emissfile,co2icefile
-real, dimension(:,:,:), allocatable :: tnsfile,unsfile,vnsfile
-real, dimension(:,:,:,:), allocatable :: tfile,ufile,vfile,rfile,hfile
-real, dimension(:,:,:,:), allocatable :: eta_gcm
-real, dimension(:,:,:,:), allocatable :: tfileorig,ufileorig,vfileorig
-real, dimension(:,:,:,:), allocatable :: tsoilfile, dsoilfile, isoilfile
-real, dimension(:,:,:,:), allocatable :: waterfile, watericefile
-real, dimension(:,:,:,:), allocatable :: dustfile,dustnfile
-
-!! Reading the parameter file 
-integer :: tmp,increment,FILES
-integer :: tmp2,tmp3,tmp4
-character*1 :: answer
-character*4 :: str
-character*2 :: str2, str3, str4
-integer, dimension(:), allocatable :: time_out
-character*13, dimension(:), allocatable :: date_out
-character*19, dimension(:), allocatable :: date_out2
-!***************************************************************************
-!***************************************************************************
-
-
-!!---------------------
-!! Open the datafile
-!!---------------------
-ierr=NF_OPEN ("input_diagfi.nc",NF_NOWRITE,nid)
-IF (ierr.NE.NF_NOERR) THEN
-   write(*,*)'**** Please create a symbolic link called input_diagfi.nc'
-   CALL ABORT
-ENDIF
-
-
-!!--------------------------
-!! Ask for data references
-!!--------------------------
-!write(*,*) "Create n files. What is n ?"
-read(*,*) FILES
-allocate(time_out(FILES))
-allocate(date_out(FILES))
-allocate(date_out2(FILES))
-
-!write(*,*) ""
-!write(*,*) "INPUT: Diagfi time"
-!write(*,*) "list of n subscripts, separated by <Return>s"
-increment=1
-do while (increment.ne.FILES+1)
-    read(*,*) tmp
-    time_out(increment)=tmp
-    increment=increment+1
-enddo
-
-!write(*,*) ""
-!write(*,*) "OUTPUT: WRF time"
-!write(*,*) "fill 4 lines indicating"
-!write(*,*) "-year (4 digits)"
-!write(*,*) "-month (2 digits)" 
-!write(*,*) "-day (2 digits)"
-!write(*,*) "-hour (2 digits)"
-increment=1
-do while (increment.ne.FILES+1)
-    read(*,*) str
-    read(*,*) str2
-    read(*,*) str3
-    read(*,*) str4
-    date_out(increment)=str//'-'//str2//'-'//str3//'_'//str4
-    date_out2(increment)=str//'-'//str2//'-'//str3//'_'//str4//':00:00'
-    !print *,date_out(increment)
-    !write(*,*) "ok? (y/n)"
-    read(*,*) answer
-    if (answer.eq.'n') then
-        !write(*,*) "please write again"
-    else 
-        !write(*,*) "next one, please"    
-        increment=increment+1
-    endif
-enddo
-
-
-!!-------------------
-!! Get array sizes
-!!-------------------
-SELECT CASE(ident)
-CASE('LMD')
-ierr=NF_INQ_DIMID(nid,"Time",timedim)
-CASE('OXF')
-ierr=NF_INQ_DIMID(nid,"time",timedim)
-END SELECT
-ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
-  write(*,*) "timelen: ",timelen
-
-SELECT CASE(ident)
-CASE('LMD')
-ierr=NF_INQ_DIMID(nid,"latitude",latdim)
-CASE('OXF')
-ierr=NF_INQ_DIMID(nid,"lat",latdim)
-END SELECT
-ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
-  write(*,*) "latlen: ",latlen
-
-SELECT CASE(ident)
-CASE('LMD')
-ierr=NF_INQ_DIMID(nid,"longitude",londim)
-CASE('OXF')
-ierr=NF_INQ_DIMID(nid,"lon",londim)
-END SELECT
-ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
-  write(*,*) "lonlen: ",lonlen
-
-SELECT CASE(ident)
-CASE('LMD')
-ierr=NF_INQ_DIMID(nid,"altitude",altdim)
-CASE('OXF')
-ierr=NF_INQ_DIMID(nid,"sigma",altdim)
-END SELECT
-ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
-  write(*,*) "altlen: ",altlen
-
-
-
-!!-------------------------
-!! Allocate local arrays
-!!-------------------------
-allocate(eta_gcm(lonlen,latlen,altlen,timelen))
-allocate(tfile(lonlen,latlen,altlen,timelen))
-allocate(tsoilfile(lonlen,latlen,altlen,timelen))
-allocate(dsoilfile(lonlen,latlen,altlen,timelen))
-allocate(isoilfile(lonlen,latlen,altlen,timelen))
-allocate(tfileorig(lonlen,latlen,altlen,timelen))
-allocate(ufile(lonlen,latlen,altlen,timelen))
-allocate(ufileorig(lonlen,latlen,altlen,timelen))
-allocate(vfile(lonlen,latlen,altlen,timelen))
-allocate(vfileorig(lonlen,latlen,altlen,timelen))
-allocate(rfile(lonlen,latlen,altlen,timelen))
-allocate(hfile(lonlen,latlen,altlen,timelen))  
-allocate(waterfile(lonlen,latlen,altlen,timelen))
-allocate(watericefile(lonlen,latlen,altlen,timelen))
-allocate(dustfile(lonlen,latlen,altlen,timelen))
-allocate(dustnfile(lonlen,latlen,altlen,timelen))
-allocate(psfile(lonlen,latlen,timelen)) 
-allocate(tsfile(lonlen,latlen,timelen))
-allocate(tnsfile(lonlen,latlen,timelen))
-allocate(unsfile(lonlen,latlen,timelen))
-allocate(vnsfile(lonlen,latlen,timelen))
-allocate(emissfile(lonlen,latlen,timelen))
-allocate(co2icefile(lonlen,latlen,timelen))
-allocate(interm(lonlen,latlen))
-allocate(gwparam(lonlen,latlen,5))
-allocate(ghtsfile(lonlen,latlen))  !! no scan axis
-allocate(vide(lonlen,latlen))
-allocate(ones(lonlen,latlen))
-allocate(lat(latlen), lon(lonlen), alt(altlen), time(timelen))
-allocate(aps(altlen),bps(altlen),levels(altlen))
-
-tfile(:,:,:,:)=0
-tsoilfile(:,:,:,:)=0 
-isoilfile(:,:,:,:)=0 
-dsoilfile(:,:,:,:)=0 
-tfileorig(:,:,:,:)=0
-ufileorig(:,:,:,:)=0
-vfileorig(:,:,:,:)=0
-ufile(:,:,:,:)=0 
-vfile(:,:,:,:)=0 
-rfile(:,:,:,:)=0 
-hfile(:,:,:,:)=0
-waterfile(:,:,:,:)=0
-watericefile(:,:,:,:)=0
-dustfile(:,:,:,:)=0
-dustnfile(:,:,:,:)=0
-psfile(:,:,:)=0 
-tsfile(:,:,:)=0
-emissfile(:,:,:)=0
-co2icefile(:,:,:)=0
-tnsfile(:,:,:)=0
-unsfile(:,:,:)=0
-vnsfile(:,:,:)=0
-interm(:,:)=0
-gwparam(:,:,:)=0
-ghtsfile(:,:)=0
-vide(:,:)=0
-ones(:,:)=0
-
-
-!!-------------------
-!! Read dimensions
-!!-------------------
-
-print *,'>>> Read dimensions !'
-
-print *,'Time'
-SELECT CASE(ident)
-CASE('LMD')
-   ierr = NF_INQ_VARID (nid, "Time",nvarid)
-CASE('OXF')
-   ierr = NF_INQ_VARID (nid, "time",nvarid)
-END SELECT
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <Time> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, time)
-#endif
-   print *,time(1),' ... to ... ',time(timelen)
-
-print *,'Latitude'
-SELECT CASE(ident)
-CASE('LMD')
-   ierr = NF_INQ_VARID (nid, "latitude",nvarid)
-CASE('OXF')
-   ierr = NF_INQ_VARID (nid, "lat",nvarid)
-END SELECT
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <latitude> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lat)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, lat)
-#endif
-   print *,lat(1),' ... to ... ',lat(latlen),' ... step: ',lat(latlen)-lat(latlen-1)
-
-print *,'Longitude'
-SELECT CASE(ident)
-CASE('LMD')
-   ierr = NF_INQ_VARID (nid, "longitude",nvarid)
-CASE('OXF')
-   ierr = NF_INQ_VARID (nid, "lon",nvarid)
-END SELECT
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <longitude> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lon)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, lon)
-#endif
-   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
-
-
-!!-------------------------------------------
-!! Reading 2D and 3D meteorological fields
-!!-------------------------------------------
-
-IF (blank .EQV. .false.) THEN
-
-print *,'>>> Read 2D optional fields !'
-
-print *,'Emissivity'
-   ierr = NF_INQ_VARID (nid,"emis",nvarid)
-IF (ierr .NE. NF_NOERR) THEN
-        PRINT *, '...warning: not found in diagfi !'
-        PRINT *, '...will be filled with a prescribed value', emiss_prescribed
-        emissfile(:,:,:)=emiss_prescribed
-ELSE
-#ifdef NC_DOUBLE
-        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, emissfile)
-#else
-        ierr = NF_GET_VAR_REAL(nid, nvarid, emissfile)
-#endif
-ENDIF   
-
-print *,'CO2 ice'
-   ierr = NF_INQ_VARID (nid,"co2ice",nvarid)
-IF (ierr .NE. NF_NOERR) THEN
-        PRINT *, '...warning: not found in diagfi !'
-        PRINT *, '...will be filled with a prescribed value', co2ice_prescribed
-        co2icefile(:,:,:)=co2ice_prescribed
-ELSE
-#ifdef NC_DOUBLE
-        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, co2icefile)
-#else
-        ierr = NF_GET_VAR_REAL(nid, nvarid, co2icefile)
-#endif
-ENDIF
-
-print *,'>>> Read 2D surface fields !'
-
-print *,'Surface Pressure'
-   ierr = NF_INQ_VARID (nid,"ps",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <ps> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, psfile)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, psfile)
-#endif
-
-print *,'Ground Temperature'
-   ierr = NF_INQ_VARID (nid,"tsurf",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <tsurf> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsfile)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, tsfile)
-#endif
-
-
-!!"atmospheric" surface temperature is taken
-!!from original diagfi.nc first level
-!!... level is ~3-5 meters
-print *,'Near-Surface Temperature'
-   ierr = NF_INQ_VARID (nid,"temp",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      ierr = NF_INQ_VARID (nid,"t",nvarid)
-        IF (ierr .NE. NF_NOERR) THEN
-           PRINT *, "Error: Readmeteo <temp> not found"
-           stop
-        ENDIF
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfileorig)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, tfileorig)
-#endif
-   tnsfile=tfileorig(:,:,1,:)
-
-!!"atmospheric" surface u is taken
-!!from original diagfi.nc first level
-!!... level is ~3-5 meters
-print *,'Near-Surface Zonal Wind'
-   ierr = NF_INQ_VARID (nid,"u",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-     PRINT *, "Error: Readmeteo <u> not found"
-     stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufileorig)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, ufileorig)
-#endif
-   unsfile=ufileorig(:,:,1,:)
-
-!!"atmospheric" surface v is taken
-!!from original diagfi.nc first level
-!!... level is ~3-5 meters
-print *,'Near-Surface Meridional Wind'
-   ierr = NF_INQ_VARID (nid,"v",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-     PRINT *, "Error: Readmeteo <v> not found"
-     stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfileorig)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, vfileorig)
-#endif
-   vnsfile=vfileorig(:,:,1,:)
-
-SELECT CASE(ident)
-CASE('LMD')
-   print *,'Geopotential height at the ground'
-   ierr = NF_INQ_VARID (nid,"phisinit",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <phisinit> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ghtsfile)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, ghtsfile)
-#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 ...'
-
-print *,'Temperature'
-   ierr = NF_INQ_VARID (nid,"temp",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      ierr = NF_INQ_VARID (nid,"t",nvarid)
-        IF (ierr .NE. NF_NOERR) THEN
-          PRINT *, "Error: Readmeteo <t> not found"
-          stop
-        ENDIF
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfile)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, tfile)
-#endif
-
-print *,'Zonal wind'   
-   ierr = NF_INQ_VARID (nid,"u",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <u> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufile)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, ufile)
-#endif
-
-print *,'Meridional wind'
-   ierr = NF_INQ_VARID (nid,"v",nvarid)
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT *, "Error: Readmeteo <v> not found"
-      stop
-   ENDIF
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfile)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, vfile)
-#endif
-
-
-!!------------------------
-!! special water stuff
-!!------------------------
-    print *,'Water vapor'
-    ierr=NF_INQ_VARID(nid,"q02",nvarid)
-    if (ierr.ne.NF_NOERR) then
-      write(*,*) "...No q02 - Water vapor set to 0"
-      waterfile(:,:,:,:)=0.  
-    endif
-    ierr=NF_GET_VAR_REAL(nid,nvarid,waterfile)
-
-    print *,'Water ice'
-    ierr=NF_INQ_VARID(nid,"q01",nvarid)
-    if (ierr.ne.NF_NOERR) then
-      write(*,*) "...No q01 - Water ice set to 0"  
-      watericefile(:,:,:,:)=0.
-    endif
-    ierr=NF_GET_VAR_REAL(nid,nvarid,watericefile)
-!!------------------------
-!! special water stuff
-!!------------------------
-
-!!------------------------
-!! special dust stuff
-!!------------------------
-    print *,'Dust mass'
-    ierr=NF_INQ_VARID(nid,"dustq",nvarid)
-    if (ierr.ne.NF_NOERR) then
-      write(*,*) "...No dustq - Dust mass set to 0"
-      dustfile(:,:,:,:)=0.
-    endif
-    ierr=NF_GET_VAR_REAL(nid,nvarid,dustfile)
-
-    print *,'Dust number'
-    ierr=NF_INQ_VARID(nid,"dustN",nvarid)
-    if (ierr.ne.NF_NOERR) then
-      write(*,*) "...No dustN - Dust number set to 0"
-      dustnfile(:,:,:,:)=0.
-    endif
-    ierr=NF_GET_VAR_REAL(nid,nvarid,dustnfile)
-!!------------------------
-!! special dust stuff
-!!------------------------
-
-
-!SELECT CASE(ident)
-!CASE('LMD')
-
-    print *,'Soil temperatures'
-    ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
-    if (ierr.ne.NF_NOERR) then
-        write(*,*) "...No tsoil - Soil temperatures set isothermal with surface temperature"
-        DO l=1,altlen
-                tsoilfile(:,:,l,:)=tsfile(:,:,:)
-        ENDDO
-    endif
-    ierr=NF_GET_VAR_REAL(nid,nvarid,tsoilfile)
-  
-!!!!!!!!
-!!!!!!!! new physics (but still compatible with old physics)
-
-    print *,'Soil depths'
-    ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
-    if (ierr.ne.NF_NOERR) then
-        write(*,*) "...No soildepth - Set to default"  !!! see soil_settings in LMD physics
-        DO l=1,altlen
-                dsoilfile(1,1,l,:)=-999.
-        ENDDO
-    endif
-    ierr=NF_GET_VAR_REAL(nid,nvarid,dsoilfile(1,1,:,1))
-    print *, dsoilfile(1,1,:,1)
-    DO m=1,lonlen
-    DO n=1,latlen
-    DO p=1,timelen
-      dsoilfile(m,n,:,p) = dsoilfile(1,1,:,1) 
-    ENDDO
-    ENDDO
-    ENDDO
-
-    print *,'Soil thermal inertia'
-    ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
-    if (ierr.ne.NF_NOERR) then
-        write(*,*) "...No soil therm. inert. - Set to -999"
-        DO l=1,altlen
-                isoilfile(:,:,l,:)=-999.
-        ENDDO
-    endif
-    ierr=NF_GET_VAR_REAL(nid,nvarid,isoilfile)
-
-!!!!!!!!
-!!!!!!!! new physics
-
-
-!END SELECT
-
-ierr=NF_CLOSE(nid)
-
-!!!!lott stuff
-!print *,'get lott'
-!ierr=NF_OPEN("init_lott.nc",NF_NOWRITE,nid3)
-!ierr=NF_INQ_VARID(nid3,"ZMEA",nvarid)
-!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
-!gwparam(:,:,1)=interm(:,:)
-!ierr=NF_INQ_VARID(nid3,"ZSTD",nvarid)
-!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
-!gwparam(:,:,2)=interm(:,:)
-!ierr=NF_INQ_VARID(nid3,"ZSIG",nvarid)
-!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
-!gwparam(:,:,3)=interm(:,:)
-!ierr=NF_INQ_VARID(nid3,"ZGAM",nvarid)
-!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
-!gwparam(:,:,4)=interm(:,:)
-!ierr=NF_INQ_VARID(nid3,"ZTHE",nvarid)
-!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
-!gwparam(:,:,5)=interm(:,:)
-!ierr=NF_CLOSE(nid3)
-
-ENDIF
-
-
-!!-----------------------------
-!! Loop on the written files
-!! >>> what follows is pretty repetitive ... 
-!!        gotta do something (one more loop?)
-!!-----------------------------
-
-!!! Equivalent eta levels for WRF interpolation in initialize_real
-!print *,'Computing equivalent eta levels'
-!	
-!	!ptop will be passed through RH(surface)
-!	ptop=MINVAL(aps(altlen)+bps(altlen)*psfile(:,:,:))
-!	print *, 'ptop', ptop
-!
-!print *, 'sample: equivalent eta levels at i=1,j=1'
-!DO k = 1,altlen
-!        levels(k)=-k
-!        eta_gcm(:,:,k,:)=(aps(k)+bps(k)*psfile(:,:,:)-ptop)/(psfile(:,:,:)-ptop)
-!        print *, eta_gcm(1,1,k,1) 
-!END DO
-
-!!---better to pass pressure values 
-!!---the eta treatment is done in module_initialize_real
-DO k = 1,altlen
-        levels(k)=-k
-        !!dummy decreasing levels
-END DO
-
-
-
-
-!! Dummy surface level is XLVL=200100.
-
-
-DO l=1,FILES
-
-
-!!---------------------------------------------
-!! Write file in intermediate format for WPS
-!! 1. Surface data
-!!---------------------------------------------
-!!
-!! a mettre pour tous sinon
-!!     WRF_DEBUG: Warning DIM             4 , NAME num_metgrid_levels REDIFINED  by
-!!            var DUSTN            26           25  in wrf_io.F90 line         2349
-!!
-
-!
-! These variables remain the same in the "loop"
-!
-HDATE=date_out(l)
-IFV=4
-XFCST=0.
-SOURCE=ident
-NX=lonlen
-NY=latlen
-ALLOCATE(SLAB(NX,NY))
-IPROJ=0
-STARTLOC='SWCORNER'
-STARTLAT=lat(1) 
-STARTLON=lon(1)
-DELTALAT=lat(latlen)-lat(latlen-1)
-DELTALON=lon(lonlen)-lon(lonlen-1)
-!
-! Open the data destination file
-!
-output="./WPSFEED/"//ident//":"//date_out2(l)       
-open(UNIT=1,			&
-	FILE=output,		&
-	STATUS='REPLACE',	&
-	FORM='UNFORMATTED',	&
-	ACCESS='SEQUENTIAL')
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='T'
-UNITS='K'
-DESC='Atmospheric temperature'
-XLVL=200100.
-!SLAB=tsfile(:,:,time_out(l))
-!SLAB=tfileorig(:,:,1,time_out(l))
-SLAB=tnsfile(:,:,time_out(l))
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='U'
-UNITS='m s{-1}'
-DESC='Zonal wind'
-XLVL=200100.
-!SLAB=ufile(:,:,1,time_out(l))
-!SLAB=ufileorig(:,:,1,time_out(l))
-SLAB=unsfile(:,:,time_out(l))
-        ! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='V'
-UNITS='m s{-1}'
-DESC='Meridional wind'
-XLVL=200100.
-!SLAB=vfile(:,:,1,time_out(l))
-!SLAB=vfileorig(:,:,1,time_out(l))
-SLAB=vnsfile(:,:,time_out(l))
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='RH'
-UNITS=''
-DESC='Customized 2D static field'
-XLVL=200100.
-!SLAB=co2icefile(:,:,time_out(l))
-SLAB=vide(:,:)
-!SLAB=vide(:,:)+ptop
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='SOILHGT'
-UNITS='m' 
-DESC='Terrain field of source analysis'  !Ground geopotential height
-XLVL=200100.
-SLAB=ghtsfile(:,:)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='PSFC'
-UNITS='Pa'
-DESC='Surface Pressure'
-XLVL=200100.
-SLAB=psfile(:,:,time_out(l))
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='ST000010'
-UNITS=''
-DESC='Emissivity'
-XLVL=200100.
-SLAB=emissfile(:,:,time_out(l))
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='ST010040'
-UNITS=''
-DESC='CO2 ice'
-XLVL=200100.
-SLAB=co2icefile(:,:,time_out(l))
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='ST040100'
-UNITS=''
-DESC='ZMEA'
-XLVL=200100.
-!SLAB=vide(:,:)
-SLAB=gwparam(:,:,1)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='ST100200'
-UNITS=''
-DESC='ZSTD'
-XLVL=200100.
-!SLAB=vide(:,:)
-SLAB=gwparam(:,:,2)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='LANDSEA'
-UNITS='0/1 Flag'
-DESC='Land/Sea flag'
-XLVL=200100.
-SLAB=ones(:,:)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='SKINTEMP'
-UNITS='K'
-DESC='Ground temperature'
-XLVL=200100.
-!SLAB=vide(:,:)
-SLAB=tsfile(:,:,time_out(l))
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='SM000010'
-UNITS=''
-DESC='ZSIG'
-XLVL=200100.
-!SLAB=vide(:,:)
-SLAB=gwparam(:,:,3)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='SM010040'
-UNITS=''
-DESC='ZGAM'
-XLVL=200100.
-!SLAB=vide(:,:)
-SLAB=gwparam(:,:,4)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='SM040100'
-UNITS=''
-DESC='ZTHE'
-XLVL=200100.
-!SLAB=vide(:,:)
-SLAB=gwparam(:,:,5)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='SM100200'
-UNITS='fraction'
-DESC='Relative humidity'
-XLVL=200100.
-SLAB=vide(:,:)
-	! And now put everything in the destination file
-	! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
-	! ... Data
-	write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-
-!------------------------!
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------!
-FIELD='HV'
-UNITS='kg/kg'
-DESC='Water vapor'
-XLVL=200100.
-SLAB=waterfile(:,:,1,time_out(l))
-        ! And now put everything in the destination file
-        ! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-        ! ... Data
-        write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------!
-FIELD='HI'
-UNITS='kg/kg'
-DESC='Water ice'
-XLVL=200100.
-SLAB=watericefile(:,:,1,time_out(l))
-        ! And now put everything in the destination file
-        ! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-        ! ... Data
-        write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------!
-FIELD='TSOIL'
-UNITS='K'
-DESC='Soil temperature'
-XLVL=200100.
-SLAB=tsoilfile(:,:,1,time_out(l))
-        ! And now put everything in the destination file
-        ! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-        ! ... Data
-        write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------!
-FIELD='DSOIL'
-UNITS='m'
-DESC='Soil depths'
-XLVL=200100.
-SLAB=dsoilfile(:,:,1,time_out(l))
-        ! And now put everything in the destination file
-        ! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-        ! ... Data
-        write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------!
-FIELD='ISOIL'
-UNITS='tiu'
-DESC='Soil thermal inertia'
-XLVL=200100.
-SLAB=isoilfile(:,:,1,time_out(l))
-        ! And now put everything in the destination file
-        ! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-        ! ... Data
-        write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------!
-FIELD='DUSTQ'
-UNITS='kg/kg'
-DESC='Dust mixing ratio'
-XLVL=200100.
-SLAB=dustfile(:,:,1,time_out(l))
-        ! And now put everything in the destination file
-        ! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-        ! ... Data
-        write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------!
-FIELD='DUSTN'
-UNITS='part/kg' 
-DESC='Dust number density'
-XLVL=200100.
-SLAB=dustnfile(:,:,1,time_out(l))
-        ! And now put everything in the destination file
-        ! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-        ! ... Data
-        write(1) SLAB
-!print *,'The field '//DESC//' was written to '//output
-
-
-!!----------------------------------------------------
-!! Write file in intermediate format for WPS
-!! 2. 3D meteorological data
-!! >>> same stuff, but handle vertical profile
-!! >>> !! XLVL must be decreasing (pressure levels) 
-!!----------------------------------------------------
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='T'
-UNITS='K'
-DESC='Atmospheric temperature'
-DO k = 1,altlen
-	XLVL=levels(k)
-	SLAB=tfile(:,:,k,time_out(l))
-		! And now put everything in the destination file
-		! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-		! ... Data
-		write(1) SLAB
-END DO
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='U'
-UNITS='m s{-1}'
-DESC='Zonal wind'
-DO k = 1,altlen
-        XLVL=levels(k)
-	SLAB=ufile(:,:,k,time_out(l))
-		! And now put everything in the destination file
-		! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-		! ... Data
-		write(1) SLAB
-END DO
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='V'
-UNITS='m s{-1}'
-DESC='Meridional wind'
-DO k = 1,altlen
-        XLVL=levels(k)
-	SLAB=vfile(:,:,k,time_out(l))
-		! And now put everything in the destination file
-		! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-		! ... Data
-		write(1) SLAB
-END DO
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='RH'
-UNITS=''
-DESC='Customized 2D static field'
-DESC='Eta-levels from the GCM'
-DESC='Pressure values from the GCM'
-DO k = 1,altlen
-        XLVL=levels(k)
-SELECT CASE(ident)
-CASE('LMD')
-        SLAB=aps(k)+bps(k)*psfile(:,:,time_out(l))
-CASE('OXF')
-	SLAB=bps(k)*psfile(:,:,time_out(l))
-END SELECT
-		! And now put everything in the destination file
-		! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-		! ... Data
-		write(1) SLAB
-END DO
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!    
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------! 
-FIELD='HGT'
-UNITS='m'
-DESC='Height'
-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))))
-CASE('OXF')
-	SLAB=10000.*alog(610./(bps(k)*psfile(:,:,time_out(l))))
-END SELECT
-		! And now put everything in the destination file
-		! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
-		! ... Data
-		write(1) SLAB
-END DO
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------!
-FIELD='HV'
-UNITS='kg/kg'
-DESC='Water vapor'
-DO k = 1,altlen
-        XLVL=levels(k)
-        SLAB=waterfile(:,:,k,time_out(l))
-                ! And now put everything in the destination file
-                ! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-                ! ... Data
-                write(1) SLAB
-END DO
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------!
-FIELD='HI'
-UNITS='kg/kg'
-DESC='Water ice'
-DO k = 1,altlen
-        XLVL=levels(k)
-        SLAB=watericefile(:,:,k,time_out(l))
-                ! And now put everything in the destination file
-                ! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-                ! ... Data
-                write(1) SLAB
-END DO
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------!
-FIELD='TSOIL'
-UNITS='K'
-DESC='Soil temperature'
-DO k = 1,altlen
-        XLVL=levels(k)
-        SLAB=tsoilfile(:,:,k,time_out(l))
-                ! And now put everything in the destination file
-                ! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-                ! ... Data
-                write(1) SLAB
-END DO
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------!
-FIELD='DSOIL'
-UNITS='m'
-DESC='Soil depths'
-DO k = 1,altlen
-        XLVL=levels(k)
-        SLAB=dsoilfile(:,:,k,time_out(l))
-                ! And now put everything in the destination file
-                ! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-                ! ... Data
-                write(1) SLAB
-END DO
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------!
-FIELD='ISOIL'
-UNITS='tiu'
-DESC='Soil thermal inertia'
-DO k = 1,altlen
-        XLVL=levels(k)
-        SLAB=isoilfile(:,:,k,time_out(l))
-                ! And now put everything in the destination file
-                ! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-                ! ... Data
-                write(1) SLAB
-END DO
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------!
-FIELD='DUSTQ'
-UNITS='kg/kg'
-DESC='Dust mixing ratio'
-DO k = 1,altlen
-        XLVL=levels(k)
-        SLAB=dustfile(:,:,k,time_out(l))
-                ! And now put everything in the destination file
-                ! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-                ! ... Data
-                write(1) SLAB
-END DO
-!print *,'The field '//DESC//' was written to '//output
-
-!------------------------!
-! >>> Write a variable   !
-!    ... Copy&Paste part !
-!------------------------!
-FIELD='DUSTN'
-UNITS='part/kg'
-DESC='Dust number density'
-DO k = 1,altlen
-        XLVL=levels(k)
-        SLAB=dustnfile(:,:,k,time_out(l))
-                ! And now put everything in the destination file
-                ! ... Header
-        write(1) IFV
-        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
-        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
-                ! ... Data
-                write(1) SLAB
-END DO
-!print *,'The field '//DESC//' was written to '//output
-
-print *,'****done file '//output, int(100.*float(l)/float(FILES)), ' % '
-close(1)
-
-DEALLOCATE(SLAB)
-
-END DO !! end of the files loop
-
-
-!!-------------------------
-!! DeAllocate local arrays
-!!-------------------------
-deallocate(eta_gcm)
-deallocate(tfile)
-deallocate(tsoilfile)
-deallocate(isoilfile)
-deallocate(dsoilfile)
-deallocate(tfileorig)
-deallocate(ufile)
-deallocate(ufileorig)
-deallocate(vfile)
-deallocate(vfileorig)
-deallocate(rfile)
-deallocate(hfile)  
-deallocate(waterfile)
-deallocate(watericefile)
-deallocate(dustfile)
-deallocate(dustnfile)
-deallocate(psfile) 
-deallocate(tsfile)
-deallocate(tnsfile)
-deallocate(unsfile)
-deallocate(vnsfile)
-deallocate(emissfile)
-deallocate(co2icefile)
-deallocate(ghtsfile)	!! no scan axis
-deallocate(vide)
-deallocate(ones)
-deallocate(lat, lon, alt, time)
-deallocate(aps,bps,levels)
-
-
-print *, '------------------------'
-print *, 'End of pre-WPS process !'
-print *, '------------------------'
-print *, 'Here is what you should set in namelist.wps:' 
-print *, " start_date = '"//date_out2(1)//"'"
-print *, " end_date = '"//date_out2(FILES)//"'"
-print *, '------------------------'
-print *, 'Any date between those is ok'
-print *, 'for example, second available date is ... '//date_out2(2)
-print *, '------------------------'
-
-end
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WPS/wps_mars/metgrid/METGRID.TBL.ARW_MarsBase
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WPS/wps_mars/metgrid/METGRID.TBL.ARW_MarsBase	(revision 70)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WPS/wps_mars/metgrid/METGRID.TBL.ARW_MarsBase	(revision 72)
@@ -192,4 +192,5 @@
         flag_in_output=FLAG_SM040100
 ========================================
+# WATER ICE ON THE SURFACE
 name=SM100200
 #        interp_option=sixteen_pt+four_pt+average_4pt+search
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WPS/wps_mars/metgrid/METGRID.TBL.ARW_MarsBase_newphys
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WPS/wps_mars/metgrid/METGRID.TBL.ARW_MarsBase_newphys	(revision 70)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WPS/wps_mars/metgrid/METGRID.TBL.ARW_MarsBase_newphys	(revision 72)
@@ -192,4 +192,5 @@
         flag_in_output=FLAG_SM040100
 ========================================
+# WATER ICE ON THE SURFACE
 name=SM100200
 #        interp_option=sixteen_pt+four_pt+average_4pt+search
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WPS/wps_mars/switch_WPS_Mars
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WPS/wps_mars/switch_WPS_Mars	(revision 70)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WPS/wps_mars/switch_WPS_Mars	(revision 72)
@@ -23,5 +23,6 @@
 cd metgrid
 \rm METGRID.TBL 2> /dev/null
-cp ../wps_mars/metgrid/METGRID.TBL.ARW_MarsBase METGRID.TBL.ARW
+#cp ../wps_mars/metgrid/METGRID.TBL.ARW_MarsBase METGRID.TBL.ARW
+cp ../wps_mars/metgrid/METGRID.TBL.ARW_MarsBase_newphys METGRID.TBL.ARW
 cd src
 rm module_date_pack.f90
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM	(revision 70)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM	(revision 72)
@@ -81,5 +81,5 @@
 ####
 state  real  TSURF     ij   misc  1  -  rhd  "TSURF"     "SURF TEMPERATURE"                "K"       #SAVEMARS2 tsurf
-state  real  CO2ICE    ij   misc  1  -  rd   "CO2ICE"    "CO2 GROUND ICE"                  "kg"      #SAVEMARS2 co2ice
+state  real  CO2ICE    ij   misc  1  -  rhd  "CO2ICE"    "CO2 GROUND ICE"                  "kg m-2"  #SAVEMARS2 co2ice
 state  real  EMISSIV   ij   misc  1  -  rd   "EMISSIV"   "EMISSIVITY"                      ""        #SAVEMARS2 emis
 state  real  THETA     ij   misc  1  -  rd   "THETA"     "SLOPE INCLINATION"               "deg"     #SAVEMARS2 theta_sl
@@ -90,11 +90,12 @@
 state  real  SWUP      ij   misc  1  -  rd   "SWUP"      "UPWARD SW FLUX AT TOP"           "W m-2"   #SAVEMARS2 fluxtop_sw_tot    
 state  real  LWUP      ij   misc  1  -  rd   "LWUP"      "UPWARD LW FLUX AT TOP"           "W m-2"   #SAVEMARS2 fluxtop_lw
-state  real  MTOT      ij   misc  1  -  rd   "MTOT"      "TOTAL MASS WATER VAPOR"          "kg m-2"  #SAVEMARS2 mtot
-state  real  ICETOT    ij   misc  1  -  rd   "ICETOT"    "TOTAL MASS WATER ICE"            "kg m-2"  #SAVEMARS2 icetot
-state  real  RAVE      ij   misc  1  -  rd   "RAVE"      "MEAN ICE RADIUS"                 "m"       #SAVEMARS2 rave
+state  real  MTOT      ij   misc  1  -  rhd  "MTOT"      "TOTAL MASS WATER VAPOR"          "kg m-2"  #SAVEMARS2 mtot
+state  real  ICETOT    ij   misc  1  -  rhd  "ICETOT"    "TOTAL MASS WATER ICE"            "kg m-2"  #SAVEMARS2 icetot
+state  real  RAVE      ij   misc  1  -  rhd  "RAVE"      "MEAN ICE RADIUS"                 "m"       #SAVEMARS2 rave
 state  real  RICE      ikj  misc  1  -  rd   "RICE"      "ICE RADIUS"                      "m"       #SAVEMARS3 rice
 state  real  HR_SW     ikj  misc  1  -  rd   "HR_SW"     "HEATING RATE SW"                 "K/s"     #SAVEMARS3 zdtsw
 state  real  HR_LW     ikj  misc  1  -  rd   "HR_LW"     "HEATING RATE LW"                 "K/s"     #SAVEMARS3 zdtlw
 state  real  HR_SH     ikj  misc  1  -  rd   "HR_SH"     "HEATING RATE sens. heat"         "K/s"     #SAVEMARS3 zdtdif
+state  real  QSURF     ij   misc  1  -  rhd  "QSURF"     "TRACER AT SURFACE"               "kg m-2"  #SAVEMARS2 qsurflast
 ####
 ####
@@ -202,5 +203,6 @@
 state  real   slpy            ij   misc  1  -  i012rd  "SLOPEY"      "slope y direction"       "rad"
 state  real   mars_emiss      ij   misc  1  -  i012rd  "MARS_EMISS"  "emissivity"              "0 - 1 fraction"
-state  real   mars_cice       ij   misc  1  -  i012rd  "MARS_CICE"   "co2 ice"                 "kg"
+state  real   mars_cice       ij   misc  1  -  i012rd  "MARS_CICE"   "co2 ice"                 "kg/m2"
+state  real   mars_wice       ij   misc  1  -  i012rd  "MARS_WICE"   "h2o ice"                 "kg/m2"
 state  real   mars_tsoil     ilj   misc  1  -  i012rd  "MARS_TSOIL"  "soil temperatures"       "K"
 state  real   mars_isoil     ilj   misc  1  -  i012rd  "MARS_ISOIL"  "soil thermal inertia"    "tiu"
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F	(revision 70)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F	(revision 72)
@@ -429,4 +429,5 @@
     grid%mars_emiss(i,j)=0.95
     grid%mars_cice(i,j)=0.
+    grid%mars_wice(i,j)=0.
     grid%slpx(i,j) = 0.
     grid%slpy(i,j) = 0.
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_real.F
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_real.F	(revision 70)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_real.F	(revision 72)
@@ -277,6 +277,8 @@
               grid%mars_emiss(i,j)=grid%st000010(i,j)
               grid%mars_cice(i,j)=grid%st010040(i,j)
+              grid%mars_wice(i,j)=grid%sm100200(i,j)
                   !! one more security ... co2ice cannot be negative
                   IF (grid%mars_cice(i,j) .lt. 0.) grid%mars_cice(i,j)=0.
+                  IF (grid%mars_wice(i,j) .lt. 0.) grid%mars_wice(i,j)=0.
 
               DO k = 1, config_flags%num_soil_layers  
@@ -484,4 +486,5 @@
 grid%mars_emiss(i,j) = 0.95
 grid%mars_cice(i,j) = 0.
+grid%mars_wice(i,j) = 0.
 grid%mars_tsoil(i,:,j)=grid%mars_tsoil(its+1,:,jts+1)
 
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/solve_em.F
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/solve_em.F	(revision 70)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/solve_em.F	(revision 72)
@@ -853,4 +853,5 @@
      &        ,MARS_ALB=grid%albedo_gcm,MARS_TI=grid%therm_inert  &
      &        ,MARS_CICE=grid%mars_cice,MARS_EMISS=grid%mars_emiss   &  
+     &        ,MARS_WICE=grid%mars_wice         &
      &        ,MARS_TSOIL=grid%mars_tsoil       &
 #ifdef NEWPHYS
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/meso_physiq.F
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/meso_physiq.F	(revision 70)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/meso_physiq.F	(revision 72)
@@ -184,4 +184,5 @@
       REAL ust(ngridmx)    !! pour LES avec isfflx!=0
       LOGICAL flag_LES     !! pour LES avec isfflx!=0
+      REAL qsurflast(ngridmx) !! pour diagnostics
 c ****WRF
       REAL pday
@@ -1496,4 +1497,12 @@
          wtsurf(ig) = tsurf(ig)    !! surface temperature
          wco2ice(ig) = co2ice(ig)  !! co2 ice 
+
+         !!! TEMP TEMP TEMP TEMP TEMP TEMP TEMP
+         !!! specific to WRF WRF WRF
+         !!! just to output water ice on surface
+         !!! [it might not be water ice on surface but OK]
+         !!! uncomment the Registry entry
+         qsurflast(ig) = qsurf(ig,nqmx)
+
       enddo
 c
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/phys/module_lmd_driver.F
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/phys/module_lmd_driver.F	(revision 70)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/phys/module_lmd_driver.F	(revision 72)
@@ -30,4 +30,5 @@
         MARS_MODE, &
         MARS_ALB,MARS_TI,MARS_CICE,MARS_EMISS, &
+        MARS_WICE, &
         MARS_TSOIL, &
 #ifdef NEWPHYS
@@ -109,4 +110,5 @@
      XLAT,XLONG,HT,  &
      MARS_ALB,MARS_TI,MARS_EMISS,MARS_CICE, &
+     MARS_WICE, &
      SLPX,SLPY
 ! 3D arrays 
@@ -440,15 +442,15 @@
 !ENDIF   
 
-print *, 'check dynamics'
-  !!! in some cases, weird values are displayed
-  !!! despite the fact that outputs are OK...
-  !print *, 'u', MAXVAL(u), MINVAL(u)
-  !print *, 'v', MAXVAL(v), MINVAL(v)
-  !print *, 'w', MAXVAL(w), MINVAL(w)
-  !print *, 't', MAXVAL(t), MINVAL(t, MASK = t > 0)
-print *, 'u', u(10,1,10), u(10,15,10)
-print *, 'v', v(10,1,10), v(10,15,10)
-print *, 'w', w(10,1,10), w(10,15,10) 
-print *, 't', t(10,1,10), t(10,15,10)
+!print *, 'check dynamics'
+!  !!! in some cases, weird values are displayed
+!  !!! despite the fact that outputs are OK...
+!  !print *, 'u', MAXVAL(u), MINVAL(u)
+!  !print *, 'v', MAXVAL(v), MINVAL(v)
+!  !print *, 'w', MAXVAL(w), MINVAL(w)
+!  !print *, 't', MAXVAL(t), MINVAL(t, MASK = t > 0)
+!print *, 'u', u(10,1,10), u(10,15,10)
+!print *, 'v', v(10,1,10), v(10,15,10)
+!print *, 'w', w(10,1,10), w(10,15,10) 
+!print *, 't', t(10,1,10), t(10,15,10)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -786,4 +788,23 @@
 
 !-------------------!
+! Tracer at surface !
+!-------------------!
+#ifdef NEWPHYS
+  !!! a faire !!!
+  PRINT *, 'WARNING WARNING no tracer at surface'
+  qsurf_val(:)=0.
+#else
+SELECT CASE (MARS_MODE)
+    CASE(0)  !! NO TRACERS (mars=0)
+    qsurf_val(:)=0.
+    CASE(1)  !! WATER CYCLE (mars=1)
+    qsurf_val(2)=MARS_WICE(i,j)    !! attention... H2O ice is tracer nqmx in qsurf in LMD physics
+    qsurf_val(1)=0.                
+    CASE(2)  !! DUST CYCLE (mars=2)
+    qsurf_val(:)=0.
+END SELECT
+#endif
+
+!-------------------!
 ! Slope inclination !
 !-------------------!
@@ -832,4 +853,5 @@
 PRINT *,'v_prof ',v_prof
 PRINT *,'tsoil ',tsoil_val 
+PRINT *,'qsurf ',qsurf_val 
 #ifdef NEWPHYS
 PRINT *,'isoil ',isoil_val
@@ -844,6 +866,5 @@
 !-------------------------!
 !-------------------------!
-q2_val(:)=0      !PBL wind variance
-qsurf_val(:)=0   !Tracer on surface
+q2_val(:)=0.      !PBL wind variance
 
 !-----------------!
@@ -1167,5 +1188,5 @@
 !!PRINT *, 'min',MINVAL(RUBLTEN, MASK=RUBLTEN/=0.),&
 !!        ' at',MINLOC(RUBLTEN, MASK=RUBLTEN/=0.)
-PRINT *, RUBLTEN(10,1,10), RUBLTEN(10,15,10)
+!PRINT *, RUBLTEN(10,1,10), RUBLTEN(10,15,10)
 PRINT *, 'v non-zero tendencies'
 !!PRINT *, 'max',MAXVAL(RVBLTEN, MASK=RVBLTEN/=0.),&
@@ -1173,5 +1194,5 @@
 !!PRINT *, 'min',MINVAL(RVBLTEN, MASK=RVBLTEN/=0.),&
 !!        ' at',MINLOC(RVBLTEN, MASK=RVBLTEN/=0.)
-PRINT *, RVBLTEN(10,1,10), RVBLTEN(10,15,10)
+!PRINT *, RVBLTEN(10,1,10), RVBLTEN(10,15,10)
 !!! STOP IF CRASH
 !IF (MAXVAL(RUBLTEN, MASK=RUBLTEN/=0.) == 0.) STOP
