Context Navigation

simu_MCS.F90 @ 2540

Last change on this file since 2540 was 2540, checked in by abierjon, 3 years ago

Mars GCM:
Make simu_MCS.F90 compatible with DYNAMICO lon-lat output files, with a check on the latitude array order when doing the interpolation

File size: 89.6 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 ! 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	! Initialize logicals
542	dustok1 = .false.
543	dustok2 = .false.
544	dustok3 = .false.
545	wiceok = .false.
546
547	! Get nbvarfile (total number of variables in the GCM file)
548	status=NF90_INQUIRE(gcmfid,nVariables=nbvarfile)
549	error_text="Error : Pb with nf90_inquire(gcmfid,nVariables=nbvarfile)"
550	call status_check(status,error_text)
551
552	! List of variables that should not be processed
553	notprocessed(1)='Time'
554	notprocessed(2)='controle'
555	notprocessed(3)='rlonu'
556	notprocessed(4)='latitude'
557	notprocessed(5)='longitude'
558	notprocessed(6)='altitude'
559	notprocessed(7)='rlatv'
560	notprocessed(8)='aps'
561	notprocessed(9)='bps'
562	notprocessed(10)='ap'
563	notprocessed(11)='bp'
564	notprocessed(12)='cu'
565	notprocessed(13)='cv'
566	notprocessed(14)='aire'
567	notprocessed(15)='phisinit'
568
569
570	! List of variables in the GCM file
571	write(,)""
572	write(,)"List of variables in the GCM file :"
573	Nnotprocessed=0
574	do v=1,nbvarfile
575	status=NF90_INQUIRE_VARIABLE(gcmfid,v,name=GCMvarname)
576	! GCMvarname now contains the "name" of variable of ID # v
577	var_ok=.true.
578	do vnot=1,15
579	if (GCMvarname.eq.notprocessed(vnot)) then
580	var_ok=.false.
581	Nnotprocessed=Nnotprocessed+1
582	endif
583	enddo
584	if (var_ok) write(,) trim(GCMvarname)
585
586	! Detect if we can compute dust and wice opacities
587	if (trim(GCMvarname).eq."dso") then
588	dustok1 = .true.
589	else if (trim(GCMvarname).eq."dsodust") then
590	dustok2 = .true.
591	else if (trim(GCMvarname).eq."dustq") then
592	dustok3 = .true.
593	else if (trim(GCMvarname).eq."h2o_ice") then
594	wiceok = .true.
595	endif
596	enddo
597
598	! Nnotprocessed: # of variables that won't be processed
599	! nbvarfile: total # of variables in file
600	! +2: the dust and wice opacities
601	allocate(gcm_vars(nbvarfile-Nnotprocessed+2),stat=status)
602	if (status.ne.0) then
603	write(,) "Error: failed allocation of gcm_vars(nbvarfile-Nnotprocessed+2)"
604	write(,) " nbvarfile=",nbvarfile
605	write(,) " Nnotprocessed=",Nnotprocessed
606	stop
607	endif
608
609	! List of variables to process
610	write(,)
611	write(,) "Which variables do you want to redistribute ?"
612	write(,) "list of <variables> (separated by <Return>s)"
613	write(,) "(an empty line , i.e: just <Return>, implies end of list)"
614	write(,) "NB: this program handles only 4D netcdf variables for now"
615	nbvar=0
616	read(*,'(a50)') GCMvarname
617	do while ((GCMvarname.ne.' ').AND.(trim(GCMvarname).ne."all"))
618	nbvar=nbvar+1
619	gcm_vars(nbvar)=GCMvarname
620	read(*,'(a50)') GCMvarname
621	enddo
622
623	if (GCMvarname.eq."all") then
624	nbvar=nbvarfile-Nnotprocessed
625	do v=Nnotprocessed+1,nbvarfile
626	status=nf90_inquire_variable(gcmfid,v,name=gcm_vars(v-Nnotprocessed))
627	enddo
628	! Variables names from the file are stored in gcm_vars()
629	nbvar=nbvarfile-Nnotprocessed
630	do v=1,nbvar
631	status=nf90_inquire_variable(gcmfid,v+Nnotprocessed,name=gcm_vars(v))
632	write(*,'(a9,1x,i2,1x,a1,1x,a64)') "variable ",v,":",gcm_vars(v)
633	enddo
634	else if(nbvar==0) then
635	write(,) "No variables to process in the GCM file... program stopped"
636	stop
637	endif ! of if (GCMvarname.eq."all")
638
639	!================================================================
640	! 2.1.2 Handle dust and wice opacities (first set-ups)
641	!================================================================
642	! 2nd part is in section 2.4.2
643	write(,)
644	! Load atmospheric density "rho"
645	if (dustok1.or.dustok2.or.dustok3.or.wiceok) then
646	! Check that the GCM file contains that variable
647	status=nf90_inq_varid(gcmfid,"rho",GCMvarid)
648	if (status.ne.nf90_noerr) then
649	write(,) "Failed to find variable rho in "//trim(gcmfile)
650	write(,) "No computation of opacities..."
651	dustok1 =.false.
652	dustok2 =.false.
653	dustok3 =.false.
654	wiceok =.false.
655	else
656	! Length allocation for each dimension of the 4D variable
657	allocate(rho(GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen))
658
659	! Load datasets
660	if (.not.is_stats) then
661	status=NF90_GET_VAR(gcmfid,GCMvarid,rho)
662	error_text="Failed to load rho"
663	call status_check(status,error_text)
664	else
665	! if it is a stats file, we load only the first sol, and then copy it to all the other sols
666	status=NF90_GET_VAR(gcmfid,GCMvarid,rho(:,:,:,1:GCMstatstimelen))
667	error_text="Failed to load rho"
668	call status_check(status,error_text)
669	! write(,)"GCMstatstimelen = ", GCMstatstimelen
670	! write(,)"GCMtimelen = ", GCMtimelen
671	do l=(GCMstatstimelen+1),GCMtimelen
672	if (modulo(l,GCMstatstimelen).ne.0) then
673	rho(:,:,:,l) = rho(:,:,:,modulo(l,GCMstatstimelen))
674	else ! if l is a multiple of GCMstatstimelen, since the index modulo(l,GCMstatstimelen)=0
675	! doesn't exist, we make a special case
676	rho(:,:,:,l) = rho(:,:,:,GCMstatstimelen)
677	endif
678	enddo
679	endif
680	write(,) "Variable rho loaded from the GCM file"
681	endif
682	endif ! dustok1.or.dustok2.or.dustok3.or.wiceok
683
684	! Dust and wice opacity booleans
685	if (dustok1.or.dustok2.or.dustok3) then
686	nbvar=nbvar+1
687	gcm_vars(nbvar)="dust"
688	endif
689
690	if (wiceok) then
691	nbvar=nbvar+1
692	gcm_vars(nbvar)="wice"
693	endif
694
695	!write(,) "gcm_vars retrieved : ",gcm_vars(1:nbvar)
696
697	!===============================================================================
698	! 2.2 Definition of LT variables from obsfile to outfile
699	!===============================================================================
700	! --> the day/night loop begins here
701
702	!****************** NOTA BENE (cf sections 2.2 and 4)***********************
703	! We execute the program a first time with the daytime values, and then a second
704	! time with the nighttime values.
705	!*******************************************************************************
706
707	dayornight = "dayside" ! we begin with daytime temperature
708	write(,)"" ; write(,) "Beginning the 1st loop, on daytime values"; write(,)""
709	DAY_OR_NIGHT: DO ! (the end of the loop is in section 4.)
710
711	SELECT CASE (dayornight)
712	CASE ("dayside")
713	OBSLTave_name = "dtimeave"
714	OBSLTmax_name = "dtimemax"
715	OBSLTmin_name = "dtimemin"
716	CASE ("nightside")
717	OBSLTave_name = "ntimeave"
718	OBSLTmax_name = "ntimemax"
719	OBSLTmin_name = "ntimemin"
720	END SELECT
721
722	!================================================================
723	! 2.2.1 Average of Local Time in the OBS bins
724	!================================================================
725	! Read the OBS file
726	!------------------
727	status=nf90_inq_varid(obsfid,trim(OBSLTave_name),LT_id)
728	error_text="Failed to find Observer local time ("//trim(OBSLTave_name)//") ID in "//trim(obsfile)
729	call status_check(status,error_text)
730	status=nf90_inquire_variable(obsfid,LT_id,dimids=LTshape)
731	error_text="Failed to get the dim shape of variable "//trim(OBSLTave_name)
732	call status_check(status,error_text)
733
734	! Length allocation for each dimension of the 3D variable
735	allocate(OBSLT(OBSlonlen,OBSlatlen,OBSLslen))
736
737	! Load datasets
738	status=NF90_GET_VAR(obsfid,LT_id,OBSLT)
739	error_text="Failed to load "//trim(OBSLTave_name)//" from the obsfile"
740	call status_check(status,error_text)
741	write(,) trim(OBSLTave_name)," loaded from the obsfile"
742
743	! Get LT missing_value attribute
744	status=nf90_get_att(obsfid,LT_id,"_FillValue",LTmiss_val)
745	error_text="Failed to load missing_value attribute"
746	call status_check(status,error_text)
747
748	! Create the variable in the outfile
749	!-----------------------------------
750	! Switch to netcdf define mode
751	status=nf90_redef(outfid)
752	error_text="Error: could not switch to define mode in the outfile"
753	call status_check(status,error_text)
754
755	! Definition of the variable
756	status=NF90_DEF_VAR(outfid,trim(OBSLTave_name),nf90_float,LTshape,LT_id)
757	error_text="Error: could not define the variable "//trim(OBSLTave_name)//" in the outfile"
758	call status_check(status,error_text)
759
760	! Write the attributes
761	select case (dayornight)
762	case ("dayside")
763	status=nf90_put_att(outfid,LT_id,"long_name","Average local time in bin - day side [6h, 18h]")
764	case ("nightside")
765	status=nf90_put_att(outfid,LT_id,"long_name","Average local time in bin - night side [18h, 6h]")
766	end select
767	status=nf90_put_att(outfid,LT_id,"units","hours")
768	status=nf90_put_att(outfid,LT_id,"_FillValue",LTmiss_val)
769
770	! End the netcdf define mode (and thus enter the "data writing" mode)
771	status=nf90_enddef(outfid)
772	error_text="Error: could not close the define mode of the outfile"
773	call status_check(status,error_text)
774
775	! Write the data in the output file
776	status = NF90_PUT_VAR(outfid, LT_id, OBSLT) ! write the MCS d/ntimeave as the output Local Time
777	error_text="Error: could not write "//trim(OBSLTave_name)//" data in the outfile"
778	call status_check(status,error_text)
779
780	write(,)"Local Time (",trim(OBSLTave_name),") has been created in the outfile"
781	write(*,'(" with missing_value attribute : ",1pe12.5)')LTmiss_val
782
783	!================================================================
784	! 2.2.2 Maximum of Local Time in the OBS bins
785	!================================================================
786	! Read the OBS file
787	!------------------
788	status=nf90_inq_varid(obsfid,trim(OBSLTmax_name),LT_id)
789	error_text="Failed to find Observer max local time ("//trim(OBSLTmax_name)//") ID in "//trim(obsfile)
790	call status_check(status,error_text)
791	status=nf90_inquire_variable(obsfid,LT_id,dimids=LTshape)
792	error_text="Failed to get the dim shape of variable "//trim(OBSLTmax_name)
793	call status_check(status,error_text)
794
795	! Length allocation for each dimension of the 3D variable
796	allocate(OBSLTmax(OBSlonlen,OBSlatlen,OBSLslen))
797
798	! Load datasets
799	status=NF90_GET_VAR(obsfid,LT_id,OBSLTmax)
800	error_text="Failed to load "//trim(OBSLTmax_name)//" from the obsfile"
801	call status_check(status,error_text)
802	write(,) trim(OBSLTmax_name)," loaded from the obsfile"
803
804	!================================================================
805	! 2.2.3 Minimum of Local Time in the OBS bins
806	!================================================================
807	! Read the OBS file
808	!------------------
809	status=nf90_inq_varid(obsfid,trim(OBSLTmin_name),LT_id)
810	error_text="Failed to find Observer min local time ("//trim(OBSLTmin_name)//") ID in "//trim(obsfile)
811	call status_check(status,error_text)
812	status=nf90_inquire_variable(obsfid,LT_id,dimids=LTshape)
813	error_text="Failed to obtain information on variable "//trim(OBSLTmin_name)
814	call status_check(status,error_text)
815
816	! Length allocation for each dimension of the 3D variable
817	allocate(OBSLTmin(OBSlonlen,OBSlatlen,OBSLslen))
818
819	! Load datasets
820	status=NF90_GET_VAR(obsfid,LT_id,OBSLTmin)
821	error_text="Failed to load "//trim(OBSLTmin_name)//" from the obsfile"
822	call status_check(status,error_text)
823	write(,) trim(OBSLTmin_name)," loaded from the obsfile"
824	write(,)""
825
826	!===============================================================================
827	! 2.3 Definition of numbin variables from obsfile to outfile
828	!===============================================================================
829	!================================================================
830	! 2.3.1 Number of values in the OBS "temp" bins
831	!================================================================
832	SELECT CASE (dayornight)
833	CASE ("dayside")
834	numbin_name = "dnumbintemp"
835	CASE ("nightside")
836	numbin_name = "nnumbintemp"
837	END SELECT
838
839	! Read the OBS file
840	!------------------
841	status=nf90_inq_varid(obsfid,trim(numbin_name),numbin_id)
842	error_text="Failed to find Observer nb of temp values in bin ("//trim(numbin_name)//")'s ID in "//trim(obsfile)
843	call status_check(status,error_text)
844	status=nf90_inquire_variable(obsfid,numbin_id,dimids=numbinshape)
845	error_text="Failed to obtain information on variable "//trim(numbin_name)
846	call status_check(status,error_text)
847
848	! Length allocation for each dimension of the 4D variable
849	allocate(numbin(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
850
851	! Load datasets
852	status=NF90_GET_VAR(obsfid,numbin_id,numbin)
853	error_text="Failed to load "//trim(numbin_name)//" from the obsfile"
854	call status_check(status,error_text)
855	write(,) trim(numbin_name)," loaded from the obsfile"
856
857	! Create the variable in the outfile
858	!-----------------------------------
859	! Switch to netcdf define mode
860	status=nf90_redef(outfid)
861	error_text="Error: could not switch to define mode in the outfile"
862	call status_check(status,error_text)
863
864	! Definition of the variable
865	status=NF90_DEF_VAR(outfid,trim(numbin_name),nf90_float,numbinshape,numbin_id)
866	error_text="Error: could not define the variable "//trim(numbin_name)//" in the outfile"
867	call status_check(status,error_text)
868
869	! Write the attributes
870	select case (dayornight)
871	case ("dayside")
872	status=nf90_put_att(outfid,numbin_id,"long_name","Number of temp values in bin - day side")
873	case ("nightside")
874	status=nf90_put_att(outfid,numbin_id,"long_name","Number of temp values in bin - night side")
875	end select
876
877	! End the netcdf define mode (and thus enter the "data writing" mode)
878	status=nf90_enddef(outfid)
879	error_text="Error: could not close the define mode of the outfile"
880	call status_check(status,error_text)
881
882	! Write the data in the output file
883	status = NF90_PUT_VAR(outfid, numbin_id, numbin)
884	error_text="Error: could not write "//trim(numbin_name)//" data in the outfile"
885	call status_check(status,error_text)
886
887	write(,)"Number of temp values in bin (",trim(numbin_name),") has been created in the outfile"
888	write(,)""
889	deallocate(numbin)
890
891	!================================================================
892	! 2.3.2 Number of values in the OBS "dust" bins
893	!================================================================
894	SELECT CASE (dayornight)
895	CASE ("dayside")
896	numbin_name = "dnumbindust"
897	CASE ("nightside")
898	numbin_name = "nnumbindust"
899	END SELECT
900
901	! Read the OBS file
902	!------------------
903	status=nf90_inq_varid(obsfid,trim(numbin_name),numbin_id)
904	error_text="Failed to find Observer nb of dust values in bin ("//trim(numbin_name)//")'s ID in "//trim(obsfile)
905	call status_check(status,error_text)
906	status=nf90_inquire_variable(obsfid,numbin_id,dimids=numbinshape)
907	error_text="Failed to obtain information on variable "//trim(numbin_name)
908	call status_check(status,error_text)
909
910	! Length allocation for each dimension of the 4D variable
911	allocate(numbin(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
912
913	! Load datasets
914	status=NF90_GET_VAR(obsfid,numbin_id,numbin)
915	error_text="Failed to load "//trim(numbin_name)//" from the obsfile"
916	call status_check(status,error_text)
917	write(,) trim(numbin_name)," loaded from the obsfile"
918
919	! Create the variable in the outfile
920	!-----------------------------------
921	! Switch to netcdf define mode
922	status=nf90_redef(outfid)
923	error_text="Error: could not switch to define mode in the outfile"
924	call status_check(status,error_text)
925
926	! Definition of the variable
927	status=NF90_DEF_VAR(outfid,trim(numbin_name),nf90_float,numbinshape,numbin_id)
928	error_text="Error: could not define the variable "//trim(numbin_name)//" in the outfile"
929	call status_check(status,error_text)
930
931	! Write the attributes
932	select case (dayornight)
933	case ("dayside")
934	status=nf90_put_att(outfid,numbin_id,"long_name","Number of dust values in bin - day side")
935	case ("nightside")
936	status=nf90_put_att(outfid,numbin_id,"long_name","Number of dust values in bin - night side")
937	end select
938
939	! End the netcdf define mode (and thus enter the "data writing" mode)
940	status=nf90_enddef(outfid)
941	error_text="Error: could not close the define mode of the outfile"
942	call status_check(status,error_text)
943
944	! Write the data in the output file
945	status = NF90_PUT_VAR(outfid, numbin_id, numbin)
946	error_text="Error: could not write "//trim(numbin_name)//" data in the outfile"
947	call status_check(status,error_text)
948
949	write(,)"Number of dust values in bin (",trim(numbin_name),") has been created in the outfile"
950	write(,)""
951	deallocate(numbin)
952
953	!================================================================
954	! 2.3.3 Number of values in the OBS "wice" bins
955	!================================================================
956	SELECT CASE (dayornight)
957	CASE ("dayside")
958	numbin_name = "dnumbinwice"
959	CASE ("nightside")
960	numbin_name = "nnumbinwice"
961	END SELECT
962
963	! Read the OBS file
964	!------------------
965	status=nf90_inq_varid(obsfid,trim(numbin_name),numbin_id)
966	error_text="Failed to find Observer nb of wice values in bin ("//trim(numbin_name)//")'s ID in "//trim(obsfile)
967	call status_check(status,error_text)
968	status=nf90_inquire_variable(obsfid,numbin_id,dimids=numbinshape)
969	error_text="Failed to obtain information on variable "//trim(numbin_name)
970	call status_check(status,error_text)
971
972	! Length allocation for each dimension of the 4D variable
973	allocate(numbin(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
974
975	! Load datasets
976	status=NF90_GET_VAR(obsfid,numbin_id,numbin)
977	error_text="Failed to load "//trim(numbin_name)//" from the obsfile"
978	call status_check(status,error_text)
979	write(,) trim(numbin_name)," loaded from the obsfile"
980
981	! Create the variable in the outfile
982	!-----------------------------------
983	! Switch to netcdf define mode
984	status=nf90_redef(outfid)
985	error_text="Error: could not switch to define mode in the outfile"
986	call status_check(status,error_text)
987
988	! Definition of the variable
989	status=NF90_DEF_VAR(outfid,trim(numbin_name),nf90_float,numbinshape,numbin_id)
990	error_text="Error: could not define the variable "//trim(numbin_name)//" in the outfile"
991	call status_check(status,error_text)
992
993	! Write the attributes
994	select case (dayornight)
995	case ("dayside")
996	status=nf90_put_att(outfid,numbin_id,"long_name","Number of wice values in bin - day side")
997	case ("nightside")
998	status=nf90_put_att(outfid,numbin_id,"long_name","Number of wice values in bin - night side")
999	end select
1000
1001	! End the netcdf define mode (and thus enter the "data writing" mode)
1002	status=nf90_enddef(outfid)
1003	error_text="Error: could not close the define mode of the outfile"
1004	call status_check(status,error_text)
1005
1006	! Write the data in the output file
1007	status = NF90_PUT_VAR(outfid, numbin_id, numbin)
1008	error_text="Error: could not write "//trim(numbin_name)//" data in the outfile"
1009	call status_check(status,error_text)
1010
1011	write(,)"Number of wice values in bin (",trim(numbin_name),") has been created in the outfile"
1012	write(,)""
1013	deallocate(numbin)
1014
1015
1016	!===============================================================================
1017	! 2.4 Opening of the GCM variables
1018	!===============================================================================
1019	! --> the var loop begins here
1020
1021	VAR: DO v=1,nbvar ! LOOP ON ALL THE GCM VARIABLES TO PROCESS
1022	! (the end of the loop is in section 3.4)
1023
1024	GCMvarname = gcm_vars(v)
1025
1026	! Detect the dust and wice opacities special cases
1027	if (trim(GCMvarname).eq."dust") then
1028	if (dustok1) then ! "dso" is detected in gcmfile
1029	GCMvarname="dso"
1030	dustok2=.false.
1031	dustok3=.false.
1032	else if (dustok2) then ! "dsodust" is detected in gcmfile
1033	GCMvarname="dsodust"
1034	dustok3=.false.
1035	else if (dustok3) then ! "dustq" is detected in gcmfile
1036	GCMvarname="dustq"
1037	endif
1038	write(,) "Computing dust opacity..."
1039	endif
1040	if (trim(GCMvarname).eq."wice") then ! "h2o_ice" detected in gcmfile
1041	GCMvarname="h2o_ice"
1042	write(,) "Computing water ice opacity..."
1043	endif
1044
1045	!================================================================
1046	! 2.4.1 Generic reading of the variable
1047	!================================================================
1048	! Check that the GCM file contains that variable
1049	status=nf90_inq_varid(gcmfid,trim(GCMvarname),GCMvarid)
1050	if (status.ne.nf90_noerr) then
1051	write(,) "Failed to find variable "//trim(GCMvarname)//" in "//trim(gcmfile)
1052	write(,) "We'll skip that variable..."
1053	CYCLE VAR ! go directly to the next variable
1054	endif
1055
1056	! Sanity checks on the variable
1057	status=nf90_inquire_variable(gcmfid,GCMvarid,ndims=nbdim,dimids=GCMvarshape)
1058	error_text="Failed to obtain information on variable "//trim(GCMvarname)
1059	call status_check(status,error_text)
1060
1061	! Check that it is a 4D variable
1062	if (nbdim.ne.4) then
1063	write(,) "Error:",trim(GCMvarname)," is not a 4D variable"
1064	write(,) "We'll skip that variable...";write(,)""
1065	CYCLE VAR ! go directly to the next variable
1066	endif
1067	! Check that its dimensions are indeed lon,lat,alt,time (in the right order)
1068	if (GCMvarshape(1).ne.lon_dimid_gcm) then
1069	write(,) "Error, expected first dimension of ",trim(GCMvarname)," to be longitude!"
1070	write(,) "We'll skip that variable..."
1071	CYCLE VAR ! go directly to the next variable
1072	endif
1073	if (GCMvarshape(2).ne.lat_dimid_gcm) then
1074	write(,) "Error, expected second dimension of ",trim(GCMvarname)," to be latitude!"
1075	write(,) "We'll skip that variable..."
1076	CYCLE VAR ! go directly to the next variable
1077	endif
1078	if (GCMvarshape(3).ne.alt_dimid_gcm) then
1079	write(,) "Error, expected third dimension of ",trim(GCMvarname)," to be altitude!"
1080	write(,) "We'll skip that variable..."
1081	CYCLE VAR ! go directly to the next variable
1082	endif
1083	if (GCMvarshape(4).ne.time_dimid_gcm) then
1084	write(,) "Error, expected fourth dimension of ",trim(GCMvarname)," to be time!"
1085	write(,) "We'll skip that variable..."
1086	CYCLE VAR ! go directly to the next variable
1087	endif
1088
1089	! Length allocation for each dimension of the 4D variable
1090	allocate(GCM_var(GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen))
1091
1092	! Load datasets
1093	if (.not.is_stats) then
1094	status=NF90_GET_VAR(gcmfid,GCMvarid,GCM_var)
1095	error_text="Failed to load "//trim(GCMvarname)
1096	call status_check(status,error_text)
1097	else
1098	! if it is a stats file, we load only the first sol, and then copy it to all the other sols
1099	status=NF90_GET_VAR(gcmfid,GCMvarid,GCM_var(:,:,:,1:GCMstatstimelen))
1100	error_text="Failed to load "//trim(GCMvarname)
1101	call status_check(status,error_text)
1102	! write(,)"GCMstatstimelen = ", GCMstatstimelen
1103	! write(,)"GCMtimelen = ", GCMtimelen
1104	do l=(GCMstatstimelen+1),GCMtimelen
1105	if (modulo(l,GCMstatstimelen).ne.0) then
1106	GCM_var(:,:,:,l) = GCM_var(:,:,:,modulo(l,GCMstatstimelen))
1107	else ! if l is a multiple of GCMstatstimelen, since the index modulo(l,GCMstatstimelen)=0
1108	! doesn't exist, we make a special case
1109	GCM_var(:,:,:,l) = GCM_var(:,:,:,GCMstatstimelen)
1110	endif
1111	enddo
1112	endif
1113	write(,) "Variable ",trim(GCMvarname)," loaded from the GCM file"
1114
1115	! Get dataset's missing_value attribute
1116	status=nf90_get_att(gcmfid,GCMvarid,"missing_value",GCMmiss_val)
1117	error_text="Failed to load missing_value attribute"
1118	call status_check(status,error_text)
1119
1120	! Get other variable's attributes
1121	status=nf90_get_att(gcmfid,GCMvarid,"long_name",long_name)
1122	if (status.ne.nf90_noerr) then
1123	! if no attribute "long_name", try "title"
1124	status=nf90_get_att(gcmfid,GCMvarid,"title",long_name)
1125	endif
1126	status=nf90_get_att(gcmfid,GCMvarid,"units",units)
1127
1128	!================================================================
1129	! 2.4.2 Handle dust and wice opacities (second part)
1130	!================================================================
1131	! DUST
1132	!-----
1133	if (trim(gcm_vars(v)).eq."dust") then
1134
1135	IF (dustok1.or.dustok2) THEN
1136	! Dust opacity computed from its density-scaled opacity
1137	! write(,)long_name(index(long_name,"(")+1:index(long_name,")")-1)
1138	do i=1,GCMlonlen
1139	do j=1,GCMlatlen
1140	do k=1,GCMaltlen
1141	do l=1,GCMtimelen
1142	if (GCM_var(i,j,k,l).ne.GCMmiss_val) then
1143	! Multiply by rho to have opacity [1/km]
1144	GCM_var(i,j,k,l) = GCM_var(i,j,k,l) * rho(i,j,k,l) *1000.
1145
1146	if (long_name(index(long_name,"(")+1:index(long_name,")")-1).eq."TES") then
1147	! The density-scaled opacity was calibrated on TES wavelength (9.3um)
1148	! so we must recalibrate it to MCS wavelength (21.6um) using recalibration
1149	! coefficients from Montabone et al. 2015, section 2.3
1150	GCM_var(i,j,k,l) = 1.3/2.7 * GCM_var(i,j,k,l)
1151	endif
1152	endif
1153	enddo
1154	enddo
1155	enddo
1156	enddo
1157
1158	long_name = "IR Dust opacity (from DSO)"
1159
1160	ELSE IF (dustok3) THEN
1161	! Dust opacity computed from its mass mixing ratio
1162	do i=1,GCMlonlen
1163	do j=1,GCMlatlen
1164	do k=1,GCMaltlen
1165	do l=1,GCMtimelen
1166	if (GCM_var(i,j,k,l).ne.GCMmiss_val) then
1167	! Opacity is computed from the equation of Heavens et al. 2014, section 2.3,
1168	! assuming a rho_dust=3000 kg/m3 and an effective radius reff=1.06 microns
1169	! + the opacity is here in 1/km and the MMR in kg/kg
1170	GCM_var(i,j,k,l) = GCM_var(i,j,k,l) * rho(i,j,k,l) / 0.012 * 1000
1171	endif
1172	enddo
1173	enddo
1174	enddo
1175	enddo
1176
1177	long_name = "IR Dust opacity (from MMR)"
1178
1179	ENDIF
1180
1181	GCMvarname = gcm_vars(v) ! reput the right name in GCMvarname
1182	units = "opacity/km"
1183	endif ! trim(gcm_vars(v)).eq."dust"
1184
1185
1186	! WICE
1187	!-----
1188	if (trim(gcm_vars(v)).eq."wice") then
1189	! Water ice opacity computed from its mass mixing ratio
1190	do i=1,GCMlonlen
1191	do j=1,GCMlatlen
1192	do k=1,GCMaltlen
1193	do l=1,GCMtimelen
1194	if (GCM_var(i,j,k,l).ne.GCMmiss_val) then
1195	! Opacity at MCS wavelength (11.9um) is computed from an equation
1196	! similar to the one of Heavens et al. 2014, section 2.3.
1197	! We assume a rho_wice=920 kg/m3, an effective radius reff=3um,
1198	! an extinction coefficient Qext(wvl,reff)=1.54471
1199	GCM_var(i,j,k,l) = 7501.54471 GCM_var(i,j,k,l) * rho(i,j,k,l) / (920*3e-6)
1200	endif
1201	enddo
1202	enddo
1203	enddo
1204	enddo
1205
1206	long_name = "IR Water ice opacity (from MMR)"
1207
1208	GCMvarname = gcm_vars(v) ! reput the right name in GCMvarname
1209	units = "opacity/km"
1210	endif ! trim(gcm_vars(v)).eq."wice"
1211
1212	!===============================================================================
1213	! 2.5 Opening of the associated MCS variables
1214	!===============================================================================
1215	! Observer variables to extract :
1216	IF ((index(GCMvarname,"dust").ne.0).or.(index(GCMvarname,"dso").ne.0)) THEN
1217	! if the variable name contains "dust" or "dso". Especially for the targeted variables :
1218	! dustq,dustN,dsodust,reffdust,opadust, and their equivalents for stormdust & topdust
1219	OBSvarname = "dust"
1220	numbin_name = "numbindust"
1221	ELSE IF ((trim(GCMvarname).eq."h2o_ice").or.(trim(GCMvarname).eq."wice") &
1222	.or.(trim(GCMvarname).eq."reffice").or.(trim(GCMvarname).eq."opawice")) THEN
1223	OBSvarname = "wice"
1224	numbin_name = "numbinwice"
1225	ELSE ! default case is temp binning, since it contains the most values
1226	OBSvarname = "temp"
1227	numbin_name = "numbintemp"
1228	ENDIF
1229
1230	SELECT CASE (dayornight)
1231	CASE ("dayside")
1232	OBSvarname = "d"//OBSvarname
1233	numbin_name = "d"//numbin_name
1234	CASE ("nightside")
1235	OBSvarname = "n"//OBSvarname
1236	numbin_name = "n"//numbin_name
1237	END SELECT
1238
1239	!================================================================
1240	! 2.5.1 MCS reference variable (for the missing values)
1241	!================================================================
1242	! Check that the observation file contains that variable
1243	status=nf90_inq_varid(obsfid,trim(OBSvarname),OBSvarid)
1244	error_text="Failed to find variable "//trim(OBSvarname)//" in "//trim(obsfile)
1245	call status_check(status,error_text)
1246
1247	! Sanity checks on the variable
1248	status=nf90_inquire_variable(obsfid,OBSvarid,ndims=nbdim,dimids=OBSvarshape)
1249	error_text="Failed to obtain information on variable "//trim(OBSvarname)
1250	call status_check(status,error_text)
1251
1252	! Check that it is a 4D variable
1253	if (nbdim.ne.4) then
1254	write(,) "Error, expected a 4D (lon-lat-alt-time) variable for ",trim(OBSvarname)
1255	stop
1256	endif
1257	! Check that its dimensions are indeed lon,lat,alt,time (in the right order)
1258	if (OBSvarshape(1).ne.lon_dimid_obs) then
1259	write(,) "Error, expected first dimension of ",trim(OBSvarname)," to be longitude!"
1260	stop
1261	endif
1262	if (OBSvarshape(2).ne.lat_dimid_obs) then
1263	write(,) "Error, expected second dimension of ",trim(OBSvarname)," to be latitude!"
1264	stop
1265	endif
1266	if (OBSvarshape(3).ne.alt_dimid_obs) then
1267	write(,) "Error, expected third dimension of ",trim(OBSvarname)," to be altitude!"
1268	stop
1269	endif
1270	if (OBSvarshape(4).ne.time_dimid_obs) then
1271	write(,) "Error, expected fourth dimension of ",trim(OBSvarname)," to be time!"
1272	stop
1273	endif
1274
1275	! Length allocation for each dimension of the 4D variable
1276	allocate(OBS_var(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
1277
1278	! Load datasets
1279	status=NF90_GET_VAR(obsfid,OBSvarid,OBS_var)
1280	error_text="Failed to load "//trim(OBSvarname)//" from the obsfile"
1281	call status_check(status,error_text)
1282	write(,) trim(OBSvarname)," loaded from the obsfile as reference variable"
1283
1284	! Get OBS_var missing_value attribute
1285	status=nf90_get_att(obsfid,OBSvarid,"_FillValue",OBSmiss_val)
1286	error_text="Failed to load missing_value attribute"
1287	call status_check(status,error_text)
1288
1289	!================================================================
1290	! 2.5.2 Number of values in the OBS bin (for the sol binning)
1291	!================================================================
1292	! Check that the observation file contains that variable
1293	status=nf90_inq_varid(obsfid,trim(numbin_name),numbin_id)
1294
1295	! Checks have already been done in section 2.3
1296
1297	! Length allocation for each dimension of the 4D variable
1298	allocate(numbin(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
1299
1300	! Load datasets
1301	status=NF90_GET_VAR(obsfid,numbin_id,numbin)
1302
1303	!===============================================================================
1304	! 2.6 Definition of GCM variables in outfile
1305	!===============================================================================
1306	! Switch to netcdf define mode
1307	status=nf90_redef(outfid)
1308	error_text="Error: could not switch to define mode in the outfile"
1309	call status_check(status,error_text)
1310
1311	! Definition of the variable
1312	SELECT CASE (dayornight)
1313	CASE ("dayside")
1314	outvarname = "d"//GCMvarname
1315	CASE ("nightside")
1316	outvarname = "n"//GCMvarname
1317	END SELECT
1318	status=NF90_DEF_VAR(outfid,trim(outvarname),nf90_float,OBSvarshape,outvarid)
1319	error_text="Error: could not define the variable "//trim(outvarname)//" in the outfile"
1320	call status_check(status,error_text)
1321
1322	! Write the attributes
1323	SELECT CASE (dayornight)
1324	CASE ("dayside")
1325	long_name = trim(long_name)//" - day side"
1326	status=nf90_put_att(outfid,outvarid,"long_name",long_name)
1327	CASE ("nightside")
1328	long_name = trim(long_name)//" - night side"
1329	status=nf90_put_att(outfid,outvarid,"long_name",long_name)
1330	END SELECT
1331	status=nf90_put_att(outfid,outvarid,"units",units)
1332	status=nf90_put_att(outfid,outvarid,"_FillValue",OBSmiss_val)
1333	comment = "Reference numbin: "//trim(numbin_name)
1334	status=nf90_put_att(outfid,outvarid,"comment",comment)
1335
1336	write(,)trim(outvarname)," has been created in the outfile"
1337	write(*,'(" with missing_value attribute : ",1pe12.5)')OBSmiss_val
1338	write(,)""
1339
1340	! End the netcdf define mode (and thus enter the "data writing" mode)
1341	status=nf90_enddef(outfid)
1342	error_text="Error: could not close the define mode of the outfile"
1343	call status_check(status,error_text)
1344
1345	allocate(outvar(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
1346
1347
1348	!===================================================================================
1349	! 3. EXTRACTION OF THE VARIABLE
1350	!===================================================================================
1351	!===============================================================================
1352	! 3.1 Do some checks and preparations before the extraction
1353	!===============================================================================
1354	Ls: do l=1,OBSLslen ! loop on all observed Solar longitudes
1355	Ls_val=OBSLs(l) ! Ls_val=center of the output bin
1356	if ((Ls_val.lt.0.).or.(Ls_val.gt.360.)) then
1357	write(,) "Unexpected value for OBSLs: ",Ls_val
1358	stop
1359	endif
1360
1361	! Convert the Ls bin into a sol interval on which the binning is done :
1362	!----------------------------------------------------------------------
1363	!-> get the index of the maximum value of GCM sol (m_maxsol) that is lower than Ls bin's superior bound (maxsol)
1364	call ls2sol(Ls_val+OBSdeltaLs/2.,maxsol)
1365	m_maxsol=1
1366	do while (((GCMtime(m_maxsol)+0.5).lt.maxsol).AND.(m_maxsol.le.(GCMtimelen-1)))
1367	!! The +0.5 is there to take into account the whole planet (lon=[-180°;180°]) and not just the lon=0° from the GCM
1368	m_maxsol=m_maxsol+1
1369	enddo
1370	!-> get the index of the minimum value of GCM sol (m_minsol) that is greater than Ls bin's inferior bound (minsol)
1371	call ls2sol(Ls_val-OBSdeltaLs/2.,minsol)
1372	m_minsol=1
1373	do while (((GCMtime(m_minsol)-0.5).le.minsol).AND.(m_minsol.le.(GCMtimelen-1)))
1374	!! Same comment for the -0.5
1375	m_minsol=m_minsol+1
1376	enddo
1377	if (m_minsol.gt.m_maxsol) then
1378	write(,) "No value in gcmfile between sol=",minsol," and sol=",maxsol," (Ls=",Ls_val,"Â°)"
1379	! Write a missing_value to output
1380	outvar(:,:,:,l) = OBSmiss_val
1381	CYCLE Ls ! go directly to the next Ls
1382	endif
1383	! Get all the integer values of GCM sols that fit in this interval
1384	solcount=floor(GCMtime(m_maxsol))-ceiling(GCMtime(m_minsol))+1
1385	! sols that are not fully in the interval are not counted
1386	allocate(int_sol_list(solcount))
1387	! write(,) "GCMminsol=", GCMtime(m_minsol)
1388	! write(,)"GCMmaxsol=", GCMtime(m_maxsol)
1389	do m=1,solcount
1390	int_sol_list(m)=ceiling(GCMtime(m_minsol)) + m-1
1391	enddo
1392	! write(,)"int_sol_list=",int_sol_list
1393
1394
1395	latitude: do j=1,OBSlatlen ! loop on all observed latitudes
1396	lat_val=OBSlat(j)
1397	if ((lat_val.lt.-90.).or.(lat_val.gt.90.)) then
1398	write(,) "Unexpected value for OBSlat: ",lat_val
1399	stop
1400	endif
1401
1402	longitude: do i=1,OBSlonlen ! loop on all observed longitudes
1403	lon_val=OBSlon(i)
1404	if ((lon_val.lt.-360.).or.(lon_val.gt.360.)) then
1405	write(,) "Unexpected value for lon_val: ",lon_val
1406	stop
1407	endif
1408	! We want lon_val in [-180:180] for the subroutine extraction
1409	if (lon_val.lt.-180.) lon_val=lon_val+360.
1410	if (lon_val.gt.180.) lon_val=lon_val-360.
1411
1412	LT_val=OBSLT(i,j,l) ! find the Observer average LT value at bin(lon_val, lat_val, Ls_val)
1413
1414	if ((LT_val.lt.0.).or.(LT_val.gt.24.)) then
1415	if (LT_val.eq.LTmiss_val) then
1416	! write(,) "Missing value in obsfile for LT_val"
1417	! Write a missing_value to output
1418	outvar(i,j,:,l) = OBSmiss_val
1419	CYCLE longitude ! go directly to the next longitude
1420	else
1421	write(,) "Unexpected value for LT_val: ",LT_val
1422	stop
1423	endif
1424	endif
1425
1426	altitude: do k=1,OBSaltlen ! loop on all observed altitudes
1427	alt_val=OBSalt(k)
1428	if (OBS_var(i,j,k,l).eq.OBSmiss_val) then
1429	! write(,) "Missing value in obsfile for ",OBSvarname
1430	! Write a missing_value to output
1431	outvar(i,j,k,l) = OBSmiss_val
1432	CYCLE altitude ! go directly to the next altitude
1433	endif
1434
1435	!===============================================================================
1436	! 3.2 Compute GCM sol date corresponding to Observer Ls (via m_(min/max)sol)
1437	! and LT (via OBSLT(min/max))
1438	!===============================================================================
1439	LTcount=floor(numbin(i,j,k,l)) ! find the Observer number of temp values
1440	! at bin(lon_val,lat_val,alt_val,Ls_val)
1441	if (LTcount.eq.0.) then
1442	! Write a missing_value to output
1443	outvar(i,j,k,l) = OBSmiss_val
1444	CYCLE altitude ! go directly to the next altitude
1445	endif
1446	if (LTcount.lt.0.) then
1447	write(,) "Unexpected value for LTcount: ",LTcount
1448	stop
1449	endif
1450
1451	! Generate the sol list for the interpolation
1452	allocate(sol_list(solcount*LTcount))
1453	call gen_sol_list(solcount,int_sol_list,LTcount,LT_val,OBSLTmax(i,j,l),&
1454	OBSLTmin(i,j,l),lon_val,LTmod,dayornight,&
1455	sol_list)
1456
1457	solerrcount=0
1458	solbinned_value=0
1459	sol_bin: do m=1,solcount*LTcount ! loop on all GCM sols of the bin
1460	sol=sol_list(m)
1461	! write(,)"sol=",sol
1462	!===============================================================================
1463	! 3.3 Do the interpolation and binning for the given location
1464	!===============================================================================
1465	call extraction(lon_val,lat_val,alt_val,sol,&
1466	GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen,&
1467	GCMlon,GCMlat,GCMalt,GCMtime,&
1468	GCM_var,GCMmiss_val,GCMalttype,GCMvarname,extr_value)
1469
1470	if (extr_value.eq.GCMmiss_val) then
1471	! write(,) "Missing value in gcmfile at lon=",lon_val,"; lat=",lat_val,"; alt=",alt_val,"; sol=",sol
1472	solerrcount=solerrcount+1
1473	CYCLE sol_bin ! go directly to the next GCM sol of the bin
1474	endif
1475	solbinned_value=solbinned_value+extr_value
1476
1477	enddo sol_bin ! end loop on all GCM sols of the bin
1478	if ((solcount*LTcount-solerrcount).ne.0) then
1479	solbinned_value=solbinned_value/(solcount*LTcount-solerrcount)
1480	else
1481	! write(,)"No GCM value in this sol bin"
1482	solbinned_value=OBSmiss_val
1483	endif
1484	! Write value to output
1485	outvar(i,j,k,l)=solbinned_value
1486
1487	errcount=errcount+solerrcount
1488
1489	deallocate(sol_list)
1490
1491	enddo altitude ! end loop on observed altitudes
1492	enddo longitude ! end loop on observed longitudes
1493	enddo latitude ! end loop on observed latitudes
1494
1495	deallocate(int_sol_list)
1496
1497	enddo Ls ! end loop on observed Solar longitudes
1498
1499	! write(,)"Nb of GCM missing values :",errcount
1500
1501	!===============================================================================
1502	! 3.4 Write the data in the netcdf output file
1503	!===============================================================================
1504	status = nf90_put_var(outfid, outvarid, outvar)
1505	error_text="Error: could not write "//trim(outvarname)//" data in the outfile"
1506	call status_check(status,error_text)
1507
1508	! Deallocations before going to the next GCM variable
1509	deallocate(GCM_var)
1510	deallocate(OBS_var)
1511	deallocate(numbin)
1512	deallocate(outvar)
1513	ENDDO VAR ! end loop on variables
1514
1515	!===================================================================================
1516	! 4. END OF THE DAY/NIGHT LOOP
1517	!===================================================================================
1518	IF (dayornight.EQ."dayside") THEN
1519	! this is the end of the first loop (on daytime values)
1520	! and we still have to loop on nighttime values
1521	dayornight="nightside"
1522	deallocate(OBSLT)
1523	deallocate(OBSLTmax)
1524	deallocate(OBSLTmin)
1525	write(,)""
1526	write(,) "Beginning the 2nd loop, on nighttime values"; write(,)""
1527	CYCLE DAY_OR_NIGHT
1528	ELSE ! i.e. dayornight="nightside"
1529	! this is the end of the second loop (on nighttime values)
1530	! and thus the end of the program
1531	EXIT DAY_OR_NIGHT
1532	ENDIF
1533	ENDDO DAY_OR_NIGHT ! end of the day/night loop that begins in section 1.3
1534
1535	!===================================================================================
1536	! 5. CLOSE THE FILES
1537	!===================================================================================
1538	status = nf90_close(gcmfid)
1539	error_text="Error: could not close file "//trim(gcmfile)
1540	call status_check(status,error_text)
1541	status = nf90_close(obsfid)
1542	error_text="Error: could not close file "//trim(obsfile)
1543	call status_check(status,error_text)
1544	status = nf90_close(outfid)
1545	error_text="Error: could not close file "//trim(outfile)
1546	call status_check(status,error_text)
1547
1548	write(,)"";write(,)"-> Program simu_MCS completed!"
1549
1550	end program simu_MCS
1551
1552
1553
1554	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1555	!===================================================================================================
1556	! Subroutines
1557	!===================================================================================================
1558
1559	subroutine extraction(lon,lat,alt,sol,&
1560	lonlen,latlen,altlen,timelen,&
1561	longitude,latitude,altitude,time,&
1562	field,missing_value,alttype,varname,value)
1563
1564	implicit none
1565	!================================================================
1566	! Arguments:
1567	!================================================================
1568	real,intent(in) :: lon ! sought longitude (deg, in [-180:180])
1569	real,intent(in) :: lat ! sought latitude (deg, in [-90:90])
1570	real,intent(in) :: alt ! sought altitude (m or Pa)
1571	real,intent(in) :: sol ! sought date (sols)
1572	integer,intent(in) :: lonlen
1573	integer,intent(in) :: latlen
1574	integer,intent(in) :: altlen
1575	integer,intent(in) :: timelen
1576	real,intent(in) :: longitude(lonlen)
1577	real,intent(in) :: latitude(latlen)
1578	real,intent(in) :: altitude(altlen)
1579	real,intent(in) :: time(timelen)
1580	real,intent(in) :: field(lonlen,latlen,altlen,timelen)
1581	real,intent(in) :: missing_value ! default value for "no data"
1582	character,intent(in) :: alttype ! altitude coord. type:'z' (altitude, m)
1583	! 'p' (pressure, Pa)
1584	character(len=*),intent(in) :: varname ! variable name (in GCM file)
1585	real,intent(out) :: value
1586
1587	!================================================================
1588	! Local variables:
1589	!================================================================
1590	real,save :: prev_lon=-99 ! previous value of 'lon' routine was called with
1591	real,save :: prev_lat=-99 ! previous value of 'lat' routine was called with
1592	real,save :: prev_alt=-9.e20 ! ! previous value of 'alt'
1593	real,save :: prev_sol=-99 ! previous value of 'sol' routine was called with
1594
1595	! encompasing indexes:
1596	integer,save :: ilon_inf=-1,ilon_sup=-1 ! along longitude
1597	integer,save :: ilat_inf=-1,ilat_sup=-1 ! along latitude
1598	integer,save :: ialt_inf=-1,ialt_sup=-1 ! along altitude
1599	integer,save :: itim_inf=-1,itim_sup=-1 ! along time
1600
1601	! intermediate interpolated values
1602	real :: t_interp(2,2,2) ! after time interpolation
1603	real :: zt_interp(2,2) ! after altitude interpolation
1604	real :: yzt_interp(2) ! after latitude interpolation
1605	real :: coeff ! interpolation coefficient
1606
1607	integer :: i
1608
1609	! By default, set value to missing_value
1610	value=missing_value
1611
1612	!================================================================
1613	! 1. Find encompassing indexes
1614	!================================================================
1615	if (lon.ne.prev_lon) then
1616	do i=1,lonlen-1
1617	if (longitude(i).le.lon) then
1618	ilon_inf=i
1619	else
1620	exit
1621	endif
1622	enddo
1623	ilon_sup=ilon_inf+1
1624	endif ! of if (lon.ne.prev_lon)
1625	!write(,) 'ilon_inf=',ilon_inf," longitude(ilon_inf)=",longitude(ilon_inf)
1626
1627	if (lat.ne.prev_lat) then
1628	if (latitude(1).gt.latitude(2)) then
1629	! decreasing latitudes, from 90N to (-)90S (LMDZ-like)
1630	do i=1,latlen-1
1631	if (latitude(i).ge.lat) then
1632	ilat_inf=i
1633	else
1634	exit
1635	endif
1636	enddo
1637	else
1638	! increasing latitudes, from (-)90S to 90N (DYNAMICO-like)
1639	do i=1,latlen-1
1640	if (latitude(i).le.lat) then
1641	ilat_inf=i
1642	else
1643	exit
1644	endif
1645	enddo
1646	endif
1647	ilat_sup=ilat_inf+1
1648	endif ! of if (lat.ne.prev_lat)
1649	!write(,) 'ilat_inf=',ilat_inf," latitude(ilat_inf)=",latitude(ilat_inf)
1650
1651	if (alt.ne.prev_alt) then
1652	if (alttype.eq.'p') then ! pressures are ordered from max to min
1653	!handle special case for alt not in the altitude(1:altlen) interval
1654	if ((alt.gt.altitude(1)).or.(alt.lt.altitude(altlen))) then
1655	ialt_inf=-1
1656	ialt_sup=-1
1657	! return to main program (with value=missing_value)
1658	! write(,)"Problem in extraction : GCM alt p"
1659	return
1660	else ! general case
1661	do i=1,altlen-1
1662	if (altitude(i).ge.alt) then
1663	ialt_inf=i
1664	else
1665	exit
1666	endif
1667	enddo
1668	ialt_sup=ialt_inf+1
1669	endif ! of if ((alt.gt.altitude(1)).or.(alt.lt.altitude(altlen)))
1670	else ! altitudes (m) are ordered from min to max
1671	!handle special case for alt not in the altitude(1:altlen) interval
1672	if ((alt.lt.altitude(1)).or.(alt.gt.altitude(altlen))) then
1673	ialt_inf=-1
1674	ialt_sup=-1
1675	! return to main program (with value=missing_value)
1676	! write(,)"Problem in extraction : GCM alt z"
1677	return
1678	else ! general case
1679	do i=1,altlen-1
1680	if (altitude(i).le.alt) then
1681	ialt_inf=i
1682	else
1683	exit
1684	endif
1685	enddo
1686	ialt_sup=ialt_inf+1
1687	endif ! of if ((alt.lt.altitude(1)).or.(alt.gt.altitude(altlen)))
1688	endif ! of if (alttype.eq.'p')
1689	endif ! of if (alt.ne.prev_alt)
1690	!write(,) 'ialt_inf=',ialt_inf," altitude(ialt_inf)=",altitude(ialt_inf)
1691
1692	if (sol.ne.prev_sol) then
1693	!handle special case for sol not in the time(1):time(timenlen) interval
1694	if ((sol.lt.time(1)).or.(sol.gt.time(timelen))) then
1695	itim_inf=-1
1696	itim_sup=-1
1697	! return to main program (with value=missing_value)
1698	! write(,)"Problem in extraction : GCM sol"
1699	return
1700	else ! general case
1701	do i=1,timelen-1
1702	if (time(i).le.sol) then
1703	itim_inf=i
1704	else
1705	exit
1706	endif
1707	enddo
1708	itim_sup=itim_inf+1
1709	endif ! of if ((sol.lt.time(1)).or.(sol.gt.time(timenlen)))
1710	endif ! of if (sol.ne.prev_sol)
1711	!write(,) 'itim_inf=',itim_inf," time(itim_inf)=",time(itim_inf)
1712	!write(,) 'itim_sup=',itim_sup," time(itim_sup)=",time(itim_sup)
1713
1714	!================================================================
1715	! 2. Interpolate
1716	!================================================================
1717	! check that there are no "missing_value" in the field() elements we need
1718	! otherwise return to main program (with value=missing_value)
1719	if (field(ilon_inf,ilat_inf,ialt_inf,itim_inf).eq.missing_value) then
1720	! write(,)"Error 1 in extraction"
1721	return
1722	endif
1723	if (field(ilon_inf,ilat_inf,ialt_inf,itim_sup).eq.missing_value) then
1724	! write(,)"Error 2 in extraction"
1725	return
1726	endif
1727	if (field(ilon_inf,ilat_inf,ialt_sup,itim_inf).eq.missing_value) then
1728	! write(,)"Error 3 in extraction"
1729	return
1730	endif
1731	if (field(ilon_inf,ilat_inf,ialt_sup,itim_sup).eq.missing_value) then
1732	! write(,)"Error 4 in extraction"
1733	return
1734	endif
1735	if (field(ilon_inf,ilat_sup,ialt_inf,itim_inf).eq.missing_value) then
1736	! write(,)"Error 5 in extraction"
1737	return
1738	endif
1739	if (field(ilon_inf,ilat_sup,ialt_inf,itim_sup).eq.missing_value) then
1740	! write(,)"Error 6 in extraction"
1741	return
1742	endif
1743	if (field(ilon_inf,ilat_sup,ialt_sup,itim_inf).eq.missing_value) then
1744	! write(,)"Error 7 in extraction"
1745	return
1746	endif
1747	if (field(ilon_inf,ilat_sup,ialt_sup,itim_sup).eq.missing_value) then
1748	! write(,)"Error 8 in extraction"
1749	return
1750	endif
1751	if (field(ilon_sup,ilat_inf,ialt_inf,itim_inf).eq.missing_value) then
1752	! write(,)"Error 9 in extraction"
1753	return
1754	endif
1755	if (field(ilon_sup,ilat_inf,ialt_inf,itim_sup).eq.missing_value) then
1756	! write(,)"Error 10 in extraction"
1757	return
1758	endif
1759	if (field(ilon_sup,ilat_inf,ialt_sup,itim_inf).eq.missing_value) then
1760	! write(,)"Error 11 in extraction"
1761	return
1762	endif
1763	if (field(ilon_sup,ilat_inf,ialt_sup,itim_sup).eq.missing_value) then
1764	! write(,)"Error 12 in extraction"
1765	return
1766	endif
1767	if (field(ilon_sup,ilat_sup,ialt_inf,itim_inf).eq.missing_value) then
1768	! write(,)"Error 13 in extraction"
1769	return
1770	endif
1771	if (field(ilon_sup,ilat_sup,ialt_inf,itim_sup).eq.missing_value) then
1772	! write(,)"Error 14 in extraction"
1773	return
1774	endif
1775	if (field(ilon_sup,ilat_sup,ialt_sup,itim_inf).eq.missing_value) then
1776	! write(,)"Error 15 in extraction"
1777	return
1778	endif
1779	if (field(ilon_sup,ilat_sup,ialt_sup,itim_sup).eq.missing_value) then
1780	! write(,)"Error 16 in extraction"
1781	return
1782	endif
1783
1784	!================================================================
1785	! 2.1 Linear interpolation in time
1786	!================================================================
1787	coeff=(sol-time(itim_inf))/(time(itim_sup)-time(itim_inf))
1788	t_interp(1,1,1)=field(ilon_inf,ilat_inf,ialt_inf,itim_inf)+ &
1789	coeff*(field(ilon_inf,ilat_inf,ialt_inf,itim_sup)- &
1790	field(ilon_inf,ilat_inf,ialt_inf,itim_inf))
1791	t_interp(1,1,2)=field(ilon_inf,ilat_inf,ialt_sup,itim_inf)+ &
1792	coeff*(field(ilon_inf,ilat_inf,ialt_sup,itim_sup)- &
1793	field(ilon_inf,ilat_inf,ialt_sup,itim_inf))
1794	t_interp(1,2,1)=field(ilon_inf,ilat_sup,ialt_inf,itim_inf)+ &
1795	coeff*(field(ilon_inf,ilat_sup,ialt_inf,itim_sup)- &
1796	field(ilon_inf,ilat_sup,ialt_inf,itim_inf))
1797	t_interp(1,2,2)=field(ilon_inf,ilat_sup,ialt_sup,itim_inf)+ &
1798	coeff*(field(ilon_inf,ilat_sup,ialt_sup,itim_sup)- &
1799	field(ilon_inf,ilat_sup,ialt_sup,itim_inf))
1800	t_interp(2,1,1)=field(ilon_sup,ilat_inf,ialt_inf,itim_inf)+ &
1801	coeff*(field(ilon_sup,ilat_inf,ialt_inf,itim_sup)- &
1802	field(ilon_sup,ilat_inf,ialt_inf,itim_inf))
1803	t_interp(2,1,2)=field(ilon_sup,ilat_inf,ialt_sup,itim_inf)+ &
1804	coeff*(field(ilon_sup,ilat_inf,ialt_sup,itim_sup)- &
1805	field(ilon_sup,ilat_inf,ialt_sup,itim_inf))
1806	t_interp(2,2,1)=field(ilon_sup,ilat_sup,ialt_inf,itim_inf)+ &
1807	coeff*(field(ilon_sup,ilat_sup,ialt_inf,itim_sup)- &
1808	field(ilon_sup,ilat_sup,ialt_inf,itim_inf))
1809	t_interp(2,2,2)=field(ilon_sup,ilat_sup,ialt_sup,itim_inf)+ &
1810	coeff*(field(ilon_sup,ilat_sup,ialt_sup,itim_sup)- &
1811	field(ilon_sup,ilat_sup,ialt_sup,itim_inf))
1812
1813	!================================================================
1814	! 2.2 Vertical interpolation
1815	!================================================================
1816	if (((varname=='rho').or.(varname=='pressure')).and.(alttype=='z')) then
1817	! do the interpolation on the log of the quantity
1818	coeff=(alt-altitude(ialt_inf))/(altitude(ialt_sup)-altitude(ialt_inf))
1819	zt_interp(1,1)=log(t_interp(1,1,1))+coeff* &
1820	(log(t_interp(1,1,2))-log(t_interp(1,1,1)))
1821	zt_interp(1,2)=log(t_interp(1,2,1))+coeff* &
1822	(log(t_interp(1,2,2))-log(t_interp(1,2,1)))
1823	zt_interp(2,1)=log(t_interp(2,1,1))+coeff* &
1824	(log(t_interp(2,1,2))-log(t_interp(2,1,1)))
1825	zt_interp(2,2)=log(t_interp(2,2,1))+coeff* &
1826	(log(t_interp(2,2,2))-log(t_interp(2,2,1)))
1827	zt_interp(1:2,1:2)=exp(zt_interp(1:2,1:2))
1828	else ! general case
1829	coeff=(alt-altitude(ialt_inf))/(altitude(ialt_sup)-altitude(ialt_inf))
1830	zt_interp(1,1)=t_interp(1,1,1)+coeff*(t_interp(1,1,2)-t_interp(1,1,1))
1831	zt_interp(1,2)=t_interp(1,2,1)+coeff*(t_interp(1,2,2)-t_interp(1,2,1))
1832	zt_interp(2,1)=t_interp(2,1,1)+coeff*(t_interp(2,1,2)-t_interp(2,1,1))
1833	zt_interp(2,2)=t_interp(2,2,1)+coeff*(t_interp(2,2,2)-t_interp(2,2,1))
1834	endif
1835
1836	!================================================================
1837	! 2.3 Latitudinal interpolation
1838	!================================================================
1839	coeff=(lat-latitude(ilat_inf))/(latitude(ilat_sup)-latitude(ilat_inf))
1840	yzt_interp(1)=zt_interp(1,1)+coeff*(zt_interp(1,2)-zt_interp(1,1))
1841	yzt_interp(2)=zt_interp(2,1)+coeff*(zt_interp(2,2)-zt_interp(2,1))
1842
1843	!================================================================
1844	! 2.4 longitudinal interpolation
1845	!================================================================
1846	coeff=(lon-longitude(ilon_inf))/(longitude(ilon_sup)-longitude(ilon_inf))
1847	value=yzt_interp(1)+coeff*(yzt_interp(2)-yzt_interp(1))
1848
1849	end subroutine extraction
1850
1851	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1852
1853	subroutine inidim(outfid,lonlen,latlen,altlen,timelen,lon,lat,alt,time,units_alt,&
1854	londimid,latdimid,altdimid,timedimid)
1855	!================================================================
1856	! Purpose:
1857	! Initialize a data file (NetCDF format)
1858	!================================================================
1859	! Remarks:
1860	! The NetCDF file remains open
1861	!================================================================
1862	use netcdf ! NetCDF definitions
1863	implicit none
1864	!================================================================
1865	! Arguments:
1866	!================================================================
1867	integer, intent(in):: outfid
1868	! outfid: [netcdf] file ID
1869	integer, intent(in):: lonlen
1870	! lonlen: longitude length (# of longitude values)
1871	integer, intent(in):: latlen
1872	! latlen: latitude length (# of latitude values)
1873	integer, intent(in):: altlen
1874	! altlen: altitude length (# of altitude values)
1875	integer, intent(in):: timelen
1876	! timelen: time length (# of time values)
1877	real, intent(in):: lon(lonlen)
1878	! lon(): longitude
1879	real, intent(in):: lat(latlen)
1880	! lat(): latitude
1881	real, intent(in):: alt(altlen)
1882	! alt(): altitude
1883	real, intent(in):: time(timelen)
1884	! time(): time (Ls)
1885	character(len=1), intent(in) :: units_alt
1886	! units_alt: altitude coord. type:'z' (altitude, m) 'p' (pressure, Pa)
1887	integer,intent(inout):: londimid
1888	! londimid: [netcdf] lon() (i.e.: longitude) ID in MCS and output files (they are the same)
1889	integer,intent(inout) :: latdimid
1890	! latdimid: [netcdf] lat() ID in MCS and output files (they are the same)
1891	integer,intent(inout):: altdimid
1892	! altdimid: [netcdf] alt() ID in MCS and output files (they are the same)
1893	integer,intent(inout):: timedimid
1894	! timedimid: [netcdf] time() ID in MCS and output files (they are the same)
1895
1896	!================================================================
1897	! Local variables:
1898	!================================================================
1899	integer :: lonvarid,latvarid,altvarid,timevarid,status
1900	! *varid: [netcdf] ID of a variable
1901	! status: [netcdf] return error code (from called subroutines)
1902	character(len=256) :: error_text
1903	integer :: d ! loop index on dimensions ID
1904
1905	!===============================================================
1906	! 1. Define/write "dimensions" in the same order than MCS file
1907	! and get their IDs
1908	!================================================================
1909	do d=1,4
1910	if (altdimid.eq.d) then
1911	status=nf90_def_dim(outfid, "altitude", altlen, altdimid)
1912	error_text="Error: could not define the altitude dimension in the outfile"
1913	call status_check(status,error_text)
1914	else if (timedimid.eq.d) then
1915	status=nf90_def_dim(outfid, "time", timelen, timedimid)
1916	error_text="Error: could not define the time dimension in the outfile"
1917	call status_check(status,error_text)
1918	else if (latdimid.eq.d) then
1919	status=nf90_def_dim(outfid, "latitude", latlen, latdimid)
1920	error_text="Error: could not define the latitude dimension in the outfile"
1921	call status_check(status,error_text)
1922	else if (londimid.eq.d) then
1923	status=nf90_def_dim(outfid, "longitude", lonlen, londimid)
1924	error_text="Error: could not define the longitude dimension in the outfile"
1925	call status_check(status,error_text)
1926	endif
1927	enddo
1928
1929	!================================================================
1930	! 2. Define "variables" and their attributes
1931	!================================================================
1932	!================================================================
1933	! 2.1 Write "longitude" (data and attributes)
1934	!================================================================
1935
1936	! Insert the definition of the variable
1937	status=nf90_def_var(outfid,"longitude",nf90_float,(/londimid/),lonvarid)
1938	error_text="Error: could not define the longitude variable in the outfile"
1939	call status_check(status,error_text)
1940
1941	! Write the attributes
1942	status=nf90_put_att(outfid,lonvarid,"long_name","longitude")
1943	error_text="Error: could not put the long_name attribute for the longitude variable in the outfile"
1944	call status_check(status,error_text)
1945
1946	status=nf90_put_att(outfid,lonvarid,"units","degrees_east")
1947	error_text="Error: could not put the units attribute for the longitude variable in the outfile"
1948	call status_check(status,error_text)
1949
1950	!================================================================
1951	! 2.2 "latitude"
1952	!================================================================
1953
1954	! Insert the definition of the variable
1955	status=nf90_def_var(outfid,"latitude",nf90_float,(/latdimid/),latvarid)
1956	error_text="Error: could not define the latitude variable in the outfile"
1957	call status_check(status,error_text)
1958
1959	! Write the attributes
1960	status=nf90_put_att(outfid,latvarid,"long_name","latitude")
1961	error_text="Error: could not put the long_name attribute for the latitude variable in the outfile"
1962	call status_check(status,error_text)
1963
1964	status=nf90_put_att(outfid,latvarid,"units","degrees_north")
1965	error_text="Error: could not put the units attribute for the latitude variable in the outfile"
1966	call status_check(status,error_text)
1967
1968	!================================================================
1969	! 2.3 Write "altitude" (data and attributes)
1970	!================================================================
1971
1972	! Insert the definition of the variable
1973	status=nf90_def_var(outfid,"altitude",nf90_float,(/altdimid/),altvarid)
1974	error_text="Error: could not define the altitude variable in the outfile"
1975	call status_check(status,error_text)
1976
1977	! Write the attributes
1978	if (units_alt.eq.'p') then ! pressure coordinate
1979	status=nf90_put_att(outfid,altvarid,"long_name","pressure")
1980	error_text="Error: could not put the long_name attribute for the altitude variable in the outfile"
1981	call status_check(status,error_text)
1982
1983	status=nf90_put_att(outfid,altvarid,'units',"Pa")
1984	error_text="Error: could not put the units attribute for the altitude variable in the outfile"
1985	call status_check(status,error_text)
1986
1987	status=nf90_put_att(outfid,altvarid,'positive',"down")
1988	error_text="Error: could not put the positive attribute for the altitude variable in the outfile"
1989	call status_check(status,error_text)
1990
1991	status=nf90_put_att(outfid,altvarid,'comment',&
1992	"The vertical variable is in fact pressure, not altitude. We just call it altitude for easy reading in Ferret and Grads.")
1993	error_text="Error: could not put the comment attribute for the altitude variable in the outfile"
1994	call status_check(status,error_text)
1995
1996	else if (units_alt.eq.'z') then ! altitude coordinate
1997	status=nf90_put_att(outfid,altvarid,"long_name","altitude")
1998	error_text="Error: could not put the long_name attribute for the altitude variable in the outfile"
1999	call status_check(status,error_text)
2000
2001	status=nf90_put_att(outfid,altvarid,'units',"m")
2002	error_text="Error: could not put the units attribute for the altitude variable in the outfile"
2003	call status_check(status,error_text)
2004
2005	status=nf90_put_att(outfid,altvarid,'positive',"up")
2006	error_text="Error: could not put the positive attribute for the altitude variable in the outfile"
2007	call status_check(status,error_text)
2008
2009	else
2010	write(,)"I do not understand this unit type ",units_alt," for outfile altitude!"
2011	stop
2012	end if
2013
2014	!================================================================
2015	! 2.4 "time"
2016	!================================================================
2017
2018	! Insert the definition of the variable
2019	status=nf90_def_var(outfid,"time",nf90_float,(/timedimid/),timevarid)
2020	error_text="Error: could not define the time variable in the outfile"
2021	call status_check(status,error_text)
2022
2023	! Write the attributes
2024	status=nf90_put_att(outfid,timevarid,"long_name","Solar longitude")
2025	error_text="Error: could not put the long_name attribute for the time variable in the outfile"
2026	call status_check(status,error_text)
2027
2028	status=nf90_put_att(outfid,timevarid,"units","days since 0000-01-1 00:00:00")
2029	error_text="Error: could not put the units attribute for the time variable in the outfile"
2030	call status_check(status,error_text)
2031
2032	status=nf90_put_att(outfid,timevarid,"comment",&
2033	"Units is in fact degrees, but set to a dummy value of days for compatibility with Ferret and Grads.")
2034	error_text="Error: could not put the comment attribute for the time variable in the outfile"
2035	call status_check(status,error_text)
2036
2037	!================================================================
2038	! 2.5 End netcdf define mode
2039	!================================================================
2040	status=nf90_enddef(outfid)
2041	error_text="Error: could not end the define mode of the outfile"
2042	call status_check(status,error_text)
2043
2044	!================================================================
2045	! 3. Write "variables" (data of the dimension variables)
2046	!================================================================
2047	! "time"
2048	status=nf90_put_var(outfid,timevarid,time)
2049	error_text="Error: could not write the time variable's data in the outfile"
2050	call status_check(status,error_text)
2051
2052	! "latitude"
2053	status=nf90_put_var(outfid,latvarid,lat)
2054	error_text="Error: could not write the latitude variable's data in the outfile"
2055	call status_check(status,error_text)
2056
2057	! "longitude"
2058	status=nf90_put_var(outfid,lonvarid,lon)
2059	error_text="Error: could not write the longitude variable's data in the outfile"
2060	call status_check(status,error_text)
2061
2062	! "altitude"
2063	status=nf90_put_var(outfid,altvarid,alt)
2064	error_text="Error: could not write the altitude variable's data in the outfile"
2065	call status_check(status,error_text)
2066
2067	end subroutine inidim
2068
2069	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2070
2071	subroutine ls2sol(ls,sol)
2072
2073	implicit none
2074	!================================================================
2075	! Arguments:
2076	!================================================================
2077	real,intent(in) :: ls
2078	real,intent(out) :: sol
2079
2080	!================================================================
2081	! Local:
2082	!================================================================
2083	double precision xref,zx0,zteta,zz
2084	!xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
2085	double precision year_day
2086	double precision peri_day,timeperi,e_elips
2087	double precision pi,degrad
2088	parameter (year_day=668.6d0) ! number of sols in a martian year
2089	parameter (peri_day=485.35d0) ! date (in sols) of perihelion
2090	!timeperi: 2pi( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
2091	parameter (timeperi=1.90258341759902d0)
2092	parameter (e_elips=0.0934d0) ! eccentricity of orbit
2093	parameter (pi=3.14159265358979d0)
2094	parameter (degrad=57.2957795130823d0)
2095
2096	if (abs(ls).lt.1.0e-5) then
2097	if (ls.ge.0.0) then
2098	sol = 0.0
2099	else
2100	sol = real(year_day)
2101	end if
2102	return
2103	end if
2104
2105	zteta = ls/degrad + timeperi
2106	zx0 = 2.0datan(dtan(0.5zteta)/dsqrt((1.+e_elips)/(1.-e_elips)))
2107	xref = zx0-e_elips*dsin(zx0)
2108	zz = xref/(2.*pi)
2109	sol = real(zz*year_day + peri_day)
2110	if (sol.lt.0.0) sol = sol + real(year_day)
2111	if (sol.ge.year_day) sol = sol - real(year_day)
2112
2113
2114	end subroutine ls2sol
2115
2116	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2117
2118	subroutine gen_sol_list(solcount,int_sol_list,LTcount,LTave,LTmax,LTmin,lon,f_LT,dayornight,&
2119	sol_list)
2120	!================================================================
2121	! Purpose:
2122	! Generate a list that is the combination of two lists :
2123	! - int_sol_list, which is a list of integer values of sol
2124	! - LT_list, which contains a number LTcount of LT values in the
2125	! interval [LTmin;LTmax] which are evenly distributed around
2126	! LTave (so that their average is LTave)
2127	!================================================================
2128
2129	implicit none
2130	!================================================================
2131	! Arguments:
2132	!================================================================
2133	integer, intent(in) :: solcount ! nb of integer values of sol
2134	real, intent(in) :: int_sol_list(solcount) ! list of the integer values of sol
2135	integer, intent(in) :: LTcount ! nb of LT per sol to be interpolated
2136	real, intent(in) :: LTave ! average of LT
2137	real, intent(in) :: LTmax, LTmin ! bounds of the LT interval for the interpolation
2138	real, intent(in) :: lon ! longitude value for the transformation into a sol value at lon=0°
2139	external f_LT ! function that changes LT interval for night LT
2140	real f_LT
2141	character (len=10), intent(in) :: dayornight ! to know if we have day or night values
2142
2143	real, intent(out) :: sol_list(solcount*LTcount) ! all the sol values at lon=0° to interpolate
2144
2145	!================================================================
2146	! Local variables:
2147	!================================================================
2148	integer :: N
2149	integer :: a,b,c ! loop iteration indexes
2150	real :: LT_list(LTcount)
2151
2152	N = floor(LTcount/2.)
2153
2154	!================================================================
2155	! 1. Filling of LT_list
2156	!================================================================
2157	SELECT CASE (dayornight)
2158	CASE ("dayside")
2159	if (abs(LTave-LTmax).le.abs(LTave-LTmin)) then ! LTave is closer to LTmax than to LTmin
2160	do a=1,N
2161	LT_list(a)=LTave+a*abs(LTave-LTmax)/(N+1)
2162	LT_list(N+a)=LTave-a*abs(LTave-LTmax)/(N+1)
2163	enddo
2164	else ! LTave is closer to LTmin than to LTmax
2165	do a=1,N
2166	LT_list(a)=LTave+a*abs(LTave-LTmin)/(N+1)
2167	LT_list(N+a)=LTave-a*abs(LTave-LTmin)/(N+1)
2168	enddo
2169	endif
2170	CASE ("nightside")
2171	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
2172	do a=1,N
2173	LT_list(a)=LTave+a*abs(f_LT(LTave)-f_LT(LTmax))/(N+1)
2174	LT_list(a)=mod(LT_list(a),24.) ! reput the LT in a [0;24[ interval if needed
2175	LT_list(N+a)=LTave-a*abs(f_LT(LTave)-f_LT(LTmax))/(N+1)
2176	LT_list(N+a)=mod(LT_list(N+a),24.)
2177	enddo
2178	else ! LTave is closer to LTmin than to LTmax
2179	do a=1,N
2180	LT_list(a)=LTave+a*abs(f_LT(LTave)-f_LT(LTmin))/(N+1)
2181	LT_list(a)=mod(LT_list(a),24.)
2182	LT_list(N+a)=LTave-a*abs(f_LT(LTave)-f_LT(LTmin))/(N+1)
2183	LT_list(N+a)=mod(LT_list(N+a),24.)
2184	enddo
2185	endif
2186	END SELECT
2187
2188	if (mod(LTcount,2).eq.1) then ! if LTcount is an odd number
2189	LT_list(LTcount)=LTave ! add LTave to the list
2190	endif
2191
2192	!================================================================
2193	! 2. Combination of int_sol_list & LT_list into sol_list
2194	!================================================================
2195	c=1
2196	do a=1,solcount
2197	do b=1,LTcount
2198	sol_list(c)=int_sol_list(a)+(LT_list(b)-lon/15.)/24.
2199	c=c+1
2200	enddo
2201	enddo
2202
2203	end subroutine gen_sol_list
2204
2205	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2206
2207	subroutine status_check(status,error_text)
2208
2209	use netcdf
2210	implicit none
2211	!================================================================
2212	! Arguments:
2213	!================================================================
2214	integer,intent(in) :: status
2215	character(len=256),intent(in) :: error_text
2216
2217	if (status.ne.nf90_noerr) then
2218	write(,)trim(error_text)
2219	write(,)trim(nf90_strerror(status))
2220	stop
2221	endif
2222
2223	end subroutine status_check
2224
2225	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2226
2227	real function LTmod(LT)
2228	!================================================================
2229	! Purpose:
2230	! For night local times management, replace hours
2231	! from a [0;24[ interval to a [-12;12[ interval in which
2232	! midnight = 0 (in order to ensure continuity when comparing
2233	! night local times)
2234	!================================================================
2235	implicit none
2236	real,intent(in) :: LT
2237
2238	LTmod = 2*mod(LT,12.)-LT
2239	return
2240	end function LTmod

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

Download in other formats:

Original Format