source: trunk/LMDZ.TITAN/Tools/energy.F90 @ 842

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

SL: tools for postprocessing (Veznus and Titan); see DOC/documentation/vt-tools.pdf
File size: 14.0 KB
Line 
1program energy
2
3! SL 12/2009:
4! This program reads 4D (lon-lat-alt-time) fields directly from LMD outputs without regrid : histmth
5!
6! it computes:
7! dmass -- 4D -- mass of each cell
8! sek   -- 4D -- specific kinetic energy
9! ek    -- 1D -- integrated kinetic energy
10! sep   -- 4D -- specific potential energy
11! ep    -- 1D -- integrated potential energy
12!
13! Minimal requirements and dependencies:
14! The dataset must include the following data:
15! - surface pressure
16! - atmospheric temperature
17! - zonal and meridional winds
18
19! VERTICAL WIND SPEED IS NEGLECTED IN KINETIC ENERGY
20
21implicit none
22
23include "netcdf.inc" ! NetCDF definitions
24
25character (len=128) :: infile ! input file name (name_P.nc)
26character (len=128) :: outfile ! output file name
27
28character (len=64) :: text ! to store some text
29integer infid ! NetCDF input file ID
30integer outfid ! NetCDF output file ID
31integer lon_dimid,lat_dimid,alt_dimid,time_dimid ! NetCDF dimension IDs
32integer lon_varid,lat_varid,alt_varid,time_varid
33integer              :: datashape1d ! shape of 1D datasets
34integer,dimension(2) :: datashape2d ! shape of 2D datasets
35integer,dimension(3) :: datashape3d ! shape of 3D datasets
36integer,dimension(4) :: datashape4d ! shape of 4D datasets
37
38real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
39real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
40real :: pi
41real,dimension(:),allocatable :: lon ! longitude
42integer lonlength ! # of grid points along longitude
43real,dimension(:),allocatable :: lat ! latitude
44real,dimension(:),allocatable :: latrad ! latitude in radian
45integer latlength ! # of grid points along latitude
46real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
47integer altlength ! # of grid point along presnivs (of input datasets)
48real,dimension(:),allocatable :: time ! time
49integer timelength ! # of points along time
50real,dimension(:,:,:),allocatable :: ps ! surface pressure
51real,dimension(:,:,:,:),allocatable :: temp ! atmospheric temperature
52real,dimension(:,:,:,:),allocatable :: vitu ! zonal wind (in m/s)
53real,dimension(:,:,:,:),allocatable :: vitv ! meridional wind (in m/s)
54
55real,dimension(:,:,:,:),allocatable :: rayon ! distance to center (m)
56real,dimension(:,:,:,:),allocatable :: grav ! gravity field (m s-2)
57real,dimension(:,:,:,:),allocatable :: dmass ! mass in cell (kg)
58real,dimension(:,:,:,:),allocatable :: sek ! specific kinetic energy
59real,dimension(:),allocatable :: ek ! total kinetic energy
60real,dimension(:,:,:,:),allocatable :: sep ! specific potential energy
61real,dimension(:),allocatable :: ep ! total potential energy
62
63integer ierr,ierr1,ierr2 ! NetCDF routines return codes
64integer i,j,ilon,ilat,ilev,itim ! for loops
65
66real :: deltalat,deltalon ! lat and lon intervals in radians
67real,dimension(:,:,:,:),allocatable :: deltap ! pressure thickness of each layer (Pa)
68real :: tmpp ! temporary pressure
69real :: dz ! altitude diff
70real :: signe ! orientation of lon axis for mountain torque computation
71logical :: lmdflag
72
73real :: cpdet
74
75include "planet.h"
76
77!===============================================================================
78! 1. Input parameters
79!===============================================================================
80
81pi = 2.*asin(1.)
82
83write(*,*) ""
84write(*,*) "You are working on the atmosphere of ",planet
85
86!===============================================================================
87! 1.1 Input file
88!===============================================================================
89
90write(*,*) ""
91write(*,*) " Program valid for Venus or Titan LMD, or Venus CAM output files"
92write(*,*) "Enter input file name:"
93
94read(*,'(a128)') infile
95write(*,*) ""
96
97! open input file
98
99ierr = NF_OPEN(infile,NF_NOWRITE,infid)
100if (ierr.ne.NF_NOERR) then
101   write(*,*) 'ERROR: Pb opening file ',trim(infile)
102   stop ""
103endif
104
105!===============================================================================
106! 1.2 Get grids in lon,lat,alt(pressure),time
107!===============================================================================
108
109call get_iddim(infid,lat_varid,latlength,lon_varid,lonlength,&
110                     alt_varid,altlength,time_varid,timelength,lmdflag )
111
112allocate(lon(lonlength))
113ierr=NF_GET_VAR_REAL(infid,lon_varid,lon)
114if (ierr.ne.NF_NOERR) stop "Error: Failed reading longitude"
115if(lon(1).gt.lon(2)) then
116  signe=-1.
117else
118  signe=1.
119endif
120
121allocate(lat(latlength))
122ierr=NF_GET_VAR_REAL(infid,lat_varid,lat)
123if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
124
125allocate(latrad(latlength))
126latrad = lat*pi/180.
127
128! Lat, lon pressure intervals
129deltalat = abs(latrad(2)-latrad(1))
130deltalon = abs(lon(2)-lon(1))*pi/180.
131
132allocate(plev(altlength))
133ierr=NF_GET_VAR_REAL(infid,alt_varid,plev)
134if (ierr.ne.NF_NOERR) stop "Error: Failed reading presnivs (ie pressure levels)"
135if(.not.lmdflag) then
136! in CAM files, pressure in mbar and reversed...
137  call reverselev(altlength,plev)
138  plev=plev*100.  ! convertion to Pa
139endif
140
141allocate(time(timelength))
142ierr=NF_GET_VAR_REAL(infid,time_varid,time)
143if (ierr.ne.NF_NOERR) stop "Error: Failed reading time"
144
145! Time axis IN PLANET DAYS
146
147if(.not.lmdflag) then
148! in CAM files, time in Earth days...
149!   => seconds
150  time=time*86400.
151endif
152time=time/localday
153
154!===============================================================================
155! 1.3 Get output file name
156!===============================================================================
157write(*,*) ""
158!write(*,*) "Enter output file name"
159!read(*,*) outfile
160outfile=infile(1:len_trim(infile)-3)//"_NRG.nc"
161write(*,*) "Output file name is: "//trim(outfile)
162
163
164
165!===============================================================================
166! 2.1 Store needed fields
167!===============================================================================
168
169!===============================================================================
170! 2.1.1 Surface pressure
171!===============================================================================
172allocate(ps(lonlength,latlength,timelength))
173
174text="ps"
175call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
176if (ierr1.ne.NF_NOERR) then
177  write(*,*) "  looking for psol instead... "
178  text="psol"
179  call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
180  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get psol ID"
181endif
182if (ierr2.ne.NF_NOERR) stop "Error: Failed reading surface pressure"
183if((.not.lmdflag).and.(planet.eq."Venus")) call reverse3d(lonlength,latlength,timelength,ps)
184
185!===============================================================================
186! 2.1.2 Atmospheric temperature
187!===============================================================================
188allocate(temp(lonlength,latlength,altlength,timelength))
189
190if(lmdflag) then
191  text="temp"
192else
193  text="T"
194endif
195call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
196if (ierr1.ne.NF_NOERR) then
197  write(*,*) "  looking for t instead... "
198  text="t"
199  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
200  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
201endif
202if (ierr2.ne.NF_NOERR) stop "Error: Failed reading temperature"
203
204if(.not.lmdflag) call reverse4d(lonlength,latlength,altlength,timelength,temp)
205
206!===============================================================================
207! 2.1.3 Winds
208!===============================================================================
209allocate(vitu(lonlength,latlength,altlength,timelength))
210allocate(vitv(lonlength,latlength,altlength,timelength))
211
212! zonal wind vitu (in m/s)
213if(lmdflag) then
214  text="vitu"
215else
216  text="U"
217endif
218call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,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
222if(.not.lmdflag) call reverse4d(lonlength,latlength,altlength,timelength,vitu)
223
224! meridional wind vitv (in m/s)
225if(lmdflag) then
226  text="vitv"
227else
228  text="V"
229endif
230call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
231if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
232if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
233
234if(.not.lmdflag) call reverse4d(lonlength,latlength,altlength,timelength,vitv)
235
236!===============================================================================
237! 2.1.4 Altitude above areoide
238!===============================================================================
239! Only needed if g(z) on Titan...
240
241!allocate(za(lonlength,latlength,altlength,timelength))
242
243!text="zareoid"
244!call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
245!if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
246!if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
247
248!===============================================================================
249! 2.2 Computations
250!===============================================================================
251
252!===============================================================================
253! 2.2.2 Mass in cells
254!===============================================================================
255allocate(rayon(lonlength,latlength,altlength,timelength))
256allocate(grav(lonlength,latlength,altlength,timelength))
257allocate(dmass(lonlength,latlength,altlength,timelength))
258
259do itim=1,timelength
260do ilon=1,lonlength
261 do ilat=1,latlength
262  do ilev=1,altlength
263! Need to be consistent with GCM computations
264!    if (za(ilon,ilat,ilev,itim).ne.miss_val) then
265     rayon(ilon,ilat,ilev,itim) = a0
266!     rayon(ilon,ilat,ilev,itim) = za(ilon,ilat,ilev,itim) + a0
267      grav(ilon,ilat,ilev,itim) = g0*a0*a0 &
268                 /(rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim))
269!    else
270!     rayon(ilon,ilat,ilev,itim) = miss_val
271!      grav(ilon,ilat,ilev,itim) = miss_val
272!    endif
273  enddo
274 enddo
275enddo
276enddo ! timelength
277
278call cellmass(infid,latlength,lonlength,altlength,timelength, &
279              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
280              dmass )
281
282!===============================================================================
283! 2.2.6 Specific energies
284!===============================================================================
285allocate(sek(lonlength,latlength,altlength,timelength))
286allocate(sep(lonlength,latlength,altlength,timelength))
287
288do itim=1,timelength
289
290do ilon=1,lonlength
291 do ilat=1,latlength
292  do ilev=1,altlength
293    if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
294      if ((vitu(ilon,ilat,ilev,itim).lt.miss_val) &
295     .and.(vitv(ilon,ilat,ilev,itim).lt.miss_val)) then
296      sek(ilon,ilat,ilev,itim) = 0.5 * &
297       ( vitu(ilon,ilat,ilev,itim)*vitu(ilon,ilat,ilev,itim) &
298       + vitv(ilon,ilat,ilev,itim)*vitv(ilon,ilat,ilev,itim) )
299       else
300      sek(ilon,ilat,ilev,itim) = miss_val
301       endif
302      if (temp(ilon,ilat,ilev,itim).ne.miss_val) then
303      sep(ilon,ilat,ilev,itim) = temp(ilon,ilat,ilev,itim) &
304                             * cpdet(temp(ilon,ilat,ilev,itim))
305       else
306      sep(ilon,ilat,ilev,itim) = miss_val
307       endif
308    else
309      sek(ilon,ilat,ilev,itim) = miss_val
310      sep(ilon,ilat,ilev,itim) = miss_val
311    endif
312  enddo
313 enddo
314enddo
315
316enddo ! timelength
317
318!===============================================================================
319! 2.2.7 Total energies
320!===============================================================================
321allocate(ek(timelength))
322allocate(ep(timelength))
323
324do itim=1,timelength
325
326ek(itim) = 0.
327ep(itim) = 0.
328do ilon=1,lonlength
329 do ilat=1,latlength
330  do ilev=1,altlength
331    if (sek(ilon,ilat,ilev,itim).ne.miss_val) then
332      ek(itim) = ek(itim) &
333       + sek(ilon,ilat,ilev,itim) * dmass(ilon,ilat,ilev,itim)
334    endif
335    if (sep(ilon,ilat,ilev,itim).ne.miss_val) then
336      ep(itim) = ep(itim) &
337       + sep(ilon,ilat,ilev,itim) * dmass(ilon,ilat,ilev,itim)
338    endif
339  enddo
340 enddo
341enddo
342if (ek(itim).eq.0.) then
343  ek(itim) = miss_val
344  ep(itim) = miss_val
345endif
346
347enddo ! timelength
348
349print*,"End of computations"
350
351!===============================================================================
352! 3. Create output file
353!===============================================================================
354
355! Create output file
356ierr=NF_CREATE(outfile,NF_CLOBBER,outfid)
357if (ierr.ne.NF_NOERR) then
358  write(*,*)"Error: could not create file ",outfile
359  stop
360endif
361
362!===============================================================================
363! 3.1. Define and write dimensions
364!===============================================================================
365
366call write_dim(outfid,lonlength,latlength,altlength,timelength, &
367    lon,lat,plev,time,lon_dimid,lat_dimid,alt_dimid,time_dimid)
368
369!===============================================================================
370! 3.2. Define and write variables
371!===============================================================================
372
373! 1D Variables
374
375datashape1d=time_dimid
376 
377call write_var1d(outfid,datashape1d,timelength,&
378                "ek        ", "Total kinetic energy","J         ",miss_val,&
379                 ek )
380
381call write_var1d(outfid,datashape1d,timelength,&
382                "ep        ", "Total pot energy    ","J         ",miss_val,&
383                 ep )
384
385! 3D variables
386
387datashape3d(1)=lon_dimid
388datashape3d(2)=lat_dimid
389datashape3d(3)=time_dimid
390
391call write_var3d(outfid,datashape3d,lonlength,latlength,timelength,&
392                "ps        ", "Surface pressure    ","Pa        ",miss_val,&
393                 ps )
394
395! 4D variables
396
397datashape4d(1)=lon_dimid
398datashape4d(2)=lat_dimid
399datashape4d(3)=alt_dimid
400datashape4d(4)=time_dimid
401
402call write_var4d(outfid,datashape4d,lonlength,latlength,altlength,timelength,&
403                "dmass     ", "Mass                ","kg        ",miss_val,&
404                 dmass )
405
406call write_var4d(outfid,datashape4d,lonlength,latlength,altlength,timelength,&
407                "sek       ", "Specific kin energy ","J/kg      ",miss_val,&
408                 sek )
409
410call write_var4d(outfid,datashape4d,lonlength,latlength,altlength,timelength,&
411                "sep       ", "Specific pot energy ","J/kg      ",miss_val,&
412                 sep )
413
414
415!!!! Close output file
416ierr=NF_CLOSE(outfid)
417if (ierr.ne.NF_NOERR) then
418  write(*,*) 'Error, failed to close output file ',outfile
419endif
420
421
422end program
Note: See TracBrowser for help on using the repository browser.