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

Last change on this file since 357 was 137, checked in by aslmd, 14 years ago

ATTENTION COMMIT MAJEUR NON-CONSERVATIF

  • CHANGEMENT ARBORESCENCE POUR SEPARATION CLAIRE DES COMPOSANTES et MODELES
  • UTILISATION DE LIENS SYMBOLIQUES POUR GARDER UNE BASE COMMUNE LMDZ.COMMON
  • VOIR LE FICHIER DOC/000-MODELS POUR EN SAVOIR PLUS !

EN CAS DE PROBLEME AVEC svn update IL EST CONSEILLE DE REVENIR A UN svn co

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