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)     !   
!------------------------------------------------------------------------!



!***************************************************************************
!***************************************************************************
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='OXF'


!! 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,nid3
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 :: taufile
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,"lat",latdim)
ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
  write(*,*) "latlen: ",latlen
ierr=NF_INQ_DIMID(nid,"lon",londim)
ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
  write(*,*) "lonlen: ",lonlen
ierr=NF_INQ_DIMID(nid,"sigma",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(taufile(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
taufile(:,:,:)=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, "lat",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, "lon",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 *, 'Sigma 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, "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 *,aps(1),' ... to ... ',aps(altlen)
   print *,bps(1),' ... to ... ',bps(altlen)




!!-------------------------------------------
!! Reading 2D and 3D meteorological fields
!!-------------------------------------------

IF (blank .eq. .false.) THEN

print *,'>>> Read 2D optional fields !'

print *,'Dust Opacity'
   ierr = NF_INQ_VARID (nid,"totoptdep",nvarid)
   IF (ierr .NE. NF_NOERR) THEN
      PRINT *, "Error: Readmeteo <totoptdep> not found"
      stop
   ENDIF
#ifdef NC_DOUBLE
   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, taufile)
#else
   ierr = NF_GET_VAR_REAL(nid, nvarid, taufile)
#endif
!-------------------
! at 700 Pa
taufile=taufile*700.
!-------------------


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,"t",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,:)


!------------------------
! special oxford stuff
!------------------------
!
! geopotential height ~ altimetry
!
print *,'Geopotential height at the ground'
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)



!   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,"t",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 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='Oxford Assim'
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='2D field'
XLVL=200100.
!SLAB=co2icefile(:,:,time_out(l))
SLAB=vide(:,:)
!SLAB=vide(:,:)+ptop
!SLAB=taufile(:,:,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='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))
!Oxford uses purely sigma coordinates
SLAB=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./(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)

END DO !! end of the files loop

print *, 'End of pre-WPS process !'


end
