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

Last change on this file since 1847 was 1847, checked in by aslmd, 7 years ago

an idea to cope with different planets in streamfunction
File size: 35.3 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  infile2="diagfi.nc"
399  write(*,*) "         Trying file ",trim(infile2)
400  ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
401  if (ierr.ne.NF_NOERR) then
402     infile2="diagfi1.nc"
403     write(*,*) "         Trying file ",trim(infile2)
404     ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
405     if (ierr.ne.NF_NOERR) then
406        write(*,*) "Problem: Could not find/open these files"
407        stop "Might as well stop here"
408     endif
409  endif
410
411
412   ! Get ID for cv
413   ierr=NF_INQ_VARID(infid2,"cv",tmpvarid)
414   if (ierr.ne.NF_NOERR) then
415      stop "Error: Failed to get cv ID"
416   endif
417   ! Get cv
418   ierr=NF_GET_VAR_REAL(infid2,tmpvarid,cv)
419   if (ierr.ne.NF_NOERR) then
420      stop "Error: Failed reading cv"
421   endif
422   ! Close file
423   write(*,*) 'OK, got cv'
424   ierr=NF_CLOSE(infid2)
425else
426   ierr=NF_GET_VAR_REAL(infid,tmpvarid,cv)
427   if (ierr.ne.NF_NOERR) then
428      stop "Error: Failed reading cv"
429   endif
430endif
431
432
433!Other variables: rho, temp
434
435ierr=NF_INQ_VARID(infid,"temp",tmpvarid)
436if (ierr.ne.NF_NOERR) then
437  stop "Error: Failed to get temp ID"
438else
439  allocate(temp(lonlength,latlength,altlength,timelength))
440  ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp)
441  if (ierr.ne.NF_NOERR) then
442    stop "Error: Failed reading temperature"
443  endif
444endif
445
446ierr=NF_INQ_VARID(infid,"rho",tmpvarid)
447if (ierr.ne.NF_NOERR) then
448  write(*,*) "Error: Failed to get rho ID"
449  write(*,*) "Let's compute it from temperature"
450  allocate(rho(lonlength,latlength,altlength,timelength))
451 
452  do it=1,timelength
453    do ilat=1,latlength
454      do ilon=1,lonlength
455        do ialt=1,altlength
456          rho(ilon,ilat,ialt,it)= (aps(ialt)+bps(ialt)*ps(ilon,ilat,it))/(Rgas*temp(ilon,ilat,ialt,it))
457        enddo
458      enddo
459    enddo
460  enddo
461else
462  allocate(rho(lonlength,latlength,altlength,timelength))
463  ierr=NF_GET_VAR_REAL(infid,tmpvarid,rho)
464  if (ierr.ne.NF_NOERR) then
465    stop "Error: Failed reading mass density"
466  endif
467endif
468
469!==============================================================================
470! 1.4. Output file name
471!==============================================================================
472filename=infile(1:len_trim(infile)-3)//"_stream.nc"
473write(*,*)'The output file has name:'
474write(*,*) filename
475
476
477!==============================================================================
478! 2.1. Temporal averages
479!==============================================================================
480
481allocate(ucum(lonlength,latlength,altlength)) !Time mean zonal wind
482allocate(vcum(lonlength,latlength,altlength)) !Time mean meridional wind
483allocate(tempcum(lonlength,latlength,altlength)) !Time mean temperature
484allocate(rhocum(lonlength,latlength,altlength)) !Time mean density
485allocate(pscum(lonlength,latlength)) !Time mean surface pressure
486do ilat=1,latlength
487   do ilon=1,lonlength
488      pscum(ilon,ilat)=0.
489      do ialt=1,altlength
490            ucum(ilon,ilat,ialt)=0.
491            vcum(ilon,ilat,ialt)=0.
492            tempcum(ilon,ilat,ialt)=0.
493            rhocum(ilon,ilat,ialt)=0.
494            do it=1,timelength
495               ucum(ilon,ilat,ialt)=ucum(ilon,ilat,ialt)+u(ilon,ilat,ialt,it)
496               vcum(ilon,ilat,ialt)=vcum(ilon,ilat,ialt)+v(ilon,ilat,ialt,it)
497               tempcum(ilon,ilat,ialt)=tempcum(ilon,ilat,ialt)+temp(ilon,ilat,ialt,it)
498               rhocum(ilon,ilat,ialt)=rhocum(ilon,ilat,ialt)+rho(ilon,ilat,ialt,it)
499            enddo
500            ucum(ilon,ilat,ialt)=ucum(ilon,ilat,ialt)/timelength
501            vcum(ilon,ilat,ialt)=vcum(ilon,ilat,ialt)/timelength
502            tempcum(ilon,ilat,ialt)=tempcum(ilon,ilat,ialt)/timelength
503            rhocum(ilon,ilat,ialt)=rhocum(ilon,ilat,ialt)/timelength
504      enddo
505      do it=1,timelength
506         pscum(ilon,ilat)=pscum(ilon,ilat)+ps(ilon,ilat,it)
507      enddo
508      pscum(ilon,ilat)=pscum(ilon,ilat)/timelength
509   enddo
510enddo
511
512!==============================================================================
513! 2.2. Calculation of the stream function
514!==============================================================================
515
516!Contravariant meridional wind
517allocate(vcont(lonlength,latlength,altlength))
518do ilon=1,lonlength
519   do ialt=1,altlength
520      do ilat=1,latlength-1
521         vcont(ilon,ilat,ialt)=0.5*&
522              (vcum(ilon,ilat,ialt)+vcum(ilon,ilat+1,ialt))/cv(ilon,ilat)
523      enddo
524      vcont(ilon,latlength,ialt)=0.!vcont(ilon,latlength-1,ialt)
525   enddo
526enddo
527
528!Pressure from hybrid levels
529allocate(Plev(lonlength,latlength,altlength))
530do ilon=1,lonlength
531   do ilat=1,latlength
532      do ialt=1,altlength
533         Plev(ilon,ilat,ialt)=aps(ialt)+bps(ialt)*pscum(ilon,ilat)
534      enddo
535   enddo
536enddo
537
538!Mass in the box
539allocate(mass(lonlength,latlength,altlength))
540do ilon=1,lonlength
541   do ilat=1,latlength
542      do ialt=1,altlength
543         mass(ilon,ilat,ialt)=aire(ilon,ilat)*&
544              Plev(ilon,ilat,ialt)/g0
545      enddo
546   enddo
547enddo
548
549!Mass variation in the box
550allocate(dmass(lonlength,latlength,altlength))
551do ilon=1,lonlength
552   do ilat=1,latlength
553      do ialt=1,altlength-1
554         dmass(ilon,ilat,ialt)=mass(ilon,ilat,ialt)-&
555              mass(ilon,ilat,ialt+1)
556      enddo
557      dmass(ilon,ilat,altlength)=0.!dmass(ilon,ilat,altlength-1)
558   enddo
559enddo
560
561!Stream function
562allocate(psi(latlength,altlength))
563allocate(vcontcum(latlength,altlength))
564do ilat=1,latlength-1
565   psi(ilat,altlength)=0.
566   do ialt=altlength-1,1,-1
567      psi(ilat,ialt)=psi(ilat,ialt+1)
568      vcontcum(ilat,ialt)=0.
569      do ilon=1,lonlength
570         vcontcum(ilat,ialt)=vcontcum(ilat,ialt)+vcont(ilon,ilat,ialt)/lonlength
571         psi(ilat,ialt)=psi(ilat,ialt)-(vcont(ilon,ilat,ialt)&
572              *(dmass(ilon,ilat,ialt)+dmass(ilon,ilat+1,ialt)))
573      enddo
574      psi(latlength,ialt)=0.!psi(latlength-1,ialt)
575   enddo
576enddo
577
578!==============================================================================
579! 2.3. Calculation of the total angular momentum
580!==============================================================================
581
582!Rotation speed
583omega=2.*3.141592/daysec
584
585!Angular momentum
586allocate(mom(lonlength,latlength,altlength))
587allocate(momave(latlength,altlength))
588do ilat=1,latlength
589   do ialt=1,altlength
590      momave(ilat,ialt)=0.
591      do ilon=1,lonlength
592!            mom(ilon,ilat,ialt,it)=dmass(ilon,ilat,ialt,it)*a0*&
593         mom(ilon,ilat,ialt)=a0*cos(lat(ilat)*3.141592/180.)*&
594              (omega*a0*cos(lat(ilat)*3.141592/180.)+ucum(ilon,ilat,ialt))
595         
596         momave(ilat,ialt)=momave(ilat,ialt)+mom(ilon,ilat,ialt)
597      enddo
598      momave(ilat,ialt)=momave(ilat,ialt)/lonlength
599   enddo
600enddo
601
602!==============================================================================
603! 2.4. Zonal mean winds, temperature and density
604!==============================================================================
605
606allocate(uzm(latlength,altlength)) !Zonal mean zonal wind
607allocate(vzm(latlength,altlength)) !Zonal mean meridional wind
608allocate(tempzm(latlength,altlength)) !Zonal mean temperature
609allocate(rhozm(latlength,altlength)) !Zonal mean density
610allocate(pszm(latlength)) !Zonal mean surface pressure
611allocate(phisinitzm(latlength)) !Zonal mean ground geopotential
612do ilat=1,latlength
613   phisinitzm(ilat)=0.
614   pszm(ilat)=0.
615   do ialt=1,altlength
616      uzm(ilat,ialt)=0.
617      vzm(ilat,ialt)=0.
618      tempzm(ilat,ialt)=0.
619      rhozm(ilat,ialt)=0.
620      do ilon=1,lonlength
621         uzm(ilat,ialt)=uzm(ilat,ialt)+ucum(ilon,ilat,ialt)
622         vzm(ilat,ialt)=vzm(ilat,ialt)+vcum(ilon,ilat,ialt)
623         tempzm(ilat,ialt)=tempzm(ilat,ialt)+tempcum(ilon,ilat,ialt)
624         rhozm(ilat,ialt)=rhozm(ilat,ialt)+rhocum(ilon,ilat,ialt)
625      enddo
626      uzm(ilat,ialt)=uzm(ilat,ialt)/lonlength
627      vzm(ilat,ialt)=vzm(ilat,ialt)/lonlength
628      tempzm(ilat,ialt)=tempzm(ilat,ialt)/lonlength
629      rhozm(ilat,ialt)=rhozm(ilat,ialt)/lonlength
630   enddo
631   do ilon=1,lonlength
632      pszm(ilat)=pszm(ilat)+pscum(ilon,ilat)
633      phisinitzm(ilat)=phisinitzm(ilat)+phisinit(ilon,ilat)
634   enddo
635   pszm(ilat)=pszm(ilat)/lonlength
636   phisinitzm(ilat)=phisinitzm(ilat)/lonlength 
637enddo
638
639
640!==============================================================================
641! 2.3. Recalculation of angular momentum using zonally averaged wind
642!==============================================================================
643
644do ilat=1,latlength
645   do ialt=1,altlength
646      momave(ilat,ialt)=a0*cos(lat(ilat)*3.141592/180.)*&
647           (omega*a0*cos(lat(ilat)*3.141592/180.)+uzm(ilat,ialt))
648   enddo
649enddo
650
651
652!==============================================================================
653! 3.1 Dimensions in output file
654!==============================================================================
655
656! 3.1.1 Initialize output file's lat,lon,alt and time dimensions
657call initiate (filename,lat,alt,time,nout,&
658     latdimout,londimout,altdimout,timedimout,lonvarout,timevarout)
659! Initialize output file's aps,bps variables
660call init2(infid,lonlength,latlength,altlength,&
661     nout,londimout,latdimout,altdimout)
662
663! 3.1.2 New longitude dimension/value in output file
664allocate(lon_fake(1))
665lon_fake(1)=0.
666#ifdef NC_DOUBLE
667        ierr= NF_PUT_VARA_DOUBLE(nout,lonvarout,1,1,lon_fake)
668#else
669        ierr= NF_PUT_VARA_REAL(nout,lonvarout,1,1,lon_fake)
670#endif
671
672! 3.1.3 New time dimension/value in output file
673#ifdef NC_DOUBLE
674        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,1,1,lon_fake)
675#else
676        ierr= NF_PUT_VARA_REAL(nout,timevarout,1,1,lon_fake)
677#endif
678
679
680!==============================================================================
681! 3.2 Write variables
682!==============================================================================
683
684!Define the dimensions of the variables to be written
685
686!!3.2.1 Stream function
687call def_var(nout,"psi","Stream function","",4,&
688             (/londimout,latdimout,altdimout,timedimout/),psivarout,ierr)
689if (ierr.ne.NF_NOERR) then
690   write(*,*) 'Error, could not define variable psi'
691   stop ""
692endif
693
694!Write in the output file
695ierr=NF_PUT_VAR_REAL(nout,psivarout,psi)
696if (ierr.ne.NF_NOERR) then
697   write(*,*)'Error, Failed to write variable psi'
698   stop
699endif
700
701!3.2.2 Momentum
702call def_var(nout,"momave","Angular momentum","",4,&
703             (/londimout,latdimout,altdimout,timedimout/),momvarout,ierr)
704if (ierr.ne.NF_NOERR) then
705   write(*,*) 'Error, could not define variable momave'
706   stop ""
707endif
708
709!Write in the output file
710ierr=NF_PUT_VAR_REAL(nout,momvarout,momave)
711if (ierr.ne.NF_NOERR) then
712   write(*,*)'Error, Failed to write variable momave'
713   stop
714endif
715
716
717!3.2.3 Zonal mean zonal wind
718call def_var(nout,"u","Zonal mean zonal wind","m/s",4,&
719             (/londimout,latdimout,altdimout,timedimout/),uvarout,ierr)
720if (ierr.ne.NF_NOERR) then
721   write(*,*) 'Error, could not define variable u'
722   stop ""
723endif
724
725!Write in the output file
726ierr=NF_PUT_VAR_REAL(nout,uvarout,uzm)
727if (ierr.ne.NF_NOERR) then
728   write(*,*)'Error, Failed to write variable u'
729   stop
730endif
731
732
733!3.2.4 Zonal mean meridional wind
734call def_var(nout,"v","Zonal mean meridional wind","m/s",4,&
735             (/londimout,latdimout,altdimout,timedimout/),vvarout,ierr)
736if (ierr.ne.NF_NOERR) then
737   write(*,*) 'Error, could not define variable v'
738   stop ""
739endif
740
741!Write in the output file
742ierr=NF_PUT_VAR_REAL(nout,vvarout,vzm)
743if (ierr.ne.NF_NOERR) then
744   write(*,*)'Error, Failed to write variable v'
745   stop
746endif
747
748!3.2.5 Zonal mean density
749call def_var(nout,"rho","Zonal mean atmospheric density","",4,&
750             (/londimout,latdimout,altdimout,timedimout/),rhovarout,ierr)
751if (ierr.ne.NF_NOERR) then
752   write(*,*) 'Error, could not define variable rho'
753   stop ""
754endif
755
756!Write in the output file
757ierr=NF_PUT_VAR_REAL(nout,rhovarout,rhozm)
758if (ierr.ne.NF_NOERR) then
759   write(*,*)'Error, Failed to write variable rho'
760   stop
761endif
762
763!3.2.6 Zonal mean temperature
764call def_var(nout,"temp","Zonal mean temperature","K",4,&
765             (/londimout,latdimout,altdimout,timedimout/),tempvarout,ierr)
766if (ierr.ne.NF_NOERR) then
767   write(*,*) 'Error, could not define variable temp'
768   stop ""
769endif
770
771!Write in the output file
772ierr=NF_PUT_VAR_REAL(nout,tempvarout,tempzm)
773if (ierr.ne.NF_NOERR) then
774   write(*,*)'Error, Failed to write variable temp'
775   stop
776endif
777
778!3.2.7 Zonal mean surface pressure
779call def_var(nout,"ps","Zonal mean surface pressure","Pa",3,&
780             (/londimout,latdimout,timedimout/),psvarout,ierr)
781if (ierr.ne.NF_NOERR) then
782   write(*,*) 'Error, could not define variable ps'
783   stop ""
784endif
785
786!Write in the output file
787ierr=NF_PUT_VAR_REAL(nout,psvarout,pszm)
788if (ierr.ne.NF_NOERR) then
789   write(*,*)'Error, Failed to write variable ps'
790   stop
791endif
792
793
794!3.2.8 Zonal mean geopotential
795call def_var(nout,"phisinit","Zonal mean initial geopotential","",2,&
796             (/londimout,latdimout/),phisinitvarout,ierr)
797if (ierr.ne.NF_NOERR) then
798   write(*,*) 'Error, could not define variable phisinit'
799   stop ""
800endif
801
802!Write in the output file
803ierr=NF_PUT_VAR_REAL(nout,phisinitvarout,phisinitzm)
804if (ierr.ne.NF_NOERR) then
805   write(*,*)'Error, Failed to write variable phisinit'
806   stop
807endif
808
809! Close input file
810ierr=nf_close(nid)
811
812! Close output file
813ierr=NF_CLOSE(nout)
814
815contains
816
817!******************************************************************************
818Subroutine initiate (filename,lat,alt,time,&
819                     nout,latdimout,londimout,altdimout,timedimout,&
820                     lonvarout,timevarout)
821!==============================================================================
822! Purpose:
823! Create and initialize a data file (NetCDF format)
824!==============================================================================
825! Remarks:
826! The NetCDF file (created in this subroutine) remains open
827!==============================================================================
828
829implicit none
830
831include "netcdf.inc" ! NetCDF definitions
832
833!==============================================================================
834! Arguments:
835!==============================================================================
836character (len=*), intent(in):: filename
837! filename(): the file's name
838real, dimension(:), intent(in):: lat
839! lat(): latitude
840real, dimension(:), intent(in):: alt
841! alt(): altitude
842real, dimension(:), intent(in):: time
843! time(): Time
844integer, intent(out):: nout
845! nout: [netcdf] file ID
846integer, intent(out):: latdimout
847! latdimout: [netcdf] lat() (i.e.: latitude)  ID
848integer, intent(out):: londimout
849! londimout: [netcdf] lon()  ID
850integer, intent(out):: altdimout
851! altdimout: [netcdf] alt()  ID
852integer, intent(out):: timedimout
853! timedimout: [netcdf] "Time"  ID
854integer, intent(out):: lonvarout
855! timevarout: [netcdf] Longiture (considered as a variable) ID
856integer, intent(out):: timevarout
857! timevarout: [netcdf] Time (considered as a variable) ID
858
859!==============================================================================
860! Local variables:
861!==============================================================================
862!integer :: latdim,londim,altdim,timedim
863integer :: nvarid,ierr
864! nvarid: [netcdf] ID of a variable
865! ierr: [netcdf]  return error code (from called subroutines)
866
867!==============================================================================
868! 1. Create (and open) output file
869!==============================================================================
870write(*,*) "creating "//trim(adjustl(filename))//'...'
871ierr = NF_CREATE(filename,NF_CLOBBER,nout)
872! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
873if (ierr.NE.NF_NOERR) then
874   WRITE(*,*)'ERROR: Impossible to create the file.'
875   stop ""
876endif
877
878!==============================================================================
879! 2. Define/write "dimensions" and get their IDs
880!==============================================================================
881
882ierr = NF_DEF_DIM(nout, "latitude", size(lat), latdimout)
883!ierr = NF_DEF_DIM(nout, "longitude", NF_UNLIMITED, londimout)
884ierr = NF_DEF_DIM(nout, "longitude", 1, londimout)
885ierr = NF_DEF_DIM(nout, "altitude", size(alt), altdimout)
886ierr = NF_DEF_DIM(nout, "Time", 1, timedimout)
887
888! End netcdf define mode
889ierr = NF_ENDDEF(nout)
890
891!==============================================================================
892! 3. Write "Time" (attributes)
893!==============================================================================
894
895call def_var(nout,"Time","Time","years since 0000-00-0 00:00:00",1,&
896             (/timedimout/),timevarout,ierr)
897
898!==============================================================================
899! 4. Write "latitude" (data and attributes)
900!==============================================================================
901
902call def_var(nout,"latitude","latitude","degrees_north",1,&
903             (/latdimout/),nvarid,ierr)
904
905#ifdef NC_DOUBLE
906ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lat)
907#else
908ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
909#endif
910
911!==============================================================================
912! 4. Write "longitude" (attributes)
913!==============================================================================
914
915call def_var(nout,"longitude","East longitude","degrees_east",1,&
916             (/londimout/),lonvarout,ierr)
917
918
919!==============================================================================
920! 4. Write "altitude" (data and attributes)
921!==============================================================================
922
923! Switch to netcdf define mode
924
925call def_var(nout,"altitude","Altitude","km",1,&
926             (/altdimout/),nvarid,ierr)
927
928#ifdef NC_DOUBLE
929ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,alt)
930#else
931ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
932#endif
933
934end Subroutine initiate
935!******************************************************************************
936subroutine init2(infid,lonlen,latlen,altlen, &
937                 outfid,londimout,latdimout,altdimout)
938!==============================================================================
939! Purpose:
940! Copy aps() and bps() from input file to outpout file
941!==============================================================================
942! Remarks:
943! The NetCDF files must be open
944!==============================================================================
945
946implicit none
947
948include "netcdf.inc" ! NetCDF definitions
949
950!==============================================================================
951! Arguments:
952!==============================================================================
953integer, intent(in) :: infid  ! NetCDF output file ID
954integer, intent(in) :: lonlen ! # of grid points along longitude
955integer, intent(in) :: latlen ! # of grid points along latitude
956integer, intent(in) :: altlen ! # of grid points along latitude
957integer, intent(in) :: outfid ! NetCDF output file ID
958integer, intent(in) :: londimout ! longitude dimension ID
959integer, intent(in) :: latdimout ! latitude dimension ID
960integer, intent(in) :: altdimout ! altitude dimension ID
961!==============================================================================
962! Local variables:
963!==============================================================================
964real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
965integer :: apsid,bpsid
966integer :: ierr
967integer :: tmpvarid ! temporary variable ID
968logical :: aps_ok, bps_ok ! are "phisinit" "aps" "bps" available ?
969
970
971!==============================================================================
972! 1. Read data from input file
973!==============================================================================
974
975! hybrid coordinate aps
976  allocate(aps(altlen))
977ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
978if (ierr.ne.NF_NOERR) then
979  write(*,*) "oops Failed to get aps ID. OK"
980  aps_ok=.false.
981else
982  ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
983  if (ierr.ne.NF_NOERR) then
984   stop "error: Failed reading aps"
985  endif
986  aps_ok=.true.
987endif
988
989! hybrid coordinate bps
990  allocate(bps(altlen))
991ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
992if (ierr.ne.NF_NOERR) then
993  write(*,*) "oops: Failed to get bps ID. OK"
994  bps_ok=.false.
995else
996  ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
997  if (ierr.ne.NF_NOERR) then
998    stop "Error: Failed reading bps"
999  endif
1000  bps_ok=.true.
1001endif
1002
1003!==============================================================================
1004! 2. Write
1005!==============================================================================
1006
1007!==============================================================================
1008! 2.2. Hybrid coordinates aps() and bps()
1009!==============================================================================
1010
1011IF (aps_ok) then
1012   ! define aps
1013   call def_var(nout,"aps","hybrid pressure at midlayers"," ",1,&
1014        (/altdimout/),apsid,ierr)
1015   if (ierr.ne.NF_NOERR) then
1016      stop "Error: Failed to def_var aps"
1017   endif
1018
1019   ! write aps
1020#ifdef NC_DOUBLE
1021   ierr=NF_PUT_VAR_DOUBLE(outfid,apsid,aps)
1022#else
1023   ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
1024#endif
1025   if (ierr.ne.NF_NOERR) then
1026      stop "Error: Failed to write aps"
1027   endif
1028ENDIF
1029
1030IF (bps_ok) then
1031   ! define bps
1032   call def_var(nout,"bps","hybrid sigma at midlayers"," ",1,&
1033        (/altdimout/),bpsid,ierr)
1034   if (ierr.ne.NF_NOERR) then
1035      stop "Error: Failed to def_var bps"
1036   endif
1037
1038   ! write bps
1039#ifdef NC_DOUBLE
1040   ierr=NF_PUT_VAR_DOUBLE(outfid,bpsid,bps)
1041#else
1042   ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
1043#endif
1044   if (ierr.ne.NF_NOERR) then
1045      stop "Error: Failed to write bps"
1046   endif
1047ENDIF
1048
1049
1050
1051! Cleanup
1052deallocate(aps)
1053deallocate(bps)
1054
1055end subroutine init2
1056!******************************************************************************
1057subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
1058!==============================================================================
1059! Purpose: Write a variable (i.e: add a variable to a dataset)
1060! called "name"; along with its attributes "title", "units"...
1061! to a file (following the NetCDF format)
1062!==============================================================================
1063! Remarks:
1064! The NetCDF file must be open
1065!==============================================================================
1066
1067implicit none
1068
1069include "netcdf.inc" ! NetCDF definitions
1070
1071!==============================================================================
1072! Arguments:
1073!==============================================================================
1074integer, intent(in) :: nid
1075! nid: [netcdf] file ID #
1076character (len=*), intent(in) :: name
1077! name(): [netcdf] variable's name
1078character (len=*), intent(in) :: title
1079! title(): [netcdf] variable's "title" attribute
1080character (len=*), intent(in) :: units
1081! unit(): [netcdf] variable's "units" attribute
1082integer, intent(in) :: nbdim
1083! nbdim: number of dimensions of the variable
1084integer, dimension(nbdim), intent(in) :: dim
1085! dim(nbdim): [netcdf] dimension(s) ID(s)
1086integer, intent(out) :: nvarid
1087! nvarid: [netcdf] ID # of the variable
1088integer, intent(out) :: ierr
1089! ierr: [netcdf] subroutines returned error code
1090
1091! Switch to netcdf define mode
1092ierr=NF_REDEF(nid)
1093
1094! Insert the definition of the variable
1095#ifdef NC_DOUBLE
1096ierr=NF_DEF_VAR(nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
1097#else
1098ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
1099#endif
1100
1101! Write the attributes
1102ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",len_trim(adjustl(title)),adjustl(title))
1103ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
1104
1105! End netcdf define mode
1106ierr=NF_ENDDEF(nid)
1107
1108end subroutine def_var
1109
1110end program streamfunction
Note: See TracBrowser for help on using the repository browser.