source: trunk/LMDZ.VENUS/Tools/localtime.F90 @ 1317

Last change on this file since 1317 was 1305, checked in by slebonnois, 11 years ago

SL: VENUS PHOTOCHEMISTRY. Needs Lapack (see arch files...)
File size: 33.4 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
13include "planet.h"
14
15character (len=50)  file
16! file(): input file(s) names(s)
17character (len=30), dimension(6) :: notconcat
18! notconcat(): names of the (15) variables that won't be concatenated
19character (len=50), dimension(:), allocatable :: var
20! var(): name(s) of variable(s) that will be concatenated
21character (len=50) :: tmpvar,title,units
22! tmpvar(): used to temporarily store a variable name
23! title(): [netcdf] title attribute
24! units(): [netcdf] units attribute
25character (len=100) :: filename,vartmp
26! filename(): output file name
27! vartmp(): temporary variable name (read from netcdf input file)
28!character (len=1) :: ccopy
29! ccpy: 'y' or 'n' answer
30integer :: nid,ierr,miss
31! nid: [netcdf] file ID #
32! ierr: [netcdf] subroutine returned error code
33! miss: [netcdf] subroutine returned error code
34integer :: i,j,k,inter
35! for various loops
36integer :: varid
37! varid: [netcdf] variable ID #
38real, dimension(:), allocatable:: lat,lon,alt,ctl,time
39! lat(): array, stores latitude coordinates
40! lon(): array, stores longitude coordinates
41! alt(): array, stores altitude coordinates
42! ctl(): array, stores controle array
43! time(): array, stores time coordinates
44integer :: nbvar,nbvarfile,ndim
45! nbvar: # of variables to concatenate
46! nbfile: # number of input file(s)
47! nbvarfile: total # of variables in an input file
48! ndim: [netcdf] # (3 or 4) of dimensions (for variables)
49integer :: latdim,londim,altdim,ctldim,timedim
50! latdim: [netcdf] "latitude" dim ID
51! londim: [netcdf] "longitude" dim ID
52! altdim: [netcdf] "altdim" dim ID
53! ctldim: [netcdf] "controle" dim ID
54! timedim: [netcdf] "timedim" dim ID
55integer :: latvar,lonvar,altvar,ctlvar,timevar
56! latvar: [netcdf] ID of "latitude" variable
57! lonvar: [netcdf] ID of "longitude" variable
58! altvar: [netcdf] ID of "altitude" variable
59! ctlvar: [netcdf] ID of "controle" variable
60! timevar: [netcdf] ID of "Time" variable
61integer :: latlen,lonlen,altlen,ctllen,timelen,timelen_lt,timelen_tot
62integer :: ilat,ilon,ialt,it
63! latlen: # of elements of lat() array
64! lonlen: # of elements of lon() array
65! altvar: # of elements of alt() array
66! ctlvar: # of elements of ctl() array
67! timelen: # of elemnets of time() array
68! timelen_tot:# =timelen or timelen+1 (if 1 more time to interpolate needed)
69! timelen_lt: # of elemnets of time() array in output
70integer :: nout,latdimout,londimout,altdimout,timedimout,timevarout
71integer :: nhour,nsol
72! nout: [netcdf] output file ID
73! latdimout: [netcdf] output latitude (dimension) ID
74! londimout: [netcdf] output longitude (dimension) ID
75! altdimout: [netcdf] output altitude (dimension) ID
76! timedimout: [netcdf] output time (dimension) ID
77! timevarout: [netcdf] ID of output "Time" variable
78integer :: varidout
79! varidout: [netcdf] variable ID # (of a variable to write to the output file)
80integer :: Nnotconcat,var_ok
81! Nnotconcat: # of (leading)variables that won't be concatenated
82! var_ok: flag (0 or 1)
83integer, dimension(4) :: corner,edges,dim
84! corner: [netcdf]
85! edges: [netcdf]
86! dim: [netcdf]
87real, dimension(:,:,:,:), allocatable :: var3d
88! var3D(,,,): 4D array to store a field
89
90real, dimension(:,:,:,:), allocatable :: var3d_lt
91! var3D_lt(,,,): 4D array to store a field in local time coordinate
92real, dimension(:), allocatable :: lt_gcm,lt_hour
93real, dimension(:), allocatable :: lt_out,lt_outc
94real, dimension(:), allocatable :: var_gcm
95
96real :: missing
97!PARAMETER(missing=1E+20)
98! missing: [netcdf] to handle "missing" values when reading/writing files
99real, dimension(2) :: valid_range
100! valid_range(2): [netcdf] interval in which a value is considered valid
101logical :: stats   ! stats=T when reading a "stats" kind of ile
102
103
104!==============================================================================
105! 1.1. Get input file name(s)
106!==============================================================================
107write(*,*)
108write(*,*) "which file do you want to use?  (type histmth.nc)"
109
110read(*,'(a50)') file
111
112if (len_trim(file).eq.0) then
113   write(*,*) "no file... game over"
114   stop ""
115endif
116
117!==============================================================================
118! 1.3. Open the first input file
119!==============================================================================
120
121ierr = NF_OPEN(file,NF_NOWRITE,nid)
122if (ierr.NE.NF_NOERR) then
123   write(*,*) 'ERROR: Pb opening file '//file
124   stop ""
125endif
126
127ierr=NF_INQ_NVARS(nid,nbvarfile)
128! nbvarfile now set to be the (total) number of variables in file
129
130!==============================================================================
131! 1.4. Ask for (output) "Time" axis type
132!==============================================================================
133
134notconcat(1)='time_counter'
135notconcat(2)='lat'
136notconcat(3)='lon'
137notconcat(4)='presnivs'
138notconcat(5)='aire'
139notconcat(6)='phis'
140
141!==============================================================================
142! 1.5. Get (and check) list of variables to concatenate
143!==============================================================================
144write(*,*)
145   Nnotconcat=0
146do i=1,nbvarfile
147   ierr=NF_INQ_VARNAME(nid,i,vartmp)
148   ! vartmp now contains the "name" of variable of ID # i
149   var_ok=0
150   do inter=1,6
151      if (vartmp.eq.notconcat(inter)) then
152         var_ok=1
153         Nnotconcat=Nnotconcat+1
154      endif
155   enddo       
156   if (var_ok.eq.0)  write(*,*) trim(vartmp)
157enddo
158
159! Nnotconcat: # of variables that won't be concatenated
160! nbvarfile: total # of variables in file
161allocate(var(nbvarfile-Nnotconcat))
162
163
164write(*,*)
165write(*,*) "which variables do you want to redistribute ?"
166write(*,*) "all / list of <variables> (separated by <Enter>s)"
167write(*,*) "(an empty line , i.e: just <Enter>, implies end of list)"
168nbvar=0
169read(*,'(a50)') tmpvar
170do while ((tmpvar/=' ').AND.(trim(tmpvar)/='all'))
171   nbvar=nbvar+1
172   var(nbvar)=tmpvar
173   read(*,'(a50)') tmpvar
174enddo
175
176if (tmpvar=="all") then
177         nbvar=nbvarfile-Nnotconcat
178         do j=Nnotconcat+1,nbvarfile
179            ierr=nf_inq_varname(nid,j,var(j-Nnotconcat))
180         enddo
181! Variables names from the file are catched
182   nbvar=nbvarfile-Nnotconcat
183   do i=1,nbvar
184      ierr=nf_inq_varname(nid,i+Nnotconcat,var(i))
185      write(*,'(a9,1x,i2,1x,a1,1x,a50)') "variable ",i,":",var(i)
186   enddo
187else if(nbvar==0) then
188   write(*,*) "no variable... game over"
189   stop ""
190endif ! of if (tmpvar=="all")
191
192!==============================================================================
193! 1.6. Get output file name
194!==============================================================================
195 filename=file(1:len_trim(file)-3)//"_LT.nc"
196 write(*,*) filename
197
198!==============================================================================
199! 2.1. Open input file
200!==============================================================================
201
202   if (i/=1) then
203      write(*,*)
204      write(*,*) "opening "//trim(file)//"..."
205      ierr = NF_OPEN(file,NF_NOWRITE,nid)
206      if (ierr.NE.NF_NOERR) then
207         write(*,*) 'ERROR: Pb opening file '//file
208         write(*,*) NF_STRERROR(ierr)
209         stop ""
210      endif
211   endif
212
213!==============================================================================
214! 2.2. Read (and check) dimensions of variables from input file
215!==============================================================================
216
217   ierr=NF_INQ_DIMID(nid,"lat",latdim)
218   ierr=NF_INQ_VARID(nid,"lat",latvar)
219   if (ierr.NE.NF_NOERR) then
220      write(*,*) 'ERROR: Field <lat> is missing in file'//file
221      stop "" 
222   endif
223   ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
224!  write(*,*) "latlen: ",latlen
225
226   ierr=NF_INQ_DIMID(nid,"lon",londim)
227   ierr=NF_INQ_VARID(nid,"lon",lonvar)
228   if (ierr.NE.NF_NOERR) then
229      write(*,*) 'ERROR: Field <lon> is missing in file'//file
230      stop ""
231   endif
232   ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
233!  write(*,*) "lonlen: ",lonlen
234
235   ierr=NF_INQ_DIMID(nid,"presnivs",altdim)
236   ierr=NF_INQ_VARID(nid,"presnivs",altvar)
237   if (ierr.NE.NF_NOERR) then
238      write(*,*) 'ERROR: Field <presnivs> is missing in file'//file
239      stop ""
240   endif
241   ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
242!  write(*,*) "altlen: ",altlen
243
244   ierr=NF_INQ_DIMID(nid,"index",ctldim)
245   ierr=NF_INQ_VARID(nid,"controle",ctlvar)
246   if (ierr.NE.NF_NOERR) then
247      write(*,*) 'Field <controle> is missing in file'//file
248      ctllen=0
249      !stop ""
250   else
251      ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
252   endif
253!  write(*,*) "controle: ",controle
254
255!==============================================================================
256! 2.3. Read (and check compatibility of) dimensions of
257!       variables from input file
258!==============================================================================
259
260! First call; initialize/allocate
261      allocate(lat(latlen))
262      allocate(lon(lonlen))
263      allocate(alt(altlen))
264      allocate(ctl(ctllen))
265#ifdef NC_DOUBLE
266      ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat)
267      ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon)
268      ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt)
269      if (ctllen .ne. -1) ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
270#else
271      ierr = NF_GET_VAR_REAL(nid,latvar,lat)
272      ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
273      ierr = NF_GET_VAR_REAL(nid,altvar,alt)
274      if (ctllen .ne. -1) ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
275#endif
276!==============================================================================
277! 2.4. Handle "Time" dimension from input file
278!==============================================================================
279
280!==============================================================================
281! 2.4.0 Read "Time" dimension from input file
282!==============================================================================
283   ierr=NF_INQ_DIMID(nid,"time_counter",timedim)
284   if (ierr.NE.NF_NOERR) then
285      write(*,*) 'ERROR: Dimension <time_counter> is missing in file'//file
286      stop ""
287   endif
288   ierr=NF_INQ_VARID(nid,"time_counter",timevar)
289   if (ierr.NE.NF_NOERR) then
290      write(*,*) 'ERROR: Field <time_counter> is missing in file'//file
291      stop ""
292   endif
293   ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
294!  write(*,*) "timelen: ",timelen
295
296   ! allocate time() array and fill it with values from input file
297   allocate(time(timelen))
298
299#ifdef NC_DOUBLE
300   ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
301#else
302   ierr = NF_GET_VAR_REAL(nid,timevar,time)
303#endif
304
305
306!==============================================================================
307! 2.4.1 Select local times to be saved "Time" in output file
308!==============================================================================
309      write(*,*) 'number of local time to be used ?'
310      read(*,*) nhour
311      allocate(lt_hour(nhour))
312      write(*,*) 'list of selected local time (0<t<24)'
313      do it=1,nhour
314         read(*,*) lt_hour(it)
315      end do
316
317      if ((timelen.eq.12).and.(time(1).eq.2).and.(time(timelen).eq.24))then
318         write(*,*) 'We detect a ""stats"" file'
319         stop
320      else   ! case Venus histmth.nc file
321        stats=.false.
322!       Number of sol in input file
323        nsol = int((time(timelen)-time(1))/localday)+1
324        write(*,*) time
325        write(*,*) "Number of days in input file:",nsol
326        timelen_lt=nhour*nsol
327        allocate(lt_out(timelen_lt))
328        i=0
329        do k=1,nsol
330          do it=1,nhour
331            i=i+1
332            lt_out(i)=int(time(1)/localday) + k-1  + lt_hour(it)/24.
333          end do
334        end do
335        end if
336
337      if (nsol.gt.1) then
338         timelen_tot=timelen
339      else
340!        if only 1 sol available, we must add 1 timestep for the interpolation
341         timelen_tot=timelen+1
342      endif     
343      allocate(lt_gcm(timelen_tot))
344      allocate(var_gcm(timelen_tot))
345
346      allocate(lt_outc(timelen_lt))
347         
348
349!==============================================================================
350! 2.4.1.5 Initiate dimension in output file
351!==============================================================================
352
353
354   ! Initialize output file's lat,lon,alt and time dimensions
355      call initiate (filename,lat,lon,alt,ctl,nout,&
356           latdimout,londimout,altdimout,timedimout,timevarout)
357   ! Initialize output file's phis variables
358     call init2(nid,lonlen,latlen,altlen,&
359                nout,londimout,latdimout,altdimout)
360
361
362!==============================================================================
363! 2.4.2 Write/extend "Time" dimension/values in output file
364!==============================================================================
365
366   if (stats) then
367
368     do it=1,timelen_lt
369#ifdef NC_DOUBLE
370        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,lt_hour(it))
371#else
372        ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,lt_hour(it))
373#endif
374     enddo
375   else
376     do it=1,timelen_lt
377#ifdef NC_DOUBLE
378        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,lt_out(it))
379#else
380        ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,lt_out(it))
381#endif
382     enddo
383  end if
384
385!==============================================================================
386! 2.5 Read/write variables
387!==============================================================================
388
389   do j=1,nbvar ! loop on variables to read/write
390
391!==============================================================================
392! 2.5.1 Check that variable to be read existe in input file
393!==============================================================================
394
395      write(*,*) "variable ",var(j)
396      ierr=nf_inq_varid(nid,var(j),varid)
397      if (ierr.NE.NF_NOERR) then
398         write(*,*) 'ERROR: Field <',var(j),'> not found in file'//file
399         stop ""
400      endif
401      ierr=nf_inq_varndims(nid,varid,ndim)
402
403!==============================================================================
404! 2.5.2 Prepare things in order to read/write the variable
405!==============================================================================
406
407      ! build dim(),corner() and edges() arrays
408      ! and allocate var3d() array
409      if (ndim==3) then
410         allocate(var3d(lonlen,latlen,timelen,1))
411         allocate(var3d_lt(lonlen,latlen,timelen_lt,1))
412         dim(1)=londimout
413         dim(2)=latdimout
414         dim(3)=timedimout
415
416         ! start indexes (where data values will be written)
417         corner(1)=1
418         corner(2)=1
419         corner(3)=1
420         corner(4)=1
421
422         ! length (along dimensions) of block of data to be written
423         edges(1)=lonlen
424         edges(2)=latlen
425         edges(3)=timelen_lt
426         edges(4)=1
427   
428      else if (ndim==4) then
429         allocate(var3d(lonlen,latlen,altlen,timelen))
430         allocate(var3d_lt(lonlen,latlen,altlen,timelen_lt))
431         dim(1)=londimout
432         dim(2)=latdimout
433         dim(3)=altdimout
434         dim(4)=timedimout
435
436         ! start indexes (where data values will be written)
437         corner(1)=1
438         corner(2)=1
439         corner(3)=1
440         corner(4)=1
441
442         ! length (along dimensions) of block of data to be written
443         edges(1)=lonlen
444         edges(2)=latlen
445         edges(3)=altlen
446         edges(4)=timelen_lt
447      endif
448
449         units="                                                    "
450         title="                                                    "
451         ierr=nf_get_att_text(nid,varid,"title",title)
452         ierr=nf_get_att_text(nid,varid,"units",units)
453         call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr)
454
455     
456
457!==============================================================================
458! 2.5.3 Read from input file
459!==============================================================================
460
461#ifdef NC_DOUBLE
462      ierr = NF_GET_VAR_DOUBLE(nid,varid,var3d)
463      miss=NF_GET_ATT_DOUBLE(nid,varid,"missing_value",missing)
464      miss=NF_GET_ATT_DOUBLE(nid,varid,"valid_range",valid_range)
465#else
466      ierr = NF_GET_VAR_REAL(nid,varid,var3d)
467      miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
468      miss=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
469#endif
470
471!==============================================================================
472! 2.5.4 Read from input file
473!==============================================================================
474! Interpolation in local time :
475
476        do ilon=1,lonlen
477!          write(*,*) 'processing longitude =', lon(ilon)
478!          Local time at each longitude :
479           do it=1,timelen  ! loop temps in
480! suppose that:
481!  - we're on Venus
482!  - the run is initiated when LT(lon=0)=0
483!  - the time points in output are located at mid-intervals
484               lt_gcm(it) = int(time(1)/localday) + &
485                            (it-0.5)*nsol/timelen - 0.5*lon(ilon)/180.
486           end do
487           if (nsol.eq.1) lt_gcm(timelen+1) = lt_gcm(1)+1
488
489           do it=1,timelen_lt ! loop time local time
490              lt_outc(it)=lt_out(it)
491!             Correction to use same local time previous or following
492!             day where data are missing :
493              if(lt_outc(it).gt.lt_gcm(timelen_tot))lt_outc(it)=lt_out(it)-1
494              if(lt_outc(it).lt.lt_gcm(1))lt_outc(it)=lt_out(it)+1
495           end do
496
497           if (ndim==3) then
498             do ilat=1,latlen
499               do it=1,timelen  ! loop temps in
500                  var_gcm(it) = var3d(ilon,ilat,it,1)
501               end do
502               if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
503               do it=1,timelen_lt ! loop time local time
504                  call interpolf(lt_outc(it),var3d_lt(ilon,ilat,it,1),&
505                                 missing,lt_gcm,var_gcm,timelen_tot)
506               end do   
507             enddo
508
509           else if  (ndim==4) then
510             do ilat=1,latlen
511               do ialt=1,altlen
512                 do it=1,timelen ! loop temps in
513                    var_gcm(it) = var3d(ilon,ilat,ialt,it)
514                 end do
515                 if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
516                 do it=1,timelen_lt ! loop time local time
517                    call interpolf(lt_outc(it),var3d_lt(ilon,ilat,ialt,it),&
518                                   missing,lt_gcm,var_gcm,timelen_tot)
519                 end do   
520               enddo
521             enddo
522           end if
523
524        end do
525
526
527!==============================================================================
528! 2.5.5 write (append) to the output file
529!==============================================================================
530
531#ifdef NC_DOUBLE
532      ierr= NF_PUT_VARA_DOUBLE(nout,varidout,corner,edges,var3d_lt)
533#else
534      ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3d_lt)
535#endif
536
537      if (ierr.ne.NF_NOERR) then
538         write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
539         stop ""
540      endif
541
542! In case there is a "valid_range" attribute
543      ! Write "valid_range" and "missing_value" attributes in output file
544      if (miss.eq.NF_NOERR) then
545         call missing_value(nout,varidout,valid_range,missing)
546      endif
547
548      ! free var3d() array
549      deallocate(var3d)
550      deallocate(var3d_lt)
551
552   enddo ! of do j=1,nbvar
553
554   deallocate(time)
555   deallocate(lt_gcm)
556   deallocate(lt_out)
557   deallocate(lt_outc)
558   deallocate(var_gcm)
559
560   ! Close input file
561   ierr=nf_close(nid)
562
563! Close output file
564ierr=NF_CLOSE(nout)
565
566contains
567
568!******************************************************************************
569Subroutine initiate (filename,lat,lon,alt,ctl,&
570                     nout,latdimout,londimout,altdimout,timedimout,timevarout)
571!==============================================================================
572! Purpose:
573! Create and initialize a data file (NetCDF format)
574!==============================================================================
575! Remarks:
576! The NetCDF file (created in this subroutine) remains open
577!==============================================================================
578
579implicit none
580
581include "netcdf.inc" ! NetCDF definitions
582
583!==============================================================================
584! Arguments:
585!==============================================================================
586character (len=*), intent(in):: filename
587! filename(): the file's name
588real, dimension(:), intent(in):: lat
589! lat(): latitude
590real, dimension(:), intent(in):: lon
591! lon(): longitude
592real, dimension(:), intent(in):: alt
593! alt(): altitude
594real, dimension(:), intent(in):: ctl
595! ctl(): controle
596integer, intent(out):: nout
597! nout: [netcdf] file ID
598integer, intent(out):: latdimout
599! latdimout: [netcdf] lat() (i.e.: latitude)  ID
600integer, intent(out):: londimout
601! londimout: [netcdf] lon()  ID
602integer, intent(out):: altdimout
603! altdimout: [netcdf] alt()  ID
604integer, intent(out):: timedimout
605! timedimout: [netcdf] "Time"  ID
606integer, intent(out):: timevarout
607! timevarout: [netcdf] "Time" (considered as a variable) ID
608
609!==============================================================================
610! Local variables:
611!==============================================================================
612!integer :: latdim,londim,altdim,timedim
613integer :: nvarid,ierr
614integer :: ctldimout
615! nvarid: [netcdf] ID of a variable
616! ierr: [netcdf]  return error code (from called subroutines)
617
618!==============================================================================
619! 1. Create (and open) output file
620!==============================================================================
621write(*,*) "creating "//trim(adjustl(filename))//'...'
622ierr = NF_CREATE(filename,IOR(NF_CLOBBER,NF_64BIT_OFFSET),nout)
623! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
624if (ierr.NE.NF_NOERR) then
625   WRITE(*,*)'ERROR: Impossible to create the file.'
626   stop ""
627endif
628
629!==============================================================================
630! 2. Define/write "dimensions" and get their IDs
631!==============================================================================
632
633ierr = NF_DEF_DIM(nout, "latitude", size(lat), latdimout)
634ierr = NF_DEF_DIM(nout, "longitude", size(lon), londimout)
635ierr = NF_DEF_DIM(nout, "presnivs", size(alt), altdimout)
636if (size(ctl).ne.0) ierr = NF_DEF_DIM(nout, "index", size(ctl), ctldimout)
637ierr = NF_DEF_DIM(nout, "Time", NF_UNLIMITED, timedimout)
638
639! End netcdf define mode
640ierr = NF_ENDDEF(nout)
641
642!==============================================================================
643! 3. Write "Time" and its attributes
644!==============================================================================
645
646call def_var(nout,"Time","Time","years since 0000-00-0 00:00:00",1,&
647             (/timedimout/),timevarout,ierr)
648
649!==============================================================================
650! 4. Write "latitude" (data and attributes)
651!==============================================================================
652
653call def_var(nout,"latitude","latitude","degrees_north",1,&
654             (/latdimout/),nvarid,ierr)
655
656#ifdef NC_DOUBLE
657ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lat)
658#else
659ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
660#endif
661
662!==============================================================================
663! 4. Write "longitude" (data and attributes)
664!==============================================================================
665
666call def_var(nout,"longitude","East longitude","degrees_east",1,&
667             (/londimout/),nvarid,ierr)
668
669#ifdef NC_DOUBLE
670ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lon)
671#else
672ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
673#endif
674
675!==============================================================================
676! 5. Write "altitude" (data and attributes)
677!==============================================================================
678
679! Switch to netcdf define mode
680ierr = NF_REDEF (nout)
681
682#ifdef NC_DOUBLE
683ierr = NF_DEF_VAR (nout,"presnivs",NF_DOUBLE,1,altdimout,nvarid)
684#else
685ierr = NF_DEF_VAR (nout,"presnivs",NF_FLOAT,1,altdimout,nvarid)
686#endif
687
688ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",8,"presnivs")
689ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',2,"Pa")
690ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',4,"down")
691
692! End netcdf define mode
693ierr = NF_ENDDEF(nout)
694
695#ifdef NC_DOUBLE
696ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,alt)
697#else
698ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
699#endif
700
701!==============================================================================
702! 6. Write "controle" (data and attributes)
703!==============================================================================
704
705if (size(ctl).ne.0) then
706   ! Switch to netcdf define mode
707   ierr = NF_REDEF (nout)
708
709   #ifdef NC_DOUBLE
710   ierr = NF_DEF_VAR (nout,"controle",NF_DOUBLE,1,ctldimout,nvarid)
711   #else
712   ierr = NF_DEF_VAR (nout,"controle",NF_FLOAT,1,ctldimout,nvarid)
713   #endif
714
715   ierr = NF_PUT_ATT_TEXT (nout,nvarid,"title",18,"Control parameters")
716
717   ! End netcdf define mode
718   ierr = NF_ENDDEF(nout)
719
720   #ifdef NC_DOUBLE
721   ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,ctl)
722   #else
723   ierr = NF_PUT_VAR_REAL (nout,nvarid,ctl)
724   #endif
725endif
726
727end Subroutine initiate
728!******************************************************************************
729subroutine init2(infid,lonlen,latlen,altlen, &
730                 outfid,londimout,latdimout,altdimout)
731!==============================================================================
732! Purpose:
733! Copy phis() from input file to outpout file
734!==============================================================================
735! Remarks:
736! The NetCDF files must be open
737!==============================================================================
738
739implicit none
740
741include "netcdf.inc" ! NetCDF definitions
742
743!==============================================================================
744! Arguments:
745!==============================================================================
746integer, intent(in) :: infid  ! NetCDF output file ID
747integer, intent(in) :: lonlen ! # of grid points along longitude
748integer, intent(in) :: latlen ! # of grid points along latitude
749integer, intent(in) :: altlen ! # of grid points along latitude
750integer, intent(in) :: outfid ! NetCDF output file ID
751integer, intent(in) :: londimout ! longitude dimension ID
752integer, intent(in) :: latdimout ! latitude dimension ID
753integer, intent(in) :: altdimout ! altitude dimension ID
754!==============================================================================
755! Local variables:
756!==============================================================================
757real,dimension(:,:),allocatable :: phis ! Ground geopotential
758integer :: phisid
759integer :: ierr
760integer :: tmpvarid ! temporary variable ID
761logical :: phis_ok ! is "phis" available ?
762
763
764!==============================================================================
765! 1. Read data from input file
766!==============================================================================
767
768! ground geopotential phis
769ierr=NF_INQ_VARID(infid,"phis",tmpvarid)
770allocate(phis(lonlen,latlen))
771if (ierr.ne.NF_NOERR) then
772  write(*,*) "Failed to get phis ID. OK"
773  phis = 0.
774  phis_ok = .false.
775else
776  ierr=NF_GET_VAR_REAL(infid,tmpvarid,phis)
777  if (ierr.ne.NF_NOERR) then
778    stop "Error: Failed reading phis"
779  endif
780  phis_ok = .true.
781endif
782
783
784!==============================================================================
785! 2. Write
786!==============================================================================
787
788!==============================================================================
789! 2.2. phis()
790!==============================================================================
791
792
793IF (phis_ok) THEN
794
795!define phis
796 call def_var(nout,"phis","Ground level geopotential"," ",2,&
797            (/londimout,latdimout/),phisid,ierr)
798 if (ierr.ne.NF_NOERR) then
799     stop "Error: Failed to def_var phis"
800  endif
801
802! write phis
803#ifdef NC_DOUBLE
804ierr=NF_PUT_VAR_DOUBLE(outfid,phisid,phis)
805#else
806ierr=NF_PUT_VAR_REAL(outfid,phisid,phis)
807#endif
808if (ierr.ne.NF_NOERR) then
809  stop "Error: Failed to write phis"
810endif
811
812END IF
813
814
815! Cleanup
816deallocate(phis)
817
818end subroutine init2
819!******************************************************************************
820subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
821!==============================================================================
822! Purpose: Write a variable (i.e: add a variable to a dataset)
823! called "name"; along with its attributes "title", "units"...
824! to a file (following the NetCDF format)
825!==============================================================================
826! Remarks:
827! The NetCDF file must be open
828!==============================================================================
829
830implicit none
831
832include "netcdf.inc" ! NetCDF definitions
833
834!==============================================================================
835! Arguments:
836!==============================================================================
837integer, intent(in) :: nid
838! nid: [netcdf] file ID #
839character (len=*), intent(in) :: name
840! name(): [netcdf] variable's name
841character (len=*), intent(in) :: title
842! title(): [netcdf] variable's "title" attribute
843character (len=*), intent(in) :: units
844! unit(): [netcdf] variable's "units" attribute
845integer, intent(in) :: nbdim
846! nbdim: number of dimensions of the variable
847integer, dimension(nbdim), intent(in) :: dim
848! dim(nbdim): [netcdf] dimension(s) ID(s)
849integer, intent(out) :: nvarid
850! nvarid: [netcdf] ID # of the variable
851integer, intent(out) :: ierr
852! ierr: [netcdf] subroutines returned error code
853
854! Switch to netcdf define mode
855ierr=NF_REDEF(nid)
856
857! Insert the definition of the variable
858#ifdef NC_DOUBLE
859ierr=NF_DEF_VAR(nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
860#else
861ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
862#endif
863
864! Write the attributes
865ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",len_trim(adjustl(title)),adjustl(title))
866ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
867
868! End netcdf define mode
869ierr=NF_ENDDEF(nid)
870
871end subroutine def_var
872!******************************************************************************
873subroutine  missing_value(nout,nvarid,valid_range,missing)
874!==============================================================================
875! Purpose:
876! Write "valid_range" and "missing_value" attributes (of a given
877! variable) to a netcdf file
878!==============================================================================
879! Remarks:
880! NetCDF file must be open
881! Variable (nvarid) ID must be set
882!==============================================================================
883
884implicit none
885
886include "netcdf.inc"  ! NetCDF definitions
887
888!==============================================================================
889! Arguments:
890!==============================================================================
891integer, intent(in) :: nout
892! nout: [netcdf] file ID #
893integer, intent(in) :: nvarid
894! varid: [netcdf] variable ID #
895real, dimension(2), intent(in) :: valid_range
896! valid_range(2): [netcdf] "valid_range" attribute (min and max)
897real, intent(in) :: missing
898! missing: [netcdf] "missing_value" attribute
899
900!==============================================================================
901! Local variables:
902!==============================================================================
903integer :: ierr
904! ierr: [netcdf] subroutine returned error code
905!      INTEGER nout,nvarid,ierr
906!      REAL missing
907!      REAL valid_range(2)
908
909! Switch to netcdf dataset define mode
910ierr = NF_REDEF (nout)
911if (ierr.ne.NF_NOERR) then
912   print*,'missing_value: '
913   print*, NF_STRERROR(ierr)
914endif
915
916! Write valid_range() attribute
917#ifdef NC_DOUBLE
918ierr=NF_PUT_ATT_DOUBLE(nout,nvarid,'valid_range',NF_DOUBLE,2,valid_range)
919#else
920ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
921#endif
922
923if (ierr.ne.NF_NOERR) then
924   print*,'missing_value: valid_range attribution failed'
925   print*, NF_STRERROR(ierr)
926   !write(*,*) 'NF_NOERR', NF_NOERR
927   !CALL abort
928   stop ""
929endif
930
931! Write "missing_value" attribute
932#ifdef NC_DOUBLE
933ierr= NF_PUT_ATT_DOUBLE(nout,nvarid,'missing_value',NF_DOUBLE,1,missing)
934#else
935ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
936#endif
937
938if (ierr.NE.NF_NOERR) then
939   print*, 'missing_value: missing value attribution failed'
940   print*, NF_STRERROR(ierr)
941!    WRITE(*,*) 'NF_NOERR', NF_NOERR
942!    CALL abort
943   stop ""
944endif
945
946! End netcdf dataset define mode
947ierr = NF_ENDDEF(nout)
948
949end subroutine  missing_value
950!******************************************************************************
951
952end program localtime
953
954!*****************************************************************************
955subroutine interpolf(x,y,missing,xd,yd,nd)
956!==============================================================================
957! Purpose:
958! Yield y=f(x), where xd() end yd() are arrays of known values,
959! using linear interpolation
960! If x is not included in the interval spaned by xd(), then y is set
961! to a default value 'missing'
962! Note:
963! Array xd() should contain ordered (either increasing or decreasing) abscissas
964!==============================================================================
965implicit none
966!==============================================================================
967! Arguments:
968!==============================================================================
969real,intent(in) :: x ! abscissa at which interpolation is to be done
970real,intent(in) :: missing ! missing value (if no interpolation is performed)
971integer :: nd ! size of arrays
972real,dimension(nd),intent(in) :: xd ! array of known absissas
973real,dimension(nd),intent(in) :: yd ! array of correponding values
974
975real,intent(out) :: y ! interpolated value
976!==============================================================================
977! Local variables:
978!==============================================================================
979integer :: i
980
981! default: set y to 'missing'
982y=missing
983
984   do i=1,nd-1
985     if (((x.ge.xd(i)).and.(x.le.xd(i+1))).or.&
986          ((x.le.xd(i)).and.(x.ge.xd(i+1)))) then
987        y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i))
988        exit
989     endif
990   enddo
991
992
993end subroutine interpolf
994
Note: See TracBrowser for help on using the repository browser.