Merge pull request #2477 from GEOS-ESM/feature/ygyu/aeronet_revise

Station sampler: add support to Global Historical Climatology Network Daily (GHCN-D)
GEOS-ESM · Dec 20, 2023 · 85d64fb · 85d64fb
2 parents b304b5a + ee6a9c8
commit 85d64fb
Show file tree

Hide file tree

Showing 4 changed files with 83 additions and 26 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- Station sampler: add support to Global Historical Climatology Network Daily (GHCN-D)
 - Add to trajectory sampler DEFINE_OBS_PLATFORM for reading multiple IODA files. To do this, we add union_platform function for observation.
 - New directory (`docs/tutorial/grid_comps/automatic_code_generator`) containing an example showing how to automatically generate the source code using the `MAPL_GridCompSpecs_ACG.py` tool.
 - Added/modified a few _ASSERT calls in ExtData, to better explain what is wrong in .yaml file

diff --git a/gridcomps/History/MAPL_HistoryCollection.F90 b/gridcomps/History/MAPL_HistoryCollection.F90
@@ -104,6 +104,7 @@ module MAPL_HistoryCollectionMod
      type(GriddedIOItemVector)          :: items
      character(len=ESMF_MAXSTR)         :: currentFile
      character(len=ESMF_MAXPATHLEN)     :: stationIdFile
+     integer                            :: stationSkipLine
      logical                            :: splitField
      logical                            :: regex
      logical                            :: timeseries_output = .false.

diff --git a/gridcomps/History/MAPL_HistoryGridComp.F90 b/gridcomps/History/MAPL_HistoryGridComp.F90
@@ -886,6 +886,8 @@ subroutine Initialize ( gc, import, dumexport, clock, rc )
             label=trim(string) // 'sampler_spec:', _RC)
        call ESMF_ConfigGetAttribute(cfg, value=list(n)%stationIdFile, default="", &
             label=trim(string) // 'station_id_file:', _RC)
+       call ESMF_ConfigGetAttribute(cfg, value=list(n)%stationSkipLine, default=0, &
+            label=trim(string) // 'station_skip_line:', _RC)
 
 ! Get an optional file containing a 1-D track for the output
        call ESMF_ConfigGetDim(cfg, nline, ncol,  label=trim(string)//'obs_files:', rc=rc)  ! here donot check rc on purpose
@@ -2399,7 +2401,7 @@ subroutine Initialize ( gc, import, dumexport, clock, rc )
              list(n)%trajectory = HistoryTrajectory(cfg,string,clock,_RC)
              call list(n)%trajectory%initialize(items=list(n)%items,bundle=list(n)%bundle,timeinfo=list(n)%timeInfo,vdata=list(n)%vdata,_RC)
           elseif (list(n)%sampler_spec == 'station') then
-             list(n)%station_sampler = StationSampler (trim(list(n)%stationIdFile),_RC)
+             list(n)%station_sampler = StationSampler (trim(list(n)%stationIdFile), nskip_line=list(n)%stationSkipLine, _RC)
              call list(n)%station_sampler%add_metadata_route_handle(list(n)%bundle,list(n)%timeInfo,vdata=list(n)%vdata,_RC)
           else
              global_attributes = list(n)%global_atts%define_collection_attributes(_RC)
@@ -3475,22 +3477,10 @@ subroutine Run ( gc, import, export, clock, rc )
          read(DateStamp( 1: 8),'(i8.8)') nymd
          read(DateStamp(10:15),'(i6.6)') nhms
 
-!         write(6,'(a)') 'bf fill_grads_template'
-!         write(6,'(10a)') 'filename(n), fntmpl=', trim(filename(n)), trim(fntmpl)
-!         write(6,'(10a)') 'trim(INTSTATE%expid)', trim(INTSTATE%expid)
-!         write(6,'(2x,a,10i20)') 'nymd, nhms', nymd, nhms
-
-
          call fill_grads_template ( filename(n), fntmpl, &
               experiment_id=trim(INTSTATE%expid), &
               nymd=nymd, nhms=nhms, _RC ) ! here is where we get the actual filename of file we will write
 
-!         write(6,'(a)') 'af fill_grads_template'
-!         write(6,'(a)') 'filename(n), fntmpl=', trim(filename(n)), trim(fntmpl)
-!         write(6,'(10a)') 'trim(INTSTATE%expid)', trim(INTSTATE%expid)
-!         write(6,'(2x,a,10i20)') 'nymd, nhms', nymd, nhms
-
-
          if(list(n)%monthly .and. list(n)%partial) then
             filename(n)=trim(filename(n)) // '-partial'
             list(n)%currentFile = filename(n)
@@ -3610,7 +3600,7 @@ subroutine Run ( gc, import, export, clock, rc )
             state_out = INTSTATE%GIM(n)
          end if
 
-         if (.not.list(n)%timeseries_output) then
+         if (.not.list(n)%timeseries_output .AND. list(n)%sampler_spec /= 'station') then
             IOTYPE: if (list(n)%unit < 0) then    ! CFIO
                call list(n)%mGriddedIO%bundlepost(list(n)%currentFile,oClients=o_Clients,_RC)
             else
@@ -3634,6 +3624,11 @@ subroutine Run ( gc, import, export, clock, rc )
             end if IOTYPE
          end if
 
+         if (list(n)%sampler_spec == 'station') then
+            call ESMF_ClockGet(clock,currTime=current_time,_RC)
+            call list(n)%station_sampler%append_file(current_time,_RC)
+         endif
+
       endif OUTTIME
 
       if( NewSeg(n) .and. list(n)%unit /= 0 .and. list(n)%duration /= 0 ) then
@@ -3698,10 +3693,6 @@ subroutine Run ( gc, import, export, clock, rc )
             call list(n)%trajectory%destroy_rh_regen_LS (_RC)
          end if
       end if
-      if (list(n)%sampler_spec == 'station') then
-         call ESMF_ClockGet(clock,currTime=current_time,_RC)
-         call list(n)%station_sampler%append_file(current_time,_RC)
-      endif
 
       if( Writing(n) .and. list(n)%unit < 0) then
 

diff --git a/gridcomps/History/MAPL_StationSamplerMod.F90 b/gridcomps/History/MAPL_StationSamplerMod.F90
@@ -24,6 +24,7 @@ module StationSamplerMod
      integer                  :: nstation
      integer, allocatable :: station_id(:)
      character(len=ESMF_MAXSTR), allocatable :: station_name(:)
+     character(len=ESMF_MAXSTR), allocatable :: station_fullname(:)
      real(kind=REAL64), allocatable :: lons(:)
      real(kind=REAL64), allocatable :: lats(:)
      real(kind=REAL64), allocatable :: elevs(:)
@@ -49,34 +50,51 @@ module StationSamplerMod
 
 contains
 
-  function new_StationSampler_readfile (filename,rc) result(sampler)
+  function new_StationSampler_readfile (filename,nskip_line, rc) result(sampler)
     use pflogger, only             :  Logger, logging
     implicit none
     type(StationSampler)           :: sampler
     character(len=*), intent(in)   :: filename
+    integer, optional, intent(in) :: nskip_line
     integer, optional, intent(out) :: rc
 
     integer :: unit, ios, nstation, status
-    integer :: i, ncount
+    integer :: i, j, k, ncount
     logical :: con1, con2
     character (len=1)     :: CH1
     character (len=5)     :: seq
-    character (len=100)   :: line
+    character (len=100)   :: line, line2
+    integer               :: nskip
     type(Logger), pointer :: lgr
 
     !__ 1. read from station_id_file: static
     !      plain text format:
-    !      [name,lat,lon,elev] or [id,name,lat,lon,elev]
+    !      ["name,lat,lon,elev"] or ["id,name,lat,lon,elev"]
+    !      ["name_short lat lon elev name_full"]
     !
+
     open(newunit=unit, file=trim(filename), form='formatted', &
          access='sequential', status='old', _IOSTAT)
     ios=0
     nstation=0
+    nskip=0
+    if (present(nskip_line)) then
+       nskip=nskip_line
+    end if
+    if (nskip>0) then
+       do i=1, nskip
+          read(unit, *)
+       end do
+    end if
     read(unit, '(a100)', IOSTAT=ios) line
     call count_substring(line, ',', ncount)
-    con1= ncount.GE.3 .AND. ncount.LE.4
+    con1= (ncount>=2 .AND. ncount<=4).OR.(ncount==0)
     _ASSERT(con1, 'string sequence in Aeronet file not supported')
-    if (ncount==3) then
+    if (ncount==0) then
+       seq='AFFFA'
+    elseif (ncount==2) then
+       seq='AFF'
+    elseif (ncount==3) then
        seq='AFFF'
     elseif (ncount==4) then
        CH1=line(1:1)
@@ -94,6 +112,11 @@ function new_StationSampler_readfile (filename,rc) result(sampler)
     end if
 
     rewind(unit)
+    if (nskip>0) then
+       do i=1, nskip
+          read(unit, *)
+       end do
+    end if
     ios=0
     do while (ios==0)
        read(unit, '(a100)', IOSTAT=ios) line
@@ -102,10 +125,17 @@ function new_StationSampler_readfile (filename,rc) result(sampler)
     sampler%nstation=nstation
     allocate(sampler%station_id(nstation))
     allocate(sampler%station_name(nstation))
+    allocate(sampler%station_fullname(nstation))
     allocate(sampler%lons(nstation))
     allocate(sampler%lats(nstation))
     allocate(sampler%elevs(nstation))
+
     rewind(unit)
+    if (nskip>0) then
+       do i=1, nskip
+          read(unit, *)
+       end do
+    end if
     do i=1, nstation
        if(seq=='IAFFF') then
           read(unit, *) &
@@ -119,12 +149,41 @@ function new_StationSampler_readfile (filename,rc) result(sampler)
                sampler%station_id(i), &
                sampler%lats(i), &
                sampler%lons(i)
-       elseif(trim(seq)=='AFFF') then
+       elseif(trim(seq)=='AFF' .OR. trim(seq)=='AFFF') then
           read(unit, *) &
                sampler%station_name(i), &
                sampler%lats(i), &
                sampler%lons(i)
-               sampler%station_id(i)=i
+          sampler%station_id(i)=i
+       elseif(trim(seq)=='AFFFA') then
+       ! Ex: 'ZI000067991 -22.2170   30.0000  457.0 BEITBRIDGE 67991'
+          read(unit, *) &
+               sampler%station_name(i), &
+               sampler%lats(i), &
+               sampler%lons(i)
+          sampler%station_id(i)=i
+          backspace(unit)
+          read(unit, '(a100)', IOSTAT=ios) line
+          j=index(line, '.', BACK=.true.)
+          line2=line(j+1:)
+          k=len(line2)
+          line=''
+          do j=1, k
+             CH1=line2(j:j)
+             con1= (CH1>='a'.AND.CH1<='z').OR.(CH1>='A'.AND.CH1<='Z')
+             if (con1) exit
+          enddo
+          read(line2(j:k), '(a100)') line
+          line2=trim(line)
+          k=len(line2)
+          line=''
+          do j=1, k
+             CH1=line2(j:j)
+             con1= (CH1>='0' .AND. CH1<='9')
+             if (con1) exit
+          enddo
+          if (j>k) j=k
+          sampler%station_fullname(i) = trim(line2(1:j-1))
        end if
     end do
     close(unit)
@@ -149,6 +208,7 @@ function new_StationSampler_readfile (filename,rc) result(sampler)
     _RETURN(_SUCCESS)
   end function new_StationSampler_readfile
 
+
   subroutine add_metadata_route_handle (this,bundle,timeInfo,vdata,rc)
     class(StationSampler),  intent(inout)       :: this
     type(ESMF_FieldBundle), intent(in)          :: bundle
@@ -206,6 +266,7 @@ subroutine add_metadata_route_handle (this,bundle,timeInfo,vdata,rc)
     v = Variable(type=pFIO_INT32, dimensions='station_index')
     call this%fmd%add_variable('station_id',v)
 
+
     !__ 2. filemetadata: extract field from bundle, add_variable
     !
     call ESMF_FieldBundleGet(bundle, fieldCount=fieldCount, _RC)
@@ -244,14 +305,17 @@ subroutine add_metadata_route_handle (this,bundle,timeInfo,vdata,rc)
     end do
     deallocate (fieldNameList)
 
+
     !__ 3. locstream route handle
     !
     call ESMF_FieldBundleGet(bundle,grid=grid,_RC)
     this%regridder = LocStreamRegridder(grid,this%esmf_ls,_RC)
 
+
     _RETURN(_SUCCESS)
   end subroutine add_metadata_route_handle
 
+
   subroutine append_file(this,current_time,rc)
     class(StationSampler), intent(inout) :: this
     type(ESMF_Time), intent(in) :: current_time