Changeset 2402 for trunk


Ignore:
Timestamp:
Jul 8, 2020, 11:37:08 AM (4 years ago)
Author:
emillour
Message:

General Utilities:
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/UTIL/zrecast.F90

    r1442 r2402  
    8989integer nbvarinfile ! # of variables in input file
    9090integer nbattr ! # of attributes of a given variable in input file
     91integer nbvar3dinfile ! # of 3D (lon,lat,time) variables in input file
    9192integer nbvar4dinfile ! # of 4D (lon,lat,alt,time) variables in input file
    9293integer outfid ! NetCDF output file ID
    … …  
    108109character (len=64), dimension(:), allocatable :: var
    109110! var(): names of variables that will be processed
     111integer, allocatable :: varndims(:) ! number of dimensions of var(:) variables
    110112integer nbvar ! # of variables to process
    111113integer,dimension(:),allocatable :: var_id ! IDs of variables var() (in outfile)
    … …  
    133135real,dimension(:,:,:,:),allocatable :: ra_gcm ! radial distance (m) to
    134136                                          ! center of Mars, (of GCM grid points)
     137real,dimension(:,:,:),allocatable :: surfdata ! to store a GCM 3D dataset
    135138real,dimension(:,:,:,:),allocatable :: indata ! to store a GCM 4D dataset
    136139real,dimension(:,:,:,:),allocatable :: outdata ! to store a 4D dataset
    … …  
    202205write(*,*)" The following variables have been found:"
    203206nbvar4dinfile=0
     207nbvar3dinfile=0
    204208do i=1,nbvarinfile
    205209  ! get name of variable # i
    … …  
    211215    write(*,*) trim(tmpvarname)
    212216  endif
     217  if (tmpndims.eq.3) then
     218    nbvar3dinfile=nbvar3dinfile+1
     219    write(*,*) trim(tmpvarname)
     220  endif
    213221enddo
    214222
    215 allocate(var(nbvar4dinfile),stat=ierr)
     223allocate(var(nbvar4dinfile+nbvar3dinfile),stat=ierr)
    216224if (ierr.ne.0) then
    217   write(*,*) "Error! failed memory allocation of var(nbvar4dinfile)"
     225  write(*,*) "Error! failed memory allocation of var(nbvar4dinfile+nbvar3dinfile)"
     226  stop
     227endif
     228allocate(varndims(nbvar4dinfile+nbvar3dinfile),stat=ierr)
     229if (ierr.ne.0) then
     230  write(*,*) "Error! failed memory allocation of varndims(nbvar4dinfile+nbvar3dinfile)"
    218231  stop
    219232endif
    … …  
    229242  ierr=NF_INQ_VARID(infid,tmpvarname,tmpvarid)
    230243  if (ierr.eq.NF_NOERR) then ! valid name
    231     ! also check that it is indeed a 4D variable
     244    ! also check that it is indeed a 3D or 4D variable
    232245    ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims)
    233     if (tmpndims.eq.4) then
     246    if ((tmpndims.eq.4).or.(tmpndims.eq.3)) then
    234247      nbvar=nbvar+1
    235       var(nbvar)=tmpvarname   
     248      var(nbvar)=tmpvarname
     249      varndims(nbvar)=tmpndims
    236250    else
    237       write(*,*) 'Error: ',trim(tmpvarname),' is not a 4D variable'
     251      write(*,*) 'Error: ',trim(tmpvarname),' is not a 3D or 4D variable'
    238252      write(*,*) '       (we''ll skip that one)'
    239253    endif
    … …  
    249263  nbvar=0
    250264  do i=1,nbvarinfile
    251     ! look for 4D variables
     265    ! look for 3D and 4D variables
    252266    ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
    253     if (tmpndims.eq.4) then
     267    if ((tmpndims.eq.4).or.(tmpndims.eq.3)) then
    254268      nbvar=nbvar+1
    255269      ! get the corresponding name
    256270      ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
    257271      var(nbvar)=tmpvarname
     272      varndims(nbvar)=tmpndims
    258273    endif
    259274  enddo
    … …  
    698713  if (.not.auto_mcd_levels) then
    699714    write(*,*) ""
    700     write(*,*) "Enter min and max of vertical coordinate (Pa or m)"
    701     write(*,*) " (in that order and on the same line)"
    702715    if (ztype.eq.1) then ! pressure coordinate
    703       read(*,*) pmin,pmax
     716      write(*,*) "Enter min and max of vertical coordinate (Pa)"
     717      write(*,*) " (in that order and on the same line)"
     718      read(*,*,iostat=ierr) pmin,pmax
     719      if ((ierr/=0).or.(pmin.gt.pmax)) then
     720        write(*,*) "Error: wrong vertical coordinate boundaries"
     721        write(*,*) "       pmin=",pmin
     722        write(*,*) "       pmax=",pmax
     723        stop
     724      endif
    704725    else ! altitude coordinate, except in MCD case
    705       read(*,*) zamin,zamax
     726      write(*,*) "Enter min and max of vertical coordinate (m)"
     727      write(*,*) " (in that order and on the same line)"
     728      read(*,*,iostat=ierr) zamin,zamax
     729      if ((ierr/=0).or.(zamin.gt.zamax)) then
     730        write(*,*) "Error: wrong vertical coordinate boundaries"
     731        write(*,*) "       zamin=",zamin
     732        write(*,*) "       zamax=",zamax
     733        stop
     734      endif
    706735    endif
    707736  endif
    … …  
    710739  if (ztype.eq.1) then ! pressure coordinate
    711740    write(*,*) "Number of levels along vertical coordinate?"
    712     read(*,*) nblev
     741    read(*,*,iostat=ierr) nblev
     742    if (ierr/=0) then
     743      write(*,*) "Error: bad input value for nblev ", nblev
     744      stop
     745    endif
    713746    allocate(plevel(nblev),stat=ierr)
    714747    if (ierr.ne.0) then
    … …  
    728761  else if (ztype.eq.2) then ! above areoid heights
    729762    write(*,*) "Number of levels along vertical coordinate?"
    730     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
    731768    allocate(zareoid(nblev),stat=ierr)
    732769    if (ierr.ne.0) then
    … …  
    756793    else
    757794      write(*,*) "Number of levels along vertical coordinate?"
    758       read(*,*) nblev
     795      read(*,*,iostat=ierr) nblev
     796      if (ierr/=0) then
     797        write(*,*) "Error: bad input value for nblev ", nblev
     798        stop
     799      endif
    759800      allocate(zsurface(nblev),stat=ierr)
    760801      if (ierr.ne.0) then
    … …  
    773814  else ! distance to center of planet
    774815   write(*,*) "Number of levels along vertical coordinate?"
    775     read(*,*) nblev
     816    read(*,*,iostat=ierr) nblev
     817    if (ierr/=0) then
     818      write(*,*) "Error: bad input value for nblev ", nblev
     819      stop
     820    endif
    776821    allocate(zradius(nblev),stat=ierr)
    777822    if (ierr.ne.0) then
    … …  
    792837  write(*,*) ""
    793838  write(*,*) "Number of levels along vertical coordinate?"
    794   read(*,*) nblev
     839  read(*,*,iostat=ierr) nblev
     840  if (ierr/=0) then
     841    write(*,*) "Error: bad input value for nblev ", nblev
     842    stop
     843  endif
    795844  if (ztype.eq.1) then ! pressure coordinate
    796845    allocate(plevel(nblev),stat=ierr)
    … …  
    804853    write(*,*) " (one value per line)"
    805854    do i=1,nblev
    806       read(*,*) plevel(i)
     855      read(*,*,iostat=ierr) plevel(i)
     856      if (ierr/=0) then
     857        write(*,*) "Error: bad input value for plevel(i):",plevel(i)
     858        stop
     859      endif
    807860    enddo
    808861  else if (ztype.eq.2) then ! above areoid heights
    … …  
    816869    write(*,*) " from min to max (one value per line)"
    817870    do i=1,nblev
    818       read(*,*) zareoid(i)
     871      read(*,*,iostat=ierr) zareoid(i)
     872      if (ierr/=0) then
     873        write(*,*) "Error: bad input value for zareoid(i):",zareoid(i)
     874        stop
     875      endif
    819876    enddo
    820877  else if (ztype.eq.3) then ! above local surface
    … …  
    828885    write(*,*) " from min to max (one value per line)"
    829886    do i=1,nblev
    830       read(*,*) zsurface(i)
     887      read(*,*,iostat=ierr) zsurface(i)
     888      if (ierr/=0) then
     889        write(*,*) "Error: bad input value for zsurface(i):",zsurface(i)
     890        stop
     891      endif
    831892    enddo
    832893  else ! distance to center of planet
    … …  
    840901    write(*,*) " from min to max (one value per line)"
    841902    do i=1,nblev
    842       read(*,*) zradius(i)
     903      read(*,*,iostat=ierr) zradius(i)
     904      if (ierr/=0) then
     905        write(*,*) "Error: bad input value for zradius(i):",zradius(i)
     906        stop
     907      endif
    843908    enddo
    844909  endif
    … …  
    14741539endif
    14751540
    1476 ! 3.3.3 Define 4D variables
     1541! 3.3.3 Define 3D and 4D variables
    14771542
    14781543! add pressure or zareoid
    … …  
    15381603  stop
    15391604endif
    1540 do i=1,nbvar
     1605do i=1,nbvar ! loop on 3D+4D variables
     1606  if (var(i)=="ps") cycle   ! skip surface pressure (already processed)
     1607  ! other variables
    15411608  write(*,*) ""
    15421609  write(*,*) "Creating ",trim(var(i))
    1543   ! define the variable
    1544   ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,4,datashape,var_id(i))
    1545   if (ierr.ne.NF_NOERR) then
    1546     write(*,*) 'Error, could not define variable ',trim(var(i))
    1547     write(*,*) NF_STRERROR(ierr)
    1548     stop
     1610  if (varndims(i).eq.3) then
     1611    ! define the variable
     1612    ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,3,surfdatashape,var_id(i))
     1613    if (ierr.ne.NF_NOERR) then
     1614      write(*,*) 'Error, could not define variable ',trim(var(i))
     1615      write(*,*) NF_STRERROR(ierr)
     1616      stop
     1617    endif
     1618  elseif (varndims(i).eq.4) then
     1619    ! define the variable
     1620    ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,4,datashape,var_id(i))
     1621    if (ierr.ne.NF_NOERR) then
     1622      write(*,*) 'Error, could not define variable ',trim(var(i))
     1623      write(*,*) NF_STRERROR(ierr)
     1624      stop
     1625    endif
     1626  else
     1627    write(*,*) "Unexpected value for varndims(i)!"
     1628    write(*,*) " i=",i," varndims(i)=",varndims(i)
     1629    stop
    15491630  endif
    15501631
    … …  
    17801861endif
    17811862
     1863! allocate surfdata() to store input 3D GCM data
     1864allocate(surfdata(lonlength,latlength,timelength),stat=ierr)
     1865if (ierr.ne.0) then
     1866  write(*,*) "Error: Failed to allocate surfdata(lonlength,latlength,timelength)"
     1867  write(*,*) "       lonlength=",lonlength," latlength=",latlength
     1868  write(*,*) "       timelength=",timelength
     1869  stop
     1870endif
     1871
     1872do i=1,nbvar ! loop on 3D+4D variables
     1873  if (varndims(i).eq.4) cycle ! skip 4D variables
     1874  if (var(i)=="ps") cycle ! skip surface pressure (already processed)
     1875  ! load the variable from input file
     1876  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
     1877  if (ierr.ne.NF_NOERR) then
     1878    write(*,*) "Error: Failed to get ",trim(var(i)),"ID"
     1879    write(*,*) NF_STRERROR(ierr)
     1880    stop
     1881  endif
     1882  ierr=NF_GET_VAR_REAL(infid,tmpvarid,surfdata)
     1883  if (ierr.ne.NF_NOERR) then
     1884    write(*,*) "Error: Failed to load ",trim(var(i))
     1885    write(*,*) NF_STRERROR(ierr)
     1886    stop
     1887  endif
     1888  ! write the variable to output file
     1889  ierr=NF_PUT_VARA_REAL(outfid,var_id(i),corner(1:3),edges(1:3),surfdata)
     1890  if (ierr.ne.NF_NOERR) then
     1891    write(*,*) "Error: Could not write ",trim(var(i))," data to output file"
     1892    write(*,*) NF_STRERROR(ierr)
     1893    stop
     1894  endif
     1895enddo ! of do i=1,nbvar
     1896
    17821897!===============================================================================
    17831898! 4. Interpolate and write 4D variables
    … …  
    18041919! 4.1 If output is in pressure coordinate
    18051920if (ztype.eq.1) then
    1806   do i=1,nbvar ! loop on 4D variable to process
     1921  do i=1,nbvar ! loop on 3D+4D variable to process
     1922    ! Skip 3D variables
     1923    if (varndims(i).eq.3) cycle
    18071924  ! identify and read a dataset
    18081925  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
    … …  
    18501967! 4.2 If output is in above areoid altitude
    18511968if (ztype.eq.2) then
    1852   do i=1,nbvar ! loop on 4D variable to process
     1969  do i=1,nbvar ! loop on 3D+4D variable to process
     1970    ! Skip 3D variables
     1971    if (varndims(i).eq.3) cycle
    18531972  ! identify and read a dataset
    18541973  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
    … …  
    19012020! 4.3 If output is in above local surface altitude
    19022021if (ztype.eq.3) then
    1903   do i=1,nbvar ! loop on 4D variable to process
     2022  do i=1,nbvar ! loop on 3D+4D variable to process
     2023    ! Skip 3D variables
     2024    if (varndims(i).eq.3) cycle
    19042025    write(*,*) " "
    19052026    write(*,*) "Processing "//trim(var(i))
    … …  
    19562077! 4.4 If output is in radial distance to center of planet
    19572078if (ztype.eq.4) then
    1958   do i=1,nbvar ! loop on 4D variable to process
     2079  do i=1,nbvar ! loop on 3D+4D variable to process
     2080    ! Skip 3D variables
     2081    if (varndims(i).eq.3) cycle
    19592082  ! identify and read a dataset
    19602083  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
    … …  
    20112134  write(*,*) NF_STRERROR(ierr)
    20122135endif
     2136
     2137write(*,*)
     2138write(*,*) "zrecast: All is well that ends well."
    20132139
    20142140end program
Note: See TracChangeset for help on using the changeset viewer.