source: trunk/LMDZ.TITAN.old/Tools/energy.F90 @ 3552

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

SL: corrections in Titan and Venus tools
File size: 14.1 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 ! 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 :: coslat ! cos of latitude
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.)
82miss_val = miss_val_def
83
84write(*,*) ""
85write(*,*) "You are working on the atmosphere of ",planet
86
87!===============================================================================
88! 1.1 Input file
89!===============================================================================
90
91write(*,*) ""
92write(*,*) " Program valid for Venus or Titan LMD, or Venus CAM output files"
93write(*,*) "Enter input file name:"
94
95read(*,'(a128)') infile
96write(*,*) ""
97
98! open input file
99
100ierr = NF_OPEN(infile,NF_NOWRITE,infid)
101if (ierr.ne.NF_NOERR) then
102   write(*,*) 'ERROR: Pb opening file ',trim(infile)
103   stop ""
104endif
105
106!===============================================================================
107! 1.2 Get grids in lon,lat,alt(pressure),time
108!===============================================================================
109
110call get_iddim(infid,lat_varid,latlength,lon_varid,lonlength,&
111                     alt_varid,altlength,time_varid,timelength,lmdflag )
112
113allocate(lon(lonlength))
114ierr=NF_GET_VAR_REAL(infid,lon_varid,lon)
115if (ierr.ne.NF_NOERR) stop "Error: Failed reading longitude"
116if(lon(1).gt.lon(2)) then
117  signe=-1.
118else
119  signe=1.
120endif
121
122allocate(lat(latlength))
123ierr=NF_GET_VAR_REAL(infid,lat_varid,lat)
124if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
125
126allocate(coslat(latlength))
127! Beware of rounding problems at poles...
128coslat(:) = max(0.,cos(lat(:)*pi/180.))
129
130! Lat, lon pressure intervals
131deltalat = abs(lat(2)-lat(1))*pi/180.
132deltalon = abs(lon(2)-lon(1))*pi/180.
133
134allocate(plev(altlength))
135ierr=NF_GET_VAR_REAL(infid,alt_varid,plev)
136if (ierr.ne.NF_NOERR) stop "Error: Failed reading presnivs (ie pressure levels)"
137if(.not.lmdflag) then
138! in CAM files, pressure in mbar and reversed...
139  call reverselev(altlength,plev)
140  plev=plev*100.  ! convertion to Pa
141endif
142
143allocate(time(timelength))
144ierr=NF_GET_VAR_REAL(infid,time_varid,time)
145if (ierr.ne.NF_NOERR) stop "Error: Failed reading time"
146
147! Time axis IN PLANET DAYS
148
149if(.not.lmdflag) then
150! in CAM files, time in Earth days...
151!   => seconds
152  time=time*86400.
153endif
154time=time/localday
155
156!===============================================================================
157! 1.3 Get output file name
158!===============================================================================
159write(*,*) ""
160!write(*,*) "Enter output file name"
161!read(*,*) outfile
162outfile=infile(1:len_trim(infile)-3)//"_NRG.nc"
163write(*,*) "Output file name is: "//trim(outfile)
164
165
166
167!===============================================================================
168! 2.1 Store needed fields
169!===============================================================================
170
171!===============================================================================
172! 2.1.1 Surface pressure
173!===============================================================================
174allocate(ps(lonlength,latlength,timelength))
175
176text="ps"
177call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
178if (ierr1.ne.NF_NOERR) then
179  write(*,*) "  looking for psol instead... "
180  text="psol"
181  call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
182  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get psol ID"
183endif
184if (ierr2.ne.NF_NOERR) stop "Error: Failed reading surface pressure"
185if((.not.lmdflag).and.(planet.eq."Venus")) call reverse3d(lonlength,latlength,timelength,ps)
186
187!===============================================================================
188! 2.1.2 Atmospheric temperature
189!===============================================================================
190allocate(temp(lonlength,latlength,altlength,timelength))
191
192if(lmdflag) then
193  text="temp"
194else
195  text="T"
196endif
197call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
198if (ierr1.ne.NF_NOERR) then
199  write(*,*) "  looking for t instead... "
200  text="t"
201  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
202  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
203endif
204if (ierr2.ne.NF_NOERR) stop "Error: Failed reading temperature"
205
206if(.not.lmdflag) call reverse4d(lonlength,latlength,altlength,timelength,temp)
207
208!===============================================================================
209! 2.1.3 Winds
210!===============================================================================
211allocate(vitu(lonlength,latlength,altlength,timelength))
212allocate(vitv(lonlength,latlength,altlength,timelength))
213
214! zonal wind vitu (in m/s)
215if(lmdflag) then
216  text="vitu"
217else
218  text="U"
219endif
220call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
221if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
222if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
223
224if(.not.lmdflag) call reverse4d(lonlength,latlength,altlength,timelength,vitu)
225
226! meridional wind vitv (in m/s)
227if(lmdflag) then
228  text="vitv"
229else
230  text="V"
231endif
232call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
233if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
234if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
235
236if(.not.lmdflag) call reverse4d(lonlength,latlength,altlength,timelength,vitv)
237
238!===============================================================================
239! 2.1.4 Altitude above areoide
240!===============================================================================
241! Only needed if g(z) on Titan...
242
243!allocate(za(lonlength,latlength,altlength,timelength))
244
245!text="zareoid"
246!call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
247!if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
248!if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
249
250!===============================================================================
251! 2.2 Computations
252!===============================================================================
253
254!===============================================================================
255! 2.2.2 Mass in cells
256!===============================================================================
257allocate(rayon(lonlength,latlength,altlength,timelength))
258allocate(grav(lonlength,latlength,altlength,timelength))
259allocate(dmass(lonlength,latlength,altlength,timelength))
260
261do itim=1,timelength
262do ilon=1,lonlength
263 do ilat=1,latlength
264  do ilev=1,altlength
265! Need to be consistent with GCM computations
266!    if (za(ilon,ilat,ilev,itim).ne.miss_val) then
267     rayon(ilon,ilat,ilev,itim) = a0
268!     rayon(ilon,ilat,ilev,itim) = za(ilon,ilat,ilev,itim) + a0
269      grav(ilon,ilat,ilev,itim) = g0*a0*a0 &
270                 /(rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim))
271!    else
272!     rayon(ilon,ilat,ilev,itim) = miss_val
273!      grav(ilon,ilat,ilev,itim) = miss_val
274!    endif
275  enddo
276 enddo
277enddo
278enddo ! timelength
279
280call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
281              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
282              dmass )
283
284!===============================================================================
285! 2.2.6 Specific energies
286!===============================================================================
287allocate(sek(lonlength,latlength,altlength,timelength))
288allocate(sep(lonlength,latlength,altlength,timelength))
289
290do itim=1,timelength
291
292do ilon=1,lonlength
293 do ilat=1,latlength
294  do ilev=1,altlength
295    if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
296      if ((vitu(ilon,ilat,ilev,itim).lt.miss_val) &
297     .and.(vitv(ilon,ilat,ilev,itim).lt.miss_val)) then
298      sek(ilon,ilat,ilev,itim) = 0.5 * &
299       ( vitu(ilon,ilat,ilev,itim)*vitu(ilon,ilat,ilev,itim) &
300       + vitv(ilon,ilat,ilev,itim)*vitv(ilon,ilat,ilev,itim) )
301       else
302      sek(ilon,ilat,ilev,itim) = miss_val
303       endif
304      if (temp(ilon,ilat,ilev,itim).ne.miss_val) then
305      sep(ilon,ilat,ilev,itim) = temp(ilon,ilat,ilev,itim) &
306                             * cpdet(temp(ilon,ilat,ilev,itim))
307       else
308      sep(ilon,ilat,ilev,itim) = miss_val
309       endif
310    else
311      sek(ilon,ilat,ilev,itim) = miss_val
312      sep(ilon,ilat,ilev,itim) = miss_val
313    endif
314  enddo
315 enddo
316enddo
317
318enddo ! timelength
319
320!===============================================================================
321! 2.2.7 Total energies
322!===============================================================================
323allocate(ek(timelength))
324allocate(ep(timelength))
325
326do itim=1,timelength
327
328ek(itim) = 0.
329ep(itim) = 0.
330do ilon=1,lonlength
331 do ilat=1,latlength
332  do ilev=1,altlength
333    if (sek(ilon,ilat,ilev,itim).ne.miss_val) then
334      ek(itim) = ek(itim) &
335       + sek(ilon,ilat,ilev,itim) * dmass(ilon,ilat,ilev,itim)
336    endif
337    if (sep(ilon,ilat,ilev,itim).ne.miss_val) then
338      ep(itim) = ep(itim) &
339       + sep(ilon,ilat,ilev,itim) * dmass(ilon,ilat,ilev,itim)
340    endif
341  enddo
342 enddo
343enddo
344if (ek(itim).eq.0.) then
345  ek(itim) = miss_val
346  ep(itim) = miss_val
347endif
348
349enddo ! timelength
350
351print*,"End of computations"
352
353!===============================================================================
354! 3. Create output file
355!===============================================================================
356
357! Create output file
358ierr=NF_CREATE(outfile,NF_CLOBBER,outfid)
359if (ierr.ne.NF_NOERR) then
360  write(*,*)"Error: could not create file ",outfile
361  stop
362endif
363
364!===============================================================================
365! 3.1. Define and write dimensions
366!===============================================================================
367
368call write_dim(outfid,lonlength,latlength,altlength,timelength, &
369    lon,lat,plev,time,lon_dimid,lat_dimid,alt_dimid,time_dimid)
370
371!===============================================================================
372! 3.2. Define and write variables
373!===============================================================================
374
375! 1D Variables
376
377datashape1d=time_dimid
378 
379call write_var1d(outfid,datashape1d,timelength,&
380                "ek        ", "Total kinetic energy","J         ",miss_val,&
381                 ek )
382
383call write_var1d(outfid,datashape1d,timelength,&
384                "ep        ", "Total pot energy    ","J         ",miss_val,&
385                 ep )
386
387! 3D variables
388
389datashape3d(1)=lon_dimid
390datashape3d(2)=lat_dimid
391datashape3d(3)=time_dimid
392
393call write_var3d(outfid,datashape3d,lonlength,latlength,timelength,&
394                "ps        ", "Surface pressure    ","Pa        ",miss_val,&
395                 ps )
396
397! 4D variables
398
399datashape4d(1)=lon_dimid
400datashape4d(2)=lat_dimid
401datashape4d(3)=alt_dimid
402datashape4d(4)=time_dimid
403
404call write_var4d(outfid,datashape4d,lonlength,latlength,altlength,timelength,&
405                "dmass     ", "Mass                ","kg        ",miss_val,&
406                 dmass )
407
408call write_var4d(outfid,datashape4d,lonlength,latlength,altlength,timelength,&
409                "sek       ", "Specific kin energy ","J/kg      ",miss_val,&
410                 sek )
411
412call write_var4d(outfid,datashape4d,lonlength,latlength,altlength,timelength,&
413                "sep       ", "Specific pot energy ","J/kg      ",miss_val,&
414                 sep )
415
416
417!!!! Close output file
418ierr=NF_CLOSE(outfid)
419if (ierr.ne.NF_NOERR) then
420  write(*,*) 'Error, failed to close output file ',outfile
421endif
422
423
424end program
Note: See TracBrowser for help on using the repository browser.