source: trunk/LMDZ.MARS/util/terminator.F90 @ 2280

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

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