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

Last change on this file since 433 was 324, checked in by aslmd, 14 years ago

MESOSCALE : Preparatory commit for the ultimate option mars=42 which

would allow mesoscale modeling with photochemistry.

[see r76 method to add 'mars' options]
Modified module_lmd_driver.F and Registry.EM and runmeso
Modified readmeteo.F90 and introduced an option -DPHOTOCHEM
...
Transparent to the casual user
Option mars=42 not yet finished though -- so do not use!
File size: 44.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
90
91!! Reading the parameter file
92integer :: tmp,increment,FILES
93integer :: tmp2,tmp3,tmp4
94character*1 :: answer
95character*4 :: str
96character*2 :: str2, str3, str4
97integer, dimension(:), allocatable :: time_out
98character*13, dimension(:), allocatable :: date_out
99character*19, dimension(:), allocatable :: date_out2
100
101#ifdef PHOTOCHEM
102real, dimension(:,:,:,:,:), allocatable :: chemtrac
103integer :: nchemtrac,i
104CHARACTER*20,DIMENSION(:),ALLOCATABLE :: wtnom
105#endif
106
107
108!***************************************************************************
109!***************************************************************************
110
111
112!!---------------------
113!! Open the datafile
114!!---------------------
115ierr=NF_OPEN ("input_diagfi.nc",NF_NOWRITE,nid)
116IF (ierr.NE.NF_NOERR) THEN
117   write(*,*)'**** Please create a symbolic link called input_diagfi.nc'
118   CALL ABORT
119ENDIF
120
121
122!!--------------------------
123!! Ask for data references
124!!--------------------------
125!write(*,*) "Create n files. What is n ?"
126read(*,*) FILES
127allocate(time_out(FILES))
128allocate(date_out(FILES))
129allocate(date_out2(FILES))
130
131!write(*,*) ""
132!write(*,*) "INPUT: Diagfi time"
133!write(*,*) "list of n subscripts, separated by <Return>s"
134increment=1
135do while (increment.ne.FILES+1)
136    read(*,*) tmp
137    time_out(increment)=tmp
138    increment=increment+1
139enddo
140
141!write(*,*) ""
142!write(*,*) "OUTPUT: WRF time"
143!write(*,*) "fill 4 lines indicating"
144!write(*,*) "-year (4 digits)"
145!write(*,*) "-month (2 digits)"
146!write(*,*) "-day (2 digits)"
147!write(*,*) "-hour (2 digits)"
148increment=1
149do while (increment.ne.FILES+1)
150    read(*,*) str
151    read(*,*) str2
152    read(*,*) str3
153    read(*,*) str4
154    date_out(increment)=str//'-'//str2//'-'//str3//'_'//str4
155    date_out2(increment)=str//'-'//str2//'-'//str3//'_'//str4//':00:00'
156    !print *,date_out(increment)
157    !write(*,*) "ok? (y/n)"
158    read(*,*) answer
159    if (answer.eq.'n') then
160        !write(*,*) "please write again"
161    else
162        !write(*,*) "next one, please"   
163        increment=increment+1
164    endif
165enddo
166
167
168!!-------------------
169!! Get array sizes
170!!-------------------
171SELECT CASE(ident)
172CASE('LMD')
173ierr=NF_INQ_DIMID(nid,"Time",timedim)
174CASE('OXF')
175ierr=NF_INQ_DIMID(nid,"time",timedim)
176END SELECT
177ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
178  write(*,*) "timelen: ",timelen
179
180SELECT CASE(ident)
181CASE('LMD')
182ierr=NF_INQ_DIMID(nid,"latitude",latdim)
183CASE('OXF')
184ierr=NF_INQ_DIMID(nid,"lat",latdim)
185END SELECT
186ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
187  write(*,*) "latlen: ",latlen
188
189SELECT CASE(ident)
190CASE('LMD')
191ierr=NF_INQ_DIMID(nid,"longitude",londim)
192CASE('OXF')
193ierr=NF_INQ_DIMID(nid,"lon",londim)
194END SELECT
195ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
196  write(*,*) "lonlen: ",lonlen
197
198SELECT CASE(ident)
199CASE('LMD')
200ierr=NF_INQ_DIMID(nid,"altitude",altdim)
201CASE('OXF')
202ierr=NF_INQ_DIMID(nid,"sigma",altdim)
203END SELECT
204ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
205  write(*,*) "altlen: ",altlen
206
207
208
209!!-------------------------
210!! Allocate local arrays
211!!-------------------------
212allocate(eta_gcm(lonlen,latlen,altlen,timelen))
213allocate(tfile(lonlen,latlen,altlen,timelen))
214allocate(tsoilfile(lonlen,latlen,altlen,timelen))
215allocate(dsoilfile(lonlen,latlen,altlen,timelen))
216allocate(isoilfile(lonlen,latlen,altlen,timelen))
217allocate(vertdsoil(altlen))
218!allocate(tfileorig(lonlen,latlen,altlen,timelen))
219allocate(ufile(lonlen,latlen,altlen,timelen))
220!allocate(ufileorig(lonlen,latlen,altlen,timelen))
221allocate(vfile(lonlen,latlen,altlen,timelen))
222!allocate(vfileorig(lonlen,latlen,altlen,timelen))
223allocate(rfile(lonlen,latlen,altlen,timelen))
224allocate(hfile(lonlen,latlen,altlen,timelen)) 
225allocate(waterfile(lonlen,latlen,altlen,timelen))
226allocate(watericefile(lonlen,latlen,altlen,timelen))
227!allocate(swaterfile(lonlen,latlen,timelen))
228allocate(swatericefile(lonlen,latlen,timelen))
229allocate(dustfile(lonlen,latlen,altlen,timelen))
230allocate(dustnfile(lonlen,latlen,altlen,timelen))
231allocate(psfile(lonlen,latlen,timelen))
232allocate(tsfile(lonlen,latlen,timelen))
233allocate(tnsfile(lonlen,latlen,timelen))
234allocate(unsfile(lonlen,latlen,timelen))
235allocate(vnsfile(lonlen,latlen,timelen))
236allocate(emissfile(lonlen,latlen,timelen))
237allocate(co2icefile(lonlen,latlen,timelen))
238allocate(interm(lonlen,latlen))
239allocate(gwparam(lonlen,latlen,5))
240allocate(ghtsfile(lonlen,latlen))    !! no scan axis
241allocate(vide(lonlen,latlen))
242allocate(ones(lonlen,latlen))
243allocate(lat(latlen), lon(lonlen), alt(altlen), time(timelen))
244allocate(aps(altlen),bps(altlen),levels(altlen))
245#ifdef PHOTOCHEM
246nchemtrac = 14
247allocate(wtnom(nchemtrac))
248wtnom(1)  = "c_co2"
249wtnom(2)  = "c_co"
250wtnom(3)  = "c_o"
251wtnom(4)  = "c_o1d"
252wtnom(5)  = "c_o2"
253wtnom(6)  = "c_o3"
254wtnom(7)  = "c_h"
255wtnom(8)  = "c_h2"
256wtnom(9)  = "c_oh"
257wtnom(10) = "c_ho2"
258wtnom(11) = "c_h2o2"
259wtnom(12) = "c_ch4"
260wtnom(13) = "c_n2"
261wtnom(14) = "c_ar"
262allocate(chemtrac(lonlen,latlen,altlen,timelen,nchemtrac))
263chemtrac(:,:,:,:,:)=0
264#endif
265
266tfile(:,:,:,:)=0
267tsoilfile(:,:,:,:)=0
268isoilfile(:,:,:,:)=0
269dsoilfile(:,:,:,:)=0
270vertdsoil(:)=0.
271!tfileorig(:,:,:,:)=0
272!ufileorig(:,:,:,:)=0
273!vfileorig(:,:,:,:)=0
274ufile(:,:,:,:)=0
275vfile(:,:,:,:)=0
276rfile(:,:,:,:)=0
277hfile(:,:,:,:)=0
278waterfile(:,:,:,:)=0
279watericefile(:,:,:,:)=0
280!swaterfile(:,:,:)=0
281swatericefile(:,:,:)=0
282dustfile(:,:,:,:)=0
283dustnfile(:,:,:,:)=0
284psfile(:,:,:)=0
285tsfile(:,:,:)=0
286emissfile(:,:,:)=0
287co2icefile(:,:,:)=0
288tnsfile(:,:,:)=0
289unsfile(:,:,:)=0
290vnsfile(:,:,:)=0
291interm(:,:)=0
292gwparam(:,:,:)=0
293ghtsfile(:,:)=0
294vide(:,:)=0
295ones(:,:)=0
296
297
298!!-------------------
299!! Read dimensions
300!!-------------------
301
302print *,'>>> Read dimensions !'
303
304print *,'Time'
305SELECT CASE(ident)
306CASE('LMD')
307   ierr = NF_INQ_VARID (nid, "Time",nvarid)
308CASE('OXF')
309   ierr = NF_INQ_VARID (nid, "time",nvarid)
310END SELECT
311   IF (ierr .NE. NF_NOERR) THEN
312      PRINT *, "Error: Readmeteo <Time> not found"
313      stop
314   ENDIF
315#ifdef NC_DOUBLE
316   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
317#else
318   ierr = NF_GET_VAR_REAL(nid, nvarid, time)
319#endif
320   print *,time(1),' ... to ... ',time(timelen)
321
322print *,'Latitude'
323SELECT CASE(ident)
324CASE('LMD')
325   ierr = NF_INQ_VARID (nid, "latitude",nvarid)
326CASE('OXF')
327   ierr = NF_INQ_VARID (nid, "lat",nvarid)
328END SELECT
329   IF (ierr .NE. NF_NOERR) THEN
330      PRINT *, "Error: Readmeteo <latitude> not found"
331      stop
332   ENDIF
333#ifdef NC_DOUBLE
334   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lat)
335#else
336   ierr = NF_GET_VAR_REAL(nid, nvarid, lat)
337#endif
338   print *,lat(1),' ... to ... ',lat(latlen),' ... step: ',lat(latlen)-lat(latlen-1)
339
340print *,'Longitude'
341SELECT CASE(ident)
342CASE('LMD')
343   ierr = NF_INQ_VARID (nid, "longitude",nvarid)
344CASE('OXF')
345   ierr = NF_INQ_VARID (nid, "lon",nvarid)
346END SELECT
347   IF (ierr .NE. NF_NOERR) THEN
348      PRINT *, "Error: Readmeteo <longitude> not found"
349      stop
350   ENDIF
351#ifdef NC_DOUBLE
352   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lon)
353#else
354   ierr = NF_GET_VAR_REAL(nid, nvarid, lon)
355#endif
356   print *,lon(1),' ... to ... ',lon(lonlen),' ... step: ',lon(lonlen)-lon(lonlen-1)
357
358SELECT CASE(ident)
359CASE('LMD')
360print *, 'Hybrid coordinates'
361   ierr = NF_INQ_VARID (nid, "aps", nvarid)
362   IF (ierr .NE. NF_NOERR) THEN
363      PRINT *, "Error: Readmeteo <aps> not found"
364      stop
365   ENDIF
366#ifdef NC_DOUBLE
367   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, aps)
368#else
369   ierr = NF_GET_VAR_REAL(nid, nvarid, aps)
370#endif
371   ierr = NF_INQ_VARID (nid, "bps", nvarid)
372   IF (ierr .NE. NF_NOERR) THEN
373      PRINT *, "Error: Readmeteo <bps> not found"
374      stop
375   ENDIF
376#ifdef NC_DOUBLE
377   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bps)
378#else
379   ierr = NF_GET_VAR_REAL(nid, nvarid, bps)
380#endif
381   print *,aps(1),' ... to ... ',aps(altlen)
382   print *,bps(1),' ... to ... ',bps(altlen)
383CASE('OXF')
384   ierr = NF_INQ_VARID (nid, "sigma", nvarid)
385   IF (ierr .NE. NF_NOERR) THEN
386      PRINT *, "Error: Readmeteo <sigma> not found"
387      stop
388   ENDIF
389#ifdef NC_DOUBLE
390   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bps)
391#else
392   ierr = NF_GET_VAR_REAL(nid, nvarid, bps)
393#endif
394   print *,bps(1),' ... to ... ',bps(altlen)
395END SELECT
396
397
398!!-------------------------------------------
399!! Reading 2D and 3D meteorological fields
400!!-------------------------------------------
401
402IF (blank .EQV. .false.) THEN
403
404print *,'>>> Read 2D optional fields !'
405
406print *,'Emissivity'
407   ierr = NF_INQ_VARID (nid,"emis",nvarid)
408IF (ierr .NE. NF_NOERR) THEN
409        PRINT *, '...warning: not found in diagfi !'
410        PRINT *, '...will be filled with a prescribed value', emiss_prescribed
411  emissfile(:,:,:)=emiss_prescribed
412ELSE
413#ifdef NC_DOUBLE
414  ierr = NF_GET_VAR_DOUBLE(nid, nvarid, emissfile)
415#else
416  ierr = NF_GET_VAR_REAL(nid, nvarid, emissfile)
417#endif
418ENDIF   
419
420print *,'CO2 ice'
421   ierr = NF_INQ_VARID (nid,"co2ice",nvarid)
422IF (ierr .NE. NF_NOERR) THEN
423  PRINT *, '...warning: not found in diagfi !'
424  PRINT *, '...will be filled with a prescribed value', co2ice_prescribed
425  co2icefile(:,:,:)=co2ice_prescribed
426ELSE
427#ifdef NC_DOUBLE
428  ierr = NF_GET_VAR_DOUBLE(nid, nvarid, co2icefile)
429#else
430  ierr = NF_GET_VAR_REAL(nid, nvarid, co2icefile)
431#endif
432ENDIF
433
434print *,'>>> Read 2D surface fields !'
435
436print *,'Surface Pressure'
437   ierr = NF_INQ_VARID (nid,"ps",nvarid)
438   IF (ierr .NE. NF_NOERR) THEN
439      PRINT *, "Error: Readmeteo <ps> not found"
440      stop
441   ENDIF
442#ifdef NC_DOUBLE
443   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, psfile)
444#else
445   ierr = NF_GET_VAR_REAL(nid, nvarid, psfile)
446#endif
447
448print *,'Ground Temperature'
449   ierr = NF_INQ_VARID (nid,"tsurf",nvarid)
450   IF (ierr .NE. NF_NOERR) THEN
451      PRINT *, "Error: Readmeteo <tsurf> not found"
452      stop
453   ENDIF
454#ifdef NC_DOUBLE
455   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsfile)
456#else
457   ierr = NF_GET_VAR_REAL(nid, nvarid, tsfile)
458#endif
459
460
461!!!"atmospheric" surface temperature is taken
462!!!from original diagfi.nc first level
463!!!... level is ~3-5 meters
464!print *,'Near-Surface Temperature'
465!   ierr = NF_INQ_VARID (nid,"temp",nvarid)
466!   IF (ierr .NE. NF_NOERR) THEN
467!      ierr = NF_INQ_VARID (nid,"t",nvarid)
468!       IF (ierr .NE. NF_NOERR) THEN
469!           PRINT *, "Error: Readmeteo <temp> not found"
470!           stop
471!       ENDIF
472!   ENDIF
473!#ifdef NC_DOUBLE
474!   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfileorig)
475!#else
476!   ierr = NF_GET_VAR_REAL(nid, nvarid, tfileorig)
477!#endif
478!   tnsfile=tfileorig(:,:,1,:)
479
480!!!"atmospheric" surface u is taken
481!!!from original diagfi.nc first level
482!!!... level is ~3-5 meters
483!print *,'Near-Surface Zonal Wind'
484!   ierr = NF_INQ_VARID (nid,"u",nvarid)
485!   IF (ierr .NE. NF_NOERR) THEN
486!     PRINT *, "Error: Readmeteo <u> not found"
487!     stop
488!   ENDIF
489!#ifdef NC_DOUBLE
490!   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufileorig)
491!#else
492!   ierr = NF_GET_VAR_REAL(nid, nvarid, ufileorig)
493!#endif
494!   unsfile=ufileorig(:,:,1,:)
495
496!!!"atmospheric" surface v is taken
497!!!from original diagfi.nc first level
498!!!... level is ~3-5 meters
499!print *,'Near-Surface Meridional Wind'
500!   ierr = NF_INQ_VARID (nid,"v",nvarid)
501!   IF (ierr .NE. NF_NOERR) THEN
502!     PRINT *, "Error: Readmeteo <v> not found"
503!     stop
504!   ENDIF
505!#ifdef NC_DOUBLE
506!   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfileorig)
507!#else
508!   ierr = NF_GET_VAR_REAL(nid, nvarid, vfileorig)
509!#endif
510!   vnsfile=vfileorig(:,:,1,:)
511
512SELECT CASE(ident)
513CASE('LMD')
514   print *,'Geopotential height at the ground'
515   ierr = NF_INQ_VARID (nid,"phisinit",nvarid)
516   IF (ierr .NE. NF_NOERR) THEN
517      PRINT *, "Error: Readmeteo <phisinit> not found"
518      stop
519   ENDIF
520#ifdef NC_DOUBLE
521   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ghtsfile)
522#else
523   ierr = NF_GET_VAR_REAL(nid, nvarid, ghtsfile)
524#endif
525!**** from geopotential to geopotential height
526ghtsfile=ghtsfile/grav
527!****
528CASE('OXF')
529!
530! geopotential height ~ altimetry
531!
532print *,'Geopotential height at the ground from file mountain_new.nc'
533ierr=NF_OPEN("mountain_new.nc",NF_NOWRITE,nid3)
534    if (ierr.ne.NF_NOERR) then
535      write(*,*) "Error: Could not open that file either"
536      stop "Might as well stop here"
537    endif
538ierr=NF_INQ_VARID(nid3,"orography",nvarid)
539ierr=NF_GET_VAR_REAL(nid3,nvarid,ghtsfile)
540if (ierr.ne.NF_NOERR) then
541    stop "Error: Failed reading phisinit"
542endif
543ierr=NF_CLOSE(nid3)
544END SELECT
545
546
547print *,'>>> Read 3D meteorological fields ! - This may take some time ...'
548
549print *,'Temperature'
550   ierr = NF_INQ_VARID (nid,"temp",nvarid)
551   IF (ierr .NE. NF_NOERR) THEN
552      ierr = NF_INQ_VARID (nid,"t",nvarid)
553        IF (ierr .NE. NF_NOERR) THEN
554          PRINT *, "Error: Readmeteo <t> not found"
555          stop
556        ENDIF
557   ENDIF
558#ifdef NC_DOUBLE
559   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfile)
560#else
561   ierr = NF_GET_VAR_REAL(nid, nvarid, tfile)
562#endif
563tnsfile=tfile(:,:,1,:)
564
565print *,'Zonal wind'   
566   ierr = NF_INQ_VARID (nid,"u",nvarid)
567   IF (ierr .NE. NF_NOERR) THEN
568      PRINT *, "Error: Readmeteo <u> not found"
569      stop
570   ENDIF
571#ifdef NC_DOUBLE
572   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufile)
573#else
574   ierr = NF_GET_VAR_REAL(nid, nvarid, ufile)
575#endif
576unsfile=ufile(:,:,1,:)
577
578print *,'Meridional wind'
579   ierr = NF_INQ_VARID (nid,"v",nvarid)
580   IF (ierr .NE. NF_NOERR) THEN
581      PRINT *, "Error: Readmeteo <v> not found"
582      stop
583   ENDIF
584#ifdef NC_DOUBLE
585   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfile)
586#else
587   ierr = NF_GET_VAR_REAL(nid, nvarid, vfile)
588#endif
589vnsfile=ufile(:,:,1,:)
590
591
592!!------------------------
593!! special water stuff
594!!------------------------
595    print *,'Water vapor'
596    ierr=NF_INQ_VARID(nid,"q02",nvarid)
597    if (ierr.ne.NF_NOERR) then
598      write(*,*) "...No q02 - Water vapor set to 0"
599      waterfile(:,:,:,:)=0. 
600    else
601      ierr=NF_GET_VAR_REAL(nid,nvarid,waterfile)
602    endif   
603
604    print *,'Water ice'
605    ierr=NF_INQ_VARID(nid,"q01",nvarid)
606    if (ierr.ne.NF_NOERR) then
607      write(*,*) "...No q01 - Water ice set to 0" 
608      watericefile(:,:,:,:)=0.
609    else
610      ierr=NF_GET_VAR_REAL(nid,nvarid,watericefile)
611    endif
612!    print *,'Surface Water vapor'
613!    ierr=NF_INQ_VARID(nid,"qsurf02",nvarid)
614!    if (ierr.ne.NF_NOERR) then
615!      write(*,*) "...No qsurf02 - surface Water vapor set to 0"
616!      swaterfile(:,:,:)=0.
617!    endif
618!    ierr=NF_GET_VAR_REAL(nid,nvarid,swaterfile)
619
620    print *,'Surface Water ice'
621!!!!!! ATTENTION ATTENTION
622!!!!!! water ice a la surface est qsurf(ig,nqmx)
623    ierr=NF_INQ_VARID(nid,"qsurf02",nvarid)
624    if (ierr.ne.NF_NOERR) then
625      write(*,*) "...No qsurf02 - surface Water ice set to 0"
626      swatericefile(:,:,:)=0.
627    else
628      ierr=NF_GET_VAR_REAL(nid,nvarid,swatericefile)
629    endif
630!!------------------------
631!! special water stuff
632!!------------------------
633
634!!------------------------
635!! special dust stuff
636!!------------------------
637    print *,'Dust mass'
638    ierr=NF_INQ_VARID(nid,"dustq",nvarid)
639    if (ierr.ne.NF_NOERR) then
640      write(*,*) "...No dustq - Dust mass set to 0"
641      dustfile(:,:,:,:)=0.
642    else
643      ierr=NF_GET_VAR_REAL(nid,nvarid,dustfile)
644    endif
645
646    print *,'Dust number'
647    ierr=NF_INQ_VARID(nid,"dustN",nvarid)
648    if (ierr.ne.NF_NOERR) then
649      write(*,*) "...No dustN - Dust number set to 0"
650      dustnfile(:,:,:,:)=0.
651    else
652      ierr=NF_GET_VAR_REAL(nid,nvarid,dustnfile)
653    endif
654!!------------------------
655!! special dust stuff
656!!------------------------
657
658!SELECT CASE(ident)
659!CASE('LMD')
660
661    print *,'Soil temperatures'
662    ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
663    if (ierr.ne.NF_NOERR) then
664        write(*,*) "...No tsoil - Soil temperatures set isothermal with surface temperature"
665        DO l=1,altlen
666                tsoilfile(:,:,l,:)=tsfile(:,:,:)
667        ENDDO
668    else
669        ierr=NF_GET_VAR_REAL(nid,nvarid,tsoilfile)
670    endif
671
672!!!!!!!!
673!!!!!!!! new physics (but still compatible with old physics)
674    print *,'Soil depths'
675    ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
676    if (ierr.ne.NF_NOERR) then
677        write(*,*) "...No soildepth - Set to -999"  !!! see soil_settings in LMD physics
678        DO l=1,altlen
679                vertdsoil(l)=-999.
680        ENDDO
681    else
682        ierr=NF_GET_VAR_REAL(nid,nvarid,vertdsoil)
683    endif
684    print *, 'wait a minute' !! AS: I know this could be better
685    DO m=1,lonlen
686     DO n=1,latlen
687      DO p=1,timelen
688       dsoilfile(m,n,:,p) = vertdsoil(:)
689      ENDDO
690     ENDDO
691    ENDDO
692    DEALLOCATE(vertdsoil)
693
694    print *,'Soil thermal inertia'
695    ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
696    if (ierr.ne.NF_NOERR) then
697        write(*,*) "...No soil therm. inert. - Set to -999"
698        DO l=1,altlen
699                isoilfile(:,:,l,:)=-999.
700        ENDDO
701    else
702        ierr=NF_GET_VAR_REAL(nid,nvarid,isoilfile)
703    endif
704!!!!!!!!
705!!!!!!!! new physics
706
707
708!!!!!!!!!!!!!!!!!!!!!!!!NEW PHYSICS + PHOTOCHEM
709!!!!!!!!!!!!!!!!!!!!!!!!NEW PHYSICS + PHOTOCHEM
710#ifdef PHOTOCHEM
711    print *,'photochem'
712    DO i=1,nchemtrac
713     print *,wtnom(i)
714     ierr=NF_INQ_VARID(nid,wtnom(i),nvarid)
715     if (ierr.ne.NF_NOERR) then
716       write(*,*) "...No ",wtnom(i), " - set to 0"
717       chemtrac(:,:,:,:,i)=0.
718     else
719       ierr=NF_GET_VAR_REAL(nid,nvarid,chemtrac(:,:,:,:,i))
720     endif
721    ENDDO
722#endif
723!!!!!!!!!!!!!!!!!!!!!!!!NEW PHYSICS + PHOTOCHEM
724!!!!!!!!!!!!!!!!!!!!!!!!NEW PHYSICS + PHOTOCHEM
725
726
727
728!END SELECT
729
730ierr=NF_CLOSE(nid)
731
732!!!!lott stuff
733!print *,'get lott'
734!ierr=NF_OPEN("init_lott.nc",NF_NOWRITE,nid3)
735!ierr=NF_INQ_VARID(nid3,"ZMEA",nvarid)
736!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
737!gwparam(:,:,1)=interm(:,:)
738!ierr=NF_INQ_VARID(nid3,"ZSTD",nvarid)
739!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
740!gwparam(:,:,2)=interm(:,:)
741!ierr=NF_INQ_VARID(nid3,"ZSIG",nvarid)
742!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
743!gwparam(:,:,3)=interm(:,:)
744!ierr=NF_INQ_VARID(nid3,"ZGAM",nvarid)
745!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
746!gwparam(:,:,4)=interm(:,:)
747!ierr=NF_INQ_VARID(nid3,"ZTHE",nvarid)
748!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
749!gwparam(:,:,5)=interm(:,:)
750!ierr=NF_CLOSE(nid3)
751
752ENDIF
753
754!!-----------------------------
755!! Loop on the written files
756!! >>> what follows is pretty repetitive ...
757!!        gotta do something (one more loop?)
758!!-----------------------------
759
760!!! Equivalent eta levels for WRF interpolation in initialize_real
761!print *,'Computing equivalent eta levels'
762!       
763!       !ptop will be passed through RH(surface)
764!       ptop=MINVAL(aps(altlen)+bps(altlen)*psfile(:,:,:))
765!       print *, 'ptop', ptop
766!
767!print *, 'sample: equivalent eta levels at i=1,j=1'
768!DO k = 1,altlen
769!        levels(k)=-k
770!        eta_gcm(:,:,k,:)=(aps(k)+bps(k)*psfile(:,:,:)-ptop)/(psfile(:,:,:)-ptop)
771!        print *, eta_gcm(1,1,k,1)
772!END DO
773
774!!---better to pass pressure values
775!!---the eta treatment is done in module_initialize_real
776DO k = 1,altlen
777        levels(k)=-k
778        !!dummy decreasing levels
779END DO
780
781
782
783
784!! Dummy surface level is XLVL=200100.
785
786
787DO l=1,FILES
788
789
790!!---------------------------------------------
791!! Write file in intermediate format for WPS
792!! 1. Surface data
793!!---------------------------------------------
794!!
795!! a mettre pour tous sinon
796!!     WRF_DEBUG: Warning DIM             4 , NAME num_metgrid_levels REDIFINED  by
797!!            var DUSTN            26           25  in wrf_io.F90 line         2349
798!!
799
800!
801! These variables remain the same in the "loop"
802!
803HDATE=date_out(l)
804IFV=4
805XFCST=0.
806SOURCE=ident
807NX=lonlen
808NY=latlen
809ALLOCATE(SLAB(NX,NY))
810IPROJ=0
811STARTLOC='SWCORNER'
812STARTLAT=lat(1)
813STARTLON=lon(1)
814DELTALAT=lat(latlen)-lat(latlen-1)
815DELTALON=lon(lonlen)-lon(lonlen-1)
816!
817! Open the data destination file
818!
819output="./WPSFEED/"//ident//":"//date_out2(l)       
820open(UNIT=1,                    &
821        FILE=output,            &
822        STATUS='REPLACE',       &
823        FORM='UNFORMATTED',     &
824        ACCESS='SEQUENTIAL')
825
826!------------------------!   
827! >>> Write a variable   !
828!    ... Copy&Paste part !
829!------------------------!
830FIELD='T'
831UNITS='K'
832DESC='Atmospheric temperature'
833XLVL=200100.
834SLAB=tnsfile(:,:,time_out(l))
835        ! And now put everything in the destination file
836        ! ... Header
837        write(1) IFV
838        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
839        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
840        ! ... Data
841        write(1) SLAB
842!print *,'The field '//DESC//' was written to '//output
843
844!------------------------!   
845! >>> Write a variable   !
846!    ... Copy&Paste part !
847!------------------------!
848FIELD='U'
849UNITS='m s{-1}'
850DESC='Zonal wind'
851XLVL=200100.
852SLAB=unsfile(:,:,time_out(l))
853        ! And now put everything in the destination file
854        ! ... Header
855        write(1) IFV
856        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
857        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
858        ! ... Data
859        write(1) SLAB
860!print *,'The field '//DESC//' was written to '//output
861
862!------------------------!   
863! >>> Write a variable   !
864!    ... Copy&Paste part !
865!------------------------!
866FIELD='V'
867UNITS='m s{-1}'
868DESC='Meridional wind'
869XLVL=200100.
870SLAB=vnsfile(:,:,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='RH'
885UNITS=''
886DESC='Customized 2D static field'
887XLVL=200100.
888SLAB=vide(:,:)
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='SOILHGT'
903UNITS='m'
904DESC='Terrain field of source analysis'  !Ground geopotential height
905XLVL=200100.
906SLAB=ghtsfile(:,:)
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='PSFC'
921UNITS='Pa'
922DESC='Surface Pressure'
923XLVL=200100.
924SLAB=psfile(:,:,time_out(l))
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='ST000010'
939UNITS=''
940DESC='Emissivity'
941XLVL=200100.
942SLAB=emissfile(:,:,time_out(l))
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='ST010040'
957UNITS=''
958DESC='CO2 ice'
959XLVL=200100.
960SLAB=co2icefile(:,:,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='ST040100'
975UNITS=''
976DESC='ZMEA'
977XLVL=200100.
978SLAB=gwparam(:,:,1)
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='ST100200'
993UNITS=''
994DESC='ZSTD'
995XLVL=200100.
996SLAB=gwparam(:,:,2)
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='LANDSEA'
1011UNITS='0/1 Flag'
1012DESC='Land/Sea flag'
1013XLVL=200100.
1014SLAB=ones(:,:)
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='SKINTEMP'
1029UNITS='K'
1030DESC='Ground temperature'
1031XLVL=200100.
1032SLAB=tsfile(:,:,time_out(l))
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='SM000010'
1047UNITS=''
1048DESC='ZSIG'
1049XLVL=200100.
1050SLAB=gwparam(:,:,3)
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='SM010040'
1065UNITS=''
1066DESC='ZGAM'
1067XLVL=200100.
1068SLAB=gwparam(:,:,4)
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='SM040100'
1083UNITS=''
1084DESC='ZTHE'
1085XLVL=200100.
1086SLAB=gwparam(:,:,5)
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='SM100200'
1101UNITS='kg/m2'
1102DESC='Surf water ice'
1103XLVL=200100.
1104SLAB=swatericefile(:,:,time_out(l))
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!------------------------!
1116! >>> Write a variable   !
1117!    ... Copy&Paste part !
1118!------------------------!
1119FIELD='HV'
1120UNITS='kg/kg'
1121DESC='Water vapor'
1122XLVL=200100.
1123SLAB=waterfile(:,:,1,time_out(l))
1124        ! And now put everything in the destination file
1125        ! ... Header
1126        write(1) IFV
1127        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1128        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1129        ! ... Data
1130        write(1) SLAB
1131!print *,'The field '//DESC//' was written to '//output
1132
1133!------------------------!
1134! >>> Write a variable   !
1135!    ... Copy&Paste part !
1136!------------------------!
1137FIELD='HI'
1138UNITS='kg/kg'
1139DESC='Water ice'
1140XLVL=200100.
1141SLAB=watericefile(:,:,1,time_out(l))
1142        ! And now put everything in the destination file
1143        ! ... Header
1144        write(1) IFV
1145        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1146        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1147        ! ... Data
1148        write(1) SLAB
1149!print *,'The field '//DESC//' was written to '//output
1150
1151!------------------------!
1152! >>> Write a variable   !
1153!    ... Copy&Paste part !
1154!------------------------!
1155FIELD='TSOIL'
1156UNITS='K'
1157DESC='Soil temperature'
1158XLVL=200100.
1159SLAB=tsoilfile(:,:,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='DSOIL'
1174UNITS='m'
1175DESC='Soil depths'
1176XLVL=200100.
1177SLAB=dsoilfile(:,:,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='ISOIL'
1192UNITS='tiu'
1193DESC='Soil thermal inertia'
1194XLVL=200100.
1195SLAB=isoilfile(:,:,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='DUSTQ'
1210UNITS='kg/kg'
1211DESC='Dust mixing ratio'
1212XLVL=200100.
1213SLAB=dustfile(:,:,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='DUSTN'
1228UNITS='part/kg'
1229DESC='Dust number density'
1230XLVL=200100.
1231SLAB=dustnfile(:,:,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!     PHOTOCHEMISTRY     !
1244!------------------------!
1245#ifdef PHOTOCHEM
1246    DO i=1,nchemtrac
1247       FIELD=wtnom(i)
1248       UNITS='units'
1249       DESC='desc'
1250       XLVL=200100.
1251       SLAB=chemtrac(:,:,1,time_out(l),i)
1252       ! And now put everything in the destination file
1253       ! ... Header
1254       write(1) IFV
1255       write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1256       write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1257       ! ... Data
1258       write(1) SLAB
1259    ENDDO
1260#endif
1261
1262!!----------------------------------------------------
1263!! Write file in intermediate format for WPS
1264!! 2. 3D meteorological data
1265!! >>> same stuff, but handle vertical profile
1266!! >>> !! XLVL must be decreasing (pressure levels)
1267!!----------------------------------------------------
1268
1269!------------------------!   
1270! >>> Write a variable   !
1271!    ... Copy&Paste part !
1272!------------------------!
1273FIELD='T'
1274UNITS='K'
1275DESC='Atmospheric temperature'
1276DO k = 1,altlen
1277        XLVL=levels(k)
1278        SLAB=tfile(:,:,k,time_out(l))
1279                ! And now put everything in the destination file
1280                ! ... Header
1281        write(1) IFV
1282        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1283        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1284                ! ... Data
1285                write(1) SLAB
1286END DO
1287!print *,'The field '//DESC//' was written to '//output
1288
1289!------------------------!   
1290! >>> Write a variable   !
1291!    ... Copy&Paste part !
1292!------------------------!
1293FIELD='U'
1294UNITS='m s{-1}'
1295DESC='Zonal wind'
1296DO k = 1,altlen
1297        XLVL=levels(k)
1298        SLAB=ufile(:,:,k,time_out(l))
1299                ! And now put everything in the destination file
1300                ! ... Header
1301        write(1) IFV
1302        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1303        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1304                ! ... Data
1305                write(1) SLAB
1306END DO
1307!print *,'The field '//DESC//' was written to '//output
1308
1309!------------------------!   
1310! >>> Write a variable   !
1311!    ... Copy&Paste part !
1312!------------------------!
1313FIELD='V'
1314UNITS='m s{-1}'
1315DESC='Meridional wind'
1316DO k = 1,altlen
1317        XLVL=levels(k)
1318        SLAB=vfile(:,:,k,time_out(l))
1319                ! And now put everything in the destination file
1320                ! ... Header
1321        write(1) IFV
1322        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1323        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1324                ! ... Data
1325                write(1) SLAB
1326END DO
1327!print *,'The field '//DESC//' was written to '//output
1328
1329!------------------------!   
1330! >>> Write a variable   !
1331!    ... Copy&Paste part !
1332!------------------------!
1333FIELD='RH'
1334UNITS=''
1335DESC='Customized 2D static field'
1336DESC='Eta-levels from the GCM'
1337DESC='Pressure values from the GCM'
1338DO k = 1,altlen
1339        XLVL=levels(k)
1340SELECT CASE(ident)
1341CASE('LMD')
1342        SLAB=aps(k)+bps(k)*psfile(:,:,time_out(l))
1343CASE('OXF')
1344        SLAB=bps(k)*psfile(:,:,time_out(l))
1345END SELECT
1346                ! And now put everything in the destination file
1347                ! ... Header
1348        write(1) IFV
1349        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1350        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1351                ! ... Data
1352                write(1) SLAB
1353END DO
1354!print *,'The field '//DESC//' was written to '//output
1355
1356!------------------------!   
1357! >>> Write a variable   !
1358!    ... Copy&Paste part !
1359!------------------------!
1360FIELD='HGT'
1361UNITS='m'
1362DESC='Height'
1363DO k = 1,altlen
1364        XLVL=levels(k)
1365!!*******
1366!! PROVISOIRE: normalement, il faudrait la hauteur geopotentielle
1367!!*******
1368!!however, not used by initialize_real
1369SELECT CASE(ident)
1370CASE('LMD')
1371        SLAB=10000.*alog(610./(aps(k)+bps(k)*psfile(:,:,time_out(l))))
1372CASE('OXF')
1373        SLAB=10000.*alog(610./(bps(k)*psfile(:,:,time_out(l))))
1374END SELECT
1375                ! And now put everything in the destination file
1376                ! ... Header
1377        write(1) IFV
1378        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1379        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
1380                ! ... Data
1381                write(1) SLAB
1382END DO
1383!print *,'The field '//DESC//' was written to '//output
1384
1385!------------------------!
1386! >>> Write a variable   !
1387!    ... Copy&Paste part !
1388!------------------------!
1389FIELD='HV'
1390UNITS='kg/kg'
1391DESC='Water vapor'
1392DO k = 1,altlen
1393        XLVL=levels(k)
1394        SLAB=waterfile(:,:,k,time_out(l))
1395                ! And now put everything in the destination file
1396                ! ... Header
1397        write(1) IFV
1398        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1399        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1400                ! ... Data
1401                write(1) SLAB
1402END DO
1403!print *,'The field '//DESC//' was written to '//output
1404
1405!------------------------!
1406! >>> Write a variable   !
1407!    ... Copy&Paste part !
1408!------------------------!
1409FIELD='HI'
1410UNITS='kg/kg'
1411DESC='Water ice'
1412DO k = 1,altlen
1413        XLVL=levels(k)
1414        SLAB=watericefile(:,:,k,time_out(l))
1415                ! And now put everything in the destination file
1416                ! ... Header
1417        write(1) IFV
1418        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1419        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1420                ! ... Data
1421                write(1) SLAB
1422END DO
1423!print *,'The field '//DESC//' was written to '//output
1424
1425!------------------------!
1426! >>> Write a variable   !
1427!    ... Copy&Paste part !
1428!------------------------!
1429FIELD='TSOIL'
1430UNITS='K'
1431DESC='Soil temperature'
1432DO k = 1,altlen
1433        XLVL=levels(k)
1434        SLAB=tsoilfile(:,:,k,time_out(l))
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='DSOIL'
1450UNITS='m'
1451DESC='Soil depths'
1452DO k = 1,altlen
1453        XLVL=levels(k)
1454        SLAB=dsoilfile(:,:,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='ISOIL'
1470UNITS='tiu'
1471DESC='Soil thermal inertia'
1472DO k = 1,altlen
1473        XLVL=levels(k)
1474        SLAB=isoilfile(:,:,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='DUSTQ'
1490UNITS='kg/kg'
1491DESC='Dust mixing ratio'
1492DO k = 1,altlen
1493        XLVL=levels(k)
1494        SLAB=dustfile(:,:,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='DUSTN'
1510UNITS='part/kg'
1511DESC='Dust number density'
1512DO k = 1,altlen
1513        XLVL=levels(k)
1514        SLAB=dustnfile(:,:,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!------------------------!
1527! >>> Write a variable   !
1528!     PHOTOCHEMISTRY     !
1529!------------------------!
1530#ifdef PHOTOCHEM
1531    DO i=1,nchemtrac
1532       FIELD=wtnom(i)
1533       UNITS='units'
1534       DESC='desc'
1535       DO k = 1,altlen
1536         XLVL=levels(k)
1537         SLAB=chemtrac(:,:,k,time_out(l),i)
1538         ! And now put everything in the destination file
1539         ! ... Header
1540         write(1) IFV
1541         write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
1542         write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
1543         ! ... Data
1544         write(1) SLAB
1545       END DO
1546    ENDDO
1547#endif
1548
1549print *,'****done file '//output, int(100.*float(l)/float(FILES)), ' % '
1550close(1)
1551
1552DEALLOCATE(SLAB)
1553
1554END DO !! end of the files loop
1555
1556
1557!!-------------------------
1558!! DeAllocate local arrays
1559!!-------------------------
1560deallocate(eta_gcm)
1561deallocate(tfile)
1562deallocate(tsoilfile)
1563deallocate(isoilfile)
1564deallocate(dsoilfile)
1565!deallocate(tfileorig)
1566deallocate(ufile)
1567!deallocate(ufileorig)
1568deallocate(vfile)
1569!deallocate(vfileorig)
1570deallocate(rfile)
1571deallocate(hfile) 
1572deallocate(waterfile)
1573deallocate(watericefile)
1574!deallocate(swaterfile)
1575deallocate(swatericefile)
1576deallocate(dustfile)
1577deallocate(dustnfile)
1578deallocate(psfile)
1579deallocate(tsfile)
1580deallocate(tnsfile)
1581deallocate(unsfile)
1582deallocate(vnsfile)
1583deallocate(emissfile)
1584deallocate(co2icefile)
1585deallocate(ghtsfile)    !! no scan axis
1586deallocate(vide)
1587deallocate(ones)
1588deallocate(lat, lon, alt, time)
1589deallocate(aps,bps,levels)
1590
1591#ifdef PHOTOCHEM
1592deallocate(chemtrac)
1593deallocate(wtnom)
1594#endif
1595
1596print *, '------------------------'
1597print *, 'End of pre-WPS process !'
1598print *, '------------------------'
1599print *, 'Here is what you should set in namelist.wps:'
1600print *, " start_date = '"//date_out2(1)//"'"
1601print *, " end_date = '"//date_out2(FILES)//"'"
1602print *, '------------------------'
1603print *, 'Any date between those is ok'
1604print *, 'for example, second available date is ... '//date_out2(2)
1605print *, '------------------------'
1606
1607end
Note: See TracBrowser for help on using the repository browser.