source: trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/readmeteo.F90_backup @ 11

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

spiga@svn-planeto:ajoute le modele meso-echelle martien

File size: 26.2 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! A. Spiga - 16/04/2007                                                  !
15!            22/05/2007 : need to run zrecast as a preliminary           !
16!            07/06/2007 : external parameter file 'readmeteo.def         !
17!            30/07/2007 : manage additional variables (tsoil, emiss,...) !     
18!            19/10/2007 : no need to run zrecast at all (eta levels)     !   
19!               12/2007 : include co2ice and emissivity                  !
20!                                                                        !
21! spiga@lmd.jussieu.fr                                                   !     
22!------------------------------------------------------------------------!
23
24
25
26!***************************************************************************
27!***************************************************************************
28REAL, PARAMETER :: emiss_prescribed=0.95
29REAL, PARAMETER :: co2ice_prescribed=0.
30
31
32
33
34!***************************************************************************
35!***************************************************************************
36
37REAL :: ptop
38REAL, PARAMETER :: grav=3.72
39LOGICAL, PARAMETER :: blank=.false.
40CHARACTER (LEN=3), PARAMETER :: ident='LMD'
41
42
43!! Variables to be written in the output file
44!! *** NB: conformity with metgrid/src/read_met_module.F90
45CHARACTER (LEN=33) :: output
46INTEGER :: IFV
47CHARACTER (LEN=24) :: HDATE
48REAL :: XFCST
49CHARACTER (LEN=32) :: SOURCE
50CHARACTER (LEN=9) :: FIELD
51CHARACTER (LEN=25) :: UNITS
52CHARACTER (LEN=46) :: DESC
53REAL :: XLVL
54INTEGER :: NX,NY,IPROJ
55CHARACTER (LEN=8) :: STARTLOC
56REAL :: STARTLAT,STARTLON,DELTALAT,DELTALON
57REAL, POINTER, DIMENSION(:,:) :: SLAB
58
59!! Variables related to NETCDF format
60integer :: nid,nid2,nvarid,ierr,ierr2
61integer :: timedim,londim,latdim,altdim,altdim2
62integer :: rlatvdim,rlonudim
63integer :: timelen,lonlen,latlen,altlen,altlen2
64integer :: rlatvlen,rlonulen
65integer :: soillen,soildim
66integer, dimension(4) :: corner,edges
67
68!! Intermediate data arrays
69integer :: k,l
70real, dimension(:), allocatable :: lat,lon,time,alt,aps,bps,levels
71real, dimension(:,:), allocatable :: vide,ones,ghtsfile
72real, dimension(:,:,:), allocatable :: psfile,tsfile
73real, dimension(:,:,:), allocatable :: emissfile,co2icefile
74real, dimension(:,:,:), allocatable :: tnsfile,unsfile,vnsfile
75real, dimension(:,:,:,:), allocatable :: tfile,ufile,vfile,rfile,hfile
76real, dimension(:,:,:,:), allocatable :: eta_gcm
77real, dimension(:,:,:,:), allocatable :: tfileorig,ufileorig,vfileorig
78real, dimension(:,:,:,:), allocatable :: tsoilfile
79
80!! Reading the parameter file
81integer :: tmp,increment,FILES
82integer :: tmp2,tmp3,tmp4
83character*1 :: answer
84character*4 :: str
85character*2 :: str2, str3, str4
86integer, dimension(:), allocatable :: time_out
87character*13, dimension(:), allocatable :: date_out
88character*19, dimension(:), allocatable :: date_out2
89!***************************************************************************
90!***************************************************************************
91
92
93!!---------------------
94!! Open the datafile
95!!---------------------
96ierr=NF_OPEN ("input_diagfi.nc",NF_NOWRITE,nid)
97IF (ierr.NE.NF_NOERR) THEN
98   write(*,*)'**** Please create a symbolic link called input_diagfi.nc'
99   CALL ABORT
100ENDIF
101
102
103!!--------------------------
104!! Ask for data references
105!!--------------------------
106write(*,*) "Create n files. What is n ?"
107read(*,*) FILES
108allocate(time_out(FILES))
109allocate(date_out(FILES))
110allocate(date_out2(FILES))
111
112write(*,*) ""
113write(*,*) "INPUT: Diagfi time"
114write(*,*) "list of n subscripts, separated by <Return>s"
115increment=1
116do while (increment.ne.FILES+1)
117    read(*,*) tmp
118    time_out(increment)=tmp
119    increment=increment+1
120enddo
121
122write(*,*) ""
123write(*,*) "OUTPUT: WRF time"
124write(*,*) "fill 4 lines indicating"
125write(*,*) "-year (4 digits)"
126write(*,*) "-month (2 digits)"
127write(*,*) "-day (2 digits)"
128write(*,*) "-hour (2 digits)"
129increment=1
130do while (increment.ne.FILES+1)
131    read(*,*) str
132    read(*,*) str2
133    read(*,*) str3
134    read(*,*) str4
135    date_out(increment)=str//'-'//str2//'-'//str3//'_'//str4
136    date_out2(increment)=str//'-'//str2//'-'//str3//'_'//str4//':00:00'
137    print *,date_out(increment)
138    write(*,*) "ok? (y/n)"
139    read(*,*) answer
140    if (answer.eq.'n') then
141        write(*,*) "please write again"
142    else
143        write(*,*) "next one, please"   
144        increment=increment+1
145    endif
146enddo
147
148
149!!-------------------
150!! Get array sizes
151!!-------------------
152ierr=NF_INQ_DIMID(nid,"Time",timedim)
153ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
154  write(*,*) "timelen: ",timelen
155ierr=NF_INQ_DIMID(nid,"latitude",latdim)
156ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
157  write(*,*) "latlen: ",latlen
158ierr=NF_INQ_DIMID(nid,"longitude",londim)
159ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
160  write(*,*) "lonlen: ",lonlen
161ierr=NF_INQ_DIMID(nid,"altitude",altdim)
162ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
163  write(*,*) "altlen: ",altlen
164
165
166
167!!-------------------------
168!! Allocate local arrays
169!!-------------------------
170allocate(eta_gcm(lonlen,latlen,altlen,timelen))
171allocate(tfile(lonlen,latlen,altlen,timelen))
172allocate(tfileorig(lonlen,latlen,altlen,timelen))
173allocate(ufile(lonlen,latlen,altlen,timelen))
174allocate(ufileorig(lonlen,latlen,altlen,timelen))
175allocate(vfile(lonlen,latlen,altlen,timelen))
176allocate(vfileorig(lonlen,latlen,altlen,timelen))
177allocate(rfile(lonlen,latlen,altlen,timelen))
178allocate(hfile(lonlen,latlen,altlen,timelen)) 
179allocate(psfile(lonlen,latlen,timelen))
180allocate(tsfile(lonlen,latlen,timelen))
181allocate(tnsfile(lonlen,latlen,timelen))
182allocate(unsfile(lonlen,latlen,timelen))
183allocate(vnsfile(lonlen,latlen,timelen))
184allocate(emissfile(lonlen,latlen,timelen))
185allocate(co2icefile(lonlen,latlen,timelen))
186allocate(ghtsfile(lonlen,latlen))       !! no scan axis
187allocate(vide(lonlen,latlen))
188allocate(ones(lonlen,latlen))
189allocate(lat(latlen), lon(lonlen), alt(altlen), time(timelen))
190allocate(aps(altlen),bps(altlen),levels(altlen))
191
192tfile(:,:,:,:)=0
193tfileorig(:,:,:,:)=0
194ufileorig(:,:,:,:)=0
195vfileorig(:,:,:,:)=0
196ufile(:,:,:,:)=0
197vfile(:,:,:,:)=0
198rfile(:,:,:,:)=0
199hfile(:,:,:,:)=0
200psfile(:,:,:)=0
201tsfile(:,:,:)=0
202emissfile(:,:,:)=0
203co2icefile(:,:,:)=0
204tnsfile(:,:,:)=0
205unsfile(:,:,:)=0
206vnsfile(:,:,:)=0
207ghtsfile(:,:)=0
208vide(:,:)=0
209ones(:,:)=0
210
211
212!!-------------------
213!! Read dimensions
214!!-------------------
215
216print *,'>>> Read dimensions !'
217
218print *,'Time'
219   ierr = NF_INQ_VARID (nid, "Time",nvarid)
220   IF (ierr .NE. NF_NOERR) THEN
221      PRINT *, "Error: Readmeteo <Time> not found"
222      stop
223   ENDIF
224#ifdef NC_DOUBLE
225   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
226#else
227   ierr = NF_GET_VAR_REAL(nid, nvarid, time)
228#endif
229   print *,time(1),' ... to ... ',time(timelen)
230
231print *,'Latitude'
232   ierr = NF_INQ_VARID (nid, "latitude",nvarid)
233   IF (ierr .NE. NF_NOERR) THEN
234      PRINT *, "Error: Readmeteo <latitude> not found"
235      stop
236   ENDIF
237#ifdef NC_DOUBLE
238   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lat)
239#else
240   ierr = NF_GET_VAR_REAL(nid, nvarid, lat)
241#endif
242   print *,lat(1),' ... to ... ',lat(latlen),' ... step: ',lat(latlen)-lat(latlen-1)
243
244print *,'Longitude'
245   ierr = NF_INQ_VARID (nid, "longitude",nvarid)
246   IF (ierr .NE. NF_NOERR) THEN
247      PRINT *, "Error: Readmeteo <longitude> not found"
248      stop
249   ENDIF
250#ifdef NC_DOUBLE
251   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lon)
252#else
253   ierr = NF_GET_VAR_REAL(nid, nvarid, lon)
254#endif
255   print *,lon(1),' ... to ... ',lon(lonlen),' ... step: ',lon(lonlen)-lon(lonlen-1)
256
257print *, 'Hybrid coordinates'
258   ierr = NF_INQ_VARID (nid, "aps", nvarid)
259   IF (ierr .NE. NF_NOERR) THEN
260      PRINT *, "Error: Readmeteo <aps> not found"
261      stop
262   ENDIF
263#ifdef NC_DOUBLE
264   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, aps)
265#else
266   ierr = NF_GET_VAR_REAL(nid, nvarid, aps)
267#endif
268   ierr = NF_INQ_VARID (nid, "bps", nvarid)
269   IF (ierr .NE. NF_NOERR) THEN
270      PRINT *, "Error: Readmeteo <bps> not found"
271      stop
272   ENDIF
273#ifdef NC_DOUBLE
274   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bps)
275#else
276   ierr = NF_GET_VAR_REAL(nid, nvarid, bps)
277#endif
278   print *,aps(1),' ... to ... ',aps(altlen)
279   print *,bps(1),' ... to ... ',bps(altlen)
280
281
282
283
284!!-------------------------------------------
285!! Reading 2D and 3D meteorological fields
286!!-------------------------------------------
287
288IF (blank .EQV. .false.) THEN
289
290print *,'>>> Read 2D optional fields !'
291
292print *,'Emissivity'
293   ierr = NF_INQ_VARID (nid,"emiss",nvarid)
294IF (ierr .NE. NF_NOERR) THEN
295        PRINT *, '...warning: not found in diagfi !'
296        PRINT *, '...will be filled with a prescribed value', emiss_prescribed
297        emissfile(:,:,:)=emiss_prescribed       
298ELSE
299#ifdef NC_DOUBLE
300        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, emissfile)
301#else
302        ierr = NF_GET_VAR_REAL(nid, nvarid, emissfile)
303#endif
304ENDIF   
305
306print *,'CO2 ice'
307   ierr = NF_INQ_VARID (nid,"co2ice",nvarid)
308IF (ierr .NE. NF_NOERR) THEN
309        PRINT *, '...warning: not found in diagfi !'
310        PRINT *, '...will be filled with a prescribed value', co2ice_prescribed
311        co2icefile(:,:,:)=co2ice_prescribed     
312ELSE
313#ifdef NC_DOUBLE
314        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, co2icefile)
315#else
316        ierr = NF_GET_VAR_REAL(nid, nvarid, co2icefile)
317#endif
318ENDIF
319
320print *,'>>> Read 2D surface fields !'
321
322print *,'Surface Pressure'
323   ierr = NF_INQ_VARID (nid,"ps",nvarid)
324   IF (ierr .NE. NF_NOERR) THEN
325      PRINT *, "Error: Readmeteo <ps> not found"
326      stop
327   ENDIF
328#ifdef NC_DOUBLE
329   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, psfile)
330#else
331   ierr = NF_GET_VAR_REAL(nid, nvarid, psfile)
332#endif
333
334print *,'Ground Temperature'
335   ierr = NF_INQ_VARID (nid,"tsurf",nvarid)
336   IF (ierr .NE. NF_NOERR) THEN
337      PRINT *, "Error: Readmeteo <tsurf> not found"
338      stop
339   ENDIF
340#ifdef NC_DOUBLE
341   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsfile)
342#else
343   ierr = NF_GET_VAR_REAL(nid, nvarid, tsfile)
344#endif
345
346
347!!"atmospheric" surface temperature is taken
348!!from original diagfi.nc first level
349!!... level is ~3-5 meters
350print *,'Near-Surface Temperature'
351   ierr = NF_INQ_VARID (nid,"temp",nvarid)
352   IF (ierr .NE. NF_NOERR) THEN
353      PRINT *, "Error: Readmeteo <temp> not found"
354      stop
355   ENDIF
356#ifdef NC_DOUBLE
357   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfileorig)
358#else
359   ierr = NF_GET_VAR_REAL(nid, nvarid, tfileorig)
360#endif
361   tnsfile=tfileorig(:,:,1,:)
362
363!!"atmospheric" surface u is taken
364!!from original diagfi.nc first level
365!!... level is ~3-5 meters
366print *,'Near-Surface Zonal Wind'
367   ierr = NF_INQ_VARID (nid,"u",nvarid)
368   IF (ierr .NE. NF_NOERR) THEN
369     PRINT *, "Error: Readmeteo <u> not found"
370     stop
371   ENDIF
372#ifdef NC_DOUBLE
373   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufileorig)
374#else
375   ierr = NF_GET_VAR_REAL(nid, nvarid, ufileorig)
376#endif
377   unsfile=ufileorig(:,:,1,:)
378
379!!"atmospheric" surface v is taken
380!!from original diagfi.nc first level
381!!... level is ~3-5 meters
382print *,'Near-Surface Meridional Wind'
383   ierr = NF_INQ_VARID (nid,"v",nvarid)
384   IF (ierr .NE. NF_NOERR) THEN
385     PRINT *, "Error: Readmeteo <v> not found"
386     stop
387   ENDIF
388#ifdef NC_DOUBLE
389   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfileorig)
390#else
391   ierr = NF_GET_VAR_REAL(nid, nvarid, vfileorig)
392#endif
393   vnsfile=vfileorig(:,:,1,:)
394
395
396   print *,'Geopotential height at the ground'
397   ierr = NF_INQ_VARID (nid,"phisinit",nvarid)
398   IF (ierr .NE. NF_NOERR) THEN
399      PRINT *, "Error: Readmeteo <phisinit> not found"
400      stop
401   ENDIF
402#ifdef NC_DOUBLE
403   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ghtsfile)
404#else
405   ierr = NF_GET_VAR_REAL(nid, nvarid, ghtsfile)
406#endif
407!**** from geopotential to geopotential height
408ghtsfile=ghtsfile/grav
409!****
410
411print *,'>>> Read 3D meteorological fields ! - This may take some time ...'
412
413print *,'Temperature'
414   ierr = NF_INQ_VARID (nid,"temp",nvarid)
415   IF (ierr .NE. NF_NOERR) THEN
416      PRINT *, "Error: Readmeteo <t> not found"
417      stop
418   ENDIF
419#ifdef NC_DOUBLE
420   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfile)
421#else
422   ierr = NF_GET_VAR_REAL(nid, nvarid, tfile)
423#endif
424
425print *,'Zonal wind'   
426   ierr = NF_INQ_VARID (nid,"u",nvarid)
427   IF (ierr .NE. NF_NOERR) THEN
428      PRINT *, "Error: Readmeteo <u> not found"
429      stop
430   ENDIF
431#ifdef NC_DOUBLE
432   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufile)
433#else
434   ierr = NF_GET_VAR_REAL(nid, nvarid, ufile)
435#endif
436
437print *,'Meridional wind'
438   ierr = NF_INQ_VARID (nid,"v",nvarid)
439   IF (ierr .NE. NF_NOERR) THEN
440      PRINT *, "Error: Readmeteo <v> not found"
441      stop
442   ENDIF
443#ifdef NC_DOUBLE
444   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfile)
445#else
446   ierr = NF_GET_VAR_REAL(nid, nvarid, vfile)
447#endif
448
449ENDIF
450
451
452!!-----------------------------
453!! Loop on the written files
454!! >>> what follows is pretty repetitive ...
455!!        gotta do something (one more loop?)
456!!-----------------------------
457
458!!! Equivalent eta levels for WRF interpolation in initialize_real
459!print *,'Computing equivalent eta levels'
460!       
461!       !ptop will be passed through RH(surface)
462!       ptop=MINVAL(aps(altlen)+bps(altlen)*psfile(:,:,:))
463!       print *, 'ptop', ptop
464!
465!print *, 'sample: equivalent eta levels at i=1,j=1'
466!DO k = 1,altlen
467!        levels(k)=-k
468!        eta_gcm(:,:,k,:)=(aps(k)+bps(k)*psfile(:,:,:)-ptop)/(psfile(:,:,:)-ptop)
469!        print *, eta_gcm(1,1,k,1)
470!END DO
471
472!!---better to pass pressure values
473!!---the eta treatment is done in module_initialize_real
474DO k = 1,altlen
475        levels(k)=-k
476        !!dummy decreasing levels
477END DO
478
479
480
481
482!! Dummy surface level is XLVL=200100.
483
484
485DO l=1,FILES
486
487
488!!---------------------------------------------
489!! Write file in intermediate format for WPS
490!! 1. Surface data
491!!---------------------------------------------
492
493!
494! These variables remain the same in the "loop"
495!
496HDATE=date_out(l)
497IFV=4
498XFCST=0.
499SOURCE='LMD GCM'
500NX=lonlen
501NY=latlen
502ALLOCATE(SLAB(NX,NY))
503IPROJ=0
504STARTLOC='SWCORNER'
505STARTLAT=lat(1)
506STARTLON=lon(1)
507DELTALAT=lat(latlen)-lat(latlen-1)
508DELTALON=lon(lonlen)-lon(lonlen-1)
509!
510! Open the data destination file
511!
512output=ident//"/"//ident//":"//date_out2(l)       
513open(UNIT=1,                    &
514        FILE=output,            &
515        STATUS='REPLACE',       &
516        FORM='UNFORMATTED',     &
517        ACCESS='SEQUENTIAL')
518
519!------------------------!   
520! >>> Write a variable   !
521!    ... Copy&Paste part !
522!------------------------!
523FIELD='T'
524UNITS='K'
525DESC='Atmospheric temperature'
526XLVL=200100.
527!SLAB=tsfile(:,:,time_out(l))
528!SLAB=tfileorig(:,:,1,time_out(l))
529SLAB=tnsfile(:,:,time_out(l))
530        ! And now put everything in the destination file
531        ! ... Header
532        write(1) IFV
533        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
534        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
535        ! ... Data
536        write(1) SLAB
537print *,'The field '//DESC//' was written to '//output
538
539!------------------------!   
540! >>> Write a variable   !
541!    ... Copy&Paste part !
542!------------------------!
543FIELD='U'
544UNITS='m s{-1}'
545DESC='Zonal wind'
546XLVL=200100.
547!SLAB=ufile(:,:,1,time_out(l))
548!SLAB=ufileorig(:,:,1,time_out(l))
549SLAB=unsfile(:,:,time_out(l))
550        ! And now put everything in the destination file
551        ! ... Header
552        write(1) IFV
553        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
554        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
555        ! ... Data
556        write(1) SLAB
557print *,'The field '//DESC//' was written to '//output
558
559!------------------------!   
560! >>> Write a variable   !
561!    ... Copy&Paste part !
562!------------------------!
563FIELD='V'
564UNITS='m s{-1}'
565DESC='Meridional wind'
566XLVL=200100.
567!SLAB=vfile(:,:,1,time_out(l))
568!SLAB=vfileorig(:,:,1,time_out(l))
569SLAB=vnsfile(:,:,time_out(l))
570        ! And now put everything in the destination file
571        ! ... Header
572        write(1) IFV
573        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
574        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
575        ! ... Data
576        write(1) SLAB
577print *,'The field '//DESC//' was written to '//output
578
579!------------------------!   
580! >>> Write a variable   !
581!    ... Copy&Paste part !
582!------------------------!
583FIELD='RH'
584UNITS=''
585DESC='Customized 2D static field'
586XLVL=200100.
587!SLAB=co2icefile(:,:,time_out(l))
588SLAB=vide(:,:)
589!SLAB=vide(:,:)+ptop
590        ! And now put everything in the destination file
591        ! ... Header
592        write(1) IFV
593        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
594        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
595        ! ... Data
596        write(1) SLAB
597print *,'The field '//DESC//' was written to '//output
598
599!------------------------!   
600! >>> Write a variable   !
601!    ... Copy&Paste part !
602!------------------------!
603FIELD='SOILHGT'
604UNITS='m'
605DESC='Terrain field of source analysis'  !Ground geopotential height
606XLVL=200100.
607SLAB=ghtsfile(:,:)
608        ! And now put everything in the destination file
609        ! ... Header
610        write(1) IFV
611        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
612        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
613        ! ... Data
614        write(1) SLAB
615print *,'The field '//DESC//' was written to '//output
616
617!------------------------!   
618! >>> Write a variable   !
619!    ... Copy&Paste part !
620!------------------------!
621FIELD='PSFC'
622UNITS='Pa'
623DESC='Surface Pressure'
624XLVL=200100.
625SLAB=psfile(:,:,time_out(l))
626        ! And now put everything in the destination file
627        ! ... Header
628        write(1) IFV
629        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
630        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
631        ! ... Data
632        write(1) SLAB
633print *,'The field '//DESC//' was written to '//output
634
635!------------------------!   
636! >>> Write a variable   !
637!    ... Copy&Paste part !
638!------------------------!
639FIELD='ST000010'
640UNITS=''
641DESC='Emissivity'
642XLVL=200100.
643!SLAB=vide(:,:)
644SLAB=emissfile(:,:,time_out(l))
645        ! And now put everything in the destination file
646        ! ... Header
647        write(1) IFV
648        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
649        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
650        ! ... Data
651        write(1) SLAB
652print *,'The field '//DESC//' was written to '//output
653
654!------------------------!   
655! >>> Write a variable   !
656!    ... Copy&Paste part !
657!------------------------!
658FIELD='ST010040'
659UNITS=''
660DESC='CO2 ice'
661XLVL=200100.
662!SLAB=vide(:,:)
663SLAB=co2icefile(:,:,time_out(l))
664        ! And now put everything in the destination file
665        ! ... Header
666        write(1) IFV
667        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
668        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
669        ! ... Data
670        write(1) SLAB
671print *,'The field '//DESC//' was written to '//output
672
673!------------------------!   
674! >>> Write a variable   !
675!    ... Copy&Paste part !
676!------------------------!
677FIELD='ST040100'
678UNITS='K'
679DESC='T of 40-100 cm ground layer'
680XLVL=200100.
681SLAB=vide(:,:)
682        ! And now put everything in the destination file
683        ! ... Header
684        write(1) IFV
685        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
686        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
687        ! ... Data
688        write(1) SLAB
689print *,'The field '//DESC//' was written to '//output
690
691!------------------------!   
692! >>> Write a variable   !
693!    ... Copy&Paste part !
694!------------------------!
695FIELD='ST100200'
696UNITS='K'
697DESC='T of 100-200 cm ground layer'
698XLVL=200100.
699SLAB=vide(:,:)
700        ! And now put everything in the destination file
701        ! ... Header
702        write(1) IFV
703        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
704        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
705        ! ... Data
706        write(1) SLAB
707print *,'The field '//DESC//' was written to '//output
708
709!------------------------!   
710! >>> Write a variable   !
711!    ... Copy&Paste part !
712!------------------------!
713FIELD='LANDSEA'
714UNITS='0/1 Flag'
715DESC='Land/Sea flag'
716XLVL=200100.
717SLAB=ones(:,:)
718        ! And now put everything in the destination file
719        ! ... Header
720        write(1) IFV
721        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
722        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
723        ! ... Data
724        write(1) SLAB
725print *,'The field '//DESC//' was written to '//output
726
727!------------------------!   
728! >>> Write a variable   !
729!    ... Copy&Paste part !
730!------------------------!
731FIELD='SKINTEMP'
732UNITS='K'
733DESC='Ground temperature'
734XLVL=200100.
735!SLAB=vide(:,:)
736SLAB=tsfile(:,:,time_out(l))
737        ! And now put everything in the destination file
738        ! ... Header
739        write(1) IFV
740        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
741        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
742        ! ... Data
743        write(1) SLAB
744print *,'The field '//DESC//' was written to '//output
745
746!------------------------!   
747! >>> Write a variable   !
748!    ... Copy&Paste part !
749!------------------------!
750FIELD='SM000010'
751UNITS='fraction'
752DESC='Relative humidity'
753XLVL=200100.
754SLAB=vide(:,:)
755        ! And now put everything in the destination file
756        ! ... Header
757        write(1) IFV
758        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
759        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
760        ! ... Data
761        write(1) SLAB
762print *,'The field '//DESC//' was written to '//output
763
764!------------------------!   
765! >>> Write a variable   !
766!    ... Copy&Paste part !
767!------------------------!
768FIELD='SM010040'
769UNITS='fraction'
770DESC='Relative humidity'
771XLVL=200100.
772SLAB=vide(:,:)
773        ! And now put everything in the destination file
774        ! ... Header
775        write(1) IFV
776        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
777        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
778        ! ... Data
779        write(1) SLAB
780print *,'The field '//DESC//' was written to '//output
781
782!------------------------!   
783! >>> Write a variable   !
784!    ... Copy&Paste part !
785!------------------------!
786FIELD='SM040100'
787UNITS='fraction'
788DESC='Relative humidity'
789XLVL=200100.
790SLAB=vide(:,:)
791        ! And now put everything in the destination file
792        ! ... Header
793        write(1) IFV
794        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
795        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
796        ! ... Data
797        write(1) SLAB
798print *,'The field '//DESC//' was written to '//output
799
800!------------------------!   
801! >>> Write a variable   !
802!    ... Copy&Paste part !
803!------------------------!
804FIELD='SM100200'
805UNITS='fraction'
806DESC='Relative humidity'
807XLVL=200100.
808SLAB=vide(:,:)
809        ! And now put everything in the destination file
810        ! ... Header
811        write(1) IFV
812        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
813        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
814        ! ... Data
815        write(1) SLAB
816print *,'The field '//DESC//' was written to '//output
817
818
819
820!!----------------------------------------------------
821!! Write file in intermediate format for WPS
822!! 2. 3D meteorological data
823!! >>> same stuff, but handle vertical profile
824!! >>> !! XLVL must be decreasing (pressure levels)
825!!----------------------------------------------------
826
827!------------------------!   
828! >>> Write a variable   !
829!    ... Copy&Paste part !
830!------------------------!
831FIELD='T'
832UNITS='K'
833DESC='Atmospheric temperature'
834DO k = 1,altlen
835        XLVL=levels(k)
836        SLAB=tfile(:,:,k,time_out(l))
837                ! And now put everything in the destination file
838                ! ... Header
839        write(1) IFV
840        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
841        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
842                ! ... Data
843                write(1) SLAB
844END DO
845print *,'The field '//DESC//' was written to '//output
846
847!------------------------!   
848! >>> Write a variable   !
849!    ... Copy&Paste part !
850!------------------------!
851FIELD='U'
852UNITS='m s{-1}'
853DESC='Zonal wind'
854DO k = 1,altlen
855        XLVL=levels(k)
856        SLAB=ufile(:,:,k,time_out(l))
857                ! And now put everything in the destination file
858                ! ... Header
859        write(1) IFV
860        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
861        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
862                ! ... Data
863                write(1) SLAB
864END DO
865print *,'The field '//DESC//' was written to '//output
866
867!------------------------!   
868! >>> Write a variable   !
869!    ... Copy&Paste part !
870!------------------------!
871FIELD='V'
872UNITS='m s{-1}'
873DESC='Meridional wind'
874DO k = 1,altlen
875        XLVL=levels(k)
876        SLAB=vfile(:,:,k,time_out(l))
877                ! And now put everything in the destination file
878                ! ... Header
879        write(1) IFV
880        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
881        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
882                ! ... Data
883                write(1) SLAB
884END DO
885print *,'The field '//DESC//' was written to '//output
886
887!------------------------!   
888! >>> Write a variable   !
889!    ... Copy&Paste part !
890!------------------------!
891FIELD='RH'
892UNITS=''
893DESC='Customized 2D static field'
894DESC='Eta-levels from the GCM'
895DESC='Pressure values from the GCM'
896DO k = 1,altlen
897        XLVL=levels(k)
898        !SLAB=rfile(:,:,k,time_out(l))
899        !SLAB=eta_gcm(:,:,k,time_out(l))
900        SLAB=aps(k)+bps(k)*psfile(:,:,time_out(l))
901                ! And now put everything in the destination file
902                ! ... Header
903        write(1) IFV
904        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
905        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
906                ! ... Data
907                write(1) SLAB
908END DO
909print *,'The field '//DESC//' was written to '//output
910
911!------------------------!   
912! >>> Write a variable   !
913!    ... Copy&Paste part !
914!------------------------!
915FIELD='HGT'
916UNITS='m'
917DESC='Height'
918DO k = 1,altlen
919        XLVL=levels(k)
920!!*******
921!! PROVISOIRE: normalement, il faudrait la hauteur geopotentielle
922!!*******
923!!SLAB=hfile(:,:,k,time_out(l))
924!!SLAB=vide(:,:)
925        SLAB=10000.*alog(610./(aps(k)+bps(k)*psfile(:,:,time_out(l))))
926!!however, not used by initialize_real
927                ! And now put everything in the destination file
928                ! ... Header
929        write(1) IFV
930        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
931        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
932                ! ... Data
933                write(1) SLAB
934END DO
935print *,'The field '//DESC//' was written to '//output
936
937print *,'****Close file '//output
938close(1)
939
940DEALLOCATE(SLAB)
941
942END DO !! end of the files loop
943
944
945!!-------------------------
946!! DeAllocate local arrays
947!!-------------------------
948deallocate(eta_gcm)
949deallocate(tfile)
950deallocate(tfileorig)
951deallocate(ufile)
952deallocate(ufileorig)
953deallocate(vfile)
954deallocate(vfileorig)
955deallocate(rfile)
956deallocate(hfile) 
957deallocate(psfile)
958deallocate(tsfile)
959deallocate(tnsfile)
960deallocate(unsfile)
961deallocate(vnsfile)
962deallocate(emissfile)
963deallocate(co2icefile)
964deallocate(ghtsfile)    !! no scan axis
965deallocate(vide)
966deallocate(ones)
967deallocate(lat, lon, alt, time)
968deallocate(aps,bps,levels)
969
970
971print *, 'End of pre-WPS process !'
972
973
974end
Note: See TracBrowser for help on using the repository browser.