diff --git a/.circleci/config.yml b/.circleci/config.yml index 0cd5bd5971a5..c1d9deaf44b9 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -16,12 +16,12 @@ parameters: # Anchors to prevent forgetting to update a version os_version: &os_version ubuntu20 -baselibs_version: &baselibs_version v7.14.0 +baselibs_version: &baselibs_version v7.17.0 bcs_version: &bcs_version v11.3.0 tag_build_arg_name: &tag_build_arg_name maplversion orbs: - ci: geos-esm/circleci-tools@1 + ci: geos-esm/circleci-tools@2 workflows: build-and-test: @@ -221,7 +221,7 @@ workflows: when: equal: [ "release", << pipeline.parameters.GHA_Event >> ] jobs: - - ci/publish-docker: + - ci/publish_docker: filters: tags: only: /^v.*$/ @@ -238,7 +238,7 @@ workflows: compiler_version: 2022.1.0 image_name: geos-env tag_build_arg_name: *tag_build_arg_name - - ci/publish-docker: + - ci/publish_docker: filters: tags: only: /^v.*$/ diff --git a/.github/workflows/validate_yaml_files.yml b/.github/workflows/validate_yaml_files.yml index 8dc81d4f2b54..449db6e674da 100644 --- a/.github/workflows/validate_yaml_files.yml +++ b/.github/workflows/validate_yaml_files.yml @@ -24,7 +24,7 @@ jobs: format: colored config_file: .yamllint.yml - - uses: actions/upload-artifact@v3 + - uses: actions/upload-artifact@v4 if: always() with: name: yamllint-logfile diff --git a/.github/workflows/workflow.yml b/.github/workflows/workflow.yml index 467e64ea8ce6..8d77a47ab3b9 100644 --- a/.github/workflows/workflow.yml +++ b/.github/workflows/workflow.yml @@ -17,7 +17,7 @@ jobs: name: Build and Test MAPL GNU runs-on: ubuntu-latest container: - image: gmao/ubuntu20-geos-env-mkl:v7.14.0-openmpi_4.1.4-gcc_12.1.0 + image: gmao/ubuntu20-geos-env-mkl:v7.17.0-openmpi_4.1.4-gcc_12.1.0 # Per https://github.com/actions/virtual-environments/issues/1445#issuecomment-713861495 # It seems like we might not need secrets on GitHub Actions which is good for forked # pull requests @@ -77,7 +77,7 @@ jobs: name: Build and Test MAPL Intel runs-on: ubuntu-latest container: - image: gmao/ubuntu20-geos-env:v7.14.0-intelmpi_2021.6.0-intel_2022.1.0 + image: gmao/ubuntu20-geos-env:v7.17.0-intelmpi_2021.6.0-intel_2022.1.0 # Per https://github.com/actions/virtual-environments/issues/1445#issuecomment-713861495 # It seems like we might not need secrets on GitHub Actions which is good for forked # pull requests diff --git a/CHANGELOG.md b/CHANGELOG.md index 46169ef62b46..c8e980c155d3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,6 +17,34 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Deprecated +## [2.43.0] - 2023-12-21 + +### 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 + +### Changed + +- Change the verification of the grid in MAPL_GetGlobalHorzIJIndex to avoid collective call +- Swath grid step 1: allow for destroying and regenerating swath grid and regenerating regridder route handle, and creating + allocatable metadata in griddedIO. Modifications are made to GriddedIO.F90, MAPL_AbstractRegridder.F90, and MAPL_EsmfRegridder.F90. +- Swath grid step 2: add control keywords for swath grid. Allow for filename template with DOY. Allow for missing obs files. User needs to specify index_name_lon/lat, var_name_lon/lat/time, tunit, obs_file_begin/end/interval, Epoch and Epoch_init. +- Update CI to Baselibs 7.17.0 (for future MAPL3 work) and the BCs v11.3.0 (to fix coupled run) +- Update `components.yaml` + - ESMA_env v4.24.0 (Baselibs 7.17.0) +- Update CI to use circleci-tools v2 +- Changed the Python MAPL `__init__.py` file to restore behavior from pre-Python3 transition where we did `from foo import *`. Also fix up other Python2 code to Python3. + +### Fixed + +- Fixed bug broken multi-step file output in History under certain template conditions +- [#2433] Implemented workarounds for gfortran-13 +- Missing TARGET in GriddedIO - exposed runtime error when using NAG + debug. +- Allow ExtData2G to be built as SHARED or STATIC + ## [2.42.4] - 2023-12-10 ### Changed diff --git a/CMakeLists.txt b/CMakeLists.txt index c04500255b85..371040a11279 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,7 +8,7 @@ endif () project ( MAPL - VERSION 2.42.4 + VERSION 2.43.0 LANGUAGES Fortran CXX C) # Note - CXX is required for ESMF # Set the possible values of build type for cmake-gui @@ -207,6 +207,11 @@ endif () add_definitions(-Dsys${CMAKE_SYSTEM_NAME}) +# Support for automated code generation +include(mapl_acg) +include(mapl_create_stub_component) +add_subdirectory (Apps) + # Special case - MAPL_cfio is built twice with two different precisions. add_subdirectory (MAPL_cfio MAPL_cfio_r4) add_subdirectory (MAPL_cfio MAPL_cfio_r8) @@ -245,12 +250,6 @@ if (PFUNIT_FOUND) add_subdirectory (pfunit EXCLUDE_FROM_ALL) endif () - -# Support for automated code generation -include(mapl_acg) -include(mapl_create_stub_component) -add_subdirectory (Apps) - add_subdirectory (Tests) # @env will exist here if MAPL is built as itself but not as part of, say, GEOSgcm diff --git a/Python/MAPL/Date.py b/Python/MAPL/Date.py index 60e73a5ec4b5..dd53820d56d6 100644 --- a/Python/MAPL/Date.py +++ b/Python/MAPL/Date.py @@ -149,7 +149,7 @@ def copy(self): def __iter__(self): return self - def next(self): + def __next__(self): #Last day of the month. if self.day == NumberDaysMonth(self.month, self.year): self.day = 1 @@ -230,7 +230,7 @@ def __add__(self, n): #Convert back to date format. return DateFromJDNumber(temp) else: - raise TypeError, "%s is not an integer." % str(n) + raise TypeError("%s is not an integer." % str(n)) def __sub__(self, date): """Returns the (signed) difference of days between the dates.""" @@ -253,7 +253,7 @@ def __sub__(self, date): ret += NumberDaysYear(year) return ret else: - raise TypeError, "%s is neither an integer nor a Date." % str(date) + raise TypeError("%s is neither an integer nor a Date." % str(date)) #Adding an integer is "commutative". def __radd__(self, n): @@ -283,7 +283,7 @@ def ToCOMTime(self): def DateFromJDNumber(n): """Returns a date corresponding to the given Julian day number.""" if not isinstance(n, int): - raise TypeError, "%s is not an integer." % str(n) + raise TypeError("%s is not an integer." % str(n)) a = n + 32044 b = (4*a + 3)//146097 @@ -320,21 +320,21 @@ def strpdate(s): temp = Date() curr_month = temp.month while temp.month == curr_month: - print temp - temp.next() + print(temp) + next(temp) - print "\n" + print("\n") #How many days until the end of the year? temp = Date() temp.day, temp.month = 1, 1 curr_year = temp.year while temp.year == curr_year: - print "%s is %d days away from the end of the year." % (str(temp), - temp.DaysToEndYear()) + print("%s is %d days away from the end of the year." % (str(temp), + temp.DaysToEndYear())) temp += NumberDaysMonth(temp.month) - print "\n" + print("\n") #Playing with __sub__. temp = Date() @@ -344,11 +344,11 @@ def strpdate(s): temp_list.append(temp) temp += NumberDaysMonth(temp.month) for elem in temp_list: - print "%s differs %d days from current date: %s" % (str(elem), + print("%s differs %d days from current date: %s" % (str(elem), elem - Date(), - str(Date())) + str(Date()))) - print "\n" + print("\n") #Swapping arguments works? - print 23 + Date() + print(23 + Date()) diff --git a/Python/MAPL/__init__.py b/Python/MAPL/__init__.py index ad6e2d0f4afa..db14e158682e 100644 --- a/Python/MAPL/__init__.py +++ b/Python/MAPL/__init__.py @@ -45,3 +45,11 @@ """ __version__ = "1.0.0" + +from .exp import * +from .job import * +from .run import * +from .config import * +from .history import * +from .Date import * +from .filelock import * diff --git a/Python/MAPL/exp.py b/Python/MAPL/exp.py index 70d7119c0a59..11c7c7065b61 100644 --- a/Python/MAPL/exp.py +++ b/Python/MAPL/exp.py @@ -35,7 +35,7 @@ def __del__(self): self.submit() # resubmit itsef def submit(self): - raise NotImplementedError, "Not implemented yet" + raise NotImplementedError("Not implemented yet") # -------------- @@ -74,7 +74,7 @@ def setup(inConfigFiles=None): os.mkdir(tmpdir) os.chdir(tmpdir) if os.system(cmd): - raise IOerror, "red_ma.pl did not complete successfully" + raise IOerror("red_ma.pl did not complete successfully") # Resources as specified by user # ------------------------------ @@ -82,7 +82,7 @@ def setup(inConfigFiles=None): # Setup directory tree # -------------------- - for dir in cf.regex('Dir$').values(): + for dir in list(cf.regex('Dir$').values()): os.mkdir(dir) # Populate Resources diff --git a/Python/MAPL/history.py b/Python/MAPL/history.py index 476da2173652..77baef7226a0 100644 --- a/Python/MAPL/history.py +++ b/Python/MAPL/history.py @@ -2,7 +2,7 @@ A special class for handling history resources. """ -from config import * +from .config import * class History(Config): @@ -35,15 +35,14 @@ def arc(self,outFile): *.temkplate resources. """ dict = self.regex('template$') - Tmpl = [str.replace("'","").replace(",","") for str in dict.values()] - Coll = [ str.split('.')[0].replace(",","") for str in dict.keys() ] + Tmpl = [str.replace("'","").replace(",","") for str in list(dict.values())] + Coll = [ str.split('.')[0].replace(",","") for str in list(dict.keys()) ] Arch = [str.replace("'","").replace(",","") \ - for str in self.regex('archive$').values()] + for str in list(self.regex('archive$').values())] if len(Tmpl) != len(Arch): - raise IOError,\ - "There are %d template resources but only %d archive resources."\ - %(len(Tmpl),len(Arch)) + raise IOError("There are %d template resources but only %d archive resources."\ + %(len(Tmpl),len(Arch))) header = '# PESTO resource for History Collections ' + \ '(automatically generated - do not edit)' diff --git a/Python/MAPL/job.py b/Python/MAPL/job.py index abfd3ce5d571..abf8d79f05de 100644 --- a/Python/MAPL/job.py +++ b/Python/MAPL/job.py @@ -11,8 +11,8 @@ """ -import Abstract -from exp import Exp +from . import Abstract +from .exp import Exp class Job(Exp): @@ -72,13 +72,13 @@ def __del__(self): # ----------------- def getResources(self): - raise NotImplementedError, "Not implemented yet" + raise NotImplementedError("Not implemented yet") def getRecyclables(self): - raise NotImplementedError, "Not implemented yet" + raise NotImplementedError("Not implemented yet") def putRecyclables(self): - raise NotImplementedError, "Not implemented yet" + raise NotImplementedError("Not implemented yet") # ---------------- diff --git a/Python/MAPL/run.py b/Python/MAPL/run.py index 4a6837097230..30056b326714 100644 --- a/Python/MAPL/run.py +++ b/Python/MAPL/run.py @@ -5,7 +5,7 @@ """ -from job import Job +from .job import Job class Run(Job): diff --git a/base/Base/Base_Base_implementation.F90 b/base/Base/Base_Base_implementation.F90 index ae7cf95efa56..151616a5ffd6 100644 --- a/base/Base/Base_Base_implementation.F90 +++ b/base/Base/Base_Base_implementation.F90 @@ -3210,7 +3210,7 @@ function grid_is_ok(grid) result(OK) type(ESMF_Grid), intent(inout) :: grid logical :: OK integer :: I1, I2, J1, J2, j - real(ESMF_KIND_R8), allocatable :: corner_lons(:,:), corner_lats(:,:) + real(ESMF_KIND_R8), pointer :: corner_lons(:,:), corner_lats(:,:) real(ESMF_KIND_R8) :: accurate_lat, accurate_lon real :: tolerance @@ -3218,9 +3218,11 @@ function grid_is_ok(grid) result(OK) call MAPL_GridGetInterior(grid,I1,I2,J1,J2) OK = .true. ! check the edge of face 1 along longitude - allocate(corner_lons(I2-I1+2, J2-J1+2)) - allocate(corner_lats(I2-I1+2, J2-J1+2)) - call MAPL_GridGetCorners(Grid,corner_lons,corner_lats) + call ESMF_GridGetCoord(grid,localDE=0,coordDim=1,staggerloc=ESMF_STAGGERLOC_CORNER, & + farrayPtr=corner_lons, rc=status) + call ESMF_GridGetCoord(grid,localDE=0,coordDim=2,staggerloc=ESMF_STAGGERLOC_CORNER, & + farrayPtr=corner_lats, rc=status) + if ( I1 ==1 .and. J2<=IM_WORLD ) then if (J1 == 1) then accurate_lon = 1.750d0*MAPL_PI_R8 - shift @@ -3233,7 +3235,7 @@ function grid_is_ok(grid) result(OK) endif endif - do j = J1, J2+1 + do j = J1+1, J2 accurate_lat = -alpha + (j-1)*dalpha if ( abs(accurate_lat - corner_lats(1,j-J1+1)) > 5.0*tolerance) then print*, "accurate_lat: ", accurate_lat diff --git a/base/CMakeLists.txt b/base/CMakeLists.txt index a08dacd1250a..268d7291f6f4 100644 --- a/base/CMakeLists.txt +++ b/base/CMakeLists.txt @@ -32,7 +32,7 @@ set (srcs MAPL_IO.F90 MAPL_LatLonGridFactory.F90 MAPL_TransposeRegridder.F90 MAPL_Comms.F90 MAPL_LatLonToLatLonRegridder.F90 MAPL_TripolarGridFactory.F90 - MAPL_LlcGridFactory.F90 + MAPL_LlcGridFactory.F90 MAPL_SwathGridFactory.F90 MAPL_Config.F90 MAPL_LocStreamMod.F90 MAPL_ConservativeRegridder.F90 MAPL_MaxMinMod.F90 MAPL_VerticalInterpMod.F90 MAPL_CubedSphereGridFactory.F90 MAPL_MemUtils.F90 MAPL_VerticalMethods.F90 @@ -55,7 +55,7 @@ set (srcs MAPL_Resource.F90 MAPL_XYGridFactory.F90 MAPL_NetCDF.F90 Plain_netCDF_Time.F90 - MAPL_DateTime_Parsing_ESMF.F90 + MAPL_DateTime_Parsing_ESMF.F90 MAPL_ObsUtil.F90 # Orphaned program: should not be in this library. # tstqsat.F90 ) diff --git a/base/MAPL_AbstractGridFactory.F90 b/base/MAPL_AbstractGridFactory.F90 index 2a422d617991..fabe78cfa803 100644 --- a/base/MAPL_AbstractGridFactory.F90 +++ b/base/MAPL_AbstractGridFactory.F90 @@ -82,6 +82,11 @@ module MAPL_AbstractGridFactoryMod procedure(get_file_format_vars), deferred :: get_file_format_vars procedure(decomps_are_equal), deferred :: decomps_are_equal procedure(physical_params_are_equal), deferred :: physical_params_are_equal + + procedure :: get_xy_subset + procedure :: get_xy_mask + procedure :: destroy + procedure :: get_obs_time end type AbstractGridFactory abstract interface @@ -238,6 +243,7 @@ function generate_file_reference3D(this,fpointer,metadata) result(ref) type(FileMetadata), intent(in), optional :: metaData end function generate_file_reference3D + end interface character(len=*), parameter :: MOD_NAME = 'MAPL_AbstractGridFactory::' @@ -1030,5 +1036,49 @@ function get_grid(this, unusable, rc) result(grid) end if end function get_grid - + + + ! This procedure should only be called for time dependent grids. + ! A default implementation is to fail for other grid types, so we do not + ! have to explicitly add methods to all of the existing subclasses. + subroutine get_xy_subset(this, interval, xy_subset, rc) + class(AbstractGridFactory), intent(in) :: this + type(ESMF_Time), intent(in) :: interval(2) + integer, intent(out) :: xy_subset(2,2) + integer, optional, intent(out) :: rc + integer :: status + + _RETURN(_FAILURE) + end subroutine get_xy_subset + + subroutine get_xy_mask(this, interval, xy_mask, rc) + class(AbstractGridFactory), intent(inout) :: this + type(ESMF_Time), intent(in) :: interval(2) + integer, allocatable, intent(out) :: xy_mask(:,:) + integer, optional, intent(out) :: rc + integer :: status + + _RETURN(_FAILURE) + end subroutine get_xy_mask + + ! Probably don't need to do anything more for subclasses unless they have + ! other objects that don't finalize well. (NetCDF, ESMF, MPI, ...) + subroutine destroy(this, rc) + class(AbstractGridFactory), intent(inout) :: this + integer, optional, intent(out) :: rc + integer :: status + + call ESMF_GridDestroy(this%grid, noGarbage=.true., _RC) + _RETURN(_SUCCESS) + end subroutine destroy + + subroutine get_obs_time(this, grid, obs_time, rc) + class(AbstractGridFactory), intent(inout) :: this + type (ESMF_Grid), intent(in) :: grid + real(ESMF_KIND_R4), intent(out) :: obs_time(:,:) + integer, optional, intent(out) :: rc + + _RETURN(_SUCCESS) + end subroutine get_obs_time + end module MAPL_AbstractGridFactoryMod diff --git a/base/MAPL_AbstractRegridder.F90 b/base/MAPL_AbstractRegridder.F90 index 52aa6364a388..86086af152d3 100644 --- a/base/MAPL_AbstractRegridder.F90 +++ b/base/MAPL_AbstractRegridder.F90 @@ -8,6 +8,7 @@ module MAPL_AbstractRegridderMod use ESMF use MAPL_MemUtilsMod use MAPL_ExceptionHandling + use MAPL_RegridderSpecRouteHandleMap use, intrinsic :: iso_fortran_env, only: REAL32, REAL64 implicit none private @@ -92,6 +93,9 @@ module MAPL_AbstractRegridderMod procedure :: has_undef_value procedure :: get_regrid_method + procedure :: destroy + procedure :: destroy_route_handle + end type AbstractRegridder @@ -1006,4 +1010,21 @@ integer function get_regrid_method(this) result(method) method = this%spec%regrid_method end function get_regrid_method + + subroutine destroy(this, rc) + class(AbstractRegridder), intent(inout) :: this + integer, optional, intent(out) :: rc + integer :: status + + _RETURN(_SUCCESS) + end subroutine destroy + + subroutine destroy_route_handle(this, kind, rc) + class(AbstractRegridder), intent(inout) :: this + type(ESMF_TypeKind_Flag), intent(in) :: kind + integer, optional, intent(out) :: rc + + _RETURN(_SUCCESS) + end subroutine destroy_route_handle + end module MAPL_AbstractRegridderMod diff --git a/base/MAPL_EsmfRegridder.F90 b/base/MAPL_EsmfRegridder.F90 index 382ec9cc2c4f..581545b41c57 100644 --- a/base/MAPL_EsmfRegridder.F90 +++ b/base/MAPL_EsmfRegridder.F90 @@ -53,6 +53,8 @@ module MAPL_EsmfRegridderMod procedure :: do_regrid procedure :: create_route_handle procedure :: select_route_handle + procedure :: destroy + procedure :: destroy_route_handle end type EsmfRegridder @@ -1600,4 +1602,54 @@ function select_route_handle(this, kind, do_transpose, rc) result(route_handle) end function select_route_handle + subroutine destroy(this, rc) + class(EsmfRegridder), intent(inout) :: this + integer, optional, intent(out) :: rc + integer :: status + + call this%destroy_route_handle(ESMF_TYPEKIND_R4, _RC) + + _RETURN(_SUCCESS) + end subroutine destroy + + + subroutine destroy_route_handle(this, kind, rc) + class(EsmfRegridder), intent(inout) :: this + type(ESMF_TypeKind_Flag), intent(in) :: kind + integer, optional, intent(out) :: rc + + type (RegridderSpec) :: spec + type(ESMF_RouteHandle) :: dummy_rh + type(RegridderSpecRouteHandleMap), pointer :: route_handles, transpose_route_handles + type(ESMF_RouteHandle) :: route_handle + type(RegridderSpecRouteHandleMapIterator) :: iter + integer :: status + + if (kind == ESMF_TYPEKIND_R4) then + route_handles => route_handles_r4 + transpose_route_handles => transpose_route_handles_r4 + else if(kind == ESMF_TYPEKIND_R8) then + route_handles => route_handles_r8 + transpose_route_handles => transpose_route_handles_r8 + else + _FAIL('unsupported type kind (must be R4 or R8)') + end if + + spec = this%get_spec() + + _ASSERT(route_handles%count(spec) == 1, 'Did not find this spec in route handle table.') + route_handle = route_handles%at(spec) + call ESMF_RouteHandleDestroy(route_handle, noGarbage=.true.,_RC) + iter = route_handles%find(spec) + call route_handles%erase(iter) + + _ASSERT(transpose_route_handles%count(spec) == 1, 'Did not find this spec in route handle table.') + route_handle = transpose_route_handles%at(spec) + call ESMF_RouteHandleDestroy(route_handle, noGarbage=.true., _RC) + iter = transpose_route_handles%find(spec) + call transpose_route_handles%erase(iter) + + _RETURN(_SUCCESS) + end subroutine destroy_route_handle + end module MAPL_EsmfRegridderMod diff --git a/base/MAPL_GridManager.F90 b/base/MAPL_GridManager.F90 index 360c784ea808..c26cf5db567e 100644 --- a/base/MAPL_GridManager.F90 +++ b/base/MAPL_GridManager.F90 @@ -32,6 +32,9 @@ module MAPL_GridManager_private type (Integer64GridFactoryMap) :: factories contains procedure :: add_prototype + procedure :: destroy_grid + generic :: destroy => destroy_grid + procedure :: delete !!$ procedure :: make_field !!$ procedure :: delete_field @@ -120,6 +123,7 @@ subroutine initialize_prototypes(this, unusable, rc) use MAPL_LlcGridFactoryMod, only: LlcGridFactory use MAPL_ExternalGridFactoryMod, only: ExternalGridFactory use MAPL_XYGridFactoryMod, only: XYGridFactory + use MAPL_SwathGridFactoryMod, only : SwathGridFactory class (GridManager), intent(inout) :: this class (KeywordEnforcer), optional, intent(in) :: unusable @@ -131,7 +135,8 @@ subroutine initialize_prototypes(this, unusable, rc) type (LlcGridFactory) :: llc_factory type (ExternalGridFactory) :: external_factory type (XYGridFactory) :: xy_factory - + type (SwathGridFactory) :: swath_factory + ! This is a local variable to prevent the subroutine from running ! initialiazation twice. Calling functions have their own local variables ! to prevent calling this subroutine twice, but the initialization status @@ -147,6 +152,7 @@ subroutine initialize_prototypes(this, unusable, rc) call this%prototypes%insert('llc', llc_factory) call this%prototypes%insert('External', external_factory) call this%prototypes%insert('XY', xy_factory) + call this%prototypes%insert('Swath', swath_factory) initialized = .true. end if @@ -397,6 +403,27 @@ function make_factory_from_distGrid(this, grid_type, dist_grid, lon_array, lat_a end function make_factory_from_distGrid + subroutine destroy_grid(this, grid, unusable, rc) + use ESMF + class (GridManager), target, intent(inout) :: this + type (ESMF_Grid), intent(inout) :: grid + class (KeywordEnforcer), optional, intent(in) :: unusable + integer, optional, intent(out) :: rc + + integer :: status + integer (kind=ESMF_KIND_I8) :: id + class(AbstractGridFactory), pointer :: factory + type(Integer64GridFactoryMapIterator) :: iter + + call ESMF_AttributeGet(grid, factory_id_attribute, id, _RC) + factory => this%factories%at(id) + call factory%destroy(_RC) + iter = this%factories%find(id) + call this%factories%erase(iter) + + _RETURN(_SUCCESS) + _UNUSED_DUMMY(unusable) + end subroutine destroy_grid ! Clients should use this procedure to release ESMF resources when a grid ! is no longer being used. @@ -413,15 +440,13 @@ subroutine delete(this, grid, unusable, rc) integer :: status character(len=*), parameter :: Iam= MOD_NAME // 'destroy_grid' - _UNUSED_DUMMY(unusable) - if (.not. this%keep_grids) then - call ESMF_GridDestroy(grid, rc=status) + call ESMF_GridDestroy(grid, noGarbage=.true., rc=status) _ASSERT(status==0,'failed to destroy grid') end if _RETURN(_SUCCESS) - + _UNUSED_DUMMY(unusable) end subroutine delete diff --git a/base/MAPL_LatLonToLatLonRegridder.F90 b/base/MAPL_LatLonToLatLonRegridder.F90 index af0a77dffa3f..56e89c29b028 100644 --- a/base/MAPL_LatLonToLatLonRegridder.F90 +++ b/base/MAPL_LatLonToLatLonRegridder.F90 @@ -169,7 +169,7 @@ subroutine compute_binning_weights(Weight,Xin,Xout,HasPoles,rc) dx = Xout(j_out+1)-Xin(j1) ff = ff + dx b(j1) = dx - b = b/ff + b(:) = b(:)/ff end if end associate diff --git a/base/MAPL_ObsUtil.F90 b/base/MAPL_ObsUtil.F90 new file mode 100644 index 000000000000..d4ed2f8de5ab --- /dev/null +++ b/base/MAPL_ObsUtil.F90 @@ -0,0 +1,624 @@ +#include "MAPL_ErrLog.h" +#include "unused_dummy.H" + +module MAPL_ObsUtilMod + use ESMF + use Plain_netCDF_Time + use netCDF + use MAPL_CommsMod, only : MAPL_AM_I_ROOT + use pFIO_FileMetadataMod, only : FileMetadata + use pFIO_NetCDF4_FileFormatterMod, only : NetCDF4_FileFormatter + use, intrinsic :: iso_fortran_env, only: REAL32, REAL64 + implicit none + integer, parameter :: mx_ngeoval = 60 +!! private + + public :: obs_unit + type :: obs_unit + integer :: nobs_epoch + integer :: ngeoval + logical :: export_all_geoval + type(FileMetadata), allocatable :: metadata + type(NetCDF4_FileFormatter), allocatable :: file_handle + character(len=ESMF_MAXSTR) :: name + character(len=ESMF_MAXSTR) :: obsFile_output + character(len=ESMF_MAXSTR) :: input_template + character(len=ESMF_MAXSTR) :: geoval_name(mx_ngeoval) + real(kind=REAL64), allocatable :: lons(:) + real(kind=REAL64), allocatable :: lats(:) + real(kind=REAL64), allocatable :: times_R8(:) + real(kind=REAL32), allocatable :: p2d(:) + real(kind=REAL32), allocatable :: p3d(:,:) + end type obs_unit + + type obs_platform + character (len=ESMF_MAXSTR) :: name='' + character (len=ESMF_MAXSTR) :: nc_index='' + character (len=ESMF_MAXSTR) :: nc_lon='' + character (len=ESMF_MAXSTR) :: nc_lat='' + character (len=ESMF_MAXSTR) :: nc_time='' + character (len=ESMF_MAXSTR) :: file_name_template='' + integer :: ngeoval=0 + integer :: nentry_name=0 + character (len=ESMF_MAXSTR), allocatable :: field_name(:,:) + !character (len=ESMF_MAXSTR), allocatable :: field_name(:) + end type obs_platform + + interface sort_multi_arrays_by_time + module procedure sort_three_arrays_by_time + module procedure sort_four_arrays_by_time + end interface sort_multi_arrays_by_time + +contains + + subroutine get_obsfile_Tbracket_from_epoch(currTime, & + obsfile_start_time, obsfile_end_time, obsfile_interval, & + epoch_frequency, obsfile_Ts_index, obsfile_Te_index, rc) + implicit none + type(ESMF_Time), intent(in) :: currTime + type(ESMF_Time), intent(in) :: obsfile_start_time, obsfile_end_time + type(ESMF_TimeInterval), intent(in) :: obsfile_interval, epoch_frequency + integer, intent(out) :: obsfile_Ts_index + integer, intent(out) :: obsfile_Te_index + integer, optional, intent(out) :: rc + + type(ESMF_Time) :: T1, Tn + type(ESMF_Time) :: cT1 + type(ESMF_Time) :: Ts, Te + type(ESMF_TimeInterval) :: dT1, dT2, dTs, dTe + real(ESMF_KIND_R8) :: dT0_s, dT1_s, dT2_s + real(ESMF_KIND_R8) :: s1, s2 + integer :: n1, n2 + integer :: status + + T1 = obsfile_start_time + Tn = obsfile_end_time + + cT1 = currTime + dT1 = currTime - T1 + dT2 = currTime + epoch_frequency - T1 + + call ESMF_TimeIntervalGet(obsfile_interval, s_r8=dT0_s, rc=status) + call ESMF_TimeIntervalGet(dT1, s_r8=dT1_s, rc=status) + call ESMF_TimeIntervalGet(dT2, s_r8=dT2_s, rc=status) + n1 = floor (dT1_s / dT0_s) + n2 = floor (dT2_s / dT0_s) + s1 = n1 * dT0_s + s2 = n2 * dT0_s + call ESMF_TimeIntervalSet(dTs, s_r8=s1, rc=status) + call ESMF_TimeIntervalSet(dTe, s_r8=s2, rc=status) + Ts = T1 + dTs + Te = T1 + dTe + + obsfile_Ts_index = n1 + if ( dT2_s - n2*dT0_s < 1 ) then + obsfile_Te_index = n2 - 1 + else + obsfile_Te_index = n2 + end if + + _RETURN(ESMF_SUCCESS) + + end subroutine get_obsfile_Tbracket_from_epoch + + + function get_filename_from_template (time, file_template, rc) result(filename) + use Plain_netCDF_Time, only : ESMF_time_to_two_integer + use MAPL_StringTemplate, only : fill_grads_template + type(ESMF_Time), intent(in) :: time + character(len=*), intent(in) :: file_template + character(len=ESMF_MAXSTR) :: filename + integer, optional, intent(out) :: rc + + integer :: itime(2) + integer :: nymd, nhms + integer :: status + + _FAIL ('DO not use get_filename_from_template') + call ESMF_time_to_two_integer(time, itime, _RC) + print*, 'two integer time, itime(:)', itime(1:2) + nymd = itime(1) + nhms = itime(2) + call fill_grads_template ( filename, file_template, & + experiment_id='', nymd=nymd, nhms=nhms, _RC ) + print*, 'ck: obsFile_T=', trim(filename) + _RETURN(ESMF_SUCCESS) + + end function get_filename_from_template + + + subroutine time_real_to_ESMF (times_R8_1d, times_esmf_1d, datetime_units, rc) + use MAPL_NetCDF, only : convert_NetCDF_DateTime_to_ESMF + + real(kind=REAL64), intent(in) :: times_R8_1d(:) + type(ESMF_Time), intent(inout) :: times_esmf_1d(:) + character(len=*), intent(in) :: datetime_units + integer, optional, intent(out) :: rc + + type(ESMF_TimeInterval) :: interval + type(ESMF_Time) :: time0 + type(ESMF_Time) :: time1 + character(len=:), allocatable :: tunit + + integer :: i, len + integer :: int_time + integer :: status + + len = size (times_R8_1d) + do i=1, len + int_time = times_R8_1d(i) + call convert_NetCDF_DateTime_to_ESMF(int_time, datetime_units, interval, & + time0, time=time1, time_unit=tunit, _RC) + times_esmf_1d(i) = time1 + enddo + + _RETURN(_SUCCESS) + end subroutine time_real_to_ESMF + + + subroutine reset_times_to_current_day(current_time, times_1d, rc) + type(ESMF_Time), intent(in) :: current_time + type(ESMF_Time), intent(inout) :: times_1d(:) + integer, optional, intent(out) :: rc + integer :: i,status,h,m,yp,mp,dp,s,ms,us,ns + integer :: year,month,day + + call ESMF_TimeGet(current_time,yy=year,mm=month,dd=day,_RC) + do i=1,size(times_1d) + call ESMF_TimeGet(times_1d(i),yy=yp,mm=mp,dd=dp,h=h,m=m,s=s,ms=ms,us=us,ns=ns,_RC) + call ESMF_TimeSet(times_1d(i),yy=year,mm=month,dd=day,h=h,m=m,s=s,ms=ms,us=us,ns=ns,_RC) + enddo + _RETURN(_SUCCESS) + end subroutine reset_times_to_current_day + + + + + ! --//-------------------------------------//-> + ! files + ! o o o o o o o o o o T: filename + ! <--- off set + ! o o o o o o o o o o T: file content start + ! | | + ! curr curr+Epoch + ! + + subroutine Find_M_files_for_currTime (currTime, & + obsfile_start_time, obsfile_end_time, obsfile_interval, & + epoch_frequency, file_template, M, filenames, & + T_offset_in_file_content, rc) + implicit none + type(ESMF_Time), intent(in) :: currTime + type(ESMF_Time), intent(in) :: obsfile_start_time, obsfile_end_time + type(ESMF_TimeInterval), intent(in) :: obsfile_interval, epoch_frequency + character(len=*), intent(in) :: file_template + integer, intent(out) :: M + character(len=ESMF_MAXSTR), intent(inout) :: filenames(:) + type(ESMF_TimeInterval), intent(in), optional :: T_offset_in_file_content + integer, optional, intent(out) :: rc + + type(ESMF_Time) :: T1, Tn + type(ESMF_Time) :: cT1 + type(ESMF_Time) :: Ts, Te + type(ESMF_TimeInterval) :: dT1, dT2, dTs, dTe + type(ESMF_TimeInterval) :: Toff + real(ESMF_KIND_R8) :: dT0_s, dT1_s, dT2_s + real(ESMF_KIND_R8) :: s1, s2 + character(len=ESMF_MAXSTR) :: test_file + + integer :: obsfile_Ts_index, obsfile_Te_index + integer :: n1, n2 + integer :: i, j + integer :: status + + !__ s1. Arithmetic index list based on s,e,interval + ! + if (present(T_offset_in_file_content)) then + Toff = T_offset_in_file_content + else + call ESMF_TimeIntervalSet(Toff, h=0, m=0, s=60, rc=status) + endif + + ! T1 = obsfile_start_time + Toff + ! Tn = obsfile_end_time + Toff + + T1 = obsfile_start_time + Tn = obsfile_end_time + + cT1 = currTime + dT1 = currTime - T1 + dT2 = currTime + epoch_frequency - T1 + + call ESMF_TimeIntervalGet(obsfile_interval, s_r8=dT0_s, rc=status) + call ESMF_TimeIntervalGet(dT1, s_r8=dT1_s, rc=status) + call ESMF_TimeIntervalGet(dT2, s_r8=dT2_s, rc=status) + + n1 = floor (dT1_s / dT0_s) + n2 = floor (dT2_s / dT0_s) + +! print*, 'ck dT0_s, dT1_s, dT2_s', dT0_s, dT1_s, dT2_s +! print*, '1st n1, n2', n1, n2 + + obsfile_Ts_index = n1 + if ( dT2_s - n2*dT0_s < 1 ) then + obsfile_Te_index = n2 - 1 + else + obsfile_Te_index = n2 + end if + + ! put back + n1 = obsfile_Ts_index + n2 = obsfile_Te_index + +! print*, __LINE__, __FILE__ +! print*, '2nd n1, n2', n1, n2 + + !__ s2. further test file existence + ! + j=0 + do i= n1, n2 + test_file = get_filename_from_template_use_index & + (obsfile_start_time, obsfile_interval, & + i, file_template, rc=rc) + if (test_file /= '') then + j=j+1 + filenames(j) = test_file + end if + end do + M=j + + _ASSERT ( M < size(filenames) , 'code crash, number of files exceeds upper bound') + _ASSERT (M/=0, 'M is zero, no files found for currTime') + + + _RETURN(_SUCCESS) + + end subroutine Find_M_files_for_currTime + + + subroutine read_M_files_4_swath ( filenames, Xdim, Ydim, & + index_name_lon, index_name_lat,& + var_name_lon, var_name_lat, var_name_time, & + lon, lat, time, rc ) + use pFlogger, only: logging, Logger + character(len=ESMF_MAXSTR), intent(in) :: filenames(:) + integer, intent(out) :: Xdim + integer, intent(out) :: Ydim + character(len=ESMF_MAXSTR), intent(in) :: index_name_lon + character(len=ESMF_MAXSTR), intent(in) :: index_name_lat + character(len=ESMF_MAXSTR), optional, intent(in) :: var_name_lon + character(len=ESMF_MAXSTR), optional, intent(in) :: var_name_lat + character(len=ESMF_MAXSTR), optional, intent(in) :: var_name_time + + real, optional, intent(inout) :: lon(:,:) + real, optional, intent(inout) :: lat(:,:) + !! real(ESMF_KIND_R8), optional, intent(inout) :: time_R8(:,:) + real, optional, intent(inout) :: time(:,:) + + integer, optional, intent(out) :: rc + + integer :: M + integer :: i, j, jx, status + integer :: nlon, nlat + integer :: ncid, ncid2 + character(len=ESMF_MAXSTR) :: grp1, grp2 + integer :: varid + logical :: found_group + + character(len=ESMF_MAXSTR) :: filename + integer, allocatable :: nlons(:), nlats(:) + real(ESMF_KIND_R8), allocatable :: time_loc_R8(:,:) + real(ESMF_KIND_R8), allocatable :: lon_loc(:,:) + real(ESMF_KIND_R8), allocatable :: lat_loc(:,:) + class(Logger), pointer :: lgr + + !__ s1. get Xdim Ydim + M = size(filenames) + _ASSERT(M/=0, 'M is zero, no files found') + lgr => logging%get_logger('MAPL.Sampler') + + allocate(nlons(M), nlats(M)) + jx=0 + do i = 1, M + filename = filenames(i) + CALL get_ncfile_dimension(filename, nlon=nlon, nlat=nlat, & + key_lon=index_name_lon, key_lat=index_name_lat, _RC) + nlons(i)=nlon + nlats(i)=nlat + jx=jx+nlat + + call lgr%debug('Input filename: %a', trim(filename)) + call lgr%debug('Input file : nlon, nlat= %i6 %i6', nlon, nlat) + end do + Xdim=nlon + Ydim=jx + + + !__ s2. get fields + jx=0 + do i = 1, M + filename = filenames(i) + nlon = nlons(i) + nlat = nlats(i) + + if (present(var_name_time).AND.present(time)) then + allocate (time_loc_R8(nlon, nlat)) + call get_var_from_name_w_group (var_name_time, time_loc_R8, filename, _RC) + time(1:nlon,jx+1:jx+nlat) = time_loc_R8(1:nlon,1:nlat) + deallocate(time_loc_R8) + end if + + if (present(var_name_lon).AND.present(lon)) then + allocate (lon_loc(nlon, nlat)) + call get_var_from_name_w_group (var_name_lon, lon_loc, filename, _RC) + lon(1:nlon,jx+1:jx+nlat) = lon_loc(1:nlon,1:nlat) + deallocate(lon_loc) + end if + + if (present(var_name_lat).AND.present(lat)) then + allocate (lat_loc(nlon, nlat)) + call get_var_from_name_w_group (var_name_lat, lat_loc, filename, _RC) + lat(1:nlon,jx+1:jx+nlat) = lat_loc(1:nlon,1:nlat) + deallocate(lat_loc) + end if + + jx = jx + nlat + + end do + + _RETURN(_SUCCESS) + end subroutine read_M_files_4_swath + + + ! + !-- caveat: note call this subr. on head node + ! because of (bash ls) command therein + ! + function get_filename_from_template_use_index (obsfile_start_time, obsfile_interval, & + f_index, file_template, rc) result(filename) + use Plain_netCDF_Time, only : ESMF_time_to_two_integer + use MAPL_StringTemplate, only : fill_grads_template + character(len=ESMF_MAXSTR) :: filename + type(ESMF_Time), intent(in) :: obsfile_start_time + type(ESMF_TimeInterval), intent(in) :: obsfile_interval + character(len=*), intent(in) :: file_template + integer, intent(in) :: f_index + integer, optional, intent(out) :: rc + + integer :: itime(2) + integer :: nymd, nhms + integer :: status + real(ESMF_KIND_R8) :: dT0_s + real(ESMF_KIND_R8) :: s + type(ESMF_TimeInterval) :: dT + type(ESMF_Time) :: time + integer :: i, j, u + logical :: EX + + character(len=ESMF_MAXSTR) :: file_template_left + character(len=ESMF_MAXSTR) :: file_template_right + character(len=ESMF_MAXSTR) :: filename_left + character(len=ESMF_MAXSTR) :: filename_full + character(len=ESMF_MAXSTR) :: filename2 + character(len=ESMF_MAXSTR) :: cmd + + call ESMF_TimeIntervalGet(obsfile_interval, s_r8=dT0_s, rc=status) + s = dT0_s * f_index + call ESMF_TimeIntervalSet(dT, s_r8=s, rc=status) + time = obsfile_start_time + dT + + call ESMF_time_to_two_integer(time, itime, _RC) + nymd = itime(1) + nhms = itime(2) + + ! parse time info + ! + call fill_grads_template ( filename, file_template, & + experiment_id='', nymd=nymd, nhms=nhms, _RC ) + inquire(file= trim(filename), EXIST = EX) + if(.not.EX) filename='' + + _RETURN(_SUCCESS) + + end function get_filename_from_template_use_index + + + + subroutine get_var_from_name_w_group (var_name, var2d, filename, rc) + character(len=ESMF_MAXSTR), intent(in) :: var_name, filename + real(ESMF_KIND_R8), intent(inout) :: var2d(:,:) + integer, optional, intent(out) :: rc + + integer :: i, j + character(len=ESMF_MAXSTR) :: grp1, grp2 + character(len=ESMF_MAXSTR) :: short_name + integer :: ncid, ncid2, varid + logical :: found_group + integer :: status + + + i=index(var_name, '/') + if (i>0) then + found_group = .true. + grp1 = var_name(1:i-1) + j=index(var_name(i+1:), '/') + if (j>0) then + grp2=var_name(i+1:i+j-1) + short_name=var_name(i+j+1:) + else + grp2='' + short_name=var_name(i+1:) + endif + i=i+j + else + found_group = .false. + grp1 = '' + grp2='' + short_name=var_name + endif + + call check_nc_status(nf90_open(filename, NF90_NOWRITE, ncid2), _RC) + if ( found_group ) then + call check_nc_status(nf90_inq_ncid(ncid2, grp1, ncid), _RC) + if (j>0) then + call check_nc_status(nf90_inq_ncid(ncid, grp2, ncid2), _RC) + ncid=ncid2 + endif + else + print*, 'no grp name' + ncid=ncid2 + endif + call check_nc_status(nf90_inq_varid(ncid, short_name, varid), _RC) + call check_nc_status(nf90_get_var(ncid, varid, var2d), _RC) +!! call check_nc_status(nf90_close(ncid), _RC) + + _RETURN(_SUCCESS) + + end subroutine get_var_from_name_w_group + + + subroutine sort_three_arrays_by_time(U,V,T,rc) + use MAPL_SortMod + real(ESMF_KIND_R8), intent(inout) :: U(:), V(:), T(:) + integer, optional, intent(out) :: rc + + integer :: i, len + integer, allocatable :: IA(:) + integer(ESMF_KIND_I8), allocatable :: IX(:) + real(ESMF_KIND_R8), allocatable :: X(:) + + _ASSERT (size(U)==size(V), 'U,V different dimension') + _ASSERT (size(U)==size(T), 'U,T different dimension') + len = size (T) + + allocate (IA(len), IX(len), X(len)) + do i=1, len + IX(i)=T(i) + IA(i)=i + enddo + call MAPL_Sort(IX,IA) + + X = U + do i=1, len + U(i) = X(IA(i)) + enddo + X = V + do i=1, len + V(i) = X(IA(i)) + enddo + X = T + do i=1, len + T(i) = X(IA(i)) + enddo + _RETURN(_SUCCESS) + end subroutine sort_three_arrays_by_time + + + subroutine sort_four_arrays_by_time(U,V,T,ID,rc) + use MAPL_SortMod + real(ESMF_KIND_R8) :: U(:), V(:), T(:) + integer :: ID(:) + integer, optional, intent(out) :: rc + + integer :: i, len + integer, allocatable :: IA(:) + integer(ESMF_KIND_I8), allocatable :: IX(:) + real(ESMF_KIND_R8), allocatable :: X(:) + integer, allocatable :: NX(:) + + _ASSERT(size(U)==size(V), 'U,V different dimension') + _ASSERT(size(U)==size(T), 'U,T different dimension') + len = size (T) + + allocate (IA(len), IX(len), X(len), NX(len)) + do i=1, len + IX(i)=T(i) + IA(i)=i + enddo + call MAPL_Sort(IX,IA) + + X = U + do i=1, len + U(i) = X(IA(i)) + enddo + X = V + do i=1, len + V(i) = X(IA(i)) + enddo + X = T + do i=1, len + T(i) = X(IA(i)) + enddo + NX = ID + do i=1, len + ID(i) = NX(IA(i)) + enddo + _RETURN(_SUCCESS) + end subroutine sort_four_arrays_by_time + + + + function copy_platform_nckeys(a, rc) + type(obs_platform) :: copy_platform_nckeys + type(obs_platform), intent(in) :: a + integer, optional, intent(out) :: rc + + copy_platform_nckeys%nc_index = a%nc_index + copy_platform_nckeys%nc_lon = a%nc_lon + copy_platform_nckeys%nc_lat = a%nc_lat + copy_platform_nckeys%nc_time = a%nc_time + copy_platform_nckeys%nentry_name = a%nentry_name + _RETURN(_SUCCESS) + + end function copy_platform_nckeys + + + function union_platform(a, b, rc) + type(obs_platform) :: union_platform + type(obs_platform), intent(in) :: a + type(obs_platform), intent(in) :: b + integer, optional, intent(out) :: rc + + character (len=ESMF_MAXSTR), allocatable :: field_name_loc(:,:) + integer :: nfield, nentry_name + integer, allocatable :: tag(:) + integer :: i, j, k + integer :: status + + union_platform = copy_platform_nckeys(a, _RC) + nfield = a%ngeoval + b%ngeoval + allocate (tag(b%ngeoval)) + + tag(:)=1 ! true + k=nfield + do j=1, b%ngeoval + do i=1, a%ngeoval + if ( trim(b%field_name(1,j)) == trim(a%field_name(1,i)) ) then + tag(j)=0 + endif + enddo + if (tag(j)==0) k=k-1 + enddo + union_platform%ngeoval=k + nfield=k + nentry_name=union_platform%nentry_name + if ( allocated (union_platform%field_name) ) deallocate(union_platform%field_name) + allocate(union_platform%field_name(nentry_name, nfield)) + do i=1, a%ngeoval + union_platform%field_name(:,i) = a%field_name(:,i) + enddo + if (nfield>a%ngeoval) then + k = a%ngeoval + do j=1, b%ngeoval + if (tag(j)==1) then + k = k + 1 + union_platform%field_name(:,k) = b%field_name(:,j) + end if + enddo + end if + _RETURN(_SUCCESS) + + end function union_platform + + +end module MAPL_ObsUtilMod diff --git a/base/MAPL_SwathGridFactory.F90 b/base/MAPL_SwathGridFactory.F90 new file mode 100644 index 000000000000..591c9eb562cc --- /dev/null +++ b/base/MAPL_SwathGridFactory.F90 @@ -0,0 +1,1445 @@ +#include "MAPL_Exceptions.h" +#include "MAPL_ErrLog.h" +#include "unused_dummy.H" + +module MAPL_SwathGridFactoryMod + use MAPL_AbstractGridFactoryMod + use MAPL_MinMaxMod + use MAPL_KeywordEnforcerMod + use MAPL_ExceptionHandling + use MAPL_ShmemMod + use mapl_ErrorHandlingMod + use MAPL_Constants + use MAPL_Base, only : MAPL_GridGetInterior + use ESMF + use pFIO + use MAPL_CommsMod + !!use netcdf + !!use Plain_netCDF_Time + use MAPL_ObsUtilMod + use pflogger, only : Logger, logging + use, intrinsic :: iso_fortran_env, only: REAL32 + use, intrinsic :: iso_fortran_env, only: REAL64 + implicit none + integer, parameter :: gridLabel_max = 20 + integer, parameter :: mx_file = 300 + private + + public :: SwathGridFactory + + type, extends(AbstractGridFactory) :: SwathGridFactory + private + character(len=:), allocatable :: grid_name + character(len=:), allocatable :: grid_file_name + character(len=ESMF_MAXSTR) :: filenames(mx_file) + integer :: M_file + + integer :: cell_across_swath + integer :: cell_along_swath + integer :: im_world = MAPL_UNDEFINED_INTEGER + integer :: jm_world = MAPL_UNDEFINED_INTEGER + integer :: lm = MAPL_UNDEFINED_INTEGER + logical :: force_decomposition = .false. + + integer :: epoch ! unit: second + integer(ESMF_KIND_I8) :: epoch_index(4) ! is,ie,js,je + real(ESMF_KIND_R8), allocatable:: t_alongtrack(:) + ! note: this var is not deallocated in swathfactory, use caution + character(len=ESMF_MAXSTR) :: tunit + character(len=ESMF_MAXSTR) :: index_name_lon + character(len=ESMF_MAXSTR) :: index_name_lat + character(len=ESMF_MAXSTR) :: var_name_lon + character(len=ESMF_MAXSTR) :: var_name_lat + character(len=ESMF_MAXSTR) :: var_name_time + character(len=ESMF_MAXSTR) :: input_template + logical :: found_group + + type(ESMF_Time) :: obsfile_start_time ! user specify + type(ESMF_Time) :: obsfile_end_time + type(ESMF_TimeInterval) :: obsfile_interval + type(ESMF_TimeInterval) :: EPOCH_FREQUENCY + integer :: obsfile_Ts_index ! for epoch + integer :: obsfile_Te_index + logical :: is_valid + + ! Domain decomposition: + integer :: nx = MAPL_UNDEFINED_INTEGER + integer :: ny = MAPL_UNDEFINED_INTEGER + integer, allocatable :: ims(:) + integer, allocatable :: jms(:) + ! Used for halo + type (ESMF_DELayout) :: layout + logical :: initialized_from_metadata = .false. + contains + procedure :: make_new_grid + procedure :: create_basic_grid + procedure :: add_horz_coordinates_from_file + procedure :: init_halo + procedure :: halo + + procedure :: initialize_from_file_metadata + procedure :: initialize_from_config_with_prefix + procedure :: initialize_from_esmf_distGrid + + procedure :: equals + procedure :: check_and_fill_consistency + procedure :: generate_grid_name + procedure :: to_string + + procedure :: append_metadata + procedure :: get_grid_vars + procedure :: get_file_format_vars + procedure :: append_variable_metadata + procedure :: check_decomposition + procedure :: generate_newnxy + procedure :: generate_file_bounds + procedure :: generate_file_corner_bounds + procedure :: generate_file_reference2D + procedure :: generate_file_reference3D + procedure :: decomps_are_equal + procedure :: physical_params_are_equal + + procedure :: get_xy_subset + procedure :: destroy + procedure :: get_obs_time + end type SwathGridFactory + + character(len=*), parameter :: MOD_NAME = 'MAPL_SwathGridFactory::' + + interface SwathGridFactory + module procedure SwathGridFactory_from_parameters + end interface SwathGridFactory + + interface set_with_default + module procedure set_with_default_integer + module procedure set_with_default_real + module procedure set_with_default_real64 + module procedure set_with_default_character + module procedure set_with_default_bounds + end interface set_with_default + +contains + + function SwathGridFactory_from_parameters(unusable, grid_name, & + & im_world, jm_world, lm, nx, ny, ims, jms, rc) result(factory) + type (SwathGridFactory) :: factory + class (KeywordEnforcer), optional, intent(in) :: unusable + character(len=*), optional, intent(in) :: grid_name + + ! grid details: + integer, optional, intent(in) :: im_world + integer, optional, intent(in) :: jm_world + integer, optional, intent(in) :: lm + + ! decomposition: + integer, optional, intent(in) :: nx + integer, optional, intent(in) :: ny + integer, optional, intent(in) :: ims(:) + integer, optional, intent(in) :: jms(:) + + integer, optional, intent(out) :: rc + + integer :: status + + _UNUSED_DUMMY(unusable) + + call set_with_default(factory%grid_name, grid_name, MAPL_GRID_NAME_DEFAULT) + call set_with_default(factory%nx, nx, MAPL_UNDEFINED_INTEGER) + call set_with_default(factory%ny, ny, MAPL_UNDEFINED_INTEGER) + call set_with_default(factory%im_world, im_world, MAPL_UNDEFINED_INTEGER) + call set_with_default(factory%jm_world, jm_world, MAPL_UNDEFINED_INTEGER) + call set_with_default(factory%lm, lm, MAPL_UNDEFINED_INTEGER) + + ! default is unallocated + if (present(ims)) factory%ims = ims + if (present(jms)) factory%jms = jms + + call factory%check_and_fill_consistency(_RC) + + _RETURN(_SUCCESS) + end function SwathGridFactory_from_parameters + + + function make_new_grid(this, unusable, rc) result(grid) + type (ESMF_Grid) :: grid + class (SwathGridFactory), intent(in) :: this + class (KeywordEnforcer), optional, intent(in) :: unusable + integer, optional, intent(out) :: rc + integer :: status + + _UNUSED_DUMMY(unusable) + grid = this%create_basic_grid(_RC) + call this%add_horz_coordinates_from_file(grid,_RC) + _RETURN(_SUCCESS) + end function make_new_grid + + + function create_basic_grid(this, unusable, rc) result(grid) + type (ESMF_Grid) :: grid + class (SwathGridFactory), intent(in) :: this + class (KeywordEnforcer), optional, intent(in) :: unusable + integer, optional, intent(out) :: rc + integer :: status + + _UNUSED_DUMMY(unusable) + + grid = ESMF_GridCreateNoPeriDim( & + & name = this%grid_name, & + & countsPerDEDim1=this%ims, & + & countsPerDEDim2=this%jms, & + & indexFlag=ESMF_INDEX_DELOCAL, & + & coordDep1=[1,2], & + & coordDep2=[1,2], & + & coordSys=ESMF_COORDSYS_SPH_RAD, & + & _RC) + + ! Allocate coords at default stagger location + call ESMF_GridAddCoord(grid, _RC) + + if (this%lm /= MAPL_UNDEFINED_INTEGER) then + call ESMF_AttributeSet(grid, name='GRID_LM', value=this%lm, _RC) + end if + call ESMF_AttributeSet(grid, 'GridType', 'LatLon', _RC) + call ESMF_AttributeSet(grid, 'Global', .false., _RC) + + _RETURN(_SUCCESS) + end function create_basic_grid + + + subroutine add_horz_coordinates_from_file(this, grid, unusable, rc) + use MAPL_BaseMod, only: MAPL_grid_interior + implicit none + class (SwathGridFactory), intent(in) :: this + type (ESMF_Grid), intent(inout) :: grid + class (KeywordEnforcer), optional, intent(in) :: unusable + integer, optional, intent(out) :: rc + integer :: status + + real(kind=ESMF_KIND_R8), pointer :: fptr(:,:) + real, pointer :: centers(:,:) + real, allocatable :: centers_full(:,:) + + integer :: i, j, k + integer :: Xdim, Ydim + integer :: Xdim_full, Ydim_full + integer :: nx, ny + + integer :: IM, JM + integer :: IM_WORLD, JM_WORLD + integer :: COUNTS(3), DIMS(3) + integer :: i_1, i_n, j_1, j_n ! regional array bounds + type(Logger), pointer :: lgr + + _UNUSED_DUMMY(unusable) + + + Xdim=this%im_world + Ydim=this%jm_world + Xdim_full=this%cell_across_swath + Ydim_full=this%cell_along_swath + + call MAPL_grid_interior(grid, i_1, i_n, j_1, j_n) + call MAPL_AllocateShared(centers,[Xdim,Ydim],transroot=.true.,_RC) + call MAPL_SyncSharedMemory(_RC) + + +! if (mapl_am_I_root()) then +! write(6,'(2x,a,10i8)') & +! 'ck: Xdim, Ydim, Xdim_full, Ydim_full', Xdim, Ydim, Xdim_full, Ydim_full +! write(6,'(2x,a,10i8)') & +! 'ck: i_1, i_n, j_1, j_n', i_1, i_n, j_1, j_n +! end if + + + ! read longitudes + if (MAPL_AmNodeRoot .or. (.not. MAPL_ShmInitialized)) then + allocate( centers_full(Xdim_full, Ydim_full)) + call read_M_files_4_swath (this%filenames(1:this%M_file), nx, ny, & + this%index_name_lon, this%index_name_lat, & + var_name_lon=this%var_name_lon, lon=centers_full, _RC) + k=0 + do j=this%epoch_index(3), this%epoch_index(4) + k=k+1 + centers(1:Xdim, k) = centers_full(1:Xdim, j) + enddo + centers=centers*MAPL_DEGREES_TO_RADIANS_R8 + deallocate (centers_full) + end if + call MAPL_SyncSharedMemory(_RC) + call ESMF_GridGetCoord(grid, coordDim=1, localDE=0, & + staggerloc=ESMF_STAGGERLOC_CENTER, farrayPtr=fptr, _RC) + fptr=real(centers(i_1:i_n,j_1:j_n), kind=ESMF_KIND_R8) + + + ! read latitudes + if (MAPL_AmNodeRoot .or. (.not. MAPL_ShmInitialized)) then + allocate( centers_full(Xdim_full, Ydim_full)) + call read_M_files_4_swath (this%filenames(1:this%M_file), nx, ny, & + this%index_name_lon, this%index_name_lat, & + var_name_lat=this%var_name_lat, lat=centers_full, _RC) + k=0 + do j=this%epoch_index(3), this%epoch_index(4) + k=k+1 + centers(1:Xdim, k) = centers_full(1:Xdim, j) + enddo + centers=centers*MAPL_DEGREES_TO_RADIANS_R8 + deallocate (centers_full) + end if + call MAPL_SyncSharedMemory(_RC) + call ESMF_GridGetCoord(grid, coordDim=2, localDE=0, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + farrayPtr=fptr, rc=status) + fptr=real(centers(i_1:i_n,j_1:j_n), kind=ESMF_KIND_R8) + + if(MAPL_ShmInitialized) then + call MAPL_DeAllocNodeArray(centers,_RC) + else + deallocate(centers) + end if + + _RETURN(_SUCCESS) + end subroutine add_horz_coordinates_from_file + + + subroutine initialize_from_file_metadata(this, file_metadata, unusable, force_file_coordinates, rc) + use MAPL_KeywordEnforcerMod + use MAPL_BaseMod, only: MAPL_DecomposeDim + + class (SwathGridFactory), intent(inout) :: this + type (FileMetadata), target, intent(in) :: file_metadata + class (KeywordEnforcer), optional, intent(in) :: unusable + logical, optional, intent(in) :: force_file_coordinates + integer, optional, intent(out) :: rc + + integer :: status + + class (CoordinateVariable), pointer :: v + class (*), pointer :: ptr(:) + + character(:), allocatable :: lon_name + character(:), allocatable :: lat_name + character(:), allocatable :: lev_name + integer :: i + logical :: hasLon, hasLat, hasLongitude, hasLatitude, hasLev,hasLevel,regLat,regLon + real(kind=REAL64) :: del12,delij + + integer :: i_min, i_max + real(kind=REAL64) :: d_lat, d_lat_temp, extrap_lat + logical :: is_valid, use_file_coords, compute_lons, compute_lats + + _UNUSED_DUMMY(unusable) + + if (present(force_file_coordinates)) then + use_file_coords = force_file_Coordinates + else + use_file_coords = .false. + end if + + ! Cannot assume that lats and lons are evenly spaced + + associate (im => this%im_world, jm => this%jm_world, lm => this%lm) + lon_name = 'lon' + hasLon = file_metadata%has_dimension(lon_name) + if (hasLon) then + im = file_metadata%get_dimension(lon_name, _RC) + else + lon_name = 'longitude' + hasLongitude = file_metadata%has_dimension(lon_name) + if (hasLongitude) then + im = file_metadata%get_dimension(lon_name, _RC) + else + _FAIL('no longitude coordinate') + end if + end if + lat_name = 'lat' + hasLat = file_metadata%has_dimension(lat_name) + if (hasLat) then + jm = file_metadata%get_dimension(lat_name, _RC) + else + lat_name = 'latitude' + hasLatitude = file_metadata%has_dimension(lat_name) + if (hasLatitude) then + jm = file_metadata%get_dimension(lat_name, _RC) + else + _FAIL('no latitude coordinate') + end if + end if + hasLev=.false. + hasLevel=.false. + lev_name = 'lev' + hasLev = file_metadata%has_dimension(lev_name) + if (hasLev) then + lm = file_metadata%get_dimension(lev_name,_RC) + else + lev_name = 'levels' + hasLevel = file_metadata%has_dimension(lev_name) + if (hasLevel) then + lm = file_metadata%get_dimension(lev_name,_RC) + end if + end if + end associate + + call this%make_arbitrary_decomposition(this%nx, this%ny, _RC) + + ! Determine IMS and JMS with constraint for ESMF that each DE has at least an extent + ! of 2. Required for ESMF_FieldRegrid(). + allocate(this%ims(0:this%nx-1)) + allocate(this%jms(0:this%ny-1)) + call MAPL_DecomposeDim(this%im_world, this%ims, this%nx, min_DE_extent=2) + call MAPL_DecomposeDim(this%jm_world, this%jms, this%ny, min_DE_extent=2) + + call this%check_and_fill_consistency(_RC) + + _RETURN(_SUCCESS) + + end subroutine initialize_from_file_metadata + + + subroutine initialize_from_config_with_prefix(this, config, prefix, unusable, rc) + use MPI + implicit none + class (SwathGridFactory), intent(inout) :: this + type (ESMF_Config), intent(inout) :: config + character(len=*), intent(in) :: prefix + class (KeywordEnforcer), optional, intent(in) :: unusable + integer, optional, intent(out) :: rc + integer :: status + + type(ESMF_VM) :: VM + integer:: mpic + integer:: irank, ierror + integer :: nlon, nlat, tdim + integer :: Xdim, Ydim, ntime + integer :: nx, ny + character(len=ESMF_MAXSTR) :: key_lon, key_lat, key_time + character(len=ESMF_MAXSTR) :: tunit, grp1, grp2 + character(len=ESMF_MAXSTR) :: filename, STR1, tmp + character(len=ESMF_MAXSTR) :: symd, shms + + + ! real(ESMF_KIND_R8), allocatable :: scanTime(:,:) + real, allocatable :: scanTime(:,:) + integer :: yy, mm, dd, h, m, s, sec, second + integer :: i, j, L + integer :: ncid, ncid2, varid + integer :: fid_s, fid_e + integer :: M_file + + type(ESMF_Time) :: currTime + integer (ESMF_KIND_I8) :: j0, j1, jt, jt1, jt2 + real(ESMF_KIND_R8) :: jx0, jx1 + real(ESMF_KIND_R8) :: x0, x1 + integer :: khi, klo, k, nstart, max_iter + type(Logger), pointer :: lgr + logical :: ispresent + + type(ESMF_TimeInterval) :: Toff + + _UNUSED_DUMMY(unusable) + lgr => logging%get_logger('HISTORY.sampler') + + call ESMF_VmGetCurrent(VM, _RC) + + ! input : config + ! output: this%epoch_index, nx, ny + ! + ! Read in specs, crop epoch_index based on scanTime + ! + + + !__ s1. read in file spec. + ! + call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'GRIDNAME:', default=MAPL_GRID_NAME_DEFAULT) + this%grid_name = trim(tmp) + call ESMF_ConfigGetAttribute(config, this%nx, label=prefix//'NX:', default=MAPL_UNDEFINED_INTEGER) + call ESMF_ConfigGetAttribute(config, this%ny, label=prefix//'NY:', default=MAPL_UNDEFINED_INTEGER) + call ESMF_ConfigGetAttribute(config, this%lm, label=prefix//'LM:', default=MAPL_UNDEFINED_INTEGER) + call ESMF_ConfigGetAttribute(config, this%input_template, label=prefix//'GRID_FILE:', default='unknown.txt', _RC) + call ESMF_ConfigGetAttribute(config, this%epoch, label=prefix//'Epoch:', default=300, _RC) + call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'Epoch_init:', default='2006', _RC) + + + call ESMF_ConfigGetAttribute(config, value=STR1, default="", & + label= prefix// 'obs_file_begin:', _RC) + + if (trim(STR1)=='') then + _FAIL('obs_file_begin missing, code crash') + else + call ESMF_TimeSet(this%obsfile_start_time, timestring=STR1, _RC) + end if + + call ESMF_ConfigGetAttribute(config, value=STR1, default="", & + label=prefix // 'obs_file_end:', _RC) + + if (trim(STR1)=='') then + _FAIL('obs_file_end missing, code crash') + else + call ESMF_TimeSet(this%obsfile_end_time, timestring=STR1, _RC) + end if + + call ESMF_ConfigGetAttribute(config, value=STR1, default="", & + label= prefix// 'obs_file_interval:', _RC) + _ASSERT(STR1/='', 'fatal error: obs_file_interval not provided in RC file') + + +! if (mapl_am_I_root()) then +! write(6,'(//2x, a)') 'SWATH initialize_from_config_with_prefix' +! print*, 'obs_file_begin: str1=', trim(STR1) +! write(6,105) 'obs_file_begin provided: ', trim(STR1) +! print*, 'obs_file_end: str1=', trim(STR1) +! write(6,105) 'obs_file_end provided:', trim(STR1) +! write(6,105) 'obs_file_interval:', trim(STR1) +! write(6,106) 'Epoch (hhmmss) :', this%epoch +! end if + + + i= index( trim(STR1), ' ' ) + if (i>0) then + symd=STR1(1:i-1) + shms=STR1(i+1:) + else + symd='' + shms=trim(STR1) + endif + call convert_twostring_2_esmfinterval (symd, shms, this%obsfile_interval, _RC) + + second = hms_2_s(this%Epoch) + call ESMF_TimeIntervalSet(this%epoch_frequency, s=second, _RC) + + if ( index(tmp, 'T') /= 0 .OR. index(tmp, '-') /= 0 ) then + call ESMF_TimeSet(currTime, timeString=tmp, _RC) + else + read(tmp,'(i4,5i2)') yy,mm,dd,h,m,s + call ESMF_Timeset(currTime, yy=yy, mm=mm, dd=dd, h=h, m=m, s=s, _RC) + endif + + call lgr%debug(' %a %a', 'input_template =', trim(this%input_template)) + !!write(6,'(2x,a,/,4i8,/,5(2x,a))') 'nx,ny,lm,epoch -- filename,tmp', & + !! this%nx,this%ny,this%lm,this%epoch,& + !! trim(filename),trim(tmp) + !!print*, 'ck: Epoch_init:', trim(tmp) + + + call ESMF_ConfigGetAttribute(config, value=this%index_name_lon, default="", & + label=prefix // 'index_name_lon:', _RC) + call ESMF_ConfigGetAttribute(config, value=this%index_name_lat, default="", & + label=prefix // 'index_name_lat:', _RC) + call ESMF_ConfigGetAttribute(config, this%var_name_lon, & + label=prefix // 'var_name_lon:', default="", _RC) + call ESMF_ConfigGetAttribute(config, this%var_name_lat, & + label=prefix // 'var_name_lat:', default="", _RC) + call ESMF_ConfigGetAttribute(config, this%var_name_time, default="", & + label=prefix//'var_name_time:', _RC) + call ESMF_ConfigGetAttribute(config, this%tunit, default="", & + label=prefix//'tunit:', _RC) + + + + !__ s2. find obsFile even if missing on disk and get array: this%t_alongtrack(:) + ! + call ESMF_VMGet(vm, mpiCommunicator=mpic, _RC) + call MPI_COMM_RANK(mpic, irank, ierror) + + if (irank==0) & + write(6,'(10(2x,a20,2x,a40,/))') & + 'index_name_lon:', trim(this%index_name_lon), & + 'index_name_lat:', trim(this%index_name_lat), & + 'var_name_lon:', trim(this%var_name_lon), & + 'var_name_lat:', trim(this%var_name_lat), & + 'var_name_time:', trim(this%var_name_time), & + 'tunit:', trim(this%tunit) + + if (irank==0) then + call ESMF_TimeIntervalSet(Toff, h=0, m=0, s=0, _RC) + call Find_M_files_for_currTime (currTime, & + this%obsfile_start_time, this%obsfile_end_time, this%obsfile_interval, & + this%epoch_frequency, this%input_template, M_file, this%filenames, & + T_offset_in_file_content = Toff, _RC) + this%M_file = M_file + write(6,'(10(2x,a20,2x,i40))') & + 'M_file:', M_file + do i=1, M_file + write(6,'(10(2x,a20,2x,a))') & + 'filenames(i):', trim(this%filenames(i)) + end do + + call read_M_files_4_swath (this%filenames(1:M_file), nx, ny, & + this%index_name_lon, this%index_name_lat, _RC) + nlon=nx + nlat=ny + allocate(scanTime(nlon, nlat)) + allocate(this%t_alongtrack(nlat)) + + call read_M_files_4_swath (this%filenames(1:M_file), nx, ny, & + this%index_name_lon, this%index_name_lat, & + var_name_time=this%var_name_time, time=scanTime, _RC) + + + do j=1, nlat + this%t_alongtrack(j)= scanTime(1,j) + enddo + nstart = 1 + ! + ! redefine nstart to skip un-defined time value + ! If the t_alongtrack contains undefined values, use this code + ! + x0 = this%t_alongtrack(1) + x1 = 1.d16 + if (x0 > x1) then + ! + ! bisect backward finding the first index arr[n] < x1 + klo=1 + khi=nlat + max_iter = int( log( real(nlat) ) / log(2.d0) ) + 2 + do i=1, max_iter + k = (klo+khi)/2 + if ( this%t_alongtrack(k) < x1 ) then + khi=k + else + nstart = khi + exit + endif + enddo + call lgr%debug('%a %i4', 'nstart', nstart) + call lgr%debug('%a %i4', 'this%t_alongtrack(nstart)', this%t_alongtrack(nstart)) + endif + + this%cell_across_swath = nlon + this%cell_along_swath = nlat + deallocate(scanTime) +!! write(6,*) 'this%t_alongtrack(j)=', this%t_alongtrack(::100) + + + ! P2. + ! determine im_world from Epoch + ! ----------------------------- + ! t_axis = t_alongtrack = t_a + ! convert currTime to j0 + ! use Epoch to find j1 + ! search j0, j1 in t_a + + call time_esmf_2_nc_int (currTime, this%tunit, j0, _RC) + sec = hms_2_s (this%Epoch) + j1= j0 + sec + jx0= j0 + jx1= j1 + call lgr%debug ('%a %i16 %i16', 'j0, j1 ', j0, j1) + + + this%epoch_index(1)= 1 + this%epoch_index(2)= this%cell_across_swath + call bisect( this%t_alongtrack, jx0, jt1, n_LB=int(nstart, ESMF_KIND_I8), n_UB=int(this%cell_along_swath, ESMF_KIND_I8), rc=rc) + call bisect( this%t_alongtrack, jx1, jt2, n_LB=int(nstart, ESMF_KIND_I8), n_UB=int(this%cell_along_swath, ESMF_KIND_I8), rc=rc) + + + if (jt1==jt2) then + _FAIL('Epoch Time is too small, empty swath grid is generated, increase Epoch') + endif + jt1 = jt1 + 1 ! (x1,x2] design + this%epoch_index(3)= jt1 + this%epoch_index(4)= jt2 + Xdim = this%cell_across_swath + Ydim = this%epoch_index(4) - this%epoch_index(3) + 1 + + call lgr%debug ('%a %i4 %i4', 'bisect for j0: rc, jt', rc, jt1) + call lgr%debug ('%a %i4 %i4', 'bisect for j1: rc, jt', rc, jt2) + call lgr%debug ('%a %i4 %i4', 'Xdim, Ydim', Xdim, Ydim) + call lgr%debug ('%a %i4 %i4 %i4 %i4', 'this%epoch_index(4)', & + this%epoch_index(1), this%epoch_index(2), & + this%epoch_index(3), this%epoch_index(4)) + + this%im_world = Xdim + this%jm_world = Ydim + end if + + call MPI_bcast(this%M_file, 1, MPI_INTEGER, 0, mpic, ierror) + do i=1, this%M_file + call MPI_bcast(this%filenames(i), ESMF_MAXSTR, MPI_CHARACTER, 0, mpic, ierror) + end do + call MPI_bcast(this%epoch_index, 4, MPI_INTEGER8, 0, mpic, ierror) + call MPI_bcast(this%im_world, 1, MPI_INTEGER, 0, mpic, ierror) + call MPI_bcast(this%jm_world, 1, MPI_INTEGER, 0, mpic, ierror) + call MPI_bcast(this%cell_across_swath, 1, MPI_INTEGER, 0, mpic, ierror) + call MPI_bcast(this%cell_along_swath, 1, MPI_INTEGER, 0, mpic, ierror) + ! donot need to bcast this%along_track (root only) + + call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'IMS_FILE:', rc=status) + if ( status == _SUCCESS ) then + call get_ims_from_file(this%ims, trim(tmp),this%nx, _RC) + else + call get_multi_integer(this%ims, 'IMS:', _RC) + endif + call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'JMS_FILE:', rc=status) + if ( status == _SUCCESS ) then + call get_ims_from_file(this%jms, trim(tmp),this%ny, _RC) + else + call get_multi_integer(this%jms, 'JMS:', _RC) + endif + ! ims is set at here + call this%check_and_fill_consistency(_RC) + + + _RETURN(_SUCCESS) + +105 format (1x,a,2x,a) +106 format (1x,a,2x,10i8) + + contains + + subroutine get_multi_integer(values, label, rc) + integer, allocatable, intent(out) :: values(:) + character(len=*) :: label + integer, optional, intent(out) :: rc + + integer :: i + integer :: n + integer :: tmp + integer :: status + logical :: isPresent + + call ESMF_ConfigFindLabel(config, label=prefix//label, isPresent=isPresent, _RC) + + if (.not. isPresent) then + _RETURN(_SUCCESS) + end if + + ! First pass: count values + n = 0 + do + call ESMF_ConfigGetAttribute(config, tmp, rc=status) + if (status /= _SUCCESS) then + exit + else + n = n + 1 + end if + end do + + + ! Second pass: allocate and fill + allocate(values(n), stat=status) ! no point in checking status + _VERIFY(status) + call ESMF_ConfigFindLabel(config, label=prefix//label,_RC) + do i = 1, n + call ESMF_ConfigGetAttribute(config, values(i), _RC) + write(6,*) 'values(i)=', values(i) + end do + + _RETURN(_SUCCESS) + + end subroutine get_multi_integer + + subroutine get_ims_from_file(values, file_name, n, rc) + integer, allocatable, intent(out) :: values(:) + character(len=*), intent(in) :: file_name + integer, intent(in) :: n + integer, optional, intent(out) :: rc + + logical :: FileExists + integer :: i, total, unit + integer :: status + + inquire(FILE = trim(file_name), EXIST=FileExists) + allocate(values(n), stat=status) ! no point in checking status + _VERIFY(status) + + _ASSERT(FileExists, "File <"//trim(file_name)//"> not found") + if (MAPL_AM_I_Root(VM)) then + open(newunit=UNIT, file=trim(file_name), form="formatted", iostat=status ) + _VERIFY(STATUS) + read(UNIT,*) total + _ASSERT(total == n, trim(file_name) // " n is different from total") + do i = 1,total + read(UNIT,*) values(i) + enddo + close(UNIT) + endif + + call MAPL_CommsBcast(VM, values, n=N, ROOT=MAPL_Root, _RC) + _RETURN(_SUCCESS) + + end subroutine get_ims_from_file + + subroutine get_range(range, label, rc) + type(RealMinMax), intent(out) :: range + character(len=*) :: label + integer, optional, intent(out) :: rc + + integer :: i + integer :: n + integer :: status + logical :: isPresent + + call ESMF_ConfigFindLabel(config, label=prefix//label,isPresent=isPresent,_RC) + if (.not. isPresent) then + _RETURN(_SUCCESS) + end if + + ! Must be 2 values: min and max + call ESMF_ConfigGetAttribute(config, range%min, _RC) + call ESMF_ConfigGetAttribute(config, range%max, _RC) + + _RETURN(_SUCCESS) + + end subroutine get_range + + + end subroutine initialize_from_config_with_prefix + + + + function to_string(this) result(string) + character(len=:), allocatable :: string + class (SwathGridFactory), intent(in) :: this + + _UNUSED_DUMMY(this) + string = 'SwathGridFactory' + + end function to_string + + + subroutine check_and_fill_consistency(this, unusable, rc) + use MAPL_BaseMod, only: MAPL_DecomposeDim + class (SwathGridFactory), intent(inout) :: this + class (KeywordEnforcer), optional, intent(in) :: unusable + integer, optional, intent(out) :: rc + + integer :: status + logical :: verify_decomp + + _UNUSED_DUMMY(unusable) + + if (.not. allocated(this%grid_name)) then + this%grid_name = MAPL_GRID_NAME_DEFAULT + end if + + ! Check decomposition/bounds + ! WY notes: should not have this assert + !_ASSERT(allocated(this%ims) .eqv. allocated(this%jms), 'inconsistent options') + call verify(this%nx, this%im_world, this%ims, rc=status) + call verify(this%ny, this%jm_world, this%jms, rc=status) + + if (.not.this%force_decomposition) then + verify_decomp = this%check_decomposition(_RC) + if ( (.not.verify_decomp) ) then + call this%generate_newnxy(_RC) + end if + end if + + _RETURN(_SUCCESS) + + contains + + subroutine verify(n, m_world, ms, rc) + integer, intent(inout) :: n + integer, intent(inout) :: m_world + integer, allocatable, intent(inout) :: ms(:) + integer, optional, intent(out) :: rc + + integer :: status + + if (allocated(ms)) then + _ASSERT(size(ms) > 0, 'degenerate topology') + + if (n == MAPL_UNDEFINED_INTEGER) then + n = size(ms) + else + _ASSERT(n == size(ms), 'inconsistent topology') + end if + + if (m_world == MAPL_UNDEFINED_INTEGER) then + m_world = sum(ms) + else + _ASSERT(m_world == sum(ms), 'inconsistent decomponsition') + end if + + else + + _ASSERT(n /= MAPL_UNDEFINED_INTEGER, 'uninitialized topology') + _ASSERT(m_world /= MAPL_UNDEFINED_INTEGER,'uninitialized dimension') + allocate(ms(n), stat=status) + _VERIFY(status) + !call MAPL_DecomposeDim(m_world, ms, n, min_DE_extent=2) + call MAPL_DecomposeDim(m_world, ms, n) + + end if + + _RETURN(_SUCCESS) + + end subroutine verify + + end subroutine check_and_fill_consistency + + + elemental subroutine set_with_default_integer(to, from, default) + integer, intent(out) :: to + integer, optional, intent(in) :: from + integer, intent(in) :: default + + if (present(from)) then + to = from + else + to = default + end if + + end subroutine set_with_default_integer + + elemental subroutine set_with_default_real64(to, from, default) + real(REAL64), intent(out) :: to + real(REAL64), optional, intent(in) :: from + real(REAL64), intent(in) :: default + + if (present(from)) then + to = from + else + to = default + end if + + end subroutine set_with_default_real64 + + elemental subroutine set_with_default_real(to, from, default) + real, intent(out) :: to + real, optional, intent(in) :: from + real, intent(in) :: default + + if (present(from)) then + to = from + else + to = default + end if + + end subroutine set_with_default_real + + subroutine set_with_default_character(to, from, default) + character(len=:), allocatable, intent(out) :: to + character(len=*), optional, intent(in) :: from + character(len=*), intent(in) :: default + + if (present(from)) then + to = from + else + to = default + end if + + end subroutine set_with_default_character + + + elemental subroutine set_with_default_bounds(to, from, default) + type (RealMinMax), intent(out) :: to + type (RealMinMax), optional, intent(in) :: from + type (RealMinMax), intent(in) :: default + + if (present(from)) then + to = from + else + to = default + end if + + end subroutine set_with_default_bounds + + + ! MAPL uses values in lon_array and lat_array only to determine the + ! general positioning. Actual coordinates are then recomputed. + ! This helps to avoid roundoff differences from slightly different + ! input files. + subroutine initialize_from_esmf_distGrid(this, dist_grid, lon_array, lat_array, unusable, rc) + use MAPL_ConfigMod + use MAPL_Constants, only: PI => MAPL_PI_R8 + class (SwathGridFactory), intent(inout) :: this + type (ESMF_DistGrid), intent(in) :: dist_grid + type (ESMF_LocalArray), intent(in) :: lon_array + type (ESMF_LocalArray), intent(in) :: lat_array + class (KeywordEnforcer), optional, intent(in) :: unusable + integer, optional, intent(out) :: rc + + integer :: dim_count, tile_count + integer, allocatable :: max_index(:,:) + integer :: status + character(len=2) :: pole ,dateline + + type (ESMF_Config) :: config + type (ESMF_VM) :: vm + integer :: nPet + real(kind=REAL32), pointer :: lon(:) + real(kind=REAL32), pointer :: lat(:) + integer :: nx_guess,nx,ny + integer :: i + + real, parameter :: tiny = 1.e-4 + + _FAIL ('stop: not implemented: subroutine initialize_from_esmf_distGrid') + + _UNUSED_DUMMY(unusable) + + call ESMF_DistGridGet(dist_grid, dimCount=dim_count, tileCount=tile_count) + allocate(max_index(dim_count, tile_count)) + call ESMF_DistGridGet(dist_grid, maxindexPTile=max_index) + + config = MAPL_ConfigCreate(_RC) + call MAPL_ConfigSetAttribute(config, max_index(1,1), 'IM_WORLD:', _RC) + call MAPL_ConfigSetAttribute(config, max_index(2,1), 'JM_WORLD:', _RC) + call MAPL_ConfigSetAttribute(config, max_index(3,1), 'LM:', _RC) + + lon => null() + lat => null() + call ESMF_LocalArrayGet(lon_array, farrayPtr=lon, _RC) + call ESMF_LocalArrayGet(lat_array, farrayPtr=lat, _RC) + + call ESMF_VMGetCurrent(vm, _RC) + call ESMF_VMGet(vm, PETcount=nPet, _RC) + + nx_guess = nint(sqrt(real(nPet))) + do nx = nx_guess,1,-1 + ny=nPet/nx + if (nx*ny==nPet) then + call MAPL_ConfigSetAttribute(config, nx, 'NX:') + call MAPL_ConfigSetAttribute(config, ny, 'NY:') + exit + end if + enddo + + call this%initialize(config, _RC) + + end subroutine initialize_from_esmf_distGrid + + + function decomps_are_equal(this,a) result(equal) + class (SwathGridFactory), intent(in) :: this + class (AbstractGridFactory), intent(in) :: a + logical :: equal + + select type (a) + class default + equal = .false. + return + class is (SwathGridFactory) + equal = .true. + + + equal = size(a%ims)==size(this%ims) .and. size(a%jms)==size(this%jms) + if (.not. equal) return + + ! same decomposition + equal = all(a%ims == this%ims) .and. all(a%jms == this%jms) + if (.not. equal) return + + end select + + end function decomps_are_equal + + + function physical_params_are_equal(this, a) result(equal) + class (SwathGridFactory), intent(in) :: this + class (AbstractGridFactory), intent(in) :: a + logical :: equal + + select type (a) + class default + equal = .false. + return + class is (SwathGridFactory) + equal = .true. + + equal = (a%im_world == this%im_world) .and. (a%jm_world == this%jm_world) + if (.not. equal) return + end select + + end function physical_params_are_equal + + logical function equals(a, b) + class (SwathGridFactory), intent(in) :: a + class (AbstractGridFactory), intent(in) :: b + + select type (b) + class default + equals = .false. + return + class is (SwathGridFactory) + equals = .true. + + equals = (a%lm == b%lm) + if (.not. equals) return + + equals = a%decomps_are_equal(b) + if (.not. equals) return + + equals = a%physical_params_are_equal(b) + if (.not. equals) return + + end select + + end function equals + + + function generate_grid_name(this) result(name) + character(len=:), allocatable :: name + class (SwathGridFactory), intent(in) :: this +! from tclune: This needs thought. I suspect we want something that indicates this is a swath grid. + character(len=4) :: im_string, jm_string + name = im_string // 'x' // jm_string + end function generate_grid_name + + + function check_decomposition(this,unusable,rc) result(can_decomp) + class (SwathGridFactory), target, intent(inout) :: this + class (KeywordEnforcer), optional, intent(in) :: unusable + integer, optional, intent(out) :: rc + logical :: can_decomp + integer :: n + _UNUSED_DUMMY(unusable) + + can_decomp = .true. + if (this%im_world==1 .and. this%jm_world==1) then + _RETURN(_SUCCESS) + end if + n = this%im_world/this%nx + if (n < 2) can_decomp = .false. + n = this%jm_world/this%ny + if (n < 2) can_decomp = .false. + _RETURN(_SUCCESS) + end function check_decomposition + + + subroutine generate_newnxy(this,unusable,rc) + use MAPL_BaseMod, only: MAPL_DecomposeDim + class (SwathGridFactory), target, intent(inout) :: this + class (KeywordEnforcer), optional, intent(in) :: unusable + integer, optional, intent(out) :: rc + integer :: n + + _UNUSED_DUMMY(unusable) + + n = this%im_world/this%nx + if (n < 2) then + this%nx = generate_new_decomp(this%im_world,this%nx) + deallocate(this%ims) + allocate(this%ims(0:this%nx-1)) + call MAPL_DecomposeDim(this%im_world, this%ims, this%nx) + end if + n = this%jm_world/this%ny + if (n < 2) then + this%ny = generate_new_decomp(this%jm_world,this%ny) + deallocate(this%jms) + allocate(this%jms(0:this%ny-1)) + call MAPL_DecomposeDim(this%jm_world, this%jms, this%ny) + end if + + _RETURN(_SUCCESS) + + end subroutine generate_newnxy + + function generate_new_decomp(im,nd) result(n) + integer, intent(in) :: im, nd + integer :: n + logical :: canNotDecomp + + canNotDecomp = .true. + n = nd + do while(canNotDecomp) + if ( (im/n) < 2) then + n = n/2 + else + canNotDecomp = .false. + end if + enddo + end function generate_new_decomp + + subroutine init_halo(this, unusable, rc) + class (SwathGridFactory), target, intent(inout) :: this + class (KeywordEnforcer), optional, intent(in) :: unusable + integer, optional, intent(out) :: rc + _FAIL('Stop: subroutine init_halo is not needed for SwathGridFactory') + end subroutine init_halo + + subroutine halo(this, array, unusable, halo_width, rc) + use MAPL_CommsMod + class (SwathGridFactory), intent(inout) :: this + real(kind=REAL32), intent(inout) :: array(:,:) + class (KeywordEnforcer), optional, intent(in) :: unusable + integer, optional, intent(in) :: halo_width + integer, optional, intent(out) :: rc + _FAIL( 'Stop: subroutine halo is not needed for SwathGridFactory') + end subroutine halo + + + subroutine append_metadata(this, metadata) + use MAPL_Constants + class (SwathGridFactory), intent(inout) :: this + type (FileMetadata), intent(inout) :: metadata + + type (Variable) :: v + real(kind=REAL64), allocatable :: temp_coords(:) + + character(len=ESMF_MAXSTR) :: key_lon + character(len=ESMF_MAXSTR) :: key_lat + + ! Horizontal grid dimensions + call metadata%add_dimension('lon', this%im_world) + call metadata%add_dimension('lat', this%jm_world) + + ! Coordinate variables + v = Variable(type=PFIO_REAL64, dimensions='lon,lat') + call v%add_attribute('long_name', 'longitude') + call v%add_attribute('units', 'degrees_east') + call metadata%add_variable('lons', v) + + v = Variable(type=PFIO_REAL64, dimensions='lon,lat') + call v%add_attribute('long_name', 'latitude') + call v%add_attribute('units', 'degrees_north') + call metadata%add_variable('lats', v) + + end subroutine append_metadata + + + function get_grid_vars(this) result(vars) + class (SwathGridFactory), intent(inout) :: this + + character(len=:), allocatable :: vars + character(len=ESMF_MAXSTR) :: key_lon + character(len=ESMF_MAXSTR) :: key_lat + _UNUSED_DUMMY(this) + + !!key_lon=trim(this%var_name_lon) + !!key_lat=trim(this%var_name_lat) + vars = 'lon,lat' + + end function get_grid_vars + + + function get_file_format_vars(this) result(vars) + class (SwathGridFactory), intent(inout) :: this + + character(len=:), allocatable :: vars + _UNUSED_DUMMY(this) + + vars = 'lon,lat' + end function get_file_format_vars + + + subroutine append_variable_metadata(this,var) + class (SwathGridFactory), intent(inout) :: this + type(Variable), intent(inout) :: var + _UNUSED_DUMMY(this) + _UNUSED_DUMMY(var) + end subroutine append_variable_metadata + + + subroutine generate_file_bounds(this,grid,local_start,global_start,global_count,metadata,rc) + use MAPL_BaseMod + class(SwathGridFactory), intent(inout) :: this + type(ESMF_Grid), intent(inout) :: grid + integer, allocatable, intent(out) :: local_start(:) + integer, allocatable, intent(out) :: global_start(:) + integer, allocatable, intent(out) :: global_count(:) + type(FileMetaData), intent(in), optional :: metaData + integer, optional, intent(out) :: rc + + integer :: status + integer :: global_dim(3), i1,j1,in,jn + + _UNUSED_DUMMY(this) + + call MAPL_GridGet(grid,globalCellCountPerDim=global_dim,_RC) + call MAPL_GridGetInterior(grid,i1,in,j1,jn) + allocate(local_start,source=[i1,j1]) + allocate(global_start,source=[1,1]) + allocate(global_count,source=[global_dim(1),global_dim(2)]) + + _RETURN(_SUCCESS) + + end subroutine generate_file_bounds + + + subroutine generate_file_corner_bounds(this,grid,local_start,global_start,global_count,rc) + use esmf + class (SwathGridFactory), intent(inout) :: this + type(ESMF_Grid), intent(inout) :: grid + integer, allocatable, intent(out) :: local_start(:) + integer, allocatable, intent(out) :: global_start(:) + integer, allocatable, intent(out) :: global_count(:) + integer, optional, intent(out) :: rc + + _UNUSED_DUMMY(this) + _UNUSED_DUMMY(grid) + _UNUSED_DUMMY(local_start) + _UNUSED_DUMMY(global_start) + _UNUSED_DUMMY(global_count) + + _FAIL('unimplemented') + _RETURN(_SUCCESS) + end subroutine generate_file_corner_bounds + + function generate_file_reference2D(this,fpointer) result(ref) + use pFIO + type(ArrayReference) :: ref + class(SwathGridFactory), intent(inout) :: this + real, pointer, intent(in) :: fpointer(:,:) + _UNUSED_DUMMY(this) + ref = ArrayReference(fpointer) + end function generate_file_reference2D + + function generate_file_reference3D(this,fpointer,metaData) result(ref) + use pFIO + type(ArrayReference) :: ref + class(SwathGridFactory), intent(inout) :: this + real, pointer, intent(in) :: fpointer(:,:,:) + type(FileMetaData), intent(in), optional :: metaData + _UNUSED_DUMMY(this) + ref = ArrayReference(fpointer) + end function generate_file_reference3D + + + subroutine get_xy_subset(this, interval, xy_subset, rc) + use MPI + class(SwathGridFactory), intent(in) :: this + type(ESMF_Time), intent(in) :: interval(2) + integer, intent(out) :: xy_subset(2,2) + integer, optional, intent(out) :: rc + + type(ESMF_VM) :: VM + integer:: mpic + integer:: irank, ierror + + integer :: status + type(ESMF_Time) :: T1, T2 + integer(ESMF_KIND_I8) :: i1, i2 + real(ESMF_KIND_R8) :: iT1, iT2 + integer(ESMF_KIND_I8) :: index1, index2 + integer :: jlo, jhi, je + + + call ESMF_VmGetCurrent(VM, _RC) + call ESMF_VMGet(vm, mpiCommunicator=mpic, _RC) + call MPI_COMM_RANK(mpic, irank, ierror) + + if (irank==0) then + ! xtrack + xy_subset(1:2,1)=this%epoch_index(1:2) + + ! atrack + T1= interval(1) + T2= interval(2) + + ! this%t_alongtrack + ! + call time_esmf_2_nc_int (T1, this%tunit, i1, _RC) + call time_esmf_2_nc_int (T2, this%tunit, i2, _RC) + iT1 = i1 ! int to real*8 + iT2 = i2 + if (this%epoch_index(3) > 2) then + jlo = this%epoch_index(3) - 2 + else + jlo = this%epoch_index(3) + end if + jhi = this%epoch_index(4) + 1 + ! + ! -- it is possible some obs files are missing + ! + call bisect( this%t_alongtrack, iT1, index1, n_LB=int(jlo, ESMF_KIND_I8), n_UB=int(jhi, ESMF_KIND_I8), rc=rc) + call bisect( this%t_alongtrack, iT2, index2, n_LB=int(jlo, ESMF_KIND_I8), n_UB=int(jhi, ESMF_KIND_I8), rc=rc) + + !! complex version + !! ! (x1, x2] design in bisect + !! if (index1==jlo-1) then + !! je = index1 + 1 + !! else + !! je = index1 + !! end if + !! xy_subset(1, 2) = je + !! if (index2==jlo-1) then + !! je = index2 + 1 + !! else + !! je = index2 + !! end if + !! xy_subset(2, 2) = je + + ! simple version + xy_subset(1, 2)=index1+1 ! atrack + xy_subset(2, 2)=index2 + + ! + !- relative + ! + xy_subset(1,2)= xy_subset(1,2) - this%epoch_index(3) + 1 + xy_subset(2,2)= xy_subset(2,2) - this%epoch_index(3) + 1 + end if + + call MPI_bcast(xy_subset, 4, MPI_INTEGER, 0, mpic, ierror) + + _RETURN(_SUCCESS) + end subroutine get_xy_subset + + + subroutine destroy(this, rc) + class(SwathGridFactory), intent(inout) :: this + integer, optional, intent(out) :: rc + integer :: i + return + end subroutine destroy + + + ! here grid == external_grid + ! because this%grid is protected in AbstractGridFactory + subroutine get_obs_time(this, grid, obs_time, rc) + use MAPL_BaseMod, only: MAPL_grid_interior + class(SwathGridFactory), intent(inout) :: this + type (ESMF_Grid), intent(in) :: grid + real(ESMF_KIND_R4), intent(out) :: obs_time(:,:) + integer, optional, intent(out) :: rc + integer :: status + + integer :: i_1, i_n, j_1, j_n ! regional array bounds + + !! shared mem + real(kind=ESMF_KIND_R8), pointer :: fptr(:,:) + real, pointer :: centers(:,:) + real, allocatable :: centers_full(:,:) + + integer :: i, j, k + integer :: Xdim, Ydim + integer :: Xdim_full, Ydim_full + integer :: nx, ny + integer :: IM_WORLD, JM_WORLD + + + !- shared mem case in MPI + ! + Xdim=this%im_world + Ydim=this%jm_world + Xdim_full=this%cell_across_swath + Ydim_full=this%cell_along_swath + + call MAPL_grid_interior(grid, i_1, i_n, j_1, j_n) + call MAPL_AllocateShared(centers,[Xdim,Ydim],transroot=.true.,_RC) + call MAPL_SyncSharedMemory(_RC) + + + ! read Time and set + if (MAPL_AmNodeRoot .or. (.not. MAPL_ShmInitialized)) then + allocate( centers_full(Xdim_full, Ydim_full)) + call read_M_files_4_swath (this%filenames(1:this%M_file), nx, ny, & + this%index_name_lon, this%index_name_lat, & + var_name_time=this%var_name_time, time=centers_full, _RC) + !!call get_v2d_netcdf(this%grid_file_name, time_name, centers_full, Xdim_full, Ydim_full) + k=0 + do j=this%epoch_index(3), this%epoch_index(4) + k=k+1 + centers(1:Xdim, k) = centers_full(1:Xdim, j) + enddo + deallocate (centers_full) + end if + call MAPL_SyncSharedMemory(_RC) + + !(Xdim, Ydim) + obs_time = centers(i_1:i_n,j_1:j_n) + + if(MAPL_ShmInitialized) then + call MAPL_DeAllocNodeArray(centers,_RC) + else + deallocate(centers) + end if + + _RETURN(_SUCCESS) + end subroutine get_obs_time + + +end module MAPL_SwathGridFactoryMod diff --git a/base/Plain_netCDF_Time.F90 b/base/Plain_netCDF_Time.F90 index 85ff1507b407..be20b3d76bb1 100644 --- a/base/Plain_netCDF_Time.F90 +++ b/base/Plain_netCDF_Time.F90 @@ -87,7 +87,6 @@ subroutine get_ncfile_dimension(filename, nlon, nlat, tdim, key_lon, key_lat, ke lat_name=trim(key_lat) call check_nc_status(nf90_inq_dimid(ncid, trim(lat_name), dimid), _RC) call check_nc_status(nf90_inquire_dimension(ncid, dimid, len=nlat), _RC) - call check_nc_status(nf90_close(ncid), _RC) endif if(present(key_time)) then @@ -142,7 +141,7 @@ subroutine get_attribute_from_group(filename, group_name, var_name, attr_name, a endif attr = str(1:i+5)//trim(str2) deallocate(str) - call check_nc_status(nf90_close(ncid), _RC) + call check_nc_status(nf90_close(ncid2), _RC) _RETURN(_SUCCESS) @@ -225,8 +224,11 @@ subroutine get_v1d_netcdf_R8(filename, name, array, Xdim, group_name, rc) end if call check_nc_status(nf90_inq_varid(ncid, name, varid), _RC) call check_nc_status(nf90_get_var(ncid, varid, array), _RC) - call check_nc_status(nf90_close(ncid), _RC) - + if(present(group_name)) then + call check_nc_status(nf90_close(ncid2), _RC) + else + call check_nc_status(nf90_close(ncid), _RC) + end if _RETURN(_SUCCESS) end subroutine get_v1d_netcdf_R8 @@ -446,10 +448,10 @@ subroutine bisect_find_LB_R8_I8(xa, x, n, n_LB, n_UB, rc) integer :: i, nmax LB=1; UB=size(xa,1) - if(present(n_LB)) LB=n_LB - if(present(n_UB)) UB=n_UB + if(present(n_LB)) LB=max(LB, n_LB) + if(present(n_UB)) UB=min(UB, n_UB) klo=LB; khi=UB; dk=1 - + if ( xa(LB ) > xa(UB) ) then klo= UB khi= LB @@ -512,3 +514,203 @@ subroutine convert_twostring_2_esmfinterval(symd, shms, interval, rc) end subroutine convert_twostring_2_esmfinterval end module Plain_NetCDF_Time + + +module MAPL_scan_pattern_in_file + +! procedure :: matchbgn +! procedure :: matches +! procedure :: scan_begin +! procedure :: scan_contain +! generic :: scan_count_matchbgn +! generic :: go_last_pattern + +contains + subroutine scan_begin (iunps, substring, rew) + implicit none + ! unit of input + integer, intent(in) :: iunps + ! Label to be matched + character (len=*), intent(in) :: substring + logical, intent(in) :: rew + ! String read from file + character (len=100) :: line + ! Flag if .true. rewind the file + !logical, external :: matchbgn +! logical :: matchbgn + integer :: ios + ! + ios = 0 + if (rew) rewind (iunps) + do while (ios==0) + read (iunps, '(a100)', iostat = ios) line + if (matchbgn (trim(line), trim(substring)) ) return + enddo + return + end subroutine scan_begin + + + subroutine scan_contain (iunps, stop_string, rew) + !--------------------------------------------------------------------- + ! + implicit none + integer, intent(in) :: iunps + character (len=*), intent(in) :: stop_string + logical, intent(in) :: rew ! if rewind + character (len=100) :: line +!! logical :: matches ! function name + integer :: ios + ! + ios = 0 + if (rew) rewind (iunps) + do while (ios==0) + read (iunps, '(a100)', iostat = ios) line + if (matches (trim(line), trim(stop_string)) ) return + enddo + return + end subroutine scan_contain + + + + subroutine scan_count_match_bgn (iunps, string, count, rew) + !--------------------------------------------------------------------- + ! + implicit none + integer, intent(in) :: iunps + character (len=*), intent(in) :: string + integer, intent(out) :: count + logical, intent(in) :: rew ! if rewind + character (len=100) :: line +!! logical :: matches ! function name + integer :: ios + ! + ios = 0 + count = 0 + if (rew) rewind (iunps) + do while (ios==0) + read (iunps, '(a100)', iostat = ios) line + if (matchbgn (line, string) ) then + count = count + 1 + endif + enddo + return + end subroutine scan_count_match_bgn + + + subroutine go_last_patn (iunps, substring, outline, rew) + !--------------------------------------------------------------------- + ! + implicit none + integer, intent(in) :: iunps + logical, intent(in) :: rew + character (len=*), intent(in) :: substring + character (len=150), intent(out) :: outline ! fixed + character (len=150) :: line + integer :: ios, nr, mx + + if (rew) rewind (iunps) + ios=0 + nr=0 + do while (ios==0) + read (iunps, '(a150)', iostat = ios, err = 300) line + if (index(line, substring).ge.1 ) then + nr=nr+1 + ! write (6,*) 'nr', nr + endif + enddo + + rewind (iunps) + ios=0 + mx=0 + do while (ios==0) + read (iunps, '(a150)', iostat = ios, err = 300) line + if (index(line, substring).ge.1 ) then + mx=mx+1 + if (mx.eq.nr) then + outline=line + return + endif + endif + enddo +300 continue + end subroutine go_last_patn + + + function matchbgn ( string, substring ) + ! only begin with + ! string: main-str + ! substring: sub-str + IMPLICIT NONE + CHARACTER (LEN=*), INTENT(IN) :: string, substring + LOGICAL :: matchbgn + if (index(string, substring).eq.1) then + matchbgn = .TRUE. + else + matchbgn = .FALSE. + endif + return + end function matchbgn + + + !----------------------------------------------------------------------- + function matches( string, substring ) + !----------------------------------------------------------------------- + ! + ! ... .TRUE. if string is contained in substring, .FALSE. otherwise + ! + IMPLICIT NONE + ! + CHARACTER (LEN=*), INTENT(IN) :: string, substring + LOGICAL :: matches + INTEGER :: l + + l=index (string, substring) + if (l.ge.1) then + matches = .TRUE. + else + matches = .FALSE. + endif + RETURN + end function matches + + + subroutine split_string_by_space (string_in, length_mx, & + mxseg, nseg, str_piece, jstatus) + integer, intent (in) :: length_mx + character (len=length_mx), intent (in) :: string_in + integer, intent (in) :: mxseg + integer, intent (out):: nseg + character (len=length_mx), intent (out):: str_piece(mxseg) + integer, intent (out):: jstatus + character (len=length_mx) :: string + character (len=1) :: mark + integer :: ios + integer :: wc + ! + ! "xxxx yy zz uu vv" + ! + ! split by space '' + mark=' ' + wc=0 + ios=0 + string = trim( adjustl(string_in) ) + do while (ios==0) + i = index (string, mark) + !!print*, 'index=', i + if (i > 1) then + wc = wc + 1 + str_piece(wc)=trim(adjustl(string(1:i))) + !!write(6,*) 'str_piece(wc)=', trim(str_piece(wc)) + string = trim(adjustl(string(i:))) + else + ios=1 + end if + if (LEN_TRIM(adjustl(string)) == 0) ios=1 + end do + nseg=wc + return + end subroutine split_string_by_space + + + +end module MAPL_scan_pattern_in_file diff --git a/base/StringTemplate.F90 b/base/StringTemplate.F90 index 1b13af15edfa..c3efbdeecece 100644 --- a/base/StringTemplate.F90 +++ b/base/StringTemplate.F90 @@ -5,13 +5,14 @@ module MAPL_StringTemplate use MAPL_ExceptionHandling use MAPL_KeywordEnforcerMod + implicit none private public fill_grads_template public StrTemplate -character(len=2), parameter :: valid_tokens(14) = ["y4","y2","m1","m2","mc","Mc","MC","d1","d2","h1","h2","h3","n2","S2"] +character(len=2), parameter :: valid_tokens(15) = ["y4","y2","m1","m2","mc","Mc","MC","d1","d2","h1","h2","h3","n2","S2","D3"] character(len=3),parameter :: mon_lc(12) = [& 'jan','feb','mar','apr','may','jun', & 'jul','aug','sep','oct','nov','dec'] @@ -165,6 +166,10 @@ function evaluate_token(token,year,month,day,hour,minute,second,preserve) result logical, intent(in) :: preserve character(len=4) :: buffer character(len=1) :: c1,c2 + type(ESMF_Time) :: time + integer(ESMF_KIND_I4) :: doy + integer :: status, rc + type(ESMF_Calendar) :: gregorianCalendar c1=token(1:1) c2=token(2:2) select case(c1) @@ -208,6 +213,22 @@ function evaluate_token(token,year,month,day,hour,minute,second,preserve) result else buffer="%"//token end if + case("D") ! dayOfYear + if (.not.skip_token(day,preserve)) then + if (c2 == "3") then + gregorianCalendar = ESMF_CalendarCreate(ESMF_CALKIND_GREGORIAN, & + name='Gregorian_obs' , rc=rc) + call ESMF_TimeSet(time, yy=year, mm=month, dd=day, & + calendar=gregorianCalendar, _RC) + call ESMF_TimeGet(time, dayOfYear=doy, _RC) + call ESMF_CalendarDestroy(gregorianCalendar) + write(buffer,'(i3.3)')doy + else + _FAIL('Day of Year must be %D3') + end if + else + buffer="%"//token + end if case("h") if (.not.skip_token(hour,preserve)) then if (c2 == "3") then diff --git a/components.yaml b/components.yaml index 970c7762769f..e0bbd84af99c 100644 --- a/components.yaml +++ b/components.yaml @@ -5,7 +5,7 @@ MAPL: ESMA_env: local: ./ESMA_env remote: ../ESMA_env.git - tag: v4.20.5 + tag: v4.24.0 develop: main ESMA_cmake: diff --git a/docs/tutorial/grid_comps/CMakeLists.txt b/docs/tutorial/grid_comps/CMakeLists.txt index 2a006d2e4feb..9cdb243357f7 100644 --- a/docs/tutorial/grid_comps/CMakeLists.txt +++ b/docs/tutorial/grid_comps/CMakeLists.txt @@ -4,3 +4,4 @@ add_subdirectory (leaf_comp_a) add_subdirectory (leaf_comp_b) add_subdirectory (parent_with_one_child) add_subdirectory (parent_with_two_children) +add_subdirectory (automatic_code_generator_example) diff --git a/docs/tutorial/grid_comps/automatic_code_generator_example/ACG_GridComp.F90 b/docs/tutorial/grid_comps/automatic_code_generator_example/ACG_GridComp.F90 new file mode 100644 index 000000000000..c3eed7ab9585 --- /dev/null +++ b/docs/tutorial/grid_comps/automatic_code_generator_example/ACG_GridComp.F90 @@ -0,0 +1,113 @@ +#include "MAPL_Generic.h" +#include "MAPL_Exceptions.h" +!------------------------------------------------------------------------------ +!> +!### MODULE: `ACG_GridComp` +! +! This module is created to show how to automatically regenerate code segments +! for the registration and access of ESMF states member variables. +! It is not meant to be executed in an application but only to be compiled. +! +module ACG_GridComp + + use ESMF + use MAPL + + implicit none + private + + public SetServices + +!------------------------------------------------------------------------------ + contains +!------------------------------------------------------------------------------ +!> +! `SetServices` uses MAPL_GenericSetServices, which sets +! the Initialize and Finalize services to generic versions. +! It also allocates our instance of a generic state and puts it in the +! gridded component (GC). Here we only set the run method and +! declare the data services. +! + subroutine SetServices(GC,rc) + + type(ESMF_GridComp), intent(inout) :: GC !! gridded component + integer, optional :: rc !! return code + + integer :: status + + call MAPL_GridCompSetEntryPoint ( gc, ESMF_METHOD_INITIALIZE, initialize, _RC) + call MAPL_GridCompSetEntryPoint ( gc, ESMF_METHOD_RUN, run, _RC) + +#include "ACG_Export___.h" +#include "ACG_Import___.h" + +! Set generic services +! ---------------------------------- + call MAPL_GenericSetServices(GC, _RC) + + _RETURN(_SUCCESS) + + end subroutine SetServices + +!------------------------------------------------------------------------------ +!> +! `initialize` is meant to initialize the `ACG` gridded component. +! It primarily creates its exports. +! + subroutine initialize(GC, import, export, clock, rc) + + type (ESMF_GridComp), intent(inout) :: GC !! Gridded component + type (ESMF_State), intent(inout) :: import !! Import state + type (ESMF_State), intent(inout) :: export !! Export state + type (ESMF_Clock), intent(inout) :: clock !! The clock + integer, optional, intent( out) :: RC !! Error code +! +! Locals + integer :: status + + call MAPL_GridCreate(GC, _RC) + +! Call Generic Initialize +! ---------------------------------------- + call MAPL_GenericInitialize(GC, import, export, clock, _RC) + + _RETURN(_SUCCESS) + + end subroutine initialize + +!------------------------------------------------------------------------------ +!> +! `run` is the Run method for `ACG`. +! + subroutine run(GC, import, export, clock, rc) + + type (ESMF_GridComp), intent(inout) :: GC !! Gridded component + type (ESMF_State), intent(inout) :: import !! Import state + type (ESMF_State), intent(inout) :: export !! Export state + type (ESMF_Clock), intent(inout) :: clock !! The clock + integer, optional, intent( out) :: RC !! Error code +! +! Locals + type (MAPL_MetaComp), pointer :: MAPL + integer :: status + +#include "ACG_DeclarePointer___.h" + +!**************************************************************************** +! Begin... + + ! Get my internal MAPL_Generic state + ! ----------------------------------- + call MAPL_GetObjectFromGC ( GC, MAPL, _RC) + +#include "ACG_GetPointer___.h" + + + _RETURN(_SUCCESS) + + _UNUSED_DUMMY(import) + _UNUSED_DUMMY(clock) + + end subroutine run + +end module ACG_GridComp diff --git a/docs/tutorial/grid_comps/automatic_code_generator_example/ACG_StateSpecs.rc b/docs/tutorial/grid_comps/automatic_code_generator_example/ACG_StateSpecs.rc new file mode 100644 index 000000000000..386d1f122034 --- /dev/null +++ b/docs/tutorial/grid_comps/automatic_code_generator_example/ACG_StateSpecs.rc @@ -0,0 +1,57 @@ +schema_version: 2.0.0 +component: ACG + +category: IMPORT +#---------------------------------------------------------------------------- +# VARIABLE | DIMENSIONS | Additional Metadata +#---------------------------------------------------------------------------- + NAME | UNITS | DIMS | VLOC | RESTART | LONG NAME +#---------------------------------------------------------------------------- + ZLE | m | xyz | E | | geopotential_height + T | K | xyz | C | OPT | air_temperature + PLE | Pa | xyz | E | OPT | air_pressure + +category: EXPORT +#--------------------------------------------------------------------------- +# VARIABLE | DIMENSIONS | Additional Metadata +#--------------------------------------------------------------------------- + NAME | UNITS | DIMS | VLOC | LONG NAME +#--------------------------------------------------------------------------- + ZPBLCN | m | xy | N | boundary_layer_depth + CNV_FRC | 1 | xy | N | convective_fraction + +category: INTERNAL +#--------------------------------------------------------------------------- +# VARIABLE | DIMENSION | Additional Metadata +#--------------------------------------------------------------------------- + NAME | UNITS | DIMS | VLOC | ADD2EXPORT | FRIENDLYTO | LONG NAME +#--------------------------------------------------------------------------- + + +#******************************************************** +# +# Legend +# +#------------------------------------------------------------------ +# Column label | MAPL keyword/interpretation | Default +#--------------|--------------------------------------------------- +# NAME | short_name | +# UNITS | units | +# DIMS | dims | +# VLOC | VLocation | MAPL_VLocationNone +# LONG NAME | long_name | +# COND | if () then | .FALSE. +# NUM_SUBTILES | num_subtiles +# ... +#------------------------------------------------------------------ +# +#-------------------------------------------- +# Entry alias | Column | MAPL keyword/interpretation +#--------------|----------------------------- +# xyz | DIMS | MAPL_HorzVert +# xy | DIMS | MAPL_HorzOnly +# z | DIMS | MAPL_VertOnly (plus ungridded) +# C | VLOC | MAPL_VlocationCenter +# E | VLOC | MAPL_VlocationEdge +# N | VLOC | MAPL_VlocationNone +#-------------------------------------------- diff --git a/docs/tutorial/grid_comps/automatic_code_generator_example/CMakeLists.txt b/docs/tutorial/grid_comps/automatic_code_generator_example/CMakeLists.txt new file mode 100644 index 000000000000..4ae20760f332 --- /dev/null +++ b/docs/tutorial/grid_comps/automatic_code_generator_example/CMakeLists.txt @@ -0,0 +1,21 @@ +esma_set_this (OVERRIDE MAPL.acg) + +set (srcs + ACG_GridComp.F90 + ) + +esma_add_library (${this} SRCS ${srcs} DEPENDENCIES MAPL.shared MAPL TYPE ${MAPL_LIBRARY_TYPE}) + +target_link_libraries(${this} PRIVATE esmf) + +target_include_directories (${this} PUBLIC $) + +set_target_properties (${this} PROPERTIES Fortran_MODULE_DIRECTORY ${include_${this}}) + +mapl_acg (${this} ACG_StateSpecs.rc + IMPORT_SPECS EXPORT_SPECS + GET_POINTERS DECLARE_POINTERS) + +if (NOT CMAKE_Fortran_COMPILER_ID MATCHES "NAG") + target_link_libraries(${this} PRIVATE OpenMP::OpenMP_Fortran) +endif () diff --git a/docs/user_guide/docs/mapl_code_generator.md b/docs/user_guide/docs/mapl_code_generator.md index d832986e37c4..fd3766fb614a 100644 --- a/docs/user_guide/docs/mapl_code_generator.md +++ b/docs/user_guide/docs/mapl_code_generator.md @@ -7,7 +7,7 @@ The number of the those variables can be large and make the declaration process MAPL has a utility tool (named [MAPL_GridCompSpecs_ACG.py ](https://github.com/GEOS-ESM/MAPL/blob/main/Apps/MAPL_GridCompSpecs_ACG.py)) that simplifies and facilitates the registration and access of member variables of the various states (Export, Import, and Internal) of gridded components. -The tool relies on a formatted ASCII file (`spec`` file) to autmatically generate, at compilation time, include files that have the necessary code segments for defining and accessing the expected state member variables. +The tool relies on a formatted ASCII file (`spec` file) to autmatically generate, at compilation time, include files that have the necessary code segments for defining and accessing the expected state member variables. In this document, we describe the [steps](https://github.com/GEOS-ESM/MAPL/wiki/Setting-Up-MAPL-Automatic-Code-Generator) to follow to use the tool. To simplify this documents, we use the words _Imports_, _Exports_ and _Internals_ to refer to member variables of the Import, Export and Internal states, respectively. @@ -138,6 +138,7 @@ Assume that we create such a file (that we name `MyComponent_StateSpecs.rc`) and ``` +schema_version: 2.0.0 component: MyComponent category: IMPORT @@ -196,7 +197,12 @@ category: INTERNAL #-------------------------------------------- ``` -Running `MAPL_GridCompSpecs_ACG.py` on the file `MyComponent_StateSpecs.rc` generates at compilation time four (4) includes files: +#### Remark +It is required to have the settings for the two variable `schema_version` (here `2.0.0`) +and `component` (here `MyComponent`) on top of the `spec` file. + + +Running `MAPL_GridCompSpecs_ACG.py` on the file `MyComponent_StateSpecs.rc` generates at compilation time four (4) include files: 1. `MyComponent_Export___.h` for the `MAPL_AddExportSpec` calls in the `SetServices` routine: @@ -307,6 +313,16 @@ mapl_acg (${this} MyComponent_StateSpecs.rc Note, if in your case, there is no Internal state, `INTERNAL_SPECS` needs not to be added in the above command. But there is no harm including it. +### Sample code +We provide a sample code (gridded component module, `spec` and `CMakeLists.txt` files) that shows +how the automatic code generator is used. The code is available at: + +``` + docs/tutorial/grid_comps/automatic_code_generator_example +``` + +The code is provided for illustration only and compiled with the entire MAPL package. + ### Future Work A future version of the tool will support a YAML specification file. diff --git a/gridcomps/ExtData2G/CMakeLists.txt b/gridcomps/ExtData2G/CMakeLists.txt index 96aff8e544f0..97d1e5d41c92 100644 --- a/gridcomps/ExtData2G/CMakeLists.txt +++ b/gridcomps/ExtData2G/CMakeLists.txt @@ -23,7 +23,7 @@ set (srcs ) -esma_add_library (${this} SRCS ${srcs} DEPENDENCIES MAPL.shared MAPL.base MAPL.generic MAPL.griddedio TYPE SHARED) +esma_add_library (${this} SRCS ${srcs} DEPENDENCIES MAPL.shared MAPL.base MAPL.generic MAPL.griddedio TYPE ${MAPL_LIBRARY_TYPE}) target_link_libraries (${this} PUBLIC GFTL::gftl GFTL_SHARED::gftl-shared esmf NetCDF::NetCDF_Fortran PRIVATE MPI::MPI_Fortran) target_include_directories (${this} PUBLIC $) diff --git a/gridcomps/ExtData2G/ExtDataGridCompNG.F90 b/gridcomps/ExtData2G/ExtDataGridCompNG.F90 index 0da681786a62..2a14407b4ced 100644 --- a/gridcomps/ExtData2G/ExtDataGridCompNG.F90 +++ b/gridcomps/ExtData2G/ExtDataGridCompNG.F90 @@ -387,14 +387,16 @@ SUBROUTINE Initialize_ ( GC, IMPORT, EXPORT, CLOCK, rc ) do j=1,num_rules num_primary=num_primary+1 write(sidx,'(I1)')j - call config_yaml%fillin_primary(current_base_name//"+"//sidx,current_base_name,self%primary%item(num_primary),time,clock,_RC) + call config_yaml%fillin_primary(current_base_name//"+"//sidx,current_base_name,self%primary%item(num_primary),time,clock,rc=status) + _ASSERT(status==0, "ExtData multi-rule problem with BASE NAME "//TRIM(current_base_name)) allocate(self%primary%item(num_primary)%start_end_time(2)) self%primary%item(num_primary)%start_end_time(1)=time_ranges(j) self%primary%item(num_primary)%start_end_time(2)=time_ranges(j+1) enddo else num_primary=num_primary+1 - call config_yaml%fillin_primary(current_base_name,current_base_name,self%primary%item(num_primary),time,clock,_RC) + call config_yaml%fillin_primary(current_base_name,current_base_name,self%primary%item(num_primary),time,clock,rc=status) + _ASSERT(status==0, "ExtData single-rule problem with BASE NAME "//TRIM(current_base_name)) end if call ESMF_StateGet(Export,current_base_name,state_item_type,_RC) if (state_item_type /= ESMF_STATEITEM_NOTFOUND) then @@ -1757,7 +1759,7 @@ function get_item_index(this,base_name,current_time,rc) result(item_index) exit end if enddo - _ASSERT(found,"no item with that basename found") + _ASSERT(found,"ExtData no item with basename '"//TRIM(base_name)//"' found") item_index = -1 if (num_rules == 1) then @@ -1771,7 +1773,7 @@ function get_item_index(this,base_name,current_time,rc) result(item_index) endif enddo end if - _ASSERT(item_index/=-1,"did not find item") + _ASSERT(item_index/=-1,"ExtData did not find item index for basename "//TRIM(base_name)) _RETURN(_SUCCESS) end function get_item_index diff --git a/gridcomps/ExtData2G/ExtDataOldTypesCreator.F90 b/gridcomps/ExtData2G/ExtDataOldTypesCreator.F90 index 1abc5720c795..af453e701f5b 100644 --- a/gridcomps/ExtData2G/ExtDataOldTypesCreator.F90 +++ b/gridcomps/ExtData2G/ExtDataOldTypesCreator.F90 @@ -142,7 +142,13 @@ subroutine fillin_primary(this,item_name,base_name,primary_item,time,clock,unusa ! file_template primary_item%isConst = .false. if (index(rule%collection,"/dev/null")==0) then - dataset => this%file_stream_map%at(trim(rule%collection)) + + if ( ASSOCIATED(this%file_stream_map%at(trim(rule%collection))) ) then + dataset => this%file_stream_map%at(trim(rule%collection)) + else + _FAIL("ExtData problem with collection "//TRIM(rule%collection)) + end if + primary_item%file_template = dataset%file_template get_range = trim(time_sample%extrap_outside) /= "none" call dataset%detect_metadata(primary_item%file_metadata,time,rule%multi_rule,get_range=get_range,_RC) diff --git a/gridcomps/History/CMakeLists.txt b/gridcomps/History/CMakeLists.txt index c00678537245..25ba48139cfe 100644 --- a/gridcomps/History/CMakeLists.txt +++ b/gridcomps/History/CMakeLists.txt @@ -5,6 +5,7 @@ set (srcs MAPL_HistoryTrajectoryMod_smod.F90 MAPL_HistoryCollection.F90 MAPL_HistoryGridComp.F90 + MAPL_EpochSwathMod.F90 MAPL_StationSamplerMod.F90 ) diff --git a/gridcomps/History/MAPL_EpochSwathMod.F90 b/gridcomps/History/MAPL_EpochSwathMod.F90 new file mode 100644 index 000000000000..62b94145df5f --- /dev/null +++ b/gridcomps/History/MAPL_EpochSwathMod.F90 @@ -0,0 +1,1234 @@ +! +! __ Analogy to GriddedIO.F90 with a twist for Epoch Swath grid +! +#include "MAPL_Generic.h" + +module MAPL_EpochSwathMod + use ESMF + use ESMFL_Mod + use MAPL_AbstractGridFactoryMod + use MAPL_AbstractRegridderMod + use MAPL_GridManagerMod + use MAPL_BaseMod + use MAPL_NewRegridderManager + use MAPL_RegridMethods + use MAPL_TimeDataMod + use MAPL_VerticalDataMod + use MAPL_Constants + use pFIO + use MAPL_GriddedIOItemVectorMod + use MAPL_GriddedIOItemMod + use MAPL_ExceptionHandling + use pFIO_ClientManagerMod + use MAPL_DataCollectionMod + use MAPL_DataCollectionManagerMod + use gFTL_StringVector + use gFTL_StringStringMap + use MAPL_StringGridMapMod + use MAPL_FileMetadataUtilsMod + use MAPL_DownbitMod + use Plain_netCDF_Time + use, intrinsic :: ISO_C_BINDING + use, intrinsic :: iso_fortran_env, only: REAL64 + use ieee_arithmetic, only: isnan => ieee_is_nan + implicit none + + private + + type, public :: samplerHQ + type(ESMF_Clock) :: clock + type(ESMF_Alarm) :: alarm + type(ESMF_Time) :: RingTime + type(ESMF_TimeInterval) :: Frequency_epoch + type(ESMF_config) :: config_grid_save + type(ESMF_grid) :: ogrid + character(len=ESMF_MAXSTR) :: grid_type + real*8 :: arr(2) + + contains + procedure :: create_grid + procedure :: regrid_accumulate => regrid_accumulate_on_xysubset + procedure :: destroy_rh_regen_ogrid + procedure :: fill_time_in_bundle + end type samplerHQ + + interface samplerHQ + module procedure new_samplerHQ + end interface samplerHQ + + type, public :: sampler + type(FileMetaData), allocatable :: metadata + type(fileMetadataUtils), pointer :: current_file_metadata + integer :: write_collection_id + integer :: read_collection_id + integer :: metadata_collection_id + class (AbstractRegridder), pointer :: regrid_handle => null() + type(ESMF_Grid) :: output_grid + logical :: doVertRegrid = .false. + type(ESMF_FieldBundle) :: output_bundle + type(ESMF_FieldBundle) :: input_bundle + type(ESMF_FieldBundle) :: acc_bundle + type(ESMF_Time) :: startTime + integer :: regrid_method = REGRID_METHOD_BILINEAR + integer :: nbits_to_keep = MAPL_NBITS_NOT_SET + real, allocatable :: lons(:,:),lats(:,:) + real, allocatable :: corner_lons(:,:),corner_lats(:,:) + real, allocatable :: times(:) + type(TimeData) :: timeInfo + type(VerticalData) :: vdata + type(GriddedIOitemVector) :: items + integer :: deflateLevel = 0 + integer :: quantizeAlgorithm = 1 + integer :: quantizeLevel = 0 + integer, allocatable :: chunking(:) + logical :: itemOrderAlphabetical = .true. + integer :: fraction + logical :: have_initalized + contains +!! procedure :: CreateFileMetaData + procedure :: Create_bundle_RH + procedure :: CreateVariable + procedure :: regridScalar + procedure :: regridVector + procedure :: set_param + procedure :: set_default_chunking + procedure :: check_chunking + procedure :: alphabatize_variables + procedure :: addVariable_to_acc_bundle + procedure :: addVariable_to_output_bundle + procedure :: interp_accumulate_fields + end type sampler + + interface sampler + module procedure new_sampler + end interface sampler + +contains + + function new_samplerHQ(clock, config, key, rc) result(hq) + implicit none + type(samplerHQ) :: hq + type(ESMF_Clock), intent(in) :: clock + type(ESMF_Config), intent(inout) :: config + character(len=*), intent(in) :: key + integer, optional, intent(out) :: rc + + character(len=ESMF_MAXSTR) :: time_string + integer :: status + integer :: time_integer + type(ESMF_Time) :: RingTime_epoch + type(ESMF_Time) :: startTime + type(ESMF_Time) :: currTime + type(ESMF_TimeInterval) :: timeStep + type(ESMF_TimeInterval) :: Frequency_epoch + + integer :: sec, second + integer :: n1 + type(ESMF_Config) :: cf + + + hq%clock= clock + hq%config_grid_save= config + + hq%arr(1:2) = -2.d0 + call ESMF_ClockGet ( clock, CurrTime=currTime, _RC ) + call ESMF_ClockGet ( clock, timestep=timestep, _RC ) + call ESMF_ClockGet ( clock, startTime=startTime, _RC ) + call ESMF_ConfigGetAttribute(config, value=time_integer, label=trim(key)//'.Epoch:', default=0, _RC) + _ASSERT(time_integer /= 0, 'Epoch value in config wrong') + second = hms_2_s (time_integer) + call ESMF_TimeIntervalSet(frequency_epoch, s=second, _RC) + hq%frequency_epoch = frequency_epoch + hq%RingTime = currTime + hq%alarm = ESMF_AlarmCreate( clock=clock, RingInterval=Frequency_epoch, & + RingTime=hq%RingTime, sticky=.false., _RC ) + + _RETURN(_SUCCESS) + + end function new_samplerHQ + + + !--------------------------------------------------! + ! __ set + ! - ogrid via grid_manager%make_grid + ! using currTime and HQ%config_grid_save + !--------------------------------------------------! + function create_grid(this, key, currTime, grid_type, rc) result(ogrid) + type (ESMF_Grid) :: ogrid + class(samplerHQ) :: this + character(len=*), intent(in) :: key + type(ESMF_Time), intent(inout) :: currTime + character(len=*), optional, intent(in) :: grid_type + integer, intent(out), optional :: rc + integer :: status + + type(ESMF_Config) :: config_grid + character(len=ESMF_MAXSTR) :: time_string + + logical :: ispresent + + if (present(grid_type)) this%grid_type = trim(grid_type) + config_grid = this%config_grid_save + call ESMF_TimeGet(currTime, timeString=time_string, _RC) + ! + ! -- the `ESMF_ConfigSetAttribute` shows a risk + ! to overwrite the nextline in config + ! + call ESMF_ConfigSetAttribute( config_grid, trim(time_string), label=trim(key)//'.Epoch_init:', _RC) + ogrid = grid_manager%make_grid(config_grid, prefix=trim(key)//'.', _RC ) + this%ogrid = ogrid + _RETURN(_SUCCESS) + + end function create_grid + + + subroutine regrid_accumulate_on_xysubset (this, sp, rc) + class(samplerHQ) :: this + class(sampler), intent(inout) :: sp + integer, intent(out), optional :: rc + integer :: status + + class(AbstractGridFactory), pointer :: factory + integer :: xy_subset(2,2) + type(ESMF_Time) :: timeset(2) + type(ESMF_Time) :: current_time + type(ESMF_TimeInterval) :: dur + character(len=ESMF_MAXSTR) :: time_string + + integer, allocatable :: global_xy_mask(:,:) + integer, allocatable :: local_xy_mask(:,:) + + integer :: counts(5) + integer :: dims(3) + integer :: m1, m2 + + ! __ s1. get xy_subset + + factory => grid_manager%get_factory(this%ogrid,_RC) + call ESMF_ClockGet(this%clock,currTime=current_time,_RC) + call ESMF_ClockGet(this%clock,timeStep=dur, _RC ) + timeset(1) = current_time - dur + timeset(2) = current_time + call factory%get_xy_subset( timeset, xy_subset, _RC) + + ! __ s2. interpolate then save data using xy_mask + + call sp%interp_accumulate_fields (xy_subset, _RC) + + _RETURN(ESMF_SUCCESS) + + end subroutine regrid_accumulate_on_xysubset + + + subroutine destroy_rh_regen_ogrid (this, key_grid_label, output_grids, sp, rc) + implicit none + class(samplerHQ) :: this + class(sampler) :: sp + type (StringGridMap), target, intent(inout) :: output_grids + character(len=*), intent(in) :: key_grid_label + integer, intent(out), optional :: rc + integer :: status + + class(AbstractGridFactory), pointer :: factory + type(ESMF_Time) :: currTime + type(ESMF_TimeInterval) :: dur + character(len=ESMF_MAXSTR) :: time_string + + type(ESMF_Grid), pointer :: pgrid + type(ESMF_Grid) :: ogrid + type(ESMF_Grid) :: input_grid + character(len=ESMF_MAXSTR) :: key_str + type (StringGridMapIterator) :: iter + character(len=:), pointer :: key + type (ESMF_Config) :: config_grid + + integer :: i, numVars + character(len=ESMF_MAXSTR), allocatable :: names(:) + type(ESMF_Field) :: field + + if ( .NOT. ESMF_AlarmIsRinging(this%alarm) ) then + write(6,*) 'ck: regen, not in alarming' + rc=0 + return + endif + + + !__ s1. destroy ogrid + regen ogrid + + key_str=trim(key_grid_label) + pgrid => output_grids%at(trim(key_grid_label)) + call grid_manager%destroy(pgrid,_RC) + + call ESMF_ClockGet (this%clock, CurrTime=currTime, _RC ) + iter = output_grids%begin() + do while (iter /= output_grids%end()) + key => iter%key() + if (trim(key)==trim(key_str)) then + ogrid = this%create_grid (key_str, currTime, _RC) + call output_grids%set(key, ogrid) + this%ogrid = ogrid + endif + call iter%next() + enddo + + + !__ s2. destroy RH + + call sp%regrid_handle%destroy(_RC) + + + !__ s3. destroy acc_bundle / output_bundle + + call ESMF_FieldBundleGet(sp%acc_bundle,fieldCount=numVars,_RC) + allocate(names(numVars),stat=status) + call ESMF_FieldBundleGet(sp%acc_bundle,fieldNameList=names,_RC) + do i=1,numVars + call ESMF_FieldBundleGet(sp%acc_bundle,trim(names(i)),field=field,_RC) + call ESMF_FieldDestroy(field,noGarbage=.true., _RC) + enddo + call ESMF_FieldBundleDestroy(sp%acc_bundle,noGarbage=.true.,_RC) + + call ESMF_FieldBundleGet(sp%output_bundle,fieldCount=numVars,_RC) + allocate(names(numVars),stat=status) + call ESMF_FieldBundleGet(sp%output_bundle,fieldNameList=names,_RC) + do i=1,numVars + call ESMF_FieldBundleGet(sp%output_bundle,trim(names(i)),field=field,_RC) + call ESMF_FieldDestroy(field,noGarbage=.true., _RC) + enddo + call ESMF_FieldBundleDestroy(sp%output_bundle,noGarbage=.true.,_RC) + + _RETURN(ESMF_SUCCESS) + + end subroutine destroy_rh_regen_ogrid + + + subroutine fill_time_in_bundle (this, xname, bundle, rc) + implicit none + class(samplerHQ) :: this + character(len=*), intent(in) :: xname + type(ESMF_FieldBundle), intent(inout) :: bundle + integer, optional, intent(out) :: rc + integer :: status + + class(AbstractGridFactory), pointer :: factory + type(ESMF_Field) :: field + real(kind=ESMF_KIND_R4), pointer :: ptr2d(:,:) + + ! __ get field xname='time' + call ESMF_FieldBundleGet (bundle, xname, field=field, _RC) + call ESMF_FieldGet (field, farrayptr=ptr2d, _RC) + + ! __ obs_time from swath factory + factory => grid_manager%get_factory(this%ogrid,_RC) + call factory%get_obs_time (this%ogrid, ptr2d, _RC) + + _RETURN(ESMF_SUCCESS) + + end subroutine fill_time_in_bundle + + + function new_sampler(metadata,input_bundle,output_bundle,write_collection_id,read_collection_id, & + metadata_collection_id,regrid_method,fraction,items,rc) result(GriddedIO) + type(sampler) :: GriddedIO + type(Filemetadata), intent(in), optional :: metadata + type(ESMF_FieldBundle), intent(in), optional :: input_bundle + type(ESMF_FieldBundle), intent(in), optional :: output_bundle + integer, intent(in), optional :: write_collection_id + integer, intent(in), optional :: read_collection_id + integer, intent(in), optional :: metadata_collection_id + integer, intent(in), optional :: regrid_method + integer, intent(in), optional :: fraction + type(GriddedIOitemVector), intent(in), optional :: items + integer, intent(out), optional :: rc + + if (present(metadata)) GriddedIO%metadata=metadata + if (present(input_bundle)) GriddedIO%input_bundle=input_bundle + if (present(output_bundle)) GriddedIO%output_bundle=output_bundle + if (present(regrid_method)) GriddedIO%regrid_method=regrid_method + if (present(write_collection_id)) GriddedIO%write_collection_id=write_collection_id + if (present(read_collection_id)) GriddedIO%read_collection_id=read_collection_id + if (present(metadata_collection_id)) GriddedIO%metadata_collection_id=metadata_collection_id + if (present(items)) GriddedIO%items=items + if (present(fraction)) GriddedIO%fraction=fraction + _RETURN(ESMF_SUCCESS) + end function new_sampler + + + subroutine Create_bundle_RH(this,items,bundle,timeInfo,vdata,ogrid,global_attributes,rc) + class (sampler), intent(inout) :: this + type(GriddedIOitemVector), target, intent(inout) :: items + type(ESMF_FieldBundle), intent(inout) :: bundle + type(TimeData), optional, intent(inout) :: timeInfo + type(VerticalData), intent(inout), optional :: vdata + type (ESMF_Grid), intent(inout), pointer, optional :: ogrid + type(StringStringMap), target, intent(in), optional :: global_attributes + integer, intent(out), optional :: rc + + type(ESMF_Grid) :: input_grid + class (AbstractGridFactory), pointer :: factory + + type(ESMF_Field) :: new_field + type(GriddedIOitemVectorIterator) :: iter + type(GriddedIOitem), pointer :: item + type(stringVector) :: order + integer :: metadataVarsSize + type(StringStringMapIterator) :: s_iter + character(len=:), pointer :: attr_name, attr_val + integer :: status + + this%items = items + this%input_bundle = bundle + this%output_bundle = ESMF_FieldBundleCreate(rc=status) + _VERIFY(status) + if(present(timeInfo)) this%timeInfo = timeInfo + call ESMF_FieldBundleGet(this%input_bundle,grid=input_grid,rc=status) + _VERIFY(status) + if (present(ogrid)) then + this%output_grid=ogrid + else + call ESMF_FieldBundleGet(this%input_bundle,grid=this%output_grid,rc=status) + _VERIFY(status) + end if + this%regrid_handle => new_regridder_manager%make_regridder(input_grid,this%output_grid,this%regrid_method,rc=status) + _VERIFY(status) + + ! We get the regrid_method here because in the case of Identity, we set it to + ! REGRID_METHOD_IDENTITY in the regridder constructor if identity. Now we need + ! to change the regrid_method in the GriddedIO object to be the same as the + ! the regridder object. + this%regrid_method = this%regrid_handle%get_regrid_method() + + call ESMF_FieldBundleSet(this%output_bundle,grid=this%output_grid,rc=status) + _VERIFY(status) + factory => get_factory(this%output_grid,rc=status) + _VERIFY(status) + + ! __ please note, metadata in this section is not used in put_var to netCDF + ! the design used mGriddedIO%metadata in MAPL_HistoryGridComp.F90 + ! In other words, factory%append_metadata appeared here and in GriddedIO.F90 + ! + if (allocated(this%metadata)) then + deallocate (this%metadata) + end if + allocate(this%metadata) + call factory%append_metadata(this%metadata) + if (present(vdata)) then + this%vdata=vdata + else + this%vdata=VerticalData(rc=status) + _VERIFY(status) + end if + + call this%vdata%append_vertical_metadata(this%metadata,this%input_bundle,rc=status) + _VERIFY(status) + this%doVertRegrid = (this%vdata%regrid_type /= VERTICAL_METHOD_NONE) + if (this%vdata%regrid_type == VERTICAL_METHOD_ETA2LEV) call this%vdata%get_interpolating_variable(this%input_bundle,rc=status) + _VERIFY(status) + + iter = this%items%begin() + do while (iter /= this%items%end()) + item => iter%get() + if (item%itemType == ItemTypeScalar) then + call this%CreateVariable(item%xname,rc=status) + _VERIFY(status) + else if (item%itemType == ItemTypeVector) then + call this%CreateVariable(item%xname,rc=status) + _VERIFY(status) + call this%CreateVariable(item%yname,rc=status) + _VERIFY(status) + end if + call iter%next() + enddo + + + ! __ add acc_bundle and output_bundle + ! + this%acc_bundle = ESMF_FieldBundleCreate(_RC) + call ESMF_FieldBundleSet(this%acc_bundle,grid=this%output_grid,_RC) + iter = this%items%begin() + do while (iter /= this%items%end()) + item => iter%get() + call this%addVariable_to_acc_bundle(item%xname,_RC) + if (item%itemType == ItemTypeVector) then + call this%addVariable_to_acc_bundle(item%yname,_RC) + end if + call iter%next() + enddo + + + ! __ add time to acc_bundle + ! + new_field = ESMF_FieldCreate(this%output_grid ,name='time', & + typekind=ESMF_TYPEKIND_R4,_RC) + call MAPL_FieldBundleAdd( this%acc_bundle, new_field, _RC ) + + + _RETURN(_SUCCESS) + end subroutine Create_Bundle_RH + + + subroutine set_param(this,deflation,quantize_algorithm,quantize_level,chunking,nbits_to_keep,regrid_method,itemOrder,write_collection_id,rc) + class (sampler), intent(inout) :: this + integer, optional, intent(in) :: deflation + integer, optional, intent(in) :: quantize_algorithm + integer, optional, intent(in) :: quantize_level + integer, optional, intent(in) :: chunking(:) + integer, optional, intent(in) :: nbits_to_keep + integer, optional, intent(in) :: regrid_method + logical, optional, intent(in) :: itemOrder + integer, optional, intent(in) :: write_collection_id + integer, optional, intent(out) :: rc + + integer :: status + + if (present(regrid_method)) this%regrid_method=regrid_method + if (present(nbits_to_keep)) this%nbits_to_keep=nbits_to_keep + if (present(deflation)) this%deflateLevel = deflation + if (present(quantize_algorithm)) this%quantizeAlgorithm = quantize_algorithm + if (present(quantize_level)) this%quantizeLevel = quantize_level + if (present(chunking)) then + allocate(this%chunking,source=chunking,stat=status) + _VERIFY(status) + end if + if (present(itemOrder)) this%itemOrderAlphabetical = itemOrder + if (present(write_collection_id)) this%write_collection_id=write_collection_id + _RETURN(ESMF_SUCCESS) + + end subroutine set_param + + subroutine set_default_chunking(this,rc) + class (sampler), intent(inout) :: this + integer, optional, intent(out) :: rc + + integer :: global_dim(3) + integer :: status + + call MAPL_GridGet(this%output_grid,globalCellCountPerDim=global_dim,rc=status) + _VERIFY(status) + if (global_dim(1)*6 == global_dim(2)) then + allocate(this%chunking(5)) + this%chunking(1) = global_dim(1) + this%chunking(2) = global_dim(1) + this%chunking(3) = 1 + this%chunking(4) = 1 + this%chunking(5) = 1 + else + allocate(this%chunking(4)) + this%chunking(1) = global_dim(1) + this%chunking(2) = global_dim(2) + this%chunking(3) = 1 + this%chunking(4) = 1 + endif + _RETURN(ESMF_SUCCESS) + + end subroutine set_default_chunking + + subroutine check_chunking(this,lev_size,rc) + class (sampler), intent(inout) :: this + integer, intent(in) :: lev_size + integer, optional, intent(out) :: rc + + integer :: global_dim(3) + integer :: status + character(len=5) :: c1,c2 + + call MAPL_GridGet(this%output_grid,globalCellCountPerDim=global_dim,rc=status) + _VERIFY(status) + if (global_dim(1)*6 == global_dim(2)) then + write(c2,'(I5)')global_dim(1) + write(c1,'(I5)')this%chunking(1) + _ASSERT(this%chunking(1) <= global_dim(1), "Chunk for Xdim "//c1//" must be less than or equal to "//c2) + write(c1,'(I5)')this%chunking(2) + _ASSERT(this%chunking(2) <= global_dim(1), "Chunk for Ydim "//c1//" must be less than or equal to "//c2) + _ASSERT(this%chunking(3) <= 6, "Chunksize for face dimension must be 6 or less") + if (lev_size > 0) then + write(c2,'(I5)')lev_size + write(c1,'(I5)')this%chunking(4) + _ASSERT(this%chunking(4) <= lev_size, "Chunk for level size "//c1//" must be less than or equal to "//c2) + end if + _ASSERT(this%chunking(5) == 1, "Time must have chunk size of 1") + else + write(c2,'(I5)')global_dim(1) + write(c1,'(I5)')this%chunking(1) + _ASSERT(this%chunking(1) <= global_dim(1), "Chunk for lon "//c1//" must be less than or equal to "//c2) + write(c2,'(I5)')global_dim(2) + write(c1,'(I5)')this%chunking(2) + _ASSERT(this%chunking(2) <= global_dim(2), "Chunk for lat "//c1//" must be less than or equal to "//c2) + if (lev_size > 0) then + write(c2,'(I5)')lev_size + write(c1,'(I5)')this%chunking(3) + _ASSERT(this%chunking(3) <= lev_size, "Chunk for level size "//c1//" must be less than or equal to "//c2) + end if + _ASSERT(this%chunking(4) == 1, "Time must have chunk size of 1") + endif + _RETURN(ESMF_SUCCESS) + + end subroutine check_chunking + + subroutine CreateVariable(this,itemName,rc) + class (sampler), intent(inout) :: this + character(len=*), intent(in) :: itemName + integer, optional, intent(out) :: rc + + integer :: status + + type(ESMF_Field) :: field,newField + class (AbstractGridFactory), pointer :: factory + integer :: fieldRank + logical :: isPresent + character(len=ESMF_MAXSTR) :: varName,longName,units + character(len=:), allocatable :: grid_dims + character(len=:), allocatable :: vdims + type(Variable) :: v + + call ESMF_FieldBundleGet(this%input_bundle,itemName,field=field,rc=status) + _VERIFY(status) + factory => get_factory(this%output_grid,rc=status) + _VERIFY(status) + + + call ESMF_FieldGet(field,rank=fieldRank,rc=status) + _VERIFY(status) + call ESMF_FieldGet(field,name=varName,rc=status) + _VERIFY(status) + call ESMF_AttributeGet(field,name="LONG_NAME",isPresent=isPresent,rc=status) + _VERIFY(status) + if ( isPresent ) then + call ESMF_AttributeGet (FIELD, NAME="LONG_NAME",VALUE=LongName, RC=STATUS) + _VERIFY(STATUS) + else + LongName = varName + endif + call ESMF_AttributeGet(field,name="UNITS",isPresent=isPresent,rc=status) + _VERIFY(status) + if ( isPresent ) then + call ESMF_AttributeGet (FIELD, NAME="UNITS",VALUE=units, RC=STATUS) + _VERIFY(STATUS) + else + units = 'unknown' + endif + + + ! finally make a new field if neccessary + if (this%doVertRegrid .and. (fieldRank ==3) ) then + newField = MAPL_FieldCreate(field,this%output_grid,lm=this%vData%lm,rc=status) + _VERIFY(status) + call MAPL_FieldBundleAdd(this%output_bundle,newField,rc=status) + _VERIFY(status) + else + newField = MAPL_FieldCreate(field,this%output_grid,rc=status) + _VERIFY(status) + call MAPL_FieldBundleAdd(this%output_bundle,newField,rc=status) + _VERIFY(status) + end if + _RETURN(_SUCCESS) + + end subroutine CreateVariable + + + subroutine RegridScalar(this,itemName,rc) + class (sampler), intent(inout) :: this + character(len=*), intent(in) :: itemName + integer, optional, intent(out) :: rc + + integer :: status + + type(ESMF_Field) :: field,outField + integer :: fieldRank + real, pointer :: ptr3d(:,:,:),outptr3d(:,:,:) + real, pointer :: ptr2d(:,:), outptr2d(:,:) + real, allocatable, target :: ptr3d_inter(:,:,:) + type(ESMF_Grid) :: gridIn,gridOut + logical :: hasDE_in, hasDE_out + logical :: first_entry + + call ESMF_FieldBundleGet(this%output_bundle,itemName,field=outField,rc=status) + _VERIFY(status) + call ESMF_FieldBundleGet(this%input_bundle,grid=gridIn,rc=status) + _VERIFY(status) + call ESMF_FieldBundleGet(this%output_bundle,grid=gridOut,rc=status) + _VERIFY(status) + hasDE_in = MAPL_GridHasDE(gridIn,rc=status) + _VERIFY(status) + hasDE_out = MAPL_GridHasDE(gridOut,rc=status) + _VERIFY(status) + first_entry = .true. + if (this%doVertRegrid) then + call ESMF_FieldBundleGet(this%input_bundle,itemName,field=field,rc=status) + _VERIFY(status) + call ESMF_FieldGet(Field,rank=fieldRank,rc=status) + _VERIFY(status) + if (fieldRank==3) then + if (hasDE_in) then + call ESMF_FieldGet(field,farrayPtr=ptr3d,rc=status) + _VERIFY(status) + else + allocate(ptr3d(0,0,0)) + end if + allocate(ptr3d_inter(size(ptr3d,1),size(ptr3d,2),this%vdata%lm),stat=status) + _VERIFY(status) + if (this%vdata%regrid_type==VERTICAL_METHOD_SELECT) then + call this%vdata%regrid_select_level(ptr3d,ptr3d_inter,rc=status) + _VERIFY(status) + else if (this%vdata%regrid_type==VERTICAL_METHOD_ETA2LEV) then + call this%vdata%regrid_eta_to_pressure(ptr3d,ptr3d_inter,rc=status) + _VERIFY(status) + else if (this%vdata%regrid_type==VERTICAL_METHOD_FLIP) then + call this%vdata%flip_levels(ptr3d,ptr3d_inter,rc=status) + _VERIFY(status) + end if + ptr3d => ptr3d_inter + end if + else + if (first_entry) then + nullify(ptr3d) + first_entry = .false. + end if + end if + + call ESMF_FieldBundleGet(this%input_bundle,itemName,field=field,rc=status) + _VERIFY(status) + call ESMF_FieldGet(field,rank=fieldRank,rc=status) + _VERIFY(status) + if (fieldRank==2) then + if (hasDE_in) then + call MAPL_FieldGetPointer(field,ptr2d,rc=status) + _VERIFY(status) + else + allocate(ptr2d(0,0)) + end if + if (hasDE_out) then + call MAPL_FieldGetPointer(OutField,outptr2d,rc=status) + _VERIFY(status) + else + allocate(outptr2d(0,0)) + end if + if (gridIn==gridOut) then + outPtr2d=ptr2d + else + if (this%regrid_method==REGRID_METHOD_FRACTION) ptr2d=ptr2d-this%fraction + call this%regrid_handle%regrid(ptr2d,outPtr2d,rc=status) + _VERIFY(status) + end if + +!! print *, maxval(ptr2d) +!! print *, minval(ptr2d) +!! print *, maxval(outptr2d) +!! print *, minval(outptr2d) + + else if (fieldRank==3) then + if (.not.associated(ptr3d)) then + if (hasDE_in) then + call ESMF_FieldGet(field,farrayPtr=ptr3d,rc=status) + _VERIFY(status) + else + allocate(ptr3d(0,0,0)) + end if + end if + if (hasDE_out) then + call MAPL_FieldGetPointer(OutField,outptr3d,rc=status) + _VERIFY(status) + else + allocate(outptr3d(0,0,0)) + end if + if (gridIn==gridOut) then + outPtr3d=Ptr3d + else + if (this%regrid_method==REGRID_METHOD_FRACTION) ptr3d=ptr3d-this%fraction + call this%regrid_handle%regrid(ptr3d,outPtr3d,rc=status) + _VERIFY(status) + end if + else + _FAIL('rank not supported') + end if + + if (allocated(ptr3d_inter)) deallocate(ptr3d_inter) + _RETURN(_SUCCESS) + + end subroutine RegridScalar + + subroutine RegridVector(this,xName,yName,rc) + class (sampler), intent(inout) :: this + character(len=*), intent(in) :: xName + character(len=*), intent(in) :: yName + integer, optional, intent(out) :: rc + + integer :: status + + type(ESMF_Field) :: xfield,xoutField + type(ESMF_Field) :: yfield,youtField + integer :: fieldRank + real, pointer :: xptr3d(:,:,:),xoutptr3d(:,:,:) + real, pointer :: xptr2d(:,:), xoutptr2d(:,:) + real, allocatable, target :: xptr3d_inter(:,:,:) + real, pointer :: yptr3d(:,:,:),youtptr3d(:,:,:) + real, pointer :: yptr2d(:,:), youtptr2d(:,:) + real, allocatable, target :: yptr3d_inter(:,:,:) + type(ESMF_Grid) :: gridIn, gridOut + logical :: hasDE_in, hasDE_out + + call ESMF_FieldBundleGet(this%output_bundle,xName,field=xoutField,rc=status) + _VERIFY(status) + call ESMF_FieldBundleGet(this%output_bundle,yName,field=youtField,rc=status) + _VERIFY(status) + call ESMF_FieldBundleGet(this%input_bundle,grid=gridIn,rc=status) + _VERIFY(status) + call ESMF_FieldBundleGet(this%output_bundle,grid=gridOut,rc=status) + _VERIFY(status) + hasDE_in = MAPL_GridHasDE(gridIn,rc=status) + _VERIFY(status) + hasDE_out = MAPL_GridHasDE(gridOut,rc=status) + _VERIFY(status) + + if (this%doVertRegrid) then + call ESMF_FieldBundleGet(this%input_bundle,xName,field=xfield,rc=status) + _VERIFY(status) + call ESMF_FieldGet(xField,rank=fieldRank,rc=status) + _VERIFY(status) + if (fieldRank==3) then + if (hasDE_in) then + call ESMF_FieldGet(xfield,farrayPtr=xptr3d,rc=status) + _VERIFY(status) + else + allocate(xptr3d(0,0,0)) + end if + allocate(xptr3d_inter(size(xptr3d,1),size(xptr3d,2),this%vdata%lm),stat=status) + _VERIFY(status) + if (this%vdata%regrid_type==VERTICAL_METHOD_SELECT) then + call this%vdata%regrid_select_level(xptr3d,xptr3d_inter,rc=status) + _VERIFY(status) + else if (this%vdata%regrid_type==VERTICAL_METHOD_ETA2LEV) then + call this%vdata%regrid_eta_to_pressure(xptr3d,xptr3d_inter,rc=status) + _VERIFY(status) + else if (this%vdata%regrid_type==VERTICAL_METHOD_FLIP) then + call this%vdata%flip_levels(xptr3d,xptr3d_inter,rc=status) + _VERIFY(status) + end if + xptr3d => xptr3d_inter + end if + call ESMF_FieldBundleGet(this%input_bundle,yName,field=yfield,rc=status) + _VERIFY(status) + call ESMF_FieldGet(yField,rank=fieldRank,rc=status) + _VERIFY(status) + if (fieldRank==3) then + if (hasDE_in) then + call ESMF_FieldGet(yfield,farrayPtr=yptr3d,rc=status) + _VERIFY(status) + else + allocate(yptr3d(0,0,0)) + end if + allocate(yptr3d_inter(size(yptr3d,1),size(yptr3d,2),this%vdata%lm),stat=status) + _VERIFY(status) + if (this%vdata%regrid_type==VERTICAL_METHOD_SELECT) then + call this%vdata%regrid_select_level(yptr3d,yptr3d_inter,rc=status) + _VERIFY(status) + else if (this%vdata%regrid_type==VERTICAL_METHOD_ETA2LEV) then + call this%vdata%regrid_eta_to_pressure(yptr3d,yptr3d_inter,rc=status) + _VERIFY(status) + else if (this%vdata%regrid_type==VERTICAL_METHOD_FLIP) then + call this%vdata%flip_levels(yptr3d,yptr3d_inter,rc=status) + _VERIFY(status) + end if + yptr3d => yptr3d_inter + end if + else + if (associated(xptr3d)) nullify(xptr3d) + if (associated(yptr3d)) nullify(yptr3d) + end if + + call ESMF_FieldBundleGet(this%input_bundle,xname,field=xfield,rc=status) + _VERIFY(status) + call ESMF_FieldBundleGet(this%input_bundle,yname,field=yfield,rc=status) + _VERIFY(status) + call ESMF_FieldGet(xfield,rank=fieldRank,rc=status) + _VERIFY(status) + if (fieldRank==2) then + if (hasDE_in) then + call MAPL_FieldGetPointer(xfield,xptr2d,rc=status) + _VERIFY(status) + call MAPL_FieldGetPointer(yfield,yptr2d,rc=status) + _VERIFY(status) + else + allocate(xptr2d(0,0)) + allocate(yptr2d(0,0)) + end if + + if (hasDE_in) then + call MAPL_FieldGetPointer(xOutField,xoutptr2d,rc=status) + _VERIFY(status) + call MAPL_FieldGetPointer(yOutField,youtptr2d,rc=status) + _VERIFY(status) + else + allocate(xoutptr2d(0,0)) + allocate(youtptr2d(0,0)) + end if + + + if (gridIn==gridOut) then + xoutPtr2d=xptr2d + youtPtr2d=yptr2d + else + call this%regrid_handle%regrid(xptr2d,yptr2d,xoutPtr2d,youtPtr2d,rc=status) + _VERIFY(status) + end if + else if (fieldRank==3) then + if (.not.associated(xptr3d)) then + if (hasDE_in) then + call MAPL_FieldGetPointer(xfield,xptr3d,rc=status) + _VERIFY(status) + else + allocate(xptr3d(0,0,0)) + end if + end if + if (.not.associated(yptr3d)) then + if (hasDE_in) then + call MAPL_FieldGetPointer(yfield,yptr3d,rc=status) + _VERIFY(status) + else + allocate(yptr3d(0,0,0)) + end if + end if + + if (hasDE_out) then + call MAPL_FieldGetPointer(xOutField,xoutptr3d,rc=status) + _VERIFY(status) + call MAPL_FieldGetPointer(yOutField,youtptr3d,rc=status) + _VERIFY(status) + else + allocate(xoutptr3d(0,0,0)) + allocate(youtptr3d(0,0,0)) + end if + + if (gridIn==gridOut) then + xoutPtr3d=xptr3d + youtPtr3d=yptr3d + else + call this%regrid_handle%regrid(xptr3d,yptr3d,xoutPtr3d,youtPtr3d,rc=status) + _VERIFY(status) + end if + end if + + if (allocated(xptr3d_inter)) deallocate(xptr3d_inter) + if (allocated(yptr3d_inter)) deallocate(yptr3d_inter) + _RETURN(_SUCCESS) + + end subroutine RegridVector + + + subroutine alphabatize_variables(this,nfixedVars,rc) + class (sampler), intent(inout) :: this + integer, intent(in) :: nFixedVars + integer, optional, intent(out) :: rc + + type(StringVector) :: order + type(StringVector) :: newOrder + character(len=:), pointer :: v1 + character(len=ESMF_MAXSTR) :: c1,c2 + character(len=ESMF_MAXSTR), allocatable :: temp(:) + logical :: swapped + integer :: n,i + integer :: status + + order = this%metadata%get_order(rc=status) + _VERIFY(status) + n = Order%size() + allocate(temp(nFixedVars+1:n)) + do i=1,n + v1 => order%at(i) + if ( i > nFixedVars) temp(i)=trim(v1) + enddo + + swapped = .true. + do while(swapped) + swapped = .false. + do i=nFixedVars+1,n-1 + c1 = temp(i) + c2 = temp(i+1) + if (c1 > c2) then + temp(i+1)=c1 + temp(i)=c2 + swapped =.true. + end if + enddo + enddo + + do i=1,nFixedVars + v1 => Order%at(i) + call newOrder%push_back(v1) + enddo + do i=nFixedVars+1,n + call newOrder%push_back(trim(temp(i))) + enddo + call this%metadata%set_order(newOrder,rc=status) + _VERIFY(status) + deallocate(temp) + + _RETURN(_SUCCESS) + + end subroutine alphabatize_variables + + + subroutine addVariable_to_acc_bundle(this,itemName,rc) + class (sampler), intent(inout) :: this + character(len=*), intent(in) :: itemName + integer, optional, intent(out) :: rc + + type(ESMF_Field) :: field,newField + type(ESMF_Array) :: array1 + real(KIND=ESMF_KIND_R4), pointer :: pt2d(:,:) + class (AbstractGridFactory), pointer :: factory + integer :: fieldRank + logical :: isPresent + integer :: status + + call ESMF_FieldBundleGet(this%input_bundle,itemName,field=field,_RC) + call ESMF_FieldGet(field,rank=fieldRank,rc=status) + if (this%doVertRegrid .and. (fieldRank ==3) ) then + newField = MAPL_FieldCreate(field,this%output_grid,lm=this%vData%lm,_RC) + else + newField = MAPL_FieldCreate(field,this%output_grid,_RC) + end if + call MAPL_FieldBundleAdd(this%acc_bundle,newField,_RC) + + _RETURN(_SUCCESS) + + end subroutine addVariable_to_acc_bundle + + + subroutine addVariable_to_output_bundle(this,itemName,rc) + class (sampler), intent(inout) :: this + character(len=*), intent(in) :: itemName + integer, optional, intent(out) :: rc + + type(ESMF_Field) :: field,newField + class (AbstractGridFactory), pointer :: factory + integer :: fieldRank + logical :: isPresent + integer :: status + + call ESMF_FieldBundleGet(this%input_bundle,itemName,field=field,_RC) + call ESMF_FieldGet(field,rank=fieldRank,rc=status) + if (this%doVertRegrid .and. (fieldRank ==3) ) then + newField = MAPL_FieldCreate(field,this%output_grid,lm=this%vData%lm,_RC) + else + newField = MAPL_FieldCreate(field,this%output_grid,_RC) + end if + call MAPL_FieldBundleAdd(this%output_bundle,newField,_RC) + + _RETURN(_SUCCESS) + end subroutine addVariable_to_output_bundle + + + + !! -- based on subroutine bundlepost(this,filename,oClients,rc) + !! + subroutine interp_accumulate_fields (this,xy_subset,rc) + implicit none + class (sampler) :: this + integer, intent(in) :: xy_subset(2,2) + !!integer, intent(in) :: xy_mask(:,:) + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_Field) :: outField + type(ESMF_Field) :: new_outField + type(ESMF_Grid) :: grid + integer :: tindex + type(ArrayReference) :: ref + + type(GriddedIOitemVectorIterator) :: iter + type(GriddedIOitem), pointer :: item + logical :: have_time + + type(ESMF_Array) :: array1, array2 + integer :: is,ie,js,je + + integer :: rank, rank1, rank2 + real(KIND=ESMF_KIND_R4), pointer :: pt2d(:,:), pt2d_(:,:) + real(KIND=ESMF_KIND_R4), pointer :: pt3d(:,:,:), pt3d_(:,:,:) + + integer :: localDe, localDECount + integer, dimension(:), allocatable :: LB, UB, exclusiveCount + integer, dimension(:), allocatable :: compLB, compUB, compCount + integer :: dimCount + integer :: y1, y2 + integer :: j, jj + integer :: ii1, iin, jj1, jjn + integer, dimension(:), allocatable :: j1, j2 + + is=xy_subset(1,1); ie=xy_subset(2,1) + js=xy_subset(1,2); je=xy_subset(2,2) + + if (this%vdata%regrid_type==VERTICAL_METHOD_ETA2LEV) then + call this%vdata%setup_eta_to_pressure(regrid_handle=this%regrid_handle,output_grid=this%output_grid,rc=status) + _VERIFY(status) + end if + + call ESMF_FieldBundleGet(this%output_bundle, grid=grid, _RC) + call ESMF_GridGet(grid, localDECount=localDECount, dimCount=dimCount, _RC) + allocate ( LB(dimCount), UB(dimCount), exclusiveCount(dimCount) ) + allocate ( compLB(dimCount), compUB(dimCount), compCount(dimCount) ) + + allocate ( j1(0:localDEcount-1) ) ! start + allocate ( j2(0:localDEcount-1) ) ! end + + _ASSERT ( localDEcount == 1, 'failed, due to localDEcount > 1') + call MAPL_GridGetInterior(grid,ii1,iin,jj1,jjn) +!! write(6,*) 'MAPL_GridGetInterior, ii1,iin,jj1,jjn', ii1,iin,jj1,jjn +!! print*, 'js,je ', js, je + + LB(1)=ii1; LB(2)=jj1 + UB(1)=iin; UB(2)=jjn + + do localDe=0, localDEcount-1 + ! + ! is/ie, js/je, [LB, UB] + ! + ! + y1=jj1; y2=jjn + if (y1 < js) then + if (y2 < js) then + j1(localDe)=-1 + j2(localDe)=-1 + elseif (y2 < je) then + j1(localDe)=js + j2(localDe)=y2 + else + j1(localDe)=js + j2(localDe)=je + endif + elseif (y1 <= je) then + j1(localDe)=y1 + if (y2 < je) then + j2(localDe)=y2 + else + j2(localDe)=je + endif + else + j1(localDe)=-1 + j2(localDe)=-1 + endif + enddo + +!! write(6,*) 'ck bundlepost_acc' +!! write(6,*) 'j1(localDe)', j1(0:localDeCount-1) +!! write(6,*) 'j2(localDe)', j2(0:localDeCount-1) + + + iter = this%items%begin() + do while (iter /= this%items%end()) + item => iter%get() + if (item%itemType == ItemTypeScalar) then + call this%RegridScalar(item%xname,rc=status) + _VERIFY(status) + call ESMF_FieldBundleGet(this%output_bundle,item%xname,field=outField, _RC) + _VERIFY(status) + if (this%vdata%regrid_type==VERTICAL_METHOD_ETA2LEV) then + call this%vdata%correct_topo(outField,rc=status) + _VERIFY(status) + end if + + ! -- mask the time interval + ! store the time interval fields into new bundle + call ESMF_FieldGet(outField, Array=array1, _RC) + call ESMF_FieldBundleGet(this%acc_bundle,item%xname,field=new_outField,_RC) + call ESMF_FieldGet(new_outField, Array=array2, _RC) + call ESMF_ArrayGet(array1, rank=rank, _RC) + if (rank==2) then + call ESMF_ArrayGet(array1, farrayptr=pt2d, _RC) + call ESMF_ArrayGet(array2, farrayptr=pt2d_, _RC) + localDe=0 + if (j1(localDe)>0) then + do j= j1(localDe), j2(localDe) + jj= j-jj1+1 ! j_local +!! write(6,*) 'j, jj', j, jj + pt2d_(:,jj) = pt2d(:,jj) + enddo + endif + elseif (rank==3) then + call ESMF_ArrayGet(array1, farrayptr=pt3d, _RC) + call ESMF_ArrayGet(array2, farrayptr=pt3d_, _RC) + do localDe=0, localDEcount-1 + if (j1(localDe)>0) then + do j= j1(localDe), j2(localDe) + jj= j-jj1+1 + pt3d_(:,jj,:) = pt3d(:,jj,:) + enddo + endif + enddo + else + _FAIL('failed interp_accumulate_fields') + endif + + else if (item%itemType == ItemTypeVector) then + _FAIL('ItemTypeVector not implemented') + end if + call iter%next() + enddo + + _RETURN(ESMF_SUCCESS) + + end subroutine interp_accumulate_fields + + + subroutine get_xy_mask(grid, xy_subset, xy_mask, rc) + implicit none + type(ESMF_Grid), intent(in) :: grid + integer, intent(in) :: xy_subset(2,2) + integer, intent(out) :: xy_mask(:,:) + integer, optional, intent(out) :: rc + + integer :: status + integer :: ii1, iin, jj1, jjn ! local box for localDE + integer :: is, ie, js, je ! global box for each time-interval + integer :: j1p, jnp ! local y-index for each time-interval + + integer :: dimCount + integer :: y1, y2 + integer :: j, jj + integer :: j1, j2 + + is=xy_subset(1,1); ie=xy_subset(2,1) + js=xy_subset(1,2); je=xy_subset(2,2) + + call MAPL_GridGetInterior(grid,ii1,iin,jj1,jjn) + write(6,*) 'MAPL_GridGetInterior, ii1,iin,jj1,jjn', ii1,iin,jj1,jjn + + y1=jj1; y2=jjn + if (y1 < js) then + if (y2 < js) then + j1=-1 + j2=-1 + elseif (y2 < je) then + j1=js + j2=y2 + else + j1=js + j2=je + endif + elseif (y1 <= je) then + j1=y1 + if (y2 < je) then + j2=y2 + else + j2=je + endif + else + j1=-1 + j2=-1 + endif + +!! write(6,*) 'get_xy_mask: j1,j2=', j1, j2 + xy_mask(:,:) = 0 + if (j1 > 0) then + do jj = j1, j2 + xy_mask(:, jj) = 1 + enddo + end if + + if(present(rc)) rc=0 + + end subroutine get_xy_mask + + +end module MAPL_EpochSwathMod diff --git a/gridcomps/History/MAPL_HistoryCollection.F90 b/gridcomps/History/MAPL_HistoryCollection.F90 index b5c0c3d35f19..0bc60881028d 100644 --- a/gridcomps/History/MAPL_HistoryCollection.F90 +++ b/gridcomps/History/MAPL_HistoryCollection.F90 @@ -11,6 +11,7 @@ module MAPL_HistoryCollectionMod use HistoryTrajectoryMod use StationSamplerMod use gFTL_StringStringMap + use MAPL_EpochSwathMod implicit none private @@ -55,6 +56,7 @@ module MAPL_HistoryCollectionMod integer,pointer :: expSTATE (:) integer :: unit type(ESMF_FieldBundle) :: bundle + type(sampler) :: xsampler type(MAPL_CFIO) :: MCFIO type(MAPL_GriddedIO) :: mGriddedIO type(VerticalData) :: vdata @@ -102,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 index d9356e9bb70c..a29439905bbd 100644 --- a/gridcomps/History/MAPL_HistoryGridComp.F90 +++ b/gridcomps/History/MAPL_HistoryGridComp.F90 @@ -58,6 +58,8 @@ module MAPL_HistoryGridCompMod use MAPL_TimeUtilsMod, only: is_valid_time, is_valid_date use gFTL_StringStringMap !use ESMF_CFIOMOD + use MAPL_EpochSwathMod + use pflogger, only: Logger, logging use mpi @@ -143,6 +145,8 @@ module MAPL_HistoryGridCompMod public HISTORY_ExchangeListWrap + type(samplerHQ) :: Hsampler + contains !===================================================================== @@ -484,7 +488,6 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) ! Read User-Supplied History Lists from Config File ! ------------------------------------------------- call ESMF_GridCompGet( gc, config=config, _RC ) - call ESMF_ConfigGetAttribute ( config, value=INTSTATE%expsrc, & label ='EXPSRC:', default='', _RC ) call ESMF_ConfigGetAttribute ( config, value=INTSTATE%expid, & @@ -577,7 +580,6 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) end if call ESMF_ConfigNextLine ( config,tableEnd=tend,_RC ) enddo - if (nlist == 0) then _RETURN(ESMF_SUCCESS) end if @@ -614,13 +616,18 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) call MAPL_ConfigSetAttribute(config, value=nx,label=trim(key)//".NX:",_RC) call MAPL_ConfigSetAttribute(config, value=ny,label=trim(key)//".NY:",_RC) end if - output_grid = grid_manager%make_grid(config, prefix=key//'.', _RC) + + if (trim(grid_type)/='Swath') then + output_grid = grid_manager%make_grid(config, prefix=key//'.', _RC) + else + Hsampler = samplerHQ(clock, config, key, _RC) + output_grid = Hsampler%create_grid(key, currTime, grid_type=grid_type, _RC) + end if call IntState%output_grids%set(key, output_grid) call iter%next() end do end block OUTPUT_GRIDS end if - if (intstate%version >= 2) then call ESMF_ConfigFindLabel(config, 'FIELD_SETS:', _RC) table_end = .false. @@ -634,7 +641,7 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) end if call ESMF_ConfigNextLine ( config,tableEnd=table_end,_RC ) enddo - + field_set_iter = intState%field_sets%begin() do while (field_set_iter /= intState%field_sets%end()) key => field_set_iter%key() @@ -645,7 +652,6 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) end if - allocate(IntState%Regrid(nlist), _STAT) allocate( Vvarn(nlist), _STAT) allocate(INTSTATE%STAMPOFFSET(nlist), _STAT) @@ -654,17 +660,14 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) ! ---------------------------------------------------------------------------- if( MAPL_AM_I_ROOT(vm) ) then - call ESMF_ConfigGetAttribute(config, value=HIST_CF, & label="HIST_CF:", default="HIST.rc", _RC ) unitr = GETFILE(HIST_CF, FORM='formatted', _RC) - ! for each collection do n = 1, nlist rewind(unitr) string = trim( list(n)%collection ) // '.' unitw = GETFILE(trim(string)//'rcx', FORM='formatted', _RC) - match = .false. contLine = .false. con3 = .false. @@ -698,7 +701,12 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) end if - +! Overwrite the above process if HISTORY.rc encounters DEFINE_OBS_PLATFORM for OSSE +! ---------------------------------------------------------------------------- + if( MAPL_AM_I_ROOT(vm) ) then + call regen_rcx_for_obs_platform (config, nlist, list, _RC) + end if + call ESMF_VMbarrier(vm, _RC) ! Initialize History Lists @@ -878,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 @@ -886,8 +896,6 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) list(n)%timeseries_output = .true. endif endif - call ESMF_ConfigGetAttribute(cfg, value=list(n)%recycle_track, default=.false., & - label=trim(string) // 'recycle_track:', _RC) ! Handle "backwards" mode: this is hidden (i.e. not documented) feature ! Defaults to .false. @@ -1379,7 +1387,6 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) end if enddo enddo - ! Get Output Export States ! ------------------------ @@ -2364,6 +2371,16 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) else list(n)%vdata = VerticalData(positive=list(n)%positive,_RC) end if + if (trim(list(n)%output_grid_label)=='SwathGrid') then + call list(n)%xsampler%set_param(deflation=list(n)%deflate,_RC) + call list(n)%xsampler%set_param(quantize_algorithm=list(n)%quantize_algorithm,_RC) + call list(n)%xsampler%set_param(quantize_level=list(n)%quantize_level,_RC) + call list(n)%xsampler%set_param(chunking=list(n)%chunkSize,_RC) + call list(n)%xsampler%set_param(nbits_to_keep=list(n)%nbits_to_keep,_RC) + call list(n)%xsampler%set_param(regrid_method=list(n)%regrid_method,_RC) + call list(n)%xsampler%set_param(itemOrder=intState%fileOrderAlphabetical,_RC) + endif + call list(n)%mGriddedIO%set_param(deflation=list(n)%deflate,_RC) call list(n)%mGriddedIO%set_param(quantize_algorithm=list(n)%quantize_algorithm,_RC) call list(n)%mGriddedIO%set_param(quantize_level=list(n)%quantize_level,_RC) @@ -2371,30 +2388,36 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) call list(n)%mGriddedIO%set_param(nbits_to_keep=list(n)%nbits_to_keep,_RC) call list(n)%mGriddedIO%set_param(regrid_method=list(n)%regrid_method,_RC) call list(n)%mGriddedIO%set_param(itemOrder=intState%fileOrderAlphabetical,_RC) + if (list(n)%monthly) then - nextMonth = currTime - oneMonth - dur = nextMonth - currTime - call ESMF_TimeIntervalGet(dur, s=sec, _RC) + nextMonth = currTime - oneMonth + dur = nextMonth - currTime + call ESMF_TimeIntervalGet(dur, s=sec, _RC) list(n)%timeInfo = TimeData(clock,tm,sec,IntState%stampoffset(n),funits='days') else list(n)%timeInfo = TimeData(clock,tm,MAPL_nsecf(list(n)%frequency),IntState%stampoffset(n),integer_time=intstate%integer_time) end if if (list(n)%timeseries_output) then list(n)%trajectory = HistoryTrajectory(cfg,string,clock,_RC) - call list(n)%trajectory%initialize(list(n)%items,list(n)%bundle,list(n)%timeInfo,vdata=list(n)%vdata,recycle_track=list(n)%recycle_track,_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) - if (trim(list(n)%output_grid_label)/='') then + if (trim(list(n)%output_grid_label)=='SwathGrid') then pgrid => IntState%output_grids%at(trim(list(n)%output_grid_label)) - call list(n)%mGriddedIO%CreateFileMetaData(list(n)%items,list(n)%bundle,list(n)%timeInfo,ogrid=pgrid,vdata=list(n)%vdata,global_attributes=global_attributes,_RC) + call list(n)%xsampler%Create_bundle_RH(list(n)%items,list(n)%bundle,ogrid=pgrid,vdata=list(n)%vdata,global_attributes=global_attributes,_RC) else - call list(n)%mGriddedIO%CreateFileMetaData(list(n)%items,list(n)%bundle,list(n)%timeInfo,vdata=list(n)%vdata,global_attributes=global_attributes,_RC) - end if - collection_id = o_Clients%add_hist_collection(list(n)%mGriddedIO%metadata, mode = create_mode) - call list(n)%mGriddedIO%set_param(write_collection_id=collection_id) + if (trim(list(n)%output_grid_label)/='') then + pgrid => IntState%output_grids%at(trim(list(n)%output_grid_label)) + call list(n)%mGriddedIO%CreateFileMetaData(list(n)%items,list(n)%bundle,list(n)%timeInfo,ogrid=pgrid,vdata=list(n)%vdata,global_attributes=global_attributes,_RC) + else + call list(n)%mGriddedIO%CreateFileMetaData(list(n)%items,list(n)%bundle,list(n)%timeInfo,vdata=list(n)%vdata,global_attributes=global_attributes,_RC) + end if + collection_id = o_Clients%add_hist_collection(list(n)%mGriddedIO%metadata, mode = create_mode) + call list(n)%mGriddedIO%set_param(write_collection_id=collection_id) + endif end if end if call ESMF_ConfigDestroy(cfg, _RC) @@ -3222,7 +3245,15 @@ subroutine Run ( gc, import, export, clock, rc ) type(ESMF_Time) :: lastMonth type(ESMF_TimeInterval) :: dur, oneMonth integer :: sec + type (StringGridMap) :: pt_output_grids + character(len=ESMF_MAXSTR) :: key_grid_label + type (ESMF_Grid), pointer :: pgrid + integer :: collection_id + integer :: create_mode + type(StringStringMap) :: global_attributes + type(timeData) :: timeinfo_uninit + type(ESMF_Grid) :: new_grid ! variables for "backwards" mode logical :: fwd logical, allocatable :: Ignore(:) @@ -3230,12 +3261,13 @@ subroutine Run ( gc, import, export, clock, rc ) ! ErrLog vars integer :: status logical :: file_exists + type(GriddedIOitem) :: item + type(Logger), pointer :: lgr !============================================================================= ! Begin... - _UNUSED_DUMMY(import) _UNUSED_DUMMY(export) @@ -3339,6 +3371,8 @@ subroutine Run ( gc, import, export, clock, rc ) Writing(n) = .false. else if (list(n)%timeseries_output) then Writing(n) = ESMF_AlarmIsRinging ( list(n)%trajectory%alarm ) + else if (trim(list(n)%output_grid_label)=='SwathGrid') then + Writing(n) = ESMF_AlarmIsRinging ( Hsampler%alarm ) else Writing(n) = ESMF_AlarmIsRinging ( list(n)%his_alarm ) endif @@ -3379,8 +3413,41 @@ subroutine Run ( gc, import, export, clock, rc ) end do + if(any(Writing)) call WRITE_PARALLEL("") + + ! swath only + epoch_swath_grid_case: do n=1,nlist + if (trim(list(n)%output_grid_label)=='SwathGrid') then + call Hsampler%regrid_accumulate(list(n)%xsampler,_RC) + + if( ESMF_AlarmIsRinging ( Hsampler%alarm ) ) then + create_mode = PFIO_NOCLOBBER ! defaut no overwrite + if (intState%allow_overwrite) create_mode = PFIO_CLOBBER + + ! add time to items + ! true metadata comes here from mGriddedIO%metadata + ! the mGriddedIO below only touches metadata, collection_id etc., it is safe. + ! + if (.NOT. list(n)%xsampler%have_initalized) then + list(n)%xsampler%have_initalized = .true. + global_attributes = list(n)%global_atts%define_collection_attributes(_RC) + endif + item%itemType = ItemTypeScalar + item%xname = 'time' + call list(n)%items%push_back(item) + call Hsampler%fill_time_in_bundle ('time', list(n)%xsampler%acc_bundle, _RC) + call list(n)%mGriddedIO%destroy(_RC) + call list(n)%mGriddedIO%CreateFileMetaData(list(n)%items,list(n)%xsampler%acc_bundle,timeinfo_uninit,vdata=list(n)%vdata,global_attributes=global_attributes,_RC) + call list(n)%items%pop_back() + + collection_id = o_Clients%add_hist_collection(list(n)%mGriddedIO%metadata, mode = create_mode) + call list(n)%mGriddedIO%set_param(write_collection_id=collection_id) + endif + end if + end do epoch_swath_grid_case + ! Write Id and time ! ----------------- @@ -3410,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) @@ -3457,7 +3512,7 @@ subroutine Run ( gc, import, export, clock, rc ) end if elseif (list(n)%sampler_spec == 'station') then if (list(n)%unit.eq.0) then - if (mapl_am_i_root()) call lgr%debug('%a %a',& + call lgr%debug('%a %a',& "Station_data output to new file:",trim(filename(n))) call list(n)%station_sampler%close_file_handle(_RC) call list(n)%station_sampler%create_file_handle(filename(n),_RC) @@ -3471,7 +3526,9 @@ subroutine Run ( gc, import, export, clock, rc ) inquire (file=trim(filename(n)),exist=file_exists) _ASSERT(.not.file_exists,trim(filename(n))//" being created for History output already exists") end if - call list(n)%mGriddedIO%modifyTime(oClients=o_Clients,_RC) + if (trim(list(n)%output_grid_label)/='SwathGrid') then + call list(n)%mGriddedIO%modifyTime(oClients=o_Clients,_RC) + endif list(n)%currentFile = filename(n) list(n)%unit = -1 else @@ -3491,6 +3548,7 @@ subroutine Run ( gc, import, export, clock, rc ) enddo OPENLOOP call MAPL_TimerOff(GENSTATE,"----IO Create") + call MAPL_TimerOn(GENSTATE,"----IO Write") call MAPL_TimerOn(GENSTATE,"-----IO Post") POSTLOOP: do n=1,nlist @@ -3542,11 +3600,9 @@ 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 if( INTSTATE%LCTL(n) ) then @@ -3568,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 @@ -3579,6 +3640,7 @@ subroutine Run ( gc, import, export, clock, rc ) enddo POSTLOOP + if (any(writing)) then call o_Clients%done_collective_stage(_RC) call o_Clients%post_wait() @@ -3589,6 +3651,23 @@ subroutine Run ( gc, import, export, clock, rc ) call MAPL_TimerOn(GENSTATE,"----IO Write") call MAPL_TimerOn(GENSTATE,"-----IO Wait") + + ! destroy ogrid/RH/acc_bundle, regenerate them + ! swath only + epoch_swath_regen_grid: do n=1,nlist + if (trim(list(n)%output_grid_label)=='SwathGrid') then + if( ESMF_AlarmIsRinging ( Hsampler%alarm ) ) then + key_grid_label = list(n)%output_grid_label + call Hsampler%destroy_rh_regen_ogrid ( key_grid_label, IntState%output_grids, list(n)%xsampler, _RC ) + pgrid => IntState%output_grids%at(trim(list(n)%output_grid_label)) + call list(n)%xsampler%Create_bundle_RH(list(n)%items,list(n)%bundle,ogrid=pgrid,& + vdata=list(n)%vdata,global_attributes=global_attributes,_RC) + if( MAPL_AM_I_ROOT() ) write(6,'(//)') + endif + end if + end do epoch_swath_regen_grid + + WAITLOOP: do n=1,nlist if( Writing(n) .and. list(n)%unit < 0) then @@ -3614,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 @@ -5165,6 +5240,310 @@ function get_acc_offset(current_time,ref_time,rc) result(acc_offset) end if _RETURN(_SUCCESS) end function + + + ! __ read data to object: obs_platform + ! __ for each collection: find union fields, write to collection.rcx + ! + subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) + use MAPL_scan_pattern_in_file + use MAPL_ObsUtilMod, only : obs_platform, union_platform + ! + ! Plan: + !- read and write schema + !- extract union of field lines, print out to rc + type(ESMF_Config), intent(inout) :: config + integer, intent(in) :: nlist + type(HistoryCollection), pointer :: list(:) + integer, intent(inout), optional :: rc + + character(len=ESMF_MAXSTR) :: HIST_CF + integer :: n, unitr, unitw + logical :: match, contLine, con3 + integer :: status + + character (len=ESMF_MAXSTR) :: fname + character (len=ESMF_MAXSTR) :: marker + character (len=ESMF_MAXSTR) :: line, line2 + character (len=ESMF_MAXSTR) :: string + character (len=ESMF_MAXSTR), allocatable :: str_piece(:) + type(obs_platform), allocatable :: PLFS(:) + type(obs_platform) :: p1 + integer :: k, i, j + integer :: ios, ngeoval, count, nplf + integer :: length_mx + integer :: mxseg + integer :: nseg + integer :: nfield, nplatform + integer :: nentry_name + logical :: obs_flag + integer, allocatable :: map(:) + type(Logger), pointer :: lgr + + lgr => logging%get_logger('HISTORY.sampler') + + ! + ! -- note: work on HEAD node + ! + call ESMF_ConfigGetAttribute(config, value=HIST_CF, & + label="HIST_CF:", default="HIST.rc", _RC ) + unitr = GETFILE(HIST_CF, FORM='formatted', _RC) + + call scan_count_match_bgn (unitr, 'PLATFORM.', count, .false.) + rewind(unitr) + call lgr%debug('%a %i8','count PLATFORM.', count) + if (count==0) then + rc = 0 + return + endif + nplf = count + allocate (PLFS(nplf)) + allocate (map(nplf)) + + ! __ s1. scan get platform name + nc_index/lat/lon/time + do k=1, count + call scan_begin(unitr, 'PLATFORM.', .false.) + backspace(unitr) + read(unitr, '(a)') line + i=index(line, '.') + j=index(line, ':') + _ASSERT(i>1 .AND. j>1, 'keyword PLATFORM.X is not found') + PLFS(k)%name = line(i+1:j-1) + marker=line(1:j) + + call lgr%debug('%a %a', 'marker=', trim(marker)) + call scan_contain(unitr, marker, .true.) + call scan_contain(unitr, 'index:', .false.) + backspace(unitr) + read(unitr, '(a)') line + i=index(line, ':') + PLFS(k)%nc_index = trim(line(i+1:)) + + call scan_contain(unitr, marker, .true.) + call scan_contain(unitr, 'longitude:', .false.) + backspace(unitr) + read(unitr, '(a)') line + i=index(line, ':') + PLFS(k)%nc_lon = trim(line(i+1:)) + + call scan_contain(unitr, marker, .true.) + call scan_contain(unitr, 'latitude:', .false.) + backspace(unitr) + read(unitr, '(a)') line + i=index(line, ':') + PLFS(k)%nc_lat = trim(line(i+1:)) + + call scan_contain(unitr, marker, .true.) + call scan_contain(unitr, 'time:', .false.) + backspace(unitr) + read(unitr, '(a)') line + i=index(line, ':') + PLFS(k)%nc_time = trim(line(i+1:)) + + call scan_contain(unitr, marker, .true.) + call scan_contain(unitr, 'file_name_template:', .false.) + backspace(unitr) + read(unitr, '(a)') line + i=index(line, ':') + PLFS(k)%file_name_template = trim(line(i+1:)) + + call lgr%debug('%a %a %a %a %a', & + trim( PLFS(k)%name ), & + trim( PLFS(k)%nc_lon ), & + trim( PLFS(k)%nc_lat ), & + trim( PLFS(k)%nc_time ), & + trim( PLFS(k)%file_name_template ) ) + + end do + + + ! __ s2.1 scan fields: get ngeoval / nentry_name = nword + length_mx = ESMF_MAXSTR + mxseg = 10 + allocate (str_piece(mxseg)) + rewind(unitr) + do k=1, count + call scan_begin(unitr, 'PLATFORM.', .false.) + backspace(unitr) + read(unitr, '(a)') line + i=index(line, 'PLATFORM.') + j=index(line, ':') + marker=line(1:j) + call scan_begin(unitr, marker, .true.) + call scan_contain(unitr, 'geovals_fields:', .false.) + ios=0 + ngeoval=0 + do while (ios == 0) + read (unitr, '(A)' ) line + i=index(line, '::') + if (i==0) then + ngeoval = ngeoval + 1 + call split_string_by_space (line, length_mx, mxseg, & + nseg, str_piece, status) + else + exit + endif + enddo + PLFS(k)%ngeoval = ngeoval + PLFS(k)%nentry_name = nseg +!! call lgr%debug('%a %i','ngeoval=', ngeoval) + + allocate ( PLFS(k)%field_name (nseg, ngeoval) ) + nentry_name = nseg ! assume the same for each field_name + end do + + + ! __ s2.2 scan fields: get splitted PLFS(k)%field_name + rewind(unitr) + do k=1, count + call scan_begin(unitr, 'PLATFORM.', .false.) + backspace(unitr) + read(unitr, '(a)') line + i=index(line, 'PLATFORM.') + j=index(line, ':') + marker=line(1:j) + ! + call scan_begin(unitr, marker, .true.) + call scan_contain(unitr, 'geovals_fields:', .false.) + ios=0 + ngeoval=0 + do while (ios == 0) + read (unitr, '(A)' ) line + i=index(line, '::') + if (i==0) then + ngeoval = ngeoval + 1 + call split_string_by_space (line, length_mx, mxseg, & + nseg, str_piece, status) + PLFS(k)%field_name (1:nseg, ngeoval) = str_piece(1:nseg) + else + exit + endif + enddo + end do + deallocate(str_piece) + rewind(unitr) + + !!do k=1, nplf + !! do i=1, ngeoval + !! write(6,*) 'PLFS(k)%field_name (1:nseg, ngeoval)=', PLFS(k)%field_name (1:nseg,i) + !! enddo + !!enddo + !!write(6,*) 'nlist=', nlist + + + ! __ s3: Add more entry: 'obs_files:' and 'fields:' to rcx + ! for each collection + obs_flag=.false. + do n = 1, nlist + rewind(unitr) + string = trim( list(n)%collection ) // '.' + unitw = GETFILE(trim(string)//'rcx', FORM='formatted', _RC) + match = .false. + contLine = .false. + obs_flag = .false. + do while (.true.) + read(unitr, '(A)', end=1236) line + j = index( adjustl(line), trim(adjustl(string)) ) + match = (j == 1) + if (match) then + j = index(line, trim(string)//'fields:') + contLine = (j > 0) + end if + if (match .or. contLine) then + write(unitw,'(A)') trim(line) + end if + if (contLine) then + if (adjustl(line) == '::') contLine = .false. + end if + if ( index(line, trim(string)//'ObsPlatforms:') > 0 ) then + obs_flag =.true. + line2 = line + endif + end do +1236 continue + + if (obs_flag) then + + ! __ write common nc_index,time,lon,lat + k=1 ! plat form # 1 + write(unitw, '(2(2x,a))') trim(string)//'nc_Index: ', trim(adjustl(PLFS(k)%nc_index)) + write(unitw, '(2(2x,a))') trim(string)//'nc_Time: ', trim(adjustl(PLFS(k)%nc_time)) + write(unitw, '(2(2x,a))') trim(string)//'nc_Longitude:', trim(adjustl(PLFS(k)%nc_lon)) + write(unitw, '(2(2x,a))') trim(string)//'nc_Latitude: ', trim(adjustl(PLFS(k)%nc_lat)) + write(unitw, '(/)') + + length_mx = ESMF_MAXSTR + mxseg = 100 + allocate (str_piece(mxseg)) + i = index(line2, ':') + line = adjustl ( line2(i+1:) ) + call split_string_by_space (line, length_mx, mxseg, & + nplatform, str_piece, status) +! write(6,*) 'nplatform=', nplatform +! write(6,*) 'str_piece=', str_piece(1:nplatform) +! do j=1, nplf +! write(6,*) 'PLFS(j)%name=', trim( PLFS(j)%name ) +! enddo + + ! + ! a) union the platform + ! + ! + ! find the index for each str_piece + map(:) = -1 + do i=1, nplatform ! loc collection + do j=1, nplf ! tot + if ( trim(str_piece(i)) == trim( PLFS(j)%name ) ) then + map(i)=j + end if + end do + end do + deallocate(str_piece) + + !!write(6,*) 'map(:)=', map(:) + do i=1, nplatform + k=map(i) + if (i==1) then + p1 = PLFS(k) + else + p1 = union_platform(p1, PLFS(k), _RC) + end if + end do + + nfield = p1%ngeoval + nentry_name = p1%nentry_name + do j=1, nfield + line='' + do i=1, nentry_name + line = trim(line)//' '//trim(p1%field_name(i,j)) + enddo + if (j==1) then + write(unitw, '(10(2x,a))') trim(string)//'fields:', trim(line) + else + write(unitw, '(12x,a)') trim(line) + end if + end do + write(unitw,'(a,/)') '::' + write(unitw,'(a)') 'geovals.obs_files: # table start from next line' + + do k=1, nplatform + write(unitw, '(a)') trim(adjustl(PLFS(k)%file_name_template)) + do j=1, PLFS(k)%ngeoval + line='' + do i=1, nentry_name + line = trim(line)//' '//trim(adjustl(PLFS(k)%field_name(i,j))) + enddo + write(unitw, '(a)') trim(adjustl(line)) + enddo + write(unitw, '(20a)') (('-'), j=1,20) + enddo + write(unitw,'(a)') '::' + end if + call free_file(unitw, _RC) + end do + call free_file(unitr, _RC) + end subroutine regen_rcx_for_obs_platform + end module MAPL_HistoryGridCompMod diff --git a/gridcomps/History/MAPL_HistoryTrajectoryMod.F90 b/gridcomps/History/MAPL_HistoryTrajectoryMod.F90 index aa5755f6aaf7..34aa412ef6c1 100644 --- a/gridcomps/History/MAPL_HistoryTrajectoryMod.F90 +++ b/gridcomps/History/MAPL_HistoryTrajectoryMod.F90 @@ -6,33 +6,19 @@ module HistoryTrajectoryMod use MAPL_VerticalDataMod use LocStreamFactoryMod use MAPL_LocstreamRegridderMod - use, intrinsic :: iso_fortran_env, only: REAL32, REAL64 + use MAPL_ObsUtilMod + use, intrinsic :: iso_fortran_env, only: REAL64 implicit none - private - public :: obs_unit - type :: obs_unit - integer :: nobs_epoch - type(FileMetadata), allocatable :: metadata - type(NetCDF4_FileFormatter), allocatable :: file_handle - character(len=ESMF_MAXSTR) :: name - character(len=ESMF_MAXSTR) :: obsFile_output - character(len=ESMF_MAXSTR) :: input_template - real(kind=REAL64), allocatable :: lons(:) - real(kind=REAL64), allocatable :: lats(:) - real(kind=REAL64), allocatable :: times_R8(:) - real(kind=REAL32), allocatable :: p2d(:) - real(kind=REAL32), allocatable :: p3d(:,:) - end type obs_unit + private public :: HistoryTrajectory - type :: HistoryTrajectory private type(ESMF_LocStream) :: LS_rt type(ESMF_LocStream) :: LS_ds type(LocStreamFactory) :: locstream_factory - type(obs_unit), allocatable :: obs(:) + type(obs_unit), allocatable :: obs(:) type(ESMF_Time), allocatable :: times(:) real(kind=REAL64), allocatable :: lons(:) real(kind=REAL64), allocatable :: lats(:) @@ -50,8 +36,7 @@ module HistoryTrajectoryMod logical :: do_vertical_regrid type(LocstreamRegridder) :: regridder - type(TimeData) :: time_info - logical :: recycle_track + type(TimeData) :: time_info type(ESMF_Clock) :: clock type(ESMF_Alarm), public :: alarm type(ESMF_Time) :: RingTime @@ -79,30 +64,21 @@ module HistoryTrajectoryMod logical :: is_valid contains procedure :: initialize - procedure :: reinitialize procedure :: create_variable => create_metadata_variable procedure :: create_file_handle procedure :: close_file_handle procedure :: append_file procedure :: create_new_bundle - procedure :: reset_times_to_current_day - procedure :: time_real_to_ESMF procedure :: create_grid procedure :: regrid_accumulate => regrid_accumulate_on_xsubset procedure :: destroy_rh_regen_LS procedure :: get_x_subset - procedure :: get_obsfile_Tbracket_from_epoch - procedure :: get_filename_from_template_use_index end type HistoryTrajectory interface HistoryTrajectory module procedure HistoryTrajectory_from_config end interface HistoryTrajectory - interface sort_multi_arrays_by_time - module procedure sort_three_arrays_by_time - module procedure sort_four_arrays_by_time - end interface sort_multi_arrays_by_time interface module function HistoryTrajectory_from_config(config,string,clock,rc) result(traj) @@ -113,21 +89,16 @@ module function HistoryTrajectory_from_config(config,string,clock,rc) result(tra integer, optional, intent(out) :: rc end function HistoryTrajectory_from_config - module subroutine initialize(this,items,bundle,timeInfo,vdata,recycle_track,rc) + module subroutine initialize(this,items,bundle,timeInfo,vdata,reinitialize,rc) class(HistoryTrajectory), intent(inout) :: this - type(GriddedIOitemVector), target, intent(inout) :: items - type(ESMF_FieldBundle), intent(inout) :: bundle - type(TimeData), intent(inout) :: timeInfo + type(GriddedIOitemVector), optional, intent(inout) :: items + type(ESMF_FieldBundle), optional, intent(inout) :: bundle + type(TimeData), optional, intent(inout) :: timeInfo type(VerticalData), optional, intent(inout) :: vdata - logical, optional, intent(inout) :: recycle_track + logical, optional, intent(in) :: reinitialize integer, optional, intent(out) :: rc end subroutine initialize - module subroutine reinitialize(this,rc) - class(HistoryTrajectory), intent(inout) :: this - integer, optional, intent(out) :: rc - end subroutine reinitialize - module subroutine create_metadata_variable(this,vname,rc) class(HistoryTrajectory), intent(inout) :: this character(len=*), intent(in) :: vname @@ -157,27 +128,6 @@ module subroutine append_file(this,current_time,rc) integer, optional, intent(out) :: rc end subroutine append_file - module subroutine reset_times_to_current_day(this,rc) - class(HistoryTrajectory), intent(Inout) :: this - integer, optional, intent(out) :: rc - end subroutine reset_times_to_current_day - - module subroutine sort_three_arrays_by_time(U,V,T,rc) - real(ESMF_KIND_R8) :: U(:), V(:), T(:) - integer, optional, intent(out) :: rc - end subroutine sort_three_arrays_by_time - - module subroutine sort_four_arrays_by_time(U,V,T,ID,rc) - real(ESMF_KIND_R8) :: U(:), V(:), T(:) - integer :: ID(:) - integer, optional, intent(out) :: rc - end subroutine sort_four_arrays_by_time - - module subroutine time_real_to_ESMF (this,rc) - class(HistoryTrajectory), intent(inout) :: this - integer, optional, intent(out) :: rc - end subroutine time_real_to_ESMF - module subroutine create_grid(this, rc) class(HistoryTrajectory), intent(inout) :: this integer, optional, intent(out) :: rc @@ -201,26 +151,5 @@ module subroutine destroy_rh_regen_LS (this, rc) integer, optional, intent(out) :: rc end subroutine destroy_rh_regen_LS - module subroutine get_obsfile_Tbracket_from_epoch(this, currTime, rc) - class(HistoryTrajectory), intent(inout) :: this - type(ESMF_Time), intent(in) :: currTime - integer, optional, intent(out) :: rc - end subroutine get_obsfile_Tbracket_from_epoch - - module function get_filename_from_template (time, file_template, rc) result(filename) - type(ESMF_Time), intent(in) :: time - character(len=*), intent(in) :: file_template - character(len=ESMF_MAXSTR) :: filename - integer, optional, intent(out) :: rc - end function get_filename_from_template - - module function get_filename_from_template_use_index (this, f_index, file_template, rc) result(filename) - class(HistoryTrajectory), intent(inout) :: this - character(len=*), intent(in) :: file_template - character(len=ESMF_MAXSTR) :: filename - integer, intent(in) :: f_index - integer, optional, intent(out) :: rc - end function get_filename_from_template_use_index - end interface end module HistoryTrajectoryMod diff --git a/gridcomps/History/MAPL_HistoryTrajectoryMod_smod.F90 b/gridcomps/History/MAPL_HistoryTrajectoryMod_smod.F90 index 012c3ba6b48d..1f13c1b6a3d0 100644 --- a/gridcomps/History/MAPL_HistoryTrajectoryMod_smod.F90 +++ b/gridcomps/History/MAPL_HistoryTrajectoryMod_smod.F90 @@ -19,7 +19,7 @@ use MAPL_NetCDF use MAPL_StringTemplate use Plain_netCDF_Time - use MAPL_ISO8601_DateTime_ESMF + use MAPL_ObsUtilMod use, intrinsic :: iso_fortran_env, only: REAL32 use, intrinsic :: iso_fortran_env, only: REAL64 implicit none @@ -27,17 +27,23 @@ contains module procedure HistoryTrajectory_from_config + use BinIOMod use pflogger, only : Logger, logging type(ESMF_Time) :: currTime type(ESMF_TimeInterval) :: epoch_frequency type(ESMF_TimeInterval) :: obs_time_span integer :: time_integer, second integer :: status - character(len=ESMF_MAXSTR) :: STR1 + character(len=ESMF_MAXSTR) :: STR1, line character(len=ESMF_MAXSTR) :: symd, shms - integer :: nline, ncol + integer :: nline, col + integer, allocatable :: ncol(:) + character(len=ESMF_MAXSTR), allocatable :: word(:) + integer :: nobs, head, jvar logical :: tend - integer :: i, j, k + integer :: i, j, k, M + integer :: count + integer :: unitr, unitw type(Logger), pointer :: lgr traj%clock=clock @@ -60,30 +66,7 @@ label=trim(string) // 'nc_Longitude:', _RC) call ESMF_ConfigGetAttribute(config, value=traj%nc_latitude, default="", & label=trim(string) // 'nc_Latitude:', _RC) - call ESMF_ConfigGetDim(config, nline, ncol, label=trim(string)//'obs_files:', rc=rc) - _ASSERT(rc==0 .AND. nline > 0, 'obs_files not found') - traj%nobs_type = nline - allocate (traj%obs(nline)) - do k=1, nline - allocate (traj%obs(k)%metadata) - if (mapl_am_i_root()) then - allocate (traj%obs(k)%file_handle) - end if - end do - call ESMF_ConfigFindLabel( config, trim(string)//'obs_files:', rc=rc) - lgr => logging%get_logger('HISTORY.sampler') - call lgr%debug('%a %i8', 'nobs_type=', nline) - do i=1, nline - call ESMF_ConfigNextLine( config, tableEnd=tend, rc=rc) - call ESMF_ConfigGetAttribute( config, traj%obs(i)%input_template, rc=rc) - call lgr%debug('%a %i4 %a %a', 'obs(', i, ') input_template =', & - trim(traj%obs(i)%input_template)) - j=index(traj%obs(i)%input_template , '%') - k=index(traj%obs(i)%input_template , '/', back=.true.) - _ASSERT(j>0, '% is not found, template is wrong') - traj%obs(i)%name = traj%obs(i)%input_template(k+1:j-1) - enddo call ESMF_ConfigGetAttribute(config, value=STR1, default="", & label=trim(string) // 'obs_file_begin:', _RC) @@ -133,100 +116,130 @@ call convert_twostring_2_esmfinterval (symd, shms, traj%obsfile_interval, _RC) traj%is_valid = .true. - _RETURN(_SUCCESS) - -105 format (1x,a,2x,a) -106 format (1x,a,2x,i8) - end procedure + ! __ s1. overall print + call ESMF_ConfigGetDim(config, nline, col, label=trim(string)//'obs_files:', rc=rc) + _ASSERT(rc==0 .AND. nline > 0, 'obs_files not found') + !! write(6,*) 'nline, col', nline, col + allocate(ncol(1:nline)) + + call ESMF_ConfigFindLabel( config, trim(string)//'obs_files:', _RC ) + do i = 1, nline + call ESMF_ConfigNextLine(config, _RC) + ncol(i) = ESMF_ConfigGetLen(config, _RC) + !!write(6,*) 'line', i, 'ncol(i)', ncol(i) + enddo - module procedure initialize - integer :: status - type(ESMF_Grid) :: grid - type(variable) :: v - type(GriddedIOitemVectorIterator) :: iter - type(GriddedIOitem), pointer :: item - type(ESMF_Time) :: currTime - integer :: k + + ! __ s2. find nobs && distinguish design with vs wo '------' + nobs=0 + call ESMF_ConfigFindLabel( config, trim(string)//'obs_files:', _RC) + do i=1, nline + call ESMF_ConfigNextLine( config, tableEnd=tend, _RC) + call ESMF_ConfigGetAttribute( config, STR1, _RC) + if ( index(trim(STR1), '-----') > 0 ) nobs=nobs+1 + enddo - this%bundle=bundle - this%items=items - if (present(vdata)) then - this%vdata=vdata + ! __ s3. retrieve template and geoval, set metadata file_handle + lgr => logging%get_logger('HISTORY.sampler') + if ( nobs == 0 ) then + ! + ! treatment-1: + ! + traj%nobs_type = nline ! here .rc format cannot have empty spaces + allocate (traj%obs(nline)) + call ESMF_ConfigFindLabel( config, trim(string)//'obs_files:', _RC) + do i=1, nline + call ESMF_ConfigNextLine( config, tableEnd=tend, _RC) + call ESMF_ConfigGetAttribute( config, traj%obs(i)%input_template, _RC) + traj%obs(i)%export_all_geoval = .true. + enddo else - this%vdata=VerticalData(_RC) - end if - do k=1, this%nobs_type - call this%vdata%append_vertical_metadata(this%obs(k)%metadata,this%bundle,_RC) - end do - this%do_vertical_regrid = (this%vdata%regrid_type /= VERTICAL_METHOD_NONE) - if (this%vdata%regrid_type == VERTICAL_METHOD_ETA2LEV) call this%vdata%get_interpolating_variable(this%bundle,_RC) + ! + !-- selectively output geovals + ! treatment-2: + ! + traj%nobs_type = nobs + allocate (traj%obs(nobs)) + ! + nobs=0 ! reuse counter + head=1 + jvar=0 - call ESMF_ClockGet (this%clock, CurrTime=currTime, _RC) - call this%get_obsfile_Tbracket_from_epoch(currTime, _RC) - if (this%obsfile_Te_index < 0) then - if (mapl_am_I_root()) then - write(6,*) "model start time is earlier than obsfile_start_time" - write(6,*) "solution: adjust obsfile_start_time and Epoch in rc file" - end if - _FAIL("obs file not found at init time") - endif - call this%create_grid(_RC) - call ESMF_FieldBundleGet(this%bundle,grid=grid,_RC) - this%regridder = LocStreamRegridder(grid,this%LS_ds,_RC) - this%output_bundle = this%create_new_bundle(_RC) - this%acc_bundle = this%create_new_bundle(_RC) - this%time_info = timeInfo + ! + ! count '------' in history.rc as special markers for ngeoval + ! + call ESMF_ConfigFindLabel(config, trim(string)//'obs_files:', _RC) + do i=1, nline + call ESMF_ConfigNextLine(config, tableEnd=tend, _RC) + M = ncol(i) + _ASSERT(M>=1, '# of columns should be >= 1') + allocate (word(M)) + count=0 + do col=1, M + call ESMF_ConfigGetAttribute(config, word(col), _RC) + if (trim(word(col))/=',') then + count=count+1 + end if + enddo + if (count ==1 .or. count==2) then + ! 1-item case: file template or one-var + ! 2-item : var1 , 'root' case + STR1=trim(word(1)) + else + ! 3-item : var1 , 'root', var1_alias case + STR1=trim(word(M)) + end if + deallocate(word) + if ( index(trim(STR1), '-----') == 0 ) then + if (head==1 .AND. trim(STR1)/='') then + nobs=nobs+1 + traj%obs(nobs)%input_template = trim(STR1) + traj%obs(nobs)%export_all_geoval = .false. + head=0 + else + if (trim(STR1)/='') then + jvar=jvar+1 + traj%obs(nobs)%geoval_name(jvar) = trim(STR1) + end if + end if + else + traj%obs(nobs)%ngeoval=jvar + head=1 + jvar=0 + endif + enddo + end if - do k=1, this%nobs_type - call this%obs(k)%metadata%add_dimension(this%nc_index, this%obs(k)%nobs_epoch) - if (this%time_info%integer_time) then - v = Variable(type=PFIO_INT32,dimensions=this%nc_index) - else - v = Variable(type=PFIO_REAL32,dimensions=this%nc_index) + do k=1, traj%nobs_type + allocate (traj%obs(k)%metadata) + if (mapl_am_i_root()) then + allocate (traj%obs(k)%file_handle) end if - call v%add_attribute('units', this%datetime_units) - call v%add_attribute('long_name', 'dateTime') - call this%obs(k)%metadata%add_variable(this%var_name_time,v) - - v = variable(type=PFIO_REAL64,dimensions=this%nc_index) - call v%add_attribute('units','degrees_east') - call v%add_attribute('long_name','longitude') - call this%obs(k)%metadata%add_variable(this%var_name_lon,v) - - v = variable(type=PFIO_REAL64,dimensions=this%nc_index) - call v%add_attribute('units','degrees_north') - call v%add_attribute('long_name','latitude') - call this%obs(k)%metadata%add_variable(this%var_name_lat,v) end do - - iter = this%items%begin() - do while (iter /= this%items%end()) - item => iter%get() - if (item%itemType == ItemTypeScalar) then - call this%create_variable(item%xname,_RC) - else if (item%itemType == ItemTypeVector) then - call this%create_variable(item%xname,_RC) - call this%create_variable(item%yname,_RC) - end if - call iter%next() - enddo - - this%recycle_track=.false. - if (present(recycle_track)) then - this%recycle_track=recycle_track - end if - if (this%recycle_track) then - call this%reset_times_to_current_day(_RC) - end if - + + call lgr%debug('%a %i8', 'nobs_type=', traj%nobs_type) + do i=1, traj%nobs_type + call lgr%debug('%a %i4 %a %a', 'obs(', i, ') input_template =', & + trim(traj%obs(i)%input_template)) + j=index(traj%obs(i)%input_template , '%') + k=index(traj%obs(i)%input_template , '/', back=.true.) + _ASSERT(j>0, '% is not found, template is wrong') + traj%obs(i)%name = traj%obs(i)%input_template(k+1:j-1) + end do + _RETURN(_SUCCESS) - end procedure initialize - +105 format (1x,a,2x,a) +106 format (1x,a,2x,i8) + end procedure HistoryTrajectory_from_config - module procedure reinitialize + + ! + !-- integrate both initialize and reinitialize + ! + module procedure initialize integer :: status type(ESMF_Grid) :: grid type(variable) :: v @@ -235,13 +248,26 @@ type(ESMF_Time) :: currTime integer :: k - do k=1, this%nobs_type - allocate (this%obs(k)%metadata) - if (mapl_am_i_root()) then - allocate (this%obs(k)%file_handle) + if (.not. present(reinitialize)) then + if(present(bundle)) this%bundle=bundle + if(present(items)) this%items=items + if(present(timeInfo)) this%time_info=timeInfo + if (present(vdata)) then + this%vdata=vdata + else + this%vdata=VerticalData(_RC) end if - end do - + else + if (reinitialize) then + do k=1, this%nobs_type + allocate (this%obs(k)%metadata) + if (mapl_am_i_root()) then + allocate (this%obs(k)%file_handle) + end if + end do + end if + end if + do k=1, this%nobs_type call this%vdata%append_vertical_metadata(this%obs(k)%metadata,this%bundle,_RC) end do @@ -249,7 +275,10 @@ if (this%vdata%regrid_type == VERTICAL_METHOD_ETA2LEV) call this%vdata%get_interpolating_variable(this%bundle,_RC) call ESMF_ClockGet (this%clock, CurrTime=currTime, _RC) - call this%get_obsfile_Tbracket_from_epoch(currTime, _RC) + call get_obsfile_Tbracket_from_epoch(currTime, & + this%obsfile_start_time, this%obsfile_end_time, & + this%obsfile_interval, this%epoch_frequency, & + this%obsfile_Ts_index, this%obsfile_Te_index, _RC) if (this%obsfile_Te_index < 0) then if (mapl_am_I_root()) then write(6,*) "model start time is earlier than obsfile_start_time" @@ -264,6 +293,7 @@ this%output_bundle = this%create_new_bundle(_RC) this%acc_bundle = this%create_new_bundle(_RC) + do k=1, this%nobs_type call this%obs(k)%metadata%add_dimension(this%nc_index, this%obs(k)%nobs_epoch) if (this%time_info%integer_time) then @@ -286,6 +316,7 @@ call this%obs(k)%metadata%add_variable(this%var_name_lat,v) end do + ! push varible names down to each obs(k); see create_metadata_variable iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() @@ -297,9 +328,11 @@ end if call iter%next() enddo + _RETURN(_SUCCESS) - end procedure reinitialize + end procedure initialize + module procedure create_metadata_variable @@ -308,7 +341,7 @@ logical :: is_present integer :: field_rank, status character(len=ESMF_MAXSTR) :: var_name,long_name,units,vdims - integer :: k + integer :: k, ig call ESMF_FieldBundleGet(this%bundle,vname,field=field,_RC) call ESMF_FieldGet(field,name=var_name,rank=field_rank,_RC) @@ -337,7 +370,11 @@ call v%add_attribute('valid_range',(/-MAPL_UNDEF,MAPL_UNDEF/)) do k = 1, this%nobs_type - call this%obs(k)%metadata%add_variable(trim(var_name),v,_RC) + do ig = 1, this%obs(k)%ngeoval + if (trim(var_name) == trim(this%obs(k)%geoval_name(ig))) then + call this%obs(k)%metadata%add_variable(trim(var_name),v,_RC) + endif + enddo enddo _RETURN(_SUCCESS) @@ -483,10 +520,16 @@ i=index(this%nc_longitude, '/') _ASSERT (i>0, 'group name not found') grp_name = this%nc_longitude(1:i-1) - this%var_name_lat = this%nc_latitude(i+1:) this%var_name_lon = this%nc_longitude(i+1:) + i=index(this%nc_latitude, '/') + this%var_name_lat = this%nc_latitude(i+1:) + i=index(this%nc_time, '/') this%var_name_time= this%nc_time(i+1:) + call lgr%debug('%a', 'grp_name,this%var_name_lat,this%var_name_lon,this%var_name_time') + call lgr%debug('%a %a %a %a', & + trim(grp_name),trim(this%var_name_lat),trim(this%var_name_lon),trim(this%var_name_time)) + L=0 fid_s=this%obsfile_Ts_index fid_e=this%obsfile_Te_index @@ -494,9 +537,9 @@ allocate(this%lons(0),this%lats(0),_STAT) allocate(this%times_R8(0),_STAT) allocate(this%obstype_id(0),_STAT) - this%epoch_index(1:2)=0 + this%epoch_index(1:2) = 0 this%nobs_epoch = 0 - rc=0 + rc = 0 return end if @@ -505,10 +548,14 @@ do k=1, this%nobs_type j = max (fid_s, L) do while (j<=fid_e) - filename = this%get_filename_from_template_use_index(j, this%obs(k)%input_template, _RC) - call lgr%debug('%a %a', 'input filename: ', trim(filename)) - call get_ncfile_dimension(filename, tdim=num_times, key_time=this%nc_index, _RC) - len = len + num_times + filename = get_filename_from_template_use_index( & + this%obsfile_start_time, this%obsfile_interval, & + j, this%obs(k)%input_template, _RC) + if (filename /= '') then + call lgr%debug('%a %a', 'true filename: ', trim(filename)) + call get_ncfile_dimension(filename, tdim=num_times, key_time=this%nc_index, _RC) + len = len + num_times + end if j=j+1 enddo enddo @@ -522,17 +569,22 @@ do k=1, this%nobs_type j = max (fid_s, L) do while (j<=fid_e) - filename = this%get_filename_from_template_use_index(j, this%obs(k)%input_template, _RC) - call get_ncfile_dimension(trim(filename), tdim=num_times, key_time=this%nc_index, _RC) - call get_v1d_netcdf_R8 (filename, this%var_name_lon, lons_full(len+1:), num_times, group_name=grp_name) - call get_v1d_netcdf_R8 (filename, this%var_name_lat, lats_full(len+1:), num_times, group_name=grp_name) - call get_v1d_netcdf_R8 (filename, this%var_name_time, times_R8_full(len+1:), num_times, group_name=grp_name) - call get_attribute_from_group (filename, grp_name, this%var_name_time, "units", timeunits_file) - obstype_id_full(len+1:len+num_times) = k - call lgr%debug('%a %f25.12, %f25.12', 'times_R8_full(1:200:100)', & - times_R8_full(1), times_R8_full(200)) - - len = len + num_times + filename = get_filename_from_template_use_index( & + this%obsfile_start_time, this%obsfile_interval, & + j, this%obs(k)%input_template, _RC) + if (filename /= '') then + call get_ncfile_dimension(trim(filename), tdim=num_times, key_time=this%nc_index, _RC) + call get_v1d_netcdf_R8 (filename, this%var_name_lon, lons_full(len+1:), num_times, group_name=grp_name) + call get_v1d_netcdf_R8 (filename, this%var_name_lat, lats_full(len+1:), num_times, group_name=grp_name) + call get_v1d_netcdf_R8 (filename, this%var_name_time, times_R8_full(len+1:), num_times, group_name=grp_name) + + call get_attribute_from_group (filename, grp_name, this%var_name_time, "units", timeunits_file) + obstype_id_full(len+1:len+num_times) = k + call lgr%debug('%a %f25.12, %f25.12', 'times_R8_full(1:200:100)', & + times_R8_full(1), times_R8_full(200)) + + len = len + num_times + end if j=j+1 enddo enddo @@ -668,6 +720,8 @@ ! defer destroy fieldB at regen_grid step ! end if + + _RETURN(_SUCCESS) end procedure create_grid @@ -689,7 +743,7 @@ integer :: lm integer :: rank integer :: status - integer :: j, k + integer :: j, k, ig integer, allocatable :: ix(:) if (.NOT. this%is_valid) then @@ -733,6 +787,8 @@ iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() + !!write(6, '(2x,a,2x,a)') 'item%xname', trim(item%xname) + if (item%itemType == ItemTypeScalar) then call ESMF_FieldBundleGet(this%acc_bundle,trim(item%xname),field=acc_field,_RC) call ESMF_FieldGet(acc_field,rank=rank,_RC) @@ -771,8 +827,12 @@ is = 1 nx = this%obs(k)%nobs_epoch if (nx>0) then - call this%obs(k)%file_handle%put_var(trim(item%xname), this%obs(k)%p2d(1:nx), & - start=[is],count=[nx]) + do ig = 1, this%obs(k)%ngeoval + if (trim(item%xname) == trim(this%obs(k)%geoval_name(ig))) then + call this%obs(k)%file_handle%put_var(trim(item%xname), this%obs(k)%p2d(1:nx), & + start=[is],count=[nx]) + end if + enddo endif enddo end if @@ -817,8 +877,12 @@ is = 1 nx = this%obs(k)%nobs_epoch if (nx>0) then - call this%obs(k)%file_handle%put_var(trim(item%xname), this%obs(k)%p3d(:,:), & - start=[is,1],count=[nx,size(p_acc_rt_3d,2)]) + do ig = 1, this%obs(k)%ngeoval + if (trim(item%xname) == trim(this%obs(k)%geoval_name(ig))) then + call this%obs(k)%file_handle%put_var(trim(item%xname), this%obs(k)%p3d(:,:), & + start=[is,1],count=[nx,size(p_acc_rt_3d,2)]) + end if + end do endif enddo !!write(6,'(10f8.2)') p_acc_rt_3d(:,:) @@ -874,7 +938,6 @@ call this%vdata%setup_eta_to_pressure(_RC) endif - iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() @@ -992,7 +1055,7 @@ this%epoch_index(1:2)=0 - call this%reinitialize(_RC) + call this%initialize(reinitialize=.true., _RC) _RETURN(ESMF_SUCCESS) @@ -1054,200 +1117,4 @@ end procedure get_x_subset - module procedure get_obsfile_Tbracket_from_epoch - implicit none - integer :: status - - type(ESMF_Time) :: T1, Tn - type(ESMF_Time) :: cT1 - type(ESMF_Time) :: Ts, Te - type(ESMF_TimeInterval) :: dT1, dT2, dTs, dTe - real(ESMF_KIND_R8) :: dT0_s, dT1_s, dT2_s - real(ESMF_KIND_R8) :: s1, s2 - integer :: n1, n2 - - T1 = this%obsfile_start_time - Tn = this%obsfile_end_time - - cT1 = currTime - dT1 = currTime - T1 - dT2 = currTime + this%epoch_frequency - T1 - - call ESMF_TimeIntervalGet(this%obsfile_interval, s_r8=dT0_s, rc=status) - call ESMF_TimeIntervalGet(dT1, s_r8=dT1_s, rc=status) - call ESMF_TimeIntervalGet(dT2, s_r8=dT2_s, rc=status) - n1 = floor (dT1_s / dT0_s) - n2 = floor (dT2_s / dT0_s) - s1 = n1 * dT0_s - s2 = n2 * dT0_s - call ESMF_TimeIntervalSet(dTs, s_r8=s1, rc=status) - call ESMF_TimeIntervalSet(dTe, s_r8=s2, rc=status) - Ts = T1 + dTs - Te = T1 + dTe - - this%obsfile_Ts_index = n1 - if ( dT2_s - n2*dT0_s < 1 ) then - this%obsfile_Te_index = n2 - 1 - else - this%obsfile_Te_index = n2 - end if - - _RETURN(ESMF_SUCCESS) - end procedure get_obsfile_Tbracket_from_epoch - - - module procedure get_filename_from_template - integer :: itime(2) - integer :: nymd, nhms - integer :: status - - stop 'DO not use get_filename_from_template' - call ESMF_time_to_two_integer(time, itime, _RC) - print*, 'two integer time, itime(:)', itime(1:2) - nymd = itime(1) - nhms = itime(2) - call fill_grads_template ( filename, file_template, & - experiment_id='', nymd=nymd, nhms=nhms, _RC ) - print*, 'ck: this%obsFile_T=', trim(filename) - _RETURN(ESMF_SUCCESS) - end procedure get_filename_from_template - - - module procedure get_filename_from_template_use_index - integer :: itime(2) - integer :: nymd, nhms - integer :: status - real(ESMF_KIND_R8) :: dT0_s - real(ESMF_KIND_R8) :: s - type(ESMF_TimeInterval) :: dT - type(ESMF_Time) :: time - - call ESMF_TimeIntervalGet(this%obsfile_interval, s_r8=dT0_s, rc=status) - s = dT0_s * f_index - call ESMF_TimeIntervalSet(dT, s_r8=s, rc=status) - time = this%obsfile_start_time + dT - - call ESMF_time_to_two_integer(time, itime, _RC) - nymd = itime(1) - nhms = itime(2) - call fill_grads_template ( filename, file_template, & - experiment_id='', nymd=nymd, nhms=nhms, _RC ) - - _RETURN(ESMF_SUCCESS) - end procedure get_filename_from_template_use_index - - - - module procedure time_real_to_ESMF - type(ESMF_TimeInterval) :: interval - type(ESMF_Time) :: time0 - type(ESMF_Time) :: time1 - character(len=:), allocatable :: tunit - character(len=ESMF_MAXSTR) :: datetime_units - integer :: i, len - integer :: int_time - integer :: status - - datetime_units = this%datetime_units - len = size (this%times_R8) - do i=1, len - int_time = this%times_R8(i) - call convert_NetCDF_DateTime_to_ESMF(int_time, datetime_units, interval, time0, time=time1, time_unit=tunit, _RC) - this%times(i) = time1 - enddo - - _RETURN(_SUCCESS) - end procedure time_real_to_ESMF - - - - - module procedure reset_times_to_current_day - - integer :: i,status,h,m,yp,mp,dp,s,ms,us,ns - type(ESMF_Clock) :: clock - type(ESMF_Time) :: current_time - integer :: year,month,day - - call this%time_info%get(clock=clock,_RC) - call ESMF_ClockGet(clock,currtime=current_time,_RC) - call ESMF_TimeGet(current_time,yy=year,mm=month,dd=day,_RC) - do i=1,size(this%times) - call ESMF_TimeGet(this%times(i),yy=yp,mm=mp,dd=dp,h=h,m=m,s=s,ms=ms,us=us,ns=ns,_RC) - call ESMF_TimeSet(this%times(i),yy=year,mm=month,dd=day,h=h,m=m,s=s,ms=ms,us=us,ns=ns,_RC) - enddo - - end procedure reset_times_to_current_day - - - module procedure sort_three_arrays_by_time - integer :: i, len - integer, allocatable :: IA(:) - integer(ESMF_KIND_I8), allocatable :: IX(:) - real(ESMF_KIND_R8), allocatable :: X(:) - - _ASSERT (size(U)==size(V), 'U,V different dimension') - _ASSERT (size(U)==size(T), 'U,T different dimension') - len = size (T) - - allocate (IA(len), IX(len), X(len)) - do i=1, len - IX(i)=T(i) - IA(i)=i - enddo - call MAPL_Sort(IX,IA) - - X = U - do i=1, len - U(i) = X(IA(i)) - enddo - X = V - do i=1, len - V(i) = X(IA(i)) - enddo - X = T - do i=1, len - T(i) = X(IA(i)) - enddo - _RETURN(_SUCCESS) - end procedure sort_three_arrays_by_time - - - module procedure sort_four_arrays_by_time - integer :: i, len - integer, allocatable :: IA(:) - integer(ESMF_KIND_I8), allocatable :: IX(:) - real(ESMF_KIND_R8), allocatable :: X(:) - integer, allocatable :: NX(:) - - _ASSERT (size(U)==size(V), 'U,V different dimension') - _ASSERT (size(U)==size(T), 'U,T different dimension') - len = size (T) - - allocate (IA(len), IX(len), X(len), NX(len)) - do i=1, len - IX(i)=T(i) - IA(i)=i - enddo - call MAPL_Sort(IX,IA) - - X = U - do i=1, len - U(i) = X(IA(i)) - enddo - X = V - do i=1, len - V(i) = X(IA(i)) - enddo - X = T - do i=1, len - T(i) = X(IA(i)) - enddo - NX = ID - do i=1, len - ID(i) = NX(IA(i)) - enddo - _RETURN(_SUCCESS) - end procedure sort_four_arrays_by_time - end submodule HistoryTrajectory_implement diff --git a/gridcomps/History/MAPL_StationSamplerMod.F90 b/gridcomps/History/MAPL_StationSamplerMod.F90 index fff4cf073a0f..7a6da62832c5 100644 --- 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 diff --git a/griddedio/GriddedIO.F90 b/griddedio/GriddedIO.F90 index 20d09063a7f8..2fa6bb522d89 100644 --- a/griddedio/GriddedIO.F90 +++ b/griddedio/GriddedIO.F90 @@ -31,7 +31,7 @@ module MAPL_GriddedIOMod private type, public :: MAPL_GriddedIO - type(FileMetaData) :: metadata + type(FileMetaData), allocatable :: metadata type(fileMetadataUtils), pointer :: current_file_metadata integer :: write_collection_id integer :: read_collection_id @@ -73,6 +73,7 @@ module MAPL_GriddedIOMod procedure :: request_data_from_file procedure :: process_data_from_file procedure :: swap_undef_value + procedure :: destroy end type MAPL_GriddedIO interface MAPL_GriddedIO @@ -114,7 +115,7 @@ subroutine CreateFileMetaData(this,items,bundle,timeInfo,vdata,ogrid,global_attr type(TimeData), intent(inout) :: timeInfo type(VerticalData), intent(inout), optional :: vdata type (ESMF_Grid), intent(inout), pointer, optional :: ogrid - type(StringStringMap), intent(in), target, optional :: global_attributes + type(StringStringMap), target, intent(in), optional :: global_attributes integer, intent(out), optional :: rc type(ESMF_Grid) :: input_grid @@ -128,6 +129,11 @@ subroutine CreateFileMetaData(this,items,bundle,timeInfo,vdata,ogrid,global_attr character(len=:), pointer :: attr_name, attr_val integer :: status + if ( allocated (this%metadata) ) deallocate(this%metadata) + allocate(this%metadata) + + call MAPL_FieldBundleDestroy(this%output_bundle, _RC) + this%items = items this%input_bundle = bundle this%output_bundle = ESMF_FieldBundleCreate(rc=status) @@ -141,9 +147,11 @@ subroutine CreateFileMetaData(this,items,bundle,timeInfo,vdata,ogrid,global_attr call ESMF_FieldBundleGet(this%input_bundle,grid=this%output_grid,rc=status) _VERIFY(status) end if + this%regrid_handle => new_regridder_manager%make_regridder(input_grid,this%output_grid,this%regrid_method,rc=status) _VERIFY(status) + ! We get the regrid_method here because in the case of Identity, we set it to ! REGRID_METHOD_IDENTITY in the regridder constructor if identity. Now we need ! to change the regrid_method in the GriddedIO object to be the same as the @@ -156,6 +164,7 @@ subroutine CreateFileMetaData(this,items,bundle,timeInfo,vdata,ogrid,global_attr _VERIFY(status) call factory%append_metadata(this%metadata) + if (present(vdata)) then this%vdata=vdata else @@ -179,6 +188,7 @@ subroutine CreateFileMetaData(this,items,bundle,timeInfo,vdata,ogrid,global_attr call this%check_chunking(this%vdata%lm,_RC) end if + order = this%metadata%get_order(rc=status) _VERIFY(status) metadataVarsSize = order%size() @@ -213,7 +223,16 @@ subroutine CreateFileMetaData(this,items,bundle,timeInfo,vdata,ogrid,global_attr end if _RETURN(_SUCCESS) - end subroutine CreateFileMetaData + end subroutine CreateFileMetaData + + + subroutine destroy(this, rc) + class (MAPL_GriddedIO), intent(inout) :: this + integer, intent(out), optional :: rc + if(allocated(this%chunking)) deallocate(this%chunking) + _RETURN(_SUCCESS) + end subroutine destroy + subroutine set_param(this,deflation,quantize_algorithm,quantize_level,chunking,nbits_to_keep,regrid_method,itemOrder,write_collection_id,rc) class (MAPL_GriddedIO), intent(inout) :: this @@ -483,6 +502,7 @@ subroutine bundlepost(this,filename,oClients,rc) end if else tindex = -1 + call this%stage2DLatLon(filename,oClients=oClients,_RC) end if if (this%vdata%regrid_type==VERTICAL_METHOD_ETA2LEV) then @@ -812,14 +832,14 @@ subroutine RegridVector(this,xName,yName,rc) end subroutine RegridVector subroutine stage2DLatLon(this, fileName, oClients, rc) - class (MAPL_GriddedIO), intent(inout) :: this + class (MAPL_GriddedIO), target, intent(inout) :: this character(len=*), intent(in) :: fileName - type (ClientManager), optional, intent(inout) :: oClients + type (ClientManager), optional, target, intent(inout) :: oClients integer, optional, intent(out) :: rc integer :: status real(REAL64), pointer :: ptr2d(:,:) - type(ArrayReference) :: ref + type(ArrayReference), target :: ref class (AbstractGridFactory), pointer :: factory integer, allocatable :: localStart(:),globalStart(:),globalCount(:) logical :: hasll @@ -840,9 +860,11 @@ subroutine stage2DLatLon(this, fileName, oClients, rc) farrayPtr=ptr2d, rc=status) _VERIFY(STATUS) this%lons=ptr2d*MAPL_RADIANS_TO_DEGREES + ref = ArrayReference(this%lons) - call oClients%collective_stage_data(this%write_collection_id,trim(filename),'lons', & - ref,start=localStart, global_start=GlobalStart, global_count=GlobalCount) + call oClients%collective_stage_data(this%write_collection_id,trim(filename),'lons', & + ref,start=localStart, global_start=GlobalStart, global_count=GlobalCount) + call ESMF_GridGetCoord(this%output_grid, localDE=0, coordDim=2, & staggerloc=ESMF_STAGGERLOC_CENTER, & farrayPtr=ptr2d, rc=status) @@ -883,6 +905,7 @@ subroutine stage2DLatLon(this, fileName, oClients, rc) call oClients%collective_stage_data(this%write_collection_id,trim(filename),'corner_lats', & ref,start=localStart, global_start=GlobalStart, global_count=GlobalCount) end if + _RETURN(_SUCCESS) end subroutine stage2DLatLon @@ -936,11 +959,12 @@ subroutine stageData(this, field, fileName, tIndex, oClients, rc) allocate(ptr2d(0,0)) end if ref = factory%generate_file_reference2D(Ptr2D) - allocate(localStart,source=[gridLocalStart,1]) if (tindex > -1) then + allocate(localStart,source=[gridLocalStart,1]) allocate(globalStart,source=[gridGlobalStart,tindex]) allocate(globalCount,source=[gridGlobalCount,1]) else + allocate(localStart,source=[gridLocalStart]) allocate(globalStart,source=gridGlobalStart) allocate(globalCount,source=gridGlobalCount) end if @@ -957,17 +981,19 @@ subroutine stageData(this, field, fileName, tIndex, oClients, rc) allocate(ptr3d(0,0,0)) end if ref = factory%generate_file_reference3D(Ptr3D) - allocate(localStart,source=[gridLocalStart,1,1]) if (tindex > -1) then + allocate(localStart,source=[gridLocalStart,1,1]) allocate(globalStart,source=[gridGlobalStart,1,tindex]) allocate(globalCount,source=[gridGlobalCount,lm,1]) else + allocate(localStart,source=[gridLocalStart,1]) allocate(globalStart,source=[gridGlobalStart,1]) allocate(globalCount,source=[gridGlobalCount,lm]) end if else _FAIL( "Rank not supported") end if + call oClients%collective_stage_data(this%write_collection_id,trim(filename),trim(fieldName), & ref,start=localStart, global_start=GlobalStart, global_count=GlobalCount) _RETURN(_SUCCESS) diff --git a/profiler/tests/test_MeterNodeIterator.pf b/profiler/tests/test_MeterNodeIterator.pf index 4e15735d2a3b..77d8253353b9 100644 --- a/profiler/tests/test_MeterNodeIterator.pf +++ b/profiler/tests/test_MeterNodeIterator.pf @@ -12,11 +12,11 @@ contains class (AbstractMeterNodeIterator), allocatable :: iter_1 class (AbstractMeterNodeIterator), allocatable :: iter_2 - node = MeterNode('all', AdvancedMeter(MpiTimerGauge())) - iter_1 = node%begin() - iter_2 = node%begin() + allocate(iter_1, source=node%begin()) + allocate(iter_2, source=node%begin()) + @assertTrue(iter_1 == iter_2) @assertFalse(iter_1 /= iter_2) @assertTrue(iter_1 /= node%end()) @@ -46,8 +46,8 @@ contains class (AbstractMeterNodeIterator), allocatable :: iter_2 node = MeterNode('all', AdvancedMeter(MpiTimerGauge())) - iter_1 = node%begin() - iter_2 = node%begin() + allocate(iter_1, source=node%begin()) + allocate(iter_2, source=node%begin()) call node%add_child('a', AdvancedMeter(MpiTimerGauge())) @@ -73,7 +73,7 @@ contains node = MeterNode('all', AdvancedMeter(MpiTimerGauge())) count = 0 - iter = node%begin() + allocate(iter, source=node%begin()) do while (iter /= node%end()) count = count + 1 call iter%next() @@ -98,7 +98,7 @@ contains call node%add_child('c', AdvancedMeter(MpiTimerGauge())) count = 0 - iter = node%begin() + allocate(iter, source=node%begin()) do while (iter /= node%end()) count = count + 1 call iter%next() @@ -162,7 +162,7 @@ contains count = 0 - iter = node%begin() + allocate(iter, source=node%begin()) do while (iter /= node%end()) count = count + 1 t => iter%get_meter()