Context Navigation

readmeteo_newphys.F90 @ 55

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

LMD_MM_MARS: element manquant pour runs traceurs avec nouvelle physique
--> avec mars=0 et non actif, le modele fonctionne toujours bien [avec starts JBM et conditions de ses runs]
--> mettre les traceurs avec mars=11 [actifs ou non] semble faire crasher le modele: pourquoi ?

M 54 mesoscale/LMD_MM_MARS/SRC/WRFV2/call_meso_physiq1.inc
M 54 mesoscale/LMD_MM_MARS/SRC/WRFV2/call_meso_physiq2.inc
M 54 mesoscale/LMD_MM_MARS/SRC/WRFV2/call_meso_physiq3.inc
M 54 mesoscale/LMD_MM_MARS/SRC/WRFV2/call_meso_physiq4.inc
M 54 mars/libf/phymars/meso_physiq.F
le tableau tnom est defini dans module_lmd_driver en fonction de config_flag%mars [MARS_MODE]

et passe dans meso_physiq [pour ensuite servir dans initracer]

--> le COMMON advtrac.h est alors necessaire

M 54 mesoscale/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM
config_flag%mars defini une serie de traceurs ordonnee dans le tableau scalar

M 54 mesoscale/LMD_MM_MARS/SRC/WRFV2/share/module_model_constants.F
on passe au cp et R du GCM martien pour une complete correspondance

M 54 mesoscale/LMD_MM_MARS/SRC/WRFV2/call_meso_physiq.inc
ce fichier est en fait inutile et pourrait etre supprime

M 54 mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_real.F
M 54 mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/solve_em.F
M 54 mesoscale/LMD_MM_MARS/SRC/WRFV2/phys/module_lmd_driver.F
M 54 mesoscale/LMD_MM_MARS/SRC/WRFV2/main/real_em.F
modifications pour prise en compte des traceurs avec la nouvelle physique
--> avantage des traceurs generiques dans la physique
--> pour l'instant dans SCALAR l'indice 1 est dummy, 2/3 water vapor/ice et dernier CO2

M 54 mesoscale/LMD_MM_MARS/SRC/PREP_MARS/readmeteo_newphys.F90
M 54 mesoscale/LMD_MM_MARS/SRC/WPS/wps_mars/metgrid/METGRID.TBL.ARW_MarsBase_newphys
modifications pour possibilites de passer les tableaux necessaires pour

initialiser et guider les traceurs dans la nouvelle physique

M mars/libf/phymars/dimradmars.h
M mars/libf/phymars/callradite.F
version traceurs non actifs pour tests basiques

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

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

Download in other formats:

Original Format