Changeset 2434 for trunk/LMDZ.MARS/util

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

Location:

trunk/LMDZ.MARS/util

Files:

• concatnc.F90 (modified) (46 diffs)
• concatnc.def (modified) (1 diff)
• localtime.F90 (modified) (44 diffs)
• solzenangle.F90 (modified) (60 diffs)
• solzenangle.def (modified) (1 diff)

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

trunk/LMDZ.MARS/util/concatnc.def

@@ -3 +3 @@
 diagfi12.nc        
 
-515              ! Starting sol to use in the output file?
-sol              ! output timescale ("sol" or "ls")
+n              ! Starting sol of the run stored in the first input file?
+sol              ! output timescale ("sol" or "Ls" or "adls")
 tsurf            
 ps
 
 
+
 -----------------------------------------------------------------------
-ABOVE is the list of inputs to be fed to "concatnc.e" if you don't want
+ABOVE is the list of inputs to be fed to "concatnc.e" if you don t want
 to reply directly to the program: 
 
 1) List of N files to be read  (N lines)
 2) blank line at the end
-3) Starting sol of the output file?
-(or a blank if you want to use the first file's
- starting sol = time(1)+controle(4) )
-4) Output timescale ("sol" or "ls")
-5) List of X variables to be kept (X lines) or 'all'
-6) blank line at the end
+3) Starting sol of the run stored in the first input file? juste write "n" if you 
+   wish to keep the date of the input files (recommended)
+4) List of X variables to be kept (X lines) or 'all'
+5) blank line at the end
 
 The use is : 

trunk/LMDZ.MARS/util/localtime.F90

@@ -14 +14 @@
 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=30), dimension(16) :: notprocessed
+! notprocessed(): names of the (16) variables that won't be processed
 character (len=50), dimension(:), allocatable :: var
-! var(): name(s) of variable(s) that will be concatenated
-character (len=50) :: tmpvar,title,units
+! var(): name(s) of variable(s) that will be processed
+character (len=100) :: tmpvar,long_name,units
 ! tmpvar(): used to temporarily store a variable name
-! title(): [netcdf] title attribute
+! long_name(): [netcdf] long_name attribute
 ! units(): [netcdf] units attribute
 character (len=100) :: filename,vartmp
@@ -58 +58 @@
 ! ctldim: [netcdf] "controle" dim ID
 ! timedim: [netcdf] "timedim" dim ID
+integer :: gcmlayerdim ! NetCDF dimension ID for # of layers in GCM
 integer :: latvar,lonvar,altvar,ctlvar,timevar
 ! latvar: [netcdf] ID of "latitude" variable
@@ -68 +69 @@
 ! latlen: # of elements of lat() array
 ! lonlen: # of elements of lon() array
-! altvar: # of elements of alt() array
-! ctlvar: # of elements of ctl() array
+! altlen: # of elements of alt() array
+! ctllen: # of elements of ctl() array
 ! timelen: # of elemnets of time() array
 ! timelen_tot:# =timelen or timelen+1 (if 1 more time to interpolate needed)
 ! timelen_lt: # of elemnets of time() array in output
+integer :: nhour,nsol
+integer :: GCM_layers ! number of GCM atmospheric layers (may not be
+! same as altlen if file is an output of zrecast)
 integer :: nout,latdimout,londimout,altdimout,timedimout,timevarout
-integer :: nhour,nsol
 ! nout: [netcdf] output file ID
 ! latdimout: [netcdf] output latitude (dimension) ID
@@ -81 +84 @@
 ! timedimout: [netcdf] output time (dimension) ID
 ! timevarout: [netcdf] ID of output "Time" variable
+integer :: layerdimout ! NetCDF dimension ID for # of layers in GCM
 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
+integer :: Nnotprocessed,var_ok
+! Nnotprocessed: # of (leading)variables that won't be concatenated
 ! var_ok: flag (0 or 1)
 integer, dimension(4) :: corner,edges,dim
@@ -143 +147 @@
 !==============================================================================
 
-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'
+notprocessed(1)='Time'
+notprocessed(2)='controle'
+notprocessed(3)='rlonu'
+notprocessed(4)='latitude'
+notprocessed(5)='longitude'
+notprocessed(6)='altitude'
+notprocessed(7)='rlatv'
+notprocessed(8)='aps'
+notprocessed(9)='bps'
+notprocessed(10)='ap'
+notprocessed(11)='bp'
+notprocessed(12)='soildepth'
+notprocessed(13)='cu'
+notprocessed(14)='cv'
+notprocessed(15)='aire'
+notprocessed(16)='phisinit'
+
 
 !==============================================================================
@@ -163 +169 @@
 !==============================================================================
 write(*,*)
-   Nnotconcat=0
+   Nnotprocessed=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
+   do inter=1,16
+      if (vartmp.eq.notprocessed(inter)) then
          var_ok=1 
-         Nnotconcat=Nnotconcat+1
+         Nnotprocessed=Nnotprocessed+1
       endif 
    enddo        
@@ -177 +183 @@
 enddo
 
-! Nnotconcat: # of variables that won't be processed
+! Nnotprocessed: # of variables that won't be processed
 ! nbvarfile: total # of variables in file
-allocate(var(nbvarfile-Nnotconcat),stat=ierr)
+allocate(var(nbvarfile-Nnotprocessed),stat=ierr)
 if (ierr.ne.0) then
-  write(*,*) "Error: failed allocation of var(nbvarfile-Nnotconcat)"
+  write(*,*) "Error: failed allocation of var(nbvarfile-Nnotprocessed)"
   write(*,*) "  nbvarfile=",nbvarfile
-  write(*,*) "  Nnotconcat=",Nnotconcat
+  write(*,*) "  Nnotprocessed=",Nnotprocessed
   stop
 endif
@@ -201 +207 @@
 
 if (tmpvar=="all") then
-         nbvar=nbvarfile-Nnotconcat
-         do j=Nnotconcat+1,nbvarfile
-            ierr=nf_inq_varname(nid,j,var(j-Nnotconcat))
+         nbvar=nbvarfile-Nnotprocessed
+         do j=Nnotprocessed+1,nbvarfile
+            ierr=nf_inq_varname(nid,j,var(j-Nnotprocessed))
          enddo
 ! Variables names from the file are stored in var()
-   nbvar=nbvarfile-Nnotconcat
+   nbvar=nbvarfile-Nnotprocessed
    do i=1,nbvar
-      ierr=nf_inq_varname(nid,i+Nnotconcat,var(i))
+      ierr=nf_inq_varname(nid,i+Nnotprocessed,var(i))
       write(*,'(a9,1x,i2,1x,a1,1x,a50)') "variable ",i,":",var(i)
    enddo
@@ -220 +226 @@
 !==============================================================================
  filename=file(1:len_trim(file)-3)//"_LT.nc"
- write(*,*) trim(filename)
+ write(*,*) "Output file :"//trim(filename)
 
 !==============================================================================
@@ -226 +232 @@
 !==============================================================================
 
-      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
+   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
 
 !==============================================================================
@@ -277 +283 @@
    ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
 !  write(*,*) "altlen: ",altlen
+
+! load size of aps() or sigma() (in case it is not altlen)
+   ! default is that GCM_layers=altlen
+   ! but for outputs of zrecast, it may be a different value
+   ierr=NF_INQ_DIMID(nid,"GCM_layers",gcmlayerdim)
+   if (ierr.ne.NF_NOERR) then
+     ! didn't find a GCM_layers dimension; therefore we have:
+     GCM_layers=altlen
+   else
+     ! load value of GCM_layers
+     ierr=NF_INQ_DIMLEN(nid,gcmlayerdim,GCM_layers)
+   endif
+!   write(*,*) "GCM_layers=",GCM_layers
 
    ierr=NF_INQ_DIMID(nid,"index",ctldim)
@@ -299 +318 @@
 !==============================================================================
 
-! 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
+   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
+
+   ierr = NF_GET_VAR_REAL(nid,latvar,lat)
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to load latitude"
+     write(*,*) NF_STRERROR(ierr)
+     stop
+   endif
+
+   ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to load longitude"
+     write(*,*) NF_STRERROR(ierr)
+     stop
+   endif
+
+   ierr = NF_GET_VAR_REAL(nid,altvar,alt)
+   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)
+   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,timevar,"title",long_name)
+   endif
+   ierr=nf_get_att_text(nid,altvar,'units',altunits)
+   ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
+
+   if (ctllen .ne. 0) then
+     ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
+     if (ierr.ne.0) then
+       write(*,*) "Error: failed to load controle"
+       write(*,*) NF_STRERROR(ierr)
+       stop
+     endif
+   endif ! of if (ctllen .ne. 0)
 
 !==============================================================================
@@ -391 +397 @@
 !  write(*,*) "timelen: ",timelen
 
-   ! Sanity Check: Does "Time" has a "title" attribute
+   ! Sanity Check: Does "Time" has a "long_name" 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
+   ierr=nf_get_att_text(nid,timevar,"long_name",long_name)
+   if (ierr.ne.nf_noerr) then
+   ! if no attribute "long_name", try "title"
+     ierr=nf_get_att_text(nid,timevar,"title",long_name)
+   endif
+   if ((ierr.EQ.NF_NOERR).and.(long_name.eq."Solar longitude")) then
      write(*,*) "ERROR: Time axis in input file is in Solar Longitude!"
      write(*,*) "       localtime requires sols as time axis!"
@@ -408 +418 @@
    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"
@@ -432 +438 @@
       write(*,*) 'list of selected local time (0<t<24)'
       do it=1,nhour
-         read(*,*) lt_hour(it)
+        read(*,*) lt_hour(it)
       end do
 
       if ((timelen.eq.12).and.(time(1).eq.2).and.(time(timelen).eq.24))then
-         write(*,*) 'We detect a ""stats"" file'
-         stats=.true.
-         timelen_lt=nhour
-         allocate(lt_out(timelen_lt),stat=ierr)
-         if (ierr.ne.0) then
-           write(*,*) "Error: failed to allocate lt_hour(nhour)"
-           stop
-         endif
-         do it=1,nhour
-           lt_out(it)=lt_hour(it)/24.
-         end do
+        write(*,*) 'We detect a ""stats"" file'
+        stats=.true.
+        timelen_lt=nhour
+        allocate(lt_out(timelen_lt),stat=ierr)
+        if (ierr.ne.0) then
+          write(*,*) "Error: failed to allocate lt_hour(nhour)"
+          stop
+        endif
+        do it=1,nhour
+          lt_out(it)=lt_hour(it)/24.
+        end do
 !        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
+        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.
@@ -469 +475 @@
           end do
         end do
-        end if
+      end if
 
       if (nsol.gt.1) then
@@ -499 +505 @@
 
    ! Initialize output file's lat,lon,alt and time dimensions
-     call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
+     call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,GCM_layers,&
            altlong_name,altunits,altpositive,&
-           nout,latdimout,londimout,altdimout,timedimout,timevarout)
+           nout,latdimout,londimout,altdimout,timedimout,&
+           layerdimout,timevarout)
 
    ! Initialize output file's aps,bps and phisinit variables
-     call init2(nid,lonlen,latlen,altlen,altdim,&
-                nout,londimout,latdimout,altdimout)
+     call init2(nid,lonlen,latlen,altlen,GCM_layers,&
+                nout,londimout,latdimout,altdimout,&
+                layerdimout)
 
 
@@ -515 +523 @@
 
      do it=1,timelen_lt
-#ifdef NC_DOUBLE
-        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,lt_hour(it))
-#else
         ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,lt_hour(it))
-#endif
      enddo
    else
      do it=1,timelen_lt
-#ifdef NC_DOUBLE
-        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,lt_out(it))
-#else
         ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,lt_out(it))
-#endif
         if (ierr.NE.NF_NOERR) then
           write(*,*) "Error , failed to write Time"
@@ -553 +553 @@
       ierr=nf_inq_varndims(nid,varid,ndim)
 
+      ! Check that it is a 3D or 4D variable
+      if (ndim.lt.3) then
+        write(*,*) "Error:",trim(var(j))," is nor a 3D neither a 4D variable"
+        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
@@ -599 +605 @@
       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)
+      units=" "
+      long_name=" "
+      ierr=nf_get_att_text(nid,timevar,"long_name",long_name)
+      if (ierr.ne.nf_noerr) then
+      ! if no attribute "long_name", try "title"
+        ierr=nf_get_att_text(nid,timevar,"title",long_name)
+      endif
+      ierr=nf_get_att_text(nid,varid,"units",units)
+      call def_var(nout,var(j),long_name,units,ndim,dim,varidout,ierr)
 
       
@@ -611 +621 @@
 !==============================================================================
 
-#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
 
 !==============================================================================
@@ -677 +681 @@
 !==============================================================================
 
-#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
@@ -712 +712 @@
 ierr=NF_CLOSE(nout)
 
+write(*,*) "Program completed !"
+
 contains
 
 !******************************************************************************
-Subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
+Subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,GCM_layers,&
                      altlong_name,altunits,altpositive,&
-                     nout,latdimout,londimout,altdimout,timedimout,timevarout)
+                     nout,latdimout,londimout,altdimout,timedimout,&
+                     layerdimout,timevarout)
 !==============================================================================
 ! Purpose:
@@ -744 +747 @@
 real, intent(in):: ctl(ctllen)
 ! ctl(): controle
+integer,intent(in) :: GCM_layers ! number of GCM layers
 character (len=*), intent(in) :: altlong_name
 ! altlong_name(): [netcdf] altitude "long_name" attribute
@@ -760 +764 @@
 integer, intent(out):: timedimout
 ! timedimout: [netcdf] "Time"  ID
+integer,intent(out) :: layerdimout
+! layerdimout: [netcdf] "GCM_layers" ID
 integer, intent(out):: timevarout
 ! timevarout: [netcdf] "Time" (considered as a variable) ID
@@ -821 +827 @@
 endif
 
+ierr = NF_DEF_DIM(nout, "GCM_layers", GCM_layers, layerdimout)
+
 ! End netcdf define mode
 ierr = NF_ENDDEF(nout)
@@ -843 +851 @@
              (/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>'
@@ -861 +865 @@
              (/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>'
@@ -879 +880 @@
 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))
@@ -892 +889 @@
 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 
+ierr = NF_PUT_VAR_REAL (nout,nvarid,alt) 
 if (ierr.NE.NF_NOERR) then
    WRITE(*,*)'initiate: error failed writing variable <altitude>'
@@ -911 +904 @@
    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
    if (ierr.NE.NF_NOERR) then
       WRITE(*,*)'initiate: error failed writing variable <controle>'
@@ -936 +921 @@
 end Subroutine initiate
 !******************************************************************************
-subroutine init2(infid,lonlen,latlen,altlen,altdimid, &
-                 outfid,londimout,latdimout,altdimout)
+subroutine init2(infid,lonlen,latlen,altlen,GCM_layers,&
+                 outfid,londimout,latdimout,altdimout,&
+                 layerdimout)
 !==============================================================================
 ! Purpose:
@@ -957 +943 @@
 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) :: GCM_layers ! # of GCM atmospheric layers
 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
+integer, intent(in) :: layerdimout ! GCM_layers dimension ID
 !==============================================================================
 ! Local variables:
@@ -971 +958 @@
 integer :: tmpdimid ! temporary dimension ID
 integer :: tmpvarid ! temporary variable ID
-logical :: phis, aps_ok, bps_ok ! are "phisinit" "aps" "bps" available ?
+logical :: phis, hybrid ! are "phisinit" and "aps"/"bps" available ?
 
 
@@ -979 +966 @@
 
 ! hybrid coordinate aps
-  allocate(aps(altlen),stat=ierr)
+!write(*,*) "aps: altlen=",altlen," GCM_layers=",GCM_layers
+allocate(aps(GCM_layers),stat=ierr)
+if (ierr.ne.0) then
+  write(*,*) "init2: failed to allocate aps!"
+  stop
+endif
+ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
+if (ierr.ne.NF_NOERR) then
+  write(*,*) "Ooops. Failed to get aps ID. OK."
+  hybrid=.false.
+else
+  hybrid=.true.
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
+  hybrid=.true.
+  if (ierr.ne.NF_NOERR) then
+    stop "init2 Error: Failed reading aps"
+  endif
+
+  ! hybrid coordinate bps
+!  write(*,*) "bps: altlen=",altlen," GCM_layers=",GCM_layers
+  allocate(bps(GCM_layers),stat=ierr)
   if (ierr.ne.0) then
-    write(*,*) "init2: failed to allocate aps(altlen)"
+    write(*,*) "init2: failed to allocate bps!"
     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.
+  ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
+  if (ierr.ne.NF_NOERR) then
+    write(*,*) "Ooops. Failed to get bps ID. OK."
+    hybrid=.false.
   else
     ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
@@ -1027 +1000 @@
       stop "init2 Error: Failed reading bps"
     endif
-    bps_ok=.true.
   endif
 endif
@@ -1059 +1031 @@
 !==============================================================================
 
-IF (aps_ok) then 
+IF (hybrid) then 
   ! define aps
   call def_var(outfid,"aps","hybrid pressure at midlayers"," ",1,&
-               (/altdimout/),apsid,ierr)
+               (/layerdimout/),apsid,ierr)
   if (ierr.ne.NF_NOERR) then
     stop "Error: Failed to def_var aps"
@@ -1068 +1040 @@
 
   ! 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)
+               (/layerdimout/),bpsid,ierr)
   if (ierr.ne.NF_NOERR) then
     stop "Error: Failed to def_var bps"
@@ -1087 +1053 @@
 
   ! 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)
+  
+  deallocate(aps)
+  deallocate(bps)
+ENDIF ! of IF (hybrid)
 
 !==============================================================================
@@ -1111 +1076 @@
 
   ! 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
 
+  deallocate(phisinit)
 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)
+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)
 !==============================================================================
@@ -1151 +1109 @@
 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
@@ -1168 +1126 @@
 
 ! 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))
 
@@ -1232 +1186 @@
 ! 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
@@ -1246 +1196 @@
   endif
 endif
-
+if (miss.eq.NF_NOERR) then
 ! 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

trunk/LMDZ.MARS/util/solzenangle.F90

@@ -15 +15 @@
 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=30), dimension(16) :: notprocessed
+! notprocessed(): names of the (16) variables that won't be processed
 character (len=50), dimension(:), allocatable :: var
-! var(): name(s) of variable(s) that will be concatenated
-character (len=50) :: tmpvar,title,units
+! var(): name(s) of variable(s) that will be processed
+character (len=100) :: tmpvar,long_name,units
 ! tmpvar(): used to temporarily store a variable name
-! title(): [netcdf] title attribute
+! long_name(): [netcdf] long_name attribute
 ! units(): [netcdf] units attribute
 character (len=100) :: filename,vartmp
@@ -62 +62 @@
 ! ctldim: [netcdf] "controle" dim ID
 ! timedim: [netcdf] "timedim" dim ID
+integer :: gcmlayerdim ! NetCDF dimension ID for # of layers in GCM
 integer :: latvar,lonvar,altvar,ctlvar,timevar
 ! latvar: [netcdf] ID of "latitude" variable
@@ -72 +73 @@
 ! latlen: # of elements of lat() array
 ! lonlen: # of elements of lon() array
-! altvar: # of elements of alt() array
-! ctlvar: # of elements of ctl() array
+! altlen: # of elements of alt() array
+! ctllen: # 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 :: GCM_layers ! number of GCM atmospheric layers (may not be
+! same as altlen if file is an output of zrecast)
 integer :: nout,latdimout,londimout,altdimout,timedimout,timevarout
 ! nout: [netcdf] output file ID
@@ -85 +88 @@
 ! timedimout: [netcdf] output time (dimension) ID
 ! timevarout: [netcdf] ID of output "Time" variable
+integer :: layerdimout ! NetCDF dimension ID for # of layers in GCM
 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
+integer :: Nnotprocessed,var_ok
+! Nnotprocessed: # of (leading)variables that won't be concatenated
 ! var_ok: flag (0 or 1)
 integer, dimension(4) :: corner,edges,dim
@@ -107 +111 @@
 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 :: startsol ! starting date (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
@@ -121 +125 @@
 ! valid_range(2): [netcdf] interval in which a value is considered valid
 logical :: stats   ! stats=T when reading a "stats" kind of ile
+real :: year_day_gcm = 669. ! number of sols in a year in GCM
+
+
 
 
@@ -160 +167 @@
 !==============================================================================
 
-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'
+notprocessed(1)='Time'
+notprocessed(2)='controle'
+notprocessed(3)='rlonu'
+notprocessed(4)='latitude'
+notprocessed(5)='longitude'
+notprocessed(6)='altitude'
+notprocessed(7)='rlatv'
+notprocessed(8)='aps'
+notprocessed(9)='bps'
+notprocessed(10)='ap'
+notprocessed(11)='bp'
+notprocessed(12)='soildepth'
+notprocessed(13)='cu'
+notprocessed(14)='cv'
+notprocessed(15)='aire'
+notprocessed(16)='phisinit'
 
 !==============================================================================
@@ -180 +188 @@
 !==============================================================================
 write(*,*)
-   Nnotconcat=0
+   Nnotprocessed=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
+   do inter=1,16
+      if (vartmp.eq.notprocessed(inter)) then
          var_ok=1 
-         Nnotconcat=Nnotconcat+1
+         Nnotprocessed=Nnotprocessed+1
       endif 
    enddo        
@@ -194 +202 @@
 enddo
 
-! Nnotconcat: # of variables that won't be processed
+! Nnotprocessed: # of variables that won't be processed
 ! nbvarfile: total # of variables in file
-allocate(var(nbvarfile-Nnotconcat),stat=ierr)
+allocate(var(nbvarfile-Nnotprocessed),stat=ierr)
 if (ierr.ne.0) then
-  write(*,*) "Error: failed allocation of var(nbvarfile-Nnotconcat)"
+  write(*,*) "Error: failed allocation of var(nbvarfile-Nnotprocessed)"
   write(*,*) "  nbvarfile=",nbvarfile
-  write(*,*) "  Nnotconcat=",Nnotconcat
+  write(*,*) "  Nnotprocessed=",Nnotprocessed
   stop
 endif
@@ -218 +226 @@
 
 if (tmpvar=="all") then
-         nbvar=nbvarfile-Nnotconcat
-         do j=Nnotconcat+1,nbvarfile
-            ierr=nf_inq_varname(nid,j,var(j-Nnotconcat))
+         nbvar=nbvarfile-Nnotprocessed
+         do j=Nnotprocessed+1,nbvarfile
+            ierr=nf_inq_varname(nid,j,var(j-Nnotprocessed))
          enddo
 ! Variables names from the file are stored in var()
-   nbvar=nbvarfile-Nnotconcat
+   nbvar=nbvarfile-Nnotprocessed
    do i=1,nbvar
-      ierr=nf_inq_varname(nid,i+Nnotconcat,var(i))
+      ierr=nf_inq_varname(nid,i+Nnotprocessed,var(i))
       write(*,'(a9,1x,i2,1x,a1,1x,a50)') "variable ",i,":",var(i)
    enddo
@@ -290 +298 @@
 !  write(*,*) "altlen: ",altlen
 
+! load size of aps() or sigma() (in case it is not altlen)
+   ! default is that GCM_layers=altlen
+   ! but for outputs of zrecast, it may be a different value
+   ierr=NF_INQ_DIMID(nid,"GCM_layers",gcmlayerdim)
+   if (ierr.ne.NF_NOERR) then
+     ! didn't find a GCM_layers dimension; therefore we have:
+     GCM_layers=altlen
+   else
+     ! load value of GCM_layers
+     ierr=NF_INQ_DIMLEN(nid,gcmlayerdim,GCM_layers)
+   endif
+!   write(*,*) "GCM_layers=",GCM_layers
+
    ierr=NF_INQ_DIMID(nid,"index",ctldim)
    if (ierr.NE.NF_NOERR) then
@@ -304 +325 @@
       ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
    endif
-!  write(*,*) "controle: ",controle
+
 
 !==============================================================================
@@ -311 +332 @@
 !==============================================================================
 
-! First call; initialize/allocate
    allocate(lat(latlen),stat=ierr)
    if (ierr.ne.0) then
@@ -333 +353 @@
    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"
@@ -344 +360 @@
    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"
@@ -355 +367 @@
    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"
@@ -367 +375 @@
 ! 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)
 
    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"
@@ -403 +411 @@
 !  write(*,*) "timelen: ",timelen
 
-   ! Sanity Check: Does "Time" has a "title" attribute
+   ! Sanity Check: Does "Time" has a "long_name" 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
+   ierr=nf_get_att_text(nid,timevar,"long_name",long_name)
+   if (ierr.ne.nf_noerr) then
+   ! if no attribute "long_name", try "title"
+     ierr=nf_get_att_text(nid,timevar,"title",long_name)
+   endif
+   if ((ierr.EQ.NF_NOERR).and.(long_name.eq."Solar longitude")) then
      write(*,*) "ERROR: Time axis in input file is in Solar Longitude!"
      write(*,*) "       localtime requires sols as time axis!"
@@ -420 +432 @@
    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"
@@ -434 +442 @@
 ! 2.4.1 Select local times to be saved "Time" in output file
 !==============================================================================
-   write(*,*) "Planet side of the Solar Zenith Angle ?"
+   write(*,*) "Planet side of the Solar Zenith Angle ? (morning or evening)"
    read(*,*) planetside
    if ((trim(planetside).ne."morning").and.(trim(planetside).ne."evening")) then
@@ -448 +456 @@
      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
+
+!  Sol and season:  
+
+   if (ctllen .ne. 0) then
+      startsol = int(int(ctl(4))+int(time(1) + ctl(27)))
+   else 
+      startsol = int(time(1))
+   end if
 
    if ((timelen.eq.12).and.(time(1).eq.2).and.(time(timelen).eq.24))then
@@ -466 +478 @@
      nsol = int(time(timelen)) - int(time(1))
    end if
-
    allocate(intsol(nsol),stat=ierr)
    if (ierr.ne.0) then
@@ -477 +488 @@
      stop
    endif
+
+   if (stats) then
+         write(*,*) "What is the mean Ls of this stats file (degree) ?" 
+         read(*,*,iostat=ierr) tmpLs
+         startsol=0.
+   else 
+      write(*,*) "Beginning date (sol since Ls=0) of this file is :",startsol 
+      write(*,*) "If you do not want to change it,  just write 'n' (recommended)"
+      write(*,*) "If you really wish to change it, write it now (not recommended!)"
+      read(*,*,iostat=ierr) startsol
+      if((ierr.eq.0).and.(ctllen.ne.0))  ctl(4)=startsol
+      write(*,*) "Beginning date (sol since Ls=0) will be :",startsol 
+   end if
+
    do it=1,nsol
      intsol(it)= int(time(1)) + it-1.
      do ilon=1,lonlen
-!         For the computation of Ls, we try to take into account the rotation time
-!         of Mars it's supposed that the morning/evening
+!         For the computation of Ls, we try to take into account the
+!         rotation time of Mars. It is 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
+         tmpsol = intsol(it) + (6.-lon(ilon)/15.)/24. + startsol
        else !if trim(planetside).eq."evening"
-         tmpsol = intsol(it) + (18.-lon(ilon)/15.)/24. + starttimeoffset            
+         tmpsol = intsol(it) + (18.-lon(ilon)/15.)/24. + startsol            
        endif
-       call sol2ls(tmpsol,tmpLs)
+       if (.not.stats) call sol2ls(tmpsol,tmpLs)
        do ilat=1,latlen
 !           Compute the Local Time of the solar zenith angle at a given longitude
@@ -541 +566 @@
 
   ! Initialize output file's lat,lon,alt and time dimensions
-   call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
+   call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,GCM_layers,&
          altlong_name,altunits,altpositive,&
-         nout,latdimout,londimout,altdimout,timedimout,timevarout)
-
+         nout,latdimout,londimout,altdimout,timedimout,&
+         layerdimout,timevarout)
+         
   ! Initialize output file's aps,bps and phisinit variables
-   call init2(nid,lonlen,latlen,altlen,altdim,&
-              nout,londimout,latdimout,altdimout)
-
+   call init2(nid,lonlen,latlen,altlen,GCM_layers,&
+              nout,londimout,latdimout,altdimout,&
+              layerdimout)
 
 !==============================================================================
@@ -557 +583 @@
 
      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
+       ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,intsol(it)*24.)
      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
+       ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,intsol(it))
+       if (ierr.NE.NF_NOERR) then
+         write(*,*) "Error , failed to write Time"
+         stop
+       endif
      enddo
    end if
@@ -584 +602 @@
 
 !==============================================================================
-! 2.5.1 Check that variable to be read existe in input file
+! 2.5.1 Check that variable to be read exists in input file
 !==============================================================================
 
@@ -595 +613 @@
       ierr=nf_inq_varndims(nid,varid,ndim)
 
+      ! Check that it is a 3D or 4D variable
+      if (ndim.lt.3) then
+        write(*,*) "Error:",trim(var(j))," is nor a 3D neither a 4D variable"
+        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
@@ -642 +667 @@
 
       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)
 
       
@@ -653 +682 @@
 !==============================================================================
 
-#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
-
-      ! If there is no missing value attribute in the input variable, create one
+
+     ! If there is no missing value attribute in the input variable, create one
       if (miss.ne.NF_NOERR) then
          missing = miss_lt
@@ -673 +696 @@
 !==============================================================================
 
-      do ilon=1,lonlen
+       do ilon=1,lonlen
 !        Recast GCM local time (sol decimal value) at each longitude :
          do it=1,timelen  ! loop time input file
@@ -688 +711 @@
                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
+             if (nsol.eq.1) then 
+!              If 1 sol file (stats), use all data from that sol.
+               if(lt_outc(it).gt.lt_gcm(timelen_tot)) lt_outc(it)=lt_outc(it)-1
+               if(lt_outc(it).lt.lt_gcm(1)) lt_outc(it)=lt_outc(it)+1
+             else
+!              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
+             end if
            enddo ! local time
 
@@ -710 +739 @@
              enddo ! local time
 
-           else if  (ndim==4) then
+           else if (ndim==4) then
              do ialt=1,altlen
                do it=1,timelen ! loop time input file
@@ -725 +754 @@
                enddo ! local time   
              enddo ! alt
-           endif
+           endif ! ndim
 
          enddo ! lat
-      end do ! lon
+       end do ! lon
 
 
@@ -735 +764 @@
 !==============================================================================
 
-#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
+       ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3d_lt)
+
+       if (ierr.ne.NF_NOERR) then
          write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
          stop ""
-      endif
+       endif
 
       ! Write "missing_value" (and "valid_range" if there is one) attributes in output file
-      call missing_value(nout,varidout,validr,miss,valid_range,missing)
+       call missing_value(nout,varidout,validr,miss,valid_range,missing)
 
       ! free var3d() array
-      deallocate(var3d)
-      deallocate(var3d_lt)
+       deallocate(var3d)
+       deallocate(var3d_lt)
 
    enddo ! of do j=1,nbvar
@@ -780 +805 @@
    write(szastr,*) sza
    if (trim(planetside).eq."morning") then
-     title="Morning-side LTST at sza="//trim(adjustl(szastr))//"deg"
+     long_name="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
+     long_name="Evening-side LTST at sza="//trim(adjustl(szastr))//"deg"            
+   endif
+   call def_var(nout,tmpvar,long_name,units,3,dim,varidout,ierr)
+
    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)
@@ -798 +819 @@
    validr=NF_NOERR+1
    miss=NF_NOERR
-   call missing_value(nout,varidout,validr,miss,(0,0),miss_lt)
+   valid_range(:)=0.
+   call missing_value(nout,varidout,validr,miss,valid_range,miss_lt)
    write(*,*)"with missing value = ",miss_lt
 !==============================================================================
@@ -837 +859 @@
 !******************************************************************************
 
-subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
+subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,GCM_layers,&
                      altlong_name,altunits,altpositive,&
-                     nout,latdimout,londimout,altdimout,timedimout,timevarout)
+                     nout,latdimout,londimout,altdimout,timedimout,&
+                     layerdimout,timevarout)
 !==============================================================================
 ! Purpose:
@@ -866 +889 @@
 real, intent(in):: ctl(ctllen)
 ! ctl(): controle
+integer,intent(in) :: GCM_layers ! number of GCM layers
 character (len=*), intent(in) :: altlong_name
 ! altlong_name(): [netcdf] altitude "long_name" attribute
@@ -882 +906 @@
 integer, intent(out):: timedimout
 ! timedimout: [netcdf] "Time"  ID
+integer,intent(out) :: layerdimout
+! layerdimout: [netcdf] "GCM_layers" ID
 integer, intent(out):: timevarout
 ! timevarout: [netcdf] "Time" (considered as a variable) ID
@@ -943 +969 @@
 endif
 
+ierr = NF_DEF_DIM(nout, "GCM_layers", GCM_layers, layerdimout)
+if (ierr.NE.NF_NOERR) then
+   WRITE(*,*)'initiate: error failed to define dimension <GCM_layers>'
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
 ! End netcdf define mode
 ierr = NF_ENDDEF(nout)
@@ -950 +983 @@
    stop ""
 endif
-
 !==============================================================================
 ! 3. Write "Time" and its attributes
@@ -960 +992 @@
 ierr = NF_REDEF (nout)
 ierr = NF_PUT_ATT_TEXT(nout,timevarout,'comment',135,&
-      "The true unit of the variable Time is the integer value of sol at lon=0deg. A false unit is put to enable reading from Grads or Ferret")
+"The true unit of the variable Time is the integer value &
+ of sol at lon=0deg. A false unit is put to enable reading from Grads or Ferret")
 ! End netcdf define mode
 ierr = NF_ENDDEF(nout)
@@ -971 +1004 @@
              (/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>'
@@ -989 +1018 @@
              (/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>'
@@ -1007 +1032 @@
 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))
@@ -1020 +1041 @@
 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 
+ierr = NF_PUT_VAR_REAL (nout,nvarid,alt) 
 if (ierr.NE.NF_NOERR) then
    WRITE(*,*)'initiate: error failed writing variable <altitude>'
@@ -1039 +1056 @@
    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
    if (ierr.NE.NF_NOERR) then
       WRITE(*,*)'initiate: error failed writing variable <controle>'
@@ -1065 +1074 @@
 
 !******************************************************************************
-subroutine init2(infid,lonlen,latlen,altlen,altdimid, &
-                 outfid,londimout,latdimout,altdimout)
+subroutine init2(infid,lonlen,latlen,altlen,GCM_layers, &
+                 outfid,londimout,latdimout,altdimout, &
+                 layerdimout)
 !==============================================================================
 ! Purpose:
@@ -1086 +1096 @@
 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) :: GCM_layers ! # of GCM atmospheric layers
 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
+integer, intent(in) :: layerdimout ! GCM_layers dimension ID
 !==============================================================================
 ! Local variables:
@@ -1100 +1111 @@
 integer :: tmpdimid ! temporary dimension ID
 integer :: tmpvarid ! temporary variable ID
-logical :: phis, aps_ok, bps_ok ! are "phisinit" "aps" "bps" available ?
+integer :: tmplen ! temporary variable length
+logical :: phis,hybrid ! are "phisinit" and "aps"/"bps" available ?
 
 
@@ -1108 +1120 @@
 
 ! hybrid coordinate aps
-  allocate(aps(altlen),stat=ierr)
+!write(*,*) "aps: altlen=",altlen," GCM_layers=",GCM_layers
+allocate(aps(GCM_layers),stat=ierr)
+if (ierr.ne.0) then
+  write(*,*) "init2: failed to allocate aps!"
+  stop
+endif
+ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
+if (ierr.ne.NF_NOERR) then
+  write(*,*) "Ooops. Failed to get aps ID. OK."
+  hybrid=.false.
+else
+  ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
+  hybrid=.true.
+  if (ierr.ne.NF_NOERR) then
+    stop "init2 Error: Failed reading aps"
+  endif
+
+  ! hybrid coordinate bps
+!  write(*,*) "bps: altlen=",altlen," GCM_layers=",GCM_layers
+  allocate(bps(GCM_layers),stat=ierr)
   if (ierr.ne.0) then
-    write(*,*) "init2: failed to allocate aps(altlen)"
+    write(*,*) "init2: failed to allocate bps!"
     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.
+  ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
+  if (ierr.ne.NF_NOERR) then
+    write(*,*) "Ooops. Failed to get bps ID. OK."
+    hybrid=.false.
   else
     ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
@@ -1156 +1153 @@
       stop "init2 Error: Failed reading bps"
     endif
-    bps_ok=.true.
   endif
 endif
@@ -1188 +1184 @@
 !==============================================================================
 
-IF (aps_ok) then 
+IF (hybrid) then 
   ! define aps
   call def_var(outfid,"aps","hybrid pressure at midlayers"," ",1,&
-               (/altdimout/),apsid,ierr)
+               (/layerdimout/),apsid,ierr)
   if (ierr.ne.NF_NOERR) then
     stop "Error: Failed to def_var aps"
@@ -1197 +1193 @@
 
   ! 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)
+               (/layerdimout/),bpsid,ierr)
   if (ierr.ne.NF_NOERR) then
     stop "Error: Failed to def_var bps"
@@ -1216 +1206 @@
 
   ! 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)
+  
+  deallocate(aps)
+  deallocate(bps)
+ENDIF ! of IF (hybrid)
 
 !==============================================================================
@@ -1240 +1229 @@
 
   ! 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
 
+deallocate(phisinit)
 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)
+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)
 !==============================================================================
@@ -1281 +1262 @@
 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
@@ -1298 +1279 @@
 
 ! 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))
 
@@ -1363 +1340 @@
 ! 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
@@ -1380 +1353 @@
 ! 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

trunk/LMDZ.MARS/util/solzenangle.def

@@ -18 +18 @@
 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
-
+6) reference time 
+    - for a diagfi/concat : 
+        - write "n" to keep the date stored in the file (recommended)      
+        - if needed:  sol value at the beginning of the run, wrt Ls=0
+    - for a stats : Ls value at the middle of the run (degree)
 USE :
 solzenangle.e < solzenangle.def