Context Navigation

simu_MCS.F90 @ 2525

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

Mars GCM:
Make a cleaner initialization of some variables in simu_MCS.F90

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

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

Download in other formats:

Original Format