Changeset 2434 for trunk


Ignore:
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
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r2432 r2434  
    32063206  concentration in mol.cm-3 etc..). Set missing value to 1.E20
    32073207  lslin.F90: fix various problems in the output.
     3208 
     3209== 17/11/2020 == FF + AB
     3210Update utilities :
     3211concat.F90
     3212- rewrite and simplify the handling of time and offset so that any file
     3213  can be concatenated, including files from different years or stats file
     3214- use modulo to shift starting sols in concat.nc to values between 0 and 669
     3215- add a "adls" option in addition to "sol" or "ls" to add a Ls 1D variable
     3216   while keeping "sol" as the Time axis
     3217- add conservation of altitude attributes (long_name,units,positive)
     3218- enable absence of both hybrid (aps,bps) and sigma levels
     3219
     3220solzenangle.F90
     3221- improve calculation for 1 sol file (stats) to use all local time data
     3222- read the starting sol in the file instead of asking it to the user
     3223  (except for stats file) ; keep the possibility for user to change it
     3224- ask mean Ls value for stats file instead of sol number
     3225- fix crash when using "all variable" option (ticket #66)
     3226- fix bug on aps,bps by using GCM_layers dimension instead of altitude
     3227
     3228localtime.F90
     3229- fix crash when using "all variable" option (ticket #66)
     3230- fix bug on aps,bps by using GCM_layers dimension instead of altitude
     3231
     3232For all 3 utilities :
     3233- handle all ndim cases for the variables (ticket #52)
     3234- change "title" attribute into "long_name" by default (ticket #48)
     3235- extend size of long_name character string (ticket #57)
     3236- remove the use of #ifdef NC_DOUBLE (ticket #67)

  • trunk/LMDZ.MARS/util/concatnc.F90

    r2395 r2434  
    3434character (len=50), dimension(:), allocatable :: var
    3535! var(): name(s) of variable(s) that will be concatenated
    36 character (len=50) :: tmpvar,tmpfile,title,units
     36character (len=100) :: tmpvar,tmpfile,long_name,units
    3737! tmpvar(): used to temporarily store a variable name
    3838! tmpfile(): used to temporarily store a file name
    39 ! title(): [netcdf] title attribute
     39! long_name(): [netcdf] long_name attribute
    4040! units(): [netcdf] units attribute
    4141character (len=100) :: filename,vartmp
    4242! filename(): output file name
    4343! vartmp(): temporary variable name (read from netcdf input file)
    44 !character (len=1) :: ccopy
    45 ! ccpy: 'y' or 'n' answer
     44character (len=50) :: altlong_name,altunits,altpositive
     45! altlong_name(): [netcdf] altitude "long_name" attribute
     46! altunits(): [netcdf] altitude "units" attribute
     47! altpositive(): [netcdf] altitude "positive" attribute
    4648character (len=4) :: axis
    47 ! axis: "ls" or "sol"
     49! axis:  "sol" or "ls" or "adls"
    4850integer :: nid,ierr,miss
    4951! nid: [netcdf] file ID #
    … …  
    115117real :: ctlsol
    116118! ctlsol: starting sol value of the file, stored in the variable ctl()
    117 real :: memotime
    118 ! memotime: (cumulative) time value, in martian days (sols)
     119real :: previous_last_output_time, output_time
     120! previous_last_output_time: last time stored in concat.nc after treating a file
     121logical :: firstsol
     122! firstsol:  true if the first file's starting day must be used in the output
     123real :: startsol
     124! startsol: sol of the firstday as requested by the users (used if firstsol="y")
    119125real :: starttimeoffset
    120 ! starttimeoffset: offset given by the user for the output time axis
    121 character (len=1) :: firstsol
    122 ! firstsol: =[y/n], to know whether the user offset or the first file's starting day must be used in the output
     126! starttimeoffset: offset given by the user for the output time axis (0 if ! firstsol="n")
    123127real :: missing
    124128!PARAMETER(missing=1E+20)
    … …  
    126130real, dimension(2) :: valid_range
    127131! valid_range(2): [netcdf] interval in which a value is considered valid
     132real :: year_day = 669.
    128133
    129134!==============================================================================
    … …  
    154159!==============================================================================
    155160
    156 write(*,*)
    157 write(*,*) "Starting day to be stored in the output file time axis?"
    158 write(*,*) "(e.g.: 100 if you want the output file to start at time=100 sols)"
    159 write(*,*) "Your answer will bypass any starting day value stored in the 1st"
    160 write(*,*) "input file. Conversely, just type return if you want this 1st file's"
    161 write(*,*) "value to be kept for the output file starting day."
    162 read(*,*,iostat=ierr) starttimeoffset
    163 if (ierr.ne.0) then
     161 write(*,*) "Starting time in the concat file ?"
     162 write(*,*) "   If you really wish to change it, write it now (not recommanded!)"
     163 write(*,*) "   If you do not want to change it,  just write 'n' (recommanded)"
     164
     165read(*,*,iostat=ierr) startsol
     166
     167firstsol= (ierr.eq.0)
     168if (firstsol) then
    164169   ! if nothing or a character that is not a float is read
    165170   write(*,*) "1st input file's starting day will serve as starting day for the outfile."
    166    firstsol='n'
    167171else ! if the user gave a number
    168172   write(*,*) "Your input day will serve as starting day for the outfile."
    169    firstsol='y'
    170173endif
    171174
    … …  
    190193
    191194write(*,*)
    192 write(*,*) "Warning: to read the result with grads, choose sol"
    193 write(*,*) "Warning: use ferret to read the non linear scale ls"
    194 write(*,*) "Which time axis should be given in the output? (sol/ls)"
     195write(*,*) "Which time axis should be given in the output? (sol/ls/adls)"
     196write(*,*) "        (Warning: to read the result with grads, choose sol)"
     197write(*,*) "  To keep sol as the time axis, and just add Ls as a variable, type adls  "
    195198read(*,*) axis
    196199! loop as long as axis is neither "sol" nor "ls"
    197 do while ((axis/="sol").AND.(axis/="ls"))
     200do while ((axis/="sol").AND.(axis/="ls").AND.(axis/="adls"))
    198201   read(*,*) axis
    199202enddo
    … …  
    255258if (tmpvar=="all") then
    256259   if (axis=="ls") then
    257 !      write(*,*) "Do you want to keep the original file? (y/n)"
    258 !      read(*,*) ccopy
    259 !      if ((ccopy=="n").or.(ccopy=="N")) then
    260 !         do i=1,nbfile
    261 !            ierr=NF_CLOSE(nid)
    262 !            ierr = NF_OPEN(file(1),NF_WRITE,nid)
    263 !            call change_time_axis(nid,ierr)
    264 !            ierr=NF_CLOSE(nid)
    265 !            STOP ""
    266 !         enddo
    267 !      else
    268260         nbvar=nbvarfile-Nnotconcat
    269261         do j=Nnotconcat+1,nbvarfile
    270262            ierr=nf_inq_varname(nid,j,var(j-Nnotconcat))
    271263         enddo
    272 !      endif ! of if ((ccopy=="n").or.(ccopy=="N"))
    273264   endif ! of if (axis=="ls")
    274265! Variables names from the file are catched
    … …  
    394385      endif
    395386   endif
    396 !  write(*,*) "controle: ",controle
    397387
    398388   if (ctllen .ne. 0) then ! variable controle
    399389      allocate(ctl(ctllen))
    400 #ifdef NC_DOUBLE
    401       ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
    402 #else
    403390      ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
    404 #endif
    405391      ctlsol = ctl(4)
    406392   endif
    … …  
    425411   allocate(time(timelen))
    426412
    427 #ifdef NC_DOUBLE
    428    ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
    429 #else
    430413   ierr = NF_GET_VAR_REAL(nid,timevar,time)
    431 #endif
    432 
    433414
    434415!==============================================================================
    … …  
    437418!==============================================================================
    438419
    439    if (i==1) then ! First call; initialize/allocate
     420   if (i==1) then ! First File ; initialize/allocate
    440421      memolatlen=latlen
    441422      memolonlen=lonlen
    … …  
    445426      allocate(lon(lonlen))
    446427      allocate(alt(altlen))
    447 #ifdef NC_DOUBLE
    448       ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat)
    449       ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon)
    450       ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt)
    451 #else
    452428      ierr = NF_GET_VAR_REAL(nid,latvar,lat)
    453429      ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
    454430      ierr = NF_GET_VAR_REAL(nid,altvar,alt)
    455 #endif
     431!   Get altitude attributes to handle files with any altitude type
     432      ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
     433      if (ierr.ne.nf_noerr) then
     434      ! if no attribute "long_name", try "title"
     435        ierr=nf_get_att_text(nid,altvar,"title",altlong_name)
     436      endif
     437      ierr=nf_get_att_text(nid,altvar,'units',altunits)
     438      ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
    456439      write(*,*)
     440     
    457441      if (ctllen .ne. 0) then
    458          if (int(starttimeoffset) .ne. (int(ctlsol)+int(time(1))) ) then
    459            write(*,*) "WARNING: Starting day of the first file is not ",&
    460                        int(starttimeoffset)," but ",int(ctlsol)+int(time(1)),"!"
     442       
     443         if (firstsol) then
     444           starttimeoffset = float(int(ctlsol)+int(time(1) + ctl(27))) - startsol
     445         else ! if firstsol.eq..false.
     446!          if there is no user chosen first sol, nevertheless use
     447!          starttimeoffset to lower the sol number (same exact season)
     448           startsol = modulo(int(int(ctlsol)+int(time(1) + ctl(27))),int(year_day))
     449           starttimeoffset = float(int(ctlsol)+int(time(1) + ctl(27))) - startsol
    461450         endif
    462        
    463          if (firstsol.eq.'y') then
    464            starttimeoffset = float(int(ctlsol)+int(time(1))) + ctl(27) - starttimeoffset
    465          else ! if firstsol.eq.'n'
    466            starttimeoffset = 0
    467          endif
    468 
    469          memotime=float(int(ctlsol)+int(time(1))) + ctl(27)
     451
     452!        artificially estimating "previous_last_output_time" for first file read:
     453         previous_last_output_time= startsol
     454
    470455         ctl(4) = 0.  ! values written in the output
    471456         ctl(27) = 0. ! file controle variable (resp. day_ini and hour_ini)
    472457         
    473458      else ! if no variable "controle"
    474          if (int(starttimeoffset) .ne. (int(time(1))) ) then
    475            write(*,*) "WARNING: Starting day of the first file is not ",&
    476                        int(starttimeoffset)," but ",int(time(1)),"!"
    477          endif
    478 
    479          if (firstsol.eq.'y') then
    480            starttimeoffset = int(time(1)) - starttimeoffset
    481          else ! if firstsol.eq.'n'
     459         if (firstsol) then
     460           starttimeoffset = int(time(1)) - startsol
     461         else ! if firstsol.eq.false
    482462           starttimeoffset = 0
    483463         endif
    484464
    485          memotime=float(int(time(1)))
     465!        artificially estimating "previous_last_output_time" for first file read:
     466         previous_last_output_time=time(1) -(time(2) - time(1)) - starttimeoffset
    486467      endif
     468      write(*,*) "Original starting day of 1st input file", previous_last_output_time + starttimeoffset
     469      write(*,*) "Starting day of the concat file : ", previous_last_output_time
     470
    487471       
    488472   ! Initialize output file's lat,lon,alt and time dimensions
    489473      write(*,*)
    490       call initiate(filename,lat,lon,alt,ctl,GCM_layers,nout,&
    491        latdimout,londimout,altdimout,timedimout,&
     474      call initiate(filename,lat,lon,alt,ctl,GCM_layers,&
     475       altlong_name,altunits,altpositive,&
     476       nout,latdimout,londimout,altdimout,timedimout,&
    492477       layerdimout,interlayerdimout,timevarout,axis)
    493478   ! Initialize output file's aps,bps,ap,bp and phisinit variables
    494      call init2(nid,lonlen,latlen,altlen,GCM_layers,&
     479      call init2(nid,lonlen,latlen,altlen,GCM_layers,&
    495480                nout,londimout,latdimout,altdimout,&
    496481                layerdimout,interlayerdimout)
    497482   
    498    else ! Not a first call,
     483   else ! Not first file to concatenate
    499484   ! Check that latitude,longitude and altitude of current input file
    500485   ! are identical to those of the output file
    … …  
    521506
    522507   rep=0
    523    write(*,*)
    524    write(*,'(a3,1x,f6.1)') "Sol",memotime-starttimeoffset
    525 !   if (ctllen.ne.0) write(*,'(a30,1x,f6.1)') "In the current file's ctl: Sol",ctlsol
    526 
    527    ! Add (memotime)/(ctlsol) offset and write "concatenated" time values to output file
    528    if (ctllen.ne.0) then
    529      do k=reptime+1,reptime+timelen
     508
     509   do k=reptime+1,reptime+timelen
    530510       rep=rep+1
    531        if ((time(rep)+ctlsol).le.memotime) then
    532          write(*,*) "ERROR : the time intervals of the files are not dissociated"
    533          stop ""
    534        endif
    535 #ifdef NC_DOUBLE
    536        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,k,1,ctlsol+time(rep)-starttimeoffset)
    537 #else
    538        ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,ctlsol+time(rep)-starttimeoffset)
    539 #endif
    540        !
    541      enddo
    542 
    543    ! Compute new time offset (for further concatenations)
    544    memotime=ctlsol+time(timelen) ! the actual last day of the current file
    545                                  ! becomes the reference time for next file
    546                                  ! leaving potential "time holes" in the output time axis
    547                                  ! at the junction of two files' time axes
    548    
    549    else ! if no variable "controle"
    550    
    551      do k=reptime+1,reptime+timelen
    552         rep=rep+1
    553 #ifdef NC_DOUBLE
    554         ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,k,1,memotime+time(rep)-int(time(1))-starttimeoffset)
    555 #else
    556         ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,memotime+time(rep)-int(time(1))-starttimeoffset)
    557 #endif
    558      enddo
    559      
    560    ! Compute new time offset (for further concatenations)
    561    memotime=memotime+time(timelen)-time(1) ! files are concatenated one after another
    562                                            ! even if their time axes normally don't
    563                                            ! directly follow each other
    564    endif ! if ctllen.ne.0
     511       if (ctllen.ne.0) then
     512          output_time=time(rep)+ctlsol-starttimeoffset
     513!         correction if the file to concatenate seems "before" previous file
     514!         (for instance to concatenate diagfi from the previous year at the enf of a year)
     515          do while (output_time.lt.previous_last_output_time) 
     516                   output_time = output_time + year_day
     517          end do
     518       else ! if no control, just add timestep after timestep
     519          output_time= previous_last_output_time + (time(2)-time(1)) &
     520         + (time(rep)-time(1))
     521       end if
     522       if (rep.eq.1) write(*,*) "Sol", int(output_time) 
     523
     524       ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,output_time)
     525   end do
     526!  use the last output_time value to update memotime   
     527   previous_last_output_time = output_time
     528
    565529
    566530!==============================================================================
    … …  
    577541      ierr=nf_inq_varid(nid,var(j),varid)
    578542      if (ierr.NE.NF_NOERR) then
    579          write(*,*) 'ERROR: Field <',var(j),'> not found in file'//file(i)
     543         write(*,*) 'ERROR: Field <',var(j),'> not found in file '//file(i)
    580544         stop ""
    581545      endif
    582546      ierr=nf_inq_varndims(nid,varid,ndim)
    583547
     548      ! Check that it is a 1D, 3D or 4D variable
     549      if ((ndim.ne.1).and.(ndim.ne.3).and.(ndim.ne.4)) then
     550        write(*,*) "Error:",trim(var(j))," is not a time-dependent variable,"
     551        write(*,*) "so it can't be concatenated."
     552        write(*,*) "We'll skip that variable..."
     553        CYCLE ! go directly to the next variable
     554      endif
     555     
    584556!==============================================================================
    585557! 2.5.2 Prepare things in order to read/write the variable
    … …  
    644616      if (i==1) then ! First call: write some definitions to output file
    645617         units="                                                    "
    646          title="                                                    "
    647          ierr=nf_get_att_text(nid,varid,"title",title)
     618         long_name="                                                    "
     619         ierr=nf_get_att_text(nid,varid,"long_name",long_name)
     620         if (ierr.ne.nf_noerr) then
     621         ! if no attribute "long_name", try "title"
     622           ierr=nf_get_att_text(nid,varid,"title",long_name)
     623         endif
    648624         ierr=nf_get_att_text(nid,varid,"units",units)
    649          call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr)
     625         call def_var(nout,var(j),long_name,units,ndim,dim,varidout,ierr)
    650626      else
    651627         ierr=NF_INQ_VARID(nout,var(j),varidout)
    … …  
    656632!==============================================================================
    657633
    658 #ifdef NC_DOUBLE
    659       ierr = NF_GET_VAR_DOUBLE(nid,varid,var3d)
    660       ierr= NF_PUT_VARA_DOUBLE(nout,varidout,corner,edges,var3d)
    661       miss=NF_GET_ATT_DOUBLE(nid,varid,"missing_value",missing)
    662       miss=NF_GET_ATT_DOUBLE(nid,varid,"valid_range",valid_range)
    663 #else
    664634      ierr = NF_GET_VAR_REAL(nid,varid,var3d)
    665635      ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3d)
    666636      miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
    667637      miss=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
    668 #endif
    669638
    670639      if (ierr.ne.NF_NOERR) then
    … …  
    700669!==============================================================================
    701670
    702 if (axis=="ls") then
    703    write(*,*);write(*,*) "Converting the time axis from sol to Ls..."
    704    call change_time_axis(nout,ierr)
     671if ((axis=="ls").or.(axis=="adls")) then
     672   call change_time_axis(nout,ierr,axis)
    705673endif
    706674
    … …  
    712680
    713681!******************************************************************************
    714 subroutine initiate (filename,lat,lon,alt,ctl,GCM_layers,nout,&
    715          latdimout,londimout,altdimout,timedimout,&
    716          layerdimout,interlayerdimout,timevarout,axis)
     682subroutine initiate (filename,lat,lon,alt,ctl,GCM_layers,&
     683                     altlong_name,altunits,altpositive,&
     684                     nout,latdimout,londimout,altdimout,timedimout,&
     685                     layerdimout,interlayerdimout,timevarout,axis)
    717686!==============================================================================
    718687! Purpose:
    … …  
    741710! ctl(): controle
    742711integer,intent(in) :: GCM_layers ! number of GCM layers
     712character (len=*), intent(in) :: altlong_name
     713! altlong_name(): [netcdf] altitude "long_name" attribute
     714character (len=*), intent(in) :: altunits
     715! altunits(): [netcdf] altitude "units" attribute
     716character (len=*), intent(in) :: altpositive
     717! altpositive(): [netcdf] altitude "positive" attribute
    743718integer, intent(out):: nout
    744719! nout: [netcdf] file ID
    … …  
    787762ierr = NF_DEF_DIM(nout, "longitude", size(lon), londimout)
    788763ierr = NF_DEF_DIM(nout, "altitude", size(alt), altdimout)
    789 if (size(ctl).ne.0) ierr = NF_DEF_DIM(nout, "index", size(ctl), ctldimout)
     764if (ctllen.ne.0) ierr = NF_DEF_DIM(nout, "index", size(ctl), ctldimout)
    790765ierr = NF_DEF_DIM(nout, "Time", NF_UNLIMITED, timedimout)
    791766ierr = NF_DEF_DIM(nout, "GCM_layers", GCM_layers, layerdimout)
    … …  
    798773! 3. Write "Time" and its attributes
    799774!==============================================================================
    800 if (axis=="sol") then
     775if (axis=="ls") then
     776  call def_var(nout,"Time","Solar longitude","degree",1,&
     777             (/timedimout/),timevarout,ierr)
     778else
    801779  call def_var(nout,"Time","Time","days since 0000-00-0 00:00:00",1,&
    802780             (/timedimout/),timevarout,ierr)
    803 else ! Ls
    804   call def_var(nout,"Time","Solar longitude","days since 0000-00-0 00:00:00",1,&
    805              (/timedimout/),timevarout,ierr)
    806781endif
    807782!==============================================================================
    … …  
    812787             (/latdimout/),nvarid,ierr)
    813788
    814 #ifdef NC_DOUBLE
    815 ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lat)
    816 #else
    817789ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
    818 #endif
    819790
    820791!==============================================================================
    … …  
    825796             (/londimout/),nvarid,ierr)
    826797
    827 #ifdef NC_DOUBLE
    828 ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lon)
    829 #else
    830798ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
    831 #endif
    832799
    833800!==============================================================================
    … …  
    838805ierr = NF_REDEF (nout)
    839806
    840 #ifdef NC_DOUBLE
    841 ierr = NF_DEF_VAR (nout,"altitude",NF_DOUBLE,1,altdimout,nvarid)
    842 #else
    843807ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,altdimout,nvarid)
    844 #endif
    845 
    846 ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",8,"altitude")
    847 ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',2,"km")
    848 ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',2,"up")
     808
     809ierr = NF_PUT_ATT_TEXT (nout,nvarid,'long_name',len_trim(adjustl(altlong_name)),adjustl(altlong_name))
     810ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',len_trim(adjustl(altunits)),adjustl(altunits))
     811ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',len_trim(adjustl(altpositive)),adjustl(altpositive))
    849812
    850813! End netcdf define mode
    851814ierr = NF_ENDDEF(nout)
    852815
    853 #ifdef NC_DOUBLE
    854 ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,alt)
    855 #else
    856816ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
    857 #endif
    858817
    859818!==============================================================================
    … …  
    861820!==============================================================================
    862821
    863 if (size(ctl).ne.0) then
     822if (ctllen.ne.0) then
    864823   ! Switch to netcdf define mode
    865824   ierr = NF_REDEF (nout)
    866825
    867 #ifdef NC_DOUBLE
    868    ierr = NF_DEF_VAR (nout,"controle",NF_DOUBLE,1,ctldimout,nvarid)
    869 #else
    870826   ierr = NF_DEF_VAR (nout,"controle",NF_FLOAT,1,ctldimout,nvarid)
    871 #endif
    872 
    873    ierr = NF_PUT_ATT_TEXT (nout,nvarid,"title",18,"Control parameters")
     827
     828   ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",18,"Control parameters")
    874829
    875830   ! End netcdf define mode
    876831   ierr = NF_ENDDEF(nout)
    877832
    878 #ifdef NC_DOUBLE
    879    ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,ctl)
    880 #else
    881833   ierr = NF_PUT_VAR_REAL (nout,nvarid,ctl)
    882 #endif
    883834endif
    884835
    … …  
    890841!==============================================================================
    891842! Purpose:
    892 ! Copy ap() , bp(), aps(), bps(), aire() and phisinit()
    893 ! from input file to outpout file
     843! Copy ap(), bp(), aps(), bps(), aire() and phisinit()
     844! from input file to output file
    894845!==============================================================================
    895846! Remarks:
    … …  
    907858integer, intent(in) :: lonlen ! # of grid points along longitude
    908859integer, intent(in) :: latlen ! # of grid points along latitude
    909 integer, intent(in) :: altlen ! # of grid points along latitude
     860integer, intent(in) :: altlen ! # of grid points along altitude
    910861integer, intent(in) :: GCM_layers ! # of GCM atmospheric layers
    911862integer, intent(in) :: outfid ! NetCDF output file ID
    … …  
    931882logical :: hybrid ! are "aps" and "bps" available ?
    932883logical :: apbp ! are "ap" and "bp" available ?
     884logical :: sig ! is "sigma" available ?
    933885
    934886!==============================================================================
    … …  
    976928  write(*,*) "init2: failed to allocate ap!"
    977929  stop
     930endif
     931ierr=NF_INQ_VARID(infid,"ap",tmpvarid)
     932if (ierr.ne.NF_NOERR) then
     933  write(*,*) "Ooops. Failed to get ap ID. OK."
     934  apbp=.false.
    978935else
    979   ierr=NF_INQ_VARID(infid,"ap",tmpvarid)
    980   if (ierr.ne.NF_NOERR) then
    981     write(*,*) "Ooops. Failed to get ap ID. OK."
    982     apbp=.false.
    983   else
    984     ierr=NF_GET_VAR_REAL(infid,tmpvarid,ap)
    985     apbp=.true.
    986     if (ierr.ne.NF_NOERR) then
    987       stop "Error: Failed reading ap"
    988     endif
     936  ierr=NF_GET_VAR_REAL(infid,tmpvarid,ap)
     937  apbp=.true.
     938  if (ierr.ne.NF_NOERR) then
     939    stop "Error: Failed reading ap"
    989940  endif
    990941endif
    … …  
    995946  write(*,*) "init2: failed to allocate bp!"
    996947  stop
     948endif
     949ierr=NF_INQ_VARID(infid,"bp",tmpvarid)
     950if (ierr.ne.NF_NOERR) then
     951  write(*,*) "Ooops. Failed to get bp ID. OK."
     952  apbp=.false.
    997953else
    998   ierr=NF_INQ_VARID(infid,"bp",tmpvarid)
    999   if (ierr.ne.NF_NOERR) then
    1000     write(*,*) "Ooops. Failed to get bp ID. OK."
    1001     apbp=.false.
    1002   else
    1003     ierr=NF_GET_VAR_REAL(infid,tmpvarid,bp)
    1004     apbp=.true.
    1005     if (ierr.ne.NF_NOERR) then
    1006       stop "Error: Failed reading bp"
    1007     endif
     954  ierr=NF_GET_VAR_REAL(infid,tmpvarid,bp)
     955  apbp=.true.
     956  if (ierr.ne.NF_NOERR) then
     957    stop "Error: Failed reading bp"
    1008958  endif
    1009959endif
    … …  
    1017967  endif
    1018968  ierr=NF_INQ_VARID(infid,"sigma",tmpvarid)
    1019   ierr=NF_GET_VAR_REAL(infid,tmpvarid,sigma)
    1020   if (ierr.ne.NF_NOERR) then
    1021     stop "init2 Error: Failed reading sigma"
     969  if (ierr.ne.NF_NOERR) then
     970    write(*,*) "Ooops. Failed to get sigma ID. OK."
     971    sig=.false.
     972  else
     973    ierr=NF_GET_VAR_REAL(infid,tmpvarid,sigma)
     974    sig=.true.
     975    if (ierr.ne.NF_NOERR) then
     976      stop "init2 Error: Failed reading sigma"
     977    endif
    1022978  endif
    1023979endif ! of if (.not.hybrid)
    … …  
    10751031
    10761032! write aps
    1077 #ifdef NC_DOUBLE
    1078   ierr=NF_PUT_VAR_DOUBLE(outfid,apsid,aps)
    1079 #else
    10801033  ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
    1081 #endif
    10821034  if (ierr.ne.NF_NOERR) then
    10831035    stop "init2 Error: Failed to write aps"
    … …  
    10921044
    10931045! write bps
    1094 #ifdef NC_DOUBLE
    1095   ierr=NF_PUT_VAR_DOUBLE(outfid,bpsid,bps)
    1096 #else
    10971046  ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
    1098 #endif
    10991047  if (ierr.ne.NF_NOERR) then
    11001048    stop "init2 Error: Failed to write bps"
    … …  
    11101058
    11111059! write ap
    1112 #ifdef NC_DOUBLE
    1113     ierr=NF_PUT_VAR_DOUBLE(outfid,apid,ap)
    1114 #else
    11151060    ierr=NF_PUT_VAR_REAL(outfid,apid,ap)
    1116 #endif
    11171061    if (ierr.ne.NF_NOERR) then
    11181062      stop "Error: Failed to write ap"
    … …  
    11271071
    11281072!   write bp
    1129 #ifdef NC_DOUBLE
    1130     ierr=NF_PUT_VAR_DOUBLE(outfid,bpid,bp)
    1131 #else
    11321073    ierr=NF_PUT_VAR_REAL(outfid,bpid,bp)
    1133 #endif
    11341074    if (ierr.ne.NF_NOERR) then
    11351075      stop "Error: Failed to write bp"
    … …  
    11371077  endif ! of if (apbp)
    11381078
    1139 else
     1079else if (sig) then
    11401080! define sigma
    11411081  call def_var(nout,"sigma","sigma at midlayers"," ",1,&
    … …  
    11451085  endif
    11461086! write sigma
    1147 #ifdef NC_DOUBLE
    1148   ierr=NF_PUT_VAR_DOUBLE(outfid,sigmaid,sigma)
    1149 #else
    11501087  ierr=NF_PUT_VAR_REAL(outfid,sigmaid,sigma)
    1151 #endif
    11521088  if (ierr.ne.NF_NOERR) then
    11531089    stop "init2 Error: Failed to write sigma"
    … …  
    11681104 
    11691105  ! write aire
    1170 #ifdef NC_DOUBLE
    1171   ierr=NF_PUT_VAR_DOUBLE(outfid,tmpvarid,aire)
    1172 #else
    11731106  ierr=NF_PUT_VAR_REAL(outfid,tmpvarid,aire)
    1174 #endif
    11751107  if (ierr.ne.NF_NOERR) then
    11761108    stop "init2 Error: Failed to write aire"
    … …  
    11881120
    11891121  ! write phisinit
    1190 #ifdef NC_DOUBLE
    1191   ierr=NF_PUT_VAR_DOUBLE(outfid,phisinitid,phisinit)
    1192 #else
    11931122  ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
    1194 #endif
    11951123  if (ierr.ne.NF_NOERR) then
    11961124    stop "init2 Error: Failed to write phisinit"
    … …  
    12111139end subroutine init2
    12121140!******************************************************************************
    1213 subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
     1141subroutine def_var(nid,name,long_name,units,nbdim,dim,nvarid,ierr)
    12141142!==============================================================================
    12151143! Purpose: Write a variable (i.e: add a variable to a dataset)
    1216 ! called "name"; along with its attributes "title", "units"...
     1144! called "name"; along with its attributes "long_name", "units"...
    12171145! to a file (following the NetCDF format)
    12181146!==============================================================================
    … …  
    12321160character (len=*), intent(in) :: name
    12331161! name(): [netcdf] variable's name
    1234 character (len=*), intent(in) :: title
    1235 ! title(): [netcdf] variable's "title" attribute
     1162character (len=*), intent(in) :: long_name
     1163! long_name(): [netcdf] variable's "long_name" attribute
    12361164character (len=*), intent(in) :: units
    12371165! unit(): [netcdf] variable's "units" attribute
    … …  
    12491177
    12501178! Insert the definition of the variable
    1251 #ifdef NC_DOUBLE
    1252 ierr=NF_DEF_VAR(nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
    1253 #else
    12541179ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
    1255 #endif
    12561180
    12571181! Write the attributes
    1258 ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",len_trim(adjustl(title)),adjustl(title))
     1182ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name",len_trim(adjustl(long_name)),adjustl(long_name))
    12591183ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
    12601184
    … …  
    12641188end subroutine def_var
    12651189!******************************************************************************
    1266 subroutine change_time_axis(nid,ierr)
     1190subroutine change_time_axis(nid,ierr,axis)
    12671191!==============================================================================
    12681192! Purpose:
    … …  
    12851209integer, intent(out) :: ierr
    12861210! ierr: [netcdf]  return error code
     1211character (len=4),intent(in) :: axis
     1212! axis: "ls" or "sol"
    12871213
    12881214!==============================================================================
    … …  
    13151241ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
    13161242allocate(time(timelen),ls(timelen))
    1317 #ifdef NC_DOUBLE
    1318 ierr = NF_GET_VAR_DOUBLE(nid,nvarid,time)
    1319 #else
    13201243ierr = NF_GET_VAR_REAL(nid,nvarid,time)
    1321 #endif
    13221244if (ierr.ne.NF_NOERR) then
    13231245  write(*,*) "ERROR in change_time_axis: Failed to load Time"
    … …  
    13531275!==============================================================================
    13541276
    1355 #ifdef NC_DOUBLE
    1356 ierr = NF_PUT_VAR_DOUBLE(nid,nvarid,ls)
    1357 #else
     1277if (axis.eq."ls") then
     1278   write(*,*) "Converting the time axis from sol to Ls..."
     1279else ! if (axis.eq."adls")
     1280   write(*,*) "Adding Ls as a 1D time variable"
     1281   call def_var(nout,"Ls","Solar Longitude","degree",1, (/timedimout/),nvarid,ierr)
     1282end if
     1283
    13581284ierr = NF_PUT_VAR_REAL(nid,nvarid,ls)
    1359 #endif
    13601285if (ierr.ne.NF_NOERR) then
    1361   write(*,*) "ERROR in change_time_axis: Failed to write Ls"
    1362   stop
     1286   write(*,*) "ERROR in change_time_axis: Failed to write Ls"
     1287   stop
    13631288endif
    13641289
    … …  
    15021427
    15031428! Write valid_range() attribute
    1504 #ifdef NC_DOUBLE
    1505 ierr=NF_PUT_ATT_DOUBLE(nout,nvarid,'valid_range',NF_DOUBLE,2,valid_range)
    1506 #else
    15071429ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
    1508 #endif
    15091430
    15101431if (ierr.ne.NF_NOERR) then
    … …  
    15171438
    15181439! Write "missing_value" attribute
    1519 #ifdef NC_DOUBLE
    1520 ierr= NF_PUT_ATT_DOUBLE(nout,nvarid,'missing_value',NF_DOUBLE,1,missing)
    1521 #else
    15221440ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
    1523 #endif
    15241441
    15251442if (ierr.NE.NF_NOERR) then

  • trunk/LMDZ.MARS/util/concatnc.def

    r2313 r2434  
    33diagfi12.nc       
    44
    5 515              ! Starting sol to use in the output file?
    6 sol              ! output timescale ("sol" or "ls")
     5n              ! Starting sol of the run stored in the first input file?
     6sol              ! output timescale ("sol" or "Ls" or "adls")
    77tsurf           
    88ps
    99
    1010
     11
    1112-----------------------------------------------------------------------
    12 ABOVE is the list of inputs to be fed to "concatnc.e" if you don't want
     13ABOVE is the list of inputs to be fed to "concatnc.e" if you don t want
    1314to reply directly to the program:
    1415
    15161) List of N files to be read  (N lines)
    16172) blank line at the end
    17 3) Starting sol of the output file?
    18 (or a blank if you want to use the first file's
    19  starting sol = time(1)+controle(4) )
    20 4) Output timescale ("sol" or "ls")
    21 5) List of X variables to be kept (X lines) or 'all'
    22 6) blank line at the end
     183) Starting sol of the run stored in the first input file? juste write "n" if you
     19   wish to keep the date of the input files (recommended)
     204) List of X variables to be kept (X lines) or 'all'
     215) blank line at the end
    2322
    2423The use is :

  • trunk/LMDZ.MARS/util/localtime.F90

    r2280 r2434  
    1414character (len=50)  file
    1515! file(): input file(s) names(s)
    16 character (len=30), dimension(15) :: notconcat
    17 ! notconcat(): names of the (15) variables that won't be concatenated
     16character (len=30), dimension(16) :: notprocessed
     17! notprocessed(): names of the (16) variables that won't be processed
    1818character (len=50), dimension(:), allocatable :: var
    19 ! var(): name(s) of variable(s) that will be concatenated
    20 character (len=50) :: tmpvar,title,units
     19! var(): name(s) of variable(s) that will be processed
     20character (len=100) :: tmpvar,long_name,units
    2121! tmpvar(): used to temporarily store a variable name
    22 ! title(): [netcdf] title attribute
     22! long_name(): [netcdf] long_name attribute
    2323! units(): [netcdf] units attribute
    2424character (len=100) :: filename,vartmp
    … …  
    5858! ctldim: [netcdf] "controle" dim ID
    5959! timedim: [netcdf] "timedim" dim ID
     60integer :: gcmlayerdim ! NetCDF dimension ID for # of layers in GCM
    6061integer :: latvar,lonvar,altvar,ctlvar,timevar
    6162! latvar: [netcdf] ID of "latitude" variable
    … …  
    6869! latlen: # of elements of lat() array
    6970! lonlen: # of elements of lon() array
    70 ! altvar: # of elements of alt() array
    71 ! ctlvar: # of elements of ctl() array
     71! altlen: # of elements of alt() array
     72! ctllen: # of elements of ctl() array
    7273! timelen: # of elemnets of time() array
    7374! timelen_tot:# =timelen or timelen+1 (if 1 more time to interpolate needed)
    7475! timelen_lt: # of elemnets of time() array in output
     76integer :: nhour,nsol
     77integer :: GCM_layers ! number of GCM atmospheric layers (may not be
     78! same as altlen if file is an output of zrecast)
    7579integer :: nout,latdimout,londimout,altdimout,timedimout,timevarout
    76 integer :: nhour,nsol
    7780! nout: [netcdf] output file ID
    7881! latdimout: [netcdf] output latitude (dimension) ID
    … …  
    8184! timedimout: [netcdf] output time (dimension) ID
    8285! timevarout: [netcdf] ID of output "Time" variable
     86integer :: layerdimout ! NetCDF dimension ID for # of layers in GCM
    8387integer :: varidout
    8488! varidout: [netcdf] variable ID # (of a variable to write to the output file)
    85 integer :: Nnotconcat,var_ok
    86 ! Nnotconcat: # of (leading)variables that won't be concatenated
     89integer :: Nnotprocessed,var_ok
     90! Nnotprocessed: # of (leading)variables that won't be concatenated
    8791! var_ok: flag (0 or 1)
    8892integer, dimension(4) :: corner,edges,dim
    … …  
    143147!==============================================================================
    144148
    145 notconcat(1)='Time'
    146 notconcat(2)='controle'
    147 notconcat(3)='rlonu'
    148 notconcat(4)='latitude'
    149 notconcat(5)='longitude'
    150 notconcat(6)='altitude'
    151 notconcat(7)='rlatv'
    152 notconcat(8)='aps'
    153 notconcat(9)='bps'
    154 notconcat(10)='ap'
    155 notconcat(11)='bp'
    156 notconcat(12)='cu'
    157 notconcat(13)='cv'
    158 notconcat(14)='aire'
    159 notconcat(15)='phisinit'
     149notprocessed(1)='Time'
     150notprocessed(2)='controle'
     151notprocessed(3)='rlonu'
     152notprocessed(4)='latitude'
     153notprocessed(5)='longitude'
     154notprocessed(6)='altitude'
     155notprocessed(7)='rlatv'
     156notprocessed(8)='aps'
     157notprocessed(9)='bps'
     158notprocessed(10)='ap'
     159notprocessed(11)='bp'
     160notprocessed(12)='soildepth'
     161notprocessed(13)='cu'
     162notprocessed(14)='cv'
     163notprocessed(15)='aire'
     164notprocessed(16)='phisinit'
     165
    160166
    161167!==============================================================================
    … …  
    163169!==============================================================================
    164170write(*,*)
    165    Nnotconcat=0
     171   Nnotprocessed=0
    166172do i=1,nbvarfile
    167173   ierr=NF_INQ_VARNAME(nid,i,vartmp)
    168174   ! vartmp now contains the "name" of variable of ID # i
    169175   var_ok=0
    170    do inter=1,15
    171       if (vartmp.eq.notconcat(inter)) then
     176   do inter=1,16
     177      if (vartmp.eq.notprocessed(inter)) then
    172178         var_ok=1
    173          Nnotconcat=Nnotconcat+1
     179         Nnotprocessed=Nnotprocessed+1
    174180      endif
    175181   enddo       
    … …  
    177183enddo
    178184
    179 ! Nnotconcat: # of variables that won't be processed
     185! Nnotprocessed: # of variables that won't be processed
    180186! nbvarfile: total # of variables in file
    181 allocate(var(nbvarfile-Nnotconcat),stat=ierr)
     187allocate(var(nbvarfile-Nnotprocessed),stat=ierr)
    182188if (ierr.ne.0) then
    183   write(*,*) "Error: failed allocation of var(nbvarfile-Nnotconcat)"
     189  write(*,*) "Error: failed allocation of var(nbvarfile-Nnotprocessed)"
    184190  write(*,*) "  nbvarfile=",nbvarfile
    185   write(*,*) "  Nnotconcat=",Nnotconcat
     191  write(*,*) "  Nnotprocessed=",Nnotprocessed
    186192  stop
    187193endif
    … …  
    201207
    202208if (tmpvar=="all") then
    203          nbvar=nbvarfile-Nnotconcat
    204          do j=Nnotconcat+1,nbvarfile
    205             ierr=nf_inq_varname(nid,j,var(j-Nnotconcat))
     209         nbvar=nbvarfile-Nnotprocessed
     210         do j=Nnotprocessed+1,nbvarfile
     211            ierr=nf_inq_varname(nid,j,var(j-Nnotprocessed))
    206212         enddo
    207213! Variables names from the file are stored in var()
    208    nbvar=nbvarfile-Nnotconcat
     214   nbvar=nbvarfile-Nnotprocessed
    209215   do i=1,nbvar
    210       ierr=nf_inq_varname(nid,i+Nnotconcat,var(i))
     216      ierr=nf_inq_varname(nid,i+Nnotprocessed,var(i))
    211217      write(*,'(a9,1x,i2,1x,a1,1x,a50)') "variable ",i,":",var(i)
    212218   enddo
    … …  
    220226!==============================================================================
    221227 filename=file(1:len_trim(file)-3)//"_LT.nc"
    222  write(*,*) trim(filename)
     228 write(*,*) "Output file :"//trim(filename)
    223229
    224230!==============================================================================
    … …  
    226232!==============================================================================
    227233
    228       write(*,*)
    229       write(*,*) "opening "//trim(file)//"..."
    230       ierr = NF_OPEN(file,NF_NOWRITE,nid)
    231       if (ierr.NE.NF_NOERR) then
    232          write(*,*) 'ERROR: Pb opening file '//trim(file)
    233          write(*,*) NF_STRERROR(ierr)
    234          stop ""
    235       endif
     234   write(*,*)
     235   write(*,*) "opening "//trim(file)//"..."
     236   ierr = NF_OPEN(file,NF_NOWRITE,nid)
     237   if (ierr.NE.NF_NOERR) then
     238      write(*,*) 'ERROR: Pb opening file '//trim(file)
     239      write(*,*) NF_STRERROR(ierr)
     240      stop ""
     241   endif
    236242
    237243!==============================================================================
    … …  
    277283   ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
    278284!  write(*,*) "altlen: ",altlen
     285
     286! load size of aps() or sigma() (in case it is not altlen)
     287   ! default is that GCM_layers=altlen
     288   ! but for outputs of zrecast, it may be a different value
     289   ierr=NF_INQ_DIMID(nid,"GCM_layers",gcmlayerdim)
     290   if (ierr.ne.NF_NOERR) then
     291     ! didn't find a GCM_layers dimension; therefore we have:
     292     GCM_layers=altlen
     293   else
     294     ! load value of GCM_layers
     295     ierr=NF_INQ_DIMLEN(nid,gcmlayerdim,GCM_layers)
     296   endif
     297!   write(*,*) "GCM_layers=",GCM_layers
    279298
    280299   ierr=NF_INQ_DIMID(nid,"index",ctldim)
    … …  
    299318!==============================================================================
    300319
    301 ! First call; initialize/allocate
    302       allocate(lat(latlen),stat=ierr)
    303       if (ierr.ne.0) then
    304         write(*,*) "Error: failed to allocate lat(latlen)"
    305         stop
    306       endif
    307       allocate(lon(lonlen),stat=ierr)
    308       if (ierr.ne.0) then
    309         write(*,*) "Error: failed to allocate lon(lonlen)"
    310         stop
    311       endif
    312       allocate(alt(altlen),stat=ierr)
    313       if (ierr.ne.0) then
    314         write(*,*) "Error: failed to allocate alt(altlen)"
    315         stop
    316       endif
    317       allocate(ctl(ctllen),stat=ierr)
    318       if (ierr.ne.0) then
    319         write(*,*) "Error: failed to allocate ctl(ctllen)"
    320         stop
    321       endif
    322 
    323 #ifdef NC_DOUBLE
    324       ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat)
    325 #else
    326       ierr = NF_GET_VAR_REAL(nid,latvar,lat)
    327 #endif
    328       if (ierr.ne.0) then
    329         write(*,*) "Error: failed to load latitude"
    330         write(*,*) NF_STRERROR(ierr)
    331         stop
    332       endif
    333 
    334 #ifdef NC_DOUBLE
    335       ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon)
    336 #else
    337       ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
    338 #endif
    339       if (ierr.ne.0) then
    340         write(*,*) "Error: failed to load longitude"
    341         write(*,*) NF_STRERROR(ierr)
    342         stop
    343       endif
    344 
    345 #ifdef NC_DOUBLE
    346       ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt)
    347 #else
    348       ierr = NF_GET_VAR_REAL(nid,altvar,alt)
    349 #endif
    350       if (ierr.ne.0) then
    351         write(*,*) "Error: failed to load altitude"
    352         write(*,*) NF_STRERROR(ierr)
    353         stop
    354       endif
     320   allocate(lat(latlen),stat=ierr)
     321   if (ierr.ne.0) then
     322     write(*,*) "Error: failed to allocate lat(latlen)"
     323     stop
     324   endif
     325   allocate(lon(lonlen),stat=ierr)
     326   if (ierr.ne.0) then
     327     write(*,*) "Error: failed to allocate lon(lonlen)"
     328     stop
     329   endif
     330   allocate(alt(altlen),stat=ierr)
     331   if (ierr.ne.0) then
     332     write(*,*) "Error: failed to allocate alt(altlen)"
     333     stop
     334   endif
     335   allocate(ctl(ctllen),stat=ierr)
     336   if (ierr.ne.0) then
     337     write(*,*) "Error: failed to allocate ctl(ctllen)"
     338     stop
     339   endif
     340
     341   ierr = NF_GET_VAR_REAL(nid,latvar,lat)
     342   if (ierr.ne.0) then
     343     write(*,*) "Error: failed to load latitude"
     344     write(*,*) NF_STRERROR(ierr)
     345     stop
     346   endif
     347
     348   ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
     349   if (ierr.ne.0) then
     350     write(*,*) "Error: failed to load longitude"
     351     write(*,*) NF_STRERROR(ierr)
     352     stop
     353   endif
     354
     355   ierr = NF_GET_VAR_REAL(nid,altvar,alt)
     356   if (ierr.ne.0) then
     357     write(*,*) "Error: failed to load altitude"
     358     write(*,*) NF_STRERROR(ierr)
     359     stop
     360   endif
    355361! Get altitude attributes to handle files with any altitude type
    356 ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
    357 ierr=nf_get_att_text(nid,altvar,'units',altunits)
    358 ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
    359 
    360       if (ctllen .ne. 0) then
    361 #ifdef NC_DOUBLE
    362         ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
    363 #else
    364         ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
    365 #endif
    366         if (ierr.ne.0) then
    367           write(*,*) "Error: failed to load controle"
    368           write(*,*) NF_STRERROR(ierr)
    369           stop
    370         endif
    371       endif ! of if (ctllen .ne. 0)
     362   ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
     363   if (ierr.ne.nf_noerr) then
     364   ! if no attribute "long_name", try "title"
     365     ierr=nf_get_att_text(nid,timevar,"title",long_name)
     366   endif
     367   ierr=nf_get_att_text(nid,altvar,'units',altunits)
     368   ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
     369
     370   if (ctllen .ne. 0) then
     371     ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
     372     if (ierr.ne.0) then
     373       write(*,*) "Error: failed to load controle"
     374       write(*,*) NF_STRERROR(ierr)
     375       stop
     376     endif
     377   endif ! of if (ctllen .ne. 0)
    372378
    373379!==============================================================================
    … …  
    391397!  write(*,*) "timelen: ",timelen
    392398
    393    ! Sanity Check: Does "Time" has a "title" attribute
     399   ! Sanity Check: Does "Time" has a "long_name" attribute
    394400   ! and is it "Solar longitude" ?
    395    ierr=nf_get_att_text(nid,timevar,"title",title)
    396    if ((ierr.EQ.NF_NOERR).and.(title.eq."Solar longitude")) then
     401   ierr=nf_get_att_text(nid,timevar,"long_name",long_name)
     402   if (ierr.ne.nf_noerr) then
     403   ! if no attribute "long_name", try "title"
     404     ierr=nf_get_att_text(nid,timevar,"title",long_name)
     405   endif
     406   if ((ierr.EQ.NF_NOERR).and.(long_name.eq."Solar longitude")) then
    397407     write(*,*) "ERROR: Time axis in input file is in Solar Longitude!"
    398408     write(*,*) "       localtime requires sols as time axis!"
    … …  
    408418   endif
    409419
    410 #ifdef NC_DOUBLE
    411    ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
    412 #else
    413420   ierr = NF_GET_VAR_REAL(nid,timevar,time)
    414 #endif
    415421   if (ierr.NE.NF_NOERR) then
    416422     write(*,*) "Error , failed to load Time"
    … …  
    432438      write(*,*) 'list of selected local time (0<t<24)'
    433439      do it=1,nhour
    434          read(*,*) lt_hour(it)
     440        read(*,*) lt_hour(it)
    435441      end do
    436442
    437443      if ((timelen.eq.12).and.(time(1).eq.2).and.(time(timelen).eq.24))then
    438          write(*,*) 'We detect a ""stats"" file'
    439          stats=.true.
    440          timelen_lt=nhour
    441          allocate(lt_out(timelen_lt),stat=ierr)
    442          if (ierr.ne.0) then
    443            write(*,*) "Error: failed to allocate lt_hour(nhour)"
    444            stop
    445          endif
    446          do it=1,nhour
    447            lt_out(it)=lt_hour(it)/24.
    448          end do
     444        write(*,*) 'We detect a ""stats"" file'
     445        stats=.true.
     446        timelen_lt=nhour
     447        allocate(lt_out(timelen_lt),stat=ierr)
     448        if (ierr.ne.0) then
     449          write(*,*) "Error: failed to allocate lt_hour(nhour)"
     450          stop
     451        endif
     452        do it=1,nhour
     453          lt_out(it)=lt_hour(it)/24.
     454        end do
    449455!        We rewrite the time from "stats" from 0 to 1 sol...
    450          do it=1,timelen  ! loop temps in
    451                time(it) = time(it)/24.
    452          end do
    453          nsol =1
     456        do it=1,timelen  ! loop temps in
     457              time(it) = time(it)/24.
     458        end do
     459        nsol =1
    454460      else   ! case of a diagfi or concatnc file
    455461        stats=.false.
    … …  
    469475          end do
    470476        end do
    471         end if
     477      end if
    472478
    473479      if (nsol.gt.1) then
    … …  
    499505
    500506   ! Initialize output file's lat,lon,alt and time dimensions
    501      call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
     507     call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,GCM_layers,&
    502508           altlong_name,altunits,altpositive,&
    503            nout,latdimout,londimout,altdimout,timedimout,timevarout)
     509           nout,latdimout,londimout,altdimout,timedimout,&
     510           layerdimout,timevarout)
    504511
    505512   ! Initialize output file's aps,bps and phisinit variables
    506      call init2(nid,lonlen,latlen,altlen,altdim,&
    507                 nout,londimout,latdimout,altdimout)
     513     call init2(nid,lonlen,latlen,altlen,GCM_layers,&
     514                nout,londimout,latdimout,altdimout,&
     515                layerdimout)
    508516
    509517
    … …  
    515523
    516524     do it=1,timelen_lt
    517 #ifdef NC_DOUBLE
    518         ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,lt_hour(it))
    519 #else
    520525        ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,lt_hour(it))
    521 #endif
    522526     enddo
    523527   else
    524528     do it=1,timelen_lt
    525 #ifdef NC_DOUBLE
    526         ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,lt_out(it))
    527 #else
    528529        ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,lt_out(it))
    529 #endif
    530530        if (ierr.NE.NF_NOERR) then
    531531          write(*,*) "Error , failed to write Time"
    … …  
    553553      ierr=nf_inq_varndims(nid,varid,ndim)
    554554
     555      ! Check that it is a 3D or 4D variable
     556      if (ndim.lt.3) then
     557        write(*,*) "Error:",trim(var(j))," is nor a 3D neither a 4D variable"
     558        write(*,*) "We'll skip that variable..."
     559        CYCLE ! go directly to the next variable
     560      endif
    555561!==============================================================================
    556562! 2.5.2 Prepare things in order to read/write the variable
    … …  
    599605      endif
    600606
    601          units=" "
    602          title=" "
    603          ierr=nf_get_att_text(nid,varid,"title",title)
    604          ierr=nf_get_att_text(nid,varid,"units",units)
    605          call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr)
     607      units=" "
     608      long_name=" "
     609      ierr=nf_get_att_text(nid,timevar,"long_name",long_name)
     610      if (ierr.ne.nf_noerr) then
     611      ! if no attribute "long_name", try "title"
     612        ierr=nf_get_att_text(nid,timevar,"title",long_name)
     613      endif
     614      ierr=nf_get_att_text(nid,varid,"units",units)
     615      call def_var(nout,var(j),long_name,units,ndim,dim,varidout,ierr)
    606616
    607617     
    … …  
    611621!==============================================================================
    612622
    613 #ifdef NC_DOUBLE
    614       ierr = NF_GET_VAR_DOUBLE(nid,varid,var3d)
    615       miss=NF_GET_ATT_DOUBLE(nid,varid,"missing_value",missing)
    616       validr=NF_GET_ATT_DOUBLE(nid,varid,"valid_range",valid_range)
    617 #else
    618623      ierr = NF_GET_VAR_REAL(nid,varid,var3d)
    619624      miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
    620625      validr=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
    621 #endif
    622626
    623627!==============================================================================
    … …  
    677681!==============================================================================
    678682
    679 #ifdef NC_DOUBLE
    680       ierr= NF_PUT_VARA_DOUBLE(nout,varidout,corner,edges,var3d_lt)
    681 #else
    682683      ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3d_lt)
    683 #endif
    684684
    685685      if (ierr.ne.NF_NOERR) then
    … …  
    712712ierr=NF_CLOSE(nout)
    713713
     714write(*,*) "Program completed !"
     715
    714716contains
    715717
    716718!******************************************************************************
    717 Subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
     719Subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,GCM_layers,&
    718720                     altlong_name,altunits,altpositive,&
    719                      nout,latdimout,londimout,altdimout,timedimout,timevarout)
     721                     nout,latdimout,londimout,altdimout,timedimout,&
     722                     layerdimout,timevarout)
    720723!==============================================================================
    721724! Purpose:
    … …  
    744747real, intent(in):: ctl(ctllen)
    745748! ctl(): controle
     749integer,intent(in) :: GCM_layers ! number of GCM layers
    746750character (len=*), intent(in) :: altlong_name
    747751! altlong_name(): [netcdf] altitude "long_name" attribute
    … …  
    760764integer, intent(out):: timedimout
    761765! timedimout: [netcdf] "Time"  ID
     766integer,intent(out) :: layerdimout
     767! layerdimout: [netcdf] "GCM_layers" ID
    762768integer, intent(out):: timevarout
    763769! timevarout: [netcdf] "Time" (considered as a variable) ID
    … …  
    821827endif
    822828
     829ierr = NF_DEF_DIM(nout, "GCM_layers", GCM_layers, layerdimout)
     830
    823831! End netcdf define mode
    824832ierr = NF_ENDDEF(nout)
    … …  
    843851             (/latdimout/),nvarid,ierr)
    844852
    845 #ifdef NC_DOUBLE
    846 ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lat)
    847 #else
    848853ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
    849 #endif
    850854if (ierr.NE.NF_NOERR) then
    851855   WRITE(*,*)'initiate: error failed writing variable <latitude>'
    … …  
    861865             (/londimout/),nvarid,ierr)
    862866
    863 #ifdef NC_DOUBLE
    864 ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lon)
    865 #else
     867
    866868ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
    867 #endif
    868869if (ierr.NE.NF_NOERR) then
    869870   WRITE(*,*)'initiate: error failed writing variable <longitude>'
    … …  
    879880ierr = NF_REDEF (nout)
    880881
    881 #ifdef NC_DOUBLE
    882 ierr = NF_DEF_VAR (nout,"altitude",NF_DOUBLE,1,altdimout,nvarid)
    883 #else
    884882ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,altdimout,nvarid)
    885 #endif
    886883
    887884ierr = NF_PUT_ATT_TEXT (nout,nvarid,'long_name',len_trim(adjustl(altlong_name)),adjustl(altlong_name))
    … …  
    892889ierr = NF_ENDDEF(nout)
    893890
    894 #ifdef NC_DOUBLE
    895 ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,alt)
    896 #else
    897 ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
    898 #endif
     891ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
    899892if (ierr.NE.NF_NOERR) then
    900893   WRITE(*,*)'initiate: error failed writing variable <altitude>'
    … …  
    911904   ierr = NF_REDEF (nout)
    912905
    913 #ifdef NC_DOUBLE
    914    ierr = NF_DEF_VAR (nout,"controle",NF_DOUBLE,1,ctldimout,nvarid)
    915 #else
    916906   ierr = NF_DEF_VAR (nout,"controle",NF_FLOAT,1,ctldimout,nvarid)
    917 #endif
    918 
    919    ierr = NF_PUT_ATT_TEXT (nout,nvarid,"title",18,"Control parameters")
     907
     908   ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",18,"Control parameters")
    920909
    921910   ! End netcdf define mode
    922911   ierr = NF_ENDDEF(nout)
    923912
    924 #ifdef NC_DOUBLE
    925    ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,ctl)
    926 #else
    927913   ierr = NF_PUT_VAR_REAL (nout,nvarid,ctl)
    928 #endif
    929914   if (ierr.NE.NF_NOERR) then
    930915      WRITE(*,*)'initiate: error failed writing variable <controle>'
    … …  
    936921end Subroutine initiate
    937922!******************************************************************************
    938 subroutine init2(infid,lonlen,latlen,altlen,altdimid, &
    939                  outfid,londimout,latdimout,altdimout)
     923subroutine init2(infid,lonlen,latlen,altlen,GCM_layers,&
     924                 outfid,londimout,latdimout,altdimout,&
     925                 layerdimout)
    940926!==============================================================================
    941927! Purpose:
    … …  
    957943integer, intent(in) :: latlen ! # of grid points along latitude
    958944integer, intent(in) :: altlen ! # of grid points along altitude
    959 integer, intent(in) :: altdimid ! ID of altitude dimension
     945integer, intent(in) :: GCM_layers ! # of GCM atmospheric layers
    960946integer, intent(in) :: outfid ! NetCDF output file ID
    961947integer, intent(in) :: londimout ! longitude dimension ID
    962948integer, intent(in) :: latdimout ! latitude dimension ID
    963949integer, intent(in) :: altdimout ! altitude dimension ID
     950integer, intent(in) :: layerdimout ! GCM_layers dimension ID
    964951!==============================================================================
    965952! Local variables:
    … …  
    971958integer :: tmpdimid ! temporary dimension ID
    972959integer :: tmpvarid ! temporary variable ID
    973 logical :: phis, aps_ok, bps_ok ! are "phisinit" "aps" "bps" available ?
     960logical :: phis, hybrid ! are "phisinit" and "aps"/"bps" available ?
    974961
    975962
    … …  
    979966
    980967! hybrid coordinate aps
    981   allocate(aps(altlen),stat=ierr)
     968!write(*,*) "aps: altlen=",altlen," GCM_layers=",GCM_layers
     969allocate(aps(GCM_layers),stat=ierr)
     970if (ierr.ne.0) then
     971  write(*,*) "init2: failed to allocate aps!"
     972  stop
     973endif
     974ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
     975if (ierr.ne.NF_NOERR) then
     976  write(*,*) "Ooops. Failed to get aps ID. OK."
     977  hybrid=.false.
     978else
     979  hybrid=.true.
     980  ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
     981  hybrid=.true.
     982  if (ierr.ne.NF_NOERR) then
     983    stop "init2 Error: Failed reading aps"
     984  endif
     985
     986  ! hybrid coordinate bps
     987!  write(*,*) "bps: altlen=",altlen," GCM_layers=",GCM_layers
     988  allocate(bps(GCM_layers),stat=ierr)
    982989  if (ierr.ne.0) then
    983     write(*,*) "init2: failed to allocate aps(altlen)"
     990    write(*,*) "init2: failed to allocate bps!"
    984991    stop
    985992  endif
    986 
    987 ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
    988 if (ierr.ne.NF_NOERR) then
    989   write(*,*) "oops Failed to get aps ID. OK"
    990   aps_ok=.false.
    991 else
    992   ! Check that aps() is the right size (it most likely won't be if
    993   ! zrecast has be used to generate the input file)
    994   ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
    995   if (tmpdimid.ne.altdimid) then
    996     write(*,*) "init2: wrong dimension size for aps(), skipping it ..."
    997     aps_ok=.false.
    998   else
    999     ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
    1000     if (ierr.ne.NF_NOERR) then
    1001      stop "init2 error: Failed reading aps"
    1002     endif
    1003     aps_ok=.true.
    1004   endif
    1005 endif
    1006 
    1007 ! hybrid coordinate bps
    1008   allocate(bps(altlen),stat=ierr)
    1009   if (ierr.ne.0) then
    1010     write(*,*) "init2: failed to allocate bps(altlen)"
    1011     stop
    1012   endif
    1013 
    1014 ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
    1015 if (ierr.ne.NF_NOERR) then
    1016   write(*,*) "oops: Failed to get bps ID. OK"
    1017   bps_ok=.false.
    1018 else
    1019   ! Check that bps() is the right size
    1020   ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
    1021   if (tmpdimid.ne.altdimid) then
    1022     write(*,*) "init2: wrong dimension size for bps(), skipping it ..."
    1023     bps_ok=.false.
     993  ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
     994  if (ierr.ne.NF_NOERR) then
     995    write(*,*) "Ooops. Failed to get bps ID. OK."
     996    hybrid=.false.
    1024997  else
    1025998    ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
    … …  
    10271000      stop "init2 Error: Failed reading bps"
    10281001    endif
    1029     bps_ok=.true.
    10301002  endif
    10311003endif
    … …  
    10591031!==============================================================================
    10601032
    1061 IF (aps_ok) then
     1033IF (hybrid) then
    10621034  ! define aps
    10631035  call def_var(outfid,"aps","hybrid pressure at midlayers"," ",1,&
    1064                (/altdimout/),apsid,ierr)
     1036               (/layerdimout/),apsid,ierr)
    10651037  if (ierr.ne.NF_NOERR) then
    10661038    stop "Error: Failed to def_var aps"
    … …  
    10681040
    10691041  ! write aps
    1070 #ifdef NC_DOUBLE
    1071   ierr=NF_PUT_VAR_DOUBLE(outfid,apsid,aps)
    1072 #else
    10731042  ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
    1074 #endif
    10751043  if (ierr.ne.NF_NOERR) then
    10761044    stop "Error: Failed to write aps"
    10771045  endif
    1078 ENDIF ! of IF (aps_ok)
    1079 
    1080 IF (bps_ok) then
     1046 
    10811047  ! define bps
    10821048  call def_var(outfid,"bps","hybrid sigma at midlayers"," ",1,&
    1083                (/altdimout/),bpsid,ierr)
     1049               (/layerdimout/),bpsid,ierr)
    10841050  if (ierr.ne.NF_NOERR) then
    10851051    stop "Error: Failed to def_var bps"
    … …  
    10871053
    10881054  ! write bps
    1089 #ifdef NC_DOUBLE
    1090   ierr=NF_PUT_VAR_DOUBLE(outfid,bpsid,bps)
    1091 #else
    10921055  ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
    1093 #endif
    10941056  if (ierr.ne.NF_NOERR) then
    10951057    stop "Error: Failed to write bps"
    10961058  endif
    1097 ENDIF ! of IF (bps_ok)
     1059 
     1060  deallocate(aps)
     1061  deallocate(bps)
     1062ENDIF ! of IF (hybrid)
    10981063
    10991064!==============================================================================
    … …  
    11111076
    11121077  ! write phisinit
    1113 #ifdef NC_DOUBLE
    1114   ierr=NF_PUT_VAR_DOUBLE(outfid,phisinitid,phisinit)
    1115 #else
    11161078  ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
    1117 #endif
    11181079  if (ierr.ne.NF_NOERR) then
    11191080    stop "Error: Failed to write phisinit"
    11201081  endif
    11211082
     1083  deallocate(phisinit)
    11221084ENDIF ! of IF (phis)
    11231085
    11241086
    1125 ! Cleanup
    1126 deallocate(aps)
    1127 deallocate(bps)
    1128 deallocate(phisinit)
    1129 
    11301087end subroutine init2
     1088
    11311089!******************************************************************************
    1132 subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
     1090subroutine def_var(nid,name,long_name,units,nbdim,dim,nvarid,ierr)
    11331091!==============================================================================
    11341092! Purpose: Write a variable (i.e: add a variable to a dataset)
    1135 ! called "name"; along with its attributes "title", "units"...
     1093! called "name"; along with its attributes "long_name", "units"...
    11361094! to a file (following the NetCDF format)
    11371095!==============================================================================
    … …  
    11511109character (len=*), intent(in) :: name
    11521110! name(): [netcdf] variable's name
    1153 character (len=*), intent(in) :: title
    1154 ! title(): [netcdf] variable's "title" attribute
     1111character (len=*), intent(in) :: long_name
     1112! long_name(): [netcdf] variable's "long_name" attribute
    11551113character (len=*), intent(in) :: units
    11561114! unit(): [netcdf] variable's "units" attribute
    … …  
    11681126
    11691127! Insert the definition of the variable
    1170 #ifdef NC_DOUBLE
    1171 ierr=NF_DEF_VAR(nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
    1172 #else
    11731128ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
    1174 #endif
    11751129
    11761130! Write the attributes
    1177 ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",len_trim(adjustl(title)),adjustl(title))
     1131ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name",len_trim(adjustl(long_name)),adjustl(long_name))
    11781132ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
    11791133
    … …  
    12321186! Write valid_range() attribute
    12331187if (validr.eq.NF_NOERR) then
    1234 #ifdef NC_DOUBLE
    1235   ierr=NF_PUT_ATT_DOUBLE(nout,nvarid,'valid_range',NF_DOUBLE,2,valid_range)
    1236 #else
    12371188  ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
    1238 #endif
    12391189
    12401190  if (ierr.ne.NF_NOERR) then
    … …  
    12461196  endif
    12471197endif
    1248 
     1198if (miss.eq.NF_NOERR) then
    12491199! Write "missing_value" attribute
    1250 if (miss.eq.NF_NOERR) then
    1251 #ifdef NC_DOUBLE
    1252   ierr= NF_PUT_ATT_DOUBLE(nout,nvarid,'missing_value',NF_DOUBLE,1,missing)
    1253 #else
    12541200  ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
    1255 #endif
    12561201
    12571202  if (ierr.NE.NF_NOERR) then

  • trunk/LMDZ.MARS/util/solzenangle.F90

    r2416 r2434  
    1515character (len=50)  file
    1616! file(): input file(s) names(s)
    17 character (len=30), dimension(15) :: notconcat
    18 ! notconcat(): names of the (15) variables that won't be concatenated
     17character (len=30), dimension(16) :: notprocessed
     18! notprocessed(): names of the (16) variables that won't be processed
    1919character (len=50), dimension(:), allocatable :: var
    20 ! var(): name(s) of variable(s) that will be concatenated
    21 character (len=50) :: tmpvar,title,units
     20! var(): name(s) of variable(s) that will be processed
     21character (len=100) :: tmpvar,long_name,units
    2222! tmpvar(): used to temporarily store a variable name
    23 ! title(): [netcdf] title attribute
     23! long_name(): [netcdf] long_name attribute
    2424! units(): [netcdf] units attribute
    2525character (len=100) :: filename,vartmp
    … …  
    6262! ctldim: [netcdf] "controle" dim ID
    6363! timedim: [netcdf] "timedim" dim ID
     64integer :: gcmlayerdim ! NetCDF dimension ID for # of layers in GCM
    6465integer :: latvar,lonvar,altvar,ctlvar,timevar
    6566! latvar: [netcdf] ID of "latitude" variable
    … …  
    7273! latlen: # of elements of lat() array
    7374! lonlen: # of elements of lon() array
    74 ! altvar: # of elements of alt() array
    75 ! ctlvar: # of elements of ctl() array
     75! altlen: # of elements of alt() array
     76! ctllen: # of elements of ctl() array
    7677! timelen: # of elements of time() array
    7778! timelen_tot:# =timelen or timelen+1 (if 1 more time to interpolate needed)
    7879integer :: nsol
    7980! nsol: # of sols in the input file AND # of elements of time() array in output
     81integer :: GCM_layers ! number of GCM atmospheric layers (may not be
     82! same as altlen if file is an output of zrecast)
    8083integer :: nout,latdimout,londimout,altdimout,timedimout,timevarout
    8184! nout: [netcdf] output file ID
    … …  
    8588! timedimout: [netcdf] output time (dimension) ID
    8689! timevarout: [netcdf] ID of output "Time" variable
     90integer :: layerdimout ! NetCDF dimension ID for # of layers in GCM
    8791integer :: varidout
    8892! varidout: [netcdf] variable ID # (of a variable to write to the output file)
    89 integer :: Nnotconcat,var_ok
    90 ! Nnotconcat: # of (leading)variables that won't be concatenated
     93integer :: Nnotprocessed,var_ok
     94! Nnotprocessed: # of (leading)variables that won't be concatenated
    9195! var_ok: flag (0 or 1)
    9296integer, dimension(4) :: corner,edges,dim
    … …  
    107111real :: sza ! solar zenith angle
    108112character (len=30) :: szastr ! to store the sza value in a string
    109 real :: starttimeoffset ! offset (in sols) of sol 0 in input file, wrt Ls=0
     113real :: startsol ! starting date (in sols) of sol 0 in input file, wrt Ls=0
    110114real :: tmpsol ! to temporarily store a sol value
    111115real :: tmpLs ! to temporarily store a Ls value
    … …  
    121125! valid_range(2): [netcdf] interval in which a value is considered valid
    122126logical :: stats   ! stats=T when reading a "stats" kind of ile
     127real :: year_day_gcm = 669. ! number of sols in a year in GCM
     128
     129
    123130
    124131
    … …  
    160167!==============================================================================
    161168
    162 notconcat(1)='Time'
    163 notconcat(2)='controle'
    164 notconcat(3)='rlonu'
    165 notconcat(4)='latitude'
    166 notconcat(5)='longitude'
    167 notconcat(6)='altitude'
    168 notconcat(7)='rlatv'
    169 notconcat(8)='aps'
    170 notconcat(9)='bps'
    171 notconcat(10)='ap'
    172 notconcat(11)='bp'
    173 notconcat(12)='cu'
    174 notconcat(13)='cv'
    175 notconcat(14)='aire'
    176 notconcat(15)='phisinit'
     169notprocessed(1)='Time'
     170notprocessed(2)='controle'
     171notprocessed(3)='rlonu'
     172notprocessed(4)='latitude'
     173notprocessed(5)='longitude'
     174notprocessed(6)='altitude'
     175notprocessed(7)='rlatv'
     176notprocessed(8)='aps'
     177notprocessed(9)='bps'
     178notprocessed(10)='ap'
     179notprocessed(11)='bp'
     180notprocessed(12)='soildepth'
     181notprocessed(13)='cu'
     182notprocessed(14)='cv'
     183notprocessed(15)='aire'
     184notprocessed(16)='phisinit'
    177185
    178186!==============================================================================
    … …  
    180188!==============================================================================
    181189write(*,*)
    182    Nnotconcat=0
     190   Nnotprocessed=0
    183191do i=1,nbvarfile
    184192   ierr=NF_INQ_VARNAME(nid,i,vartmp)
    185193   ! vartmp now contains the "name" of variable of ID # i
    186194   var_ok=0
    187    do inter=1,15
    188       if (vartmp.eq.notconcat(inter)) then
     195   do inter=1,16
     196      if (vartmp.eq.notprocessed(inter)) then
    189197         var_ok=1
    190          Nnotconcat=Nnotconcat+1
     198         Nnotprocessed=Nnotprocessed+1
    191199      endif
    192200   enddo       
    … …  
    194202enddo
    195203
    196 ! Nnotconcat: # of variables that won't be processed
     204! Nnotprocessed: # of variables that won't be processed
    197205! nbvarfile: total # of variables in file
    198 allocate(var(nbvarfile-Nnotconcat),stat=ierr)
     206allocate(var(nbvarfile-Nnotprocessed),stat=ierr)
    199207if (ierr.ne.0) then
    200   write(*,*) "Error: failed allocation of var(nbvarfile-Nnotconcat)"
     208  write(*,*) "Error: failed allocation of var(nbvarfile-Nnotprocessed)"
    201209  write(*,*) "  nbvarfile=",nbvarfile
    202   write(*,*) "  Nnotconcat=",Nnotconcat
     210  write(*,*) "  Nnotprocessed=",Nnotprocessed
    203211  stop
    204212endif
    … …  
    218226
    219227if (tmpvar=="all") then
    220          nbvar=nbvarfile-Nnotconcat
    221          do j=Nnotconcat+1,nbvarfile
    222             ierr=nf_inq_varname(nid,j,var(j-Nnotconcat))
     228         nbvar=nbvarfile-Nnotprocessed
     229         do j=Nnotprocessed+1,nbvarfile
     230            ierr=nf_inq_varname(nid,j,var(j-Nnotprocessed))
    223231         enddo
    224232! Variables names from the file are stored in var()
    225    nbvar=nbvarfile-Nnotconcat
     233   nbvar=nbvarfile-Nnotprocessed
    226234   do i=1,nbvar
    227       ierr=nf_inq_varname(nid,i+Nnotconcat,var(i))
     235      ierr=nf_inq_varname(nid,i+Nnotprocessed,var(i))
    228236      write(*,'(a9,1x,i2,1x,a1,1x,a50)') "variable ",i,":",var(i)
    229237   enddo
    … …  
    290298!  write(*,*) "altlen: ",altlen
    291299
     300! load size of aps() or sigma() (in case it is not altlen)
     301   ! default is that GCM_layers=altlen
     302   ! but for outputs of zrecast, it may be a different value
     303   ierr=NF_INQ_DIMID(nid,"GCM_layers",gcmlayerdim)
     304   if (ierr.ne.NF_NOERR) then
     305     ! didn't find a GCM_layers dimension; therefore we have:
     306     GCM_layers=altlen
     307   else
     308     ! load value of GCM_layers
     309     ierr=NF_INQ_DIMLEN(nid,gcmlayerdim,GCM_layers)
     310   endif
     311!   write(*,*) "GCM_layers=",GCM_layers
     312
    292313   ierr=NF_INQ_DIMID(nid,"index",ctldim)
    293314   if (ierr.NE.NF_NOERR) then
    … …  
    304325      ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
    305326   endif
    306 !  write(*,*) "controle: ",controle
     327
    307328
    308329!==============================================================================
    … …  
    311332!==============================================================================
    312333
    313 ! First call; initialize/allocate
    314334   allocate(lat(latlen),stat=ierr)
    315335   if (ierr.ne.0) then
    … …  
    333353   endif
    334354
    335 #ifdef NC_DOUBLE
    336    ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat)
    337 #else
    338355   ierr = NF_GET_VAR_REAL(nid,latvar,lat)
    339 #endif
    340356   if (ierr.ne.0) then
    341357     write(*,*) "Error: failed to load latitude"
    … …  
    344360   endif
    345361
    346 #ifdef NC_DOUBLE
    347    ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon)
    348 #else
    349362   ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
    350 #endif
    351363   if (ierr.ne.0) then
    352364     write(*,*) "Error: failed to load longitude"
    … …  
    355367   endif
    356368
    357 #ifdef NC_DOUBLE
    358    ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt)
    359 #else
    360369   ierr = NF_GET_VAR_REAL(nid,altvar,alt)
    361 #endif
    362370   if (ierr.ne.0) then
    363371     write(*,*) "Error: failed to load altitude"
    … …  
    367375! Get altitude attributes to handle files with any altitude type
    368376   ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
     377   if (ierr.ne.nf_noerr) then
     378   ! if no attribute "long_name", try "title"
     379     ierr=nf_get_att_text(nid,altvar,"title",altlong_name)
     380   endif
    369381   ierr=nf_get_att_text(nid,altvar,'units',altunits)
    370382   ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
    371383
    372384   if (ctllen .ne. 0) then
    373 #ifdef NC_DOUBLE
    374      ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
    375 #else
    376385     ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
    377 #endif
    378386     if (ierr.ne.0) then
    379387       write(*,*) "Error: failed to load controle"
    … …  
    403411!  write(*,*) "timelen: ",timelen
    404412
    405    ! Sanity Check: Does "Time" has a "title" attribute
     413   ! Sanity Check: Does "Time" has a "long_name" attribute
    406414   ! and is it "Solar longitude" ?
    407    ierr=nf_get_att_text(nid,timevar,"title",title)
    408    if ((ierr.EQ.NF_NOERR).and.(title.eq."Solar longitude")) then
     415   ierr=nf_get_att_text(nid,timevar,"long_name",long_name)
     416   if (ierr.ne.nf_noerr) then
     417   ! if no attribute "long_name", try "title"
     418     ierr=nf_get_att_text(nid,timevar,"title",long_name)
     419   endif
     420   if ((ierr.EQ.NF_NOERR).and.(long_name.eq."Solar longitude")) then
    409421     write(*,*) "ERROR: Time axis in input file is in Solar Longitude!"
    410422     write(*,*) "       localtime requires sols as time axis!"
    … …  
    420432   endif
    421433
    422 #ifdef NC_DOUBLE
    423    ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
    424 #else
    425434   ierr = NF_GET_VAR_REAL(nid,timevar,time)
    426 #endif
    427435   if (ierr.NE.NF_NOERR) then
    428436     write(*,*) "Error , failed to load Time"
    … …  
    434442! 2.4.1 Select local times to be saved "Time" in output file
    435443!==============================================================================
    436    write(*,*) "Planet side of the Solar Zenith Angle ?"
     444   write(*,*) "Planet side of the Solar Zenith Angle ? (morning or evening)"
    437445   read(*,*) planetside
    438446   if ((trim(planetside).ne."morning").and.(trim(planetside).ne."evening")) then
    … …  
    448456     stop
    449457   endif
    450    write(*,*) "Beginning date of the simulation file "
    451    write(*,*) "(or median date of the simulation for a stats)?"
    452    write(*,*) "(i.e. number of sols since Ls=0 at the Time=0.0 in the input file)"
    453    read(*,*) starttimeoffset
     458
     459!  Sol and season: 
     460
     461   if (ctllen .ne. 0) then
     462      startsol = int(int(ctl(4))+int(time(1) + ctl(27)))
     463   else
     464      startsol = int(time(1))
     465   end if
    454466
    455467   if ((timelen.eq.12).and.(time(1).eq.2).and.(time(timelen).eq.24))then
    … …  
    466478     nsol = int(time(timelen)) - int(time(1))
    467479   end if
    468 
    469480   allocate(intsol(nsol),stat=ierr)
    470481   if (ierr.ne.0) then
    … …  
    477488     stop
    478489   endif
     490
     491   if (stats) then
     492         write(*,*) "What is the mean Ls of this stats file (degree) ?"
     493         read(*,*,iostat=ierr) tmpLs
     494         startsol=0.
     495   else
     496      write(*,*) "Beginning date (sol since Ls=0) of this file is :",startsol
     497      write(*,*) "If you do not want to change it,  just write 'n' (recommended)"
     498      write(*,*) "If you really wish to change it, write it now (not recommended!)"
     499      read(*,*,iostat=ierr) startsol
     500      if((ierr.eq.0).and.(ctllen.ne.0))  ctl(4)=startsol
     501      write(*,*) "Beginning date (sol since Ls=0) will be :",startsol
     502   end if
     503
    479504   do it=1,nsol
    480505     intsol(it)= int(time(1)) + it-1.
    481506     do ilon=1,lonlen
    482 !         For the computation of Ls, we try to take into account the rotation time
    483 !         of Mars it's supposed that the morning/evening
     507!         For the computation of Ls, we try to take into account the
     508!         rotation time of Mars. It is supposed that the morning/evening
    484509!         correspond to LT=6h/18h at the given longitude, which is
    485510!         then reported to a sol value at longitude 0
    486511       if (trim(planetside).eq."morning") then
    487          tmpsol = intsol(it) + (6.-lon(ilon)/15.)/24. + starttimeoffset
     512         tmpsol = intsol(it) + (6.-lon(ilon)/15.)/24. + startsol
    488513       else !if trim(planetside).eq."evening"
    489          tmpsol = intsol(it) + (18.-lon(ilon)/15.)/24. + starttimeoffset           
     514         tmpsol = intsol(it) + (18.-lon(ilon)/15.)/24. + startsol           
    490515       endif
    491        call sol2ls(tmpsol,tmpLs)
     516       if (.not.stats) call sol2ls(tmpsol,tmpLs)
    492517       do ilat=1,latlen
    493518!           Compute the Local Time of the solar zenith angle at a given longitude
    … …  
    541566
    542567  ! Initialize output file's lat,lon,alt and time dimensions
    543    call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
     568   call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,GCM_layers,&
    544569         altlong_name,altunits,altpositive,&
    545          nout,latdimout,londimout,altdimout,timedimout,timevarout)
    546 
     570         nout,latdimout,londimout,altdimout,timedimout,&
     571         layerdimout,timevarout)
     572         
    547573  ! Initialize output file's aps,bps and phisinit variables
    548    call init2(nid,lonlen,latlen,altlen,altdim,&
    549               nout,londimout,latdimout,altdimout)
    550 
     574   call init2(nid,lonlen,latlen,altlen,GCM_layers,&
     575              nout,londimout,latdimout,altdimout,&
     576              layerdimout)
    551577
    552578!==============================================================================
    … …  
    557583
    558584     do it=1,nsol
    559 #ifdef NC_DOUBLE
    560         ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,intsol(it)*24.)
    561 #else
    562         ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,intsol(it)*24.)
    563 #endif
     585       ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,intsol(it)*24.)
    564586     enddo
    565587   else
    566588     do it=1,nsol
    567 #ifdef NC_DOUBLE
    568         ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,intsol(it))
    569 #else
    570         ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,intsol(it))
    571 #endif
    572         if (ierr.NE.NF_NOERR) then
    573           write(*,*) "Error , failed to write Time"
    574           stop
    575         endif
     589       ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,intsol(it))
     590       if (ierr.NE.NF_NOERR) then
     591         write(*,*) "Error , failed to write Time"
     592         stop
     593       endif
    576594     enddo
    577595   end if
    … …  
    584602
    585603!==============================================================================
    586 ! 2.5.1 Check that variable to be read existe in input file
     604! 2.5.1 Check that variable to be read exists in input file
    587605!==============================================================================
    588606
    … …  
    595613      ierr=nf_inq_varndims(nid,varid,ndim)
    596614
     615      ! Check that it is a 3D or 4D variable
     616      if (ndim.lt.3) then
     617        write(*,*) "Error:",trim(var(j))," is nor a 3D neither a 4D variable"
     618        write(*,*) "We'll skip that variable..."
     619        CYCLE ! go directly to the next variable
     620      endif
     621     
    597622!==============================================================================
    598623! 2.5.2 Prepare things in order to read/write the variable
    … …  
    642667
    643668      units=" "
    644       title=" "
    645       ierr=nf_get_att_text(nid,varid,"title",title)
     669      long_name=" "
     670      ierr=nf_get_att_text(nid,varid,"long_name",long_name)
     671      if (ierr.ne.nf_noerr) then
     672      ! if no attribute "long_name", try "title"
     673        ierr=nf_get_att_text(nid,varid,"title",long_name)
     674      endif
    646675      ierr=nf_get_att_text(nid,varid,"units",units)
    647       call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr)
     676      call def_var(nout,var(j),long_name,units,ndim,dim,varidout,ierr)
    648677
    649678     
    … …  
    653682!==============================================================================
    654683
    655 #ifdef NC_DOUBLE
    656       ierr = NF_GET_VAR_DOUBLE(nid,varid,var3d)
    657       miss=NF_GET_ATT_DOUBLE(nid,varid,"missing_value",missing)
    658       validr=NF_GET_ATT_DOUBLE(nid,varid,"valid_range",valid_range)
    659 #else
    660684      ierr = NF_GET_VAR_REAL(nid,varid,var3d)
    661685      miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
    662686      validr=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
    663 #endif
    664 
    665       ! If there is no missing value attribute in the input variable, create one
     687
     688     ! If there is no missing value attribute in the input variable, create one
    666689      if (miss.ne.NF_NOERR) then
    667690         missing = miss_lt
    … …  
    673696!==============================================================================
    674697
    675       do ilon=1,lonlen
     698       do ilon=1,lonlen
    676699!        Recast GCM local time (sol decimal value) at each longitude :
    677700         do it=1,timelen  ! loop time input file
    … …  
    688711               lt_outc(it)=intsol(it) + lt_out(ilon,ilat,it)/24.
    689712             endif
    690 !            If the computed local time exceeds the GCM input file time bounds,
    691 !            put a missing value in the local time variable
    692              if ((lt_outc(it).gt.lt_gcm(timelen_tot)).OR.(lt_outc(it).lt.lt_gcm(1))) then
    693                lt_out(ilon,ilat,it)=miss_lt
    694                lt_outc(it)=miss_lt
    695              endif
     713             if (nsol.eq.1) then
     714!              If 1 sol file (stats), use all data from that sol.
     715               if(lt_outc(it).gt.lt_gcm(timelen_tot)) lt_outc(it)=lt_outc(it)-1
     716               if(lt_outc(it).lt.lt_gcm(1)) lt_outc(it)=lt_outc(it)+1
     717             else
     718!              If the computed local time exceeds the GCM input file time bounds,
     719!              put a missing value in the local time variable
     720               if ((lt_outc(it).gt.lt_gcm(timelen_tot)).OR.(lt_outc(it).lt.lt_gcm(1))) then
     721                 lt_out(ilon,ilat,it)=miss_lt
     722                 lt_outc(it)=miss_lt
     723               endif
     724             end if
    696725           enddo ! local time
    697726
    … …  
    710739             enddo ! local time
    711740
    712            else if  (ndim==4) then
     741           else if (ndim==4) then
    713742             do ialt=1,altlen
    714743               do it=1,timelen ! loop time input file
    … …  
    725754               enddo ! local time   
    726755             enddo ! alt
    727            endif
     756           endif ! ndim
    728757
    729758         enddo ! lat
    730       end do ! lon
     759       end do ! lon
    731760
    732761
    … …  
    735764!==============================================================================
    736765
    737 #ifdef NC_DOUBLE
    738       ierr= NF_PUT_VARA_DOUBLE(nout,varidout,corner,edges,var3d_lt)
    739 #else
    740       ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3d_lt)
    741 #endif
    742 
    743       if (ierr.ne.NF_NOERR) then
     766       ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3d_lt)
     767
     768       if (ierr.ne.NF_NOERR) then
    744769         write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
    745770         stop ""
    746       endif
     771       endif
    747772
    748773      ! Write "missing_value" (and "valid_range" if there is one) attributes in output file
    749       call missing_value(nout,varidout,validr,miss,valid_range,missing)
     774       call missing_value(nout,varidout,validr,miss,valid_range,missing)
    750775
    751776      ! free var3d() array
    752       deallocate(var3d)
    753       deallocate(var3d_lt)
     777       deallocate(var3d)
     778       deallocate(var3d_lt)
    754779
    755780   enddo ! of do j=1,nbvar
    … …  
    780805   write(szastr,*) sza
    781806   if (trim(planetside).eq."morning") then
    782      title="Morning-side LTST at sza="//trim(adjustl(szastr))//"deg"
     807     long_name="Morning-side LTST at sza="//trim(adjustl(szastr))//"deg"
    783808   else !if trim(planetside).eq."evening"
    784      title="Evening-side LTST at sza="//trim(adjustl(szastr))//"deg"           
    785    endif
    786    call def_var(nout,tmpvar,title,units,3,dim,varidout,ierr)
    787 
    788 #ifdef NC_DOUBLE
    789    ierr= NF_PUT_VARA_DOUBLE(nout,varidout,corner,edges,lt_out)
    790 #else
     809     long_name="Evening-side LTST at sza="//trim(adjustl(szastr))//"deg"           
     810   endif
     811   call def_var(nout,tmpvar,long_name,units,3,dim,varidout,ierr)
     812
    791813   ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,lt_out)
    792 #endif
    793814   if (ierr.ne.NF_NOERR) then
    794815      write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
    … …  
    798819   validr=NF_NOERR+1
    799820   miss=NF_NOERR
    800    call missing_value(nout,varidout,validr,miss,(0,0),miss_lt)
     821   valid_range(:)=0.
     822   call missing_value(nout,varidout,validr,miss,valid_range,miss_lt)
    801823   write(*,*)"with missing value = ",miss_lt
    802824!==============================================================================
    … …  
    837859!******************************************************************************
    838860
    839 subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
     861subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,GCM_layers,&
    840862                     altlong_name,altunits,altpositive,&
    841                      nout,latdimout,londimout,altdimout,timedimout,timevarout)
     863                     nout,latdimout,londimout,altdimout,timedimout,&
     864                     layerdimout,timevarout)
    842865!==============================================================================
    843866! Purpose:
    … …  
    866889real, intent(in):: ctl(ctllen)
    867890! ctl(): controle
     891integer,intent(in) :: GCM_layers ! number of GCM layers
    868892character (len=*), intent(in) :: altlong_name
    869893! altlong_name(): [netcdf] altitude "long_name" attribute
    … …  
    882906integer, intent(out):: timedimout
    883907! timedimout: [netcdf] "Time"  ID
     908integer,intent(out) :: layerdimout
     909! layerdimout: [netcdf] "GCM_layers" ID
    884910integer, intent(out):: timevarout
    885911! timevarout: [netcdf] "Time" (considered as a variable) ID
    … …  
    943969endif
    944970
     971ierr = NF_DEF_DIM(nout, "GCM_layers", GCM_layers, layerdimout)
     972if (ierr.NE.NF_NOERR) then
     973   WRITE(*,*)'initiate: error failed to define dimension <GCM_layers>'
     974   write(*,*) NF_STRERROR(ierr)
     975   stop ""
     976endif
     977
    945978! End netcdf define mode
    946979ierr = NF_ENDDEF(nout)
    … …  
    950983   stop ""
    951984endif
    952 
    953985!==============================================================================
    954986! 3. Write "Time" and its attributes
    … …  
    960992ierr = NF_REDEF (nout)
    961993ierr = NF_PUT_ATT_TEXT(nout,timevarout,'comment',135,&
    962       "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")
     994"The true unit of the variable Time is the integer value &
     995 of sol at lon=0deg. A false unit is put to enable reading from Grads or Ferret")
    963996! End netcdf define mode
    964997ierr = NF_ENDDEF(nout)
    … …  
    9711004             (/latdimout/),nvarid,ierr)
    9721005
    973 #ifdef NC_DOUBLE
    974 ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lat)
    975 #else
    9761006ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
    977 #endif
    9781007if (ierr.NE.NF_NOERR) then
    9791008   WRITE(*,*)'initiate: error failed writing variable <latitude>'
    … …  
    9891018             (/londimout/),nvarid,ierr)
    9901019
    991 #ifdef NC_DOUBLE
    992 ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lon)
    993 #else
    9941020ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
    995 #endif
    9961021if (ierr.NE.NF_NOERR) then
    9971022   WRITE(*,*)'initiate: error failed writing variable <longitude>'
    … …  
    10071032ierr = NF_REDEF (nout)
    10081033
    1009 #ifdef NC_DOUBLE
    1010 ierr = NF_DEF_VAR (nout,"altitude",NF_DOUBLE,1,altdimout,nvarid)
    1011 #else
    10121034ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,altdimout,nvarid)
    1013 #endif
    10141035
    10151036ierr = NF_PUT_ATT_TEXT (nout,nvarid,'long_name',len_trim(adjustl(altlong_name)),adjustl(altlong_name))
    … …  
    10201041ierr = NF_ENDDEF(nout)
    10211042
    1022 #ifdef NC_DOUBLE
    1023 ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,alt)
    1024 #else
    1025 ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
    1026 #endif
     1043ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
    10271044if (ierr.NE.NF_NOERR) then
    10281045   WRITE(*,*)'initiate: error failed writing variable <altitude>'
    … …  
    10391056   ierr = NF_REDEF (nout)
    10401057
    1041 #ifdef NC_DOUBLE
    1042    ierr = NF_DEF_VAR (nout,"controle",NF_DOUBLE,1,ctldimout,nvarid)
    1043 #else
    10441058   ierr = NF_DEF_VAR (nout,"controle",NF_FLOAT,1,ctldimout,nvarid)
    1045 #endif
    1046 
    1047    ierr = NF_PUT_ATT_TEXT (nout,nvarid,"title",18,"Control parameters")
     1059
     1060   ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",18,"Control parameters")
    10481061
    10491062   ! End netcdf define mode
    10501063   ierr = NF_ENDDEF(nout)
    10511064
    1052 #ifdef NC_DOUBLE
    1053    ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,ctl)
    1054 #else
    10551065   ierr = NF_PUT_VAR_REAL (nout,nvarid,ctl)
    1056 #endif
    10571066   if (ierr.NE.NF_NOERR) then
    10581067      WRITE(*,*)'initiate: error failed writing variable <controle>'
    … …  
    10651074
    10661075!******************************************************************************
    1067 subroutine init2(infid,lonlen,latlen,altlen,altdimid, &
    1068                  outfid,londimout,latdimout,altdimout)
     1076subroutine init2(infid,lonlen,latlen,altlen,GCM_layers, &
     1077                 outfid,londimout,latdimout,altdimout, &
     1078                 layerdimout)
    10691079!==============================================================================
    10701080! Purpose:
    … …  
    10861096integer, intent(in) :: latlen ! # of grid points along latitude
    10871097integer, intent(in) :: altlen ! # of grid points along altitude
    1088 integer, intent(in) :: altdimid ! ID of altitude dimension
     1098integer, intent(in) :: GCM_layers ! # of GCM atmospheric layers
    10891099integer, intent(in) :: outfid ! NetCDF output file ID
    10901100integer, intent(in) :: londimout ! longitude dimension ID
    10911101integer, intent(in) :: latdimout ! latitude dimension ID
    10921102integer, intent(in) :: altdimout ! altitude dimension ID
     1103integer, intent(in) :: layerdimout ! GCM_layers dimension ID
    10931104!==============================================================================
    10941105! Local variables:
    … …  
    11001111integer :: tmpdimid ! temporary dimension ID
    11011112integer :: tmpvarid ! temporary variable ID
    1102 logical :: phis, aps_ok, bps_ok ! are "phisinit" "aps" "bps" available ?
     1113integer :: tmplen ! temporary variable length
     1114logical :: phis,hybrid ! are "phisinit" and "aps"/"bps" available ?
    11031115
    11041116
    … …  
    11081120
    11091121! hybrid coordinate aps
    1110   allocate(aps(altlen),stat=ierr)
     1122!write(*,*) "aps: altlen=",altlen," GCM_layers=",GCM_layers
     1123allocate(aps(GCM_layers),stat=ierr)
     1124if (ierr.ne.0) then
     1125  write(*,*) "init2: failed to allocate aps!"
     1126  stop
     1127endif
     1128ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
     1129if (ierr.ne.NF_NOERR) then
     1130  write(*,*) "Ooops. Failed to get aps ID. OK."
     1131  hybrid=.false.
     1132else
     1133  ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
     1134  hybrid=.true.
     1135  if (ierr.ne.NF_NOERR) then
     1136    stop "init2 Error: Failed reading aps"
     1137  endif
     1138
     1139  ! hybrid coordinate bps
     1140!  write(*,*) "bps: altlen=",altlen," GCM_layers=",GCM_layers
     1141  allocate(bps(GCM_layers),stat=ierr)
    11111142  if (ierr.ne.0) then
    1112     write(*,*) "init2: failed to allocate aps(altlen)"
     1143    write(*,*) "init2: failed to allocate bps!"
    11131144    stop
    11141145  endif
    1115 
    1116 ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
    1117 if (ierr.ne.NF_NOERR) then
    1118   write(*,*) "oops Failed to get aps ID. OK"
    1119   aps_ok=.false.
    1120 else
    1121   ! Check that aps() is the right size (it most likely won't be if
    1122   ! zrecast has be used to generate the input file)
    1123   ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
    1124   if (tmpdimid.ne.altdimid) then
    1125     write(*,*) "init2: wrong dimension size for aps(), skipping it ..."
    1126     aps_ok=.false.
    1127   else
    1128     ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
    1129     if (ierr.ne.NF_NOERR) then
    1130      stop "init2 error: Failed reading aps"
    1131     endif
    1132     aps_ok=.true.
    1133   endif
    1134 endif
    1135 
    1136 ! hybrid coordinate bps
    1137   allocate(bps(altlen),stat=ierr)
    1138   if (ierr.ne.0) then
    1139     write(*,*) "init2: failed to allocate bps(altlen)"
    1140     stop
    1141   endif
    1142 
    1143 ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
    1144 if (ierr.ne.NF_NOERR) then
    1145   write(*,*) "oops: Failed to get bps ID. OK"
    1146   bps_ok=.false.
    1147 else
    1148   ! Check that bps() is the right size
    1149   ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
    1150   if (tmpdimid.ne.altdimid) then
    1151     write(*,*) "init2: wrong dimension size for bps(), skipping it ..."
    1152     bps_ok=.false.
     1146  ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
     1147  if (ierr.ne.NF_NOERR) then
     1148    write(*,*) "Ooops. Failed to get bps ID. OK."
     1149    hybrid=.false.
    11531150  else
    11541151    ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
    … …  
    11561153      stop "init2 Error: Failed reading bps"
    11571154    endif
    1158     bps_ok=.true.
    11591155  endif
    11601156endif
    … …  
    11881184!==============================================================================
    11891185
    1190 IF (aps_ok) then
     1186IF (hybrid) then
    11911187  ! define aps
    11921188  call def_var(outfid,"aps","hybrid pressure at midlayers"," ",1,&
    1193                (/altdimout/),apsid,ierr)
     1189               (/layerdimout/),apsid,ierr)
    11941190  if (ierr.ne.NF_NOERR) then
    11951191    stop "Error: Failed to def_var aps"
    … …  
    11971193
    11981194  ! write aps
    1199 #ifdef NC_DOUBLE
    1200   ierr=NF_PUT_VAR_DOUBLE(outfid,apsid,aps)
    1201 #else
    12021195  ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
    1203 #endif
    12041196  if (ierr.ne.NF_NOERR) then
    12051197    stop "Error: Failed to write aps"
    12061198  endif
    1207 ENDIF ! of IF (aps_ok)
    1208 
    1209 IF (bps_ok) then
     1199 
    12101200  ! define bps
    12111201  call def_var(outfid,"bps","hybrid sigma at midlayers"," ",1,&
    1212                (/altdimout/),bpsid,ierr)
     1202               (/layerdimout/),bpsid,ierr)
    12131203  if (ierr.ne.NF_NOERR) then
    12141204    stop "Error: Failed to def_var bps"
    … …  
    12161206
    12171207  ! write bps
    1218 #ifdef NC_DOUBLE
    1219   ierr=NF_PUT_VAR_DOUBLE(outfid,bpsid,bps)
    1220 #else
    12211208  ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
    1222 #endif
    12231209  if (ierr.ne.NF_NOERR) then
    12241210    stop "Error: Failed to write bps"
    12251211  endif
    1226 ENDIF ! of IF (bps_ok)
     1212 
     1213  deallocate(aps)
     1214  deallocate(bps)
     1215ENDIF ! of IF (hybrid)
    12271216
    12281217!==============================================================================
    … …  
    12401229
    12411230  ! write phisinit
    1242 #ifdef NC_DOUBLE
    1243   ierr=NF_PUT_VAR_DOUBLE(outfid,phisinitid,phisinit)
    1244 #else
    12451231  ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
    1246 #endif
    12471232  if (ierr.ne.NF_NOERR) then
    12481233    stop "Error: Failed to write phisinit"
    12491234  endif
    12501235
     1236deallocate(phisinit)
    12511237ENDIF ! of IF (phis)
    12521238
    12531239
    1254 ! Cleanup
    1255 deallocate(aps)
    1256 deallocate(bps)
    1257 deallocate(phisinit)
    1258 
    12591240end subroutine init2
    12601241
    12611242!******************************************************************************
    1262 subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
     1243subroutine def_var(nid,name,long_name,units,nbdim,dim,nvarid,ierr)
    12631244!==============================================================================
    12641245! Purpose: Write a variable (i.e: add a variable to a dataset)
    1265 ! called "name"; along with its attributes "title", "units"...
     1246! called "name"; along with its attributes "long_name", "units"...
    12661247! to a file (following the NetCDF format)
    12671248!==============================================================================
    … …  
    12811262character (len=*), intent(in) :: name
    12821263! name(): [netcdf] variable's name
    1283 character (len=*), intent(in) :: title
    1284 ! title(): [netcdf] variable's "title" attribute
     1264character (len=*), intent(in) :: long_name
     1265! long_name(): [netcdf] variable's "long_name" attribute
    12851266character (len=*), intent(in) :: units
    12861267! unit(): [netcdf] variable's "units" attribute
    … …  
    12981279
    12991280! Insert the definition of the variable
    1300 #ifdef NC_DOUBLE
    1301 ierr=NF_DEF_VAR(nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
    1302 #else
    13031281ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
    1304 #endif
    13051282
    13061283! Write the attributes
    1307 ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",len_trim(adjustl(title)),adjustl(title))
     1284ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name",len_trim(adjustl(long_name)),adjustl(long_name))
    13081285ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
    13091286
    … …  
    13631340! Write valid_range() attribute
    13641341if (validr.eq.NF_NOERR) then
    1365 #ifdef NC_DOUBLE
    1366   ierr=NF_PUT_ATT_DOUBLE(nout,nvarid,'valid_range',NF_DOUBLE,2,valid_range)
    1367 #else
    13681342  ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
    1369 #endif
    13701343
    13711344  if (ierr.ne.NF_NOERR) then
    … …  
    13801353! Write "missing_value" attribute
    13811354if (miss.eq.NF_NOERR) then
    1382 #ifdef NC_DOUBLE
    1383   ierr= NF_PUT_ATT_DOUBLE(nout,nvarid,'missing_value',NF_DOUBLE,1,missing)
    1384 #else
    13851355  ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
    1386 #endif
    13871356
    13881357  if (ierr.NE.NF_NOERR) then

  • trunk/LMDZ.MARS/util/solzenangle.def

    r2289 r2434  
    18185) solar zenith angle (within [0;180[ deg)
    19190deg = zenith ; 90deg = terminator at the surface ; >90deg = terminator in altitude
    20 6) starttimeoffset :
    21     - for a diagfi/concat : sol value at the beginning of the run, wrt Ls=0
    22     - for a stats : sol value at the middle of the run, wrt Ls=0
    23 
     206) reference time
     21    - for a diagfi/concat :
     22        - write "n" to keep the date stored in the file (recommended)     
     23        - if needed:  sol value at the beginning of the run, wrt Ls=0
     24    - for a stats : Ls value at the middle of the run (degree)
    2425USE :
    2526solzenangle.e < solzenangle.def
Note: See TracChangeset for help on using the changeset viewer.