source: trunk/LMDZ.GENERIC/utilities/localtime.F90 @ 2613

Last change on this file since 2613 was 1905, checked in by mturbet, 7 years ago

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