source: trunk/LMDZ.TITAN.old/Tools/stability.F90 @ 3026

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

SL: corrections in Titan and Venus tools
File size: 16.4 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 ! 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 :: coslat ! cos of latitude
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.)
93miss_val = miss_val_def
94
95write(*,*) ""
96write(*,*) "You are working on the atmosphere of ",planet
97
98!===============================================================================
99! 1.1 Input file
100!===============================================================================
101
102write(*,*) ""
103write(*,*) "Program valid for files with pressure axis (*_P.nc)"
104write(*,*) "Enter input file name:"
105
106read(*,'(a128)') infile
107write(*,*) ""
108
109! open input file
110
111ierr = NF_OPEN(infile,NF_NOWRITE,infid)
112if (ierr.ne.NF_NOERR) then
113   write(*,*) 'ERROR: Pb opening file ',trim(infile)
114   stop ""
115endif
116
117!===============================================================================
118! 1.2 Get grids in lon,lat,alt(pressure),time
119!===============================================================================
120
121call get_iddim(infid,lat_varid,latlength,lon_varid,lonlength,&
122                     alt_varid,altlength,time_varid,timelength,lmdflag )
123
124allocate(lon(lonlength))
125ierr=NF_GET_VAR_REAL(infid,lon_varid,lon)
126if (ierr.ne.NF_NOERR) stop "Error: Failed reading longitude"
127
128allocate(lat(latlength))
129ierr=NF_GET_VAR_REAL(infid,lat_varid,lat)
130if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
131
132allocate(coslat(latlength))
133! Beware of rounding problems at poles...
134coslat(:) = max(0.,cos(lat(:)*pi/180.))
135
136! Lat, lon pressure intervals
137deltalat = abs(lat(2)-lat(1))*pi/180.
138deltalon = abs(lon(2)-lon(1))*pi/180.
139
140allocate(plev(altlength))
141ierr=NF_GET_VAR_REAL(infid,alt_varid,plev)
142if (ierr.ne.NF_NOERR) stop "Error: Failed reading altitude (ie pressure levels)"
143
144allocate(time(timelength))
145ierr=NF_GET_VAR_REAL(infid,time_varid,time)
146if (ierr.ne.NF_NOERR) stop "Error: Failed reading time"
147
148!===============================================================================
149! 1.3 Get output file name
150!===============================================================================
151write(*,*) ""
152!write(*,*) "Enter output file name"
153!read(*,*) outfile
154outfile=infile(1:len_trim(infile)-3)//"_STA.nc"
155write(*,*) "Output file name is: "//trim(outfile)
156
157
158
159!===============================================================================
160! 2.1 Store needed fields
161!===============================================================================
162
163!===============================================================================
164! 2.1.1 Surface pressure
165!===============================================================================
166allocate(ps(lonlength,latlength,timelength))
167
168text="ps"
169call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
170if (ierr1.ne.NF_NOERR) then
171  write(*,*) "  looking for psol instead... "
172  text="psol"
173  call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
174  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get psol ID"
175endif
176if (ierr2.ne.NF_NOERR) stop "Error: Failed reading surface pressure"
177
178!===============================================================================
179! 2.1.2 Atmospheric temperature
180!===============================================================================
181allocate(temp(lonlength,latlength,altlength,timelength))
182
183text="temp"
184call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
185if (ierr1.ne.NF_NOERR) then
186  write(*,*) "  looking for t instead... "
187  text="t"
188  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
189  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
190endif
191if (ierr2.ne.NF_NOERR) stop "Error: Failed reading temperature"
192
193!===============================================================================
194! 2.1.3 Winds
195!===============================================================================
196allocate(vitu(lonlength,latlength,altlength,timelength))
197allocate(vitv(lonlength,latlength,altlength,timelength))
198
199! zonal wind vitu (in m/s)
200text="vitu"
201call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
202if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
203if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
204
205! meridional wind vitv (in m/s)
206text="vitv"
207call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
208if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
209if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
210
211!===============================================================================
212! 2.1.4 Altitude above areoide
213!===============================================================================
214! Only needed if g(z) on Titan...
215
216! allocate(za(lonlength,latlength,altlength,timelength))
217
218! text="zareoid"
219! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
220! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
221! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
222
223!===============================================================================
224!!! Allocations before timeloop
225!===============================================================================
226
227! latlength correspond a jjm+1
228! mais lonlength correspond a iim
229! pour boucler en longitude, on a besoin du point iim+1 (= 1)
230
231allocate(rayon(lonlength+1,latlength,altlength,timelength))
232allocate(grav(lonlength+1,latlength,altlength,timelength))
233allocate(dmass(lonlength+1,latlength,altlength,timelength))
234
235allocate(t3d(lonlength+1,latlength,altlength))
236allocate(u3d(lonlength+1,latlength,altlength))
237allocate(v3d(lonlength+1,latlength,altlength))
238
239allocate(t2d(latlength,altlength))
240allocate(u2d(latlength,altlength))
241allocate(v2d(latlength,altlength))
242allocate(dtsdp(latlength,altlength))
243allocate(dusdp(latlength,altlength))
244allocate(dvsdp(latlength,altlength))
245
246allocate(stab(lonlength,latlength,altlength,timelength))
247allocate(Ri(lonlength,latlength,altlength,timelength))
248allocate(deqc(latlength,altlength,timelength))
249
250allocate(stab3d(lonlength+1,latlength,altlength))
251allocate(Ri3d(lonlength+1,latlength,altlength))
252allocate(deqc2d(latlength,altlength))
253
254!===============================================================================
255! 2.2.2 Mass in cells
256!===============================================================================
257
258do itim=1,timelength
259do ilon=1,lonlength
260 do ilat=1,latlength
261  do ilev=1,altlength
262! Need to be consistent with GCM computations
263!    if (za(ilon,ilat,ilev,itim).ne.miss_val) then
264     rayon(ilon,ilat,ilev,itim) = a0
265!     rayon(ilon,ilat,ilev,itim) = za(ilon,ilat,ilev,itim) + a0
266      grav(ilon,ilat,ilev,itim) = g0*a0*a0 &
267                 /(rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim))
268!    else
269!     rayon(ilon,ilat,ilev,itim) = miss_val
270!      grav(ilon,ilat,ilev,itim) = miss_val
271!    endif
272  enddo
273 enddo
274enddo
275enddo ! timelength
276
277rayon(lonlength+1,:,:,:) = rayon(1,:,:,:)
278 grav(lonlength+1,:,:,:) =  grav(1,:,:,:)
279
280call cellmass(infid,latlength,lonlength+1,altlength,timelength,lmdflag, &
281              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
282              dmass )
283
284!===============================================================================
285!!! GLOBAL TIME LOOP !!!
286!===============================================================================
287do itim=1,timelength
288
289!===============================================================================
290! 2.2 Computations
291!===============================================================================
292
293!===============================================================================
294! 2.2.3 Init of 3D variables
295!===============================================================================
296
297do ilon=1,lonlength
298 do ilat=1,latlength
299  do ilev=1,altlength
300      t3d(ilon,ilat,ilev) = temp(ilon,ilat,ilev,itim)
301      u3d(ilon,ilat,ilev) = vitu(ilon,ilat,ilev,itim)
302      v3d(ilon,ilat,ilev) = vitv(ilon,ilat,ilev,itim)
303  enddo
304 enddo
305enddo
306
307 t3d(lonlength+1,:,:) =  t3d(1,:,:)
308 u3d(lonlength+1,:,:) =  u3d(1,:,:)
309 v3d(lonlength+1,:,:) =  v3d(1,:,:)
310
311!===============================================================================
312! 2.2.4 Stability
313!===============================================================================
314
315do ilon=1,lonlength+1
316 t2d(:,:) =  t3d(ilon,:,:)
317 call dx_dp(latlength,altlength,miss_val,plev,t2d,dtsdp)
318 do ilat=1,latlength
319  do ilev=1,altlength
320    if ((grav(ilon,ilat,ilev,itim).ne.miss_val).and. &
321        ( t3d(ilon,ilat,ilev).ne.miss_val) ) then
322      stab3d(ilon,ilat,ilev) = grav(ilon,ilat,ilev,itim)* &
323                (  1./cpdet(t2d(ilat,ilev))          &
324                 - plev(ilev)*dtsdp(ilat,ilev)/(R0*t2d(ilat,ilev)) )
325      stab3d(ilon,ilat,ilev) = stab3d(ilon,ilat,ilev)*1000. ! passage en K/km
326    else
327      stab3d(ilon,ilat,ilev) = miss_val
328    endif
329  enddo
330 enddo
331enddo
332
333!===============================================================================
334! 2.2.5 Richardson number
335!===============================================================================
336
337do ilon=1,lonlength+1
338 u2d(:,:) =  u3d(ilon,:,:)
339 v2d(:,:) =  v3d(ilon,:,:)
340 call dx_dp(latlength,altlength,miss_val,plev,u2d,dusdp)
341 call dx_dp(latlength,altlength,miss_val,plev,v2d,dvsdp)
342 do ilat=1,latlength
343  do ilev=1,altlength
344    if ((grav(ilon,ilat,ilev,itim).ne.miss_val).and. &
345        ( u3d(ilon,ilat,ilev).ne.miss_val).and. &
346        ( v3d(ilon,ilat,ilev).ne.miss_val).and. &
347        ( t3d(ilon,ilat,ilev).ne.miss_val) ) then
348      Ri3d(ilon,ilat,ilev) =  & ! attention, transfo a cause de du/dp au lieu de du/dz
349          stab3d(ilon,ilat,ilev)*t3d(ilon,ilat,ilev)*R0*R0  &
350      / (grav(ilon,ilat,ilev,itim)*plev(ilev)*plev(ilev))        &
351      / (dusdp(ilat,ilev)*dusdp(ilat,ilev)+dvsdp(ilat,ilev)*dvsdp(ilat,ilev))
352    else
353      Ri3d(ilon,ilat,ilev) = miss_val
354    endif
355  enddo
356 enddo
357enddo
358
359!===============================================================================
360! 2.2.6 Distance to cyclostrophic equilibrium
361!===============================================================================
362
363call moyzon(lonlength,latlength,altlength,miss_val,u3d,u2d)
364call moyzon(lonlength,latlength,altlength,miss_val,t3d,t2d)
365call dx_dp(latlength,altlength,miss_val,plev,u2d,dusdp)
366
367do ilat=2,latlength-1
368   if (tan(lat(ilat)*pi/180.).ne.0.) then
369     fac1 = R0/tan(lat(ilat)*pi/180.)
370   else
371     fac1 = miss_val
372   endif
373  do ilev=1,altlength
374   if ((dusdp(ilat,ilev).ne.miss_val).and. &
375       (  u2d(ilat,ilev).ne.miss_val).and. &
376       (            fac1.ne.miss_val).and. &
377       (  t2d(ilat,ilev).ne.miss_val) ) then
378    ecden = dusdp(ilat,ilev)*(2.*u2d(ilat,ilev)*plev(ilev))
379    ecnum = ((t2d(ilat+1,ilev)-t2d(ilat-1,ilev))/(2.*deltalat)*fac1-ecden)*100.
380    deqc2d(ilat,ilev) = ecnum/ecden
381   else
382    deqc2d(ilat,ilev) = miss_val
383   endif
384  enddo
385enddo
386do ilev=1,altlength
387    deqc2d(1,ilev)         = miss_val
388    deqc2d(latlength,ilev) = miss_val
389enddo
390
391!===============================================================================
392! 2.2.7 Building 3D+time variables
393!===============================================================================
394
395    deqc(:,:,itim)   = deqc2d(:,:)
396    stab(:,:,:,itim) = stab3d(1:lonlength,:,:)
397      Ri(:,:,:,itim) =   Ri3d(1:lonlength,:,:)
398
399
400enddo ! timelength
401!===============================================================================
402!!! END GLOBAL TIME LOOP !!!
403!===============================================================================
404
405print*,"End of computations"
406
407!===============================================================================
408! 3. Create output file
409!===============================================================================
410
411! Create output file
412ierr=NF_CREATE(outfile,NF_CLOBBER,outfid)
413if (ierr.ne.NF_NOERR) then
414  write(*,*)"Error: could not create file ",outfile
415  stop
416endif
417
418!===============================================================================
419! 3.1. Define and write dimensions
420!===============================================================================
421
422call write_dim(outfid,lonlength,latlength,altlength,timelength, &
423    lon,lat,plev,time,lon_dimid,lat_dimid,alt_dimid,time_dimid)
424
425!===============================================================================
426! 3.2. Define and write variables
427!===============================================================================
428
429! 3D Variables
430
431datashape3d(1)=lat_dimid
432datashape3d(2)=alt_dimid
433datashape3d(3)=time_dimid
434
435call write_var3d(outfid,datashape3d,latlength,altlength,timelength,&
436                "deqc      ", "Distance to cyclo eq","per cent  ",miss_val,&
437                 deqc )
438
439! 4D Variables
440
441datashape4d(1)=lon_dimid
442datashape4d(2)=lat_dimid
443datashape4d(3)=alt_dimid
444datashape4d(4)=time_dimid
445
446call write_var4d(outfid,datashape4d,lonlength,latlength,altlength,timelength,&
447                "stab      ", "Stability           ","K/km      ",miss_val,&
448                 stab )
449
450call write_var4d(outfid,datashape4d,lonlength,latlength,altlength,timelength,&
451                "Ri        ", "Richardson number   ","          ",miss_val,&
452                 Ri )
453
454
455!!!! Close output file
456ierr=NF_CLOSE(outfid)
457if (ierr.ne.NF_NOERR) then
458  write(*,*) 'Error, failed to close output file ',outfile
459endif
460
461
462end program
Note: See TracBrowser for help on using the repository browser.