Ignore:
Timestamp:
Apr 28, 2020, 11:44:08 AM (5 years ago)
Author:
abierjon
Message:

Mars GCM:
Resolved Ticket #46 : it is now possible to concatenate concat files with a coherent
time axis. NB : the program puts now a higher importance in respecting the controle
variable of all the input files (when there is one).
Example : concat[diagfi2 (66sols from sol 61); diagfi4 (65sols from sol 193)]
Old output Time: [61.08,..127,|127.08,..192]; New output Time: [61.08,..127,|193.08,..258]

AB

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2147 r2303  
    1919!    (case of stats file) FF, November 2011
    2020! + read and write controle field, if available. TN, October 2013
     21! + possibility to concatenate concat files with a coherent
     22!   time axis in the output. AB, April 2020
    2123! ********************************************************
    2224
    … …  
    5153integer :: varid
    5254! varid: [netcdf] variable ID #
    53 integer :: memolatlen=0,memolonlen=0,memoaltlen=0
     55integer :: memolatlen=0,memolonlen=0,memoaltlen=0,memoctllen=0
    5456! memolatlen: # of elements of lat(), read from the first input file
    5557! memolonlen: # of elements of lon(), read from the first input file
    5658! memoaltlen: # of elements of alt(), read from the first input file
     59! memoctllen: # of elements of ctl(), read from the first input file
    5760real, dimension(:), allocatable:: lat,lon,alt,ctl,time
    5861! lat(): array, stores latitude coordinates
    … …  
    109112real, dimension(:,:,:,:), allocatable :: var3d
    110113! var3D(,,,): 4D array to store a field
     114real :: ctlsol
     115! ctlsol: starting sol value of the file, stored in the variable ctl()
    111116real :: memotime
    112117! memotime: (cumulative) time value, in martian days (sols)
    … …  
    121126!==============================================================================
    122127write(*,*)
    123 write(*,*) "-- Program written specifically for NetCDF 'diagfi' files from the martian GCM --"
     128write(*,*) "-- Program written specifically for NetCDF 'diagfi'/'concat' files from the martian GCM --"
    124129write(*,*)
    125130write(*,*) "which files do you want to use?"
    … …  
    154159!==============================================================================
    155160
     161write(*,*)
     162write(*,*) "opening "//trim(file(1))//"..."
    156163ierr = NF_OPEN(file(1),NF_NOWRITE,nid)
    157164if (ierr.NE.NF_NOERR) then
    … …  
    319326   ierr=NF_INQ_VARID(nid,"latitude",latvar)
    320327   if (ierr.NE.NF_NOERR) then
    321       write(*,*) 'ERROR: Field <latitude> is missing in file'//file(i)
     328      write(*,*) 'ERROR: Field <latitude> is missing in file '//file(i)
    322329      stop "" 
    323330   endif
    … …  
    328335   ierr=NF_INQ_VARID(nid,"longitude",lonvar)
    329336   if (ierr.NE.NF_NOERR) then
    330       write(*,*) 'ERROR: Field <longitude> is missing in file'//file(i)
     337      write(*,*) 'ERROR: Field <longitude> is missing in file '//file(i)
    331338      stop ""
    332339   endif
    … …  
    337344   ierr=NF_INQ_VARID(nid,"altitude",altvar)
    338345   if (ierr.NE.NF_NOERR) then
    339       write(*,*) 'ERROR: Field <altitude> is missing in file'//file(i)
     346      write(*,*) 'ERROR: Field <altitude> is missing in file '//file(i)
    340347      stop ""
    341348   endif
    … …  
    346353   ierr=NF_INQ_VARID(nid,"controle",ctlvar)
    347354   if (ierr.NE.NF_NOERR) then
    348       write(*,*) 'Field <controle> is missing in file'//file
     355      write(*,*) 'Field <controle> is missing in file '//file(i)
    349356      ctllen=0
    350357      !stop ""
    … …  
    372379!==============================================================================
    373380
     381   if (ctllen .ne. 0) then ! variable controle
     382      allocate(ctl(ctllen))
     383#ifdef NC_DOUBLE
     384      ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
     385#else
     386      ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
     387#endif
     388      ctlsol = ctl(4)
     389   endif
     390
    374391   if (i==1) then ! First call; initialize/allocate
    375392      memolatlen=latlen
    376393      memolonlen=lonlen
    377394      memoaltlen=altlen
     395      memoctllen=ctllen
    378396      allocate(lat(latlen))
    379397      allocate(lon(lonlen))
    380398      allocate(alt(altlen))
    381       allocate(ctl(ctllen))
    382399#ifdef NC_DOUBLE
    383400      ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat)
    384401      ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon)
    385402      ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt)
    386       if (ctllen .ne. -1) ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
    387403#else
    388404      ierr = NF_GET_VAR_REAL(nid,latvar,lat)
    389405      ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
    390406      ierr = NF_GET_VAR_REAL(nid,altvar,alt)
    391       if (ctllen .ne. -1) ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
    392 #endif
    393       if (ctllen .ne. -1) then
    394          if (modulo(int(memotime),669)/=modulo(int(ctl(4)),669)) then
    395            write(*,'(2(A,I4),A)') "WARNING: Starting day of the run is not ",&
    396                                 modulo(int(memotime),669)," but ",modulo(int(ctl(4)),669),"!!"
     407#endif
     408      if (ctllen .ne. 0) then
     409         if (modulo(int(memotime),669)/=modulo(int(ctlsol),669)) then
     410           write(*,*) "WARNING: Starting day of the run is not ",&
     411                                modulo(int(memotime),669)," but ",modulo(int(ctlsol),669),"!!"
    397412           write(*,*) "Starting day of the run has been corrected."
    398            memotime=float(modulo(int(ctl(4)),669)) + ctl(27)
    399            ctl(4) = 0.
    400            ctl(27) = 0.
     413!           memotime=float(modulo(int(ctl(4)),669)) + ctl(27)
     414           memotime=float(int(ctlsol)) + ctl(27)
     415           ctl(4) = 0.  ! values written in the output
     416           ctl(27) = 0. ! file controle variable (resp. day_ini and hour_ini)
    401417         endif
    402418      endif
    … …  
    416432           write(*,*) "ERROR: Not the same latitude axis"
    417433           stop ""
    418         else if (memolonlen/=lonlen) then
     434       else if (memolonlen/=lonlen) then
    419435           write(*,*) "ERROR: Not the same longitude axis"
    420436           stop ""
    421              else if (memoaltlen/=altlen) then
    422                 write(*,*) "ERROR: Not the same altitude axis"
    423                 stop ""
     437       else if (memoaltlen/=altlen) then
     438           write(*,*) "ERROR: Not the same altitude axis"
     439           stop ""
     440       else if (memoctllen/=ctllen) then
     441           write(*,*) "ERROR: Not the same controle axis"
     442           stop ""
    424443       endif
    425444   endif ! of if (i==1)
    … …  
    434453   ierr=NF_INQ_DIMID(nid,"Time",timedim)
    435454   if (ierr.NE.NF_NOERR) then
    436       write(*,*) 'ERROR: Dimension <Time> is missing in file'//file(i)
     455      write(*,*) 'ERROR: Dimension <Time> is missing in file '//file(i)
    437456      stop ""
    438457   endif
    439458   ierr=NF_INQ_VARID(nid,"Time",timevar)
    440459   if (ierr.NE.NF_NOERR) then
    441       write(*,*) 'ERROR: Field <Time> is missing in file'//file(i)
     460      write(*,*) 'ERROR: Field <Time> is missing in file '//file(i)
    442461      stop ""
    443462   endif
    … …  
    461480   write(*,*)
    462481   write(*,'(a3,1x,f6.1)') "Sol",memotime
    463 
    464    ! Add (memotime) offset and write "concatenated" time values to output file
    465    do k=reptime+1,reptime+timelen
    466       rep=rep+1
    467 #ifdef NC_DOUBLE
    468       ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,k,1,memotime+time(rep))
    469 #else
    470       ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,memotime+time(rep))
    471 #endif
    472    enddo
     482!   if (ctllen.ne.0) write(*,'(a30,1x,f6.1)') "In the current file's ctl: Sol",ctlsol
     483
     484   ! Add (memotime)/(ctlsol) offset and write "concatenated" time values to output file
     485   if (ctllen.ne.0) then
     486     do k=reptime+1,reptime+timelen
     487       rep=rep+1
     488       if ((time(rep)+ctlsol).le.memotime)then
     489         write(*,*) "ERROR : the time intervals of the files are not dissociated"
     490         stop ""
     491       endif
     492#ifdef NC_DOUBLE
     493       ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,k,1,ctlsol+time(rep))
     494#else
     495       ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,ctlsol+time(rep))
     496#endif
     497     enddo
     498
     499   else
     500     do k=reptime+1,reptime+timelen
     501        rep=rep+1
     502#ifdef NC_DOUBLE
     503        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,k,1,memotime+time(rep))
     504#else
     505        ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,memotime+time(rep))
     506#endif
     507     enddo
     508   endif ! if ctllen.ne.0
     509
    473510   ! Compute new time offset (for further concatenations)
    474511   memotime=memotime+time(timelen)
    … …  
    597634   deallocate(time)
    598635   reptime=reptime+timelen
     636
     637   ! Deallocate controle if needed
     638   if (ctllen.ne.0) deallocate(ctl)
    599639
    600640   ! Close input file
Note: See TracChangeset for help on using the changeset viewer.