source: trunk/LMDZ.MARS/util/localtime.F90 @ 2215

Last change on this file since 2215 was 2215, checked in by emillour, 5 years ago

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