Index: trunk/LMDZ.MARS/util/zrecast.F90
===================================================================
--- trunk/LMDZ.MARS/util/zrecast.F90	(revision 2432)
+++ trunk/LMDZ.MARS/util/zrecast.F90	(revision 2528)
@@ -62,4 +62,8 @@
 ! TN 10/2013 : Read and write controle field, if available
 ! FF 10/2020 : Force interpolation in log space for variable starting by "rho"
+! FF 05/2021 : Force interpolation in log space for variable starting by "num_"
+! FF 05/2021 : "phisinit" added to the putput file
+! FF 05/2021 : Conversion to pressure coordinate possible even if temperature
+!              not available
 !
 implicit none
@@ -84,5 +88,5 @@
 integer outfid ! NetCDF output file ID
 integer lon_dimid,lat_dimid,alt_dimid,time_dimid,ctl_dimid ! NetCDF dimension IDs
-integer lon_varid,lat_varid,alt_varid,time_varid,ctl_varid
+integer lon_varid,lat_varid,alt_varid,time_varid,ctl_varid,phis_varid
 integer gcm_layers_dimid ! NetCDF dimension ID for # of layers in GCM
 integer sigma_varid,aps_varid,bps_varid
@@ -94,4 +98,5 @@
 integer,dimension(4) :: datashape ! shape of 4D datasets
 integer,dimension(3) :: surfdatashape ! shape of 3D (surface+time) datasets
+integer,dimension(2) :: phisdatashape ! shape of 2D (surface) datasets
 
 real :: miss_val= 1.E+20 ! special "missing value" to specify missing data
@@ -140,4 +145,5 @@
 real,dimension(:),allocatable :: zradius ! distance to center of planet
 logical :: have_rho ! Flag: true if density 'rho' is available
+logical :: have_temp ! Flag: true if temperature is available in some way
 logical :: have_sigma ! Flag: true if sigma levels are known (false if hybrid
                       !       coordinates are used)
@@ -937,9 +943,12 @@
      ierr=NF_INQ_VARID(infid,"teta",tmpvarid)
      if (ierr.ne.NF_NOERR) then
-       stop "Error: Failed to get t or teta ID"
+       write(*,*) "Error: Failed to get temp, t or teta ID"
+       have_temp=.false.
+       if (ztype.ne.1) stop
      endif
        ierr=NF_GET_VAR_REAL(infid,tmpvarid,teta)
        if (ierr.ne.NF_NOERR) then
-         stop "Error: Failed reading atmospheric temperature"
+         write(*,*) "Failed reading atmospheric temperature"
+         if (ztype.ne.1) stop
        endif
      
@@ -953,4 +962,5 @@
     enddo
   enddo
+  have_temp=.true.
 
   endif
@@ -960,4 +970,5 @@
     stop "Error: Failed reading atmospheric temperature"
   endif
+  have_temp=.true.
 endif
 
@@ -1358,4 +1369,21 @@
 endif
 
+! Add phisinit in outputfile
+phisdatashape(1)=lon_dimid
+phisdatashape(2)=lat_dimid
+ierr=NF_DEF_VAR(outfid,"phisinit",NF_REAL,2,phisdatashape,phis_varid)
+text='Geopotential at the surface'
+ierr=NF_PUT_ATT_TEXT(outfid,phis_varid,'long_name',len_trim(text),text)
+if (ierr.ne.NF_NOERR) then
+  stop "Error: Problem writing long_name for phisinit"
+endif
+text='m2.s2'
+ierr=NF_PUT_ATT_TEXT(outfid,phis_varid,'units',len_trim(text),text)
+if (ierr.ne.NF_NOERR) then
+  stop "Error: Problem writing units for geopotential phisinit"
+endif
+
+
+
 ! 3.3.2 Define 3D variables (ie: surface+time variables)
 
@@ -1386,5 +1414,5 @@
 
 ! add pressure or zareoid
-if (ztype.eq.1) then ! pressure vertical coordinate
+if ((ztype.eq.1).and.(have_temp)) then ! pressure vertical coordinate
   ! zareoid dataset
   ierr=NF_DEF_VAR(outfid,"zareoid",NF_REAL,4,datashape,za_varid)
@@ -1671,4 +1699,12 @@
 endif
 
+! write phisinit
+  ierr=NF_PUT_VAR_REAL(outfid,phis_varid,phisinit)
+  if (ierr.ne.NF_NOERR) then
+    write(*,*) "Error: Could not write phisinit to output file"
+    write(*,*) NF_STRERROR(ierr)
+    stop
+  endif
+
 !===============================================================================
 ! 3.5 Write 3D variables
@@ -1810,5 +1846,5 @@
   ! interpolate dataset onto new grid
   ! check if variable is "rho" (to set flag for interpolation below)
-  if (var(i)(1:3).eq.'rho') then
+  if ((var(i)(1:3).eq.'rho').or.(var(i)(1:4).eq.'num_')) then
     j=1
     write(*,*) trim(var(i)), ' is interpolated in log scale' 
@@ -1866,5 +1902,5 @@
     ! interpolate dataset onto new grid
     ! check if variable is "rho" (to set flag for interpolation below)
-    if (var(i)(1:3).eq.'rho') then
+    if ((var(i)(1:3).eq.'rho').or.(var(i)(1:4).eq.'num_')) then
       j=1
       write(*,*) trim(var(i)), ' is interpolated in log scale' 
@@ -1922,5 +1958,5 @@
   ! interpolate dataset onto new grid
   ! check if variable is "rho" (to set flag for interpolation below)
-  if (var(i)(1:3).eq.'rho') then
+  if ((var(i)(1:3).eq.'rho').or.(var(i)(1:4).eq.'num_')) then
     j=1
     write(*,*) trim(var(i)), ' is interpolated in log scale' 
@@ -2839,5 +2875,6 @@
         else if (z_new(kloop).gt.z_gcm(iloop,jloop,altlen,tloop)) then
           ! z_new(kloop) is above highest GCM level
-          newdata(iloop,jloop,kloop,tloop)=gcmdata(iloop,jloop,altlen,tloop)
+        !  newdata(iloop,jloop,kloop,tloop)=gcmdata(iloop,jloop,altlen,tloop)
+          newdata(iloop,jloop,kloop,tloop)=missing
         else ! z_new(kloop) is within range
           x=z_new(kloop)
