Context Navigation

concatnc.F90 @ 2908

Last change on this file since 2908 was 2590, checked in by lrossi, 3 years ago

MARS GCM

Utils: fixed a bug in concatnc. The program now properly carries over missing
values from the input files into the output file.

Property svn:executable set to *

File size: 52.3 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=100), 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=100), 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,validr
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(*,'(a100)') tmpvar
253	do while ((tmpvar/=' ').AND.(trim(tmpvar)/='all'))
254	nbvar=nbvar+1
255	var(nbvar)=tmpvar
256	read(*,'(a100)') 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	validr=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
697	miss=NF_GET_ATT_REAL(nid,varid,"missing_value",[missing])
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 "missing_value" 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,validr,miss,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	write(,) "init2 Error: Failed reading aps"
966	stop
967	endif
968
969	! hybrid coordinate bps
970	! write(,) "bps: altlen=",altlen," GCM_layers=",GCM_layers
971	allocate(bps(GCM_layers),stat=ierr)
972	if (ierr.ne.0) then
973	write(,) "init2: failed to allocate bps!"
974	stop
975	endif
976	ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
977	if (ierr.ne.NF_NOERR) then
978	write(,) "init2 Error: Failed to get bps ID."
979	stop
980	endif
981	ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
982	if (ierr.ne.NF_NOERR) then
983	write(,) "init2 Error: Failed reading bps"
984	stop
985	endif
986	endif
987
988	! hybrid coordinate ap
989	allocate(ap(GCM_layers+1),stat=ierr)
990	if (ierr.ne.0) then
991	write(,) "init2: failed to allocate ap!"
992	stop
993	endif
994	ierr=NF_INQ_VARID(infid,"ap",tmpvarid)
995	if (ierr.ne.NF_NOERR) then
996	write(,) "Ooops. Failed to get ap ID. OK."
997	apbp=.false.
998	else
999	ierr=NF_GET_VAR_REAL(infid,tmpvarid,ap)
1000	apbp=.true.
1001	if (ierr.ne.NF_NOERR) then
1002	write(,) "Error: Failed reading ap"
1003	stop
1004	endif
1005	endif
1006
1007	! hybrid coordinate bp
1008	allocate(bp(GCM_layers+1),stat=ierr)
1009	if (ierr.ne.0) then
1010	write(,) "init2: failed to allocate bp!"
1011	stop
1012	endif
1013	ierr=NF_INQ_VARID(infid,"bp",tmpvarid)
1014	if (ierr.ne.NF_NOERR) then
1015	write(,) "Ooops. Failed to get bp ID. OK."
1016	apbp=.false.
1017	else
1018	ierr=NF_GET_VAR_REAL(infid,tmpvarid,bp)
1019	apbp=.true.
1020	if (ierr.ne.NF_NOERR) then
1021	write(,) "Error: Failed reading bp"
1022	stop
1023	endif
1024	endif
1025
1026	! sigma levels (if any)
1027	if (.not.hybrid) then
1028	allocate(sigma(GCM_layers),stat=ierr)
1029	if (ierr.ne.0) then
1030	write(,) "init2: failed to allocate sigma"
1031	stop
1032	endif
1033	ierr=NF_INQ_VARID(infid,"sigma",tmpvarid)
1034	if (ierr.ne.NF_NOERR) then
1035	write(,) "Ooops. Failed to get sigma ID. OK."
1036	sig=.false.
1037	else
1038	ierr=NF_GET_VAR_REAL(infid,tmpvarid,sigma)
1039	sig=.true.
1040	if (ierr.ne.NF_NOERR) then
1041	write(,) "init2 Error: Failed reading sigma"
1042	stop
1043	endif
1044	endif
1045	endif ! of if (.not.hybrid)
1046
1047	! mesh area
1048	allocate(aire(lonlen,latlen),stat=ierr)
1049	if (ierr.ne.0) then
1050	write(,) "init2: failed to allocate aire!"
1051	stop
1052	endif
1053	ierr=NF_INQ_VARID(infid,"aire",tmpvarid)
1054	if (ierr.ne.NF_NOERR) then
1055	write(,)"init2 warning: Failed to get aire ID."
1056	area = .false.
1057	else
1058	ierr=NF_GET_VAR_REAL(infid,tmpvarid,aire)
1059	if (ierr.ne.NF_NOERR) then
1060	write(,) "init2 Error: Failed reading aire"
1061	stop
1062	endif
1063	area = .true.
1064	endif
1065
1066	! ground geopotential phisinit
1067	allocate(phisinit(lonlen,latlen),stat=ierr)
1068	if (ierr.ne.0) then
1069	write(,) "init2: failed to allocate phisinit!"
1070	stop
1071	endif
1072	ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
1073	if (ierr.ne.NF_NOERR) then
1074	write(,)"init2 warning: Failed to get phisinit ID."
1075	phis = .false.
1076	else
1077	ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
1078	if (ierr.ne.NF_NOERR) then
1079	write(,) "init2 Error: Failed reading phisinit"
1080	stop
1081	endif
1082	phis = .true.
1083	endif
1084
1085	!==============================================================================
1086	! 2. Write
1087	!==============================================================================
1088
1089	!==============================================================================
1090	! 2.2. Hybrid coordinates ap() , bp(), aps() and bps()
1091	!==============================================================================
1092	if(hybrid) then
1093	! define aps
1094	call def_var(nout,"aps","hybrid pressure at midlayers"," ",1,&
1095	(/layerdimout/),apsid,ierr)
1096	if (ierr.ne.NF_NOERR) then
1097	write(,) "init2 Error: Failed to def_var aps"
1098	stop
1099	endif
1100
1101	! write aps
1102	ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
1103	if (ierr.ne.NF_NOERR) then
1104	write(,) "init2 Error: Failed to write aps"
1105	stop
1106	endif
1107
1108	! define bps
1109	call def_var(nout,"bps","hybrid sigma at midlayers"," ",1,&
1110	(/layerdimout/),bpsid,ierr)
1111	if (ierr.ne.NF_NOERR) then
1112	write(,) "init2 Error: Failed to def_var bps"
1113	stop
1114	endif
1115
1116	! write bps
1117	ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
1118	if (ierr.ne.NF_NOERR) then
1119	write(,) "init2 Error: Failed to write bps"
1120	stop
1121	endif
1122
1123	if (apbp) then
1124	! define ap
1125	call def_var(nout,"ap","hybrid sigma at interlayers"," ",1,&
1126	(/interlayerdimout/),apid,ierr)
1127	if (ierr.ne.NF_NOERR) then
1128	write(,) "Error: Failed to def_var ap"
1129	stop
1130	endif
1131
1132	! write ap
1133	ierr=NF_PUT_VAR_REAL(outfid,apid,ap)
1134	if (ierr.ne.NF_NOERR) then
1135	write(,) "Error: Failed to write ap"
1136	stop
1137	endif
1138
1139	! define bp
1140	call def_var(nout,"bp","hybrid sigma at interlayers"," ",1,&
1141	(/interlayerdimout/),bpid,ierr)
1142	if (ierr.ne.NF_NOERR) then
1143	write(,) "Error: Failed to def_var bp"
1144	stop
1145	endif
1146
1147	! write bp
1148	ierr=NF_PUT_VAR_REAL(outfid,bpid,bp)
1149	if (ierr.ne.NF_NOERR) then
1150	write(,) "Error: Failed to write bp"
1151	stop
1152	endif
1153	endif ! of if (apbp)
1154
1155	else if (sig) then
1156	! define sigma
1157	call def_var(nout,"sigma","sigma at midlayers"," ",1,&
1158	(/layerdimout/),sigmaid,ierr)
1159	if (ierr.ne.NF_NOERR) then
1160	write(,) "init2 Error: Failed to def_var sigma"
1161	stop
1162	endif
1163	! write sigma
1164	ierr=NF_PUT_VAR_REAL(outfid,sigmaid,sigma)
1165	if (ierr.ne.NF_NOERR) then
1166	write(,) "init2 Error: Failed to write sigma"
1167	stop
1168	endif
1169	endif ! of if (hybrid)
1170
1171	!==============================================================================
1172	! 2.2. aire() and phisinit()
1173	!==============================================================================
1174
1175	if (area) then
1176	! define aire
1177	call def_var(nout,"aire","Mesh area","m2",2,&
1178	(/londimout,latdimout/),tmpvarid,ierr)
1179	if (ierr.ne.NF_NOERR) then
1180	write(,) "init2 Error: Failed to def_var aire"
1181	stop
1182	endif
1183
1184	! write aire
1185	ierr=NF_PUT_VAR_REAL(outfid,tmpvarid,aire)
1186	if (ierr.ne.NF_NOERR) then
1187	write(,) "init2 Error: Failed to write aire"
1188	stop
1189	endif
1190	endif ! of if (area)
1191
1192	IF (phis) THEN
1193
1194	!define phisinit
1195	call def_var(nout,"phisinit","Ground level geopotential"," ",2,&
1196	(/londimout,latdimout/),phisinitid,ierr)
1197	if (ierr.ne.NF_NOERR) then
1198	write(,) "init2 Error: Failed to def_var phisinit"
1199	stop
1200	endif
1201
1202	! write phisinit
1203	ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
1204	if (ierr.ne.NF_NOERR) then
1205	write(,) "init2 Error: Failed to write phisinit"
1206	stop
1207	endif
1208
1209	ENDIF ! of IF (phis)
1210
1211
1212	! Cleanup
1213	if (allocated(aps)) deallocate(aps)
1214	if (allocated(bps)) deallocate(bps)
1215	if (allocated(ap)) deallocate(ap)
1216	if (allocated(bp)) deallocate(bp)
1217	if (allocated(sigma)) deallocate(sigma)
1218	if (allocated(phisinit)) deallocate(phisinit)
1219	if (allocated(aire)) deallocate(aire)
1220
1221	end subroutine init2
1222	!******************************************************************************
1223	subroutine def_var(nid,name,long_name,units,nbdim,dim,nvarid,ierr)
1224	!==============================================================================
1225	! Purpose: Write a variable (i.e: add a variable to a dataset)
1226	! called "name"; along with its attributes "long_name", "units"...
1227	! to a file (following the NetCDF format)
1228	!==============================================================================
1229	! Remarks:
1230	! The NetCDF file must be open
1231	!==============================================================================
1232
1233	implicit none
1234
1235	include "netcdf.inc" ! NetCDF definitions
1236
1237	!==============================================================================
1238	! Arguments:
1239	!==============================================================================
1240	integer, intent(in) :: nid
1241	! nid: [netcdf] file ID #
1242	character (len=*), intent(in) :: name
1243	! name(): [netcdf] variable's name
1244	character (len=*), intent(in) :: long_name
1245	! long_name(): [netcdf] variable's "long_name" attribute
1246	character (len=*), intent(in) :: units
1247	! unit(): [netcdf] variable's "units" attribute
1248	integer, intent(in) :: nbdim
1249	! nbdim: number of dimensions of the variable
1250	integer, dimension(nbdim), intent(in) :: dim
1251	! dim(nbdim): [netcdf] dimension(s) ID(s)
1252	integer, intent(out) :: nvarid
1253	! nvarid: [netcdf] ID # of the variable
1254	integer, intent(out) :: ierr
1255	! ierr: [netcdf] subroutines returned error code
1256
1257	! Switch to netcdf define mode
1258	ierr=NF_REDEF(nid)
1259
1260	! Insert the definition of the variable
1261	ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
1262
1263	! Write the attributes
1264	ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name",len_trim(adjustl(long_name)),adjustl(long_name))
1265	ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
1266
1267	! End netcdf define mode
1268	ierr=NF_ENDDEF(nid)
1269
1270	end subroutine def_var
1271	!******************************************************************************
1272	subroutine change_time_axis(nid,ierr,axis)
1273	!==============================================================================
1274	! Purpose:
1275	! Read "time" variable from a dataset, convert it from "sol" (martian
1276	! days) to "Ls" (solar longitude) and write the result back in the file
1277	!==============================================================================
1278	! Remarks:
1279	! The NetCDF file must be opened before this subroutine is called
1280	!==============================================================================
1281
1282	implicit none
1283
1284	include "netcdf.inc" ! NetCDF definitions
1285
1286	!==============================================================================
1287	! Arguments:
1288	!==============================================================================
1289	integer, intent(in) :: nid
1290	! nid: [netcdf] file ID
1291	integer, intent(out) :: ierr
1292	! ierr: [netcdf] return error code
1293	character (len=4),intent(in) :: axis
1294	! axis: "ls" or "sol"
1295
1296	!==============================================================================
1297	! Local variables:
1298	!==============================================================================
1299	integer :: nvarid
1300	! nvarid: ID of the "Time" variable
1301	integer :: timelen
1302	! timelen: size of the arrays
1303	integer :: timedim
1304	! timedim: ID of the "Time" dimension
1305	integer i
1306
1307	real, dimension(:), allocatable :: time,ls
1308	! time(): time, given in sols
1309	! ls(): time, given in Ls (solar longitude)
1310
1311	!==============================================================================
1312	! 1. Read
1313	!==============================================================================
1314
1315	ierr=NF_INQ_DIMID(nid,"Time",timedim)
1316	ierr=NF_INQ_VARID(nid,"Time",nvarid)
1317	if (ierr.NE.NF_NOERR) then
1318	write(,) 'ERROR in change_time_axis: Field <Time> not found'
1319	print*, NF_STRERROR(ierr)
1320	stop
1321	endif
1322
1323	ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
1324	allocate(time(timelen),ls(timelen))
1325	ierr = NF_GET_VAR_REAL(nid,nvarid,time)
1326	if (ierr.ne.NF_NOERR) then
1327	write(,) "ERROR in change_time_axis: Failed to load Time"
1328	stop
1329	endif
1330
1331	!==============================================================================
1332	! 2. Convert sols to Ls
1333	!==============================================================================
1334
1335	do i=1,timelen
1336	call sol2ls(time(i),ls(i))
1337	enddo
1338
1339	!==============================================================================
1340	! 2.1 Check if there are not jumps in Ls (rounding problems in sol2ls)
1341	!==============================================================================
1342
1343	do i=1,timelen-1
1344	if ((ls(i+1)-ls(i)) > 350) then
1345	write(,) "+ 360 deg Ls jump solved:", ls(i), ls(i+1), "at timestep", i
1346	ls(i+1) = ls(i+1) - 360
1347	write(,) " corrected to now be:", ls(i), ls(i+1)
1348	else if ((ls(i)-ls(i+1)) > 350) then
1349	write(,) "- 360 deg Ls jump solved:", ls(i), ls(i+1), "at timestep", i
1350	ls(i+1) = ls(i+1) + 360
1351	write(,) " corrected to now be:", ls(i), ls(i+1)
1352	endif
1353	enddo
1354
1355	!==============================================================================
1356	! 3. Write
1357	!==============================================================================
1358
1359	if (axis.eq."ls") then
1360	write(,) "Converting the time axis from sol to Ls..."
1361	else ! if (axis.eq."adls")
1362	write(,) "Adding Ls as a 1D time variable"
1363	call def_var(nout,"Ls","Ls (Solar Longitude)","degrees",1, (/timedimout/),nvarid,ierr)
1364	end if
1365
1366	ierr = NF_PUT_VAR_REAL(nid,nvarid,ls)
1367	if (ierr.ne.NF_NOERR) then
1368	write(,) "ERROR in change_time_axis: Failed to write Ls"
1369	stop
1370	endif
1371
1372	end subroutine change_time_axis
1373	!******************************************************************************
1374	subroutine sol2ls(sol,Ls)
1375	!==============================================================================
1376	! Purpose:
1377	! Convert a date/time, given in sol (martian day),
1378	! into solar longitude date/time, in Ls (in degrees),
1379	! where sol=0 is (by definition) the northern hemisphere
1380	! spring equinox (where Ls=0).
1381	!==============================================================================
1382	! Notes:
1383	! Even though "Ls" is cyclic, if "sol" is greater than N (martian) year,
1384	! "Ls" will be increased by N*360
1385	! Won't work as expected if sol is negative (then again,
1386	! why would that ever happen?)
1387	!==============================================================================
1388
1389	implicit none
1390
1391	!==============================================================================
1392	! Arguments:
1393	!==============================================================================
1394	real,intent(in) :: sol
1395	real,intent(out) :: Ls
1396
1397	!==============================================================================
1398	! Local variables:
1399	!==============================================================================
1400	real year_day,peri_day,timeperi,e_elips,twopi,degrad
1401	data year_day /669./ ! # of sols in a martian year
1402	data peri_day /485.0/
1403	data timeperi /1.9082314/
1404	data e_elips /0.093358/
1405	data twopi /6.2831853/ ! 2.*pi
1406	data degrad /57.2957795/ ! pi/180
1407
1408	real zanom,xref,zx0,zdx,zteta,zz
1409
1410	integer count_years
1411	integer iter
1412
1413	!==============================================================================
1414	! 1. Compute Ls
1415	!==============================================================================
1416
1417	zz=(sol-peri_day)/year_day
1418	zanom=twopi*(zz-nint(zz))
1419	xref=abs(zanom)
1420
1421	! The equation zx0 - e * sin (zx0) = xref, solved by Newton
1422	zx0=xref+e_elips*sin(xref)
1423	do iter=1,20 ! typically, 2 or 3 iterations are enough
1424	zdx=-(zx0-e_elipssin(zx0)-xref)/(1.-e_elipscos(zx0))
1425	zx0=zx0+zdx
1426	if(abs(zdx).le.(1.e-7)) then
1427	! write(,)'iter:',iter,' \|zdx\|:',abs(zdx)
1428	exit
1429	endif
1430	enddo
1431
1432	if(zanom.lt.0.) zx0=-zx0
1433
1434	zteta=2.atan(sqrt((1.+e_elips)/(1.-e_elips))tan(zx0/2.))
1435	Ls=zteta-timeperi
1436
1437	if(Ls.lt.0.) then
1438	Ls=Ls+twopi
1439	else
1440	if(Ls.gt.twopi) then
1441	Ls=Ls-twopi
1442	endif
1443	endif
1444
1445	Ls=degrad*Ls
1446	! Ls is now in degrees
1447
1448	!==============================================================================
1449	! 1. Account for (eventual) years included in input date/time sol
1450	!==============================================================================
1451
1452	count_years=0 ! initialize
1453	zz=sol ! use "zz" to store (and work on) the value of sol
1454	do while (zz.ge.year_day)
1455	count_years=count_years+1
1456	zz=zz-year_day
1457	enddo
1458
1459	! Add 360 degrees to Ls for every year
1460	if (count_years.ne.0) then
1461	Ls=Ls+360.*count_years
1462	endif
1463
1464
1465	end subroutine sol2ls
1466	!******************************************************************************
1467	subroutine missing_value(nout,nvarid,validr,miss,valid_range,missing)
1468	!==============================================================================
1469	! Purpose:
1470	! Write "valid_range" and "missing_value" attributes (of a given
1471	! variable) to a netcdf file
1472	!==============================================================================
1473	! Remarks:
1474	! NetCDF file must be open
1475	! Variable (nvarid) ID must be set
1476	!==============================================================================
1477
1478	implicit none
1479
1480	include "netcdf.inc" ! NetCDF definitions
1481
1482	!==============================================================================
1483	! Arguments:
1484	!==============================================================================
1485	integer, intent(in) :: nout
1486	! nout: [netcdf] file ID #
1487	integer, intent(in) :: nvarid
1488	! varid: [netcdf] variable ID #
1489	integer, intent(in) :: validr
1490	! validr : [netcdf] routines return code for "valid_range" attribute
1491	integer, intent(in) :: miss
1492	! miss : [netcdf] routines return code for "missing_value" attribute
1493	real, dimension(2), intent(in) :: valid_range
1494	! valid_range(2): [netcdf] "valid_range" attribute (min and max)
1495	real, intent(in) :: missing
1496	! missing: [netcdf] "missing_value" attribute
1497
1498	!==============================================================================
1499	! Local variables:
1500	!==============================================================================
1501	integer :: ierr
1502	! ierr: [netcdf] subroutine returned error code
1503	! INTEGER nout,nvarid,ierr
1504	! REAL missing
1505	! REAL valid_range(2)
1506
1507	! Switch to netcdf dataset define mode
1508	ierr = NF_REDEF (nout)
1509	if (ierr.ne.NF_NOERR) then
1510	print*,'missing_value: '
1511	print*, NF_STRERROR(ierr)
1512	endif
1513
1514	! Write valid_range() attribute
1515	if (validr.eq.NF_NOERR) then
1516	ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
1517
1518	if (ierr.ne.NF_NOERR) then
1519	print*,'missing_value: valid_range attribution failed'
1520	print*, NF_STRERROR(ierr)
1521	!write(,) 'NF_NOERR', NF_NOERR
1522	!CALL abort
1523	stop
1524	endif
1525	endif
1526
1527	! Write "missing_value" attribute
1528	if (miss.eq.NF_NOERR) then
1529	ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,[missing])
1530
1531	if (ierr.NE.NF_NOERR) then
1532	print*, 'missing_value: missing value attribution failed'
1533	print*, NF_STRERROR(ierr)
1534	! WRITE(,) 'NF_NOERR', NF_NOERR
1535	! CALL abort
1536	stop
1537	endif
1538	endif
1539
1540	! End netcdf dataset define mode
1541	ierr = NF_ENDDEF(nout)
1542
1543	end subroutine missing_value
1544	!******************************************************************************
1545
1546	end program concatnc

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

Download in other formats:

Original Format