diff --git a/.circleci/config.yml b/.circleci/config.yml index 88afe2ca124e..6d6cc4e4c934 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -16,7 +16,7 @@ parameters: # Anchors to prevent forgetting to update a version os_version: &os_version ubuntu24 -baselibs_version: &baselibs_version v7.27.0 +baselibs_version: &baselibs_version v7.29.0 bcs_version: &bcs_version v11.6.0 tag_build_arg_name: &tag_build_arg_name maplversion diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 3e4740e794dc..1dae02b37a48 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -17,6 +17,9 @@ jobs: steps: - name: Checkout uses: actions/checkout@v4 + with: + fetch-depth: 0 + filter: blob:none - name: Build and Deploy Docs uses: ./.github/actions/deploy-ford-docs @@ -34,6 +37,9 @@ jobs: steps: - name: Checkout uses: actions/checkout@v4 + with: + fetch-depth: 0 + filter: blob:none - name: Build and Deploy Dev Docs uses: ./.github/actions/deploy-ford-docs diff --git a/.github/workflows/push-to-develop.yml b/.github/workflows/push-to-develop.yml index 449820dd235a..59f89da2710d 100644 --- a/.github/workflows/push-to-develop.yml +++ b/.github/workflows/push-to-develop.yml @@ -14,6 +14,7 @@ jobs: uses: actions/checkout@v4 with: fetch-depth: 0 + filter: blob:none - name: Run the action uses: devops-infra/action-pull-request@v0.5.5 with: diff --git a/.github/workflows/push-to-main.yml b/.github/workflows/push-to-main.yml index 6c8a112a92e9..e558b7db0ce2 100644 --- a/.github/workflows/push-to-main.yml +++ b/.github/workflows/push-to-main.yml @@ -14,6 +14,7 @@ jobs: uses: actions/checkout@v4 with: fetch-depth: 0 + filter: blob:none - name: Run the action uses: devops-infra/action-pull-request@v0.5.5 with: diff --git a/.github/workflows/validate_yaml_files.yml b/.github/workflows/validate_yaml_files.yml index 449db6e674da..4b4c11672284 100644 --- a/.github/workflows/validate_yaml_files.yml +++ b/.github/workflows/validate_yaml_files.yml @@ -15,7 +15,11 @@ jobs: validate-YAML: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - name: Checkout repo + uses: actions/checkout@v4 + with: + fetch-depth: 0 + filter: blob:none - id: yaml-lint name: yaml-lint uses: ibiqlik/action-yamllint@v3 diff --git a/.github/workflows/workflow.yml b/.github/workflows/workflow.yml index e5f266150bc7..865a7d3d7136 100644 --- a/.github/workflows/workflow.yml +++ b/.github/workflows/workflow.yml @@ -19,6 +19,9 @@ jobs: steps: - name: Checkout uses: actions/checkout@v4 + with: + fetch-depth: 0 + filter: blob:none - name: Build and Deploy Docs uses: ./.github/actions/deploy-ford-docs @@ -35,7 +38,7 @@ jobs: name: Build and Test MAPL GNU runs-on: ubuntu-latest container: - image: gmao/ubuntu24-geos-env-mkl:v7.27.0-openmpi_5.0.5-gcc_14.2.0 + image: gmao/ubuntu24-geos-env-mkl:v7.29.0-openmpi_5.0.5-gcc_14.2.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 @@ -86,7 +89,7 @@ jobs: name: Build and Test MAPL Intel runs-on: ubuntu-latest container: - image: gmao/ubuntu24-geos-env:v7.27.0-intelmpi_2021.13-ifort_2021.13 + image: gmao/ubuntu24-geos-env:v7.29.0-intelmpi_2021.13-ifort_2021.13 # 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 @@ -102,6 +105,7 @@ jobs: uses: actions/checkout@v4 with: fetch-depth: 1 + filter: blob:none - name: Set all directories as git safe run: | git config --global --add safe.directory '*' diff --git a/CHANGELOG.md b/CHANGELOG.md index ad6411dbd6fa..d30421e5361c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,6 +17,38 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Deprecated +## [2.52.0] - 2025-01-17 + +### Added + +- Added subroutine to read nc4 tile file +- Added optional `start_date` and `start_time` to control the output window for each History collection. No output will be written before then. If not specified, these default to the beginning of the experiment. +- Added utility to prepare inputs for `ExtDataDriver.x` so that ExtData can simulate a real GEOS run +- Added loggers when writing or reading weight files +- Added new option to AGCM.rc `overwrite_checkpoint` to allow checkpoint files to be overwritten. By default still will not overwrite checkpoints +- The trajectory sampler netCDF output variable `location_index_in_iodafile` can be turned off, after we add two control variables: `use_NWP_1_file` and `restore_2_obs_vector` for users. When set to true, the two options will select only one obs file at each Epoch interval, and will rotate the output field index back to the location vector inthe obs file before generating netCDF output. +- Support `splitfield: 1` in HISTORY.rc for trajectory sampler + +### Changed + +- Changed `MAPL_ESMFRegridder` to require the dstMaskValues to be added as grid attribute to use fixed masking, fixes UFS issue +- Increased formatting width of time index in ExtData2G diagnostic print +- Updated GitHub checkout action to use blobless clones +- Update CI to use Baselibs 7.29.0 by default + - This provides ESMF 8.8.0 +- Update `components.yaml` + - `ESMA_env` v4.34.0 + - Update to MPT 2.30 at NAS + - Update to Baselibs 7.29.0 (ESMF 8.8.0) + - `ESMA_cmake` v3.56.0 + - Use `LOCATION` Python `FIND_STRATEGY` + +### Fixed + +- Free MPI communicators after reading and/or writing of restarts +- Fixed the behavior of `MAPL_MaxMin` in presence of NaN +- Fixed bug with return codes and macros in udunits2f + ## [2.51.2] - 2024-12-19 ### Changed diff --git a/CMakeLists.txt b/CMakeLists.txt index f7b5c3ffac4d..83113960be62 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,7 +8,7 @@ endif () project ( MAPL - VERSION 2.51.2 + VERSION 2.52.0 LANGUAGES Fortran CXX C) # Note - CXX is required for ESMF # Set the possible values of build type for cmake-gui diff --git a/Tests/generate_extdatadriver_input.md b/Tests/generate_extdatadriver_input.md new file mode 100755 index 000000000000..8208b0ec9c74 --- /dev/null +++ b/Tests/generate_extdatadriver_input.md @@ -0,0 +1,58 @@ +# Introduction +This is a simple utility to generate inputs for ExtDataDriver.x so that ExtData can simulate a real GEOS run. It requires 3 things that are passed in +- A list of items ExtData needs to fill. This can be found by looking at the GEOS log +- The import spec of the GCM component from using the printspec option to GEOS +- all the needed yaml files in a directory, name of directory is passed + +This will generate the import and optionally the export (as well as History) defition for ExtDataDriver.x to spare the human the tedious work. + +# Example Inputs +To get the list of ExtDate items, just grab all the lines that look like this to a file: +``` + EXTDATA: INFO: ---- 00001: BC_AIRCRAFT + EXTDATA: INFO: ---- 00002: BC_ANTEBC1 + EXTDATA: INFO: ---- 00003: BC_ANTEBC2 + EXTDATA: INFO: ---- 00004: BC_AVIATION_CDS + EXTDATA: INFO: ---- 00005: BC_AVIATION_CRS + EXTDATA: INFO: ---- 00006: BC_AVIATION_LTO + EXTDATA: INFO: ---- 00007: BC_BIOFUEL + EXTDATA: INFO: ---- 00008: BC_BIOMASS + EXTDATA: INFO: ---- 00009: BC_SHIP + EXTDATA: INFO: ---- 00010: BRC_AIRCRAFT + EXTDATA: INFO: ---- 00011: BRC_ANTEBRC1 + EXTDATA: INFO: ---- 00012: BRC_ANTEBRC2 + EXTDATA: INFO: ---- 00013: BRC_AVIATION_CDS + EXTDATA: INFO: ---- 00014: BRC_AVIATION_CRS + EXTDATA: INFO: ---- 00015: BRC_AVIATION_LTO + EXTDATA: INFO: ---- 00016: BRC_BIOFUEL + EXTDATA: INFO: ---- 00017: BRC_BIOMASS + EXTDATA: INFO: ---- 00018: BRC_SHIP + EXTDATA: INFO: ---- 00019: BRC_TERPENE +``` + +To get the GCM component spec, run with `PRINTSPEC: 1` in the `CAP.rc` and copy lines out that look like this: +``` + #IMPORT spec for GCM + #COMPONENT, SHORT_NAME, LONG_NAME, UNIT, DIMS, CONTAINER_TYPE + GENERIC: INFO: GCM, WSUB_CLIM, stdev in vertical velocity, m s-1, 3, esmf_field + GENERIC: INFO: GCM, MEGAN_ORVC, MEGAN_ORVC, kgC/m2/s, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_CROP, CLM4_PFT_CROP, 1, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_C4_GRSS, CLM4_PFT_C4_GRSS, 1, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_C3_NARC_GRSS, CLM4_PFT_C3_NARC_GRSS, 1, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_C3_ARCT_GRSS, CLM4_PFT_C3_ARCT_GRSS, 1, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_BDLF_DECD_BORL_SHRB, CLM4_PFT_BDLF_DECD_BORL_SHRB, 1, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_BDLF_DECD_TMPT_SHRB, CLM4_PFT_BDLF_DECD_TMPT_SHRB, 1, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_BDLF_EVGN_SHRB, CLM4_PFT_BDLF_EVGN_SHRB, 1, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_BDLF_DECD_BORL_TREE, CLM4_PFT_BDLF_DECD_BORL_TREE, 1, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_BDLF_DECD_TMPT_TREE, CLM4_PFT_BDLF_DECD_TMPT_TREE, 1, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_BDLF_DECD_TROP_TREE, CLM4_PFT_BDLF_DECD_TROP_TREE, 1, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_BDLF_EVGN_TMPT_TREE, CLM4_PFT_BDLF_EVGN_TMPT_TREE, 1, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_BDLF_EVGN_TROP_TREE, CLM4_PFT_BDLF_EVGN_TROP_TREE, 1, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_NDLF_DECD_BORL_TREE, CLM4_PFT_NDLF_DECD_BORL_TREE, 1, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_NDLF_EVGN_BORL_TREE, CLM4_PFT_NDLF_EVGN_BORL_TREE, 1, 2, esmf_field + GENERIC: INFO: GCM, CLM4_PFT_NDLF_EVGN_TMPT_TREE, CLM4_PFT_NDLF_EVGN_TMPT_TREE, 1, 2, esmf_field +``` + +Finally just grab the right yaml files for ExtData. + +To run you will of course need to do some further editing of the produced files and link in the acutal data using the same convention `gcm_run.j` does. diff --git a/Tests/generate_extdatadriver_input.py b/Tests/generate_extdatadriver_input.py new file mode 100755 index 000000000000..e1487390b5d4 --- /dev/null +++ b/Tests/generate_extdatadriver_input.py @@ -0,0 +1,167 @@ +#!/usr/bin/env python3 +from yaml import load,dump +import argparse +import os +import yaml +import glob + +dims_dict = {"2":"xy", "3":"xyz"} + +def get_dims(component_map, name): + for comp in component_map: + if name in component_map[comp]: + dims = component_map[comp][name]["dims"] + return dims + +def get_vars_needed(input_file): + + output_list = [] + f = open(input_file,"r") + lines = f.readlines() + f.close() + for line in lines: + temp = line.split() + output_list.append(temp[4]) + + return output_list + +def get_extdata_map(input_dir): + + input_files = glob.glob(input_dir+"/*.yaml") + export_list = {} + for input_file in input_files: + f = open(input_file,'r') + extdata_def = yaml.safe_load(f) + f.close() + for short_name in extdata_def["Exports"]: + export_list.update({short_name:extdata_def["Exports"][short_name]}) + return export_list + +def get_block(cvs_file): + + temp = cvs_file[0].split(' ') + state_type = temp[1][1:].strip() + component = temp[4].strip() + i=2 + for line in cvs_file[2:]: + if "spec for" in line: + break + i=i+1 + return component,state_type,i + + +def get_component_map(cvs_file): + + i_start = 0 + i_end = 0 + n_lines = len(cvs_file) + components = {} + while i_start < n_lines-1: + + comp_name,state_type,i_end = get_block(cvs_file[i_start:]) + comp_map = {} + for i in range(i_end-2): + line = cvs_file[i_start+2+i] + values = line.split(',') + short_name = values[1].strip() + long_name = values[2].strip() + units = values[3].strip() + dims = values[4].strip() + item_type = values[5].strip() + comp_map.update({short_name:{"long_name":long_name,"units":units,"item_type":item_type,"dims":dims}}) + components.update({comp_name+"_"+state_type:comp_map}) + i_start = i_start + i_end + + return components + +def parse_args(): + p = argparse.ArgumentParser(description='Generarte input files for ExtDataDriver to simulate GEOS') + p.add_argument('extdata_provided',type=str,help='a list of items ExtData should fill',default=None) + p.add_argument('spec_def',type=str,help='the GEOS gcm import state from the printspec',default=None) + p.add_argument('extdata_dir',type=str,help='diretory with all the yaml imputs for extdata',default=None) + p.add_argument('-e','--export',action='store_true',help='also include exports for corresponding imports') + + return vars(p.parse_args()) + +if __name__ == '__main__': + + args = parse_args() + + extdata_list = args['extdata_provided'] + do_exports = args['export'] + input_file = args['spec_def'] + f = open(input_file,"r") + input_rc = f.readlines() + f.close() + + extdata_directory = args['extdata_dir'] + extdata_def = get_extdata_map(extdata_directory) + + f_agcm = open("AGCM.rc",'w') +# component + component_map = {} + component_map = get_component_map(input_rc) + + vars_needed = get_vars_needed(extdata_list) + + nl = "\n" + cm = " , " + +# Import state + written = [] + f_agcm.write("IMPORT_STATE::"+nl) + + for item in vars_needed: + if item in extdata_def: + long_name = "NA" + units = "NA" + dims = get_dims(component_map, item) + cdims = dims_dict[dims] + + if item not in written: + f_agcm.write(item+cm+long_name+cm+units+cm+cdims+cm+"c"+nl) + written.append(item) + + f_agcm.write("::"+nl) + +# Export state + if do_exports: + written = [] + f_agcm.write("EXPORT_STATE::"+nl) + for item in vars_needed: + if item in extdata_def: + long_name = "NA" + units = "NA" + dims = get_dims(component_map, item) + cdims = dims_dict[dims] + + if item not in written: + f_agcm.write(item+cm+long_name+cm+units+cm+cdims+cm+"c"+nl) + written.append(item) + + f_agcm.write("::"+nl) + + f_hist = open("HISTORY.rc",'w') + f_hist.write("GRID_LABELS:"+nl) + f_hist.write("::"+nl) + f_hist.write("COLLECTIONS: my_collection"+nl) + f_hist.write("::"+nl) + f_hist.write("my_collection.template: 'nc4'"+nl) + f_hist.write("my_collection.format: 'CFIO'"+nl) + f_hist.write("my_collection.frequency: '240000'"+nl) + first = True + written = [] + for item in vars_needed: + if item in extdata_def: + if item not in written: + if first: + first = False + f_hist.write("my_collection.fields:'"+item+"' , 'Root',"+"\n") + else: + f_hist.write("'"+item+"' , 'Root',"+"\n") + written.append(item) + f_hist.write("::") + + f_hist.close() + f_agcm.close() + diff --git a/base/FileIOShared.F90 b/base/FileIOShared.F90 index 3854e48ddbdb..657b677650a9 100644 --- a/base/FileIOShared.F90 +++ b/base/FileIOShared.F90 @@ -44,6 +44,7 @@ module FileIOSharedMod public ArrDescrInit public ArrDescrCreateReaderComm public ArrDescrCreateWriterComm + public ArrDescrCommFree ! Global vars: ! ------------ @@ -89,9 +90,13 @@ module FileIOSharedMod type ArrDescr integer(kind=MPI_OFFSET_KIND) :: offset character(len=MPI_MAX_INFO_VAL) :: romio_cb_read,cb_buffer_size,romio_cb_write - integer :: Xcomm, Ycomm, NX0, NY0 - integer :: readers_comm, IOscattercomm - integer :: writers_comm, IOgathercomm + integer :: NX0, NY0 + integer :: Xcomm = MPI_COMM_NULL + integer :: Ycomm = MPI_COMM_NULL + integer :: readers_comm = MPI_COMM_NULL + integer :: IOscattercomm = MPI_COMM_NULL + integer :: writers_comm = MPI_COMM_NULL + integer :: IOgathercomm = MPI_COMM_NULL integer :: myrow logical :: split_restart = .false. logical :: split_checkpoint = .false. @@ -539,12 +544,12 @@ subroutine ArrDescrInit(ArrDes,comm,im_world,jm_world,lm_world,nx,ny,num_readers ArrDes%jm_world=jm_world ArrDes%lm_world=lm_world - ArrDes%readers_comm = readers_comm - ArrDes%ioscattercomm = ioscattercomm - ArrDes%writers_comm = writers_comm - ArrDes%iogathercomm = iogathercomm - ArrDes%xcomm = xcomm - ArrDes%ycomm = ycomm + call MAPL_Comm_dup(readers_comm, ArrDes%readers_comm, status) + call MAPL_Comm_dup(ioscattercomm, ArrDes%IOscattercomm, status) + call MAPL_Comm_dup(writers_comm, ArrDes%writers_comm, status) + call MAPL_Comm_dup(iogathercomm, ArrDes%IOgathercomm, status) + call MAPL_Comm_dup(xcomm, ArrDes%Xcomm, status) + call MAPL_Comm_dup(ycomm, ArrDes%Ycomm, status) call mpi_comm_rank(arrdes%ycomm,arrdes%myrow,status) _VERIFY(status) @@ -586,12 +591,21 @@ subroutine ArrDescrSet(ArrDes, offset, & integer, optional, intent(IN ) :: writers_comm, iogathercomm integer, optional, target :: i1(:), in(:), j1(:), jn(:) integer, optional, intent(IN ) :: im_world, jm_world, lm_world - + integer :: status + if(present(offset )) ArrDes%offset = offset - if(present(readers_comm )) ArrDes%readers_comm = readers_comm - if(present(ioscattercomm)) ArrDes%ioscattercomm = ioscattercomm - if(present(writers_comm )) ArrDes%writers_comm = writers_comm - if(present(iogathercomm )) ArrDes%iogathercomm = iogathercomm + if(present(readers_comm )) then + call MAPL_Comm_dup(readers_comm, ArrDes%readers_comm, status) + end if + if(present(ioscattercomm)) then + call MAPL_Comm_dup(IOscattercomm, ArrDes%IOscattercomm, status) + end if + if(present(writers_comm )) then + call MAPL_Comm_dup(writers_comm, ArrDes%writers_comm, status) + end if + if(present(iogathercomm)) then + call MAPL_Comm_dup(IOgathercomm, ArrDes%IOgathercomm, status) + end if if(present(i1 )) ArrDes%i1 => i1 if(present(in )) ArrDes%in => in if(present(j1 )) ArrDes%j1 => j1 @@ -631,7 +645,7 @@ subroutine ArrDescrCreateWriterComm(arrdes, full_comm, num_writers, rc) call MPI_COMM_SPLIT(full_comm, color, myid, arrdes%writers_comm, status) _VERIFY(status) if (num_writers==ny) then - arrdes%IOgathercomm = arrdes%Xcomm + call MAPL_Comm_dup(arrdes%Xcomm, ArrDes%IOgathercomm, status) else j = arrdes%NY0 - mod(arrdes%NY0-1,ny_by_writers) call MPI_COMM_SPLIT(full_comm, j, myid, arrdes%IOgathercomm, status) @@ -679,7 +693,7 @@ subroutine ArrDescrCreateReaderComm(arrdes, full_comm, num_readers, rc) call MPI_COMM_SPLIT(full_comm, color, MYID, arrdes%readers_comm, status) _VERIFY(status) if (num_readers==ny) then - arrdes%IOscattercomm = arrdes%Xcomm + call MAPL_Comm_dup(arrdes%Xcomm, ArrDes%IOscattercomm, status) else j = arrdes%NY0 - mod(arrdes%NY0-1,ny_by_readers) call MPI_COMM_SPLIT(full_comm, j, MYID, arrdes%IOscattercomm, status) @@ -874,4 +888,48 @@ subroutine ArrayScatterShmR4D1(local_array, global_array, grid, mask, rc) _RETURN(ESMF_SUCCESS) end subroutine ArrayScatterShmR4D1 + subroutine ArrDescrCommFree(arrdes, rc) + type(ArrDescr), intent(inout) :: arrdes + integer, optional, intent(out) :: rc + + integer :: status + + call MAPL_CommFree(arrdes%Xcomm, _RC) + call MAPL_CommFree(arrdes%Ycomm, _RC) + call MAPL_CommFree(arrdes%readers_comm, _RC) + call MAPL_CommFree(arrdes%writers_comm, _RC) + call MAPL_CommFree(arrdes%IOgathercomm, _RC) + call MAPL_CommFree(arrdes%IOscattercomm, _RC) + + _RETURN(ESMF_SUCCESS) + end subroutine ArrDescrCommFree + + subroutine MAPL_CommFree(comm, rc) + integer, intent(inout) :: comm + integer, optional, intent(out) :: rc + + integer :: status + + if(comm /= MPI_COMM_NULL) then + call MPI_Comm_Free(comm, status) + _VERIFY(status) + comm = MPI_COMM_NULL + end if + _RETURN(ESMF_SUCCESS) + end subroutine MAPL_CommFree + + subroutine MAPL_Comm_Dup(comm, newcomm, rc) + integer, intent(in) :: comm + integer, intent(out) :: newcomm + integer, optional, intent(out) :: rc + + integer :: status + + newcomm = MPI_COMM_NULL + if(comm /= MPI_COMM_NULL) then + call MPI_Comm_Dup(comm, newcomm,status) + _VERIFY(status) + end if + _RETURN(ESMF_SUCCESS) + end subroutine MAPL_Comm_Dup end module FileIOSharedMod diff --git a/base/MAPL_EsmfRegridder.F90 b/base/MAPL_EsmfRegridder.F90 index 97587cf6421d..fe41619bcd23 100644 --- a/base/MAPL_EsmfRegridder.F90 +++ b/base/MAPL_EsmfRegridder.F90 @@ -14,6 +14,7 @@ module MAPL_EsmfRegridderMod use MAPL_RegridderSpecRouteHandleMap use MAPL_MAPLGrid use MAPL_ConstantsMod + use pFlogger, only: logging, Logger implicit none private @@ -1430,18 +1431,21 @@ subroutine create_route_handle(this, kind, rc) real(real64), pointer :: src_dummy_r8(:,:), dst_dummy_r8(:,:) type (ESMF_Field) :: src_field, dst_field - integer :: srcTermProcessing + integer :: srcTermProcessing, num_mask_values integer, pointer :: factorIndexList(:,:) integer, allocatable :: dstMaskValues(:) real(ESMF_KIND_R8), pointer :: factorList(:) type(ESMF_RouteHandle) :: dummy_rh type(ESMF_UnmappedAction_Flag) :: unmappedaction - logical :: global, isPresent, has_mask + logical :: global, isPresent, has_mask, has_dstMaskValues type(RegridderSpecRouteHandleMap), pointer :: route_handles, transpose_route_handles type(ESMF_RouteHandle) :: route_handle, transpose_route_handle character(len=ESMF_MAXPATHLEN) :: rh_file,rh_trans_file logical :: rh_file_exists, file_weights, compute_transpose + type(Logger), pointer :: lgr + lgr => logging%get_logger('MAPL') + if (kind == ESMF_TYPEKIND_R4) then route_handles => route_handles_r4 transpose_route_handles => transpose_route_handles_r4 @@ -1468,6 +1472,7 @@ subroutine create_route_handle(this, kind, rc) rh_file_exists = .false. end if if (rh_file_exists) then + call lgr%info('Reading weight file: %a', trim(rh_file)) route_handle = ESMF_RouteHandleCreate(rh_file,_RC) call route_handles%insert(spec, route_handle) if (compute_transpose) then @@ -1507,6 +1512,13 @@ subroutine create_route_handle(this, kind, rc) end if call ESMF_GridGetItem(spec%grid_out,itemflag=ESMF_GRIDITEM_MASK, & staggerloc=ESMF_STAGGERLOC_CENTER, isPresent = has_mask, _RC) + call ESMF_AttributeGet(spec%grid_out, name=MAPL_DESTINATIONMASK, isPresent=has_dstMaskValues, _RC) + if (has_dstMaskValues) then + _ASSERT(has_mask, "masking destination values when no masks is present") + call ESMF_AttributeGet(spec%grid_out, name=MAPL_DESTINATIONMASK, itemcount=num_mask_values, _RC) + allocate(dstMaskValues(num_mask_values), _STAT) + call ESMF_AttributeGet(spec%grid_out, name=MAPL_DESTINATIONMASK, valuelist=dstMaskValues, _RC) + end if counter = counter + 1 @@ -1516,7 +1528,6 @@ subroutine create_route_handle(this, kind, rc) call ESMF_AttributeGet(spec%grid_in, name='Global',value=global,rc=status) if (.not.global) unmappedaction=ESMF_UNMAPPEDACTION_IGNORE end if - if (has_mask) dstMaskValues = [MAPL_MASK_OUT] ! otherwise unallocated select case (spec%regrid_method) case (REGRID_METHOD_BILINEAR, REGRID_METHOD_BILINEAR_MONOTONIC) @@ -1586,8 +1597,10 @@ subroutine create_route_handle(this, kind, rc) call ESMF_FieldDestroy(dst_field, rc=status) _VERIFY(status) if (file_weights) then + call lgr%info('Writing weight file: %a', trim(rh_file)) call ESMF_RouteHandleWrite(route_handle,rh_file,_RC) if (compute_transpose) then + call lgr%info('Writing transpose weight file: %a', trim(rh_trans_file)) call ESMF_RouteHandleWrite(transpose_route_handle,rh_trans_file,_RC) end if end if diff --git a/base/MAPL_LocStreamMod.F90 b/base/MAPL_LocStreamMod.F90 index d0f654f73900..37a8963cc2d2 100644 --- a/base/MAPL_LocStreamMod.F90 +++ b/base/MAPL_LocStreamMod.F90 @@ -69,8 +69,8 @@ module MAPL_LocStreamMod !EOP -integer, parameter :: NumGlobalVars=4 -integer, parameter :: NumLocalVars =4 +!integer, parameter :: NumGlobalVars=4 +!integer, parameter :: NumLocalVars =4 type MAPL_GeoLocation @@ -343,18 +343,19 @@ subroutine MAPL_LocStreamCreateFromFile(LocStream, LAYOUT, FILENAME, NAME, MASK, integer :: UNIT integer :: N, I, K, L, NT type(MAPL_LocStreamType), pointer :: STREAM - real, pointer :: AVR(:,:), AVR_transpose(:,:) + real, pointer :: AVR(:,:) logical, pointer :: MSK(:) real :: X, Y, X0, Y0, XE, DX, DY integer :: II, JJ logical :: DoCoeffs character(len=MAPL_TileNameLength):: gname - + character(len=MAPL_TileNameLength):: GridNames(2) + integer :: IMs(2), JMs(2) integer :: irec integer(kind=1) :: byte(4) - integer :: I1, IN, J1, JN - integer :: iostat - logical :: isascii + integer :: I1, IN, J1, JN, N_Grids, N_PfafCat + integer :: iostat, filetype + logical :: isascii, isnc4, isbinary logical :: read_always logical, pointer :: ISMINE(:) type(MAPL_Tiling ), pointer :: TILING @@ -404,94 +405,63 @@ subroutine MAPL_LocStreamCreateFromFile(LocStream, LAYOUT, FILENAME, NAME, MASK, ! Use some heuristics to determine filetype (choices are BINARY and ASCII) !------------------------------------------------------------------------ - ! just get the unit - UNIT = GETFILE(FILENAME, DO_OPEN=0, ALL_PES=.true., RC=STATUS) - _VERIFY(STATUS) + call MAPL_NCIOGetFileType(FILENAME, filetype, _RC) - INQUIRE(IOLENGTH=IREC) BYTE - open (UNIT=UNIT, FILE=FILENAME, FORM='unformatted', ACCESS='DIRECT', RECL=IREC, IOSTAT=status) - _VERIFY(STATUS) - read (UNIT, REC=1, ERR=100) BYTE - call FREE_FILE(UNIT) - - isascii = .true. - do i=1,size(byte) - if (BYTE(I) < 7) then - isascii = .false. - exit - end if - end do + isnc4 = (filetype == MAPL_FILETYPE_NC4) + isascii = (filetype == MAPL_FILETYPE_TXT) + isbinary= (filetype == MAPL_FILETYPE_BIN) - if (isascii) then + if ( .not. isbinary) then ! Open file and read header info !------------------------------- - UNIT = GETFILE(FILENAME, form='FORMATTED', RC=status) - _VERIFY(STATUS) - -! Total number of tiles in exchange grid -!--------------------------------------- + if (isnc4) then + if (MAPL_AM_I_root()) then + call MAPL_ReadTilingNC4(FILENAME, GridName=GridNames, IM=IMs, JM=JMs, N_Grids=N_Grids, N_PfafCat=N_PfafCat, AVR=AVR, _RC) + NT = size(AVR,1) + endif + call MAPL_CommsBcast(layout, NT, 1, MAPL_Root, status) + call MAPL_CommsBcast(layout, N_Grids, 1, MAPL_Root, status) + call MAPL_CommsBcast(layout, IMs, 2, MAPL_Root, status) + call MAPL_CommsBcast(layout, JMs, 2, MAPL_Root, status) - if (use_pfaf_) then - call READ_PARALLEL(layout, hdr, UNIT=UNIT, rc=status) - _VERIFY(STATUS) - nt=hdr(1) - stream%pfafstetter_catchments=hdr(2) - else - call READ_PARALLEL(layout, nt, UNIT=UNIT, rc=status) - _VERIFY(STATUS) - end if + if (use_pfaf_) then + call MAPL_CommsBcast(layout, N_PfafCat, 1, MAPL_Root, status) + endif -! Number of grids that can be attached -!------------------------------------- + do N = 1, N_Grids + call MAPL_CommsBcast(layout, GridNames(N), MAPL_TileNameLength, MAPL_Root, status) + enddo - call READ_PARALLEL(layout, STREAM%N_GRIDS, unit=UNIT, rc=status) - _VERIFY(STATUS) + if (.not. MAPL_AM_I_root()) then + allocate(AVR(NT, NumGlobalVars+NumLocalVars*N_GRIDS)) + endif + call MAPL_CommsBcast(layout, AVR, NT*(NumGlobalVars+NumLocalVars*N_GRIDS), MAPL_Root, status) + endif -! The exchange grid is used to tile each attached grid -!----------------------------------------------------- + if (isascii) then + call MAPL_ReadTilingASCII(layout, FILENAME, GridNames, NT, IMs, JMs, N_Grids, N_PfafCat, AVR, _RC) + endif + STREAM%N_GRIDS = N_Grids allocate(STREAM%TILING(STREAM%N_GRIDS), STAT=STATUS) _VERIFY(STATUS) -! The names and sizes of the grids to be tiled -!--------------------------------------------- - - do N=1,STREAM%N_GRIDS - call READ_PARALLEL(layout, STREAM%TILING(N)%NAME, unit=UNIT, rc=status) - _VERIFY(STATUS) + do N = 1, N_Grids + STREAM%TILING(N)%NAME = GridNames(N) if (NewGridNames_) then - call GenOldGridName_(STREAM%TILING(N)%NAME) + call GenOldGridName_(STREAM%TILING(N)%NAME) end if - call READ_PARALLEL(layout, STREAM%TILING(N)%IM, unit=UNIT, rc=status) - _VERIFY(STATUS) - call READ_PARALLEL(layout, STREAM%TILING(N)%JM, unit=UNIT, rc=status) - _VERIFY(STATUS) + STREAM%TILING(N)%IM = IMs(N) + STREAM%TILING(N)%JM = JMs(N) enddo if (use_pfaf_) then + stream%pfafstetter_catchments = N_PfafCat STREAM%TILING(2)%IM = stream%pfafstetter_catchments STREAM%TILING(2)%JM = 1 STREAM%TILING(2)%name = "CATCHMENT_GRID" end if - -! Read location stream file into AVR -!--------------------------------------- - if(index(STREAM%TILING(1)%NAME,'EASE') /=0 ) then - allocate(AVR(NT,9), STAT=STATUS) ! 9 columns for EASE grid - _VERIFY(STATUS) - allocate(AVR_transpose(9,NT)) - else - allocate(AVR(NT,NumGlobalVars+NumLocalVars*STREAM%N_GRIDS), STAT=STATUS) - _VERIFY(STATUS) - allocate(AVR_transpose(NumGlobalVars+NumLocalVars*STREAM%N_GRIDS,NT), STAT=STATUS) - _VERIFY(STATUS) - endif - - call READ_PARALLEL(layout, AVR_transpose(:,:), unit=UNIT, rc=status) - AVR= transpose(AVR_transpose) - deallocate(AVR_transpose) - ! adjust EASE grid starting index. Internally, the starting index is 1 instead of 0. do N=1,STREAM%N_GRIDS if(index(STREAM%TILING(N)%NAME,'EASE') /=0 ) then @@ -499,12 +469,11 @@ subroutine MAPL_LocStreamCreateFromFile(LocStream, LAYOUT, FILENAME, NAME, MASK, AVR(:,NumGlobalVars+2+NumLocalVars*(N-1)) = AVR(:,NumGlobalVars+2+NumLocalVars*(N-1))+1 endif enddo - + !if (use_pfaf_) then !AVR(:,NumGlobalVars+2+NumLocalVars) = 1 !end if - call FREE_FILE(UNIT) ! Allocate msk for which tiles to include in the stream being created. !-------------------------------------------------------------------- @@ -880,7 +849,7 @@ subroutine MAPL_LocStreamCreateFromFile(LocStream, LAYOUT, FILENAME, NAME, MASK, ! and are expensive to reread. !------------------------------------------------------------------ - if(.not.isascii) then + if( isbinary ) then if ( MAPL_am_I_root() ) then rewind(UNIT) diff --git a/base/MAPL_MaxMinMod.F90 b/base/MAPL_MaxMinMod.F90 index 8aa0bef53b4e..e63f45866191 100644 --- a/base/MAPL_MaxMinMod.F90 +++ b/base/MAPL_MaxMinMod.F90 @@ -8,7 +8,7 @@ ! ! Author: GMAO SI-Team ! -! `MAPL_MaxMinMo` --- Global Max/Min of Arrays +! `MAPL_MaxMinMod` --- Global Max/Min of Arrays ! ! This module implements functions for calculating/printing out the global min/max ! of fortran arrays. Derived from GEOS-4 pmaxmin() functions. @@ -20,6 +20,7 @@ module MAPL_MaxMinMod Use ESMF Use MAPL_CommsMod + Use mpi implicit None @@ -66,21 +67,19 @@ subroutine pmaxmin2d_r4 ( qname, a, pmin_, pmax_, fac_ ) real(ESMF_KIND_R4) :: pmax, pmin, fac - integer :: im, jt - integer :: i, j, two=2 - real, allocatable :: qmin(:), qmax(:) real pm1(2) real pm_res(2) type(ESMF_VM) :: vm character(len=16) :: name +!NOTE: the current version does not trap error conditions returned in status integer :: status - - im = size(a,1) - jt = size(a,2) - allocate(qmin(jt),qmax(jt)) + integer :: comm + logical :: has_nans + logical :: has_nans_local + character(len=:), allocatable :: buf if ( present(fac_) ) then fac = fac_ @@ -90,25 +89,10 @@ subroutine pmaxmin2d_r4 ( qname, a, pmin_, pmax_, fac_ ) call ESMF_VmGetCurrent(vm=vm, rc=status) - do j=1,jt - pmax = a(1,j) - pmin = a(1,j) - do i=2,im - pmax = max(pmax, a(i,j)) - pmin = min(pmin, a(i,j)) - enddo - qmax(j) = pmax - qmin(j) = pmin - enddo -! -! Now find max/min of amax/amin -! - pmax = qmax(1) - pmin = qmin(1) - do j=2,jt - pmax = max(pmax, qmax(j)) - pmin = min(pmin, qmin(j)) - enddo + has_nans_local = any(a /= a) + + pmin = minval(a, mask=(a==a)) + pmax = maxval(a, mask=(a==a)) pm1(1) = pmax pm1(2) = -pmin @@ -118,13 +102,19 @@ subroutine pmaxmin2d_r4 ( qname, a, pmin_, pmax_, fac_ ) if ( present(pmax_) ) pmax_ = pmax if ( present(pmin_) ) pmin_ = pmin - deallocate(qmax,qmin) + + call ESMF_VmGet(VM, mpicommunicator=comm, rc=status) + call MPI_Reduce(has_nans_local, has_nans, 1, MPI_LOGICAL, MPI_LOR, 0, comm, status) if ( fac /= 0.0 ) then ! trick to prevent printing - if ( MAPL_am_I_root() ) then + if ( MAPL_am_I_root(vm) ) then name = ' ' name(1:len(qname)) = qname - write(*,*) name, ' max = ', pmax*fac, ' min = ', pmin*fac + buf = "" + if (has_nans) then + buf = " has NaN" + end if + write(*,*) name, ' max = ', pmax*fac, ' min = ', pmin*fac, buf return end if end if diff --git a/base/MAPL_ObsUtil.F90 b/base/MAPL_ObsUtil.F90 index 1b96fad92ae8..aed4c9adcd59 100644 --- a/base/MAPL_ObsUtil.F90 +++ b/base/MAPL_ObsUtil.F90 @@ -26,6 +26,8 @@ module MAPL_ObsUtilMod type :: obs_unit integer :: nobs_epoch integer :: ngeoval + integer :: count_location_until_matching_file + integer :: count_location_in_matching_file logical :: export_all_geoval type(FileMetadata), allocatable :: metadata type(NetCDF4_FileFormatter), allocatable :: file_handle @@ -38,6 +40,7 @@ module MAPL_ObsUtilMod real(kind=REAL64), allocatable :: lats(:) real(kind=REAL64), allocatable :: times_R8(:) integer, allocatable :: location_index_ioda(:) + integer, allocatable :: restore_index(:) real(kind=REAL32), allocatable :: p2d(:) real(kind=REAL32), allocatable :: p3d(:,:) end type obs_unit @@ -50,7 +53,7 @@ module MAPL_ObsUtilMod character (len=ESMF_MAXSTR) :: var_name_time='' character (len=ESMF_MAXSTR) :: file_name_template='' integer :: ngeoval=0 - integer :: nentry_name=0 + integer :: nfield_name_mx=12 character (len=ESMF_MAXSTR), allocatable :: field_name(:,:) !character (len=ESMF_MAXSTR), allocatable :: field_name(:) end type obs_platform @@ -768,7 +771,7 @@ function copy_platform_nckeys(a, rc) copy_platform_nckeys%var_name_lon = a%var_name_lon copy_platform_nckeys%var_name_lat = a%var_name_lat copy_platform_nckeys%var_name_time = a%var_name_time - copy_platform_nckeys%nentry_name = a%nentry_name + copy_platform_nckeys%nfield_name_mx = a%nfield_name_mx _RETURN(_SUCCESS) end function copy_platform_nckeys @@ -781,7 +784,7 @@ function union_platform(a, b, rc) integer, optional, intent(out) :: rc character (len=ESMF_MAXSTR), allocatable :: field_name_loc(:,:) - integer :: nfield, nentry_name + integer :: nfield, nfield_name_mx integer, allocatable :: tag(:) integer :: i, j, k integer :: status @@ -802,9 +805,9 @@ function union_platform(a, b, rc) enddo union_platform%ngeoval=k nfield=k - nentry_name=union_platform%nentry_name + nfield_name_mx=union_platform%nfield_name_mx if ( allocated (union_platform%field_name) ) deallocate(union_platform%field_name) - allocate(union_platform%field_name(nentry_name, nfield)) + allocate(union_platform%field_name(nfield_name_mx, nfield)) do i=1, a%ngeoval union_platform%field_name(:,i) = a%field_name(:,i) enddo @@ -950,6 +953,7 @@ subroutine fglob(search_name, filename, rc) ! give the last name if (lenmax < slen) then if (MAPL_AM_I_ROOT()) write(6,*) 'pathlen vs filename_max_char_len: ', slen, lenmax _FAIL ('PATHLEN is greater than filename_max_char_len') + STOP 'lenmax < slen' end if if (slen>0) filename(1:slen)=c_filename(1:slen) diff --git a/base/MAPL_XYGridFactory.F90 b/base/MAPL_XYGridFactory.F90 index 1bcbc57ea3ce..a451bbcd7ee0 100644 --- a/base/MAPL_XYGridFactory.F90 +++ b/base/MAPL_XYGridFactory.F90 @@ -1037,6 +1037,9 @@ subroutine add_mask(this,grid,rc) mask = MAPL_MASK_IN where(fptr==MAPL_UNDEF) mask = MAPL_MASK_OUT + call ESMF_AttributeSet(grid, name=MAPL_DESTINATIONMASK, & + itemCount=1, valueList=[MAPL_MASK_OUT], _RC) + _RETURN(_SUCCESS) end subroutine add_mask diff --git a/base/NCIO.F90 b/base/NCIO.F90 index 928566cadc4e..5ad1fc5bbd3c 100644 --- a/base/NCIO.F90 +++ b/base/NCIO.F90 @@ -22,6 +22,8 @@ module NCIOMod use MAPL_ExceptionHandling use netcdf use pFIO + use BinIOMod, only: GETFILE, READ_PARALLEL, FREE_FILE + use MAPL_Constants !use pFIO_ClientManagerMod use gFTL_StringIntegerMap use gFTL_StringVector @@ -43,6 +45,8 @@ module NCIOMod public MAPL_NCIOGetFileType public MAPL_VarReadNCPar public MAPL_VarWriteNCPar + public MAPL_ReadTilingNC4 + public MAPL_ReadTilingASCII interface MAPL_VarReadNCPar module procedure MAPL_StateVarReadNCPar @@ -78,7 +82,10 @@ module NCIOMod module procedure MAPL_VarWriteNCpar_R8_4d end interface - contains + integer, parameter, public :: NumGlobalVars=4 + integer, parameter, public :: NumLocalVars =4 + +contains subroutine MAPL_FieldReadNCPar(formatter,name,FIELD, ARRDES, HomePE, RC) @@ -3512,11 +3519,12 @@ subroutine MAPL_ArrayReadNCpar_3d(varn,filename,farrayPtr,arrDes,rc) _RETURN(ESMF_SUCCESS) end subroutine MAPL_ArrayReadNCpar_3d - subroutine MAPL_BundleWriteNCPar(Bundle, arrdes, CLOCK, filename, oClients, rc) + subroutine MAPL_BundleWriteNCPar(Bundle, arrdes, CLOCK, filename, clobber, oClients, rc) type(ESMF_FieldBundle), intent(inout) :: Bundle type(ArrDescr), intent(inout) :: arrdes type(ESMF_Clock), intent(in) :: CLOCK character(len=*), intent(in ) :: filename + logical, intent(in) :: clobber type (ClientManager), optional, intent(inout) :: oClients integer, optional, intent(out) :: rc @@ -3582,6 +3590,8 @@ subroutine MAPL_BundleWriteNCPar(Bundle, arrdes, CLOCK, filename, oClients, rc) type(ESMF_Field) :: lons_field, lats_field logical :: isGridCapture, have_oclients real(kind=ESMF_KIND_R8), pointer :: grid_lons(:,:), grid_lats(:,:), lons_field_ptr(:,:), lats_field_ptr(:,:) + integer :: pfio_mode + have_oclients = present(oClients) call ESMF_FieldBundleGet(Bundle,FieldCount=nVars, name=BundleName, rc=STATUS) @@ -4174,9 +4184,11 @@ subroutine MAPL_BundleWriteNCPar(Bundle, arrdes, CLOCK, filename, oClients, rc) else + pfio_mode = PFIO_NOCLOBBER + if (clobber) pfio_mode = PFIO_CLOBBER if (arrdes%writers_comm /= mpi_comm_null) then if (arrdes%num_writers == 1) then - call formatter%create(trim(filename), rc=status) + call formatter%create(trim(filename), mode=pfio_mode, rc=status) _VERIFY(status) call formatter%write(cf,rc=status) _VERIFY(STATUS) @@ -4189,7 +4201,7 @@ subroutine MAPL_BundleWriteNCPar(Bundle, arrdes, CLOCK, filename, oClients, rc) _VERIFY(status) call cf%add_attribute("Split_Cubed_Sphere", writer_rank, _RC) else - call formatter%create_par(trim(filename),comm=arrdes%writers_comm,info=info,rc=status) + call formatter%create_par(trim(filename),mode=pfio_mode,comm=arrdes%writers_comm,info=info,rc=status) _VERIFY(status) endif call formatter%write(cf,rc=status) @@ -4307,13 +4319,14 @@ end subroutine add_fvar end subroutine MAPL_BundleWriteNCPar - subroutine MAPL_StateVarWriteNCPar(filename, STATE, ARRDES, CLOCK, NAME, forceWriteNoRestart, oClients, RC) + subroutine MAPL_StateVarWriteNCPar(filename, STATE, ARRDES, CLOCK, NAME, forceWriteNoRestart, clobber, oClients, RC) character(len=*) , intent(IN ) :: filename type (ESMF_State) , intent(IN ) :: STATE type(ArrDescr) , intent(INOUT) :: ARRDES type(ESMF_Clock) , intent(IN ) :: CLOCK character(len=*), optional, intent(IN ) :: NAME logical, optional, intent(IN ) :: forceWriteNoRestart + logical, optional, intent(in ) :: clobber type (ClientManager), optional, intent(inout) :: oClients integer, optional, intent( OUT) :: RC @@ -4339,6 +4352,7 @@ subroutine MAPL_StateVarWriteNCPar(filename, STATE, ARRDES, CLOCK, NAME, forceWr logical :: is_test_framework, isGridCapture integer :: fieldIsValid type(ESMF_Array) :: array + logical :: local_clobber call ESMF_StateGet(STATE,ITEMCOUNT=ITEMCOUNT,RC=STATUS) _VERIFY(STATUS) @@ -4358,9 +4372,9 @@ subroutine MAPL_StateVarWriteNCPar(filename, STATE, ARRDES, CLOCK, NAME, forceWr _VERIFY(STATUS) forceWriteNoRestart_ = .false. - if(present(forceWriteNoRestart)) then - forceWriteNoRestart_ = forceWriteNoRestart - endif + if(present(forceWriteNoRestart)) forceWriteNoRestart_ = forceWriteNoRestart + local_clobber = .false. + if (present(clobber)) local_clobber = clobber if(present(NAME)) then DOIT = ITEMNAMES==NAME @@ -4494,7 +4508,7 @@ subroutine MAPL_StateVarWriteNCPar(filename, STATE, ARRDES, CLOCK, NAME, forceWr call ESMF_AttributeSet(bundle_write, name="MAPL_GridCapture", value=isGridCapture, _RC) end if - call MAPL_BundleWriteNCPar(Bundle_Write, arrdes, CLOCK, filename, oClients=oClients, rc=status) + call MAPL_BundleWriteNCPar(Bundle_Write, arrdes, CLOCK, filename, clobber=local_clobber, oClients=oClients, rc=status) _VERIFY(STATUS) _RETURN(ESMF_SUCCESS) @@ -4502,8 +4516,12 @@ subroutine MAPL_StateVarWriteNCPar(filename, STATE, ARRDES, CLOCK, NAME, forceWr end subroutine MAPL_StateVarWriteNCPar subroutine MAPL_NCIOGetFileType(filename,filetype,rc) - implicit none + ! filetype = 0, hdf5 + ! filetype = 1, ascii + ! filetype = 2, binary or unknown + ! filetype = -1, unknown + implicit none ! Arguments !---------- character(len=*), intent(IN ) :: filename @@ -4521,7 +4539,7 @@ subroutine MAPL_NCIOGetFileType(filename,filetype,rc) integer :: irec integer :: unit integer :: i, cwrd - logical :: typehdf5 + logical :: typehdf5, isascii character(len=12) :: fmt INQUIRE(IOLENGTH=IREC) WORD @@ -4533,9 +4551,9 @@ subroutine MAPL_NCIOGetFileType(filename,filetype,rc) read (UNIT, REC=2, ERR=100) TwoWords(5:8) close(UNIT) - typehdf5 = .true. - filetype = -1 ! Unknown + filetype = MAPL_FILETYPE_UNK ! Unknown + typehdf5 = .true. do i = 1, 8 if (iachar(TwoWords(i)) /= hdf5(i)) then typehdf5 = .false. @@ -4543,20 +4561,23 @@ subroutine MAPL_NCIOGetFileType(filename,filetype,rc) end if end do if (typehdf5) then - filetype = 0 ! HDF5 + filetype = MAPL_FILETYPE_NC4 _RETURN(ESMF_SUCCESS) - - end if - ! Attempt to identify as fortran binary - cwrd = transfer(TwoWords(1:4), irec) - ! check if divisible by 4 - irec = cwrd/4 - filetype = irec - if (cwrd /= 4*irec) then - _RETURN(ESMF_FAILURE) end if - filetype = -1 + isascii = .true. + do i = 1, 4 + if (iachar(TwoWords(i)) < 7) then + isascii = .false. + exit + end if + end do + if (isascii) then + filetype = MAPL_FILETYPE_TXT + _RETURN(ESMF_SUCCESS) + endif + + filetype = MAPL_FILETYPE_BIN _RETURN(ESMF_SUCCESS) 100 continue @@ -5057,5 +5078,241 @@ function create_flipped_field(field,rc) result(flipped_field) _RETURN(_SUCCESS) end function create_flipped_field + subroutine MAPL_ReadTilingNC4(File, GridName, im, jm, nx, ny, n_Grids, iTable, rTable, N_PfafCat, AVR,rc) + character(*), intent(IN) :: File + character(*), optional, intent(out) :: GridName(:) + integer, optional, intent(out) :: IM(:), JM(:) + integer, optional, intent(out) :: nx, ny, n_Grids + integer, optional, allocatable, intent(out) :: iTable(:,:) + real(kind=REAL64), optional, allocatable, intent(out) :: rTable(:,:) + integer, optional, intent(out) :: N_PfafCat + real, optional, pointer, intent(out) :: AVR(:,:) ! used by GEOSgcm + integer, optional, intent(out) :: rc + + type (Attribute), pointer :: ref + character(len=:), allocatable :: attr + type (NetCDF4_FileFormatter) :: formatter + type (FileMetadata) :: meta + character(len=1) :: str_num + integer :: ng, ntile, status, ll + class(*), pointer :: attr_val(:) + class(*), pointer :: char_val + integer, allocatable :: tmp_int(:) + real(kind=REAL64),allocatable :: fr(:) + + integer :: NumCol + integer, allocatable :: iTable_(:,:) + real(kind=REAL64), allocatable :: rTable_(:,:) + + call formatter%open(File, pFIO_READ, rc=status) + meta = formatter%read(rc=status) + + ng = 1 + if (meta%has_attribute("Grid1_Name")) ng = 2 ! for ng=1 (i.e., EASE), attribute is just "Grid_Name" + ntile = meta%get_dimension('tile') + + if (present(n_Grids)) then + n_Grids = ng + endif + + if (present(nx)) then + ref => meta%get_attribute('raster_nx') + attr_val => ref%get_values() + select type(attr_val) + type is (integer(INT32)) + nx = attr_val(1) + endselect + endif + if (present(ny)) then + ref => meta%get_attribute('raster_ny') + attr_val => ref%get_values() + select type (attr_val) + type is (integer(INT32)) + ny = attr_val(1) + endselect + endif + + if (present(N_PfafCat)) then + ref => meta%get_attribute('N_PfafCat') + attr_val => ref%get_values() + select type (attr_val) + type is (integer(INT32)) + N_PfafCat = attr_val(1) + endselect + endif + do ll = 1, ng + if (ng == 1) then + str_num = '' + else + write(str_num, '(i0)') ll + endif + + if (present(GridName)) then + attr = 'Grid'//trim(str_num)//'_Name' + ref =>meta%get_attribute(attr) + char_val => ref%get_value() + select type(char_val) + type is(character(*)) + GridName(ll) = char_val + class default + print('unsupported subclass (not string) of attribute named '//attr) + end select + endif + if (present(IM)) then + attr = 'IM'//trim(str_num) + ref =>meta%get_attribute(attr) + attr_val => ref%get_values() + select type(attr_val) + type is( integer(INT32)) + IM(ll) = attr_val(1) + end select + endif + if (present(JM)) then + attr = 'JM'//trim(str_num) + ref =>meta%get_attribute(attr) + attr_val => ref%get_values() + select type(attr_val) + type is(integer(INT32)) + JM(ll) = attr_val(1) + end select + endif + enddo + + if (present(iTable) .or. present(AVR) ) then + allocate(iTable_(ntile,0:7)) + allocate(tmp_int(ntile)) + call formatter%get_var('typ', iTable_(:,0)) + do ll = 1, ng + if (ng == 1) then + str_num='' + else + write(str_num, '(i0)') ll + endif + + call formatter%get_var('i_indg' //trim(str_num), tmp_int, rc=status) + iTable_(:,ll*2) = tmp_int + call formatter%get_var('j_indg' //trim(str_num), tmp_int, rc=status) + iTable_(:,ll*2+1) = tmp_int + call formatter%get_var('pfaf_index'//trim(str_num), tmp_int, rc=status) + if ( ng == 1) then + iTable_(:,4) = tmp_int + else + iTable_(:,5+ll) = tmp_int + endif + enddo + endif + if (present(rTable) .or. present(AVR) ) then + allocate(rTable_(ntile,10)) + call formatter%get_var('com_lon', rTable_(:,1), rc=status) + call formatter%get_var('com_lat', rTable_(:,2), rc=status) + call formatter%get_var('area', rTable_(:,3), rc=status) + do ll = 1, ng + if (ng == 1) then + str_num='' + else + write(str_num, '(i0)') ll + endif + call formatter%get_var('frac_cell' //trim(str_num), rTable_(:,3+ll), rc=status) + enddo + call formatter%get_var('min_lon', rTable_(:, 6), rc=status) + call formatter%get_var('max_lon', rTable_(:, 7), rc=status) + call formatter%get_var('min_lat', rTable_(:, 8), rc=status) + call formatter%get_var('max_lat', rTable_(:, 9), rc=status) + call formatter%get_var('elev', rTable_(:,10), rc=status) + endif + if (present(AVR)) then + ! In GEOSgcm, it already assumes ng = 2, so NumCol = 10 + NumCol = NumGlobalVars+NumLocalVars*ng + allocate(AVR(ntile, NumCol)) + AVR(:, 1) = iTable_(:,0) + ! for EASE grid, the second collum is replaced by the area + AVR(:, 2) = rTable_(:,3) + AVR(:, 3) = rTable_(:,1) + AVR(:, 4) = rTable_(:,2) + + AVR(:, 5) = iTable_(:,2) + AVR(:, 6) = iTable_(:,3) + AVR(:, 7) = rTable_(:,4) + if (ng == 1) then + AVR(:,8) = iTable_(:,4) + else + AVR(:, 8) = iTable_(:,6) + + AVR(:, 9) = iTable_(:,4) + AVR(:, 10) = iTable_(:,5) + AVR(:, 11) = rTable_(:,5) + AVR(:, 12) = iTable_(:,7) + endif + endif + + if (present(iTable)) then + call move_alloc(iTable_, iTable) + endif + + if (present(rTable)) then + call move_alloc(rTable_, rTable) + do ll = 1, ng + where ( rTable(:,3+ll) /=0.0 ) rTable(:,3+ll) = rTable(:,3)/rTable(:,3+ll) + enddo + endif + _RETURN(ESMF_SUCCESS) + end subroutine MAPL_ReadTilingNC4 + + subroutine MAPL_ReadTilingASCII(layout, FileName, GridName, NT, im, jm, n_Grids, N_PfafCat, AVR,rc) + type(ESMF_DELayout), intent(IN) :: LAYOUT + character(*), intent(IN) :: FileName + character(*), intent(out) :: GridName(:) + integer, intent(out) :: NT + integer, intent(out) :: IM(:), JM(:) + integer, intent(out) :: n_Grids + integer, intent(out) :: N_PfafCat + real, pointer, intent(out) :: AVR(:,:) ! used by GEOSgcm + integer, optional, intent(out) :: rc + + integer :: unit, status, hdr(2), N + real, pointer :: AVR_Transpose(:,:) + + UNIT = GETFILE(FILENAME, form='FORMATTED', RC=status) + _VERIFY(STATUS) + +! Total number of tiles in exchange grid +!--------------------------------------- + + call READ_PARALLEL(layout, hdr, UNIT=UNIT, rc=status) + _VERIFY(STATUS) + NT = hdr(1) + N_PfafCat = hdr(2) + + call READ_PARALLEL(layout, N_GRIDS, unit=UNIT, rc=status) + _VERIFY(STATUS) + + do N = 1, N_GRIDS + call READ_PARALLEL(layout, GridName(N), unit=UNIT, rc=status) + _VERIFY(STATUS) + call READ_PARALLEL(layout, IM(N), unit=UNIT, rc=status) + _VERIFY(STATUS) + call READ_PARALLEL(layout, JM(N), unit=UNIT, rc=status) + _VERIFY(STATUS) + enddo + + if(index(GridNAME(1),'EASE') /=0 ) then + allocate(AVR(NT,9), STAT=STATUS) ! 9 columns for EASE grid + _VERIFY(STATUS) + allocate(AVR_transpose(9,NT)) + else + allocate(AVR(NT,NumGlobalVars+NumLocalVars*N_GRIDS), STAT=STATUS) + _VERIFY(STATUS) + allocate(AVR_transpose(NumGlobalVars+NumLocalVars*N_GRIDS,NT), STAT=STATUS) + _VERIFY(STATUS) + endif + + call READ_PARALLEL(layout, AVR_transpose(:,:), unit=UNIT, rc=status) + AVR = transpose(AVR_transpose) + deallocate(AVR_transpose) + call FREE_FILE(UNIT) + + _RETURN(ESMF_SUCCESS) + + end subroutine MAPL_ReadTilingASCII end module NCIOMod diff --git a/base/Plain_netCDF_Time.F90 b/base/Plain_netCDF_Time.F90 index b9c163816647..924ba9323978 100644 --- a/base/Plain_netCDF_Time.F90 +++ b/base/Plain_netCDF_Time.F90 @@ -840,6 +840,13 @@ subroutine split_string_by_space (string_in, length_mx, & wc=0 ios=0 string = trim( adjustl(string_in) ) + str_piece(:)='' + i = index (string, mark) + if (i==0) then + nseg=1 + str_piece(1)=string + return + end if do while (ios==0) i = index (string, mark) !!print*, 'index=', i @@ -858,5 +865,47 @@ subroutine split_string_by_space (string_in, length_mx, & end subroutine split_string_by_space + subroutine split_string_by_seperator (string_in, length_mx, seperator_in, & + mxseg, nseg, str_piece, jstatus) + character (len=length_mx), intent (in) :: string_in + integer, intent (in) :: length_mx + character (len=length_mx), intent (in) :: seperator_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_sc, string_oper, string_aux + character (len=1) :: mark, CH + integer :: ios + integer :: wc + integer :: len1, len2 + ! + ! "xxxx; yy; zz; uu, vv," + ! seperator = ";," + ! + + !__ s1. replace seperator by space + ! + string_sc = trim( adjustl(string_in) ) + string_oper = trim( adjustl(seperator_in) ) + len1 = len_trim(string_sc) + len2 = len_trim(string_oper) + string_aux=string_sc + do i = 1, len1 + CH = string_sc(i:i) + do j = 1, len2 + mark = string_oper(j:j) + if (CH==mark) then + string_aux(i:i)=' ' + end if + end do + end do + + !__ s2. split by space + call split_string_by_space (string_aux, length_mx, & + mxseg, nseg, str_piece, jstatus) + + return + end subroutine split_string_by_seperator end module MAPL_scan_pattern_in_file diff --git a/base/tests/utCFIO_Array.F90 b/base/tests/utCFIO_Array.F90 index 8f866b0b16ca..81537aba0533 100644 --- a/base/tests/utCFIO_Array.F90 +++ b/base/tests/utCFIO_Array.F90 @@ -6,6 +6,7 @@ Program utCFIO + use, intrinsic :: iso_fortran_env, only: REAL64 use ESMF use MAPL_BaseMod @@ -332,7 +333,7 @@ subroutine MAPL_GridGetLatLons ( grid, lons, lats ) type(ESMF_Array) :: eARRAY(2) - real(KIND=8), pointer :: R8D2(:,:) + real(KIND=REAL64), pointer :: R8D2(:,:) real, pointer :: lons2d(:,:), lats2d(:,:) real, pointer :: LONSLocal(:,:), LATSlocal(:,:) integer :: IM_WORLD, JM_WORLD, dims(3) diff --git a/components.yaml b/components.yaml index 3d65a1df9aa4..bedab2379e88 100644 --- a/components.yaml +++ b/components.yaml @@ -5,13 +5,13 @@ MAPL: ESMA_env: local: ./ESMA_env remote: ../ESMA_env.git - tag: v4.32.0 + tag: v4.34.0 develop: main ESMA_cmake: local: ./ESMA_cmake remote: ../ESMA_cmake.git - tag: v3.55.0 + tag: v3.56.0 develop: develop ecbuild: diff --git a/generic/MAPL_Generic.F90 b/generic/MAPL_Generic.F90 index 28c371198197..1df5c7986940 100644 --- a/generic/MAPL_Generic.F90 +++ b/generic/MAPL_Generic.F90 @@ -1956,10 +1956,13 @@ subroutine capture(POS, PHASE, GC, IMPORT, EXPORT, CLOCK, RC) integer :: hdr type(ESMF_Time) :: start_time, curr_time, target_time character(len=1) :: phase_ + logical :: clobber_file + call ESMF_GridCompGet(GC, NAME=comp_name, _RC) call MAPL_InternalStateGet (GC, STATE, _RC) + call MAPL_GetResource(state, clobber_file, LABEL="overwrite_checkpoint:", default = .false., _RC) call ESMF_ClockGet(clock, startTime=start_time, currTime=curr_time, _RC) call MAPL_GetResource(STATE, time_label, label='TARGET_TIME:', default='') @@ -1979,12 +1982,15 @@ subroutine capture(POS, PHASE, GC, IMPORT, EXPORT, CLOCK, RC) write(phase_, '(i1)') phase call MAPL_ESMFStateWriteToFile(import, CLOCK, trim(FILENAME)//"import_"//trim(POS)//"_runPhase"//phase_, & - FILETYPE, STATE, .false., state%grid%write_restart_by_oserver, _RC) + FILETYPE, STATE, .false., clobber=clobber_file, & + write_with_oserver=state%grid%write_restart_by_oserver, _RC) call MAPL_ESMFStateWriteToFile(export, CLOCK, trim(FILENAME)//"export_"//trim(POS)//"_runPhase"//phase_, & - FILETYPE, STATE, .false., state%grid%write_restart_by_oserver, _RC) + FILETYPE, STATE, .false., clobber=clobber_file, & + write_with_oserver=state%grid%write_restart_by_oserver, _RC) call MAPL_GetResource(STATE, hdr, default=0, LABEL="INTERNAL_HEADER:", _RC) call MAPL_ESMFStateWriteToFile(internal, CLOCK, trim(FILENAME)//"internal_"//trim(POS)//"_runPhase"//phase_, & - FILETYPE, STATE, hdr/=0, state%grid%write_restart_by_oserver, _RC) + FILETYPE, STATE, hdr/=0, clobber=clobber_file, & + write_with_oserver=state%grid%write_restart_by_oserver, _RC) end if _RETURN(_SUCCESS) end subroutine capture @@ -2247,6 +2253,7 @@ recursive subroutine MAPL_GenericFinalize ( GC, IMPORT, EXPORT, CLOCK, RC ) type(ESMF_State), pointer :: child_import_state type(ESMF_State), pointer :: child_export_state type(ESMF_State), pointer :: internal_state + logical :: clobber_file !============================================================================= ! Begin... @@ -2264,6 +2271,8 @@ recursive subroutine MAPL_GenericFinalize ( GC, IMPORT, EXPORT, CLOCK, RC ) call MAPL_InternalStateRetrieve(GC, STATE, RC=status) _VERIFY(status) + call MAPL_GetResource(state, clobber_file, LABEL="overwrite_checkpoint:", default = .false., _RC) + ! Finalize the children ! --------------------- @@ -2365,7 +2374,8 @@ recursive subroutine MAPL_GenericFinalize ( GC, IMPORT, EXPORT, CLOCK, RC ) _VERIFY(status) internal_state => state%get_internal_state() call MAPL_ESMFStateWriteToFile(internal_state,CLOCK,FILENAME, & - FILETYPE, STATE, hdr/=0, state%grid%write_restart_by_oserver, RC=status) + FILETYPE, STATE, hdr/=0, clobber=clobber_file, & + write_with_oserver=state%grid%write_restart_by_oserver, RC=status) _VERIFY(status) endif @@ -2389,7 +2399,8 @@ recursive subroutine MAPL_GenericFinalize ( GC, IMPORT, EXPORT, CLOCK, RC ) endif #endif call MAPL_ESMFStateWriteToFile(IMPORT,CLOCK,FILENAME, & - FILETYPE, STATE, .FALSE., state%grid%write_restart_by_oserver, RC=status) + FILETYPE, STATE, .FALSE., clobber=clobber_file, & + write_with_oserver=state%grid%write_restart_by_oserver, RC=status) _VERIFY(status) endif @@ -2444,7 +2455,8 @@ subroutine checkpoint_export_state(rc) endif #endif call MAPL_ESMFStateWriteToFile(EXPORT,CLOCK,FILENAME, & - FILETYPE, STATE, .FALSE., state%grid%write_restart_by_oserver, RC=status) + FILETYPE, STATE, .FALSE., clobber=clobber_file, & + write_with_oserver=state%grid%write_restart_by_oserver, RC=status) _VERIFY(status) endif _RETURN(_SUCCESS) @@ -2700,6 +2712,7 @@ subroutine MAPL_StateRecord( GC, IMPORT, EXPORT, CLOCK, RC ) integer :: hdr character(len=ESMF_MAXSTR) :: FILETYPE type(ESMF_State), pointer :: internal_state + logical :: clobber_file !============================================================================= ! Begin... @@ -2718,6 +2731,7 @@ subroutine MAPL_StateRecord( GC, IMPORT, EXPORT, CLOCK, RC ) call MAPL_InternalStateRetrieve(GC, STATE, RC=status) _VERIFY(status) + call MAPL_GetResource(state, clobber_file, LABEL="overwrite_checkpoint:", default = .false., _RC) if (.not.associated(STATE%RECORD)) then _RETURN(ESMF_SUCCESS) end if @@ -2730,7 +2744,8 @@ subroutine MAPL_StateRecord( GC, IMPORT, EXPORT, CLOCK, RC ) end if call MAPL_ESMFStateWriteToFile(IMPORT, CLOCK, & STATE%RECORD%IMP_FNAME, & - FILETYPE, STATE, .FALSE., state%grid%write_restart_by_oserver, & + FILETYPE, STATE, .FALSE., clobber=clobber_file, & + write_with_oserver=state%grid%write_restart_by_oserver, & RC=status) _VERIFY(status) end if @@ -2747,7 +2762,8 @@ subroutine MAPL_StateRecord( GC, IMPORT, EXPORT, CLOCK, RC ) internal_state => STATE%get_internal_state() call MAPL_ESMFStateWriteToFile(internal_state, CLOCK, & STATE%RECORD%INT_FNAME, & - FILETYPE, STATE, hdr/=0, state%grid%write_restart_by_oserver, & + FILETYPE, STATE, hdr/=0, clobber=clobber_file, & + write_with_oserver=state%grid%write_restart_by_oserver, & RC=status) _VERIFY(status) end if @@ -5714,7 +5730,7 @@ end subroutine MAPL_GenericStateClockAdd !============================================================================= !============================================================================= - subroutine MAPL_ESMFStateWriteToFile(STATE,CLOCK,FILENAME,FILETYPE,MPL,HDR, write_with_oserver,RC) + subroutine MAPL_ESMFStateWriteToFile(STATE,CLOCK,FILENAME,FILETYPE,MPL,HDR, write_with_oserver,clobber,RC) type(ESMF_State), intent(INOUT) :: STATE type(ESMF_Clock), intent(IN ) :: CLOCK character(len=*), intent(IN ) :: FILENAME @@ -5722,6 +5738,7 @@ subroutine MAPL_ESMFStateWriteToFile(STATE,CLOCK,FILENAME,FILETYPE,MPL,HDR, writ type(MAPL_MetaComp), intent(INOUT) :: MPL logical, intent(IN ) :: HDR logical, optional, intent(in ) :: write_with_oserver + logical, optional, intent(in ) :: clobber integer, optional, intent( OUT) :: RC character(len=ESMF_MAXSTR), parameter :: IAm="MAPL_ESMFStateWriteToFile" @@ -5742,10 +5759,12 @@ subroutine MAPL_ESMFStateWriteToFile(STATE,CLOCK,FILENAME,FILETYPE,MPL,HDR, writ integer :: attr character(len=MPI_MAX_INFO_VAL ) :: romio_cb_write logical :: nwrgt1 - logical :: empty, local_write_with_oserver + logical :: empty, local_write_with_oserver, local_clobber local_write_with_oserver=.false. if (present(write_with_oserver)) local_write_with_oserver = write_with_oserver + local_clobber = .false. + if (present(clobber)) local_clobber = clobber ! Check if state is empty. If "yes", simply return empty = MAPL_IsStateEmpty(state, _RC) @@ -5913,7 +5932,7 @@ subroutine MAPL_ESMFStateWriteToFile(STATE,CLOCK,FILENAME,FILETYPE,MPL,HDR, writ if(filetype=='pbinary' ) then - arrdes%ycomm = mpl%grid%Ycomm + call MPI_Comm_Dup(mpl%grid%Ycomm,ArrDes%Ycomm, status) call MAPL_VarWrite(UNIT=UNIT, STATE=STATE, arrdes=arrdes, rc=status) _VERIFY(status) @@ -5926,9 +5945,9 @@ subroutine MAPL_ESMFStateWriteToFile(STATE,CLOCK,FILENAME,FILETYPE,MPL,HDR, writ elseif(filetype=='pnc4') then if (local_write_with_oserver) then - call MAPL_VarWriteNCPar(filename,STATE,ArrDes,CLOCK, oClients=o_clients, _RC) + call MAPL_VarWriteNCPar(filename,STATE,ArrDes,CLOCK, clobber=local_clobber, oClients=o_clients, _RC) else - call MAPL_VarWriteNCPar(filename,STATE,ArrDes,CLOCK, _RC) + call MAPL_VarWriteNCPar(filename,STATE,ArrDes,CLOCK, clobber=local_clobber, _RC) end if elseif(UNIT/=0) then @@ -5943,6 +5962,7 @@ subroutine MAPL_ESMFStateWriteToFile(STATE,CLOCK,FILENAME,FILETYPE,MPL,HDR, writ _VERIFY(status) endif + call ArrDescrCommFree(arrdes, _RC) _RETURN(ESMF_SUCCESS) end subroutine MAPL_ESMFStateWriteToFile @@ -6041,7 +6061,7 @@ subroutine MAPL_ESMFStateReadFromFile(STATE,CLOCK,FILENAME,MPL,HDR,RC) nwrgt1 = (mpl%grid%num_readers > 1) - isNC4 = -100 + isNC4 = MAPL_FILETYPE_UNK if (on_tiles) mpl%grid%split_restart = .false. if(INDEX(FNAME,'*') == 0) then if (AmIRoot) then @@ -6091,7 +6111,7 @@ subroutine MAPL_ESMFStateReadFromFile(STATE,CLOCK,FILENAME,MPL,HDR,RC) !end if if (FileExists) then - if (isNC4 == 0) then + if (isNC4 == MAPL_FILETYPE_NC4) then filetype = 'pnc4' else if (.not.nwrgt1) then @@ -6235,7 +6255,7 @@ subroutine MAPL_ESMFStateReadFromFile(STATE,CLOCK,FILENAME,MPL,HDR,RC) if(filetype=='pbinary') then call ArrDescrSet(arrdes, offset) - arrdes%Ycomm = mpl%grid%Ycomm + call MPI_Comm_Dup(mpl%grid%Ycomm,ArrDes%Ycomm, status) call MAPL_VarRead(UNIT=UNIT, STATE=STATE, arrdes=arrdes, RC=status) _VERIFY(status) if (AmReader) then @@ -6266,6 +6286,7 @@ subroutine MAPL_ESMFStateReadFromFile(STATE,CLOCK,FILENAME,MPL,HDR,RC) call MAPL_AttributeSet(STATE, NAME="MAPL_InitStatus", VALUE=MAPL_InitialRestart, RC=status) _VERIFY(status) + call ArrDescrCommFree(arrdes, _RC) _RETURN(ESMF_SUCCESS) contains @@ -9837,7 +9858,7 @@ subroutine READIT(WHICH) im_world = mpl%grid%im_world, & jm_world = mpl%grid%jm_world) endif - arrdes%Ycomm = mpl%grid%Ycomm + call MPI_Comm_Dup(mpl%grid%Ycomm,ArrDes%Ycomm, status) if (AmReader) then call MPI_COMM_RANK(mpl%grid%readers_comm, io_rank, status) @@ -10318,11 +10339,14 @@ recursive subroutine MAPL_GenericStateSave( GC, IMPORT, EXPORT, CLOCK, RC ) type(ESMF_State), pointer :: child_import_state type(ESMF_State), pointer :: child_export_state type (ESMF_State), pointer :: internal_state + logical :: clobber_file _UNUSED_DUMMY(EXPORT) call MAPL_InternalStateRetrieve(GC, STATE, RC=status) _VERIFY(status) + call MAPL_GetResource(state, clobber_file, LABEL="overwrite_checkpoint:", default = .false., _RC) + call MAPL_GetResource( STATE, FILENAME, & LABEL="IMPORT_CHECKPOINT_FILE:", & RC=status) @@ -10406,7 +10430,7 @@ recursive subroutine MAPL_GenericStateSave( GC, IMPORT, EXPORT, CLOCK, RC ) end if call MAPL_ESMFStateWriteToFile(IMPORT, CLOCK, & STATE%initial_state%IMP_FNAME, & - CFILETYPE, STATE, .FALSE., write_with_oserver = state%grid%write_restart_by_oserver, & + CFILETYPE, STATE, .FALSE., clobber=clobber_file, write_with_oserver = state%grid%write_restart_by_oserver, & RC=status) _VERIFY(status) end if @@ -10422,7 +10446,7 @@ recursive subroutine MAPL_GenericStateSave( GC, IMPORT, EXPORT, CLOCK, RC ) internal_state => STATE%get_internal_state() call MAPL_ESMFStateWriteToFile(internal_state, CLOCK, & STATE%initial_state%INT_FNAME, & - CFILETYPE, STATE, hdr/=0, write_with_oserver = state%grid%write_restart_by_oserver, & + CFILETYPE, STATE, hdr/=0, clobber=clobber_file, write_with_oserver = state%grid%write_restart_by_oserver, & RC=status) _VERIFY(status) end if @@ -11017,10 +11041,14 @@ subroutine ArrDescrSetNCPar(ArrDes, MPL, tile, offset, num_readers, num_writers, arrdes%NX0 = mpl%grid%NX0 arrdes%tile=.true. arrdes%grid=tilegrid - call ArrDescrCreateWriterComm(arrdes,mpl%grid%comm,mpl%grid%num_writers,_RC) - call ArrDescrCreateReaderComm(arrdes,mpl%grid%comm,mpl%grid%num_readers,_RC) - arrdes%iogathercomm = mpl%grid%comm - arrdes%ioscattercomm = mpl%grid%comm + if (present(num_writers)) then + call ArrDescrCreateWriterComm(arrdes,mpl%grid%comm,mpl%grid%num_writers,_RC) + end if + if (present(num_readers)) then + call ArrDescrCreateReaderComm(arrdes,mpl%grid%comm,mpl%grid%num_readers,_RC) + end if + call MPI_Comm_Dup(mpl%grid%comm,ArrDes%iogathercomm, status) + call MPI_Comm_Dup(mpl%grid%comm,ArrDes%ioscattercomm, status) arrdes%split_restart = .false. arrdes%split_checkpoint = .false. else @@ -11045,8 +11073,12 @@ subroutine ArrDescrSetNCPar(ArrDes, MPL, tile, offset, num_readers, num_writers, arrdes%NX0 = mpl%grid%NX0 arrdes%tile=.false. arrdes%grid=MPL%GRID%ESMFGRID - call ArrDescrCreateWriterComm(arrdes,mpl%grid%comm,mpl%grid%num_writers,_RC) - call ArrDescrCreateReaderComm(arrdes,mpl%grid%comm,mpl%grid%num_readers,_RC) + if (present(num_writers)) then + call ArrDescrCreateWriterComm(arrdes,mpl%grid%comm,mpl%grid%num_writers,_RC) + end if + if (present(num_readers)) then + call ArrDescrCreateReaderComm(arrdes,mpl%grid%comm,mpl%grid%num_readers,_RC) + end if call mpi_comm_rank(arrdes%ycomm,arrdes%myrow,status) _VERIFY(status) arrdes%split_restart = mpl%grid%split_restart diff --git a/gridcomps/ExtData2G/ExtDataGridCompNG.F90 b/gridcomps/ExtData2G/ExtDataGridCompNG.F90 index 9abc4b2d1fb5..71de43fd8ec7 100644 --- a/gridcomps/ExtData2G/ExtDataGridCompNG.F90 +++ b/gridcomps/ExtData2G/ExtDataGridCompNG.F90 @@ -1360,7 +1360,7 @@ subroutine IOBundle_Add_Entry(IOBundles,item,entry_num,rc) io_bundle = ExtDataNG_IOBundle(MAPL_ExtDataLeft, entry_num, current_file, time_index, item%trans, item%fracval, item%file_template, & item%pfioCollection_id,item%iclient_collection_id,itemsL,on_tiles,_RC) call IOBundles%push_back(io_bundle) - call extdata_lgr%info('%a updated L bracket with: %a at time index %i3 ',item%name, current_file, time_index) + call extdata_lgr%info('%a updated L bracket with: %a at time index %i0 ',item%name, current_file, time_index) end if end if call item%modelGridFields%comp1%get_parameters('R',update=update,file=current_file,time_index=time_index) @@ -1370,7 +1370,7 @@ subroutine IOBundle_Add_Entry(IOBundles,item,entry_num,rc) io_bundle = ExtDataNG_IOBundle(MAPL_ExtDataRight, entry_num, current_file, time_index, item%trans, item%fracval, item%file_template, & item%pfioCollection_id,item%iclient_collection_id,itemsR,on_tiles,_RC) call IOBundles%push_back(io_bundle) - call extdata_lgr%info('%a updated R bracket with: %a at time index %i3 ',item%name,current_file, time_index) + call extdata_lgr%info('%a updated R bracket with: %a at time index %i0 ',item%name,current_file, time_index) end if end if diff --git a/gridcomps/History/MAPL_HistoryCollection.F90 b/gridcomps/History/MAPL_HistoryCollection.F90 index 336773e0400d..1a6827e7c985 100644 --- a/gridcomps/History/MAPL_HistoryCollection.F90 +++ b/gridcomps/History/MAPL_HistoryCollection.F90 @@ -47,12 +47,15 @@ module MAPL_HistoryCollectionMod integer :: acc_offset integer :: ref_date integer :: ref_time + integer :: start_date + integer :: start_time integer :: end_date integer :: end_time integer :: duration type(ESMF_Alarm) :: his_alarm ! when to write file type(ESMF_Alarm) :: seg_alarm ! segment alarm controls when to write to new file type(ESMF_Alarm) :: mon_alarm + type(ESMF_Alarm) :: start_alarm type(ESMF_Alarm) :: end_alarm integer,pointer :: expSTATE (:) integer :: unit @@ -70,6 +73,7 @@ module MAPL_HistoryCollectionMod integer :: verbose integer :: xyoffset logical :: disabled + logical :: skipWriting logical :: subVm logical :: backwards ! Adds support for clock running in reverse direction logical :: useNewFormat diff --git a/gridcomps/History/MAPL_HistoryGridComp.F90 b/gridcomps/History/MAPL_HistoryGridComp.F90 index 17deb2db1299..a6e95338cd1e 100644 --- a/gridcomps/History/MAPL_HistoryGridComp.F90 +++ b/gridcomps/History/MAPL_HistoryGridComp.F90 @@ -56,7 +56,7 @@ module MAPL_HistoryGridCompMod use MaskSamplerGeosatMod use MAPL_StringTemplate use regex_module - use MAPL_TimeUtilsMod, only: is_valid_time, is_valid_date + use MAPL_TimeUtilsMod, only: is_valid_time, is_valid_date, MAPL_UndefInt use gFTL_StringStringMap !use ESMF_CFIOMOD use MAPL_EpochSwathMod @@ -811,16 +811,19 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) label=trim(string) // 'ref_time:',_RC ) _ASSERT(is_valid_time(list(n)%ref_time),'Invalid ref_time') - call ESMF_ConfigGetAttribute ( cfg, list(n)%end_date, default=-999, & + call ESMF_ConfigGetAttribute ( cfg, list(n)%start_date, default=MAPL_UndefInt, & + label=trim(string) // 'start_date:',_RC ) + _ASSERT(is_valid_date(list(n)%start_date),'Invalid start_date') + call ESMF_ConfigGetAttribute ( cfg, list(n)%start_time, default=MAPL_UndefInt, & + label=trim(string) // 'start_time:',_RC ) + _ASSERT(is_valid_time(list(n)%start_time),'Invalid start_time') + + call ESMF_ConfigGetAttribute ( cfg, list(n)%end_date, default=MAPL_UndefInt, & label=trim(string) // 'end_date:',_RC ) - if (list(n)%end_date /= -999) then - _ASSERT(is_valid_date(list(n)%end_date),'Invalid end_date') - end if - call ESMF_ConfigGetAttribute ( cfg, list(n)%end_time, default=-999, & + _ASSERT(is_valid_date(list(n)%end_date),'Invalid end_date') + call ESMF_ConfigGetAttribute ( cfg, list(n)%end_time, default=MAPL_UndefInt, & label=trim(string) // 'end_time:',_RC ) - if (list(n)%end_time /= -999) then - _ASSERT(is_valid_time(list(n)%end_time),'Invalid end_time') - end if + _ASSERT(is_valid_time(list(n)%end_time),'Invalid end_time') call ESMF_ConfigGetAttribute ( cfg, list(n)%duration, default=list(n)%frequency, & label=trim(string) // 'duration:' ,_RC ) @@ -1349,9 +1352,41 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) intState%stampOffset(n) = Frequency ! we go to the beginning of the month end if +! End Alarm based on start_date and start_time +! ---------------------------------------- + if( list(n)%start_date.ne.MAPL_UndefInt .and. list(n)%start_time.ne.MAPL_UndefInt ) then + REF_TIME(1) = list(n)%start_date/10000 + REF_TIME(2) = mod(list(n)%start_date,10000)/100 + REF_TIME(3) = mod(list(n)%start_date,100) + REF_TIME(4) = list(n)%start_time/10000 + REF_TIME(5) = mod(list(n)%start_time,10000)/100 + REF_TIME(6) = mod(list(n)%start_time,100) + + call ESMF_TimeSet( RingTime, YY = REF_TIME(1), & + MM = REF_TIME(2), & + DD = REF_TIME(3), & + H = REF_TIME(4), & + M = REF_TIME(5), & + S = REF_TIME(6), calendar=cal, rc=rc ) + else + RingTime = CurrTime + end if + list(n)%start_alarm = ESMF_AlarmCreate( clock=clock, RingTime=RingTime, sticky=.false., _RC ) + + list(n)%skipWriting = .true. + if (RingTime == CurrTime) then + call ESMF_AlarmRingerOn(list(n)%start_alarm, _RC ) + list(n)%skipWriting = .false. + else + if (RingTime < CurrTime .NEQV. list(n)%backwards) then + list(n)%skipWriting = .false. + endif + end if + + ! End Alarm based on end_date and end_time ! ---------------------------------------- - if( list(n)%end_date.ne.-999 .and. list(n)%end_time.ne.-999 ) then + if( list(n)%end_date.ne.MAPL_UndefInt .and. list(n)%end_time.ne.MAPL_UndefInt ) then REF_TIME(1) = list(n)%end_date/10000 REF_TIME(2) = mod(list(n)%end_date,10000)/100 REF_TIME(3) = mod(list(n)%end_date,100) @@ -2527,9 +2562,13 @@ subroutine Initialize ( gc, import, dumexport, clock, rc ) else print *, ' Duration: ', list(n)%duration end if - if( list(n)%end_date.ne.-999 ) then - print *, ' End_Date: ', list(n)%end_date - print *, ' End_Time: ', list(n)%end_time + if( list(n)%start_date.ne.MAPL_UndefInt ) then + print *, ' Start_Date: ', list(n)%start_date + print *, ' Start_Time: ', list(n)%start_time + endif + if( list(n)%end_date.ne.MAPL_UndefInt ) then + print *, ' End_Date: ', list(n)%end_date + print *, ' End_Time: ', list(n)%end_time endif if (trim(list(n)%output_grid_label)/='') then print *, ' Regrid Mthd: ', regrid_method_int_to_string(list(n)%regrid_method) @@ -3427,6 +3466,14 @@ subroutine Run ( gc, import, export, clock, rc ) ! decide if we are writing based on alarms + do n=1,nlist + if (list(n)%skipWriting) then + if (ESMF_AlarmIsRinging(list(n)%start_alarm)) then + list(n)%skipWriting = .false. + endif + endif + end do + do n=1,nlist if (list(n)%disabled .or. ESMF_AlarmIsRinging(list(n)%end_alarm) ) then list(n)%disabled = .true. @@ -3454,6 +3501,8 @@ subroutine Run ( gc, import, export, clock, rc ) end if end if + if (list(n)%skipWriting) writing(n) = .false. + if (writing(n) .and. .not.IntState%average(n)) then ! R8 to R4 copy (if needed!) do m=1,list(n)%field_set%nfields @@ -5384,8 +5433,6 @@ function get_acc_offset(current_time,ref_time,rc) result(acc_offset) ! __ for each collection: find union fields, write to collection.rcx ! __ note: this subroutine is called by MPI root only ! - ! __ note: this subroutine is called by MPI root only - ! subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) use MAPL_scan_pattern_in_file use MAPL_ObsUtilMod, only : obs_platform, union_platform @@ -5417,7 +5464,7 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) integer :: nseg integer :: nseg_ub integer :: nfield, nplatform - integer :: nentry_name + integer :: nfield_name_max logical :: obs_flag integer, allocatable :: map(:) type(Logger), pointer :: lgr @@ -5444,7 +5491,6 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) length_mx = ESMF_MAXSTR2 mxseg = 100 - ! __ s1. scan get platform name + index_name_x var_name_lat/lon/time do k=1, nplf call scan_begin(unitr, 'PLATFORM.', .false.) @@ -5508,7 +5554,7 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) - ! __ s2.1 scan fields: only determine ngeoval / nentry_name = nword + ! __ s2.1 scan fields: only determine ngeoval / nfield_name_max = nword allocate (str_piece(mxseg)) rewind(unitr) do k=1, nplf @@ -5532,10 +5578,10 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) end if enddo PLFS(k)%ngeoval = ngeoval - PLFS(k)%nentry_name = nseg_ub + nseg_ub = PLFS(k)%nfield_name_mx allocate ( PLFS(k)%field_name (nseg_ub, ngeoval) ) PLFS(k)%field_name = '' - !! print*, 'k, ngeoval, nentry_name', k, ngeoval, nseg_ub + !! print*, 'k, ngeoval, nfield_name_max', k, ngeoval, nseg_ub end do @@ -5633,12 +5679,12 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) call split_string_by_space (line, length_mx, mxseg, & nplatform, str_piece, status) - !! to do: add debug - !write(6,*) 'line for obsplatforms=', trim(line) - !write(6,*) 'split string, nplatform=', nplatform - !write(6,*) 'nplf=', nplf - !write(6,*) 'str_piece=', str_piece(1:nplatform) - + call lgr%debug('%a %a', 'line for obsplatforms=', trim(line)) + call lgr%debug('%a %i6', 'split string, nplatform=', nplatform) + call lgr%debug('%a %i6', 'nplf=', nplf) + !if (mapl_am_i_root()) then + ! write(6,*) ' str_piece=', str_piece(1:nplatform) + !end if ! ! a) union the platform @@ -5654,7 +5700,10 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) end do end do deallocate(str_piece) - !! write(6,*) 'collection n=',n, 'map(:)=', map(:) + !if (mapl_am_i_root()) then + ! write(6,*) 'collection n=',n, 'map(:)=', map(:) + !end if + ! __ write common nc_index,time,lon,lat k=map(1) ! plat form # 1 @@ -5673,10 +5722,10 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) end do nfield = p1%ngeoval - nentry_name = p1%nentry_name + nfield_name_max = p1%nfield_name_mx do j=1, nfield line='' - do i=1, nentry_name + do i=1, nfield_name_max line = trim(line)//' '//trim(p1%field_name(i,j)) enddo if (j==1) then @@ -5695,7 +5744,7 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) write(unitw, '(a)') trim(adjustl(PLFS(k)%file_name_template)) do j=1, PLFS(k)%ngeoval line='' - do i=1, nentry_name + do i=1, nfield_name_max line = trim(line)//' '//trim(adjustl(PLFS(k)%field_name(i,j))) enddo write(unitw, '(a)') trim(adjustl(line)) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 index 49c2eff96b38..7cfccd4b6d40 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 @@ -48,11 +48,6 @@ module HistoryTrajectoryMod type(ESMF_TimeInterval), public :: epoch_frequency integer :: nobs_type -! character(len=ESMF_MAXSTR) :: nc_index -! character(len=ESMF_MAXSTR) :: nc_time -! character(len=ESMF_MAXSTR) :: nc_latitude -! character(len=ESMF_MAXSTR) :: nc_longitude - character(len=ESMF_MAXSTR) :: index_name_x character(len=ESMF_MAXSTR) :: var_name_time character(len=ESMF_MAXSTR) :: var_name_lat @@ -62,6 +57,8 @@ module HistoryTrajectoryMod character(len=ESMF_MAXSTR) :: var_name_lon_full character(len=ESMF_MAXSTR) :: datetime_units character(len=ESMF_MAXSTR) :: Location_index_name + logical :: use_NWP_1_file = .false. + logical :: restore_2_obs_vector = .false. integer :: epoch ! unit: second integer(kind=ESMF_KIND_I8) :: epoch_index(2) real(kind=ESMF_KIND_R8), pointer:: obsTime(:) @@ -74,6 +71,7 @@ module HistoryTrajectoryMod integer :: obsfile_Te_index logical :: active ! case: when no obs. exist logical :: level_by_level = .false. + ! ! note ! for MPI_GATHERV of 3D data in procedure :: append_file ! we have choice LEVEL_BY_LEVEL or ALL_AT_ONCE (timing in sec below for extdata) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 index 4ae3193970c9..d751474c3483 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 @@ -25,33 +25,36 @@ use, intrinsic :: iso_fortran_env, only: REAL64 use, intrinsic :: iso_fortran_env, only: INT64 implicit none - contains module procedure HistoryTrajectory_from_config use BinIOMod + use MAPL_scan_pattern_in_file 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, line + character(len=ESMF_MAXSTR) :: STR1, line, splitter character(len=ESMF_MAXSTR) :: symd, shms integer :: nline, col integer, allocatable :: ncol(:) character(len=ESMF_MAXSTR), allocatable :: word(:) + character(len=ESMF_MAXSTR), allocatable :: str_piece(:) integer :: nobs, head, jvar logical :: tend integer :: i, j, k, k2, M integer :: count, idx integer :: unitr, unitw + integer :: length_mx, mxseg, nseg type(GriddedIOitem) :: item type(Logger), pointer :: lgr traj%clock=clock if (present(GENSTATE)) traj%GENSTATE => GENSTATE + call ESMF_ClockGet ( clock, CurrTime=currTime, _RC ) call ESMF_ConfigGetAttribute(config, value=time_integer, label=trim(string)//'Epoch:', default=0, _RC) _ASSERT(time_integer /= 0, 'Epoch value in config wrong') @@ -71,6 +74,22 @@ label=trim(string) // 'var_name_lat:', _RC) call ESMF_ConfigGetAttribute(config, value=traj%var_name_time_full, default="", & label=trim(string) // 'var_name_time:', _RC) + call ESMF_ConfigGetAttribute(config, value=traj%use_NWP_1_file, default=.false., & + label=trim(string)//'use_NWP_1_file:', _RC) + call ESMF_ConfigGetAttribute(config, value=traj%restore_2_obs_vector, default=.false., & + label=trim(string)//'restore_2_obs_vector:', _RC) + if (mapl_am_I_root()) then + if (traj%use_NWP_1_file) then + write(6,105) 'WARNING: Traj sampler: use_NWP_1_file is true' + write(6,105) 'WARNING: USER needs to check if observation file is fetched correctly' + end if + if (traj%restore_2_obs_vector) then + write(6,105) 'WARNING: Traj sampler: restore_2_obs_vector is true' + end if + end if + if (.NOT. traj%use_NWP_1_file .AND. traj%restore_2_obs_vector) then + _FAIL('use_NWP_1_file=.false. and restore_2_obs_vector=.true. is not allowed') + end if call ESMF_ConfigGetAttribute(config, value=STR1, default="", & label=trim(string) // 'obs_file_begin:', _RC) @@ -147,6 +166,9 @@ ! __ s3. retrieve template and geoval, set metadata file_handle lgr => logging%get_logger('HISTORY.sampler') + length_mx = ESMF_MAXSTR + mxseg = 100 + allocate (str_piece(mxseg)) if ( nobs == 0 ) then ! ! treatment-1: @@ -192,11 +214,18 @@ ! 1-item case: file template or one-var ! 2-item : var1 , 'root' case STR1=trim(word(1)) + elseif ( count == 3 ) then + ! the Alias case + the splitField case + ! 3-item : var1 , 'root', var1_alias case + ! 3-item : var1 , 'root', 'TOTEXTTAU470;TOTEXTTAU550;TOTEXTTAU870', + ! 3-item : 'u;v' vector interpolation is not handled + STR1=trim(word(3)) else - ! 3-item : var1 , 'root', var1_alias case STR1=trim(word(3)) + call lgr%debug('%a %i8', 'WARNING: there are more than 3 field_names in platform rcx' ) end if deallocate(word, _STAT) + if ( index(trim(STR1), '-----') == 0 ) then if (head==1 .AND. trim(STR1)/='') then nobs=nobs+1 @@ -205,13 +234,29 @@ head=0 else if (trim(STR1)/='') then - jvar=jvar+1 - idx = index(STR1,";") - if (idx==0) then - traj%obs(nobs)%geoval_xname(jvar) = STR1 + splitter=';,' + call split_string_by_seperator (STR1, length_mx, splitter, mxseg, & + nseg, str_piece, status) + if (count < 3) then + ! case + ! 'var1' + ! 'var1' , 'ROOT' + ! 'u;v' , 'ROOT' + jvar=jvar+1 + if (nseg==1) then + traj%obs(nobs)%geoval_xname(jvar) = STR1 + else + traj%obs(nobs)%geoval_xname(jvar) = trim(str_piece(1)) + traj%obs(nobs)%geoval_yname(jvar) = trim(str_piece(2)) + end if else - traj%obs(nobs)%geoval_xname(jvar) = trim(STR1(1:idx-1)) - traj%obs(nobs)%geoval_yname(jvar) = trim(STR1(idx+1:)) + ! case + ! 'var1' , 'ROOT' , alias + ! 'var1' , 'ROOT' , 'TOTEXTTAU470;TOTEXTTAU550;TOTEXTTAU870,' split_field + do j=1, nseg + jvar=jvar+1 + traj%obs(nobs)%geoval_xname(jvar) = trim(str_piece(j)) + end do end if end if end if @@ -223,6 +268,15 @@ enddo end if + !!if (mapl_am_i_root()) then + !! do k=1, nobs + !! do j=1, traj%obs(k)%ngeoval + !! write(6, '(2x,a,2x,2i10,2x,a)') & + !! 'traj%obs(k)%geoval_xname(j), k, j, xname ', k, j, trim(traj%obs(k)%geoval_xname(j)) + !! end do + !! end do + !!end if + do k=1, traj%nobs_type allocate (traj%obs(k)%metadata, _STAT) @@ -335,10 +389,12 @@ call v%add_attribute('long_name', 'dateTime') call this%obs(k)%metadata%add_variable(this%var_name_time,v) - v = Variable(type=PFIO_INT32,dimensions=this%index_name_x) - call v%add_attribute('units', '1') - call v%add_attribute('long_name', 'Location index in corresponding IODA file') - call this%obs(k)%metadata%add_variable(this%location_index_name,v) + if (.NOT. this%restore_2_obs_vector) then + v = Variable(type=PFIO_INT32,dimensions=this%index_name_x) + call v%add_attribute('units', '1') + call v%add_attribute('long_name', 'Location index in corresponding IODA file') + call this%obs(k)%metadata%add_variable(this%location_index_name,v) + end if v = variable(type=PFIO_REAL64,dimensions=this%index_name_x) call v%add_attribute('units','degrees_east') @@ -436,6 +492,7 @@ iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() + !!if (mapl_am_I_root()) print*, 'create new bundle, this%items%xname= ', trim(item%xname) if (item%itemType == ItemTypeScalar) then call ESMF_FieldBundleGet(this%bundle,trim(item%xname),field=src_field,_RC) call ESMF_FieldGet(src_field,rank=rank,_RC) @@ -486,8 +543,7 @@ do k=1, this%nobs_type if (this%obs(k)%nobs_epoch > 0) then filename=trim(this%obs(k)%name)//trim(filename_suffix) - call lgr%debug('%a %a', & - "Sampling to new file : ",trim(filename)) + write(6,'(1x,a,2x,a)') "Sampling to new file :",trim(filename) call this%obs(k)%file_handle%create(trim(filename),_RC) call this%obs(k)%file_handle%write(this%obs(k)%metadata,_RC) end if @@ -540,8 +596,10 @@ real(REAL64) :: t_shift integer, allocatable :: obstype_id_full(:) integer, allocatable :: location_index_ioda_full(:) + integer, allocatable :: location_index_ioda_full_aux(:) integer, allocatable :: IA_full(:) + real(ESMF_KIND_R8), pointer :: ptAT(:) type(ESMF_routehandle) :: RH type(ESMF_Time) :: timeset(2) @@ -553,7 +611,7 @@ type(ESMF_VM) :: vm integer :: mypet, petcount, mpic - integer :: i, j, k, L, ii, jj + integer :: i, j, k, L, ii, jj, jj2 integer :: fid_s, fid_e integer(kind=ESMF_KIND_I8) :: j0, j1 integer(kind=ESMF_KIND_I8) :: jt1, jt2 @@ -612,17 +670,29 @@ trim(this%var_name_lat),trim(this%var_name_time)) L=0 - fid_s=this%obsfile_Ts_index - fid_e=this%obsfile_Te_index + if (this%use_NWP_1_file) then + ! NWP IODA 1 file case + fid_s=this%obsfile_Ts_index+1 ! index is downshifted by 1 in MAPL_ObsUtil.F90 + fid_e=fid_s + else + ! regular case for any trajectory + fid_s=this%obsfile_Ts_index ! index is downshifted by 1 in MAPL_ObsUtil.F90 + fid_e=this%obsfile_Te_index + end if call lgr%debug('%a %i10 %i10', & 'fid_s, fid_e', fid_s, fid_e) arr(1)=0 ! len_full if (mapl_am_I_root()) then + ! + ! __ s1. scan 1st time: determine len len = 0 do k=1, this%nobs_type j = max (fid_s, L) + jj2 = 0 ! obs(k) count location + this%obs(k)%count_location_until_matching_file = 0 ! init + this%obs(k)%count_location_in_matching_file = 0 ! init do while (j<=fid_e) filename = get_filename_from_template_use_index( & this%obsfile_start_time, this%obsfile_interval, & @@ -632,6 +702,12 @@ call lgr%debug('%a %a', 'exist: true filename :', trim(filename)) call get_ncfile_dimension(filename, tdim=num_times, key_time=this%index_name_x, _RC) len = len + num_times + jj2 = jj2 + num_times + if (j==this%obsfile_Ts_index) then + this%obs(k)%count_location_until_matching_file = jj2 + elseif (j==this%obsfile_Ts_index+1) then + this%obs(k)%count_location_in_matching_file = num_times + end if else call lgr%debug('%a %i10', 'non-exist: filename fid j :', j) call lgr%debug('%a %a', 'non-exist: missing filename:', trim(filename)) @@ -641,6 +717,9 @@ enddo arr(1)=len + ! + ! __ s2. scan 2nd time: read time ect. into array, + ! set location_index starting with the 1st element of the closest matched obs file if (len>0) then allocate(lons_full(len),lats_full(len),_STAT) allocate(times_R8_full(len),_STAT) @@ -652,6 +731,7 @@ ii = 0 do k=1, this%nobs_type j = max (fid_s, L) + jj2 = 0 ! obs(k) count location do while (j<=fid_e) filename = get_filename_from_template_use_index( & this%obsfile_start_time, this%obsfile_interval, & @@ -671,7 +751,14 @@ times_R8_full(len+1:len+num_times) = times_R8_full(len+1:len+num_times) + t_shift obstype_id_full(len+1:len+num_times) = k do jj = 1, num_times - location_index_ioda_full(len+jj) = jj + jj2 = jj2 + 1 + ! for each obs type use the correct starting point + ! if: use_nwp_1_file: index_ioda = [ 1, Nobs ] : restore_2_obs_vector is exact + ! else index_ioda = [ -M, 0 ] + [1, Nobs1] + [Nob1+1, Nobs2] : restore_2_obs_vector may fail + ! File1 File(center) File3 + ! why: bc we have no restriction on observation file name vs content, hence unexpected things can happen + ! use use_nwp_1_file + restore_2_obs_vector only when filename and content are systematic + location_index_ioda_full(len+jj) = jj2 - this%obs(k)%count_location_until_matching_file end do len = len + num_times end if @@ -681,7 +768,6 @@ end if end if - call ESMF_VMAllFullReduce(vm, sendData=arr, recvData=nx_sum, & count=1, reduceflag=ESMF_REDUCE_SUM, rc=rc) if (nx_sum == 0) then @@ -708,10 +794,15 @@ call MAPL_CommsBcast(vm, this%datetime_units, N=ESMF_MAXSTR, ROOT=MAPL_Root, _RC) - if (mapl_am_I_root()) then call sort_index (times_R8_full, IA_full, _RC) - location_index_ioda_full(:) = IA_full(:) + !! use index to sort togehter multiple arrays + allocate(location_index_ioda_full_aux, source=location_index_ioda_full, _STAT) + do jj = 1, nx_sum + ii = IA_full(jj) + location_index_ioda_full(jj) = location_index_ioda_full_aux(ii) + end do + deallocate(location_index_ioda_full_aux) ! NVHPC dies with NVFORTRAN-S-0155-Could not resolve generic procedure sort_multi_arrays_by_time call sort_four_arrays_by_time(lons_full, lats_full, times_R8_full, obstype_id_full, _RC) call ESMF_ClockGet(this%clock,currTime=current_time,_RC) @@ -721,7 +812,15 @@ sec = hms_2_s(this%Epoch) j1 = j0 + int(sec, kind=ESMF_KIND_I8) jx0 = real ( j0, kind=ESMF_KIND_R8) - jx1 = real ( j1, kind=ESMF_KIND_R8) + if (this%use_NWP_1_file) then + ! IODA case: + ! Upper bound time is set at 'Epoch + 1 second' to get the right index from bisect + ! + jx1 = real ( j1 + 1, kind=ESMF_KIND_R8) + else + ! normal case: + jx1 = real ( j1, kind=ESMF_KIND_R8) + end if nstart=1; nend=size(times_R8_full) call bisect( times_R8_full, jx0, jt1, n_LB=int(nstart, ESMF_KIND_I8), n_UB=int(nend, ESMF_KIND_I8), rc=rc) @@ -753,15 +852,9 @@ arr(1)=nx else !! doulbe check - ! (x1, x2] design in bisect + ! (x1, x2] design in bisect : y(n) < x <= y(n+1), n is intercept index this%epoch_index(1)= jt1 + 1 -!! ! (x1, x2] design in bisect -!! if (jt1==0) then -!! this%epoch_index(1)= 1 -!! else -!! this%epoch_index(1)= jt1 -!! endif _ASSERT(jt2<=len, 'bisect index for this%epoch_index(2) failed') if (jt2==0) then this%epoch_index(2)= 1 @@ -803,6 +896,9 @@ allocate (this%obs(k)%lats(nx2), _STAT) allocate (this%obs(k)%times_R8(nx2), _STAT) allocate (this%obs(k)%location_index_ioda(nx2), _STAT) + if (this%use_NWP_1_file) then + allocate (this%obs(k)%restore_index(nx2), _STAT) + end if enddo allocate(ix(this%nobs_type), _STAT) @@ -815,6 +911,11 @@ this%obs(k)%lats(ix(k)) = lats_full(j) this%obs(k)%times_R8(ix(k)) = times_R8_full(j) this%obs(k)%location_index_ioda(ix(k)) = location_index_ioda_full(j) + if (this%use_NWP_1_file) then + ! only this case, we have exact obs in 1_file <-> sampling match + this%obs(k)%restore_index(location_index_ioda_full(j)) = ix(k) + end if + ! !if (mod(k,10**8)==1) then ! write(6,*) 'this%obs(k)%times_R8(ix(k))', this%obs(k)%times_R8(ix(k)) !endif @@ -826,10 +927,14 @@ call lgr%debug('%a %i12 %i12 %i12', & 'epoch_index(1:2), nx', this%epoch_index(1), & this%epoch_index(2), this%nobs_epoch) - do k=1, this%nobs_type - call lgr%debug('%a %i4 %a %i12', & - 'obs(', k, ')%nobs_epoch', this%obs(k)%nobs_epoch ) - enddo + ! + ! Note: the time boundary issue can appear when we use python convention [T1, T2) but obs files donot + ! ioda file split [1/2data 15Z : 1/2data 21Z ] [ 1/2data 21Z : 1/2data 3Z] (aircraft) + ! ___x x x x x ___ ---------------------------------- o --o ---o -- o -- + ! debug: negative index (extra) at Tmin missing Tmax + ! debug shows: overcount at Tmin and missing points at Tmax + ! use_NPW_1_file=.true. solves this issue, due to special treatment at time boundaries + ! end if else allocate(this%lons(0),this%lats(0),_STAT) @@ -941,6 +1046,9 @@ real(REAL32), pointer :: p_src(:,:),p_dst(:,:), p_dst_t(:,:) ! _t: transpose real(REAL32), pointer :: p_dst_rt(:,:), p_acc_rt_3d(:,:) real(REAL32), pointer :: pt1(:), pt2(:) + real(REAL64), allocatable :: aux_R8(:) + real(REAL64), allocatable :: aux_R4(:) + integer, allocatable :: vec(:) type(ESMF_Field) :: acc_field_2d_chunk, acc_field_3d_chunk, chunk_field real(REAL32), pointer :: p_acc_chunk_3d(:,:),p_acc_chunk_2d(:) @@ -949,7 +1057,7 @@ integer :: lm integer :: rank integer :: status - integer :: j, k, ig + integer :: j, j2, k, kz, ig integer, allocatable :: ix(:) type(ESMF_VM) :: vm integer :: mypet, petcount, mpic, iroot @@ -977,14 +1085,31 @@ nx=this%obs(k)%nobs_epoch if (nx >0) then if (mapl_am_i_root()) then - call this%obs(k)%file_handle%put_var(this%var_name_time, real(this%obs(k)%times_R8), & - start=[is], count=[nx], _RC) - call this%obs(k)%file_handle%put_var(this%var_name_lon, this%obs(k)%lons, & - start=[is], count=[nx], _RC) - call this%obs(k)%file_handle%put_var(this%var_name_lat, this%obs(k)%lats, & - start=[is], count=[nx], _RC) - call this%obs(k)%file_handle%put_var(this%location_index_name, this%obs(k)%location_index_ioda, & - start=[is], count=[nx], _RC) + if (this%restore_2_obs_vector) then + ! restore back to obs vector + allocate (aux_R8(nx), vec(nx)) + vec(1:nx) = this%obs(k)%restore_index(1:nx) + aux_R8(1:nx) = this%obs(k)%times_R8(vec(1:nx)) + call this%obs(k)%file_handle%put_var(this%var_name_time, aux_R8, & + start=[is], count=[nx], _RC) + aux_R8(1:nx) = this%obs(k)%lons(vec(1:nx)) + call this%obs(k)%file_handle%put_var(this%var_name_lon, aux_R8, & + start=[is], count=[nx], _RC) + aux_R8(1:nx) = this%obs(k)%lats(vec(1:nx)) + call this%obs(k)%file_handle%put_var(this%var_name_lat, aux_R8, & + start=[is], count=[nx], _RC) + deallocate (aux_R8, vec) + else + ! default: location in time sequence + call this%obs(k)%file_handle%put_var(this%var_name_time, this%obs(k)%times_R8, & + start=[is], count=[nx], _RC) + call this%obs(k)%file_handle%put_var(this%var_name_lon, this%obs(k)%lons, & + start=[is], count=[nx], _RC) + call this%obs(k)%file_handle%put_var(this%var_name_lat, this%obs(k)%lats, & + start=[is], count=[nx], _RC) + call this%obs(k)%file_handle%put_var(this%location_index_name, this%obs(k)%location_index_ioda, & + start=[is], count=[nx], _RC) + end if end if end if enddo @@ -1053,9 +1178,12 @@ allocate ( p_dst_rt(lm, 1) ) end if + iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() + !!if (mapl_am_I_root()) print*, 'regrid: 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) @@ -1095,16 +1223,33 @@ endif enddo deallocate(ix, _STAT) + + ! rotate 2d field + if (this%restore_2_obs_vector) then + do k=1, this%nobs_type + nx = this%obs(k)%nobs_epoch + if (nx>0) then + allocate (aux_R4(nx), vec(nx)) + vec(1:nx) = this%obs(k)%restore_index(1:nx) + aux_R4(1:nx) = this%obs(k)%p2d(vec(1:nx)) + this%obs(k)%p2d(1:nx) = aux_R4(1:nx) + deallocate (aux_R4, vec) + end if + end do + end if + do k=1, this%nobs_type is = 1 nx = this%obs(k)%nobs_epoch if (nx>0) then do ig = 1, this%obs(k)%ngeoval + !! print*, 'this%obs(k)%geoval_xname(ig)= ', this%obs(k)%geoval_xname(ig) + if (trim(item%xname) == trim(this%obs(k)%geoval_xname(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 + end do endif enddo do k=1, this%nobs_type @@ -1145,7 +1290,6 @@ call ESMF_FieldDestroy(dst_field,noGarbage=.true.,_RC) call ESMF_FieldDestroy(src_field,noGarbage=.true.,_RC) - if (mapl_am_i_root()) then ! !-- pack fields to obs(k)%p3d and put_var @@ -1164,6 +1308,23 @@ this%obs(k)%p3d(ix(k),:) = p_acc_rt_3d(j,:) enddo deallocate(ix, _STAT) + + ! rotate 3d field + if (this%restore_2_obs_vector) then + do k=1, this%nobs_type + nx = this%obs(k)%nobs_epoch + if (nx>0) then + allocate (aux_R4(nx), vec(nx)) + vec(1:nx) = this%obs(k)%restore_index(1:nx) + do kz=1, lm + aux_R4(1:nx) = this%obs(k)%p3d(vec(1:nx), kz) + this%obs(k)%p3d(1:nx, kz) = aux_R4(1:nx) + end do + deallocate (aux_R4, vec) + end if + end do + end if + do k=1, this%nobs_type is = 1 nx = this%obs(k)%nobs_epoch @@ -1199,7 +1360,7 @@ integer :: x_subset(2) type(ESMF_Time) :: timeset(2) type(ESMF_Time) :: current_time - type(ESMF_TimeInterval) :: dur + type(ESMF_TimeInterval) :: dur, delT type(GriddedIOitemVectorIterator) :: iter type(GriddedIOitem), pointer :: item type(ESMF_Field) :: src_field,dst_field,acc_field @@ -1232,11 +1393,17 @@ call ESMF_ClockGet(this%clock,timeStep=dur, _RC ) timeset(1) = current_time - dur timeset(2) = current_time + if (this%use_NWP_1_file) then + ! + ! change UB to Epoch + 1 s to be inclusive for IODA + if ( ESMF_AlarmIsRinging (this%alarm) ) then + call ESMF_TimeIntervalSet(delT, s=1, _RC) + timeset(2) = current_time + delT + end if + end if call this%get_x_subset(timeset, x_subset, _RC) is=x_subset(1) ie=x_subset(2) - !! write(6,'(2x,a,4i10)') 'in regrid_accumulate is, ie=', is, ie - ! ! __ I designed a method to return from regridding if no valid points exist @@ -1351,6 +1518,9 @@ if (allocated (this%obs(k)%location_index_ioda)) then deallocate (this%obs(k)%location_index_ioda) end if + if (allocated (this%obs(k)%restore_index)) then + deallocate (this%obs(k)%restore_index) + end if if (allocated(this%obs(k)%p2d)) then deallocate (this%obs(k)%p2d) endif @@ -1427,8 +1597,7 @@ call bisect( this%obstime, rT1, index1, n_LB=lb, n_UB=ub, rc=rc) call bisect( this%obstime, rT2, index2, n_LB=lb, n_UB=ub, rc=rc) - ! (x1, x2] design in bisect - ! simple version + ! (x1, x2] design in bisect: y(n) < x <= y(n+1), n is output index x_subset(1) = index1+1 x_subset(2) = index2 diff --git a/shared/Constants/InternalConstants.F90 b/shared/Constants/InternalConstants.F90 index 3cad2914e86f..5b9fc763a405 100644 --- a/shared/Constants/InternalConstants.F90 +++ b/shared/Constants/InternalConstants.F90 @@ -23,6 +23,7 @@ module MAPL_InternalConstantsMod character(len=*), parameter :: MAPL_GRID_NAME_DEFAULT = 'UNKNOWN' character(len=*), parameter :: MAPL_GRID_FILE_NAME_DEFAULT = 'UNKNOWN' character(len=*), parameter :: MAPL_CF_COMPONENT_SEPARATOR = '.' + character(len=*), parameter :: MAPL_DESTINATIONMASK = "MAPL_DestinationMask" enum, bind(c) enumerator MAPL_TimerModeOld @@ -182,6 +183,11 @@ module MAPL_InternalConstantsMod enumerator MAPL_MASK_OUT enumerator MAPL_MASK_IN endenum + + integer, parameter :: MAPL_FILETYPE_NC4 = 0 + integer, parameter :: MAPL_FILETYPE_TXT = 1 + integer, parameter :: MAPL_FILETYPE_BIN = 2 + integer, parameter :: MAPL_FILETYPE_UNK = -1 !EOP end module MAPL_InternalConstantsMod diff --git a/shared/TimeUtils.F90 b/shared/TimeUtils.F90 index 57fe47f15301..460670adcdde 100644 --- a/shared/TimeUtils.F90 +++ b/shared/TimeUtils.F90 @@ -7,6 +7,8 @@ module MAPL_TimeUtilsMod public :: is_valid_time public :: is_valid_datetime + integer, public, parameter :: MAPL_UndefInt = -999 + contains logical function is_valid_date(date) result(is_valid) @@ -19,6 +21,11 @@ logical function is_valid_date(date) result(is_valid) integer :: year, month, day logical :: is_leap_year + if (date == MAPL_UndefInt) then + is_valid = .true. + return + end if + year = date/10000 month = mod(date,10000)/100 day = mod(date,100) @@ -66,6 +73,11 @@ logical function is_valid_time(time) result(is_valid) integer :: hours, minutes, seconds + if (time == MAPL_UndefInt) then + is_valid = .true. + return + end if + hours = time/10000 minutes = mod(time,10000)/100 seconds = mod(time,100) diff --git a/udunits2f/CMakeLists.txt b/udunits2f/CMakeLists.txt index 9ddd633fc535..149a5060dbcf 100644 --- a/udunits2f/CMakeLists.txt +++ b/udunits2f/CMakeLists.txt @@ -13,9 +13,13 @@ list (APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_LIST_DIR}") esma_add_library(${this} SRCS ${srcs} + DEPENDENCIES MAPL.shared TYPE SHARED ) +target_include_directories (${this} PUBLIC + $) + find_package(udunits REQUIRED) find_package(EXPAT REQUIRED) diff --git a/udunits2f/UDSystem.F90 b/udunits2f/UDSystem.F90 index 0fe1386978ed..b9860496f10d 100644 --- a/udunits2f/UDSystem.F90 +++ b/udunits2f/UDSystem.F90 @@ -1,10 +1,12 @@ -#include "error_handling.h" +#include "MAPL_ErrLog.h" +#include "MAPL_Exceptions.h" module ud2f_UDSystem use ud2f_CptrWrapper use ud2f_interfaces use ud2f_encoding use ud2f_status_codes + use MAPL_ExceptionHandling use iso_c_binding, only: c_ptr, c_associated, c_null_ptr, c_null_char use iso_c_binding, only: c_char, c_int, c_float, c_double, c_loc implicit none @@ -94,13 +96,13 @@ module ud2f_UDSystem contains - ! Check the status for the last udunits2 call - logical function success(utstatus) + ! Convert utstatus (UDUNITS2 status) to logical + logical function is_ut_success(utstatus) integer(ut_status) :: utstatus - success = (utstatus == UT_SUCCESS) + is_ut_success = (utstatus == UT_SUCCESS) - end function success + end function is_ut_success function construct_system(path, encoding) result(instance) type(UDsystem) :: instance @@ -112,7 +114,7 @@ function construct_system(path, encoding) result(instance) ! Read in unit system from path call read_xml(path, utsystem, status) - if(success(status)) then + if(is_ut_success(status)) then call instance%set_cptr(utsystem) if(present(encoding)) instance%encoding = encoding return @@ -135,7 +137,7 @@ function construct_unit(identifier) result(instance) cchar_identifier = cstring(identifier) utunit1 = ut_parse(SYSTEM_INSTANCE%get_cptr(), cchar_identifier, SYSTEM_INSTANCE%encoding) - if(success(ut_get_status())) then + if(is_ut_success(ut_get_status())) then call instance%set_cptr(utunit1) else ! Free memory in the case of failure @@ -149,7 +151,6 @@ function construct_converter(from_unit, to_unit) result(conv) type(UDUnit), intent(in) :: from_unit type(UDUnit), intent(in) :: to_unit type(c_ptr) :: cvconverter1 - logical :: convertible ! Must supply units that are initialized and convertible if(from_unit%is_free() .or. to_unit%is_free()) return @@ -157,7 +158,7 @@ function construct_converter(from_unit, to_unit) result(conv) cvconverter1 = ut_get_converter(from_unit%get_cptr(), to_unit%get_cptr()) - if(success(ut_get_status())) then + if(is_ut_success(ut_get_status())) then call conv%set_cptr(cvconverter1) else ! Free memory in the case of failure @@ -174,9 +175,10 @@ subroutine get_converter(conv, from, to, rc) integer(ut_status) :: status conv = get_converter_function(from, to) - _ASSERT(.not. conv%is_free(), UTF_CONVERTER_NOT_INITIALIZED) + _ASSERT_RC(.not. conv%is_free(), 'Failed to get converter function', UTF_CONVERTER_NOT_INITIALIZED) + _RETURN(_SUCCESS) + _UNUSED_DUMMY(status) - _RETURN(UT_SUCCESS) end subroutine get_converter ! Get converter object @@ -315,21 +317,20 @@ subroutine initialize(path, encoding, rc) integer, optional, intent(out) :: rc integer :: status - _RETURN_UNLESS(instance_is_uninitialized()) - ! System must be once and only once. - _ASSERT(instance_is_uninitialized(), UTF_DUPLICATE_INITIALIZATION) + ! System must be initialized once and only once. + _ASSERT_RC(instance_is_uninitialized(), 'UDSystem is initialized already.', UTF_INITIALIZATION_FAILURE) ! Disable error messages from udunits2 call disable_ut_error_message_handler() call initialize_system(SYSTEM_INSTANCE, path, encoding, rc=status) - if(status /= UT_SUCCESS) then + if(.not. is_ut_success(status)) then ! On failure, free memory call finalize() - _RETURN(UTF_INITIALIZATION_FAILURE) + _RETURN(_FAILURE) end if - _ASSERT(.not. SYSTEM_INSTANCE%is_free(), UTF_NOT_INITIALIZED) - _RETURN(UT_SUCCESS) + _ASSERT_RC(.not. SYSTEM_INSTANCE%is_free(), 'Failed to initialize UDSystem', UTF_NOT_INITIALIZED) + _RETURN(_SUCCESS) end subroutine initialize @@ -339,13 +340,13 @@ subroutine initialize_system(system, path, encoding, rc) integer(ut_encoding), optional, intent(in) :: encoding integer, optional, intent(out) :: rc integer :: status - type(c_ptr) :: utsystem ! A system can be initialized only once. - _ASSERT(system%is_free(), UTF_DUPLICATE_INITIALIZATION) - + _ASSERT_RC(system%is_free(), 'UDSystem is initialized already.', UTF_INITIALIZATION_FAILURE) system = UDSystem(path, encoding) - _RETURN(UT_SUCCESS) + _RETURN(_SUCCESS) + _UNUSED_DUMMY(status) + end subroutine initialize_system ! Is the instance of the unit system initialized? @@ -376,7 +377,6 @@ end subroutine free_ut_unit ! Free memory for converter subroutine free_cv_converter(this) class(Converter), intent(in) :: this - type(c_ptr) :: cvconverter1 if(this%is_free()) return call cv_free(this%get_cptr()) @@ -401,9 +401,9 @@ function are_convertible_udunit(unit1, unit2, rc) result(convertible) convertible = (ut_are_convertible(unit1%get_cptr(), unit2%get_cptr()) /= ZERO) status = ut_get_status() - _ASSERT(success(status), status) + _ASSERT_RC(is_ut_success(status), 'Unable to check are_convertible', status) - _RETURN(UT_SUCCESS) + _RETURN(_SUCCESS) end function are_convertible_udunit ! Check if units are convertible @@ -417,9 +417,9 @@ function are_convertible_str(from, to, rc) result(convertible) unit1 = UDUnit(from) unit2 = UDUnit(to) - convertible = are_convertible_udunit(unit1, unit2, _RC) + convertible = are_convertible(unit1, unit2, _RC) - _RETURN(UT_SUCCESS) + _RETURN(_SUCCESS) end function are_convertible_str ! Create C string from Fortran string diff --git a/udunits2f/tests/Test_UDSystem.pf b/udunits2f/tests/Test_UDSystem.pf index 14f8979a656d..4ce2c8083fa7 100644 --- a/udunits2f/tests/Test_UDSystem.pf +++ b/udunits2f/tests/Test_UDSystem.pf @@ -1,3 +1,4 @@ +#include "MAPL_ErrLog.h" module Test_UDsystem use funit @@ -21,9 +22,9 @@ contains integer(ut_status) :: status call initialize_udunits_system(rc=status) - @assertEqual(UT_SUCCESS, status, 'Failed to initialize') + @assertEqual(_SUCCESS, status, 'Failed to initialize') call get_converter(conv, KM, M, rc=status) - @assertEqual(UT_SUCCESS, status, 'Failed to get converter') + @assertEqual(_SUCCESS, status, 'Failed to get converter') @assertFalse(conv%is_free(), 'cv_converter is not set') cptr = conv%get_cptr() @assertTrue(c_associated(cptr), 'c_ptr is not associated') @@ -44,9 +45,9 @@ contains character(len=*), parameter :: TO_STRING = M call initialize_udunits_system(rc=status) - @assertEqual(UT_SUCCESS, status, 'Failed to initialize') + @assertEqual(_SUCCESS, status, 'Failed to initialize') call get_converter(conv, FROM_STRING, TO_STRING, rc=status) - @assertEqual(UT_SUCCESS, status, 'Failed to get converter') + @assertEqual(_SUCCESS, status, 'Failed to get converter') actual = conv%convert(FROM) @assertEqual(actual, EXPECTED, 'Actual does not equal expected.') call conv%free() @@ -65,9 +66,9 @@ contains character(len=*), parameter :: TO_STRING = M call initialize_udunits_system(rc=status) - @assertEqual(UT_SUCCESS, status, 'Failed to initialize') + @assertEqual(_SUCCESS, status, 'Failed to initialize') call get_converter(conv, FROM_STRING, TO_STRING, rc=status) - @assertEqual(UT_SUCCESS, status, 'Failed to get converter') + @assertEqual(_SUCCESS, status, 'Failed to get converter') actual = conv%convert(FROM) @assertEqual(actual, EXPECTED, 'Actual does not equal expected.') call conv%free() @@ -86,9 +87,9 @@ contains character(len=*), parameter :: TO_STRING = M call initialize_udunits_system(rc=status) - @assertEqual(UT_SUCCESS, status, 'Failed to initialize') + @assertEqual(_SUCCESS, status, 'Failed to initialize') call get_converter(conv, FROM_STRING, TO_STRING, rc=status) - @assertEqual(UT_SUCCESS, status, 'Failed to get converter') + @assertEqual(_SUCCESS, status, 'Failed to get converter') actual = conv%convert(FROM) @assertEqual(actual, EXPECTED, 'Actual does not equal expected.') call conv%free() @@ -107,9 +108,9 @@ contains character(len=*), parameter :: TO_STRING = M call initialize_udunits_system(rc=status) - @assertEqual(UT_SUCCESS, status, 'Failed to initialize') + @assertEqual(_SUCCESS, status, 'Failed to initialize') call get_converter(conv, FROM_STRING, TO_STRING, rc=status) - @assertEqual(UT_SUCCESS, status, 'Failed to get converter') + @assertEqual(_SUCCESS, status, 'Failed to get converter') actual = conv%convert(FROM) @assertEqual(actual, EXPECTED, 'Actual does not equal expected.') call conv%free() diff --git a/udunits2f/tests/Test_udunits2f.pf b/udunits2f/tests/Test_udunits2f.pf index ec51c125b14c..d2549e9865ac 100644 --- a/udunits2f/tests/Test_udunits2f.pf +++ b/udunits2f/tests/Test_udunits2f.pf @@ -1,3 +1,4 @@ +#include "MAPL_ErrLog.h" module Test_udunits2f use funit @@ -51,7 +52,7 @@ contains integer(ut_status) :: status call initialize_udunits_system(rc=status) - @assertEqual(UT_SUCCESS, status, 'Failed to initialize') + @assertEqual(_SUCCESS, status, 'Failed to initialize') unit1 = UDUnit(KM) @assertFalse(unit1%is_free(), 'ut_unit is not set (default encoding)') @@ -68,7 +69,7 @@ contains integer(ut_status) :: status call initialize_udunits_system(rc=status) - @assertEqual(UT_SUCCESS, status, 'Failed to initialize') + @assertEqual(_SUCCESS, status, 'Failed to initialize') unit1 = UDUnit(KM) unit2 = UDUnit(M) conv = Converter(unit1, unit2) @@ -124,7 +125,7 @@ contains logical :: convertible call initialize_udunits_system(rc=status) - @assertEqual(UT_SUCCESS, status, 'Failed to initialize') + @assertEqual(_SUCCESS, status, 'Failed to initialize') unit1 = UDUnit(KM) unit2 = UDUnit(M) convertible = are_convertible(unit1, unit2, rc=status) @@ -147,7 +148,7 @@ contains logical :: convertible call initialize_udunits_system(rc=status) - @assertEqual(UT_SUCCESS, status, 'Failed to initialize') + @assertEqual(_SUCCESS, status, 'Failed to initialize') unit1 = UDUnit(KM) unit2 = UDUnit(S) convertible = are_convertible(unit1, unit2, rc=status) @@ -155,7 +156,7 @@ contains if(.not. convertible) then @assertFalse(status == UT_BAD_ARG, 'One of the units is null.') @assertFalse(status == UT_NOT_SAME_SYSTEM, 'Units belong to different systems.') - @assertTrue(status == UT_SUCCESS, 'Units are not convertible.') + @assertTrue(status == _SUCCESS, 'Units are not convertible.') end if call unit1%free()