Changeset 2434 for trunk/LMDZ.MARS/util/solzenangle.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/solzenangle.F90 (modified) (60 diffs)

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