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

Last change on this file since 779 was 665, checked in by aslmd, 13 years ago

MESOSCALE : output more physical variables in Registry. update GCM for storm
simulations. bug fix in readmeteo with ccn. UTIL PYTHON : bug fix for
operation cat (time axis)

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