Changeset 1454


Ignore:
Timestamp:
Jun 16, 2015, 7:01:30 PM (10 years ago)
Author:
slebonnois
Message:

SL: corrections in Titan and Venus tools

Location:
trunk
Files:
31 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/Tools/angmom.F90

    r1356 r1454  
    4747integer,dimension(4) :: datashape4d ! shape of 4D datasets
    4848
    49 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    50 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
     49real :: miss_val ! special "missing value" to specify missing data
     50real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
    5151real :: pi
    5252real,dimension(:),allocatable :: lon ! longitude
    5353integer lonlength ! # of grid points along longitude
    5454real,dimension(:),allocatable :: lat ! latitude
    55 real,dimension(:),allocatable :: latrad ! latitude in radian
     55real,dimension(:),allocatable :: coslat ! cos of latitude
    5656integer latlength ! # of grid points along latitude
    5757real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
    … …  
    117117
    118118pi = 2.*asin(1.)
     119miss_val = miss_val_def
    119120
    120121write(*,*) ""
    … …  
    160161if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
    161162
    162 allocate(latrad(latlength))
    163 latrad = lat*pi/180.
     163allocate(coslat(latlength))
     164! Beware of rounding problems at poles...
     165coslat(:) = max(0.,cos(lat(:)*pi/180.))
    164166
    165167! Lat, lon pressure intervals
    166 deltalat = abs(latrad(2)-latrad(1))
     168deltalat = abs(lat(2)-lat(1))*pi/180.
    167169deltalon = abs(lon(2)-lon(1))*pi/180.
    168170
    … …  
    281283endif
    282284
    283 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
     285call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
    284286if (ierr1.ne.NF_NOERR) stop "Error: Failed to get U ID"
    285287if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
    … …  
    295297
    296298!text="zareoid"
    297 !call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
     299!call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
    298300!if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
    299301!if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
    … …  
    309311  text="DUVDF"
    310312endif
    311 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duvdf,ierr1,ierr2)
     313call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duvdf,miss_val,ierr1,ierr2)
    312314if (ierr1.ne.NF_NOERR) then
    313315  write(*,*) "Failed to get duvdf ID"
    … …  
    357359
    358360text="dugwo"
    359 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwo,ierr1,ierr2)
     361call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwo,miss_val,ierr1,ierr2)
    360362if (ierr1.ne.NF_NOERR) then
    361363  write(*,*) "Failed to get dugwo ID"
    … …  
    367369
    368370text="dugwno"
    369 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwno,ierr1,ierr2)
     371call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwno,miss_val,ierr1,ierr2)
    370372if (ierr1.ne.NF_NOERR) then
    371373  write(*,*) "Failed to get dugwno ID"
    … …  
    390392
    391393text="duajs"
    392 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duajs,ierr1,ierr2)
     394call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duajs,miss_val,ierr1,ierr2)
    393395if (ierr1.ne.NF_NOERR) then
    394396  write(*,*) "Failed to get duajs ID"
    … …  
    414416  text="DUDYN"
    415417endif
    416 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dudyn,ierr1,ierr2)
     418call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dudyn,miss_val,ierr1,ierr2)
    417419if (ierr1.ne.NF_NOERR) then
    418420  write(*,*) "Failed to get dudyn ID"
    … …  
    486488enddo ! timelength
    487489
    488 call cellmass(infid,latlength,lonlength,altlength,timelength, &
    489               lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
     490call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
     491              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
    490492              dmass )
    491493
    … …  
    504506      osam(ilon,ilat,ilev,itim) = &
    505507         rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim) &
    506        * cos(latrad(ilat))*cos(latrad(ilat)) &
     508       * coslat(ilat)*coslat(ilat) &
    507509       * omega
    508510      rsam(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim) &
    509        * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))
     511       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)
    510512    else
    511513      osam(ilon,ilat,ilev,itim) = miss_val
    … …  
    569571       tmpp = (ps(ilon+1,ilat,itim)  +ps(ilon,ilat,itim))/2.
    570572      endif
    571       tmou(itim) = tmou(itim) + a0*a0* deltalat*cos(latrad(ilat)) &
     573      tmou(itim) = tmou(itim) + a0*a0* deltalat*coslat(ilat) &
    572574       * signe*dz*tmpp &
    573575       / hadley
    … …  
    582584    if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
    583585      tdyn(itim) = tdyn(itim) + dudyn(ilon,ilat,ilev,itim) &
    584        * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
     586       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
    585587       * dmass(ilon,ilat,ilev,itim) &
    586588       / hadley
    … …  
    599601    if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
    600602      tajs(itim) = tajs(itim) + duajs(ilon,ilat,ilev,itim) &
    601        * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
     603       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
    602604       * dmass(ilon,ilat,ilev,itim) &
    603605       / hadley
    … …  
    616618    if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
    617619      tbls(itim) = tbls(itim) + duvdf(ilon,ilat,ilev,itim) &
    618        * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
     620       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
    619621       * dmass(ilon,ilat,ilev,itim) &
    620622       / hadley
    … …  
    633635    if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
    634636      tgwo(itim) = tgwo(itim) + dugwo(ilon,ilat,ilev,itim) &
    635        * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
     637       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
    636638       * dmass(ilon,ilat,ilev,itim) &
    637639       / hadley
    … …  
    650652    if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
    651653      tgwno(itim) = tgwno(itim) + dugwno(ilon,ilat,ilev,itim) &
    652          * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))       &
     654         * rayon(ilon,ilat,ilev,itim)*coslat(ilat)       &
    653655         * dmass(ilon,ilat,ilev,itim)  &
    654656       / hadley

  • trunk/LMDZ.TITAN/Tools/compile_pgf

    r1356 r1454  
    11# path for netcdf should be adapted to your configuration !
    22
    3 \cp -f $2.h planet.h
     3#\cp -f $2.h planet.h
    44 pgf95 -Bstatic cpdet.F90 moyzon.F moyzon2.F moytim.F dx_dp.F epflux.F90 \
    55    io.F90 dmass.F90 reverse.F90 $1.F90 \

  • trunk/LMDZ.TITAN/Tools/compilefft_pgf

    r816 r1454  
    22# path for fftw3  should be adapted to your configuration !
    33
    4 \rm planet.h
    5 ln -s $2.h planet.h
     4#\rm planet.h
     5#ln -s $2.h planet.h
    66 pgf95 -Bstatic cpdet.F90 moyzon.F moyzon2.F moytim.F dx_dp.F epflux.F90 \
    77    io.F90 dmass.F90 reverse.F90 $1.F90 \

  • trunk/LMDZ.TITAN/Tools/dmass.F90

    r816 r1454  
    1 subroutine cellmass(infid,latlength,lonlength,altlength,timelength, &
    2                      lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
     1subroutine cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
     2                     miss_val,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
    33                     dmass )
    44
    … …  
    1414integer timelength ! # of points along time
    1515logical lmdflag ! true=LMD, false=CAM
     16real :: miss_val ! missing value
    1617real :: deltalat,deltalon ! lat and lon intervals in radians
    1718real,dimension(latlength) :: latrad ! latitudes in radians
    … …  
    3233real :: p0 ! GCM reference pressure
    3334integer tmpvarid
    34 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    35 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
    3635
    3736include "planet.h"
    … …  
    4645if(lmdflag) then  ! LMD vs CAM
    4746  text="pres"
    48   call get_var4d(infid,lonlength,latlength,altlength,timelength,text,press,ierr1,ierr2)
     47  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,press,miss_val,ierr1,ierr2)
    4948  if (ierr1.ne.NF_NOERR) then
    5049! assume we are using _P files

  • trunk/LMDZ.TITAN/Tools/dx_dp.F

    r816 r1454  
    3030c-----------------------------------------------------------------------
    3131      do j=1,jjp1
    32          if ((x(j,1).lt.indefini).and.(x(j,2).lt.indefini)) then 
     32         if ((x(j,1).ne.indefini).and.(x(j,2).ne.indefini)) then 
    3333           dxdp(j,1) = (x(j,2)-x(j,1))/(pniv(2) - pniv(1))
    3434         else
    … …  
    3737       
    3838         do l=2,llm-1
    39            if ((x(j,l-1).lt.indefini).and.(x(j,l+1).lt.indefini))then
     39           if ((x(j,l-1).ne.indefini).and.(x(j,l+1).ne.indefini))then
    4040             dxdp(j,l)= (x(j,l+1)-x(j,l-1))/(pniv(l+1) - pniv(l-1))
    41          else if((x(j,l+1).lt.indefini).and.(x(j,l).lt.indefini))then
     41         else if((x(j,l+1).ne.indefini).and.(x(j,l).ne.indefini))then
    4242             dxdp(j,l)= (x(j,l+1)-x(j,l))  /(pniv(l+1) - pniv(l))
    43          else if((x(j,l-1).lt.indefini).and.(x(j,l).lt.indefini))then
     43         else if((x(j,l-1).ne.indefini).and.(x(j,l).ne.indefini))then
    4444             dxdp(j,l)= (x(j,l)-x(j,l-1))  /(pniv(l)   - pniv(l-1))
    4545           else
    … …  
    4747           end if
    4848         end do
    49          if ((x(j,llm).lt.indefini).and.(x(j,llm-1).lt.indefini)) then
     49         if ((x(j,llm).ne.indefini).and.(x(j,llm-1).ne.indefini)) then
    5050         dxdp(j,llm)= (x(j,llm)-x(j,llm-1))/(pniv(llm)-pniv(llm-1))
    5151         else
    … …  
    5858      do j=1,jjp1
    5959         do l = 1, llm
    60            if ( (abs(dxdp(j,l)).gt.(indefini/100.)).and.
     60           if ( (abs(dxdp(j,l)).gt.(abs(indefini)/100.)).and.
    6161     .          (dxdp(j,l).ne.indefini)) then
    6262              write(*,*) '----> j= ', j , ' l= ' ,  l

  • trunk/LMDZ.TITAN/Tools/energy.F90

    r816 r1454  
    3636integer,dimension(4) :: datashape4d ! shape of 4D datasets
    3737
    38 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    39 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
     38real :: miss_val ! special "missing value" to specify missing data
     39real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
    4040real :: pi
    4141real,dimension(:),allocatable :: lon ! longitude
    4242integer lonlength ! # of grid points along longitude
    4343real,dimension(:),allocatable :: lat ! latitude
    44 real,dimension(:),allocatable :: latrad ! latitude in radian
     44real,dimension(:),allocatable :: coslat ! cos of latitude
    4545integer latlength ! # of grid points along latitude
    4646real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
    … …  
    8080
    8181pi = 2.*asin(1.)
     82miss_val = miss_val_def
    8283
    8384write(*,*) ""
    … …  
    123124if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
    124125
    125 allocate(latrad(latlength))
    126 latrad = lat*pi/180.
     126allocate(coslat(latlength))
     127! Beware of rounding problems at poles...
     128coslat(:) = max(0.,cos(lat(:)*pi/180.))
    127129
    128130! Lat, lon pressure intervals
    129 deltalat = abs(latrad(2)-latrad(1))
     131deltalat = abs(lat(2)-lat(1))*pi/180.
    130132deltalon = abs(lon(2)-lon(1))*pi/180.
    131133
    … …  
    193195  text="T"
    194196endif
    195 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     197call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    196198if (ierr1.ne.NF_NOERR) then
    197199  write(*,*) "  looking for t instead... "
    198200  text="t"
    199   call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     201  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    200202  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
    201203endif
    … …  
    216218  text="U"
    217219endif
    218 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
     220call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
    219221if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
    220222if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
    … …  
    228230  text="V"
    229231endif
    230 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
     232call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
    231233if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
    232234if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
    … …  
    242244
    243245!text="zareoid"
    244 !call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
     246!call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
    245247!if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
    246248!if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
    … …  
    276278enddo ! timelength
    277279
    278 call cellmass(infid,latlength,lonlength,altlength,timelength, &
    279               lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
     280call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
     281              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
    280282              dmass )
    281283

  • trunk/LMDZ.TITAN/Tools/epflux.F90

    r816 r1454  
    1 subroutine  epflux(iip1,jjp1,llm,indefini,rlatu,rbar &
     1subroutine  epflux(iip1,jjp1,llm,indefini,latdeg,rbar &
    22                  ,teta,u3d,v3d,w3d,press &
    33                  ,epy2d,epz2d,divep2d,vtem2d,wtem2d,acc_meridien2d &
    … …  
    4141     
    4242      integer :: iip1,jjp1,llm
    43       real :: indefini,rlatu(jjp1),rbar(jjp1,llm)
     43      real :: indefini,latdeg(jjp1),rbar(jjp1,llm)
    4444      REAL :: teta(iip1,jjp1,llm)
    4545      REAL :: u3d(iip1,jjp1,llm),v3d(iip1,jjp1,llm)
    … …  
    7575!  Autres
    7676!  ------
     77      real :: rlatu(jjp1) ! lat in radians (beware of rounding effects at poles)
    7778      real :: pi
    7879      REAL :: tetas(llm)
    … …  
    8889      jjm = jjp1-1
    8990      pi = 2.*asin(1.)
     91
     92! To avoid rounding effects at poles
     93      rlatu(:)=latdeg(:)*pi/180.00001
    9094
    9195!     Calcul de moyennes zonales
    … …  
    108112         n= 0
    109113         do j=2,jjm
    110             if (tetabar(j,l).lt.indefini) then
     114            if (tetabar(j,l).ne.indefini) then
    111115               tetas(l) = tetas(l) + tetabar(j,l)
    112116               n = n+1
    … …  
    237241      do l=1,llm
    238242
    239        if ((ubar(2,l).lt.indefini).and.(ubar(1,l).lt.indefini))then
     243       if ((ubar(2,l).ne.indefini).and.(ubar(1,l).ne.indefini))then
    240244      ducosdlat(1,l) = ( ubar(2,l)*cos(rlatu(2))-ubar(1,l)*cos(rlatu(1)) ) &
    241245                      /( rlatu(2) - rlatu(1))
    … …  
    244248       end if
    245249        do j=2,jjm
    246        if ((ubar(j+1,l).lt.indefini).and.(ubar(j-1,l).lt.indefini))then
     250       if ((ubar(j+1,l).ne.indefini).and.(ubar(j-1,l).ne.indefini))then
    247251      ducosdlat(j,l) = ( ubar(j+1,l)*cos(rlatu(j+1))-ubar(j-1,l)*cos(rlatu(j-1))) &
    248252                      /( rlatu(j+1) - rlatu(j-1))
    … …  
    251255       end if
    252256        enddo
    253        if ((ubar(jjp1,l).lt.indefini).and.(ubar(jjm,l).lt.indefini))then
     257       if ((ubar(jjp1,l).ne.indefini).and.(ubar(jjm,l).ne.indefini))then
    254258      ducosdlat(jjp1,l) = ( ubar(jjp1,l)*cos(rlatu(jjp1))  &
    255259                               -ubar(jjm,l)*cos(rlatu(jjm))  ) &
    … …  
    260264           
    261265        do j=1,jjp1
    262        if ((vptetapbar(j,l).lt.indefini).and. &
    263            (dtetabardp(j,l).lt.indefini)) then
     266       if ((vptetapbar(j,l).ne.indefini).and. &
     267           (dtetabardp(j,l).ne.indefini)) then
    264268           vz(j,l) = vptetapbar(j,l)/dtetabardp(j,l)
    265269         wlat(j,l) = cos(rlatu(j))*vz(j,l)
    … …  
    284288!     ------------
    285289
    286             if (   (rbar(j,l).lt.indefini).and. &
    287                 (vpupbar(j,l).lt.indefini).and. &
    288                 (dubardp(j,l).lt.indefini).and. &
    289                      (vz(j,l).lt.indefini)) then
     290            if (   (rbar(j,l).ne.indefini).and. &
     291                (vpupbar(j,l).ne.indefini).and. &
     292                (dubardp(j,l).ne.indefini).and. &
     293                     (vz(j,l).ne.indefini)) then
    290294           
    291295     epy2d(j,l) = rbar(j,l) * cos(rlatu(j))* ( vpupbar(j,l) &
    … …  
    303307!     ------------
    304308
    305             if (   (rbar(j,l).lt.indefini).and. &
    306                 (wpupbar(j,l).lt.indefini).and. &
    307               (ducosdlat(j,l).lt.indefini).and. &
    308                      (vz(j,l).lt.indefini)) then
     309            if (   (rbar(j,l).ne.indefini).and. &
     310                (wpupbar(j,l).ne.indefini).and. &
     311              (ducosdlat(j,l).ne.indefini).and. &
     312                     (vz(j,l).ne.indefini)) then
    309313
    310314     epz2d(j,l) = rbar(j,l) * cos(rlatu(j))*                    &
    … …  
    330334      do l=1,llm
    331335
    332           if ( (wlat(2,l).lt.indefini) .and. &
    333                (wlat(1,l).lt.indefini) ) then
     336          if ( (wlat(2,l).ne.indefini) .and. &
     337               (wlat(1,l).ne.indefini) ) then
    334338               dwlatdlat(1,l)=(wlat(2,l)-wlat(1,l)) &
    335339                     /(rlatu(2)-rlatu(1))
    … …  
    338342          end if
    339343        do j=2,jjm
    340           if ( (wlat(j+1,l).lt.indefini) .and. &
    341                (wlat(j-1,l).lt.indefini) ) then
     344          if ( (wlat(j+1,l).ne.indefini) .and. &
     345               (wlat(j-1,l).ne.indefini) ) then
    342346               dwlatdlat(j,l)=(wlat(j+1,l)-wlat(j-1,l)) &
    343347                     /(rlatu(j+1)-rlatu(j-1))
    … …  
    346350          end if
    347351        enddo
    348           if ( (wlat(jjp1,l).lt.indefini) .and. &
    349                (wlat(jjm, l).lt.indefini) ) then
     352          if ( (wlat(jjp1,l).ne.indefini) .and. &
     353               (wlat(jjm, l).ne.indefini) ) then
    350354               dwlatdlat(jjp1,l)=(wlat(jjp1,l)-wlat(jjm,l)) &
    351355                     /(rlatu(jjp1)-rlatu(jjm))
    … …  
    355359
    356360        do j=1,jjp1
    357           if ( (dvzdp(j,l).lt.indefini) &
    358            .and.(vbar(j,l).lt.indefini)) then
     361          if ( (dvzdp(j,l).ne.indefini) &
     362           .and.(vbar(j,l).ne.indefini)) then
    359363            vtem2d(j,l) = vbar(j,l)-dvzdp(j,l)
    360364          else
    … …  
    364368
    365369        do j=1,jjp1
    366           if (   (rbar(j,l).lt.indefini).and. &
    367             (dwlatdlat(j,l).lt.indefini).and. &
    368                  (wbar(j,l).lt.indefini)) then
     370          if (   (rbar(j,l).ne.indefini).and. &
     371            (dwlatdlat(j,l).ne.indefini).and. &
     372                 (wbar(j,l).ne.indefini)) then
    369373           wtem2d(j,l) = wbar(j,l)+dwlatdlat(j,l)/(rbar(j,l)*cos(rlatu(j)))
    370374          else
    … …  
    385389      do l=1,llm
    386390        do j=1,jjp1
    387           if (   (rbar(j,l).lt.indefini) &
    388           .and.(vtem2d(j,l).lt.indefini) &
    389           .and.(wtem2d(j,l).lt.indefini) &
    390        .and.(ducosdlat(j,l).lt.indefini) &
    391          .and.(dubardp(j,l).lt.indefini)) then
     391          if (   (rbar(j,l).ne.indefini) &
     392          .and.(vtem2d(j,l).ne.indefini) &
     393          .and.(wtem2d(j,l).ne.indefini) &
     394       .and.(ducosdlat(j,l).ne.indefini) &
     395         .and.(dubardp(j,l).ne.indefini)) then
    392396       acc_meridien2d(j,l) =  &
    393397           vtem2d(j,l)*(ducosdlat(j,l)/(rbar(j,l)*cos(rlatu(j))) &
    … …  
    410414        divep2d(1,l) =0
    411415        do j=2,jjm
    412           if (    (rbar(j,l).lt.indefini) .and. &
    413                (epy2d(j+1,l).lt.indefini) .and. &
    414                (epy2d(j-1,l).lt.indefini) .and. &
    415                 (depzdp(j,l).lt.indefini)    ) then
     416          if (    (rbar(j,l).ne.indefini) .and. &
     417               (epy2d(j+1,l).ne.indefini) .and. &
     418               (epy2d(j-1,l).ne.indefini) .and. &
     419                (depzdp(j,l).ne.indefini)    ) then
    416420
    417421           depycosdlat=(epy2d(j+1,l)*cos(rlatu(j+1))   &

  • trunk/LMDZ.TITAN/Tools/fft.F90

    r888 r1454  
    6363integer,dimension(4) :: datashape4d ! shape of 4D datasets
    6464
    65 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    66 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
     65real :: miss_val ! special "missing value" to specify missing data
     66real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
    6767real :: pi
    6868real,dimension(:),allocatable :: lon ! longitude
    6969integer lonlength ! # of grid points along longitude
    7070real,dimension(:),allocatable :: lat ! latitude
    71 real,dimension(:),allocatable :: latrad ! latitude in radian
    7271integer latlength ! # of grid points along latitude
    7372real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
    … …  
    125124logical :: lmdflag
    126125
     126! Tuning parameters
     127real :: fcoup1,fcoup2,width
     128real :: fcoup1tmp,fcoup2tmp,widthtmp
     129logical,dimension(4) :: ok_out
     130character (len=1) :: ok_outtmp
     131
    127132include "planet.h"
    128 include "filter.h"
    129133
    130134#include <fftw3.f>
    … …  
    135139
    136140pi = 2.*asin(1.)
     141miss_val = miss_val_def
    137142
    138143write(*,*) ""
    139144write(*,*) "You are working on the atmosphere of ",planet
     145
     146! initialisation
     147!----------------
     148
     149! Par defaut
     150
     151! Define the filters
     152! Low  cutting frequency, in Hz    : fcoup1
     153fcoup1=2.5e-6
     154! High cutting frequency, in Hz    : fcoup2
     155fcoup2=6.5e-6
     156! Half-width of the filters, in Hz : width
     157width=4.e-7
     158! Outputs (U,     V,      W,     T)
     159ok_out=(/.true.,.true.,.false.,.true./)
     160
     161print*,"Low  cutting frequency, in Hz ?"
     162print*,"between 1e-5 and 1e-7, 0 for default => 2.5e-6"
     163read(*,*) fcoup1tmp
     164if ((fcoup1tmp.lt.1e-5).and.(fcoup1tmp.gt.1e-7)) fcoup1=fcoup1tmp
     165print*,"=",fcoup1
     166print*,"High cutting frequency, in Hz ?"
     167print*,"between 1e-5 and 1e-7, 0 for default => 6.5e-6"
     168read(*,*) fcoup2tmp
     169if ((fcoup2tmp.lt.1e-5).and.(fcoup2tmp.gt.1e-7)) fcoup2=fcoup2tmp
     170print*,"=",fcoup2
     171print*,"Half-width of the filters, in Hz ?"
     172print*,"between 1e-6 and 1e-8, 0 for default => 4e-7)"
     173read(*,*) widthtmp
     174if ((widthtmp.lt.1e-6).and.(widthtmp.gt.1e-8)) width=widthtmp
     175print*,"=",width
     176!width = 3./time(timelength)
     177
     178! Outputs
     179print*,"Output of zonal wind ? (y or n, default is y)"
     180read(*,'(a1)') ok_outtmp
     181if (ok_outtmp.eq."n") ok_out(1)=.false.
     182print*,"=",ok_out(1)
     183print*,"Output of meridional wind ? (y or n, default is y)"
     184read(*,'(a1)') ok_outtmp
     185if (ok_outtmp.eq."n") ok_out(2)=.false.
     186print*,"=",ok_out(2)
     187print*,"Output of vertical wind ? (y or n, default is n)"
     188read(*,'(a1)') ok_outtmp
     189if (ok_outtmp.eq."y") ok_out(3)=.true.
     190print*,"=",ok_out(3)
     191print*,"Output of temperature ? (y or n, default is y)"
     192read(*,'(a1)') ok_outtmp
     193if (ok_outtmp.eq."n") ok_out(4)=.false.
     194print*,"=",ok_out(4)
    140195
    141196!===============================================================================
    … …  
    172227ierr=NF_GET_VAR_REAL(infid,lat_varid,lat)
    173228if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
    174 
    175 allocate(latrad(latlength))
    176 latrad = lat*pi/180.
    177229
    178230allocate(plev(altlength))
    … …  
    212264
    213265text="temp"
    214 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     266call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    215267if (ierr1.ne.NF_NOERR) then
    216268  write(*,*) "  looking for t instead... "
    217269  text="t"
    218   call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     270  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    219271  if (ierr1.ne.NF_NOERR) then
    220272    print*,"Error: Failed to get temperature ID"
    … …  
    236288
    237289text="vitu"
    238 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
     290call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
    239291if (ierr1.ne.NF_NOERR) then
    240292  print*,"Error: Failed to get vitu ID"
    … …  
    252304
    253305text="vitv"
    254 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
     306call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
    255307if (ierr1.ne.NF_NOERR) then
    256308  print*,"Error: Failed to get vitv ID"
    … …  
    268320
    269321text="vitw"
    270 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
     322call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
    271323if (ierr1.ne.NF_NOERR) then
    272324  print*,"Error: Failed to get vitw ID"
    … …  
    288340! 2.2.1 FFT and filtering
    289341!===============================================================================
     342
    290343! allocations
    291344!-------------
    … …  
    405458enddo
    406459
    407 ! Define the filters
    408 
    409 ! IN filter.h :
    410 ! Low  cutting frequency, in Hz    : fcoup1
    411 ! High cutting frequency, in Hz    : fcoup2
    412 ! Half-width of the filters, in Hz : width
    413 
    414 !print*,"Low  cutting frequency, in Hz ?"
    415 !read(*,*) fcoup1
    416 !print*,"High cutting frequency, in Hz ?"
    417 !read(*,*) fcoup2
    418 !print*,"Half-width of the filters, in Hz ?"
    419 !read(*,*) width
    420 !width = 3./time(timelength)
    421 
    422460do itim=1,M_fft+1
    423461  if (freq(itim).lt.(fcoup1-width)) then

  • trunk/LMDZ.TITAN/Tools/io.F90

    r816 r1454  
    218218!===========================================================================
    219219
    220 subroutine get_var4d(infid,dim1,dim2,dim3,dim4,text,var,ierr1,ierr2)
     220subroutine get_var4d(infid,dim1,dim2,dim3,dim4,text,var,missing,ierr1,ierr2)
    221221
    222222implicit none
    … …  
    229229character (len=64) :: text ! name of variable to read
    230230real,dimension(dim1,dim2,dim3,dim4) :: var ! variable to read
    231 integer :: ierr1,ierr2 ! NetCDF routines return code
     231real :: missing ! missing value
     232integer :: ierr1,ierr2,miss ! NetCDF routines return code
    232233
    233234! local
    … …  
    240241  write(*,*) "ID ok for ",trim(text)
    241242  ierr2=NF_GET_VAR_REAL(infid,tmpvarid,var)
     243  miss=NF_GET_ATT_REAL(infid,tmpvarid,"missing_value",missing)
    242244endif
    243245

  • trunk/LMDZ.TITAN/Tools/moytim.F

    r816 r1454  
    3333           n = 0
    3434           do t=1,nt
    35               if (x(i,j,l,t).lt.indefini) then
     35              if (x(i,j,l,t).ne.indefini) then
    3636                 xmean(i,j,l) = xmean(i,j,l) + x(i,j,l,t)
    3737                 n = n+1

  • trunk/LMDZ.TITAN/Tools/moyzon.F

    r816 r1454  
    3232           n = 0
    3333           do i=1,iim
    34               if (x(i,j,l).lt.indefini) then
     34              if (x(i,j,l).ne.indefini) then
    3535                 xbar(j,l) = xbar(j,l) + x(i,j,l)
    3636                 n = n+1

  • trunk/LMDZ.TITAN/Tools/moyzon2.F

    r816 r1454  
    3434           n = 0
    3535           do i=1,iim
    36               if (x(i,j,l).lt.indefini) then
     36              if (x(i,j,l).ne.indefini) then
    3737                 xbar(j,l) = xbar(j,l) + x(i,j,l)
    3838                 n = n+1

  • trunk/LMDZ.TITAN/Tools/psi.F90

    r1120 r1454  
    3131integer,dimension(4) :: datashape4d ! shape of 4D datasets
    3232
    33 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    34 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
     33real :: miss_val ! special "missing value" to specify missing data
     34real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
    3535real :: pi
    3636real,dimension(:),allocatable :: lon ! longitude
    3737integer lonlength ! # of grid points along longitude
    3838real,dimension(:),allocatable :: lat ! latitude
    39 real,dimension(:),allocatable :: latrad ! latitude in radian
     39real,dimension(:),allocatable :: coslat ! cos of latitude
    4040integer latlength ! # of grid points along latitude
    4141real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
    … …  
    6666
    6767pi = 2.*asin(1.)
     68miss_val = miss_val_def
    6869
    6970write(*,*) ""
    … …  
    104105if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
    105106
    106 allocate(latrad(latlength))
    107 latrad = lat*pi/180.
     107allocate(coslat(latlength))
     108! Beware of rounding problems at poles...
     109coslat(:) = max(0.,cos(lat(:)*pi/180.))
    108110
    109111! Lat, lon pressure intervals
    110 deltalat = abs(latrad(2)-latrad(1))
     112deltalat = abs(lat(2)-lat(1))*pi/180.
    111113deltalon = abs(lon(2)-lon(1))*pi/180.
    112114
    … …  
    160162! meridional wind vitv (in m/s)
    161163text="vitv"
    162 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
     164call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
    163165if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
    164166if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
    … …  
    173175
    174176text="zareoid"
    175 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
     177call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
    176178if (ierr1.ne.NF_NOERR) then
    177179  write(*,*) "  looking for geop instead... "
    178180  text="geop"
    179   call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
     181  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
    180182  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get geop ID"
    181183do itim=1,timelength
    … …  
    225227enddo ! timelength
    226228
    227 call cellmass(infid,latlength,lonlength,altlength,timelength, &
    228               lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
     229call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
     230              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
    229231              dmass )
    230232

  • trunk/LMDZ.TITAN/Tools/stability.F90

    r816 r1454  
    3131integer,dimension(4) :: datashape4d ! shape of 4D datasets
    3232
    33 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    34 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
     33real :: miss_val ! special "missing value" to specify missing data
     34real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
    3535real :: pi
    3636real,dimension(:),allocatable :: lon ! longitude
    3737integer lonlength ! # of grid points along longitude
    3838real,dimension(:),allocatable :: lat ! latitude
    39 real,dimension(:),allocatable :: latrad ! latitude in radian
     39real,dimension(:),allocatable :: coslat ! cos of latitude
    4040integer latlength ! # of grid points along latitude
    4141real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
    … …  
    9191
    9292pi = 2.*asin(1.)
     93miss_val = miss_val_def
    9394
    9495write(*,*) ""
    … …  
    129130if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
    130131
    131 allocate(latrad(latlength))
    132 latrad = lat*pi/180.
     132allocate(coslat(latlength))
     133! Beware of rounding problems at poles...
     134coslat(:) = max(0.,cos(lat(:)*pi/180.))
    133135
    134136! Lat, lon pressure intervals
    135 deltalat = abs(latrad(2)-latrad(1))
     137deltalat = abs(lat(2)-lat(1))*pi/180.
    136138deltalon = abs(lon(2)-lon(1))*pi/180.
    137139
    … …  
    180182
    181183text="temp"
    182 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     184call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    183185if (ierr1.ne.NF_NOERR) then
    184186  write(*,*) "  looking for t instead... "
    185187  text="t"
    186   call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     188  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    187189  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
    188190endif
    … …  
    197199! zonal wind vitu (in m/s)
    198200text="vitu"
    199 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
     201call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
    200202if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
    201203if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
    … …  
    203205! meridional wind vitv (in m/s)
    204206text="vitv"
    205 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
     207call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
    206208if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
    207209if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
    … …  
    215217
    216218! text="zareoid"
    217 ! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
     219! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
    218220! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
    219221! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
    … …  
    276278 grav(lonlength+1,:,:,:) =  grav(1,:,:,:)
    277279
    278 call cellmass(infid,latlength,lonlength+1,altlength,timelength, &
    279               lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
     280call cellmass(infid,latlength,lonlength+1,altlength,timelength,lmdflag, &
     281              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
    280282              dmass )
    281283
    … …  
    364366
    365367do ilat=2,latlength-1
    366    if (tan(latrad(ilat)).ne.0.) then
    367      fac1 = R0/tan(latrad(ilat))
     368   if (tan(lat(ilat)*pi/180.).ne.0.) then
     369     fac1 = R0/tan(lat(ilat)*pi/180.)
    368370   else
    369371     fac1 = miss_val

  • trunk/LMDZ.TITAN/Tools/tem.F90

    r816 r1454  
    3737integer,dimension(3) :: datashape3d ! shape of 3D datasets
    3838
    39 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    40 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
     39real :: miss_val ! special "missing value" to specify missing data
     40real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
    4141real :: pi
    4242real,dimension(:),allocatable :: lon ! longitude
    4343integer lonlength ! # of grid points along longitude
    4444real,dimension(:),allocatable :: lat ! latitude
    45 real,dimension(:),allocatable :: latrad ! latitude in radian
     45real,dimension(:),allocatable :: coslat ! cos of latitude
    4646integer latlength ! # of grid points along latitude
    4747real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
    … …  
    107107
    108108pi = 2.*asin(1.)
     109miss_val = miss_val_def
    109110
    110111write(*,*) ""
    … …  
    145146if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
    146147
    147 allocate(latrad(latlength))
    148 latrad = lat*pi/180.
     148allocate(coslat(latlength))
     149! Beware of rounding problems at poles...
     150coslat(:) = max(0.,cos(lat(:)*pi/180.))
    149151
    150152! Lat, lon pressure intervals
    151 deltalat = abs(latrad(2)-latrad(1))
     153deltalat = abs(lat(2)-lat(1))*pi/180.
    152154deltalon = abs(lon(2)-lon(1))*pi/180.
    153155
    … …  
    196198
    197199text="temp"
    198 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     200call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    199201if (ierr1.ne.NF_NOERR) then
    200202  write(*,*) "  looking for t instead... "
    201203  text="t"
    202   call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     204  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    203205  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
    204206endif
    … …  
    214216! zonal wind vitu (in m/s)
    215217text="vitu"
    216 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
     218call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
    217219if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
    218220if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
    … …  
    220222! meridional wind vitv (in m/s)
    221223text="vitv"
    222 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
     224call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
    223225if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
    224226if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
    … …  
    226228! vertical wind vitw (in Pa/s)
    227229text="vitw"
    228 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
     230call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
    229231if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitw ID"
    230232if (ierr2.ne.NF_NOERR) stop "Error: Failed reading vertical wind"
    … …  
    238240
    239241! text="zareoid"
    240 ! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
     242! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
    241243! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
    242244! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
    … …  
    312314 grav(lonlength+1,:,:,:) =  grav(1,:,:,:)
    313315
    314 call cellmass(infid,latlength,lonlength+1,altlength,timelength, &
    315               lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
     316call cellmass(infid,latlength,lonlength+1,altlength,timelength,lmdflag, &
     317              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
    316318              dmass )
    317319
    … …  
    359361call moyzon(lonlength,latlength,altlength,miss_val,u3d,ubar)
    360362
    361 call epflux(lonlength+1,latlength,altlength,miss_val,latrad,rbar &
     363call epflux(lonlength+1,latlength,altlength,miss_val,lat,rbar &
    362364            ,teta,u3d,v3d,w3d,plev &
    363365            ,epy2d,epz2d,div2d,vtem2d,wtem2d,ammc2d &
    … …  
    376378! rsurg: r*dp/g = dm/(r cos(lat) dlon dlat) !!!
    377379     rsurg(ilon,ilat,ilev) = dmass(ilon,ilat,ilev,itim) &
    378           / (r3d(ilon,ilat,ilev)*cos(latrad(ilat))*deltalon*deltalat)
     380          / (r3d(ilon,ilat,ilev)*coslat(ilat)*deltalon*deltalat)
    379381    else
    380382     rsurg(ilon,ilat,ilev) = miss_val
    … …  
    391393        (rsurgbar(ilat,ilev).ne.miss_val) ) then
    392394      vm(ilat,ilev) = vtem2d(ilat,ilev) &
    393             * 2.*pi*rsurgbar(ilat,ilev)*cos(latrad(ilat))
     395            * 2.*pi*rsurgbar(ilat,ilev)*coslat(ilat)
    394396    else
    395397      vm(ilat,ilev) = miss_val

  • trunk/LMDZ.TITAN/Tools/tmc.F90

    r816 r1454  
    4949integer,dimension(4) :: datashape4d ! shape of 4D datasets
    5050
    51 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    52 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
     51real :: miss_val ! special "missing value" to specify missing data
     52real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
    5353real :: pi
    5454real,dimension(:),allocatable :: lon ! longitude
    5555integer lonlength ! # of grid points along longitude
    5656real,dimension(:),allocatable :: lat ! latitude
    57 real,dimension(:),allocatable :: latrad ! latitude in radian
     57real,dimension(:),allocatable :: coslat ! cos of latitude
    5858integer latlength ! # of grid points along latitude
    5959real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
    … …  
    141141
    142142pi = 2.*asin(1.)
     143miss_val = miss_val_def
    143144
    144145write(*,*) ""
    … …  
    179180if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
    180181
    181 allocate(latrad(latlength))
    182 latrad = lat*pi/180.
     182allocate(coslat(latlength))
     183! Beware of rounding problems at poles...
     184coslat(:) = max(0.,cos(lat(:)*pi/180.))
    183185
    184186! Lat, lon pressure intervals
    185 deltalat = abs(latrad(2)-latrad(1))
     187deltalat = abs(lat(2)-lat(1))*pi/180.
    186188deltalon = abs(lon(2)-lon(1))*pi/180.
    187189
    … …  
    233235! zonal wind vitu (in m/s)
    234236text="vitu"
    235 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
     237call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
    236238if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
    237239if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
    … …  
    239241! meridional wind vitv (in m/s)
    240242text="vitv"
    241 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
     243call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
    242244if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
    243245if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
    … …  
    245247! vertical wind vitw (in Pa/s)
    246248text="vitw"
    247 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
     249call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
    248250if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitw ID"
    249251if (ierr2.ne.NF_NOERR) stop "Error: Failed reading vertical wind"
    … …  
    257259
    258260! text="zareoid"
    259 ! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
     261! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
    260262! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
    261263! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
    … …  
    293295enddo ! timelength
    294296
    295 call cellmass(infid,latlength,lonlength,altlength,timelength, &
    296               lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
     297call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
     298              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
    297299              dmass )
    298300
    … …  
    316318      osam(ilon,ilat,ilev,itim) = &
    317319         rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim) &
    318        * cos(latrad(ilat))*cos(latrad(ilat)) &
     320       * coslat(ilat)*coslat(ilat) &
    319321       * omega
    320322      rsam(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim) &
    321        * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))
     323       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)
    322324      tsam(ilon,ilat,ilev,itim) = osam(ilon,ilat,ilev,itim)&
    323325                                + rsam(ilon,ilat,ilev,itim)

  • trunk/LMDZ.VENUS/Tools/angmom.F90

    r1356 r1454  
    4747integer,dimension(4) :: datashape4d ! shape of 4D datasets
    4848
    49 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    50 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
     49real :: miss_val ! special "missing value" to specify missing data
     50real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
    5151real :: pi
    5252real,dimension(:),allocatable :: lon ! longitude
    5353integer lonlength ! # of grid points along longitude
    5454real,dimension(:),allocatable :: lat ! latitude
    55 real,dimension(:),allocatable :: latrad ! latitude in radian
     55real,dimension(:),allocatable :: coslat ! cos of latitude
    5656integer latlength ! # of grid points along latitude
    5757real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
    … …  
    117117
    118118pi = 2.*asin(1.)
     119miss_val = miss_val_def
    119120
    120121write(*,*) ""
    … …  
    160161if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
    161162
    162 allocate(latrad(latlength))
    163 latrad = lat*pi/180.
     163allocate(coslat(latlength))
     164! Beware of rounding problems at poles...
     165coslat(:) = max(0.,cos(lat(:)*pi/180.))
    164166
    165167! Lat, lon pressure intervals
    166 deltalat = abs(latrad(2)-latrad(1))
     168deltalat = abs(lat(2)-lat(1))*pi/180.
    167169deltalon = abs(lon(2)-lon(1))*pi/180.
    168170
    … …  
    281283endif
    282284
    283 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
     285call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
    284286if (ierr1.ne.NF_NOERR) stop "Error: Failed to get U ID"
    285287if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
    … …  
    295297
    296298!text="zareoid"
    297 !call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
     299!call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
    298300!if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
    299301!if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
    … …  
    309311  text="DUVDF"
    310312endif
    311 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duvdf,ierr1,ierr2)
     313call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duvdf,miss_val,ierr1,ierr2)
    312314if (ierr1.ne.NF_NOERR) then
    313315  write(*,*) "Failed to get duvdf ID"
    … …  
    357359
    358360text="dugwo"
    359 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwo,ierr1,ierr2)
     361call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwo,miss_val,ierr1,ierr2)
    360362if (ierr1.ne.NF_NOERR) then
    361363  write(*,*) "Failed to get dugwo ID"
    … …  
    367369
    368370text="dugwno"
    369 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwno,ierr1,ierr2)
     371call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwno,miss_val,ierr1,ierr2)
    370372if (ierr1.ne.NF_NOERR) then
    371373  write(*,*) "Failed to get dugwno ID"
    … …  
    390392
    391393text="duajs"
    392 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duajs,ierr1,ierr2)
     394call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duajs,miss_val,ierr1,ierr2)
    393395if (ierr1.ne.NF_NOERR) then
    394396  write(*,*) "Failed to get duajs ID"
    … …  
    414416  text="DUDYN"
    415417endif
    416 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dudyn,ierr1,ierr2)
     418call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dudyn,miss_val,ierr1,ierr2)
    417419if (ierr1.ne.NF_NOERR) then
    418420  write(*,*) "Failed to get dudyn ID"
    … …  
    486488enddo ! timelength
    487489
    488 call cellmass(infid,latlength,lonlength,altlength,timelength, &
    489               lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
     490call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
     491              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
    490492              dmass )
    491493
    … …  
    504506      osam(ilon,ilat,ilev,itim) = &
    505507         rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim) &
    506        * cos(latrad(ilat))*cos(latrad(ilat)) &
     508       * coslat(ilat)*coslat(ilat) &
    507509       * omega
    508510      rsam(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim) &
    509        * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))
     511       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)
    510512    else
    511513      osam(ilon,ilat,ilev,itim) = miss_val
    … …  
    569571       tmpp = (ps(ilon+1,ilat,itim)  +ps(ilon,ilat,itim))/2.
    570572      endif
    571       tmou(itim) = tmou(itim) + a0*a0* deltalat*cos(latrad(ilat)) &
     573      tmou(itim) = tmou(itim) + a0*a0* deltalat*coslat(ilat) &
    572574       * signe*dz*tmpp &
    573575       / hadley
    … …  
    582584    if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
    583585      tdyn(itim) = tdyn(itim) + dudyn(ilon,ilat,ilev,itim) &
    584        * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
     586       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
    585587       * dmass(ilon,ilat,ilev,itim) &
    586588       / hadley
    … …  
    599601    if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
    600602      tajs(itim) = tajs(itim) + duajs(ilon,ilat,ilev,itim) &
    601        * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
     603       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
    602604       * dmass(ilon,ilat,ilev,itim) &
    603605       / hadley
    … …  
    616618    if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
    617619      tbls(itim) = tbls(itim) + duvdf(ilon,ilat,ilev,itim) &
    618        * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
     620       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
    619621       * dmass(ilon,ilat,ilev,itim) &
    620622       / hadley
    … …  
    633635    if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
    634636      tgwo(itim) = tgwo(itim) + dugwo(ilon,ilat,ilev,itim) &
    635        * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
     637       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
    636638       * dmass(ilon,ilat,ilev,itim) &
    637639       / hadley
    … …  
    650652    if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
    651653      tgwno(itim) = tgwno(itim) + dugwno(ilon,ilat,ilev,itim) &
    652          * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))       &
     654         * rayon(ilon,ilat,ilev,itim)*coslat(ilat)       &
    653655         * dmass(ilon,ilat,ilev,itim)  &
    654656       / hadley

  • trunk/LMDZ.VENUS/Tools/compilefft_pgf

    r816 r1454  
    22# path for fftw3  should be adapted to your configuration !
    33
    4 \rm planet.h
    5 ln -s $2.h planet.h
     4#\rm planet.h
     5#ln -s $2.h planet.h
    66 pgf95 -Bstatic cpdet.F90 moyzon.F moyzon2.F moytim.F dx_dp.F epflux.F90 \
    77    io.F90 dmass.F90 reverse.F90 $1.F90 \

  • trunk/LMDZ.VENUS/Tools/dmass.F90

    r816 r1454  
    1 subroutine cellmass(infid,latlength,lonlength,altlength,timelength, &
    2                      lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
     1subroutine cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
     2                     miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
    33                     dmass )
    44
    … …  
    1414integer timelength ! # of points along time
    1515logical lmdflag ! true=LMD, false=CAM
     16real :: miss_val ! missing value
    1617real :: deltalat,deltalon ! lat and lon intervals in radians
    17 real,dimension(latlength) :: latrad ! latitudes in radians
     18real,dimension(latlength) :: coslat ! cos of latitudes
    1819real,dimension(altlength) :: plev ! Pressure levels (Pa)
    1920real,dimension(lonlength,latlength,timelength),intent(in) :: ps ! surface pressure
    … …  
    3233real :: p0 ! GCM reference pressure
    3334integer tmpvarid
    34 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    35 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
    3635
    3736include "planet.h"
    … …  
    4645if(lmdflag) then  ! LMD vs CAM
    4746  text="pres"
    48   call get_var4d(infid,lonlength,latlength,altlength,timelength,text,press,ierr1,ierr2)
     47  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,press,miss_val,ierr1,ierr2)
    4948  if (ierr1.ne.NF_NOERR) then
    5049! assume we are using _P files
    … …  
    119118 do ilat=1,latlength
    120119  tmpp=ps(ilon,ilat,itim)
    121   pp=press(ilon,ilat,:,itim)
     120! Beware when miss_val is negative...
     121  do ilev=1,altlength
     122    if (press(ilon,ilat,ilev,itim).eq.miss_val) then
     123      pp(ilev)=1.e33
     124    else
     125      pp(ilev)=press(ilon,ilat,ilev,itim)
     126    endif
     127  enddo
    122128if (tmpp.ge.pp(1)) then
    123129  deltap(ilon,ilat,1,itim)= tmpp - exp((log(pp(1))+log(pp(2)))/2.) ! initialization, rectified with ps below
    … …  
    136142 do ilat=1,latlength
    137143  tmpp=ps(ilon,ilat,itim)
    138   pp=press(ilon,ilat,:,itim)
     144! Beware when miss_val is negative...
     145  do ilev=1,altlength
     146    if (press(ilon,ilat,ilev,itim).eq.miss_val) then
     147      pp(ilev)=1.e33
     148    else
     149      pp(ilev)=press(ilon,ilat,ilev,itim)
     150    endif
     151  enddo
    139152  do ilev=altlength,2,-1
    140153    if ((pp(ilev).le.tmpp).and.(pp(ilev-1).gt.tmpp)) then
    … …  
    163176   .and.  (grav(ilon,ilat,ilev,itim).ne.miss_val)) then
    164177     dmass(ilon,ilat,ilev,itim) = &
    165                rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))*deltalon &
     178               rayon(ilon,ilat,ilev,itim)*coslat(ilat)*deltalon &
    166179             * rayon(ilon,ilat,ilev,itim)*deltalat &
    167180             * deltap(ilon,ilat,ilev,itim)/grav(ilon,ilat,ilev,itim)

  • trunk/LMDZ.VENUS/Tools/dx_dp.F

    r816 r1454  
    3030c-----------------------------------------------------------------------
    3131      do j=1,jjp1
    32          if ((x(j,1).lt.indefini).and.(x(j,2).lt.indefini)) then 
     32         if ((x(j,1).ne.indefini).and.(x(j,2).ne.indefini)) then 
    3333           dxdp(j,1) = (x(j,2)-x(j,1))/(pniv(2) - pniv(1))
    3434         else
    … …  
    3737       
    3838         do l=2,llm-1
    39            if ((x(j,l-1).lt.indefini).and.(x(j,l+1).lt.indefini))then
     39           if ((x(j,l-1).ne.indefini).and.(x(j,l+1).ne.indefini))then
    4040             dxdp(j,l)= (x(j,l+1)-x(j,l-1))/(pniv(l+1) - pniv(l-1))
    41          else if((x(j,l+1).lt.indefini).and.(x(j,l).lt.indefini))then
     41         else if((x(j,l+1).ne.indefini).and.(x(j,l).ne.indefini))then
    4242             dxdp(j,l)= (x(j,l+1)-x(j,l))  /(pniv(l+1) - pniv(l))
    43          else if((x(j,l-1).lt.indefini).and.(x(j,l).lt.indefini))then
     43         else if((x(j,l-1).ne.indefini).and.(x(j,l).ne.indefini))then
    4444             dxdp(j,l)= (x(j,l)-x(j,l-1))  /(pniv(l)   - pniv(l-1))
    4545           else
    … …  
    4747           end if
    4848         end do
    49          if ((x(j,llm).lt.indefini).and.(x(j,llm-1).lt.indefini)) then
     49         if ((x(j,llm).ne.indefini).and.(x(j,llm-1).ne.indefini)) then
    5050         dxdp(j,llm)= (x(j,llm)-x(j,llm-1))/(pniv(llm)-pniv(llm-1))
    5151         else
    … …  
    5858      do j=1,jjp1
    5959         do l = 1, llm
    60            if ( (abs(dxdp(j,l)).gt.(indefini/100.)).and.
     60           if ( (abs(dxdp(j,l)).gt.(abs(indefini)/100.)).and.
    6161     .          (dxdp(j,l).ne.indefini)) then
    6262              write(*,*) '----> j= ', j , ' l= ' ,  l

  • trunk/LMDZ.VENUS/Tools/energy.F90

    r816 r1454  
    3636integer,dimension(4) :: datashape4d ! shape of 4D datasets
    3737
    38 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    39 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
     38real :: miss_val ! special "missing value" to specify missing data
     39real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
    4040real :: pi
    4141real,dimension(:),allocatable :: lon ! longitude
    4242integer lonlength ! # of grid points along longitude
    4343real,dimension(:),allocatable :: lat ! latitude
    44 real,dimension(:),allocatable :: latrad ! latitude in radian
     44real,dimension(:),allocatable :: coslat ! cos of latitude
    4545integer latlength ! # of grid points along latitude
    4646real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
    … …  
    8080
    8181pi = 2.*asin(1.)
     82miss_val = miss_val_def
    8283
    8384write(*,*) ""
    … …  
    123124if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
    124125
    125 allocate(latrad(latlength))
    126 latrad = lat*pi/180.
     126allocate(coslat(latlength))
     127! Beware of rounding problems at poles...
     128coslat(:) = max(0.,cos(lat(:)*pi/180.))
    127129
    128130! Lat, lon pressure intervals
    129 deltalat = abs(latrad(2)-latrad(1))
     131deltalat = abs(lat(2)-lat(1))*pi/180.
    130132deltalon = abs(lon(2)-lon(1))*pi/180.
    131133
    … …  
    193195  text="T"
    194196endif
    195 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     197call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    196198if (ierr1.ne.NF_NOERR) then
    197199  write(*,*) "  looking for t instead... "
    198200  text="t"
    199   call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     201  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    200202  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
    201203endif
    … …  
    216218  text="U"
    217219endif
    218 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
     220call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
    219221if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
    220222if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
    … …  
    228230  text="V"
    229231endif
    230 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
     232call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
    231233if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
    232234if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
    … …  
    242244
    243245!text="zareoid"
    244 !call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
     246!call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
    245247!if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
    246248!if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
    … …  
    276278enddo ! timelength
    277279
    278 call cellmass(infid,latlength,lonlength,altlength,timelength, &
    279               lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
     280call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
     281              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
    280282              dmass )
    281283

  • trunk/LMDZ.VENUS/Tools/epflux.F90

    r816 r1454  
    1 subroutine  epflux(iip1,jjp1,llm,indefini,rlatu,rbar &
     1subroutine  epflux(iip1,jjp1,llm,indefini,latdeg,rbar &
    22                  ,teta,u3d,v3d,w3d,press &
    33                  ,epy2d,epz2d,divep2d,vtem2d,wtem2d,acc_meridien2d &
    … …  
    4141     
    4242      integer :: iip1,jjp1,llm
    43       real :: indefini,rlatu(jjp1),rbar(jjp1,llm)
     43      real :: indefini,latdeg(jjp1),rbar(jjp1,llm)
    4444      REAL :: teta(iip1,jjp1,llm)
    4545      REAL :: u3d(iip1,jjp1,llm),v3d(iip1,jjp1,llm)
    … …  
    7575!  Autres
    7676!  ------
     77      real :: rlatu(jjp1) ! lat in radians (beware of rounding effects at poles)
    7778      real :: pi
    7879      REAL :: tetas(llm)
    … …  
    8889      jjm = jjp1-1
    8990      pi = 2.*asin(1.)
     91
     92! To avoid rounding effects at poles
     93      rlatu(:)=latdeg(:)*pi/180.00001
    9094
    9195!     Calcul de moyennes zonales
    … …  
    108112         n= 0
    109113         do j=2,jjm
    110             if (tetabar(j,l).lt.indefini) then
     114            if (tetabar(j,l).ne.indefini) then
    111115               tetas(l) = tetas(l) + tetabar(j,l)
    112116               n = n+1
    … …  
    237241      do l=1,llm
    238242
    239        if ((ubar(2,l).lt.indefini).and.(ubar(1,l).lt.indefini))then
     243       if ((ubar(2,l).ne.indefini).and.(ubar(1,l).ne.indefini))then
    240244      ducosdlat(1,l) = ( ubar(2,l)*cos(rlatu(2))-ubar(1,l)*cos(rlatu(1)) ) &
    241245                      /( rlatu(2) - rlatu(1))
    … …  
    244248       end if
    245249        do j=2,jjm
    246        if ((ubar(j+1,l).lt.indefini).and.(ubar(j-1,l).lt.indefini))then
     250       if ((ubar(j+1,l).ne.indefini).and.(ubar(j-1,l).ne.indefini))then
    247251      ducosdlat(j,l) = ( ubar(j+1,l)*cos(rlatu(j+1))-ubar(j-1,l)*cos(rlatu(j-1))) &
    248252                      /( rlatu(j+1) - rlatu(j-1))
    … …  
    251255       end if
    252256        enddo
    253        if ((ubar(jjp1,l).lt.indefini).and.(ubar(jjm,l).lt.indefini))then
     257       if ((ubar(jjp1,l).ne.indefini).and.(ubar(jjm,l).ne.indefini))then
    254258      ducosdlat(jjp1,l) = ( ubar(jjp1,l)*cos(rlatu(jjp1))  &
    255259                               -ubar(jjm,l)*cos(rlatu(jjm))  ) &
    … …  
    260264           
    261265        do j=1,jjp1
    262        if ((vptetapbar(j,l).lt.indefini).and. &
    263            (dtetabardp(j,l).lt.indefini)) then
     266       if ((vptetapbar(j,l).ne.indefini).and. &
     267           (dtetabardp(j,l).ne.indefini)) then
    264268           vz(j,l) = vptetapbar(j,l)/dtetabardp(j,l)
    265269         wlat(j,l) = cos(rlatu(j))*vz(j,l)
    … …  
    284288!     ------------
    285289
    286             if (   (rbar(j,l).lt.indefini).and. &
    287                 (vpupbar(j,l).lt.indefini).and. &
    288                 (dubardp(j,l).lt.indefini).and. &
    289                      (vz(j,l).lt.indefini)) then
     290            if (   (rbar(j,l).ne.indefini).and. &
     291                (vpupbar(j,l).ne.indefini).and. &
     292                (dubardp(j,l).ne.indefini).and. &
     293                     (vz(j,l).ne.indefini)) then
    290294           
    291295     epy2d(j,l) = rbar(j,l) * cos(rlatu(j))* ( vpupbar(j,l) &
    … …  
    303307!     ------------
    304308
    305             if (   (rbar(j,l).lt.indefini).and. &
    306                 (wpupbar(j,l).lt.indefini).and. &
    307               (ducosdlat(j,l).lt.indefini).and. &
    308                      (vz(j,l).lt.indefini)) then
     309            if (   (rbar(j,l).ne.indefini).and. &
     310                (wpupbar(j,l).ne.indefini).and. &
     311              (ducosdlat(j,l).ne.indefini).and. &
     312                     (vz(j,l).ne.indefini)) then
    309313
    310314     epz2d(j,l) = rbar(j,l) * cos(rlatu(j))*                    &
    … …  
    330334      do l=1,llm
    331335
    332           if ( (wlat(2,l).lt.indefini) .and. &
    333                (wlat(1,l).lt.indefini) ) then
     336          if ( (wlat(2,l).ne.indefini) .and. &
     337               (wlat(1,l).ne.indefini) ) then
    334338               dwlatdlat(1,l)=(wlat(2,l)-wlat(1,l)) &
    335339                     /(rlatu(2)-rlatu(1))
    … …  
    338342          end if
    339343        do j=2,jjm
    340           if ( (wlat(j+1,l).lt.indefini) .and. &
    341                (wlat(j-1,l).lt.indefini) ) then
     344          if ( (wlat(j+1,l).ne.indefini) .and. &
     345               (wlat(j-1,l).ne.indefini) ) then
    342346               dwlatdlat(j,l)=(wlat(j+1,l)-wlat(j-1,l)) &
    343347                     /(rlatu(j+1)-rlatu(j-1))
    … …  
    346350          end if
    347351        enddo
    348           if ( (wlat(jjp1,l).lt.indefini) .and. &
    349                (wlat(jjm, l).lt.indefini) ) then
     352          if ( (wlat(jjp1,l).ne.indefini) .and. &
     353               (wlat(jjm, l).ne.indefini) ) then
    350354               dwlatdlat(jjp1,l)=(wlat(jjp1,l)-wlat(jjm,l)) &
    351355                     /(rlatu(jjp1)-rlatu(jjm))
    … …  
    355359
    356360        do j=1,jjp1
    357           if ( (dvzdp(j,l).lt.indefini) &
    358            .and.(vbar(j,l).lt.indefini)) then
     361          if ( (dvzdp(j,l).ne.indefini) &
     362           .and.(vbar(j,l).ne.indefini)) then
    359363            vtem2d(j,l) = vbar(j,l)-dvzdp(j,l)
    360364          else
    … …  
    364368
    365369        do j=1,jjp1
    366           if (   (rbar(j,l).lt.indefini).and. &
    367             (dwlatdlat(j,l).lt.indefini).and. &
    368                  (wbar(j,l).lt.indefini)) then
     370          if (   (rbar(j,l).ne.indefini).and. &
     371            (dwlatdlat(j,l).ne.indefini).and. &
     372                 (wbar(j,l).ne.indefini)) then
    369373           wtem2d(j,l) = wbar(j,l)+dwlatdlat(j,l)/(rbar(j,l)*cos(rlatu(j)))
    370374          else
    … …  
    385389      do l=1,llm
    386390        do j=1,jjp1
    387           if (   (rbar(j,l).lt.indefini) &
    388           .and.(vtem2d(j,l).lt.indefini) &
    389           .and.(wtem2d(j,l).lt.indefini) &
    390        .and.(ducosdlat(j,l).lt.indefini) &
    391          .and.(dubardp(j,l).lt.indefini)) then
     391          if (   (rbar(j,l).ne.indefini) &
     392          .and.(vtem2d(j,l).ne.indefini) &
     393          .and.(wtem2d(j,l).ne.indefini) &
     394       .and.(ducosdlat(j,l).ne.indefini) &
     395         .and.(dubardp(j,l).ne.indefini)) then
    392396       acc_meridien2d(j,l) =  &
    393397           vtem2d(j,l)*(ducosdlat(j,l)/(rbar(j,l)*cos(rlatu(j))) &
    … …  
    410414        divep2d(1,l) =0
    411415        do j=2,jjm
    412           if (    (rbar(j,l).lt.indefini) .and. &
    413                (epy2d(j+1,l).lt.indefini) .and. &
    414                (epy2d(j-1,l).lt.indefini) .and. &
    415                 (depzdp(j,l).lt.indefini)    ) then
     416          if (    (rbar(j,l).ne.indefini) .and. &
     417               (epy2d(j+1,l).ne.indefini) .and. &
     418               (epy2d(j-1,l).ne.indefini) .and. &
     419                (depzdp(j,l).ne.indefini)    ) then
    416420
    417421           depycosdlat=(epy2d(j+1,l)*cos(rlatu(j+1))   &

  • trunk/LMDZ.VENUS/Tools/fft.F90

    r888 r1454  
    6363integer,dimension(4) :: datashape4d ! shape of 4D datasets
    6464
    65 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    66 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
     65real :: miss_val ! special "missing value" to specify missing data
     66real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
    6767real :: pi
    6868real,dimension(:),allocatable :: lon ! longitude
    6969integer lonlength ! # of grid points along longitude
    7070real,dimension(:),allocatable :: lat ! latitude
    71 real,dimension(:),allocatable :: latrad ! latitude in radian
    7271integer latlength ! # of grid points along latitude
    7372real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
    … …  
    125124logical :: lmdflag
    126125
     126! Tuning parameters
     127real :: fcoup1,fcoup2,width
     128real :: fcoup1tmp,fcoup2tmp,widthtmp
     129logical,dimension(4) :: ok_out
     130character (len=1) :: ok_outtmp
     131
    127132include "planet.h"
    128 include "filter.h"
    129133
    130134#include <fftw3.f>
    … …  
    135139
    136140pi = 2.*asin(1.)
     141miss_val = miss_val_def
    137142
    138143write(*,*) ""
    139144write(*,*) "You are working on the atmosphere of ",planet
     145
     146! initialisation
     147!----------------
     148
     149! Par defaut
     150
     151! Define the filters
     152! Low  cutting frequency, in Hz    : fcoup1
     153fcoup1=6.e-7
     154! High cutting frequency, in Hz    : fcoup2
     155fcoup2=1.4e-6
     156! Half-width of the filters, in Hz : width
     157width=2.e-7
     158! Outputs (U,     V,      W,     T)
     159ok_out=(/.true.,.true.,.false.,.true./)
     160
     161print*,"Low  cutting frequency, in Hz ?"
     162print*,"between 1e-5 and 1e-7, 0 for default => 6.e-7"
     163read(*,*) fcoup1tmp
     164if ((fcoup1tmp.lt.1e-5).and.(fcoup1tmp.gt.1e-7)) fcoup1=fcoup1tmp
     165print*,"=",fcoup1
     166print*,"High cutting frequency, in Hz ?"
     167print*,"between 1e-5 and 1e-7, 0 for default => 1.4e-6"
     168read(*,*) fcoup2tmp
     169if ((fcoup2tmp.lt.1e-5).and.(fcoup2tmp.gt.1e-7)) fcoup2=fcoup2tmp
     170print*,"=",fcoup2
     171print*,"Half-width of the filters, in Hz ?"
     172print*,"between 1e-6 and 1e-8, 0 for default => 2e-7)"
     173read(*,*) widthtmp
     174if ((widthtmp.lt.1e-6).and.(widthtmp.gt.1e-8)) width=widthtmp
     175print*,"=",width
     176!width = 3./time(timelength)
     177
     178! Outputs
     179print*,"Output of zonal wind ? (y or n, default is y)"
     180read(*,'(a1)') ok_outtmp
     181if (ok_outtmp.eq."n") ok_out(1)=.false.
     182print*,"=",ok_out(1)
     183print*,"Output of meridional wind ? (y or n, default is y)"
     184read(*,'(a1)') ok_outtmp
     185if (ok_outtmp.eq."n") ok_out(2)=.false.
     186print*,"=",ok_out(2)
     187print*,"Output of vertical wind ? (y or n, default is n)"
     188read(*,'(a1)') ok_outtmp
     189if (ok_outtmp.eq."y") ok_out(3)=.true.
     190print*,"=",ok_out(3)
     191print*,"Output of temperature ? (y or n, default is y)"
     192read(*,'(a1)') ok_outtmp
     193if (ok_outtmp.eq."n") ok_out(4)=.false.
     194print*,"=",ok_out(4)
    140195
    141196!===============================================================================
    … …  
    172227ierr=NF_GET_VAR_REAL(infid,lat_varid,lat)
    173228if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
    174 
    175 allocate(latrad(latlength))
    176 latrad = lat*pi/180.
    177229
    178230allocate(plev(altlength))
    … …  
    212264
    213265text="temp"
    214 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     266call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    215267if (ierr1.ne.NF_NOERR) then
    216268  write(*,*) "  looking for t instead... "
    217269  text="t"
    218   call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     270  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    219271  if (ierr1.ne.NF_NOERR) then
    220272    print*,"Error: Failed to get temperature ID"
    … …  
    236288
    237289text="vitu"
    238 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
     290call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
    239291if (ierr1.ne.NF_NOERR) then
    240292  print*,"Error: Failed to get vitu ID"
    … …  
    252304
    253305text="vitv"
    254 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
     306call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
    255307if (ierr1.ne.NF_NOERR) then
    256308  print*,"Error: Failed to get vitv ID"
    … …  
    268320
    269321text="vitw"
    270 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
     322call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
    271323if (ierr1.ne.NF_NOERR) then
    272324  print*,"Error: Failed to get vitw ID"
    … …  
    288340! 2.2.1 FFT and filtering
    289341!===============================================================================
     342
    290343! allocations
    291344!-------------
    … …  
    405458enddo
    406459
    407 ! Define the filters
    408 
    409 ! IN filter.h :
    410 ! Low  cutting frequency, in Hz    : fcoup1
    411 ! High cutting frequency, in Hz    : fcoup2
    412 ! Half-width of the filters, in Hz : width
    413 
    414 !print*,"Low  cutting frequency, in Hz ?"
    415 !read(*,*) fcoup1
    416 !print*,"High cutting frequency, in Hz ?"
    417 !read(*,*) fcoup2
    418 !print*,"Half-width of the filters, in Hz ?"
    419 !read(*,*) width
    420 !width = 3./time(timelength)
    421 
    422460do itim=1,M_fft+1
    423461  if (freq(itim).lt.(fcoup1-width)) then

  • trunk/LMDZ.VENUS/Tools/io.F90

    r816 r1454  
    218218!===========================================================================
    219219
    220 subroutine get_var4d(infid,dim1,dim2,dim3,dim4,text,var,ierr1,ierr2)
     220subroutine get_var4d(infid,dim1,dim2,dim3,dim4,text,var,missing,ierr1,ierr2)
    221221
    222222implicit none
    … …  
    229229character (len=64) :: text ! name of variable to read
    230230real,dimension(dim1,dim2,dim3,dim4) :: var ! variable to read
    231 integer :: ierr1,ierr2 ! NetCDF routines return code
     231real :: missing ! missing value
     232integer :: ierr1,ierr2,miss ! NetCDF routines return code
    232233
    233234! local
    … …  
    240241  write(*,*) "ID ok for ",trim(text)
    241242  ierr2=NF_GET_VAR_REAL(infid,tmpvarid,var)
     243  miss=NF_GET_ATT_REAL(infid,tmpvarid,"missing_value",missing)
    242244endif
    243245

  • trunk/LMDZ.VENUS/Tools/moytim.F

    r816 r1454  
    3333           n = 0
    3434           do t=1,nt
    35               if (x(i,j,l,t).lt.indefini) then
     35              if (x(i,j,l,t).ne.indefini) then
    3636                 xmean(i,j,l) = xmean(i,j,l) + x(i,j,l,t)
    3737                 n = n+1

  • trunk/LMDZ.VENUS/Tools/moyzon.F

    r816 r1454  
    3232           n = 0
    3333           do i=1,iim
    34               if (x(i,j,l).lt.indefini) then
     34              if (x(i,j,l).ne.indefini) then
    3535                 xbar(j,l) = xbar(j,l) + x(i,j,l)
    3636                 n = n+1

  • trunk/LMDZ.VENUS/Tools/moyzon2.F

    r816 r1454  
    3434           n = 0
    3535           do i=1,iim
    36               if (x(i,j,l).lt.indefini) then
     36              if (x(i,j,l).ne.indefini) then
    3737                 xbar(j,l) = xbar(j,l) + x(i,j,l)
    3838                 n = n+1

  • trunk/LMDZ.VENUS/Tools/psi.F90

    r980 r1454  
    3131integer,dimension(4) :: datashape4d ! shape of 4D datasets
    3232
    33 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    34 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
     33real :: miss_val ! special "missing value" to specify missing data
     34real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
    3535real :: pi
    3636real,dimension(:),allocatable :: lon ! longitude
    3737integer lonlength ! # of grid points along longitude
    3838real,dimension(:),allocatable :: lat ! latitude
    39 real,dimension(:),allocatable :: latrad ! latitude in radian
     39real,dimension(:),allocatable :: coslat ! cos of latitude
    4040integer latlength ! # of grid points along latitude
    4141real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
    … …  
    6666
    6767pi = 2.*asin(1.)
     68miss_val = miss_val_def
    6869
    6970write(*,*) ""
    … …  
    104105if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
    105106
    106 allocate(latrad(latlength))
    107 latrad = lat*pi/180.
     107allocate(coslat(latlength))
     108! Beware of rounding problems at poles...
     109coslat(:) = max(0.,cos(lat(:)*pi/180.))
    108110
    109111! Lat, lon pressure intervals
    110 deltalat = abs(latrad(2)-latrad(1))
     112deltalat = abs(lat(2)-lat(1))*pi/180.
    111113deltalon = abs(lon(2)-lon(1))*pi/180.
    112114
    … …  
    160162! meridional wind vitv (in m/s)
    161163text="vitv"
    162 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
     164call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
    163165if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
    164166if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
    … …  
    173175
    174176text="zareoid"
    175 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
     177call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
    176178if (ierr1.ne.NF_NOERR) then
    177179  write(*,*) "  looking for geop instead... "
    178180  text="geop"
    179   call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
     181  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
    180182  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get geop ID"
    181183do itim=1,timelength
    … …  
    226228
    227229call cellmass(infid,latlength,lonlength,altlength,timelength, &
    228               lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
     230              lmdflag,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
    229231              dmass )
    230232

  • trunk/LMDZ.VENUS/Tools/stability.F90

    r816 r1454  
    3131integer,dimension(4) :: datashape4d ! shape of 4D datasets
    3232
    33 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    34 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
     33real :: miss_val ! special "missing value" to specify missing data
     34real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
    3535real :: pi
    3636real,dimension(:),allocatable :: lon ! longitude
    3737integer lonlength ! # of grid points along longitude
    3838real,dimension(:),allocatable :: lat ! latitude
    39 real,dimension(:),allocatable :: latrad ! latitude in radian
     39real,dimension(:),allocatable :: coslat ! cos of latitude
    4040integer latlength ! # of grid points along latitude
    4141real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
    … …  
    9191
    9292pi = 2.*asin(1.)
     93miss_val = miss_val_def
    9394
    9495write(*,*) ""
    … …  
    129130if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
    130131
    131 allocate(latrad(latlength))
    132 latrad = lat*pi/180.
     132allocate(coslat(latlength))
     133! Beware of rounding problems at poles...
     134coslat(:) = max(0.,cos(lat(:)*pi/180.))
    133135
    134136! Lat, lon pressure intervals
    135 deltalat = abs(latrad(2)-latrad(1))
     137deltalat = abs(lat(2)-lat(1))*pi/180.
    136138deltalon = abs(lon(2)-lon(1))*pi/180.
    137139
    … …  
    180182
    181183text="temp"
    182 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     184call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    183185if (ierr1.ne.NF_NOERR) then
    184186  write(*,*) "  looking for t instead... "
    185187  text="t"
    186   call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     188  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    187189  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
    188190endif
    … …  
    197199! zonal wind vitu (in m/s)
    198200text="vitu"
    199 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
     201call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
    200202if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
    201203if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
    … …  
    203205! meridional wind vitv (in m/s)
    204206text="vitv"
    205 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
     207call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
    206208if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
    207209if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
    … …  
    215217
    216218! text="zareoid"
    217 ! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
     219! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
    218220! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
    219221! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
    … …  
    276278 grav(lonlength+1,:,:,:) =  grav(1,:,:,:)
    277279
    278 call cellmass(infid,latlength,lonlength+1,altlength,timelength, &
    279               lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
     280call cellmass(infid,latlength,lonlength+1,altlength,timelength,lmdflag, &
     281              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
    280282              dmass )
    281283
    … …  
    364366
    365367do ilat=2,latlength-1
    366    if (tan(latrad(ilat)).ne.0.) then
    367      fac1 = R0/tan(latrad(ilat))
     368   if (tan(lat(ilat)*pi/180.).ne.0.) then
     369     fac1 = R0/tan(lat(ilat)*pi/180.)
    368370   else
    369371     fac1 = miss_val

  • trunk/LMDZ.VENUS/Tools/tem.F90

    r816 r1454  
    3737integer,dimension(3) :: datashape3d ! shape of 3D datasets
    3838
    39 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    40 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
     39real :: miss_val ! special "missing value" to specify missing data
     40real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
    4141real :: pi
    4242real,dimension(:),allocatable :: lon ! longitude
    4343integer lonlength ! # of grid points along longitude
    4444real,dimension(:),allocatable :: lat ! latitude
    45 real,dimension(:),allocatable :: latrad ! latitude in radian
     45real,dimension(:),allocatable :: coslat ! cos of latitude
    4646integer latlength ! # of grid points along latitude
    4747real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
    … …  
    107107
    108108pi = 2.*asin(1.)
     109miss_val = miss_val_def
    109110
    110111write(*,*) ""
    … …  
    145146if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
    146147
    147 allocate(latrad(latlength))
    148 latrad = lat*pi/180.
     148allocate(coslat(latlength))
     149! Beware of rounding problems at poles...
     150coslat(:) = max(0.,cos(lat(:)*pi/180.))
    149151
    150152! Lat, lon pressure intervals
    151 deltalat = abs(latrad(2)-latrad(1))
     153deltalat = abs(lat(2)-lat(1))*pi/180.
    152154deltalon = abs(lon(2)-lon(1))*pi/180.
    153155
    … …  
    196198
    197199text="temp"
    198 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     200call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    199201if (ierr1.ne.NF_NOERR) then
    200202  write(*,*) "  looking for t instead... "
    201203  text="t"
    202   call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
     204  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
    203205  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
    204206endif
    … …  
    214216! zonal wind vitu (in m/s)
    215217text="vitu"
    216 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
     218call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
    217219if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
    218220if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
    … …  
    220222! meridional wind vitv (in m/s)
    221223text="vitv"
    222 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
     224call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
    223225if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
    224226if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
    … …  
    226228! vertical wind vitw (in Pa/s)
    227229text="vitw"
    228 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
     230call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
    229231if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitw ID"
    230232if (ierr2.ne.NF_NOERR) stop "Error: Failed reading vertical wind"
    … …  
    238240
    239241! text="zareoid"
    240 ! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
     242! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
    241243! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
    242244! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
    … …  
    312314 grav(lonlength+1,:,:,:) =  grav(1,:,:,:)
    313315
    314 call cellmass(infid,latlength,lonlength+1,altlength,timelength, &
    315               lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
     316call cellmass(infid,latlength,lonlength+1,altlength,timelength,lmdflag, &
     317              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
    316318              dmass )
    317319
    … …  
    359361call moyzon(lonlength,latlength,altlength,miss_val,u3d,ubar)
    360362
    361 call epflux(lonlength+1,latlength,altlength,miss_val,latrad,rbar &
     363call epflux(lonlength+1,latlength,altlength,miss_val,lat,rbar &
    362364            ,teta,u3d,v3d,w3d,plev &
    363365            ,epy2d,epz2d,div2d,vtem2d,wtem2d,ammc2d &
    … …  
    376378! rsurg: r*dp/g = dm/(r cos(lat) dlon dlat) !!!
    377379     rsurg(ilon,ilat,ilev) = dmass(ilon,ilat,ilev,itim) &
    378           / (r3d(ilon,ilat,ilev)*cos(latrad(ilat))*deltalon*deltalat)
     380          / (r3d(ilon,ilat,ilev)*coslat(ilat)*deltalon*deltalat)
    379381    else
    380382     rsurg(ilon,ilat,ilev) = miss_val
    … …  
    391393        (rsurgbar(ilat,ilev).ne.miss_val) ) then
    392394      vm(ilat,ilev) = vtem2d(ilat,ilev) &
    393             * 2.*pi*rsurgbar(ilat,ilev)*cos(latrad(ilat))
     395            * 2.*pi*rsurgbar(ilat,ilev)*coslat(ilat)
    394396    else
    395397      vm(ilat,ilev) = miss_val

  • trunk/LMDZ.VENUS/Tools/tmc.F90

    r816 r1454  
    4949integer,dimension(4) :: datashape4d ! shape of 4D datasets
    5050
    51 real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
    52 real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
     51real :: miss_val ! special "missing value" to specify missing data
     52real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
    5353real :: pi
    5454real,dimension(:),allocatable :: lon ! longitude
    5555integer lonlength ! # of grid points along longitude
    5656real,dimension(:),allocatable :: lat ! latitude
    57 real,dimension(:),allocatable :: latrad ! latitude in radian
     57real,dimension(:),allocatable :: coslat ! cos of latitude
    5858integer latlength ! # of grid points along latitude
    5959real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
    … …  
    141141
    142142pi = 2.*asin(1.)
     143miss_val = miss_val_def
    143144
    144145write(*,*) ""
    … …  
    179180if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
    180181
    181 allocate(latrad(latlength))
    182 latrad = lat*pi/180.
     182allocate(coslat(latlength))
     183! Beware of rounding problems at poles...
     184coslat(:) = max(0.,cos(lat(:)*pi/180.))
    183185
    184186! Lat, lon pressure intervals
    185 deltalat = abs(latrad(2)-latrad(1))
     187deltalat = abs(lat(2)-lat(1))*pi/180.
    186188deltalon = abs(lon(2)-lon(1))*pi/180.
    187189
    … …  
    233235! zonal wind vitu (in m/s)
    234236text="vitu"
    235 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
     237call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
    236238if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
    237239if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
    … …  
    239241! meridional wind vitv (in m/s)
    240242text="vitv"
    241 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
     243call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
    242244if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
    243245if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
    … …  
    245247! vertical wind vitw (in Pa/s)
    246248text="vitw"
    247 call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
     249call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
    248250if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitw ID"
    249251if (ierr2.ne.NF_NOERR) stop "Error: Failed reading vertical wind"
    … …  
    257259
    258260! text="zareoid"
    259 ! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
     261! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
    260262! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
    261263! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
    … …  
    293295enddo ! timelength
    294296
    295 call cellmass(infid,latlength,lonlength,altlength,timelength, &
    296               lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
     297call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
     298              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
    297299              dmass )
    298300
    … …  
    316318      osam(ilon,ilat,ilev,itim) = &
    317319         rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim) &
    318        * cos(latrad(ilat))*cos(latrad(ilat)) &
     320       * coslat(ilat)*coslat(ilat) &
    319321       * omega
    320322      rsam(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim) &
    321        * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))
     323       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)
    322324      tsam(ilon,ilat,ilev,itim) = osam(ilon,ilat,ilev,itim)&
    323325                                + rsam(ilon,ilat,ilev,itim)
Note: See TracChangeset for help on using the changeset viewer.