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

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

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