source: trunk/LMDZ.MARS/util/streamfunction.F90 @ 1733

Last change on this file since 1733 was 644, checked in by tnavarro, 13 years ago

Added streamfunction utility in the svn repository. File from /u/forget/WWW/datagcm/Utilities/SOURCES/streamfunction.F90 with a small bug correction from EM.
File size: 35.0 KB
Line 
1program streamfunction
2
3! ----------------------------------------------------------------------------
4! Program to calculate the stream function and the total angular momentum
5! from stats and diagfi files
6! input : diagfi.nc  / concat.nc / stats.nc kind of files
7! author: F. Gonzalez-Galindo
8! ----------------------------------------------------------------------------
9
10implicit none
11
12include "netcdf.inc" ! NetCDF definitions
13
14character (len=128) :: infile ! input file name (diagfi.nc or stats.nc format)
15character (len=128) :: infile2 ! second input file (may be needed for 'aire' and 'cv')
16character (len=64) :: text ! to store some text
17
18character (len=50), dimension(:), allocatable :: var
19! var(): name(s) of variable(s) that will be concatenated
20character (len=100) :: filename
21! filename(): output file name
22integer :: nid,ierr,miss
23! nid: [netcdf] file ID #
24! ierr: [netcdf] subroutine returned error code
25! miss: [netcdf] subroutine returned error code
26integer :: tmpvarid
27! varid: temporary store a variable ID #
28integer tmpdimid ! temporarily store a dimension ID
29integer infid ! NetCDF input file ID (of diagfi.nc or stats.nc format)
30integer infid2 ! NetCDF input file which contains 'phisini' dataset (diagfi.nc)
31integer nbvarinfile ! # of variables in input file
32real, dimension(:), allocatable:: lat,lon,alt,time
33! lat(): array, stores latitude coordinates
34! lon(): array, stores longitude coordinates
35! alt(): array, stores altitude coordinates
36! time(): array, stores time coordinates
37integer :: ilat,ilon,ialt,it,i
38! ilat: index for loops on latitude
39! ilon: index for loops on longitude
40! ialt: index for loops on altitude
41! it: # index for loops on time
42! i: index for loops
43integer :: nout,latdimout,londimout,altdimout,timedimout,lonvarout,timevarout
44integer :: psivarout,momvarout,uvarout,vvarout,rhovarout,tempvarout
45integer :: psvarout, phisinitvarout
46! nout: [netcdf] output file ID
47! latdimout: [netcdf] output latitude (dimension) ID
48! londimout: [netcdf] output longitude (dimension) ID
49! altdimout: [netcdf] output altitude (dimension) ID
50! timedimout: [netcdf] output time (dimension) ID
51! lonvarout: [netcdf] ID of output "Time" variable
52integer lonlength ! # of grid points along longitude
53integer latlength ! # of grid points along latitude
54integer altlength ! # of grid point along altitude (of input datasets)
55integer timelength ! # of points along time
56real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
57real,dimension(:),allocatable :: sigma ! sigma levels
58real,dimension(:),allocatable :: lon_fake ! Fake longitude to be written
59real,dimension(:,:),allocatable :: aire ! Area of the box
60real,dimension(:,:),allocatable :: cv ! Conversion between natural and covariant v
61real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
62real,dimension(:,:,:),allocatable :: ps ! GCM surface pressure
63real,dimension(:,:,:,:),allocatable :: press ! GCM atmospheric pressure
64real,dimension(:,:,:,:),allocatable :: temp ! GCM atmospheric temperature
65real,dimension(:,:,:,:),allocatable :: u ! GCM zonal wind
66real,dimension(:,:,:,:),allocatable :: v ! GCM meridional wind
67real,dimension(:,:,:,:),allocatable :: rho ! GCM atmospheric density
68real,dimension(:,:,:),allocatable :: vcont ! Covariant meridional wind
69real,dimension(:,:),allocatable :: vcontcum ! Zonal mean covariant meridional wind
70real,dimension(:,:,:),allocatable :: plev ! Pressure from hybrid coordinates
71real,dimension(:,:,:),allocatable :: mass ! Mass in a grid box
72real,dimension(:,:,:),allocatable :: dmass ! Mass variation
73real,dimension(:,:),allocatable :: psi ! Stream function
74real,dimension(:,:,:),allocatable :: mom ! Angular momentum
75real,dimension(:,:),allocatable :: momave ! Zonally averaged angular momentum
76real,dimension(:,:,:),allocatable :: ucum ! Temporally averaged zonal wind
77real,dimension(:,:,:),allocatable :: vcum ! Temporally averaged meridional wind
78real,dimension(:,:,:),allocatable :: rhocum ! Temporally averaged density
79real,dimension(:,:,:),allocatable :: tempcum ! Temporally averaged zonal wind
80real,dimension(:,:),allocatable :: pscum ! Temporally averaged zonal wind
81real,dimension(:,:),allocatable :: uzm ! Zonally averaged zonal wind
82real,dimension(:,:),allocatable :: vzm ! Zonally averaged meridional wind
83real,dimension(:,:),allocatable :: rhozm ! Zonally averaged density
84real,dimension(:,:),allocatable :: tempzm ! Zonally averaged temperature
85real,dimension(:),allocatable :: pszm ! Zonally averaged surface pressure
86real,dimension(:),allocatable :: phisinitzm ! Zonally averaged ground geopotential
87
88!Parameters to calculate the stream function
89real :: g0 ! exact mean gravity at radius=3396.km (Lemoine et al. 2001)
90data    g0 /3.7257964/
91real :: a0 
92data    a0 /3396000/  ! radius=3396.km
93real :: daysec
94data    daysec /88775/ ! Duration of day=88775 s
95real :: omega
96real :: Rgas   ! gas constant
97data    Rgas /190./
98
99
100!===============================================================================
101! 1.1 Input file
102!===============================================================================
103
104write(*,*) ""
105write(*,*) " Program valid for diagfi.nc, concatnc.nc and stats.nc files"
106write(*,*) "Enter input file name:"
107
108read(*,'(a128)') infile
109write(*,*) ""
110
111! open input file
112
113ierr = NF_OPEN(infile,NF_NOWRITE,infid)
114if (ierr.ne.NF_NOERR) then
115   write(*,*) 'ERROR: Pb opening input file'
116   stop ""
117endif
118
119!===============================================================================
120! 1.2 Get grids in lon,lat,alt,time,
121!     as well as hybrid coordinates aps() and bps() (or sigma levels sigma())
122!     aire() and cv() from input file
123!===============================================================================
124
125! 1.2.1 longitude, latitude, altitude and time
126
127! latitude
128ierr=NF_INQ_DIMID(infid,"latitude",tmpdimid)
129if (ierr.ne.NF_NOERR) then
130  stop "Error: Failed to get latitude dimension ID"
131else
132  ierr=NF_INQ_VARID(infid,"latitude",tmpvarid)
133  if (ierr.ne.NF_NOERR) then
134    stop "Error: Failed to get latitude ID"
135  else
136    ierr=NF_INQ_DIMLEN(infid,tmpdimid,latlength)
137    if (ierr.ne.NF_NOERR) then
138      stop "Error: Failed to get latitude length"
139    else
140      allocate(lat(latlength))
141      ierr=NF_GET_VAR_REAL(infid,tmpvarid,lat)
142      if (ierr.ne.NF_NOERR) then
143        stop "Error: Failed reading latitude"
144      endif
145    endif
146  endif
147endif
148
149! longitude
150ierr=NF_INQ_DIMID(infid,"longitude",tmpdimid)
151if (ierr.ne.NF_NOERR) then
152  stop "Error: Failed to get longitude dimension ID"
153else
154  ierr=NF_INQ_VARID(infid,"longitude",tmpvarid)
155  if (ierr.ne.NF_NOERR) then
156    stop "Error: Failed to get longitude ID"
157  else
158    ierr=NF_INQ_DIMLEN(infid,tmpdimid,lonlength)
159    if (ierr.ne.NF_NOERR) then
160      stop "Error: Failed to get longitude length"
161    else
162      allocate(lon(lonlength))
163      ierr=NF_GET_VAR_REAL(infid,tmpvarid,lon)
164      if (ierr.ne.NF_NOERR) then
165        stop "Error: Failed reading longitude"
166      endif
167    endif
168  endif
169endif
170
171! time
172ierr=NF_INQ_DIMID(infid,"Time",tmpdimid)
173if (ierr.ne.NF_NOERR) then
174  stop "Error: Failed to get Time dimension ID"
175else
176  ierr=NF_INQ_VARID(infid,"Time",tmpvarid)
177  if (ierr.ne.NF_NOERR) then
178    stop "Error: Failed to get Time ID"
179  else
180    ierr=NF_INQ_DIMLEN(infid,tmpdimid,timelength)
181    if (ierr.ne.NF_NOERR) then
182      stop "Error: Failed to get Time length"
183    else
184      allocate(time(timelength))
185      ierr=NF_GET_VAR_REAL(infid,tmpvarid,time)
186      if (ierr.ne.NF_NOERR) then
187        stop "Error: Failed reading Time"
188      endif
189    endif
190  endif
191endif
192
193! altlength
194ierr=NF_INQ_DIMID(infid,"altitude",tmpdimid)
195if (ierr.ne.NF_NOERR) then
196  stop "Error: Failed to get altitude dimension ID"
197else
198  ierr=NF_INQ_VARID(infid,"altitude",tmpvarid)
199  if (ierr.ne.NF_NOERR) then
200     stop "Error: Failed to get altitude length"
201  else
202     ierr=NF_INQ_DIMLEN(infid,tmpdimid,altlength)
203     if (ierr.ne.NF_NOERR) then
204        stop "Error: Failed to get altitude length"
205     else   
206        allocate(alt(altlength))
207        ierr=NF_GET_VAR_REAL(infid,tmpvarid,alt)
208        if (ierr.ne.NF_NOERR) then
209           stop "Error: Failed reading Altitude"
210        endif
211     endif
212  endif
213endif
214
215! 1.2.2 Get hybrid coordinates
216
217! look for hybrid coordinates
218
219! hybrid coordinate aps
220ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
221if (ierr.ne.NF_NOERR) then
222   stop "Error: Failed to get aps ID"
223else
224   allocate(aps(altlength))
225   ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
226   if (ierr.ne.NF_NOERR) then
227      stop "Error: Failed reading aps"
228   endif
229endif
230
231! hybrid coordinate bps
232ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
233if (ierr.ne.NF_NOERR) then
234   stop "Error: Failed to get bps ID"
235else
236   allocate(bps(altlength))
237   ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
238   if (ierr.ne.NF_NOERR) then
239      stop "Error: Failed reading bps"
240   endif
241endif
242
243!surface pressure
244ierr=NF_INQ_VARID(infid,"ps",tmpvarid)
245if (ierr.ne.NF_NOERR) then
246   stop "Error: Failed to get ps ID"
247else
248   allocate(ps(lonlength,latlength,timelength))
249   ierr=NF_GET_VAR_REAL(infid,tmpvarid,ps)
250   if (ierr.ne.NF_NOERR) then
251      stop "Error: Failed reading ps"
252   endif
253endif
254
255
256!===============================================================================
257! 1.3 Reads variables in input file
258!===============================================================================
259
260ierr=NF_INQ_NVARS(infid,nbvarinfile)
261if (ierr.ne.NF_NOERR) then
262  write(*,*) 'ERROR: Failed geting number of variables from file'
263  stop
264endif
265
266!1.3.1 Zonal wind, meridional wind
267
268!Zonal wind
269ierr=NF_INQ_VARID(infid,"u",tmpvarid)
270if (ierr.ne.NF_NOERR) then
271  stop "Error: Failed to get u ID"
272else
273  allocate(u(lonlength,latlength,altlength,timelength))
274  ierr=NF_GET_VAR_REAL(infid,tmpvarid,u)
275  if (ierr.ne.NF_NOERR) then
276    stop "Error: Failed reading zonal wind"
277  endif
278endif
279
280!Meridional wind
281ierr=NF_INQ_VARID(infid,"v",tmpvarid)
282if (ierr.ne.NF_NOERR) then
283  stop "Error: Failed to get v ID"
284else
285  allocate(v(lonlength,latlength,altlength,timelength))
286  ierr=NF_GET_VAR_REAL(infid,tmpvarid,v)
287  if (ierr.ne.NF_NOERR) then
288    stop "Error: Failed reading zonal wind"
289  endif
290endif
291
292
293! 1.3.2 aire
294
295allocate(aire(lonlength,latlength))
296! look for 'aire' in current file
297ierr=NF_INQ_VARID(infid,"aire",tmpvarid)
298if (ierr.ne.NF_NOERR) then
299  write(*,*) "Warning: Failed to get aire ID from file ",trim(infile)
300  infile2="diagfi.nc"
301  write(*,*) "         Trying file ",trim(infile2)
302  ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
303  if (ierr.ne.NF_NOERR) then
304     infile2="diagfi1.nc"
305     write(*,*) "         Trying file ",trim(infile2)
306     ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
307     if (ierr.ne.NF_NOERR) then
308        write(*,*) "Problem: Could not find/open these files"
309        stop "Might as well stop here"
310     endif
311  endif
312
313   ! Get ID for aire
314   ierr=NF_INQ_VARID(infid2,"aire",tmpvarid)
315   if (ierr.ne.NF_NOERR) then
316      stop "Error: Failed to get aire ID"
317   endif
318   ! Get aire
319   ierr=NF_GET_VAR_REAL(infid2,tmpvarid,aire)
320   if (ierr.ne.NF_NOERR) then
321      stop "Error: Failed reading aire"
322   endif
323   ! Close file
324   write(*,*) 'OK, got aire'
325   ierr=NF_CLOSE(infid2)
326else
327   ierr=NF_GET_VAR_REAL(infid,tmpvarid,aire)
328   if (ierr.ne.NF_NOERR) then
329      stop "Error: Failed reading aire"
330   endif
331endif
332
333
334! 1.3.3 phisinit
335
336allocate(phisinit(lonlength,latlength))
337! look for '' in current file
338ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
339if (ierr.ne.NF_NOERR) then
340  write(*,*) "Warning: Failed to get phisinit ID from file ",trim(infile)
341  infile2="diagfi.nc"
342  write(*,*) "         Trying file ",trim(infile2)
343  ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
344  if (ierr.ne.NF_NOERR) then
345     infile2="diagfi1.nc"
346     write(*,*) "         Trying file ",trim(infile2)
347     ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
348     if (ierr.ne.NF_NOERR) then
349        write(*,*) "Problem: Could not find/open these files"
350        stop "Might as well stop here"
351     endif
352  endif
353
354
355   ! Get ID for phisinit
356   ierr=NF_INQ_VARID(infid2,"phisinit",tmpvarid)
357   if (ierr.ne.NF_NOERR) then
358      stop "Error: Failed to get phisinit ID"
359   endif
360   ! Get phisinit
361   ierr=NF_GET_VAR_REAL(infid2,tmpvarid,phisinit)
362   if (ierr.ne.NF_NOERR) then
363      stop "Error: Failed reading phisinit"
364   endif
365   ! Close file
366   write(*,*) 'OK, got phisinit'
367   ierr=NF_CLOSE(infid2)
368else
369   ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
370   if (ierr.ne.NF_NOERR) then
371      stop "Error: Failed reading phisinit"
372   endif
373endif
374
375
376! 1.3.4 cv
377
378allocate(cv(lonlength,latlength))
379! look for 'cv' in current file
380ierr=NF_INQ_VARID(infid,"cv",tmpvarid)
381if (ierr.ne.NF_NOERR) then
382  write(*,*) "Warning: Failed to get cv ID from file ",trim(infile)
383  infile2="diagfi.nc"
384  write(*,*) "         Trying file ",trim(infile2)
385  ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
386  if (ierr.ne.NF_NOERR) then
387     infile2="diagfi1.nc"
388     write(*,*) "         Trying file ",trim(infile2)
389     ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
390     if (ierr.ne.NF_NOERR) then
391        write(*,*) "Problem: Could not find/open these files"
392        stop "Might as well stop here"
393     endif
394  endif
395
396
397   ! Get ID for cv
398   ierr=NF_INQ_VARID(infid2,"cv",tmpvarid)
399   if (ierr.ne.NF_NOERR) then
400      stop "Error: Failed to get cv ID"
401   endif
402   ! Get cv
403   ierr=NF_GET_VAR_REAL(infid2,tmpvarid,cv)
404   if (ierr.ne.NF_NOERR) then
405      stop "Error: Failed reading cv"
406   endif
407   ! Close file
408   write(*,*) 'OK, got cv'
409   ierr=NF_CLOSE(infid2)
410else
411   ierr=NF_GET_VAR_REAL(infid,tmpvarid,cv)
412   if (ierr.ne.NF_NOERR) then
413      stop "Error: Failed reading cv"
414   endif
415endif
416
417
418!Other variables: rho, temp
419
420ierr=NF_INQ_VARID(infid,"temp",tmpvarid)
421if (ierr.ne.NF_NOERR) then
422  stop "Error: Failed to get temp ID"
423else
424  allocate(temp(lonlength,latlength,altlength,timelength))
425  ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp)
426  if (ierr.ne.NF_NOERR) then
427    stop "Error: Failed reading temperature"
428  endif
429endif
430
431ierr=NF_INQ_VARID(infid,"rho",tmpvarid)
432if (ierr.ne.NF_NOERR) then
433  write(*,*) "Error: Failed to get rho ID"
434  write(*,*) "Let's compute it from temperature"
435  allocate(rho(lonlength,latlength,altlength,timelength))
436 
437  do it=1,timelength
438    do ilat=1,latlength
439      do ilon=1,lonlength
440        do ialt=1,altlength
441          rho(ilon,ilat,ialt,it)= (aps(ialt)+bps(ialt)*ps(ilon,ilat,it))/(Rgas*temp(ilon,ilat,ialt,it))
442        enddo
443      enddo
444    enddo
445  enddo
446else
447  allocate(rho(lonlength,latlength,altlength,timelength))
448  ierr=NF_GET_VAR_REAL(infid,tmpvarid,rho)
449  if (ierr.ne.NF_NOERR) then
450    stop "Error: Failed reading mass density"
451  endif
452endif
453
454!==============================================================================
455! 1.4. Output file name
456!==============================================================================
457filename=infile(1:len_trim(infile)-3)//"_stream.nc"
458write(*,*)'The output file has name:'
459write(*,*) filename
460
461
462!==============================================================================
463! 2.1. Temporal averages
464!==============================================================================
465
466allocate(ucum(lonlength,latlength,altlength)) !Time mean zonal wind
467allocate(vcum(lonlength,latlength,altlength)) !Time mean meridional wind
468allocate(tempcum(lonlength,latlength,altlength)) !Time mean temperature
469allocate(rhocum(lonlength,latlength,altlength)) !Time mean density
470allocate(pscum(lonlength,latlength)) !Time mean surface pressure
471do ilat=1,latlength
472   do ilon=1,lonlength
473      pscum(ilon,ilat)=0.
474      do ialt=1,altlength
475            ucum(ilon,ilat,ialt)=0.
476            vcum(ilon,ilat,ialt)=0.
477            tempcum(ilon,ilat,ialt)=0.
478            rhocum(ilon,ilat,ialt)=0.
479            do it=1,timelength
480               ucum(ilon,ilat,ialt)=ucum(ilon,ilat,ialt)+u(ilon,ilat,ialt,it)
481               vcum(ilon,ilat,ialt)=vcum(ilon,ilat,ialt)+v(ilon,ilat,ialt,it)
482               tempcum(ilon,ilat,ialt)=tempcum(ilon,ilat,ialt)+temp(ilon,ilat,ialt,it)
483               rhocum(ilon,ilat,ialt)=rhocum(ilon,ilat,ialt)+rho(ilon,ilat,ialt,it)
484            enddo
485            ucum(ilon,ilat,ialt)=ucum(ilon,ilat,ialt)/timelength
486            vcum(ilon,ilat,ialt)=vcum(ilon,ilat,ialt)/timelength
487            tempcum(ilon,ilat,ialt)=tempcum(ilon,ilat,ialt)/timelength
488            rhocum(ilon,ilat,ialt)=rhocum(ilon,ilat,ialt)/timelength
489      enddo
490      do it=1,timelength
491         pscum(ilon,ilat)=pscum(ilon,ilat)+ps(ilon,ilat,it)
492      enddo
493      pscum(ilon,ilat)=pscum(ilon,ilat)/timelength
494   enddo
495enddo
496
497!==============================================================================
498! 2.2. Calculation of the stream function
499!==============================================================================
500
501!Contravariant meridional wind
502allocate(vcont(lonlength,latlength,altlength))
503do ilon=1,lonlength
504   do ialt=1,altlength
505      do ilat=1,latlength-1
506         vcont(ilon,ilat,ialt)=0.5*&
507              (vcum(ilon,ilat,ialt)+vcum(ilon,ilat+1,ialt))/cv(ilon,ilat)
508      enddo
509      vcont(ilon,latlength,ialt)=0.!vcont(ilon,latlength-1,ialt)
510   enddo
511enddo
512
513!Pressure from hybrid levels
514allocate(Plev(lonlength,latlength,altlength))
515do ilon=1,lonlength
516   do ilat=1,latlength
517      do ialt=1,altlength
518         Plev(ilon,ilat,ialt)=aps(ialt)+bps(ialt)*pscum(ilon,ilat)
519      enddo
520   enddo
521enddo
522
523!Mass in the box
524allocate(mass(lonlength,latlength,altlength))
525do ilon=1,lonlength
526   do ilat=1,latlength
527      do ialt=1,altlength
528         mass(ilon,ilat,ialt)=aire(ilon,ilat)*&
529              Plev(ilon,ilat,ialt)/g0
530      enddo
531   enddo
532enddo
533
534!Mass variation in the box
535allocate(dmass(lonlength,latlength,altlength))
536do ilon=1,lonlength
537   do ilat=1,latlength
538      do ialt=1,altlength-1
539         dmass(ilon,ilat,ialt)=mass(ilon,ilat,ialt)-&
540              mass(ilon,ilat,ialt+1)
541      enddo
542      dmass(ilon,ilat,altlength)=0.!dmass(ilon,ilat,altlength-1)
543   enddo
544enddo
545
546!Stream function
547allocate(psi(latlength,altlength))
548allocate(vcontcum(latlength,altlength))
549do ilat=1,latlength-1
550   psi(ilat,altlength)=0.
551   do ialt=altlength-1,1,-1
552      psi(ilat,ialt)=psi(ilat,ialt+1)
553      vcontcum(ilat,ialt)=0.
554      do ilon=1,lonlength
555         vcontcum(ilat,ialt)=vcontcum(ilat,ialt)+vcont(ilon,ilat,ialt)/lonlength
556         psi(ilat,ialt)=psi(ilat,ialt)-(vcont(ilon,ilat,ialt)&
557              *(dmass(ilon,ilat,ialt)+dmass(ilon,ilat+1,ialt)))
558      enddo
559      psi(latlength,ialt)=0.!psi(latlength-1,ialt)
560   enddo
561enddo
562
563!==============================================================================
564! 2.3. Calculation of the total angular momentum
565!==============================================================================
566
567!Rotation speed
568omega=2.*3.141592/daysec
569
570!Angular momentum
571allocate(mom(lonlength,latlength,altlength))
572allocate(momave(latlength,altlength))
573do ilat=1,latlength
574   do ialt=1,altlength
575      momave(ilat,ialt)=0.
576      do ilon=1,lonlength
577!            mom(ilon,ilat,ialt,it)=dmass(ilon,ilat,ialt,it)*a0*&
578         mom(ilon,ilat,ialt)=a0*cos(lat(ilat)*3.141592/180.)*&
579              (omega*a0*cos(lat(ilat)*3.141592/180.)+ucum(ilon,ilat,ialt))
580         
581         momave(ilat,ialt)=momave(ilat,ialt)+mom(ilon,ilat,ialt)
582      enddo
583      momave(ilat,ialt)=momave(ilat,ialt)/lonlength
584   enddo
585enddo
586
587!==============================================================================
588! 2.4. Zonal mean winds, temperature and density
589!==============================================================================
590
591allocate(uzm(latlength,altlength)) !Zonal mean zonal wind
592allocate(vzm(latlength,altlength)) !Zonal mean meridional wind
593allocate(tempzm(latlength,altlength)) !Zonal mean temperature
594allocate(rhozm(latlength,altlength)) !Zonal mean density
595allocate(pszm(latlength)) !Zonal mean surface pressure
596allocate(phisinitzm(latlength)) !Zonal mean ground geopotential
597do ilat=1,latlength
598   phisinitzm(ilat)=0.
599   pszm(ilat)=0.
600   do ialt=1,altlength
601      uzm(ilat,ialt)=0.
602      vzm(ilat,ialt)=0.
603      tempzm(ilat,ialt)=0.
604      rhozm(ilat,ialt)=0.
605      do ilon=1,lonlength
606         uzm(ilat,ialt)=uzm(ilat,ialt)+ucum(ilon,ilat,ialt)
607         vzm(ilat,ialt)=vzm(ilat,ialt)+vcum(ilon,ilat,ialt)
608         tempzm(ilat,ialt)=tempzm(ilat,ialt)+tempcum(ilon,ilat,ialt)
609         rhozm(ilat,ialt)=rhozm(ilat,ialt)+rhocum(ilon,ilat,ialt)
610      enddo
611      uzm(ilat,ialt)=uzm(ilat,ialt)/lonlength
612      vzm(ilat,ialt)=vzm(ilat,ialt)/lonlength
613      tempzm(ilat,ialt)=tempzm(ilat,ialt)/lonlength
614      rhozm(ilat,ialt)=rhozm(ilat,ialt)/lonlength
615   enddo
616   do ilon=1,lonlength
617      pszm(ilat)=pszm(ilat)+pscum(ilon,ilat)
618      phisinitzm(ilat)=phisinitzm(ilat)+phisinit(ilon,ilat)
619   enddo
620   pszm(ilat)=pszm(ilat)/lonlength
621   phisinitzm(ilat)=phisinitzm(ilat)/lonlength 
622enddo
623
624
625!==============================================================================
626! 2.3. Recalculation of angular momentum using zonally averaged wind
627!==============================================================================
628
629do ilat=1,latlength
630   do ialt=1,altlength
631      momave(ilat,ialt)=a0*cos(lat(ilat)*3.141592/180.)*&
632           (omega*a0*cos(lat(ilat)*3.141592/180.)+uzm(ilat,ialt))
633   enddo
634enddo
635
636
637!==============================================================================
638! 3.1 Dimensions in output file
639!==============================================================================
640
641! 3.1.1 Initialize output file's lat,lon,alt and time dimensions
642call initiate (filename,lat,alt,time,nout,&
643     latdimout,londimout,altdimout,timedimout,lonvarout,timevarout)
644! Initialize output file's aps,bps variables
645call init2(infid,lonlength,latlength,altlength,&
646     nout,londimout,latdimout,altdimout)
647
648! 3.1.2 New longitude dimension/value in output file
649allocate(lon_fake(1))
650lon_fake(1)=0.
651#ifdef NC_DOUBLE
652        ierr= NF_PUT_VARA_DOUBLE(nout,lonvarout,1,1,lon_fake)
653#else
654        ierr= NF_PUT_VARA_REAL(nout,lonvarout,1,1,lon_fake)
655#endif
656
657! 3.1.3 New time dimension/value in output file
658#ifdef NC_DOUBLE
659        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,1,1,lon_fake)
660#else
661        ierr= NF_PUT_VARA_REAL(nout,timevarout,1,1,lon_fake)
662#endif
663
664
665!==============================================================================
666! 3.2 Write variables
667!==============================================================================
668
669!Define the dimensions of the variables to be written
670
671!!3.2.1 Stream function
672call def_var(nout,"psi","Stream function","",4,&
673             (/londimout,latdimout,altdimout,timedimout/),psivarout,ierr)
674if (ierr.ne.NF_NOERR) then
675   write(*,*) 'Error, could not define variable psi'
676   stop ""
677endif
678
679!Write in the output file
680ierr=NF_PUT_VAR_REAL(nout,psivarout,psi)
681if (ierr.ne.NF_NOERR) then
682   write(*,*)'Error, Failed to write variable psi'
683   stop
684endif
685
686!3.2.2 Momentum
687call def_var(nout,"momave","Angular momentum","",4,&
688             (/londimout,latdimout,altdimout,timedimout/),momvarout,ierr)
689if (ierr.ne.NF_NOERR) then
690   write(*,*) 'Error, could not define variable momave'
691   stop ""
692endif
693
694!Write in the output file
695ierr=NF_PUT_VAR_REAL(nout,momvarout,momave)
696if (ierr.ne.NF_NOERR) then
697   write(*,*)'Error, Failed to write variable momave'
698   stop
699endif
700
701
702!3.2.3 Zonal mean zonal wind
703call def_var(nout,"u","Zonal mean zonal wind","m/s",4,&
704             (/londimout,latdimout,altdimout,timedimout/),uvarout,ierr)
705if (ierr.ne.NF_NOERR) then
706   write(*,*) 'Error, could not define variable u'
707   stop ""
708endif
709
710!Write in the output file
711ierr=NF_PUT_VAR_REAL(nout,uvarout,uzm)
712if (ierr.ne.NF_NOERR) then
713   write(*,*)'Error, Failed to write variable u'
714   stop
715endif
716
717
718!3.2.4 Zonal mean meridional wind
719call def_var(nout,"v","Zonal mean meridional wind","m/s",4,&
720             (/londimout,latdimout,altdimout,timedimout/),vvarout,ierr)
721if (ierr.ne.NF_NOERR) then
722   write(*,*) 'Error, could not define variable v'
723   stop ""
724endif
725
726!Write in the output file
727ierr=NF_PUT_VAR_REAL(nout,vvarout,vzm)
728if (ierr.ne.NF_NOERR) then
729   write(*,*)'Error, Failed to write variable v'
730   stop
731endif
732
733!3.2.5 Zonal mean density
734call def_var(nout,"rho","Zonal mean atmospheric density","",4,&
735             (/londimout,latdimout,altdimout,timedimout/),rhovarout,ierr)
736if (ierr.ne.NF_NOERR) then
737   write(*,*) 'Error, could not define variable rho'
738   stop ""
739endif
740
741!Write in the output file
742ierr=NF_PUT_VAR_REAL(nout,rhovarout,rhozm)
743if (ierr.ne.NF_NOERR) then
744   write(*,*)'Error, Failed to write variable rho'
745   stop
746endif
747
748!3.2.6 Zonal mean temperature
749call def_var(nout,"temp","Zonal mean temperature","K",4,&
750             (/londimout,latdimout,altdimout,timedimout/),tempvarout,ierr)
751if (ierr.ne.NF_NOERR) then
752   write(*,*) 'Error, could not define variable temp'
753   stop ""
754endif
755
756!Write in the output file
757ierr=NF_PUT_VAR_REAL(nout,tempvarout,tempzm)
758if (ierr.ne.NF_NOERR) then
759   write(*,*)'Error, Failed to write variable temp'
760   stop
761endif
762
763!3.2.7 Zonal mean surface pressure
764call def_var(nout,"ps","Zonal mean surface pressure","Pa",3,&
765             (/londimout,latdimout,timedimout/),psvarout,ierr)
766if (ierr.ne.NF_NOERR) then
767   write(*,*) 'Error, could not define variable ps'
768   stop ""
769endif
770
771!Write in the output file
772ierr=NF_PUT_VAR_REAL(nout,psvarout,pszm)
773if (ierr.ne.NF_NOERR) then
774   write(*,*)'Error, Failed to write variable ps'
775   stop
776endif
777
778
779!3.2.8 Zonal mean geopotential
780call def_var(nout,"phisinit","Zonal mean initial geopotential","",2,&
781             (/londimout,latdimout/),phisinitvarout,ierr)
782if (ierr.ne.NF_NOERR) then
783   write(*,*) 'Error, could not define variable phisinit'
784   stop ""
785endif
786
787!Write in the output file
788ierr=NF_PUT_VAR_REAL(nout,phisinitvarout,phisinitzm)
789if (ierr.ne.NF_NOERR) then
790   write(*,*)'Error, Failed to write variable phisinit'
791   stop
792endif
793
794! Close input file
795ierr=nf_close(nid)
796
797! Close output file
798ierr=NF_CLOSE(nout)
799
800contains
801
802!******************************************************************************
803Subroutine initiate (filename,lat,alt,time,&
804                     nout,latdimout,londimout,altdimout,timedimout,&
805                     lonvarout,timevarout)
806!==============================================================================
807! Purpose:
808! Create and initialize a data file (NetCDF format)
809!==============================================================================
810! Remarks:
811! The NetCDF file (created in this subroutine) remains open
812!==============================================================================
813
814implicit none
815
816include "netcdf.inc" ! NetCDF definitions
817
818!==============================================================================
819! Arguments:
820!==============================================================================
821character (len=*), intent(in):: filename
822! filename(): the file's name
823real, dimension(:), intent(in):: lat
824! lat(): latitude
825real, dimension(:), intent(in):: alt
826! alt(): altitude
827real, dimension(:), intent(in):: time
828! time(): Time
829integer, intent(out):: nout
830! nout: [netcdf] file ID
831integer, intent(out):: latdimout
832! latdimout: [netcdf] lat() (i.e.: latitude)  ID
833integer, intent(out):: londimout
834! londimout: [netcdf] lon()  ID
835integer, intent(out):: altdimout
836! altdimout: [netcdf] alt()  ID
837integer, intent(out):: timedimout
838! timedimout: [netcdf] "Time"  ID
839integer, intent(out):: lonvarout
840! timevarout: [netcdf] Longiture (considered as a variable) ID
841integer, intent(out):: timevarout
842! timevarout: [netcdf] Time (considered as a variable) ID
843
844!==============================================================================
845! Local variables:
846!==============================================================================
847!integer :: latdim,londim,altdim,timedim
848integer :: nvarid,ierr
849! nvarid: [netcdf] ID of a variable
850! ierr: [netcdf]  return error code (from called subroutines)
851
852!==============================================================================
853! 1. Create (and open) output file
854!==============================================================================
855write(*,*) "creating "//trim(adjustl(filename))//'...'
856ierr = NF_CREATE(filename,NF_CLOBBER,nout)
857! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
858if (ierr.NE.NF_NOERR) then
859   WRITE(*,*)'ERROR: Impossible to create the file.'
860   stop ""
861endif
862
863!==============================================================================
864! 2. Define/write "dimensions" and get their IDs
865!==============================================================================
866
867ierr = NF_DEF_DIM(nout, "latitude", size(lat), latdimout)
868!ierr = NF_DEF_DIM(nout, "longitude", NF_UNLIMITED, londimout)
869ierr = NF_DEF_DIM(nout, "longitude", 1, londimout)
870ierr = NF_DEF_DIM(nout, "altitude", size(alt), altdimout)
871ierr = NF_DEF_DIM(nout, "Time", 1, timedimout)
872
873! End netcdf define mode
874ierr = NF_ENDDEF(nout)
875
876!==============================================================================
877! 3. Write "Time" (attributes)
878!==============================================================================
879
880call def_var(nout,"Time","Time","years since 0000-00-0 00:00:00",1,&
881             (/timedimout/),timevarout,ierr)
882
883!==============================================================================
884! 4. Write "latitude" (data and attributes)
885!==============================================================================
886
887call def_var(nout,"latitude","latitude","degrees_north",1,&
888             (/latdimout/),nvarid,ierr)
889
890#ifdef NC_DOUBLE
891ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lat)
892#else
893ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
894#endif
895
896!==============================================================================
897! 4. Write "longitude" (attributes)
898!==============================================================================
899
900call def_var(nout,"longitude","East longitude","degrees_east",1,&
901             (/londimout/),lonvarout,ierr)
902
903
904!==============================================================================
905! 4. Write "altitude" (data and attributes)
906!==============================================================================
907
908! Switch to netcdf define mode
909
910call def_var(nout,"altitude","Altitude","km",1,&
911             (/altdimout/),nvarid,ierr)
912
913#ifdef NC_DOUBLE
914ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,alt)
915#else
916ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
917#endif
918
919end Subroutine initiate
920!******************************************************************************
921subroutine init2(infid,lonlen,latlen,altlen, &
922                 outfid,londimout,latdimout,altdimout)
923!==============================================================================
924! Purpose:
925! Copy aps() and bps() from input file to outpout file
926!==============================================================================
927! Remarks:
928! The NetCDF files must be open
929!==============================================================================
930
931implicit none
932
933include "netcdf.inc" ! NetCDF definitions
934
935!==============================================================================
936! Arguments:
937!==============================================================================
938integer, intent(in) :: infid  ! NetCDF output file ID
939integer, intent(in) :: lonlen ! # of grid points along longitude
940integer, intent(in) :: latlen ! # of grid points along latitude
941integer, intent(in) :: altlen ! # of grid points along latitude
942integer, intent(in) :: outfid ! NetCDF output file ID
943integer, intent(in) :: londimout ! longitude dimension ID
944integer, intent(in) :: latdimout ! latitude dimension ID
945integer, intent(in) :: altdimout ! altitude dimension ID
946!==============================================================================
947! Local variables:
948!==============================================================================
949real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
950integer :: apsid,bpsid
951integer :: ierr
952integer :: tmpvarid ! temporary variable ID
953logical :: aps_ok, bps_ok ! are "phisinit" "aps" "bps" available ?
954
955
956!==============================================================================
957! 1. Read data from input file
958!==============================================================================
959
960! hybrid coordinate aps
961  allocate(aps(altlen))
962ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
963if (ierr.ne.NF_NOERR) then
964  write(*,*) "oops Failed to get aps ID. OK"
965  aps_ok=.false.
966else
967  ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
968  if (ierr.ne.NF_NOERR) then
969   stop "error: Failed reading aps"
970  endif
971  aps_ok=.true.
972endif
973
974! hybrid coordinate bps
975  allocate(bps(altlen))
976ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
977if (ierr.ne.NF_NOERR) then
978  write(*,*) "oops: Failed to get bps ID. OK"
979  bps_ok=.false.
980else
981  ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
982  if (ierr.ne.NF_NOERR) then
983    stop "Error: Failed reading bps"
984  endif
985  bps_ok=.true.
986endif
987
988!==============================================================================
989! 2. Write
990!==============================================================================
991
992!==============================================================================
993! 2.2. Hybrid coordinates aps() and bps()
994!==============================================================================
995
996IF (aps_ok) then
997   ! define aps
998   call def_var(nout,"aps","hybrid pressure at midlayers"," ",1,&
999        (/altdimout/),apsid,ierr)
1000   if (ierr.ne.NF_NOERR) then
1001      stop "Error: Failed to def_var aps"
1002   endif
1003
1004   ! write aps
1005#ifdef NC_DOUBLE
1006   ierr=NF_PUT_VAR_DOUBLE(outfid,apsid,aps)
1007#else
1008   ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
1009#endif
1010   if (ierr.ne.NF_NOERR) then
1011      stop "Error: Failed to write aps"
1012   endif
1013ENDIF
1014
1015IF (bps_ok) then
1016   ! define bps
1017   call def_var(nout,"bps","hybrid sigma at midlayers"," ",1,&
1018        (/altdimout/),bpsid,ierr)
1019   if (ierr.ne.NF_NOERR) then
1020      stop "Error: Failed to def_var bps"
1021   endif
1022
1023   ! write bps
1024#ifdef NC_DOUBLE
1025   ierr=NF_PUT_VAR_DOUBLE(outfid,bpsid,bps)
1026#else
1027   ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
1028#endif
1029   if (ierr.ne.NF_NOERR) then
1030      stop "Error: Failed to write bps"
1031   endif
1032ENDIF
1033
1034
1035
1036! Cleanup
1037deallocate(aps)
1038deallocate(bps)
1039
1040end subroutine init2
1041!******************************************************************************
1042subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
1043!==============================================================================
1044! Purpose: Write a variable (i.e: add a variable to a dataset)
1045! called "name"; along with its attributes "title", "units"...
1046! to a file (following the NetCDF format)
1047!==============================================================================
1048! Remarks:
1049! The NetCDF file must be open
1050!==============================================================================
1051
1052implicit none
1053
1054include "netcdf.inc" ! NetCDF definitions
1055
1056!==============================================================================
1057! Arguments:
1058!==============================================================================
1059integer, intent(in) :: nid
1060! nid: [netcdf] file ID #
1061character (len=*), intent(in) :: name
1062! name(): [netcdf] variable's name
1063character (len=*), intent(in) :: title
1064! title(): [netcdf] variable's "title" attribute
1065character (len=*), intent(in) :: units
1066! unit(): [netcdf] variable's "units" attribute
1067integer, intent(in) :: nbdim
1068! nbdim: number of dimensions of the variable
1069integer, dimension(nbdim), intent(in) :: dim
1070! dim(nbdim): [netcdf] dimension(s) ID(s)
1071integer, intent(out) :: nvarid
1072! nvarid: [netcdf] ID # of the variable
1073integer, intent(out) :: ierr
1074! ierr: [netcdf] subroutines returned error code
1075
1076! Switch to netcdf define mode
1077ierr=NF_REDEF(nid)
1078
1079! Insert the definition of the variable
1080#ifdef NC_DOUBLE
1081ierr=NF_DEF_VAR(nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
1082#else
1083ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
1084#endif
1085
1086! Write the attributes
1087ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",len_trim(adjustl(title)),adjustl(title))
1088ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
1089
1090! End netcdf define mode
1091ierr=NF_ENDDEF(nid)
1092
1093end subroutine def_var
1094
1095end program streamfunction
Note: See TracBrowser for help on using the repository browser.