Changeset 2434 for trunk/LMDZ.MARS/util/concatnc.F90

Timestamp:

Nov 17, 2020, 1:36:34 PM (4 years ago)

Author:

abierjon

Message:

Mars GCM:
Update utilities :
concat.F90

• rewrite and simplify the handling of time and offset so that any file can be concatenated, including files from different years or stats file
• use modulo to shift starting sols in concat.nc to values between 0 and 669
• add a "adls" option in addition to "sol" or "ls" to add a Ls 1D variable

while keeping "sol" as the Time axis

• add conservation of altitude attributes (long_name,units,positive)
• enable absence of both hybrid (aps,bps) and sigma levels

solzenangle.F90

• improve calculation for 1 sol file (stats) to use all local time data
• read the starting sol in the file instead of asking it to the user (except for stats file) ; keep the possibility for user to change it
• ask mean Ls value for stats file instead of sol number
• fix crash when using "all variable" option (ticket #66)
• fix bug on aps,bps by using GCM_layers dimension instead of altitude

localtime.F90

• fix crash when using "all variable" option (ticket #66)
• fix bug on aps,bps by using GCM_layers dimension instead of altitude

For all 3 utilities :

• handle all ndim cases for the variables (ticket #52)
• change "title" attribute into "long_name" by default (ticket #48)
• extend size of long_name character string (ticket #57)
• remove the use of #ifdef NC_DOUBLE (ticket #67)

AB

File:

• trunk/LMDZ.MARS/util/concatnc.F90 (modified) (46 diffs)

trunk/LMDZ.MARS/util/concatnc.F90

@@ -34 +34 @@
 character (len=50), dimension(:), allocatable :: var
 ! var(): name(s) of variable(s) that will be concatenated
-character (len=50) :: tmpvar,tmpfile,title,units
+character (len=100) :: tmpvar,tmpfile,long_name,units
 ! tmpvar(): used to temporarily store a variable name
 ! tmpfile(): used to temporarily store a file name
-! title(): [netcdf] title attribute
+! long_name(): [netcdf] long_name 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=4) :: axis
-! axis: "ls" or "sol"
+! axis:  "sol" or "ls" or "adls"
 integer :: nid,ierr,miss
 ! nid: [netcdf] file ID #
@@ -115 +117 @@
 real :: ctlsol
 ! ctlsol: starting sol value of the file, stored in the variable ctl()
-real :: memotime
-! memotime: (cumulative) time value, in martian days (sols)
+real :: previous_last_output_time, output_time
+! previous_last_output_time: last time stored in concat.nc after treating a file 
+logical :: firstsol
+! firstsol:  true if the first file's starting day must be used in the output
+real :: startsol
+! startsol: sol of the firstday as requested by the users (used if firstsol="y")
 real :: starttimeoffset
-! starttimeoffset: offset given by the user for the output time axis
-character (len=1) :: firstsol
-! firstsol: =[y/n], to know whether the user offset or the first file's starting day must be used in the output
+! starttimeoffset: offset given by the user for the output time axis (0 if ! firstsol="n")
 real :: missing
 !PARAMETER(missing=1E+20)
@@ -126 +130 @@
 real, dimension(2) :: valid_range
 ! valid_range(2): [netcdf] interval in which a value is considered valid
+real :: year_day = 669.
 
 !==============================================================================
@@ -154 +159 @@
 !==============================================================================
 
-write(*,*)
-write(*,*) "Starting day to be stored in the output file time axis?"
-write(*,*) "(e.g.: 100 if you want the output file to start at time=100 sols)"
-write(*,*) "Your answer will bypass any starting day value stored in the 1st"
-write(*,*) "input file. Conversely, just type return if you want this 1st file's"
-write(*,*) "value to be kept for the output file starting day."
-read(*,*,iostat=ierr) starttimeoffset
-if (ierr.ne.0) then
+ write(*,*) "Starting time in the concat file ?"
+ write(*,*) "   If you really wish to change it, write it now (not recommanded!)"
+ write(*,*) "   If you do not want to change it,  just write 'n' (recommanded)"
+
+read(*,*,iostat=ierr) startsol 
+
+firstsol= (ierr.eq.0)
+if (firstsol) then
    ! if nothing or a character that is not a float is read
    write(*,*) "1st input file's starting day will serve as starting day for the outfile."
-   firstsol='n'
 else ! if the user gave a number
    write(*,*) "Your input day will serve as starting day for the outfile."
-   firstsol='y'
 endif
 
@@ -190 +193 @@
 
 write(*,*)
-write(*,*) "Warning: to read the result with grads, choose sol"
-write(*,*) "Warning: use ferret to read the non linear scale ls"
-write(*,*) "Which time axis should be given in the output? (sol/ls)"
+write(*,*) "Which time axis should be given in the output? (sol/ls/adls)"
+write(*,*) "        (Warning: to read the result with grads, choose sol)"
+write(*,*) "  To keep sol as the time axis, and just add Ls as a variable, type adls  "
 read(*,*) axis
 ! loop as long as axis is neither "sol" nor "ls"
-do while ((axis/="sol").AND.(axis/="ls"))
+do while ((axis/="sol").AND.(axis/="ls").AND.(axis/="adls"))
    read(*,*) axis
 enddo
@@ -255 +258 @@
 if (tmpvar=="all") then
    if (axis=="ls") then
-!      write(*,*) "Do you want to keep the original file? (y/n)"
-!      read(*,*) ccopy
-!      if ((ccopy=="n").or.(ccopy=="N")) then
-!         do i=1,nbfile
-!            ierr=NF_CLOSE(nid)
-!            ierr = NF_OPEN(file(1),NF_WRITE,nid)
-!            call change_time_axis(nid,ierr)
-!            ierr=NF_CLOSE(nid)
-!            STOP ""
-!         enddo 
-!      else
          nbvar=nbvarfile-Nnotconcat
          do j=Nnotconcat+1,nbvarfile
             ierr=nf_inq_varname(nid,j,var(j-Nnotconcat))
          enddo
-!      endif ! of if ((ccopy=="n").or.(ccopy=="N"))
    endif ! of if (axis=="ls")
 ! Variables names from the file are catched
@@ -394 +385 @@
       endif
    endif
-!  write(*,*) "controle: ",controle
 
    if (ctllen .ne. 0) then ! variable controle
       allocate(ctl(ctllen))
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
-#else
       ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
-#endif
       ctlsol = ctl(4)
    endif
@@ -425 +411 @@
    allocate(time(timelen))
 
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
-#else
    ierr = NF_GET_VAR_REAL(nid,timevar,time)
-#endif
-
 
 !==============================================================================
@@ -437 +418 @@
 !==============================================================================
 
-   if (i==1) then ! First call; initialize/allocate
+   if (i==1) then ! First File ; initialize/allocate
       memolatlen=latlen
       memolonlen=lonlen
@@ -445 +426 @@
       allocate(lon(lonlen))
       allocate(alt(altlen))
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat)
-      ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon)
-      ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt)
-#else
       ierr = NF_GET_VAR_REAL(nid,latvar,lat)
       ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
       ierr = NF_GET_VAR_REAL(nid,altvar,alt)
-#endif
+!   Get altitude attributes to handle files with any altitude type
+      ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
+      if (ierr.ne.nf_noerr) then
+      ! if no attribute "long_name", try "title"
+        ierr=nf_get_att_text(nid,altvar,"title",altlong_name)
+      endif
+      ierr=nf_get_att_text(nid,altvar,'units',altunits)
+      ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
       write(*,*)
+      
       if (ctllen .ne. 0) then
-         if (int(starttimeoffset) .ne. (int(ctlsol)+int(time(1))) ) then
-           write(*,*) "WARNING: Starting day of the first file is not ",&
-                       int(starttimeoffset)," but ",int(ctlsol)+int(time(1)),"!"
+        
+         if (firstsol) then
+           starttimeoffset = float(int(ctlsol)+int(time(1) + ctl(27))) - startsol
+         else ! if firstsol.eq..false.
+!          if there is no user chosen first sol, nevertheless use
+!          starttimeoffset to lower the sol number (same exact season)
+           startsol = modulo(int(int(ctlsol)+int(time(1) + ctl(27))),int(year_day)) 
+           starttimeoffset = float(int(ctlsol)+int(time(1) + ctl(27))) - startsol
          endif
-        
-         if (firstsol.eq.'y') then
-           starttimeoffset = float(int(ctlsol)+int(time(1))) + ctl(27) - starttimeoffset
-         else ! if firstsol.eq.'n'
-           starttimeoffset = 0
-         endif
-
-         memotime=float(int(ctlsol)+int(time(1))) + ctl(27)
+
+!        artificially estimating "previous_last_output_time" for first file read: 
+         previous_last_output_time= startsol
+
          ctl(4) = 0.  ! values written in the output
          ctl(27) = 0. ! file controle variable (resp. day_ini and hour_ini)
          
       else ! if no variable "controle"
-         if (int(starttimeoffset) .ne. (int(time(1))) ) then
-           write(*,*) "WARNING: Starting day of the first file is not ",&
-                       int(starttimeoffset)," but ",int(time(1)),"!"
-         endif
-
-         if (firstsol.eq.'y') then
-           starttimeoffset = int(time(1)) - starttimeoffset
-         else ! if firstsol.eq.'n'
+         if (firstsol) then
+           starttimeoffset = int(time(1)) - startsol
+         else ! if firstsol.eq.false
            starttimeoffset = 0
          endif
 
-         memotime=float(int(time(1)))
+!        artificially estimating "previous_last_output_time" for first file read: 
+         previous_last_output_time=time(1) -(time(2) - time(1)) - starttimeoffset 
       endif
+      write(*,*) "Original starting day of 1st input file", previous_last_output_time + starttimeoffset
+      write(*,*) "Starting day of the concat file : ", previous_last_output_time 
+
         
    ! Initialize output file's lat,lon,alt and time dimensions
       write(*,*)
-      call initiate(filename,lat,lon,alt,ctl,GCM_layers,nout,&
-       latdimout,londimout,altdimout,timedimout,&
+      call initiate(filename,lat,lon,alt,ctl,GCM_layers,&
+       altlong_name,altunits,altpositive,&
+       nout,latdimout,londimout,altdimout,timedimout,&
        layerdimout,interlayerdimout,timevarout,axis)
    ! Initialize output file's aps,bps,ap,bp and phisinit variables
-     call init2(nid,lonlen,latlen,altlen,GCM_layers,&
+      call init2(nid,lonlen,latlen,altlen,GCM_layers,&
                 nout,londimout,latdimout,altdimout,&
                 layerdimout,interlayerdimout)
    
-   else ! Not a first call,
+   else ! Not first file to concatenate
    ! Check that latitude,longitude and altitude of current input file
    ! are identical to those of the output file
@@ -521 +506 @@
 
    rep=0
-   write(*,*)
-   write(*,'(a3,1x,f6.1)') "Sol",memotime-starttimeoffset
-!   if (ctllen.ne.0) write(*,'(a30,1x,f6.1)') "In the current file's ctl: Sol",ctlsol
-
-   ! Add (memotime)/(ctlsol) offset and write "concatenated" time values to output file
-   if (ctllen.ne.0) then
-     do k=reptime+1,reptime+timelen
+
+   do k=reptime+1,reptime+timelen
        rep=rep+1
-       if ((time(rep)+ctlsol).le.memotime) then
-         write(*,*) "ERROR : the time intervals of the files are not dissociated"
-         stop ""
-       endif
-#ifdef NC_DOUBLE
-       ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,k,1,ctlsol+time(rep)-starttimeoffset)
-#else
-       ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,ctlsol+time(rep)-starttimeoffset)
-#endif
-       ! 
-     enddo
-
-   ! Compute new time offset (for further concatenations)
-   memotime=ctlsol+time(timelen) ! the actual last day of the current file
-                                 ! becomes the reference time for next file
-                                 ! leaving potential "time holes" in the output time axis
-                                 ! at the junction of two files' time axes
-   
-   else ! if no variable "controle"
-   
-     do k=reptime+1,reptime+timelen
-        rep=rep+1
-#ifdef NC_DOUBLE
-        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,k,1,memotime+time(rep)-int(time(1))-starttimeoffset)
-#else
-        ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,memotime+time(rep)-int(time(1))-starttimeoffset)
-#endif
-     enddo
-     
-   ! Compute new time offset (for further concatenations)
-   memotime=memotime+time(timelen)-time(1) ! files are concatenated one after another
-                                           ! even if their time axes normally don't 
-                                           ! directly follow each other
-   endif ! if ctllen.ne.0
+       if (ctllen.ne.0) then
+          output_time=time(rep)+ctlsol-starttimeoffset
+!         correction if the file to concatenate seems "before" previous file
+!         (for instance to concatenate diagfi from the previous year at the enf of a year)
+          do while (output_time.lt.previous_last_output_time)  
+                   output_time = output_time + year_day
+          end do
+       else ! if no control, just add timestep after timestep
+          output_time= previous_last_output_time + (time(2)-time(1)) &
+         + (time(rep)-time(1))
+       end if
+       if (rep.eq.1) write(*,*) "Sol", int(output_time)  
+
+       ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,output_time)
+   end do
+!  use the last output_time value to update memotime   
+   previous_last_output_time = output_time 
+
 
 !==============================================================================
@@ -577 +541 @@
       ierr=nf_inq_varid(nid,var(j),varid)
       if (ierr.NE.NF_NOERR) then
-         write(*,*) 'ERROR: Field <',var(j),'> not found in file'//file(i)
+         write(*,*) 'ERROR: Field <',var(j),'> not found in file '//file(i)
          stop ""
       endif
       ierr=nf_inq_varndims(nid,varid,ndim)
 
+      ! Check that it is a 1D, 3D or 4D variable
+      if ((ndim.ne.1).and.(ndim.ne.3).and.(ndim.ne.4)) then
+        write(*,*) "Error:",trim(var(j))," is not a time-dependent variable,"
+        write(*,*) "so it can't be concatenated."
+        write(*,*) "We'll skip that variable..."
+        CYCLE ! go directly to the next variable
+      endif
+      
 !==============================================================================
 ! 2.5.2 Prepare things in order to read/write the variable
@@ -644 +616 @@
       if (i==1) then ! First call: write some definitions to output file
          units="                                                    "
-         title="                                                    "
-         ierr=nf_get_att_text(nid,varid,"title",title)
+         long_name="                                                    "
+         ierr=nf_get_att_text(nid,varid,"long_name",long_name)
+         if (ierr.ne.nf_noerr) then
+         ! if no attribute "long_name", try "title"
+           ierr=nf_get_att_text(nid,varid,"title",long_name)
+         endif
          ierr=nf_get_att_text(nid,varid,"units",units)
-         call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr)
+         call def_var(nout,var(j),long_name,units,ndim,dim,varidout,ierr)
       else
          ierr=NF_INQ_VARID(nout,var(j),varidout)
@@ -656 +632 @@
 !==============================================================================
 
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid,varid,var3d)
-      ierr= NF_PUT_VARA_DOUBLE(nout,varidout,corner,edges,var3d)
-      miss=NF_GET_ATT_DOUBLE(nid,varid,"missing_value",missing)
-      miss=NF_GET_ATT_DOUBLE(nid,varid,"valid_range",valid_range)
-#else
       ierr = NF_GET_VAR_REAL(nid,varid,var3d)
       ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3d)
       miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
       miss=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
-#endif
 
       if (ierr.ne.NF_NOERR) then
@@ -700 +669 @@
 !==============================================================================
 
-if (axis=="ls") then
-   write(*,*);write(*,*) "Converting the time axis from sol to Ls..."
-   call change_time_axis(nout,ierr)
+if ((axis=="ls").or.(axis=="adls")) then
+   call change_time_axis(nout,ierr,axis)
 endif
 
@@ -712 +680 @@
 
 !******************************************************************************
-subroutine initiate (filename,lat,lon,alt,ctl,GCM_layers,nout,&
-         latdimout,londimout,altdimout,timedimout,&
-         layerdimout,interlayerdimout,timevarout,axis)
+subroutine initiate (filename,lat,lon,alt,ctl,GCM_layers,&
+                     altlong_name,altunits,altpositive,&
+                     nout,latdimout,londimout,altdimout,timedimout,&
+                     layerdimout,interlayerdimout,timevarout,axis)
 !==============================================================================
 ! Purpose:
@@ -741 +710 @@
 ! ctl(): controle
 integer,intent(in) :: GCM_layers ! number of GCM layers
+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
@@ -787 +762 @@
 ierr = NF_DEF_DIM(nout, "longitude", size(lon), londimout)
 ierr = NF_DEF_DIM(nout, "altitude", size(alt), altdimout)
-if (size(ctl).ne.0) ierr = NF_DEF_DIM(nout, "index", size(ctl), ctldimout)
+if (ctllen.ne.0) ierr = NF_DEF_DIM(nout, "index", size(ctl), ctldimout)
 ierr = NF_DEF_DIM(nout, "Time", NF_UNLIMITED, timedimout)
 ierr = NF_DEF_DIM(nout, "GCM_layers", GCM_layers, layerdimout)
@@ -798 +773 @@
 ! 3. Write "Time" and its attributes
 !==============================================================================
-if (axis=="sol") then
+if (axis=="ls") then
+  call def_var(nout,"Time","Solar longitude","degree",1,&
+             (/timedimout/),timevarout,ierr)
+else 
   call def_var(nout,"Time","Time","days since 0000-00-0 00:00:00",1,&
              (/timedimout/),timevarout,ierr)
-else ! Ls
-  call def_var(nout,"Time","Solar longitude","days since 0000-00-0 00:00:00",1,&
-             (/timedimout/),timevarout,ierr)
 endif
 !==============================================================================
@@ -812 +787 @@
              (/latdimout/),nvarid,ierr)
 
-#ifdef NC_DOUBLE
-ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lat)
-#else
 ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
-#endif
 
 !==============================================================================
@@ -825 +796 @@
              (/londimout/),nvarid,ierr)
 
-#ifdef NC_DOUBLE
-ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lon)
-#else
 ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
-#endif
 
 !==============================================================================
@@ -838 +805 @@
 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",8,"altitude")
-ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',2,"km")
-ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',2,"up")
+
+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
 
 !==============================================================================
@@ -861 +820 @@
 !==============================================================================
 
-if (size(ctl).ne.0) then
+if (ctllen.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")
+
+   ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",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
 endif
 
@@ -890 +841 @@
 !==============================================================================
 ! Purpose:
-! Copy ap() , bp(), aps(), bps(), aire() and phisinit()
-! from input file to outpout file
+! Copy ap(), bp(), aps(), bps(), aire() and phisinit()
+! from input file to output file
 !==============================================================================
 ! Remarks:
@@ -907 +858 @@
 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 latitude
+integer, intent(in) :: altlen ! # of grid points along altitude
 integer, intent(in) :: GCM_layers ! # of GCM atmospheric layers
 integer, intent(in) :: outfid ! NetCDF output file ID
@@ -931 +882 @@
 logical :: hybrid ! are "aps" and "bps" available ?
 logical :: apbp ! are "ap" and "bp" available ?
+logical :: sig ! is "sigma" available ?
 
 !==============================================================================
@@ -976 +928 @@
   write(*,*) "init2: failed to allocate ap!"
   stop
+endif
+ierr=NF_INQ_VARID(infid,"ap",tmpvarid)
+if (ierr.ne.NF_NOERR) then
+  write(*,*) "Ooops. Failed to get ap ID. OK."
+  apbp=.false.
 else
-  ierr=NF_INQ_VARID(infid,"ap",tmpvarid)
-  if (ierr.ne.NF_NOERR) then
-    write(*,*) "Ooops. Failed to get ap ID. OK."
-    apbp=.false.
-  else
-    ierr=NF_GET_VAR_REAL(infid,tmpvarid,ap)
-    apbp=.true.
-    if (ierr.ne.NF_NOERR) then
-      stop "Error: Failed reading ap"
-    endif
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,ap)
+  apbp=.true.
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed reading ap"
   endif
 endif
@@ -995 +946 @@
   write(*,*) "init2: failed to allocate bp!"
   stop
+endif
+ierr=NF_INQ_VARID(infid,"bp",tmpvarid)
+if (ierr.ne.NF_NOERR) then
+  write(*,*) "Ooops. Failed to get bp ID. OK."
+  apbp=.false.
 else
-  ierr=NF_INQ_VARID(infid,"bp",tmpvarid)
-  if (ierr.ne.NF_NOERR) then
-    write(*,*) "Ooops. Failed to get bp ID. OK."
-    apbp=.false.
-  else
-    ierr=NF_GET_VAR_REAL(infid,tmpvarid,bp)
-    apbp=.true.
-    if (ierr.ne.NF_NOERR) then
-      stop "Error: Failed reading bp"
-    endif
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,bp)
+  apbp=.true.
+  if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed reading bp"
   endif
 endif
@@ -1017 +967 @@
   endif
   ierr=NF_INQ_VARID(infid,"sigma",tmpvarid)
-  ierr=NF_GET_VAR_REAL(infid,tmpvarid,sigma)
-  if (ierr.ne.NF_NOERR) then
-    stop "init2 Error: Failed reading sigma"
+  if (ierr.ne.NF_NOERR) then
+    write(*,*) "Ooops. Failed to get sigma ID. OK."
+    sig=.false.
+  else
+    ierr=NF_GET_VAR_REAL(infid,tmpvarid,sigma)
+    sig=.true.
+    if (ierr.ne.NF_NOERR) then
+      stop "init2 Error: Failed reading sigma"
+    endif
   endif
 endif ! of if (.not.hybrid)
@@ -1075 +1031 @@
 
 ! 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 "init2 Error: Failed to write aps"
@@ -1092 +1044 @@
 
 ! 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 "init2 Error: Failed to write bps"
@@ -1110 +1058 @@
 
 ! write ap
-#ifdef NC_DOUBLE
-    ierr=NF_PUT_VAR_DOUBLE(outfid,apid,ap)
-#else
     ierr=NF_PUT_VAR_REAL(outfid,apid,ap)
-#endif
     if (ierr.ne.NF_NOERR) then
       stop "Error: Failed to write ap"
@@ -1127 +1071 @@
 
 !   write bp
-#ifdef NC_DOUBLE
-    ierr=NF_PUT_VAR_DOUBLE(outfid,bpid,bp)
-#else
     ierr=NF_PUT_VAR_REAL(outfid,bpid,bp)
-#endif
     if (ierr.ne.NF_NOERR) then
       stop "Error: Failed to write bp"
@@ -1137 +1077 @@
   endif ! of if (apbp)
 
-else
+else if (sig) then
 ! define sigma
   call def_var(nout,"sigma","sigma at midlayers"," ",1,&
@@ -1145 +1085 @@
   endif
 ! write sigma
-#ifdef NC_DOUBLE
-  ierr=NF_PUT_VAR_DOUBLE(outfid,sigmaid,sigma)
-#else
   ierr=NF_PUT_VAR_REAL(outfid,sigmaid,sigma)
-#endif
   if (ierr.ne.NF_NOERR) then
     stop "init2 Error: Failed to write sigma"
@@ -1168 +1104 @@
   
   ! write aire
-#ifdef NC_DOUBLE
-  ierr=NF_PUT_VAR_DOUBLE(outfid,tmpvarid,aire)
-#else
   ierr=NF_PUT_VAR_REAL(outfid,tmpvarid,aire)
-#endif
   if (ierr.ne.NF_NOERR) then
     stop "init2 Error: Failed to write aire"
@@ -1188 +1120 @@
 
   ! 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 "init2 Error: Failed to write phisinit"
@@ -1211 +1139 @@
 end subroutine init2
 !******************************************************************************
-subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
+subroutine def_var(nid,name,long_name,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"...
+! called "name"; along with its attributes "long_name", "units"...
 ! to a file (following the NetCDF format)
 !==============================================================================
@@ -1232 +1160 @@
 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) :: long_name
+! long_name(): [netcdf] variable's "long_name" attribute
 character (len=*), intent(in) :: units
 ! unit(): [netcdf] variable's "units" attribute
@@ -1249 +1177 @@
 
 ! 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,"long_name",len_trim(adjustl(long_name)),adjustl(long_name))
 ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
 
@@ -1264 +1188 @@
 end subroutine def_var
 !******************************************************************************
-subroutine change_time_axis(nid,ierr)
+subroutine change_time_axis(nid,ierr,axis)
 !==============================================================================
 ! Purpose: 
@@ -1285 +1209 @@
 integer, intent(out) :: ierr
 ! ierr: [netcdf]  return error code
+character (len=4),intent(in) :: axis
+! axis: "ls" or "sol"
 
 !==============================================================================
@@ -1315 +1241 @@
 ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
 allocate(time(timelen),ls(timelen))
-#ifdef NC_DOUBLE
-ierr = NF_GET_VAR_DOUBLE(nid,nvarid,time)
-#else
 ierr = NF_GET_VAR_REAL(nid,nvarid,time)
-#endif
 if (ierr.ne.NF_NOERR) then
   write(*,*) "ERROR in change_time_axis: Failed to load Time"
@@ -1353 +1275 @@
 !==============================================================================
 
-#ifdef NC_DOUBLE
-ierr = NF_PUT_VAR_DOUBLE(nid,nvarid,ls)
-#else
+if (axis.eq."ls") then 
+   write(*,*) "Converting the time axis from sol to Ls..."
+else ! if (axis.eq."adls")
+   write(*,*) "Adding Ls as a 1D time variable"
+   call def_var(nout,"Ls","Solar Longitude","degree",1, (/timedimout/),nvarid,ierr)
+end if
+
 ierr = NF_PUT_VAR_REAL(nid,nvarid,ls)
-#endif
 if (ierr.ne.NF_NOERR) then
-  write(*,*) "ERROR in change_time_axis: Failed to write Ls"
-  stop
+   write(*,*) "ERROR in change_time_axis: Failed to write Ls"
+   stop
 endif
 
@@ -1502 +1427 @@
 
 ! Write valid_range() attribute
-#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
@@ -1517 +1438 @@
 
 ! Write "missing_value" attribute
-#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