source: trunk/LMDZ.VENUS/Tools/tmc.F90 @ 3315

Last change on this file since 3315 was 1468, checked in by slebonnois, 9 years ago

SL: bug corrections
File size: 24.1 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 ! 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 :: coslat ! cos of latitude
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.)
143miss_val = miss_val_def
144
145write(*,*) ""
146write(*,*) "You are working on the atmosphere of ",planet
147
148!===============================================================================
149! 1.1 Input file
150!===============================================================================
151
152write(*,*) ""
153write(*,*) "Program valid for files with pressure axis (*_P.nc)"
154write(*,*) "Enter input file name:"
155
156read(*,'(a128)') infile
157write(*,*) ""
158
159! open input file
160
161ierr = NF_OPEN(infile,NF_NOWRITE,infid)
162if (ierr.ne.NF_NOERR) then
163   write(*,*) 'ERROR: Pb opening file ',trim(infile)
164   stop ""
165endif
166
167!===============================================================================
168! 1.2 Get grids in lon,lat,alt(pressure),time
169!===============================================================================
170
171call get_iddim(infid,lat_varid,latlength,lon_varid,lonlength,&
172                     alt_varid,altlength,time_varid,timelength,lmdflag )
173
174allocate(lon(lonlength))
175ierr=NF_GET_VAR_REAL(infid,lon_varid,lon)
176if (ierr.ne.NF_NOERR) stop "Error: Failed reading longitude"
177
178allocate(lat(latlength))
179ierr=NF_GET_VAR_REAL(infid,lat_varid,lat)
180if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
181
182allocate(coslat(latlength))
183! Beware of rounding problems at poles...
184coslat(:) = max(0.,cos(lat(:)*pi/180.))
185
186! Lat, lon pressure intervals
187deltalat = abs(lat(2)-lat(1))*pi/180.
188deltalon = abs(lon(2)-lon(1))*pi/180.
189
190allocate(plev(altlength))
191ierr=NF_GET_VAR_REAL(infid,alt_varid,plev)
192if (ierr.ne.NF_NOERR) stop "Error: Failed reading altitude (ie pressure levels)"
193
194allocate(time(timelength))
195ierr=NF_GET_VAR_REAL(infid,time_varid,time)
196if (ierr.ne.NF_NOERR) stop "Error: Failed reading time"
197
198!===============================================================================
199! 1.3 Get output file name
200!===============================================================================
201write(*,*) ""
202!write(*,*) "Enter output file name"
203!read(*,*) outfile
204outfile=infile(1:len_trim(infile)-3)//"_TMC.nc"
205write(*,*) "Output file name is: "//trim(outfile)
206
207
208
209!===============================================================================
210! 2.1 Store needed fields
211!===============================================================================
212
213!===============================================================================
214! 2.1.1 Surface pressure
215!===============================================================================
216allocate(ps(lonlength,latlength,timelength))
217
218text="ps"
219call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
220if (ierr1.ne.NF_NOERR) then
221  write(*,*) "  looking for psol instead... "
222  text="psol"
223  call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
224  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get psol ID"
225endif
226if (ierr2.ne.NF_NOERR) stop "Error: Failed reading surface pressure"
227
228!===============================================================================
229! 2.1.3 Winds
230!===============================================================================
231allocate(vitu(lonlength,latlength,altlength,timelength))
232allocate(vitv(lonlength,latlength,altlength,timelength))
233allocate(vitw(lonlength,latlength,altlength,timelength))
234
235! zonal wind vitu (in m/s)
236text="vitu"
237call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
238if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
239if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
240
241! meridional wind vitv (in m/s)
242text="vitv"
243call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
244if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
245if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
246
247! vertical wind vitw (in Pa/s)
248text="vitw"
249call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
250if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitw ID"
251if (ierr2.ne.NF_NOERR) stop "Error: Failed reading vertical wind"
252
253!===============================================================================
254! 2.1.4 Altitude above areoide
255!===============================================================================
256! Only needed if g(z) on Titan...
257
258! allocate(za(lonlength,latlength,altlength,timelength))
259
260! text="zareoid"
261! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
262! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
263! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
264
265!===============================================================================
266! 2.2 Computations
267!===============================================================================
268
269!===============================================================================
270! 2.2.2 Mass in cells
271!===============================================================================
272allocate(rayon(lonlength,latlength,altlength,timelength))
273allocate( grav(lonlength,latlength,altlength,timelength))
274allocate(dmass(lonlength,latlength,altlength,timelength))
275allocate(dmassmean(lonlength,latlength,altlength))
276allocate(dmassmeanbar(latlength,altlength))
277
278do itim=1,timelength
279do ilon=1,lonlength
280 do ilat=1,latlength
281  do ilev=1,altlength
282! Need to be consistent with GCM computations
283    if (vitu(ilon,ilat,ilev,itim).ne.miss_val) then
284     rayon(ilon,ilat,ilev,itim) = a0
285!     rayon(ilon,ilat,ilev,itim) = za(ilon,ilat,ilev,itim) + a0
286      grav(ilon,ilat,ilev,itim) = g0*a0*a0 &
287                 /(rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim))
288    else
289     rayon(ilon,ilat,ilev,itim) = miss_val
290      grav(ilon,ilat,ilev,itim) = miss_val
291    endif
292  enddo
293 enddo
294enddo
295enddo ! timelength
296
297call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
298              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
299              dmass )
300
301call moytim(lonlength,latlength,altlength,timelength,miss_val,dmass,dmassmean)
302call moyzon2(lonlength,latlength,altlength,miss_val,dmassmean,dmassmeanbar)
303
304print*,"OK dmass"
305
306!===============================================================================
307! 2.2.3 Specific angular momentum
308!===============================================================================
309allocate(osam(lonlength,latlength,altlength,timelength))
310allocate(rsam(lonlength,latlength,altlength,timelength))
311allocate(tsam(lonlength,latlength,altlength,timelength))
312
313do itim=1,timelength
314do ilon=1,lonlength
315 do ilat=1,latlength
316  do ilev=1,altlength
317    if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
318      osam(ilon,ilat,ilev,itim) = &
319         rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim) &
320       * coslat(ilat)*coslat(ilat) &
321       * omega
322      rsam(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim) &
323       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)
324      tsam(ilon,ilat,ilev,itim) = osam(ilon,ilat,ilev,itim)&
325                                + rsam(ilon,ilat,ilev,itim)
326    else
327      osam(ilon,ilat,ilev,itim) = miss_val
328      rsam(ilon,ilat,ilev,itim) = miss_val
329      tsam(ilon,ilat,ilev,itim) = miss_val
330    endif
331  enddo
332 enddo
333enddo
334enddo ! timelength
335
336print*,"debut tprt ang"
337
338!===============================================================================
339! 2.2.4 Angular momentum transport
340!===============================================================================
341! allocations
342!-------------
343allocate(vang(lonlength,latlength,altlength,timelength))
344allocate(wang(lonlength,latlength,altlength,timelength))
345allocate(  vstar(lonlength,latlength,altlength,timelength))
346allocate(  wstar(lonlength,latlength,altlength,timelength))
347allocate(angstar(lonlength,latlength,altlength,timelength))
348allocate(  vprim(lonlength,latlength,altlength,timelength))
349allocate(  wprim(lonlength,latlength,altlength,timelength))
350allocate(angprim(lonlength,latlength,altlength,timelength))
351allocate(angprimvprim(lonlength,latlength,altlength,timelength))
352allocate(angprimwprim(lonlength,latlength,altlength,timelength))
353
354! lon,lat,alt
355allocate(vangmean(lonlength,latlength,altlength))
356allocate(wangmean(lonlength,latlength,altlength))
357allocate(  vmean(lonlength,latlength,altlength))
358allocate(  wmean(lonlength,latlength,altlength))
359allocate(angmean(lonlength,latlength,altlength))
360allocate(angprimvprimmean(lonlength,latlength,altlength))
361allocate(angprimwprimmean(lonlength,latlength,altlength))
362allocate(  vstarmean(lonlength,latlength,altlength))
363allocate(  wstarmean(lonlength,latlength,altlength))
364allocate(angstarmean(lonlength,latlength,altlength))
365allocate(angstarmeanvstarmean(lonlength,latlength,altlength))
366allocate(angstarmeanwstarmean(lonlength,latlength,altlength))
367allocate(  v3d(lonlength,latlength,altlength))
368allocate(  w3d(lonlength,latlength,altlength))
369allocate(ang3d(lonlength,latlength,altlength))
370
371! lat,alt,time
372allocate(  vbar(latlength,altlength,timelength))
373allocate(  wbar(latlength,altlength,timelength))
374allocate(angbar(latlength,altlength,timelength))
375
376!lat,alt
377allocate(  vmeanbar(latlength,altlength))
378allocate(  wmeanbar(latlength,altlength))
379allocate(angmeanbar(latlength,altlength))
380allocate(  vbar2d(latlength,altlength))
381allocate(  wbar2d(latlength,altlength))
382allocate(angbar2d(latlength,altlength))
383
384allocate(totvang(latlength,altlength))
385allocate(totwang(latlength,altlength))
386allocate(mmcvang(latlength,altlength))
387allocate(mmcwang(latlength,altlength))
388allocate(trsvang(latlength,altlength))
389allocate(trswang(latlength,altlength))
390allocate(stnvang(latlength,altlength))
391allocate(stnwang(latlength,altlength))
392
393! intermediates
394!-----------------
395
396do itim=1,timelength
397   v3d(:,:,:) = vitv(:,:,:,itim)
398   w3d(:,:,:) = vitw(:,:,:,itim)
399 ang3d(:,:,:) = tsam(:,:,:,itim)
400 call moyzon2(lonlength,latlength,altlength,miss_val,  v3d,  vbar2d)
401 call moyzon2(lonlength,latlength,altlength,miss_val,  w3d,  wbar2d)
402 call moyzon2(lonlength,latlength,altlength,miss_val,ang3d,angbar2d)
403   vbar(:,:,itim) =   vbar2d(:,:)
404   wbar(:,:,itim) =   wbar2d(:,:)
405 angbar(:,:,itim) = angbar2d(:,:)
406enddo ! timelength
407
408do ilon=1,lonlength
409 do ilat=1,latlength
410  do ilev=1,altlength
411   do itim=1,timelength
412    if ((vitv(ilon,ilat,ilev,itim).ne.miss_val).and. &
413        (vbar(ilat,ilev,itim)     .ne.miss_val)) then
414   vstar(ilon,ilat,ilev,itim) = vitv(ilon,ilat,ilev,itim)-  vbar(ilat,ilev,itim)
415    else
416   vstar(ilon,ilat,ilev,itim) = miss_val
417    endif
418    if ((vitw(ilon,ilat,ilev,itim).ne.miss_val).and. &
419        (wbar(ilat,ilev,itim)     .ne.miss_val)) then
420   wstar(ilon,ilat,ilev,itim) = vitw(ilon,ilat,ilev,itim)-  wbar(ilat,ilev,itim)
421    else
422   wstar(ilon,ilat,ilev,itim) = miss_val
423    endif
424    if ((  tsam(ilon,ilat,ilev,itim).ne.miss_val).and. &
425        (angbar(ilat,ilev,itim)     .ne.miss_val)) then
426 angstar(ilon,ilat,ilev,itim) = tsam(ilon,ilat,ilev,itim)-angbar(ilat,ilev,itim)
427    else
428 angstar(ilon,ilat,ilev,itim) = miss_val
429    endif
430   enddo
431  enddo
432 enddo
433enddo ! lonlength
434call moytim(lonlength,latlength,altlength,timelength,miss_val,  vstar,  vstarmean)
435call moytim(lonlength,latlength,altlength,timelength,miss_val,  wstar,  wstarmean)
436call moytim(lonlength,latlength,altlength,timelength,miss_val,angstar,angstarmean)
437do ilon=1,lonlength
438 do ilat=1,latlength
439  do ilev=1,altlength
440    if ((angstarmean(ilon,ilat,ilev).ne.miss_val).and. &
441        (  vstarmean(ilon,ilat,ilev).ne.miss_val)) then
442angstarmeanvstarmean(ilon,ilat,ilev) = angstarmean(ilon,ilat,ilev)*vstarmean(ilon,ilat,ilev)
443    else
444angstarmeanvstarmean(ilon,ilat,ilev) = miss_val
445    endif
446    if ((angstarmean(ilon,ilat,ilev).ne.miss_val).and. &
447        (  wstarmean(ilon,ilat,ilev).ne.miss_val)) then
448angstarmeanwstarmean(ilon,ilat,ilev) = angstarmean(ilon,ilat,ilev)*wstarmean(ilon,ilat,ilev)
449    else
450angstarmeanwstarmean(ilon,ilat,ilev) = miss_val
451    endif
452  enddo
453 enddo
454enddo ! lonlength
455
456call moytim(lonlength,latlength,altlength,timelength,miss_val,vitv,  vmean)
457call moytim(lonlength,latlength,altlength,timelength,miss_val,vitw,  wmean)
458call moytim(lonlength,latlength,altlength,timelength,miss_val,tsam,angmean)
459call moyzon2(lonlength,latlength,altlength,miss_val,  vmean,  vmeanbar)
460call moyzon2(lonlength,latlength,altlength,miss_val,  wmean,  wmeanbar)
461call moyzon2(lonlength,latlength,altlength,miss_val,angmean,angmeanbar)
462
463do ilon=1,lonlength
464 do ilat=1,latlength
465  do ilev=1,altlength
466   do itim=1,timelength
467    if ((vitv(ilon,ilat,ilev,itim).ne.miss_val).and. &
468        (tsam(ilon,ilat,ilev,itim).ne.miss_val)) then
469vang(ilon,ilat,ilev,itim) = vitv(ilon,ilat,ilev,itim)*tsam(ilon,ilat,ilev,itim)
470    else
471vang(ilon,ilat,ilev,itim) = miss_val
472    endif
473    if ((vitw(ilon,ilat,ilev,itim).ne.miss_val).and. &
474        (tsam(ilon,ilat,ilev,itim).ne.miss_val)) then
475wang(ilon,ilat,ilev,itim) = vitw(ilon,ilat,ilev,itim)*tsam(ilon,ilat,ilev,itim)
476    else
477wang(ilon,ilat,ilev,itim) = miss_val
478    endif
479   enddo
480  enddo
481 enddo
482enddo ! lonlength
483call moytim(lonlength,latlength,altlength,timelength,miss_val,vang,vangmean)
484call moytim(lonlength,latlength,altlength,timelength,miss_val,wang,wangmean)
485
486do ilon=1,lonlength
487 do ilat=1,latlength
488  do ilev=1,altlength
489   do itim=1,timelength
490    if ((vitv(ilon,ilat,ilev,itim).ne.miss_val).and. &
491        (vmean(ilon,ilat,ilev)    .ne.miss_val)) then
492  vprim(ilon,ilat,ilev,itim) = vitv(ilon,ilat,ilev,itim)-  vmean(ilon,ilat,ilev)
493    else
494  vprim(ilon,ilat,ilev,itim) = miss_val
495    endif
496    if ((vitw(ilon,ilat,ilev,itim).ne.miss_val).and. &
497        (wmean(ilon,ilat,ilev)    .ne.miss_val)) then
498  wprim(ilon,ilat,ilev,itim) = vitw(ilon,ilat,ilev,itim)-  wmean(ilon,ilat,ilev)
499    else
500  wprim(ilon,ilat,ilev,itim) = miss_val
501    endif
502    if ((tsam(ilon,ilat,ilev,itim).ne.miss_val).and. &
503        (angmean(ilon,ilat,ilev)  .ne.miss_val)) then
504angprim(ilon,ilat,ilev,itim) = tsam(ilon,ilat,ilev,itim)-angmean(ilon,ilat,ilev)
505    else
506angprim(ilon,ilat,ilev,itim) = miss_val
507    endif
508    if ((angprim(ilon,ilat,ilev,itim).ne.miss_val).and. &
509        (  vprim(ilon,ilat,ilev,itim).ne.miss_val)) then
510angprimvprim(ilon,ilat,ilev,itim) = angprim(ilon,ilat,ilev,itim)*vprim(ilon,ilat,ilev,itim)
511    else
512angprimvprim(ilon,ilat,ilev,itim) = miss_val
513    endif
514    if ((angprim(ilon,ilat,ilev,itim).ne.miss_val).and. &
515        (  wprim(ilon,ilat,ilev,itim).ne.miss_val)) then
516angprimwprim(ilon,ilat,ilev,itim) = angprim(ilon,ilat,ilev,itim)*wprim(ilon,ilat,ilev,itim)
517    else
518angprimwprim(ilon,ilat,ilev,itim) = miss_val
519    endif
520   enddo
521  enddo
522 enddo
523enddo ! lonlength
524call moytim(lonlength,latlength,altlength,timelength,miss_val,&
525                  angprimvprim,angprimvprimmean)
526call moytim(lonlength,latlength,altlength,timelength,miss_val,&
527                  angprimwprim,angprimwprimmean)
528
529! ang transport terms
530!----------------------
531
532call moyzon2(lonlength,latlength,altlength,miss_val,vangmean,totvang)
533call moyzon2(lonlength,latlength,altlength,miss_val,wangmean,totwang)
534
535do ilat=1,latlength
536 do ilev=1,altlength
537    if ((angmeanbar(ilat,ilev).ne.miss_val).and. &
538        (  vmeanbar(ilat,ilev).ne.miss_val)) then
539mmcvang(ilat,ilev) = angmeanbar(ilat,ilev)*vmeanbar(ilat,ilev)
540    else
541mmcvang(ilat,ilev) = miss_val
542    endif
543    if ((angmeanbar(ilat,ilev).ne.miss_val).and. &
544        (  wmeanbar(ilat,ilev).ne.miss_val)) then
545mmcwang(ilat,ilev) = angmeanbar(ilat,ilev)*wmeanbar(ilat,ilev)
546    else
547mmcwang(ilat,ilev) = miss_val
548    endif
549 enddo
550enddo
551
552call moyzon2(lonlength,latlength,altlength,miss_val,angprimvprimmean,trsvang)
553call moyzon2(lonlength,latlength,altlength,miss_val,angprimwprimmean,trswang)
554
555call moyzon2(lonlength,latlength,altlength,miss_val,angstarmeanvstarmean,stnvang)
556call moyzon2(lonlength,latlength,altlength,miss_val,angstarmeanwstarmean,stnwang)
557
558
559print*,"End of computations"
560
561!===============================================================================
562! 3. Create output file
563!===============================================================================
564
565! Create output file
566ierr=NF_CREATE(outfile,NF_CLOBBER,outfid)
567if (ierr.ne.NF_NOERR) then
568  write(*,*)"Error: could not create file ",outfile
569  stop
570endif
571
572!===============================================================================
573! 3.1. Define and write dimensions
574!===============================================================================
575
576call write_dim(outfid,lonlength,latlength,altlength,timelength, &
577    lon,lat,plev,time,lon_dimid,lat_dimid,alt_dimid,time_dimid)
578
579!===============================================================================
580! 3.2. Define and write variables
581!===============================================================================
582
583datashape2d(1)=lat_dimid
584datashape2d(2)=alt_dimid
585
586call write_var2d(outfid,datashape2d,latlength,altlength,&
587                "totvang   ", "tot hor trpt of ang ","m3 s-2    ",miss_val,&
588                 totvang )
589
590call write_var2d(outfid,datashape2d,latlength,altlength,&
591                "totwang   ", "tot ver trpt of ang ","m3 s-2    ",miss_val,&
592                 totwang )
593
594call write_var2d(outfid,datashape2d,latlength,altlength,&
595                "mmcvang   ", "MMC hor trpt of ang ","m3 s-2    ",miss_val,&
596                 mmcvang )
597
598call write_var2d(outfid,datashape2d,latlength,altlength,&
599                "mmcwang   ", "MMC ver trpt of ang ","m3 s-2    ",miss_val,&
600                 mmcwang )
601
602call write_var2d(outfid,datashape2d,latlength,altlength,&
603                "trsvang   ", "trs hor trpt of ang ","m3 s-2    ",miss_val,&
604                 trsvang )
605
606call write_var2d(outfid,datashape2d,latlength,altlength,&
607                "trswang   ", "trs ver trpt of ang ","m3 s-2    ",miss_val,&
608                 trswang )
609
610call write_var2d(outfid,datashape2d,latlength,altlength,&
611                "stnvang   ", "stn hor trpt of ang ","m3 s-2    ",miss_val,&
612                 stnvang )
613
614call write_var2d(outfid,datashape2d,latlength,altlength,&
615                "stnwang   ", "stn ver trpt of ang ","m3 s-2    ",miss_val,&
616                 stnwang )
617
618call write_var2d(outfid,datashape2d,latlength,altlength,&
619                "dmass     ", "mass in each cell   ","kg        ",miss_val,&
620                 dmassmeanbar )
621
622
623!!!! Close output file
624ierr=NF_CLOSE(outfid)
625if (ierr.ne.NF_NOERR) write(*,*) 'Error, failed to close output file ',outfile
626
627
628end program
Note: See TracBrowser for help on using the repository browser.