Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 2288)
+++ /trunk/LMDZ.MARS/README	(revision 2289)
@@ -2950,2 +2950,8 @@
 Update reference run.def.64x48x49.MCD5 to keep up with recent (r2181) changes
 dissip_fac_mid and dissip_fac_up are now read in from the run.def file.
+
+== 16/04/2020 == AB
+Add the program solzenangle.F90 to the LMDZ.MARS/util/ directory, as well as a solzenangle.def
+This program actually replaces the former program terminator.F90, and can be used to interpolate GCM
+files at one given solar zenith angle (on the morning-side or the evening-side) all around the globe
+(the terminator corresponding to sza = 90deg)
Index: /trunk/LMDZ.MARS/util/solzenangle.F90
===================================================================
--- /trunk/LMDZ.MARS/util/solzenangle.F90	(revision 2289)
+++ /trunk/LMDZ.MARS/util/solzenangle.F90	(revision 2289)
@@ -0,0 +1,1541 @@
+program solzenangle
+
+! ----------------------------------------------------------------------------
+! Program to redistribute and interpolate the variables at the local time
+! corresponding to one solar zenith angle (morning or evening) everywhere
+! It can especially be used for the terminator (sza=90deg)
+! input : diagfi.nc  / concat.nc / stats.nc kind of files
+! authors: A.Bierjon,F.Forget,2020
+! ----------------------------------------------------------------------------
+
+implicit none
+
+include "netcdf.inc" ! NetCDF definitions
+
+character (len=50)  file
+! file(): input file(s) names(s)
+character (len=30), dimension(15) :: notconcat
+! notconcat(): names of the (15) variables that won't be concatenated
+character (len=50), dimension(:), allocatable :: var
+! var(): name(s) of variable(s) that will be concatenated
+character (len=50) :: tmpvar,title,units
+! tmpvar(): used to temporarily store a variable name
+! title(): [netcdf] title attribute
+! units(): [netcdf] units attribute
+character (len=100) :: filename,vartmp
+! filename(): output file name
+! vartmp(): temporary variable name (read from netcdf input file)
+!character (len=1) :: ccopy
+! ccpy: 'y' or 'n' answer
+
+character (len=50) :: altlong_name,altunits,altpositive
+! altlong_name(): [netcdf] altitude "long_name" attribute
+! altunits(): [netcdf] altitude "units" attribute
+! altpositive(): [netcdf] altitude "positive" attribute
+
+character (len=7) :: planetside
+! planetside(): planet side of the solar zenith angle ("morning" or "evening")
+
+integer :: nid,ierr,miss,validr
+! nid: [netcdf] file ID #
+! ierr: [netcdf] subroutine returned error code
+! miss: [netcdf] subroutine returned error code
+integer :: i,j,k,inter
+! for various loops
+integer :: varid
+! varid: [netcdf] variable ID #
+real, dimension(:), allocatable:: lat,lon,alt,ctl,time
+! lat(): array, stores latitude coordinates
+! lon(): array, stores longitude coordinates
+! alt(): array, stores altitude coordinates
+! ctl(): array, stores controle array
+! time(): array, stores time coordinates
+integer :: nbvar,nbvarfile,ndim
+! nbvar: # of variables to concatenate
+! nbfile: # number of input file(s)
+! nbvarfile: total # of variables in an input file
+! ndim: [netcdf] # (3 or 4) of dimensions (for variables)
+integer :: latdim,londim,altdim,ctldim,timedim
+! latdim: [netcdf] "latitude" dim ID
+! londim: [netcdf] "longitude" dim ID
+! altdim: [netcdf] "altdim" dim ID
+! ctldim: [netcdf] "controle" dim ID
+! timedim: [netcdf] "timedim" dim ID
+integer :: latvar,lonvar,altvar,ctlvar,timevar
+! latvar: [netcdf] ID of "latitude" variable
+! lonvar: [netcdf] ID of "longitude" variable
+! altvar: [netcdf] ID of "altitude" variable
+! ctlvar: [netcdf] ID of "controle" variable
+! timevar: [netcdf] ID of "Time" variable
+integer :: latlen,lonlen,altlen,ctllen,timelen,timelen_tot
+integer :: ilat,ilon,ialt,it
+! latlen: # of elements of lat() array
+! lonlen: # of elements of lon() array
+! altvar: # of elements of alt() array
+! ctlvar: # of elements of ctl() array
+! timelen: # of elements of time() array
+! timelen_tot:# =timelen or timelen+1 (if 1 more time to interpolate needed)
+integer :: nsol
+! nsol: # of sols in the input file AND # of elements of time() array in output
+integer :: nout,latdimout,londimout,altdimout,timedimout,timevarout
+! nout: [netcdf] output file ID
+! latdimout: [netcdf] output latitude (dimension) ID
+! londimout: [netcdf] output longitude (dimension) ID
+! altdimout: [netcdf] output altitude (dimension) ID
+! timedimout: [netcdf] output time (dimension) ID
+! timevarout: [netcdf] ID of output "Time" variable
+integer :: varidout 
+! varidout: [netcdf] variable ID # (of a variable to write to the output file)
+integer :: Nnotconcat,var_ok
+! Nnotconcat: # of (leading)variables that won't be concatenated
+! var_ok: flag (0 or 1)
+integer, dimension(4) :: corner,edges,dim
+! corner: [netcdf]
+! edges: [netcdf]
+! dim: [netcdf]
+real, dimension(:,:,:,:), allocatable :: var3d
+! var3D(,,,): 4D array to store a field
+
+real, dimension(:,:,:,:), allocatable :: var3d_lt
+! var3D_lt(,,,): 4D array to store a field in local time coordinate
+real, dimension(:), allocatable :: lt_gcm
+real, dimension(:), allocatable :: lt_outc
+real, dimension(:), allocatable :: var_gcm
+real, dimension(:), allocatable :: intsol
+real, dimension(:,:,:), allocatable :: lt_out
+real :: sza ! solar zenith angle
+character (len=30) :: szastr ! to store the sza value in a string
+real :: starttimeoffset ! offset (in sols) of sol 0 in input file, wrt Ls=0
+real :: tmpsol ! to temporarily store a sol value
+real :: tmpLs ! to temporarily store a Ls value
+real :: tmpLT ! to temporarily store a Local Time value
+
+real :: missing
+!PARAMETER(missing=1E+20)
+! missing: [netcdf] to handle "missing" values when reading/writing files
+real :: miss_lt
+PARAMETER(miss_lt=1E+20)
+! miss_lt: [netcdf] missing value to write for the Local Time in the output file
+real, dimension(2) :: valid_range
+! valid_range(2): [netcdf] interval in which a value is considered valid
+logical :: stats   ! stats=T when reading a "stats" kind of ile
+
+
+!==============================================================================
+! 1.1. Get input file name(s)
+!==============================================================================
+write(*,*) "Welcome in the solar zenith angle interpolator program !"
+write(*,*) 
+write(*,*) "which file do you want to use?  (diagfi... stats...  concat...)"
+
+read(*,'(a50)') file
+
+if (len_trim(file).eq.0) then
+   write(*,*) "no file... game over"
+   stop ""
+endif
+
+!==============================================================================
+! 1.2. Open the input file
+!==============================================================================
+
+ierr = NF_OPEN(file,NF_NOWRITE,nid)
+if (ierr.NE.NF_NOERR) then
+   write(*,*) 'ERROR: Pb opening file '//trim(file)
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+ierr=NF_INQ_NVARS(nid,nbvarfile)
+! nbvarfile now set to be the (total) number of variables in file
+if (ierr.NE.NF_NOERR) then
+   write(*,*) 'ERROR: Pb with NF_INQ_NVARS'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+!==============================================================================
+! 1.3. List of variables that should not be processed
+!==============================================================================
+
+notconcat(1)='Time'
+notconcat(2)='controle'
+notconcat(3)='rlonu'
+notconcat(4)='latitude'
+notconcat(5)='longitude'
+notconcat(6)='altitude'
+notconcat(7)='rlatv'
+notconcat(8)='aps'
+notconcat(9)='bps'
+notconcat(10)='ap'
+notconcat(11)='bp'
+notconcat(12)='cu'
+notconcat(13)='cv'
+notconcat(14)='aire'
+notconcat(15)='phisinit'
+
+!==============================================================================
+! 1.4. Get (and check) list of variables to process
+!==============================================================================
+write(*,*)
+   Nnotconcat=0
+do i=1,nbvarfile
+   ierr=NF_INQ_VARNAME(nid,i,vartmp)
+   ! vartmp now contains the "name" of variable of ID # i
+   var_ok=0
+   do inter=1,15
+      if (vartmp.eq.notconcat(inter)) then
+         var_ok=1 
+         Nnotconcat=Nnotconcat+1
+      endif 
+   enddo        
+   if (var_ok.eq.0)  write(*,*) trim(vartmp)
+enddo
+
+! Nnotconcat: # of variables that won't be processed
+! nbvarfile: total # of variables in file
+allocate(var(nbvarfile-Nnotconcat),stat=ierr)
+if (ierr.ne.0) then
+  write(*,*) "Error: failed allocation of var(nbvarfile-Nnotconcat)"
+  write(*,*) "  nbvarfile=",nbvarfile
+  write(*,*) "  Nnotconcat=",Nnotconcat
+  stop
+endif
+
+
+write(*,*)
+write(*,*) "which variables do you want to redistribute ?"
+write(*,*) "all / list of <variables> (separated by <Enter>s)"
+write(*,*) "(an empty line , i.e: just <Enter>, implies end of list)"
+nbvar=0
+read(*,'(a50)') tmpvar
+do while ((tmpvar/=' ').AND.(trim(tmpvar)/='all'))
+   nbvar=nbvar+1
+   var(nbvar)=tmpvar
+   read(*,'(a50)') tmpvar
+enddo
+
+if (tmpvar=="all") then
+         nbvar=nbvarfile-Nnotconcat
+         do j=Nnotconcat+1,nbvarfile
+            ierr=nf_inq_varname(nid,j,var(j-Nnotconcat))
+         enddo
+! Variables names from the file are stored in var()
+   nbvar=nbvarfile-Nnotconcat
+   do i=1,nbvar
+      ierr=nf_inq_varname(nid,i+Nnotconcat,var(i))
+      write(*,'(a9,1x,i2,1x,a1,1x,a50)') "variable ",i,":",var(i)
+   enddo
+else if(nbvar==0) then
+   write(*,*) "no variable... game over"
+   stop ""
+endif ! of if (tmpvar=="all")
+
+
+!==============================================================================
+! 2.1. Open input file
+!==============================================================================
+
+      write(*,*) 
+      write(*,*) "opening "//trim(file)//"..."
+      ierr = NF_OPEN(file,NF_NOWRITE,nid)
+      if (ierr.NE.NF_NOERR) then
+         write(*,*) 'ERROR: Pb opening file '//trim(file)
+         write(*,*) NF_STRERROR(ierr)
+         stop ""
+      endif
+
+!==============================================================================
+! 2.2. Read (and check) dimensions of variables from input file
+!==============================================================================
+
+   ierr=NF_INQ_DIMID(nid,"latitude",latdim)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'ERROR: Dimension <latitude> is missing in file '//trim(file)
+      stop ""  
+   endif
+   ierr=NF_INQ_VARID(nid,"latitude",latvar)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'ERROR: Field <latitude> is missing in file '//trim(file)
+      stop ""  
+   endif
+   ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
+!  write(*,*) "latlen: ",latlen
+
+   ierr=NF_INQ_DIMID(nid,"longitude",londim)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'ERROR: Dimension <longitude> is missing in file '//trim(file)
+      stop ""  
+   endif
+   ierr=NF_INQ_VARID(nid,"longitude",lonvar)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'ERROR: Field <longitude> is missing in file'//trim(file)
+      stop "" 
+   endif
+   ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
+!  write(*,*) "lonlen: ",lonlen
+
+   ierr=NF_INQ_DIMID(nid,"altitude",altdim)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'ERROR: Dimension <altitude> is missing in file '//trim(file)
+      stop ""  
+   endif
+   ierr=NF_INQ_VARID(nid,"altitude",altvar)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'ERROR: Field <altitude> is missing in file'//trim(file)
+      stop ""
+   endif
+   ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
+!  write(*,*) "altlen: ",altlen
+
+   ierr=NF_INQ_DIMID(nid,"index",ctldim)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) ' Dimension <index> is missing in file '//trim(file)
+      ctllen=0
+      !stop ""  
+   endif
+   ierr=NF_INQ_VARID(nid,"controle",ctlvar)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'Field <controle> is missing in file '//trim(file)
+      ctllen=0
+      !stop ""
+   else
+      ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
+   endif
+!  write(*,*) "controle: ",controle
+
+!==============================================================================
+! 2.3. Read (and check compatibility of) dimensions of
+!       variables from input file
+!==============================================================================
+
+! First call; initialize/allocate
+      allocate(lat(latlen),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate lat(latlen)"
+        stop
+      endif
+      allocate(lon(lonlen),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate lon(lonlen)"
+        stop
+      endif
+      allocate(alt(altlen),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate alt(altlen)"
+        stop
+      endif
+      allocate(ctl(ctllen),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate ctl(ctllen)"
+        stop
+      endif
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat)
+#else
+      ierr = NF_GET_VAR_REAL(nid,latvar,lat)
+#endif
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to load latitude"
+        write(*,*) NF_STRERROR(ierr)
+        stop
+      endif
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon)
+#else
+      ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
+#endif
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to load longitude"
+        write(*,*) NF_STRERROR(ierr)
+        stop
+      endif
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt)
+#else
+      ierr = NF_GET_VAR_REAL(nid,altvar,alt)
+#endif
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to load altitude"
+        write(*,*) NF_STRERROR(ierr)
+        stop
+      endif
+! Get altitude attributes to handle files with any altitude type
+ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
+ierr=nf_get_att_text(nid,altvar,'units',altunits)
+ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
+
+      if (ctllen .ne. 0) then
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
+#else
+        ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
+#endif
+        if (ierr.ne.0) then
+          write(*,*) "Error: failed to load controle"
+          write(*,*) NF_STRERROR(ierr)
+          stop
+        endif
+      endif ! of if (ctllen .ne. 0)
+
+!==============================================================================
+! 2.4. Handle "Time" dimension from input file
+!==============================================================================
+
+!==============================================================================
+! 2.4.0 Read "Time" dimension from input file
+!==============================================================================
+   ierr=NF_INQ_DIMID(nid,"Time",timedim)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'ERROR: Dimension <Time> is missing in file'//trim(file)
+      stop ""
+   endif
+   ierr=NF_INQ_VARID(nid,"Time",timevar)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'ERROR: Field <Time> is missing in file'//trim(file)
+      stop ""
+   endif
+   ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
+!  write(*,*) "timelen: ",timelen
+
+   ! Sanity Check: Does "Time" has a "title" attribute
+   ! and is it "Solar longitude" ?
+   ierr=nf_get_att_text(nid,timevar,"title",title)
+   if ((ierr.EQ.NF_NOERR).and.(title.eq."Solar longitude")) then
+     write(*,*) "ERROR: Time axis in input file is in Solar Longitude!"
+     write(*,*) "       localtime requires sols as time axis!"
+     write(*,*) " Might as well stop here."
+     stop
+   endif
+
+   ! allocate time() array and fill it with values from input file
+   allocate(time(timelen),stat=ierr)
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to allocate time(timelen)"
+     stop
+   endif
+
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
+#else
+   ierr = NF_GET_VAR_REAL(nid,timevar,time)
+#endif
+   if (ierr.NE.NF_NOERR) then
+     write(*,*) "Error , failed to load Time"
+     write(*,*) NF_STRERROR(ierr)
+     stop
+   endif
+
+
+!==============================================================================
+! 2.4.1 Select local times to be saved "Time" in output file
+!==============================================================================
+      write(*,*) "Planet side of the Solar Zenith Angle ?"
+      read(*,*) planetside
+      if ((trim(planetside).ne."morning").and.(trim(planetside).ne."evening")) then
+        write(*,*) "planetside = "//planetside
+        write(*,*) "Error: the planet side must be morning or evening"
+        stop
+      endif
+      write(*,*) "Solar Zenith angle to interpolate ?"
+      write(*,*) "(between 0 and 180 deg, 90deg for the terminator at the surface)"
+      read(*,*) sza
+      if ((sza.lt.0).or.(sza.ge.180)) then
+        write(*,*) "ERROR: the solar zenith angle must lie within [0;180[ deg"
+        stop
+      endif
+      write(*,*) "Beginning date of the simulation file "
+      write(*,*) "(or median date of the simulation for a stats)?"
+      write(*,*) "(i.e. number of sols since Ls=0 at the Time=0.0 in the input file)"
+      read(*,*) starttimeoffset
+
+      if ((timelen.eq.12).and.(time(1).eq.2).and.(time(timelen).eq.24))then
+         write(*,*) 'We detect a ""stats"" file'
+         stats=.true.
+!        We rewrite the time from "stats" from 0 to 1 sol...
+         do it=1,timelen  ! loop temps in
+               time(it) = time(it)/24.
+         end do
+         nsol =1
+      else   ! case of a diagfi or concatnc file 
+        stats=.false.
+!       Number of sol in input file
+        nsol = int(time(timelen)) - int(time(1))
+      end if
+      
+      allocate(intsol(nsol),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate intsol(nsol)"
+        stop
+      endif
+      allocate(lt_out(lonlen,latlen,nsol),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate lt_out(lonlen,latlen,nsol)"
+        stop
+      endif
+      do it=1,nsol
+        intsol(it)= int(time(1)) + it-1.
+        do ilon=1,lonlen
+!         For the computation of Ls, it's supposed that the morning/evening
+!         correspond to LT=6h/18h at the given longitude, which is
+!         then reported to a sol value at longitude 0
+          if (trim(planetside).eq."morning") then
+            tmpsol = intsol(it) + (6.-lon(ilon)/15.)/24. + starttimeoffset
+          else !if trim(planetside).eq."evening"
+            tmpsol = intsol(it) + (18.-lon(ilon)/15.)/24. + starttimeoffset            
+          endif
+          call sol2ls(tmpsol,tmpLs)
+          do ilat=1,latlen
+!           Compute the Local Time of the solar zenith angle at a given longitude
+            call LTsolzenangle(lat(ilat),tmpLs,planetside,sza,miss_lt,tmpLT)
+!           Deduce the sol value at this longitude
+            if (tmpLT.eq.miss_lt) then 
+              ! if there is no LT corresponding to the sza, put a missing value in lt_out
+              lt_out(ilon,ilat,it)=miss_lt
+            else
+              lt_out(ilon,ilat,it)=tmpLT
+            endif
+          enddo
+        enddo
+      enddo
+
+      if (nsol.gt.1) then
+         timelen_tot=timelen
+      else
+!        if only 1 sol available, we must add 1 timestep for the interpolation
+         timelen_tot=timelen+1
+      endif      
+      allocate(lt_gcm(timelen_tot),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate lt_gcm(timelen_tot)"
+        stop
+      endif
+      allocate(var_gcm(timelen_tot),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate var_gcm(timelen_tot)"
+        stop
+      endif
+      allocate(lt_outc(nsol),stat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) "Error: failed to allocate lt_outc(nsol)"
+        stop
+      endif
+
+!==============================================================================
+! 2.4.2. Get output file name
+!==============================================================================
+    if (trim(planetside).eq."morning") then 
+      filename=file(1:len_trim(file)-3)//"_MO.nc"
+    else ! if trim(planetside).eq."evening"
+      filename=file(1:len_trim(file)-3)//"_EV.nc"
+    endif
+
+!==============================================================================
+! 2.4.3. Initiate dimension in output file
+!==============================================================================
+
+
+   ! Initialize output file's lat,lon,alt and time dimensions
+    call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
+          altlong_name,altunits,altpositive,&
+          nout,latdimout,londimout,altdimout,timedimout,timevarout)
+
+   ! Initialize output file's aps,bps and phisinit variables
+    call init2(nid,lonlen,latlen,altlen,altdim,&
+               nout,londimout,latdimout,altdimout)
+
+
+!==============================================================================
+! 2.4.4. Write/extend "Time" dimension/values in output file
+!==============================================================================
+
+   if (stats) then 
+
+     do it=1,nsol
+#ifdef NC_DOUBLE
+        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,intsol(it)*24.)
+#else
+        ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,intsol(it)*24.)
+#endif
+     enddo
+   else
+     do it=1,nsol
+#ifdef NC_DOUBLE
+        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,intsol(it))
+#else
+        ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,intsol(it))
+#endif
+        if (ierr.NE.NF_NOERR) then
+          write(*,*) "Error , failed to write Time"
+          stop
+        endif
+     enddo
+   end if
+
+!==============================================================================
+! 2.5 Read/write variables
+!==============================================================================
+
+   do j=1,nbvar ! loop on variables to read/write
+
+!==============================================================================
+! 2.5.1 Check that variable to be read existe in input file
+!==============================================================================
+
+      write(*,*) "variable ",var(j)
+      ierr=nf_inq_varid(nid,var(j),varid)
+      if (ierr.NE.NF_NOERR) then
+         write(*,*) 'ERROR: Field <',var(j),'> not found in file'//file
+         stop ""
+      endif
+      ierr=nf_inq_varndims(nid,varid,ndim)
+
+!==============================================================================
+! 2.5.2 Prepare things in order to read/write the variable
+!==============================================================================
+
+      ! build dim(),corner() and edges() arrays
+      ! and allocate var3d() array
+      if (ndim==3) then
+         allocate(var3d(lonlen,latlen,timelen,1))
+         allocate(var3d_lt(lonlen,latlen,nsol,1))
+         dim(1)=londimout
+         dim(2)=latdimout
+         dim(3)=timedimout
+
+         ! start indexes (where data values will be written)
+         corner(1)=1
+         corner(2)=1
+         corner(3)=1
+         corner(4)=1
+
+         ! length (along dimensions) of block of data to be written
+         edges(1)=lonlen 
+         edges(2)=latlen 
+         edges(3)=nsol
+         edges(4)=1
+   
+      else if (ndim==4) then
+         allocate(var3d(lonlen,latlen,altlen,timelen))
+         allocate(var3d_lt(lonlen,latlen,altlen,nsol))
+         dim(1)=londimout
+         dim(2)=latdimout
+         dim(3)=altdimout
+         dim(4)=timedimout
+
+         ! start indexes (where data values will be written)
+         corner(1)=1
+         corner(2)=1
+         corner(3)=1
+         corner(4)=1
+
+         ! length (along dimensions) of block of data to be written
+         edges(1)=lonlen
+         edges(2)=latlen
+         edges(3)=altlen
+         edges(4)=nsol
+      endif
+
+         units=" "
+         title=" "
+         ierr=nf_get_att_text(nid,varid,"title",title)
+         ierr=nf_get_att_text(nid,varid,"units",units)
+         call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr)
+
+      
+
+!==============================================================================
+! 2.5.3 Read from input file 
+!==============================================================================
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid,varid,var3d)
+      miss=NF_GET_ATT_DOUBLE(nid,varid,"missing_value",missing)
+      validr=NF_GET_ATT_DOUBLE(nid,varid,"valid_range",valid_range)
+#else
+      ierr = NF_GET_VAR_REAL(nid,varid,var3d)
+      miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
+      validr=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
+#endif
+
+!==============================================================================
+! 2.5.4 Interpolation in local time 
+!==============================================================================
+
+        do ilon=1,lonlen
+!          write(*,*) 'processing longitude =', lon(ilon)
+!          Local time at each longitude :
+           do it=1,timelen  ! loop temps in
+               lt_gcm(it) = time(it) + lon(ilon)/360.
+           enddo
+           if (nsol.eq.1) lt_gcm(timelen+1) = lt_gcm(1)+1
+
+           do ilat=1,latlen
+             do it=1,nsol ! loop time local time
+               if (lt_out(ilon,ilat,it).eq.miss_lt) then
+                 lt_outc(it)=miss_lt
+               else
+                 lt_outc(it)=intsol(it) + lt_out(ilon,ilat,it)/24.
+               endif
+!              If the computed local time exceeds the GCM input file time bounds,
+!              put a missing value in the local time variable
+               if ((lt_outc(it).gt.lt_gcm(timelen_tot)).OR.(lt_outc(it).lt.lt_gcm(1))) then
+                 lt_out(ilon,ilat,it)=miss_lt
+                 lt_outc(it)=miss_lt
+               endif
+             enddo ! local time
+
+             if (ndim==3) then
+               do it=1,timelen  ! loop temps in
+                  var_gcm(it) = var3d(ilon,ilat,it,1)
+               enddo
+               if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
+               do it=1,nsol ! loop time local time
+                 if (lt_outc(it).eq.miss_lt) then
+                   var3d_lt(ilon,ilat,it,1)=missing
+                 else
+                   call interpolf(lt_outc(it),var3d_lt(ilon,ilat,it,1),&
+                                 missing,lt_gcm,var_gcm,timelen_tot)
+                 endif
+               enddo ! local time
+
+             else if  (ndim==4) then
+               do ialt=1,altlen
+                 do it=1,timelen ! loop temps in
+                    var_gcm(it) = var3d(ilon,ilat,ialt,it)
+                 enddo
+                 if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
+                 do it=1,nsol ! loop time local time
+                   if (lt_outc(it).eq.miss_lt) then
+                     var3d_lt(ilon,ilat,ialt,it)=missing
+                   else
+                     call interpolf(lt_outc(it),var3d_lt(ilon,ilat,ialt,it),&
+                                   missing,lt_gcm,var_gcm,timelen_tot)
+                   endif
+                 enddo ! local time   
+               enddo ! alt
+             endif
+             
+           enddo ! lat
+        end do ! lon
+
+
+!==============================================================================
+! 2.5.5 write (append) to the output file
+!==============================================================================
+
+#ifdef NC_DOUBLE
+      ierr= NF_PUT_VARA_DOUBLE(nout,varidout,corner,edges,var3d_lt)
+#else
+      ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3d_lt)
+#endif
+
+      if (ierr.ne.NF_NOERR) then
+         write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
+         stop ""
+      endif
+
+! In case there is a "valid_range" attribute
+      ! Write "valid_range" and "missing_value" attributes in output file
+      if (miss.eq.NF_NOERR) then
+         call missing_value(nout,varidout,validr,miss,valid_range,missing)
+      endif
+
+      ! free var3d() array
+      deallocate(var3d)
+      deallocate(var3d_lt)
+
+   enddo ! of do j=1,nbvar
+
+!==============================================================================
+! 2.5.6 Write the solar zenith angle Local Time variable in output file
+!==============================================================================
+   write(*,*) "variable LTsolzenangle"
+   ! build dim(),corner() and edges() arrays
+   dim(1)=londimout
+   dim(2)=latdimout
+   dim(3)=timedimout
+
+   ! start indexes (where data values will be written)
+   corner(1)=1
+   corner(2)=1
+   corner(3)=1
+   corner(4)=1
+
+   ! length (along dimensions) of block of data to be written
+   edges(1)=lonlen 
+   edges(2)=latlen 
+   edges(3)=nsol
+   edges(4)=1
+
+   tmpvar="LTsolzenangle"
+   units="hours"
+   write(szastr,*) sza
+   if (trim(planetside).eq."morning") then
+     title="Morning-side LTST at sza="//trim(adjustl(szastr))//"deg"
+   else !if trim(planetside).eq."evening"
+     title="Evening-side LTST at sza="//trim(adjustl(szastr))//"deg"            
+   endif
+   call def_var(nout,tmpvar,title,units,3,dim,varidout,ierr)
+
+#ifdef NC_DOUBLE
+   ierr= NF_PUT_VARA_DOUBLE(nout,varidout,corner,edges,lt_out)
+#else
+   ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,lt_out)
+#endif
+   if (ierr.ne.NF_NOERR) then
+      write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
+      stop ""
+   endif
+
+   validr=NF_NOERR+1
+   miss=NF_NOERR
+   call missing_value(nout,varidout,validr,miss,(0,0),miss_lt)
+   write(*,*)"with missing value = ",miss_lt
+!==============================================================================
+! 3. END
+!==============================================================================
+   deallocate(time)
+   deallocate(intsol)
+   deallocate(lt_gcm)
+   deallocate(lt_out)
+   deallocate(lt_outc)
+   deallocate(var_gcm)
+
+   ! Close input file
+   ierr=nf_close(nid)
+
+! Close output file
+ierr=NF_CLOSE(nout)
+
+write(*,*) "Program completed !"
+
+end program solzenangle
+
+
+
+!******************************************************************************
+!******************************************************************************
+
+subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
+                     altlong_name,altunits,altpositive,&
+                     nout,latdimout,londimout,altdimout,timedimout,timevarout)
+!==============================================================================
+! Purpose:
+! Create and initialize a data file (NetCDF format)
+!==============================================================================
+! Remarks:
+! The NetCDF file (created in this subroutine) remains open
+!==============================================================================
+
+implicit none
+
+include "netcdf.inc" ! NetCDF definitions
+
+!==============================================================================
+! Arguments:
+!==============================================================================
+character (len=*), intent(in):: filename
+! filename(): the file's name
+integer,intent(in) ::latlen,lonlen,altlen,ctllen
+real, intent(in):: lat(latlen)
+! lat(): latitude
+real, intent(in):: lon(lonlen)
+! lon(): longitude
+real, intent(in):: alt(altlen)
+! alt(): altitude
+real, intent(in):: ctl(ctllen)
+! ctl(): controle
+character (len=*), intent(in) :: altlong_name
+! altlong_name(): [netcdf] altitude "long_name" attribute
+character (len=*), intent(in) :: altunits
+! altunits(): [netcdf] altitude "units" attribute
+character (len=*), intent(in) :: altpositive
+! altpositive(): [netcdf] altitude "positive" attribute
+integer, intent(out):: nout
+! nout: [netcdf] file ID
+integer, intent(out):: latdimout
+! latdimout: [netcdf] lat() (i.e.: latitude)  ID
+integer, intent(out):: londimout
+! londimout: [netcdf] lon()  ID
+integer, intent(out):: altdimout
+! altdimout: [netcdf] alt()  ID
+integer, intent(out):: timedimout
+! timedimout: [netcdf] "Time"  ID
+integer, intent(out):: timevarout
+! timevarout: [netcdf] "Time" (considered as a variable) ID
+
+!==============================================================================
+! Local variables:
+!==============================================================================
+!integer :: latdim,londim,altdim,timedim
+integer :: nvarid,ierr
+integer :: ctldimout
+! nvarid: [netcdf] ID of a variable
+! ierr: [netcdf]  return error code (from called subroutines)
+
+!==============================================================================
+! 1. Create (and open) output file
+!==============================================================================
+write(*,*) "creating "//trim(adjustl(filename))//'...'
+ierr = NF_CREATE(filename,IOR(NF_CLOBBER,NF_64BIT_OFFSET),nout)
+! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'ERROR: Impossible to create the file ',trim(filename)
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+!==============================================================================
+! 2. Define/write "dimensions" and get their IDs
+!==============================================================================
+
+ierr = NF_DEF_DIM(nout, "latitude", latlen, latdimout)
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed to define dimension <latitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+ierr = NF_DEF_DIM(nout, "longitude", lonlen, londimout)
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed to define dimension <longitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+ierr = NF_DEF_DIM(nout, "altitude", altlen, altdimout)
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed to define dimension <altitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+if (size(ctl).ne.0) then
+  ierr = NF_DEF_DIM(nout, "index", ctllen, ctldimout)
+  if (ierr.NE.NF_NOERR) then
+    WRITE(*,*)'initiate: error failed to define dimension <index>'
+    write(*,*) NF_STRERROR(ierr)
+    stop ""
+  endif
+endif
+ierr = NF_DEF_DIM(nout, "Time", NF_UNLIMITED, timedimout)
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed to define dimension <Time>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+! End netcdf define mode
+ierr = NF_ENDDEF(nout)
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed to switch out of define mode'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+!==============================================================================
+! 3. Write "Time" and its attributes
+!==============================================================================
+
+call def_var(nout,"Time","Time","years since 0000-00-0 00:00:00",1,&
+             (/timedimout/),timevarout,ierr)
+
+!==============================================================================
+! 4. Write "latitude" (data and attributes)
+!==============================================================================
+
+call def_var(nout,"latitude","latitude","degrees_north",1,&
+             (/latdimout/),nvarid,ierr)
+
+#ifdef NC_DOUBLE
+ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lat)
+#else
+ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
+#endif
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed writing variable <latitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+!==============================================================================
+! 4. Write "longitude" (data and attributes)
+!==============================================================================
+
+call def_var(nout,"longitude","East longitude","degrees_east",1,&
+             (/londimout/),nvarid,ierr)
+
+#ifdef NC_DOUBLE
+ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lon)
+#else
+ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
+#endif
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed writing variable <longitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+!==============================================================================
+! 5. Write "altitude" (data and attributes)
+!==============================================================================
+
+! Switch to netcdf define mode
+ierr = NF_REDEF (nout)
+
+#ifdef NC_DOUBLE
+ierr = NF_DEF_VAR (nout,"altitude",NF_DOUBLE,1,altdimout,nvarid)
+#else
+ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,altdimout,nvarid)
+#endif
+
+ierr = NF_PUT_ATT_TEXT (nout,nvarid,'long_name',len_trim(adjustl(altlong_name)),adjustl(altlong_name))
+ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',len_trim(adjustl(altunits)),adjustl(altunits))
+ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',len_trim(adjustl(altpositive)),adjustl(altpositive))
+
+! End netcdf define mode
+ierr = NF_ENDDEF(nout)
+
+#ifdef NC_DOUBLE
+ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,alt)
+#else
+ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
+#endif 
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed writing variable <altitude>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+!==============================================================================
+! 6. Write "controle" (data and attributes)
+!==============================================================================
+
+if (size(ctl).ne.0) then
+   ! Switch to netcdf define mode
+   ierr = NF_REDEF (nout)
+
+#ifdef NC_DOUBLE
+   ierr = NF_DEF_VAR (nout,"controle",NF_DOUBLE,1,ctldimout,nvarid)
+#else
+   ierr = NF_DEF_VAR (nout,"controle",NF_FLOAT,1,ctldimout,nvarid)
+#endif
+
+   ierr = NF_PUT_ATT_TEXT (nout,nvarid,"title",18,"Control parameters")
+
+   ! End netcdf define mode
+   ierr = NF_ENDDEF(nout)
+
+#ifdef NC_DOUBLE
+   ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,ctl)
+#else
+   ierr = NF_PUT_VAR_REAL (nout,nvarid,ctl)
+#endif
+   if (ierr.NE.NF_NOERR) then
+      WRITE(*,*)'initiate: error failed writing variable <controle>'
+      write(*,*) NF_STRERROR(ierr)
+      stop ""
+   endif
+endif
+
+end Subroutine initiate
+
+!******************************************************************************
+subroutine init2(infid,lonlen,latlen,altlen,altdimid, &
+                 outfid,londimout,latdimout,altdimout)
+!==============================================================================
+! Purpose:
+! Copy aps(), bps() and phisinit() from input file to output file
+!==============================================================================
+! Remarks:
+! The NetCDF files must be open
+!==============================================================================
+
+implicit none
+
+include "netcdf.inc" ! NetCDF definitions
+
+!==============================================================================
+! Arguments:
+!==============================================================================
+integer, intent(in) :: infid  ! NetCDF output file ID
+integer, intent(in) :: lonlen ! # of grid points along longitude
+integer, intent(in) :: latlen ! # of grid points along latitude
+integer, intent(in) :: altlen ! # of grid points along altitude
+integer, intent(in) :: altdimid ! ID of altitude dimension
+integer, intent(in) :: outfid ! NetCDF output file ID
+integer, intent(in) :: londimout ! longitude dimension ID
+integer, intent(in) :: latdimout ! latitude dimension ID
+integer, intent(in) :: altdimout ! altitude dimension ID
+!==============================================================================
+! Local variables:
+!==============================================================================
+real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
+real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
+integer :: apsid,bpsid,phisinitid
+integer :: ierr
+integer :: tmpdimid ! temporary dimension ID
+integer :: tmpvarid ! temporary variable ID
+logical :: phis, aps_ok, bps_ok ! are "phisinit" "aps" "bps" available ?
+
+
+!==============================================================================
+! 1. Read data from input file
+!==============================================================================
+
+! hybrid coordinate aps
+  allocate(aps(altlen),stat=ierr)
+  if (ierr.ne.0) then
+    write(*,*) "init2: failed to allocate aps(altlen)"
+    stop
+  endif
+
+ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
+if (ierr.ne.NF_NOERR) then
+  write(*,*) "oops Failed to get aps ID. OK"
+  aps_ok=.false.
+else
+  ! Check that aps() is the right size (it most likely won't be if
+  ! zrecast has be used to generate the input file)
+  ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
+  if (tmpdimid.ne.altdimid) then
+    write(*,*) "init2: wrong dimension size for aps(), skipping it ..."
+    aps_ok=.false.
+  else
+    ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
+    if (ierr.ne.NF_NOERR) then
+     stop "init2 error: Failed reading aps"
+    endif
+    aps_ok=.true.
+  endif
+endif
+
+! hybrid coordinate bps
+  allocate(bps(altlen),stat=ierr)
+  if (ierr.ne.0) then
+    write(*,*) "init2: failed to allocate bps(altlen)"
+    stop
+  endif
+
+ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
+if (ierr.ne.NF_NOERR) then
+  write(*,*) "oops: Failed to get bps ID. OK"
+  bps_ok=.false.
+else
+  ! Check that bps() is the right size
+  ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
+  if (tmpdimid.ne.altdimid) then
+    write(*,*) "init2: wrong dimension size for bps(), skipping it ..."
+    bps_ok=.false.
+  else
+    ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
+    if (ierr.ne.NF_NOERR) then
+      stop "init2 Error: Failed reading bps"
+    endif
+    bps_ok=.true.
+  endif
+endif
+
+! ground geopotential phisinit
+allocate(phisinit(lonlen,latlen),stat=ierr)
+if (ierr.ne.0) then
+  write(*,*) "init2: failed to allocate phisinit(lonlen,latlen)"
+  stop
+endif
+ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
+if (ierr.ne.NF_NOERR) then
+  write(*,*) "Failed to get phisinit ID. OK"
+  phisinit = 0.
+  phis = .false.
+else
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
+  if (ierr.ne.NF_NOERR) then
+    stop "init2 Error: Failed reading phisinit"
+  endif
+  phis = .true.
+endif
+
+
+!==============================================================================
+! 2. Write
+!==============================================================================
+
+!==============================================================================
+! 2.2. Hybrid coordinates aps() and bps()
+!==============================================================================
+
+IF (aps_ok) then 
+  ! define aps
+  call def_var(outfid,"aps","hybrid pressure at midlayers"," ",1,&
+               (/altdimout/),apsid,ierr)
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed to def_var aps"
+  endif
+
+  ! write aps
+#ifdef NC_DOUBLE
+  ierr=NF_PUT_VAR_DOUBLE(outfid,apsid,aps)
+#else
+  ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
+#endif
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed to write aps"
+  endif
+ENDIF ! of IF (aps_ok)
+
+IF (bps_ok) then 
+  ! define bps
+  call def_var(outfid,"bps","hybrid sigma at midlayers"," ",1,&
+               (/altdimout/),bpsid,ierr)
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed to def_var bps"
+  endif
+
+  ! write bps
+#ifdef NC_DOUBLE
+  ierr=NF_PUT_VAR_DOUBLE(outfid,bpsid,bps)
+#else
+  ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
+#endif
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed to write bps"
+  endif
+ENDIF ! of IF (bps_ok)
+
+!==============================================================================
+! 2.2. phisinit()
+!==============================================================================
+
+IF (phis) THEN
+
+  !define phisinit
+  call def_var(outfid,"phisinit","Ground level geopotential"," ",2,&
+              (/londimout,latdimout/),phisinitid,ierr)
+  if (ierr.ne.NF_NOERR) then
+     stop "Error: Failed to def_var phisinit"
+  endif
+
+  ! write phisinit
+#ifdef NC_DOUBLE
+  ierr=NF_PUT_VAR_DOUBLE(outfid,phisinitid,phisinit)
+#else
+  ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
+#endif
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed to write phisinit"
+  endif
+
+ENDIF ! of IF (phis)
+
+
+! Cleanup
+deallocate(aps)
+deallocate(bps)
+deallocate(phisinit)
+
+end subroutine init2
+
+!******************************************************************************
+subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
+!==============================================================================
+! Purpose: Write a variable (i.e: add a variable to a dataset)
+! called "name"; along with its attributes "title", "units"...
+! to a file (following the NetCDF format)
+!==============================================================================
+! Remarks:
+! The NetCDF file must be open
+!==============================================================================
+
+implicit none
+
+include "netcdf.inc" ! NetCDF definitions
+
+!==============================================================================
+! Arguments:
+!==============================================================================
+integer, intent(in) :: nid
+! nid: [netcdf] file ID #
+character (len=*), intent(in) :: name
+! name(): [netcdf] variable's name
+character (len=*), intent(in) :: title
+! title(): [netcdf] variable's "title" attribute
+character (len=*), intent(in) :: units
+! unit(): [netcdf] variable's "units" attribute
+integer, intent(in) :: nbdim
+! nbdim: number of dimensions of the variable
+integer, dimension(nbdim), intent(in) :: dim
+! dim(nbdim): [netcdf] dimension(s) ID(s)
+integer, intent(out) :: nvarid
+! nvarid: [netcdf] ID # of the variable
+integer, intent(out) :: ierr
+! ierr: [netcdf] subroutines returned error code
+
+! Switch to netcdf define mode
+ierr=NF_REDEF(nid)
+
+! Insert the definition of the variable
+#ifdef NC_DOUBLE
+ierr=NF_DEF_VAR(nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
+#else
+ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
+#endif
+
+! Write the attributes
+ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",len_trim(adjustl(title)),adjustl(title))
+ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
+
+! End netcdf define mode
+ierr=NF_ENDDEF(nid)
+
+end subroutine def_var
+
+!******************************************************************************
+subroutine  missing_value(nout,nvarid,validr,miss,valid_range,missing)
+!==============================================================================
+! Purpose:
+! Write "valid_range" and "missing_value" attributes (of a given
+! variable) to a netcdf file
+!==============================================================================
+! Remarks:
+! NetCDF file must be open
+! Variable (nvarid) ID must be set
+!==============================================================================
+
+implicit none
+
+include "netcdf.inc"  ! NetCDF definitions
+
+!==============================================================================
+! Arguments:
+!==============================================================================
+integer, intent(in) :: nout
+! nout: [netcdf] file ID #
+integer, intent(in) :: nvarid
+! varid: [netcdf] variable ID #
+integer, intent(in) :: validr
+! validr : [netcdf] routines return code for "valid_range" attribute
+integer, intent(in) :: miss
+! miss : [netcdf] routines return code for "missing_value" attribute
+real, dimension(2), intent(in) :: valid_range
+! valid_range(2): [netcdf] "valid_range" attribute (min and max)
+real, intent(in) :: missing
+! missing: [netcdf] "missing_value" attribute
+
+!==============================================================================
+! Local variables:
+!==============================================================================
+integer :: ierr
+! ierr: [netcdf] subroutine returned error code
+!      INTEGER nout,nvarid,ierr
+!      REAL missing
+!      REAL valid_range(2)
+
+! Switch to netcdf dataset define mode
+ierr = NF_REDEF (nout)
+if (ierr.ne.NF_NOERR) then
+   print*,'missing_value: '
+   print*, NF_STRERROR(ierr)
+endif
+
+! Write valid_range() attribute
+if (validr.eq.NF_NOERR) then
+#ifdef NC_DOUBLE
+  ierr=NF_PUT_ATT_DOUBLE(nout,nvarid,'valid_range',NF_DOUBLE,2,valid_range)
+#else
+  ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
+#endif
+
+  if (ierr.ne.NF_NOERR) then
+     print*,'missing_value: valid_range attribution failed'
+     print*, NF_STRERROR(ierr)
+     !write(*,*) 'NF_NOERR', NF_NOERR
+     !CALL abort
+     stop ""
+  endif
+endif
+
+! Write "missing_value" attribute
+if (miss.eq.NF_NOERR) then
+#ifdef NC_DOUBLE
+  ierr= NF_PUT_ATT_DOUBLE(nout,nvarid,'missing_value',NF_DOUBLE,1,missing)
+#else
+  ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
+#endif
+
+  if (ierr.NE.NF_NOERR) then
+     print*, 'missing_value: missing value attribution failed'
+     print*, NF_STRERROR(ierr)
+  !    WRITE(*,*) 'NF_NOERR', NF_NOERR
+  !    CALL abort
+     stop ""
+  endif
+endif
+
+! End netcdf dataset define mode
+ierr = NF_ENDDEF(nout)
+
+end subroutine  missing_value
+
+!*****************************************************************************
+subroutine interpolf(x,y,missing,xd,yd,nd)
+!==============================================================================
+! Purpose:
+! Yield y=f(x), where xd() end yd() are arrays of known values,
+! using linear interpolation
+! If x is not included in the interval spaned by xd(), then y is set
+! to a default value 'missing'
+! Note:
+! Array xd() should contain ordered (either increasing or decreasing) abscissas
+!==============================================================================
+implicit none
+!==============================================================================
+! Arguments:
+!==============================================================================
+real,intent(in) :: x ! abscissa at which interpolation is to be done
+real,intent(in) :: missing ! missing value (if no interpolation is performed)
+integer,intent(in) :: nd ! size of arrays
+real,dimension(nd),intent(in) :: xd ! array of known absissas
+real,dimension(nd),intent(in) :: yd ! array of correponding values
+
+real,intent(out) :: y ! interpolated value
+!==============================================================================
+! Local variables:
+!==============================================================================
+integer :: i
+
+! default: set y to 'missing'
+y=missing
+
+   do i=1,nd-1
+     if (((x.ge.xd(i)).and.(x.le.xd(i+1))).or.&
+          ((x.le.xd(i)).and.(x.ge.xd(i+1)))) then
+        if ((yd(i).eq.missing).or.(yd(i+1).eq.missing)) then
+          ! cannot perform the interpolation if an encompassing value
+          ! is already set to 'missing'
+        else
+          !linear interpolation based on encompassing values
+          y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i))
+        endif
+        exit
+     endif
+   enddo
+
+
+end subroutine interpolf
+
+!******************************************************************************
+subroutine LTsolzenangle(lat,ls,planetside,sza,missval,lt_sza)
+! compute the localtime by inversing
+! the solar zenith angle equation :
+! cos(sza)=sin(lat)*sin(declin)
+!         +cos(lat)*cos(declin)*cos(2.*pi*(localtime/24.-.5))
+
+
+implicit none
+
+real,intent(in) :: lat ! latitude (deg)
+real,intent(in) :: ls ! solar longitude (deg)
+character (len=7), intent(in) :: planetside ! planet side : "morning" or "evening"
+real,intent(in) :: sza ! solar zenith angle (deg) (within [0;180[)
+real,intent(in) :: missval ! missing value for when there is no associated LT (polar day/night)
+real,intent(out) :: lt_sza ! local true solar time (hours) corresponding to the sza
+
+! local variables:
+double precision :: declin ! declination of the Sun (rad)
+double precision :: tmpcos ! to temporarily store a cosine value
+double precision,parameter :: obliquity=25.1919d0
+double precision,parameter :: pi=3.14159265358979d0
+double precision,parameter :: degtorad=pi/180.0d0
+double precision,parameter :: radtodeg=180.0d0/pi
+
+! Compute Sun's declination
+declin=ASIN(sin(ls*degtorad)*sin(obliquity*degtorad))
+
+! Small checks
+if ((cos(lat*degtorad).eq.0).or.(cos(declin).eq.0)) then
+  write(*,*) "ERROR in LTsolzenangle : lat and declin can't have their cosine equal to zero"
+  stop
+endif
+
+! Compute Local Time
+tmpcos = cos(sza*degtorad)/(cos(lat*degtorad)*cos(declin)) - tan(lat*degtorad)*tan(declin)
+
+if ((tmpcos.gt.1.).or.(tmpcos.lt.(-1.))) then
+  ! Detect Polar Day and Polar Night (values out of the arccos definition domain)
+  lt_sza = missval
+  
+else
+
+  if (trim(planetside).eq."morning") then
+  ! Local Time of the Morning-side Solar Zenith Angle
+    lt_sza = 24*(0.5 - ACOS(tmpcos)/(2.*pi))
+  else ! if trim(planetside).eq."evening"
+  ! Local Time of the Evening-side Solar Zenith Angle
+    lt_sza = 24*(0.5 + ACOS(tmpcos)/(2.*pi))
+  endif
+  
+endif
+
+end subroutine LTsolzenangle
+
+!******************************************************************************
+subroutine sol2ls(sol,ls)
+!  convert a given martian day number (sol)
+!  into corresponding solar longitude, Ls (in degr.),
+!  where sol=0=Ls=0 is the
+!  northern hemisphere spring equinox.
+
+implicit none
+
+!input 
+real,intent(in) :: sol
+!output 
+real,intent(out) :: ls
+
+!local variables
+double precision :: year_day,peri_day,timeperi,e_elips
+parameter (year_day=668.6d0) ! number of martian days (sols) in a martian year
+parameter (peri_day=485.35d0) ! perihelion date (in sols)
+parameter (e_elips=0.09340d0) ! orbital eccentricity
+parameter (timeperi=1.90258341759902d0) ! timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
+
+double precision :: pi,radtodeg
+parameter (pi=3.14159265358979d0)
+parameter (radtodeg=57.2957795130823d0) !  radtodeg: 180/pi
+
+double precision :: zanom,xref,zx0,zdx,zteta,zz
+!  xref: mean anomaly, zx0: eccentric anomaly, zteta: true anomaly        
+integer :: iter
+
+zz=(sol-peri_day)/year_day
+zanom=2.*pi*(zz-nint(zz))
+xref=dabs(zanom)
+
+!  The equation zx0 - e * sin (zx0) = xref, solved by Newton
+zx0=xref+e_elips*dsin(xref)
+iter=1
+iteration:do while (iter.le.10)
+   zdx=-(zx0-e_elips*dsin(zx0)-xref)/(1.-e_elips*dcos(zx0))
+   if(dabs(zdx).le.(1.d-7)) then ! typically, 2 or 3 iterations are enough
+     EXIT iteration 
+   endif
+   zx0=zx0+zdx
+   iter=iter+1
+enddo iteration
+zx0=zx0+zdx
+if(zanom.lt.0.) zx0=-zx0
+
+! compute true anomaly zteta, now that eccentric anomaly zx0 is known
+zteta=2.*datan(dsqrt((1.+e_elips)/(1.-e_elips))*dtan(zx0/2.))
+
+! compute Ls
+ls=real(zteta-timeperi)
+if(ls.lt.0.) ls=ls+2.*real(pi)
+if(ls.gt.2.*pi) ls=ls-2.*real(pi)
+! convert Ls in deg.
+ls=real(radtodeg)*ls
+
+end subroutine sol2ls
+
+!******************************************************************************
Index: /trunk/LMDZ.MARS/util/solzenangle.def
===================================================================
--- /trunk/LMDZ.MARS/util/solzenangle.def	(revision 2289)
+++ /trunk/LMDZ.MARS/util/solzenangle.def	(revision 2289)
@@ -0,0 +1,27 @@
+diagfi.nc
+tsurf
+ps
+temp
+
+morning
+90
+0
+
+-----------------------------------------------------------------------
+ABOVE is the list of inputs to be fed to "solzenangle.e" if you don t want
+to reply directly to the program:
+
+1) complete name of file to be read
+2) List of X variables to be kept (X lines) or 'all'
+3) blank line at the end
+4) planet side wrt to noon (morning or evening)
+5) solar zenith angle (within [0;180[ deg)
+0deg = zenith ; 90deg = terminator at the surface ; >90deg = terminator in altitude
+6) starttimeoffset :
+    - for a diagfi/concat : sol value at the beginning of the run, wrt Ls=0
+    - for a stats : sol value at the middle of the run, wrt Ls=0
+
+USE :
+solzenangle.e < solzenangle.def
+
+
Index: unk/LMDZ.MARS/util/terminator.F90
===================================================================
--- /trunk/LMDZ.MARS/util/terminator.F90	(revision 2288)
+++ 	(revision )
@@ -1,1519 +1,0 @@
-program terminator
-
-! ----------------------------------------------------------------------------
-! Program to redistribute and interpolate the variables at the local time
-! corresponding to one terminator (morning or evening) everywhere
-! input : diagfi.nc  / concat.nc / stats.nc kind of files
-! authors: A.Bierjon,F.Forget,2020
-! ----------------------------------------------------------------------------
-
-implicit none
-
-include "netcdf.inc" ! NetCDF definitions
-
-character (len=50)  file
-! file(): input file(s) names(s)
-character (len=30), dimension(15) :: notconcat
-! notconcat(): names of the (15) variables that won't be concatenated
-character (len=50), dimension(:), allocatable :: var
-! var(): name(s) of variable(s) that will be concatenated
-character (len=50) :: tmpvar,title,units
-! tmpvar(): used to temporarily store a variable name
-! title(): [netcdf] title attribute
-! units(): [netcdf] units attribute
-character (len=100) :: filename,vartmp
-! filename(): output file name
-! vartmp(): temporary variable name (read from netcdf input file)
-!character (len=1) :: ccopy
-! ccpy: 'y' or 'n' answer
-
-character (len=50) :: altlong_name,altunits,altpositive
-! altlong_name(): [netcdf] altitude "long_name" attribute
-! altunits(): [netcdf] altitude "units" attribute
-! altpositive(): [netcdf] altitude "positive" attribute
-
-character (len=8) :: sunterm
-! sunterm(): type of terminator ("sunrise" for the morning, "sunset" for the evening)
-
-integer :: nid,ierr,miss,validr
-! nid: [netcdf] file ID #
-! ierr: [netcdf] subroutine returned error code
-! miss: [netcdf] subroutine returned error code
-integer :: i,j,k,inter
-! for various loops
-integer :: varid
-! varid: [netcdf] variable ID #
-real, dimension(:), allocatable:: lat,lon,alt,ctl,time
-! lat(): array, stores latitude coordinates
-! lon(): array, stores longitude coordinates
-! alt(): array, stores altitude coordinates
-! ctl(): array, stores controle array
-! time(): array, stores time coordinates
-integer :: nbvar,nbvarfile,ndim
-! nbvar: # of variables to concatenate
-! nbfile: # number of input file(s)
-! nbvarfile: total # of variables in an input file
-! ndim: [netcdf] # (3 or 4) of dimensions (for variables)
-integer :: latdim,londim,altdim,ctldim,timedim
-! latdim: [netcdf] "latitude" dim ID
-! londim: [netcdf] "longitude" dim ID
-! altdim: [netcdf] "altdim" dim ID
-! ctldim: [netcdf] "controle" dim ID
-! timedim: [netcdf] "timedim" dim ID
-integer :: latvar,lonvar,altvar,ctlvar,timevar
-! latvar: [netcdf] ID of "latitude" variable
-! lonvar: [netcdf] ID of "longitude" variable
-! altvar: [netcdf] ID of "altitude" variable
-! ctlvar: [netcdf] ID of "controle" variable
-! timevar: [netcdf] ID of "Time" variable
-integer :: latlen,lonlen,altlen,ctllen,timelen,timelen_tot
-integer :: ilat,ilon,ialt,it
-! latlen: # of elements of lat() array
-! lonlen: # of elements of lon() array
-! altvar: # of elements of alt() array
-! ctlvar: # of elements of ctl() array
-! timelen: # of elements of time() array
-! timelen_tot:# =timelen or timelen+1 (if 1 more time to interpolate needed)
-integer :: nsol
-! nsol: # of sols in the input file AND # of elements of time() array in output
-integer :: nout,latdimout,londimout,altdimout,timedimout,timevarout
-! nout: [netcdf] output file ID
-! latdimout: [netcdf] output latitude (dimension) ID
-! londimout: [netcdf] output longitude (dimension) ID
-! altdimout: [netcdf] output altitude (dimension) ID
-! timedimout: [netcdf] output time (dimension) ID
-! timevarout: [netcdf] ID of output "Time" variable
-integer :: varidout 
-! varidout: [netcdf] variable ID # (of a variable to write to the output file)
-integer :: Nnotconcat,var_ok
-! Nnotconcat: # of (leading)variables that won't be concatenated
-! var_ok: flag (0 or 1)
-integer, dimension(4) :: corner,edges,dim
-! corner: [netcdf]
-! edges: [netcdf]
-! dim: [netcdf]
-real, dimension(:,:,:,:), allocatable :: var3d
-! var3D(,,,): 4D array to store a field
-
-real, dimension(:,:,:,:), allocatable :: var3d_lt
-! var3D_lt(,,,): 4D array to store a field in local time coordinate
-real, dimension(:), allocatable :: lt_gcm
-real, dimension(:), allocatable :: lt_outc
-real, dimension(:), allocatable :: var_gcm
-real, dimension(:), allocatable :: intsol
-real, dimension(:,:,:), allocatable :: lt_out
-real :: starttimeoffset ! offset (in sols) of sol 0 in input file, wrt Ls=0
-real :: tmpsol ! to temporarily store a sol value
-real :: tmpLs ! to temporarily store a Ls value
-real :: tmpLT ! to temporarily store a Local Time value
-
-real :: missing
-!PARAMETER(missing=1E+20)
-! missing: [netcdf] to handle "missing" values when reading/writing files
-real :: miss_lt
-PARAMETER(miss_lt=1E+20)
-! miss_lt: [netcdf] missing value to write for the terminator Local Time in the output file
-real, dimension(2) :: valid_range
-! valid_range(2): [netcdf] interval in which a value is considered valid
-logical :: stats   ! stats=T when reading a "stats" kind of ile
-
-
-!==============================================================================
-! 1.1. Get input file name(s)
-!==============================================================================
-write(*,*) 
-write(*,*) "which file do you want to use?  (diagfi... stats...  concat...)"
-
-read(*,'(a50)') file
-
-if (len_trim(file).eq.0) then
-   write(*,*) "no file... game over"
-   stop ""
-endif
-
-!==============================================================================
-! 1.2. Open the input file
-!==============================================================================
-
-ierr = NF_OPEN(file,NF_NOWRITE,nid)
-if (ierr.NE.NF_NOERR) then
-   write(*,*) 'ERROR: Pb opening file '//trim(file)
-   write(*,*) NF_STRERROR(ierr)
-   stop ""
-endif
-
-ierr=NF_INQ_NVARS(nid,nbvarfile)
-! nbvarfile now set to be the (total) number of variables in file
-if (ierr.NE.NF_NOERR) then
-   write(*,*) 'ERROR: Pb with NF_INQ_NVARS'
-   write(*,*) NF_STRERROR(ierr)
-   stop ""
-endif
-
-!==============================================================================
-! 1.3. List of variables that should not be processed
-!==============================================================================
-
-notconcat(1)='Time'
-notconcat(2)='controle'
-notconcat(3)='rlonu'
-notconcat(4)='latitude'
-notconcat(5)='longitude'
-notconcat(6)='altitude'
-notconcat(7)='rlatv'
-notconcat(8)='aps'
-notconcat(9)='bps'
-notconcat(10)='ap'
-notconcat(11)='bp'
-notconcat(12)='cu'
-notconcat(13)='cv'
-notconcat(14)='aire'
-notconcat(15)='phisinit'
-
-!==============================================================================
-! 1.4. Get (and check) list of variables to process
-!==============================================================================
-write(*,*)
-   Nnotconcat=0
-do i=1,nbvarfile
-   ierr=NF_INQ_VARNAME(nid,i,vartmp)
-   ! vartmp now contains the "name" of variable of ID # i
-   var_ok=0
-   do inter=1,15
-      if (vartmp.eq.notconcat(inter)) then
-         var_ok=1 
-         Nnotconcat=Nnotconcat+1
-      endif 
-   enddo        
-   if (var_ok.eq.0)  write(*,*) trim(vartmp)
-enddo
-
-! Nnotconcat: # of variables that won't be processed
-! nbvarfile: total # of variables in file
-allocate(var(nbvarfile-Nnotconcat),stat=ierr)
-if (ierr.ne.0) then
-  write(*,*) "Error: failed allocation of var(nbvarfile-Nnotconcat)"
-  write(*,*) "  nbvarfile=",nbvarfile
-  write(*,*) "  Nnotconcat=",Nnotconcat
-  stop
-endif
-
-
-write(*,*)
-write(*,*) "which variables do you want to redistribute ?"
-write(*,*) "all / list of <variables> (separated by <Enter>s)"
-write(*,*) "(an empty line , i.e: just <Enter>, implies end of list)"
-nbvar=0
-read(*,'(a50)') tmpvar
-do while ((tmpvar/=' ').AND.(trim(tmpvar)/='all'))
-   nbvar=nbvar+1
-   var(nbvar)=tmpvar
-   read(*,'(a50)') tmpvar
-enddo
-
-if (tmpvar=="all") then
-         nbvar=nbvarfile-Nnotconcat
-         do j=Nnotconcat+1,nbvarfile
-            ierr=nf_inq_varname(nid,j,var(j-Nnotconcat))
-         enddo
-! Variables names from the file are stored in var()
-   nbvar=nbvarfile-Nnotconcat
-   do i=1,nbvar
-      ierr=nf_inq_varname(nid,i+Nnotconcat,var(i))
-      write(*,'(a9,1x,i2,1x,a1,1x,a50)') "variable ",i,":",var(i)
-   enddo
-else if(nbvar==0) then
-   write(*,*) "no variable... game over"
-   stop ""
-endif ! of if (tmpvar=="all")
-
-
-!==============================================================================
-! 2.1. Open input file
-!==============================================================================
-
-      write(*,*) 
-      write(*,*) "opening "//trim(file)//"..."
-      ierr = NF_OPEN(file,NF_NOWRITE,nid)
-      if (ierr.NE.NF_NOERR) then
-         write(*,*) 'ERROR: Pb opening file '//trim(file)
-         write(*,*) NF_STRERROR(ierr)
-         stop ""
-      endif
-
-!==============================================================================
-! 2.2. Read (and check) dimensions of variables from input file
-!==============================================================================
-
-   ierr=NF_INQ_DIMID(nid,"latitude",latdim)
-   if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Dimension <latitude> is missing in file '//trim(file)
-      stop ""  
-   endif
-   ierr=NF_INQ_VARID(nid,"latitude",latvar)
-   if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Field <latitude> is missing in file '//trim(file)
-      stop ""  
-   endif
-   ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
-!  write(*,*) "latlen: ",latlen
-
-   ierr=NF_INQ_DIMID(nid,"longitude",londim)
-   if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Dimension <longitude> is missing in file '//trim(file)
-      stop ""  
-   endif
-   ierr=NF_INQ_VARID(nid,"longitude",lonvar)
-   if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Field <longitude> is missing in file'//trim(file)
-      stop "" 
-   endif
-   ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
-!  write(*,*) "lonlen: ",lonlen
-
-   ierr=NF_INQ_DIMID(nid,"altitude",altdim)
-   if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Dimension <altitude> is missing in file '//trim(file)
-      stop ""  
-   endif
-   ierr=NF_INQ_VARID(nid,"altitude",altvar)
-   if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Field <altitude> is missing in file'//trim(file)
-      stop ""
-   endif
-   ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
-!  write(*,*) "altlen: ",altlen
-
-   ierr=NF_INQ_DIMID(nid,"index",ctldim)
-   if (ierr.NE.NF_NOERR) then
-      write(*,*) ' Dimension <index> is missing in file '//trim(file)
-      ctllen=0
-      !stop ""  
-   endif
-   ierr=NF_INQ_VARID(nid,"controle",ctlvar)
-   if (ierr.NE.NF_NOERR) then
-      write(*,*) 'Field <controle> is missing in file '//trim(file)
-      ctllen=0
-      !stop ""
-   else
-      ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
-   endif
-!  write(*,*) "controle: ",controle
-
-!==============================================================================
-! 2.3. Read (and check compatibility of) dimensions of
-!       variables from input file
-!==============================================================================
-
-! First call; initialize/allocate
-      allocate(lat(latlen),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate lat(latlen)"
-        stop
-      endif
-      allocate(lon(lonlen),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate lon(lonlen)"
-        stop
-      endif
-      allocate(alt(altlen),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate alt(altlen)"
-        stop
-      endif
-      allocate(ctl(ctllen),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate ctl(ctllen)"
-        stop
-      endif
-
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat)
-#else
-      ierr = NF_GET_VAR_REAL(nid,latvar,lat)
-#endif
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to load latitude"
-        write(*,*) NF_STRERROR(ierr)
-        stop
-      endif
-
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon)
-#else
-      ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
-#endif
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to load longitude"
-        write(*,*) NF_STRERROR(ierr)
-        stop
-      endif
-
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt)
-#else
-      ierr = NF_GET_VAR_REAL(nid,altvar,alt)
-#endif
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to load altitude"
-        write(*,*) NF_STRERROR(ierr)
-        stop
-      endif
-! Get altitude attributes to handle files with any altitude type
-ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
-ierr=nf_get_att_text(nid,altvar,'units',altunits)
-ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
-
-      if (ctllen .ne. 0) then
-#ifdef NC_DOUBLE
-        ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
-#else
-        ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
-#endif
-        if (ierr.ne.0) then
-          write(*,*) "Error: failed to load controle"
-          write(*,*) NF_STRERROR(ierr)
-          stop
-        endif
-      endif ! of if (ctllen .ne. 0)
-
-!==============================================================================
-! 2.4. Handle "Time" dimension from input file
-!==============================================================================
-
-!==============================================================================
-! 2.4.0 Read "Time" dimension from input file
-!==============================================================================
-   ierr=NF_INQ_DIMID(nid,"Time",timedim)
-   if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Dimension <Time> is missing in file'//trim(file)
-      stop ""
-   endif
-   ierr=NF_INQ_VARID(nid,"Time",timevar)
-   if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Field <Time> is missing in file'//trim(file)
-      stop ""
-   endif
-   ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
-!  write(*,*) "timelen: ",timelen
-
-   ! Sanity Check: Does "Time" has a "title" attribute
-   ! and is it "Solar longitude" ?
-   ierr=nf_get_att_text(nid,timevar,"title",title)
-   if ((ierr.EQ.NF_NOERR).and.(title.eq."Solar longitude")) then
-     write(*,*) "ERROR: Time axis in input file is in Solar Longitude!"
-     write(*,*) "       localtime requires sols as time axis!"
-     write(*,*) " Might as well stop here."
-     stop
-   endif
-
-   ! allocate time() array and fill it with values from input file
-   allocate(time(timelen),stat=ierr)
-   if (ierr.ne.0) then
-     write(*,*) "Error: failed to allocate time(timelen)"
-     stop
-   endif
-
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
-#else
-   ierr = NF_GET_VAR_REAL(nid,timevar,time)
-#endif
-   if (ierr.NE.NF_NOERR) then
-     write(*,*) "Error , failed to load Time"
-     write(*,*) NF_STRERROR(ierr)
-     stop
-   endif
-
-
-!==============================================================================
-! 2.4.1 Select local times to be saved "Time" in output file
-!==============================================================================
-      write(*,*) 'Type of terminator?'
-      read(*,*) sunterm
-      if ((trim(sunterm).ne."sunrise").and.(trim(sunterm).ne."sunset")) then
-        write(*,*) "sunterm = "//sunterm
-        write(*,*) "Error: terminator type must be sunrise or sunset"
-        stop
-      endif
-      write(*,*) "Beginning date of the simulation file "
-      write(*,*) "(or median date of the simulation for a stats)?"
-      write(*,*) "(i.e. number of sols since Ls=0 at the Time=0.0 in the input file)"
-      read(*,*) starttimeoffset
-
-      if ((timelen.eq.12).and.(time(1).eq.2).and.(time(timelen).eq.24))then
-         write(*,*) 'We detect a ""stats"" file'
-         stats=.true.
-!        We rewrite the time from "stats" from 0 to 1 sol...
-         do it=1,timelen  ! loop temps in
-               time(it) = time(it)/24.
-         end do
-         nsol =1
-      else   ! case of a diagfi or concatnc file 
-        stats=.false.
-!       Number of sol in input file
-        nsol = int(time(timelen)) - int(time(1))
-      end if
-      
-      allocate(intsol(nsol),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate intsol(nsol)"
-        stop
-      endif
-      allocate(lt_out(lonlen,latlen,nsol),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate lt_out(lonlen,latlen,nsol)"
-        stop
-      endif
-      do it=1,nsol
-        intsol(it)= int(time(1)) + it-1.
-        do ilon=1,lonlen
-!         For the computation of Ls, it's supposed that the morning/evening
-!         terminator happens at LT=6h/18h at the given longitude, which is
-!         then reported to a sol value at longitude 0
-          if (trim(sunterm).eq."sunrise") then
-            tmpsol = intsol(it) + (6.-lon(ilon)/15.)/24. + starttimeoffset
-          else !if trim(sunterm).eq."sunset"
-            tmpsol = intsol(it) + (18.-lon(ilon)/15.)/24. + starttimeoffset            
-          endif
-          call sol2ls(tmpsol,tmpLs)
-          do ilat=1,latlen
-!           Compute the Local Time of the terminator at a given longitude
-            call LTterminator(lat(ilat),tmpLs,sunterm,miss_lt,tmpLT)
-!           Deduce the sol value at this longitude
-            if (tmpLT.eq.miss_lt) then 
-              ! if there is no terminator, put a missing value in lt_out
-              lt_out(ilon,ilat,it)=miss_lt
-            else
-              lt_out(ilon,ilat,it)=tmpLT
-            endif
-          enddo
-        enddo
-      enddo
-
-      if (nsol.gt.1) then
-         timelen_tot=timelen
-      else
-!        if only 1 sol available, we must add 1 timestep for the interpolation
-         timelen_tot=timelen+1
-      endif      
-      allocate(lt_gcm(timelen_tot),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate lt_gcm(timelen_tot)"
-        stop
-      endif
-      allocate(var_gcm(timelen_tot),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate var_gcm(timelen_tot)"
-        stop
-      endif
-      allocate(lt_outc(nsol),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate lt_outc(nsol)"
-        stop
-      endif
-
-!==============================================================================
-! 2.4.2. Get output file name
-!==============================================================================
-    if (trim(sunterm).eq."sunrise") then 
-      filename=file(1:len_trim(file)-3)//"_TMO.nc"
-    else ! if trim(sunterm).eq."sunset"
-      filename=file(1:len_trim(file)-3)//"_TEV.nc"
-    endif
-
-!==============================================================================
-! 2.4.3. Initiate dimension in output file
-!==============================================================================
-
-
-   ! Initialize output file's lat,lon,alt and time dimensions
-    call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
-          altlong_name,altunits,altpositive,&
-          nout,latdimout,londimout,altdimout,timedimout,timevarout)
-
-   ! Initialize output file's aps,bps and phisinit variables
-    call init2(nid,lonlen,latlen,altlen,altdim,&
-               nout,londimout,latdimout,altdimout)
-
-
-!==============================================================================
-! 2.4.4. Write/extend "Time" dimension/values in output file
-!==============================================================================
-
-   if (stats) then 
-
-     do it=1,nsol
-#ifdef NC_DOUBLE
-        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,intsol(it)*24.)
-#else
-        ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,intsol(it)*24.)
-#endif
-     enddo
-   else
-     do it=1,nsol
-#ifdef NC_DOUBLE
-        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,intsol(it))
-#else
-        ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,intsol(it))
-#endif
-        if (ierr.NE.NF_NOERR) then
-          write(*,*) "Error , failed to write Time"
-          stop
-        endif
-     enddo
-   end if
-
-!==============================================================================
-! 2.5 Read/write variables
-!==============================================================================
-
-   do j=1,nbvar ! loop on variables to read/write
-
-!==============================================================================
-! 2.5.1 Check that variable to be read existe in input file
-!==============================================================================
-
-      write(*,*) "variable ",var(j)
-      ierr=nf_inq_varid(nid,var(j),varid)
-      if (ierr.NE.NF_NOERR) then
-         write(*,*) 'ERROR: Field <',var(j),'> not found in file'//file
-         stop ""
-      endif
-      ierr=nf_inq_varndims(nid,varid,ndim)
-
-!==============================================================================
-! 2.5.2 Prepare things in order to read/write the variable
-!==============================================================================
-
-      ! build dim(),corner() and edges() arrays
-      ! and allocate var3d() array
-      if (ndim==3) then
-         allocate(var3d(lonlen,latlen,timelen,1))
-         allocate(var3d_lt(lonlen,latlen,nsol,1))
-         dim(1)=londimout
-         dim(2)=latdimout
-         dim(3)=timedimout
-
-         ! start indexes (where data values will be written)
-         corner(1)=1
-         corner(2)=1
-         corner(3)=1
-         corner(4)=1
-
-         ! length (along dimensions) of block of data to be written
-         edges(1)=lonlen 
-         edges(2)=latlen 
-         edges(3)=nsol
-         edges(4)=1
-   
-      else if (ndim==4) then
-         allocate(var3d(lonlen,latlen,altlen,timelen))
-         allocate(var3d_lt(lonlen,latlen,altlen,nsol))
-         dim(1)=londimout
-         dim(2)=latdimout
-         dim(3)=altdimout
-         dim(4)=timedimout
-
-         ! start indexes (where data values will be written)
-         corner(1)=1
-         corner(2)=1
-         corner(3)=1
-         corner(4)=1
-
-         ! length (along dimensions) of block of data to be written
-         edges(1)=lonlen
-         edges(2)=latlen
-         edges(3)=altlen
-         edges(4)=nsol
-      endif
-
-         units=" "
-         title=" "
-         ierr=nf_get_att_text(nid,varid,"title",title)
-         ierr=nf_get_att_text(nid,varid,"units",units)
-         call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr)
-
-      
-
-!==============================================================================
-! 2.5.3 Read from input file 
-!==============================================================================
-
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid,varid,var3d)
-      miss=NF_GET_ATT_DOUBLE(nid,varid,"missing_value",missing)
-      validr=NF_GET_ATT_DOUBLE(nid,varid,"valid_range",valid_range)
-#else
-      ierr = NF_GET_VAR_REAL(nid,varid,var3d)
-      miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
-      validr=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
-#endif
-
-!==============================================================================
-! 2.5.4 Interpolation in local time 
-!==============================================================================
-
-        do ilon=1,lonlen
-!          write(*,*) 'processing longitude =', lon(ilon)
-!          Local time at each longitude :
-           do it=1,timelen  ! loop temps in
-               lt_gcm(it) = time(it) + lon(ilon)/360.
-           enddo
-           if (nsol.eq.1) lt_gcm(timelen+1) = lt_gcm(1)+1
-
-           do ilat=1,latlen
-             do it=1,nsol ! loop time local time
-               if (lt_out(ilon,ilat,it).eq.miss_lt) then
-                 lt_outc(it)=miss_lt
-               else
-                 lt_outc(it)=intsol(it) + lt_out(ilon,ilat,it)/24.
-               endif
-!              If the computed local time exceeds the GCM input file time bounds,
-!              put a missing value in the local time variable
-               if ((lt_outc(it).gt.lt_gcm(timelen_tot)).OR.(lt_outc(it).lt.lt_gcm(1))) then
-                 lt_out(ilon,ilat,it)=miss_lt
-                 lt_outc(it)=miss_lt
-               endif
-             enddo ! local time
-
-             if (ndim==3) then
-               do it=1,timelen  ! loop temps in
-                  var_gcm(it) = var3d(ilon,ilat,it,1)
-               enddo
-               if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
-               do it=1,nsol ! loop time local time
-                 if (lt_outc(it).eq.miss_lt) then
-                   var3d_lt(ilon,ilat,it,1)=missing
-                 else
-                   call interpolf(lt_outc(it),var3d_lt(ilon,ilat,it,1),&
-                                 missing,lt_gcm,var_gcm,timelen_tot)
-                 endif
-               enddo ! local time
-
-             else if  (ndim==4) then
-               do ialt=1,altlen
-                 do it=1,timelen ! loop temps in
-                    var_gcm(it) = var3d(ilon,ilat,ialt,it)
-                 enddo
-                 if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
-                 do it=1,nsol ! loop time local time
-                   if (lt_outc(it).eq.miss_lt) then
-                     var3d_lt(ilon,ilat,ialt,it)=missing
-                   else
-                     call interpolf(lt_outc(it),var3d_lt(ilon,ilat,ialt,it),&
-                                   missing,lt_gcm,var_gcm,timelen_tot)
-                   endif
-                 enddo ! local time   
-               enddo ! alt
-             endif
-             
-           enddo ! lat
-        end do ! lon
-
-
-!==============================================================================
-! 2.5.5 write (append) to the output file
-!==============================================================================
-
-#ifdef NC_DOUBLE
-      ierr= NF_PUT_VARA_DOUBLE(nout,varidout,corner,edges,var3d_lt)
-#else
-      ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3d_lt)
-#endif
-
-      if (ierr.ne.NF_NOERR) then
-         write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
-         stop ""
-      endif
-
-! In case there is a "valid_range" attribute
-      ! Write "valid_range" and "missing_value" attributes in output file
-      if (miss.eq.NF_NOERR) then
-         call missing_value(nout,varidout,validr,miss,valid_range,missing)
-      endif
-
-      ! free var3d() array
-      deallocate(var3d)
-      deallocate(var3d_lt)
-
-   enddo ! of do j=1,nbvar
-
-!==============================================================================
-! 2.5.6 Write the terminator Local Time variable in output file
-!==============================================================================
-   write(*,*) "variable LTterminator"
-   ! build dim(),corner() and edges() arrays
-   dim(1)=londimout
-   dim(2)=latdimout
-   dim(3)=timedimout
-
-   ! start indexes (where data values will be written)
-   corner(1)=1
-   corner(2)=1
-   corner(3)=1
-   corner(4)=1
-
-   ! length (along dimensions) of block of data to be written
-   edges(1)=lonlen 
-   edges(2)=latlen 
-   edges(3)=nsol
-   edges(4)=1
-
-   tmpvar="LTterminator"
-   units="hour(LTST)"
-   if (trim(sunterm).eq."sunrise") then
-     title="Morning Terminator Local Time"
-   else !if trim(sunterm).eq."sunset"
-     title="Evening Terminator Local Time"            
-   endif
-   call def_var(nout,tmpvar,title,units,3,dim,varidout,ierr)
-
-#ifdef NC_DOUBLE
-   ierr= NF_PUT_VARA_DOUBLE(nout,varidout,corner,edges,lt_out)
-#else
-   ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,lt_out)
-#endif
-   if (ierr.ne.NF_NOERR) then
-      write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
-      stop ""
-   endif
-
-   validr=NF_NOERR+1
-   miss=NF_NOERR
-   call missing_value(nout,varidout,validr,miss,(0,0),miss_lt)
-   write(*,*)"with missing value = ",miss_lt
-!==============================================================================
-! 3. END
-!==============================================================================
-   deallocate(time)
-   deallocate(intsol)
-   deallocate(lt_gcm)
-   deallocate(lt_out)
-   deallocate(lt_outc)
-   deallocate(var_gcm)
-
-   ! Close input file
-   ierr=nf_close(nid)
-
-! Close output file
-ierr=NF_CLOSE(nout)
-
-write(*,*) "Program completed !"
-
-end program terminator
-
-
-
-!******************************************************************************
-!******************************************************************************
-
-subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
-                     altlong_name,altunits,altpositive,&
-                     nout,latdimout,londimout,altdimout,timedimout,timevarout)
-!==============================================================================
-! Purpose:
-! Create and initialize a data file (NetCDF format)
-!==============================================================================
-! Remarks:
-! The NetCDF file (created in this subroutine) remains open
-!==============================================================================
-
-implicit none
-
-include "netcdf.inc" ! NetCDF definitions
-
-!==============================================================================
-! Arguments:
-!==============================================================================
-character (len=*), intent(in):: filename
-! filename(): the file's name
-integer,intent(in) ::latlen,lonlen,altlen,ctllen
-real, intent(in):: lat(latlen)
-! lat(): latitude
-real, intent(in):: lon(lonlen)
-! lon(): longitude
-real, intent(in):: alt(altlen)
-! alt(): altitude
-real, intent(in):: ctl(ctllen)
-! ctl(): controle
-character (len=*), intent(in) :: altlong_name
-! altlong_name(): [netcdf] altitude "long_name" attribute
-character (len=*), intent(in) :: altunits
-! altunits(): [netcdf] altitude "units" attribute
-character (len=*), intent(in) :: altpositive
-! altpositive(): [netcdf] altitude "positive" attribute
-integer, intent(out):: nout
-! nout: [netcdf] file ID
-integer, intent(out):: latdimout
-! latdimout: [netcdf] lat() (i.e.: latitude)  ID
-integer, intent(out):: londimout
-! londimout: [netcdf] lon()  ID
-integer, intent(out):: altdimout
-! altdimout: [netcdf] alt()  ID
-integer, intent(out):: timedimout
-! timedimout: [netcdf] "Time"  ID
-integer, intent(out):: timevarout
-! timevarout: [netcdf] "Time" (considered as a variable) ID
-
-!==============================================================================
-! Local variables:
-!==============================================================================
-!integer :: latdim,londim,altdim,timedim
-integer :: nvarid,ierr
-integer :: ctldimout
-! nvarid: [netcdf] ID of a variable
-! ierr: [netcdf]  return error code (from called subroutines)
-
-!==============================================================================
-! 1. Create (and open) output file
-!==============================================================================
-write(*,*) "creating "//trim(adjustl(filename))//'...'
-ierr = NF_CREATE(filename,IOR(NF_CLOBBER,NF_64BIT_OFFSET),nout)
-! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
-if (ierr.NE.NF_NOERR) then
-   WRITE(*,*)'ERROR: Impossible to create the file ',trim(filename)
-   write(*,*) NF_STRERROR(ierr)
-   stop ""
-endif
-
-!==============================================================================
-! 2. Define/write "dimensions" and get their IDs
-!==============================================================================
-
-ierr = NF_DEF_DIM(nout, "latitude", latlen, latdimout)
-if (ierr.NE.NF_NOERR) then
-   WRITE(*,*)'initiate: error failed to define dimension <latitude>'
-   write(*,*) NF_STRERROR(ierr)
-   stop ""
-endif
-ierr = NF_DEF_DIM(nout, "longitude", lonlen, londimout)
-if (ierr.NE.NF_NOERR) then
-   WRITE(*,*)'initiate: error failed to define dimension <longitude>'
-   write(*,*) NF_STRERROR(ierr)
-   stop ""
-endif
-ierr = NF_DEF_DIM(nout, "altitude", altlen, altdimout)
-if (ierr.NE.NF_NOERR) then
-   WRITE(*,*)'initiate: error failed to define dimension <altitude>'
-   write(*,*) NF_STRERROR(ierr)
-   stop ""
-endif
-if (size(ctl).ne.0) then
-  ierr = NF_DEF_DIM(nout, "index", ctllen, ctldimout)
-  if (ierr.NE.NF_NOERR) then
-    WRITE(*,*)'initiate: error failed to define dimension <index>'
-    write(*,*) NF_STRERROR(ierr)
-    stop ""
-  endif
-endif
-ierr = NF_DEF_DIM(nout, "Time", NF_UNLIMITED, timedimout)
-if (ierr.NE.NF_NOERR) then
-   WRITE(*,*)'initiate: error failed to define dimension <Time>'
-   write(*,*) NF_STRERROR(ierr)
-   stop ""
-endif
-
-! End netcdf define mode
-ierr = NF_ENDDEF(nout)
-if (ierr.NE.NF_NOERR) then
-   WRITE(*,*)'initiate: error failed to switch out of define mode'
-   write(*,*) NF_STRERROR(ierr)
-   stop ""
-endif
-
-!==============================================================================
-! 3. Write "Time" and its attributes
-!==============================================================================
-
-call def_var(nout,"Time","Time","years since 0000-00-0 00:00:00",1,&
-             (/timedimout/),timevarout,ierr)
-
-!==============================================================================
-! 4. Write "latitude" (data and attributes)
-!==============================================================================
-
-call def_var(nout,"latitude","latitude","degrees_north",1,&
-             (/latdimout/),nvarid,ierr)
-
-#ifdef NC_DOUBLE
-ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lat)
-#else
-ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
-#endif
-if (ierr.NE.NF_NOERR) then
-   WRITE(*,*)'initiate: error failed writing variable <latitude>'
-   write(*,*) NF_STRERROR(ierr)
-   stop ""
-endif
-
-!==============================================================================
-! 4. Write "longitude" (data and attributes)
-!==============================================================================
-
-call def_var(nout,"longitude","East longitude","degrees_east",1,&
-             (/londimout/),nvarid,ierr)
-
-#ifdef NC_DOUBLE
-ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lon)
-#else
-ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
-#endif
-if (ierr.NE.NF_NOERR) then
-   WRITE(*,*)'initiate: error failed writing variable <longitude>'
-   write(*,*) NF_STRERROR(ierr)
-   stop ""
-endif
-
-!==============================================================================
-! 5. Write "altitude" (data and attributes)
-!==============================================================================
-
-! Switch to netcdf define mode
-ierr = NF_REDEF (nout)
-
-#ifdef NC_DOUBLE
-ierr = NF_DEF_VAR (nout,"altitude",NF_DOUBLE,1,altdimout,nvarid)
-#else
-ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,altdimout,nvarid)
-#endif
-
-ierr = NF_PUT_ATT_TEXT (nout,nvarid,'long_name',len_trim(adjustl(altlong_name)),adjustl(altlong_name))
-ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',len_trim(adjustl(altunits)),adjustl(altunits))
-ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',len_trim(adjustl(altpositive)),adjustl(altpositive))
-
-! End netcdf define mode
-ierr = NF_ENDDEF(nout)
-
-#ifdef NC_DOUBLE
-ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,alt)
-#else
-ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
-#endif 
-if (ierr.NE.NF_NOERR) then
-   WRITE(*,*)'initiate: error failed writing variable <altitude>'
-   write(*,*) NF_STRERROR(ierr)
-   stop ""
-endif
-
-!==============================================================================
-! 6. Write "controle" (data and attributes)
-!==============================================================================
-
-if (size(ctl).ne.0) then
-   ! Switch to netcdf define mode
-   ierr = NF_REDEF (nout)
-
-#ifdef NC_DOUBLE
-   ierr = NF_DEF_VAR (nout,"controle",NF_DOUBLE,1,ctldimout,nvarid)
-#else
-   ierr = NF_DEF_VAR (nout,"controle",NF_FLOAT,1,ctldimout,nvarid)
-#endif
-
-   ierr = NF_PUT_ATT_TEXT (nout,nvarid,"title",18,"Control parameters")
-
-   ! End netcdf define mode
-   ierr = NF_ENDDEF(nout)
-
-#ifdef NC_DOUBLE
-   ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,ctl)
-#else
-   ierr = NF_PUT_VAR_REAL (nout,nvarid,ctl)
-#endif
-   if (ierr.NE.NF_NOERR) then
-      WRITE(*,*)'initiate: error failed writing variable <controle>'
-      write(*,*) NF_STRERROR(ierr)
-      stop ""
-   endif
-endif
-
-end Subroutine initiate
-
-!******************************************************************************
-subroutine init2(infid,lonlen,latlen,altlen,altdimid, &
-                 outfid,londimout,latdimout,altdimout)
-!==============================================================================
-! Purpose:
-! Copy aps(), bps() and phisinit() from input file to output file
-!==============================================================================
-! Remarks:
-! The NetCDF files must be open
-!==============================================================================
-
-implicit none
-
-include "netcdf.inc" ! NetCDF definitions
-
-!==============================================================================
-! Arguments:
-!==============================================================================
-integer, intent(in) :: infid  ! NetCDF output file ID
-integer, intent(in) :: lonlen ! # of grid points along longitude
-integer, intent(in) :: latlen ! # of grid points along latitude
-integer, intent(in) :: altlen ! # of grid points along altitude
-integer, intent(in) :: altdimid ! ID of altitude dimension
-integer, intent(in) :: outfid ! NetCDF output file ID
-integer, intent(in) :: londimout ! longitude dimension ID
-integer, intent(in) :: latdimout ! latitude dimension ID
-integer, intent(in) :: altdimout ! altitude dimension ID
-!==============================================================================
-! Local variables:
-!==============================================================================
-real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
-real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
-integer :: apsid,bpsid,phisinitid
-integer :: ierr
-integer :: tmpdimid ! temporary dimension ID
-integer :: tmpvarid ! temporary variable ID
-logical :: phis, aps_ok, bps_ok ! are "phisinit" "aps" "bps" available ?
-
-
-!==============================================================================
-! 1. Read data from input file
-!==============================================================================
-
-! hybrid coordinate aps
-  allocate(aps(altlen),stat=ierr)
-  if (ierr.ne.0) then
-    write(*,*) "init2: failed to allocate aps(altlen)"
-    stop
-  endif
-
-ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
-if (ierr.ne.NF_NOERR) then
-  write(*,*) "oops Failed to get aps ID. OK"
-  aps_ok=.false.
-else
-  ! Check that aps() is the right size (it most likely won't be if
-  ! zrecast has be used to generate the input file)
-  ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
-  if (tmpdimid.ne.altdimid) then
-    write(*,*) "init2: wrong dimension size for aps(), skipping it ..."
-    aps_ok=.false.
-  else
-    ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
-    if (ierr.ne.NF_NOERR) then
-     stop "init2 error: Failed reading aps"
-    endif
-    aps_ok=.true.
-  endif
-endif
-
-! hybrid coordinate bps
-  allocate(bps(altlen),stat=ierr)
-  if (ierr.ne.0) then
-    write(*,*) "init2: failed to allocate bps(altlen)"
-    stop
-  endif
-
-ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
-if (ierr.ne.NF_NOERR) then
-  write(*,*) "oops: Failed to get bps ID. OK"
-  bps_ok=.false.
-else
-  ! Check that bps() is the right size
-  ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
-  if (tmpdimid.ne.altdimid) then
-    write(*,*) "init2: wrong dimension size for bps(), skipping it ..."
-    bps_ok=.false.
-  else
-    ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
-    if (ierr.ne.NF_NOERR) then
-      stop "init2 Error: Failed reading bps"
-    endif
-    bps_ok=.true.
-  endif
-endif
-
-! ground geopotential phisinit
-allocate(phisinit(lonlen,latlen),stat=ierr)
-if (ierr.ne.0) then
-  write(*,*) "init2: failed to allocate phisinit(lonlen,latlen)"
-  stop
-endif
-ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
-if (ierr.ne.NF_NOERR) then
-  write(*,*) "Failed to get phisinit ID. OK"
-  phisinit = 0.
-  phis = .false.
-else
-  ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
-  if (ierr.ne.NF_NOERR) then
-    stop "init2 Error: Failed reading phisinit"
-  endif
-  phis = .true.
-endif
-
-
-!==============================================================================
-! 2. Write
-!==============================================================================
-
-!==============================================================================
-! 2.2. Hybrid coordinates aps() and bps()
-!==============================================================================
-
-IF (aps_ok) then 
-  ! define aps
-  call def_var(outfid,"aps","hybrid pressure at midlayers"," ",1,&
-               (/altdimout/),apsid,ierr)
-  if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed to def_var aps"
-  endif
-
-  ! write aps
-#ifdef NC_DOUBLE
-  ierr=NF_PUT_VAR_DOUBLE(outfid,apsid,aps)
-#else
-  ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
-#endif
-  if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed to write aps"
-  endif
-ENDIF ! of IF (aps_ok)
-
-IF (bps_ok) then 
-  ! define bps
-  call def_var(outfid,"bps","hybrid sigma at midlayers"," ",1,&
-               (/altdimout/),bpsid,ierr)
-  if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed to def_var bps"
-  endif
-
-  ! write bps
-#ifdef NC_DOUBLE
-  ierr=NF_PUT_VAR_DOUBLE(outfid,bpsid,bps)
-#else
-  ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
-#endif
-  if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed to write bps"
-  endif
-ENDIF ! of IF (bps_ok)
-
-!==============================================================================
-! 2.2. phisinit()
-!==============================================================================
-
-IF (phis) THEN
-
-  !define phisinit
-  call def_var(outfid,"phisinit","Ground level geopotential"," ",2,&
-              (/londimout,latdimout/),phisinitid,ierr)
-  if (ierr.ne.NF_NOERR) then
-     stop "Error: Failed to def_var phisinit"
-  endif
-
-  ! write phisinit
-#ifdef NC_DOUBLE
-  ierr=NF_PUT_VAR_DOUBLE(outfid,phisinitid,phisinit)
-#else
-  ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
-#endif
-  if (ierr.ne.NF_NOERR) then
-    stop "Error: Failed to write phisinit"
-  endif
-
-ENDIF ! of IF (phis)
-
-
-! Cleanup
-deallocate(aps)
-deallocate(bps)
-deallocate(phisinit)
-
-end subroutine init2
-
-!******************************************************************************
-subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
-!==============================================================================
-! Purpose: Write a variable (i.e: add a variable to a dataset)
-! called "name"; along with its attributes "title", "units"...
-! to a file (following the NetCDF format)
-!==============================================================================
-! Remarks:
-! The NetCDF file must be open
-!==============================================================================
-
-implicit none
-
-include "netcdf.inc" ! NetCDF definitions
-
-!==============================================================================
-! Arguments:
-!==============================================================================
-integer, intent(in) :: nid
-! nid: [netcdf] file ID #
-character (len=*), intent(in) :: name
-! name(): [netcdf] variable's name
-character (len=*), intent(in) :: title
-! title(): [netcdf] variable's "title" attribute
-character (len=*), intent(in) :: units
-! unit(): [netcdf] variable's "units" attribute
-integer, intent(in) :: nbdim
-! nbdim: number of dimensions of the variable
-integer, dimension(nbdim), intent(in) :: dim
-! dim(nbdim): [netcdf] dimension(s) ID(s)
-integer, intent(out) :: nvarid
-! nvarid: [netcdf] ID # of the variable
-integer, intent(out) :: ierr
-! ierr: [netcdf] subroutines returned error code
-
-! Switch to netcdf define mode
-ierr=NF_REDEF(nid)
-
-! Insert the definition of the variable
-#ifdef NC_DOUBLE
-ierr=NF_DEF_VAR(nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
-#else
-ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
-#endif
-
-! Write the attributes
-ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",len_trim(adjustl(title)),adjustl(title))
-ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
-
-! End netcdf define mode
-ierr=NF_ENDDEF(nid)
-
-end subroutine def_var
-
-!******************************************************************************
-subroutine  missing_value(nout,nvarid,validr,miss,valid_range,missing)
-!==============================================================================
-! Purpose:
-! Write "valid_range" and "missing_value" attributes (of a given
-! variable) to a netcdf file
-!==============================================================================
-! Remarks:
-! NetCDF file must be open
-! Variable (nvarid) ID must be set
-!==============================================================================
-
-implicit none
-
-include "netcdf.inc"  ! NetCDF definitions
-
-!==============================================================================
-! Arguments:
-!==============================================================================
-integer, intent(in) :: nout
-! nout: [netcdf] file ID #
-integer, intent(in) :: nvarid
-! varid: [netcdf] variable ID #
-integer, intent(in) :: validr
-! validr : [netcdf] routines return code for "valid_range" attribute
-integer, intent(in) :: miss
-! miss : [netcdf] routines return code for "missing_value" attribute
-real, dimension(2), intent(in) :: valid_range
-! valid_range(2): [netcdf] "valid_range" attribute (min and max)
-real, intent(in) :: missing
-! missing: [netcdf] "missing_value" attribute
-
-!==============================================================================
-! Local variables:
-!==============================================================================
-integer :: ierr
-! ierr: [netcdf] subroutine returned error code
-!      INTEGER nout,nvarid,ierr
-!      REAL missing
-!      REAL valid_range(2)
-
-! Switch to netcdf dataset define mode
-ierr = NF_REDEF (nout)
-if (ierr.ne.NF_NOERR) then
-   print*,'missing_value: '
-   print*, NF_STRERROR(ierr)
-endif
-
-! Write valid_range() attribute
-if (validr.eq.NF_NOERR) then
-#ifdef NC_DOUBLE
-  ierr=NF_PUT_ATT_DOUBLE(nout,nvarid,'valid_range',NF_DOUBLE,2,valid_range)
-#else
-  ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
-#endif
-
-  if (ierr.ne.NF_NOERR) then
-     print*,'missing_value: valid_range attribution failed'
-     print*, NF_STRERROR(ierr)
-     !write(*,*) 'NF_NOERR', NF_NOERR
-     !CALL abort
-     stop ""
-  endif
-endif
-
-! Write "missing_value" attribute
-if (miss.eq.NF_NOERR) then
-#ifdef NC_DOUBLE
-  ierr= NF_PUT_ATT_DOUBLE(nout,nvarid,'missing_value',NF_DOUBLE,1,missing)
-#else
-  ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
-#endif
-
-  if (ierr.NE.NF_NOERR) then
-     print*, 'missing_value: missing value attribution failed'
-     print*, NF_STRERROR(ierr)
-  !    WRITE(*,*) 'NF_NOERR', NF_NOERR
-  !    CALL abort
-     stop ""
-  endif
-endif
-
-! End netcdf dataset define mode
-ierr = NF_ENDDEF(nout)
-
-end subroutine  missing_value
-
-!*****************************************************************************
-subroutine interpolf(x,y,missing,xd,yd,nd)
-!==============================================================================
-! Purpose:
-! Yield y=f(x), where xd() end yd() are arrays of known values,
-! using linear interpolation
-! If x is not included in the interval spaned by xd(), then y is set
-! to a default value 'missing'
-! Note:
-! Array xd() should contain ordered (either increasing or decreasing) abscissas
-!==============================================================================
-implicit none
-!==============================================================================
-! Arguments:
-!==============================================================================
-real,intent(in) :: x ! abscissa at which interpolation is to be done
-real,intent(in) :: missing ! missing value (if no interpolation is performed)
-integer,intent(in) :: nd ! size of arrays
-real,dimension(nd),intent(in) :: xd ! array of known absissas
-real,dimension(nd),intent(in) :: yd ! array of correponding values
-
-real,intent(out) :: y ! interpolated value
-!==============================================================================
-! Local variables:
-!==============================================================================
-integer :: i
-
-! default: set y to 'missing'
-y=missing
-
-   do i=1,nd-1
-     if (((x.ge.xd(i)).and.(x.le.xd(i+1))).or.&
-          ((x.le.xd(i)).and.(x.ge.xd(i+1)))) then
-        if ((yd(i).eq.missing).or.(yd(i+1).eq.missing)) then
-          ! cannot perform the interpolation if an encompassing value
-          ! is already set to 'missing'
-        else
-          !linear interpolation based on encompassing values
-          y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i))
-        endif
-        exit
-     endif
-   enddo
-
-
-end subroutine interpolf
-
-!******************************************************************************
-subroutine LTterminator(lat,ls,sunterm,missval,lt_term)
-! compute the terminator localtime by inversing
-! the solar zenith angle equation :
-! mu0=sin(lat)*sin(declin)+
-!     cos(lat)*cos(declin)*cos(2.*pi*(localtime/24.-.5))
-! for mu0=cos(sza)=0 (since sza=90deg at the terminator)
-
-implicit none
-
-real,intent(in) :: lat ! latitude (deg)
-real,intent(in) :: ls ! solar longitude (deg)
-character (len=8), intent(in) :: sunterm ! terminator type : "sunrise" for morning,"sunset" for evening
-real,intent(in) :: missval ! missing value for when there is no terminator (polar day/night)
-real,intent(out) :: lt_term ! local true solar time (hours) of the terminator
-
-! local variables:
-double precision :: declin ! declination of the Sun (rad)
-double precision,parameter :: obliquity=25.1919d0
-double precision,parameter :: pi=3.14159265358979d0
-double precision,parameter :: degtorad=pi/180.0d0
-double precision,parameter :: radtodeg=180.0d0/pi
-
-! Compute Sun's declination
-declin=ASIN(sin(ls*degtorad)*sin(obliquity*degtorad))
-
-! Compute Local Time
-if (((-tan(lat*degtorad)*tan(declin)).gt.1.).or.(((-tan(lat*degtorad)*tan(declin)).lt.(-1.)))) then
-  ! Detect Polar Day and Polar Night (values out of the arccos definition domain)
-  lt_term = missval
-  
-else
-
-  if (trim(sunterm).eq."sunrise") then
-  ! Local Time of the Morning Terminator
-    lt_term = 24*(0.5 - ACOS(-tan(lat*degtorad)*tan(declin))/(2.*pi))
-  else ! if trim(sunterm).eq."sunset"
-  ! Local Time of the Evening Terminator
-    lt_term = 24*(0.5 + ACOS(-tan(lat*degtorad)*tan(declin))/(2.*pi))
-  endif
-  
-endif
-
-end subroutine LTterminator
-
-!******************************************************************************
-subroutine sol2ls(sol,ls)
-!  convert a given martian day number (sol)
-!  into corresponding solar longitude, Ls (in degr.),
-!  where sol=0=Ls=0 is the
-!  northern hemisphere spring equinox.
-
-implicit none
-
-!input 
-real,intent(in) :: sol
-!output 
-real,intent(out) :: ls
-
-!local variables
-double precision :: year_day,peri_day,timeperi,e_elips
-parameter (year_day=668.6d0) ! number of martian days (sols) in a martian year
-parameter (peri_day=485.35d0) ! perihelion date (in sols)
-parameter (e_elips=0.09340d0) ! orbital eccentricity
-parameter (timeperi=1.90258341759902d0) ! timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
-
-double precision :: pi,radtodeg
-parameter (pi=3.14159265358979d0)
-parameter (radtodeg=57.2957795130823d0) !  radtodeg: 180/pi
-
-double precision :: zanom,xref,zx0,zdx,zteta,zz
-!  xref: mean anomaly, zx0: eccentric anomaly, zteta: true anomaly        
-integer :: iter
-
-zz=(sol-peri_day)/year_day
-zanom=2.*pi*(zz-nint(zz))
-xref=dabs(zanom)
-
-!  The equation zx0 - e * sin (zx0) = xref, solved by Newton
-zx0=xref+e_elips*dsin(xref)
-iter=1
-iteration:do while (iter.le.10)
-   zdx=-(zx0-e_elips*dsin(zx0)-xref)/(1.-e_elips*dcos(zx0))
-   if(dabs(zdx).le.(1.d-7)) then ! typically, 2 or 3 iterations are enough
-     EXIT iteration 
-   endif
-   zx0=zx0+zdx
-   iter=iter+1
-enddo iteration
-zx0=zx0+zdx
-if(zanom.lt.0.) zx0=-zx0
-
-! compute true anomaly zteta, now that eccentric anomaly zx0 is known
-zteta=2.*datan(dsqrt((1.+e_elips)/(1.-e_elips))*dtan(zx0/2.))
-
-! compute Ls
-ls=real(zteta-timeperi)
-if(ls.lt.0.) ls=ls+2.*real(pi)
-if(ls.gt.2.*pi) ls=ls-2.*real(pi)
-! convert Ls in deg.
-ls=real(radtodeg)*ls
-
-end subroutine sol2ls
-
-!******************************************************************************
Index: unk/LMDZ.MARS/util/terminator.def
===================================================================
--- /trunk/LMDZ.MARS/util/terminator.def	(revision 2288)
+++ 	(revision )
@@ -1,24 +1,0 @@
-diagfi.nc
-tsurf
-ps
-temp
-
-sunrise
-0
-
------------------------------------------------------------------------
-ABOVE is the list of inputs to be fed to "terminator.e" if you don t want
-to reply directly to the program:
-
-1) complete name of file to be read
-2) List of X variables to be kept (X lines) or 'all'
-3) blank line at the end
-4) terminator type (sunrise or sunset)
-5) starttimeoffset :
-    - for a diagfi/concat : sol value at the beginning of the run, wrt Ls=0
-    - for a stats : sol value at the middle of the run, wrt Ls=0
-
-USE :
-terminator.e < terminator.def
-
-