source: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/Generic/readmeteo.F90

Last change on this file was 2070, checked in by mlefevre, 6 years ago

Add necessary files to create start state for generic mesoscale model
File size: 48.5 KB
RevLine 
[2070]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=10.5
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, "ap", 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, "bp", 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,"h2o_vap",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,"h2o_ice",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      co2file(:,:,:,:)=0.0
648    else
649      ierr=NF_GET_VAR_REAL(nid,nvarid,co2file)
650    endif
651
652!!------------------------
653!! special dust stuff
654!!------------------------
655    print *,'Dust mass'
656    ierr=NF_INQ_VARID(nid,"dustq",nvarid)
657    if (ierr.ne.NF_NOERR) then
658      write(*,*) "...No dustq - Dust mass set to 0"
659      dustfile(:,:,:,:)=0.
660    else
661      ierr=NF_GET_VAR_REAL(nid,nvarid,dustfile)
662    endif
663
664    print *,'Dust number'
665    ierr=NF_INQ_VARID(nid,"dustN",nvarid)
666    if (ierr.ne.NF_NOERR) then
667      write(*,*) "...No dustN - Dust number set to 0"
668      dustnfile(:,:,:,:)=0.
669    else
670      ierr=NF_GET_VAR_REAL(nid,nvarid,dustnfile)
671    endif
672
673    print *,'CCN mass'
674    ierr=NF_INQ_VARID(nid,"ccn",nvarid)
675    if (ierr.ne.NF_NOERR) then
676      write(*,*) "...No ccn - CCN mass set to 0"
677      ccnfile(:,:,:,:)=0.
678    else
679      ierr=NF_GET_VAR_REAL(nid,nvarid,ccnfile)
680    endif
681
682    print *,'CCN number'
683    ierr=NF_INQ_VARID(nid,"ccnN",nvarid)
684    if (ierr.ne.NF_NOERR) then
685      write(*,*) "...No ccnN - CCN number set to 0"
686      ccnnfile(:,:,:,:)=0.
687    else
688      ierr=NF_GET_VAR_REAL(nid,nvarid,ccnnfile)
689    endif
690!!------------------------
691!! special dust stuff
692!!------------------------
693
694!SELECT CASE(ident)
695!CASE('LMD')
696
697    print *,'Soil temperatures'
698    ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
699    if (ierr.ne.NF_NOERR) then
700        write(*,*) "...No tsoil - Soil temperatures set isothermal with surface temperature"
701        DO l=1,altlen
702                tsoilfile(:,:,l,:)=tsfile(:,:,:)
703        ENDDO
704    else
705        ierr=NF_GET_VAR_REAL(nid,nvarid,tsoilfile)
706    endif
707
708!!!!!!!!
709!!!!!!!! new physics (but still compatible with old physics)
710    print *,'Soil depths'
711    ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
712    if (ierr.ne.NF_NOERR) then
713        write(*,*) "...No soildepth - Set to -999"  !!! see soil_settings in LMD physics
714        DO l=1,altlen
715                vertdsoil(l)=-999.
716        ENDDO
717    else
718        ierr=NF_GET_VAR_REAL(nid,nvarid,vertdsoil)
719    endif
720    print *, 'wait a minute' !! AS: I know this could be better
721    DO m=1,lonlen
722     DO n=1,latlen
723      DO p=1,timelen
724       dsoilfile(m,n,:,p) = vertdsoil(:)
725      ENDDO
726     ENDDO
727    ENDDO
728    DEALLOCATE(vertdsoil)
729
730    print *,'Soil thermal inertia'
731    ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
732    if (ierr.ne.NF_NOERR) then
733        write(*,*) "...No soil therm. inert. - Set to -999"
734        DO l=1,altlen
735                isoilfile(:,:,l,:)=-999.
736        ENDDO
737    else
738        ierr=NF_GET_VAR_REAL(nid,nvarid,isoilfile)
739    endif
740!!!!!!!!
741!!!!!!!! new physics
742
743
744!!!!!!!!!!!!!!!!!!!!!!!!NEW PHYSICS + PHOTOCHEM
745!!!!!!!!!!!!!!!!!!!!!!!!NEW PHYSICS + PHOTOCHEM
746#ifdef PHOTOCHEM
747    print *,'photochem'
748    DO i=1,nchemtrac
749     print *,wtnom(i)
750     ierr=NF_INQ_VARID(nid,wtnom(i),nvarid)
751     if (ierr.ne.NF_NOERR) then
752       write(*,*) "...No ",wtnom(i), " - set to 0"
753       chemtrac(:,:,:,:,i)=0.
754     else
755       ierr=NF_GET_VAR_REAL(nid,nvarid,chemtrac(:,:,:,:,i))
756     endif
757    ENDDO
758#endif
759!!!!!!!!!!!!!!!!!!!!!!!!NEW PHYSICS + PHOTOCHEM
760!!!!!!!!!!!!!!!!!!!!!!!!NEW PHYSICS + PHOTOCHEM
761
762
763
764!END SELECT
765
766ierr=NF_CLOSE(nid)
767
768!!!!lott stuff
769!print *,'get lott'
770!ierr=NF_OPEN("init_lott.nc",NF_NOWRITE,nid3)
771!ierr=NF_INQ_VARID(nid3,"ZMEA",nvarid)
772!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
773!gwparam(:,:,1)=interm(:,:)
774!ierr=NF_INQ_VARID(nid3,"ZSTD",nvarid)
775!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
776!gwparam(:,:,2)=interm(:,:)
777!ierr=NF_INQ_VARID(nid3,"ZSIG",nvarid)
778!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
779!gwparam(:,:,3)=interm(:,:)
780!ierr=NF_INQ_VARID(nid3,"ZGAM",nvarid)
781!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
782!gwparam(:,:,4)=interm(:,:)
783!ierr=NF_INQ_VARID(nid3,"ZTHE",nvarid)
784!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
785!gwparam(:,:,5)=interm(:,:)
786!ierr=NF_CLOSE(nid3)
787
788ENDIF
789
790!!-----------------------------
791!! Loop on the written files
792!! >>> what follows is pretty repetitive ...
793!!        gotta do something (one more loop?)
794!!-----------------------------
795
796!!! Equivalent eta levels for WRF interpolation in initialize_real
797!print *,'Computing equivalent eta levels'
798!       
799!       !ptop will be passed through RH(surface)
800!       ptop=MINVAL(aps(altlen)+bps(altlen)*psfile(:,:,:))
801!       print *, 'ptop', ptop
802!
803!print *, 'sample: equivalent eta levels at i=1,j=1'
804!DO k = 1,altlen
805!        levels(k)=-k
806!        eta_gcm(:,:,k,:)=(aps(k)+bps(k)*psfile(:,:,:)-ptop)/(psfile(:,:,:)-ptop)
807!        print *, eta_gcm(1,1,k,1)
808!END DO
809
810!!---better to pass pressure values
811!!---the eta treatment is done in module_initialize_real
812DO k = 1,altlen
813        levels(k)=-k
814        !!dummy decreasing levels
815END DO
816
817
818
819
820!! Dummy surface level is XLVL=200100.
821
822
823DO l=1,FILES
824
825
826!!---------------------------------------------
827!! Write file in intermediate format for WPS
828!! 1. Surface data
829!!---------------------------------------------
830!!
831!! a mettre pour tous sinon
832!!     WRF_DEBUG: Warning DIM             4 , NAME num_metgrid_levels REDIFINED  by
833!!            var DUSTN            26           25  in wrf_io.F90 line         2349
834!!
835
836!
837! These variables remain the same in the "loop"
838!
839HDATE=date_out(l)
840IFV=4
841XFCST=0.
842SOURCE=ident
843NX=lonlen
844NY=latlen
845ALLOCATE(SLAB(NX,NY))
846IPROJ=0
847STARTLOC='SWCORNER'
848STARTLAT=lat(1)
849STARTLON=lon(1)
850DELTALAT=lat(latlen)-lat(latlen-1)
851DELTALON=lon(lonlen)-lon(lonlen-1)
852!
853! Open the data destination file
854!
855output="./WPSFEED/"//ident//":"//date_out2(l)       
856open(UNIT=1,                    &
857        FILE=output,            &
858        STATUS='REPLACE',       &
859        FORM='UNFORMATTED',     &
860        ACCESS='SEQUENTIAL')
861
862!------------------------!   
863! >>> Write a variable   !
864!    ... Copy&Paste part !
865!------------------------!
866FIELD='T'
867UNITS='K'
868DESC='Atmospheric temperature'
869XLVL=200100.
870SLAB=tnsfile(:,:,time_out(l))
871        ! And now put everything in the destination file
872        ! ... Header
873        write(1) IFV
874        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
875        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
876        ! ... Data
877        write(1) SLAB
878!print *,'The field '//DESC//' was written to '//output
879
880!------------------------!   
881! >>> Write a variable   !
882!    ... Copy&Paste part !
883!------------------------!
884FIELD='U'
885UNITS='m s{-1}'
886DESC='Zonal wind'
887XLVL=200100.
888SLAB=unsfile(:,:,time_out(l))
889        ! And now put everything in the destination file
890        ! ... Header
891        write(1) IFV
892        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
893        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
894        ! ... Data
895        write(1) SLAB
896!print *,'The field '//DESC//' was written to '//output
897
898!------------------------!   
899! >>> Write a variable   !
900!    ... Copy&Paste part !
901!------------------------!
902FIELD='V'
903UNITS='m s{-1}'
904DESC='Meridional wind'
905XLVL=200100.
906SLAB=vnsfile(:,:,time_out(l))
907        ! And now put everything in the destination file
908        ! ... Header
909        write(1) IFV
910        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
911        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
912        ! ... Data
913        write(1) SLAB
914!print *,'The field '//DESC//' was written to '//output
915
916!------------------------!   
917! >>> Write a variable   !
918!    ... Copy&Paste part !
919!------------------------!
920FIELD='RH'
921UNITS=''
922DESC='Customized 2D static field'
923XLVL=200100.
924SLAB=vide(:,:)
925        ! And now put everything in the destination file
926        ! ... Header
927        write(1) IFV
928        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
929        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
930        ! ... Data
931        write(1) SLAB
932!print *,'The field '//DESC//' was written to '//output
933
934!------------------------!   
935! >>> Write a variable   !
936!    ... Copy&Paste part !
937!------------------------!
938FIELD='SOILHGT'
939UNITS='m'
940DESC='Terrain field of source analysis'  !Ground geopotential height
941XLVL=200100.
942SLAB=ghtsfile(:,:)
943        ! And now put everything in the destination file
944        ! ... Header
945        write(1) IFV
946        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
947        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
948        ! ... Data
949        write(1) SLAB
950!print *,'The field '//DESC//' was written to '//output
951
952!------------------------!   
953! >>> Write a variable   !
954!    ... Copy&Paste part !
955!------------------------!
956FIELD='PSFC'
957UNITS='Pa'
958DESC='Surface Pressure'
959XLVL=200100.
960SLAB=psfile(:,:,time_out(l))
961        ! And now put everything in the destination file
962        ! ... Header
963        write(1) IFV
964        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
965        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
966        ! ... Data
967        write(1) SLAB
968!print *,'The field '//DESC//' was written to '//output
969
970!------------------------!   
971! >>> Write a variable   !
972!    ... Copy&Paste part !
973!------------------------!
974FIELD='ST000010'
975UNITS=''
976DESC='Emissivity'
977XLVL=200100.
978SLAB=emissfile(:,:,time_out(l))
979        ! And now put everything in the destination file
980        ! ... Header
981        write(1) IFV
982        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
983        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
984        ! ... Data
985        write(1) SLAB
986!print *,'The field '//DESC//' was written to '//output
987
988!------------------------!   
989! >>> Write a variable   !
990!    ... Copy&Paste part !
991!------------------------!
992FIELD='ST010040'
993UNITS=''
994DESC='CO2 ice'
995XLVL=200100.
996SLAB=co2icefile(:,:,time_out(l))
997        ! And now put everything in the destination file
998        ! ... Header
999        write(1) IFV
1000        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1001        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1002        ! ... Data
1003        write(1) SLAB
1004!print *,'The field '//DESC//' was written to '//output
1005
1006!------------------------!   
1007! >>> Write a variable   !
1008!    ... Copy&Paste part !
1009!------------------------!
1010FIELD='ST040100'
1011UNITS=''
1012DESC='ZMEA'
1013XLVL=200100.
1014SLAB=gwparam(:,:,1)
1015        ! And now put everything in the destination file
1016        ! ... Header
1017        write(1) IFV
1018        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1019        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1020        ! ... Data
1021        write(1) SLAB
1022!print *,'The field '//DESC//' was written to '//output
1023
1024!------------------------!   
1025! >>> Write a variable   !
1026!    ... Copy&Paste part !
1027!------------------------!
1028FIELD='ST100200'
1029UNITS=''
1030DESC='ZSTD'
1031XLVL=200100.
1032SLAB=gwparam(:,:,2)
1033        ! And now put everything in the destination file
1034        ! ... Header
1035        write(1) IFV
1036        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1037        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1038        ! ... Data
1039        write(1) SLAB
1040!print *,'The field '//DESC//' was written to '//output
1041
1042!------------------------!   
1043! >>> Write a variable   !
1044!    ... Copy&Paste part !
1045!------------------------!
1046FIELD='LANDSEA'
1047UNITS='0/1 Flag'
1048DESC='Land/Sea flag'
1049XLVL=200100.
1050SLAB=ones(:,:)
1051        ! And now put everything in the destination file
1052        ! ... Header
1053        write(1) IFV
1054        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1055        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1056        ! ... Data
1057        write(1) SLAB
1058!print *,'The field '//DESC//' was written to '//output
1059
1060!------------------------!   
1061! >>> Write a variable   !
1062!    ... Copy&Paste part !
1063!------------------------!
1064FIELD='SKINTEMP'
1065UNITS='K'
1066DESC='Ground temperature'
1067XLVL=200100.
1068SLAB=tsfile(:,:,time_out(l))
1069        ! And now put everything in the destination file
1070        ! ... Header
1071        write(1) IFV
1072        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1073        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1074        ! ... Data
1075        write(1) SLAB
1076!print *,'The field '//DESC//' was written to '//output
1077
1078!------------------------!   
1079! >>> Write a variable   !
1080!    ... Copy&Paste part !
1081!------------------------!
1082FIELD='SM000010'
1083UNITS=''
1084DESC='ZSIG'
1085XLVL=200100.
1086SLAB=gwparam(:,:,3)
1087        ! And now put everything in the destination file
1088        ! ... Header
1089        write(1) IFV
1090        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1091        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1092        ! ... Data
1093        write(1) SLAB
1094!print *,'The field '//DESC//' was written to '//output
1095
1096!------------------------!   
1097! >>> Write a variable   !
1098!    ... Copy&Paste part !
1099!------------------------!
1100FIELD='SM010040'
1101UNITS=''
1102DESC='ZGAM'
1103XLVL=200100.
1104SLAB=gwparam(:,:,4)
1105        ! And now put everything in the destination file
1106        ! ... Header
1107        write(1) IFV
1108        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1109        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1110        ! ... Data
1111        write(1) SLAB
1112!print *,'The field '//DESC//' was written to '//output
1113
1114!------------------------!   
1115! >>> Write a variable   !
1116!    ... Copy&Paste part !
1117!------------------------!
1118FIELD='SM040100'
1119UNITS=''
1120DESC='ZTHE'
1121XLVL=200100.
1122SLAB=gwparam(:,:,5)
1123        ! And now put everything in the destination file
1124        ! ... Header
1125        write(1) IFV
1126        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1127        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1128        ! ... Data
1129        write(1) SLAB
1130!print *,'The field '//DESC//' was written to '//output
1131
1132!------------------------!   
1133! >>> Write a variable   !
1134!    ... Copy&Paste part !
1135!------------------------!
1136FIELD='SM100200'
1137UNITS='kg/m2'
1138DESC='Surf water ice'
1139XLVL=200100.
1140SLAB=swatericefile(:,:,time_out(l))
1141        ! And now put everything in the destination file
1142        ! ... Header
1143        write(1) IFV
1144        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1145        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1146        ! ... Data
1147        write(1) SLAB
1148!print *,'The field '//DESC//' was written to '//output
1149
1150
1151!------------------------!
1152! >>> Write a variable   !
1153!    ... Copy&Paste part !
1154!------------------------!
1155FIELD='HV'
1156UNITS='kg/kg'
1157DESC='Water vapor'
1158XLVL=200100.
1159SLAB=waterfile(:,:,1,time_out(l))
1160        ! And now put everything in the destination file
1161        ! ... Header
1162        write(1) IFV
1163        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1164        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1165        ! ... Data
1166        write(1) SLAB
1167!print *,'The field '//DESC//' was written to '//output
1168
1169!------------------------!
1170! >>> Write a variable   !
1171!    ... Copy&Paste part !
1172!------------------------!
1173FIELD='HI'
1174UNITS='kg/kg'
1175DESC='Water ice'
1176XLVL=200100.
1177SLAB=watericefile(:,:,1,time_out(l))
1178        ! And now put everything in the destination file
1179        ! ... Header
1180        write(1) IFV
1181        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1182        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1183        ! ... Data
1184        write(1) SLAB
1185!print *,'The field '//DESC//' was written to '//output
1186
1187!------------------------!
1188! >>> Write a variable   !
1189!    ... Copy&Paste part !
1190!------------------------!
1191FIELD='TSOIL'
1192UNITS='K'
1193DESC='Soil temperature'
1194XLVL=200100.
1195SLAB=tsoilfile(:,:,1,time_out(l))
1196        ! And now put everything in the destination file
1197        ! ... Header
1198        write(1) IFV
1199        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1200        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1201        ! ... Data
1202        write(1) SLAB
1203!print *,'The field '//DESC//' was written to '//output
1204
1205!------------------------!
1206! >>> Write a variable   !
1207!    ... Copy&Paste part !
1208!------------------------!
1209FIELD='DSOIL'
1210UNITS='m'
1211DESC='Soil depths'
1212XLVL=200100.
1213SLAB=dsoilfile(:,:,1,time_out(l))
1214        ! And now put everything in the destination file
1215        ! ... Header
1216        write(1) IFV
1217        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1218        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1219        ! ... Data
1220        write(1) SLAB
1221!print *,'The field '//DESC//' was written to '//output
1222
1223!------------------------!
1224! >>> Write a variable   !
1225!    ... Copy&Paste part !
1226!------------------------!
1227FIELD='ISOIL'
1228UNITS='tiu'
1229DESC='Soil thermal inertia'
1230XLVL=200100.
1231SLAB=isoilfile(:,:,1,time_out(l))
1232        ! And now put everything in the destination file
1233        ! ... Header
1234        write(1) IFV
1235        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1236        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1237        ! ... Data
1238        write(1) SLAB
1239!print *,'The field '//DESC//' was written to '//output
1240
1241!------------------------!
1242! >>> Write a variable   !
1243!    ... Copy&Paste part !
1244!------------------------!
1245FIELD='CO2'
1246UNITS='kg/kg'
1247DESC='CO2 mixing ratio'
1248XLVL=200100.
1249SLAB=co2file(:,:,1,time_out(l))
1250        ! And now put everything in the destination file
1251        ! ... Header
1252        write(1) IFV
1253        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1254        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1255        ! ... Data
1256        write(1) SLAB
1257!print *,'The field '//DESC//' was written to '//output
1258
1259!------------------------!
1260! >>> Write a variable   !
1261!    ... Copy&Paste part !
1262!------------------------!
1263FIELD='DUSTQ'
1264UNITS='kg/kg'
1265DESC='Dust mixing ratio'
1266XLVL=200100.
1267SLAB=dustfile(:,:,1,time_out(l))
1268        ! And now put everything in the destination file
1269        ! ... Header
1270        write(1) IFV
1271        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1272        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1273        ! ... Data
1274        write(1) SLAB
1275!print *,'The field '//DESC//' was written to '//output
1276
1277!------------------------!
1278! >>> Write a variable   !
1279!    ... Copy&Paste part !
1280!------------------------!
1281FIELD='DUSTN'
1282UNITS='part/kg'
1283DESC='Dust number density'
1284XLVL=200100.
1285SLAB=dustnfile(:,:,1,time_out(l))
1286        ! And now put everything in the destination file
1287        ! ... Header
1288        write(1) IFV
1289        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1290        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1291        ! ... Data
1292        write(1) SLAB
1293!print *,'The field '//DESC//' was written to '//output
1294
1295!------------------------!
1296! >>> Write a variable   !
1297!    ... Copy&Paste part !
1298!------------------------!
1299FIELD='CCNQ'
1300UNITS='kg/kg'
1301DESC='CCN mixing ratio'
1302XLVL=200100.
1303SLAB=ccnfile(:,:,1,time_out(l))
1304        ! And now put everything in the destination file
1305        ! ... Header
1306        write(1) IFV
1307        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1308        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1309        ! ... Data
1310        write(1) SLAB
1311!print *,'The field '//DESC//' was written to '//output
1312
1313!------------------------!
1314! >>> Write a variable   !
1315!    ... Copy&Paste part !
1316!------------------------!
1317FIELD='CCNN'
1318UNITS='part/kg'
1319DESC='CCN number density'
1320XLVL=200100.
1321SLAB=ccnnfile(:,:,1,time_out(l))
1322        ! And now put everything in the destination file
1323        ! ... Header
1324        write(1) IFV
1325        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1326        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1327        ! ... Data
1328        write(1) SLAB
1329!print *,'The field '//DESC//' was written to '//output
1330
1331
1332!------------------------!
1333! >>> Write a variable   !
1334!     PHOTOCHEMISTRY     !
1335!------------------------!
1336#ifdef PHOTOCHEM
1337    DO i=1,nchemtrac
1338       FIELD=wtnom(i)
1339       UNITS='units'
1340       DESC='desc'
1341       XLVL=200100.
1342       SLAB=chemtrac(:,:,1,time_out(l),i)
1343       ! And now put everything in the destination file
1344       ! ... Header
1345       write(1) IFV
1346       write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1347       write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1348       ! ... Data
1349       write(1) SLAB
1350    ENDDO
1351#endif
1352
1353!!----------------------------------------------------
1354!! Write file in intermediate format for WPS
1355!! 2. 3D meteorological data
1356!! >>> same stuff, but handle vertical profile
1357!! >>> !! XLVL must be decreasing (pressure levels)
1358!!----------------------------------------------------
1359
1360!------------------------!   
1361! >>> Write a variable   !
1362!    ... Copy&Paste part !
1363!------------------------!
1364FIELD='T'
1365UNITS='K'
1366DESC='Atmospheric temperature'
1367DO k = 1,altlen
1368        XLVL=levels(k)
1369        SLAB=tfile(:,:,k,time_out(l))
1370                ! And now put everything in the destination file
1371                ! ... Header
1372        write(1) IFV
1373        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1374        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1375                ! ... Data
1376                write(1) SLAB
1377END DO
1378!print *,'The field '//DESC//' was written to '//output
1379
1380!------------------------!   
1381! >>> Write a variable   !
1382!    ... Copy&Paste part !
1383!------------------------!
1384FIELD='U'
1385UNITS='m s{-1}'
1386DESC='Zonal wind'
1387DO k = 1,altlen
1388        XLVL=levels(k)
1389        SLAB=ufile(:,:,k,time_out(l))
1390                ! And now put everything in the destination file
1391                ! ... Header
1392        write(1) IFV
1393        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1394        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1395                ! ... Data
1396                write(1) SLAB
1397END DO
1398!print *,'The field '//DESC//' was written to '//output
1399
1400!------------------------!   
1401! >>> Write a variable   !
1402!    ... Copy&Paste part !
1403!------------------------!
1404FIELD='V'
1405UNITS='m s{-1}'
1406DESC='Meridional wind'
1407DO k = 1,altlen
1408        XLVL=levels(k)
1409        SLAB=vfile(:,:,k,time_out(l))
1410                ! And now put everything in the destination file
1411                ! ... Header
1412        write(1) IFV
1413        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1414        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1415                ! ... Data
1416                write(1) SLAB
1417END DO
1418!print *,'The field '//DESC//' was written to '//output
1419
1420!------------------------!   
1421! >>> Write a variable   !
1422!    ... Copy&Paste part !
1423!------------------------!
1424FIELD='RH'
1425UNITS=''
1426DESC='Customized 2D static field'
1427DESC='Eta-levels from the GCM'
1428DESC='Pressure values from the GCM'
1429DO k = 1,altlen
1430        XLVL=levels(k)
1431SELECT CASE(ident)
1432CASE('LMD')
1433        SLAB=aps(k)+bps(k)*psfile(:,:,time_out(l))
1434CASE('OXF')
1435        SLAB=bps(k)*psfile(:,:,time_out(l))
1436END SELECT
1437                ! And now put everything in the destination file
1438                ! ... Header
1439        write(1) IFV
1440        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1441        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1442                ! ... Data
1443                write(1) SLAB
1444END DO
1445!print *,'The field '//DESC//' was written to '//output
1446
1447!------------------------!   
1448! >>> Write a variable   !
1449!    ... Copy&Paste part !
1450!------------------------!
1451FIELD='HGT'
1452UNITS='m'
1453DESC='Height'
1454DO k = 1,altlen
1455        XLVL=levels(k)
1456!!*******
1457!! PROVISOIRE: normalement, il faudrait la hauteur geopotentielle
1458!!*******
1459!!however, not used by initialize_real
1460SELECT CASE(ident)
1461CASE('LMD')
1462        SLAB=10000.*alog(610./(aps(k)+bps(k)*psfile(:,:,time_out(l))))
1463CASE('OXF')
1464        SLAB=10000.*alog(610./(bps(k)*psfile(:,:,time_out(l))))
1465END SELECT
1466                ! And now put everything in the destination file
1467                ! ... Header
1468        write(1) IFV
1469        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1470        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1471                ! ... Data
1472                write(1) SLAB
1473END DO
1474!print *,'The field '//DESC//' was written to '//output
1475
1476!------------------------!
1477! >>> Write a variable   !
1478!    ... Copy&Paste part !
1479!------------------------!
1480FIELD='HV'
1481UNITS='kg/kg'
1482DESC='Water vapor'
1483DO k = 1,altlen
1484        XLVL=levels(k)
1485        SLAB=waterfile(:,:,k,time_out(l))
1486                ! And now put everything in the destination file
1487                ! ... Header
1488        write(1) IFV
1489        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1490        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1491                ! ... Data
1492                write(1) SLAB
1493END DO
1494!print *,'The field '//DESC//' was written to '//output
1495
1496!------------------------!
1497! >>> Write a variable   !
1498!    ... Copy&Paste part !
1499!------------------------!
1500FIELD='HI'
1501UNITS='kg/kg'
1502DESC='Water ice'
1503DO k = 1,altlen
1504        XLVL=levels(k)
1505        SLAB=watericefile(:,:,k,time_out(l))
1506                ! And now put everything in the destination file
1507                ! ... Header
1508        write(1) IFV
1509        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1510        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1511                ! ... Data
1512                write(1) SLAB
1513END DO
1514!print *,'The field '//DESC//' was written to '//output
1515
1516!------------------------!
1517! >>> Write a variable   !
1518!    ... Copy&Paste part !
1519!------------------------!
1520FIELD='TSOIL'
1521UNITS='K'
1522DESC='Soil temperature'
1523DO k = 1,altlen
1524        XLVL=levels(k)
1525        SLAB=tsoilfile(:,:,k,time_out(l))
1526                ! And now put everything in the destination file
1527                ! ... Header
1528        write(1) IFV
1529        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1530        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1531                ! ... Data
1532                write(1) SLAB
1533END DO
1534!print *,'The field '//DESC//' was written to '//output
1535
1536!------------------------!
1537! >>> Write a variable   !
1538!    ... Copy&Paste part !
1539!------------------------!
1540FIELD='DSOIL'
1541UNITS='m'
1542DESC='Soil depths'
1543DO k = 1,altlen
1544        XLVL=levels(k)
1545        SLAB=dsoilfile(:,:,k,time_out(l))
1546                ! And now put everything in the destination file
1547                ! ... Header
1548        write(1) IFV
1549        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1550        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1551                ! ... Data
1552                write(1) SLAB
1553END DO
1554!print *,'The field '//DESC//' was written to '//output
1555
1556!------------------------!
1557! >>> Write a variable   !
1558!    ... Copy&Paste part !
1559!------------------------!
1560FIELD='ISOIL'
1561UNITS='tiu'
1562DESC='Soil thermal inertia'
1563DO k = 1,altlen
1564        XLVL=levels(k)
1565        SLAB=isoilfile(:,:,k,time_out(l))
1566                ! And now put everything in the destination file
1567                ! ... Header
1568        write(1) IFV
1569        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1570        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1571                ! ... Data
1572                write(1) SLAB
1573END DO
1574!print *,'The field '//DESC//' was written to '//output
1575
1576!------------------------!
1577! >>> Write a variable   !
1578!    ... Copy&Paste part !
1579!------------------------!
1580FIELD='CO2'
1581UNITS='kg/kg'
1582DESC='CO2 mixing ratio'
1583DO k = 1,altlen
1584        XLVL=levels(k)
1585        SLAB=co2file(:,:,k,time_out(l))
1586                ! And now put everything in the destination file
1587                ! ... Header
1588        write(1) IFV
1589        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1590        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1591                ! ... Data
1592                write(1) SLAB
1593END DO
1594!print *,'The field '//DESC//' was written to '//output
1595
1596!------------------------!
1597! >>> Write a variable   !
1598!    ... Copy&Paste part !
1599!------------------------!
1600FIELD='DUSTQ'
1601UNITS='kg/kg'
1602DESC='Dust mixing ratio'
1603DO k = 1,altlen
1604        XLVL=levels(k)
1605        SLAB=dustfile(:,:,k,time_out(l))
1606                ! And now put everything in the destination file
1607                ! ... Header
1608        write(1) IFV
1609        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1610        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1611                ! ... Data
1612                write(1) SLAB
1613END DO
1614!print *,'The field '//DESC//' was written to '//output
1615
1616!------------------------!
1617! >>> Write a variable   !
1618!    ... Copy&Paste part !
1619!------------------------!
1620FIELD='DUSTN'
1621UNITS='part/kg'
1622DESC='Dust number density'
1623DO k = 1,altlen
1624        XLVL=levels(k)
1625        SLAB=dustnfile(:,:,k,time_out(l))
1626                ! And now put everything in the destination file
1627                ! ... Header
1628        write(1) IFV
1629        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1630        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1631                ! ... Data
1632                write(1) SLAB
1633END DO
1634!print *,'The field '//DESC//' was written to '//output
1635
1636!------------------------!
1637! >>> Write a variable   !
1638!    ... Copy&Paste part !
1639!------------------------!
1640FIELD='CCNQ'
1641UNITS='kg/kg'
1642DESC='CCN mixing ratio'
1643DO k = 1,altlen
1644        XLVL=levels(k)
1645        SLAB=ccnfile(:,:,k,time_out(l))
1646                ! And now put everything in the destination file
1647                ! ... Header
1648        write(1) IFV
1649        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1650        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1651                ! ... Data
1652                write(1) SLAB
1653END DO
1654!print *,'The field '//DESC//' was written to '//output
1655
1656!------------------------!
1657! >>> Write a variable   !
1658!    ... Copy&Paste part !
1659!------------------------!
1660FIELD='CCNN'
1661UNITS='part/kg'
1662DESC='CCN number density'
1663DO k = 1,altlen
1664        XLVL=levels(k)
1665        SLAB=ccnnfile(:,:,k,time_out(l))
1666                ! And now put everything in the destination file
1667                ! ... Header
1668        write(1) IFV
1669        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1670        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1671                ! ... Data
1672                write(1) SLAB
1673END DO
1674!print *,'The field '//DESC//' was written to '//output
1675
1676
1677!------------------------!
1678! >>> Write a variable   !
1679!     PHOTOCHEMISTRY     !
1680!------------------------!
1681#ifdef PHOTOCHEM
1682    DO i=1,nchemtrac
1683       FIELD=wtnom(i)
1684       UNITS='units'
1685       DESC='desc'
1686       DO k = 1,altlen
1687         XLVL=levels(k)
1688         SLAB=chemtrac(:,:,k,time_out(l),i)
1689         ! And now put everything in the destination file
1690         ! ... Header
1691         write(1) IFV
1692         write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1693         write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1694         ! ... Data
1695         write(1) SLAB
1696       END DO
1697    ENDDO
1698#endif
1699
1700print *,'****done file '//output, int(100.*float(l)/float(FILES)), ' % '
1701close(1)
1702
1703DEALLOCATE(SLAB)
1704
1705END DO !! end of the files loop
1706
1707
1708!!-------------------------
1709!! DeAllocate local arrays
1710!!-------------------------
1711deallocate(eta_gcm)
1712deallocate(tfile)
1713deallocate(tsoilfile)
1714deallocate(isoilfile)
1715deallocate(dsoilfile)
1716!deallocate(tfileorig)
1717deallocate(ufile)
1718!deallocate(ufileorig)
1719deallocate(vfile)
1720!deallocate(vfileorig)
1721deallocate(rfile)
1722deallocate(hfile) 
1723deallocate(waterfile)
1724deallocate(watericefile)
1725!deallocate(swaterfile)
1726deallocate(swatericefile)
1727deallocate(dustfile)
1728deallocate(dustnfile)
1729deallocate(co2file)
1730deallocate(ccnfile)
1731deallocate(ccnnfile)
1732deallocate(psfile)
1733deallocate(tsfile)
1734deallocate(tnsfile)
1735deallocate(unsfile)
1736deallocate(vnsfile)
1737deallocate(emissfile)
1738deallocate(co2icefile)
1739deallocate(ghtsfile)    !! no scan axis
1740deallocate(vide)
1741deallocate(ones)
1742deallocate(lat, lon, alt, time)
1743deallocate(aps,bps,levels)
1744
1745#ifdef PHOTOCHEM
1746deallocate(chemtrac)
1747deallocate(wtnom)
1748#endif
1749
1750print *, '------------------------'
1751print *, 'End of pre-WPS process !'
1752print *, '------------------------'
1753print *, 'Here is what you should set in namelist.wps:'
1754print *, " start_date = '"//date_out2(1)//"'"
1755print *, " end_date = '"//date_out2(FILES)//"'"
1756print *, '------------------------'
1757print *, 'Any date between those is ok'
1758print *, 'for example, second available date is ... '//date_out2(2)
1759print *, '------------------------'
1760
1761end
Note: See TracBrowser for help on using the repository browser.