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

Timestamp:

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

Author:

slebonnois

Message:

SL: corrections in Titan and Venus tools

File:

• trunk/LMDZ.TITAN/Tools/tem.F90 (modified) (12 diffs)

trunk/LMDZ.TITAN/Tools/tem.F90

@@ -37 +37 @@
 integer,dimension(3) :: datashape3d ! shape of 3D 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)
@@ -107 +107 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -145 +146 @@
 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.
 
@@ -196 +198 @@
 
 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
@@ -214 +216 @@
 ! 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"
@@ -220 +222 @@
 ! 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"
@@ -226 +228 @@
 ! vertical wind vitw (in Pa/s)
 text="vitw"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitw ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading vertical wind"
@@ -238 +240 @@
 
 ! 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"
@@ -312 +314 @@
  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 )
 
@@ -359 +361 @@
 call moyzon(lonlength,latlength,altlength,miss_val,u3d,ubar)
 
-call epflux(lonlength+1,latlength,altlength,miss_val,latrad,rbar &
+call epflux(lonlength+1,latlength,altlength,miss_val,lat,rbar &
             ,teta,u3d,v3d,w3d,plev &
             ,epy2d,epz2d,div2d,vtem2d,wtem2d,ammc2d &
@@ -376 +378 @@
 ! rsurg: r*dp/g = dm/(r cos(lat) dlon dlat) !!!
      rsurg(ilon,ilat,ilev) = dmass(ilon,ilat,ilev,itim) &
-          / (r3d(ilon,ilat,ilev)*cos(latrad(ilat))*deltalon*deltalat)
+          / (r3d(ilon,ilat,ilev)*coslat(ilat)*deltalon*deltalat)
     else
      rsurg(ilon,ilat,ilev) = miss_val
@@ -391 +393 @@
         (rsurgbar(ilat,ilev).ne.miss_val) ) then
       vm(ilat,ilev) = vtem2d(ilat,ilev) &
-            * 2.*pi*rsurgbar(ilat,ilev)*cos(latrad(ilat))
+            * 2.*pi*rsurgbar(ilat,ilev)*coslat(ilat)
     else
       vm(ilat,ilev) = miss_val