Context Navigation

readmeteo_newphys.F90 @ 39

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

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

Download in other formats:

Original Format