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

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

Mars GCM:
Update utilities :
concat.F90

  • rewrite and simplify the handling of time and offset so that any file can be concatenated, including files from different years or stats file
  • use modulo to shift starting sols in concat.nc to values between 0 and 669
  • add a "adls" option in addition to "sol" or "ls" to add a Ls 1D variable

while keeping "sol" as the Time axis

  • add conservation of altitude attributes (long_name,units,positive)
  • enable absence of both hybrid (aps,bps) and sigma levels

solzenangle.F90

  • improve calculation for 1 sol file (stats) to use all local time data
  • read the starting sol in the file instead of asking it to the user (except for stats file) ; keep the possibility for user to change it
  • ask mean Ls value for stats file instead of sol number
  • fix crash when using "all variable" option (ticket #66)
  • fix bug on aps,bps by using GCM_layers dimension instead of altitude

localtime.F90

  • fix crash when using "all variable" option (ticket #66)
  • fix bug on aps,bps by using GCM_layers dimension instead of altitude

For all 3 utilities :

  • handle all ndim cases for the variables (ticket #52)
  • change "title" attribute into "long_name" by default (ticket #48)
  • extend size of long_name character string (ticket #57)
  • remove the use of #ifdef NC_DOUBLE (ticket #67)

AB
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=50), 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(*,'(a50)') tmpvar
221do while ((tmpvar/=' ').AND.(trim(tmpvar)/='all'))
222   nbvar=nbvar+1
223   var(nbvar)=tmpvar
224   read(*,'(a50)') 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    stop "init2 Error: Failed reading aps"
1137  endif
1138
1139  ! hybrid coordinate bps
1140!  write(*,*) "bps: altlen=",altlen," GCM_layers=",GCM_layers
1141  allocate(bps(GCM_layers),stat=ierr)
1142  if (ierr.ne.0) then
1143    write(*,*) "init2: failed to allocate bps!"
1144    stop
1145  endif
1146  ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
1147  if (ierr.ne.NF_NOERR) then
1148    write(*,*) "Ooops. Failed to get bps ID. OK."
1149    hybrid=.false.
1150  else
1151    ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
1152    if (ierr.ne.NF_NOERR) then
1153      stop "init2 Error: Failed reading bps"
1154    endif
1155  endif
1156endif
1157
1158! ground geopotential phisinit
1159allocate(phisinit(lonlen,latlen),stat=ierr)
1160if (ierr.ne.0) then
1161  write(*,*) "init2: failed to allocate phisinit(lonlen,latlen)"
1162  stop
1163endif
1164ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
1165if (ierr.ne.NF_NOERR) then
1166  write(*,*) "Failed to get phisinit ID. OK"
1167  phisinit = 0.
1168  phis = .false.
1169else
1170  ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
1171  if (ierr.ne.NF_NOERR) then
1172    stop "init2 Error: Failed reading phisinit"
1173  endif
1174  phis = .true.
1175endif
1176
1177
1178!==============================================================================
1179! 2. Write
1180!==============================================================================
1181
1182!==============================================================================
1183! 2.2. Hybrid coordinates aps() and bps()
1184!==============================================================================
1185
1186IF (hybrid) then
1187  ! define aps
1188  call def_var(outfid,"aps","hybrid pressure at midlayers"," ",1,&
1189               (/layerdimout/),apsid,ierr)
1190  if (ierr.ne.NF_NOERR) then
1191    stop "Error: Failed to def_var aps"
1192  endif
1193
1194  ! write aps
1195  ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
1196  if (ierr.ne.NF_NOERR) then
1197    stop "Error: Failed to write aps"
1198  endif
1199 
1200  ! define bps
1201  call def_var(outfid,"bps","hybrid sigma at midlayers"," ",1,&
1202               (/layerdimout/),bpsid,ierr)
1203  if (ierr.ne.NF_NOERR) then
1204    stop "Error: Failed to def_var bps"
1205  endif
1206
1207  ! write bps
1208  ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
1209  if (ierr.ne.NF_NOERR) then
1210    stop "Error: Failed to write bps"
1211  endif
1212 
1213  deallocate(aps)
1214  deallocate(bps)
1215ENDIF ! of IF (hybrid)
1216
1217!==============================================================================
1218! 2.2. phisinit()
1219!==============================================================================
1220
1221IF (phis) THEN
1222
1223  !define phisinit
1224  call def_var(outfid,"phisinit","Ground level geopotential"," ",2,&
1225              (/londimout,latdimout/),phisinitid,ierr)
1226  if (ierr.ne.NF_NOERR) then
1227     stop "Error: Failed to def_var phisinit"
1228  endif
1229
1230  ! write phisinit
1231  ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
1232  if (ierr.ne.NF_NOERR) then
1233    stop "Error: Failed to write phisinit"
1234  endif
1235
1236deallocate(phisinit)
1237ENDIF ! of IF (phis)
1238
1239
1240end subroutine init2
1241
1242!******************************************************************************
1243subroutine def_var(nid,name,long_name,units,nbdim,dim,nvarid,ierr)
1244!==============================================================================
1245! Purpose: Write a variable (i.e: add a variable to a dataset)
1246! called "name"; along with its attributes "long_name", "units"...
1247! to a file (following the NetCDF format)
1248!==============================================================================
1249! Remarks:
1250! The NetCDF file must be open
1251!==============================================================================
1252
1253implicit none
1254
1255include "netcdf.inc" ! NetCDF definitions
1256
1257!==============================================================================
1258! Arguments:
1259!==============================================================================
1260integer, intent(in) :: nid
1261! nid: [netcdf] file ID #
1262character (len=*), intent(in) :: name
1263! name(): [netcdf] variable's name
1264character (len=*), intent(in) :: long_name
1265! long_name(): [netcdf] variable's "long_name" attribute
1266character (len=*), intent(in) :: units
1267! unit(): [netcdf] variable's "units" attribute
1268integer, intent(in) :: nbdim
1269! nbdim: number of dimensions of the variable
1270integer, dimension(nbdim), intent(in) :: dim
1271! dim(nbdim): [netcdf] dimension(s) ID(s)
1272integer, intent(out) :: nvarid
1273! nvarid: [netcdf] ID # of the variable
1274integer, intent(out) :: ierr
1275! ierr: [netcdf] subroutines returned error code
1276
1277! Switch to netcdf define mode
1278ierr=NF_REDEF(nid)
1279
1280! Insert the definition of the variable
1281ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
1282
1283! Write the attributes
1284ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name",len_trim(adjustl(long_name)),adjustl(long_name))
1285ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
1286
1287! End netcdf define mode
1288ierr=NF_ENDDEF(nid)
1289
1290end subroutine def_var
1291
1292!******************************************************************************
1293subroutine  missing_value(nout,nvarid,validr,miss,valid_range,missing)
1294!==============================================================================
1295! Purpose:
1296! Write "valid_range" and "missing_value" attributes (of a given
1297! variable) to a netcdf file
1298!==============================================================================
1299! Remarks:
1300! NetCDF file must be open
1301! Variable (nvarid) ID must be set
1302!==============================================================================
1303
1304implicit none
1305
1306include "netcdf.inc"  ! NetCDF definitions
1307
1308!==============================================================================
1309! Arguments:
1310!==============================================================================
1311integer, intent(in) :: nout
1312! nout: [netcdf] file ID #
1313integer, intent(in) :: nvarid
1314! varid: [netcdf] variable ID #
1315integer, intent(in) :: validr
1316! validr : [netcdf] routines return code for "valid_range" attribute
1317integer, intent(in) :: miss
1318! miss : [netcdf] routines return code for "missing_value" attribute
1319real, dimension(2), intent(in) :: valid_range
1320! valid_range(2): [netcdf] "valid_range" attribute (min and max)
1321real, intent(in) :: missing
1322! missing: [netcdf] "missing_value" attribute
1323
1324!==============================================================================
1325! Local variables:
1326!==============================================================================
1327integer :: ierr
1328! ierr: [netcdf] subroutine returned error code
1329!      INTEGER nout,nvarid,ierr
1330!      REAL missing
1331!      REAL valid_range(2)
1332
1333! Switch to netcdf dataset define mode
1334ierr = NF_REDEF (nout)
1335if (ierr.ne.NF_NOERR) then
1336   print*,'missing_value: '
1337   print*, NF_STRERROR(ierr)
1338endif
1339
1340! Write valid_range() attribute
1341if (validr.eq.NF_NOERR) then
1342  ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
1343
1344  if (ierr.ne.NF_NOERR) then
1345     print*,'missing_value: valid_range attribution failed'
1346     print*, NF_STRERROR(ierr)
1347     !write(*,*) 'NF_NOERR', NF_NOERR
1348     !CALL abort
1349     stop ""
1350  endif
1351endif
1352
1353! Write "missing_value" attribute
1354if (miss.eq.NF_NOERR) then
1355  ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
1356
1357  if (ierr.NE.NF_NOERR) then
1358     print*, 'missing_value: missing value attribution failed'
1359     print*, NF_STRERROR(ierr)
1360  !    WRITE(*,*) 'NF_NOERR', NF_NOERR
1361  !    CALL abort
1362     stop ""
1363  endif
1364endif
1365
1366! End netcdf dataset define mode
1367ierr = NF_ENDDEF(nout)
1368
1369end subroutine  missing_value
1370
1371!*****************************************************************************
1372subroutine interpolf(x,y,missing,xd,yd,nd)
1373!==============================================================================
1374! Purpose:
1375! Yield y=f(x), where xd() end yd() are arrays of known values,
1376! using linear interpolation
1377! If x is not included in the interval spaned by xd(), then y is set
1378! to a default value 'missing'
1379! Note:
1380! Array xd() should contain ordered (either increasing or decreasing) abscissas
1381!==============================================================================
1382implicit none
1383!==============================================================================
1384! Arguments:
1385!==============================================================================
1386real,intent(in) :: x ! abscissa at which interpolation is to be done
1387real,intent(in) :: missing ! missing value (if no interpolation is performed)
1388integer,intent(in) :: nd ! size of arrays
1389real,dimension(nd),intent(in) :: xd ! array of known absissas
1390real,dimension(nd),intent(in) :: yd ! array of correponding values
1391
1392real,intent(out) :: y ! interpolated value
1393!==============================================================================
1394! Local variables:
1395!==============================================================================
1396integer :: i
1397
1398! default: set y to 'missing'
1399y=missing
1400
1401   do i=1,nd-1
1402     if (((x.ge.xd(i)).and.(x.le.xd(i+1))).or.&
1403          ((x.le.xd(i)).and.(x.ge.xd(i+1)))) then
1404        if ((yd(i).eq.missing).or.(yd(i+1).eq.missing)) then
1405          ! cannot perform the interpolation if an encompassing value
1406          ! is already set to 'missing'
1407        else
1408          !linear interpolation based on encompassing values
1409          y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i))
1410        endif
1411        exit
1412     endif
1413   enddo
1414
1415
1416end subroutine interpolf
1417
1418!******************************************************************************
1419subroutine LTsolzenangle(lat,ls,planetside,sza,missval,lt_sza)
1420! compute the localtime by inversing
1421! the solar zenith angle equation :
1422! cos(sza)=sin(lat)*sin(declin)
1423!         +cos(lat)*cos(declin)*cos(2.*pi*(localtime/24.-.5))
1424
1425
1426implicit none
1427
1428real,intent(in) :: lat ! latitude (deg)
1429real,intent(in) :: ls ! solar longitude (deg)
1430character (len=7), intent(in) :: planetside ! planet side : "morning" or "evening"
1431real,intent(in) :: sza ! solar zenith angle (deg) (within [0;180[)
1432real,intent(in) :: missval ! missing value for when there is no associated LT (polar day/night)
1433real,intent(out) :: lt_sza ! local true solar time (hours) corresponding to the sza
1434
1435! local variables:
1436double precision :: declin ! declination of the Sun (rad)
1437double precision :: tmpcos ! to temporarily store a cosine value
1438double precision,parameter :: obliquity=25.1919d0
1439double precision,parameter :: pi=3.14159265358979d0
1440double precision,parameter :: degtorad=pi/180.0d0
1441double precision,parameter :: radtodeg=180.0d0/pi
1442
1443! Compute Sun's declination
1444declin=ASIN(sin(ls*degtorad)*sin(obliquity*degtorad))
1445
1446! Small checks
1447if ((cos(lat*degtorad).eq.0).or.(cos(declin).eq.0)) then
1448  write(*,*) "ERROR in LTsolzenangle : lat and declin can't have their cosine equal to zero"
1449  stop
1450endif
1451
1452! Compute Local Time
1453tmpcos = cos(sza*degtorad)/(cos(lat*degtorad)*cos(declin)) - tan(lat*degtorad)*tan(declin)
1454
1455if ((tmpcos.gt.1.).or.(tmpcos.lt.(-1.))) then
1456  ! Detect Polar Day and Polar Night (values out of the arccos definition domain)
1457  lt_sza = missval
1458 
1459else
1460
1461  if (trim(planetside).eq."morning") then
1462  ! Local Time of the Morning-side Solar Zenith Angle
1463    lt_sza = 24*(0.5 - ACOS(tmpcos)/(2.*pi))
1464  else ! if trim(planetside).eq."evening"
1465  ! Local Time of the Evening-side Solar Zenith Angle
1466    lt_sza = 24*(0.5 + ACOS(tmpcos)/(2.*pi))
1467  endif
1468 
1469endif
1470
1471end subroutine LTsolzenangle
1472
1473!******************************************************************************
1474subroutine sol2ls(sol,ls)
1475!  convert a given martian day number (sol)
1476!  into corresponding solar longitude, Ls (in degr.),
1477!  where sol=0=Ls=0 is the
1478!  northern hemisphere spring equinox.
1479
1480implicit none
1481
1482!input
1483real,intent(in) :: sol
1484!output
1485real,intent(out) :: ls
1486
1487!local variables
1488double precision :: year_day,peri_day,timeperi,e_elips
1489parameter (year_day=668.6d0) ! number of martian days (sols) in a martian year
1490parameter (peri_day=485.35d0) ! perihelion date (in sols)
1491parameter (e_elips=0.09340d0) ! orbital eccentricity
1492parameter (timeperi=1.90258341759902d0) ! timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
1493
1494double precision :: pi,radtodeg
1495parameter (pi=3.14159265358979d0)
1496parameter (radtodeg=57.2957795130823d0) !  radtodeg: 180/pi
1497
1498double precision :: zanom,xref,zx0,zdx,zteta,zz
1499!  xref: mean anomaly, zx0: eccentric anomaly, zteta: true anomaly       
1500integer :: iter
1501
1502zz=(sol-peri_day)/year_day
1503zanom=2.*pi*(zz-nint(zz))
1504xref=dabs(zanom)
1505
1506!  The equation zx0 - e * sin (zx0) = xref, solved by Newton
1507zx0=xref+e_elips*dsin(xref)
1508iter=1
1509iteration:do while (iter.le.10)
1510   zdx=-(zx0-e_elips*dsin(zx0)-xref)/(1.-e_elips*dcos(zx0))
1511   if(dabs(zdx).le.(1.d-7)) then ! typically, 2 or 3 iterations are enough
1512     EXIT iteration
1513   endif
1514   zx0=zx0+zdx
1515   iter=iter+1
1516enddo iteration
1517zx0=zx0+zdx
1518if(zanom.lt.0.) zx0=-zx0
1519
1520! compute true anomaly zteta, now that eccentric anomaly zx0 is known
1521zteta=2.*datan(dsqrt((1.+e_elips)/(1.-e_elips))*dtan(zx0/2.))
1522
1523! compute Ls
1524ls=real(zteta-timeperi)
1525if(ls.lt.0.) ls=ls+2.*real(pi)
1526if(ls.gt.2.*pi) ls=ls-2.*real(pi)
1527! convert Ls in deg.
1528ls=real(radtodeg)*ls
1529
1530end subroutine sol2ls
1531
1532!******************************************************************************
Note: See TracBrowser for help on using the repository browser.