Context Navigation

solzenangle.F90 @ 2289

Last change on this file since 2289 was 2289, checked in by abierjon, 5 years ago
Mars GCM: Add the program solzenangle.F90 to the LMDZ.MARS/util/ directory, as well as a solzenangle.def This program actually replaces and seals the dark fate of the former program terminator.F90 It can be used to interpolate GCM files at one given solar zenith angle (on the morning-side or the evening-side) all around the globe (the terminator corresponding to sza = 90deg) By Alwont Beebak
File size: 51.0 KB

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

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

Download in other formats:

Original Format