Ignore:
Timestamp:
Jul 7, 2020, 3:49:34 PM (4 years ago)
Author:
emillour
Message:

Mars GCM:
Improve "zrecast" utility: add possibility to transfer 3D (lon-lat-time)
variables to output file and add some extra sanity checks of user inputs.
EM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/util/zrecast.F90

    r2140 r2401  
    7979integer nbvarinfile ! # of variables in input file
    8080integer nbattr ! # of attributes of a given variable in input file
     81integer nbvar3dinfile ! # of 3D (lon,lat,time) variables in input file
    8182integer nbvar4dinfile ! # of 4D (lon,lat,alt,time) variables in input file
    8283integer outfid ! NetCDF output file ID
    … …  
    9798character (len=64), dimension(:), allocatable :: var
    9899! var(): names of variables that will be processed
     100integer, allocatable :: varndims(:) ! number of dimensions of var(:) variables
    99101integer nbvar ! # of variables to process
    100102integer,dimension(:),allocatable :: var_id ! IDs of variables var() (in outfile)
    … …  
    120122real,dimension(:,:,:,:),allocatable :: ra_gcm ! radial distance (m) to
    121123                                          ! center of Mars, (of GCM grid points)
     124real,dimension(:,:,:),allocatable :: surfdata ! to store a GCM 3D dataset
    122125real,dimension(:,:,:,:),allocatable :: indata ! to store a GCM 4D dataset
    123126real,dimension(:,:,:,:),allocatable :: outdata ! to store a 4D dataset
    … …  
    186189write(*,*)" The following variables have been found:"
    187190nbvar4dinfile=0
     191nbvar3dinfile=0
    188192do i=1,nbvarinfile
    189193  ! get name of variable # i
    … …  
    195199    write(*,*) trim(tmpvarname)
    196200  endif
     201  if (tmpndims.eq.3) then
     202    nbvar3dinfile=nbvar3dinfile+1
     203    write(*,*) trim(tmpvarname)
     204  endif
    197205enddo
    198206
    199 allocate(var(nbvar4dinfile),stat=ierr)
     207allocate(var(nbvar4dinfile+nbvar3dinfile),stat=ierr)
    200208if (ierr.ne.0) then
    201   write(*,*) "Error! failed memory allocation of var(nbvar4dinfile)"
     209  write(*,*) "Error! failed memory allocation of var(nbvar4dinfile+nbvar3dinfile)"
     210  stop
     211endif
     212allocate(varndims(nbvar4dinfile+nbvar3dinfile),stat=ierr)
     213if (ierr.ne.0) then
     214  write(*,*) "Error! failed memory allocation of varndims(nbvar4dinfile+nbvar3dinfile)"
    202215  stop
    203216endif
    … …  
    213226  ierr=NF_INQ_VARID(infid,tmpvarname,tmpvarid)
    214227  if (ierr.eq.NF_NOERR) then ! valid name
    215     ! also check that it is indeed a 4D variable
     228    ! also check that it is indeed a 3D or 4D variable
    216229    ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims)
    217     if (tmpndims.eq.4) then
     230    if ((tmpndims.eq.4).or.(tmpndims.eq.3)) then
    218231      nbvar=nbvar+1
    219       var(nbvar)=tmpvarname   
     232      var(nbvar)=tmpvarname
     233      varndims(nbvar)=tmpndims
    220234    else
    221       write(*,*) 'Error: ',trim(tmpvarname),' is not a 4D variable'
     235      write(*,*) 'Error: ',trim(tmpvarname),' is not a 3D or 4D variable'
    222236      write(*,*) '       (we''ll skip that one)'
    223237    endif
    … …  
    233247  nbvar=0
    234248  do i=1,nbvarinfile
    235     ! look for 4D variables
     249    ! look for 3D and 4D variables
    236250    ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
    237     if (tmpndims.eq.4) then
     251    if ((tmpndims.eq.4).or.(tmpndims.eq.3)) then
    238252      nbvar=nbvar+1
    239253      ! get the corresponding name
    240254      ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
    241255      var(nbvar)=tmpvarname
     256      varndims(nbvar)=tmpndims
    242257    endif
    243258  enddo
    … …  
    622637  if (.not.auto_mcd_levels) then
    623638    write(*,*) ""
    624     write(*,*) "Enter min and max of vertical coordinate (Pa or m)"
    625     write(*,*) " (in that order and on the same line)"
    626639    if (ztype.eq.1) then ! pressure coordinate
    627       read(*,*) pmin,pmax
     640      write(*,*) "Enter min and max of vertical coordinate (Pa)"
     641      write(*,*) " (in that order and on the same line)"
     642      read(*,*,iostat=ierr) pmin,pmax
     643      if ((ierr/=0).or.(pmin.gt.pmax)) then
     644        write(*,*) "Error: wrong vertical coordinate boundaries"
     645        write(*,*) "       pmin=",pmin
     646        write(*,*) "       pmax=",pmax
     647        stop
     648      endif
    628649    else ! altitude coordinate, except in MCD case
    629       read(*,*) zamin,zamax
     650      write(*,*) "Enter min and max of vertical coordinate (m)"
     651      write(*,*) " (in that order and on the same line)"
     652      read(*,*,iostat=ierr) zamin,zamax
     653      if ((ierr/=0).or.(zamin.gt.zamax)) then
     654        write(*,*) "Error: wrong vertical coordinate boundaries"
     655        write(*,*) "       zamin=",zamin
     656        write(*,*) "       zamax=",zamax
     657        stop
     658      endif
    630659    endif
    631660  endif
    … …  
    634663  if (ztype.eq.1) then ! pressure coordinate
    635664    write(*,*) "Number of levels along vertical coordinate?"
    636     read(*,*) nblev
     665    read(*,*,iostat=ierr) nblev
     666    if (ierr/=0) then
     667      write(*,*) "Error: bad input value for nblev ", nblev
     668      stop
     669    endif
    637670    allocate(plevel(nblev),stat=ierr)
    638671    if (ierr.ne.0) then
    … …  
    652685  else if (ztype.eq.2) then ! above areoid heights
    653686    write(*,*) "Number of levels along vertical coordinate?"
    654     read(*,*) nblev
     687    read(*,*,iostat=ierr) nblev
     688    if (ierr/=0) then
     689      write(*,*) "Error: bad input value for nblev ", nblev
     690      stop
     691    endif
    655692    allocate(zareoid(nblev),stat=ierr)
    656693    if (ierr.ne.0) then
    … …  
    680717    else
    681718      write(*,*) "Number of levels along vertical coordinate?"
    682       read(*,*) nblev
     719      read(*,*,iostat=ierr) nblev
     720      if (ierr/=0) then
     721        write(*,*) "Error: bad input value for nblev ", nblev
     722        stop
     723      endif
    683724      allocate(zsurface(nblev),stat=ierr)
    684725      if (ierr.ne.0) then
    … …  
    697738  else ! distance to center of planet
    698739   write(*,*) "Number of levels along vertical coordinate?"
    699     read(*,*) nblev
     740    read(*,*,iostat=ierr) nblev
     741    if (ierr/=0) then
     742      write(*,*) "Error: bad input value for nblev ", nblev
     743      stop
     744    endif
    700745    allocate(zradius(nblev),stat=ierr)
    701746    if (ierr.ne.0) then
    … …  
    716761  write(*,*) ""
    717762  write(*,*) "Number of levels along vertical coordinate?"
    718   read(*,*) nblev
     763  read(*,*,iostat=ierr) nblev
     764  if (ierr/=0) then
     765    write(*,*) "Error: bad input value for nblev ", nblev
     766    stop
     767  endif
    719768  if (ztype.eq.1) then ! pressure coordinate
    720769    allocate(plevel(nblev),stat=ierr)
    … …  
    728777    write(*,*) " (one value per line)"
    729778    do i=1,nblev
    730       read(*,*) plevel(i)
     779      read(*,*,iostat=ierr) plevel(i)
     780      if (ierr/=0) then
     781        write(*,*) "Error: bad input value for plevel(i):",plevel(i)
     782        stop
     783      endif
    731784    enddo
    732785  else if (ztype.eq.2) then ! above areoid heights
    … …  
    740793    write(*,*) " from min to max (one value per line)"
    741794    do i=1,nblev
    742       read(*,*) zareoid(i)
     795      read(*,*,iostat=ierr) zareoid(i)
     796      if (ierr/=0) then
     797        write(*,*) "Error: bad input value for zareoid(i):",zareoid(i)
     798        stop
     799      endif
    743800    enddo
    744801  else if (ztype.eq.3) then ! above local surface
    … …  
    752809    write(*,*) " from min to max (one value per line)"
    753810    do i=1,nblev
    754       read(*,*) zsurface(i)
     811      read(*,*,iostat=ierr) zsurface(i)
     812      if (ierr/=0) then
     813        write(*,*) "Error: bad input value for zsurface(i):",zsurface(i)
     814        stop
     815      endif
    755816    enddo
    756817  else ! distance to center of planet
    … …  
    764825    write(*,*) " from min to max (one value per line)"
    765826    do i=1,nblev
    766       read(*,*) zradius(i)
     827      read(*,*,iostat=ierr) zradius(i)
     828      if (ierr/=0) then
     829        write(*,*) "Error: bad input value for zradius(i):",zradius(i)
     830        stop
     831      endif
    767832    enddo
    768833  endif
    … …  
    13171382endif
    13181383
    1319 ! 3.3.3 Define 4D variables
     1384! 3.3.3 Define 3D and 4D variables
    13201385
    13211386! add pressure or zareoid
    … …  
    13811446  stop
    13821447endif
    1383 do i=1,nbvar
     1448do i=1,nbvar ! loop on 3D+4D variables
     1449  if (var(i)=="ps") cycle   ! skip surface pressure (already processed)
     1450  ! other variables
    13841451  write(*,*) ""
    13851452  write(*,*) "Creating ",trim(var(i))
    1386   ! define the variable
    1387   ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,4,datashape,var_id(i))
    1388   if (ierr.ne.NF_NOERR) then
    1389     write(*,*) 'Error, could not define variable ',trim(var(i))
    1390     write(*,*) NF_STRERROR(ierr)
    1391     stop
     1453  if (varndims(i).eq.3) then
     1454    ! define the variable
     1455    ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,3,surfdatashape,var_id(i))
     1456    if (ierr.ne.NF_NOERR) then
     1457      write(*,*) 'Error, could not define variable ',trim(var(i))
     1458      write(*,*) NF_STRERROR(ierr)
     1459      stop
     1460    endif
     1461  elseif (varndims(i).eq.4) then
     1462    ! define the variable
     1463    ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,4,datashape,var_id(i))
     1464    if (ierr.ne.NF_NOERR) then
     1465      write(*,*) 'Error, could not define variable ',trim(var(i))
     1466      write(*,*) NF_STRERROR(ierr)
     1467      stop
     1468    endif
     1469  else
     1470    write(*,*) "Unexpected value for varndims(i)!"
     1471    write(*,*) " i=",i," varndims(i)=",varndims(i)
     1472    stop
    13921473  endif
    13931474
    … …  
    16051686endif
    16061687
     1688! allocate surfdata() to store input 3D GCM data
     1689allocate(surfdata(lonlength,latlength,timelength),stat=ierr)
     1690if (ierr.ne.0) then
     1691  write(*,*) "Error: Failed to allocate surfdata(lonlength,latlength,timelength)"
     1692  write(*,*) "       lonlength=",lonlength," latlength=",latlength
     1693  write(*,*) "       timelength=",timelength
     1694  stop
     1695endif
     1696
     1697do i=1,nbvar ! loop on 3D+4D variables
     1698  if (varndims(i).eq.4) cycle ! skip 4D variables
     1699  if (var(i)=="ps") cycle ! skip surface pressure (already processed)
     1700  ! load the variable from input file
     1701  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
     1702  if (ierr.ne.NF_NOERR) then
     1703    write(*,*) "Error: Failed to get ",trim(var(i)),"ID"
     1704    write(*,*) NF_STRERROR(ierr)
     1705    stop
     1706  endif
     1707  ierr=NF_GET_VAR_REAL(infid,tmpvarid,surfdata)
     1708  if (ierr.ne.NF_NOERR) then
     1709    write(*,*) "Error: Failed to load ",trim(var(i))
     1710    write(*,*) NF_STRERROR(ierr)
     1711    stop
     1712  endif
     1713  ! write the variable to output file
     1714  ierr=NF_PUT_VARA_REAL(outfid,var_id(i),corner(1:3),edges(1:3),surfdata)
     1715  if (ierr.ne.NF_NOERR) then
     1716    write(*,*) "Error: Could not write ",trim(var(i))," data to output file"
     1717    write(*,*) NF_STRERROR(ierr)
     1718    stop
     1719  endif
     1720enddo ! of do i=1,nbvar
     1721
    16071722!===============================================================================
    16081723! 4. Interpolate and write 4D variables
    … …  
    16291744! 4.1 If output is in pressure coordinate
    16301745if (ztype.eq.1) then
    1631   do i=1,nbvar ! loop on 4D variable to process
     1746  do i=1,nbvar ! loop on 3D+4D variable to process
     1747    ! Skip 3D variables
     1748    if (varndims(i).eq.3) cycle
    16321749  ! identify and read a dataset
    16331750  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
    … …  
    16751792! 4.2 If output is in above areoid altitude
    16761793if (ztype.eq.2) then
    1677   do i=1,nbvar ! loop on 4D variable to process
     1794  do i=1,nbvar ! loop on 3D+4D variable to process
     1795    ! Skip 3D variables
     1796    if (varndims(i).eq.3) cycle
    16781797  ! identify and read a dataset
    16791798  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
    … …  
    17261845! 4.3 If output is in above local surface altitude
    17271846if (ztype.eq.3) then
    1728   do i=1,nbvar ! loop on 4D variable to process
     1847  do i=1,nbvar ! loop on 3D+4D variable to process
     1848    ! Skip 3D variables
     1849    if (varndims(i).eq.3) cycle
    17291850    write(*,*) " "
    17301851    write(*,*) "Processing "//trim(var(i))
    … …  
    17811902! 4.4 If output is in radial distance to center of planet
    17821903if (ztype.eq.4) then
    1783   do i=1,nbvar ! loop on 4D variable to process
     1904  do i=1,nbvar ! loop on 3D+4D variable to process
     1905    ! Skip 3D variables
     1906    if (varndims(i).eq.3) cycle
    17841907  ! identify and read a dataset
    17851908  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
    … …  
    18361959  write(*,*) NF_STRERROR(ierr)
    18371960endif
     1961
     1962write(*,*)
     1963write(*,*) "zrecast: All is well that ends well."
    18381964
    18391965end program
    … …  
    28462972! data structure of gravity field
    28472973        double precision v0,omega,ae,gm
    2848         double precision clm(0:ndeg,0:ndeg),slm(0:ndeg,0:ndeg)
    2849         common /gmm1/v0,omega,ae,gm,clm,slm
     2974       double precision clm(0:ndeg,0:ndeg),slm(0:ndeg,0:ndeg)
     2975       common /gmm1/v0,omega,ae,gm,clm,slm
    28502976        integer lmin,lmax
    28512977        double precision root(nd2p3)
    28522978        double precision requator
    2853         common /sqr/ lmin,lmax,root,requator
     2979       common /sqr/ lmin,lmax,root,requator
    28542980
    28552981!===============================================================================
    … …  
    28652991   dlat=lat(jloop) ! double precision version of lat
    28662992   call geoid(dlon,dlat,dareoid)
    2867    rareoid(iloop,jloop)=dareoid
     2993   rareoid(iloop,jloop)=real(dareoid)
    28682994 enddo
    28692995enddo
    … …  
    29013027! data structure of gravity field
    29023028         integer ndeg,ndeg2,nd2p3
    2903           parameter (ndeg=90,ndeg2=2*ndeg,nd2p3=ndeg2+3)
    2904           double precision c(0:ndeg,0:ndeg),s(0:ndeg,0:ndeg)
    2905           double precision  v0,omega,ae,gm !,clm,slm
    2906           double precision  plm(ndeg)
    2907           common /gmm1/ v0,omega,ae,gm,c,s
    2908           integer lmin,lmax
     3029         parameter (ndeg=90,ndeg2=2*ndeg,nd2p3=ndeg2+3)
     3030         double precision c(0:ndeg,0:ndeg),s(0:ndeg,0:ndeg)
     3031         double precision  v0,omega,ae,gm !,clm,slm
     3032         double precision  plm(ndeg)
     3033         common /gmm1/ v0,omega,ae,gm,c,s
     3034         integer lmin,lmax
    29093035          double precision root(nd2p3)
    29103036          double precision r
    … …  
    29163042        double precision x,xi,sum
    29173043       
    2918           ae= 3396000.d0
    2919           gm =42828.37d9
    2920           omega=0.70882181d-4
     3044         ae= 3396000.d0
     3045         gm =42828.37d9
     3046         omega=0.70882181d-4
    29213047! this value makes the equatorial mean radius equal to 3396 km.
    2922         do k=1,nd2p3
    2923         root(k)=sqrt(dble(k))
    2924         enddo
     3048       do k=1,nd2p3
     3049        root(k)=sqrt(dble(k))
     3050       enddo
    29253051!  initialize coefficients c and s
    29263052c(:,:)=0
    … …  
    94889614!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    94899615
    9490         r = 3396000.d0 !  MOLA potential surface
    9491         xi = ae/r
    9492         m=0
    9493         x=0.
    9494         sum=0.
    9495         call lgndr(lmax,m,x,plm,root)
    9496         do l=2, lmax
    9497           sum = sum + c(l,m) *xi**l *plm(l-m+1)
    9498         enddo
    9499         sum = sum+1.
    9500 ! centrifugal potential 
    9501         v0=(gm*sum/r + 0.5*omega**2 * r**2)
    9502 !       write(out,*)'v0=', v0
    9503         return
     9616       r = 3396000.d0 !  MOLA potential surface
     9617       xi = ae/r
     9618       m=0
     9619       x=0.
     9620       sum=0.
     9621       call lgndr(lmax,m,x,plm,root)
     9622       do l=2, lmax
     9623         sum = sum + c(l,m) *xi**l *plm(l-m+1)
     9624       enddo
     9625       sum = sum+1.
     9626! centrifugal potential       
     9627       v0=(gm*sum/r + 0.5*omega**2 * r**2)
     9628!       write(out,*)'v0=', v0
     9629       return
    95049630
    95059631end subroutine areoid_ini
    … …  
    95089634
    95099635
    9510         subroutine geoid(dlon,dlat,rg)
     9636       subroutine geoid(dlon,dlat,rg)
    95119637! marsgeoid.f
    95129638! calculate radii of surface of constant gravity potential on rotating body
    … …  
    95219647       
    95229648        double precision pi,d2r
    9523         parameter (pi=3.141592653589792D0, d2r=pi/180.d0)
     9649       parameter (pi=3.141592653589792D0, d2r=pi/180.d0)
    95249650! data structure of gravity field
    95259651        integer ndeg,ndeg1,ndeg2,nd2p3
    9526         parameter (ndeg=90,ndeg1=ndeg+1,ndeg2=2*ndeg,nd2p3=ndeg2+3)
     9652       parameter (ndeg=90,ndeg1=ndeg+1,ndeg2=2*ndeg,nd2p3=ndeg2+3)
    95279653        double precision v0,omega,ae,gm
    9528         double precision c(0:ndeg,0:ndeg),s(0:ndeg,0:ndeg)
    9529         common /gmm1/v0,omega,ae,gm,c,s
     9654       double precision c(0:ndeg,0:ndeg),s(0:ndeg,0:ndeg)
     9655       common /gmm1/v0,omega,ae,gm,c,s
    95309656        integer lmin,lmax
    95319657        double precision root(nd2p3)
    95329658        double precision r
    9533         common /sqr/ lmin,lmax,root,r
    9534         double precision plm(ndeg1)
    9535         double precision tol
     9659       common /sqr/ lmin,lmax,root,r
     9660       double precision plm(ndeg1)
     9661       double precision tol
    95369662        data tol /0.125d0/  ! tolerence on computed value of rg
    95379663
    … …  
    95419667        double precision diff
    95429668! save r
    9543         rg = r
    9544 
    9545         rlon=dlon*d2r
    9546         rlat=dlat*d2r
    9547         x = sin(rlat)
    9548         cslt= cos(rlat)
    9549 
    9550         do i=1,8 ! usually 3 iterations suffice
    9551 
    9552           xi = ae/r
    9553           sum = 0.
    9554           do m=0, lmax
    9555            call lgndr(lmax,m,x,plm,root)
    9556            do l=m, lmax
    9557              cslm =  c(l,m)*cos(m*rlon) + s(l,m)*sin(m*rlon)
    9558              sum = sum + cslm*xi**l *plm(l-m+1)
    9559            enddo
    9560           enddo
    9561           sum=sum+1.
    9562 ! centrifugal potential 
    9563           rg=(gm*sum+0.5*(omega**2)*(r**3)*(cslt**2))/v0
     9669       rg = r
     9670
     9671       rlon=dlon*d2r
     9672       rlat=dlat*d2r
     9673       x = sin(rlat)
     9674       cslt= cos(rlat)
     9675
     9676        do i=1,8 ! usually 3 iterations suffice
     9677
     9678         xi = ae/r
     9679         sum = 0.
     9680         do m=0, lmax
     9681          call lgndr(lmax,m,x,plm,root)
     9682          do l=m, lmax
     9683            cslm =  c(l,m)*cos(m*rlon) + s(l,m)*sin(m*rlon)
     9684            sum = sum + cslm*xi**l *plm(l-m+1)
     9685          enddo
     9686         enddo
     9687         sum=sum+1.
     9688! centrifugal potential       
     9689         rg=(gm*sum+0.5*(omega**2)*(r**3)*(cslt**2))/v0
    95649690!          write(out,*) i,r,rg
    95659691          diff = r-rg
    95669692          r = rg
    95679693          if( abs(diff).lt. tol) goto 400
    9568         enddo ! do i=1,8
    9569 
    9570  400    continue
    9571 
    9572         end subroutine geoid
     9694       enddo ! do i=1,8
     9695
     9696 400       continue
     9697
     9698       end subroutine geoid
    95739699
    95749700!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    96089734      if (LMAX.gt.M+1) then
    96099735         DO LL=M+2,LMAX
    9610            PLL=SQR(2*LL+1)*SQR(2*LL-1)/(SQR(LL+M)*SQR(LL-M)) &
     9736          PLL=SQR(2*LL+1)*SQR(2*LL-1)/(SQR(LL+M)*SQR(LL-M)) &
    96119737     & * (X*PMMP1-PMM*SQR(LL+M-1)*SQR(LL-M-1)/(SQR(2*LL-1)*SQR(2*LL-3)))
    96129738           PLM(LL-M+1)=OM*PLL
Note: See TracChangeset for help on using the changeset viewer.