source: trunk/LMDZ.MARS/util/solzenangle.F90 @ 3026

Last change on this file since 3026 was 2577, checked in by emillour, 3 years ago

Mars GCM utilities:
Fixes in the utilities for the picky gfortran 10+ compiler.
JL+EM
File size: 52.6 KB
Line 
1program solzenangle
2
3! ----------------------------------------------------------------------------
4! Program to redistribute and interpolate the variables at the local time
5! corresponding to one solar zenith angle (morning or evening) everywhere
6! It can especially be used for the terminator (sza=90deg)
7! input : diagfi.nc  / concat.nc / stats.nc kind of files
8! authors: A.Bierjon,F.Forget,2020
9! ----------------------------------------------------------------------------
10
11implicit none
12
13include "netcdf.inc" ! NetCDF definitions
14
15character (len=50)  file
16! file(): input file(s) names(s)
17character (len=30), dimension(16) :: notprocessed
18! notprocessed(): names of the (16) variables that won't be processed
19character (len=100), dimension(:), allocatable :: var
20! var(): name(s) of variable(s) that will be processed
21character (len=100) :: tmpvar,long_name,units
22! tmpvar(): used to temporarily store a variable name
23! long_name(): [netcdf] long_name attribute
24! units(): [netcdf] units attribute
25character (len=100) :: filename,vartmp
26! filename(): output file name
27! vartmp(): temporary variable name (read from netcdf input file)
28character (len=500) :: globatt
29! globatt: [netcdf] global attribute of the output file
30
31character (len=50) :: altlong_name,altunits,altpositive
32! altlong_name(): [netcdf] altitude "long_name" attribute
33! altunits(): [netcdf] altitude "units" attribute
34! altpositive(): [netcdf] altitude "positive" attribute
35
36character (len=7) :: planetside
37! planetside(): planet side of the solar zenith angle ("morning" or "evening")
38
39integer :: nid,ierr,miss,validr
40! nid: [netcdf] file ID #
41! ierr: [netcdf] subroutine returned error code
42! miss: [netcdf] subroutine returned error code
43integer :: i,j,k,inter
44! for various loops
45integer :: varid
46! varid: [netcdf] variable ID #
47real, dimension(:), allocatable:: lat,lon,alt,ctl,time
48! lat(): array, stores latitude coordinates
49! lon(): array, stores longitude coordinates
50! alt(): array, stores altitude coordinates
51! ctl(): array, stores controle array
52! time(): array, stores time coordinates
53integer :: nbvar,nbvarfile,ndim
54! nbvar: # of variables to concatenate
55! nbfile: # number of input file(s)
56! nbvarfile: total # of variables in an input file
57! ndim: [netcdf] # (3 or 4) of dimensions (for variables)
58integer :: latdim,londim,altdim,ctldim,timedim
59! latdim: [netcdf] "latitude" dim ID
60! londim: [netcdf] "longitude" dim ID
61! altdim: [netcdf] "altdim" dim ID
62! ctldim: [netcdf] "controle" dim ID
63! timedim: [netcdf] "timedim" dim ID
64integer :: gcmlayerdim ! NetCDF dimension ID for # of layers in GCM
65integer :: latvar,lonvar,altvar,ctlvar,timevar
66! latvar: [netcdf] ID of "latitude" variable
67! lonvar: [netcdf] ID of "longitude" variable
68! altvar: [netcdf] ID of "altitude" variable
69! ctlvar: [netcdf] ID of "controle" variable
70! timevar: [netcdf] ID of "Time" variable
71integer :: latlen,lonlen,altlen,ctllen,timelen,timelen_tot
72integer :: ilat,ilon,ialt,it
73! latlen: # of elements of lat() array
74! lonlen: # of elements of lon() array
75! altlen: # of elements of alt() array
76! ctllen: # of elements of ctl() array
77! timelen: # of elements of time() array
78! timelen_tot:# =timelen or timelen+1 (if 1 more time to interpolate needed)
79integer :: nsol
80! nsol: # of sols in the input file AND # of elements of time() array in output
81integer :: GCM_layers ! number of GCM atmospheric layers (may not be
82! same as altlen if file is an output of zrecast)
83integer :: nout,latdimout,londimout,altdimout,timedimout,timevarout
84! nout: [netcdf] output file ID
85! latdimout: [netcdf] output latitude (dimension) ID
86! londimout: [netcdf] output longitude (dimension) ID
87! altdimout: [netcdf] output altitude (dimension) ID
88! timedimout: [netcdf] output time (dimension) ID
89! timevarout: [netcdf] ID of output "Time" variable
90integer :: layerdimout ! NetCDF dimension ID for # of layers in GCM
91integer :: varidout
92! varidout: [netcdf] variable ID # (of a variable to write to the output file)
93integer :: Nnotprocessed,var_ok
94! Nnotprocessed: # of (leading)variables that won't be concatenated
95! var_ok: flag (0 or 1)
96integer, dimension(4) :: corner,edges,dim
97! corner: [netcdf]
98! edges: [netcdf]
99! dim: [netcdf]
100real, dimension(:,:,:,:), allocatable :: var3d
101! var3D(,,,): 4D array to store a field
102
103real, dimension(:,:,:,:), allocatable :: var3d_lt
104! var3D_lt(,,,): 4D array to store a field in local time coordinate
105real, dimension(:), allocatable :: lt_gcm
106real, dimension(:), allocatable :: lt_outc ! local lt_out, put in sol decimal value
107real, dimension(:), allocatable :: var_gcm
108real, dimension(:), allocatable :: intsol
109real, dimension(:,:,:), allocatable :: lt_out
110! lt_out: computed local time (hour) corresponding to sza, put in output file
111real :: sza ! solar zenith angle
112character (len=30) :: szastr ! to store the sza value in a string
113real :: startsol ! starting date (in sols) of sol 0 in input file, wrt Ls=0
114real :: tmpsol ! to temporarily store a sol value
115real :: tmpLs ! to temporarily store a Ls value
116real :: tmpLT ! to temporarily store a Local Time value
117
118real :: missing
119!PARAMETER(missing=1E+20)
120! missing: [netcdf] to handle "missing" values when reading/writing files
121real :: miss_lt
122PARAMETER(miss_lt=1E+20)
123! miss_lt: [netcdf] missing value to write for the Local Time in the output file
124real, dimension(2) :: valid_range
125! valid_range(2): [netcdf] interval in which a value is considered valid
126logical :: stats   ! stats=T when reading a "stats" kind of ile
127real :: year_day_gcm = 669. ! number of sols in a year in GCM
128
129
130
131
132!==============================================================================
133! 1.1. Get input file name(s)
134!==============================================================================
135write(*,*) "Welcome in the solar zenith angle interpolator program !"
136write(*,*)
137write(*,*) "which file do you want to use?  (diagfi... stats...  concat...)"
138
139read(*,'(a50)') file
140
141if (len_trim(file).eq.0) then
142   write(*,*) "no file... game over"
143   stop
144endif
145
146!==============================================================================
147! 1.2. Open the input file
148!==============================================================================
149
150ierr = NF_OPEN(file,NF_NOWRITE,nid)
151if (ierr.NE.NF_NOERR) then
152   write(*,*) 'ERROR: Pb opening file '//trim(file)
153   write(*,*) NF_STRERROR(ierr)
154   stop
155endif
156
157ierr=NF_INQ_NVARS(nid,nbvarfile)
158! nbvarfile now set to be the (total) number of variables in file
159if (ierr.NE.NF_NOERR) then
160   write(*,*) 'ERROR: Pb with NF_INQ_NVARS'
161   write(*,*) NF_STRERROR(ierr)
162   stop
163endif
164
165!==============================================================================
166! 1.3. List of variables that should not be processed
167!==============================================================================
168
169notprocessed(1)='Time'
170notprocessed(2)='controle'
171notprocessed(3)='rlonu'
172notprocessed(4)='latitude'
173notprocessed(5)='longitude'
174notprocessed(6)='altitude'
175notprocessed(7)='rlatv'
176notprocessed(8)='aps'
177notprocessed(9)='bps'
178notprocessed(10)='ap'
179notprocessed(11)='bp'
180notprocessed(12)='soildepth'
181notprocessed(13)='cu'
182notprocessed(14)='cv'
183notprocessed(15)='aire'
184notprocessed(16)='phisinit'
185
186!==============================================================================
187! 1.4. Get (and check) list of variables to process
188!==============================================================================
189write(*,*)
190   Nnotprocessed=0
191do i=1,nbvarfile
192   ierr=NF_INQ_VARNAME(nid,i,vartmp)
193   ! vartmp now contains the "name" of variable of ID # i
194   var_ok=0
195   do inter=1,16
196      if (vartmp.eq.notprocessed(inter)) then
197         var_ok=1
198         Nnotprocessed=Nnotprocessed+1
199      endif
200   enddo       
201   if (var_ok.eq.0)  write(*,*) trim(vartmp)
202enddo
203
204! Nnotprocessed: # of variables that won't be processed
205! nbvarfile: total # of variables in file
206allocate(var(nbvarfile-Nnotprocessed),stat=ierr)
207if (ierr.ne.0) then
208  write(*,*) "Error: failed allocation of var(nbvarfile-Nnotprocessed)"
209  write(*,*) "  nbvarfile=",nbvarfile
210  write(*,*) "  Nnotprocessed=",Nnotprocessed
211  stop
212endif
213
214
215write(*,*)
216write(*,*) "which variables do you want to redistribute ?"
217write(*,*) "all / list of <variables> (separated by <Enter>s)"
218write(*,*) "(an empty line , i.e: just <Enter>, implies end of list)"
219nbvar=0
220read(*,'(a100)') tmpvar
221do while ((tmpvar/=' ').AND.(trim(tmpvar)/='all'))
222   nbvar=nbvar+1
223   var(nbvar)=tmpvar
224   read(*,'(a100)') tmpvar
225enddo
226
227if (tmpvar=="all") then
228         nbvar=nbvarfile-Nnotprocessed
229         do j=Nnotprocessed+1,nbvarfile
230            ierr=nf_inq_varname(nid,j,var(j-Nnotprocessed))
231         enddo
232! Variables names from the file are stored in var()
233   nbvar=nbvarfile-Nnotprocessed
234   do i=1,nbvar
235      ierr=nf_inq_varname(nid,i+Nnotprocessed,var(i))
236      write(*,'(a9,1x,i2,1x,a1,1x,a50)') "variable ",i,":",var(i)
237   enddo
238else if(nbvar==0) then
239   write(*,*) "no variable... game over"
240   stop
241endif ! of if (tmpvar=="all")
242
243
244!==============================================================================
245! 2.1. Open input file
246!==============================================================================
247
248write(*,*)
249write(*,*) "opening "//trim(file)//"..."
250ierr = NF_OPEN(file,NF_NOWRITE,nid)
251if (ierr.NE.NF_NOERR) then
252   write(*,*) 'ERROR: Pb opening file '//trim(file)
253   write(*,*) NF_STRERROR(ierr)
254   stop
255endif
256
257!==============================================================================
258! 2.2. Read (and check) dimensions of variables from input file
259!==============================================================================
260
261   ierr=NF_INQ_DIMID(nid,"latitude",latdim)
262   if (ierr.NE.NF_NOERR) then
263      write(*,*) 'ERROR: Dimension <latitude> is missing in file '//trim(file)
264      stop 
265   endif
266   ierr=NF_INQ_VARID(nid,"latitude",latvar)
267   if (ierr.NE.NF_NOERR) then
268      write(*,*) 'ERROR: Field <latitude> is missing in file '//trim(file)
269      stop 
270   endif
271   ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
272!  write(*,*) "latlen: ",latlen
273
274   ierr=NF_INQ_DIMID(nid,"longitude",londim)
275   if (ierr.NE.NF_NOERR) then
276      write(*,*) 'ERROR: Dimension <longitude> is missing in file '//trim(file)
277      stop 
278   endif
279   ierr=NF_INQ_VARID(nid,"longitude",lonvar)
280   if (ierr.NE.NF_NOERR) then
281      write(*,*) 'ERROR: Field <longitude> is missing in file'//trim(file)
282      stop
283   endif
284   ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
285!  write(*,*) "lonlen: ",lonlen
286
287   ierr=NF_INQ_DIMID(nid,"altitude",altdim)
288   if (ierr.NE.NF_NOERR) then
289      write(*,*) 'ERROR: Dimension <altitude> is missing in file '//trim(file)
290      stop 
291   endif
292   ierr=NF_INQ_VARID(nid,"altitude",altvar)
293   if (ierr.NE.NF_NOERR) then
294      write(*,*) 'ERROR: Field <altitude> is missing in file'//trim(file)
295      stop
296   endif
297   ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
298!  write(*,*) "altlen: ",altlen
299
300! load size of aps() or sigma() (in case it is not altlen)
301   ! default is that GCM_layers=altlen
302   ! but for outputs of zrecast, it may be a different value
303   ierr=NF_INQ_DIMID(nid,"GCM_layers",gcmlayerdim)
304   if (ierr.ne.NF_NOERR) then
305     ! didn't find a GCM_layers dimension; therefore we have:
306     GCM_layers=altlen
307   else
308     ! load value of GCM_layers
309     ierr=NF_INQ_DIMLEN(nid,gcmlayerdim,GCM_layers)
310   endif
311!   write(*,*) "GCM_layers=",GCM_layers
312
313   ierr=NF_INQ_DIMID(nid,"index",ctldim)
314   if (ierr.NE.NF_NOERR) then
315      write(*,*) ' Dimension <index> is missing in file '//trim(file)
316      ctllen=0
317      !stop 
318   endif
319   ierr=NF_INQ_VARID(nid,"controle",ctlvar)
320   if (ierr.NE.NF_NOERR) then
321      write(*,*) 'Field <controle> is missing in file '//trim(file)
322      ctllen=0
323      !stop
324   else
325      ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
326   endif
327
328
329!==============================================================================
330! 2.3. Read (and check compatibility of) dimensions of
331!       variables from input file
332!==============================================================================
333
334   allocate(lat(latlen),stat=ierr)
335   if (ierr.ne.0) then
336     write(*,*) "Error: failed to allocate lat(latlen)"
337     stop
338   endif
339   allocate(lon(lonlen),stat=ierr)
340   if (ierr.ne.0) then
341     write(*,*) "Error: failed to allocate lon(lonlen)"
342     stop
343   endif
344   allocate(alt(altlen),stat=ierr)
345   if (ierr.ne.0) then
346     write(*,*) "Error: failed to allocate alt(altlen)"
347     stop
348   endif
349   allocate(ctl(ctllen),stat=ierr)
350   if (ierr.ne.0) then
351     write(*,*) "Error: failed to allocate ctl(ctllen)"
352     stop
353   endif
354
355   ierr = NF_GET_VAR_REAL(nid,latvar,lat)
356   if (ierr.ne.0) then
357     write(*,*) "Error: failed to load latitude"
358     write(*,*) NF_STRERROR(ierr)
359     stop
360   endif
361
362   ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
363   if (ierr.ne.0) then
364     write(*,*) "Error: failed to load longitude"
365     write(*,*) NF_STRERROR(ierr)
366     stop
367   endif
368
369   ierr = NF_GET_VAR_REAL(nid,altvar,alt)
370   if (ierr.ne.0) then
371     write(*,*) "Error: failed to load altitude"
372     write(*,*) NF_STRERROR(ierr)
373     stop
374   endif
375! Get altitude attributes to handle files with any altitude type
376   ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
377   if (ierr.ne.nf_noerr) then
378   ! if no attribute "long_name", try "title"
379     ierr=nf_get_att_text(nid,altvar,"title",altlong_name)
380   endif
381   ierr=nf_get_att_text(nid,altvar,'units',altunits)
382   ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
383
384   if (ctllen .ne. 0) then
385     ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
386     if (ierr.ne.0) then
387       write(*,*) "Error: failed to load controle"
388       write(*,*) NF_STRERROR(ierr)
389       stop
390     endif
391   endif ! of if (ctllen .ne. 0)
392
393!==============================================================================
394! 2.4. Handle "Time" dimension from input file
395!==============================================================================
396
397!==============================================================================
398! 2.4.0 Read "Time" dimension from input file
399!==============================================================================
400   ierr=NF_INQ_DIMID(nid,"Time",timedim)
401   if (ierr.NE.NF_NOERR) then
402      write(*,*) 'ERROR: Dimension <Time> is missing in file'//trim(file)
403      stop
404   endif
405   ierr=NF_INQ_VARID(nid,"Time",timevar)
406   if (ierr.NE.NF_NOERR) then
407      write(*,*) 'ERROR: Field <Time> is missing in file'//trim(file)
408      stop
409   endif
410   ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
411!  write(*,*) "timelen: ",timelen
412
413   ! Sanity Check: Does "Time" has a "long_name" attribute
414   ! and is it "Solar longitude" ?
415   ierr=nf_get_att_text(nid,timevar,"long_name",long_name)
416   if (ierr.ne.nf_noerr) then
417   ! if no attribute "long_name", try "title"
418     ierr=nf_get_att_text(nid,timevar,"title",long_name)
419   endif
420   if ((ierr.EQ.NF_NOERR).and.(long_name.eq."Solar longitude")) then
421     write(*,*) "ERROR: Time axis in input file is in Solar Longitude!"
422     write(*,*) "       localtime requires sols as time axis!"
423     write(*,*) " Might as well stop here."
424     stop
425   endif
426
427   ! allocate time() array and fill it with values from input file
428   allocate(time(timelen),stat=ierr)
429   if (ierr.ne.0) then
430     write(*,*) "Error: failed to allocate time(timelen)"
431     stop
432   endif
433
434   ierr = NF_GET_VAR_REAL(nid,timevar,time)
435   if (ierr.NE.NF_NOERR) then
436     write(*,*) "Error , failed to load Time"
437     write(*,*) NF_STRERROR(ierr)
438     stop
439   endif
440
441!==============================================================================
442! 2.4.1 Select local times to be saved "Time" in output file
443!==============================================================================
444   write(*,*) "Planet side of the Solar Zenith Angle ? (morning or evening)"
445   read(*,*) planetside
446   if ((trim(planetside).ne."morning").and.(trim(planetside).ne."evening")) then
447     write(*,*) "planetside = "//planetside
448     write(*,*) "Error: the planet side must be morning or evening"
449     stop
450   endif
451   write(*,*) "Solar Zenith angle to interpolate ?"
452   write(*,*) "(between 0 and 180 deg, 90deg for the terminator at the surface)"
453   read(*,*) sza
454   if ((sza.lt.0).or.(sza.ge.180)) then
455     write(*,*) "ERROR: the solar zenith angle must lie within [0;180[ deg"
456     stop
457   endif
458
459!  Sol and season: 
460
461   if (ctllen .ne. 0) then
462      startsol = int(int(ctl(4))+int(time(1) + ctl(27)))
463   else
464      startsol = int(time(1))
465   end if
466
467   if ((timelen.eq.12).and.(time(1).eq.2).and.(time(timelen).eq.24))then
468      write(*,*) 'We detect a ""stats"" file'
469      stats=.true.
470!        We rewrite the time from "stats" from 0 to 1 sol...
471      do it=1,timelen  ! loop time input file
472            time(it) = time(it)/24.
473      end do
474      nsol =1
475   else   ! case of a diagfi or concatnc file
476     stats=.false.
477!       Number of sol in input file
478     nsol = int(time(timelen)) - int(time(1))
479   end if
480   allocate(intsol(nsol),stat=ierr)
481   if (ierr.ne.0) then
482     write(*,*) "Error: failed to allocate intsol(nsol)"
483     stop
484   endif
485   allocate(lt_out(lonlen,latlen,nsol),stat=ierr)
486   if (ierr.ne.0) then
487     write(*,*) "Error: failed to allocate lt_out(lonlen,latlen,nsol)"
488     stop
489   endif
490
491   if (stats) then
492         write(*,*) "What is the mean Ls of this stats file (degree) ?"
493         read(*,*,iostat=ierr) tmpLs
494         startsol=0.
495   else
496      write(*,*) "Beginning date (sol since Ls=0) of this file is :",startsol
497      write(*,*) "If you do not want to change it,  just write 'n' (recommended)"
498      write(*,*) "If you really wish to change it, write it now (not recommended!)"
499      read(*,*,iostat=ierr) startsol
500      if((ierr.eq.0).and.(ctllen.ne.0))  ctl(4)=startsol
501      write(*,*) "Beginning date (sol since Ls=0) will be :",startsol
502   end if
503
504   do it=1,nsol
505     intsol(it)= int(time(1)) + it-1.
506     do ilon=1,lonlen
507!         For the computation of Ls, we try to take into account the
508!         rotation time of Mars. It is supposed that the morning/evening
509!         correspond to LT=6h/18h at the given longitude, which is
510!         then reported to a sol value at longitude 0
511       if (trim(planetside).eq."morning") then
512         tmpsol = intsol(it) + (6.-lon(ilon)/15.)/24. + startsol
513       else !if trim(planetside).eq."evening"
514         tmpsol = intsol(it) + (18.-lon(ilon)/15.)/24. + startsol           
515       endif
516       if (.not.stats) call sol2ls(tmpsol,tmpLs)
517       do ilat=1,latlen
518!           Compute the Local Time of the solar zenith angle at a given longitude
519         call LTsolzenangle(lat(ilat),tmpLs,planetside,sza,miss_lt,tmpLT)
520!           Deduce the sol value at this longitude
521         if (tmpLT.eq.miss_lt) then
522           ! if there is no LT corresponding to the sza, put a missing value in lt_out
523           lt_out(ilon,ilat,it)=miss_lt
524         else
525           lt_out(ilon,ilat,it)=tmpLT
526         endif
527       enddo
528     enddo
529   enddo
530
531   if (nsol.gt.1) then
532      timelen_tot=timelen
533   else
534!        if only 1 sol available, we must add 1 timestep for the interpolation
535      timelen_tot=timelen+1
536   endif     
537   allocate(lt_gcm(timelen_tot),stat=ierr)
538   if (ierr.ne.0) then
539     write(*,*) "Error: failed to allocate lt_gcm(timelen_tot)"
540     stop
541   endif
542   allocate(var_gcm(timelen_tot),stat=ierr)
543   if (ierr.ne.0) then
544     write(*,*) "Error: failed to allocate var_gcm(timelen_tot)"
545     stop
546   endif
547   allocate(lt_outc(nsol),stat=ierr)
548   if (ierr.ne.0) then
549     write(*,*) "Error: failed to allocate lt_outc(nsol)"
550     stop
551   endif
552
553!==============================================================================
554! 2.4.2. Get output file name
555!==============================================================================
556   if (trim(planetside).eq."morning") then
557     filename=file(1:len_trim(file)-3)//"_MO.nc"
558   else ! if trim(planetside).eq."evening"
559     filename=file(1:len_trim(file)-3)//"_EV.nc"
560   endif
561
562!==============================================================================
563! 2.4.3. Initiate dimension in output file
564!==============================================================================
565
566
567  ! Initialize output file's lat,lon,alt and time dimensions
568   call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,GCM_layers,&
569         altlong_name,altunits,altpositive,&
570         nout,latdimout,londimout,altdimout,timedimout,&
571         layerdimout,timevarout)
572         
573  ! Initialize output file's aps,bps and phisinit variables
574   call init2(nid,lonlen,latlen,altlen,GCM_layers,&
575              nout,londimout,latdimout,altdimout,&
576              layerdimout)
577
578!==============================================================================
579! 2.4.4. Write/extend "Time" dimension/values in output file
580!==============================================================================
581
582   if (stats) then
583
584     do it=1,nsol
585       ierr= NF_PUT_VARA_REAL(nout,timevarout,(/it/),(/1/),[intsol(it)*24.])
586     enddo
587   else
588     do it=1,nsol
589       ierr= NF_PUT_VARA_REAL(nout,timevarout,(/it/),(/1/),intsol(it))
590       if (ierr.NE.NF_NOERR) then
591         write(*,*) "Error , failed to write Time"
592         stop
593       endif
594     enddo
595   end if
596
597!==============================================================================
598! 2.5 Read/write variables
599!==============================================================================
600
601   do j=1,nbvar ! loop on variables to read/write
602
603!==============================================================================
604! 2.5.1 Check that variable to be read exists in input file
605!==============================================================================
606
607      write(*,*) "variable ",var(j)
608      ierr=nf_inq_varid(nid,var(j),varid)
609      if (ierr.NE.NF_NOERR) then
610         write(*,*) 'ERROR: Field <',var(j),'> not found in file'//file
611         stop
612      endif
613      ierr=nf_inq_varndims(nid,varid,ndim)
614
615      ! Check that it is a 3D or 4D variable
616      if (ndim.lt.3) then
617        write(*,*) "Error:",trim(var(j))," is nor a 3D neither a 4D variable"
618        write(*,*) "We'll skip that variable..."
619        CYCLE ! go directly to the next variable
620      endif
621     
622!==============================================================================
623! 2.5.2 Prepare things in order to read/write the variable
624!==============================================================================
625
626      ! build dim(),corner() and edges() arrays
627      ! and allocate var3d() array
628      if (ndim==3) then
629         allocate(var3d(lonlen,latlen,timelen,1))
630         allocate(var3d_lt(lonlen,latlen,nsol,1))
631         dim(1)=londimout
632         dim(2)=latdimout
633         dim(3)=timedimout
634
635         ! start indexes (where data values will be written)
636         corner(1)=1
637         corner(2)=1
638         corner(3)=1
639         corner(4)=1
640
641         ! length (along dimensions) of block of data to be written
642         edges(1)=lonlen
643         edges(2)=latlen
644         edges(3)=nsol
645         edges(4)=1
646   
647      else if (ndim==4) then
648         allocate(var3d(lonlen,latlen,altlen,timelen))
649         allocate(var3d_lt(lonlen,latlen,altlen,nsol))
650         dim(1)=londimout
651         dim(2)=latdimout
652         dim(3)=altdimout
653         dim(4)=timedimout
654
655         ! start indexes (where data values will be written)
656         corner(1)=1
657         corner(2)=1
658         corner(3)=1
659         corner(4)=1
660
661         ! length (along dimensions) of block of data to be written
662         edges(1)=lonlen
663         edges(2)=latlen
664         edges(3)=altlen
665         edges(4)=nsol
666      endif
667
668      units=" "
669      long_name=" "
670      ierr=nf_get_att_text(nid,varid,"long_name",long_name)
671      if (ierr.ne.nf_noerr) then
672      ! if no attribute "long_name", try "title"
673        ierr=nf_get_att_text(nid,varid,"title",long_name)
674      endif
675      ierr=nf_get_att_text(nid,varid,"units",units)
676      call def_var(nout,var(j),long_name,units,ndim,dim,varidout,ierr)
677
678     
679
680!==============================================================================
681! 2.5.3 Read from input file
682!==============================================================================
683
684      ierr = NF_GET_VAR_REAL(nid,varid,var3d)
685      miss=NF_GET_ATT_REAL(nid,varid,"missing_value",[missing])
686      validr=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
687
688     ! If there is no missing value attribute in the input variable, create one
689      if (miss.ne.NF_NOERR) then
690         missing = miss_lt
691         miss = NF_NOERR
692      endif
693
694!==============================================================================
695! 2.5.4 Interpolation in local time
696!==============================================================================
697
698       do ilon=1,lonlen
699!        Recast GCM local time (sol decimal value) at each longitude :
700         do it=1,timelen  ! loop time input file
701             lt_gcm(it) = time(it) + lon(ilon)/360.
702         enddo
703         if (nsol.eq.1) lt_gcm(timelen+1) = lt_gcm(1)+1
704
705         do ilat=1,latlen
706           do it=1,nsol ! loop time local time
707             if (lt_out(ilon,ilat,it).eq.miss_lt) then
708               lt_outc(it)=miss_lt
709             else
710!              Put computed sza local time (hour) into a sol decimal value :
711               lt_outc(it)=intsol(it) + lt_out(ilon,ilat,it)/24.
712             endif
713             if (nsol.eq.1) then
714!              If 1 sol file (stats), use all data from that sol.
715               if(lt_outc(it).gt.lt_gcm(timelen_tot)) lt_outc(it)=lt_outc(it)-1
716               if(lt_outc(it).lt.lt_gcm(1)) lt_outc(it)=lt_outc(it)+1
717             else
718!              If the computed local time exceeds the GCM input file time bounds,
719!              put a missing value in the local time variable
720               if ((lt_outc(it).gt.lt_gcm(timelen_tot)).OR.(lt_outc(it).lt.lt_gcm(1))) then
721                 lt_out(ilon,ilat,it)=miss_lt
722                 lt_outc(it)=miss_lt
723               endif
724             end if
725           enddo ! local time
726
727           if (ndim==3) then
728             do it=1,timelen  ! loop time input file
729                var_gcm(it) = var3d(ilon,ilat,it,1)
730             enddo
731             if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
732             do it=1,nsol ! loop time local time
733               if (lt_outc(it).eq.miss_lt) then
734                 var3d_lt(ilon,ilat,it,1)=missing
735               else
736                 call interpolf(lt_outc(it),var3d_lt(ilon,ilat,it,1),&
737                               missing,lt_gcm,var_gcm,timelen_tot)
738               endif
739             enddo ! local time
740
741           else if (ndim==4) then
742             do ialt=1,altlen
743               do it=1,timelen ! loop time input file
744                  var_gcm(it) = var3d(ilon,ilat,ialt,it)
745               enddo
746               if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
747               do it=1,nsol ! loop time local time
748                 if (lt_outc(it).eq.miss_lt) then
749                   var3d_lt(ilon,ilat,ialt,it)=missing
750                 else
751                   call interpolf(lt_outc(it),var3d_lt(ilon,ilat,ialt,it),&
752                                 missing,lt_gcm,var_gcm,timelen_tot)
753                 endif
754               enddo ! local time   
755             enddo ! alt
756           endif ! ndim
757
758         enddo ! lat
759       end do ! lon
760
761
762!==============================================================================
763! 2.5.5 write (append) to the output file
764!==============================================================================
765
766       ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3d_lt)
767
768       if (ierr.ne.NF_NOERR) then
769         write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
770         stop
771       endif
772
773      ! Write "missing_value" (and "valid_range" if there is one) attributes in output file
774       call missing_value(nout,varidout,validr,miss,valid_range,missing)
775
776      ! free var3d() array
777       deallocate(var3d)
778       deallocate(var3d_lt)
779
780   enddo ! of do j=1,nbvar
781
782!==============================================================================
783! 2.5.6 Write the solar zenith angle Local Time variable in output file
784!==============================================================================
785   write(*,*) "variable LTsolzenangle"
786   ! build dim(),corner() and edges() arrays
787   dim(1)=londimout
788   dim(2)=latdimout
789   dim(3)=timedimout
790
791   ! start indexes (where data values will be written)
792   corner(1)=1
793   corner(2)=1
794   corner(3)=1
795   corner(4)=1
796
797   ! length (along dimensions) of block of data to be written
798   edges(1)=lonlen
799   edges(2)=latlen
800   edges(3)=nsol
801   edges(4)=1
802
803   tmpvar="LTsolzenangle"
804   units="hours"
805   write(szastr,*) sza
806   if (trim(planetside).eq."morning") then
807     long_name="Morning-side LTST at sza="//trim(adjustl(szastr))//"deg"
808   else !if trim(planetside).eq."evening"
809     long_name="Evening-side LTST at sza="//trim(adjustl(szastr))//"deg"           
810   endif
811   call def_var(nout,tmpvar,long_name,units,3,dim,varidout,ierr)
812
813   ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,lt_out)
814   if (ierr.ne.NF_NOERR) then
815      write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
816      stop
817   endif
818
819   validr=NF_NOERR+1
820   miss=NF_NOERR
821   valid_range(:)=0.
822   call missing_value(nout,varidout,validr,miss,valid_range,miss_lt)
823   write(*,*)"with missing value = ",miss_lt
824!==============================================================================
825! 3. END
826!==============================================================================
827   deallocate(time)
828   deallocate(intsol)
829   deallocate(lt_gcm)
830   deallocate(lt_out)
831   deallocate(lt_outc)
832   deallocate(var_gcm)
833
834   ! Close input file
835   ierr=nf_close(nid)
836
837! Write output file's global attribute
838globatt = "GCM simulation that has been interpolated at the prescribed "//trim(planetside)//&
839          "-side solar zenith angle of "//trim(adjustl(szastr))//"deg all around the globe. "//&
840          "Note that the variable LTsolzenangle stands for the Local (at lon,lat) True Solar Time corresponding "//&
841          "to the SZA. You can derive the corresponding sol at lon=0deg by combining LTsolzenangle, longitude and "//&
842          "Time (integer sol value at lon=0deg) variables. For any question on the program, please contact "//&
843          "antoine.bierjon@lmd.jussieu.fr"
844! Switch to netcdf define mode
845ierr = NF_REDEF (nout)
846ierr = NF_PUT_ATT_TEXT(nout,NF_GLOBAL,'comment',len_trim(globatt),trim(globatt))
847! End netcdf define mode
848ierr = NF_ENDDEF(nout)
849! Close output file
850ierr=NF_CLOSE(nout)
851
852write(*,*) "Program completed !"
853
854end program solzenangle
855
856
857
858!******************************************************************************
859!******************************************************************************
860
861subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,GCM_layers,&
862                     altlong_name,altunits,altpositive,&
863                     nout,latdimout,londimout,altdimout,timedimout,&
864                     layerdimout,timevarout)
865!==============================================================================
866! Purpose:
867! Create and initialize a data file (NetCDF format)
868!==============================================================================
869! Remarks:
870! The NetCDF file (created in this subroutine) remains open
871!==============================================================================
872
873implicit none
874
875include "netcdf.inc" ! NetCDF definitions
876
877!==============================================================================
878! Arguments:
879!==============================================================================
880character (len=*), intent(in):: filename
881! filename(): the file's name
882integer,intent(in) ::latlen,lonlen,altlen,ctllen
883real, intent(in):: lat(latlen)
884! lat(): latitude
885real, intent(in):: lon(lonlen)
886! lon(): longitude
887real, intent(in):: alt(altlen)
888! alt(): altitude
889real, intent(in):: ctl(ctllen)
890! ctl(): controle
891integer,intent(in) :: GCM_layers ! number of GCM layers
892character (len=*), intent(in) :: altlong_name
893! altlong_name(): [netcdf] altitude "long_name" attribute
894character (len=*), intent(in) :: altunits
895! altunits(): [netcdf] altitude "units" attribute
896character (len=*), intent(in) :: altpositive
897! altpositive(): [netcdf] altitude "positive" attribute
898integer, intent(out):: nout
899! nout: [netcdf] file ID
900integer, intent(out):: latdimout
901! latdimout: [netcdf] lat() (i.e.: latitude)  ID
902integer, intent(out):: londimout
903! londimout: [netcdf] lon()  ID
904integer, intent(out):: altdimout
905! altdimout: [netcdf] alt()  ID
906integer, intent(out):: timedimout
907! timedimout: [netcdf] "Time"  ID
908integer,intent(out) :: layerdimout
909! layerdimout: [netcdf] "GCM_layers" ID
910integer, intent(out):: timevarout
911! timevarout: [netcdf] "Time" (considered as a variable) ID
912
913!==============================================================================
914! Local variables:
915!==============================================================================
916!integer :: latdim,londim,altdim,timedim
917integer :: nvarid,ierr
918integer :: ctldimout
919! nvarid: [netcdf] ID of a variable
920! ierr: [netcdf]  return error code (from called subroutines)
921
922!==============================================================================
923! 1. Create (and open) output file
924!==============================================================================
925write(*,*) "creating "//trim(adjustl(filename))//'...'
926ierr = NF_CREATE(filename,IOR(NF_CLOBBER,NF_64BIT_OFFSET),nout)
927! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
928if (ierr.NE.NF_NOERR) then
929   WRITE(*,*)'ERROR: Impossible to create the file ',trim(filename)
930   write(*,*) NF_STRERROR(ierr)
931   stop
932endif
933
934!==============================================================================
935! 2. Define/write "dimensions" and get their IDs
936!==============================================================================
937
938ierr = NF_DEF_DIM(nout, "latitude", latlen, latdimout)
939if (ierr.NE.NF_NOERR) then
940   WRITE(*,*)'initiate: error failed to define dimension <latitude>'
941   write(*,*) NF_STRERROR(ierr)
942   stop
943endif
944ierr = NF_DEF_DIM(nout, "longitude", lonlen, londimout)
945if (ierr.NE.NF_NOERR) then
946   WRITE(*,*)'initiate: error failed to define dimension <longitude>'
947   write(*,*) NF_STRERROR(ierr)
948   stop
949endif
950ierr = NF_DEF_DIM(nout, "altitude", altlen, altdimout)
951if (ierr.NE.NF_NOERR) then
952   WRITE(*,*)'initiate: error failed to define dimension <altitude>'
953   write(*,*) NF_STRERROR(ierr)
954   stop
955endif
956if (size(ctl).ne.0) then
957  ierr = NF_DEF_DIM(nout, "index", ctllen, ctldimout)
958  if (ierr.NE.NF_NOERR) then
959    WRITE(*,*)'initiate: error failed to define dimension <index>'
960    write(*,*) NF_STRERROR(ierr)
961    stop
962  endif
963endif
964ierr = NF_DEF_DIM(nout, "Time", NF_UNLIMITED, timedimout)
965if (ierr.NE.NF_NOERR) then
966   WRITE(*,*)'initiate: error failed to define dimension <Time>'
967   write(*,*) NF_STRERROR(ierr)
968   stop
969endif
970
971ierr = NF_DEF_DIM(nout, "GCM_layers", GCM_layers, layerdimout)
972if (ierr.NE.NF_NOERR) then
973   WRITE(*,*)'initiate: error failed to define dimension <GCM_layers>'
974   write(*,*) NF_STRERROR(ierr)
975   stop
976endif
977
978! End netcdf define mode
979ierr = NF_ENDDEF(nout)
980if (ierr.NE.NF_NOERR) then
981   WRITE(*,*)'initiate: error failed to switch out of define mode'
982   write(*,*) NF_STRERROR(ierr)
983   stop
984endif
985!==============================================================================
986! 3. Write "Time" and its attributes
987!==============================================================================
988
989call def_var(nout,"Time","Time","years since 0000-00-0 00:00:00",1,&
990             (/timedimout/),timevarout,ierr)
991! Switch to netcdf define mode
992ierr = NF_REDEF (nout)
993ierr = NF_PUT_ATT_TEXT(nout,timevarout,'comment',135,&
994"The true unit of the variable Time is the integer value "//&
995"of sol at lon=0deg. A false unit is put to enable reading from Grads or Ferret")
996! End netcdf define mode
997ierr = NF_ENDDEF(nout)
998
999!==============================================================================
1000! 4. Write "latitude" (data and attributes)
1001!==============================================================================
1002
1003call def_var(nout,"latitude","latitude","degrees_north",1,&
1004             (/latdimout/),nvarid,ierr)
1005
1006ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
1007if (ierr.NE.NF_NOERR) then
1008   WRITE(*,*)'initiate: error failed writing variable <latitude>'
1009   write(*,*) NF_STRERROR(ierr)
1010   stop
1011endif
1012
1013!==============================================================================
1014! 4. Write "longitude" (data and attributes)
1015!==============================================================================
1016
1017call def_var(nout,"longitude","East longitude","degrees_east",1,&
1018             (/londimout/),nvarid,ierr)
1019
1020ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
1021if (ierr.NE.NF_NOERR) then
1022   WRITE(*,*)'initiate: error failed writing variable <longitude>'
1023   write(*,*) NF_STRERROR(ierr)
1024   stop
1025endif
1026
1027!==============================================================================
1028! 5. Write "altitude" (data and attributes)
1029!==============================================================================
1030
1031! Switch to netcdf define mode
1032ierr = NF_REDEF (nout)
1033
1034ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,[altdimout],nvarid)
1035
1036ierr = NF_PUT_ATT_TEXT (nout,nvarid,'long_name',len_trim(adjustl(altlong_name)),adjustl(altlong_name))
1037ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',len_trim(adjustl(altunits)),adjustl(altunits))
1038ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',len_trim(adjustl(altpositive)),adjustl(altpositive))
1039
1040! End netcdf define mode
1041ierr = NF_ENDDEF(nout)
1042
1043ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
1044if (ierr.NE.NF_NOERR) then
1045   WRITE(*,*)'initiate: error failed writing variable <altitude>'
1046   write(*,*) NF_STRERROR(ierr)
1047   stop
1048endif
1049
1050!==============================================================================
1051! 6. Write "controle" (data and attributes)
1052!==============================================================================
1053
1054if (size(ctl).ne.0) then
1055   ! Switch to netcdf define mode
1056   ierr = NF_REDEF (nout)
1057
1058   ierr = NF_DEF_VAR (nout,"controle",NF_FLOAT,1,[ctldimout],nvarid)
1059
1060   ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",18,"Control parameters")
1061
1062   ! End netcdf define mode
1063   ierr = NF_ENDDEF(nout)
1064
1065   ierr = NF_PUT_VAR_REAL (nout,nvarid,ctl)
1066   if (ierr.NE.NF_NOERR) then
1067      WRITE(*,*)'initiate: error failed writing variable <controle>'
1068      write(*,*) NF_STRERROR(ierr)
1069      stop
1070   endif
1071endif
1072
1073end Subroutine initiate
1074
1075!******************************************************************************
1076subroutine init2(infid,lonlen,latlen,altlen,GCM_layers, &
1077                 outfid,londimout,latdimout,altdimout, &
1078                 layerdimout)
1079!==============================================================================
1080! Purpose:
1081! Copy aps(), bps() and phisinit() from input file to output file
1082!==============================================================================
1083! Remarks:
1084! The NetCDF files must be open
1085!==============================================================================
1086
1087implicit none
1088
1089include "netcdf.inc" ! NetCDF definitions
1090
1091!==============================================================================
1092! Arguments:
1093!==============================================================================
1094integer, intent(in) :: infid  ! NetCDF output file ID
1095integer, intent(in) :: lonlen ! # of grid points along longitude
1096integer, intent(in) :: latlen ! # of grid points along latitude
1097integer, intent(in) :: altlen ! # of grid points along altitude
1098integer, intent(in) :: GCM_layers ! # of GCM atmospheric layers
1099integer, intent(in) :: outfid ! NetCDF output file ID
1100integer, intent(in) :: londimout ! longitude dimension ID
1101integer, intent(in) :: latdimout ! latitude dimension ID
1102integer, intent(in) :: altdimout ! altitude dimension ID
1103integer, intent(in) :: layerdimout ! GCM_layers dimension ID
1104!==============================================================================
1105! Local variables:
1106!==============================================================================
1107real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
1108real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
1109integer :: apsid,bpsid,phisinitid
1110integer :: ierr
1111integer :: tmpdimid ! temporary dimension ID
1112integer :: tmpvarid ! temporary variable ID
1113integer :: tmplen ! temporary variable length
1114logical :: phis,hybrid ! are "phisinit" and "aps"/"bps" available ?
1115
1116
1117!==============================================================================
1118! 1. Read data from input file
1119!==============================================================================
1120
1121! hybrid coordinate aps
1122!write(*,*) "aps: altlen=",altlen," GCM_layers=",GCM_layers
1123allocate(aps(GCM_layers),stat=ierr)
1124if (ierr.ne.0) then
1125  write(*,*) "init2: failed to allocate aps!"
1126  stop
1127endif
1128ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
1129if (ierr.ne.NF_NOERR) then
1130  write(*,*) "Ooops. Failed to get aps ID. OK."
1131  hybrid=.false.
1132else
1133  ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
1134  hybrid=.true.
1135  if (ierr.ne.NF_NOERR) then
1136    write(*,*) "init2 Error: Failed reading aps"
1137    stop
1138  endif
1139
1140  ! hybrid coordinate bps
1141!  write(*,*) "bps: altlen=",altlen," GCM_layers=",GCM_layers
1142  allocate(bps(GCM_layers),stat=ierr)
1143  if (ierr.ne.0) then
1144    write(*,*) "init2: failed to allocate bps!"
1145    stop
1146  endif
1147  ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
1148  if (ierr.ne.NF_NOERR) then
1149    write(*,*) "Ooops. Failed to get bps ID. OK."
1150    hybrid=.false.
1151  else
1152    ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
1153    if (ierr.ne.NF_NOERR) then
1154      write(*,*) "init2 Error: Failed reading bps"
1155      stop
1156    endif
1157  endif
1158endif
1159
1160! ground geopotential phisinit
1161allocate(phisinit(lonlen,latlen),stat=ierr)
1162if (ierr.ne.0) then
1163  write(*,*) "init2: failed to allocate phisinit(lonlen,latlen)"
1164  stop
1165endif
1166ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
1167if (ierr.ne.NF_NOERR) then
1168  write(*,*) "Failed to get phisinit ID. OK"
1169  phisinit = 0.
1170  phis = .false.
1171else
1172  ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
1173  if (ierr.ne.NF_NOERR) then
1174    write(*,*) "init2 Error: Failed reading phisinit"
1175    stop
1176  endif
1177  phis = .true.
1178endif
1179
1180
1181!==============================================================================
1182! 2. Write
1183!==============================================================================
1184
1185!==============================================================================
1186! 2.2. Hybrid coordinates aps() and bps()
1187!==============================================================================
1188
1189IF (hybrid) then
1190  ! define aps
1191  call def_var(outfid,"aps","hybrid pressure at midlayers"," ",1,&
1192               (/layerdimout/),apsid,ierr)
1193  if (ierr.ne.NF_NOERR) then
1194    write(*,*) "Error: Failed to def_var aps"
1195    stop
1196  endif
1197
1198  ! write aps
1199  ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
1200  if (ierr.ne.NF_NOERR) then
1201    write(*,*) "Error: Failed to write aps"
1202    stop
1203  endif
1204 
1205  ! define bps
1206  call def_var(outfid,"bps","hybrid sigma at midlayers"," ",1,&
1207               (/layerdimout/),bpsid,ierr)
1208  if (ierr.ne.NF_NOERR) then
1209    write(*,*) "Error: Failed to def_var bps"
1210    stop
1211  endif
1212
1213  ! write bps
1214  ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
1215  if (ierr.ne.NF_NOERR) then
1216    write(*,*) "Error: Failed to write bps"
1217    stop
1218  endif
1219 
1220  deallocate(aps)
1221  deallocate(bps)
1222ENDIF ! of IF (hybrid)
1223
1224!==============================================================================
1225! 2.2. phisinit()
1226!==============================================================================
1227
1228IF (phis) THEN
1229
1230  !define phisinit
1231  call def_var(outfid,"phisinit","Ground level geopotential"," ",2,&
1232              (/londimout,latdimout/),phisinitid,ierr)
1233  if (ierr.ne.NF_NOERR) then
1234     write(*,*) "Error: Failed to def_var phisinit"
1235     stop
1236  endif
1237
1238  ! write phisinit
1239  ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
1240  if (ierr.ne.NF_NOERR) then
1241    write(*,*) "Error: Failed to write phisinit"
1242    stop
1243  endif
1244
1245deallocate(phisinit)
1246ENDIF ! of IF (phis)
1247
1248
1249end subroutine init2
1250
1251!******************************************************************************
1252subroutine def_var(nid,name,long_name,units,nbdim,dim,nvarid,ierr)
1253!==============================================================================
1254! Purpose: Write a variable (i.e: add a variable to a dataset)
1255! called "name"; along with its attributes "long_name", "units"...
1256! to a file (following the NetCDF format)
1257!==============================================================================
1258! Remarks:
1259! The NetCDF file must be open
1260!==============================================================================
1261
1262implicit none
1263
1264include "netcdf.inc" ! NetCDF definitions
1265
1266!==============================================================================
1267! Arguments:
1268!==============================================================================
1269integer, intent(in) :: nid
1270! nid: [netcdf] file ID #
1271character (len=*), intent(in) :: name
1272! name(): [netcdf] variable's name
1273character (len=*), intent(in) :: long_name
1274! long_name(): [netcdf] variable's "long_name" attribute
1275character (len=*), intent(in) :: units
1276! unit(): [netcdf] variable's "units" attribute
1277integer, intent(in) :: nbdim
1278! nbdim: number of dimensions of the variable
1279integer, dimension(nbdim), intent(in) :: dim
1280! dim(nbdim): [netcdf] dimension(s) ID(s)
1281integer, intent(out) :: nvarid
1282! nvarid: [netcdf] ID # of the variable
1283integer, intent(out) :: ierr
1284! ierr: [netcdf] subroutines returned error code
1285
1286! Switch to netcdf define mode
1287ierr=NF_REDEF(nid)
1288
1289! Insert the definition of the variable
1290ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
1291
1292! Write the attributes
1293ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name",len_trim(adjustl(long_name)),adjustl(long_name))
1294ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
1295
1296! End netcdf define mode
1297ierr=NF_ENDDEF(nid)
1298
1299end subroutine def_var
1300
1301!******************************************************************************
1302subroutine  missing_value(nout,nvarid,validr,miss,valid_range,missing)
1303!==============================================================================
1304! Purpose:
1305! Write "valid_range" and "missing_value" attributes (of a given
1306! variable) to a netcdf file
1307!==============================================================================
1308! Remarks:
1309! NetCDF file must be open
1310! Variable (nvarid) ID must be set
1311!==============================================================================
1312
1313implicit none
1314
1315include "netcdf.inc"  ! NetCDF definitions
1316
1317!==============================================================================
1318! Arguments:
1319!==============================================================================
1320integer, intent(in) :: nout
1321! nout: [netcdf] file ID #
1322integer, intent(in) :: nvarid
1323! varid: [netcdf] variable ID #
1324integer, intent(in) :: validr
1325! validr : [netcdf] routines return code for "valid_range" attribute
1326integer, intent(in) :: miss
1327! miss : [netcdf] routines return code for "missing_value" attribute
1328real, dimension(2), intent(in) :: valid_range
1329! valid_range(2): [netcdf] "valid_range" attribute (min and max)
1330real, intent(in) :: missing
1331! missing: [netcdf] "missing_value" attribute
1332
1333!==============================================================================
1334! Local variables:
1335!==============================================================================
1336integer :: ierr
1337! ierr: [netcdf] subroutine returned error code
1338!      INTEGER nout,nvarid,ierr
1339!      REAL missing
1340!      REAL valid_range(2)
1341
1342! Switch to netcdf dataset define mode
1343ierr = NF_REDEF (nout)
1344if (ierr.ne.NF_NOERR) then
1345   print*,'missing_value: '
1346   print*, NF_STRERROR(ierr)
1347endif
1348
1349! Write valid_range() attribute
1350if (validr.eq.NF_NOERR) then
1351  ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
1352
1353  if (ierr.ne.NF_NOERR) then
1354     print*,'missing_value: valid_range attribution failed'
1355     print*, NF_STRERROR(ierr)
1356     !write(*,*) 'NF_NOERR', NF_NOERR
1357     !CALL abort
1358     stop
1359  endif
1360endif
1361
1362! Write "missing_value" attribute
1363if (miss.eq.NF_NOERR) then
1364  ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,[missing])
1365
1366  if (ierr.NE.NF_NOERR) then
1367     print*, 'missing_value: missing value attribution failed'
1368     print*, NF_STRERROR(ierr)
1369  !    WRITE(*,*) 'NF_NOERR', NF_NOERR
1370  !    CALL abort
1371     stop
1372  endif
1373endif
1374
1375! End netcdf dataset define mode
1376ierr = NF_ENDDEF(nout)
1377
1378end subroutine  missing_value
1379
1380!*****************************************************************************
1381subroutine interpolf(x,y,missing,xd,yd,nd)
1382!==============================================================================
1383! Purpose:
1384! Yield y=f(x), where xd() end yd() are arrays of known values,
1385! using linear interpolation
1386! If x is not included in the interval spaned by xd(), then y is set
1387! to a default value 'missing'
1388! Note:
1389! Array xd() should contain ordered (either increasing or decreasing) abscissas
1390!==============================================================================
1391implicit none
1392!==============================================================================
1393! Arguments:
1394!==============================================================================
1395real,intent(in) :: x ! abscissa at which interpolation is to be done
1396real,intent(in) :: missing ! missing value (if no interpolation is performed)
1397integer,intent(in) :: nd ! size of arrays
1398real,dimension(nd),intent(in) :: xd ! array of known absissas
1399real,dimension(nd),intent(in) :: yd ! array of correponding values
1400
1401real,intent(out) :: y ! interpolated value
1402!==============================================================================
1403! Local variables:
1404!==============================================================================
1405integer :: i
1406
1407! default: set y to 'missing'
1408y=missing
1409
1410   do i=1,nd-1
1411     if (((x.ge.xd(i)).and.(x.le.xd(i+1))).or.&
1412          ((x.le.xd(i)).and.(x.ge.xd(i+1)))) then
1413        if ((yd(i).eq.missing).or.(yd(i+1).eq.missing)) then
1414          ! cannot perform the interpolation if an encompassing value
1415          ! is already set to 'missing'
1416        else
1417          !linear interpolation based on encompassing values
1418          y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i))
1419        endif
1420        exit
1421     endif
1422   enddo
1423
1424
1425end subroutine interpolf
1426
1427!******************************************************************************
1428subroutine LTsolzenangle(lat,ls,planetside,sza,missval,lt_sza)
1429! compute the localtime by inversing
1430! the solar zenith angle equation :
1431! cos(sza)=sin(lat)*sin(declin)
1432!         +cos(lat)*cos(declin)*cos(2.*pi*(localtime/24.-.5))
1433
1434
1435implicit none
1436
1437real,intent(in) :: lat ! latitude (deg)
1438real,intent(in) :: ls ! solar longitude (deg)
1439character (len=7), intent(in) :: planetside ! planet side : "morning" or "evening"
1440real,intent(in) :: sza ! solar zenith angle (deg) (within [0;180[)
1441real,intent(in) :: missval ! missing value for when there is no associated LT (polar day/night)
1442real,intent(out) :: lt_sza ! local true solar time (hours) corresponding to the sza
1443
1444! local variables:
1445double precision :: declin ! declination of the Sun (rad)
1446double precision :: tmpcos ! to temporarily store a cosine value
1447double precision,parameter :: obliquity=25.1919d0
1448double precision,parameter :: pi=3.14159265358979d0
1449double precision,parameter :: degtorad=pi/180.0d0
1450double precision,parameter :: radtodeg=180.0d0/pi
1451
1452! Compute Sun's declination
1453declin=ASIN(sin(ls*degtorad)*sin(obliquity*degtorad))
1454
1455! Small checks
1456if ((cos(lat*degtorad).eq.0).or.(cos(declin).eq.0)) then
1457  write(*,*) "ERROR in LTsolzenangle : lat and declin can't have their cosine equal to zero"
1458  stop
1459endif
1460
1461! Compute Local Time
1462tmpcos = cos(sza*degtorad)/(cos(lat*degtorad)*cos(declin)) - tan(lat*degtorad)*tan(declin)
1463
1464if ((tmpcos.gt.1.).or.(tmpcos.lt.(-1.))) then
1465  ! Detect Polar Day and Polar Night (values out of the arccos definition domain)
1466  lt_sza = missval
1467 
1468else
1469
1470  if (trim(planetside).eq."morning") then
1471  ! Local Time of the Morning-side Solar Zenith Angle
1472    lt_sza = 24*(0.5 - ACOS(tmpcos)/(2.*pi))
1473  else ! if trim(planetside).eq."evening"
1474  ! Local Time of the Evening-side Solar Zenith Angle
1475    lt_sza = 24*(0.5 + ACOS(tmpcos)/(2.*pi))
1476  endif
1477 
1478endif
1479
1480end subroutine LTsolzenangle
1481
1482!******************************************************************************
1483subroutine sol2ls(sol,ls)
1484!  convert a given martian day number (sol)
1485!  into corresponding solar longitude, Ls (in degr.),
1486!  where sol=0=Ls=0 is the
1487!  northern hemisphere spring equinox.
1488
1489implicit none
1490
1491!input
1492real,intent(in) :: sol
1493!output
1494real,intent(out) :: ls
1495
1496!local variables
1497double precision :: year_day,peri_day,timeperi,e_elips
1498parameter (year_day=668.6d0) ! number of martian days (sols) in a martian year
1499parameter (peri_day=485.35d0) ! perihelion date (in sols)
1500parameter (e_elips=0.09340d0) ! orbital eccentricity
1501parameter (timeperi=1.90258341759902d0) ! timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
1502
1503double precision :: pi,radtodeg
1504parameter (pi=3.14159265358979d0)
1505parameter (radtodeg=57.2957795130823d0) !  radtodeg: 180/pi
1506
1507double precision :: zanom,xref,zx0,zdx,zteta,zz
1508!  xref: mean anomaly, zx0: eccentric anomaly, zteta: true anomaly       
1509integer :: iter
1510
1511zz=(sol-peri_day)/year_day
1512zanom=2.*pi*(zz-nint(zz))
1513xref=dabs(zanom)
1514
1515!  The equation zx0 - e * sin (zx0) = xref, solved by Newton
1516zx0=xref+e_elips*dsin(xref)
1517iter=1
1518iteration:do while (iter.le.10)
1519   zdx=-(zx0-e_elips*dsin(zx0)-xref)/(1.-e_elips*dcos(zx0))
1520   if(dabs(zdx).le.(1.d-7)) then ! typically, 2 or 3 iterations are enough
1521     EXIT iteration
1522   endif
1523   zx0=zx0+zdx
1524   iter=iter+1
1525enddo iteration
1526zx0=zx0+zdx
1527if(zanom.lt.0.) zx0=-zx0
1528
1529! compute true anomaly zteta, now that eccentric anomaly zx0 is known
1530zteta=2.*datan(dsqrt((1.+e_elips)/(1.-e_elips))*dtan(zx0/2.))
1531
1532! compute Ls
1533ls=real(zteta-timeperi)
1534if(ls.lt.0.) ls=ls+2.*real(pi)
1535if(ls.gt.2.*pi) ls=ls-2.*real(pi)
1536! convert Ls in deg.
1537ls=real(radtodeg)*ls
1538
1539end subroutine sol2ls
1540
1541!******************************************************************************
Note: See TracBrowser for help on using the repository browser.