Ignore:
Timestamp:
Oct 9, 2023, 9:05:12 AM (14 months ago)
Author:
emillour
Message:

Mars PCM utility:
Improve concatnc's handling of XIOS output files (handle change in some
variable names, e.g. "area" or "phisfi").
EM

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.MARS/util/concatnc.F90

    • Property svn:executable deleted
    r2590 r3077  
    3030character (len=100), dimension(1000) :: file
    3131! file(): input file(s) names(s)
    32 character (len=30), dimension(16) :: notconcat
    33 ! notconcat(): names of the (16) variables that won't be concatenated
     32character (len=30), dimension(24) :: notconcat
     33! notconcat(): names of the (24) variables that won't be concatenated
    3434character (len=100), dimension(:), allocatable :: var
    3535! var(): name(s) of variable(s) that will be concatenated
     
    131131! valid_range(2): [netcdf] interval in which a value is considered valid
    132132real :: year_day = 669.
    133 integer :: XIOS =0
    134133
    135134!==============================================================================
     
    220219notconcat(15)='phisinit'
    221220notconcat(16)='soildepth'
    222 
     221! for XIOS files:
     222notconcat(17)='lat'
     223notconcat(18)='lon'
     224notconcat(19)='time_instant'
     225notconcat(20)='time_instant_bounds'
     226notconcat(21)='time_counter'
     227notconcat(22)='time_counter_bounds'
     228notconcat(23)='area'
     229notconcat(24)='phisfi'
    223230
    224231!==============================================================================
     
    314321      ierr = NF_OPEN(file(i),NF_NOWRITE,nid)
    315322      if (ierr.NE.NF_NOERR) then
    316          write(*,*) 'ERROR: Pb opening file '//file(i)
     323         write(*,*) 'ERROR: Pb opening file '//trim(file(i))
    317324         write(*,*) NF_STRERROR(ierr)
    318325         stop
     
    329336   ierr=NF_INQ_VARID(nid,"latitude",latvar)
    330337   if (ierr.NE.NF_NOERR) then
     338     !failed to find "latitude"; look for "lat" instead
    331339     ierr=NF_INQ_DIMID(nid,"lat",latdim)
    332340     ierr=NF_INQ_VARID(nid,"lat",latvar)
    333      if (ierr.EQ.NF_NOERR) then
    334         write(*,*) 'XIOS concat spotted'
    335         XIOS=1
    336      endif
    337341   endif 
    338342   if (ierr.NE.NF_NOERR) then
    339       write(*,*) 'ERROR: Field <latitude> or <lat> is missing in file '//file(i)
    340       stop 
     343     write(*,*) 'ERROR: Neither <latitude> nor <lat> in file '//trim(file(i))
     344     stop 
    341345   endif
    342346   ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
    343347!  write(*,*) "latlen: ",latlen
    344348
    345    if (XIOS.EQ.0) then
    346        ierr=NF_INQ_DIMID(nid,"longitude",londim)
    347        ierr=NF_INQ_VARID(nid,"longitude",lonvar)
    348        if (ierr.NE.NF_NOERR) then
    349           write(*,*) 'ERROR: Field <longitude> is missing in file '//file(i)
    350           stop
    351        endif
    352        ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
    353        !  write(*,*) "lonlen: ",lonlen
    354    else
    355        ierr=NF_INQ_DIMID(nid,"lon",londim)
    356        ierr=NF_INQ_VARID(nid,"lon",lonvar)
    357        if (ierr.NE.NF_NOERR) then
    358           write(*,*) 'ERROR: Field <lon> is missing in file '//file(i)
    359           stop
    360        endif
    361        ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
    362        !  write(*,*) "lonlen: ",lonlen
     349   ierr=NF_INQ_DIMID(nid,"longitude",londim)
     350   ierr=NF_INQ_VARID(nid,"longitude",lonvar)
     351   if (ierr.NE.NF_NOERR) then
     352     ! failed to find "longitude"; look for "lon" instead
     353     ierr=NF_INQ_DIMID(nid,"lon",londim)
     354     ierr=NF_INQ_VARID(nid,"lon",lonvar)
    363355   endif
     356   if (ierr.NE.NF_NOERR) then
     357     write(*,*) 'ERROR: Neither <longitude> nor <lon> in file '//trim(file(i))
     358     stop
     359   endif
     360   ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
     361   !  write(*,*) "lonlen: ",lonlen
    364362
    365363   ierr=NF_INQ_DIMID(nid,"altitude",altdim)
    366364   ierr=NF_INQ_VARID(nid,"altitude",altvar)
    367365   if (ierr.NE.NF_NOERR) then
    368       write(*,*) 'ERROR: Field <altitude> is missing in file '//file(i)
     366      write(*,*) 'ERROR: Field <altitude> is missing in file '//trim(file(i))
    369367      stop
    370368   endif
     
    388386! 2.2.1 Check presence and read "controle" dimension & variable from input file
    389387!==============================================================================
    390    if (XIOS.EQ.0) then
    391         ierr=NF_INQ_DIMID(nid,"index",ctldim)
    392         if (ierr.NE.NF_NOERR) then
    393           write(*,*) 'Dimension <index> is missing in file '//file(i)
    394           write(*,*) "The program continues..."
    395           ctllen=0
    396           !stop
    397        else
    398           ierr=NF_INQ_VARID(nid,"controle",ctlvar)
    399           if (ierr.NE.NF_NOERR) then
    400             write(*,*) 'Field <controle> is missing in file '//file(i)
    401             write(*,*) "The program continues..."
    402             ctllen=0
    403             !stop
    404           else
    405             ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
    406           endif
    407        endif
    408    else
     388
     389   ctllen=-666 ! dummy initialization
     390   ierr=NF_INQ_DIMID(nid,"index",ctldim)
     391   if (ierr.NE.NF_NOERR) then
     392     ctllen=0
     393     ! no "index", look for "controle_axe" instead
    409394     ierr=NF_INQ_DIMID(nid,"controle_axe",ctldim)
    410395     if (ierr.NE.NF_NOERR) then
    411        write(*,*) 'Dimension <controle_axe> is missing in file '//file(i)
    412        write(*,*) "The program continues..."
     396       write(*,*) 'Neither <index> nor <controle_axe> in file '//trim(file(i))
     397       write(*,*) "The program continues without <controle>..."
    413398       ctllen=0
    414        !stop
    415399     else
    416        ierr=NF_INQ_VARID(nid,"controle",ctlvar)
    417        if (ierr.NE.NF_NOERR) then
    418           write(*,*) 'Field <controle> is missing in file '//file(i)
    419            write(*,*) "The program continues..."
    420           ctllen=0
    421           !stop
    422         else
    423           ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
    424        endif
     400       ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
     401     endif
     402   else
     403     ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
     404   endif
     405
     406   ! allocate ctl (even if ctllen==0; required to not trigger errors in debug
     407   !  mode when ctl is passed as an argument to routines)
     408   allocate(ctl(ctllen))
     409
     410   if (ctllen .ne. 0) then ! variable controle available
     411     ierr=NF_INQ_VARID(nid,"controle",ctlvar)
     412     if (ierr.NE.NF_NOERR) then
     413       write(*,*) 'Field <controle> is missing in file '//trim(file(i))
     414       write(*,*) "The program continues without <controle>..."
     415       ctllen=0 ! seting that implies there is no "controle" later on
     416     else
     417      ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
     418      ctlsol = ctl(4)
    425419     endif
    426420   endif
    427421
    428    if (ctllen .ne. 0) then ! variable controle
    429       allocate(ctl(ctllen))
    430       ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
    431       ctlsol = ctl(4)
    432    endif
    433 
    434422!==============================================================================
    435423! 2.3 Read "Time" dimension & variable from input file
    436424!==============================================================================
    437    if (XIOS.EQ.0) then
     425
    438426   ierr=NF_INQ_DIMID(nid,"Time",timedim)
    439427   if (ierr.NE.NF_NOERR) then
    440       write(*,*) 'ERROR: Dimension <Time> is missing in file '//file(i)
    441       stop
     428     ! Time not found, look for "time_counter" instead
     429     ierr=NF_INQ_DIMID(nid,"time_counter",timedim)
     430     if (ierr.NE.NF_NOERR) then
     431       write(*,*) 'ERROR: Neither <Time> nor <time_counter> dim in file '//trim(file(i))
     432       stop
     433     endif
    442434   endif
     435
     436   ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
     437!  write(*,*) "timelen: ",timelen
     438
    443439   ierr=NF_INQ_VARID(nid,"Time",timevar)
    444440   if (ierr.NE.NF_NOERR) then
    445       write(*,*) 'ERROR: Field <Time> is missing in file '//file(i)
    446       stop
    447    endif
    448    ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
    449 !  write(*,*) "timelen: ",timelen
    450    else
    451    ierr=NF_INQ_DIMID(nid,"time_counter",timedim)
    452    if (ierr.NE.NF_NOERR) then
    453       write(*,*) 'ERROR: Dimension <time_counter> is missing in file '//file(i)
    454       stop
    455    endif
    456    ierr=NF_INQ_VARID(nid,"time_counter",timevar)
    457    if (ierr.NE.NF_NOERR) then
    458       write(*,*) 'ERROR: Field <time_counter> is missing in file '//file(i)
    459       stop
    460    endif
    461    ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
    462 !  write(*,*) "timelen: ",timelen
     441     ! "Time" not found, look for "time_counter" instead
     442     ierr=NF_INQ_VARID(nid,"time_counter",timevar)
     443     if (ierr.NE.NF_NOERR) then
     444       write(*,*) 'ERROR: Neither <Time> nor <time_counter> in file '//trim(file(i))
     445       stop
     446     endif
    463447   endif
    464448
     
    568552!         correction if the file to concatenate seems "before" previous file
    569553!         (for instance to concatenate diagfi from the previous year at the enf of a year)
    570           print *, "time(rep)", time(rep), rep
    571           print *, "ctlsol", ctlsol
    572           print *, "starttimeoffset", starttimeoffset
    573           print *, "output_time", output_time
    574           print *, "previous_last_output_time", previous_last_output_time
     554          if (rep.eq.1) then ! only print these checks once
     555            print *, "time(rep=1)", time(rep)
     556            print *, "ctlsol", ctlsol
     557            print *, "starttimeoffset", starttimeoffset
     558            print *, "output_time", output_time
     559            print *, "previous_last_output_time", previous_last_output_time
     560            write(*,*)
     561          endif
    575562          do while (output_time.lt.previous_last_output_time) 
    576563                   output_time = output_time + year_day
     
    717704   reptime=reptime+timelen
    718705
    719    ! Deallocate controle if needed
    720    if (ctllen.ne.0) deallocate(ctl)
     706   ! Deallocate controle (is allocated and possibly read for each input file)
     707   deallocate(ctl)
    721708
    722709   ! Close input file
     
    10541041if (ierr.ne.NF_NOERR) then
    10551042  write(*,*)"init2 warning: Failed to get aire ID."
    1056   area = .false.
     1043  ! try "area" instead of "aire"
     1044  ierr=NF_INQ_VARID(infid,"area",tmpvarid)
     1045  if (ierr.ne.NF_NOERR) then
     1046    write(*,*)"init2 warning: Failed to get area ID."
     1047    area = .false.
     1048  else
     1049    ierr=NF_GET_VAR_REAL(infid,tmpvarid,aire)
     1050    if (ierr.ne.NF_NOERR) then
     1051      write(*,*) "init2 Error: Failed reading area"
     1052      stop
     1053    endif
     1054    area = .true.   
     1055  endif
    10571056else
    10581057  ierr=NF_GET_VAR_REAL(infid,tmpvarid,aire)
     
    10731072if (ierr.ne.NF_NOERR) then
    10741073  write(*,*)"init2 warning: Failed to get phisinit ID."
    1075   phis = .false.
     1074  ! try "phisfi" instead of "phisinit"
     1075  ierr=NF_INQ_VARID(infid,"phisfi",tmpvarid)
     1076  if (ierr.ne.NF_NOERR) then
     1077    write(*,*)"init2 warning: Failed to get phisfi ID."
     1078    phis = .false.
     1079  else
     1080    ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
     1081    if (ierr.ne.NF_NOERR) then
     1082      write(*,*) "init2 Error: Failed reading phisfi"
     1083      stop
     1084    endif
     1085  endif
     1086  phis = .true. 
    10761087else
    10771088  ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
Note: See TracChangeset for help on using the changeset viewer.