source: trunk/LMDZ.MARS/util/simu_MCS.F90 @ 2404

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

Mars GCM:
Add the new version of the MCS observer simulator program simu_MCS.F90, replacing the old version simu_MCS_temp.F90.
The program can now handle (ie, interpolate and bin like the MCS file) any 4D GCM variable, by linking them to
a MCS variable among temp (default), dust and water ice, serving as reference for the binning.
More details in the preface of simu_MCS.F90 and in simu_MCS.def

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