source: trunk/LMDZ.TITAN.old/Tools/tem.F90 @ 3094

Last change on this file since 3094 was 1454, checked in by slebonnois, 10 years ago

SL: corrections in Titan and Venus tools
File size: 18.6 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 ! 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 :: coslat ! cos of latitude
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.)
109miss_val = miss_val_def
110
111write(*,*) ""
112write(*,*) "You are working on the atmosphere of ",planet
113
114!===============================================================================
115! 1.1 Input file
116!===============================================================================
117
118write(*,*) ""
119write(*,*) "Program valid for files with pressure axis (*_P.nc)"
120write(*,*) "Enter input file name:"
121
122read(*,'(a128)') infile
123write(*,*) ""
124
125! open input file
126
127ierr = NF_OPEN(infile,NF_NOWRITE,infid)
128if (ierr.ne.NF_NOERR) then
129   write(*,*) 'ERROR: Pb opening file ',trim(infile)
130   stop ""
131endif
132
133!===============================================================================
134! 1.2 Get grids in lon,lat,alt(pressure),time
135!===============================================================================
136
137call get_iddim(infid,lat_varid,latlength,lon_varid,lonlength,&
138                     alt_varid,altlength,time_varid,timelength,lmdflag )
139
140allocate(lon(lonlength))
141ierr=NF_GET_VAR_REAL(infid,lon_varid,lon)
142if (ierr.ne.NF_NOERR) stop "Error: Failed reading longitude"
143
144allocate(lat(latlength))
145ierr=NF_GET_VAR_REAL(infid,lat_varid,lat)
146if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
147
148allocate(coslat(latlength))
149! Beware of rounding problems at poles...
150coslat(:) = max(0.,cos(lat(:)*pi/180.))
151
152! Lat, lon pressure intervals
153deltalat = abs(lat(2)-lat(1))*pi/180.
154deltalon = abs(lon(2)-lon(1))*pi/180.
155
156allocate(plev(altlength))
157ierr=NF_GET_VAR_REAL(infid,alt_varid,plev)
158if (ierr.ne.NF_NOERR) stop "Error: Failed reading altitude (ie pressure levels)"
159
160allocate(time(timelength))
161ierr=NF_GET_VAR_REAL(infid,time_varid,time)
162if (ierr.ne.NF_NOERR) stop "Error: Failed reading time"
163
164!===============================================================================
165! 1.3 Get output file name
166!===============================================================================
167write(*,*) ""
168!write(*,*) "Enter output file name"
169!read(*,*) outfile
170outfile=infile(1:len_trim(infile)-3)//"_TEM.nc"
171write(*,*) "Output file name is: "//trim(outfile)
172
173
174
175!===============================================================================
176! 2.1 Store needed fields
177!===============================================================================
178
179!===============================================================================
180! 2.1.1 Surface pressure
181!===============================================================================
182allocate(ps(lonlength,latlength,timelength))
183
184text="ps"
185call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
186if (ierr1.ne.NF_NOERR) then
187  write(*,*) "  looking for psol instead... "
188  text="psol"
189  call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
190  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get psol ID"
191endif
192if (ierr2.ne.NF_NOERR) stop "Error: Failed reading surface pressure"
193
194!===============================================================================
195! 2.1.2 Atmospheric temperature
196!===============================================================================
197allocate(temp(lonlength,latlength,altlength,timelength))
198
199text="temp"
200call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
201if (ierr1.ne.NF_NOERR) then
202  write(*,*) "  looking for t instead... "
203  text="t"
204  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
205  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
206endif
207if (ierr2.ne.NF_NOERR) stop "Error: Failed reading temperature"
208
209!===============================================================================
210! 2.1.3 Winds
211!===============================================================================
212allocate(vitu(lonlength,latlength,altlength,timelength))
213allocate(vitv(lonlength,latlength,altlength,timelength))
214allocate(vitw(lonlength,latlength,altlength,timelength))
215
216! zonal wind vitu (in m/s)
217text="vitu"
218call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
219if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
220if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
221
222! meridional wind vitv (in m/s)
223text="vitv"
224call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
225if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
226if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
227
228! vertical wind vitw (in Pa/s)
229text="vitw"
230call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
231if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitw ID"
232if (ierr2.ne.NF_NOERR) stop "Error: Failed reading vertical wind"
233
234!===============================================================================
235! 2.1.4 Altitude above areoide
236!===============================================================================
237! Only needed if g(z) on Titan...
238
239! allocate(za(lonlength,latlength,altlength,timelength))
240
241! text="zareoid"
242! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
243! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
244! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
245
246!===============================================================================
247!!! Allocations before timeloop
248!===============================================================================
249
250! latlength correspond a jjm+1
251! mais lonlength correspond a iim
252! pour boucler en longitude, on a besoin du point iim+1 (= 1)
253
254allocate(rayon(lonlength+1,latlength,altlength,timelength))
255allocate(grav(lonlength+1,latlength,altlength,timelength))
256allocate(dmass(lonlength+1,latlength,altlength,timelength))
257
258allocate(r3d(lonlength+1,latlength,altlength))
259allocate(rsurg(lonlength+1,latlength,altlength))
260allocate(rbar(latlength,altlength))
261allocate(rsurgbar(latlength,altlength))
262
263allocate(t3d(lonlength+1,latlength,altlength))
264allocate(u3d(lonlength+1,latlength,altlength))
265allocate(v3d(lonlength+1,latlength,altlength))
266allocate(w3d(lonlength+1,latlength,altlength))
267allocate(pk3d(lonlength+1,latlength,altlength))
268allocate(teta(lonlength+1,latlength,altlength))
269
270allocate(epy(latlength,altlength,timelength))
271allocate(epz(latlength,altlength,timelength))
272allocate(divep(latlength,altlength,timelength))
273allocate(ammctem(latlength,altlength,timelength))
274allocate(uzon(latlength,altlength,timelength))
275allocate(vtem(latlength,altlength,timelength))
276allocate(wtem(latlength,altlength,timelength))
277
278allocate(epy2d(latlength,altlength))
279allocate(epz2d(latlength,altlength))
280allocate(div2d(latlength,altlength))
281allocate(ammc2d(latlength,altlength))
282allocate(ubar(latlength,altlength))
283allocate(vtem2d(latlength,altlength))
284allocate(wtem2d(latlength,altlength))
285
286allocate(vm(latlength,altlength))
287allocate(psi(latlength,altlength))
288allocate(psitem(latlength,altlength,timelength))
289
290!===============================================================================
291! 2.2.2 Mass in cells
292!===============================================================================
293
294do itim=1,timelength
295do ilon=1,lonlength
296 do ilat=1,latlength
297  do ilev=1,altlength
298! Need to be consistent with GCM computations
299!    if (za(ilon,ilat,ilev,itim).ne.miss_val) then
300     rayon(ilon,ilat,ilev,itim) = a0
301!     rayon(ilon,ilat,ilev,itim) = za(ilon,ilat,ilev,itim) + a0
302      grav(ilon,ilat,ilev,itim) = g0*a0*a0 &
303                 /(rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim))
304!    else
305!     rayon(ilon,ilat,ilev,itim) = miss_val
306!      grav(ilon,ilat,ilev,itim) = miss_val
307!    endif
308  enddo
309 enddo
310enddo
311enddo ! timelength
312
313rayon(lonlength+1,:,:,:) = rayon(1,:,:,:)
314 grav(lonlength+1,:,:,:) =  grav(1,:,:,:)
315
316call cellmass(infid,latlength,lonlength+1,altlength,timelength,lmdflag, &
317              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
318              dmass )
319
320!===============================================================================
321!!! GLOBAL TIME LOOP !!!
322!===============================================================================
323do itim=1,timelength
324
325!===============================================================================
326! 2.2 Computations
327!===============================================================================
328
329!===============================================================================
330! 2.2.3 Init of 3D variables
331!===============================================================================
332
333do ilon=1,lonlength
334 do ilat=1,latlength
335  do ilev=1,altlength
336      r3d(ilon,ilat,ilev) = rayon(ilon,ilat,ilev,itim)
337      t3d(ilon,ilat,ilev) = temp(ilon,ilat,ilev,itim)
338      u3d(ilon,ilat,ilev) = vitu(ilon,ilat,ilev,itim)
339      v3d(ilon,ilat,ilev) = vitv(ilon,ilat,ilev,itim)
340      w3d(ilon,ilat,ilev) = vitw(ilon,ilat,ilev,itim)
341     pk3d(ilon,ilat,ilev) = cp0*(plev(ilev)/psref)**(R0/cp0)
342  enddo
343 enddo
344enddo
345
346 t3d(lonlength+1,:,:) =  t3d(1,:,:)
347 u3d(lonlength+1,:,:) =  u3d(1,:,:)
348 v3d(lonlength+1,:,:) =  v3d(1,:,:)
349 w3d(lonlength+1,:,:) =  w3d(1,:,:)
350pk3d(lonlength+1,:,:) = pk3d(1,:,:)
351
352call t2tpot((lonlength+1)*latlength*altlength,t3d,teta,pk3d)
353
354!===============================================================================
355! 2.2.4 TEM and Eliassen-Palm
356!===============================================================================
357
358print*,"eliasflu_meridien",itim
359
360call moyzon(lonlength,latlength,altlength,miss_val,r3d,rbar)
361call moyzon(lonlength,latlength,altlength,miss_val,u3d,ubar)
362
363call epflux(lonlength+1,latlength,altlength,miss_val,lat,rbar &
364            ,teta,u3d,v3d,w3d,plev &
365            ,epy2d,epz2d,div2d,vtem2d,wtem2d,ammc2d &
366!           ,vpupbar2d,wpupbar2d,vpvpbar2d,wpvpbar2d &
367!           ,vptetapbar2d,wptetapbar2d &
368           )
369
370!===============================================================================
371! 2.2.5 Stream function
372!===============================================================================
373
374do ilon=1,lonlength+1
375 do ilat=1,latlength
376  do ilev=1,altlength
377    if (dmass(ilon,ilat,ilev,itim).ne.miss_val) then
378! rsurg: r*dp/g = dm/(r cos(lat) dlon dlat) !!!
379     rsurg(ilon,ilat,ilev) = dmass(ilon,ilat,ilev,itim) &
380          / (r3d(ilon,ilat,ilev)*coslat(ilat)*deltalon*deltalat)
381    else
382     rsurg(ilon,ilat,ilev) = miss_val
383    endif
384  enddo
385 enddo
386enddo
387
388call moyzon(lonlength,latlength,altlength,miss_val,rsurg,rsurgbar)
389
390do ilat=1,latlength
391 do ilev=1,altlength
392    if (  (vtem2d(ilat,ilev).ne.miss_val).and. &
393        (rsurgbar(ilat,ilev).ne.miss_val) ) then
394      vm(ilat,ilev) = vtem2d(ilat,ilev) &
395            * 2.*pi*rsurgbar(ilat,ilev)*coslat(ilat)
396    else
397      vm(ilat,ilev) = miss_val
398    endif
399 enddo
400enddo
401
402
403do ilat=1,latlength
404  psi(ilat,altlength) = 0.
405    if (vm(ilat,altlength).ne.miss_val) then
406      psi(ilat,altlength) = psi(ilat,altlength) &
407           + vm(ilat,altlength)
408    endif
409 do ilev=altlength-1,1,-1
410  psi(ilat,ilev) = psi(ilat,ilev+1)
411    if (vm(ilat,ilev).ne.miss_val) then
412      psi(ilat,ilev) = psi(ilat,ilev) &
413           + vm(ilat,ilev)
414    endif
415 enddo
416enddo
417
418!===============================================================================
419! 2.2.6 Building 2D+time variables
420!===============================================================================
421
422    epy(:,:,itim) = epy2d(:,:)
423    epz(:,:,itim) = epz2d(:,:)
424  divep(:,:,itim) = div2d(:,:)
425ammctem(:,:,itim) = ammc2d(:,:)
426   uzon(:,:,itim) =   ubar(:,:)
427   vtem(:,:,itim) = vtem2d(:,:)
428   wtem(:,:,itim) = wtem2d(:,:)
429 psitem(:,:,itim) = psi(:,:)
430
431
432enddo ! timelength
433!===============================================================================
434!!! END GLOBAL TIME LOOP !!!
435!===============================================================================
436
437print*,"End of computations"
438
439!===============================================================================
440! 3. Create output file
441!===============================================================================
442
443! Create output file
444ierr=NF_CREATE(outfile,NF_CLOBBER,outfid)
445if (ierr.ne.NF_NOERR) then
446  write(*,*)"Error: could not create file ",outfile
447  stop
448endif
449
450!===============================================================================
451! 3.1. Define and write dimensions
452!===============================================================================
453
454call write_dim(outfid,lonlength,latlength,altlength,timelength, &
455    lon,lat,plev,time,lon_dimid,lat_dimid,alt_dimid,time_dimid)
456
457!===============================================================================
458! 3.2. Define and write variables
459!===============================================================================
460
461datashape3d(1)=lat_dimid
462datashape3d(2)=alt_dimid
463datashape3d(3)=time_dimid
464
465call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
466                "epy       ", "EP flux on lat      ","m3 s-2    ",miss_val,&
467                 epy )
468
469call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
470                "epz       ", "EP flux on press    ","m3 s-2    ",miss_val,&
471                 epz )
472
473call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
474                "divep     ", "Div of EP flux      ","m s-2     ",miss_val,&
475                 divep )
476
477call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
478                "ammctem   ", "acc by residual mmc ","m s-2     ",miss_val,&
479                 ammctem )
480
481call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
482                "uzon      ", "Mean zonal wind     ","m s-1     ",miss_val,&
483                 uzon )
484
485call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
486                "vtem      ", "Resid TEM merid wind","m s-1     ",miss_val,&
487                 vtem )
488
489call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
490                "wtem      ", "Resid TEM verti wind","Pa s-1    ",miss_val,&
491                 wtem )
492
493call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
494                "psitem    ", "Resid stream funct  ","kg s-1    ",miss_val,&
495                 psitem )
496
497
498!!!! Close output file
499ierr=NF_CLOSE(outfid)
500if (ierr.ne.NF_NOERR) write(*,*) 'Error, failed to close output file ',outfile
501
502
503end program
Note: See TracBrowser for help on using the repository browser.