Context Navigation

solzenangle.F90 @ 2567

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