source: trunk/LMDZ.TITAN/Tools/tmc.F90 @ 1379

Last change on this file since 1379 was 816, checked in by slebonnois, 12 years ago

SL: tools for postprocessing (Veznus and Titan); see DOC/documentation/vt-tools.pdf
File size: 24.0 KB
Line 
1program tmc
2
3! SL 01/2010:
4! This program reads 4D (lon-lat-alt-time) fields recast in log P coordinates
5!
6! it computes angular momentum transport from high-frequency outputs:
7!
8! totvang -- 2D -- Meridional transport of angular momentum, total (m3 s-2)
9! totwang -- 2D --   Vertical transport of angular momentum, total (m3 s-2)
10! mmcvang -- 2D -- Meridional transport of angular momentum, by MMC (m3 s-2)
11! mmcwang -- 2D --   Vertical transport of angular momentum, by MMC (m3 s-2)
12! trsvang -- 2D -- Meridional transport of angular momentum, transients (m3 s-2)
13! trswang -- 2D --   Vertical transport of angular momentum, transients (m3 s-2)
14! stnvang -- 2D -- Meridional transport of angular momentum, stationaries (m3 s-2)
15! stnwang -- 2D --   Vertical transport of angular momentum, stationaries (m3 s-2)
16! dmass   -- 2D -- Mass in each cell (dmassmeanbar)
17!
18! Minimal requirements and dependencies:
19! The dataset must include the following data:
20! - pressure vertical coordinate
21! - surface pressure
22! - zonal, meridional and vertical (Pa/s) winds
23! - altitude above areoid
24!
25! Convention: qbar  <=> zonal average    / qstar = q - qbar
26!             qmean <=> temporal average / qprim = q - qmean
27!
28!  Therefore: ((qv)mean)bar                         (total)
29!                          =  qmeanbar *  vmeanbar  (mmc)
30!                      + (qstarmean * vstarmean)bar (stn)
31!                         + (qprim * vprim)meanbar  (trs)
32!
33!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
34
35implicit none
36
37include "netcdf.inc" ! NetCDF definitions
38
39character (len=128) :: infile ! input file name (name_P.nc)
40character (len=128) :: outfile ! output file name
41
42character (len=64) :: text ! to store some text
43integer infid ! NetCDF input file ID
44integer outfid ! NetCDF output file ID
45integer lon_dimid,lat_dimid,alt_dimid,time_dimid ! NetCDF dimension IDs
46integer lon_varid,lat_varid,alt_varid,time_varid
47integer,dimension(2) :: datashape2d ! shape of 3D datasets
48integer,dimension(3) :: datashape3d ! shape of 3D datasets
49integer,dimension(4) :: datashape4d ! shape of 4D datasets
50
51real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
52real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
53real :: pi
54real,dimension(:),allocatable :: lon ! longitude
55integer lonlength ! # of grid points along longitude
56real,dimension(:),allocatable :: lat ! latitude
57real,dimension(:),allocatable :: latrad ! latitude in radian
58integer latlength ! # of grid points along latitude
59real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
60integer altlength ! # of grid point along altitude (of input datasets)
61real,dimension(:),allocatable :: time ! time
62integer timelength ! # of points along time
63real,dimension(:,:,:),allocatable :: ps ! surface pressure
64real,dimension(:,:,:,:),allocatable :: vitu ! zonal wind (in m/s)
65real,dimension(:,:,:,:),allocatable :: vitv ! meridional wind (in m/s)
66real,dimension(:,:,:,:),allocatable :: vitw ! vertical wind (in Pa/s, then converted in m/s)
67real,dimension(:,:,:,:),allocatable :: za ! above areoid levels (m)
68
69!!! output variables
70real,dimension(:,:),allocatable :: totvang ! merid transport of ang momentum
71real,dimension(:,:),allocatable :: totwang ! verti transport of ang momentum
72real,dimension(:,:),allocatable :: mmcvang ! merid transport of ang momentum
73real,dimension(:,:),allocatable :: mmcwang ! verti transport of ang momentum
74real,dimension(:,:),allocatable :: trsvang ! merid transport of ang momentum
75real,dimension(:,:),allocatable :: trswang ! verti transport of ang momentum
76real,dimension(:,:),allocatable :: stnvang ! merid transport of ang momentum
77real,dimension(:,:),allocatable :: stnwang ! verti transport of ang momentum
78real,dimension(:,:),allocatable :: dmassmeanbar
79
80! local variables
81real :: deltalat,deltalon ! lat and lon intervals in radians
82real,dimension(:,:,:,:),allocatable :: rayon ! distance to center (m)
83real,dimension(:,:,:,:),allocatable :: grav ! gravity field (m s-2)
84real,dimension(:,:,:,:),allocatable :: dmass ! mass in cell (kg)
85real,dimension(:,:,:),allocatable   :: dmassmean
86
87real,dimension(:,:,:,:),allocatable :: osam ! planetary rotation specific ang (m2/s)
88real,dimension(:,:,:,:),allocatable :: rsam ! zonal wind specific ang (m2/s)
89real,dimension(:,:,:,:),allocatable :: tsam ! total specific ang (m2/s)
90
91real,dimension(:,:,:,:),allocatable :: vang ! v * specific ang (m3/s)
92real,dimension(:,:,:,:),allocatable :: wang ! w * specific ang (m3/s)
93real,dimension(:,:,:,:),allocatable :: vstar
94real,dimension(:,:,:,:),allocatable :: wstar
95real,dimension(:,:,:,:),allocatable :: angstar
96real,dimension(:,:,:,:),allocatable :: vprim
97real,dimension(:,:,:,:),allocatable :: wprim
98real,dimension(:,:,:,:),allocatable :: angprim
99real,dimension(:,:,:,:),allocatable :: angprimvprim
100real,dimension(:,:,:,:),allocatable :: angprimwprim
101
102! lon,lat,alt
103real,dimension(:,:,:),allocatable :: vangmean
104real,dimension(:,:,:),allocatable :: wangmean
105real,dimension(:,:,:),allocatable :: vmean
106real,dimension(:,:,:),allocatable :: wmean
107real,dimension(:,:,:),allocatable :: angmean
108real,dimension(:,:,:),allocatable :: angprimvprimmean
109real,dimension(:,:,:),allocatable :: angprimwprimmean
110real,dimension(:,:,:),allocatable :: vstarmean
111real,dimension(:,:,:),allocatable :: wstarmean
112real,dimension(:,:,:),allocatable :: angstarmean
113real,dimension(:,:,:),allocatable :: angstarmeanvstarmean
114real,dimension(:,:,:),allocatable :: angstarmeanwstarmean
115real,dimension(:,:,:),allocatable :: v3d    ! intermediate for vbar
116real,dimension(:,:,:),allocatable :: w3d    ! intermediate for wbar
117real,dimension(:,:,:),allocatable :: ang3d  ! intermediate for angbar
118
119! lat,alt,time
120real,dimension(:,:,:),allocatable :: vbar
121real,dimension(:,:,:),allocatable :: wbar
122real,dimension(:,:,:),allocatable :: angbar
123
124!lat,alt
125real,dimension(:,:),allocatable :: vmeanbar
126real,dimension(:,:),allocatable :: wmeanbar
127real,dimension(:,:),allocatable :: angmeanbar
128real,dimension(:,:),allocatable :: vbar2d     ! intermediate for vbar
129real,dimension(:,:),allocatable :: wbar2d     ! intermediate for wbar
130real,dimension(:,:),allocatable :: angbar2d   ! intermediate for qbar
131
132integer ierr,ierr1,ierr2 ! NetCDF routines return codes
133integer i,j,ilon,ilat,ilev,itim ! for loops
134logical :: lmdflag
135
136include "planet.h"
137
138!===============================================================================
139! 1. Input parameters
140!===============================================================================
141
142pi = 2.*asin(1.)
143
144write(*,*) ""
145write(*,*) "You are working on the atmosphere of ",planet
146
147!===============================================================================
148! 1.1 Input file
149!===============================================================================
150
151write(*,*) ""
152write(*,*) "Program valid for files with pressure axis (*_P.nc)"
153write(*,*) "Enter input file name:"
154
155read(*,'(a128)') infile
156write(*,*) ""
157
158! open input file
159
160ierr = NF_OPEN(infile,NF_NOWRITE,infid)
161if (ierr.ne.NF_NOERR) then
162   write(*,*) 'ERROR: Pb opening file ',trim(infile)
163   stop ""
164endif
165
166!===============================================================================
167! 1.2 Get grids in lon,lat,alt(pressure),time
168!===============================================================================
169
170call get_iddim(infid,lat_varid,latlength,lon_varid,lonlength,&
171                     alt_varid,altlength,time_varid,timelength,lmdflag )
172
173allocate(lon(lonlength))
174ierr=NF_GET_VAR_REAL(infid,lon_varid,lon)
175if (ierr.ne.NF_NOERR) stop "Error: Failed reading longitude"
176
177allocate(lat(latlength))
178ierr=NF_GET_VAR_REAL(infid,lat_varid,lat)
179if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
180
181allocate(latrad(latlength))
182latrad = lat*pi/180.
183
184! Lat, lon pressure intervals
185deltalat = abs(latrad(2)-latrad(1))
186deltalon = abs(lon(2)-lon(1))*pi/180.
187
188allocate(plev(altlength))
189ierr=NF_GET_VAR_REAL(infid,alt_varid,plev)
190if (ierr.ne.NF_NOERR) stop "Error: Failed reading altitude (ie pressure levels)"
191
192allocate(time(timelength))
193ierr=NF_GET_VAR_REAL(infid,time_varid,time)
194if (ierr.ne.NF_NOERR) stop "Error: Failed reading time"
195
196!===============================================================================
197! 1.3 Get output file name
198!===============================================================================
199write(*,*) ""
200!write(*,*) "Enter output file name"
201!read(*,*) outfile
202outfile=infile(1:len_trim(infile)-3)//"_TMC.nc"
203write(*,*) "Output file name is: "//trim(outfile)
204
205
206
207!===============================================================================
208! 2.1 Store needed fields
209!===============================================================================
210
211!===============================================================================
212! 2.1.1 Surface pressure
213!===============================================================================
214allocate(ps(lonlength,latlength,timelength))
215
216text="ps"
217call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
218if (ierr1.ne.NF_NOERR) then
219  write(*,*) "  looking for psol instead... "
220  text="psol"
221  call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
222  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get psol ID"
223endif
224if (ierr2.ne.NF_NOERR) stop "Error: Failed reading surface pressure"
225
226!===============================================================================
227! 2.1.3 Winds
228!===============================================================================
229allocate(vitu(lonlength,latlength,altlength,timelength))
230allocate(vitv(lonlength,latlength,altlength,timelength))
231allocate(vitw(lonlength,latlength,altlength,timelength))
232
233! zonal wind vitu (in m/s)
234text="vitu"
235call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
236if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
237if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
238
239! meridional wind vitv (in m/s)
240text="vitv"
241call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
242if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
243if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
244
245! vertical wind vitw (in Pa/s)
246text="vitw"
247call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
248if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitw ID"
249if (ierr2.ne.NF_NOERR) stop "Error: Failed reading vertical wind"
250
251!===============================================================================
252! 2.1.4 Altitude above areoide
253!===============================================================================
254! Only needed if g(z) on Titan...
255
256! allocate(za(lonlength,latlength,altlength,timelength))
257
258! text="zareoid"
259! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
260! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
261! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
262
263!===============================================================================
264! 2.2 Computations
265!===============================================================================
266
267!===============================================================================
268! 2.2.2 Mass in cells
269!===============================================================================
270allocate(rayon(lonlength,latlength,altlength,timelength))
271allocate( grav(lonlength,latlength,altlength,timelength))
272allocate(dmass(lonlength,latlength,altlength,timelength))
273allocate(dmassmean(lonlength,latlength,altlength))
274allocate(dmassmeanbar(latlength,altlength))
275
276do itim=1,timelength
277do ilon=1,lonlength
278 do ilat=1,latlength
279  do ilev=1,altlength
280! Need to be consistent with GCM computations
281!    if (za(ilon,ilat,ilev,itim).ne.miss_val) then
282     rayon(ilon,ilat,ilev,itim) = a0
283!     rayon(ilon,ilat,ilev,itim) = za(ilon,ilat,ilev,itim) + a0
284      grav(ilon,ilat,ilev,itim) = g0*a0*a0 &
285                 /(rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim))
286!    else
287!     rayon(ilon,ilat,ilev,itim) = miss_val
288!      grav(ilon,ilat,ilev,itim) = miss_val
289!    endif
290  enddo
291 enddo
292enddo
293enddo ! timelength
294
295call cellmass(infid,latlength,lonlength,altlength,timelength, &
296              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
297              dmass )
298
299call moytim(lonlength,latlength,altlength,timelength,miss_val,dmass,dmassmean)
300call moyzon2(lonlength,latlength,altlength,miss_val,dmassmean,dmassmeanbar)
301
302print*,"OK dmass"
303
304!===============================================================================
305! 2.2.3 Specific angular momentum
306!===============================================================================
307allocate(osam(lonlength,latlength,altlength,timelength))
308allocate(rsam(lonlength,latlength,altlength,timelength))
309allocate(tsam(lonlength,latlength,altlength,timelength))
310
311do itim=1,timelength
312do ilon=1,lonlength
313 do ilat=1,latlength
314  do ilev=1,altlength
315    if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
316      osam(ilon,ilat,ilev,itim) = &
317         rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim) &
318       * cos(latrad(ilat))*cos(latrad(ilat)) &
319       * omega
320      rsam(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim) &
321       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))
322      tsam(ilon,ilat,ilev,itim) = osam(ilon,ilat,ilev,itim)&
323                                + rsam(ilon,ilat,ilev,itim)
324    else
325      osam(ilon,ilat,ilev,itim) = miss_val
326      rsam(ilon,ilat,ilev,itim) = miss_val
327      tsam(ilon,ilat,ilev,itim) = miss_val
328    endif
329  enddo
330 enddo
331enddo
332enddo ! timelength
333
334print*,"debut tprt ang"
335
336!===============================================================================
337! 2.2.4 Angular momentum transport
338!===============================================================================
339! allocations
340!-------------
341allocate(vang(lonlength,latlength,altlength,timelength))
342allocate(wang(lonlength,latlength,altlength,timelength))
343allocate(  vstar(lonlength,latlength,altlength,timelength))
344allocate(  wstar(lonlength,latlength,altlength,timelength))
345allocate(angstar(lonlength,latlength,altlength,timelength))
346allocate(  vprim(lonlength,latlength,altlength,timelength))
347allocate(  wprim(lonlength,latlength,altlength,timelength))
348allocate(angprim(lonlength,latlength,altlength,timelength))
349allocate(angprimvprim(lonlength,latlength,altlength,timelength))
350allocate(angprimwprim(lonlength,latlength,altlength,timelength))
351
352! lon,lat,alt
353allocate(vangmean(lonlength,latlength,altlength))
354allocate(wangmean(lonlength,latlength,altlength))
355allocate(  vmean(lonlength,latlength,altlength))
356allocate(  wmean(lonlength,latlength,altlength))
357allocate(angmean(lonlength,latlength,altlength))
358allocate(angprimvprimmean(lonlength,latlength,altlength))
359allocate(angprimwprimmean(lonlength,latlength,altlength))
360allocate(  vstarmean(lonlength,latlength,altlength))
361allocate(  wstarmean(lonlength,latlength,altlength))
362allocate(angstarmean(lonlength,latlength,altlength))
363allocate(angstarmeanvstarmean(lonlength,latlength,altlength))
364allocate(angstarmeanwstarmean(lonlength,latlength,altlength))
365allocate(  v3d(lonlength,latlength,altlength))
366allocate(  w3d(lonlength,latlength,altlength))
367allocate(ang3d(lonlength,latlength,altlength))
368
369! lat,alt,time
370allocate(  vbar(latlength,altlength,timelength))
371allocate(  wbar(latlength,altlength,timelength))
372allocate(angbar(latlength,altlength,timelength))
373
374!lat,alt
375allocate(  vmeanbar(latlength,altlength))
376allocate(  wmeanbar(latlength,altlength))
377allocate(angmeanbar(latlength,altlength))
378allocate(  vbar2d(latlength,altlength))
379allocate(  wbar2d(latlength,altlength))
380allocate(angbar2d(latlength,altlength))
381
382allocate(totvang(latlength,altlength))
383allocate(totwang(latlength,altlength))
384allocate(mmcvang(latlength,altlength))
385allocate(mmcwang(latlength,altlength))
386allocate(trsvang(latlength,altlength))
387allocate(trswang(latlength,altlength))
388allocate(stnvang(latlength,altlength))
389allocate(stnwang(latlength,altlength))
390
391! intermediates
392!-----------------
393
394do itim=1,timelength
395   v3d(:,:,:) = vitv(:,:,:,itim)
396   w3d(:,:,:) = vitw(:,:,:,itim)
397 ang3d(:,:,:) = tsam(:,:,:,itim)
398 call moyzon2(lonlength,latlength,altlength,miss_val,  v3d,  vbar2d)
399 call moyzon2(lonlength,latlength,altlength,miss_val,  w3d,  wbar2d)
400 call moyzon2(lonlength,latlength,altlength,miss_val,ang3d,angbar2d)
401   vbar(:,:,itim) =   vbar2d(:,:)
402   wbar(:,:,itim) =   wbar2d(:,:)
403 angbar(:,:,itim) = angbar2d(:,:)
404enddo ! timelength
405
406do ilon=1,lonlength
407 do ilat=1,latlength
408  do ilev=1,altlength
409   do itim=1,timelength
410    if ((vitv(ilon,ilat,ilev,itim).ne.miss_val).and. &
411        (vbar(ilat,ilev,itim)     .ne.miss_val)) then
412   vstar(ilon,ilat,ilev,itim) = vitv(ilon,ilat,ilev,itim)-  vbar(ilat,ilev,itim)
413    else
414   vstar(ilon,ilat,ilev,itim) = miss_val
415    endif
416    if ((vitw(ilon,ilat,ilev,itim).ne.miss_val).and. &
417        (wbar(ilat,ilev,itim)     .ne.miss_val)) then
418   wstar(ilon,ilat,ilev,itim) = vitw(ilon,ilat,ilev,itim)-  wbar(ilat,ilev,itim)
419    else
420   wstar(ilon,ilat,ilev,itim) = miss_val
421    endif
422    if ((  tsam(ilon,ilat,ilev,itim).ne.miss_val).and. &
423        (angbar(ilat,ilev,itim)     .ne.miss_val)) then
424 angstar(ilon,ilat,ilev,itim) = tsam(ilon,ilat,ilev,itim)-angbar(ilat,ilev,itim)
425    else
426 angstar(ilon,ilat,ilev,itim) = miss_val
427    endif
428   enddo
429  enddo
430 enddo
431enddo ! lonlength
432call moytim(lonlength,latlength,altlength,timelength,miss_val,  vstar,  vstarmean)
433call moytim(lonlength,latlength,altlength,timelength,miss_val,  wstar,  wstarmean)
434call moytim(lonlength,latlength,altlength,timelength,miss_val,angstar,angstarmean)
435do ilon=1,lonlength
436 do ilat=1,latlength
437  do ilev=1,altlength
438    if ((angstarmean(ilon,ilat,ilev).ne.miss_val).and. &
439        (  vstarmean(ilon,ilat,ilev).ne.miss_val)) then
440angstarmeanvstarmean(ilon,ilat,ilev) = angstarmean(ilon,ilat,ilev)*vstarmean(ilon,ilat,ilev)
441    else
442angstarmeanvstarmean(ilon,ilat,ilev) = miss_val
443    endif
444    if ((angstarmean(ilon,ilat,ilev).ne.miss_val).and. &
445        (  wstarmean(ilon,ilat,ilev).ne.miss_val)) then
446angstarmeanwstarmean(ilon,ilat,ilev) = angstarmean(ilon,ilat,ilev)*wstarmean(ilon,ilat,ilev)
447    else
448angstarmeanwstarmean(ilon,ilat,ilev) = miss_val
449    endif
450  enddo
451 enddo
452enddo ! lonlength
453
454call moytim(lonlength,latlength,altlength,timelength,miss_val,vitv,  vmean)
455call moytim(lonlength,latlength,altlength,timelength,miss_val,vitw,  wmean)
456call moytim(lonlength,latlength,altlength,timelength,miss_val,tsam,angmean)
457call moyzon2(lonlength,latlength,altlength,miss_val,  vmean,  vmeanbar)
458call moyzon2(lonlength,latlength,altlength,miss_val,  wmean,  wmeanbar)
459call moyzon2(lonlength,latlength,altlength,miss_val,angmean,angmeanbar)
460
461do ilon=1,lonlength
462 do ilat=1,latlength
463  do ilev=1,altlength
464   do itim=1,timelength
465    if ((vitv(ilon,ilat,ilev,itim).ne.miss_val).and. &
466        (tsam(ilon,ilat,ilev,itim).ne.miss_val)) then
467vang(ilon,ilat,ilev,itim) = vitv(ilon,ilat,ilev,itim)*tsam(ilon,ilat,ilev,itim)
468    else
469vang(ilon,ilat,ilev,itim) = miss_val
470    endif
471    if ((vitw(ilon,ilat,ilev,itim).ne.miss_val).and. &
472        (tsam(ilon,ilat,ilev,itim).ne.miss_val)) then
473wang(ilon,ilat,ilev,itim) = vitw(ilon,ilat,ilev,itim)*tsam(ilon,ilat,ilev,itim)
474    else
475wang(ilon,ilat,ilev,itim) = miss_val
476    endif
477   enddo
478  enddo
479 enddo
480enddo ! lonlength
481call moytim(lonlength,latlength,altlength,timelength,miss_val,vang,vangmean)
482call moytim(lonlength,latlength,altlength,timelength,miss_val,wang,wangmean)
483
484do ilon=1,lonlength
485 do ilat=1,latlength
486  do ilev=1,altlength
487   do itim=1,timelength
488    if ((vitv(ilon,ilat,ilev,itim).ne.miss_val).and. &
489        (vmean(ilon,ilat,ilev)    .ne.miss_val)) then
490  vprim(ilon,ilat,ilev,itim) = vitv(ilon,ilat,ilev,itim)-  vmean(ilon,ilat,ilev)
491    else
492  vprim(ilon,ilat,ilev,itim) = miss_val
493    endif
494    if ((vitw(ilon,ilat,ilev,itim).ne.miss_val).and. &
495        (wmean(ilon,ilat,ilev)    .ne.miss_val)) then
496  wprim(ilon,ilat,ilev,itim) = vitw(ilon,ilat,ilev,itim)-  wmean(ilon,ilat,ilev)
497    else
498  wprim(ilon,ilat,ilev,itim) = miss_val
499    endif
500    if ((tsam(ilon,ilat,ilev,itim).ne.miss_val).and. &
501        (angmean(ilon,ilat,ilev)  .ne.miss_val)) then
502angprim(ilon,ilat,ilev,itim) = tsam(ilon,ilat,ilev,itim)-angmean(ilon,ilat,ilev)
503    else
504angprim(ilon,ilat,ilev,itim) = miss_val
505    endif
506    if ((angprim(ilon,ilat,ilev,itim).ne.miss_val).and. &
507        (  vprim(ilon,ilat,ilev,itim).ne.miss_val)) then
508angprimvprim(ilon,ilat,ilev,itim) = angprim(ilon,ilat,ilev,itim)*vprim(ilon,ilat,ilev,itim)
509    else
510angprimvprim(ilon,ilat,ilev,itim) = miss_val
511    endif
512    if ((angprim(ilon,ilat,ilev,itim).ne.miss_val).and. &
513        (  wprim(ilon,ilat,ilev,itim).ne.miss_val)) then
514angprimwprim(ilon,ilat,ilev,itim) = angprim(ilon,ilat,ilev,itim)*wprim(ilon,ilat,ilev,itim)
515    else
516angprimwprim(ilon,ilat,ilev,itim) = miss_val
517    endif
518   enddo
519  enddo
520 enddo
521enddo ! lonlength
522call moytim(lonlength,latlength,altlength,timelength,miss_val,&
523                  angprimvprim,angprimvprimmean)
524call moytim(lonlength,latlength,altlength,timelength,miss_val,&
525                  angprimwprim,angprimwprimmean)
526
527! ang transport terms
528!----------------------
529
530call moyzon2(lonlength,latlength,altlength,miss_val,vangmean,totvang)
531call moyzon2(lonlength,latlength,altlength,miss_val,wangmean,totwang)
532
533do ilat=1,latlength
534 do ilev=1,altlength
535    if ((angmeanbar(ilat,ilev).ne.miss_val).and. &
536        (  vmeanbar(ilat,ilev).ne.miss_val)) then
537mmcvang(ilat,ilev) = angmeanbar(ilat,ilev)*vmeanbar(ilat,ilev)
538    else
539mmcvang(ilat,ilev) = miss_val
540    endif
541    if ((angmeanbar(ilat,ilev).ne.miss_val).and. &
542        (  wmeanbar(ilat,ilev).ne.miss_val)) then
543mmcwang(ilat,ilev) = angmeanbar(ilat,ilev)*wmeanbar(ilat,ilev)
544    else
545mmcwang(ilat,ilev) = miss_val
546    endif
547 enddo
548enddo
549
550call moyzon2(lonlength,latlength,altlength,miss_val,angprimvprimmean,trsvang)
551call moyzon2(lonlength,latlength,altlength,miss_val,angprimwprimmean,trswang)
552
553call moyzon2(lonlength,latlength,altlength,miss_val,angstarmeanvstarmean,stnvang)
554call moyzon2(lonlength,latlength,altlength,miss_val,angstarmeanwstarmean,stnwang)
555
556
557print*,"End of computations"
558
559!===============================================================================
560! 3. Create output file
561!===============================================================================
562
563! Create output file
564ierr=NF_CREATE(outfile,NF_CLOBBER,outfid)
565if (ierr.ne.NF_NOERR) then
566  write(*,*)"Error: could not create file ",outfile
567  stop
568endif
569
570!===============================================================================
571! 3.1. Define and write dimensions
572!===============================================================================
573
574call write_dim(outfid,lonlength,latlength,altlength,timelength, &
575    lon,lat,plev,time,lon_dimid,lat_dimid,alt_dimid,time_dimid)
576
577!===============================================================================
578! 3.2. Define and write variables
579!===============================================================================
580
581datashape2d(1)=lat_dimid
582datashape2d(2)=alt_dimid
583
584call write_var2d(outfid,datashape2d,latlength,altlength,&
585                "totvang   ", "tot hor trpt of ang ","m3 s-2    ",miss_val,&
586                 totvang )
587
588call write_var2d(outfid,datashape2d,latlength,altlength,&
589                "totwang   ", "tot ver trpt of ang ","m3 s-2    ",miss_val,&
590                 totwang )
591
592call write_var2d(outfid,datashape2d,latlength,altlength,&
593                "mmcvang   ", "MMC hor trpt of ang ","m3 s-2    ",miss_val,&
594                 mmcvang )
595
596call write_var2d(outfid,datashape2d,latlength,altlength,&
597                "mmcwang   ", "MMC ver trpt of ang ","m3 s-2    ",miss_val,&
598                 mmcwang )
599
600call write_var2d(outfid,datashape2d,latlength,altlength,&
601                "trsvang   ", "trs hor trpt of ang ","m3 s-2    ",miss_val,&
602                 trsvang )
603
604call write_var2d(outfid,datashape2d,latlength,altlength,&
605                "trswang   ", "trs ver trpt of ang ","m3 s-2    ",miss_val,&
606                 trswang )
607
608call write_var2d(outfid,datashape2d,latlength,altlength,&
609                "stnvang   ", "stn hor trpt of ang ","m3 s-2    ",miss_val,&
610                 stnvang )
611
612call write_var2d(outfid,datashape2d,latlength,altlength,&
613                "stnwang   ", "stn ver trpt of ang ","m3 s-2    ",miss_val,&
614                 stnwang )
615
616call write_var2d(outfid,datashape2d,latlength,altlength,&
617                "dmass     ", "mass in each cell   ","kg        ",miss_val,&
618                 dmassmeanbar )
619
620
621!!!! Close output file
622ierr=NF_CLOSE(outfid)
623if (ierr.ne.NF_NOERR) write(*,*) 'Error, failed to close output file ',outfile
624
625
626end program
Note: See TracBrowser for help on using the repository browser.