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

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

Mars GCM:
Correction in solzenangle.F90 : the missing value was not initialized when the attribute was not
present in the input file.
+ Resolved Ticket #61 : add comments in the code and in the output file to make the program more
transparent for the users

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