Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 2400)
+++ /trunk/LMDZ.MARS/README	(revision 2401)
@@ -3124,2 +3124,6 @@
 == 07/07/2020 == DB
 Add kstar parameter to control kmin value. kmi=1/lambda_max, to ensure the "subgrid scale" characteristic, we have to constrain the GW's wavelength by the size of the mesh
+
+== 07/07/2020 == EM
+Improve "zrecast" utility: add possibility to transfer 3D (lon-lat-time)
+variables to output file and add some extra sanity checks of user inputs.
Index: /trunk/LMDZ.MARS/util/zrecast.F90
===================================================================
--- /trunk/LMDZ.MARS/util/zrecast.F90	(revision 2400)
+++ /trunk/LMDZ.MARS/util/zrecast.F90	(revision 2401)
@@ -79,4 +79,5 @@
 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
@@ -97,4 +98,5 @@
 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)
@@ -120,4 +122,5 @@
 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
@@ -186,4 +189,5 @@
 write(*,*)" The following variables have been found:"
 nbvar4dinfile=0
+nbvar3dinfile=0
 do i=1,nbvarinfile
   ! get name of variable # i
@@ -195,9 +199,18 @@
     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
@@ -213,11 +226,12 @@
   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
@@ -233,11 +247,12 @@
   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
@@ -622,10 +637,24 @@
   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
@@ -634,5 +663,9 @@
   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
@@ -652,5 +685,9 @@
   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
@@ -680,5 +717,9 @@
     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
@@ -697,5 +738,9 @@
   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
@@ -716,5 +761,9 @@
   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)
@@ -728,5 +777,9 @@
     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
@@ -740,5 +793,9 @@
     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
@@ -752,5 +809,9 @@
     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
@@ -764,5 +825,9 @@
     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
@@ -1317,5 +1382,5 @@
 endif
 
-! 3.3.3 Define 4D variables
+! 3.3.3 Define 3D and 4D variables
 
 ! add pressure or zareoid
@@ -1381,13 +1446,29 @@
   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
 
@@ -1605,4 +1686,38 @@
 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
@@ -1629,5 +1744,7 @@
 ! 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)
@@ -1675,5 +1792,7 @@
 ! 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)
@@ -1726,5 +1845,7 @@
 ! 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))
@@ -1781,5 +1902,7 @@
 ! 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)
@@ -1836,4 +1959,7 @@
   write(*,*) NF_STRERROR(ierr)
 endif
+
+write(*,*) 
+write(*,*) "zrecast: All is well that ends well."
 
 end program
@@ -2846,10 +2972,10 @@
 ! data structure of gravity field
         double precision v0,omega,ae,gm
-	double precision clm(0:ndeg,0:ndeg),slm(0:ndeg,0:ndeg)
-	common /gmm1/v0,omega,ae,gm,clm,slm
+       double precision clm(0:ndeg,0:ndeg),slm(0:ndeg,0:ndeg)
+       common /gmm1/v0,omega,ae,gm,clm,slm
         integer lmin,lmax
         double precision root(nd2p3)
         double precision requator
-	common /sqr/ lmin,lmax,root,requator
+       common /sqr/ lmin,lmax,root,requator
 
 !===============================================================================
@@ -2865,5 +2991,5 @@
    dlat=lat(jloop) ! double precision version of lat
    call geoid(dlon,dlat,dareoid)
-   rareoid(iloop,jloop)=dareoid
+   rareoid(iloop,jloop)=real(dareoid)
  enddo
 enddo
@@ -2901,10 +3027,10 @@
 ! data structure of gravity field
          integer ndeg,ndeg2,nd2p3
-	  parameter (ndeg=90,ndeg2=2*ndeg,nd2p3=ndeg2+3)
-	  double precision c(0:ndeg,0:ndeg),s(0:ndeg,0:ndeg)
-	  double precision  v0,omega,ae,gm !,clm,slm
-	  double precision  plm(ndeg)
-	  common /gmm1/ v0,omega,ae,gm,c,s
-	  integer lmin,lmax
+         parameter (ndeg=90,ndeg2=2*ndeg,nd2p3=ndeg2+3)
+         double precision c(0:ndeg,0:ndeg),s(0:ndeg,0:ndeg)
+         double precision  v0,omega,ae,gm !,clm,slm
+         double precision  plm(ndeg)
+         common /gmm1/ v0,omega,ae,gm,c,s
+         integer lmin,lmax
           double precision root(nd2p3)
           double precision r
@@ -2916,11 +3042,11 @@
         double precision x,xi,sum
         
-	  ae= 3396000.d0
-	  gm =42828.37d9
-	  omega=0.70882181d-4
+         ae= 3396000.d0
+         gm =42828.37d9
+         omega=0.70882181d-4
 ! this value makes the equatorial mean radius equal to 3396 km.
-	do k=1,nd2p3
-	 root(k)=sqrt(dble(k))
-	enddo
+       do k=1,nd2p3
+        root(k)=sqrt(dble(k))
+       enddo
 !  initialize coefficients c and s
 c(:,:)=0
@@ -9488,18 +9614,18 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-	r = 3396000.d0 !  MOLA potential surface
-	xi = ae/r
-	m=0
-	x=0.
-	sum=0.
-	call lgndr(lmax,m,x,plm,root)
-	do l=2, lmax
-	  sum = sum + c(l,m) *xi**l *plm(l-m+1)
-	enddo
-	sum = sum+1.
-! centrifugal potential	
-	v0=(gm*sum/r + 0.5*omega**2 * r**2)
-!	write(out,*)'v0=', v0
-	return
+       r = 3396000.d0 !  MOLA potential surface
+       xi = ae/r
+       m=0
+       x=0.
+       sum=0.
+       call lgndr(lmax,m,x,plm,root)
+       do l=2, lmax
+         sum = sum + c(l,m) *xi**l *plm(l-m+1)
+       enddo
+       sum = sum+1.
+! centrifugal potential       
+       v0=(gm*sum/r + 0.5*omega**2 * r**2)
+!       write(out,*)'v0=', v0
+       return
 
 end subroutine areoid_ini
@@ -9508,5 +9634,5 @@
 
 
-	subroutine geoid(dlon,dlat,rg)
+       subroutine geoid(dlon,dlat,rg)
 ! marsgeoid.f
 ! calculate radii of surface of constant gravity potential on rotating body
@@ -9521,17 +9647,17 @@
         
         double precision pi,d2r
-	parameter (pi=3.141592653589792D0, d2r=pi/180.d0)
+       parameter (pi=3.141592653589792D0, d2r=pi/180.d0)
 ! data structure of gravity field
         integer ndeg,ndeg1,ndeg2,nd2p3
-	parameter (ndeg=90,ndeg1=ndeg+1,ndeg2=2*ndeg,nd2p3=ndeg2+3)
+       parameter (ndeg=90,ndeg1=ndeg+1,ndeg2=2*ndeg,nd2p3=ndeg2+3)
         double precision v0,omega,ae,gm
-	double precision c(0:ndeg,0:ndeg),s(0:ndeg,0:ndeg)
-	common /gmm1/v0,omega,ae,gm,c,s
+       double precision c(0:ndeg,0:ndeg),s(0:ndeg,0:ndeg)
+       common /gmm1/v0,omega,ae,gm,c,s
         integer lmin,lmax
         double precision root(nd2p3)
         double precision r
-	common /sqr/ lmin,lmax,root,r
-	double precision plm(ndeg1)
-	double precision tol
+       common /sqr/ lmin,lmax,root,r
+       double precision plm(ndeg1)
+       double precision tol
         data tol /0.125d0/  ! tolerence on computed value of rg
 
@@ -9541,34 +9667,34 @@
         double precision diff
 ! save r
-	rg = r
-
-	rlon=dlon*d2r
-	rlat=dlat*d2r
-	x = sin(rlat)
-	cslt= cos(rlat)
-
- 	do i=1,8 ! usually 3 iterations suffice
-
-	  xi = ae/r
-	  sum = 0.
-	  do m=0, lmax
-	   call lgndr(lmax,m,x,plm,root)
-	   do l=m, lmax
-	     cslm =  c(l,m)*cos(m*rlon) + s(l,m)*sin(m*rlon)
-	     sum = sum + cslm*xi**l *plm(l-m+1)
-	   enddo
-	  enddo
-	  sum=sum+1.
-! centrifugal potential	
-	  rg=(gm*sum+0.5*(omega**2)*(r**3)*(cslt**2))/v0
+       rg = r
+
+       rlon=dlon*d2r
+       rlat=dlat*d2r
+       x = sin(rlat)
+       cslt= cos(rlat)
+
+        do i=1,8 ! usually 3 iterations suffice
+
+         xi = ae/r
+         sum = 0.
+         do m=0, lmax
+          call lgndr(lmax,m,x,plm,root)
+          do l=m, lmax
+            cslm =  c(l,m)*cos(m*rlon) + s(l,m)*sin(m*rlon)
+            sum = sum + cslm*xi**l *plm(l-m+1)
+          enddo
+         enddo
+         sum=sum+1.
+! centrifugal potential       
+         rg=(gm*sum+0.5*(omega**2)*(r**3)*(cslt**2))/v0
 !          write(out,*) i,r,rg
           diff = r-rg
           r = rg
           if( abs(diff).lt. tol) goto 400
-	enddo ! do i=1,8
-
- 400	continue
-
-	end subroutine geoid
+       enddo ! do i=1,8
+
+ 400       continue
+
+       end subroutine geoid
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -9608,5 +9734,5 @@
       if (LMAX.gt.M+1) then
          DO LL=M+2,LMAX
-	   PLL=SQR(2*LL+1)*SQR(2*LL-1)/(SQR(LL+M)*SQR(LL-M)) &
+          PLL=SQR(2*LL+1)*SQR(2*LL-1)/(SQR(LL+M)*SQR(LL-M)) &
      & * (X*PMMP1-PMM*SQR(LL+M-1)*SQR(LL-M-1)/(SQR(2*LL-1)*SQR(2*LL-3)))
            PLM(LL-M+1)=OM*PLL
