Context Navigation

concatnc.F90 @ 2546

Last change on this file since 2546 was 2546, checked in by romain.vande, 3 years ago

MARS Dynamico:

Utilitaire concatnc, localtime, zrecast can process Xios output
RV

Property svn:executable set to *

File size: 51.6 KB

Line
1	program concatnc
2
3
4	! ********************************************************
5	! Program to concatenate data from Netcdf "diagfi" files,
6	! outputs from the martian GCM
7	! input : diagfi.nc / concat.nc / stats.nc kind of files
8	! author: Y. Wanherdrick
9	! + aps(), bps() and phisinit() are now also written
10	! to output file (E. Millour, september 2006)
11	! if available (F. Forget, october 2006)
12	! + handle 1D data (EM, January 2007)
13	! + ap(), bp() (F. Forget, February 2008)
14	! + handle the possibility that number of GCM layers (aps, bps
15	! or sigma) may be different from number of vertical levels
16	! of data (which is the case for outputs from 'zrecast')
17	! (EM, April 2010)
18	! + handle absence of ap() and bp() if aps and bps are available
19	! (case of stats file) FF, November 2011
20	! + read and write controle field, if available. TN, October 2013
21	! + possibility to concatenate concat files with a coherent
22	! time axis in the output, and to redistribute the time axis
23	! from an offset given by the user. AB, April 2020
24	! ********************************************************
25
26	implicit none
27
28	include "netcdf.inc" ! NetCDF definitions
29
30	character (len=80), dimension(1000) :: file
31	! file(): input file(s) names(s)
32	character (len=30), dimension(16) :: notconcat
33	! notconcat(): names of the (16) variables that won't be concatenated
34	character (len=50), dimension(:), allocatable :: var
35	! var(): name(s) of variable(s) that will be concatenated
36	character (len=100) :: tmpvar,tmpfile,long_name,units
37	! tmpvar(): used to temporarily store a variable name
38	! tmpfile(): used to temporarily store a file name
39	! long_name(): [netcdf] long_name attribute
40	! units(): [netcdf] units attribute
41	character (len=100) :: filename,vartmp
42	! filename(): output file name
43	! vartmp(): temporary variable name (read from netcdf input file)
44	character (len=50) :: altlong_name,altunits,altpositive
45	! altlong_name(): [netcdf] altitude "long_name" attribute
46	! altunits(): [netcdf] altitude "units" attribute
47	! altpositive(): [netcdf] altitude "positive" attribute
48	character (len=4) :: axis
49	! axis: "sol" or "ls" or "adls"
50	integer :: nid,ierr,miss
51	! nid: [netcdf] file ID #
52	! ierr: [netcdf] subroutine returned error code
53	! miss: [netcdf] subroutine returned error code
54	integer :: i,j,k,inter
55	! for various loops
56	integer :: varid
57	! varid: [netcdf] variable ID #
58	integer :: memolatlen=0,memolonlen=0,memoaltlen=0,memoctllen=0
59	! memolatlen: # of elements of lat(), read from the first input file
60	! memolonlen: # of elements of lon(), read from the first input file
61	! memoaltlen: # of elements of alt(), read from the first input file
62	! memoctllen: # of elements of ctl(), read from the first input file
63	real, dimension(:), allocatable:: lat,lon,alt,ctl,time
64	! lat(): array, stores latitude coordinates
65	! lon(): array, stores longitude coordinates
66	! alt(): array, stores altitude coordinates
67	! ctl(): array, stores controle coordinates
68	! time(): array, stores time coordinates
69	integer :: nbvar,nbfile,nbvarfile,ndim
70	! nbvar: # of variables to concatenate
71	! nbfile: # number of input file(s)
72	! nbvarfile: total # of variables in an input file
73	! ndim: [netcdf] # (3 or 4) of dimensions (for variables)
74	integer :: latdim,londim,altdim,ctldim,timedim
75	! latdim: [netcdf] "latitude" dim ID
76	! londim: [netcdf] "longitude" dim ID
77	! altdim: [netcdf] "altdim" dim ID
78	! ctldim: [netcdf] "ctldim" dim ID
79	! timedim: [netcdf] "timedim" dim ID
80	integer :: gcmlayerdim ! NetCDF dimension ID for # of layers in GCM
81	integer :: latvar,lonvar,altvar,ctlvar,timevar
82	! latvar: [netcdf] ID of "latitude" variable
83	! lonvar: [netcdf] ID of "longitude" variable
84	! altvar: [netcdf] ID of "altitude" variable
85	! ctlvar: [netcdf] ID of "controle" variable
86	! timevar: [netcdf] ID of "Time" variable
87	integer :: latlen,lonlen,altlen,ctllen,timelen
88	! latlen: # of elements of lat() array
89	! lonlen: # of elements of lon() array
90	! altlen: # of elements of alt() array
91	! ctllen: # of elements of ctl() array
92	! timelen: # of elemnets of time() array
93	integer :: GCM_layers ! number of GCM atmospheric layers (may not be
94	! same as altlen if file is an output of zrecast)
95	integer :: nout,latdimout,londimout,altdimout,timedimout,timevarout
96	! nout: [netcdf] output file ID
97	! latdimout: [netcdf] output latitude (dimension) ID
98	! londimout: [netcdf] output longitude (dimension) ID
99	! altdimout: [netcdf] output altitude (dimension) ID
100	! timedimout: [netcdf] output time (dimension) ID
101	! timevarout: [netcdf] ID of output "Time" variable
102	integer :: layerdimout ! NetCDF dimension ID for # of layers in GCM
103	integer :: interlayerdimout ! dimension ID for # of interlayers in GCM
104	integer :: reptime,rep,varidout
105	! reptime: total length of concatenated time() arrays
106	! rep: # number of elements of a time() array to write to the output file
107	! varidout: [netcdf] variable ID # (of a variable to write to the output file)
108	integer :: Nnotconcat,var_ok
109	! Nnotconcat: # of (leading)variables that won't be concatenated
110	! var_ok: flag (0 or 1)
111	integer, dimension(4) :: corner,edges,dim
112	! corner: [netcdf]
113	! edges: [netcdf]
114	! dim: [netcdf]
115	real, dimension(:,:,:,:), allocatable :: var3d
116	! var3D(,,,): 4D array to store a field
117	real :: ctlsol
118	! ctlsol: starting sol value of the file, stored in the variable ctl()
119	real :: previous_last_output_time, output_time
120	! previous_last_output_time: last time stored in concat.nc after treating a file
121	logical :: firstsol
122	! firstsol: true if the first file's starting day must be used in the output
123	real :: startsol
124	! startsol: sol of the firstday as requested by the users (used if firstsol="y")
125	real :: starttimeoffset
126	! starttimeoffset: offset given by the user for the output time axis (0 if ! firstsol="n")
127	real :: missing
128	!PARAMETER(missing=1E+20)
129	! missing: [netcdf] to handle "missing" values when reading/writing files
130	real, dimension(2) :: valid_range
131	! valid_range(2): [netcdf] interval in which a value is considered valid
132	real :: year_day = 669.
133	integer :: XIOS =0
134
135	!==============================================================================
136	! 1.1. Get input file name(s)
137	!==============================================================================
138	write(,)
139	write(,) "-- Program written specifically for NetCDF 'diagfi'/'concat' files"
140	write(,) " from the martian GCM --"
141	write(,)
142	write(,) "which files do you want to use?"
143	write(,) "<Enter> when list ok"
144
145	nbfile=0
146	read(*,'(a50)') tmpfile
147	do While (len_trim(tmpfile).ne.0)
148	nbfile=nbfile+1
149	file(nbfile)=tmpfile
150	read(*,'(a50)') tmpfile
151	enddo
152
153	if(nbfile==0) then
154	write(,) "no file... game over"
155	stop ""
156	endif
157
158	!==============================================================================
159	! 1.2. Ask for starting day value (starttimeoffset)
160	!==============================================================================
161
162	write(,) "Starting time in the concat file ?"
163	write(,) " If you really wish to change it, write it now (not recommanded!)"
164	write(,) " If you do not want to change it, just write 'n' (recommanded)"
165
166	read(,,iostat=ierr) startsol
167
168	firstsol= (ierr.eq.0)
169	if (firstsol) then
170	! if nothing or a character that is not a float is read
171	write(,) "1st input file's starting day will serve as starting day for the outfile."
172	else ! if the user gave a number
173	write(,) "Your input day will serve as starting day for the outfile."
174	endif
175
176	!==============================================================================
177	! 1.3. Open the first input file
178	!==============================================================================
179
180	write(,)
181	write(,) "opening "//trim(file(1))//"..."
182	ierr = NF_OPEN(file(1),NF_NOWRITE,nid)
183	if (ierr.NE.NF_NOERR) then
184	write(,) 'ERROR: Pb opening file '//file(1)
185	stop ""
186	endif
187
188	ierr=NF_INQ_NVARS(nid,nbvarfile)
189	! nbvarfile now set to be the (total) number of variables in file
190
191	!==============================================================================
192	! 1.4. Ask for (output) "Time" axis type
193	!==============================================================================
194
195	write(,)
196	write(,) "Which time axis should be given in the output? (sol/ls/adls)"
197	write(,) " (Warning: to read the result with grads, choose sol)"
198	write(,) " To keep sol as the time axis, and just add Ls as a variable, type adls "
199	read(,) axis
200	! loop as long as axis is neither "sol" nor "ls"
201	do while ((axis/="sol").AND.(axis/="ls").AND.(axis/="adls"))
202	read(,) axis
203	enddo
204
205	! The following variables don't need to be concatenated
206	notconcat(1)='Time'
207	notconcat(2)='controle'
208	notconcat(3)='rlonu'
209	notconcat(4)='latitude'
210	notconcat(5)='longitude'
211	notconcat(6)='altitude'
212	notconcat(7)='rlatv'
213	notconcat(8)='aps'
214	notconcat(9)='bps'
215	notconcat(10)='ap'
216	notconcat(11)='bp'
217	notconcat(12)='cu'
218	notconcat(13)='cv'
219	notconcat(14)='aire'
220	notconcat(15)='phisinit'
221	notconcat(16)='soildepth'
222
223
224	!==============================================================================
225	! 1.5. Get (and check) list of variables to concatenate
226	!==============================================================================
227	write(,)
228	Nnotconcat=0
229	do i=1,nbvarfile
230	ierr=NF_INQ_VARNAME(nid,i,vartmp)
231	! vartmp now contains the "name" of variable of ID # i
232	var_ok=0
233	do inter=1,size(notconcat)
234	if (vartmp.eq.notconcat(inter)) then
235	var_ok=1
236	Nnotconcat=Nnotconcat+1
237	endif
238	enddo
239	if (var_ok.eq.0) write(,) trim(vartmp)
240	enddo
241
242	! Nnotconcat: # of variables that won't be concatenated
243	! nbvarfile: total # of variables in file
244	allocate(var(nbvarfile-Nnotconcat))
245
246
247	write(,)
248	write(,) "which variables do you want to concatenate?"
249	write(,) "all / list of <variables> (separated by <Enter>s)"
250	write(,) "(an empty line , i.e: just <Enter>, implies end of list)"
251	nbvar=0
252	read(*,'(a50)') tmpvar
253	do while ((tmpvar/=' ').AND.(trim(tmpvar)/='all'))
254	nbvar=nbvar+1
255	var(nbvar)=tmpvar
256	read(*,'(a50)') tmpvar
257	enddo
258
259	if (tmpvar=="all") then
260	if (axis=="ls") then
261	nbvar=nbvarfile-Nnotconcat
262	do j=Nnotconcat+1,nbvarfile
263	ierr=nf_inq_varname(nid,j,var(j-Nnotconcat))
264	enddo
265	endif ! of if (axis=="ls")
266	! Variables names from the file are catched
267	nbvar=nbvarfile-Nnotconcat
268	j=1
269	do i=1,nbvarfile
270	ierr=nf_inq_varname(nid,i,vartmp)
271	var_ok=0
272	do inter=1,size(notconcat)
273	if (vartmp.eq.notconcat(inter)) then
274	var_ok=1
275	endif
276	enddo
277	if (var_ok.eq.0) then
278	if (j .gt. nbvar) then
279	write(,) "PROBLEM HERE !", var
280	stop
281	endif
282	var(j) = vartmp
283	write(*,'(a9,1x,i2,1x,a1,1x,a50)') "variable ",j,":",var(j)
284	j=j+1
285	endif
286	enddo
287	else if(nbvar==0) then
288	write(,) "no variable... game over"
289	stop ""
290	endif ! of if (tmpvar=="all")
291
292
293	!==============================================================================
294	! 1.6. Get output file name
295	!==============================================================================
296	filename="concat.nc"
297
298
299	!==============================================================================
300	! 2. Concatenate input file(s) into output file
301	!==============================================================================
302
303	reptime=0
304
305	do i=1,nbfile
306
307	!==============================================================================
308	! 2.1. Open input file
309	!==============================================================================
310
311	if (i/=1) then
312	write(,)
313	write(,) "opening "//trim(file(i))//"..."
314	ierr = NF_OPEN(file(i),NF_NOWRITE,nid)
315	if (ierr.NE.NF_NOERR) then
316	write(,) 'ERROR: Pb opening file '//file(i)
317	write(,) NF_STRERROR(ierr)
318	stop ""
319	endif
320	endif
321
322	!==============================================================================
323	! 2.2. Read (and check) dimensions of variables from input file
324	!==============================================================================
325	!==============================================================================
326	! 2.2.1 Lat, lon, alt, GCM_layers
327	!==============================================================================
328	ierr=NF_INQ_DIMID(nid,"latitude",latdim)
329	ierr=NF_INQ_VARID(nid,"latitude",latvar)
330	if (ierr.NE.NF_NOERR) then
331	ierr=NF_INQ_DIMID(nid,"lat",latdim)
332	ierr=NF_INQ_VARID(nid,"lat",latvar)
333	if (ierr.EQ.NF_NOERR) then
334	write(,) 'XIOS concat spotted'
335	XIOS=1
336	endif
337	endif
338	if (ierr.NE.NF_NOERR) then
339	write(,) 'ERROR: Field <latitude> or <lat> is missing in file '//file(i)
340	stop ""
341	endif
342	ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
343	! write(,) "latlen: ",latlen
344
345	if (XIOS.EQ.0) then
346	ierr=NF_INQ_DIMID(nid,"longitude",londim)
347	ierr=NF_INQ_VARID(nid,"longitude",lonvar)
348	if (ierr.NE.NF_NOERR) then
349	write(,) 'ERROR: Field <longitude> is missing in file '//file(i)
350	stop ""
351	endif
352	ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
353	! write(,) "lonlen: ",lonlen
354	else
355	ierr=NF_INQ_DIMID(nid,"lon",londim)
356	ierr=NF_INQ_VARID(nid,"lon",lonvar)
357	if (ierr.NE.NF_NOERR) then
358	write(,) 'ERROR: Field <lon> is missing in file '//file(i)
359	stop ""
360	endif
361	ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
362	! write(,) "lonlen: ",lonlen
363	endif
364
365	ierr=NF_INQ_DIMID(nid,"altitude",altdim)
366	ierr=NF_INQ_VARID(nid,"altitude",altvar)
367	if (ierr.NE.NF_NOERR) then
368	write(,) 'ERROR: Field <altitude> is missing in file '//file(i)
369	stop ""
370	endif
371	ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
372	! write(,) "altlen: ",altlen
373
374	! load size of aps() or sigma() (in case it is not altlen)
375	! default is that GCM_layers=altlen
376	! but for outputs of zrecast, it may be a different value
377	ierr=NF_INQ_DIMID(nid,"GCM_layers",gcmlayerdim)
378	if (ierr.ne.NF_NOERR) then
379	! didn't find a GCM_layers dimension; therefore we have:
380	GCM_layers=altlen
381	else
382	! load value of GCM_layers
383	ierr=NF_INQ_DIMLEN(nid,gcmlayerdim,GCM_layers)
384	endif
385	! write(,) "GCM_layers=",GCM_layers
386
387	!==============================================================================
388	! 2.2.1 Check presence and read "controle" dimension & variable from input file
389	!==============================================================================
390	if (XIOS.EQ.0) then
391	ierr=NF_INQ_DIMID(nid,"index",ctldim)
392	if (ierr.NE.NF_NOERR) then
393	write(,) 'Dimension <index> is missing in file '//file(i)
394	write(,) "The program continues..."
395	ctllen=0
396	!stop ""
397	else
398	ierr=NF_INQ_VARID(nid,"controle",ctlvar)
399	if (ierr.NE.NF_NOERR) then
400	write(,) 'Field <controle> is missing in file '//file(i)
401	write(,) "The program continues..."
402	ctllen=0
403	!stop ""
404	else
405	ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
406	endif
407	endif
408	else
409	ierr=NF_INQ_DIMID(nid,"controle_axe",ctldim)
410	if (ierr.NE.NF_NOERR) then
411	write(,) 'Dimension <controle_axe> is missing in file '//file(i)
412	write(,) "The program continues..."
413	ctllen=0
414	!stop ""
415	else
416	ierr=NF_INQ_VARID(nid,"controle",ctlvar)
417	if (ierr.NE.NF_NOERR) then
418	write(,) 'Field <controle> is missing in file '//file(i)
419	write(,) "The program continues..."
420	ctllen=0
421	!stop ""
422	else
423	ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
424	endif
425	endif
426	endif
427
428	if (ctllen .ne. 0) then ! variable controle
429	allocate(ctl(ctllen))
430	ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
431	ctlsol = ctl(4)
432	endif
433
434	!==============================================================================
435	! 2.3 Read "Time" dimension & variable from input file
436	!==============================================================================
437	if (XIOS.EQ.0) then
438	ierr=NF_INQ_DIMID(nid,"Time",timedim)
439	if (ierr.NE.NF_NOERR) then
440	write(,) 'ERROR: Dimension <Time> is missing in file '//file(i)
441	stop ""
442	endif
443	ierr=NF_INQ_VARID(nid,"Time",timevar)
444	if (ierr.NE.NF_NOERR) then
445	write(,) 'ERROR: Field <Time> is missing in file '//file(i)
446	stop ""
447	endif
448	ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
449	! write(,) "timelen: ",timelen
450	else
451	ierr=NF_INQ_DIMID(nid,"time_counter",timedim)
452	if (ierr.NE.NF_NOERR) then
453	write(,) 'ERROR: Dimension <time_counter> is missing in file '//file(i)
454	stop ""
455	endif
456	ierr=NF_INQ_VARID(nid,"time_counter",timevar)
457	if (ierr.NE.NF_NOERR) then
458	write(,) 'ERROR: Field <time_counter> is missing in file '//file(i)
459	stop ""
460	endif
461	ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
462	! write(,) "timelen: ",timelen
463	endif
464
465	! allocate time() array and fill it with values from input file
466	allocate(time(timelen))
467
468	ierr = NF_GET_VAR_REAL(nid,timevar,time)
469
470	!==============================================================================
471	! 2.3. Read (and check compatibility of) dimensions of
472	! variables from input file
473	!==============================================================================
474
475	if (i==1) then ! First File ; initialize/allocate
476	memolatlen=latlen
477	memolonlen=lonlen
478	memoaltlen=altlen
479	memoctllen=ctllen
480	allocate(lat(latlen))
481	allocate(lon(lonlen))
482	allocate(alt(altlen))
483	ierr = NF_GET_VAR_REAL(nid,latvar,lat)
484	ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
485	ierr = NF_GET_VAR_REAL(nid,altvar,alt)
486	! Get altitude attributes to handle files with any altitude type
487	ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
488	if (ierr.ne.nf_noerr) then
489	! if no attribute "long_name", try "title"
490	ierr=nf_get_att_text(nid,altvar,"title",altlong_name)
491	endif
492	ierr=nf_get_att_text(nid,altvar,'units',altunits)
493	ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
494	write(,)
495
496	if (ctllen .ne. 0) then
497
498	if (firstsol) then
499	starttimeoffset = float(int(ctlsol)+int(time(1) + ctl(27))) - startsol
500	else ! if firstsol.eq..false.
501	! if there is no user chosen first sol, nevertheless use
502	! starttimeoffset to lower the sol number (same exact season)
503	startsol = modulo(int(int(ctlsol)+int(time(1) + ctl(27))),int(year_day))
504	starttimeoffset = float(int(ctlsol)+int(time(1) + ctl(27))) - startsol
505	endif
506
507	! artificially estimating "previous_last_output_time" for first file read:
508	previous_last_output_time= startsol
509
510	ctl(4) = 0. ! values written in the output
511	ctl(27) = 0. ! file controle variable (resp. day_ini and hour_ini)
512
513	else ! if no variable "controle"
514	if (firstsol) then
515	starttimeoffset = int(time(1)) - startsol
516	else ! if firstsol.eq.false
517	starttimeoffset = 0
518	endif
519
520	! artificially estimating "previous_last_output_time" for first file read:
521	previous_last_output_time=time(1) -(time(2) - time(1)) - starttimeoffset
522	endif
523	write(,) "Original starting day of 1st input file", previous_last_output_time + starttimeoffset
524	write(,) "Starting day of the concat file : ", previous_last_output_time
525
526
527	! Initialize output file's lat,lon,alt and time dimensions
528	write(,)
529	call initiate(filename,lat,lon,alt,ctl,GCM_layers,&
530	altlong_name,altunits,altpositive,&
531	nout,latdimout,londimout,altdimout,timedimout,&
532	layerdimout,interlayerdimout,timevarout,axis)
533	! Initialize output file's aps,bps,ap,bp and phisinit variables
534	call init2(nid,lonlen,latlen,altlen,GCM_layers,&
535	nout,londimout,latdimout,altdimout,&
536	layerdimout,interlayerdimout)
537
538	else ! Not first file to concatenate
539	! Check that latitude,longitude and altitude of current input file
540	! are identical to those of the output file
541	! and that either all the files or no file contain the controle variable
542	if (memolatlen/=latlen) then
543	write(,) "ERROR: Not the same latitude axis"
544	stop ""
545	else if (memolonlen/=lonlen) then
546	write(,) "ERROR: Not the same longitude axis"
547	stop ""
548	else if (memoaltlen/=altlen) then
549	write(,) "ERROR: Not the same altitude axis"
550	stop ""
551	else if (memoctllen/=ctllen) then
552	write(,) "ERROR: Not the same controle axis"
553	stop ""
554	endif
555	endif ! of if (i==1)
556
557
558	!==============================================================================
559	! 2.4 Write/extend "Time" dimension/values in output file
560	!==============================================================================
561
562	rep=0
563
564	do k=reptime+1,reptime+timelen
565	rep=rep+1
566	if (ctllen.ne.0) then
567	output_time=time(rep)+ctlsol-starttimeoffset
568	! correction if the file to concatenate seems "before" previous file
569	! (for instance to concatenate diagfi from the previous year at the enf of a year)
570	print *, "time(rep)", time(rep), rep
571	print *, "ctlsol", ctlsol
572	print *, "starttimeoffset", starttimeoffset
573	print *, "output_time", output_time
574	print *, "previous_last_output_time", previous_last_output_time
575	do while (output_time.lt.previous_last_output_time)
576	output_time = output_time + year_day
577	end do
578	else ! if no control, just add timestep after timestep
579	output_time= previous_last_output_time + (time(2)-time(1)) &
580	+ (time(rep)-time(1))
581	end if
582	if (rep.eq.1) write(,) "Sol", int(output_time)
583
584	ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,output_time)
585	end do
586	! use the last output_time value to update memotime
587	previous_last_output_time = output_time
588
589
590	!==============================================================================
591	! 2.5 Read/write variables
592	!==============================================================================
593
594	do j=1,nbvar ! loop on variables to read/write
595
596	!==============================================================================
597	! 2.5.1 Check that variable to be read existe in input file
598	!==============================================================================
599
600	write(,) "variable ",var(j)
601	ierr=nf_inq_varid(nid,var(j),varid)
602	if (ierr.NE.NF_NOERR) then
603	write(,) 'ERROR: Field <',var(j),'> not found in file '//file(i)
604	stop ""
605	endif
606	ierr=nf_inq_varndims(nid,varid,ndim)
607
608	! Check that it is a 1D, 3D or 4D variable
609	if ((ndim.ne.1).and.(ndim.ne.3).and.(ndim.ne.4)) then
610	write(,) "Error:",trim(var(j))," is not a time-dependent variable,"
611	write(,) "so it can't be concatenated."
612	write(,) "We'll skip that variable..."
613	CYCLE ! go directly to the next variable
614	endif
615
616	!==============================================================================
617	! 2.5.2 Prepare things in order to read/write the variable
618	!==============================================================================
619
620	! build dim(),corner() and edges() arrays
621	! and allocate var3d() array
622	if (ndim==1) then
623	allocate(var3d(timelen,1,1,1))
624	dim(1)=timedimout
625
626	! start indexes (where data values will be written)
627	corner(1)=reptime+1
628	corner(2)=1
629	corner(3)=1
630	corner(4)=1
631
632	! length (along dimensions) of block of data to be written
633	edges(1)=timelen
634	edges(2)=1
635	edges(3)=1
636	edges(4)=1
637
638	else if (ndim==3) then
639	allocate(var3d(lonlen,latlen,timelen,1))
640	dim(1)=londimout
641	dim(2)=latdimout
642	dim(3)=timedimout
643
644	! start indexes (where data values will be written)
645	corner(1)=1
646	corner(2)=1
647	corner(3)=reptime+1
648	corner(4)=1
649
650	! length (along dimensions) of block of data to be written
651	edges(1)=lonlen
652	edges(2)=latlen
653	edges(3)=timelen
654	edges(4)=1
655
656	else if (ndim==4) then
657	allocate(var3d(lonlen,latlen,altlen,timelen))
658	dim(1)=londimout
659	dim(2)=latdimout
660	dim(3)=altdimout
661	dim(4)=timedimout
662
663	! start indexes (where data values will be written)
664	corner(1)=1
665	corner(2)=1
666	corner(3)=1
667	corner(4)=reptime+1
668
669	! length (along dimensions) of block of data to be written
670	edges(1)=lonlen
671	edges(2)=latlen
672	edges(3)=altlen
673	edges(4)=timelen
674	endif
675
676	if (i==1) then ! First call: write some definitions to output file
677	units=" "
678	long_name=" "
679	ierr=nf_get_att_text(nid,varid,"long_name",long_name)
680	if (ierr.ne.nf_noerr) then
681	! if no attribute "long_name", try "title"
682	ierr=nf_get_att_text(nid,varid,"title",long_name)
683	endif
684	ierr=nf_get_att_text(nid,varid,"units",units)
685	call def_var(nout,var(j),long_name,units,ndim,dim,varidout,ierr)
686	else
687	ierr=NF_INQ_VARID(nout,var(j),varidout)
688	endif
689
690	!==============================================================================
691	! 2.5.3 Read from input file and write (append) to the output file
692	!==============================================================================
693
694	ierr = NF_GET_VAR_REAL(nid,varid,var3d)
695	ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3d)
696	miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
697	miss=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
698
699	if (ierr.ne.NF_NOERR) then
700	write(,) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
701	stop ""
702	endif
703
704	! In case there is a "valid_range" attribute
705	! Write "valid_range" and "missing_value" attributes in output file
706	if (miss.eq.NF_NOERR) then
707	call missing_value(nout,varidout,valid_range,missing)
708	endif
709
710	! free var3d() array
711	deallocate(var3d)
712
713	enddo ! of do j=1,nbvar
714
715	! Free time() and compute/store array length (for further concatenations)
716	deallocate(time)
717	reptime=reptime+timelen
718
719	! Deallocate controle if needed
720	if (ctllen.ne.0) deallocate(ctl)
721
722	! Close input file
723	ierr=nf_close(nid)
724
725	enddo ! of i=1,nbfile
726
727	!==============================================================================
728	! 3. If required, change time axis (from sols to Ls)
729	!==============================================================================
730
731	if ((axis=="ls").or.(axis=="adls")) then
732	call change_time_axis(nout,ierr,axis)
733	endif
734
735	! Close output file
736	ierr=NF_CLOSE(nout)
737
738	write(,);write(,) "Program concatnc completed !"
739	contains
740
741	!******************************************************************************
742	subroutine initiate (filename,lat,lon,alt,ctl,GCM_layers,&
743	altlong_name,altunits,altpositive,&
744	nout,latdimout,londimout,altdimout,timedimout,&
745	layerdimout,interlayerdimout,timevarout,axis)
746	!==============================================================================
747	! Purpose:
748	! Create and initialize a data file (NetCDF format)
749	!==============================================================================
750	! Remarks:
751	! The NetCDF file (created in this subroutine) remains open
752	!==============================================================================
753
754	implicit none
755
756	include "netcdf.inc" ! NetCDF definitions
757
758	!==============================================================================
759	! Arguments:
760	!==============================================================================
761	character (len=*), intent(in):: filename
762	! filename(): the file's name
763	real, dimension(:), intent(in):: lat
764	! lat(): latitude
765	real, dimension(:), intent(in):: lon
766	! lon(): longitude
767	real, dimension(:), intent(in):: alt
768	! alt(): altitude
769	real, dimension(:), intent(in):: ctl
770	! ctl(): controle
771	integer,intent(in) :: GCM_layers ! number of GCM layers
772	character (len=*), intent(in) :: altlong_name
773	! altlong_name(): [netcdf] altitude "long_name" attribute
774	character (len=*), intent(in) :: altunits
775	! altunits(): [netcdf] altitude "units" attribute
776	character (len=*), intent(in) :: altpositive
777	! altpositive(): [netcdf] altitude "positive" attribute
778	integer, intent(out):: nout
779	! nout: [netcdf] file ID
780	integer, intent(out):: latdimout
781	! latdimout: [netcdf] lat() (i.e.: latitude) ID
782	integer, intent(out):: londimout
783	! londimout: [netcdf] lon() ID
784	integer, intent(out):: altdimout
785	! altdimout: [netcdf] alt() ID
786	integer, intent(out):: timedimout
787	! timedimout: [netcdf] "Time" ID
788	integer,intent(out) :: layerdimout
789	! layerdimout: [netcdf] "GCM_layers" ID
790	integer,intent(out) :: interlayerdimout
791	! layerdimout: [netcdf] "GCM_layers+1" ID
792	integer, intent(out):: timevarout
793	! timevarout: [netcdf] "Time" (considered as a variable) ID
794	character (len=4),intent(in) :: axis
795	! axis: "ls" or "sol"
796
797	!==============================================================================
798	! Local variables:
799	!==============================================================================
800	!integer :: latdim,londim,altdim,timedim
801	integer :: nvarid,ierr
802	integer :: ctldimout
803	! nvarid: [netcdf] ID of a variable
804	! ierr: [netcdf] return error code (from called subroutines)
805
806	!==============================================================================
807	! 1. Create (and open) output file
808	!==============================================================================
809	write(,) "creating "//trim(adjustl(filename))//'...'
810	ierr = NF_CREATE(filename,IOR(NF_CLOBBER,NF_64BIT_OFFSET),nout)
811	! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
812	if (ierr.NE.NF_NOERR) then
813	WRITE(,)'ERROR: Impossible to create the file.'
814	stop ""
815	endif
816
817	!==============================================================================
818	! 2. Define/write "dimensions" and get their IDs
819	!==============================================================================
820
821	ierr = NF_DEF_DIM(nout, "latitude", size(lat), latdimout)
822	ierr = NF_DEF_DIM(nout, "longitude", size(lon), londimout)
823	ierr = NF_DEF_DIM(nout, "altitude", size(alt), altdimout)
824	if (ctllen.ne.0) ierr = NF_DEF_DIM(nout, "index", size(ctl), ctldimout)
825	ierr = NF_DEF_DIM(nout, "Time", NF_UNLIMITED, timedimout)
826	ierr = NF_DEF_DIM(nout, "GCM_layers", GCM_layers, layerdimout)
827	ierr = NF_DEF_DIM(nout, "GCM_interlayers",GCM_layers+1,interlayerdimout)
828
829	! End netcdf define mode
830	ierr = NF_ENDDEF(nout)
831
832	!==============================================================================
833	! 3. Write "Time" and its attributes
834	!==============================================================================
835	if (axis=="ls") then
836	call def_var(nout,"Time","Ls (Solar Longitude)","degrees",1,&
837	(/timedimout/),timevarout,ierr)
838	else
839	call def_var(nout,"Time","Time","days since 0000-00-0 00:00:00",1,&
840	(/timedimout/),timevarout,ierr)
841	endif
842	!==============================================================================
843	! 4. Write "latitude" (data and attributes)
844	!==============================================================================
845
846	call def_var(nout,"latitude","latitude","degrees_north",1,&
847	(/latdimout/),nvarid,ierr)
848
849	ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
850
851	!==============================================================================
852	! 4. Write "longitude" (data and attributes)
853	!==============================================================================
854
855	call def_var(nout,"longitude","East longitude","degrees_east",1,&
856	(/londimout/),nvarid,ierr)
857
858	ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
859
860	!==============================================================================
861	! 5. Write "altitude" (data and attributes)
862	!==============================================================================
863
864	! Switch to netcdf define mode
865	ierr = NF_REDEF (nout)
866
867	ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,altdimout,nvarid)
868
869	ierr = NF_PUT_ATT_TEXT (nout,nvarid,'long_name',len_trim(adjustl(altlong_name)),adjustl(altlong_name))
870	ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',len_trim(adjustl(altunits)),adjustl(altunits))
871	ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',len_trim(adjustl(altpositive)),adjustl(altpositive))
872
873	! End netcdf define mode
874	ierr = NF_ENDDEF(nout)
875
876	ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
877
878	!==============================================================================
879	! 6. Write "controle" (data and attributes)
880	!==============================================================================
881
882	if (ctllen.ne.0) then
883	! Switch to netcdf define mode
884	ierr = NF_REDEF (nout)
885
886	ierr = NF_DEF_VAR (nout,"controle",NF_FLOAT,1,ctldimout,nvarid)
887
888	ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",18,"Control parameters")
889
890	! End netcdf define mode
891	ierr = NF_ENDDEF(nout)
892
893	ierr = NF_PUT_VAR_REAL (nout,nvarid,ctl)
894	endif
895
896	end Subroutine initiate
897	!******************************************************************************
898	subroutine init2(infid,lonlen,latlen,altlen,GCM_layers, &
899	outfid,londimout,latdimout,altdimout, &
900	layerdimout,interlayerdimout)
901	!==============================================================================
902	! Purpose:
903	! Copy ap(), bp(), aps(), bps(), aire() and phisinit()
904	! from input file to output file
905	!==============================================================================
906	! Remarks:
907	! The NetCDF files must be open
908	!==============================================================================
909
910	implicit none
911
912	include "netcdf.inc" ! NetCDF definitions
913
914	!==============================================================================
915	! Arguments:
916	!==============================================================================
917	integer, intent(in) :: infid ! NetCDF output file ID
918	integer, intent(in) :: lonlen ! # of grid points along longitude
919	integer, intent(in) :: latlen ! # of grid points along latitude
920	integer, intent(in) :: altlen ! # of grid points along altitude
921	integer, intent(in) :: GCM_layers ! # of GCM atmospheric layers
922	integer, intent(in) :: outfid ! NetCDF output file ID
923	integer, intent(in) :: londimout ! longitude dimension ID
924	integer, intent(in) :: latdimout ! latitude dimension ID
925	integer, intent(in) :: altdimout ! altitude dimension ID
926	integer, intent(in) :: layerdimout ! GCM_layers dimension ID
927	integer, intent(in) :: interlayerdimout ! GCM_layers+1 dimension ID
928	!==============================================================================
929	! Local variables:
930	!==============================================================================
931	real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
932	real,dimension(:),allocatable :: ap,bp ! hybrid vertical coordinates
933	real,dimension(:),allocatable :: sigma ! sigma levels
934	real,dimension(:,:),allocatable :: aire ! mesh areas
935	real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
936	integer :: apsid,bpsid,sigmaid,phisinitid
937	integer :: apid,bpid
938	integer :: ierr
939	integer :: tmpvarid ! temporary variable ID
940	logical :: area ! is "aire" available ?
941	logical :: phis ! is "phisinit" available ?
942	logical :: hybrid ! are "aps" and "bps" available ?
943	logical :: apbp ! are "ap" and "bp" available ?
944	logical :: sig ! is "sigma" available ?
945
946	!==============================================================================
947	! 1. Read data from input file
948	!==============================================================================
949
950	! hybrid coordinate aps
951	!write(,) "aps: altlen=",altlen," GCM_layers=",GCM_layers
952	allocate(aps(GCM_layers),stat=ierr)
953	if (ierr.ne.0) then
954	write(,) "init2: failed to allocate aps!"
955	stop
956	endif
957	ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
958	if (ierr.ne.NF_NOERR) then
959	write(,) "Ooops. Failed to get aps ID. OK, will look for sigma coord."
960	hybrid=.false.
961	else
962	ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
963	hybrid=.true.
964	if (ierr.ne.NF_NOERR) then
965	stop "init2 Error: Failed reading aps"
966	endif
967
968	! hybrid coordinate bps
969	! write(,) "bps: altlen=",altlen," GCM_layers=",GCM_layers
970	allocate(bps(GCM_layers),stat=ierr)
971	if (ierr.ne.0) then
972	write(,) "init2: failed to allocate bps!"
973	stop
974	endif
975	ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
976	if (ierr.ne.NF_NOERR) then
977	stop "init2 Error: Failed to get bps ID."
978	endif
979	ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
980	if (ierr.ne.NF_NOERR) then
981	stop "init2 Error: Failed reading bps"
982	endif
983	endif
984
985	! hybrid coordinate ap
986	allocate(ap(GCM_layers+1),stat=ierr)
987	if (ierr.ne.0) then
988	write(,) "init2: failed to allocate ap!"
989	stop
990	endif
991	ierr=NF_INQ_VARID(infid,"ap",tmpvarid)
992	if (ierr.ne.NF_NOERR) then
993	write(,) "Ooops. Failed to get ap ID. OK."
994	apbp=.false.
995	else
996	ierr=NF_GET_VAR_REAL(infid,tmpvarid,ap)
997	apbp=.true.
998	if (ierr.ne.NF_NOERR) then
999	stop "Error: Failed reading ap"
1000	endif
1001	endif
1002
1003	! hybrid coordinate bp
1004	allocate(bp(GCM_layers+1),stat=ierr)
1005	if (ierr.ne.0) then
1006	write(,) "init2: failed to allocate bp!"
1007	stop
1008	endif
1009	ierr=NF_INQ_VARID(infid,"bp",tmpvarid)
1010	if (ierr.ne.NF_NOERR) then
1011	write(,) "Ooops. Failed to get bp ID. OK."
1012	apbp=.false.
1013	else
1014	ierr=NF_GET_VAR_REAL(infid,tmpvarid,bp)
1015	apbp=.true.
1016	if (ierr.ne.NF_NOERR) then
1017	stop "Error: Failed reading bp"
1018	endif
1019	endif
1020
1021	! sigma levels (if any)
1022	if (.not.hybrid) then
1023	allocate(sigma(GCM_layers),stat=ierr)
1024	if (ierr.ne.0) then
1025	write(,) "init2: failed to allocate sigma"
1026	stop
1027	endif
1028	ierr=NF_INQ_VARID(infid,"sigma",tmpvarid)
1029	if (ierr.ne.NF_NOERR) then
1030	write(,) "Ooops. Failed to get sigma ID. OK."
1031	sig=.false.
1032	else
1033	ierr=NF_GET_VAR_REAL(infid,tmpvarid,sigma)
1034	sig=.true.
1035	if (ierr.ne.NF_NOERR) then
1036	stop "init2 Error: Failed reading sigma"
1037	endif
1038	endif
1039	endif ! of if (.not.hybrid)
1040
1041	! mesh area
1042	allocate(aire(lonlen,latlen),stat=ierr)
1043	if (ierr.ne.0) then
1044	write(,) "init2: failed to allocate aire!"
1045	stop
1046	endif
1047	ierr=NF_INQ_VARID(infid,"aire",tmpvarid)
1048	if (ierr.ne.NF_NOERR) then
1049	write(,)"init2 warning: Failed to get aire ID."
1050	area = .false.
1051	else
1052	ierr=NF_GET_VAR_REAL(infid,tmpvarid,aire)
1053	if (ierr.ne.NF_NOERR) then
1054	stop "init2 Error: Failed reading aire"
1055	endif
1056	area = .true.
1057	endif
1058
1059	! ground geopotential phisinit
1060	allocate(phisinit(lonlen,latlen),stat=ierr)
1061	if (ierr.ne.0) then
1062	write(,) "init2: failed to allocate phisinit!"
1063	stop
1064	endif
1065	ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
1066	if (ierr.ne.NF_NOERR) then
1067	write(,)"init2 warning: Failed to get phisinit ID."
1068	phis = .false.
1069	else
1070	ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
1071	if (ierr.ne.NF_NOERR) then
1072	stop "init2 Error: Failed reading phisinit"
1073	endif
1074	phis = .true.
1075	endif
1076
1077	!==============================================================================
1078	! 2. Write
1079	!==============================================================================
1080
1081	!==============================================================================
1082	! 2.2. Hybrid coordinates ap() , bp(), aps() and bps()
1083	!==============================================================================
1084	if(hybrid) then
1085	! define aps
1086	call def_var(nout,"aps","hybrid pressure at midlayers"," ",1,&
1087	(/layerdimout/),apsid,ierr)
1088	if (ierr.ne.NF_NOERR) then
1089	stop "init2 Error: Failed to def_var aps"
1090	endif
1091
1092	! write aps
1093	ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
1094	if (ierr.ne.NF_NOERR) then
1095	stop "init2 Error: Failed to write aps"
1096	endif
1097
1098	! define bps
1099	call def_var(nout,"bps","hybrid sigma at midlayers"," ",1,&
1100	(/layerdimout/),bpsid,ierr)
1101	if (ierr.ne.NF_NOERR) then
1102	stop "init2 Error: Failed to def_var bps"
1103	endif
1104
1105	! write bps
1106	ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
1107	if (ierr.ne.NF_NOERR) then
1108	stop "init2 Error: Failed to write bps"
1109	endif
1110
1111	if (apbp) then
1112	! define ap
1113	call def_var(nout,"ap","hybrid sigma at interlayers"," ",1,&
1114	(/interlayerdimout/),apid,ierr)
1115	if (ierr.ne.NF_NOERR) then
1116	stop "Error: Failed to def_var ap"
1117	endif
1118
1119	! write ap
1120	ierr=NF_PUT_VAR_REAL(outfid,apid,ap)
1121	if (ierr.ne.NF_NOERR) then
1122	stop "Error: Failed to write ap"
1123	endif
1124
1125	! define bp
1126	call def_var(nout,"bp","hybrid sigma at interlayers"," ",1,&
1127	(/interlayerdimout/),bpid,ierr)
1128	if (ierr.ne.NF_NOERR) then
1129	stop "Error: Failed to def_var bp"
1130	endif
1131
1132	! write bp
1133	ierr=NF_PUT_VAR_REAL(outfid,bpid,bp)
1134	if (ierr.ne.NF_NOERR) then
1135	stop "Error: Failed to write bp"
1136	endif
1137	endif ! of if (apbp)
1138
1139	else if (sig) then
1140	! define sigma
1141	call def_var(nout,"sigma","sigma at midlayers"," ",1,&
1142	(/layerdimout/),sigmaid,ierr)
1143	if (ierr.ne.NF_NOERR) then
1144	stop "init2 Error: Failed to def_var sigma"
1145	endif
1146	! write sigma
1147	ierr=NF_PUT_VAR_REAL(outfid,sigmaid,sigma)
1148	if (ierr.ne.NF_NOERR) then
1149	stop "init2 Error: Failed to write sigma"
1150	endif
1151	endif ! of if (hybrid)
1152
1153	!==============================================================================
1154	! 2.2. aire() and phisinit()
1155	!==============================================================================
1156
1157	if (area) then
1158	! define aire
1159	call def_var(nout,"aire","Mesh area","m2",2,&
1160	(/londimout,latdimout/),tmpvarid,ierr)
1161	if (ierr.ne.NF_NOERR) then
1162	stop "init2 Error: Failed to def_var aire"
1163	endif
1164
1165	! write aire
1166	ierr=NF_PUT_VAR_REAL(outfid,tmpvarid,aire)
1167	if (ierr.ne.NF_NOERR) then
1168	stop "init2 Error: Failed to write aire"
1169	endif
1170	endif ! of if (area)
1171
1172	IF (phis) THEN
1173
1174	!define phisinit
1175	call def_var(nout,"phisinit","Ground level geopotential"," ",2,&
1176	(/londimout,latdimout/),phisinitid,ierr)
1177	if (ierr.ne.NF_NOERR) then
1178	stop "init2 Error: Failed to def_var phisinit"
1179	endif
1180
1181	! write phisinit
1182	ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
1183	if (ierr.ne.NF_NOERR) then
1184	stop "init2 Error: Failed to write phisinit"
1185	endif
1186
1187	ENDIF ! of IF (phis)
1188
1189
1190	! Cleanup
1191	if (allocated(aps)) deallocate(aps)
1192	if (allocated(bps)) deallocate(bps)
1193	if (allocated(ap)) deallocate(ap)
1194	if (allocated(bp)) deallocate(bp)
1195	if (allocated(sigma)) deallocate(sigma)
1196	if (allocated(phisinit)) deallocate(phisinit)
1197	if (allocated(aire)) deallocate(aire)
1198
1199	end subroutine init2
1200	!******************************************************************************
1201	subroutine def_var(nid,name,long_name,units,nbdim,dim,nvarid,ierr)
1202	!==============================================================================
1203	! Purpose: Write a variable (i.e: add a variable to a dataset)
1204	! called "name"; along with its attributes "long_name", "units"...
1205	! to a file (following the NetCDF format)
1206	!==============================================================================
1207	! Remarks:
1208	! The NetCDF file must be open
1209	!==============================================================================
1210
1211	implicit none
1212
1213	include "netcdf.inc" ! NetCDF definitions
1214
1215	!==============================================================================
1216	! Arguments:
1217	!==============================================================================
1218	integer, intent(in) :: nid
1219	! nid: [netcdf] file ID #
1220	character (len=*), intent(in) :: name
1221	! name(): [netcdf] variable's name
1222	character (len=*), intent(in) :: long_name
1223	! long_name(): [netcdf] variable's "long_name" attribute
1224	character (len=*), intent(in) :: units
1225	! unit(): [netcdf] variable's "units" attribute
1226	integer, intent(in) :: nbdim
1227	! nbdim: number of dimensions of the variable
1228	integer, dimension(nbdim), intent(in) :: dim
1229	! dim(nbdim): [netcdf] dimension(s) ID(s)
1230	integer, intent(out) :: nvarid
1231	! nvarid: [netcdf] ID # of the variable
1232	integer, intent(out) :: ierr
1233	! ierr: [netcdf] subroutines returned error code
1234
1235	! Switch to netcdf define mode
1236	ierr=NF_REDEF(nid)
1237
1238	! Insert the definition of the variable
1239	ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
1240
1241	! Write the attributes
1242	ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name",len_trim(adjustl(long_name)),adjustl(long_name))
1243	ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
1244
1245	! End netcdf define mode
1246	ierr=NF_ENDDEF(nid)
1247
1248	end subroutine def_var
1249	!******************************************************************************
1250	subroutine change_time_axis(nid,ierr,axis)
1251	!==============================================================================
1252	! Purpose:
1253	! Read "time" variable from a dataset, convert it from "sol" (martian
1254	! days) to "Ls" (solar longitude) and write the result back in the file
1255	!==============================================================================
1256	! Remarks:
1257	! The NetCDF file must be opened before this subroutine is called
1258	!==============================================================================
1259
1260	implicit none
1261
1262	include "netcdf.inc" ! NetCDF definitions
1263
1264	!==============================================================================
1265	! Arguments:
1266	!==============================================================================
1267	integer, intent(in) :: nid
1268	! nid: [netcdf] file ID
1269	integer, intent(out) :: ierr
1270	! ierr: [netcdf] return error code
1271	character (len=4),intent(in) :: axis
1272	! axis: "ls" or "sol"
1273
1274	!==============================================================================
1275	! Local variables:
1276	!==============================================================================
1277	integer :: nvarid
1278	! nvarid: ID of the "Time" variable
1279	integer :: timelen
1280	! timelen: size of the arrays
1281	integer :: timedim
1282	! timedim: ID of the "Time" dimension
1283	integer i
1284
1285	real, dimension(:), allocatable :: time,ls
1286	! time(): time, given in sols
1287	! ls(): time, given in Ls (solar longitude)
1288
1289	!==============================================================================
1290	! 1. Read
1291	!==============================================================================
1292
1293	ierr=NF_INQ_DIMID(nid,"Time",timedim)
1294	ierr=NF_INQ_VARID(nid,"Time",nvarid)
1295	if (ierr.NE.NF_NOERR) then
1296	write(,) 'ERROR in change_time_axis: Field <Time> not found'
1297	print*, NF_STRERROR(ierr)
1298	stop ""
1299	endif
1300
1301	ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
1302	allocate(time(timelen),ls(timelen))
1303	ierr = NF_GET_VAR_REAL(nid,nvarid,time)
1304	if (ierr.ne.NF_NOERR) then
1305	write(,) "ERROR in change_time_axis: Failed to load Time"
1306	stop
1307	endif
1308
1309	!==============================================================================
1310	! 2. Convert sols to Ls
1311	!==============================================================================
1312
1313	do i=1,timelen
1314	call sol2ls(time(i),ls(i))
1315	enddo
1316
1317	!==============================================================================
1318	! 2.1 Check if there are not jumps in Ls (rounding problems in sol2ls)
1319	!==============================================================================
1320
1321	do i=1,timelen-1
1322	if ((ls(i+1)-ls(i)) > 350) then
1323	write(,) "+ 360 deg Ls jump solved:", ls(i), ls(i+1), "at timestep", i
1324	ls(i+1) = ls(i+1) - 360
1325	write(,) " corrected to now be:", ls(i), ls(i+1)
1326	else if ((ls(i)-ls(i+1)) > 350) then
1327	write(,) "- 360 deg Ls jump solved:", ls(i), ls(i+1), "at timestep", i
1328	ls(i+1) = ls(i+1) + 360
1329	write(,) " corrected to now be:", ls(i), ls(i+1)
1330	endif
1331	enddo
1332
1333	!==============================================================================
1334	! 3. Write
1335	!==============================================================================
1336
1337	if (axis.eq."ls") then
1338	write(,) "Converting the time axis from sol to Ls..."
1339	else ! if (axis.eq."adls")
1340	write(,) "Adding Ls as a 1D time variable"
1341	call def_var(nout,"Ls","Ls (Solar Longitude)","degrees",1, (/timedimout/),nvarid,ierr)
1342	end if
1343
1344	ierr = NF_PUT_VAR_REAL(nid,nvarid,ls)
1345	if (ierr.ne.NF_NOERR) then
1346	write(,) "ERROR in change_time_axis: Failed to write Ls"
1347	stop
1348	endif
1349
1350	end subroutine change_time_axis
1351	!******************************************************************************
1352	subroutine sol2ls(sol,Ls)
1353	!==============================================================================
1354	! Purpose:
1355	! Convert a date/time, given in sol (martian day),
1356	! into solar longitude date/time, in Ls (in degrees),
1357	! where sol=0 is (by definition) the northern hemisphere
1358	! spring equinox (where Ls=0).
1359	!==============================================================================
1360	! Notes:
1361	! Even though "Ls" is cyclic, if "sol" is greater than N (martian) year,
1362	! "Ls" will be increased by N*360
1363	! Won't work as expected if sol is negative (then again,
1364	! why would that ever happen?)
1365	!==============================================================================
1366
1367	implicit none
1368
1369	!==============================================================================
1370	! Arguments:
1371	!==============================================================================
1372	real,intent(in) :: sol
1373	real,intent(out) :: Ls
1374
1375	!==============================================================================
1376	! Local variables:
1377	!==============================================================================
1378	real year_day,peri_day,timeperi,e_elips,twopi,degrad
1379	data year_day /669./ ! # of sols in a martian year
1380	data peri_day /485.0/
1381	data timeperi /1.9082314/
1382	data e_elips /0.093358/
1383	data twopi /6.2831853/ ! 2.*pi
1384	data degrad /57.2957795/ ! pi/180
1385
1386	real zanom,xref,zx0,zdx,zteta,zz
1387
1388	integer count_years
1389	integer iter
1390
1391	!==============================================================================
1392	! 1. Compute Ls
1393	!==============================================================================
1394
1395	zz=(sol-peri_day)/year_day
1396	zanom=twopi*(zz-nint(zz))
1397	xref=abs(zanom)
1398
1399	! The equation zx0 - e * sin (zx0) = xref, solved by Newton
1400	zx0=xref+e_elips*sin(xref)
1401	do iter=1,20 ! typically, 2 or 3 iterations are enough
1402	zdx=-(zx0-e_elipssin(zx0)-xref)/(1.-e_elipscos(zx0))
1403	zx0=zx0+zdx
1404	if(abs(zdx).le.(1.e-7)) then
1405	! write(,)'iter:',iter,' \|zdx\|:',abs(zdx)
1406	exit
1407	endif
1408	enddo
1409
1410	if(zanom.lt.0.) zx0=-zx0
1411
1412	zteta=2.atan(sqrt((1.+e_elips)/(1.-e_elips))tan(zx0/2.))
1413	Ls=zteta-timeperi
1414
1415	if(Ls.lt.0.) then
1416	Ls=Ls+twopi
1417	else
1418	if(Ls.gt.twopi) then
1419	Ls=Ls-twopi
1420	endif
1421	endif
1422
1423	Ls=degrad*Ls
1424	! Ls is now in degrees
1425
1426	!==============================================================================
1427	! 1. Account for (eventual) years included in input date/time sol
1428	!==============================================================================
1429
1430	count_years=0 ! initialize
1431	zz=sol ! use "zz" to store (and work on) the value of sol
1432	do while (zz.ge.year_day)
1433	count_years=count_years+1
1434	zz=zz-year_day
1435	enddo
1436
1437	! Add 360 degrees to Ls for every year
1438	if (count_years.ne.0) then
1439	Ls=Ls+360.*count_years
1440	endif
1441
1442
1443	end subroutine sol2ls
1444	!******************************************************************************
1445	subroutine missing_value(nout,nvarid,valid_range,missing)
1446	!==============================================================================
1447	! Purpose:
1448	! Write "valid_range" and "missing_value" attributes (of a given
1449	! variable) to a netcdf file
1450	!==============================================================================
1451	! Remarks:
1452	! NetCDF file must be open
1453	! Variable (nvarid) ID must be set
1454	!==============================================================================
1455
1456	implicit none
1457
1458	include "netcdf.inc" ! NetCDF definitions
1459
1460	!==============================================================================
1461	! Arguments:
1462	!==============================================================================
1463	integer, intent(in) :: nout
1464	! nout: [netcdf] file ID #
1465	integer, intent(in) :: nvarid
1466	! varid: [netcdf] variable ID #
1467	real, dimension(2), intent(in) :: valid_range
1468	! valid_range(2): [netcdf] "valid_range" attribute (min and max)
1469	real, intent(in) :: missing
1470	! missing: [netcdf] "missing_value" attribute
1471
1472	!==============================================================================
1473	! Local variables:
1474	!==============================================================================
1475	integer :: ierr
1476	! ierr: [netcdf] subroutine returned error code
1477	! INTEGER nout,nvarid,ierr
1478	! REAL missing
1479	! REAL valid_range(2)
1480
1481	! Switch to netcdf dataset define mode
1482	ierr = NF_REDEF (nout)
1483	if (ierr.ne.NF_NOERR) then
1484	print*,'missing_value: '
1485	print*, NF_STRERROR(ierr)
1486	endif
1487
1488	! Write valid_range() attribute
1489	ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
1490
1491	if (ierr.ne.NF_NOERR) then
1492	print*,'missing_value: valid_range attribution failed'
1493	print*, NF_STRERROR(ierr)
1494	!write(,) 'NF_NOERR', NF_NOERR
1495	!CALL abort
1496	stop ""
1497	endif
1498
1499	! Write "missing_value" attribute
1500	ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
1501
1502	if (ierr.NE.NF_NOERR) then
1503	print*, 'missing_value: missing value attribution failed'
1504	print*, NF_STRERROR(ierr)
1505	! WRITE(,) 'NF_NOERR', NF_NOERR
1506	! CALL abort
1507	stop ""
1508	endif
1509
1510	! End netcdf dataset define mode
1511	ierr = NF_ENDDEF(nout)
1512
1513	end subroutine missing_value
1514	!******************************************************************************
1515
1516	end program concatnc

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

Download in other formats:

Original Format