Context Navigation

localtime.F90 @ 2440

Last change on this file since 2440 was 2440, checked in by abierjon, 4 years ago

Mars GCM:
Harmonization of the utility programs lslin, concatnc and localtime when dealing
with Solar longitude as the Time variable (for localtime, the check on long_name
attribute is still compliant with files from previous versions)

File size: 42.6 KB

Line
1	program localtime
2
3	! ----------------------------------------------------------------------------
4	! Program to redistribute and interpolate the variable a the same
5	! local times everywhere
6	! input : diagfi.nc / concat.nc / stats.nc kind of files
7	! author: F. Forget
8	! ----------------------------------------------------------------------------
9
10	implicit none
11
12	include "netcdf.inc" ! NetCDF definitions
13
14	character (len=50) file
15	! file(): input file(s) names(s)
16	character (len=30), dimension(16) :: notprocessed
17	! notprocessed(): names of the (16) variables that won't be processed
18	character (len=50), dimension(:), allocatable :: var
19	! var(): name(s) of variable(s) that will be processed
20	character (len=100) :: tmpvar,long_name,units
21	! tmpvar(): used to temporarily store a variable name
22	! long_name(): [netcdf] long_name attribute
23	! units(): [netcdf] units attribute
24	character (len=100) :: filename,vartmp
25	! filename(): output file name
26	! vartmp(): temporary variable name (read from netcdf input file)
27	!character (len=1) :: ccopy
28	! ccpy: 'y' or 'n' answer
29
30	character (len=50) :: altlong_name,altunits,altpositive
31	! altlong_name(): [netcdf] altitude "long_name" attribute
32	! altunits(): [netcdf] altitude "units" attribute
33	! altpositive(): [netcdf] altitude "positive" attribute
34
35	integer :: nid,ierr,miss,validr
36	! nid: [netcdf] file ID #
37	! ierr: [netcdf] subroutine returned error code
38	! miss: [netcdf] subroutine returned error code
39	integer :: i,j,k,inter
40	! for various loops
41	integer :: varid
42	! varid: [netcdf] variable ID #
43	real, dimension(:), allocatable:: lat,lon,alt,ctl,time
44	! lat(): array, stores latitude coordinates
45	! lon(): array, stores longitude coordinates
46	! alt(): array, stores altitude coordinates
47	! ctl(): array, stores controle array
48	! time(): array, stores time coordinates
49	integer :: nbvar,nbvarfile,ndim
50	! nbvar: # of variables to concatenate
51	! nbfile: # number of input file(s)
52	! nbvarfile: total # of variables in an input file
53	! ndim: [netcdf] # (3 or 4) of dimensions (for variables)
54	integer :: latdim,londim,altdim,ctldim,timedim
55	! latdim: [netcdf] "latitude" dim ID
56	! londim: [netcdf] "longitude" dim ID
57	! altdim: [netcdf] "altdim" dim ID
58	! ctldim: [netcdf] "controle" dim ID
59	! timedim: [netcdf] "timedim" dim ID
60	integer :: gcmlayerdim ! NetCDF dimension ID for # of layers in GCM
61	integer :: latvar,lonvar,altvar,ctlvar,timevar
62	! latvar: [netcdf] ID of "latitude" variable
63	! lonvar: [netcdf] ID of "longitude" variable
64	! altvar: [netcdf] ID of "altitude" variable
65	! ctlvar: [netcdf] ID of "controle" variable
66	! timevar: [netcdf] ID of "Time" variable
67	integer :: latlen,lonlen,altlen,ctllen,timelen,timelen_lt,timelen_tot
68	integer :: ilat,ilon,ialt,it
69	! latlen: # of elements of lat() array
70	! lonlen: # of elements of lon() array
71	! altlen: # of elements of alt() array
72	! ctllen: # of elements of ctl() array
73	! timelen: # of elemnets of time() array
74	! timelen_tot:# =timelen or timelen+1 (if 1 more time to interpolate needed)
75	! timelen_lt: # of elemnets of time() array in output
76	integer :: nhour,nsol
77	integer :: GCM_layers ! number of GCM atmospheric layers (may not be
78	! same as altlen if file is an output of zrecast)
79	integer :: nout,latdimout,londimout,altdimout,timedimout,timevarout
80	! nout: [netcdf] output file ID
81	! latdimout: [netcdf] output latitude (dimension) ID
82	! londimout: [netcdf] output longitude (dimension) ID
83	! altdimout: [netcdf] output altitude (dimension) ID
84	! timedimout: [netcdf] output time (dimension) ID
85	! timevarout: [netcdf] ID of output "Time" variable
86	integer :: layerdimout ! NetCDF dimension ID for # of layers in GCM
87	integer :: varidout
88	! varidout: [netcdf] variable ID # (of a variable to write to the output file)
89	integer :: Nnotprocessed,var_ok
90	! Nnotprocessed: # of (leading)variables that won't be concatenated
91	! var_ok: flag (0 or 1)
92	integer, dimension(4) :: corner,edges,dim
93	! corner: [netcdf]
94	! edges: [netcdf]
95	! dim: [netcdf]
96	real, dimension(:,:,:,:), allocatable :: var3d
97	! var3D(,,,): 4D array to store a field
98
99	real, dimension(:,:,:,:), allocatable :: var3d_lt
100	! var3D_lt(,,,): 4D array to store a field in local time coordinate
101	real, dimension(:), allocatable :: lt_gcm,lt_hour
102	real, dimension(:), allocatable :: lt_out,lt_outc
103	real, dimension(:), allocatable :: var_gcm
104
105	real :: missing
106	!PARAMETER(missing=1E+20)
107	! missing: [netcdf] to handle "missing" values when reading/writing files
108	real, dimension(2) :: valid_range
109	! valid_range(2): [netcdf] interval in which a value is considered valid
110	logical :: stats ! stats=T when reading a "stats" kind of ile
111
112
113	!==============================================================================
114	! 1.1. Get input file name(s)
115	!==============================================================================
116	write(,)
117	write(,) "which file do you want to use? (diagfi... stats... concat...)"
118
119	read(*,'(a50)') file
120
121	if (len_trim(file).eq.0) then
122	write(,) "no file... game over"
123	stop ""
124	endif
125
126	!==============================================================================
127	! 1.3. Open the input file
128	!==============================================================================
129
130	ierr = NF_OPEN(file,NF_NOWRITE,nid)
131	if (ierr.NE.NF_NOERR) then
132	write(,) 'ERROR: Pb opening file '//trim(file)
133	write(,) NF_STRERROR(ierr)
134	stop ""
135	endif
136
137	ierr=NF_INQ_NVARS(nid,nbvarfile)
138	! nbvarfile now set to be the (total) number of variables in file
139	if (ierr.NE.NF_NOERR) then
140	write(,) 'ERROR: Pb with NF_INQ_NVARS'
141	write(,) NF_STRERROR(ierr)
142	stop ""
143	endif
144
145	!==============================================================================
146	! 1.4. List of variables that should not be processed
147	!==============================================================================
148
149	notprocessed(1)='Time'
150	notprocessed(2)='controle'
151	notprocessed(3)='rlonu'
152	notprocessed(4)='latitude'
153	notprocessed(5)='longitude'
154	notprocessed(6)='altitude'
155	notprocessed(7)='rlatv'
156	notprocessed(8)='aps'
157	notprocessed(9)='bps'
158	notprocessed(10)='ap'
159	notprocessed(11)='bp'
160	notprocessed(12)='soildepth'
161	notprocessed(13)='cu'
162	notprocessed(14)='cv'
163	notprocessed(15)='aire'
164	notprocessed(16)='phisinit'
165
166
167	!==============================================================================
168	! 1.5. Get (and check) list of variables to process
169	!==============================================================================
170	write(,)
171	Nnotprocessed=0
172	do i=1,nbvarfile
173	ierr=NF_INQ_VARNAME(nid,i,vartmp)
174	! vartmp now contains the "name" of variable of ID # i
175	var_ok=0
176	do inter=1,16
177	if (vartmp.eq.notprocessed(inter)) then
178	var_ok=1
179	Nnotprocessed=Nnotprocessed+1
180	endif
181	enddo
182	if (var_ok.eq.0) write(,) trim(vartmp)
183	enddo
184
185	! Nnotprocessed: # of variables that won't be processed
186	! nbvarfile: total # of variables in file
187	allocate(var(nbvarfile-Nnotprocessed),stat=ierr)
188	if (ierr.ne.0) then
189	write(,) "Error: failed allocation of var(nbvarfile-Nnotprocessed)"
190	write(,) " nbvarfile=",nbvarfile
191	write(,) " Nnotprocessed=",Nnotprocessed
192	stop
193	endif
194
195
196	write(,)
197	write(,) "which variables do you want to redistribute ?"
198	write(,) "all / list of <variables> (separated by <Enter>s)"
199	write(,) "(an empty line , i.e: just <Enter>, implies end of list)"
200	nbvar=0
201	read(*,'(a50)') tmpvar
202	do while ((tmpvar/=' ').AND.(trim(tmpvar)/='all'))
203	nbvar=nbvar+1
204	var(nbvar)=tmpvar
205	read(*,'(a50)') tmpvar
206	enddo
207
208	if (tmpvar=="all") then
209	nbvar=nbvarfile-Nnotprocessed
210	do j=Nnotprocessed+1,nbvarfile
211	ierr=nf_inq_varname(nid,j,var(j-Nnotprocessed))
212	enddo
213	! Variables names from the file are stored in var()
214	nbvar=nbvarfile-Nnotprocessed
215	do i=1,nbvar
216	ierr=nf_inq_varname(nid,i+Nnotprocessed,var(i))
217	write(*,'(a9,1x,i2,1x,a1,1x,a50)') "variable ",i,":",var(i)
218	enddo
219	else if(nbvar==0) then
220	write(,) "no variable... game over"
221	stop ""
222	endif ! of if (tmpvar=="all")
223
224	!==============================================================================
225	! 1.6. Get output file name
226	!==============================================================================
227	filename=file(1:len_trim(file)-3)//"_LT.nc"
228	write(,) "Output file :"//trim(filename)
229
230	!==============================================================================
231	! 2.1. Open input file
232	!==============================================================================
233
234	write(,)
235	write(,) "opening "//trim(file)//"..."
236	ierr = NF_OPEN(file,NF_NOWRITE,nid)
237	if (ierr.NE.NF_NOERR) then
238	write(,) 'ERROR: Pb opening file '//trim(file)
239	write(,) NF_STRERROR(ierr)
240	stop ""
241	endif
242
243	!==============================================================================
244	! 2.2. Read (and check) dimensions of variables from input file
245	!==============================================================================
246
247	ierr=NF_INQ_DIMID(nid,"latitude",latdim)
248	if (ierr.NE.NF_NOERR) then
249	write(,) 'ERROR: Dimension <latitude> is missing in file '//trim(file)
250	stop ""
251	endif
252	ierr=NF_INQ_VARID(nid,"latitude",latvar)
253	if (ierr.NE.NF_NOERR) then
254	write(,) 'ERROR: Field <latitude> is missing in file '//trim(file)
255	stop ""
256	endif
257	ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
258	! write(,) "latlen: ",latlen
259
260	ierr=NF_INQ_DIMID(nid,"longitude",londim)
261	if (ierr.NE.NF_NOERR) then
262	write(,) 'ERROR: Dimension <longitude> is missing in file '//trim(file)
263	stop ""
264	endif
265	ierr=NF_INQ_VARID(nid,"longitude",lonvar)
266	if (ierr.NE.NF_NOERR) then
267	write(,) 'ERROR: Field <longitude> is missing in file '//trim(file)
268	stop ""
269	endif
270	ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
271	! write(,) "lonlen: ",lonlen
272
273	ierr=NF_INQ_DIMID(nid,"altitude",altdim)
274	if (ierr.NE.NF_NOERR) then
275	write(,) 'ERROR: Dimension <altitude> is missing in file '//trim(file)
276	stop ""
277	endif
278	ierr=NF_INQ_VARID(nid,"altitude",altvar)
279	if (ierr.NE.NF_NOERR) then
280	write(,) 'ERROR: Field <altitude> is missing in file '//trim(file)
281	stop ""
282	endif
283	ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
284	! write(,) "altlen: ",altlen
285
286	! load size of aps() or sigma() (in case it is not altlen)
287	! default is that GCM_layers=altlen
288	! but for outputs of zrecast, it may be a different value
289	ierr=NF_INQ_DIMID(nid,"GCM_layers",gcmlayerdim)
290	if (ierr.ne.NF_NOERR) then
291	! didn't find a GCM_layers dimension; therefore we have:
292	GCM_layers=altlen
293	else
294	! load value of GCM_layers
295	ierr=NF_INQ_DIMLEN(nid,gcmlayerdim,GCM_layers)
296	endif
297	! write(,) "GCM_layers=",GCM_layers
298
299	ierr=NF_INQ_DIMID(nid,"index",ctldim)
300	if (ierr.NE.NF_NOERR) then
301	write(,) 'Dimension <index> is missing in file '//trim(file)
302	ctllen=0
303	!stop ""
304	endif
305	ierr=NF_INQ_VARID(nid,"controle",ctlvar)
306	if (ierr.NE.NF_NOERR) then
307	write(,) 'Field <controle> is missing in file '//trim(file)
308	ctllen=0
309	!stop ""
310	else
311	ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
312	endif
313	! write(,) "controle: ",controle
314
315	!==============================================================================
316	! 2.3. Read (and check compatibility of) dimensions of
317	! variables from input file
318	!==============================================================================
319
320	allocate(lat(latlen),stat=ierr)
321	if (ierr.ne.0) then
322	write(,) "Error: failed to allocate lat(latlen)"
323	stop
324	endif
325	allocate(lon(lonlen),stat=ierr)
326	if (ierr.ne.0) then
327	write(,) "Error: failed to allocate lon(lonlen)"
328	stop
329	endif
330	allocate(alt(altlen),stat=ierr)
331	if (ierr.ne.0) then
332	write(,) "Error: failed to allocate alt(altlen)"
333	stop
334	endif
335	allocate(ctl(ctllen),stat=ierr)
336	if (ierr.ne.0) then
337	write(,) "Error: failed to allocate ctl(ctllen)"
338	stop
339	endif
340
341	ierr = NF_GET_VAR_REAL(nid,latvar,lat)
342	if (ierr.ne.0) then
343	write(,) "Error: failed to load latitude"
344	write(,) NF_STRERROR(ierr)
345	stop
346	endif
347
348	ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
349	if (ierr.ne.0) then
350	write(,) "Error: failed to load longitude"
351	write(,) NF_STRERROR(ierr)
352	stop
353	endif
354
355	ierr = NF_GET_VAR_REAL(nid,altvar,alt)
356	if (ierr.ne.0) then
357	write(,) "Error: failed to load altitude"
358	write(,) NF_STRERROR(ierr)
359	stop
360	endif
361	! Get altitude attributes to handle files with any altitude type
362	ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
363	if (ierr.ne.nf_noerr) then
364	! if no attribute "long_name", try "title"
365	ierr=nf_get_att_text(nid,timevar,"title",long_name)
366	endif
367	ierr=nf_get_att_text(nid,altvar,'units',altunits)
368	ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
369
370	if (ctllen .ne. 0) then
371	ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
372	if (ierr.ne.0) then
373	write(,) "Error: failed to load controle"
374	write(,) NF_STRERROR(ierr)
375	stop
376	endif
377	endif ! of if (ctllen .ne. 0)
378
379	!==============================================================================
380	! 2.4. Handle "Time" dimension from input file
381	!==============================================================================
382
383	!==============================================================================
384	! 2.4.0 Read "Time" dimension from input file
385	!==============================================================================
386	ierr=NF_INQ_DIMID(nid,"Time",timedim)
387	if (ierr.NE.NF_NOERR) then
388	write(,) 'ERROR: Dimension <Time> is missing in file'//trim(file)
389	stop ""
390	endif
391	ierr=NF_INQ_VARID(nid,"Time",timevar)
392	if (ierr.NE.NF_NOERR) then
393	write(,) 'ERROR: Field <Time> is missing in file'//trim(file)
394	stop ""
395	endif
396	ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
397	! write(,) "timelen: ",timelen
398
399	! Sanity Check: Does "Time" has a "long_name" attribute
400	! and is it "Solar longitude" ?
401	ierr=nf_get_att_text(nid,timevar,"long_name",long_name)
402	if (ierr.ne.nf_noerr) then
403	! if no attribute "long_name", try "title"
404	ierr=nf_get_att_text(nid,timevar,"title",long_name)
405	endif
406	if ((ierr.EQ.NF_NOERR).and.(index(long_name,"Solar").ne.0)) then
407	! if long_name contains "Solar"
408	write(,) "ERROR: Time axis in input file is in Solar Longitude!"
409	write(,) " localtime requires sols as time axis!"
410	write(,) " Might as well stop here."
411	stop
412	endif
413
414	! allocate time() array and fill it with values from input file
415	allocate(time(timelen),stat=ierr)
416	if (ierr.ne.0) then
417	write(,) "Error: failed to allocate time(timelen)"
418	stop
419	endif
420
421	ierr = NF_GET_VAR_REAL(nid,timevar,time)
422	if (ierr.NE.NF_NOERR) then
423	write(,) "Error , failed to load Time"
424	write(,) NF_STRERROR(ierr)
425	stop
426	endif
427
428
429	!==============================================================================
430	! 2.4.1 Select local times to be saved "Time" in output file
431	!==============================================================================
432	write(,) 'number of local time to be used ?'
433	read(,) nhour
434	allocate(lt_hour(nhour),stat=ierr)
435	if (ierr.ne.0) then
436	write(,) "Error: failed to allocate lt_hour(nhour)"
437	stop
438	endif
439	write(,) 'list of selected local time (0<t<24)'
440	do it=1,nhour
441	read(,) lt_hour(it)
442	end do
443
444	if ((timelen.eq.12).and.(time(1).eq.2).and.(time(timelen).eq.24))then
445	write(,) 'We detect a ""stats"" file'
446	stats=.true.
447	timelen_lt=nhour
448	allocate(lt_out(timelen_lt),stat=ierr)
449	if (ierr.ne.0) then
450	write(,) "Error: failed to allocate lt_hour(nhour)"
451	stop
452	endif
453	do it=1,nhour
454	lt_out(it)=lt_hour(it)/24.
455	end do
456	! We rewrite the time from "stats" from 0 to 1 sol...
457	do it=1,timelen ! loop temps in
458	time(it) = time(it)/24.
459	end do
460	nsol =1
461	else ! case of a diagfi or concatnc file
462	stats=.false.
463	! Number of sol in input file
464	nsol = int(time(timelen)) - int(time(1))
465	timelen_lt=nhour*nsol
466	allocate(lt_out(timelen_lt),stat=ierr)
467	if (ierr.ne.0) then
468	write(,) "Error: failed to allocate lt_hour(nhour)"
469	stop
470	endif
471	i=0
472	do k=1,nsol
473	do it=1,nhour
474	i=i+1
475	lt_out(i)=int(time(1)) + k-1 + lt_hour(it)/24.
476	end do
477	end do
478	end if
479
480	if (nsol.gt.1) then
481	timelen_tot=timelen
482	else
483	! if only 1 sol available, we must add 1 timestep for the interpolation
484	timelen_tot=timelen+1
485	endif
486	allocate(lt_gcm(timelen_tot),stat=ierr)
487	if (ierr.ne.0) then
488	write(,) "Error: failed to allocate lt_gcm(timelen_tot)"
489	stop
490	endif
491	allocate(var_gcm(timelen_tot),stat=ierr)
492	if (ierr.ne.0) then
493	write(,) "Error: failed to allocate var_gcm(timelen_tot)"
494	stop
495	endif
496	allocate(lt_outc(timelen_lt),stat=ierr)
497	if (ierr.ne.0) then
498	write(,) "Error: failed to allocate lt_outc(timelen_lt)"
499	stop
500	endif
501
502	!==============================================================================
503	! 2.4.1.5 Initiate dimension in output file
504	!==============================================================================
505
506
507	! Initialize output file's lat,lon,alt and time dimensions
508	call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,GCM_layers,&
509	altlong_name,altunits,altpositive,&
510	nout,latdimout,londimout,altdimout,timedimout,&
511	layerdimout,timevarout)
512
513	! Initialize output file's aps,bps and phisinit variables
514	call init2(nid,lonlen,latlen,altlen,GCM_layers,&
515	nout,londimout,latdimout,altdimout,&
516	layerdimout)
517
518
519	!==============================================================================
520	! 2.4.2 Write/extend "Time" dimension/values in output file
521	!==============================================================================
522
523	if (stats) then
524
525	do it=1,timelen_lt
526	ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,lt_hour(it))
527	enddo
528	else
529	do it=1,timelen_lt
530	ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,lt_out(it))
531	if (ierr.NE.NF_NOERR) then
532	write(,) "Error , failed to write Time"
533	stop
534	endif
535	enddo
536	end if
537
538	!==============================================================================
539	! 2.5 Read/write variables
540	!==============================================================================
541
542	do j=1,nbvar ! loop on variables to read/write
543
544	!==============================================================================
545	! 2.5.1 Check that variable to be read existe in input file
546	!==============================================================================
547
548	write(,) "variable ",var(j)
549	ierr=nf_inq_varid(nid,var(j),varid)
550	if (ierr.NE.NF_NOERR) then
551	write(,) 'ERROR: Field <',var(j),'> not found in file'//file
552	stop ""
553	endif
554	ierr=nf_inq_varndims(nid,varid,ndim)
555
556	! Check that it is a 3D or 4D variable
557	if (ndim.lt.3) then
558	write(,) "Error:",trim(var(j))," is nor a 3D neither a 4D variable"
559	write(,) "We'll skip that variable..."
560	CYCLE ! go directly to the next variable
561	endif
562	!==============================================================================
563	! 2.5.2 Prepare things in order to read/write the variable
564	!==============================================================================
565
566	! build dim(),corner() and edges() arrays
567	! and allocate var3d() array
568	if (ndim==3) then
569	allocate(var3d(lonlen,latlen,timelen,1))
570	allocate(var3d_lt(lonlen,latlen,timelen_lt,1))
571	dim(1)=londimout
572	dim(2)=latdimout
573	dim(3)=timedimout
574
575	! start indexes (where data values will be written)
576	corner(1)=1
577	corner(2)=1
578	corner(3)=1
579	corner(4)=1
580
581	! length (along dimensions) of block of data to be written
582	edges(1)=lonlen
583	edges(2)=latlen
584	edges(3)=timelen_lt
585	edges(4)=1
586
587	else if (ndim==4) then
588	allocate(var3d(lonlen,latlen,altlen,timelen))
589	allocate(var3d_lt(lonlen,latlen,altlen,timelen_lt))
590	dim(1)=londimout
591	dim(2)=latdimout
592	dim(3)=altdimout
593	dim(4)=timedimout
594
595	! start indexes (where data values will be written)
596	corner(1)=1
597	corner(2)=1
598	corner(3)=1
599	corner(4)=1
600
601	! length (along dimensions) of block of data to be written
602	edges(1)=lonlen
603	edges(2)=latlen
604	edges(3)=altlen
605	edges(4)=timelen_lt
606	endif
607
608	units=" "
609	long_name=" "
610	ierr=nf_get_att_text(nid,varid,"long_name",long_name)
611	if (ierr.ne.nf_noerr) then
612	! if no attribute "long_name", try "title"
613	ierr=nf_get_att_text(nid,varid,"title",long_name)
614	endif
615	ierr=nf_get_att_text(nid,varid,"units",units)
616	call def_var(nout,var(j),long_name,units,ndim,dim,varidout,ierr)
617
618
619
620	!==============================================================================
621	! 2.5.3 Read from input file
622	!==============================================================================
623
624	ierr = NF_GET_VAR_REAL(nid,varid,var3d)
625	miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
626	validr=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
627
628	!==============================================================================
629	! 2.5.4 Read from input file
630	!==============================================================================
631	! Interpolation in local time :
632
633	do ilon=1,lonlen
634	! write(,) 'processing longitude =', lon(ilon)
635	! Local time at each longitude :
636	do it=1,timelen ! loop temps in
637	lt_gcm(it) = time(it) +0.5*lon(ilon)/180.
638	end do
639	if (nsol.eq.1) lt_gcm(timelen+1) = lt_gcm(1)+1
640
641	do it=1,timelen_lt ! loop time local time
642	lt_outc(it)=lt_out(it)
643	! Correction to use same local time previous or following
644	! day where data are missing :
645
646	if(lt_outc(it).gt.lt_gcm(timelen_tot))lt_outc(it)=lt_out(it)-1
647	if(lt_outc(it).lt.lt_gcm(1))lt_outc(it)=lt_out(it)+1
648	end do
649
650	if (ndim==3) then
651	do ilat=1,latlen
652	do it=1,timelen ! loop temps in
653	var_gcm(it) = var3d(ilon,ilat,it,1)
654	end do
655	if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
656	do it=1,timelen_lt ! loop time local time
657	call interpolf(lt_outc(it),var3d_lt(ilon,ilat,it,1),&
658	missing,lt_gcm,var_gcm,timelen_tot)
659	end do
660	enddo
661
662	else if (ndim==4) then
663	do ilat=1,latlen
664	do ialt=1,altlen
665	do it=1,timelen ! loop temps in
666	var_gcm(it) = var3d(ilon,ilat,ialt,it)
667	end do
668	if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
669	do it=1,timelen_lt ! loop time local time
670	call interpolf(lt_outc(it),var3d_lt(ilon,ilat,ialt,it),&
671	missing,lt_gcm,var_gcm,timelen_tot)
672	end do
673	enddo
674	enddo
675	end if
676
677	end do
678
679
680	!==============================================================================
681	! 2.5.5 write (append) to the output file
682	!==============================================================================
683
684	ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3d_lt)
685
686	if (ierr.ne.NF_NOERR) then
687	write(,) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
688	stop ""
689	endif
690
691	! In case there is a "valid_range" attribute
692	! Write "valid_range" and "missing_value" attributes in output file
693	if (miss.eq.NF_NOERR) then
694	call missing_value(nout,varidout,validr,miss,valid_range,missing)
695	endif
696
697	! free var3d() array
698	deallocate(var3d)
699	deallocate(var3d_lt)
700
701	enddo ! of do j=1,nbvar
702
703	deallocate(time)
704	deallocate(lt_gcm)
705	deallocate(lt_out)
706	deallocate(lt_outc)
707	deallocate(var_gcm)
708
709	! Close input file
710	ierr=nf_close(nid)
711
712	! Close output file
713	ierr=NF_CLOSE(nout)
714
715	write(,) "Program completed !"
716
717	contains
718
719	!******************************************************************************
720	Subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,GCM_layers,&
721	altlong_name,altunits,altpositive,&
722	nout,latdimout,londimout,altdimout,timedimout,&
723	layerdimout,timevarout)
724	!==============================================================================
725	! Purpose:
726	! Create and initialize a data file (NetCDF format)
727	!==============================================================================
728	! Remarks:
729	! The NetCDF file (created in this subroutine) remains open
730	!==============================================================================
731
732	implicit none
733
734	include "netcdf.inc" ! NetCDF definitions
735
736	!==============================================================================
737	! Arguments:
738	!==============================================================================
739	character (len=*), intent(in):: filename
740	! filename(): the file's name
741	integer,intent(in) ::latlen,lonlen,altlen,ctllen
742	real, intent(in):: lat(latlen)
743	! lat(): latitude
744	real, intent(in):: lon(lonlen)
745	! lon(): longitude
746	real, intent(in):: alt(altlen)
747	! alt(): altitude
748	real, intent(in):: ctl(ctllen)
749	! ctl(): controle
750	integer,intent(in) :: GCM_layers ! number of GCM layers
751	character (len=*), intent(in) :: altlong_name
752	! altlong_name(): [netcdf] altitude "long_name" attribute
753	character (len=*), intent(in) :: altunits
754	! altunits(): [netcdf] altitude "units" attribute
755	character (len=*), intent(in) :: altpositive
756	! altpositive(): [netcdf] altitude "positive" attribute
757	integer, intent(out):: nout
758	! nout: [netcdf] file ID
759	integer, intent(out):: latdimout
760	! latdimout: [netcdf] lat() (i.e.: latitude) ID
761	integer, intent(out):: londimout
762	! londimout: [netcdf] lon() ID
763	integer, intent(out):: altdimout
764	! altdimout: [netcdf] alt() ID
765	integer, intent(out):: timedimout
766	! timedimout: [netcdf] "Time" ID
767	integer,intent(out) :: layerdimout
768	! layerdimout: [netcdf] "GCM_layers" ID
769	integer, intent(out):: timevarout
770	! timevarout: [netcdf] "Time" (considered as a variable) ID
771
772	!==============================================================================
773	! Local variables:
774	!==============================================================================
775	!integer :: latdim,londim,altdim,timedim
776	integer :: nvarid,ierr
777	integer :: ctldimout
778	! nvarid: [netcdf] ID of a variable
779	! ierr: [netcdf] return error code (from called subroutines)
780
781	!==============================================================================
782	! 1. Create (and open) output file
783	!==============================================================================
784	write(,) "creating "//trim(adjustl(filename))//'...'
785	ierr = NF_CREATE(filename,IOR(NF_CLOBBER,NF_64BIT_OFFSET),nout)
786	! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
787	if (ierr.NE.NF_NOERR) then
788	WRITE(,)'ERROR: Impossible to create the file ',trim(filename)
789	write(,) NF_STRERROR(ierr)
790	stop ""
791	endif
792
793	!==============================================================================
794	! 2. Define/write "dimensions" and get their IDs
795	!==============================================================================
796
797	ierr = NF_DEF_DIM(nout, "latitude", latlen, latdimout)
798	if (ierr.NE.NF_NOERR) then
799	WRITE(,)'initiate: error failed to define dimension <latitude>'
800	write(,) NF_STRERROR(ierr)
801	stop ""
802	endif
803	ierr = NF_DEF_DIM(nout, "longitude", lonlen, londimout)
804	if (ierr.NE.NF_NOERR) then
805	WRITE(,)'initiate: error failed to define dimension <longitude>'
806	write(,) NF_STRERROR(ierr)
807	stop ""
808	endif
809	ierr = NF_DEF_DIM(nout, "altitude", altlen, altdimout)
810	if (ierr.NE.NF_NOERR) then
811	WRITE(,)'initiate: error failed to define dimension <altitude>'
812	write(,) NF_STRERROR(ierr)
813	stop ""
814	endif
815	if (size(ctl).ne.0) then
816	ierr = NF_DEF_DIM(nout, "index", ctllen, ctldimout)
817	if (ierr.NE.NF_NOERR) then
818	WRITE(,)'initiate: error failed to define dimension <index>'
819	write(,) NF_STRERROR(ierr)
820	stop ""
821	endif
822	endif
823	ierr = NF_DEF_DIM(nout, "Time", NF_UNLIMITED, timedimout)
824	if (ierr.NE.NF_NOERR) then
825	WRITE(,)'initiate: error failed to define dimension <Time>'
826	write(,) NF_STRERROR(ierr)
827	stop ""
828	endif
829
830	ierr = NF_DEF_DIM(nout, "GCM_layers", GCM_layers, layerdimout)
831
832	! End netcdf define mode
833	ierr = NF_ENDDEF(nout)
834	if (ierr.NE.NF_NOERR) then
835	WRITE(,)'initiate: error failed to switch out of define mode'
836	write(,) NF_STRERROR(ierr)
837	stop ""
838	endif
839
840	!==============================================================================
841	! 3. Write "Time" and its attributes
842	!==============================================================================
843
844	call def_var(nout,"Time","Time","years since 0000-00-0 00:00:00",1,&
845	(/timedimout/),timevarout,ierr)
846
847	!==============================================================================
848	! 4. Write "latitude" (data and attributes)
849	!==============================================================================
850
851	call def_var(nout,"latitude","latitude","degrees_north",1,&
852	(/latdimout/),nvarid,ierr)
853
854	ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
855	if (ierr.NE.NF_NOERR) then
856	WRITE(,)'initiate: error failed writing variable <latitude>'
857	write(,) NF_STRERROR(ierr)
858	stop ""
859	endif
860
861	!==============================================================================
862	! 4. Write "longitude" (data and attributes)
863	!==============================================================================
864
865	call def_var(nout,"longitude","East longitude","degrees_east",1,&
866	(/londimout/),nvarid,ierr)
867
868
869	ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
870	if (ierr.NE.NF_NOERR) then
871	WRITE(,)'initiate: error failed writing variable <longitude>'
872	write(,) NF_STRERROR(ierr)
873	stop ""
874	endif
875
876	!==============================================================================
877	! 5. Write "altitude" (data and attributes)
878	!==============================================================================
879
880	! Switch to netcdf define mode
881	ierr = NF_REDEF (nout)
882
883	ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,altdimout,nvarid)
884
885	ierr = NF_PUT_ATT_TEXT (nout,nvarid,'long_name',len_trim(adjustl(altlong_name)),adjustl(altlong_name))
886	ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',len_trim(adjustl(altunits)),adjustl(altunits))
887	ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',len_trim(adjustl(altpositive)),adjustl(altpositive))
888
889	! End netcdf define mode
890	ierr = NF_ENDDEF(nout)
891
892	ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
893	if (ierr.NE.NF_NOERR) then
894	WRITE(,)'initiate: error failed writing variable <altitude>'
895	write(,) NF_STRERROR(ierr)
896	stop ""
897	endif
898
899	!==============================================================================
900	! 6. Write "controle" (data and attributes)
901	!==============================================================================
902
903	if (size(ctl).ne.0) then
904	! Switch to netcdf define mode
905	ierr = NF_REDEF (nout)
906
907	ierr = NF_DEF_VAR (nout,"controle",NF_FLOAT,1,ctldimout,nvarid)
908
909	ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",18,"Control parameters")
910
911	! End netcdf define mode
912	ierr = NF_ENDDEF(nout)
913
914	ierr = NF_PUT_VAR_REAL (nout,nvarid,ctl)
915	if (ierr.NE.NF_NOERR) then
916	WRITE(,)'initiate: error failed writing variable <controle>'
917	write(,) NF_STRERROR(ierr)
918	stop ""
919	endif
920	endif
921
922	end Subroutine initiate
923	!******************************************************************************
924	subroutine init2(infid,lonlen,latlen,altlen,GCM_layers,&
925	outfid,londimout,latdimout,altdimout,&
926	layerdimout)
927	!==============================================================================
928	! Purpose:
929	! Copy aps(), bps() and phisinit() from input file to output file
930	!==============================================================================
931	! Remarks:
932	! The NetCDF files must be open
933	!==============================================================================
934
935	implicit none
936
937	include "netcdf.inc" ! NetCDF definitions
938
939	!==============================================================================
940	! Arguments:
941	!==============================================================================
942	integer, intent(in) :: infid ! NetCDF output file ID
943	integer, intent(in) :: lonlen ! # of grid points along longitude
944	integer, intent(in) :: latlen ! # of grid points along latitude
945	integer, intent(in) :: altlen ! # of grid points along altitude
946	integer, intent(in) :: GCM_layers ! # of GCM atmospheric layers
947	integer, intent(in) :: outfid ! NetCDF output file ID
948	integer, intent(in) :: londimout ! longitude dimension ID
949	integer, intent(in) :: latdimout ! latitude dimension ID
950	integer, intent(in) :: altdimout ! altitude dimension ID
951	integer, intent(in) :: layerdimout ! GCM_layers dimension ID
952	!==============================================================================
953	! Local variables:
954	!==============================================================================
955	real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
956	real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
957	integer :: apsid,bpsid,phisinitid
958	integer :: ierr
959	integer :: tmpdimid ! temporary dimension ID
960	integer :: tmpvarid ! temporary variable ID
961	logical :: phis, hybrid ! are "phisinit" and "aps"/"bps" available ?
962
963
964	!==============================================================================
965	! 1. Read data from input file
966	!==============================================================================
967
968	! hybrid coordinate aps
969	!write(,) "aps: altlen=",altlen," GCM_layers=",GCM_layers
970	allocate(aps(GCM_layers),stat=ierr)
971	if (ierr.ne.0) then
972	write(,) "init2: failed to allocate aps!"
973	stop
974	endif
975	ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
976	if (ierr.ne.NF_NOERR) then
977	write(,) "Ooops. Failed to get aps ID. OK."
978	hybrid=.false.
979	else
980	hybrid=.true.
981	ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
982	hybrid=.true.
983	if (ierr.ne.NF_NOERR) then
984	stop "init2 Error: Failed reading aps"
985	endif
986
987	! hybrid coordinate bps
988	! write(,) "bps: altlen=",altlen," GCM_layers=",GCM_layers
989	allocate(bps(GCM_layers),stat=ierr)
990	if (ierr.ne.0) then
991	write(,) "init2: failed to allocate bps!"
992	stop
993	endif
994	ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
995	if (ierr.ne.NF_NOERR) then
996	write(,) "Ooops. Failed to get bps ID. OK."
997	hybrid=.false.
998	else
999	ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
1000	if (ierr.ne.NF_NOERR) then
1001	stop "init2 Error: Failed reading bps"
1002	endif
1003	endif
1004	endif
1005
1006	! ground geopotential phisinit
1007	allocate(phisinit(lonlen,latlen),stat=ierr)
1008	if (ierr.ne.0) then
1009	write(,) "init2: failed to allocate phisinit(lonlen,latlen)"
1010	stop
1011	endif
1012	ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
1013	if (ierr.ne.NF_NOERR) then
1014	write(,) "Failed to get phisinit ID. OK"
1015	phisinit = 0.
1016	phis = .false.
1017	else
1018	ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
1019	if (ierr.ne.NF_NOERR) then
1020	stop "init2 Error: Failed reading phisinit"
1021	endif
1022	phis = .true.
1023	endif
1024
1025
1026	!==============================================================================
1027	! 2. Write
1028	!==============================================================================
1029
1030	!==============================================================================
1031	! 2.2. Hybrid coordinates aps() and bps()
1032	!==============================================================================
1033
1034	IF (hybrid) then
1035	! define aps
1036	call def_var(outfid,"aps","hybrid pressure at midlayers"," ",1,&
1037	(/layerdimout/),apsid,ierr)
1038	if (ierr.ne.NF_NOERR) then
1039	stop "Error: Failed to def_var aps"
1040	endif
1041
1042	! write aps
1043	ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
1044	if (ierr.ne.NF_NOERR) then
1045	stop "Error: Failed to write aps"
1046	endif
1047
1048	! define bps
1049	call def_var(outfid,"bps","hybrid sigma at midlayers"," ",1,&
1050	(/layerdimout/),bpsid,ierr)
1051	if (ierr.ne.NF_NOERR) then
1052	stop "Error: Failed to def_var bps"
1053	endif
1054
1055	! write bps
1056	ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
1057	if (ierr.ne.NF_NOERR) then
1058	stop "Error: Failed to write bps"
1059	endif
1060
1061	deallocate(aps)
1062	deallocate(bps)
1063	ENDIF ! of IF (hybrid)
1064
1065	!==============================================================================
1066	! 2.2. phisinit()
1067	!==============================================================================
1068
1069	IF (phis) THEN
1070
1071	!define phisinit
1072	call def_var(outfid,"phisinit","Ground level geopotential"," ",2,&
1073	(/londimout,latdimout/),phisinitid,ierr)
1074	if (ierr.ne.NF_NOERR) then
1075	stop "Error: Failed to def_var phisinit"
1076	endif
1077
1078	! write phisinit
1079	ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
1080	if (ierr.ne.NF_NOERR) then
1081	stop "Error: Failed to write phisinit"
1082	endif
1083
1084	deallocate(phisinit)
1085	ENDIF ! of IF (phis)
1086
1087
1088	end subroutine init2
1089
1090	!******************************************************************************
1091	subroutine def_var(nid,name,long_name,units,nbdim,dim,nvarid,ierr)
1092	!==============================================================================
1093	! Purpose: Write a variable (i.e: add a variable to a dataset)
1094	! called "name"; along with its attributes "long_name", "units"...
1095	! to a file (following the NetCDF format)
1096	!==============================================================================
1097	! Remarks:
1098	! The NetCDF file must be open
1099	!==============================================================================
1100
1101	implicit none
1102
1103	include "netcdf.inc" ! NetCDF definitions
1104
1105	!==============================================================================
1106	! Arguments:
1107	!==============================================================================
1108	integer, intent(in) :: nid
1109	! nid: [netcdf] file ID #
1110	character (len=*), intent(in) :: name
1111	! name(): [netcdf] variable's name
1112	character (len=*), intent(in) :: long_name
1113	! long_name(): [netcdf] variable's "long_name" attribute
1114	character (len=*), intent(in) :: units
1115	! unit(): [netcdf] variable's "units" attribute
1116	integer, intent(in) :: nbdim
1117	! nbdim: number of dimensions of the variable
1118	integer, dimension(nbdim), intent(in) :: dim
1119	! dim(nbdim): [netcdf] dimension(s) ID(s)
1120	integer, intent(out) :: nvarid
1121	! nvarid: [netcdf] ID # of the variable
1122	integer, intent(out) :: ierr
1123	! ierr: [netcdf] subroutines returned error code
1124
1125	! Switch to netcdf define mode
1126	ierr=NF_REDEF(nid)
1127
1128	! Insert the definition of the variable
1129	ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
1130
1131	! Write the attributes
1132	ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name",len_trim(adjustl(long_name)),adjustl(long_name))
1133	ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
1134
1135	! End netcdf define mode
1136	ierr=NF_ENDDEF(nid)
1137
1138	end subroutine def_var
1139	!******************************************************************************
1140	subroutine missing_value(nout,nvarid,validr,miss,valid_range,missing)
1141	!==============================================================================
1142	! Purpose:
1143	! Write "valid_range" and "missing_value" attributes (of a given
1144	! variable) to a netcdf file
1145	!==============================================================================
1146	! Remarks:
1147	! NetCDF file must be open
1148	! Variable (nvarid) ID must be set
1149	!==============================================================================
1150
1151	implicit none
1152
1153	include "netcdf.inc" ! NetCDF definitions
1154
1155	!==============================================================================
1156	! Arguments:
1157	!==============================================================================
1158	integer, intent(in) :: nout
1159	! nout: [netcdf] file ID #
1160	integer, intent(in) :: nvarid
1161	! varid: [netcdf] variable ID #
1162	integer, intent(in) :: validr
1163	! validr : [netcdf] routines return code for "valid_range" attribute
1164	integer, intent(in) :: miss
1165	! miss : [netcdf] routines return code for "missing_value" attribute
1166	real, dimension(2), intent(in) :: valid_range
1167	! valid_range(2): [netcdf] "valid_range" attribute (min and max)
1168	real, intent(in) :: missing
1169	! missing: [netcdf] "missing_value" attribute
1170
1171	!==============================================================================
1172	! Local variables:
1173	!==============================================================================
1174	integer :: ierr
1175	! ierr: [netcdf] subroutine returned error code
1176	! INTEGER nout,nvarid,ierr
1177	! REAL missing
1178	! REAL valid_range(2)
1179
1180	! Switch to netcdf dataset define mode
1181	ierr = NF_REDEF (nout)
1182	if (ierr.ne.NF_NOERR) then
1183	print*,'missing_value: '
1184	print*, NF_STRERROR(ierr)
1185	endif
1186
1187	! Write valid_range() attribute
1188	if (validr.eq.NF_NOERR) then
1189	ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
1190
1191	if (ierr.ne.NF_NOERR) then
1192	print*,'missing_value: valid_range attribution failed'
1193	print*, NF_STRERROR(ierr)
1194	!write(,) 'NF_NOERR', NF_NOERR
1195	!CALL abort
1196	stop ""
1197	endif
1198	endif
1199	if (miss.eq.NF_NOERR) then
1200	! Write "missing_value" attribute
1201	ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
1202
1203	if (ierr.NE.NF_NOERR) then
1204	print*, 'missing_value: missing value attribution failed'
1205	print*, NF_STRERROR(ierr)
1206	! WRITE(,) 'NF_NOERR', NF_NOERR
1207	! CALL abort
1208	stop ""
1209	endif
1210	endif
1211
1212	! End netcdf dataset define mode
1213	ierr = NF_ENDDEF(nout)
1214
1215	end subroutine missing_value
1216
1217	!*****************************************************************************
1218	subroutine interpolf(x,y,missing,xd,yd,nd)
1219	!==============================================================================
1220	! Purpose:
1221	! Yield y=f(x), where xd() end yd() are arrays of known values,
1222	! using linear interpolation
1223	! If x is not included in the interval spaned by xd(), then y is set
1224	! to a default value 'missing'
1225	! Note:
1226	! Array xd() should contain ordered (either increasing or decreasing) abscissas
1227	!==============================================================================
1228	implicit none
1229	!==============================================================================
1230	! Arguments:
1231	!==============================================================================
1232	real,intent(in) :: x ! abscissa at which interpolation is to be done
1233	real,intent(in) :: missing ! missing value (if no interpolation is performed)
1234	integer :: nd ! size of arrays
1235	real,dimension(nd),intent(in) :: xd ! array of known absissas
1236	real,dimension(nd),intent(in) :: yd ! array of correponding values
1237
1238	real,intent(out) :: y ! interpolated value
1239	!==============================================================================
1240	! Local variables:
1241	!==============================================================================
1242	integer :: i
1243
1244	! default: set y to 'missing'
1245	y=missing
1246
1247	do i=1,nd-1
1248	if (((x.ge.xd(i)).and.(x.le.xd(i+1))).or.&
1249	((x.le.xd(i)).and.(x.ge.xd(i+1)))) then
1250	if ((yd(i).eq.missing).or.(yd(i+1).eq.missing)) then
1251	! cannot perform the interpolation if an encompasing value
1252	! is already set to 'missing'
1253	else
1254	!linear interpolation based on encompassing values
1255	y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i))
1256	endif
1257	exit
1258	endif
1259	enddo
1260
1261
1262	end subroutine interpolf
1263
1264	!******************************************************************************
1265
1266	end program localtime
1267

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

Download in other formats:

Original Format