source: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/readmeteo.F90 @ 3582

Last change on this file since 3582 was 1732, checked in by aslmd, 7 years ago

MESOSCALE. restored PREP_MARS to r1723 before readmeteo.F90 and create_readmeteo.F90 were deleted -- and history lost
File size: 48.5 KB
Line 
1program readmeteo
2
3implicit none
4include "netcdf.inc"
5
6!------------------------------------------------------------------------!
7! Readmeteo prepares an initial state for the WPS pre-processor of WRF   !
8!                                                                        !
9! Input is a diagfi.nc NETCDF file from a LMD/GCM simulation             !
10!                                                                        !
11! Outputs are WRFSI intermediate format files ready for metgrid.exe      !
12!       http://www.mmm.ucar.edu/wrf/OnLineTutorial/WPS/IM_si.htm         !
13!                                                                        !
14! **** Please create a WPSFEED folder (or a symbolic link) ****          !       
15!                                                                        !
16! A. Spiga - 16/04/2007                                                  !
17!            22/05/2007 : need to run zrecast as a preliminary           !
18!            07/06/2007 : external parameter file 'readmeteo.def         !
19!            30/07/2007 : manage additional variables (tsoil, emiss,...) !     
20!            19/10/2007 : no need to run zrecast at all (eta levels)     !   
21!               12/2007 : include co2ice and emissivity                  !
22!               02/2008 : include water vapor and ice                    !
23!               01/2010 : possible use of diagsoil for new physics       !
24!               02/2011 : add dust tracers, correct surface              !
25!                                                                        !
26! spiga@lmd.jussieu.fr                                                   !     
27!------------------------------------------------------------------------!
28
29
30
31!***************************************************************************
32!***************************************************************************
33REAL, PARAMETER :: emiss_prescribed=0.95
34REAL, PARAMETER :: co2ice_prescribed=0.
35CHARACTER (LEN=3), PARAMETER :: ident='LMD'
36
37
38
39!***************************************************************************
40!***************************************************************************
41
42REAL :: ptop
43REAL, PARAMETER :: grav=3.72
44LOGICAL, PARAMETER :: blank=.false.
45
46
47
48!! Variables to be written in the output file
49!! *** NB: conformity with metgrid/src/read_met_module.F90
50CHARACTER (LEN=33) :: output
51INTEGER :: IFV
52CHARACTER (LEN=24) :: HDATE
53REAL :: XFCST
54CHARACTER (LEN=32) :: SOURCE
55CHARACTER (LEN=9) :: FIELD
56CHARACTER (LEN=25) :: UNITS
57CHARACTER (LEN=46) :: DESC
58REAL :: XLVL
59INTEGER :: NX,NY,IPROJ
60CHARACTER (LEN=8) :: STARTLOC
61REAL :: STARTLAT,STARTLON,DELTALAT,DELTALON
62REAL, POINTER, DIMENSION(:,:) :: SLAB
63
64!! Variables related to NETCDF format
65integer :: nid,nid2,nid3,nvarid,ierr,ierr2
66integer :: timedim,londim,latdim,altdim,altdim2
67integer :: rlatvdim,rlonudim
68integer :: timelen,lonlen,latlen,altlen,altlen2
69integer :: rlatvlen,rlonulen
70integer :: soillen,soildim
71integer :: nphys
72integer, dimension(4) :: corner,edges
73
74!! Intermediate data arrays
75integer :: k,l,m,n,p
76real, dimension(:), allocatable :: lat,lon,time,alt,aps,bps,levels,vertdsoil
77real, dimension(:,:), allocatable :: vide,ones,ghtsfile
78real, dimension(:,:), allocatable :: interm
79real, dimension(:,:,:), allocatable :: gwparam
80real, dimension(:,:,:), allocatable :: psfile,tsfile
81real, dimension(:,:,:), allocatable :: emissfile,co2icefile
82real, dimension(:,:,:), allocatable :: tnsfile,unsfile,vnsfile
83real, dimension(:,:,:,:), allocatable :: tfile,ufile,vfile,rfile,hfile
84real, dimension(:,:,:,:), allocatable :: eta_gcm
85!real, dimension(:,:,:,:), allocatable :: tfileorig,ufileorig,vfileorig
86real, dimension(:,:,:,:), allocatable :: tsoilfile, dsoilfile, isoilfile
87real, dimension(:,:,:,:), allocatable :: waterfile, watericefile
88real, dimension(:,:,:), allocatable :: swatericefile!, swaterfile
89real, dimension(:,:,:,:), allocatable :: dustfile,dustnfile
90real, dimension(:,:,:,:), allocatable :: co2file
91real, dimension(:,:,:,:), allocatable :: ccnfile,ccnnfile
92
93!! Reading the parameter file
94integer :: tmp,increment,FILES
95integer :: tmp2,tmp3,tmp4
96character*1 :: answer
97character*4 :: str
98character*2 :: str2, str3, str4
99integer, dimension(:), allocatable :: time_out
100character*13, dimension(:), allocatable :: date_out
101character*19, dimension(:), allocatable :: date_out2
102
103#ifdef PHOTOCHEM
104real, dimension(:,:,:,:,:), allocatable :: chemtrac
105integer :: nchemtrac,i
106CHARACTER*20,DIMENSION(:),ALLOCATABLE :: wtnom
107#endif
108
109
110!***************************************************************************
111!***************************************************************************
112
113
114!!---------------------
115!! Open the datafile
116!!---------------------
117ierr=NF_OPEN ("input_diagfi.nc",NF_NOWRITE,nid)
118IF (ierr.NE.NF_NOERR) THEN
119   write(*,*)'**** Please create a symbolic link called input_diagfi.nc'
120   CALL ABORT
121ENDIF
122
123
124!!--------------------------
125!! Ask for data references
126!!--------------------------
127!write(*,*) "Create n files. What is n ?"
128read(*,*) FILES
129allocate(time_out(FILES))
130allocate(date_out(FILES))
131allocate(date_out2(FILES))
132
133!write(*,*) ""
134!write(*,*) "INPUT: Diagfi time"
135!write(*,*) "list of n subscripts, separated by <Return>s"
136increment=1
137do while (increment.ne.FILES+1)
138    read(*,*) tmp
139    time_out(increment)=tmp
140    increment=increment+1
141enddo
142
143!write(*,*) ""
144!write(*,*) "OUTPUT: WRF time"
145!write(*,*) "fill 4 lines indicating"
146!write(*,*) "-year (4 digits)"
147!write(*,*) "-month (2 digits)"
148!write(*,*) "-day (2 digits)"
149!write(*,*) "-hour (2 digits)"
150increment=1
151do while (increment.ne.FILES+1)
152    read(*,*) str
153    read(*,*) str2
154    read(*,*) str3
155    read(*,*) str4
156    date_out(increment)=str//'-'//str2//'-'//str3//'_'//str4
157    date_out2(increment)=str//'-'//str2//'-'//str3//'_'//str4//':00:00'
158    !print *,date_out(increment)
159    !write(*,*) "ok? (y/n)"
160    read(*,*) answer
161    if (answer.eq.'n') then
162        !write(*,*) "please write again"
163    else
164        !write(*,*) "next one, please"   
165        increment=increment+1
166    endif
167enddo
168
169
170!!-------------------
171!! Get array sizes
172!!-------------------
173SELECT CASE(ident)
174CASE('LMD')
175ierr=NF_INQ_DIMID(nid,"Time",timedim)
176CASE('OXF')
177ierr=NF_INQ_DIMID(nid,"time",timedim)
178END SELECT
179ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
180  write(*,*) "timelen: ",timelen
181
182SELECT CASE(ident)
183CASE('LMD')
184ierr=NF_INQ_DIMID(nid,"latitude",latdim)
185CASE('OXF')
186ierr=NF_INQ_DIMID(nid,"lat",latdim)
187END SELECT
188ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
189  write(*,*) "latlen: ",latlen
190
191SELECT CASE(ident)
192CASE('LMD')
193ierr=NF_INQ_DIMID(nid,"longitude",londim)
194CASE('OXF')
195ierr=NF_INQ_DIMID(nid,"lon",londim)
196END SELECT
197ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
198  write(*,*) "lonlen: ",lonlen
199
200SELECT CASE(ident)
201CASE('LMD')
202ierr=NF_INQ_DIMID(nid,"altitude",altdim)
203CASE('OXF')
204ierr=NF_INQ_DIMID(nid,"sigma",altdim)
205END SELECT
206ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
207  write(*,*) "altlen: ",altlen
208
209
210
211!!-------------------------
212!! Allocate local arrays
213!!-------------------------
214allocate(eta_gcm(lonlen,latlen,altlen,timelen))
215allocate(tfile(lonlen,latlen,altlen,timelen))
216allocate(tsoilfile(lonlen,latlen,altlen,timelen))
217allocate(dsoilfile(lonlen,latlen,altlen,timelen))
218allocate(isoilfile(lonlen,latlen,altlen,timelen))
219allocate(vertdsoil(altlen))
220!allocate(tfileorig(lonlen,latlen,altlen,timelen))
221allocate(ufile(lonlen,latlen,altlen,timelen))
222!allocate(ufileorig(lonlen,latlen,altlen,timelen))
223allocate(vfile(lonlen,latlen,altlen,timelen))
224!allocate(vfileorig(lonlen,latlen,altlen,timelen))
225allocate(rfile(lonlen,latlen,altlen,timelen))
226allocate(hfile(lonlen,latlen,altlen,timelen)) 
227allocate(waterfile(lonlen,latlen,altlen,timelen))
228allocate(watericefile(lonlen,latlen,altlen,timelen))
229!allocate(swaterfile(lonlen,latlen,timelen))
230allocate(swatericefile(lonlen,latlen,timelen))
231allocate(dustfile(lonlen,latlen,altlen,timelen))
232allocate(dustnfile(lonlen,latlen,altlen,timelen))
233allocate(co2file(lonlen,latlen,altlen,timelen))
234allocate(ccnfile(lonlen,latlen,altlen,timelen))
235allocate(ccnnfile(lonlen,latlen,altlen,timelen))
236allocate(psfile(lonlen,latlen,timelen))
237allocate(tsfile(lonlen,latlen,timelen))
238allocate(tnsfile(lonlen,latlen,timelen))
239allocate(unsfile(lonlen,latlen,timelen))
240allocate(vnsfile(lonlen,latlen,timelen))
241allocate(emissfile(lonlen,latlen,timelen))
242allocate(co2icefile(lonlen,latlen,timelen))
243allocate(interm(lonlen,latlen))
244allocate(gwparam(lonlen,latlen,5))
245allocate(ghtsfile(lonlen,latlen))    !! no scan axis
246allocate(vide(lonlen,latlen))
247allocate(ones(lonlen,latlen))
248allocate(lat(latlen), lon(lonlen), alt(altlen), time(timelen))
249allocate(aps(altlen),bps(altlen),levels(altlen))
250#ifdef PHOTOCHEM
251nchemtrac = 14
252allocate(wtnom(nchemtrac))
253wtnom(1)  = "c_co2"
254wtnom(2)  = "c_co"
255wtnom(3)  = "c_o"
256wtnom(4)  = "c_o1d"
257wtnom(5)  = "c_o2"
258wtnom(6)  = "c_o3"
259wtnom(7)  = "c_h"
260wtnom(8)  = "c_h2"
261wtnom(9)  = "c_oh"
262wtnom(10) = "c_ho2"
263wtnom(11) = "c_h2o2"
264wtnom(12) = "c_ch4"
265wtnom(13) = "c_n2"
266wtnom(14) = "c_ar"
267allocate(chemtrac(lonlen,latlen,altlen,timelen,nchemtrac))
268chemtrac(:,:,:,:,:)=0
269#endif
270
271tfile(:,:,:,:)=0
272tsoilfile(:,:,:,:)=0
273isoilfile(:,:,:,:)=0
274dsoilfile(:,:,:,:)=0
275vertdsoil(:)=0.
276!tfileorig(:,:,:,:)=0
277!ufileorig(:,:,:,:)=0
278!vfileorig(:,:,:,:)=0
279ufile(:,:,:,:)=0
280vfile(:,:,:,:)=0
281rfile(:,:,:,:)=0
282hfile(:,:,:,:)=0
283waterfile(:,:,:,:)=0
284watericefile(:,:,:,:)=0
285!swaterfile(:,:,:)=0
286swatericefile(:,:,:)=0
287dustfile(:,:,:,:)=0
288dustnfile(:,:,:,:)=0
289co2file(:,:,:,:)=0
290ccnfile(:,:,:,:)=0
291ccnnfile(:,:,:,:)=0
292psfile(:,:,:)=0
293tsfile(:,:,:)=0
294emissfile(:,:,:)=0
295co2icefile(:,:,:)=0
296tnsfile(:,:,:)=0
297unsfile(:,:,:)=0
298vnsfile(:,:,:)=0
299interm(:,:)=0
300gwparam(:,:,:)=0
301ghtsfile(:,:)=0
302vide(:,:)=0
303ones(:,:)=0
304
305
306!!-------------------
307!! Read dimensions
308!!-------------------
309
310print *,'>>> Read dimensions !'
311
312print *,'Time'
313SELECT CASE(ident)
314CASE('LMD')
315   ierr = NF_INQ_VARID (nid, "Time",nvarid)
316CASE('OXF')
317   ierr = NF_INQ_VARID (nid, "time",nvarid)
318END SELECT
319   IF (ierr .NE. NF_NOERR) THEN
320      PRINT *, "Error: Readmeteo <Time> not found"
321      stop
322   ENDIF
323#ifdef NC_DOUBLE
324   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
325#else
326   ierr = NF_GET_VAR_REAL(nid, nvarid, time)
327#endif
328   print *,time(1),' ... to ... ',time(timelen)
329
330print *,'Latitude'
331SELECT CASE(ident)
332CASE('LMD')
333   ierr = NF_INQ_VARID (nid, "latitude",nvarid)
334CASE('OXF')
335   ierr = NF_INQ_VARID (nid, "lat",nvarid)
336END SELECT
337   IF (ierr .NE. NF_NOERR) THEN
338      PRINT *, "Error: Readmeteo <latitude> not found"
339      stop
340   ENDIF
341#ifdef NC_DOUBLE
342   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lat)
343#else
344   ierr = NF_GET_VAR_REAL(nid, nvarid, lat)
345#endif
346   print *,lat(1),' ... to ... ',lat(latlen),' ... step: ',lat(latlen)-lat(latlen-1)
347
348print *,'Longitude'
349SELECT CASE(ident)
350CASE('LMD')
351   ierr = NF_INQ_VARID (nid, "longitude",nvarid)
352CASE('OXF')
353   ierr = NF_INQ_VARID (nid, "lon",nvarid)
354END SELECT
355   IF (ierr .NE. NF_NOERR) THEN
356      PRINT *, "Error: Readmeteo <longitude> not found"
357      stop
358   ENDIF
359#ifdef NC_DOUBLE
360   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lon)
361#else
362   ierr = NF_GET_VAR_REAL(nid, nvarid, lon)
363#endif
364   print *,lon(1),' ... to ... ',lon(lonlen),' ... step: ',lon(lonlen)-lon(lonlen-1)
365
366SELECT CASE(ident)
367CASE('LMD')
368print *, 'Hybrid coordinates'
369   ierr = NF_INQ_VARID (nid, "aps", nvarid)
370   IF (ierr .NE. NF_NOERR) THEN
371      PRINT *, "Error: Readmeteo <aps> not found"
372      stop
373   ENDIF
374#ifdef NC_DOUBLE
375   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, aps)
376#else
377   ierr = NF_GET_VAR_REAL(nid, nvarid, aps)
378#endif
379   ierr = NF_INQ_VARID (nid, "bps", nvarid)
380   IF (ierr .NE. NF_NOERR) THEN
381      PRINT *, "Error: Readmeteo <bps> not found"
382      stop
383   ENDIF
384#ifdef NC_DOUBLE
385   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bps)
386#else
387   ierr = NF_GET_VAR_REAL(nid, nvarid, bps)
388#endif
389   print *,aps(1),' ... to ... ',aps(altlen)
390   print *,bps(1),' ... to ... ',bps(altlen)
391CASE('OXF')
392   ierr = NF_INQ_VARID (nid, "sigma", nvarid)
393   IF (ierr .NE. NF_NOERR) THEN
394      PRINT *, "Error: Readmeteo <sigma> not found"
395      stop
396   ENDIF
397#ifdef NC_DOUBLE
398   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bps)
399#else
400   ierr = NF_GET_VAR_REAL(nid, nvarid, bps)
401#endif
402   print *,bps(1),' ... to ... ',bps(altlen)
403END SELECT
404
405
406!!-------------------------------------------
407!! Reading 2D and 3D meteorological fields
408!!-------------------------------------------
409
410IF (blank .EQV. .false.) THEN
411
412print *,'>>> Read 2D optional fields !'
413
414print *,'Emissivity'
415   ierr = NF_INQ_VARID (nid,"emis",nvarid)
416IF (ierr .NE. NF_NOERR) THEN
417        PRINT *, '...warning: not found in diagfi !'
418        PRINT *, '...will be filled with a prescribed value', emiss_prescribed
419  emissfile(:,:,:)=emiss_prescribed
420ELSE
421#ifdef NC_DOUBLE
422  ierr = NF_GET_VAR_DOUBLE(nid, nvarid, emissfile)
423#else
424  ierr = NF_GET_VAR_REAL(nid, nvarid, emissfile)
425#endif
426ENDIF   
427
428print *,'CO2 ice'
429   ierr = NF_INQ_VARID (nid,"co2ice",nvarid)
430IF (ierr .NE. NF_NOERR) THEN
431  PRINT *, '...warning: not found in diagfi !'
432  PRINT *, '...will be filled with a prescribed value', co2ice_prescribed
433  co2icefile(:,:,:)=co2ice_prescribed
434ELSE
435#ifdef NC_DOUBLE
436  ierr = NF_GET_VAR_DOUBLE(nid, nvarid, co2icefile)
437#else
438  ierr = NF_GET_VAR_REAL(nid, nvarid, co2icefile)
439#endif
440ENDIF
441
442print *,'>>> Read 2D surface fields !'
443
444print *,'Surface Pressure'
445   ierr = NF_INQ_VARID (nid,"ps",nvarid)
446   IF (ierr .NE. NF_NOERR) THEN
447      PRINT *, "Error: Readmeteo <ps> not found"
448      stop
449   ENDIF
450#ifdef NC_DOUBLE
451   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, psfile)
452#else
453   ierr = NF_GET_VAR_REAL(nid, nvarid, psfile)
454#endif
455
456print *,'Ground Temperature'
457   ierr = NF_INQ_VARID (nid,"tsurf",nvarid)
458   IF (ierr .NE. NF_NOERR) THEN
459      PRINT *, "Error: Readmeteo <tsurf> not found"
460      stop
461   ENDIF
462#ifdef NC_DOUBLE
463   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsfile)
464#else
465   ierr = NF_GET_VAR_REAL(nid, nvarid, tsfile)
466#endif
467
468
469!!!"atmospheric" surface temperature is taken
470!!!from original diagfi.nc first level
471!!!... level is ~3-5 meters
472!print *,'Near-Surface Temperature'
473!   ierr = NF_INQ_VARID (nid,"temp",nvarid)
474!   IF (ierr .NE. NF_NOERR) THEN
475!      ierr = NF_INQ_VARID (nid,"t",nvarid)
476!       IF (ierr .NE. NF_NOERR) THEN
477!           PRINT *, "Error: Readmeteo <temp> not found"
478!           stop
479!       ENDIF
480!   ENDIF
481!#ifdef NC_DOUBLE
482!   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfileorig)
483!#else
484!   ierr = NF_GET_VAR_REAL(nid, nvarid, tfileorig)
485!#endif
486!   tnsfile=tfileorig(:,:,1,:)
487
488!!!"atmospheric" surface u is taken
489!!!from original diagfi.nc first level
490!!!... level is ~3-5 meters
491!print *,'Near-Surface Zonal Wind'
492!   ierr = NF_INQ_VARID (nid,"u",nvarid)
493!   IF (ierr .NE. NF_NOERR) THEN
494!     PRINT *, "Error: Readmeteo <u> not found"
495!     stop
496!   ENDIF
497!#ifdef NC_DOUBLE
498!   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufileorig)
499!#else
500!   ierr = NF_GET_VAR_REAL(nid, nvarid, ufileorig)
501!#endif
502!   unsfile=ufileorig(:,:,1,:)
503
504!!!"atmospheric" surface v is taken
505!!!from original diagfi.nc first level
506!!!... level is ~3-5 meters
507!print *,'Near-Surface Meridional Wind'
508!   ierr = NF_INQ_VARID (nid,"v",nvarid)
509!   IF (ierr .NE. NF_NOERR) THEN
510!     PRINT *, "Error: Readmeteo <v> not found"
511!     stop
512!   ENDIF
513!#ifdef NC_DOUBLE
514!   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfileorig)
515!#else
516!   ierr = NF_GET_VAR_REAL(nid, nvarid, vfileorig)
517!#endif
518!   vnsfile=vfileorig(:,:,1,:)
519
520SELECT CASE(ident)
521CASE('LMD')
522   print *,'Geopotential height at the ground'
523   ierr = NF_INQ_VARID (nid,"phisinit",nvarid)
524   IF (ierr .NE. NF_NOERR) THEN
525      PRINT *, "Error: Readmeteo <phisinit> not found"
526      stop
527   ENDIF
528#ifdef NC_DOUBLE
529   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ghtsfile)
530#else
531   ierr = NF_GET_VAR_REAL(nid, nvarid, ghtsfile)
532#endif
533!**** from geopotential to geopotential height
534ghtsfile=ghtsfile/grav
535!****
536CASE('OXF')
537!
538! geopotential height ~ altimetry
539!
540print *,'Geopotential height at the ground from file mountain_new.nc'
541ierr=NF_OPEN("mountain_new.nc",NF_NOWRITE,nid3)
542    if (ierr.ne.NF_NOERR) then
543      write(*,*) "Error: Could not open that file either"
544      stop "Might as well stop here"
545    endif
546ierr=NF_INQ_VARID(nid3,"orography",nvarid)
547ierr=NF_GET_VAR_REAL(nid3,nvarid,ghtsfile)
548if (ierr.ne.NF_NOERR) then
549    stop "Error: Failed reading phisinit"
550endif
551ierr=NF_CLOSE(nid3)
552END SELECT
553
554
555print *,'>>> Read 3D meteorological fields ! - This may take some time ...'
556
557print *,'Temperature'
558   ierr = NF_INQ_VARID (nid,"temp",nvarid)
559   IF (ierr .NE. NF_NOERR) THEN
560      ierr = NF_INQ_VARID (nid,"t",nvarid)
561        IF (ierr .NE. NF_NOERR) THEN
562          PRINT *, "Error: Readmeteo <t> not found"
563          stop
564        ENDIF
565   ENDIF
566#ifdef NC_DOUBLE
567   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfile)
568#else
569   ierr = NF_GET_VAR_REAL(nid, nvarid, tfile)
570#endif
571tnsfile=tfile(:,:,1,:)
572
573print *,'Zonal wind'   
574   ierr = NF_INQ_VARID (nid,"u",nvarid)
575   IF (ierr .NE. NF_NOERR) THEN
576      PRINT *, "Error: Readmeteo <u> not found"
577      stop
578   ENDIF
579#ifdef NC_DOUBLE
580   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufile)
581#else
582   ierr = NF_GET_VAR_REAL(nid, nvarid, ufile)
583#endif
584unsfile=ufile(:,:,1,:)
585
586print *,'Meridional wind'
587   ierr = NF_INQ_VARID (nid,"v",nvarid)
588   IF (ierr .NE. NF_NOERR) THEN
589      PRINT *, "Error: Readmeteo <v> not found"
590      stop
591   ENDIF
592#ifdef NC_DOUBLE
593   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfile)
594#else
595   ierr = NF_GET_VAR_REAL(nid, nvarid, vfile)
596#endif
597vnsfile=ufile(:,:,1,:)
598
599
600!!------------------------
601!! special water stuff
602!!------------------------
603    print *,'Water vapor'
604    ierr=NF_INQ_VARID(nid,"q02",nvarid)
605    if (ierr.ne.NF_NOERR) then
606      write(*,*) "...No q02 - Water vapor set to 0"
607      waterfile(:,:,:,:)=0. 
608    else
609      ierr=NF_GET_VAR_REAL(nid,nvarid,waterfile)
610    endif   
611
612    print *,'Water ice'
613    ierr=NF_INQ_VARID(nid,"q01",nvarid)
614    if (ierr.ne.NF_NOERR) then
615      write(*,*) "...No q01 - Water ice set to 0" 
616      watericefile(:,:,:,:)=0.
617    else
618      ierr=NF_GET_VAR_REAL(nid,nvarid,watericefile)
619    endif
620!    print *,'Surface Water vapor'
621!    ierr=NF_INQ_VARID(nid,"qsurf02",nvarid)
622!    if (ierr.ne.NF_NOERR) then
623!      write(*,*) "...No qsurf02 - surface Water vapor set to 0"
624!      swaterfile(:,:,:)=0.
625!    endif
626!    ierr=NF_GET_VAR_REAL(nid,nvarid,swaterfile)
627
628    print *,'Surface Water ice'
629!!!!!! ATTENTION ATTENTION
630!!!!!! water ice a la surface est qsurf(ig,nqmx)
631    ierr=NF_INQ_VARID(nid,"qsurf02",nvarid)
632    if (ierr.ne.NF_NOERR) then
633      write(*,*) "...No qsurf02 - surface Water ice set to 0"
634      swatericefile(:,:,:)=0.
635    else
636      ierr=NF_GET_VAR_REAL(nid,nvarid,swatericefile)
637    endif
638!!------------------------
639!! special water stuff
640!!------------------------
641
642    print *,'CO2 mass mixing ratio'
643    ierr=NF_INQ_VARID(nid,"co2",nvarid)
644    if (ierr.ne.NF_NOERR) then
645      write(*,*) "...No co2 - co2 mixing ratio set to 0.95"
646      co2file(:,:,:,:)=0.95
647    else
648      ierr=NF_GET_VAR_REAL(nid,nvarid,co2file)
649    endif
650
651!!------------------------
652!! special dust stuff
653!!------------------------
654    print *,'Dust mass'
655    ierr=NF_INQ_VARID(nid,"dustq",nvarid)
656    if (ierr.ne.NF_NOERR) then
657      write(*,*) "...No dustq - Dust mass set to 0"
658      dustfile(:,:,:,:)=0.
659    else
660      ierr=NF_GET_VAR_REAL(nid,nvarid,dustfile)
661    endif
662
663    print *,'Dust number'
664    ierr=NF_INQ_VARID(nid,"dustN",nvarid)
665    if (ierr.ne.NF_NOERR) then
666      write(*,*) "...No dustN - Dust number set to 0"
667      dustnfile(:,:,:,:)=0.
668    else
669      ierr=NF_GET_VAR_REAL(nid,nvarid,dustnfile)
670    endif
671
672    print *,'CCN mass'
673    ierr=NF_INQ_VARID(nid,"ccn",nvarid)
674    if (ierr.ne.NF_NOERR) then
675      write(*,*) "...No ccn - CCN mass set to 0"
676      ccnfile(:,:,:,:)=0.
677    else
678      ierr=NF_GET_VAR_REAL(nid,nvarid,ccnfile)
679    endif
680
681    print *,'CCN number'
682    ierr=NF_INQ_VARID(nid,"ccnN",nvarid)
683    if (ierr.ne.NF_NOERR) then
684      write(*,*) "...No ccnN - CCN number set to 0"
685      ccnnfile(:,:,:,:)=0.
686    else
687      ierr=NF_GET_VAR_REAL(nid,nvarid,ccnnfile)
688    endif
689!!------------------------
690!! special dust stuff
691!!------------------------
692
693!SELECT CASE(ident)
694!CASE('LMD')
695
696    print *,'Soil temperatures'
697    ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
698    if (ierr.ne.NF_NOERR) then
699        write(*,*) "...No tsoil - Soil temperatures set isothermal with surface temperature"
700        DO l=1,altlen
701                tsoilfile(:,:,l,:)=tsfile(:,:,:)
702        ENDDO
703    else
704        ierr=NF_GET_VAR_REAL(nid,nvarid,tsoilfile)
705    endif
706
707!!!!!!!!
708!!!!!!!! new physics (but still compatible with old physics)
709    print *,'Soil depths'
710    ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
711    if (ierr.ne.NF_NOERR) then
712        write(*,*) "...No soildepth - Set to -999"  !!! see soil_settings in LMD physics
713        DO l=1,altlen
714                vertdsoil(l)=-999.
715        ENDDO
716    else
717        ierr=NF_GET_VAR_REAL(nid,nvarid,vertdsoil)
718    endif
719    print *, 'wait a minute' !! AS: I know this could be better
720    DO m=1,lonlen
721     DO n=1,latlen
722      DO p=1,timelen
723       dsoilfile(m,n,:,p) = vertdsoil(:)
724      ENDDO
725     ENDDO
726    ENDDO
727    DEALLOCATE(vertdsoil)
728
729    print *,'Soil thermal inertia'
730    ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
731    if (ierr.ne.NF_NOERR) then
732        write(*,*) "...No soil therm. inert. - Set to -999"
733        DO l=1,altlen
734                isoilfile(:,:,l,:)=-999.
735        ENDDO
736    else
737        ierr=NF_GET_VAR_REAL(nid,nvarid,isoilfile)
738    endif
739!!!!!!!!
740!!!!!!!! new physics
741
742
743!!!!!!!!!!!!!!!!!!!!!!!!NEW PHYSICS + PHOTOCHEM
744!!!!!!!!!!!!!!!!!!!!!!!!NEW PHYSICS + PHOTOCHEM
745#ifdef PHOTOCHEM
746    print *,'photochem'
747    DO i=1,nchemtrac
748     print *,wtnom(i)
749     ierr=NF_INQ_VARID(nid,wtnom(i),nvarid)
750     if (ierr.ne.NF_NOERR) then
751       write(*,*) "...No ",wtnom(i), " - set to 0"
752       chemtrac(:,:,:,:,i)=0.
753     else
754       ierr=NF_GET_VAR_REAL(nid,nvarid,chemtrac(:,:,:,:,i))
755     endif
756    ENDDO
757#endif
758!!!!!!!!!!!!!!!!!!!!!!!!NEW PHYSICS + PHOTOCHEM
759!!!!!!!!!!!!!!!!!!!!!!!!NEW PHYSICS + PHOTOCHEM
760
761
762
763!END SELECT
764
765ierr=NF_CLOSE(nid)
766
767!!!!lott stuff
768!print *,'get lott'
769!ierr=NF_OPEN("init_lott.nc",NF_NOWRITE,nid3)
770!ierr=NF_INQ_VARID(nid3,"ZMEA",nvarid)
771!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
772!gwparam(:,:,1)=interm(:,:)
773!ierr=NF_INQ_VARID(nid3,"ZSTD",nvarid)
774!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
775!gwparam(:,:,2)=interm(:,:)
776!ierr=NF_INQ_VARID(nid3,"ZSIG",nvarid)
777!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
778!gwparam(:,:,3)=interm(:,:)
779!ierr=NF_INQ_VARID(nid3,"ZGAM",nvarid)
780!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
781!gwparam(:,:,4)=interm(:,:)
782!ierr=NF_INQ_VARID(nid3,"ZTHE",nvarid)
783!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
784!gwparam(:,:,5)=interm(:,:)
785!ierr=NF_CLOSE(nid3)
786
787ENDIF
788
789!!-----------------------------
790!! Loop on the written files
791!! >>> what follows is pretty repetitive ...
792!!        gotta do something (one more loop?)
793!!-----------------------------
794
795!!! Equivalent eta levels for WRF interpolation in initialize_real
796!print *,'Computing equivalent eta levels'
797!       
798!       !ptop will be passed through RH(surface)
799!       ptop=MINVAL(aps(altlen)+bps(altlen)*psfile(:,:,:))
800!       print *, 'ptop', ptop
801!
802!print *, 'sample: equivalent eta levels at i=1,j=1'
803!DO k = 1,altlen
804!        levels(k)=-k
805!        eta_gcm(:,:,k,:)=(aps(k)+bps(k)*psfile(:,:,:)-ptop)/(psfile(:,:,:)-ptop)
806!        print *, eta_gcm(1,1,k,1)
807!END DO
808
809!!---better to pass pressure values
810!!---the eta treatment is done in module_initialize_real
811DO k = 1,altlen
812        levels(k)=-k
813        !!dummy decreasing levels
814END DO
815
816
817
818
819!! Dummy surface level is XLVL=200100.
820
821
822DO l=1,FILES
823
824
825!!---------------------------------------------
826!! Write file in intermediate format for WPS
827!! 1. Surface data
828!!---------------------------------------------
829!!
830!! a mettre pour tous sinon
831!!     WRF_DEBUG: Warning DIM             4 , NAME num_metgrid_levels REDIFINED  by
832!!            var DUSTN            26           25  in wrf_io.F90 line         2349
833!!
834
835!
836! These variables remain the same in the "loop"
837!
838HDATE=date_out(l)
839IFV=4
840XFCST=0.
841SOURCE=ident
842NX=lonlen
843NY=latlen
844ALLOCATE(SLAB(NX,NY))
845IPROJ=0
846STARTLOC='SWCORNER'
847STARTLAT=lat(1)
848STARTLON=lon(1)
849DELTALAT=lat(latlen)-lat(latlen-1)
850DELTALON=lon(lonlen)-lon(lonlen-1)
851!
852! Open the data destination file
853!
854output="./WPSFEED/"//ident//":"//date_out2(l)       
855open(UNIT=1,                    &
856        FILE=output,            &
857        STATUS='REPLACE',       &
858        FORM='UNFORMATTED',     &
859        ACCESS='SEQUENTIAL')
860
861!------------------------!   
862! >>> Write a variable   !
863!    ... Copy&Paste part !
864!------------------------!
865FIELD='T'
866UNITS='K'
867DESC='Atmospheric temperature'
868XLVL=200100.
869SLAB=tnsfile(:,:,time_out(l))
870        ! And now put everything in the destination file
871        ! ... Header
872        write(1) IFV
873        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
874        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
875        ! ... Data
876        write(1) SLAB
877!print *,'The field '//DESC//' was written to '//output
878
879!------------------------!   
880! >>> Write a variable   !
881!    ... Copy&Paste part !
882!------------------------!
883FIELD='U'
884UNITS='m s{-1}'
885DESC='Zonal wind'
886XLVL=200100.
887SLAB=unsfile(:,:,time_out(l))
888        ! And now put everything in the destination file
889        ! ... Header
890        write(1) IFV
891        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
892        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
893        ! ... Data
894        write(1) SLAB
895!print *,'The field '//DESC//' was written to '//output
896
897!------------------------!   
898! >>> Write a variable   !
899!    ... Copy&Paste part !
900!------------------------!
901FIELD='V'
902UNITS='m s{-1}'
903DESC='Meridional wind'
904XLVL=200100.
905SLAB=vnsfile(:,:,time_out(l))
906        ! And now put everything in the destination file
907        ! ... Header
908        write(1) IFV
909        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
910        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
911        ! ... Data
912        write(1) SLAB
913!print *,'The field '//DESC//' was written to '//output
914
915!------------------------!   
916! >>> Write a variable   !
917!    ... Copy&Paste part !
918!------------------------!
919FIELD='RH'
920UNITS=''
921DESC='Customized 2D static field'
922XLVL=200100.
923SLAB=vide(:,:)
924        ! And now put everything in the destination file
925        ! ... Header
926        write(1) IFV
927        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
928        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
929        ! ... Data
930        write(1) SLAB
931!print *,'The field '//DESC//' was written to '//output
932
933!------------------------!   
934! >>> Write a variable   !
935!    ... Copy&Paste part !
936!------------------------!
937FIELD='SOILHGT'
938UNITS='m'
939DESC='Terrain field of source analysis'  !Ground geopotential height
940XLVL=200100.
941SLAB=ghtsfile(:,:)
942        ! And now put everything in the destination file
943        ! ... Header
944        write(1) IFV
945        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
946        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
947        ! ... Data
948        write(1) SLAB
949!print *,'The field '//DESC//' was written to '//output
950
951!------------------------!   
952! >>> Write a variable   !
953!    ... Copy&Paste part !
954!------------------------!
955FIELD='PSFC'
956UNITS='Pa'
957DESC='Surface Pressure'
958XLVL=200100.
959SLAB=psfile(:,:,time_out(l))
960        ! And now put everything in the destination file
961        ! ... Header
962        write(1) IFV
963        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
964        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
965        ! ... Data
966        write(1) SLAB
967!print *,'The field '//DESC//' was written to '//output
968
969!------------------------!   
970! >>> Write a variable   !
971!    ... Copy&Paste part !
972!------------------------!
973FIELD='ST000010'
974UNITS=''
975DESC='Emissivity'
976XLVL=200100.
977SLAB=emissfile(:,:,time_out(l))
978        ! And now put everything in the destination file
979        ! ... Header
980        write(1) IFV
981        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
982        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
983        ! ... Data
984        write(1) SLAB
985!print *,'The field '//DESC//' was written to '//output
986
987!------------------------!   
988! >>> Write a variable   !
989!    ... Copy&Paste part !
990!------------------------!
991FIELD='ST010040'
992UNITS=''
993DESC='CO2 ice'
994XLVL=200100.
995SLAB=co2icefile(:,:,time_out(l))
996        ! And now put everything in the destination file
997        ! ... Header
998        write(1) IFV
999        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1000        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1001        ! ... Data
1002        write(1) SLAB
1003!print *,'The field '//DESC//' was written to '//output
1004
1005!------------------------!   
1006! >>> Write a variable   !
1007!    ... Copy&Paste part !
1008!------------------------!
1009FIELD='ST040100'
1010UNITS=''
1011DESC='ZMEA'
1012XLVL=200100.
1013SLAB=gwparam(:,:,1)
1014        ! And now put everything in the destination file
1015        ! ... Header
1016        write(1) IFV
1017        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1018        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1019        ! ... Data
1020        write(1) SLAB
1021!print *,'The field '//DESC//' was written to '//output
1022
1023!------------------------!   
1024! >>> Write a variable   !
1025!    ... Copy&Paste part !
1026!------------------------!
1027FIELD='ST100200'
1028UNITS=''
1029DESC='ZSTD'
1030XLVL=200100.
1031SLAB=gwparam(:,:,2)
1032        ! And now put everything in the destination file
1033        ! ... Header
1034        write(1) IFV
1035        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1036        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1037        ! ... Data
1038        write(1) SLAB
1039!print *,'The field '//DESC//' was written to '//output
1040
1041!------------------------!   
1042! >>> Write a variable   !
1043!    ... Copy&Paste part !
1044!------------------------!
1045FIELD='LANDSEA'
1046UNITS='0/1 Flag'
1047DESC='Land/Sea flag'
1048XLVL=200100.
1049SLAB=ones(:,:)
1050        ! And now put everything in the destination file
1051        ! ... Header
1052        write(1) IFV
1053        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1054        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1055        ! ... Data
1056        write(1) SLAB
1057!print *,'The field '//DESC//' was written to '//output
1058
1059!------------------------!   
1060! >>> Write a variable   !
1061!    ... Copy&Paste part !
1062!------------------------!
1063FIELD='SKINTEMP'
1064UNITS='K'
1065DESC='Ground temperature'
1066XLVL=200100.
1067SLAB=tsfile(:,:,time_out(l))
1068        ! And now put everything in the destination file
1069        ! ... Header
1070        write(1) IFV
1071        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1072        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1073        ! ... Data
1074        write(1) SLAB
1075!print *,'The field '//DESC//' was written to '//output
1076
1077!------------------------!   
1078! >>> Write a variable   !
1079!    ... Copy&Paste part !
1080!------------------------!
1081FIELD='SM000010'
1082UNITS=''
1083DESC='ZSIG'
1084XLVL=200100.
1085SLAB=gwparam(:,:,3)
1086        ! And now put everything in the destination file
1087        ! ... Header
1088        write(1) IFV
1089        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1090        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1091        ! ... Data
1092        write(1) SLAB
1093!print *,'The field '//DESC//' was written to '//output
1094
1095!------------------------!   
1096! >>> Write a variable   !
1097!    ... Copy&Paste part !
1098!------------------------!
1099FIELD='SM010040'
1100UNITS=''
1101DESC='ZGAM'
1102XLVL=200100.
1103SLAB=gwparam(:,:,4)
1104        ! And now put everything in the destination file
1105        ! ... Header
1106        write(1) IFV
1107        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1108        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1109        ! ... Data
1110        write(1) SLAB
1111!print *,'The field '//DESC//' was written to '//output
1112
1113!------------------------!   
1114! >>> Write a variable   !
1115!    ... Copy&Paste part !
1116!------------------------!
1117FIELD='SM040100'
1118UNITS=''
1119DESC='ZTHE'
1120XLVL=200100.
1121SLAB=gwparam(:,:,5)
1122        ! And now put everything in the destination file
1123        ! ... Header
1124        write(1) IFV
1125        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1126        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1127        ! ... Data
1128        write(1) SLAB
1129!print *,'The field '//DESC//' was written to '//output
1130
1131!------------------------!   
1132! >>> Write a variable   !
1133!    ... Copy&Paste part !
1134!------------------------!
1135FIELD='SM100200'
1136UNITS='kg/m2'
1137DESC='Surf water ice'
1138XLVL=200100.
1139SLAB=swatericefile(:,:,time_out(l))
1140        ! And now put everything in the destination file
1141        ! ... Header
1142        write(1) IFV
1143        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1144        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1145        ! ... Data
1146        write(1) SLAB
1147!print *,'The field '//DESC//' was written to '//output
1148
1149
1150!------------------------!
1151! >>> Write a variable   !
1152!    ... Copy&Paste part !
1153!------------------------!
1154FIELD='HV'
1155UNITS='kg/kg'
1156DESC='Water vapor'
1157XLVL=200100.
1158SLAB=waterfile(:,:,1,time_out(l))
1159        ! And now put everything in the destination file
1160        ! ... Header
1161        write(1) IFV
1162        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1163        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1164        ! ... Data
1165        write(1) SLAB
1166!print *,'The field '//DESC//' was written to '//output
1167
1168!------------------------!
1169! >>> Write a variable   !
1170!    ... Copy&Paste part !
1171!------------------------!
1172FIELD='HI'
1173UNITS='kg/kg'
1174DESC='Water ice'
1175XLVL=200100.
1176SLAB=watericefile(:,:,1,time_out(l))
1177        ! And now put everything in the destination file
1178        ! ... Header
1179        write(1) IFV
1180        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1181        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1182        ! ... Data
1183        write(1) SLAB
1184!print *,'The field '//DESC//' was written to '//output
1185
1186!------------------------!
1187! >>> Write a variable   !
1188!    ... Copy&Paste part !
1189!------------------------!
1190FIELD='TSOIL'
1191UNITS='K'
1192DESC='Soil temperature'
1193XLVL=200100.
1194SLAB=tsoilfile(:,:,1,time_out(l))
1195        ! And now put everything in the destination file
1196        ! ... Header
1197        write(1) IFV
1198        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1199        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1200        ! ... Data
1201        write(1) SLAB
1202!print *,'The field '//DESC//' was written to '//output
1203
1204!------------------------!
1205! >>> Write a variable   !
1206!    ... Copy&Paste part !
1207!------------------------!
1208FIELD='DSOIL'
1209UNITS='m'
1210DESC='Soil depths'
1211XLVL=200100.
1212SLAB=dsoilfile(:,:,1,time_out(l))
1213        ! And now put everything in the destination file
1214        ! ... Header
1215        write(1) IFV
1216        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1217        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1218        ! ... Data
1219        write(1) SLAB
1220!print *,'The field '//DESC//' was written to '//output
1221
1222!------------------------!
1223! >>> Write a variable   !
1224!    ... Copy&Paste part !
1225!------------------------!
1226FIELD='ISOIL'
1227UNITS='tiu'
1228DESC='Soil thermal inertia'
1229XLVL=200100.
1230SLAB=isoilfile(:,:,1,time_out(l))
1231        ! And now put everything in the destination file
1232        ! ... Header
1233        write(1) IFV
1234        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1235        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1236        ! ... Data
1237        write(1) SLAB
1238!print *,'The field '//DESC//' was written to '//output
1239
1240!------------------------!
1241! >>> Write a variable   !
1242!    ... Copy&Paste part !
1243!------------------------!
1244FIELD='CO2'
1245UNITS='kg/kg'
1246DESC='CO2 mixing ratio'
1247XLVL=200100.
1248SLAB=co2file(:,:,1,time_out(l))
1249        ! And now put everything in the destination file
1250        ! ... Header
1251        write(1) IFV
1252        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1253        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1254        ! ... Data
1255        write(1) SLAB
1256!print *,'The field '//DESC//' was written to '//output
1257
1258!------------------------!
1259! >>> Write a variable   !
1260!    ... Copy&Paste part !
1261!------------------------!
1262FIELD='DUSTQ'
1263UNITS='kg/kg'
1264DESC='Dust mixing ratio'
1265XLVL=200100.
1266SLAB=dustfile(:,:,1,time_out(l))
1267        ! And now put everything in the destination file
1268        ! ... Header
1269        write(1) IFV
1270        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1271        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1272        ! ... Data
1273        write(1) SLAB
1274!print *,'The field '//DESC//' was written to '//output
1275
1276!------------------------!
1277! >>> Write a variable   !
1278!    ... Copy&Paste part !
1279!------------------------!
1280FIELD='DUSTN'
1281UNITS='part/kg'
1282DESC='Dust number density'
1283XLVL=200100.
1284SLAB=dustnfile(:,:,1,time_out(l))
1285        ! And now put everything in the destination file
1286        ! ... Header
1287        write(1) IFV
1288        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1289        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1290        ! ... Data
1291        write(1) SLAB
1292!print *,'The field '//DESC//' was written to '//output
1293
1294!------------------------!
1295! >>> Write a variable   !
1296!    ... Copy&Paste part !
1297!------------------------!
1298FIELD='CCNQ'
1299UNITS='kg/kg'
1300DESC='CCN mixing ratio'
1301XLVL=200100.
1302SLAB=ccnfile(:,:,1,time_out(l))
1303        ! And now put everything in the destination file
1304        ! ... Header
1305        write(1) IFV
1306        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1307        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1308        ! ... Data
1309        write(1) SLAB
1310!print *,'The field '//DESC//' was written to '//output
1311
1312!------------------------!
1313! >>> Write a variable   !
1314!    ... Copy&Paste part !
1315!------------------------!
1316FIELD='CCNN'
1317UNITS='part/kg'
1318DESC='CCN number density'
1319XLVL=200100.
1320SLAB=ccnnfile(:,:,1,time_out(l))
1321        ! And now put everything in the destination file
1322        ! ... Header
1323        write(1) IFV
1324        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1325        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1326        ! ... Data
1327        write(1) SLAB
1328!print *,'The field '//DESC//' was written to '//output
1329
1330
1331!------------------------!
1332! >>> Write a variable   !
1333!     PHOTOCHEMISTRY     !
1334!------------------------!
1335#ifdef PHOTOCHEM
1336    DO i=1,nchemtrac
1337       FIELD=wtnom(i)
1338       UNITS='units'
1339       DESC='desc'
1340       XLVL=200100.
1341       SLAB=chemtrac(:,:,1,time_out(l),i)
1342       ! And now put everything in the destination file
1343       ! ... Header
1344       write(1) IFV
1345       write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1346       write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1347       ! ... Data
1348       write(1) SLAB
1349    ENDDO
1350#endif
1351
1352!!----------------------------------------------------
1353!! Write file in intermediate format for WPS
1354!! 2. 3D meteorological data
1355!! >>> same stuff, but handle vertical profile
1356!! >>> !! XLVL must be decreasing (pressure levels)
1357!!----------------------------------------------------
1358
1359!------------------------!   
1360! >>> Write a variable   !
1361!    ... Copy&Paste part !
1362!------------------------!
1363FIELD='T'
1364UNITS='K'
1365DESC='Atmospheric temperature'
1366DO k = 1,altlen
1367        XLVL=levels(k)
1368        SLAB=tfile(:,:,k,time_out(l))
1369                ! And now put everything in the destination file
1370                ! ... Header
1371        write(1) IFV
1372        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1373        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1374                ! ... Data
1375                write(1) SLAB
1376END DO
1377!print *,'The field '//DESC//' was written to '//output
1378
1379!------------------------!   
1380! >>> Write a variable   !
1381!    ... Copy&Paste part !
1382!------------------------!
1383FIELD='U'
1384UNITS='m s{-1}'
1385DESC='Zonal wind'
1386DO k = 1,altlen
1387        XLVL=levels(k)
1388        SLAB=ufile(:,:,k,time_out(l))
1389                ! And now put everything in the destination file
1390                ! ... Header
1391        write(1) IFV
1392        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1393        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1394                ! ... Data
1395                write(1) SLAB
1396END DO
1397!print *,'The field '//DESC//' was written to '//output
1398
1399!------------------------!   
1400! >>> Write a variable   !
1401!    ... Copy&Paste part !
1402!------------------------!
1403FIELD='V'
1404UNITS='m s{-1}'
1405DESC='Meridional wind'
1406DO k = 1,altlen
1407        XLVL=levels(k)
1408        SLAB=vfile(:,:,k,time_out(l))
1409                ! And now put everything in the destination file
1410                ! ... Header
1411        write(1) IFV
1412        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1413        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1414                ! ... Data
1415                write(1) SLAB
1416END DO
1417!print *,'The field '//DESC//' was written to '//output
1418
1419!------------------------!   
1420! >>> Write a variable   !
1421!    ... Copy&Paste part !
1422!------------------------!
1423FIELD='RH'
1424UNITS=''
1425DESC='Customized 2D static field'
1426DESC='Eta-levels from the GCM'
1427DESC='Pressure values from the GCM'
1428DO k = 1,altlen
1429        XLVL=levels(k)
1430SELECT CASE(ident)
1431CASE('LMD')
1432        SLAB=aps(k)+bps(k)*psfile(:,:,time_out(l))
1433CASE('OXF')
1434        SLAB=bps(k)*psfile(:,:,time_out(l))
1435END SELECT
1436                ! And now put everything in the destination file
1437                ! ... Header
1438        write(1) IFV
1439        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1440        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1441                ! ... Data
1442                write(1) SLAB
1443END DO
1444!print *,'The field '//DESC//' was written to '//output
1445
1446!------------------------!   
1447! >>> Write a variable   !
1448!    ... Copy&Paste part !
1449!------------------------!
1450FIELD='HGT'
1451UNITS='m'
1452DESC='Height'
1453DO k = 1,altlen
1454        XLVL=levels(k)
1455!!*******
1456!! PROVISOIRE: normalement, il faudrait la hauteur geopotentielle
1457!!*******
1458!!however, not used by initialize_real
1459SELECT CASE(ident)
1460CASE('LMD')
1461        SLAB=10000.*alog(610./(aps(k)+bps(k)*psfile(:,:,time_out(l))))
1462CASE('OXF')
1463        SLAB=10000.*alog(610./(bps(k)*psfile(:,:,time_out(l))))
1464END SELECT
1465                ! And now put everything in the destination file
1466                ! ... Header
1467        write(1) IFV
1468        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1469        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1470                ! ... Data
1471                write(1) SLAB
1472END DO
1473!print *,'The field '//DESC//' was written to '//output
1474
1475!------------------------!
1476! >>> Write a variable   !
1477!    ... Copy&Paste part !
1478!------------------------!
1479FIELD='HV'
1480UNITS='kg/kg'
1481DESC='Water vapor'
1482DO k = 1,altlen
1483        XLVL=levels(k)
1484        SLAB=waterfile(:,:,k,time_out(l))
1485                ! And now put everything in the destination file
1486                ! ... Header
1487        write(1) IFV
1488        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1489        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1490                ! ... Data
1491                write(1) SLAB
1492END DO
1493!print *,'The field '//DESC//' was written to '//output
1494
1495!------------------------!
1496! >>> Write a variable   !
1497!    ... Copy&Paste part !
1498!------------------------!
1499FIELD='HI'
1500UNITS='kg/kg'
1501DESC='Water ice'
1502DO k = 1,altlen
1503        XLVL=levels(k)
1504        SLAB=watericefile(:,:,k,time_out(l))
1505                ! And now put everything in the destination file
1506                ! ... Header
1507        write(1) IFV
1508        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1509        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1510                ! ... Data
1511                write(1) SLAB
1512END DO
1513!print *,'The field '//DESC//' was written to '//output
1514
1515!------------------------!
1516! >>> Write a variable   !
1517!    ... Copy&Paste part !
1518!------------------------!
1519FIELD='TSOIL'
1520UNITS='K'
1521DESC='Soil temperature'
1522DO k = 1,altlen
1523        XLVL=levels(k)
1524        SLAB=tsoilfile(:,:,k,time_out(l))
1525                ! And now put everything in the destination file
1526                ! ... Header
1527        write(1) IFV
1528        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1529        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1530                ! ... Data
1531                write(1) SLAB
1532END DO
1533!print *,'The field '//DESC//' was written to '//output
1534
1535!------------------------!
1536! >>> Write a variable   !
1537!    ... Copy&Paste part !
1538!------------------------!
1539FIELD='DSOIL'
1540UNITS='m'
1541DESC='Soil depths'
1542DO k = 1,altlen
1543        XLVL=levels(k)
1544        SLAB=dsoilfile(:,:,k,time_out(l))
1545                ! And now put everything in the destination file
1546                ! ... Header
1547        write(1) IFV
1548        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1549        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1550                ! ... Data
1551                write(1) SLAB
1552END DO
1553!print *,'The field '//DESC//' was written to '//output
1554
1555!------------------------!
1556! >>> Write a variable   !
1557!    ... Copy&Paste part !
1558!------------------------!
1559FIELD='ISOIL'
1560UNITS='tiu'
1561DESC='Soil thermal inertia'
1562DO k = 1,altlen
1563        XLVL=levels(k)
1564        SLAB=isoilfile(:,:,k,time_out(l))
1565                ! And now put everything in the destination file
1566                ! ... Header
1567        write(1) IFV
1568        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1569        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1570                ! ... Data
1571                write(1) SLAB
1572END DO
1573!print *,'The field '//DESC//' was written to '//output
1574
1575!------------------------!
1576! >>> Write a variable   !
1577!    ... Copy&Paste part !
1578!------------------------!
1579FIELD='CO2'
1580UNITS='kg/kg'
1581DESC='CO2 mixing ratio'
1582DO k = 1,altlen
1583        XLVL=levels(k)
1584        SLAB=co2file(:,:,k,time_out(l))
1585                ! And now put everything in the destination file
1586                ! ... Header
1587        write(1) IFV
1588        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1589        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1590                ! ... Data
1591                write(1) SLAB
1592END DO
1593!print *,'The field '//DESC//' was written to '//output
1594
1595!------------------------!
1596! >>> Write a variable   !
1597!    ... Copy&Paste part !
1598!------------------------!
1599FIELD='DUSTQ'
1600UNITS='kg/kg'
1601DESC='Dust mixing ratio'
1602DO k = 1,altlen
1603        XLVL=levels(k)
1604        SLAB=dustfile(:,:,k,time_out(l))
1605                ! And now put everything in the destination file
1606                ! ... Header
1607        write(1) IFV
1608        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1609        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1610                ! ... Data
1611                write(1) SLAB
1612END DO
1613!print *,'The field '//DESC//' was written to '//output
1614
1615!------------------------!
1616! >>> Write a variable   !
1617!    ... Copy&Paste part !
1618!------------------------!
1619FIELD='DUSTN'
1620UNITS='part/kg'
1621DESC='Dust number density'
1622DO k = 1,altlen
1623        XLVL=levels(k)
1624        SLAB=dustnfile(:,:,k,time_out(l))
1625                ! And now put everything in the destination file
1626                ! ... Header
1627        write(1) IFV
1628        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1629        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1630                ! ... Data
1631                write(1) SLAB
1632END DO
1633!print *,'The field '//DESC//' was written to '//output
1634
1635!------------------------!
1636! >>> Write a variable   !
1637!    ... Copy&Paste part !
1638!------------------------!
1639FIELD='CCNQ'
1640UNITS='kg/kg'
1641DESC='CCN mixing ratio'
1642DO k = 1,altlen
1643        XLVL=levels(k)
1644        SLAB=ccnfile(:,:,k,time_out(l))
1645                ! And now put everything in the destination file
1646                ! ... Header
1647        write(1) IFV
1648        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1649        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1650                ! ... Data
1651                write(1) SLAB
1652END DO
1653!print *,'The field '//DESC//' was written to '//output
1654
1655!------------------------!
1656! >>> Write a variable   !
1657!    ... Copy&Paste part !
1658!------------------------!
1659FIELD='CCNN'
1660UNITS='part/kg'
1661DESC='CCN number density'
1662DO k = 1,altlen
1663        XLVL=levels(k)
1664        SLAB=ccnnfile(:,:,k,time_out(l))
1665                ! And now put everything in the destination file
1666                ! ... Header
1667        write(1) IFV
1668        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1669        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1670                ! ... Data
1671                write(1) SLAB
1672END DO
1673!print *,'The field '//DESC//' was written to '//output
1674
1675
1676!------------------------!
1677! >>> Write a variable   !
1678!     PHOTOCHEMISTRY     !
1679!------------------------!
1680#ifdef PHOTOCHEM
1681    DO i=1,nchemtrac
1682       FIELD=wtnom(i)
1683       UNITS='units'
1684       DESC='desc'
1685       DO k = 1,altlen
1686         XLVL=levels(k)
1687         SLAB=chemtrac(:,:,k,time_out(l),i)
1688         ! And now put everything in the destination file
1689         ! ... Header
1690         write(1) IFV
1691         write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1692         write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1693         ! ... Data
1694         write(1) SLAB
1695       END DO
1696    ENDDO
1697#endif
1698
1699print *,'****done file '//output, int(100.*float(l)/float(FILES)), ' % '
1700close(1)
1701
1702DEALLOCATE(SLAB)
1703
1704END DO !! end of the files loop
1705
1706
1707!!-------------------------
1708!! DeAllocate local arrays
1709!!-------------------------
1710deallocate(eta_gcm)
1711deallocate(tfile)
1712deallocate(tsoilfile)
1713deallocate(isoilfile)
1714deallocate(dsoilfile)
1715!deallocate(tfileorig)
1716deallocate(ufile)
1717!deallocate(ufileorig)
1718deallocate(vfile)
1719!deallocate(vfileorig)
1720deallocate(rfile)
1721deallocate(hfile) 
1722deallocate(waterfile)
1723deallocate(watericefile)
1724!deallocate(swaterfile)
1725deallocate(swatericefile)
1726deallocate(dustfile)
1727deallocate(dustnfile)
1728deallocate(co2file)
1729deallocate(ccnfile)
1730deallocate(ccnnfile)
1731deallocate(psfile)
1732deallocate(tsfile)
1733deallocate(tnsfile)
1734deallocate(unsfile)
1735deallocate(vnsfile)
1736deallocate(emissfile)
1737deallocate(co2icefile)
1738deallocate(ghtsfile)    !! no scan axis
1739deallocate(vide)
1740deallocate(ones)
1741deallocate(lat, lon, alt, time)
1742deallocate(aps,bps,levels)
1743
1744#ifdef PHOTOCHEM
1745deallocate(chemtrac)
1746deallocate(wtnom)
1747#endif
1748
1749print *, '------------------------'
1750print *, 'End of pre-WPS process !'
1751print *, '------------------------'
1752print *, 'Here is what you should set in namelist.wps:'
1753print *, " start_date = '"//date_out2(1)//"'"
1754print *, " end_date = '"//date_out2(FILES)//"'"
1755print *, '------------------------'
1756print *, 'Any date between those is ok'
1757print *, 'for example, second available date is ... '//date_out2(2)
1758print *, '------------------------'
1759
1760end
Note: See TracBrowser for help on using the repository browser.