Context Navigation

simu_MCS.F90 @ 2443

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

Mars GCM:
Changes in aeroptical.F90 :

add the flag nf90_64bit_offset when creating the output file. Without it, the program couldn't deal with large files (concat) with multiple aerosols
modify the code structure to make it a bit less memory-intensive
add topdust in computed aerosol opacities
add the possibility to compute absorption opacities, instead of extinction : the user must specify 'ext' or 'abs' in aeroptical.def
update the description of aeroptical in util/README

+ adapt simu_MCS to take into account aeroptical output variables when choosing

the binning reference

+ add outputs for topdust in physiq_mod

File size: 89.2 KB

Line
1	program simu_MCS
2	! program for a MCS observer simulator of GCM data
3	! author : Antoine Bierjon, 2019-2020
4	! contact : antoine.bierjon@lmd.jussieu.fr
5	!
6	!===================================================================================================
7	! PREFACE
8	!===================================================================================================
9	! This program loads on one hand, a GCM file, which can be :
10	! - a GCM diagfi.nc or concat.nc that covers the observation period,
11	! - a GCM stats file (stats.nc) which gets extended to cover the observation period
12	!
13	! On the other hand, it loads a MRO/MCS data file (binned by Luca Montabone).
14	! This MCS file serves then as a reference for the interpolation and the binning of the
15	! GCM variables contained in the GCM file and given as inputs in the simu_MCS.def
16	!
17	! Since the MCS data is binned by intervals of Ls, the GCM simulation's data is interpolated at
18	! the MCS spatial coordinates, at every GCM sol value contained in each interval of Ls, before
19	! being averaged to create output bins for each interval. The binning in sols also takes into
20	! account the variability of Local Time in the bin and tries to represent it (with a number of LT
21	! values in each sol bin equal to the number of values in each associated MCS bin : numbintemp/dust/wice,
22	! and centered around the MCS bin LT average : timeave).
23	!
24	! There are also specific GCM variables that the program looks for in order to make them
25	! comparable with their equivalent in MCS files. These variables are :
26	! GCM MCS
27	! # temp ---> dtemp,ntemp
28	! # dso(/dsodust/qdust)+rho ---> ddust,ndust (dust opacity/km)
29	! # h2o_ice+rho ---> dwice,nwice (water ice opacity/km)
30	!
31	! Eventually, the program outputs a netcdf file, filled with the GCM binned data for dayside
32	! and nightside and edited with the same format than the MCS file (and with the same
33	! dimensions' declaration order).
34	!
35	! Minimal requirements :
36	! - the MCS file must contain the variables :
37	! dtimeave,dtimemax,dtimemin,
38	! dtemp,dnumbintemp,
39	! ddust,dnumbindust,
40	! dwice,dnumbinwice,
41	! ntimeave,ntimemax,ntimemin,
42	! ntemp,nnumbintemp,
43	! ndust,nnumbindust,
44	! nwice,nnumbinwice
45	! - the GCM file must :
46	! # have the same altitude type as the MCS file ;
47	! # cover completely the observation period
48	!
49	! See also NOTA BENE in section 2.2
50	!
51	!
52	!
53	! Algorithm :
54	! 0. Variable declarations
55	! 1. OPENING OF THE FILES AND INITIALIZATION OF THE DIMENSIONS
56	! 1.1 MCS data file : obsfile
57	! 1.1.1 Open the Observer data file
58	! 1.1.2 Get dimensions lon,lat,alt,time from the observation file
59	! 1.2. GCM simulation file : gcmfile
60	! 1.2.1 Open the GCM simulation file
61	! 1.2.2 Get dimensions lon,lat,alt,time from the GCM file
62	! 1.3 Create the output file & initialize the coordinates
63	! 2. VARIABLES MANAGEMENT
64	! 2.1 List of the GCM variables to be interpolated
65	! 2.1.1 Read the GCM variables
66	! 2.1.2 Handle dust and wice opacities (first set-ups)
67	! 2.2 Definition of LT variables from obsfile to outfile
68	![day/night loop begins...
69	! 2.2.1 Average of Local Time in the OBS bins
70	! 2.2.2 Maximum of Local Time in the OBS bins
71	! 2.2.3 Minimum of Local Time in the OBS bins
72	! 2.3 Definition of numbin variables from obsfile to outfile
73	! 2.3.1 Number of values in the OBS "temp" bins
74	! 2.3.2 Number of values in the OBS "dust" bins
75	! 2.3.3 Number of values in the OBS "wice" bins
76	! 2.4 Opening of the GCM variables
77	! [var loop begins...
78	! 2.4.1 Generic reading of the variable
79	! 2.4.2 Handle dust and wice opacities (second part)
80	! 2.5 Opening of the associated MCS variables
81	! 2.5.1 MCS reference variable (for the missing values)
82	! 2.5.2 Number of values in the OBS bin (for the sol binning)
83	! 2.6 Definition of GCM variables in outfile
84	! 3. EXTRACTION OF THE VARIABLE
85	! [coordinates loop begins...
86	! 3.1 Do some checks and preparations before the extraction
87	! 3.2 Compute GCM sol date corresponding to Observer Ls (via m_(min/max)sol)
88	! and LT (via OBSLT(min/max))
89	! 3.3 Do the interpolation and binning for the given location
90	! ..coordinates loop ends]
91	! 3.4 Write the data in the netcdf output file
92	! ..var loop ends]
93	! 4. END OF THE DAY/NIGHT LOOP
94	!..day/night loop ends]
95	! 5. CLOSE THE FILES
96	!
97	! Subroutines
98	! extraction
99	! inidim
100	! ls2sol
101	! gen_sol_list
102	! status_check
103	! LTmod
104	!===================================================================================================
105
106
107	use netcdf
108
109
110	!===================================================================================
111	! 0. Variable declarations
112	!===================================================================================
113
114	implicit none ! for no implicitly typed variables
115
116	!------------------------
117	! Files:
118	character(len=256) :: gcmfile ! GCM simulation file
119	logical :: is_stats = .false. ! to check if the GCM file is a stats.nc or a diagfi.nc file
120	character(len=256) :: obsfile ! observation file = MCS/Observer data file
121	character(len=256) :: outfile ! output file
122
123	!------------------------
124	! NetCDF stuff
125	integer :: status ! NetCDF routines return code
126	character (len=256) :: error_text ! to store the error text
127	integer :: gcmfid ! NetCDF gcm file ID
128	integer :: obsfid ! NetCDF observation file ID
129	integer :: outfid ! NetCDF output file ID
130	integer :: GCMvarid, OBSvarid, LT_id, numbin_id, outvarid ! to store the ID of a variable
131	integer :: lat_dimid_obs,lon_dimid_obs,alt_dimid_obs,time_dimid_obs ! dimensions' ID in OBS and output files
132	integer :: lat_dimid_gcm,lon_dimid_gcm,alt_dimid_gcm,time_dimid_gcm ! dimensions' ID in GCM file
133	integer :: GCMvarshape(4), OBSvarshape(4), LTshape(3),numbinshape(4) ! to store a variable's coordinates order
134
135	!------------------------
136	! Dimensions
137	real,dimension(:),allocatable :: GCMlon, OBSlon ! longitude in the GCM & the Observer files
138	integer GCMlonlen, OBSlonlen ! # of grid points along GCMlon & OBSlon
139	real,dimension(:),allocatable :: GCMlat, OBSlat ! latitude in the GCM & the Observer files
140	integer GCMlatlen, OBSlatlen ! # of grid points along GCMlat & OBSlat
141	real,dimension(:),allocatable :: GCMalt, OBSalt ! altitude/pressure in the GCM & the Observer files
142	integer GCMaltlen, OBSaltlen ! # of grid point along GCMalt & OBSalt
143	character(len=1) :: GCMalttype, OBSalttype ! altitude coord. type:'z' (altitude, m) 'p' (pressure, Pa)
144	real,dimension(:),allocatable :: GCMtime, OBSLs ! time in the GCM diagfi (sols) & the Observer files (Ls)
145	real,dimension(:),allocatable :: GCMstatstime ! time in the GCM stats file (LT at lon 0°)
146	integer :: GCMtimelen, GCMstatstimelen, OBSLslen ! # of points along GCMtime, GCMstatstime, OBSLs
147	real :: starttimeoffset=0. ! offset (in sols) wrt Ls=0 of sol 0 in GCM file
148
149	!------------------------
150	! Variables
151	character(len=10) :: dayornight ! are we in the "dayside" or "nightside" loop?
152	character(len=64),dimension(:),allocatable :: gcm_vars ! list of GCM variables to interpolate
153	character(len=10),dimension(15) :: notprocessed ! names of the (15) variables that won't be processed
154	integer :: nbvarfile,Nnotprocessed ! nbs of variables to deal with the non-processed ones
155	integer :: nbvar ! nb of variables that will be processed
156	logical :: var_ok ! is this variable to be processed?
157	logical :: dustok1,dustok2,dustok3,wiceok = .false. ! is it possible to compute opacities and how?
158	character(len=64) :: GCMvarname,OBSvarname,outvarname ! name of the variables
159	integer :: nbdim ! nb of dimensions of a variable
160	real,dimension(:,:,:,:),allocatable :: GCM_var,OBS_var ! the 4D variable extracted from GCM & OBS files
161	real,dimension(:,:,:,:),allocatable :: rho ! atmospheric density for opacities' computation
162	character(len=64) :: long_name,units,comment ! netcdf attributes of the variable
163
164	character(len=64) :: OBSLTave_name,OBSLTmax_name,OBSLTmin_name
165	! names of the average, max and min of LT
166	real,dimension(:,:,:),allocatable :: OBSLT,OBSLTmax,OBSLTmin
167	! 3D variables extracted from obsfile (ave, max and min of LT in a bin)
168
169	character (len=64) :: numbin_name ! name of the nb of values in an OBS bin
170	real,dimension(:,:,:,:),allocatable :: numbin ! nb of values in an OBS temp/dust/wice bin
171
172	real :: GCMmiss_val, OBSmiss_val, LTmiss_val ! value to denote non-existant data
173	real :: extr_value ! result of the extraction subroutine
174	real, dimension(:,:,:,:), allocatable :: outvar ! outvar(,,,): 4D array to store the output variable's data
175
176	!------------------------
177	! Time binning management
178	real :: OBSdeltaLs ! difference of Ls between each observation bin
179	real :: sol, maxsol, minsol ! sol date corresponding to Observed Ls and LT
180	integer :: m_maxsol, m_minsol ! indexes of the maximum and minimum GCM sol taken in a bin for interpolation
181
182	external LTmod ! declaration of the function LTmod
183	real :: LTmod ! declaration of the type of the function LTmod
184	integer :: LTcount ! nb of LT samples on which the interpolation is performed,
185	! for every LT interval (= numbin[lon,lat,alt,Ls])
186
187	integer :: solcount ! number of GCM sol integer values in one Ls interval
188	real,dimension(:),allocatable :: int_sol_list ! list of the integer values of GCM sol
189	real,dimension(:),allocatable :: sol_list ! list of the sol values used for the interpolation
190	integer :: solerrcount ! nb of GCM missing values during interpolation, removed for the binning
191	integer :: errcount = 0 ! total number of GCM missing values
192	real :: solbinned_value ! extracted value averaged on sol samples, which is finally put in the output bin
193
194	!------------------------
195	! Extraction & loop indices
196	real :: lon_val, lat_val, alt_val, Ls_val, LT_val ! where and when the output file is written at
197
198	integer :: i,j,k,l ! loop iteration indices
199	integer :: m ! sol binning loops iteration index
200	integer :: v,vnot ! variable loops indices
201
202
203	!===================================================================================
204	! 1. OPENING OF THE FILES AND INITIALIZATION OF THE DIMENSIONS
205	!===================================================================================
206	write(,) "Welcome in the MRO/MCS Observer Simulator program !"
207	!===============================================================================
208	! 1.1 MCS data file : obsfile
209	!===============================================================================
210	!================================================================
211	! 1.1.1 Open the Observer data file
212	!================================================================
213	! Ask the user to give a netcdf observation file
214	WRITE(,) "-> Enter observation file name :"
215	READ(,) obsfile
216
217	status=NF90_OPEN(obsfile,nf90_nowrite,obsfid) ! nowrite mode=the program can only read the opened file
218	error_text="Error: could not open file "//trim(obsfile)
219	call status_check(status,error_text)
220
221	!================================================================
222	! 1.1.2 Get dimensions lon,lat,alt,time from the observation file
223	!================================================================
224	! OBS Latitude
225	!--------------
226	status=nf90_inq_dimid(obsfid,"latitude",lat_dimid_obs)
227	error_text="Failed to find Observer latitude dimension"
228	call status_check(status,error_text)
229
230	status=nf90_inquire_dimension(obsfid,lat_dimid_obs,len=OBSlatlen)
231	error_text="Failed to find Observer latitude length"
232	call status_check(status,error_text)
233	allocate(OBSlat(OBSlatlen))
234
235	status=nf90_inq_varid(obsfid,"latitude",OBSvarid)
236	error_text="Failed to find Observer latitude ID"
237	call status_check(status,error_text)
238
239	! Read OBSlat
240	status=NF90_GET_VAR(obsfid,OBSvarid,OBSlat)
241	error_text="Failed to load OBSlat"
242	call status_check(status,error_text)
243
244
245	! OBS Longitude
246	!--------------
247	status=nf90_inq_dimid(obsfid,"longitude",lon_dimid_obs)
248	error_text="Failed to find Observer longitude dimension"
249	call status_check(status,error_text)
250
251	status=nf90_inquire_dimension(obsfid,lon_dimid_obs,len=OBSlonlen)
252	error_text="Failed to find Observer longitude length"
253	call status_check(status,error_text)
254	allocate(OBSlon(OBSlonlen))
255
256	status=nf90_inq_varid(obsfid,"longitude",OBSvarid)
257	error_text="Failed to find Observer longitude ID"
258	call status_check(status,error_text)
259
260	! Read OBSlon
261	status=NF90_GET_VAR(obsfid,OBSvarid,OBSlon)
262	error_text="Failed to load OBSlon"
263	call status_check(status,error_text)
264
265
266	! OBS Time (Ls)
267	!--------------
268	status=nf90_inq_dimid(obsfid,"time",time_dimid_obs)
269	error_text="Failed to find Observer time (Ls) dimension"
270	call status_check(status,error_text)
271
272	status=nf90_inquire_dimension(obsfid,time_dimid_obs,len=OBSLslen)
273	error_text="Failed to find Observer time (Ls) length"
274	call status_check(status,error_text)
275	allocate(OBSLs(OBSLslen))
276
277	status=nf90_inq_varid(obsfid,"time",OBSvarid)
278	error_text="Failed to find Observer time (Ls) ID"
279	call status_check(status,error_text)
280
281	! Read OBSLs
282	status=NF90_GET_VAR(obsfid,OBSvarid,OBSLs)
283	error_text="Failed to load OBSLs"
284	call status_check(status,error_text)
285
286	! Get the observation timestep between bins
287	OBSdeltaLs=OBSLs(2)-OBSLs(1)
288
289
290	! OBS Altitude
291	!--------------
292	status=nf90_inq_dimid(obsfid,"altitude",alt_dimid_obs)
293	error_text="Failed to find Observer altitude dimension"
294	call status_check(status,error_text)
295
296	status=nf90_inquire_dimension(obsfid,alt_dimid_obs,len=OBSaltlen)
297	error_text="Failed to find Observer altitude length"
298	call status_check(status,error_text)
299	allocate(OBSalt(OBSaltlen))
300
301	status=nf90_inq_varid(obsfid,"altitude",OBSvarid)
302	error_text="Failed to find Observer altitude ID"
303	call status_check(status,error_text)
304
305	! Read OBSalt
306	status=NF90_GET_VAR(obsfid,OBSvarid,OBSalt)
307	error_text="Failed to load OBSalt"
308	call status_check(status,error_text)
309
310	! Check altitude attribute "units" to find out altitude type and compare with the GCM file
311	status=nf90_get_att(obsfid,OBSvarid,"units",units)
312	error_text="Failed to load Observer altitude units attribute"
313	call status_check(status,error_text)
314	! an unknown and invisible character is placed just after the unit's
315	! characters in the Observer file so we only take the first characters
316	! corresponding to the sought unit
317	if (trim(units(1:2)).eq."Pa") then
318	units="Pa"
319	OBSalttype='p' ! pressure coordinate
320	else if (trim(units(1:2)).eq."m") then
321	units="m"
322	OBSalttype='z' ! altitude coordinate
323	else
324	write(,)" I do not understand this unit ",trim(units)," for Observer altitude!"
325	stop
326	endif
327
328	!===============================================================================
329	! 1.2. GCM simulation file : gcmfile
330	!===============================================================================
331	!================================================================
332	! 1.2.1 Open the GCM simulation file
333	!================================================================
334	! Ask the user to give a netcdf input file
335	write(,)"";WRITE(,) "-> Enter input file name (GCM simulation) :"
336	READ(,) gcmfile
337
338	! Open GCM file
339	status=NF90_OPEN(gcmfile,nf90_nowrite,gcmfid)
340	! nowrite mode=the program can only read the opened file
341	error_text="Failed to open datafile "//trim(gcmfile)
342	call status_check(status,error_text)
343
344	!================================================================
345	! 1.2.2 Get dimensions lon,lat,alt,time from the GCM file
346	!================================================================
347	! GCM Latitude
348	!--------------
349	status=nf90_inq_dimid(gcmfid,"latitude",lat_dimid_gcm)
350	error_text="Failed to find GCM latitude dimension"
351	call status_check(status,error_text)
352
353	status=nf90_inquire_dimension(gcmfid,lat_dimid_gcm,len=GCMlatlen)
354	error_text="Failed to find GCM latitude length"
355	call status_check(status,error_text)
356	allocate(GCMlat(GCMlatlen))
357
358	status=nf90_inq_varid(gcmfid,"latitude",GCMvarid)
359	error_text="Failed to find GCM latitude ID"
360	call status_check(status,error_text)
361
362	! Read GCMlat
363	status=NF90_GET_VAR(gcmfid,GCMvarid,GCMlat)
364	error_text="Failed to load GCMlat"
365	call status_check(status,error_text)
366
367
368	! GCM Longitude
369	!--------------
370	status=nf90_inq_dimid(gcmfid,"longitude",lon_dimid_gcm)
371	error_text="Failed to find GCM longitude dimension"
372	call status_check(status,error_text)
373
374	status=nf90_inquire_dimension(gcmfid,lon_dimid_gcm,len=GCMlonlen)
375	error_text="Failed to find GCM longitude length"
376	call status_check(status,error_text)
377	allocate(GCMlon(GCMlonlen))
378
379	status=nf90_inq_varid(gcmfid,"longitude",GCMvarid)
380	error_text="Failed to find GCM longitude ID"
381	call status_check(status,error_text)
382
383	! Read GCMlon
384	status=NF90_GET_VAR(gcmfid,GCMvarid,GCMlon)
385	error_text="Failed to load GCMlon"
386	call status_check(status,error_text)
387
388
389	! GCM Time
390	!--------------
391	status=nf90_inq_dimid(gcmfid,"Time",time_dimid_gcm)
392	error_text="Failed to find GCM time dimension"
393	call status_check(status,error_text)
394
395	status=nf90_inquire_dimension(gcmfid,time_dimid_gcm,len=GCMtimelen)
396	error_text="Failed to find GCM time length"
397	call status_check(status,error_text)
398	allocate(GCMtime(GCMtimelen))
399
400	status=nf90_inq_varid(gcmfid,"Time",GCMvarid)
401	error_text="Failed to find GCM time ID"
402	call status_check(status,error_text)
403
404	status=NF90_GET_VAR(gcmfid,GCMvarid,GCMtime)
405	error_text="Failed to load GCMtime"
406	call status_check(status,error_text)
407
408	! is_stats ?
409	IF ((GCMtimelen.eq.12).and.(GCMtime(1).eq.2.).and.(GCMtime(GCMtimelen).eq.24.)) then
410	! if GCM file is a stats, time is in LT at longitude 0° and not in sols at longitude 0°
411	write(,)"The GCM file is recognized as a stats file."
412	is_stats = .true.
413	deallocate(GCMtime)
414	GCMstatstimelen = GCMtimelen
415	allocate(GCMstatstime(GCMstatstimelen))
416	status=NF90_GET_VAR(gcmfid,GCMvarid,GCMstatstime)
417	error_text="Failed to load GCMstatstime (LT at lon 0)"
418	call status_check(status,error_text)
419	ELSE
420	write(,)"The GCM file is recognized as a diagfi/concatnc file."
421	ENDIF
422
423	! Simulation time offset management
424	WRITE(,) "Beginning date of the simulation file?"
425	WRITE(,) "(i.e. number of sols since Ls=0 at the Time=0.0 in the GCM file)"
426	READ(,) starttimeoffset
427	if (.not.is_stats) then
428	! Add the offset to GCMtime(:) if the file is not a stats file
429	GCMtime(:)=GCMtime(:)+starttimeoffset
430	endif
431
432	! Check of temporal coherence between gcmfile & obsfile
433	call ls2sol(OBSLs(OBSLslen),maxsol) ! maximum date considered
434	call ls2sol(OBSLs(1),minsol) ! minimum date considered
435
436	IF (.not.is_stats) then ! if it is a diagfi, we check the time coherence between the 2 files
437	if ((maxsol.gt.maxval(GCMtime)).or.(minsol.lt.minval(GCMtime))) then
438	write(,)"Error : obsfile temporal bounds exceed the GCM simulation bounds."
439	write(,)"Please use a GCM file whose time interval contains the observation period."
440	stop
441	else
442	write(,)"Both files are temporally coherent. The program continues..."
443	endif
444
445	ELSE ! if it is a stats, we create the array GCMtime array (in sols) covering the observation period
446	! and filled with the mean GCM day stored in stats.nc
447
448	GCMtimelen = ((ceiling(maxsol)-floor(minsol)+1)+2) ! we add 2 days in the beginning and the end
449	! to be sure we cover the observation period
450	allocate(GCMtime(GCMstatstimelen * GCMtimelen))
451	do l=1,GCMtimelen
452	do m=1,GCMstatstimelen
453	GCMtime(m+(l-1)*GCMstatstimelen) = (floor(minsol)-1) + (l-1) + GCMstatstime(m)/24.
454	enddo
455	enddo
456	GCMtimelen = GCMstatstimelen * GCMtimelen
457	write(,)"GCMtime has been created from the stats.nc time and the observation period. The program continues..."
458	ENDIF
459
460
461	! GCM Altitude
462	!--------------
463	status=nf90_inq_dimid(gcmfid,"altitude",alt_dimid_gcm)
464	error_text="Failed to find GCM altitude dimension"
465	call status_check(status,error_text)
466
467	status=nf90_inquire_dimension(gcmfid,alt_dimid_gcm,len=GCMaltlen)
468	error_text="Failed to find GCM altitude length"
469	call status_check(status,error_text)
470	allocate(GCMalt(GCMaltlen))
471
472	status=nf90_inq_varid(gcmfid,"altitude",GCMvarid)
473	error_text="Failed to find GCM altitude ID"
474	call status_check(status,error_text)
475
476	! Read GCMalt
477	status=NF90_GET_VAR(gcmfid,GCMvarid,GCMalt)
478	error_text="Failed to load GCMalt"
479	call status_check(status,error_text)
480
481	! Check altitude attribute "units" to find out altitude type
482	status=nf90_get_att(gcmfid,GCMvarid,"units",units)
483	error_text="Failed to load GCM altitude units attribute"
484	call status_check(status,error_text)
485	if (trim(units).eq."Pa") then
486	GCMalttype='p' ! pressure coordinate
487	else if (trim(units).eq."m") then
488	GCMalttype='z' ! altitude coordinate
489	else
490	write(,)"I do not understand this unit ",trim(units)," for GCM altitude!"
491	if (OBSalttype.eq.'p') then
492	write(,)"Please use zrecast to put the altitude in the same type as the MCS file (pressure in Pa)"
493	else if (OBSalttype.eq.'z') then
494	write(,)"Please use zrecast to put the altitude in the same type as the MCS file (altitude in m)"
495	endif
496	stop
497	endif
498	IF(OBSalttype.ne.GCMalttype) then
499	write(,)"Observer altitude type (", OBSalttype,") and ", &
500	"GCM altitude type (",GCMalttype,") don't match!"
501	stop
502	ENDIF
503
504	!===============================================================================
505	! 1.3 Create the output file & initialize the coordinates
506	!===============================================================================
507	! Name of the outfile
508	IF (.not.is_stats) then
509	outfile=obsfile(1:index(obsfile, ".nc")-1)//"_GCMdiagfi.nc"
510	ELSE
511	outfile=obsfile(1:index(obsfile, ".nc")-1)//"_GCMstats.nc"
512	ENDIF
513
514	! Creation of the outfile
515	status=NF90_CREATE(outfile,nf90_clobber,outfid)!NB: clobber mode=overwrite any dataset with the same file name !
516	error_text="Error: could not create file "//trim(outfile)
517	call status_check(status,error_text)
518	write(,)"";WRITE(,)"-> Output file is: ",trim(outfile)
519
520	! Creation of the dimensions
521	call inidim(outfid,OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen,OBSlon,OBSlat,OBSalt,OBSLs,OBSalttype,&
522	lon_dimid_obs,lat_dimid_obs,alt_dimid_obs,time_dimid_obs)
523
524	!write(,)"Dimensions ID in the outfile are :"
525	!write(,)"lon_dimid=",lon_dimid_obs
526	!write(,)"lat_dimid=",lat_dimid_obs
527	!write(,)"alt_dimid=",alt_dimid_obs
528	!write(,)"time_dimid=",time_dimid_obs
529
530
531
532	!===================================================================================
533	! 2. VARIABLES MANAGEMENT
534	!===================================================================================
535	!===============================================================================
536	! 2.1 List of the GCM variables to be interpolated
537	!===============================================================================
538	!================================================================
539	! 2.1.1 Read the GCM variables
540	!================================================================
541	! Get nbvarfile (total number of variables in the GCM file)
542	status=NF90_INQUIRE(gcmfid,nVariables=nbvarfile)
543	error_text="Error : Pb with nf90_inquire(gcmfid,nVariables=nbvarfile)"
544	call status_check(status,error_text)
545
546	! List of variables that should not be processed
547	notprocessed(1)='Time'
548	notprocessed(2)='controle'
549	notprocessed(3)='rlonu'
550	notprocessed(4)='latitude'
551	notprocessed(5)='longitude'
552	notprocessed(6)='altitude'
553	notprocessed(7)='rlatv'
554	notprocessed(8)='aps'
555	notprocessed(9)='bps'
556	notprocessed(10)='ap'
557	notprocessed(11)='bp'
558	notprocessed(12)='cu'
559	notprocessed(13)='cv'
560	notprocessed(14)='aire'
561	notprocessed(15)='phisinit'
562
563
564	! List of variables in the GCM file
565	write(,)""
566	write(,)"List of variables in the GCM file :"
567	Nnotprocessed=0
568	do v=1,nbvarfile
569	status=NF90_INQUIRE_VARIABLE(gcmfid,v,name=GCMvarname)
570	! GCMvarname now contains the "name" of variable of ID # v
571	var_ok=.true.
572	do vnot=1,15
573	if (GCMvarname.eq.notprocessed(vnot)) then
574	var_ok=.false.
575	Nnotprocessed=Nnotprocessed+1
576	endif
577	enddo
578	if (var_ok) write(,) trim(GCMvarname)
579
580	! Detect if we can compute dust and wice opacities
581	if (trim(GCMvarname).eq."dso") then
582	dustok1 = .true.
583	else if (trim(GCMvarname).eq."dsodust") then
584	dustok2 = .true.
585	else if (trim(GCMvarname).eq."dustq") then
586	dustok3 = .true.
587	else if (trim(GCMvarname).eq."h2o_ice") then
588	wiceok = .true.
589	endif
590	enddo
591
592	! Nnotprocessed: # of variables that won't be processed
593	! nbvarfile: total # of variables in file
594	! +2: the dust and wice opacities
595	allocate(gcm_vars(nbvarfile-Nnotprocessed+2),stat=status)
596	if (status.ne.0) then
597	write(,) "Error: failed allocation of gcm_vars(nbvarfile-Nnotprocessed+2)"
598	write(,) " nbvarfile=",nbvarfile
599	write(,) " Nnotprocessed=",Nnotprocessed
600	stop
601	endif
602
603	! List of variables to process
604	write(,)
605	write(,) "Which variables do you want to redistribute ?"
606	write(,) "list of <variables> (separated by <Return>s)"
607	write(,) "(an empty line , i.e: just <Return>, implies end of list)"
608	write(,) "NB: this program handles only 4D netcdf variables for now"
609	nbvar=0
610	read(*,'(a50)') GCMvarname
611	do while ((GCMvarname.ne.' ').AND.(trim(GCMvarname).ne."all"))
612	nbvar=nbvar+1
613	gcm_vars(nbvar)=GCMvarname
614	read(*,'(a50)') GCMvarname
615	enddo
616
617	if (GCMvarname.eq."all") then
618	nbvar=nbvarfile-Nnotprocessed
619	do v=Nnotprocessed+1,nbvarfile
620	status=nf90_inquire_variable(gcmfid,v,name=gcm_vars(v-Nnotprocessed))
621	enddo
622	! Variables names from the file are stored in gcm_vars()
623	nbvar=nbvarfile-Nnotprocessed
624	do v=1,nbvar
625	status=nf90_inquire_variable(gcmfid,v+Nnotprocessed,name=gcm_vars(v))
626	write(*,'(a9,1x,i2,1x,a1,1x,a64)') "variable ",v,":",gcm_vars(v)
627	enddo
628	else if(nbvar==0) then
629	write(,) "No variables to process in the GCM file... program stopped"
630	stop
631	endif ! of if (GCMvarname.eq."all")
632
633	!================================================================
634	! 2.1.2 Handle dust and wice opacities (first set-ups)
635	!================================================================
636	! 2nd part is in section 2.4.2
637	write(,)
638	! Load atmospheric density "rho"
639	if (dustok1.or.dustok2.or.dustok3.or.wiceok) then
640	! Check that the GCM file contains that variable
641	status=nf90_inq_varid(gcmfid,"rho",GCMvarid)
642	if (status.ne.nf90_noerr) then
643	write(,) "Failed to find variable rho in "//trim(gcmfile)
644	write(,) "No computation of opacities..."
645	dustok1 =.false.
646	dustok2 =.false.
647	dustok3 =.false.
648	wiceok =.false.
649	else
650	! Length allocation for each dimension of the 4D variable
651	allocate(rho(GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen))
652
653	! Load datasets
654	if (.not.is_stats) then
655	status=NF90_GET_VAR(gcmfid,GCMvarid,rho)
656	error_text="Failed to load rho"
657	call status_check(status,error_text)
658	else
659	! if it is a stats file, we load only the first sol, and then copy it to all the other sols
660	status=NF90_GET_VAR(gcmfid,GCMvarid,rho(:,:,:,1:GCMstatstimelen))
661	error_text="Failed to load rho"
662	call status_check(status,error_text)
663	! write(,)"GCMstatstimelen = ", GCMstatstimelen
664	! write(,)"GCMtimelen = ", GCMtimelen
665	do l=(GCMstatstimelen+1),GCMtimelen
666	if (modulo(l,GCMstatstimelen).ne.0) then
667	rho(:,:,:,l) = rho(:,:,:,modulo(l,GCMstatstimelen))
668	else ! if l is a multiple of GCMstatstimelen, since the index modulo(l,GCMstatstimelen)=0
669	! doesn't exist, we make a special case
670	rho(:,:,:,l) = rho(:,:,:,GCMstatstimelen)
671	endif
672	enddo
673	endif
674	write(,) "Variable rho loaded from the GCM file"
675	endif
676	endif ! dustok1.or.dustok2.or.dustok3.or.wiceok
677
678	! Dust and wice opacity booleans
679	if (dustok1.or.dustok2.or.dustok3) then
680	nbvar=nbvar+1
681	gcm_vars(nbvar)="dust"
682	endif
683
684	if (wiceok) then
685	nbvar=nbvar+1
686	gcm_vars(nbvar)="wice"
687	endif
688
689	!===============================================================================
690	! 2.2 Definition of LT variables from obsfile to outfile
691	!===============================================================================
692	! --> the day/night loop begins here
693
694	!****************** NOTA BENE (cf sections 2.2 and 4)***********************
695	! We execute the program a first time with the daytime values, and then a second
696	! time with the nighttime values.
697	!*******************************************************************************
698
699	dayornight = "dayside" ! we begin with daytime temperature
700	write(,)"" ; write(,) "Beginning the 1st loop, on daytime values"; write(,)""
701	DAY_OR_NIGHT: DO ! (the end of the loop is in section 4.)
702
703	SELECT CASE (dayornight)
704	CASE ("dayside")
705	OBSLTave_name = "dtimeave"
706	OBSLTmax_name = "dtimemax"
707	OBSLTmin_name = "dtimemin"
708	CASE ("nightside")
709	OBSLTave_name = "ntimeave"
710	OBSLTmax_name = "ntimemax"
711	OBSLTmin_name = "ntimemin"
712	END SELECT
713
714	!================================================================
715	! 2.2.1 Average of Local Time in the OBS bins
716	!================================================================
717	! Read the OBS file
718	!------------------
719	status=nf90_inq_varid(obsfid,trim(OBSLTave_name),LT_id)
720	error_text="Failed to find Observer local time ("//trim(OBSLTave_name)//") ID in "//trim(obsfile)
721	call status_check(status,error_text)
722	status=nf90_inquire_variable(obsfid,LT_id,dimids=LTshape)
723	error_text="Failed to get the dim shape of variable "//trim(OBSLTave_name)
724	call status_check(status,error_text)
725
726	! Length allocation for each dimension of the 3D variable
727	allocate(OBSLT(OBSlonlen,OBSlatlen,OBSLslen))
728
729	! Load datasets
730	status=NF90_GET_VAR(obsfid,LT_id,OBSLT)
731	error_text="Failed to load "//trim(OBSLTave_name)//" from the obsfile"
732	call status_check(status,error_text)
733	write(,) trim(OBSLTave_name)," loaded from the obsfile"
734
735	! Get LT missing_value attribute
736	status=nf90_get_att(obsfid,LT_id,"_FillValue",LTmiss_val)
737	error_text="Failed to load missing_value attribute"
738	call status_check(status,error_text)
739
740	! Create the variable in the outfile
741	!-----------------------------------
742	! Switch to netcdf define mode
743	status=nf90_redef(outfid)
744	error_text="Error: could not switch to define mode in the outfile"
745	call status_check(status,error_text)
746
747	! Definition of the variable
748	status=NF90_DEF_VAR(outfid,trim(OBSLTave_name),nf90_float,LTshape,LT_id)
749	error_text="Error: could not define the variable "//trim(OBSLTave_name)//" in the outfile"
750	call status_check(status,error_text)
751
752	! Write the attributes
753	select case (dayornight)
754	case ("dayside")
755	status=nf90_put_att(outfid,LT_id,"long_name","Average local time in bin - day side [6h, 18h]")
756	case ("nightside")
757	status=nf90_put_att(outfid,LT_id,"long_name","Average local time in bin - night side [18h, 6h]")
758	end select
759	status=nf90_put_att(outfid,LT_id,"units","hours")
760	status=nf90_put_att(outfid,LT_id,"_FillValue",LTmiss_val)
761
762	! End the netcdf define mode (and thus enter the "data writing" mode)
763	status=nf90_enddef(outfid)
764	error_text="Error: could not close the define mode of the outfile"
765	call status_check(status,error_text)
766
767	! Write the data in the output file
768	status = NF90_PUT_VAR(outfid, LT_id, OBSLT) ! write the MCS d/ntimeave as the output Local Time
769	error_text="Error: could not write "//trim(OBSLTave_name)//" data in the outfile"
770	call status_check(status,error_text)
771
772	write(,)"Local Time (",trim(OBSLTave_name),") has been created in the outfile"
773	write(*,'(" with missing_value attribute : ",1pe12.5)')LTmiss_val
774
775	!================================================================
776	! 2.2.2 Maximum of Local Time in the OBS bins
777	!================================================================
778	! Read the OBS file
779	!------------------
780	status=nf90_inq_varid(obsfid,trim(OBSLTmax_name),LT_id)
781	error_text="Failed to find Observer max local time ("//trim(OBSLTmax_name)//") ID in "//trim(obsfile)
782	call status_check(status,error_text)
783	status=nf90_inquire_variable(obsfid,LT_id,dimids=LTshape)
784	error_text="Failed to get the dim shape of variable "//trim(OBSLTmax_name)
785	call status_check(status,error_text)
786
787	! Length allocation for each dimension of the 3D variable
788	allocate(OBSLTmax(OBSlonlen,OBSlatlen,OBSLslen))
789
790	! Load datasets
791	status=NF90_GET_VAR(obsfid,LT_id,OBSLTmax)
792	error_text="Failed to load "//trim(OBSLTmax_name)//" from the obsfile"
793	call status_check(status,error_text)
794	write(,) trim(OBSLTmax_name)," loaded from the obsfile"
795
796	!================================================================
797	! 2.2.3 Minimum of Local Time in the OBS bins
798	!================================================================
799	! Read the OBS file
800	!------------------
801	status=nf90_inq_varid(obsfid,trim(OBSLTmin_name),LT_id)
802	error_text="Failed to find Observer min local time ("//trim(OBSLTmin_name)//") ID in "//trim(obsfile)
803	call status_check(status,error_text)
804	status=nf90_inquire_variable(obsfid,LT_id,dimids=LTshape)
805	error_text="Failed to obtain information on variable "//trim(OBSLTmin_name)
806	call status_check(status,error_text)
807
808	! Length allocation for each dimension of the 3D variable
809	allocate(OBSLTmin(OBSlonlen,OBSlatlen,OBSLslen))
810
811	! Load datasets
812	status=NF90_GET_VAR(obsfid,LT_id,OBSLTmin)
813	error_text="Failed to load "//trim(OBSLTmin_name)//" from the obsfile"
814	call status_check(status,error_text)
815	write(,) trim(OBSLTmin_name)," loaded from the obsfile"
816	write(,)""
817
818	!===============================================================================
819	! 2.3 Definition of numbin variables from obsfile to outfile
820	!===============================================================================
821	!================================================================
822	! 2.3.1 Number of values in the OBS "temp" bins
823	!================================================================
824	SELECT CASE (dayornight)
825	CASE ("dayside")
826	numbin_name = "dnumbintemp"
827	CASE ("nightside")
828	numbin_name = "nnumbintemp"
829	END SELECT
830
831	! Read the OBS file
832	!------------------
833	status=nf90_inq_varid(obsfid,trim(numbin_name),numbin_id)
834	error_text="Failed to find Observer nb of temp values in bin ("//trim(numbin_name)//")'s ID in "//trim(obsfile)
835	call status_check(status,error_text)
836	status=nf90_inquire_variable(obsfid,numbin_id,dimids=numbinshape)
837	error_text="Failed to obtain information on variable "//trim(numbin_name)
838	call status_check(status,error_text)
839
840	! Length allocation for each dimension of the 4D variable
841	allocate(numbin(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
842
843	! Load datasets
844	status=NF90_GET_VAR(obsfid,numbin_id,numbin)
845	error_text="Failed to load "//trim(numbin_name)//" from the obsfile"
846	call status_check(status,error_text)
847	write(,) trim(numbin_name)," loaded from the obsfile"
848
849	! Create the variable in the outfile
850	!-----------------------------------
851	! Switch to netcdf define mode
852	status=nf90_redef(outfid)
853	error_text="Error: could not switch to define mode in the outfile"
854	call status_check(status,error_text)
855
856	! Definition of the variable
857	status=NF90_DEF_VAR(outfid,trim(numbin_name),nf90_float,numbinshape,numbin_id)
858	error_text="Error: could not define the variable "//trim(numbin_name)//" in the outfile"
859	call status_check(status,error_text)
860
861	! Write the attributes
862	select case (dayornight)
863	case ("dayside")
864	status=nf90_put_att(outfid,numbin_id,"long_name","Number of temp values in bin - day side")
865	case ("nightside")
866	status=nf90_put_att(outfid,numbin_id,"long_name","Number of temp values in bin - night side")
867	end select
868
869	! End the netcdf define mode (and thus enter the "data writing" mode)
870	status=nf90_enddef(outfid)
871	error_text="Error: could not close the define mode of the outfile"
872	call status_check(status,error_text)
873
874	! Write the data in the output file
875	status = NF90_PUT_VAR(outfid, numbin_id, numbin)
876	error_text="Error: could not write "//trim(numbin_name)//" data in the outfile"
877	call status_check(status,error_text)
878
879	write(,)"Number of temp values in bin (",trim(numbin_name),") has been created in the outfile"
880	write(,)""
881	deallocate(numbin)
882
883	!================================================================
884	! 2.3.2 Number of values in the OBS "dust" bins
885	!================================================================
886	SELECT CASE (dayornight)
887	CASE ("dayside")
888	numbin_name = "dnumbindust"
889	CASE ("nightside")
890	numbin_name = "nnumbindust"
891	END SELECT
892
893	! Read the OBS file
894	!------------------
895	status=nf90_inq_varid(obsfid,trim(numbin_name),numbin_id)
896	error_text="Failed to find Observer nb of dust values in bin ("//trim(numbin_name)//")'s ID in "//trim(obsfile)
897	call status_check(status,error_text)
898	status=nf90_inquire_variable(obsfid,numbin_id,dimids=numbinshape)
899	error_text="Failed to obtain information on variable "//trim(numbin_name)
900	call status_check(status,error_text)
901
902	! Length allocation for each dimension of the 4D variable
903	allocate(numbin(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
904
905	! Load datasets
906	status=NF90_GET_VAR(obsfid,numbin_id,numbin)
907	error_text="Failed to load "//trim(numbin_name)//" from the obsfile"
908	call status_check(status,error_text)
909	write(,) trim(numbin_name)," loaded from the obsfile"
910
911	! Create the variable in the outfile
912	!-----------------------------------
913	! Switch to netcdf define mode
914	status=nf90_redef(outfid)
915	error_text="Error: could not switch to define mode in the outfile"
916	call status_check(status,error_text)
917
918	! Definition of the variable
919	status=NF90_DEF_VAR(outfid,trim(numbin_name),nf90_float,numbinshape,numbin_id)
920	error_text="Error: could not define the variable "//trim(numbin_name)//" in the outfile"
921	call status_check(status,error_text)
922
923	! Write the attributes
924	select case (dayornight)
925	case ("dayside")
926	status=nf90_put_att(outfid,numbin_id,"long_name","Number of dust values in bin - day side")
927	case ("nightside")
928	status=nf90_put_att(outfid,numbin_id,"long_name","Number of dust values in bin - night side")
929	end select
930
931	! End the netcdf define mode (and thus enter the "data writing" mode)
932	status=nf90_enddef(outfid)
933	error_text="Error: could not close the define mode of the outfile"
934	call status_check(status,error_text)
935
936	! Write the data in the output file
937	status = NF90_PUT_VAR(outfid, numbin_id, numbin)
938	error_text="Error: could not write "//trim(numbin_name)//" data in the outfile"
939	call status_check(status,error_text)
940
941	write(,)"Number of dust values in bin (",trim(numbin_name),") has been created in the outfile"
942	write(,)""
943	deallocate(numbin)
944
945	!================================================================
946	! 2.3.3 Number of values in the OBS "wice" bins
947	!================================================================
948	SELECT CASE (dayornight)
949	CASE ("dayside")
950	numbin_name = "dnumbinwice"
951	CASE ("nightside")
952	numbin_name = "nnumbinwice"
953	END SELECT
954
955	! Read the OBS file
956	!------------------
957	status=nf90_inq_varid(obsfid,trim(numbin_name),numbin_id)
958	error_text="Failed to find Observer nb of wice values in bin ("//trim(numbin_name)//")'s ID in "//trim(obsfile)
959	call status_check(status,error_text)
960	status=nf90_inquire_variable(obsfid,numbin_id,dimids=numbinshape)
961	error_text="Failed to obtain information on variable "//trim(numbin_name)
962	call status_check(status,error_text)
963
964	! Length allocation for each dimension of the 4D variable
965	allocate(numbin(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
966
967	! Load datasets
968	status=NF90_GET_VAR(obsfid,numbin_id,numbin)
969	error_text="Failed to load "//trim(numbin_name)//" from the obsfile"
970	call status_check(status,error_text)
971	write(,) trim(numbin_name)," loaded from the obsfile"
972
973	! Create the variable in the outfile
974	!-----------------------------------
975	! Switch to netcdf define mode
976	status=nf90_redef(outfid)
977	error_text="Error: could not switch to define mode in the outfile"
978	call status_check(status,error_text)
979
980	! Definition of the variable
981	status=NF90_DEF_VAR(outfid,trim(numbin_name),nf90_float,numbinshape,numbin_id)
982	error_text="Error: could not define the variable "//trim(numbin_name)//" in the outfile"
983	call status_check(status,error_text)
984
985	! Write the attributes
986	select case (dayornight)
987	case ("dayside")
988	status=nf90_put_att(outfid,numbin_id,"long_name","Number of wice values in bin - day side")
989	case ("nightside")
990	status=nf90_put_att(outfid,numbin_id,"long_name","Number of wice values in bin - night side")
991	end select
992
993	! End the netcdf define mode (and thus enter the "data writing" mode)
994	status=nf90_enddef(outfid)
995	error_text="Error: could not close the define mode of the outfile"
996	call status_check(status,error_text)
997
998	! Write the data in the output file
999	status = NF90_PUT_VAR(outfid, numbin_id, numbin)
1000	error_text="Error: could not write "//trim(numbin_name)//" data in the outfile"
1001	call status_check(status,error_text)
1002
1003	write(,)"Number of wice values in bin (",trim(numbin_name),") has been created in the outfile"
1004	write(,)""
1005	deallocate(numbin)
1006
1007
1008	!===============================================================================
1009	! 2.4 Opening of the GCM variables
1010	!===============================================================================
1011	! --> the var loop begins here
1012
1013	VAR: DO v=1,nbvar ! LOOP ON ALL THE GCM VARIABLES TO PROCESS
1014	! (the end of the loop is in section 3.4)
1015
1016	GCMvarname = gcm_vars(v)
1017
1018	! Detect the dust and wice opacities special cases
1019	if (trim(GCMvarname).eq."dust") then
1020	if (dustok1) then ! "dso" is detected in gcmfile
1021	GCMvarname="dso"
1022	dustok2=.false.
1023	dustok3=.false.
1024	else if (dustok2) then ! "dsodust" is detected in gcmfile
1025	GCMvarname="dsodust"
1026	dustok3=.false.
1027	else if (dustok3) then ! "dustq" is detected in gcmfile
1028	GCMvarname="dustq"
1029	endif
1030	write(,) "Computing dust opacity..."
1031	endif
1032	if (trim(GCMvarname).eq."wice") then ! "h2o_ice" detected in gcmfile
1033	GCMvarname="h2o_ice"
1034	write(,) "Computing water ice opacity..."
1035	endif
1036
1037	!================================================================
1038	! 2.4.1 Generic reading of the variable
1039	!================================================================
1040	! Check that the GCM file contains that variable
1041	status=nf90_inq_varid(gcmfid,trim(GCMvarname),GCMvarid)
1042	if (status.ne.nf90_noerr) then
1043	write(,) "Failed to find variable "//trim(GCMvarname)//" in "//trim(gcmfile)
1044	write(,) "We'll skip that variable..."
1045	CYCLE VAR ! go directly to the next variable
1046	endif
1047
1048	! Sanity checks on the variable
1049	status=nf90_inquire_variable(gcmfid,GCMvarid,ndims=nbdim,dimids=GCMvarshape)
1050	error_text="Failed to obtain information on variable "//trim(GCMvarname)
1051	call status_check(status,error_text)
1052
1053	! Check that it is a 4D variable
1054	if (nbdim.ne.4) then
1055	write(,) "Error:",trim(GCMvarname)," is not a 4D variable"
1056	write(,) "We'll skip that variable...";write(,)""
1057	CYCLE VAR ! go directly to the next variable
1058	endif
1059	! Check that its dimensions are indeed lon,lat,alt,time (in the right order)
1060	if (GCMvarshape(1).ne.lon_dimid_gcm) then
1061	write(,) "Error, expected first dimension of ",trim(GCMvarname)," to be longitude!"
1062	write(,) "We'll skip that variable..."
1063	CYCLE VAR ! go directly to the next variable
1064	endif
1065	if (GCMvarshape(2).ne.lat_dimid_gcm) then
1066	write(,) "Error, expected second dimension of ",trim(GCMvarname)," to be latitude!"
1067	write(,) "We'll skip that variable..."
1068	CYCLE VAR ! go directly to the next variable
1069	endif
1070	if (GCMvarshape(3).ne.alt_dimid_gcm) then
1071	write(,) "Error, expected third dimension of ",trim(GCMvarname)," to be altitude!"
1072	write(,) "We'll skip that variable..."
1073	CYCLE VAR ! go directly to the next variable
1074	endif
1075	if (GCMvarshape(4).ne.time_dimid_gcm) then
1076	write(,) "Error, expected fourth dimension of ",trim(GCMvarname)," to be time!"
1077	write(,) "We'll skip that variable..."
1078	CYCLE VAR ! go directly to the next variable
1079	endif
1080
1081	! Length allocation for each dimension of the 4D variable
1082	allocate(GCM_var(GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen))
1083
1084	! Load datasets
1085	if (.not.is_stats) then
1086	status=NF90_GET_VAR(gcmfid,GCMvarid,GCM_var)
1087	error_text="Failed to load "//trim(GCMvarname)
1088	call status_check(status,error_text)
1089	else
1090	! if it is a stats file, we load only the first sol, and then copy it to all the other sols
1091	status=NF90_GET_VAR(gcmfid,GCMvarid,GCM_var(:,:,:,1:GCMstatstimelen))
1092	error_text="Failed to load "//trim(GCMvarname)
1093	call status_check(status,error_text)
1094	! write(,)"GCMstatstimelen = ", GCMstatstimelen
1095	! write(,)"GCMtimelen = ", GCMtimelen
1096	do l=(GCMstatstimelen+1),GCMtimelen
1097	if (modulo(l,GCMstatstimelen).ne.0) then
1098	GCM_var(:,:,:,l) = GCM_var(:,:,:,modulo(l,GCMstatstimelen))
1099	else ! if l is a multiple of GCMstatstimelen, since the index modulo(l,GCMstatstimelen)=0
1100	! doesn't exist, we make a special case
1101	GCM_var(:,:,:,l) = GCM_var(:,:,:,GCMstatstimelen)
1102	endif
1103	enddo
1104	endif
1105	write(,) "Variable ",trim(GCMvarname)," loaded from the GCM file"
1106
1107	! Get dataset's missing_value attribute
1108	status=nf90_get_att(gcmfid,GCMvarid,"missing_value",GCMmiss_val)
1109	error_text="Failed to load missing_value attribute"
1110	call status_check(status,error_text)
1111
1112	! Get other variable's attributes
1113	status=nf90_get_att(gcmfid,GCMvarid,"long_name",long_name)
1114	if (status.ne.nf90_noerr) then
1115	! if no attribute "long_name", try "title"
1116	status=nf90_get_att(gcmfid,GCMvarid,"title",long_name)
1117	endif
1118	status=nf90_get_att(gcmfid,GCMvarid,"units",units)
1119
1120	!================================================================
1121	! 2.4.2 Handle dust and wice opacities (second part)
1122	!================================================================
1123	! DUST
1124	!-----
1125	if (trim(gcm_vars(v)).eq."dust") then
1126
1127	IF (dustok1.or.dustok2) THEN
1128	! Dust opacity computed from its density-scaled opacity
1129	! write(,)long_name(index(long_name,"(")+1:index(long_name,")")-1)
1130	do i=1,GCMlonlen
1131	do j=1,GCMlatlen
1132	do k=1,GCMaltlen
1133	do l=1,GCMtimelen
1134	if (GCM_var(i,j,k,l).ne.GCMmiss_val) then
1135	! Multiply by rho to have opacity [1/km]
1136	GCM_var(i,j,k,l) = GCM_var(i,j,k,l) * rho(i,j,k,l) *1000.
1137
1138	if (long_name(index(long_name,"(")+1:index(long_name,")")-1).eq."TES") then
1139	! The density-scaled opacity was calibrated on TES wavelength (9.3um)
1140	! so we must recalibrate it to MCS wavelength (21.6um) using recalibration
1141	! coefficients from Montabone et al. 2015, section 2.3
1142	GCM_var(i,j,k,l) = 1.3/2.7 * GCM_var(i,j,k,l)
1143	endif
1144	endif
1145	enddo
1146	enddo
1147	enddo
1148	enddo
1149
1150	long_name = "IR Dust opacity (from DSO)"
1151
1152	ELSE IF (dustok3) THEN
1153	! Dust opacity computed from its mass mixing ratio
1154	do i=1,GCMlonlen
1155	do j=1,GCMlatlen
1156	do k=1,GCMaltlen
1157	do l=1,GCMtimelen
1158	if (GCM_var(i,j,k,l).ne.GCMmiss_val) then
1159	! Opacity is computed from the equation of Heavens et al. 2014, section 2.3,
1160	! assuming a rho_dust=3000 kg/m3 and an effective radius reff=1.06 microns
1161	! + the opacity is here in 1/km and the MMR in kg/kg
1162	GCM_var(i,j,k,l) = GCM_var(i,j,k,l) * rho(i,j,k,l) / 0.012 * 1000
1163	endif
1164	enddo
1165	enddo
1166	enddo
1167	enddo
1168
1169	long_name = "IR Dust opacity (from MMR)"
1170
1171	ENDIF
1172
1173	GCMvarname = gcm_vars(v) ! reput the right name in GCMvarname
1174	units = "opacity/km"
1175	endif ! trim(gcm_vars(v)).eq."dust"
1176
1177
1178	! WICE
1179	!-----
1180	if (trim(gcm_vars(v)).eq."wice") then
1181	! Water ice opacity computed from its mass mixing ratio
1182	do i=1,GCMlonlen
1183	do j=1,GCMlatlen
1184	do k=1,GCMaltlen
1185	do l=1,GCMtimelen
1186	if (GCM_var(i,j,k,l).ne.GCMmiss_val) then
1187	! Opacity at MCS wavelength (11.9um) is computed from an equation
1188	! similar to the one of Heavens et al. 2014, section 2.3.
1189	! We assume a rho_wice=920 kg/m3, an effective radius reff=3um,
1190	! an extinction coefficient Qext(wvl,reff)=1.54471
1191	GCM_var(i,j,k,l) = 7501.54471 GCM_var(i,j,k,l) * rho(i,j,k,l) / (920*3e-6)
1192	endif
1193	enddo
1194	enddo
1195	enddo
1196	enddo
1197
1198	long_name = "IR Water ice opacity (from MMR)"
1199
1200	GCMvarname = gcm_vars(v) ! reput the right name in GCMvarname
1201	units = "opacity/km"
1202	endif ! trim(gcm_vars(v)).eq."wice"
1203
1204	!===============================================================================
1205	! 2.5 Opening of the associated MCS variables
1206	!===============================================================================
1207	! Observer variables to extract :
1208	IF ((index(GCMvarname,"dust").ne.0).or.(index(GCMvarname,"dso").ne.0)) THEN
1209	! if the variable name contains "dust" or "dso". Especially for the targeted variables :
1210	! dustq,dustN,dsodust,reffdust,opadust, and their equivalents for stormdust & topdust
1211	OBSvarname = "dust"
1212	numbin_name = "numbindust"
1213	ELSE IF ((trim(GCMvarname).eq."h2o_ice").or.(trim(GCMvarname).eq."wice") &
1214	.or.(trim(GCMvarname).eq."reffice").or.(trim(GCMvarname).eq."opawice")) THEN
1215	OBSvarname = "wice"
1216	numbin_name = "numbinwice"
1217	ELSE ! default case is temp binning, since it contains the most values
1218	OBSvarname = "temp"
1219	numbin_name = "numbintemp"
1220	ENDIF
1221
1222	SELECT CASE (dayornight)
1223	CASE ("dayside")
1224	OBSvarname = "d"//OBSvarname
1225	numbin_name = "d"//numbin_name
1226	CASE ("nightside")
1227	OBSvarname = "n"//OBSvarname
1228	numbin_name = "n"//numbin_name
1229	END SELECT
1230
1231	!================================================================
1232	! 2.5.1 MCS reference variable (for the missing values)
1233	!================================================================
1234	! Check that the observation file contains that variable
1235	status=nf90_inq_varid(obsfid,trim(OBSvarname),OBSvarid)
1236	error_text="Failed to find variable "//trim(OBSvarname)//" in "//trim(obsfile)
1237	call status_check(status,error_text)
1238
1239	! Sanity checks on the variable
1240	status=nf90_inquire_variable(obsfid,OBSvarid,ndims=nbdim,dimids=OBSvarshape)
1241	error_text="Failed to obtain information on variable "//trim(OBSvarname)
1242	call status_check(status,error_text)
1243
1244	! Check that it is a 4D variable
1245	if (nbdim.ne.4) then
1246	write(,) "Error, expected a 4D (lon-lat-alt-time) variable for ",trim(OBSvarname)
1247	stop
1248	endif
1249	! Check that its dimensions are indeed lon,lat,alt,time (in the right order)
1250	if (OBSvarshape(1).ne.lon_dimid_obs) then
1251	write(,) "Error, expected first dimension of ",trim(OBSvarname)," to be longitude!"
1252	stop
1253	endif
1254	if (OBSvarshape(2).ne.lat_dimid_obs) then
1255	write(,) "Error, expected second dimension of ",trim(OBSvarname)," to be latitude!"
1256	stop
1257	endif
1258	if (OBSvarshape(3).ne.alt_dimid_obs) then
1259	write(,) "Error, expected third dimension of ",trim(OBSvarname)," to be altitude!"
1260	stop
1261	endif
1262	if (OBSvarshape(4).ne.time_dimid_obs) then
1263	write(,) "Error, expected fourth dimension of ",trim(OBSvarname)," to be time!"
1264	stop
1265	endif
1266
1267	! Length allocation for each dimension of the 4D variable
1268	allocate(OBS_var(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
1269
1270	! Load datasets
1271	status=NF90_GET_VAR(obsfid,OBSvarid,OBS_var)
1272	error_text="Failed to load "//trim(OBSvarname)//" from the obsfile"
1273	call status_check(status,error_text)
1274	write(,) trim(OBSvarname)," loaded from the obsfile as reference variable"
1275
1276	! Get OBS_var missing_value attribute
1277	status=nf90_get_att(obsfid,OBSvarid,"_FillValue",OBSmiss_val)
1278	error_text="Failed to load missing_value attribute"
1279	call status_check(status,error_text)
1280
1281	!================================================================
1282	! 2.5.2 Number of values in the OBS bin (for the sol binning)
1283	!================================================================
1284	! Check that the observation file contains that variable
1285	status=nf90_inq_varid(obsfid,trim(numbin_name),numbin_id)
1286
1287	! Checks have already been done in section 2.3
1288
1289	! Length allocation for each dimension of the 4D variable
1290	allocate(numbin(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
1291
1292	! Load datasets
1293	status=NF90_GET_VAR(obsfid,numbin_id,numbin)
1294
1295	!===============================================================================
1296	! 2.6 Definition of GCM variables in outfile
1297	!===============================================================================
1298	! Switch to netcdf define mode
1299	status=nf90_redef(outfid)
1300	error_text="Error: could not switch to define mode in the outfile"
1301	call status_check(status,error_text)
1302
1303	! Definition of the variable
1304	SELECT CASE (dayornight)
1305	CASE ("dayside")
1306	outvarname = "d"//GCMvarname
1307	CASE ("nightside")
1308	outvarname = "n"//GCMvarname
1309	END SELECT
1310	status=NF90_DEF_VAR(outfid,trim(outvarname),nf90_float,OBSvarshape,outvarid)
1311	error_text="Error: could not define the variable "//trim(outvarname)//" in the outfile"
1312	call status_check(status,error_text)
1313
1314	! Write the attributes
1315	SELECT CASE (dayornight)
1316	CASE ("dayside")
1317	long_name = trim(long_name)//" - day side"
1318	status=nf90_put_att(outfid,outvarid,"long_name",long_name)
1319	CASE ("nightside")
1320	long_name = trim(long_name)//" - night side"
1321	status=nf90_put_att(outfid,outvarid,"long_name",long_name)
1322	END SELECT
1323	status=nf90_put_att(outfid,outvarid,"units",units)
1324	status=nf90_put_att(outfid,outvarid,"_FillValue",OBSmiss_val)
1325	comment = "Reference numbin: "//trim(numbin_name)
1326	status=nf90_put_att(outfid,outvarid,"comment",comment)
1327
1328	write(,)trim(outvarname)," has been created in the outfile"
1329	write(*,'(" with missing_value attribute : ",1pe12.5)')OBSmiss_val
1330	write(,)""
1331
1332	! End the netcdf define mode (and thus enter the "data writing" mode)
1333	status=nf90_enddef(outfid)
1334	error_text="Error: could not close the define mode of the outfile"
1335	call status_check(status,error_text)
1336
1337	allocate(outvar(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
1338
1339
1340	!===================================================================================
1341	! 3. EXTRACTION OF THE VARIABLE
1342	!===================================================================================
1343	!===============================================================================
1344	! 3.1 Do some checks and preparations before the extraction
1345	!===============================================================================
1346	Ls: do l=1,OBSLslen ! loop on all observed Solar longitudes
1347	Ls_val=OBSLs(l) ! Ls_val=center of the output bin
1348	if ((Ls_val.lt.0.).or.(Ls_val.gt.360.)) then
1349	write(,) "Unexpected value for OBSLs: ",Ls_val
1350	stop
1351	endif
1352
1353	! Convert the Ls bin into a sol interval on which the binning is done :
1354	!----------------------------------------------------------------------
1355	!-> get the index of the maximum value of GCM sol (m_maxsol) that is lower than Ls bin's superior bound (maxsol)
1356	call ls2sol(Ls_val+OBSdeltaLs/2.,maxsol)
1357	m_maxsol=1
1358	do while (((GCMtime(m_maxsol)+0.5).lt.maxsol).AND.(m_maxsol.le.(GCMtimelen-1)))
1359	!! The +0.5 is there to take into account the whole planet (lon=[-180°;180°]) and not just the lon=0° from the GCM
1360	m_maxsol=m_maxsol+1
1361	enddo
1362	!-> get the index of the minimum value of GCM sol (m_minsol) that is greater than Ls bin's inferior bound (minsol)
1363	call ls2sol(Ls_val-OBSdeltaLs/2.,minsol)
1364	m_minsol=1
1365	do while (((GCMtime(m_minsol)-0.5).le.minsol).AND.(m_minsol.le.(GCMtimelen-1)))
1366	!! Same comment for the -0.5
1367	m_minsol=m_minsol+1
1368	enddo
1369	if (m_minsol.gt.m_maxsol) then
1370	write(,) "No value in gcmfile between sol=",minsol," and sol=",maxsol," (Ls=",Ls_val,"Â°)"
1371	! Write a missing_value to output
1372	outvar(:,:,:,l) = OBSmiss_val
1373	CYCLE Ls ! go directly to the next Ls
1374	endif
1375	! Get all the integer values of GCM sols that fit in this interval
1376	solcount=floor(GCMtime(m_maxsol))-ceiling(GCMtime(m_minsol))+1
1377	! sols that are not fully in the interval are not counted
1378	allocate(int_sol_list(solcount))
1379	! write(,) "GCMminsol=", GCMtime(m_minsol)
1380	! write(,)"GCMmaxsol=", GCMtime(m_maxsol)
1381	do m=1,solcount
1382	int_sol_list(m)=ceiling(GCMtime(m_minsol)) + m-1
1383	enddo
1384	! write(,)"int_sol_list=",int_sol_list
1385
1386
1387	latitude: do j=1,OBSlatlen ! loop on all observed latitudes
1388	lat_val=OBSlat(j)
1389	if ((lat_val.lt.-90.).or.(lat_val.gt.90.)) then
1390	write(,) "Unexpected value for OBSlat: ",lat_val
1391	stop
1392	endif
1393
1394	longitude: do i=1,OBSlonlen ! loop on all observed longitudes
1395	lon_val=OBSlon(i)
1396	if ((lon_val.lt.-360.).or.(lon_val.gt.360.)) then
1397	write(,) "Unexpected value for lon_val: ",lon_val
1398	stop
1399	endif
1400	! We want lon_val in [-180:180] for the subroutine extraction
1401	if (lon_val.lt.-180.) lon_val=lon_val+360.
1402	if (lon_val.gt.180.) lon_val=lon_val-360.
1403
1404	LT_val=OBSLT(i,j,l) ! find the Observer average LT value at bin(lon_val, lat_val, Ls_val)
1405
1406	if ((LT_val.lt.0.).or.(LT_val.gt.24.)) then
1407	if (LT_val.eq.LTmiss_val) then
1408	! write(,) "Missing value in obsfile for LT_val"
1409	! Write a missing_value to output
1410	outvar(i,j,:,l) = OBSmiss_val
1411	CYCLE longitude ! go directly to the next longitude
1412	else
1413	write(,) "Unexpected value for LT_val: ",LT_val
1414	stop
1415	endif
1416	endif
1417
1418	altitude: do k=1,OBSaltlen ! loop on all observed altitudes
1419	alt_val=OBSalt(k)
1420	if (OBS_var(i,j,k,l).eq.OBSmiss_val) then
1421	! write(,) "Missing value in obsfile for ",OBSvarname
1422	! Write a missing_value to output
1423	outvar(i,j,k,l) = OBSmiss_val
1424	CYCLE altitude ! go directly to the next altitude
1425	endif
1426
1427	!===============================================================================
1428	! 3.2 Compute GCM sol date corresponding to Observer Ls (via m_(min/max)sol)
1429	! and LT (via OBSLT(min/max))
1430	!===============================================================================
1431	LTcount=floor(numbin(i,j,k,l)) ! find the Observer number of temp values
1432	! at bin(lon_val,lat_val,alt_val,Ls_val)
1433	if (LTcount.eq.0.) then
1434	! Write a missing_value to output
1435	outvar(i,j,k,l) = OBSmiss_val
1436	CYCLE altitude ! go directly to the next altitude
1437	endif
1438	if (LTcount.lt.0.) then
1439	write(,) "Unexpected value for LTcount: ",LTcount
1440	stop
1441	endif
1442
1443	! Generate the sol list for the interpolation
1444	allocate(sol_list(solcount*LTcount))
1445	call gen_sol_list(solcount,int_sol_list,LTcount,LT_val,OBSLTmax(i,j,l),&
1446	OBSLTmin(i,j,l),lon_val,LTmod,dayornight,&
1447	sol_list)
1448
1449	solerrcount=0
1450	solbinned_value=0
1451	sol_bin: do m=1,solcount*LTcount ! loop on all GCM sols of the bin
1452	sol=sol_list(m)
1453	! write(,)"sol=",sol
1454	!===============================================================================
1455	! 3.3 Do the interpolation and binning for the given location
1456	!===============================================================================
1457	call extraction(lon_val,lat_val,alt_val,sol,&
1458	GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen,&
1459	GCMlon,GCMlat,GCMalt,GCMtime,&
1460	GCM_var,GCMmiss_val,GCMalttype,GCMvarname,extr_value)
1461
1462	if (extr_value.eq.GCMmiss_val) then
1463	! write(,) "Missing value in gcmfile at lon=",lon_val,"; lat=",lat_val,"; alt=",alt_val,"; sol=",sol
1464	solerrcount=solerrcount+1
1465	CYCLE sol_bin ! go directly to the next GCM sol of the bin
1466	endif
1467	solbinned_value=solbinned_value+extr_value
1468
1469	enddo sol_bin ! end loop on all GCM sols of the bin
1470	if ((solcount*LTcount-solerrcount).ne.0) then
1471	solbinned_value=solbinned_value/(solcount*LTcount-solerrcount)
1472	else
1473	! write(,)"No GCM value in this sol bin"
1474	solbinned_value=OBSmiss_val
1475	endif
1476	! Write value to output
1477	outvar(i,j,k,l)=solbinned_value
1478
1479	errcount=errcount+solerrcount
1480
1481	deallocate(sol_list)
1482
1483	enddo altitude ! end loop on observed altitudes
1484	enddo longitude ! end loop on observed longitudes
1485	enddo latitude ! end loop on observed latitudes
1486
1487	deallocate(int_sol_list)
1488
1489	enddo Ls ! end loop on observed Solar longitudes
1490
1491	! write(,)"Nb of GCM missing values :",errcount
1492
1493	!===============================================================================
1494	! 3.4 Write the data in the netcdf output file
1495	!===============================================================================
1496	status = nf90_put_var(outfid, outvarid, outvar)
1497	error_text="Error: could not write "//trim(outvarname)//" data in the outfile"
1498	call status_check(status,error_text)
1499
1500	! Deallocations before going to the next GCM variable
1501	deallocate(GCM_var)
1502	deallocate(OBS_var)
1503	deallocate(numbin)
1504	deallocate(outvar)
1505	ENDDO VAR ! end loop on variables
1506
1507	!===================================================================================
1508	! 4. END OF THE DAY/NIGHT LOOP
1509	!===================================================================================
1510	IF (dayornight.EQ."dayside") THEN
1511	! this is the end of the first loop (on daytime values)
1512	! and we still have to loop on nighttime values
1513	dayornight="nightside"
1514	deallocate(OBSLT)
1515	deallocate(OBSLTmax)
1516	deallocate(OBSLTmin)
1517	write(,)""
1518	write(,) "Beginning the 2nd loop, on nighttime values"; write(,)""
1519	CYCLE DAY_OR_NIGHT
1520	ELSE ! i.e. dayornight="nightside"
1521	! this is the end of the second loop (on nighttime values)
1522	! and thus the end of the program
1523	EXIT DAY_OR_NIGHT
1524	ENDIF
1525	ENDDO DAY_OR_NIGHT ! end of the day/night loop that begins in section 1.3
1526
1527	!===================================================================================
1528	! 5. CLOSE THE FILES
1529	!===================================================================================
1530	status = nf90_close(gcmfid)
1531	error_text="Error: could not close file "//trim(gcmfile)
1532	call status_check(status,error_text)
1533	status = nf90_close(obsfid)
1534	error_text="Error: could not close file "//trim(obsfile)
1535	call status_check(status,error_text)
1536	status = nf90_close(outfid)
1537	error_text="Error: could not close file "//trim(outfile)
1538	call status_check(status,error_text)
1539
1540	write(,)"";write(,)"-> Program simu_MCS completed!"
1541
1542	end program simu_MCS
1543
1544
1545
1546	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1547	!===================================================================================================
1548	! Subroutines
1549	!===================================================================================================
1550
1551	subroutine extraction(lon,lat,alt,sol,&
1552	lonlen,latlen,altlen,timelen,&
1553	longitude,latitude,altitude,time,&
1554	field,missing_value,alttype,varname,value)
1555
1556	implicit none
1557	!================================================================
1558	! Arguments:
1559	!================================================================
1560	real,intent(in) :: lon ! sought longitude (deg, in [-180:180])
1561	real,intent(in) :: lat ! sought latitude (deg, in [-90:90])
1562	real,intent(in) :: alt ! sought altitude (m or Pa)
1563	real,intent(in) :: sol ! sought date (sols)
1564	integer,intent(in) :: lonlen
1565	integer,intent(in) :: latlen
1566	integer,intent(in) :: altlen
1567	integer,intent(in) :: timelen
1568	real,intent(in) :: longitude(lonlen)
1569	real,intent(in) :: latitude(latlen)
1570	real,intent(in) :: altitude(altlen)
1571	real,intent(in) :: time(timelen)
1572	real,intent(in) :: field(lonlen,latlen,altlen,timelen)
1573	real,intent(in) :: missing_value ! default value for "no data"
1574	character,intent(in) :: alttype ! altitude coord. type:'z' (altitude, m)
1575	! 'p' (pressure, Pa)
1576	character(len=*),intent(in) :: varname ! variable name (in GCM file)
1577	real,intent(out) :: value
1578
1579	!================================================================
1580	! Local variables:
1581	!================================================================
1582	real,save :: prev_lon=-99 ! previous value of 'lon' routine was called with
1583	real,save :: prev_lat=-99 ! previous value of 'lat' routine was called with
1584	real,save :: prev_alt=-9.e20 ! ! previous value of 'alt'
1585	real,save :: prev_sol=-99 ! previous value of 'sol' routine was called with
1586
1587	! encompasing indexes:
1588	integer,save :: ilon_inf=-1,ilon_sup=-1 ! along longitude
1589	integer,save :: ilat_inf=-1,ilat_sup=-1 ! along latitude
1590	integer,save :: ialt_inf=-1,ialt_sup=-1 ! along altitude
1591	integer,save :: itim_inf=-1,itim_sup=-1 ! along time
1592
1593	! intermediate interpolated values
1594	real :: t_interp(2,2,2) ! after time interpolation
1595	real :: zt_interp(2,2) ! after altitude interpolation
1596	real :: yzt_interp(2) ! after latitude interpolation
1597	real :: coeff ! interpolation coefficient
1598
1599	integer :: i
1600
1601	! By default, set value to missing_value
1602	value=missing_value
1603
1604	!================================================================
1605	! 1. Find encompassing indexes
1606	!================================================================
1607	if (lon.ne.prev_lon) then
1608	do i=1,lonlen-1
1609	if (longitude(i).le.lon) then
1610	ilon_inf=i
1611	else
1612	exit
1613	endif
1614	enddo
1615	ilon_sup=ilon_inf+1
1616	endif ! of if (lon.ne.prev_lon)
1617	!write(,) 'ilon_inf=',ilon_inf," longitude(ilon_inf)=",longitude(ilon_inf)
1618
1619	if (lat.ne.prev_lat) then
1620	! recall that latitudes start from north pole
1621	do i=1,latlen-1
1622	if (latitude(i).ge.lat) then
1623	ilat_inf=i
1624	else
1625	exit
1626	endif
1627	enddo
1628	ilat_sup=ilat_inf+1
1629	endif ! of if (lat.ne.prev_lat)
1630	!write(,) 'ilat_inf=',ilat_inf," latitude(ilat_inf)=",latitude(ilat_inf)
1631
1632	if (alt.ne.prev_alt) then
1633	if (alttype.eq.'p') then ! pressures are ordered from max to min
1634	!handle special case for alt not in the altitude(1:altlen) interval
1635	if ((alt.gt.altitude(1)).or.(alt.lt.altitude(altlen))) then
1636	ialt_inf=-1
1637	ialt_sup=-1
1638	! return to main program (with value=missing_value)
1639	! write(,)"Problem in extraction : GCM alt p"
1640	return
1641	else ! general case
1642	do i=1,altlen-1
1643	if (altitude(i).ge.alt) then
1644	ialt_inf=i
1645	else
1646	exit
1647	endif
1648	enddo
1649	ialt_sup=ialt_inf+1
1650	endif ! of if ((alt.gt.altitude(1)).or.(alt.lt.altitude(altlen)))
1651	else ! altitudes (m) are ordered from min to max
1652	!handle special case for alt not in the altitude(1:altlen) interval
1653	if ((alt.lt.altitude(1)).or.(alt.gt.altitude(altlen))) then
1654	ialt_inf=-1
1655	ialt_sup=-1
1656	! return to main program (with value=missing_value)
1657	! write(,)"Problem in extraction : GCM alt z"
1658	return
1659	else ! general case
1660	do i=1,altlen-1
1661	if (altitude(i).le.alt) then
1662	ialt_inf=i
1663	else
1664	exit
1665	endif
1666	enddo
1667	ialt_sup=ialt_inf+1
1668	endif ! of if ((alt.lt.altitude(1)).or.(alt.gt.altitude(altlen)))
1669	endif ! of if (alttype.eq.'p')
1670	endif ! of if (alt.ne.prev_alt)
1671	!write(,) 'ialt_inf=',ialt_inf," altitude(ialt_inf)=",altitude(ialt_inf)
1672
1673	if (sol.ne.prev_sol) then
1674	!handle special case for sol not in the time(1):time(timenlen) interval
1675	if ((sol.lt.time(1)).or.(sol.gt.time(timelen))) then
1676	itim_inf=-1
1677	itim_sup=-1
1678	! return to main program (with value=missing_value)
1679	! write(,)"Problem in extraction : GCM sol"
1680	return
1681	else ! general case
1682	do i=1,timelen-1
1683	if (time(i).le.sol) then
1684	itim_inf=i
1685	else
1686	exit
1687	endif
1688	enddo
1689	itim_sup=itim_inf+1
1690	endif ! of if ((sol.lt.time(1)).or.(sol.gt.time(timenlen)))
1691	endif ! of if (sol.ne.prev_sol)
1692	!write(,) 'itim_inf=',itim_inf," time(itim_inf)=",time(itim_inf)
1693	!write(,) 'itim_sup=',itim_sup," time(itim_sup)=",time(itim_sup)
1694
1695	!================================================================
1696	! 2. Interpolate
1697	!================================================================
1698	! check that there are no "missing_value" in the field() elements we need
1699	! otherwise return to main program (with value=missing_value)
1700	if (field(ilon_inf,ilat_inf,ialt_inf,itim_inf).eq.missing_value) then
1701	! write(,)"Error 1 in extraction"
1702	return
1703	endif
1704	if (field(ilon_inf,ilat_inf,ialt_inf,itim_sup).eq.missing_value) then
1705	! write(,)"Error 2 in extraction"
1706	return
1707	endif
1708	if (field(ilon_inf,ilat_inf,ialt_sup,itim_inf).eq.missing_value) then
1709	! write(,)"Error 3 in extraction"
1710	return
1711	endif
1712	if (field(ilon_inf,ilat_inf,ialt_sup,itim_sup).eq.missing_value) then
1713	! write(,)"Error 4 in extraction"
1714	return
1715	endif
1716	if (field(ilon_inf,ilat_sup,ialt_inf,itim_inf).eq.missing_value) then
1717	! write(,)"Error 5 in extraction"
1718	return
1719	endif
1720	if (field(ilon_inf,ilat_sup,ialt_inf,itim_sup).eq.missing_value) then
1721	! write(,)"Error 6 in extraction"
1722	return
1723	endif
1724	if (field(ilon_inf,ilat_sup,ialt_sup,itim_inf).eq.missing_value) then
1725	! write(,)"Error 7 in extraction"
1726	return
1727	endif
1728	if (field(ilon_inf,ilat_sup,ialt_sup,itim_sup).eq.missing_value) then
1729	! write(,)"Error 8 in extraction"
1730	return
1731	endif
1732	if (field(ilon_sup,ilat_inf,ialt_inf,itim_inf).eq.missing_value) then
1733	! write(,)"Error 9 in extraction"
1734	return
1735	endif
1736	if (field(ilon_sup,ilat_inf,ialt_inf,itim_sup).eq.missing_value) then
1737	! write(,)"Error 10 in extraction"
1738	return
1739	endif
1740	if (field(ilon_sup,ilat_inf,ialt_sup,itim_inf).eq.missing_value) then
1741	! write(,)"Error 11 in extraction"
1742	return
1743	endif
1744	if (field(ilon_sup,ilat_inf,ialt_sup,itim_sup).eq.missing_value) then
1745	! write(,)"Error 12 in extraction"
1746	return
1747	endif
1748	if (field(ilon_sup,ilat_sup,ialt_inf,itim_inf).eq.missing_value) then
1749	! write(,)"Error 13 in extraction"
1750	return
1751	endif
1752	if (field(ilon_sup,ilat_sup,ialt_inf,itim_sup).eq.missing_value) then
1753	! write(,)"Error 14 in extraction"
1754	return
1755	endif
1756	if (field(ilon_sup,ilat_sup,ialt_sup,itim_inf).eq.missing_value) then
1757	! write(,)"Error 15 in extraction"
1758	return
1759	endif
1760	if (field(ilon_sup,ilat_sup,ialt_sup,itim_sup).eq.missing_value) then
1761	! write(,)"Error 16 in extraction"
1762	return
1763	endif
1764
1765	!================================================================
1766	! 2.1 Linear interpolation in time
1767	!================================================================
1768	coeff=(sol-time(itim_inf))/(time(itim_sup)-time(itim_inf))
1769	t_interp(1,1,1)=field(ilon_inf,ilat_inf,ialt_inf,itim_inf)+ &
1770	coeff*(field(ilon_inf,ilat_inf,ialt_inf,itim_sup)- &
1771	field(ilon_inf,ilat_inf,ialt_inf,itim_inf))
1772	t_interp(1,1,2)=field(ilon_inf,ilat_inf,ialt_sup,itim_inf)+ &
1773	coeff*(field(ilon_inf,ilat_inf,ialt_sup,itim_sup)- &
1774	field(ilon_inf,ilat_inf,ialt_sup,itim_inf))
1775	t_interp(1,2,1)=field(ilon_inf,ilat_sup,ialt_inf,itim_inf)+ &
1776	coeff*(field(ilon_inf,ilat_sup,ialt_inf,itim_sup)- &
1777	field(ilon_inf,ilat_sup,ialt_inf,itim_inf))
1778	t_interp(1,2,2)=field(ilon_inf,ilat_sup,ialt_sup,itim_inf)+ &
1779	coeff*(field(ilon_inf,ilat_sup,ialt_sup,itim_sup)- &
1780	field(ilon_inf,ilat_sup,ialt_sup,itim_inf))
1781	t_interp(2,1,1)=field(ilon_sup,ilat_inf,ialt_inf,itim_inf)+ &
1782	coeff*(field(ilon_sup,ilat_inf,ialt_inf,itim_sup)- &
1783	field(ilon_sup,ilat_inf,ialt_inf,itim_inf))
1784	t_interp(2,1,2)=field(ilon_sup,ilat_inf,ialt_sup,itim_inf)+ &
1785	coeff*(field(ilon_sup,ilat_inf,ialt_sup,itim_sup)- &
1786	field(ilon_sup,ilat_inf,ialt_sup,itim_inf))
1787	t_interp(2,2,1)=field(ilon_sup,ilat_sup,ialt_inf,itim_inf)+ &
1788	coeff*(field(ilon_sup,ilat_sup,ialt_inf,itim_sup)- &
1789	field(ilon_sup,ilat_sup,ialt_inf,itim_inf))
1790	t_interp(2,2,2)=field(ilon_sup,ilat_sup,ialt_sup,itim_inf)+ &
1791	coeff*(field(ilon_sup,ilat_sup,ialt_sup,itim_sup)- &
1792	field(ilon_sup,ilat_sup,ialt_sup,itim_inf))
1793
1794	!================================================================
1795	! 2.2 Vertical interpolation
1796	!================================================================
1797	if (((varname=='rho').or.(varname=='pressure')).and.(alttype=='z')) then
1798	! do the interpolation on the log of the quantity
1799	coeff=(alt-altitude(ialt_inf))/(altitude(ialt_sup)-altitude(ialt_inf))
1800	zt_interp(1,1)=log(t_interp(1,1,1))+coeff* &
1801	(log(t_interp(1,1,2))-log(t_interp(1,1,1)))
1802	zt_interp(1,2)=log(t_interp(1,2,1))+coeff* &
1803	(log(t_interp(1,2,2))-log(t_interp(1,2,1)))
1804	zt_interp(2,1)=log(t_interp(2,1,1))+coeff* &
1805	(log(t_interp(2,1,2))-log(t_interp(2,1,1)))
1806	zt_interp(2,2)=log(t_interp(2,2,1))+coeff* &
1807	(log(t_interp(2,2,2))-log(t_interp(2,2,1)))
1808	zt_interp(1:2,1:2)=exp(zt_interp(1:2,1:2))
1809	else ! general case
1810	coeff=(alt-altitude(ialt_inf))/(altitude(ialt_sup)-altitude(ialt_inf))
1811	zt_interp(1,1)=t_interp(1,1,1)+coeff*(t_interp(1,1,2)-t_interp(1,1,1))
1812	zt_interp(1,2)=t_interp(1,2,1)+coeff*(t_interp(1,2,2)-t_interp(1,2,1))
1813	zt_interp(2,1)=t_interp(2,1,1)+coeff*(t_interp(2,1,2)-t_interp(2,1,1))
1814	zt_interp(2,2)=t_interp(2,2,1)+coeff*(t_interp(2,2,2)-t_interp(2,2,1))
1815	endif
1816
1817	!================================================================
1818	! 2.3 Latitudinal interpolation
1819	!================================================================
1820	coeff=(lat-latitude(ilat_inf))/(latitude(ilat_sup)-latitude(ilat_inf))
1821	yzt_interp(1)=zt_interp(1,1)+coeff*(zt_interp(1,2)-zt_interp(1,1))
1822	yzt_interp(2)=zt_interp(2,1)+coeff*(zt_interp(2,2)-zt_interp(2,1))
1823
1824	!================================================================
1825	! 2.4 longitudinal interpolation
1826	!================================================================
1827	coeff=(lon-longitude(ilon_inf))/(longitude(ilon_sup)-longitude(ilon_inf))
1828	value=yzt_interp(1)+coeff*(yzt_interp(2)-yzt_interp(1))
1829
1830	end subroutine extraction
1831
1832	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1833
1834	subroutine inidim(outfid,lonlen,latlen,altlen,timelen,lon,lat,alt,time,units_alt,&
1835	londimid,latdimid,altdimid,timedimid)
1836	!================================================================
1837	! Purpose:
1838	! Initialize a data file (NetCDF format)
1839	!================================================================
1840	! Remarks:
1841	! The NetCDF file remains open
1842	!================================================================
1843	use netcdf ! NetCDF definitions
1844	implicit none
1845	!================================================================
1846	! Arguments:
1847	!================================================================
1848	integer, intent(in):: outfid
1849	! outfid: [netcdf] file ID
1850	integer, intent(in):: lonlen
1851	! lonlen: longitude length (# of longitude values)
1852	integer, intent(in):: latlen
1853	! latlen: latitude length (# of latitude values)
1854	integer, intent(in):: altlen
1855	! altlen: altitude length (# of altitude values)
1856	integer, intent(in):: timelen
1857	! timelen: time length (# of time values)
1858	real, intent(in):: lon(lonlen)
1859	! lon(): longitude
1860	real, intent(in):: lat(latlen)
1861	! lat(): latitude
1862	real, intent(in):: alt(altlen)
1863	! alt(): altitude
1864	real, intent(in):: time(timelen)
1865	! time(): time (Ls)
1866	character(len=1), intent(in) :: units_alt
1867	! units_alt: altitude coord. type:'z' (altitude, m) 'p' (pressure, Pa)
1868	integer,intent(inout):: londimid
1869	! londimid: [netcdf] lon() (i.e.: longitude) ID in MCS and output files (they are the same)
1870	integer,intent(inout) :: latdimid
1871	! latdimid: [netcdf] lat() ID in MCS and output files (they are the same)
1872	integer,intent(inout):: altdimid
1873	! altdimid: [netcdf] alt() ID in MCS and output files (they are the same)
1874	integer,intent(inout):: timedimid
1875	! timedimid: [netcdf] time() ID in MCS and output files (they are the same)
1876
1877	!================================================================
1878	! Local variables:
1879	!================================================================
1880	integer :: lonvarid,latvarid,altvarid,timevarid,status
1881	! *varid: [netcdf] ID of a variable
1882	! status: [netcdf] return error code (from called subroutines)
1883	character(len=256) :: error_text
1884	integer :: d ! loop index on dimensions ID
1885
1886	!===============================================================
1887	! 1. Define/write "dimensions" in the same order than MCS file
1888	! and get their IDs
1889	!================================================================
1890	do d=1,4
1891	if (altdimid.eq.d) then
1892	status=nf90_def_dim(outfid, "altitude", altlen, altdimid)
1893	error_text="Error: could not define the altitude dimension in the outfile"
1894	call status_check(status,error_text)
1895	else if (timedimid.eq.d) then
1896	status=nf90_def_dim(outfid, "time", timelen, timedimid)
1897	error_text="Error: could not define the time dimension in the outfile"
1898	call status_check(status,error_text)
1899	else if (latdimid.eq.d) then
1900	status=nf90_def_dim(outfid, "latitude", latlen, latdimid)
1901	error_text="Error: could not define the latitude dimension in the outfile"
1902	call status_check(status,error_text)
1903	else if (londimid.eq.d) then
1904	status=nf90_def_dim(outfid, "longitude", lonlen, londimid)
1905	error_text="Error: could not define the longitude dimension in the outfile"
1906	call status_check(status,error_text)
1907	endif
1908	enddo
1909
1910	!================================================================
1911	! 2. Define "variables" and their attributes
1912	!================================================================
1913	!================================================================
1914	! 2.1 Write "longitude" (data and attributes)
1915	!================================================================
1916
1917	! Insert the definition of the variable
1918	status=nf90_def_var(outfid,"longitude",nf90_float,(/londimid/),lonvarid)
1919	error_text="Error: could not define the longitude variable in the outfile"
1920	call status_check(status,error_text)
1921
1922	! Write the attributes
1923	status=nf90_put_att(outfid,lonvarid,"long_name","longitude")
1924	error_text="Error: could not put the long_name attribute for the longitude variable in the outfile"
1925	call status_check(status,error_text)
1926
1927	status=nf90_put_att(outfid,lonvarid,"units","degrees_east")
1928	error_text="Error: could not put the units attribute for the longitude variable in the outfile"
1929	call status_check(status,error_text)
1930
1931	!================================================================
1932	! 2.2 "latitude"
1933	!================================================================
1934
1935	! Insert the definition of the variable
1936	status=nf90_def_var(outfid,"latitude",nf90_float,(/latdimid/),latvarid)
1937	error_text="Error: could not define the latitude variable in the outfile"
1938	call status_check(status,error_text)
1939
1940	! Write the attributes
1941	status=nf90_put_att(outfid,latvarid,"long_name","latitude")
1942	error_text="Error: could not put the long_name attribute for the latitude variable in the outfile"
1943	call status_check(status,error_text)
1944
1945	status=nf90_put_att(outfid,latvarid,"units","degrees_north")
1946	error_text="Error: could not put the units attribute for the latitude variable in the outfile"
1947	call status_check(status,error_text)
1948
1949	!================================================================
1950	! 2.3 Write "altitude" (data and attributes)
1951	!================================================================
1952
1953	! Insert the definition of the variable
1954	status=nf90_def_var(outfid,"altitude",nf90_float,(/altdimid/),altvarid)
1955	error_text="Error: could not define the altitude variable in the outfile"
1956	call status_check(status,error_text)
1957
1958	! Write the attributes
1959	if (units_alt.eq.'p') then ! pressure coordinate
1960	status=nf90_put_att(outfid,altvarid,"long_name","pressure")
1961	error_text="Error: could not put the long_name attribute for the altitude variable in the outfile"
1962	call status_check(status,error_text)
1963
1964	status=nf90_put_att(outfid,altvarid,'units',"Pa")
1965	error_text="Error: could not put the units attribute for the altitude variable in the outfile"
1966	call status_check(status,error_text)
1967
1968	status=nf90_put_att(outfid,altvarid,'positive',"down")
1969	error_text="Error: could not put the positive attribute for the altitude variable in the outfile"
1970	call status_check(status,error_text)
1971
1972	status=nf90_put_att(outfid,altvarid,'comment',&
1973	"The vertical variable is in fact pressure, not altitude. We just call it altitude for easy reading in Ferret and Grads.")
1974	error_text="Error: could not put the comment attribute for the altitude variable in the outfile"
1975	call status_check(status,error_text)
1976
1977	else if (units_alt.eq.'z') then ! altitude coordinate
1978	status=nf90_put_att(outfid,altvarid,"long_name","altitude")
1979	error_text="Error: could not put the long_name attribute for the altitude variable in the outfile"
1980	call status_check(status,error_text)
1981
1982	status=nf90_put_att(outfid,altvarid,'units',"m")
1983	error_text="Error: could not put the units attribute for the altitude variable in the outfile"
1984	call status_check(status,error_text)
1985
1986	status=nf90_put_att(outfid,altvarid,'positive',"up")
1987	error_text="Error: could not put the positive attribute for the altitude variable in the outfile"
1988	call status_check(status,error_text)
1989
1990	else
1991	write(,)"I do not understand this unit type ",units_alt," for outfile altitude!"
1992	stop
1993	end if
1994
1995	!================================================================
1996	! 2.4 "time"
1997	!================================================================
1998
1999	! Insert the definition of the variable
2000	status=nf90_def_var(outfid,"time",nf90_float,(/timedimid/),timevarid)
2001	error_text="Error: could not define the time variable in the outfile"
2002	call status_check(status,error_text)
2003
2004	! Write the attributes
2005	status=nf90_put_att(outfid,timevarid,"long_name","Solar longitude")
2006	error_text="Error: could not put the long_name attribute for the time variable in the outfile"
2007	call status_check(status,error_text)
2008
2009	status=nf90_put_att(outfid,timevarid,"units","days since 0000-01-1 00:00:00")
2010	error_text="Error: could not put the units attribute for the time variable in the outfile"
2011	call status_check(status,error_text)
2012
2013	status=nf90_put_att(outfid,timevarid,"comment",&
2014	"Units is in fact degrees, but set to a dummy value of days for compatibility with Ferret and Grads.")
2015	error_text="Error: could not put the comment attribute for the time variable in the outfile"
2016	call status_check(status,error_text)
2017
2018	!================================================================
2019	! 2.5 End netcdf define mode
2020	!================================================================
2021	status=nf90_enddef(outfid)
2022	error_text="Error: could not end the define mode of the outfile"
2023	call status_check(status,error_text)
2024
2025	!================================================================
2026	! 3. Write "variables" (data of the dimension variables)
2027	!================================================================
2028	! "time"
2029	status=nf90_put_var(outfid,timevarid,time)
2030	error_text="Error: could not write the time variable's data in the outfile"
2031	call status_check(status,error_text)
2032
2033	! "latitude"
2034	status=nf90_put_var(outfid,latvarid,lat)
2035	error_text="Error: could not write the latitude variable's data in the outfile"
2036	call status_check(status,error_text)
2037
2038	! "longitude"
2039	status=nf90_put_var(outfid,lonvarid,lon)
2040	error_text="Error: could not write the longitude variable's data in the outfile"
2041	call status_check(status,error_text)
2042
2043	! "altitude"
2044	status=nf90_put_var(outfid,altvarid,alt)
2045	error_text="Error: could not write the altitude variable's data in the outfile"
2046	call status_check(status,error_text)
2047
2048	end subroutine inidim
2049
2050	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2051
2052	subroutine ls2sol(ls,sol)
2053
2054	implicit none
2055	!================================================================
2056	! Arguments:
2057	!================================================================
2058	real,intent(in) :: ls
2059	real,intent(out) :: sol
2060
2061	!================================================================
2062	! Local:
2063	!================================================================
2064	double precision xref,zx0,zteta,zz
2065	!xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
2066	double precision year_day
2067	double precision peri_day,timeperi,e_elips
2068	double precision pi,degrad
2069	parameter (year_day=668.6d0) ! number of sols in a martian year
2070	parameter (peri_day=485.35d0) ! date (in sols) of perihelion
2071	!timeperi: 2pi( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
2072	parameter (timeperi=1.90258341759902d0)
2073	parameter (e_elips=0.0934d0) ! eccentricity of orbit
2074	parameter (pi=3.14159265358979d0)
2075	parameter (degrad=57.2957795130823d0)
2076
2077	if (abs(ls).lt.1.0e-5) then
2078	if (ls.ge.0.0) then
2079	sol = 0.0
2080	else
2081	sol = real(year_day)
2082	end if
2083	return
2084	end if
2085
2086	zteta = ls/degrad + timeperi
2087	zx0 = 2.0datan(dtan(0.5zteta)/dsqrt((1.+e_elips)/(1.-e_elips)))
2088	xref = zx0-e_elips*dsin(zx0)
2089	zz = xref/(2.*pi)
2090	sol = real(zz*year_day + peri_day)
2091	if (sol.lt.0.0) sol = sol + real(year_day)
2092	if (sol.ge.year_day) sol = sol - real(year_day)
2093
2094
2095	end subroutine ls2sol
2096
2097	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2098
2099	subroutine gen_sol_list(solcount,int_sol_list,LTcount,LTave,LTmax,LTmin,lon,f_LT,dayornight,&
2100	sol_list)
2101	!================================================================
2102	! Purpose:
2103	! Generate a list that is the combination of two lists :
2104	! - int_sol_list, which is a list of integer values of sol
2105	! - LT_list, which contains a number LTcount of LT values in the
2106	! interval [LTmin;LTmax] which are evenly distributed around
2107	! LTave (so that their average is LTave)
2108	!================================================================
2109
2110	implicit none
2111	!================================================================
2112	! Arguments:
2113	!================================================================
2114	integer, intent(in) :: solcount ! nb of integer values of sol
2115	real, intent(in) :: int_sol_list(solcount) ! list of the integer values of sol
2116	integer, intent(in) :: LTcount ! nb of LT per sol to be interpolated
2117	real, intent(in) :: LTave ! average of LT
2118	real, intent(in) :: LTmax, LTmin ! bounds of the LT interval for the interpolation
2119	real, intent(in) :: lon ! longitude value for the transformation into a sol value at lon=0°
2120	external f_LT ! function that changes LT interval for night LT
2121	real f_LT
2122	character (len=10), intent(in) :: dayornight ! to know if we have day or night values
2123
2124	real, intent(out) :: sol_list(solcount*LTcount) ! all the sol values at lon=0° to interpolate
2125
2126	!================================================================
2127	! Local variables:
2128	!================================================================
2129	integer :: N
2130	integer :: a,b,c ! loop iteration indexes
2131	real :: LT_list(LTcount)
2132
2133	N = floor(LTcount/2.)
2134
2135	!================================================================
2136	! 1. Filling of LT_list
2137	!================================================================
2138	SELECT CASE (dayornight)
2139	CASE ("dayside")
2140	if (abs(LTave-LTmax).le.abs(LTave-LTmin)) then ! LTave is closer to LTmax than to LTmin
2141	do a=1,N
2142	LT_list(a)=LTave+a*abs(LTave-LTmax)/(N+1)
2143	LT_list(N+a)=LTave-a*abs(LTave-LTmax)/(N+1)
2144	enddo
2145	else ! LTave is closer to LTmin than to LTmax
2146	do a=1,N
2147	LT_list(a)=LTave+a*abs(LTave-LTmin)/(N+1)
2148	LT_list(N+a)=LTave-a*abs(LTave-LTmin)/(N+1)
2149	enddo
2150	endif
2151	CASE ("nightside")
2152	if (abs(f_LT(LTave)-f_LT(LTmax)).le.abs(f_LT(LTave)-f_LT(LTmin))) then ! LTave is closer to LTmax than to LTmin
2153	do a=1,N
2154	LT_list(a)=LTave+a*abs(f_LT(LTave)-f_LT(LTmax))/(N+1)
2155	LT_list(a)=mod(LT_list(a),24.) ! reput the LT in a [0;24[ interval if needed
2156	LT_list(N+a)=LTave-a*abs(f_LT(LTave)-f_LT(LTmax))/(N+1)
2157	LT_list(N+a)=mod(LT_list(N+a),24.)
2158	enddo
2159	else ! LTave is closer to LTmin than to LTmax
2160	do a=1,N
2161	LT_list(a)=LTave+a*abs(f_LT(LTave)-f_LT(LTmin))/(N+1)
2162	LT_list(a)=mod(LT_list(a),24.)
2163	LT_list(N+a)=LTave-a*abs(f_LT(LTave)-f_LT(LTmin))/(N+1)
2164	LT_list(N+a)=mod(LT_list(N+a),24.)
2165	enddo
2166	endif
2167	END SELECT
2168
2169	if (mod(LTcount,2).eq.1) then ! if LTcount is an odd number
2170	LT_list(LTcount)=LTave ! add LTave to the list
2171	endif
2172
2173	!================================================================
2174	! 2. Combination of int_sol_list & LT_list into sol_list
2175	!================================================================
2176	c=1
2177	do a=1,solcount
2178	do b=1,LTcount
2179	sol_list(c)=int_sol_list(a)+(LT_list(b)-lon/15.)/24.
2180	c=c+1
2181	enddo
2182	enddo
2183
2184	end subroutine gen_sol_list
2185
2186	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2187
2188	subroutine status_check(status,error_text)
2189
2190	use netcdf
2191	implicit none
2192	!================================================================
2193	! Arguments:
2194	!================================================================
2195	integer,intent(in) :: status
2196	character(len=256),intent(in) :: error_text
2197
2198	if (status.ne.nf90_noerr) then
2199	write(,)trim(error_text)
2200	write(,)trim(nf90_strerror(status))
2201	stop
2202	endif
2203
2204	end subroutine status_check
2205
2206	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2207
2208	real function LTmod(LT)
2209	!================================================================
2210	! Purpose:
2211	! For night local times management, replace hours
2212	! from a [0;24[ interval to a [-12;12[ interval in which
2213	! midnight = 0 (in order to ensure continuity when comparing
2214	! night local times)
2215	!================================================================
2216	implicit none
2217	real,intent(in) :: LT
2218
2219	LTmod = 2*mod(LT,12.)-LT
2220	return
2221	end function LTmod

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

Download in other formats:

Original Format