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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.