source: trunk/LMDZ.VENUS/Tools/tem.F90 @ 937

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

SL: tools for postprocessing (Veznus and Titan); see DOC/documentation/vt-tools.pdf
File size: 18.5 KB
Line 
1program tem
2
3! SL 01/2010:
4! This program reads 4D (lon-lat-alt-time) fields recast in log P coordinates
5! Developed from the tool built by Audrey Crespin during her PhD.
6!
7! it computes TransEulerianMean variables:
8!
9! vtem   -- 3D -- Residual meridional speed (m s-1)
10! wtem   -- 3D -- Residual   vertical speed (Pa s-1)
11! psitem -- 3D -- Residual stream function (kg s-1)
12! epfy   -- 3D -- meridional component of Eliassen-Palm flux
13! epfz   -- 3D -- vertical component of Eliassen-Palm flux
14! divepf -- 3D -- Divergence of Eliassen-Palm flux
15! ammctem - 3D -- Acc due to residual MMC
16!
17! Minimal requirements and dependencies:
18! The dataset must include the following data:
19! - pressure vertical coordinate
20! - surface pressure
21! - atmospheric temperature
22! - zonal, meridional and vertical (Pa/s) winds
23! - altitude above areoid
24
25implicit none
26
27include "netcdf.inc" ! NetCDF definitions
28
29character (len=128) :: infile ! input file name (name_P.nc)
30character (len=128) :: outfile ! output file name
31
32character (len=64) :: text ! to store some text
33integer infid ! NetCDF input file ID
34integer outfid ! NetCDF output file ID
35integer lon_dimid,lat_dimid,alt_dimid,time_dimid ! NetCDF dimension IDs
36integer lon_varid,lat_varid,alt_varid,time_varid
37integer,dimension(3) :: datashape3d ! shape of 3D datasets
38
39real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
40real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
41real :: pi
42real,dimension(:),allocatable :: lon ! longitude
43integer lonlength ! # of grid points along longitude
44real,dimension(:),allocatable :: lat ! latitude
45real,dimension(:),allocatable :: latrad ! latitude in radian
46integer latlength ! # of grid points along latitude
47real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
48integer altlength ! # of grid point along altitude (of input datasets)
49real,dimension(:),allocatable :: time ! time
50integer timelength ! # of points along time
51real,dimension(:,:,:),allocatable :: ps ! surface pressure
52real,dimension(:,:,:,:),allocatable :: temp ! atmospheric temperature
53real,dimension(:,:,:,:),allocatable :: vitu ! zonal wind (in m/s)
54real,dimension(:,:,:,:),allocatable :: vitv ! meridional wind (in m/s)
55real,dimension(:,:,:,:),allocatable :: vitw ! vertical wind (in Pa/s, then converted in m/s)
56real,dimension(:,:,:,:),allocatable :: za ! above areoid levels (m)
57
58!!! output variables
59real,dimension(:,:,:),allocatable :: epy ! merid component of EP flux
60real,dimension(:,:,:),allocatable :: epz ! verti component of EP flux
61real,dimension(:,:,:),allocatable :: divep ! divergence of EP flux
62real,dimension(:,:,:),allocatable :: ammctem ! acc by residual mmc
63real,dimension(:,:,:),allocatable :: uzon ! mean zonal wind
64real,dimension(:,:,:),allocatable :: vtem ! residual merid wind
65real,dimension(:,:,:),allocatable :: wtem ! residual verti wind
66real,dimension(:,:,:),allocatable :: psitem ! residual stream function
67
68! variables prepared for computation (4D)
69real,dimension(:,:,:,:),allocatable :: rayon ! distance to center (m)
70real,dimension(:,:,:,:),allocatable :: grav ! gravity field (m s-2)
71real,dimension(:,:,:,:),allocatable :: dmass ! mass in cell (kg)
72
73! variables prepared for computation inside timeloop
74real,dimension(:,:,:),allocatable :: r3d ! distance to center (m)
75real,dimension(:,:,:),allocatable :: rsurg ! rayon/grav
76real,dimension(:,:,:),allocatable :: t3d ! temp
77real,dimension(:,:,:),allocatable :: u3d ! zonal wind
78real,dimension(:,:,:),allocatable :: v3d ! merid wind
79real,dimension(:,:,:),allocatable :: w3d ! verti wind
80real,dimension(:,:,:),allocatable :: pk3d ! Exner function
81real,dimension(:,:,:),allocatable :: teta ! potential temp
82! variables obtained from computation inside timeloop
83real,dimension(:,:),allocatable :: epy2d ! merid component of EP flux
84real,dimension(:,:),allocatable :: epz2d ! verti component of EP flux
85real,dimension(:,:),allocatable :: div2d ! divergence of EP flux
86real,dimension(:,:),allocatable :: ammc2d ! acc by residual mmc
87real,dimension(:,:),allocatable :: ubar   ! mean zonal wind
88real,dimension(:,:),allocatable :: vtem2d ! residual merid wind
89real,dimension(:,:),allocatable :: wtem2d ! residual verti wind
90
91real,dimension(:,:),allocatable :: rbar   ! distance to center (zonal ave)
92real,dimension(:,:),allocatable :: rsurgbar ! rayon/grav
93real,dimension(:,:),allocatable :: vm   ! merid mass flux (zonal ave)
94real,dimension(:,:),allocatable :: psi  ! residual stream function
95real :: deltalat,deltalon ! lat and lon intervals in radians
96real,dimension(:,:,:),allocatable :: deltap ! pressure thickness of each layer (Pa)
97
98integer ierr,ierr1,ierr2 ! NetCDF routines return codes
99integer i,j,ilon,ilat,ilev,itim ! for loops
100logical :: lmdflag
101
102include "planet.h"
103
104!===============================================================================
105! 1. Input parameters
106!===============================================================================
107
108pi = 2.*asin(1.)
109
110write(*,*) ""
111write(*,*) "You are working on the atmosphere of ",planet
112
113!===============================================================================
114! 1.1 Input file
115!===============================================================================
116
117write(*,*) ""
118write(*,*) "Program valid for files with pressure axis (*_P.nc)"
119write(*,*) "Enter input file name:"
120
121read(*,'(a128)') infile
122write(*,*) ""
123
124! open input file
125
126ierr = NF_OPEN(infile,NF_NOWRITE,infid)
127if (ierr.ne.NF_NOERR) then
128   write(*,*) 'ERROR: Pb opening file ',trim(infile)
129   stop ""
130endif
131
132!===============================================================================
133! 1.2 Get grids in lon,lat,alt(pressure),time
134!===============================================================================
135
136call get_iddim(infid,lat_varid,latlength,lon_varid,lonlength,&
137                     alt_varid,altlength,time_varid,timelength,lmdflag )
138
139allocate(lon(lonlength))
140ierr=NF_GET_VAR_REAL(infid,lon_varid,lon)
141if (ierr.ne.NF_NOERR) stop "Error: Failed reading longitude"
142
143allocate(lat(latlength))
144ierr=NF_GET_VAR_REAL(infid,lat_varid,lat)
145if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
146
147allocate(latrad(latlength))
148latrad = lat*pi/180.
149
150! Lat, lon pressure intervals
151deltalat = abs(latrad(2)-latrad(1))
152deltalon = abs(lon(2)-lon(1))*pi/180.
153
154allocate(plev(altlength))
155ierr=NF_GET_VAR_REAL(infid,alt_varid,plev)
156if (ierr.ne.NF_NOERR) stop "Error: Failed reading altitude (ie pressure levels)"
157
158allocate(time(timelength))
159ierr=NF_GET_VAR_REAL(infid,time_varid,time)
160if (ierr.ne.NF_NOERR) stop "Error: Failed reading time"
161
162!===============================================================================
163! 1.3 Get output file name
164!===============================================================================
165write(*,*) ""
166!write(*,*) "Enter output file name"
167!read(*,*) outfile
168outfile=infile(1:len_trim(infile)-3)//"_TEM.nc"
169write(*,*) "Output file name is: "//trim(outfile)
170
171
172
173!===============================================================================
174! 2.1 Store needed fields
175!===============================================================================
176
177!===============================================================================
178! 2.1.1 Surface pressure
179!===============================================================================
180allocate(ps(lonlength,latlength,timelength))
181
182text="ps"
183call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
184if (ierr1.ne.NF_NOERR) then
185  write(*,*) "  looking for psol instead... "
186  text="psol"
187  call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
188  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get psol ID"
189endif
190if (ierr2.ne.NF_NOERR) stop "Error: Failed reading surface pressure"
191
192!===============================================================================
193! 2.1.2 Atmospheric temperature
194!===============================================================================
195allocate(temp(lonlength,latlength,altlength,timelength))
196
197text="temp"
198call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
199if (ierr1.ne.NF_NOERR) then
200  write(*,*) "  looking for t instead... "
201  text="t"
202  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
203  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
204endif
205if (ierr2.ne.NF_NOERR) stop "Error: Failed reading temperature"
206
207!===============================================================================
208! 2.1.3 Winds
209!===============================================================================
210allocate(vitu(lonlength,latlength,altlength,timelength))
211allocate(vitv(lonlength,latlength,altlength,timelength))
212allocate(vitw(lonlength,latlength,altlength,timelength))
213
214! zonal wind vitu (in m/s)
215text="vitu"
216call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
217if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
218if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
219
220! meridional wind vitv (in m/s)
221text="vitv"
222call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
223if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
224if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
225
226! vertical wind vitw (in Pa/s)
227text="vitw"
228call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
229if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitw ID"
230if (ierr2.ne.NF_NOERR) stop "Error: Failed reading vertical wind"
231
232!===============================================================================
233! 2.1.4 Altitude above areoide
234!===============================================================================
235! Only needed if g(z) on Titan...
236
237! allocate(za(lonlength,latlength,altlength,timelength))
238
239! text="zareoid"
240! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
241! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
242! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
243
244!===============================================================================
245!!! Allocations before timeloop
246!===============================================================================
247
248! latlength correspond a jjm+1
249! mais lonlength correspond a iim
250! pour boucler en longitude, on a besoin du point iim+1 (= 1)
251
252allocate(rayon(lonlength+1,latlength,altlength,timelength))
253allocate(grav(lonlength+1,latlength,altlength,timelength))
254allocate(dmass(lonlength+1,latlength,altlength,timelength))
255
256allocate(r3d(lonlength+1,latlength,altlength))
257allocate(rsurg(lonlength+1,latlength,altlength))
258allocate(rbar(latlength,altlength))
259allocate(rsurgbar(latlength,altlength))
260
261allocate(t3d(lonlength+1,latlength,altlength))
262allocate(u3d(lonlength+1,latlength,altlength))
263allocate(v3d(lonlength+1,latlength,altlength))
264allocate(w3d(lonlength+1,latlength,altlength))
265allocate(pk3d(lonlength+1,latlength,altlength))
266allocate(teta(lonlength+1,latlength,altlength))
267
268allocate(epy(latlength,altlength,timelength))
269allocate(epz(latlength,altlength,timelength))
270allocate(divep(latlength,altlength,timelength))
271allocate(ammctem(latlength,altlength,timelength))
272allocate(uzon(latlength,altlength,timelength))
273allocate(vtem(latlength,altlength,timelength))
274allocate(wtem(latlength,altlength,timelength))
275
276allocate(epy2d(latlength,altlength))
277allocate(epz2d(latlength,altlength))
278allocate(div2d(latlength,altlength))
279allocate(ammc2d(latlength,altlength))
280allocate(ubar(latlength,altlength))
281allocate(vtem2d(latlength,altlength))
282allocate(wtem2d(latlength,altlength))
283
284allocate(vm(latlength,altlength))
285allocate(psi(latlength,altlength))
286allocate(psitem(latlength,altlength,timelength))
287
288!===============================================================================
289! 2.2.2 Mass in cells
290!===============================================================================
291
292do itim=1,timelength
293do ilon=1,lonlength
294 do ilat=1,latlength
295  do ilev=1,altlength
296! Need to be consistent with GCM computations
297!    if (za(ilon,ilat,ilev,itim).ne.miss_val) then
298     rayon(ilon,ilat,ilev,itim) = a0
299!     rayon(ilon,ilat,ilev,itim) = za(ilon,ilat,ilev,itim) + a0
300      grav(ilon,ilat,ilev,itim) = g0*a0*a0 &
301                 /(rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim))
302!    else
303!     rayon(ilon,ilat,ilev,itim) = miss_val
304!      grav(ilon,ilat,ilev,itim) = miss_val
305!    endif
306  enddo
307 enddo
308enddo
309enddo ! timelength
310
311rayon(lonlength+1,:,:,:) = rayon(1,:,:,:)
312 grav(lonlength+1,:,:,:) =  grav(1,:,:,:)
313
314call cellmass(infid,latlength,lonlength+1,altlength,timelength, &
315              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
316              dmass )
317
318!===============================================================================
319!!! GLOBAL TIME LOOP !!!
320!===============================================================================
321do itim=1,timelength
322
323!===============================================================================
324! 2.2 Computations
325!===============================================================================
326
327!===============================================================================
328! 2.2.3 Init of 3D variables
329!===============================================================================
330
331do ilon=1,lonlength
332 do ilat=1,latlength
333  do ilev=1,altlength
334      r3d(ilon,ilat,ilev) = rayon(ilon,ilat,ilev,itim)
335      t3d(ilon,ilat,ilev) = temp(ilon,ilat,ilev,itim)
336      u3d(ilon,ilat,ilev) = vitu(ilon,ilat,ilev,itim)
337      v3d(ilon,ilat,ilev) = vitv(ilon,ilat,ilev,itim)
338      w3d(ilon,ilat,ilev) = vitw(ilon,ilat,ilev,itim)
339     pk3d(ilon,ilat,ilev) = cp0*(plev(ilev)/psref)**(R0/cp0)
340  enddo
341 enddo
342enddo
343
344 t3d(lonlength+1,:,:) =  t3d(1,:,:)
345 u3d(lonlength+1,:,:) =  u3d(1,:,:)
346 v3d(lonlength+1,:,:) =  v3d(1,:,:)
347 w3d(lonlength+1,:,:) =  w3d(1,:,:)
348pk3d(lonlength+1,:,:) = pk3d(1,:,:)
349
350call t2tpot((lonlength+1)*latlength*altlength,t3d,teta,pk3d)
351
352!===============================================================================
353! 2.2.4 TEM and Eliassen-Palm
354!===============================================================================
355
356print*,"eliasflu_meridien",itim
357
358call moyzon(lonlength,latlength,altlength,miss_val,r3d,rbar)
359call moyzon(lonlength,latlength,altlength,miss_val,u3d,ubar)
360
361call epflux(lonlength+1,latlength,altlength,miss_val,latrad,rbar &
362            ,teta,u3d,v3d,w3d,plev &
363            ,epy2d,epz2d,div2d,vtem2d,wtem2d,ammc2d &
364!           ,vpupbar2d,wpupbar2d,vpvpbar2d,wpvpbar2d &
365!           ,vptetapbar2d,wptetapbar2d &
366           )
367
368!===============================================================================
369! 2.2.5 Stream function
370!===============================================================================
371
372do ilon=1,lonlength+1
373 do ilat=1,latlength
374  do ilev=1,altlength
375    if (dmass(ilon,ilat,ilev,itim).ne.miss_val) then
376! rsurg: r*dp/g = dm/(r cos(lat) dlon dlat) !!!
377     rsurg(ilon,ilat,ilev) = dmass(ilon,ilat,ilev,itim) &
378          / (r3d(ilon,ilat,ilev)*cos(latrad(ilat))*deltalon*deltalat)
379    else
380     rsurg(ilon,ilat,ilev) = miss_val
381    endif
382  enddo
383 enddo
384enddo
385
386call moyzon(lonlength,latlength,altlength,miss_val,rsurg,rsurgbar)
387
388do ilat=1,latlength
389 do ilev=1,altlength
390    if (  (vtem2d(ilat,ilev).ne.miss_val).and. &
391        (rsurgbar(ilat,ilev).ne.miss_val) ) then
392      vm(ilat,ilev) = vtem2d(ilat,ilev) &
393            * 2.*pi*rsurgbar(ilat,ilev)*cos(latrad(ilat))
394    else
395      vm(ilat,ilev) = miss_val
396    endif
397 enddo
398enddo
399
400
401do ilat=1,latlength
402  psi(ilat,altlength) = 0.
403    if (vm(ilat,altlength).ne.miss_val) then
404      psi(ilat,altlength) = psi(ilat,altlength) &
405           + vm(ilat,altlength)
406    endif
407 do ilev=altlength-1,1,-1
408  psi(ilat,ilev) = psi(ilat,ilev+1)
409    if (vm(ilat,ilev).ne.miss_val) then
410      psi(ilat,ilev) = psi(ilat,ilev) &
411           + vm(ilat,ilev)
412    endif
413 enddo
414enddo
415
416!===============================================================================
417! 2.2.6 Building 2D+time variables
418!===============================================================================
419
420    epy(:,:,itim) = epy2d(:,:)
421    epz(:,:,itim) = epz2d(:,:)
422  divep(:,:,itim) = div2d(:,:)
423ammctem(:,:,itim) = ammc2d(:,:)
424   uzon(:,:,itim) =   ubar(:,:)
425   vtem(:,:,itim) = vtem2d(:,:)
426   wtem(:,:,itim) = wtem2d(:,:)
427 psitem(:,:,itim) = psi(:,:)
428
429
430enddo ! timelength
431!===============================================================================
432!!! END GLOBAL TIME LOOP !!!
433!===============================================================================
434
435print*,"End of computations"
436
437!===============================================================================
438! 3. Create output file
439!===============================================================================
440
441! Create output file
442ierr=NF_CREATE(outfile,NF_CLOBBER,outfid)
443if (ierr.ne.NF_NOERR) then
444  write(*,*)"Error: could not create file ",outfile
445  stop
446endif
447
448!===============================================================================
449! 3.1. Define and write dimensions
450!===============================================================================
451
452call write_dim(outfid,lonlength,latlength,altlength,timelength, &
453    lon,lat,plev,time,lon_dimid,lat_dimid,alt_dimid,time_dimid)
454
455!===============================================================================
456! 3.2. Define and write variables
457!===============================================================================
458
459datashape3d(1)=lat_dimid
460datashape3d(2)=alt_dimid
461datashape3d(3)=time_dimid
462
463call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
464                "epy       ", "EP flux on lat      ","m3 s-2    ",miss_val,&
465                 epy )
466
467call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
468                "epz       ", "EP flux on press    ","m3 s-2    ",miss_val,&
469                 epz )
470
471call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
472                "divep     ", "Div of EP flux      ","m s-2     ",miss_val,&
473                 divep )
474
475call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
476                "ammctem   ", "acc by residual mmc ","m s-2     ",miss_val,&
477                 ammctem )
478
479call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
480                "uzon      ", "Mean zonal wind     ","m s-1     ",miss_val,&
481                 uzon )
482
483call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
484                "vtem      ", "Resid TEM merid wind","m s-1     ",miss_val,&
485                 vtem )
486
487call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
488                "wtem      ", "Resid TEM verti wind","Pa s-1    ",miss_val,&
489                 wtem )
490
491call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
492                "psitem    ", "Resid stream funct  ","kg s-1    ",miss_val,&
493                 psitem )
494
495
496!!!! Close output file
497ierr=NF_CLOSE(outfid)
498if (ierr.ne.NF_NOERR) write(*,*) 'Error, failed to close output file ',outfile
499
500
501end program
Note: See TracBrowser for help on using the repository browser.