Context Navigation

readmeteo.F90 @ 87

Last change on this file since 87 was 73, checked in by aslmd, 14 years ago
LMD_MM_MARS: corrections cycle de l'eau propagees a la nouvelle physique. + corrections readmeteo.F90 [version synchronisee precedemment n etait pas la plus a jour] + corrections api.F90 pour avoir cp, R comme GCM
File size: 41.8 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format