source: trunk/LMDZ.TITAN/Tools/stability.F90 @ 819

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

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