source: trunk/LMDZ.VENUS/Tools/localtime_mean_and_std.F90 @ 3094

Last change on this file since 3094 was 1442, checked in by slebonnois, 10 years ago

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