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

Last change on this file was 2567, checked in by emillour, 3 years ago

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