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

Last change on this file since 2513 was 2513, checked in by lrossi, 4 years ago

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