Changeset 2402

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:

• trunk/UTIL/zrecast.F90 (modified) (25 diffs)

trunk/UTIL/zrecast.F90

@@ -89 +89 @@
 integer nbvarinfile ! # of variables in input file
 integer nbattr ! # of attributes of a given variable in input file
+integer nbvar3dinfile ! # of 3D (lon,lat,time) variables in input file
 integer nbvar4dinfile ! # of 4D (lon,lat,alt,time) variables in input file
 integer outfid ! NetCDF output file ID
@@ -108 +109 @@
 character (len=64), dimension(:), allocatable :: var
 ! var(): names of variables that will be processed
+integer, allocatable :: varndims(:) ! number of dimensions of var(:) variables
 integer nbvar ! # of variables to process
 integer,dimension(:),allocatable :: var_id ! IDs of variables var() (in outfile)
@@ -133 +135 @@
 real,dimension(:,:,:,:),allocatable :: ra_gcm ! radial distance (m) to
                                           ! center of Mars, (of GCM grid points)
+real,dimension(:,:,:),allocatable :: surfdata ! to store a GCM 3D dataset
 real,dimension(:,:,:,:),allocatable :: indata ! to store a GCM 4D dataset
 real,dimension(:,:,:,:),allocatable :: outdata ! to store a 4D dataset
@@ -202 +205 @@
 write(*,*)" The following variables have been found:"
 nbvar4dinfile=0
+nbvar3dinfile=0
 do i=1,nbvarinfile
   ! get name of variable # i
@@ -211 +215 @@
     write(*,*) trim(tmpvarname)
   endif
+  if (tmpndims.eq.3) then
+    nbvar3dinfile=nbvar3dinfile+1
+    write(*,*) trim(tmpvarname)
+  endif
 enddo
 
-allocate(var(nbvar4dinfile),stat=ierr)
+allocate(var(nbvar4dinfile+nbvar3dinfile),stat=ierr)
 if (ierr.ne.0) then
-  write(*,*) "Error! failed memory allocation of var(nbvar4dinfile)"
+  write(*,*) "Error! failed memory allocation of var(nbvar4dinfile+nbvar3dinfile)"
+  stop
+endif
+allocate(varndims(nbvar4dinfile+nbvar3dinfile),stat=ierr)
+if (ierr.ne.0) then
+  write(*,*) "Error! failed memory allocation of varndims(nbvar4dinfile+nbvar3dinfile)"
   stop
 endif
@@ -229 +242 @@
   ierr=NF_INQ_VARID(infid,tmpvarname,tmpvarid)
   if (ierr.eq.NF_NOERR) then ! valid name
-    ! also check that it is indeed a 4D variable
+    ! also check that it is indeed a 3D or 4D variable
     ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims)
-    if (tmpndims.eq.4) then
+    if ((tmpndims.eq.4).or.(tmpndims.eq.3)) then
       nbvar=nbvar+1
-      var(nbvar)=tmpvarname    
+      var(nbvar)=tmpvarname
+      varndims(nbvar)=tmpndims
     else
-      write(*,*) 'Error: ',trim(tmpvarname),' is not a 4D variable'
+      write(*,*) 'Error: ',trim(tmpvarname),' is not a 3D or 4D variable'
       write(*,*) '       (we''ll skip that one)'
     endif
@@ -249 +263 @@
   nbvar=0
   do i=1,nbvarinfile
-    ! look for 4D variables
+    ! look for 3D and 4D variables
     ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
-    if (tmpndims.eq.4) then
+    if ((tmpndims.eq.4).or.(tmpndims.eq.3)) then
       nbvar=nbvar+1
       ! get the corresponding name
       ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
       var(nbvar)=tmpvarname
+      varndims(nbvar)=tmpndims
     endif
   enddo
@@ -698 +713 @@
   if (.not.auto_mcd_levels) then
     write(*,*) ""
-    write(*,*) "Enter min and max of vertical coordinate (Pa or m)"
-    write(*,*) " (in that order and on the same line)"
     if (ztype.eq.1) then ! pressure coordinate
-      read(*,*) pmin,pmax
+      write(*,*) "Enter min and max of vertical coordinate (Pa)"
+      write(*,*) " (in that order and on the same line)"
+      read(*,*,iostat=ierr) pmin,pmax
+      if ((ierr/=0).or.(pmin.gt.pmax)) then
+        write(*,*) "Error: wrong vertical coordinate boundaries"
+        write(*,*) "       pmin=",pmin
+        write(*,*) "       pmax=",pmax
+        stop
+      endif
     else ! altitude coordinate, except in MCD case
-      read(*,*) zamin,zamax
+      write(*,*) "Enter min and max of vertical coordinate (m)"
+      write(*,*) " (in that order and on the same line)"
+      read(*,*,iostat=ierr) zamin,zamax
+      if ((ierr/=0).or.(zamin.gt.zamax)) then
+        write(*,*) "Error: wrong vertical coordinate boundaries"
+        write(*,*) "       zamin=",zamin
+        write(*,*) "       zamax=",zamax
+        stop
+      endif
     endif
   endif
@@ -710 +739 @@
   if (ztype.eq.1) then ! pressure coordinate
     write(*,*) "Number of levels along vertical coordinate?"
-    read(*,*) nblev
+    read(*,*,iostat=ierr) nblev
+    if (ierr/=0) then
+      write(*,*) "Error: bad input value for nblev ", nblev
+      stop
+    endif
     allocate(plevel(nblev),stat=ierr)
     if (ierr.ne.0) then
@@ -728 +761 @@
   else if (ztype.eq.2) then ! above areoid heights
     write(*,*) "Number of levels along vertical coordinate?"
-    read(*,*) nblev
+    read(*,*,iostat=ierr) nblev
+    if (ierr/=0) then
+      write(*,*) "Error: bad input value for nblev ", nblev
+      stop
+    endif
     allocate(zareoid(nblev),stat=ierr)
     if (ierr.ne.0) then
@@ -756 +793 @@
     else
       write(*,*) "Number of levels along vertical coordinate?"
-      read(*,*) nblev
+      read(*,*,iostat=ierr) nblev
+      if (ierr/=0) then
+        write(*,*) "Error: bad input value for nblev ", nblev
+        stop
+      endif
       allocate(zsurface(nblev),stat=ierr)
       if (ierr.ne.0) then
@@ -773 +814 @@
   else ! distance to center of planet
    write(*,*) "Number of levels along vertical coordinate?"
-    read(*,*) nblev
+    read(*,*,iostat=ierr) nblev
+    if (ierr/=0) then
+      write(*,*) "Error: bad input value for nblev ", nblev
+      stop
+    endif
     allocate(zradius(nblev),stat=ierr)
     if (ierr.ne.0) then
@@ -792 +837 @@
   write(*,*) ""
   write(*,*) "Number of levels along vertical coordinate?"
-  read(*,*) nblev
+  read(*,*,iostat=ierr) nblev
+  if (ierr/=0) then
+    write(*,*) "Error: bad input value for nblev ", nblev
+    stop
+  endif
   if (ztype.eq.1) then ! pressure coordinate
     allocate(plevel(nblev),stat=ierr)
@@ -804 +853 @@
     write(*,*) " (one value per line)"
     do i=1,nblev
-      read(*,*) plevel(i)
+      read(*,*,iostat=ierr) plevel(i)
+      if (ierr/=0) then
+        write(*,*) "Error: bad input value for plevel(i):",plevel(i)
+        stop
+      endif
     enddo
   else if (ztype.eq.2) then ! above areoid heights
@@ -816 +869 @@
     write(*,*) " from min to max (one value per line)"
     do i=1,nblev
-      read(*,*) zareoid(i)
+      read(*,*,iostat=ierr) zareoid(i)
+      if (ierr/=0) then
+        write(*,*) "Error: bad input value for zareoid(i):",zareoid(i)
+        stop
+      endif
     enddo
   else if (ztype.eq.3) then ! above local surface
@@ -828 +885 @@
     write(*,*) " from min to max (one value per line)"
     do i=1,nblev
-      read(*,*) zsurface(i)
+      read(*,*,iostat=ierr) zsurface(i)
+      if (ierr/=0) then
+        write(*,*) "Error: bad input value for zsurface(i):",zsurface(i)
+        stop
+      endif
     enddo
   else ! distance to center of planet
@@ -840 +901 @@
     write(*,*) " from min to max (one value per line)"
     do i=1,nblev
-      read(*,*) zradius(i)
+      read(*,*,iostat=ierr) zradius(i)
+      if (ierr/=0) then
+        write(*,*) "Error: bad input value for zradius(i):",zradius(i)
+        stop
+      endif
     enddo
   endif
@@ -1474 +1539 @@
 endif
 
-! 3.3.3 Define 4D variables
+! 3.3.3 Define 3D and 4D variables
 
 ! add pressure or zareoid
@@ -1538 +1603 @@
   stop
 endif
-do i=1,nbvar
+do i=1,nbvar ! loop on 3D+4D variables
+  if (var(i)=="ps") cycle   ! skip surface pressure (already processed)
+  ! other variables
   write(*,*) ""
   write(*,*) "Creating ",trim(var(i))
-  ! define the variable
-  ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,4,datashape,var_id(i))
-  if (ierr.ne.NF_NOERR) then
-    write(*,*) 'Error, could not define variable ',trim(var(i))
-    write(*,*) NF_STRERROR(ierr)
-    stop 
+  if (varndims(i).eq.3) then
+    ! define the variable
+    ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,3,surfdatashape,var_id(i))
+    if (ierr.ne.NF_NOERR) then
+      write(*,*) 'Error, could not define variable ',trim(var(i))
+      write(*,*) NF_STRERROR(ierr)
+      stop 
+    endif
+  elseif (varndims(i).eq.4) then
+    ! define the variable
+    ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,4,datashape,var_id(i))
+    if (ierr.ne.NF_NOERR) then
+      write(*,*) 'Error, could not define variable ',trim(var(i))
+      write(*,*) NF_STRERROR(ierr)
+      stop 
+    endif
+  else
+    write(*,*) "Unexpected value for varndims(i)!"
+    write(*,*) " i=",i," varndims(i)=",varndims(i)
+    stop
   endif
 
@@ -1780 +1861 @@
 endif
 
+! allocate surfdata() to store input 3D GCM data
+allocate(surfdata(lonlength,latlength,timelength),stat=ierr)
+if (ierr.ne.0) then
+  write(*,*) "Error: Failed to allocate surfdata(lonlength,latlength,timelength)"
+  write(*,*) "       lonlength=",lonlength," latlength=",latlength
+  write(*,*) "       timelength=",timelength
+  stop
+endif
+
+do i=1,nbvar ! loop on 3D+4D variables
+  if (varndims(i).eq.4) cycle ! skip 4D variables
+  if (var(i)=="ps") cycle ! skip surface pressure (already processed)
+  ! load the variable from input file
+  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
+  if (ierr.ne.NF_NOERR) then
+    write(*,*) "Error: Failed to get ",trim(var(i)),"ID"
+    write(*,*) NF_STRERROR(ierr)
+    stop
+  endif
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,surfdata)
+  if (ierr.ne.NF_NOERR) then
+    write(*,*) "Error: Failed to load ",trim(var(i))
+    write(*,*) NF_STRERROR(ierr)
+    stop
+  endif
+  ! write the variable to output file
+  ierr=NF_PUT_VARA_REAL(outfid,var_id(i),corner(1:3),edges(1:3),surfdata)
+  if (ierr.ne.NF_NOERR) then
+    write(*,*) "Error: Could not write ",trim(var(i))," data to output file"
+    write(*,*) NF_STRERROR(ierr)
+    stop
+  endif
+enddo ! of do i=1,nbvar
+
 !===============================================================================
 ! 4. Interpolate and write 4D variables
@@ -1804 +1919 @@
 ! 4.1 If output is in pressure coordinate
 if (ztype.eq.1) then
-  do i=1,nbvar ! loop on 4D variable to process
+  do i=1,nbvar ! loop on 3D+4D variable to process
+    ! Skip 3D variables
+    if (varndims(i).eq.3) cycle
   ! identify and read a dataset
   ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
@@ -1850 +1967 @@
 ! 4.2 If output is in above areoid altitude
 if (ztype.eq.2) then
-  do i=1,nbvar ! loop on 4D variable to process
+  do i=1,nbvar ! loop on 3D+4D variable to process
+    ! Skip 3D variables
+    if (varndims(i).eq.3) cycle
   ! identify and read a dataset
   ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
@@ -1901 +2020 @@
 ! 4.3 If output is in above local surface altitude
 if (ztype.eq.3) then
-  do i=1,nbvar ! loop on 4D variable to process
+  do i=1,nbvar ! loop on 3D+4D variable to process
+    ! Skip 3D variables
+    if (varndims(i).eq.3) cycle
     write(*,*) " "
     write(*,*) "Processing "//trim(var(i))
@@ -1956 +2077 @@
 ! 4.4 If output is in radial distance to center of planet
 if (ztype.eq.4) then
-  do i=1,nbvar ! loop on 4D variable to process
+  do i=1,nbvar ! loop on 3D+4D variable to process
+    ! Skip 3D variables
+    if (varndims(i).eq.3) cycle
   ! identify and read a dataset
   ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
@@ -2011 +2134 @@
   write(*,*) NF_STRERROR(ierr)
 endif
+
+write(*,*) 
+write(*,*) "zrecast: All is well that ends well."
 
 end program