source: trunk/LMDZ.MARS/util/simu_MCS_temp.F90 @ 2202

Last change on this file since 2202 was 2202, checked in by abierjon, 5 years ago

Mars GCM:
Add the observer simulator program simu_MCS_temp.F90 to the LMDZ.MARS/util/ directory, as well as a simu_MCS_temp.def. This program reads a MCS file (binned by Luca Montabone), extracts the daytime and nighttime temperatures from a GCM simulation file, and outputs a netcdf file in the same format as the MCS file.
More details in simu_MCS_temp.def and in the preface of simu_MCS_temp.F90.
AB
File size: 69.2 KB
Line 
1program simu_MCS_temp
2! program for a MCS observer simulator of the temperature data
3! author : Antoine Bierjon, June-December 2019
4! PLEASE GIVE ME YOUR FEEDBACKS ON THIS PROGRAM ! :D
5!
6!===================================================================================================
7!     PREFACE
8!===================================================================================================
9! This program loads Observer's data file (MCS-type, binned by Luca Montabone), serving as a
10! reference, then extract the dayside and nightside temperature data from either :
11!       - a GCM simulation (diagfi.nc) that covers the observation period,
12!       - a GCM stats file (stats.nc) after extending it to cover the observation period.
13!
14! Since the MCS data is binned by 5°-intervals of Ls, the GCM simulation's data is interpolated at
15! the Observer's spatial coordinates, at every GCM sol value in each 5°-interval of Ls, before
16! being averaged to create output bins for each interval. The binning in sols also takes into
17! account the variability of Local Time in the bin and try to represent it too.
18!
19! Eventually, the program outputs a netcdf file, filled with the output bin's data and edited with
20! the same format than the MCS file (and with the same dimensions' declaration order).
21!
22! Minimal requirements :
23!     - the MCS file must contain the variables dtemp,dtimeave,dtimemax,dtimemin,
24!                                               ntemp,ntimeave,ntimemax,ntimemin
25!     - the GCM file must : contain the variable temp (atmospheric temperature) ;
26!                           have the same altitude type as the MCS file ;
27!                           cover completely the NON-BINNED observation period
28!     - if the GCM file is a stats file, please name it with "stats" in it (to be changed in the future)
29!
30! See also NOTA BENE in section 1.3
31! More details can also be found in my internship report (in French) :
32!           /planeto/abierjon/rapportstage2019/Rapport_de_PFE_Master_Bierjon_2019.pdf
33! or directly in my office !
34!===================================================================================================
35
36
37use netcdf
38
39!===================================================================================================
40! 0. Variable declarations
41!===================================================================================================
42
43implicit none ! for no implicitly typed variables
44
45! Files:
46character(len=256) :: gcmfile ! GCM simulation file
47logical :: is_stats = .false. ! to check if the GCM file is a stats.nc or a diagfi.nc file
48character(len=256) :: obsfile ! observation file = MCS/Observer data file
49character(len=256) :: outfile ! output file
50
51! NetCDF stuff
52integer :: status ! NetCDF routines return code
53character (len=256) :: error_text ! to store the error text
54integer :: gcmfid ! NetCDF gcm file ID
55integer :: obsfid ! NetCDF observation file ID
56integer :: outfid ! NetCDF output file ID
57integer :: GCMvarid, OBSvarid, LT_id, outvarid ! to store the ID of a variable
58integer :: lat_dimid_obs,lon_dimid_obs,alt_dimid_obs,time_dimid_obs ! dimensions' ID in OBS and output files
59integer :: lat_dimid_gcm,lon_dimid_gcm,alt_dimid_gcm,time_dimid_gcm ! dimensions' ID in GCM file
60integer :: GCMvarshape(4), OBSvarshape(4), LTshape(3) ! to store a variable's coordinates order
61integer :: nb
62
63real,dimension(:),allocatable :: GCMlon, OBSlon ! longitude in the GCM & the Observer files
64integer GCMlonlen, OBSlonlen ! # of grid points along GCMlon & OBSlon
65
66real,dimension(:),allocatable :: GCMlat, OBSlat ! latitude in the GCM & the Observer files
67integer GCMlatlen, OBSlatlen ! # of grid points along GCMlat & OBSlat
68
69real,dimension(:),allocatable :: GCMalt, OBSalt ! can be geometric heights or pressure
70integer GCMaltlen, OBSaltlen ! # of grid point along GCMalt & OBSalt
71character (len=256) :: units ! to store the units attribute of altitude
72character(len=1) :: GCMalttype, OBSalttype ! altitude coord. type:'z' (altitude, m) 'p' (pressure, Pa)
73
74real,dimension(:),allocatable :: GCMtime, OBSLs ! time in the GCM diagfi (sols) & the Observer files (Ls)
75real,dimension(:),allocatable :: GCMstatstime ! time in the GCM stats file (LT at lon 0°)
76integer GCMtimelen, GCMstatstimelen, OBSLslen ! # of points along GCMtime, GCMstatstime, OBSLs
77real :: starttimeoffset=0 ! offset (in sols) wrt Ls=0 of sol 0 in GCM file
78
79character (len=64) :: GCMvarname, OBSvarname ! to store the name of the variable to retrieve
80real,dimension(:,:,:,:),allocatable :: GCM_var,OBS_var ! the 4D variable extracted from GCM & OBS files
81character (len=64) :: OBSLTave_name, OBSLTmax_name, OBSLTmin_name ! to store the name of the average, max and min of LT (ex : dtimeave, ntimemax...)
82real,dimension(:,:,:),allocatable :: OBSLT, OBSLTmax, OBSLTmin ! 3D variables extracted from obsfile (average, max and min of LT in a bin)
83real :: GCMmiss_val, OBSmiss_val, LTmiss_val ! value to denote non-existant data
84real :: extr_value ! to store the result of the extraction subroutine
85
86real :: OBSdeltaLs ! difference of Ls between each observation bin
87real :: sol, maxsol, minsol ! sol date corresponding to Observed Ls and LT
88integer :: m_maxsol, m_minsol ! indexes of the maximum and minimum GCM sol taken in a bin for interpolation
89
90real,dimension(:,:,:),allocatable :: LTmod_arraymax, LTmod_arraymin ! declaration of the type of the function LTmod_array
91real :: LTmod ! declaration of the type of the function LTmod
92real :: maxLTvar ! maximum variability of LT within one MCS bin
93integer :: LTcount ! nb of LT samples on which the interpolation is performed, for every LT interval
94
95integer :: solcount ! number of GCM sol integer values in one Ls interval
96real,dimension(:),allocatable :: int_sol_list ! list of the integer values of GCM sol
97real,dimension(:),allocatable :: sol_list ! list of the sol values used for the interpolation
98integer :: solerrcount ! nb of GCM missing values during interpolation, which are removed for the binning
99integer :: errcount = 0 ! total number of GCM missing values
100real :: solbinned_value ! to store the extracted value averaged on sol samples, which is finally put in the output bin
101real :: GCMdeltat ! timestep of the GCM simulation (in Martian hours)
102
103real :: lon_val, lat_val, alt_val, Ls_val, LT_val ! where and when the output file is written at
104
105integer :: i,j,k,l ! loop iteration indexes
106integer :: m,m_int,n ! sol binning loops iteration index
107
108real, dimension(:,:,:,:), allocatable :: outvar ! outvar(,,,): 4D array to store the output variable's data
109
110!===================================================================================================
111! 1. Read the observation file (Observer data)
112!===================================================================================================
113!===============================================================================
114! 1.1 Open the Observer data file obsfile
115!===============================================================================
116! Ask the user to give a netcdf observation file
117write(*,*) "-> Enter observation file name :"
118read(*,*) obsfile
119
120status=nf90_open(obsfile,nf90_nowrite,obsfid) ! nowrite mode=the program can only read the opened file
121error_text="Error: could not open file "//trim(obsfile)
122call status_check(status,error_text)
123
124!===============================================================================
125! 1.2 Get grids for dimensions lon,lat,alt,time from the observation file
126!===============================================================================
127! OBS Latitude
128status=nf90_inq_dimid(obsfid,"latitude",lat_dimid_obs)
129error_text="Failed to find Observer latitude dimension"
130call status_check(status,error_text)
131
132status=nf90_inquire_dimension(obsfid,lat_dimid_obs,len=OBSlatlen)
133error_text="Failed to find Observer latitude length"
134call status_check(status,error_text)
135allocate(OBSlat(OBSlatlen))
136
137status=nf90_inq_varid(obsfid,"latitude",OBSvarid)
138error_text="Failed to find Observer latitude ID"
139call status_check(status,error_text)
140
141! Read OBSlat
142status=nf90_get_var(obsfid,OBSvarid,OBSlat)
143error_text="Failed to load OBSlat"
144call status_check(status,error_text)
145
146
147! OBS Longitude
148status=nf90_inq_dimid(obsfid,"longitude",lon_dimid_obs)
149error_text="Failed to find Observer longitude dimension"
150call status_check(status,error_text)
151
152status=nf90_inquire_dimension(obsfid,lon_dimid_obs,len=OBSlonlen)
153error_text="Failed to find Observer longitude length"
154call status_check(status,error_text)
155allocate(OBSlon(OBSlonlen))
156
157status=nf90_inq_varid(obsfid,"longitude",OBSvarid)
158error_text="Failed to find Observer longitude ID"
159call status_check(status,error_text)
160
161! Read GCMlon
162status=nf90_get_var(obsfid,OBSvarid,OBSlon)
163error_text="Failed to load OBSlon"
164call status_check(status,error_text)
165
166
167! OBS Time (Ls)
168status=nf90_inq_dimid(obsfid,"time",time_dimid_obs)
169error_text="Failed to find Observer time (Ls) dimension"
170call status_check(status,error_text)
171
172status=nf90_inquire_dimension(obsfid,time_dimid_obs,len=OBSLslen)
173error_text="Failed to find Observer time (Ls) length"
174call status_check(status,error_text)
175allocate(OBSLs(OBSLslen))
176
177status=nf90_inq_varid(obsfid,"time",OBSvarid)
178error_text="Failed to find Observer time (Ls) ID"
179call status_check(status,error_text)
180
181! Read OBSLs
182status=nf90_get_var(obsfid,OBSvarid,OBSLs)
183error_text="Failed to load OBSLs"
184call status_check(status,error_text)
185
186! Get the observation timestep between bins
187OBSdeltaLs=OBSLs(2)-OBSLs(1)
188
189! OBS Altitude
190status=nf90_inq_dimid(obsfid,"altitude",alt_dimid_obs)
191error_text="Failed to find Observer altitude dimension"
192call status_check(status,error_text)
193
194status=nf90_inquire_dimension(obsfid,alt_dimid_obs,len=OBSaltlen)
195error_text="Failed to find Observer altitude length"
196call status_check(status,error_text)
197allocate(OBSalt(OBSaltlen))
198
199status=nf90_inq_varid(obsfid,"altitude",OBSvarid)
200error_text="Failed to find Observer altitude ID"
201call status_check(status,error_text)
202
203! Read OBSalt
204status=nf90_get_var(obsfid,OBSvarid,OBSalt)
205error_text="Failed to load OBSalt"
206call status_check(status,error_text)
207
208! Check altitude attribute "units" to find out altitude type and compare with the GCM file
209status=nf90_get_att(obsfid,OBSvarid,"units",units)
210error_text="Failed to load Observer altitude units attribute"
211call status_check(status,error_text)
212! an unknown and invisible character is placed just after the unit's
213! characters in the Observer file so we only take the first characters
214! corresponding to the sought unit
215if (trim(units(1:2)).eq."Pa") then
216  units="Pa"
217  OBSalttype='p' ! pressure coordinate
218else if (trim(units(1:2)).eq."m") then
219  units="m"
220  OBSalttype='z' ! altitude coordinate
221else
222  write(*,*)" I do not understand this unit ",trim(units)," for Observer altitude!"
223  stop
224endif
225
226!===============================================================================
227! 1.3 Extraction of the wanted variables from the observation file
228!===============================================================================
229
230! Observer variables to extract :
231!******************** NOTA BENE (cf sections 1.3 and 4)*************************
232! We execute the program a first time with the daytime values, and then a second
233! time with the nighttime values. However, a lot of instructions below are
234! independent of this distinction, hence we execute them only in the first loop,
235! when OBSvarname="dtemp". It also implies a lot of IF (when it's only executed in
236! the first loop) or CASE (when it's executed in both loops but depends on wether
237! it is daytime or nighttime).
238! This should be changed one day with a cleaner version of the code...
239!*******************************************************************************
240OBSvarname="dtemp" ! we begin with daytime temperature
241write(*,*)"" ; write(*,*) "Beginning the 1st loop, on daytime values"; write(*,*)""
242DAY_OR_NIGHT: DO ! (the end of the loop is in section 4.)
243
244  ! Check that the observation file contains that variable
245  status=nf90_inq_varid(obsfid,trim(OBSvarname),OBSvarid)
246  error_text="Failed to find variable "//trim(OBSvarname)//" in "//trim(obsfile)
247  call status_check(status,error_text)
248
249  ! Sanity checks on the variable
250  status=nf90_inquire_variable(obsfid,OBSvarid,ndims=nb,dimids=OBSvarshape)
251  error_text="Failed to obtain information on variable "//trim(OBSvarname)
252  call status_check(status,error_text)
253
254  ! Check that it is a 4D variable
255  if (nb.ne.4) then
256    write(*,*) "Error, expected a 4D (lon-lat-alt-time) variable for ", trim(OBSvarname)
257    stop
258  endif
259  ! Check that its dimensions are indeed lon,lat,alt,time (in the right order)
260  if (OBSvarshape(1).ne.lon_dimid_obs) then
261    write(*,*) "Error, expected first dimension of ", trim(OBSvarname), " to be longitude!"
262    stop
263  endif
264  if (OBSvarshape(2).ne.lat_dimid_obs) then
265    write(*,*) "Error, expected second dimension of ", trim(OBSvarname), " to be latitude!"
266    stop
267  endif
268  if (OBSvarshape(3).ne.alt_dimid_obs) then
269    write(*,*) "Error, expected third dimension of ", trim(OBSvarname), " to be altitude!"
270    stop
271  endif
272  if (OBSvarshape(4).ne.time_dimid_obs) then
273    write(*,*) "Error, expected fourth dimension of ", trim(OBSvarname), " to be time!"
274    stop
275  endif
276
277  !write(*,*)"Dimensions ID in the obsfile are :" 
278  !write(*,*)"lon_dimid=",lon_dimid_obs
279  !write(*,*)"lat_dimid=",lat_dimid_obs
280  !write(*,*)"alt_dimid=",alt_dimid_obs
281  !write(*,*)"time_dimid=",time_dimid_obs
282
283  ! Length allocation for each dimension of the 4D variable
284  allocate(OBS_var(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
285
286  ! Load datasets
287  status=nf90_get_var(obsfid,OBSvarid,OBS_var)
288  error_text="Failed to load "//trim(OBSvarname)//" from the obsfile"
289  call status_check(status,error_text)
290  write(*,*) "Loaded ",trim(OBSvarname)," from the obsfile"
291
292  ! Get OBS_var missing_value attribute
293  status=nf90_get_att(obsfid,OBSvarid,"_FillValue",OBSmiss_val)
294  error_text="Failed to load missing_value attribute"
295  call status_check(status,error_text)
296  write(*,'("  with missing_value attribute : ",1pe12.5)')OBSmiss_val
297
298!===============================================================================
299! 1.4 Extraction of the OBS Local Time (LT)
300!===============================================================================
301  SELECT CASE (OBSvarname)
302  CASE ("dtemp")
303    OBSLTave_name = "dtimeave"
304    OBSLTmax_name = "dtimemax"
305    OBSLTmin_name = "dtimemin"
306  CASE ("ntemp")
307    OBSLTave_name = "ntimeave"
308    OBSLTmax_name = "ntimemax"
309    OBSLTmin_name = "ntimemin"
310  END SELECT
311 
312  status=nf90_inq_varid(obsfid,trim(OBSLTave_name),LT_id)
313  error_text="Failed to find Observer local time ("//trim(OBSLTave_name)//") dimension"
314  call status_check(status,error_text)
315  ! Sanity checks on the variable
316  status=nf90_inquire_variable(obsfid,LT_id,ndims=nb,dimids=LTshape)
317  error_text="Failed to obtain information on variable "//trim(OBSLTave_name)
318  call status_check(status,error_text)
319  ! Check that it is a 3D variable
320  if (nb.ne.3) then
321    write(*,*) "Error, expected a 3D (lon-lat-time) variable for of "//trim(OBSLTave_name)//"!"
322    stop
323  endif
324  ! Check that its dimensions are indeed lon,lat,time (in the right order)
325  if (LTshape(1).ne.lon_dimid_obs) then
326    write(*,*) "Error, expected first dimension of ", trim(OBSLTave_name), " to be longitude!"
327    stop
328  endif
329  if (LTshape(2).ne.lat_dimid_obs) then
330    write(*,*) "Error, expected second dimension of ", trim(OBSLTave_name), " to be latitude!"
331    stop
332  endif
333  if (LTshape(3).ne.time_dimid_obs) then
334    write(*,*) "Error, expected third dimension of ", trim(OBSLTave_name), " to be time!"
335    stop
336  endif
337  ! Length allocation for each dimension of the 3D variable
338  allocate(OBSLT(OBSlonlen,OBSlatlen,OBSLslen))
339  ! Load datasets
340  status=nf90_get_var(obsfid,LT_id,OBSLT)
341  error_text="Failed to load "//trim(OBSLTave_name)//" from the obsfile"
342  call status_check(status,error_text)
343  write(*,*) "Loaded ", trim(OBSLTave_name), " from the obsfile"
344  ! Get LT missing_value attribute
345  status=nf90_get_att(obsfid,LT_id,"_FillValue",LTmiss_val)
346  error_text="Failed to load missing_value attribute"
347  call status_check(status,error_text)
348  write(*,'("  with missing_value attribute : ",1pe12.5)')LTmiss_val
349
350
351  ! Get the maximum values of LT in the OBS bins
352  status=nf90_inq_varid(obsfid,trim(OBSLTmax_name),LT_id)
353  error_text="Failed to find Observer max local time ("//trim(OBSLTmax_name)//") dimension"
354  call status_check(status,error_text)
355  ! Sanity checks on the variable
356  status=nf90_inquire_variable(obsfid,LT_id,ndims=nb,dimids=LTshape)
357  error_text="Failed to obtain information on variable "//trim(OBSLTmax_name)
358  call status_check(status,error_text)
359  ! Check that it is a 3D variable
360  if (nb.ne.3) then
361    write(*,*) "Error, expected a 3D (lon-lat-time) variable for of ", trim(OBSLTmax_name), "!"
362    stop
363  endif
364  ! Check that its dimensions are indeed lon,lat,time (in the right order)
365  if (LTshape(1).ne.lon_dimid_obs) then
366    write(*,*) "Error, expected first dimension of ", trim(OBSLTmax_name), " to be longitude!"
367    stop
368  endif
369  if (LTshape(2).ne.lat_dimid_obs) then
370    write(*,*) "Error, expected second dimension of ", trim(OBSLTmax_name), " to be latitude!"
371    stop
372  endif
373  if (LTshape(3).ne.time_dimid_obs) then
374    write(*,*) "Error, expected third dimension of ", trim(OBSLTmax_name), " to be time!"
375    stop
376  endif
377  ! Length allocation for each dimension of the 3D variable
378  allocate(OBSLTmax(OBSlonlen,OBSlatlen,OBSLslen))
379  ! Load datasets
380  status=nf90_get_var(obsfid,LT_id,OBSLTmax)
381  error_text="Failed to load "//trim(OBSLTmax_name)//" from the obsfile"
382  call status_check(status,error_text)
383  write(*,*) "Loaded ", trim(OBSLTmax_name), " from the obsfile"
384
385
386  ! Get the minimum values of LT in the OBS bins
387  status=nf90_inq_varid(obsfid,trim(OBSLTmin_name),LT_id)
388  error_text="Failed to find Observer min local time ("//trim(OBSLTmin_name)//") dimension"
389  call status_check(status,error_text)
390  ! Sanity checks on the variable
391  status=nf90_inquire_variable(obsfid,LT_id,ndims=nb,dimids=LTshape)
392  error_text="Failed to obtain information on variable "//trim(OBSLTmin_name)
393  call status_check(status,error_text)
394  ! Check that it is a 3D variable
395  if (nb.ne.3) then
396    write(*,*) "Error, expected a 3D (lon-lat-time) variable for of ", trim(OBSLTmin_name), "!"
397    stop
398  endif
399  ! Check that its dimensions are indeed lon,lat,time (in the right order)
400  if (LTshape(1).ne.lon_dimid_obs) then
401    write(*,*) "Error, expected first dimension of ", trim(OBSLTmin_name), " to be longitude!"
402    stop
403  endif
404  if (LTshape(2).ne.lat_dimid_obs) then
405    write(*,*) "Error, expected second dimension of ", trim(OBSLTmin_name), " to be latitude!"
406    stop
407  endif
408  if (LTshape(3).ne.time_dimid_obs) then
409    write(*,*) "Error, expected third dimension of ", trim(OBSLTmin_name), " to be time!"
410    stop
411  endif
412  ! Length allocation for each dimension of the 3D variable
413  allocate(OBSLTmin(OBSlonlen,OBSlatlen,OBSLslen))
414  ! Load datasets
415  status=nf90_get_var(obsfid,LT_id,OBSLTmin)
416  error_text="Failed to load "//trim(OBSLTmin_name)//" from the obsfile"
417  call status_check(status,error_text)
418  write(*,*) "Loaded ", trim(OBSLTmin_name), " from the obsfile"
419
420 
421  SELECT CASE (OBSvarname)
422  CASE ("dtemp")
423    maxLTvar=maxval((OBSLTmax(:,:,:)-OBSLTmin(:,:,:)))
424    write(*,*)"The max daytime LT variability is ",maxLTvar,"hours (within a same bin)"
425  CASE ("ntemp")
426    ! in case of night local times, we replace hours from a [0;24[ interval to a
427    ! [-12;12[ interval in which midnight = 0, with the subroutine LTmod
428    allocate(LTmod_arraymax(OBSlonlen,OBSlatlen,OBSLslen))
429    allocate(LTmod_arraymin(OBSlonlen,OBSlatlen,OBSLslen))
430    do i=1,OBSlonlen
431      do j=1,OBSlatlen
432        do l=1,OBSLslen
433          LTmod_arraymax(i,j,l) = LTmod(OBSLTmax(i,j,l))
434          LTmod_arraymin(i,j,l) = LTmod(OBSLTmin(i,j,l))
435        enddo
436      enddo
437    enddo
438    maxLTvar=maxval(LTmod_arraymax(:,:,:) - LTmod_arraymin(:,:,:))
439    write(*,*)"The max nighttime LT variability is ",maxLTvar,"hours (within a same bin)"
440  END SELECT
441 
442  LTcount=floor(maxLTvar) ! integer nb of hours fitting in the broadest interval of LT
443
444!===================================================================================================
445! 2. Read the input file (GCM simulation)
446!===================================================================================================
447  IF (OBSvarname.EQ."dtemp") THEN ! reading the GCM file is done only in the first loop
448 
449!===============================================================================
450! 2.1 Open the GCM simulation file gcmfile
451!===============================================================================
452    ! Ask the user to give a netcdf input file
453    write(*,*)"";write(*,*) "-> Enter input file name (GCM simulation) :"
454    read(*,*) gcmfile
455    if (index(gcmfile,"stats").ne.0) then ! if the GCM file name contains "stats", we assume it is a stats file
456    ! => A REMPLACER PAR UN TEST SUR LA DIMENSION TIME ? REGARDER DANS AUTRES UTILITAIRES
457      is_stats = .true.
458      write(*,*)"The GCM file is recognized as a stats file."
459    else
460      write(*,*)"The GCM file is recognized as a diagfi/concatnc file."
461    endif
462
463    ! Open GCM file
464    status=nf90_open(gcmfile,nf90_nowrite,gcmfid)
465    ! nowrite mode=the program can only read the opened file
466    error_text="Failed to open datafile "//trim(gcmfile)
467    call status_check(status,error_text)
468
469!===============================================================================
470! 2.2 Get grids for dimensions lon,lat,alt,time from the GCM file
471!===============================================================================
472    ! GCM Latitude
473    status=nf90_inq_dimid(gcmfid,"latitude",lat_dimid_gcm)
474    error_text="Failed to find GCM latitude dimension"
475    call status_check(status,error_text)
476
477    status=nf90_inquire_dimension(gcmfid,lat_dimid_gcm,len=GCMlatlen)
478    error_text="Failed to find GCM latitude length"
479    call status_check(status,error_text)
480    allocate(GCMlat(GCMlatlen))
481
482    status=nf90_inq_varid(gcmfid,"latitude",GCMvarid)
483    error_text="Failed to find GCM latitude ID"
484    call status_check(status,error_text)
485
486    ! Read GCMlat
487    status=nf90_get_var(gcmfid,GCMvarid,GCMlat)
488    error_text="Failed to load GCMlat"
489    call status_check(status,error_text)
490
491
492    ! GCM Longitude
493    status=nf90_inq_dimid(gcmfid,"longitude",lon_dimid_gcm)
494    error_text="Failed to find GCM longitude dimension"
495    call status_check(status,error_text)
496
497    status=nf90_inquire_dimension(gcmfid,lon_dimid_gcm,len=GCMlonlen)
498    error_text="Failed to find GCM longitude length"
499    call status_check(status,error_text)
500    allocate(GCMlon(GCMlonlen))
501
502    status=nf90_inq_varid(gcmfid,"longitude",GCMvarid)
503    error_text="Failed to find GCM longitude ID"
504    call status_check(status,error_text)
505
506    ! Read GCMlon
507    status=nf90_get_var(gcmfid,GCMvarid,GCMlon)
508    error_text="Failed to load GCMlon"
509    call status_check(status,error_text)
510
511
512    ! GCM Time
513    status=nf90_inq_dimid(gcmfid,"Time",time_dimid_gcm)
514    error_text="Failed to find GCM time (sols) dimension"
515    call status_check(status,error_text)
516
517    if (.not.is_stats) then
518      status=nf90_inquire_dimension(gcmfid,time_dimid_gcm,len=GCMtimelen)
519      error_text="Failed to find GCM time length"
520      call status_check(status,error_text)
521      allocate(GCMtime(GCMtimelen))
522    else
523      status=nf90_inquire_dimension(gcmfid,time_dimid_gcm,len=GCMstatstimelen)
524      error_text="Failed to find GCM stats time length"
525      call status_check(status,error_text)
526      allocate(GCMstatstime(GCMstatstimelen))
527    endif
528
529    status=nf90_inq_varid(gcmfid,"Time",GCMvarid)
530    error_text="Failed to find GCM time ID"
531    call status_check(status,error_text)
532
533    ! Read GCMtime
534    if (.not.is_stats) then ! if GCM file is a diagfi, time is in sols at longitude 0°
535      status=nf90_get_var(gcmfid,GCMvarid,GCMtime)
536      error_text="Failed to load GCMtime (sols)"
537      call status_check(status,error_text)
538    else ! if GCM file is a stats, time is in LT at longitude 0°
539      status=nf90_get_var(gcmfid,GCMvarid,GCMstatstime)
540      error_text="Failed to load GCMstatstime (LT at lon 0)"
541      call status_check(status,error_text)
542    endif
543
544    ! Simulation time offset management
545    if (.not.is_stats) then
546      write(*,*) "Beginning date of the simulation file?"
547      write(*,*) "(i.e. number of sols since Ls=0 at the Time=0.0 in the GCM file)"
548      read(*,*) starttimeoffset
549      ! Add the offset to GCMtime(:)
550      GCMtime(:)=GCMtime(:)+starttimeoffset
551    endif
552
553    ! Get the timestep of the simulation (in Martian hours)
554    if (.not.is_stats) then
555      GCMdeltat = (GCMtime(2)-GCMtime(1))*24.
556    else
557      GCMdeltat = (GCMstatstime(2)-GCMstatstime(1))
558    endif
559    write(*,*)"GCM simulation's timestep is ", GCMdeltat, " hours."
560
561    ! Check of temporal coherence between gcmfile & obsfile
562    call ls2sol(OBSLs(OBSLslen),maxsol) ! maximum date considered
563    call ls2sol(OBSLs(1),minsol) ! minimum date considered
564    if (.not.is_stats) then ! if it is a diagfi, we check the time coherence between the 2 files
565      if ((maxsol.gt.maxval(GCMtime)).or.(minsol.lt.minval(GCMtime))) then
566        write(*,*)"Error : obsfile temporal bounds exceed the GCM simulation bounds."
567        write(*,*)"Please use a GCM file whose time interval contains the observation period."
568        stop
569      else
570        write(*,*)"Both files are temporally coherent. The program continues..."
571      endif
572
573    else ! if it is a stats, we create the array GCMtime array (in sols) covering the observation period
574         ! and filled with the mean GCM day stored in stats.nc
575
576      GCMtimelen = ((ceiling(maxsol)-floor(minsol)+1)+2) ! we add 2 days in the beginning and the end
577                                                         ! to be sure we cover the observation period
578      allocate(GCMtime(GCMstatstimelen * GCMtimelen))
579      do l=1,GCMtimelen
580        do m=1,GCMstatstimelen
581          GCMtime(m+(l-1)*GCMstatstimelen) = (floor(minsol)-1) + (l-1) + GCMstatstime(m)/24.
582        enddo
583      enddo
584      GCMtimelen = GCMstatstimelen * GCMtimelen
585      write(*,*)"GCMtime has been created from the stats.nc time and the observation period. The program continues..."
586    endif
587
588
589    ! GCM Altitude
590    status=nf90_inq_dimid(gcmfid,"altitude",alt_dimid_gcm)
591    error_text="Failed to find GCM altitude dimension"
592    call status_check(status,error_text)
593
594    status=nf90_inquire_dimension(gcmfid,alt_dimid_gcm,len=GCMaltlen)
595    error_text="Failed to find GCM altitude length"
596    call status_check(status,error_text)
597    allocate(GCMalt(GCMaltlen))
598
599    status=nf90_inq_varid(gcmfid,"altitude",GCMvarid)
600    error_text="Failed to find GCM altitude ID"
601    call status_check(status,error_text)
602
603    ! Read GCMalt
604    status=nf90_get_var(gcmfid,GCMvarid,GCMalt)
605    error_text="Failed to load GCMalt"
606    call status_check(status,error_text)
607
608    ! Check altitude attribute "units" to find out altitude type
609    status=nf90_get_att(gcmfid,GCMvarid,"units",units)
610    error_text="Failed to load GCM altitude units attribute"
611    call status_check(status,error_text)
612    if (trim(units).eq."Pa") then
613      GCMalttype='p' ! pressure coordinate
614    else if (trim(units).eq."m") then
615      GCMalttype='z' ! altitude coordinate
616    else
617      write(*,*)"I do not understand this unit ",trim(units)," for GCM altitude!"
618      if (OBSalttype.eq.'p') then
619        write(*,*)"Please use zrecast to put the altitude in the same type as the MCS file (pressure in Pa)"
620      else if (OBSalttype.eq.'z') then
621        write(*,*)"Please use zrecast to put the altitude in the same type as the MCS file (altitude in m)"
622      endif
623      stop
624    endif
625    if(OBSalttype.ne.GCMalttype) then
626      write(*,*)"Observer altitude type (", OBSalttype,") and ", &
627                "GCM altitude type (",GCMalttype,") don't match!"
628      stop
629    endif
630!===============================================================================
631! 2.3 Extraction of the wanted variables (the ones to compare with observer data)
632!===============================================================================
633    ! GCM variable to extract for both MCS dtemp & ntemp is : temp
634    GCMvarname="temp" ! temperature
635
636    ! Check that the GCM file contains that variable
637    status=nf90_inq_varid(gcmfid,trim(GCMvarname),GCMvarid)
638    error_text="Failed to find variable "//trim(GCMvarname)//" in "//trim(gcmfile)
639    call status_check(status,error_text)
640
641    ! Sanity checks on the variable
642    status=nf90_inquire_variable(gcmfid,GCMvarid,ndims=nb,dimids=GCMvarshape)
643    error_text="Failed to obtain information on variable "//trim(GCMvarname)
644    call status_check(status,error_text)
645
646    ! Check that it is a 4D variable
647    if (nb.ne.4) then
648      write(*,*) "Error, expected a 4D (lon-lat-alt-time) variable for ", trim(GCMvarname)
649      stop
650    endif
651    ! Check that its dimensions are indeed lon,lat,alt,time (in the right order)
652    if (GCMvarshape(1).ne.lon_dimid_gcm) then
653      write(*,*) "Error, expected first dimension of ", trim(GCMvarname), " to be longitude!"
654      stop
655    endif
656    if (GCMvarshape(2).ne.lat_dimid_gcm) then
657      write(*,*) "Error, expected second dimension of ", trim(GCMvarname), " to be latitude!"
658      stop
659    endif
660    if (GCMvarshape(3).ne.alt_dimid_gcm) then
661      write(*,*) "Error, expected third dimension of ", trim(GCMvarname), " to be altitude!"
662      stop
663    endif
664    if (GCMvarshape(4).ne.time_dimid_gcm) then
665      write(*,*) "Error, expected fourth dimension of ", trim(GCMvarname), " to be time!"
666      stop
667    endif
668
669    ! Length allocation for each dimension of the 4D variable
670    allocate(GCM_var(GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen))
671
672    ! Load datasets
673    if (.not.is_stats) then
674      status=nf90_get_var(gcmfid,GCMvarid,GCM_var)
675      error_text="Failed to load "//trim(GCMvarname)
676      call status_check(status,error_text)
677    else
678      ! if it is a stats file, we load only the first sol, and then copy it to all the other sols
679      status=nf90_get_var(gcmfid,GCMvarid,GCM_var(:,:,:,1:GCMstatstimelen))
680      error_text="Failed to load "//trim(GCMvarname)
681      call status_check(status,error_text)
682    !  write(*,*)"GCMstatstimelen = ", GCMstatstimelen
683    !  write(*,*)"GCMtimelen = ", GCMtimelen
684      do l=(GCMstatstimelen+1),GCMtimelen
685        if (modulo(l,GCMstatstimelen).ne.0) then
686          GCM_var(:,:,:,l) = GCM_var(:,:,:,modulo(l,GCMstatstimelen))
687        else ! if l is a multiple of GCMstatstimelen, since index modulo(l,GCMstatstimelen)=0 doesn't exist,
688             ! we make a special case
689          GCM_var(:,:,:,l) = GCM_var(:,:,:,GCMstatstimelen)
690        endif
691      enddo
692    endif
693    write(*,*) "Loaded ",trim(GCMvarname), " from the GCM file"
694
695    ! Get dataset's missing_value attribute
696    status=nf90_get_att(gcmfid,GCMvarid,"missing_value",GCMmiss_val)
697    error_text="Failed to load missing_value attribute"
698    call status_check(status,error_text)
699    write(*,'("  with missing_value attribute : ",1pe12.5)')GCMmiss_val
700
701  ENDIF ! end of IF (OBSvarname.EQ."dtemp")
702 
703!===================================================================================================
704! 3. Output file
705!===================================================================================================
706!===============================================================================
707! 3.1 Create the output file & initialize the coordinates
708!===============================================================================
709  IF (OBSvarname.EQ."dtemp") THEN ! creating the output file is done only in the first loop
710    ! Name of the outfile
711    if (.not.is_stats) then
712      outfile=obsfile(1:index(obsfile, ".nc")-1)//"_GCMdiagfi.nc"
713    else
714      outfile=obsfile(1:index(obsfile, ".nc")-1)//"_GCMstats.nc"
715    endif
716
717    ! Creation of the outfile
718    status=nf90_create(outfile,nf90_clobber,outfid) ! NB: clobber mode=overwrite any dataset with the same file name !
719    error_text="Error: could not create file "//trim(outfile)
720    call status_check(status,error_text)
721    write(*,*)"";write(*,*)"-> Output file is: ",trim(outfile)
722   
723    ! Creation of the dimensions
724    call inidim(outfid,OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen,OBSlon,OBSlat,OBSalt,OBSLs,OBSalttype,&
725                lon_dimid_obs,lat_dimid_obs,alt_dimid_obs,time_dimid_obs)
726
727    !write(*,*)"Dimensions ID in the outfile are :"           
728    !write(*,*)"lon_dimid=",lon_dimid_obs
729    !write(*,*)"lat_dimid=",lat_dimid_obs
730    !write(*,*)"alt_dimid=",alt_dimid_obs
731    !write(*,*)"time_dimid=",time_dimid_obs
732   
733  ENDIF ! end of IF (OBSvarname.EQ."dtemp")
734 
735!===============================================================================
736! 3.2 Create the outfile variables (other variables than the coordinates)
737!===============================================================================
738  ! Switch to netcdf define mode
739  status=nf90_redef(outfid)
740  error_text="Error: could not switch to define mode in the outfile"
741  call status_check(status,error_text)
742
743  ! Insert the definition of the variables
744  status=nf90_def_var(outfid,trim(OBSvarname),nf90_float,OBSvarshape,outvarid)
745  error_text="Error: could not define the variable "//trim(OBSvarname)//" in the outfile"
746  call status_check(status,error_text)
747
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  !! Temperature
754  SELECT CASE (OBSvarname)
755  CASE ("dtemp")
756    status=nf90_put_att(outfid,outvarid,"long_name","Temperature - day side (interpolated from GCM)")
757  CASE ("ntemp")
758    status=nf90_put_att(outfid,outvarid,"long_name","Temperature - night side (interpolated from GCM)")
759  END SELECT
760  error_text="Error: could not put the long_name attribute for the variable "//trim(OBSvarname)//" in the outfile"
761  call status_check(status,error_text)
762 
763  status=nf90_put_att(outfid,outvarid,"units","K")
764  error_text="Error: could not put the units attribute for the variable "//trim(OBSvarname)//" in the outfile"
765  call status_check(status,error_text)
766 
767  status=nf90_put_att(outfid,outvarid,"_FillValue",OBSmiss_val)
768  error_text="Error: could not put the _FillValue attribute for the variable "//trim(OBSvarname)//" in the outfile"
769  call status_check(status,error_text)
770
771  write(*,*)trim(OBSvarname)," has been created in the outfile"
772  write(*,'("  with missing_value attribute : ",1pe12.5)')OBSmiss_val
773
774  !! Local Time (average in bin)
775  SELECT CASE (OBSvarname)
776  CASE ("dtemp")
777    status=nf90_put_att(outfid,LT_id,"long_name","Local Time, evaluated as average local time in bin - day side [6h, 18h]")
778  CASE ("ntemp")
779    status=nf90_put_att(outfid,LT_id,"long_name","Local Time, evaluated as average local time in bin - night side [18h, 6h]")
780  END SELECT
781  error_text="Error: could not put the long_name attribute for the variable "//trim(OBSLTave_name)//" in the outfile"
782  call status_check(status,error_text)
783 
784  status=nf90_put_att(outfid,LT_id,"units","hours")
785  error_text="Error: could not put the units attribute for the variable "//trim(OBSLTave_name)//" in the outfile"
786  call status_check(status,error_text)
787 
788  status=nf90_put_att(outfid,LT_id,"_FillValue",LTmiss_val)
789  error_text="Error: could not put the _FillValue attribute for the variable "//trim(OBSLTave_name)//" in the outfile"
790  call status_check(status,error_text)
791
792  write(*,*)"Local Time (", trim(OBSLTave_name), ") has been created in the outfile"
793  write(*,'("  with missing_value attribute : ",1pe12.5)')OBSmiss_val
794
795  ! End the netcdf define mode (and thus enter the "data writing" mode)
796  status=nf90_enddef(outfid)
797  error_text="Error: could not close the define mode of the outfile"
798  call status_check(status,error_text)
799
800!===============================================================================
801! 3.3 Fill the output file
802!===============================================================================
803  allocate(outvar(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
804
805!!===============================================================
806!! 3.3.1 Do some checks
807!!===============================================================
808  Ls: do l=1,OBSLslen ! loop on all observed Solar longitudes
809    Ls_val=OBSLs(l) ! Ls_val=center of the output bin
810    if ((Ls_val.lt.0.).or.(Ls_val.gt.360.)) then
811      write(*,*) "Unexpected value for OBSLs: ",Ls_val
812      stop
813    endif
814
815    ! Convert the Ls bin into a sol interval on which the binning is done :
816    !-> get the index of the maximum value of GCM sol (m_maxsol) that is lower than Ls bin's superior bound (maxsol)
817    call ls2sol(Ls_val+OBSdeltaLs/2.,maxsol)
818    m_maxsol=1
819    do while (((GCMtime(m_maxsol)+0.5).lt.maxsol).AND.(m_maxsol.le.(GCMtimelen-1)))
820  !! The +0.5 is there to take into account the whole planet (lon=[-180°;180°]) and not just the lon=0° from the GCM
821      m_maxsol=m_maxsol+1
822    enddo
823    !-> get the index of the minimum value of GCM sol (m_minsol) that is greater than Ls bin's inferior bound (minsol)
824    call ls2sol(Ls_val-OBSdeltaLs/2.,minsol)
825    m_minsol=1
826    do while (((GCMtime(m_minsol)-0.5).le.minsol).AND.(m_minsol.le.(GCMtimelen-1)))
827  !! Same comment for the -0.5
828      m_minsol=m_minsol+1
829    enddo
830    if (m_minsol.gt.m_maxsol) then
831      write(*,*) "No value in gcmfile between sol=",minsol," and sol=",maxsol," (Ls=",Ls_val,"°)"
832      ! Write a missing_value to output
833      outvar(:,:,:,l)=OBSmiss_val
834      CYCLE Ls ! go directly to the next Ls
835    endif 
836    ! Get all the integer values of GCM sols that fit in this interval
837    solcount=floor(GCMtime(m_maxsol))-ceiling(GCMtime(m_minsol))+1 ! sols that are not fully in the interval are not counted
838    allocate(int_sol_list(solcount))
839    m_int=1
840  !  write(*,*) "GCMminsol=", GCMtime(m_minsol)
841  !  write(*,*)"GCMmaxsol=", GCMtime(m_maxsol)
842    do m=m_minsol,m_maxsol ! loop on all GCM sols of the bin
843      sol=GCMtime(m)
844      if (sol.eq.floor(sol)) then
845  !      write(*,*)"sol=",sol
846        int_sol_list(m_int)=sol
847        m_int=m_int+1
848      endif
849    enddo
850  !  write(*,*)"int_sol_list=",int_sol_list
851
852    latitude: do j=1,OBSlatlen ! loop on all observed latitudes
853      lat_val=OBSlat(j)
854      if ((lat_val.lt.-90.).or.(lat_val.gt.90.)) then
855        write(*,*) "Unexpected value for OBSlat: ",lat_val
856        stop
857      endif
858
859      longitude: do i=1,OBSlonlen ! loop on all observed longitudes
860        lon_val=OBSlon(i)
861        if ((lon_val.lt.-360.).or.(lon_val.gt.360.)) then
862          write(*,*) "Unexpected value for lon_val: ",lon_val
863          stop
864        endif
865        ! We want lon_val in [-180:180] for the subroutine extraction
866        if (lon_val.lt.-180.) lon_val=lon_val+360.
867        if (lon_val.gt.180.) lon_val=lon_val-360.
868
869        LT_val=OBSLT(i,j,l) ! find the Observer corresponding LT value at lon_val, lat_val, Ls_val
870       
871        if ((LT_val.lt.0.).or.(LT_val.gt.24.)) then
872          if (LT_val.eq.LTmiss_val) then
873  !            write(*,*) "Missing value in obsfile for LT_val"
874            ! Write a missing_value to output
875            outvar(i,j,:,l)=OBSmiss_val
876            CYCLE longitude ! go directly to the next longitude
877          else
878            write(*,*) "Unexpected value for LT_val: ",LT_val
879            stop
880          endif
881        endif
882        ! Generate the sol list for the interpolation
883        allocate(sol_list(solcount*LTcount))
884        call gen_sol_list(solcount,int_sol_list,LTcount,LT_val,OBSLTmax(i,j,l),OBSLTmin(i,j,l),lon_val,LTmod,OBSvarname,&
885                          sol_list)
886
887        altitude: do k=1,OBSaltlen ! loop on all observed altitudes
888          alt_val=OBSalt(k)
889          if (OBS_var(i,j,k,l).eq.OBSmiss_val) then
890  !          write(*,*) "Missing value in obsfile for ",OBSvarname
891            ! Write a missing_value to output
892            outvar(i,j,k,l)=OBSmiss_val
893            CYCLE altitude ! go directly to the next altitude
894          endif
895
896!!===============================================================
897!! 3.3.2 Compute GCM sol date corresponding to Observer Ls
898!!       (via m_(min/max)sol) and LT (via OBSLT(min/max))
899!!===============================================================       
900          solerrcount=0
901          solbinned_value=0
902          sol_bin: do n=1,solcount*LTcount ! loop on all GCM sols of the bin
903            sol=sol_list(n)
904  !          write(*,*)"sol=",sol
905!!=============================================================== 
906!! 3.3.3 Do the interpolation and binning for the given location
907!!===============================================================
908            call extraction(lon_val,lat_val,alt_val,sol,&
909                            GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen,&
910                            GCMlon,GCMlat,GCMalt,GCMtime,&
911                            GCM_var,GCMmiss_val,GCMalttype,GCMvarname,extr_value)
912
913            if (extr_value.eq.GCMmiss_val) then
914  !            write(*,*) "Missing value in gcmfile at lon=",lon_val,"; lat=",lat_val,"; alt=",alt_val,"; sol=",sol
915              solerrcount=solerrcount+1 
916              CYCLE sol_bin ! go directly to the next GCM sol of the bin
917            endif
918            solbinned_value=solbinned_value+extr_value
919
920          enddo sol_bin ! end loop on all GCM sols of the bin
921          if ((solcount*LTcount-solerrcount).ne.0) then
922            solbinned_value=solbinned_value/(solcount*LTcount-solerrcount)
923          else
924  !          write(*,*)"No GCM value in this sol bin"
925            solbinned_value=OBSmiss_val
926          endif
927          ! Write value to output
928          outvar(i,j,k,l)=solbinned_value
929
930          errcount=errcount+solerrcount
931
932        enddo altitude ! end loop on observed altitudes
933
934        deallocate(sol_list)
935
936      enddo longitude ! end loop on observed longitudes
937    enddo latitude ! end loop on observed latitudes
938
939    deallocate(int_sol_list)
940
941  enddo Ls ! end loop on observed Solar longitudes
942
943  !write(*,*)"Nb of GCM missing values :",errcount
944
945!!=============================================================== 
946!! 3.3.4 Write the data in the netcdf output file
947!!===============================================================
948  status = nf90_put_var(outfid, outvarid, outvar)
949  error_text="Error: could not write "//trim(OBSvarname)//" data in the outfile"
950  call status_check(status,error_text)
951
952  status = nf90_put_var(outfid, LT_id, OBSLT) ! write the MCS d/ntimeave as the output Local Time
953  error_text="Error: could not write Local Time data in the outfile"
954  call status_check(status,error_text)
955
956!!============================================================================== 
957!! 4. End of the Day/Night loop
958!!==============================================================================
959  IF (OBSvarname.EQ."dtemp") THEN
960    ! this is the end of the first loop (on daytime values)
961    ! and we still have to loop on nighttime values
962    OBSvarname="ntemp"
963    deallocate(OBS_var)
964    deallocate(OBSLT)
965    deallocate(OBSLTmax)
966    deallocate(OBSLTmin)
967    deallocate(outvar)
968    write(*,*)"";write(*,*)""; write(*,*) "Beginning the 2nd loop, on nighttime values" ; write(*,*)""
969    CYCLE DAY_OR_NIGHT
970  ELSE ! i.e. OBSvarname="ntemp"
971    ! this is the end of the second loop (on nighttime values)
972    ! and thus the end of the program
973    EXIT DAY_OR_NIGHT
974  ENDIF
975ENDDO DAY_OR_NIGHT ! end of the day/night loop that begins in section 1.3
976
977!===============================================================================
978! 5. Close output file
979!===============================================================================
980status = nf90_close(outfid)
981error_text="Error: could not close file "//trim(outfile)
982call status_check(status,error_text)
983
984write(*,*)"";write(*,*)"-> Program completed!"
985
986end program simu_MCS_temp
987
988
989!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
990!===================================================================================================
991! Subroutines
992!===================================================================================================
993
994subroutine extraction(lon,lat,alt,sol,&
995                  lonlen,latlen,altlen,timelen,&
996                  longitude,latitude,altitude,time,&
997                  field,missing_value,alttype,varname,value)
998
999implicit none
1000!================================================================
1001! Arguments:
1002!================================================================
1003real,intent(in) :: lon  ! sought longitude (deg, in [-180:180])
1004real,intent(in) :: lat  ! sought latitude (deg, in [-90:90])
1005real,intent(in) :: alt  ! sought altitude (m or Pa)
1006real,intent(in) :: sol  ! sought date (sols)
1007integer,intent(in) :: lonlen
1008integer,intent(in) :: latlen
1009integer,intent(in) :: altlen
1010integer,intent(in) :: timelen
1011real,intent(in) :: longitude(lonlen)
1012real,intent(in) :: latitude(latlen)
1013real,intent(in) :: altitude(altlen)
1014real,intent(in) :: time(timelen)
1015real,intent(in) :: field(lonlen,latlen,altlen,timelen)
1016real,intent(in) :: missing_value ! default value for "no data"
1017character,intent(in) :: alttype ! altitude coord. type:'z' (altitude, m)
1018                                !                      'p' (pressure, Pa)
1019character(len=*),intent(in) :: varname ! variable name (in GCM file)
1020real,intent(out) :: value
1021
1022!================================================================
1023! Local variables:
1024!================================================================
1025real,save :: prev_lon=-99 ! previous value of 'lon' routine was called with
1026real,save :: prev_lat=-99 ! previous value of 'lat' routine was called with
1027real,save :: prev_alt=-9.e20 ! ! previous value of 'alt'
1028real,save :: prev_sol=-99 ! previous value of 'sol' routine was called with
1029
1030! encompasing indexes:
1031integer,save :: ilon_inf=-1,ilon_sup=-1 ! along longitude
1032integer,save :: ilat_inf=-1,ilat_sup=-1 ! along latitude
1033integer,save :: ialt_inf=-1,ialt_sup=-1 ! along altitude
1034integer,save :: itim_inf=-1,itim_sup=-1 ! along time
1035
1036! intermediate interpolated values
1037real :: t_interp(2,2,2) ! after time interpolation
1038real :: zt_interp(2,2) ! after altitude interpolation
1039real :: yzt_interp(2) ! after latitude interpolation
1040real :: coeff ! interpolation coefficient
1041
1042integer :: i
1043
1044! By default, set value to missing_value
1045value=missing_value
1046
1047!================================================================
1048! 1. Find encompassing indexes
1049!================================================================
1050if (lon.ne.prev_lon) then
1051  do i=1,lonlen-1
1052    if (longitude(i).le.lon) then
1053      ilon_inf=i
1054    else
1055      exit
1056    endif
1057  enddo
1058  ilon_sup=ilon_inf+1
1059endif ! of if (lon.ne.prev_lon)
1060!write(*,*) 'ilon_inf=',ilon_inf," longitude(ilon_inf)=",longitude(ilon_inf)
1061
1062if (lat.ne.prev_lat) then
1063  ! recall that latitudes start from north pole
1064  do i=1,latlen-1
1065    if (latitude(i).ge.lat) then
1066      ilat_inf=i
1067    else
1068      exit
1069    endif
1070  enddo
1071  ilat_sup=ilat_inf+1
1072endif ! of if (lat.ne.prev_lat)
1073!write(*,*) 'ilat_inf=',ilat_inf," latitude(ilat_inf)=",latitude(ilat_inf)
1074
1075if (alt.ne.prev_alt) then
1076  if (alttype.eq.'p') then ! pressures are ordered from max to min
1077    !handle special case for alt not in the altitude(1:altlen) interval
1078    if ((alt.gt.altitude(1)).or.(alt.lt.altitude(altlen))) then
1079      ialt_inf=-1
1080      ialt_sup=-1
1081      ! return to main program (with value=missing_value)
1082!      write(*,*)"Problem in extraction : GCM alt p"
1083      return
1084    else ! general case
1085      do i=1,altlen-1
1086        if (altitude(i).ge.alt) then
1087          ialt_inf=i
1088        else
1089          exit
1090        endif
1091      enddo
1092      ialt_sup=ialt_inf+1
1093    endif ! of if ((alt.gt.altitude(1)).or.(alt.lt.altitude(altlen)))
1094  else ! altitudes (m) are ordered from min to max
1095    !handle special case for alt not in the altitude(1:altlen) interval
1096    if ((alt.lt.altitude(1)).or.(alt.gt.altitude(altlen))) then
1097      ialt_inf=-1
1098      ialt_sup=-1
1099      ! return to main program (with value=missing_value)
1100!      write(*,*)"Problem in extraction : GCM alt z"
1101      return
1102    else ! general case
1103      do i=1,altlen-1
1104        if (altitude(i).le.alt) then
1105          ialt_inf=i
1106        else
1107          exit
1108        endif
1109      enddo
1110      ialt_sup=ialt_inf+1
1111    endif ! of if ((alt.lt.altitude(1)).or.(alt.gt.altitude(altlen)))
1112  endif ! of if (alttype.eq.'p')
1113endif ! of if (alt.ne.prev_alt)
1114!write(*,*) 'ialt_inf=',ialt_inf," altitude(ialt_inf)=",altitude(ialt_inf)
1115
1116if (sol.ne.prev_sol) then
1117  !handle special case for sol not in the time(1):time(timenlen) interval
1118  if ((sol.lt.time(1)).or.(sol.gt.time(timelen))) then
1119    itim_inf=-1
1120    itim_sup=-1
1121    ! return to main program (with value=missing_value)
1122!    write(*,*)"Problem in extraction : GCM sol"
1123    return
1124  else ! general case
1125    do i=1,timelen-1
1126      if (time(i).le.sol) then
1127        itim_inf=i
1128      else
1129        exit
1130      endif
1131    enddo
1132    itim_sup=itim_inf+1
1133  endif ! of if ((sol.lt.time(1)).or.(sol.gt.time(timenlen)))
1134endif ! of if (sol.ne.prev_sol)
1135!write(*,*) 'itim_inf=',itim_inf," time(itim_inf)=",time(itim_inf)
1136!write(*,*) 'itim_sup=',itim_sup," time(itim_sup)=",time(itim_sup)
1137
1138!================================================================
1139! 2. Interpolate
1140!================================================================
1141! check that there are no "missing_value" in the field() elements we need
1142! otherwise return to main program (with value=missing_value)
1143if (field(ilon_inf,ilat_inf,ialt_inf,itim_inf).eq.missing_value) then
1144!  write(*,*)"Error 1 in extraction"
1145  return
1146 endif
1147if (field(ilon_inf,ilat_inf,ialt_inf,itim_sup).eq.missing_value) then
1148!  write(*,*)"Error 2 in extraction"
1149  return
1150 endif
1151if (field(ilon_inf,ilat_inf,ialt_sup,itim_inf).eq.missing_value) then
1152!  write(*,*)"Error 3 in extraction"
1153  return
1154 endif
1155if (field(ilon_inf,ilat_inf,ialt_sup,itim_sup).eq.missing_value) then
1156!  write(*,*)"Error 4 in extraction"
1157  return
1158 endif
1159if (field(ilon_inf,ilat_sup,ialt_inf,itim_inf).eq.missing_value) then
1160!  write(*,*)"Error 5 in extraction"
1161  return
1162 endif
1163if (field(ilon_inf,ilat_sup,ialt_inf,itim_sup).eq.missing_value) then
1164!  write(*,*)"Error 6 in extraction"
1165  return
1166 endif
1167if (field(ilon_inf,ilat_sup,ialt_sup,itim_inf).eq.missing_value) then
1168!  write(*,*)"Error 7 in extraction"
1169  return
1170 endif
1171if (field(ilon_inf,ilat_sup,ialt_sup,itim_sup).eq.missing_value) then
1172!  write(*,*)"Error 8 in extraction"
1173  return
1174 endif
1175if (field(ilon_sup,ilat_inf,ialt_inf,itim_inf).eq.missing_value) then
1176!  write(*,*)"Error 9 in extraction"
1177  return
1178 endif
1179if (field(ilon_sup,ilat_inf,ialt_inf,itim_sup).eq.missing_value) then
1180!  write(*,*)"Error 10 in extraction"
1181  return
1182 endif
1183if (field(ilon_sup,ilat_inf,ialt_sup,itim_inf).eq.missing_value) then
1184!  write(*,*)"Error 11 in extraction"
1185  return
1186 endif
1187if (field(ilon_sup,ilat_inf,ialt_sup,itim_sup).eq.missing_value) then
1188!  write(*,*)"Error 12 in extraction"
1189  return
1190 endif
1191if (field(ilon_sup,ilat_sup,ialt_inf,itim_inf).eq.missing_value) then
1192!  write(*,*)"Error 13 in extraction"
1193  return
1194 endif
1195if (field(ilon_sup,ilat_sup,ialt_inf,itim_sup).eq.missing_value) then
1196!  write(*,*)"Error 14 in extraction"
1197  return
1198 endif
1199if (field(ilon_sup,ilat_sup,ialt_sup,itim_inf).eq.missing_value) then
1200!  write(*,*)"Error 15 in extraction"
1201  return
1202 endif
1203if (field(ilon_sup,ilat_sup,ialt_sup,itim_sup).eq.missing_value) then
1204!  write(*,*)"Error 16 in extraction"
1205  return
1206 endif
1207
1208!================================================================
1209! 2.1 Linear interpolation in time
1210!================================================================
1211coeff=(sol-time(itim_inf))/(time(itim_sup)-time(itim_inf))
1212t_interp(1,1,1)=field(ilon_inf,ilat_inf,ialt_inf,itim_inf)+ &
1213                coeff*(field(ilon_inf,ilat_inf,ialt_inf,itim_sup)- &
1214                       field(ilon_inf,ilat_inf,ialt_inf,itim_inf))
1215t_interp(1,1,2)=field(ilon_inf,ilat_inf,ialt_sup,itim_inf)+ &
1216                coeff*(field(ilon_inf,ilat_inf,ialt_sup,itim_sup)- &
1217                       field(ilon_inf,ilat_inf,ialt_sup,itim_inf))
1218t_interp(1,2,1)=field(ilon_inf,ilat_sup,ialt_inf,itim_inf)+ &
1219                coeff*(field(ilon_inf,ilat_sup,ialt_inf,itim_sup)- &
1220                       field(ilon_inf,ilat_sup,ialt_inf,itim_inf))
1221t_interp(1,2,2)=field(ilon_inf,ilat_sup,ialt_sup,itim_inf)+ &
1222                coeff*(field(ilon_inf,ilat_sup,ialt_sup,itim_sup)- &
1223                       field(ilon_inf,ilat_sup,ialt_sup,itim_inf))
1224t_interp(2,1,1)=field(ilon_sup,ilat_inf,ialt_inf,itim_inf)+ &
1225                coeff*(field(ilon_sup,ilat_inf,ialt_inf,itim_sup)- &
1226                       field(ilon_sup,ilat_inf,ialt_inf,itim_inf))
1227t_interp(2,1,2)=field(ilon_sup,ilat_inf,ialt_sup,itim_inf)+ &
1228                coeff*(field(ilon_sup,ilat_inf,ialt_sup,itim_sup)- &
1229                       field(ilon_sup,ilat_inf,ialt_sup,itim_inf))
1230t_interp(2,2,1)=field(ilon_sup,ilat_sup,ialt_inf,itim_inf)+ &
1231                coeff*(field(ilon_sup,ilat_sup,ialt_inf,itim_sup)- &
1232                       field(ilon_sup,ilat_sup,ialt_inf,itim_inf))
1233t_interp(2,2,2)=field(ilon_sup,ilat_sup,ialt_sup,itim_inf)+ &
1234                coeff*(field(ilon_sup,ilat_sup,ialt_sup,itim_sup)- &
1235                       field(ilon_sup,ilat_sup,ialt_sup,itim_inf))
1236
1237!================================================================
1238! 2.2 Vertical interpolation
1239!================================================================
1240if (((varname=='rho').or.(varname=='pressure')).and.(alttype=='z')) then
1241  ! do the interpolation on the log of the quantity
1242  coeff=(alt-altitude(ialt_inf))/(altitude(ialt_sup)-altitude(ialt_inf))
1243  zt_interp(1,1)=log(t_interp(1,1,1))+coeff* &
1244                             (log(t_interp(1,1,2))-log(t_interp(1,1,1)))
1245  zt_interp(1,2)=log(t_interp(1,2,1))+coeff* &
1246                             (log(t_interp(1,2,2))-log(t_interp(1,2,1)))
1247  zt_interp(2,1)=log(t_interp(2,1,1))+coeff* &
1248                             (log(t_interp(2,1,2))-log(t_interp(2,1,1)))
1249  zt_interp(2,2)=log(t_interp(2,2,1))+coeff* &
1250                             (log(t_interp(2,2,2))-log(t_interp(2,2,1)))
1251  zt_interp(1:2,1:2)=exp(zt_interp(1:2,1:2))
1252else ! general case
1253  coeff=(alt-altitude(ialt_inf))/(altitude(ialt_sup)-altitude(ialt_inf))
1254  zt_interp(1,1)=t_interp(1,1,1)+coeff*(t_interp(1,1,2)-t_interp(1,1,1))
1255  zt_interp(1,2)=t_interp(1,2,1)+coeff*(t_interp(1,2,2)-t_interp(1,2,1))
1256  zt_interp(2,1)=t_interp(2,1,1)+coeff*(t_interp(2,1,2)-t_interp(2,1,1))
1257  zt_interp(2,2)=t_interp(2,2,1)+coeff*(t_interp(2,2,2)-t_interp(2,2,1))
1258endif
1259
1260!================================================================
1261! 2.3 Latitudinal interpolation
1262!================================================================
1263coeff=(lat-latitude(ilat_inf))/(latitude(ilat_sup)-latitude(ilat_inf))
1264yzt_interp(1)=zt_interp(1,1)+coeff*(zt_interp(1,2)-zt_interp(1,1))
1265yzt_interp(2)=zt_interp(2,1)+coeff*(zt_interp(2,2)-zt_interp(2,1))
1266
1267!================================================================
1268! 2.4 longitudinal interpolation
1269!================================================================
1270coeff=(lon-longitude(ilon_inf))/(longitude(ilon_sup)-longitude(ilon_inf))
1271value=yzt_interp(1)+coeff*(yzt_interp(2)-yzt_interp(1))
1272
1273end subroutine extraction
1274
1275!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1276
1277subroutine inidim(outfid,lonlen,latlen,altlen,timelen,lon,lat,alt,time,units_alt,&
1278                     londimid,latdimid,altdimid,timedimid)
1279!================================================================
1280! Purpose:
1281! Initialize a data file (NetCDF format)
1282!================================================================
1283! Remarks:
1284! The NetCDF file remains open
1285!================================================================
1286use netcdf ! NetCDF definitions
1287implicit none
1288!================================================================
1289! Arguments:
1290!================================================================
1291integer, intent(in):: outfid
1292! outfid: [netcdf] file ID
1293integer, intent(in):: lonlen
1294! lonlen: longitude length (# of longitude values)
1295integer, intent(in):: latlen
1296! latlen: latitude length (# of latitude values)
1297integer, intent(in):: altlen
1298! altlen: altitude length (# of altitude values)
1299integer, intent(in):: timelen
1300! timelen: time length (# of time values)
1301real, intent(in):: lon(lonlen)
1302! lon(): longitude
1303real, intent(in):: lat(latlen)
1304! lat(): latitude
1305real, intent(in):: alt(altlen)
1306! alt(): altitude
1307real, intent(in):: time(timelen)
1308! time(): time (Ls)
1309character(len=1), intent(in) :: units_alt
1310! units_alt: altitude coord. type:'z' (altitude, m) 'p' (pressure, Pa)
1311integer,intent(inout):: londimid
1312! londimid: [netcdf] lon() (i.e.: longitude) ID in MCS and output files (they are the same)
1313integer,intent(inout) :: latdimid
1314! latdimid: [netcdf] lat() ID in MCS and output files (they are the same)
1315integer,intent(inout):: altdimid
1316! altdimid: [netcdf] alt() ID in MCS and output files (they are the same)
1317integer,intent(inout):: timedimid
1318! timedimid: [netcdf] time() ID in MCS and output files (they are the same)
1319
1320!================================================================
1321! Local variables:
1322!================================================================
1323integer :: lonvarid,latvarid,altvarid,timevarid,status
1324! *varid: [netcdf] ID of a variable
1325! status: [netcdf]  return error code (from called subroutines)
1326character(len=256) :: error_text
1327integer :: d ! loop index on dimensions ID
1328
1329!===============================================================
1330! 1. Define/write "dimensions" in the same order than MCS file
1331!    and get their IDs
1332!================================================================
1333do d=1,4
1334  if (altdimid.eq.d) then
1335    status=nf90_def_dim(outfid, "altitude", altlen, altdimid)
1336    error_text="Error: could not define the altitude dimension in the outfile"
1337    call status_check(status,error_text)
1338  else if (timedimid.eq.d) then
1339    status=nf90_def_dim(outfid, "time", timelen, timedimid)
1340    error_text="Error: could not define the time dimension in the outfile"
1341    call status_check(status,error_text)
1342  else if (latdimid.eq.d) then
1343    status=nf90_def_dim(outfid, "latitude", latlen, latdimid)
1344    error_text="Error: could not define the latitude dimension in the outfile"
1345    call status_check(status,error_text)
1346  else if (londimid.eq.d) then
1347    status=nf90_def_dim(outfid, "longitude", lonlen, londimid)
1348    error_text="Error: could not define the longitude dimension in the outfile"
1349    call status_check(status,error_text)
1350  endif
1351enddo
1352
1353!================================================================
1354! 2. Define "variables" and their attributes
1355!================================================================
1356!================================================================
1357! 2.1 Write "longitude" (data and attributes)
1358!================================================================
1359
1360! Insert the definition of the variable
1361status=nf90_def_var(outfid,"longitude",nf90_float,(/londimid/),lonvarid)
1362error_text="Error: could not define the longitude variable in the outfile"
1363call status_check(status,error_text)
1364
1365! Write the attributes
1366status=nf90_put_att(outfid,lonvarid,"long_name","longitude")
1367error_text="Error: could not put the long_name attribute for the longitude variable in the outfile"
1368call status_check(status,error_text)
1369
1370status=nf90_put_att(outfid,lonvarid,"units","degrees_east")
1371error_text="Error: could not put the units attribute for the longitude variable in the outfile"
1372call status_check(status,error_text)
1373
1374!================================================================
1375! 2.2 "latitude"
1376!================================================================
1377
1378! Insert the definition of the variable
1379status=nf90_def_var(outfid,"latitude",nf90_float,(/latdimid/),latvarid)
1380error_text="Error: could not define the latitude variable in the outfile"
1381call status_check(status,error_text)
1382
1383! Write the attributes
1384status=nf90_put_att(outfid,latvarid,"long_name","latitude")
1385error_text="Error: could not put the long_name attribute for the latitude variable in the outfile"
1386call status_check(status,error_text)
1387
1388status=nf90_put_att(outfid,latvarid,"units","degrees_north")
1389error_text="Error: could not put the units attribute for the latitude variable in the outfile"
1390call status_check(status,error_text)
1391
1392!================================================================
1393! 2.3 Write "altitude" (data and attributes)
1394!================================================================
1395
1396! Insert the definition of the variable
1397status=nf90_def_var(outfid,"altitude",nf90_float,(/altdimid/),altvarid)
1398error_text="Error: could not define the altitude variable in the outfile"
1399call status_check(status,error_text)
1400
1401! Write the attributes
1402if (units_alt.eq.'p') then ! pressure coordinate
1403  status=nf90_put_att(outfid,altvarid,"long_name","pressure")
1404  error_text="Error: could not put the long_name attribute for the altitude variable in the outfile"
1405  call status_check(status,error_text)
1406
1407  status=nf90_put_att(outfid,altvarid,'units',"Pa")
1408  error_text="Error: could not put the units attribute for the altitude variable in the outfile"
1409  call status_check(status,error_text)
1410
1411  status=nf90_put_att(outfid,altvarid,'positive',"down")
1412  error_text="Error: could not put the positive attribute for the altitude variable in the outfile"
1413  call status_check(status,error_text)
1414
1415  status=nf90_put_att(outfid,altvarid,'comment',&
1416  "The vertical variable is in fact pressure, not altitude. We just call it altitude for easy reading in Ferret and Grads.")
1417  error_text="Error: could not put the comment attribute for the altitude variable in the outfile"
1418  call status_check(status,error_text)
1419
1420else if (units_alt.eq.'z') then ! altitude coordinate
1421  status=nf90_put_att(outfid,altvarid,"long_name","altitude")
1422  error_text="Error: could not put the long_name attribute for the altitude variable in the outfile"
1423  call status_check(status,error_text)
1424
1425  status=nf90_put_att(outfid,altvarid,'units',"m")
1426  error_text="Error: could not put the units attribute for the altitude variable in the outfile"
1427  call status_check(status,error_text)
1428
1429  status=nf90_put_att(outfid,altvarid,'positive',"up")
1430  error_text="Error: could not put the positive attribute for the altitude variable in the outfile"
1431  call status_check(status,error_text)
1432
1433else
1434  write(*,*)"I do not understand this unit type ",units_alt," for outfile altitude!"
1435  stop
1436end if
1437
1438!================================================================
1439! 2.4 "time"
1440!================================================================
1441
1442! Insert the definition of the variable
1443status=nf90_def_var(outfid,"time",nf90_float,(/timedimid/),timevarid)
1444error_text="Error: could not define the time variable in the outfile"
1445call status_check(status,error_text)
1446
1447! Write the attributes
1448status=nf90_put_att(outfid,timevarid,"long_name","Solar longitude")
1449error_text="Error: could not put the long_name attribute for the time variable in the outfile"
1450call status_check(status,error_text)
1451
1452status=nf90_put_att(outfid,timevarid,"units","days since 0000-01-1 00:00:00")
1453error_text="Error: could not put the units attribute for the time variable in the outfile"
1454call status_check(status,error_text)
1455
1456status=nf90_put_att(outfid,timevarid,"comment",&
1457"Units is in fact degrees, but set to a dummy value of days for compatibility with Ferret and Grads.")
1458error_text="Error: could not put the comment attribute for the time variable in the outfile"
1459call status_check(status,error_text)
1460
1461!================================================================
1462! 2.5 End netcdf define mode
1463!================================================================
1464status=nf90_enddef(outfid)
1465error_text="Error: could not end the define mode of the outfile"
1466call status_check(status,error_text)
1467
1468!================================================================
1469! 3. Write "variables" (data of the dimension variables)
1470!================================================================
1471! "time"
1472status=nf90_put_var(outfid,timevarid,time)
1473error_text="Error: could not write the time variable's data in the outfile"
1474call status_check(status,error_text)
1475
1476! "latitude"
1477status=nf90_put_var(outfid,latvarid,lat)
1478error_text="Error: could not write the latitude variable's data in the outfile"
1479call status_check(status,error_text)
1480
1481! "longitude"
1482status=nf90_put_var(outfid,lonvarid,lon)
1483error_text="Error: could not write the longitude variable's data in the outfile"
1484call status_check(status,error_text)
1485
1486! "altitude"
1487status=nf90_put_var(outfid,altvarid,alt)
1488error_text="Error: could not write the altitude variable's data in the outfile"
1489call status_check(status,error_text)
1490
1491end subroutine inidim
1492
1493!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1494
1495subroutine ls2sol(ls,sol)
1496
1497implicit none
1498!================================================================
1499!  Arguments:
1500!================================================================
1501real,intent(in) :: ls
1502real,intent(out) :: sol
1503
1504!================================================================
1505!  Local:
1506!================================================================
1507double precision xref,zx0,zteta,zz
1508!xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
1509double precision year_day
1510double precision peri_day,timeperi,e_elips
1511double precision pi,degrad
1512parameter (year_day=668.6d0) ! number of sols in a martian year
1513parameter (peri_day=485.35d0) ! date (in sols) of perihelion
1514!timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
1515parameter (timeperi=1.90258341759902d0)
1516parameter (e_elips=0.0934d0)  ! eccentricity of orbit
1517parameter (pi=3.14159265358979d0)
1518parameter (degrad=57.2957795130823d0)
1519
1520      if (abs(ls).lt.1.0e-5) then
1521         if (ls.ge.0.0) then
1522            sol = 0.0
1523         else
1524            sol = real(year_day)
1525         end if
1526         return
1527      end if
1528
1529      zteta = ls/degrad + timeperi
1530      zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips)))
1531      xref = zx0-e_elips*dsin(zx0)
1532      zz = xref/(2.*pi)
1533      sol = real(zz*year_day + peri_day)
1534      if (sol.lt.0.0) sol = sol + real(year_day)
1535      if (sol.ge.year_day) sol = sol - real(year_day)
1536
1537
1538end subroutine ls2sol
1539
1540!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1541
1542subroutine gen_sol_list(solcount,int_sol_list,LTcount,LTave,LTmax,LTmin,lon,f_LT,LTname,&
1543                        sol_list)
1544!================================================================
1545! Purpose:
1546! Generate a list that is the combination of two lists :
1547! - int_sol_list, which is a list of integer values of sol
1548! - LT_list, which contains a number LTcount of LT values in the
1549!   interval [LTmin;LTmax] which are evenly distributed around
1550!   LTave (so that their average is LTave)
1551!================================================================
1552
1553implicit none
1554!================================================================
1555! Arguments:
1556!================================================================
1557integer, intent(in) :: solcount ! nb of integer values of sol
1558real, intent(in) :: int_sol_list(solcount) ! list of the integer values of sol
1559integer, intent(in) :: LTcount ! nb of LT per sol to be interpolated
1560real, intent(in) :: LTave ! average of LT
1561real, intent(in) :: LTmax, LTmin ! bounds of the LT interval for the interpolation
1562real, intent(in) :: lon ! longitude value for the transformation into a sol value at lon=0°
1563external f_LT ! function that changes LT interval for night LT
1564real :: f_LT
1565character (len=64), intent(in) :: LTname ! to know if we have day or night values
1566
1567real, intent(out) :: sol_list(solcount*LTcount) ! all the sol values at lon=0° to interpolate
1568
1569!================================================================
1570! Local variables:
1571!================================================================
1572integer :: N
1573integer :: a,b,c ! loop iteration indexes
1574real :: LT_list(LTcount)
1575
1576N = floor(LTcount/2.)
1577
1578!================================================================
1579! 1. Filling of LT_list
1580!================================================================
1581SELECT CASE (LTname)
1582CASE ("dtemp")
1583  if (abs(LTave-LTmax).le.abs(LTave-LTmin)) then ! LTave is closer to LTmax than to LTmin
1584  do a=1,N
1585    LT_list(a)=LTave+a*abs(LTave-LTmax)/(N+1)
1586    LT_list(N+a)=LTave-a*abs(LTave-LTmax)/(N+1)
1587  enddo
1588else ! LTave is closer to LTmin than to LTmax
1589  do a=1,N
1590    LT_list(a)=LTave+a*abs(LTave-LTmin)/(N+1)
1591    LT_list(N+a)=LTave-a*abs(LTave-LTmin)/(N+1)
1592  enddo
1593endif
1594CASE ("ntemp")
1595  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
1596    do a=1,N
1597      LT_list(a)=LTave+a*abs(f_LT(LTave)-f_LT(LTmax))/(N+1)
1598      LT_list(a)=mod(LT_list(a),24.) ! reput the LT in a [0;24[ interval if needed
1599      LT_list(N+a)=LTave-a*abs(f_LT(LTave)-f_LT(LTmax))/(N+1)
1600      LT_list(N+a)=mod(LT_list(N+a),24.)
1601    enddo
1602  else ! LTave is closer to LTmin than to LTmax
1603    do a=1,N
1604      LT_list(a)=LTave+a*abs(f_LT(LTave)-f_LT(LTmin))/(N+1)
1605      LT_list(a)=mod(LT_list(a),24.)
1606      LT_list(N+a)=LTave-a*abs(f_LT(LTave)-f_LT(LTmin))/(N+1)
1607      LT_list(N+a)=mod(LT_list(N+a),24.)
1608    enddo
1609  endif
1610END SELECT
1611
1612if (mod(LTcount,2).eq.1) then ! if LTcount is an odd number
1613  LT_list(LTcount)=LTave ! add LTave to the list
1614endif
1615
1616!================================================================
1617! 2. Combination of int_sol_list & LT_list into sol_list
1618!================================================================
1619c=1
1620do a=1,solcount
1621  do b=1,LTcount
1622    sol_list(c)=int_sol_list(a)+(LT_list(b)-lon/15.)/24.
1623    c=c+1
1624  enddo
1625enddo
1626
1627end subroutine gen_sol_list
1628
1629!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1630
1631subroutine status_check(status,error_text)
1632
1633use netcdf
1634implicit none
1635!================================================================
1636!  Arguments:
1637!================================================================
1638integer,intent(in) :: status
1639character(len=256),intent(in) :: error_text
1640
1641if (status.ne.nf90_noerr) then
1642  write(*,*)trim(error_text)
1643  write(*,*)trim(nf90_strerror(status))
1644  stop
1645endif
1646
1647end subroutine status_check
1648
1649!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1650
1651real function LTmod(LT)
1652!================================================================
1653! Purpose:
1654! For night local times management, replace hours 
1655! from a [0;24[ interval to a [-12;12[ interval in which
1656! midnight = 0 (in order to ensure continuity when comparing
1657! night local times)
1658!================================================================
1659real :: LT
1660
1661LTmod = 2*mod(LT,12.)-LT
1662return
1663end function LTmod
Note: See TracBrowser for help on using the repository browser.