Changeset 1454 for trunk/LMDZ.TITAN/Tools/stability.F90

Timestamp:

Jun 16, 2015, 7:01:30 PM (10 years ago)

Author:

slebonnois

Message:

SL: corrections in Titan and Venus tools

File:

• trunk/LMDZ.TITAN/Tools/stability.F90 (modified) (9 diffs)

trunk/LMDZ.TITAN/Tools/stability.F90

@@ -31 +31 @@
 integer,dimension(4) :: datashape4d ! shape of 4D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
+real,dimension(:),allocatable :: coslat ! cos of latitude
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -91 +91 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -129 +130 @@
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
 
-allocate(latrad(latlength))
-latrad = lat*pi/180.
+allocate(coslat(latlength))
+! Beware of rounding problems at poles...
+coslat(:) = max(0.,cos(lat(:)*pi/180.))
 
 ! Lat, lon pressure intervals
-deltalat = abs(latrad(2)-latrad(1))
+deltalat = abs(lat(2)-lat(1))*pi/180.
 deltalon = abs(lon(2)-lon(1))*pi/180.
 
@@ -180 +182 @@
 
 text="temp"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "  looking for t instead... "
   text="t"
-  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
   if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
 endif
@@ -197 +199 @@
 ! zonal wind vitu (in m/s)
 text="vitu"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
@@ -203 +205 @@
 ! meridional wind vitv (in m/s)
 text="vitv"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
@@ -215 +217 @@
 
 ! text="zareoid"
-! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
 ! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
 ! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
@@ -276 +278 @@
  grav(lonlength+1,:,:,:) =  grav(1,:,:,:)
 
-call cellmass(infid,latlength,lonlength+1,altlength,timelength, &
-              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+call cellmass(infid,latlength,lonlength+1,altlength,timelength,lmdflag, &
+              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
               dmass )
 
@@ -364 +366 @@
 
 do ilat=2,latlength-1
-   if (tan(latrad(ilat)).ne.0.) then
-     fac1 = R0/tan(latrad(ilat))
+   if (tan(lat(ilat)*pi/180.).ne.0.) then
+     fac1 = R0/tan(lat(ilat)*pi/180.)
    else
      fac1 = miss_val