Changeset 1454 for trunk/LMDZ.VENUS/Tools/angmom.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.VENUS/Tools/angmom.F90 (modified) (18 diffs)

trunk/LMDZ.VENUS/Tools/angmom.F90

@@ -47 +47 @@
 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)
@@ -117 +117 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -160 +161 @@
 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.
 
@@ -281 +283 @@
 endif
 
-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 U ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
@@ -295 +297 @@
 
 !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"
@@ -309 +311 @@
   text="DUVDF"
 endif
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duvdf,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duvdf,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get duvdf ID"
@@ -357 +359 @@
 
 text="dugwo"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwo,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwo,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get dugwo ID"
@@ -367 +369 @@
 
 text="dugwno"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwno,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwno,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get dugwno ID"
@@ -390 +392 @@
 
 text="duajs"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duajs,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duajs,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get duajs ID"
@@ -414 +416 @@
   text="DUDYN"
 endif
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dudyn,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dudyn,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get dudyn ID"
@@ -486 +488 @@
 enddo ! timelength
 
-call cellmass(infid,latlength,lonlength,altlength,timelength, &
-              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
+              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
               dmass )
 
@@ -504 +506 @@
       osam(ilon,ilat,ilev,itim) = &
          rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim) &
-       * cos(latrad(ilat))*cos(latrad(ilat)) &
+       * coslat(ilat)*coslat(ilat) &
        * omega
       rsam(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)
     else
       osam(ilon,ilat,ilev,itim) = miss_val
@@ -569 +571 @@
        tmpp = (ps(ilon+1,ilat,itim)  +ps(ilon,ilat,itim))/2.
       endif
-      tmou(itim) = tmou(itim) + a0*a0* deltalat*cos(latrad(ilat)) &
+      tmou(itim) = tmou(itim) + a0*a0* deltalat*coslat(ilat) &
        * signe*dz*tmpp &
        / hadley
@@ -582 +584 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tdyn(itim) = tdyn(itim) + dudyn(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
        * dmass(ilon,ilat,ilev,itim) &
        / hadley
@@ -599 +601 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tajs(itim) = tajs(itim) + duajs(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
        * dmass(ilon,ilat,ilev,itim) &
        / hadley
@@ -616 +618 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tbls(itim) = tbls(itim) + duvdf(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
        * dmass(ilon,ilat,ilev,itim) &
        / hadley
@@ -633 +635 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tgwo(itim) = tgwo(itim) + dugwo(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
        * dmass(ilon,ilat,ilev,itim) &
        / hadley
@@ -650 +652 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tgwno(itim) = tgwno(itim) + dugwno(ilon,ilat,ilev,itim) &
-         * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))       &
+         * rayon(ilon,ilat,ilev,itim)*coslat(ilat)       &
          * dmass(ilon,ilat,ilev,itim)  &
        / hadley