diff --git a/.github/workflows/extbuild.yml b/.github/workflows/extbuild.yml index 52afb7b81..185f6d85d 100644 --- a/.github/workflows/extbuild.yml +++ b/.github/workflows/extbuild.yml @@ -19,11 +19,11 @@ jobs: CPPFLAGS: "-I/usr/include -I/usr/local/include " LDFLAGS: "-L/usr/lib/x86_64-linux-gnu " # Versions of all dependencies can be updated here - these match tag names in the github repo - ESMF_VERSION: v8.5.0 - ParallelIO_VERSION: pio2_6_0 + ESMF_VERSION: v8.6.1 + ParallelIO_VERSION: pio2_6_2 steps: - id: checkout-CDEPS - uses: actions/checkout@v3 + uses: actions/checkout@v4 with: submodules: recursive - id: load-env @@ -37,7 +37,7 @@ jobs: sudo apt-get install autotools-dev autoconf - name: Cache PARALLELIO id: cache-PARALLELIO - uses: actions/cache@v3 + uses: actions/cache@v4 with: path: ${GITHUB_WORKSPACE}/pio key: ${{ runner.os }}-${{ env.ParallelIO_VERSION }}-parallelio2 diff --git a/.gitmodules b/.gitmodules index e69de29bb..45cad096c 100644 --- a/.gitmodules +++ b/.gitmodules @@ -0,0 +1,3 @@ +# This is a git-fleximod adapted .gitmodules file. Any field with a name starting in fx is a git-fleximod +# specific field. See https://github.com/ESMCI/git-fleximod for details. + diff --git a/CMakeLists.txt b/CMakeLists.txt index 28fb2a269..a4e0aa82e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -71,23 +71,6 @@ add_subdirectory(streams) add_subdirectory(dshr) if(NOT DISABLE_FoX) - if(IS_DIRECTORY "${FOX_ROOT}") - message(STATUS "FoX library is already checked out!") - message(STATUS "FoX source dir: ${FOX_ROOT}") - else() - FetchContent_Declare(fox - GIT_REPOSITORY https://github.com/ESMCI/fox.git - GIT_TAG 4.1.2.1 - SOURCE_DIR ${FOX_ROOT} - BINARY_DIR ${FOX_ROOT}/.. - ) - FetchContent_GetProperties(fox) - if(NOT fox_POPULATED) - FetchContent_Populate(fox) - message(STATUS "FoX source dir: ${fox_SOURCE_DIR}") - message(STATUS "FoX binary dir: ${fox_BINARY_DIR}") - endif() - endif() add_subdirectory(fox) target_include_directories(streams PUBLIC $ @@ -99,7 +82,7 @@ endif() target_include_directories(dshr PUBLIC $ $) -foreach(COMP datm dice dlnd docn drof dwav) +foreach(COMP datm dice dglc dlnd docn drof dwav) add_subdirectory("${COMP}") if(BLD_STANDALONE) target_include_directories(${COMP} PUBLIC $ diff --git a/cime_config/buildlib b/cime_config/buildlib index a692a4e50..536bdae66 100755 --- a/cime_config/buildlib +++ b/cime_config/buildlib @@ -73,8 +73,10 @@ def buildlib(bldroot, libroot, case): strthread = "nothreads" mpilib = case.get_value("MPILIB") compiler = case.get_value("COMPILER") - sharedpath = os.path.join(compiler, mpilib, strdebug, strthread, "nuopc") - + sharedpath = os.path.join(compiler, mpilib, strdebug, strthread) + sharedroot = case.get_value("SHAREDLIBROOT") + cdepsblddir = os.path.join(sharedroot, sharedpath, "CDEPS") + logger.info("Running cmake for CDEPS") srcpath = os.path.abspath(os.path.join(os.path.dirname(__file__), os.pardir)) cmake_flags = get_standard_cmake_args(case, os.path.join(sharedpath, "cdeps")) @@ -147,31 +149,47 @@ def buildlib(bldroot, libroot, case): else: bld_time = src_time - 1 - # if any file in src is newer than CmakeFiles in the build directory, rerun cmake + # Make sure that no other process is currently trying to build this library, done with a simple lockfile + if os.path.exists(cdepsblddir): + logger.info("{} already exists, checking for lockfile".format(cdepsblddir)) + while os.path.exists(os.path.join(cdepsblddir,"lockfile")): + logger.info("Waiting for lockfile in {}".format(cdepsblddir)) + time.sleep(10) + else: + logger.info("{} does not exist, creating lockfile".format(cdepsblddir)) + os.makedirs(cdepsblddir) + with open(os.path.join(cdepsblddir,"lockfile"),"w") as fd: + fd.write(str(os.getpid())) + + try: + # if any file in src is newer than CmakeFiles in the build directory, rerun cmake + if src_time > bld_time: + logger.info("cmake_flags {}".format(cmake_flags)) + s, o, e = run_cmd( + "cmake {} ".format(cmake_flags), from_dir=bldroot, verbose=True + ) + expect(not s, "ERROR from cmake output={}, error={}".format(o, e)) + else: + # The dwav_lib is the last file built in cdeps, wait for it to be built + dwav_lib = os.path.join(bldroot, "dwav", "libdwav.a") + time_to_wait = 600 + time_counter = 0 + while not os.path.exists(dwav_lib): + time.sleep(1) + time_counter += 1 + if time_counter > time_to_wait: + break + expect(time_counter <= time_to_wait, " Timeout waiting for {}".format(dwav_lib)) - if src_time > bld_time: - logger.info("cmake_flags {}".format(cmake_flags)) s, o, e = run_cmd( - "cmake {} ".format(cmake_flags), from_dir=bldroot, verbose=True + "make install VERBOSE=1 DESTDIR={}".format(libroot), + from_dir=bldroot, + verbose=True, ) - expect(not s, "ERROR from cmake output={}, error={}".format(o, e)) - else: - # The dwav_lib is the last file built in cdeps, wait for it to be built - dwav_lib = os.path.join(bldroot, "dwav", "libdwav.a") - time_to_wait = 300 - time_counter = 0 - while not os.path.exists(dwav_lib): - time.sleep(1) - time_counter += 1 - if time_counter > time_to_wait: - break - expect(time_counter <= time_to_wait, " Timeout waiting for {}".format(dwav_lib)) - - s, o, e = run_cmd( - "make install VERBOSE=1 DESTDIR={}".format(libroot), - from_dir=bldroot, - verbose=True, - ) + finally: + if os.path.exists(os.path.join(cdepsblddir,"lockfile")): + os.remove(os.path.join(cdepsblddir,"lockfile")) + expect(not s, "ERROR from make output={}, error={}".format(o, e)) logger.info("make output={}\nerror={}".format(o, e)) if compiler == "gnu" and case.get_value("DEBUG"): @@ -184,7 +202,7 @@ def buildlib(bldroot, libroot, case): expect(False, nextline) # Link the CDEPS component directories to the location expected by cime - for comp in ("atm", "lnd", "ice", "ocn", "rof", "wav"): + for comp in ("atm", "glc", "lnd", "ice", "ocn", "rof", "wav"): compname = case.get_value("COMP_{}".format(comp.upper())) comppath = os.path.join(case.get_value("EXEROOT"), comp, "obj") if compname == "d" + comp: diff --git a/cime_config/stream_cdeps.py b/cime_config/stream_cdeps.py index 1f702a13c..b7af2f28c 100644 --- a/cime_config/stream_cdeps.py +++ b/cime_config/stream_cdeps.py @@ -76,6 +76,7 @@ def create_stream_xml( data_list_file, user_mods_file, available_neon_data=None, + available_plumber_data=None ): """ Create the stream xml file and append the required stream input data to the input data list file @@ -187,7 +188,17 @@ def create_stream_xml( {"name": "NEON.NEON_PRECIP.$NEONSITE"}, err_msg="No stream_entry {} found".format(stream_name), ) + elif stream_name.startswith("PLUMBER2"): + self.stream_nodes = super(StreamCDEPS, self).get_child( + "stream_entry", + {"name": "PLUMBER2.$PLUMBER2SITE"}, + err_msg="No stream_entry {} found".format(stream_name), + ) elif stream_name.startswith("CLM_USRDAT."): + if 'PLUMBER2' in stream_name: + # if PLUMBER2 is in the stream name + # we want to use PLUMBER2.PLUMBER2SITE instead of CLM_USRDAT.PLUMBER2 + continue self.stream_nodes = super(StreamCDEPS, self).get_child( "stream_entry", {"name": "CLM_USRDAT.$CLM_USRDAT_NAME"}, @@ -244,6 +255,13 @@ def create_stream_xml( os.path.join(rundir, "inputdata", "atm", neon) + "\n" ) + elif available_plumber_data and stream_name.startswith("PLUMBER2"): + rundir = case.get_value("RUNDIR") + for plumber in available_plumber_data: + stream_datafiles += ( + os.path.join(rundir, "inputdata", "atm", plumber) + + "\n" + ) else: stream_datafiles = child.xml_element.text stream_datafiles = self._resolve_values( diff --git a/cime_config/stream_definition_v2.0.xsd b/cime_config/stream_definition_v2.0.xsd index de3189e35..d6ec74768 100644 --- a/cime_config/stream_definition_v2.0.xsd +++ b/cime_config/stream_definition_v2.0.xsd @@ -93,7 +93,8 @@ - + + diff --git a/cime_config/testdefs/testlist_cdeps.xml b/cime_config/testdefs/testlist_cdeps.xml deleted file mode 100644 index 03cbfb0c7..000000000 --- a/cime_config/testdefs/testlist_cdeps.xml +++ /dev/null @@ -1,150 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/datm/atm_comp_nuopc.F90 b/datm/atm_comp_nuopc.F90 index 4fae638c2..6259dbc4d 100644 --- a/datm/atm_comp_nuopc.F90 +++ b/datm/atm_comp_nuopc.F90 @@ -29,7 +29,7 @@ module cdeps_datm_comp use shr_const_mod , only : shr_const_cday use shr_sys_mod , only : shr_sys_abort use shr_cal_mod , only : shr_cal_ymd2date - use shr_log_mod , only : shr_log_setLogUnit + use shr_log_mod , only : shr_log_setLogUnit use dshr_methods_mod , only : dshr_state_diagnose, chkerr, memcheck use dshr_strdata_mod , only : shr_strdata_type, shr_strdata_init_from_config, shr_strdata_advance use dshr_strdata_mod , only : shr_strdata_get_stream_pointer, shr_strdata_setOrbs @@ -88,6 +88,12 @@ module cdeps_datm_comp use datm_datamode_simple_mod , only : datm_datamode_simple_restart_write use datm_datamode_simple_mod , only : datm_datamode_simple_restart_read + use datm_datamode_simple_mod , only : datm_datamode_simple_advertise + use datm_datamode_simple_mod , only : datm_datamode_simple_init_pointers + use datm_datamode_simple_mod , only : datm_datamode_simple_advance + use datm_datamode_simple_mod , only : datm_datamode_simple_restart_write + use datm_datamode_simple_mod , only : datm_datamode_simple_restart_read + implicit none private ! except @@ -264,7 +270,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ! Obtain flds_scalar values, mpi values, multi-instance values and ! set logunit and set shr logging to my log file - call dshr_init(gcomp, 'ATM', sdat, mpicom, my_task, inst_index, inst_suffix, & + call dshr_init(gcomp, 'ATM', mpicom, my_task, inst_index, inst_suffix, & flds_scalar_name, flds_scalar_num, flds_scalar_index_nx, flds_scalar_index_ny, & logunit, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -395,6 +401,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call datm_datamode_simple_advertise(exportState, fldsExport, flds_scalar_name, & nlfilename, my_task, vm, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + case ('SIMPLE') + call datm_datamode_simple_advertise(exportState, fldsExport, flds_scalar_name, & + nlfilename, my_task, vm, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return end select end subroutine InitializeAdvertise @@ -638,6 +648,9 @@ subroutine datm_comp_run(importState, exportState, target_ymd, target_tod, targe case('SIMPLE') call datm_datamode_simple_init_pointers(exportState, sdat, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + case('SIMPLE') + call datm_datamode_simple_init_pointers(exportState, sdat, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return end select ! Read restart if needed @@ -712,6 +725,10 @@ subroutine datm_comp_run(importState, exportState, target_ymd, target_tod, targe call datm_datamode_simple_advance(target_ymd, target_tod, target_mon, & sdat%model_calendar, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + case('SIMPLE') + call datm_datamode_simple_advance(target_ymd, target_tod, target_mon, & + sdat%model_calendar, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return end select ! Write restarts if needed @@ -739,6 +756,7 @@ subroutine datm_comp_run(importState, exportState, target_ymd, target_tod, targe case('SIMPLE') call datm_datamode_simple_restart_write(case_name, inst_suffix, target_ymd, target_tod, & logunit, my_task, sdat) + if (ChkErr(rc,__LINE__,u_FILE_u)) return end select end if diff --git a/datm/cime_config/buildnml b/datm/cime_config/buildnml index afee50344..28e341c02 100755 --- a/datm/cime_config/buildnml +++ b/datm/cime_config/buildnml @@ -141,11 +141,16 @@ def _create_namelists(case, confdir, inst_string, infile, nmlgen, data_list_path available_neon_data = None if 'CLM_USRDAT' in model_grid: config['model_grid'] = 'CLM_USRDAT' - neonsite = case.get_value("NEONSITE") - if neonsite: - config['neon'] = "True" - # download and use the listing.csv file to determine data availablity - available_neon_data = _get_neon_data_availability(case, neonsite) + if 'NEON' in clm_usrdat_name: + neonsite = case.get_value("NEONSITE") + if neonsite: + config['neon'] = "True" + # download and use the listing.csv file to determine data availablity + available_neon_data = _get_neon_data_availability(case, neonsite) + if 'PLUMBER2' in clm_usrdat_name: + plumber2site = case.get_value('PLUMBER2SITE') + if plumber2site: + config['plumber'] = "True" else: config['model_grid'] = model_grid @@ -194,6 +199,8 @@ def _create_namelists(case, confdir, inst_string, infile, nmlgen, data_list_path streamlist.append(clm_usrdat_name+"_PRECIP."+neonsite) if clm_usrdat_name == 'NEON': streamlist.append(clm_usrdat_name+".NEON_PRECIP."+neonsite) + if clm_usrdat_name == 'PLUMBER2': + streamlist.append(clm_usrdat_name+"."+plumber2site) bias_correct = nmlgen.get_value("bias_correct") if bias_correct is not None: diff --git a/datm/cime_config/config_component.xml b/datm/cime_config/config_component.xml index e0a04797d..f2698e9d4 100644 --- a/datm/cime_config/config_component.xml +++ b/datm/cime_config/config_component.xml @@ -296,6 +296,15 @@ starting year to loop data over + + integer + + 9999 + run_component_datm + env_run.xml + Start year listed in PLUMBER2 filenames for certain datm_modes. Currently only used in PLUMBER2; leave as the default value (9999) for other cases. + + integer diff --git a/datm/cime_config/namelist_definition_datm.xml b/datm/cime_config/namelist_definition_datm.xml index 491b6f6d1..3c2be3995 100644 --- a/datm/cime_config/namelist_definition_datm.xml +++ b/datm/cime_config/namelist_definition_datm.xml @@ -40,6 +40,9 @@ NEON.$NEONSITE + + PLUMBER2.$PLUMBER2SITE + CLM_USRDAT.$CLM_USRDAT_NAME diff --git a/datm/cime_config/stream_definition_datm.xml b/datm/cime_config/stream_definition_datm.xml index 310c5097e..bfe6fc4a0 100644 --- a/datm/cime_config/stream_definition_datm.xml +++ b/datm/cime_config/stream_definition_datm.xml @@ -41,6 +41,8 @@ CORE_RYF9091_JRA = JRA55 repeat year forcing, v1.3, 1990-1991 (for forcing POP and CICE) CORE_RYF0304_JRA = JRA55 repeat year forcing, v1.3, 2003-2004 (for forcing POP and CICE) ERA5 = ERA5 intra-annual year forcing + NEON = Run with forcing from NEON tower data + PLUMBER2 = Run with forcing from PLUMBER2 tower data SIMPLE = Namelist-configurable, constant datm forcing for simple experiments CPLHIST = Streams for lnd or ocn/ice forcing used for spinup @@ -341,6 +343,50 @@ single + + + + + + + + none + + + $DIN_LOC_ROOT/atm/datm7/CLM1PT_data/PLUMBER2/${PLUMBER2SITE}/CLM1PT_data/CTSM_DATM_${PLUMBER2SITE}_${DATM_YR_START_FILENAME}-${DATM_YR_END}.nc + + + ZBOT Sa_z + TBOT Sa_tbot + QBOT Sa_shum + WIND Sa_wind + PRECTmms Faxa_precn + FSDS Faxa_swdn + PSRF Sa_pbot + FLDS Faxa_lwdn + + null + + none + + null + $DATM_YR_ALIGN + $DATM_YR_START + $DATM_YR_END + 0 + + linear + + + cycle + limit + + + 1.5 + + single + + @@ -3578,7 +3624,7 @@ none - $DIN_LOC_ROOT/atm/datm7/CO2/fco2_datm_globalSSP3-7.0__simyr_2014-2501_CMIP6_c190506.nc + $DIN_LOC_ROOT/atm/datm7/CO2/fco2_datm_globalSSP3-7.0_simyr_1750-2501_CMIP6_c201101.nc CO2 Sa_co2diag @@ -4001,7 +4047,7 @@ $DIN_LOC_ROOT/share/meshes/fv0.9x1.25_141008_polemod_ESMFmesh.nc - $DIN_LOC_ROOT/atm/cam/chem/trop_mozart_aero/aero/aerodep_clm_SSP370_b.e21.BWSSP370cmip6.f09_g17.CMIP6-SSP3-7.0-WACCM.001_2014-2101_monthly_0.9x1.25_c190402.nc + $DIN_LOC_ROOT/atm/cam/chem/trop_mozart_aero/aero/aerodep_clm_SSP370_b.e21.BWSSP370cmip6.f09_g17.CMIP6-SSP3-7.0-WACCM.001_1849-2101_monthly_0.9x1.25_c201103.nc BCDEPWET Faxa_bcphiwet @@ -4036,6 +4082,7 @@ 1.5 + 30 single @@ -4575,6 +4622,7 @@ 1.5 + 30 single diff --git a/datm/cime_config/testdefs/testlist_datm.xml b/datm/cime_config/testdefs/testlist_datm.xml index ce9efb45a..34bbe3fcf 100644 --- a/datm/cime_config/testdefs/testlist_datm.xml +++ b/datm/cime_config/testdefs/testlist_datm.xml @@ -1,110 +1,115 @@ - + - + + - + - + + - + - + + - + - + + - + - + + - + - + + - + - + + - + - + + - + - + + - + - + + - + - + + - + - + + - + - - - - - - - - - + + diff --git a/datm/datm_datamode_core2_mod.F90 b/datm/datm_datamode_core2_mod.F90 index 552b35dae..41d2caae8 100644 --- a/datm/datm_datamode_core2_mod.F90 +++ b/datm/datm_datamode_core2_mod.F90 @@ -40,6 +40,8 @@ module datm_datamode_core2_mod ! export state pointers real(r8), pointer :: Sa_u(:) => null() real(r8), pointer :: Sa_v(:) => null() + real(r8), pointer :: Sa_u10m(:) => null() + real(r8), pointer :: Sa_v10m(:) => null() real(r8), pointer :: Sa_z(:) => null() real(r8), pointer :: Sa_tbot(:) => null() real(r8), pointer :: Sa_ptem(:) => null() @@ -115,6 +117,8 @@ subroutine datm_datamode_core2_advertise(exportState, fldsexport, flds_scalar_na call dshr_fldList_add(fldsExport, 'Sa_z' ) call dshr_fldList_add(fldsExport, 'Sa_u' ) call dshr_fldList_add(fldsExport, 'Sa_v' ) + call dshr_fldList_add(fldsExport, 'Sa_u10m' ) + call dshr_fldList_add(fldsExport, 'Sa_v10m' ) call dshr_fldList_add(fldsExport, 'Sa_ptem' ) call dshr_fldList_add(fldsExport, 'Sa_dens' ) call dshr_fldList_add(fldsExport, 'Sa_pslv' ) @@ -219,6 +223,10 @@ subroutine datm_datamode_core2_init_pointers(exportState, sdat, datamode, factor if (ChkErr(rc,__LINE__,u_FILE_u)) return call dshr_state_getfldptr(exportState, 'Sa_v' , fldptr1=Sa_v , rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'Sa_u10m' , fldptr1=Sa_u10m , rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'Sa_v10m' , fldptr1=Sa_v10m , rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return call dshr_state_getfldptr(exportState, 'Sa_tbot' , fldptr1=Sa_tbot , rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call dshr_state_getfldptr(exportState, 'Sa_pbot' , fldptr1=Sa_pbot , rc=rc) @@ -315,6 +323,10 @@ subroutine datm_datamode_core2_advance(datamode, target_ymd, target_tod, target_ Sa_u(n) = uprime*cos(winddFactor(n)*degtorad) - vprime*sin(winddFactor(n)*degtorad) Sa_v(n) = uprime*sin(winddFactor(n)*degtorad) + vprime*cos(winddFactor(n)*degtorad) + ! Set Sa_u10m and Sa_v10m to Sa_u and Sa_v + Sa_u10m(n) = Sa_u(n) + Sa_v10m(n) = Sa_v(n) + !--- density and pslv taken directly from input stream, set pbot --- Sa_pbot(n) = Sa_pslv(n) diff --git a/datm/datm_datamode_jra_mod.F90 b/datm/datm_datamode_jra_mod.F90 index d280bf860..9eb815150 100644 --- a/datm/datm_datamode_jra_mod.F90 +++ b/datm/datm_datamode_jra_mod.F90 @@ -25,6 +25,10 @@ module datm_datamode_jra_mod ! export state pointers real(r8), pointer :: Sa_z(:) => null() + real(r8), pointer :: Sa_u(:) => null() + real(r8), pointer :: Sa_v(:) => null() + real(r8), pointer :: Sa_u10m(:) => null() + real(r8), pointer :: Sa_v10m(:) => null() real(r8), pointer :: Sa_tbot(:) => null() real(r8), pointer :: Sa_ptem(:) => null() real(r8), pointer :: Sa_shum(:) => null() @@ -88,6 +92,8 @@ subroutine datm_datamode_jra_advertise(exportState, fldsexport, flds_scalar_name call dshr_fldList_add(fldsExport, 'Sa_z' ) call dshr_fldList_add(fldsExport, 'Sa_u' ) call dshr_fldList_add(fldsExport, 'Sa_v' ) + call dshr_fldList_add(fldsExport, 'Sa_u10m' ) + call dshr_fldList_add(fldsExport, 'Sa_v10m' ) call dshr_fldList_add(fldsExport, 'Sa_ptem' ) call dshr_fldList_add(fldsExport, 'Sa_dens' ) call dshr_fldList_add(fldsExport, 'Sa_pslv' ) @@ -174,6 +180,14 @@ subroutine datm_datamode_jra_init_pointers(exportState, sdat, rc) call shr_strdata_get_stream_pointer( sdat, 'Faxa_swdn' , strm_swdn , rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'Sa_u' , fldptr1=Sa_u , rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'Sa_v' , fldptr1=Sa_v , rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'Sa_u10m' , fldptr1=Sa_u10m , rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'Sa_v10m' , fldptr1=Sa_v10m , rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return call dshr_state_getfldptr(exportState, 'Sa_z' , fldptr1=Sa_z , rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call dshr_state_getfldptr(exportState, 'Sa_tbot' , fldptr1=Sa_tbot , rc=rc) @@ -253,6 +267,10 @@ subroutine datm_datamode_jra_advance(exportstate, target_ymd, target_tod, model_ Sa_pbot(n) = Sa_pslv(n) Sa_ptem(n) = Sa_tbot(n) + ! Set Sa_u10m and Sa_v10m to Sa_u and Sa_v + Sa_u10m(n) = Sa_u(n) + Sa_v10m(n) = Sa_v(n) + ! density computation for JRA55 forcing Sa_dens(n) = Sa_pbot(n)/(rdair*Sa_tbot(n)*(1 + 0.608*Sa_shum(n))) diff --git a/dglc/CMakeLists.txt b/dglc/CMakeLists.txt new file mode 100644 index 000000000..2c62d04b0 --- /dev/null +++ b/dglc/CMakeLists.txt @@ -0,0 +1,34 @@ +project(dglc Fortran) +set(SRCFILES glc_comp_nuopc.F90 + dglc_datamode_noevolve_mod.F90) + +foreach(FILE ${SRCFILES}) + if(EXISTS "${CASEROOT}/SourceMods/src.dglc/${FILE}") + list(REMOVE_ITEM SRCFILES ${FILE}) + list(APPEND SRCFILES "${CASEROOT}/SourceMods/src.dglc/${FILE}") + message("Using ${FILE} from ${CASEROOT}/SourceMods/src.dglc") + endif() +endforeach() + +message("DGLC srcfiles are ${SRCFILES}") + +add_library(dglc ${SRCFILES}) + +add_dependencies(dglc dshr streams) +target_include_directories (dglc PRIVATE ${ESMF_F90COMPILEPATHS}) +target_include_directories (dglc PRIVATE ${CMAKE_SOURCE_DIR}) +target_include_directories (dglc PRIVATE ${PIO_Fortran_INCLUDE_DIR}) +if(NOT DISABLE_FoX) + target_include_directories (dglc PRIVATE ${CMAKE_CURRENT_BINARY_DIR}/../fox/include) +endif() + +if(BLD_STANDALONE) + # ESMX requires mod files + foreach (SRC ${SRCFILES}) + string(REGEX REPLACE "[.]F90$" ".mod" MOD ${SRC}) + if (NOT DEFINED CIMEROOT AND MOD STREQUAL glc_comp_nuopc.mod) + set(MOD cdeps_dglc_comp.mod) + endif() + install(FILES "${CMAKE_CURRENT_BINARY_DIR}/${MOD}" DESTINATION include) + endforeach () +endif() diff --git a/dglc/cime_config/buildlib b/dglc/cime_config/buildlib new file mode 120000 index 000000000..0c5e984ac --- /dev/null +++ b/dglc/cime_config/buildlib @@ -0,0 +1 @@ +../../cime_config/buildlib_comps \ No newline at end of file diff --git a/dglc/cime_config/buildnml b/dglc/cime_config/buildnml new file mode 100755 index 000000000..b68516429 --- /dev/null +++ b/dglc/cime_config/buildnml @@ -0,0 +1,197 @@ +#!/usr/bin/env python3 + +"""Namelist creator for CIME's data glc model. +""" + +# Typically ignore this. +# pylint: disable=invalid-name + +# Disable these because this is our standard setup +# pylint: disable=wildcard-import,unused-wildcard-import,wrong-import-position + +import os, sys + +_CDEPS_CONFIG = os.path.join(os.path.dirname(os.path.abspath(__file__)), os.pardir,os.pardir,"cime_config") +_CIMEROOT = os.environ.get("CIMEROOT") +if _CIMEROOT is None: + raise SystemExit("ERROR: must set CIMEROOT environment variable") +_LIBDIR = os.path.join(_CIMEROOT, "CIME", "Tools") +sys.path.append(_LIBDIR) +sys.path.append(_CDEPS_CONFIG) + +from standard_script_setup import * +from CIME.case import Case +from CIME.nmlgen import NamelistGenerator +from CIME.utils import expect, safe_copy +from CIME.XML.files import Files +from CIME.buildnml import create_namelist_infile, parse_input, copy_inputs_to_rundir +from stream_cdeps import StreamCDEPS + +logger = logging.getLogger(__name__) + +# pylint: disable=too-many-arguments,too-many-locals,too-many-branches,too-many-statements +#################################################################################### +def _create_namelists(case, confdir, inst_string, infile, nmlgen, data_list_path): +#################################################################################### + """Write out the namelist and stream xml files for this component. + + Most arguments are the same as those for `NamelistGenerator`. The + `inst_string` argument is used as a suffix to distinguish files for + different instances. The `confdir` argument is used to specify the directory + in which output files will be placed. + """ + #---------------------------------------------------- + # Write out dglc_in and dglc.streams.xml + #---------------------------------------------------- + caseroot = case.get_value("CASEROOT") + dglc_mode = case.get_value("DGLC_MODE") + glc_grid = case.get_value("GLC_GRID") + + # Check for incompatible options. + expect(glc_grid != "null", + "DGLC_GRID cannot be null") + + # Do not allow single column mode for dglc + if case.get_value('PTS_LON'): + scol_lon = float(case.get_value('PTS_LON')) + else: + scol_lon = -999. + if case.get_value('PTS_LAT'): + scol_lat = float(case.get_value('PTS_LAT')) + else: + scol_lat = -999. + if (scol_lon > -999. or scol_lat > -999.): + expect(False, + "single column mode for DGLC is not currently allowed") + + # Log some settings. + logger.debug("DGLC mode is {}".format(dglc_mode)) + + # Initialize namelist defaults + config = {} + config['dglc_mode'] = dglc_mode + config['glc_grid'] = glc_grid + + # Initialize nmlgen + nmlgen.init_defaults(infile, config) + + # Generate dglc_in + namelist_file = os.path.join(confdir, "dglc_in") + nmlgen.write_output_file(namelist_file, data_list_path, groups=['dglc_nml']) + + # Generate dglc.streams.xml + logger.debug("dglc_mode is {}".format(dglc_mode)) + + if 'noevolve' in dglc_mode: + generate_stream_file = False + else: + generate_stream_file = True + #endif + if generate_stream_file: + # Determine streams + streamlist = nmlgen.get_streams() + if type(streamlist) == type(str()): + streamlist = [] + outfile = os.path.join(confdir, "dglc.streams"+inst_string+".xml" ) + schema_file = os.path.join(_CDEPS_CONFIG,"stream_definition_v2.0.xsd") + stream_file = os.path.join(_CDEPS_CONFIG,os.pardir, "dglc","cime_config","stream_definition_dglc.xml") + streams = StreamCDEPS(stream_file, schema_file) + streams.create_stream_xml(streamlist, case, outfile, data_list_path, + os.path.join(caseroot,'user_nl_dglc_streams'+inst_string)) + #endif + +############################################################################### +def buildnml(case, caseroot, compname): +############################################################################### + rundir = case.get_value("RUNDIR") + inst_name = compname.upper()[1:] + ninst = case.get_value("NINST_"+inst_name) + if ninst is None: + ninst = case.get_value("NINST") + confdir = os.path.join(caseroot,"Buildconf",compname + "conf") + if not os.path.isdir(confdir): + os.makedirs(confdir) + + #---------------------------------------------------- + # Construct the namelist generator + #---------------------------------------------------- + # determine directory for user modified namelist_definitions.xml + user_xml_dir = os.path.join(caseroot, "SourceMods", "src." + compname) + expect (os.path.isdir(user_xml_dir), + "user_xml_dir {} does not exist ".format(user_xml_dir)) + + # NOTE: User definition *replaces* existing definition. + files = Files(comp_interface="nuopc") + definition_file = [files.get_value("NAMELIST_DEFINITION_FILE", {"component":compname})] + + user_definition = os.path.join(user_xml_dir, "namelist_definition_"+compname+".xml") + if os.path.isfile(user_definition): + definition_file = [user_definition] + for file_ in definition_file: + expect(os.path.isfile(file_), "Namelist XML file {} not found!".format(file_)) + + # Create the namelist generator object - independent of instance + nmlgen = NamelistGenerator(case, definition_file, files=files) + + #---------------------------------------------------- + # Clear out old data. + #---------------------------------------------------- + + data_list_path = os.path.join(caseroot, "Buildconf", compname+".input_data_list") + if os.path.exists(data_list_path): + os.remove(data_list_path) + + #---------------------------------------------------- + # Loop over instances + #---------------------------------------------------- + for inst_counter in range(1, ninst+1): + # determine instance string + inst_string = "" + if ninst > 1: + inst_string = '_' + "{:04d}".format(inst_counter) + + # If multi-instance case does not have restart file, use + # single-case restart for each instance + rpointer = "rpointer." + compname + if (os.path.isfile(os.path.join(rundir,rpointer)) and + (not os.path.isfile(os.path.join(rundir,rpointer + inst_string)))): + safe_copy(os.path.join(rundir, rpointer), + os.path.join(rundir, rpointer + inst_string)) + + inst_string_label = inst_string + if not inst_string_label: + inst_string_label = "\"\"" + + # create namelist output infile using user_nl_file as input + user_nl_file = os.path.join(caseroot, "user_nl_" + compname + inst_string) + expect(os.path.isfile(user_nl_file), + "Missing required user_nl_file {} ".format(user_nl_file)) + infile = os.path.join(confdir, "namelist_infile") + create_namelist_infile(case, user_nl_file, infile) + namelist_infile = [infile] + + # create namelist and xml stream file(s) + _create_namelists(case, confdir, inst_string, namelist_infile, nmlgen, data_list_path) + + # copy namelist files and stream text files, to rundir + copy_inputs_to_rundir(caseroot, compname, confdir, rundir, inst_string) + +############################################################################### +def get_user_nl_list(case): +############################################################################### + """Returns a list of user_nl_dglc* files needed in this case + This function is called by CIME to stage the user_nl_dglc* files in the case + directory. + """ + user_nl_list = ["user_nl_dglc", "user_nl_dglc_streams"] + return user_nl_list + +############################################################################### +def _main_func(): + # Build the component namelist and required stream xml files + caseroot = parse_input(sys.argv) + with Case(caseroot) as case: + buildnml(case, caseroot, "dglc") + +if __name__ == "__main__": + _main_func() diff --git a/dglc/cime_config/config_archive.xml b/dglc/cime_config/config_archive.xml new file mode 100644 index 000000000..f61daec79 --- /dev/null +++ b/dglc/cime_config/config_archive.xml @@ -0,0 +1,11 @@ + + + r + rs1 + unset + + rpointer.glc$NINST_STRING + $CASE.dglc$NINST_STRING.r.$DATENAME.nc,$CASE.dglc$NINST_STRING.rs1.$DATENAME.bin + + + diff --git a/dglc/cime_config/config_component.xml b/dglc/cime_config/config_component.xml new file mode 100644 index 000000000..a47b9d5e9 --- /dev/null +++ b/dglc/cime_config/config_component.xml @@ -0,0 +1,92 @@ + + + + + + + + Data glc model (DGLC) + no evolve mode + + + + char + dglc + dglc + case_comp + env_case.xml + Name of land component + + + + char + noevolve + noevolve + + noevolve + + run_component_dglc + env_run.xml + + NOEVOLVE mode is used in CESM as follows. + In typical runs, CISM is not evolving; CLM computes the surface mass + balance (SMB) and sends it to CISM, but CISM’s ice sheet geometry + remains fixed over the course of the run. In these runs, CISM serves + two roles in the system: + + + + + logical + FALSE + run_component_dglc + env_run.xml + + TRUE + + Whether to include the Greenland Ice Sheet in this DGLC simulation + This should generally be set at create_newcase time (via the compset). In principle it + can be changed later, but great care is needed to change a number of other variables + to be consistent (GLC_GRID, GLC_DOMAIN_MESH and possibly others). + + + + + logical + FALSE + run_component_dglc + env_run.xml + + TRUE + + Whether to include the Antarctica Ice Sheet in this DGLC simulation + This should generally be set at create_newcase time (via the compset). In principle it + can be changed later, but great care is needed to change a number of other variables + to be consistent (GLC_GRID, GLC_DOMAIN_MESH and possibly others). + + + + + logical + TRUE,FALSE + FALSE + run_component_dglc + env_run.xml + If set to true, than dglc restarts will not be read on a continuation run. + + + + + ========================================= + DGLC naming conventions + ========================================= + + + diff --git a/dglc/cime_config/namelist_definition_dglc.xml b/dglc/cime_config/namelist_definition_dglc.xml new file mode 100644 index 000000000..8bc4e0c0e --- /dev/null +++ b/dglc/cime_config/namelist_definition_dglc.xml @@ -0,0 +1,140 @@ + + + + + + + + + + char(100) + streams + streams_file + List of streams used for each supported dglc_mode + + none + + + + + char + dglc + dglc_nml + noevolve + + Copies all fields directly from the input data streams Any required + fields not found on an input stream will be set to zero. + + + noevolve + + + + + char + streams + abs + dglc_nml + + $DIN_LOC_ROOT/glc/cism/Antarctica/ISMIP6_Antarctica_8km.init.c210908.nc:$DIN_LOC_ROOT/glc/cism/Greenland/greenland_4km_epsg3413_c171126.nc + $DIN_LOC_ROOT/glc/cism/Greenland/greenland_4km_epsg3413_c171126.nc + $DIN_LOC_ROOT/glc/cism/Antarctica/ISMIP6_Antarctica_8km.init.c210908.nc + + + colon deliminted string of inputdata files + + + + + char + streams + abs + dglc_nml + + $GLC_DOMAIN_MESH + + + colon deliminted string of file(s) specifying model mesh + for more than one ice sheets this will contain an array of mesh file names + + + + + real(20) + streams + dglc_nml + + 8000,4000 + 4000 + 8000 + + + model internal grid size(s) in m which is used to compute the internal model + model area in radians**2 using the formula (model_internal_gridsize/shr_const_rearth)**2 - + for more than one ice sheet this will contain an array of model_internal_gridsize values. + + + + + integer(20) + streams + dglc_nml + + 704,416 + 704 + 416 + 76 + + + global size(s) of nx where for more than one ice sheet this + will contain an array of nx values + + + + + integer(20) + streams + dglc_nml + + 576,704 + 576 + 704 + 141 + + + global size(s) of ny where for more than one ice sheet this + will contain an array of ny values + + + + + char + dglc + dglc_nml + + main restart file name for dglc model + + + null + + + + + logical + dglc + dglc_nml + + If set to true, than dglc restarts will not be read on a continuation run. + This capability is used, for example, in CTSM spinup runs. + + + $DGLC_SKIP_RESTART_READ + + + + diff --git a/dglc/cime_config/stream_definition_dglc.xml b/dglc/cime_config/stream_definition_dglc.xml new file mode 100644 index 000000000..80ba75c5b --- /dev/null +++ b/dglc/cime_config/stream_definition_dglc.xml @@ -0,0 +1,49 @@ + + + + + + + + + + + + + unset + + + unset + + + unset unset + + null + + bilinear + + null + 0 + 0 + 0 + 0 + + lower + + + cycle + + + 1.5 + + single + + + diff --git a/dglc/cime_config/testdefs/testlist_dglc.xml b/dglc/cime_config/testdefs/testlist_dglc.xml new file mode 100644 index 000000000..8383469c7 --- /dev/null +++ b/dglc/cime_config/testdefs/testlist_dglc.xml @@ -0,0 +1,32 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dglc/cime_config/user_nl_dglc b/dglc/cime_config/user_nl_dglc new file mode 100644 index 000000000..2cbcf3997 --- /dev/null +++ b/dglc/cime_config/user_nl_dglc @@ -0,0 +1,13 @@ +!------------------------------------------------------------------------ +! Users should ONLY USE user_nl_dglc to change namelists variables +! Users should add all user specific namelist changes below in the form of +! namelist_var = new_namelist_value +! Note that any namelist variable from shr_strdata_nml and dglc_nml can +! be modified below using the above syntax +! User preview_namelists to view (not modify) the output namelist in the +! directory $CASEROOT/CaseDocs +! To modify the contents of a stream txt file, first use preview_namelists +! to obtain the contents of the stream txt files in CaseDocs, and then +! place a copy of the modified stream txt file in $CASEROOT with the string +! user_ prepended. +!------------------------------------------------------------------------ diff --git a/dglc/cime_config/user_nl_dglc_streams b/dglc/cime_config/user_nl_dglc_streams new file mode 100644 index 000000000..856ee8944 --- /dev/null +++ b/dglc/cime_config/user_nl_dglc_streams @@ -0,0 +1,33 @@ +!------------------------------------------------------------------------ +! This file is used to modify datm.streams.xml generated in $RUNDIR +! Entries should have the form +! :<= new stream_value> +! The following are accepted values for an assume streamname of foo +! foo:meshfile = character string +! foo:datafiles = comma separated string of full pathnames (e.g. file1,file2,file3...) +! foo:datavars = comma separated string of field pairs (e.g. foo foobar,foo2 foobar2...) +! foo:taxmode = one of [cycle, extend, limit] +! foo:tintalgo = one of [lower,upper,nearest,linear,coszen] +! foo:readmode = single (only suported mode right now) +! foo:mapalgo = one of [bilinear,redist,nn,consf,consd,none] +! foo:dtlimit = real (1.5 is default) +! foo:year_first = integer +! foo:year_last = integer +! foo:year_align = integer +! foo:vectors = one of [none,u:v] +! foo:lev_dimname: = one of [null,name of level dimenion name] +! foo:offset = integer +! As an example: +! foo:year_first = 1950 +! would change the stream year_first stream_entry to 1950 for the foo stream block +! NOTE: multi-line inputs are enabled by adding a \ at the end of the line +! As an emaple: +! foo:datafiles=foo1,foo2, \ +! foo3 +! Will yield the following new entry for datafiles in stream foo +! +! foo1 +! foo2 +! foo3 +! +!------------------------------------------------------------------------ diff --git a/dglc/dglc_datamode_noevolve_mod.F90 b/dglc/dglc_datamode_noevolve_mod.F90 new file mode 100644 index 000000000..6a7caa58a --- /dev/null +++ b/dglc/dglc_datamode_noevolve_mod.F90 @@ -0,0 +1,702 @@ +module dglc_datamode_noevolve_mod + + use ESMF , only : ESMF_State, ESMF_LOGMSG_INFO, ESMF_LogWrite, ESMF_SUCCESS + use ESMF , only : ESMF_Mesh, ESMF_DistGrid, ESMF_FieldBundle, ESMF_Field + use ESMF , only : ESMF_FieldBundleCreate, ESMF_FieldCreate, ESMF_MeshLoc_Element + use ESMF , only : ESMF_FieldBundleAdd, ESMF_MeshGet, ESMF_DistGridGet, ESMF_Typekind_R8 + use ESMF , only : ESMF_GridComp, ESMF_GridCompGet + use ESMF , only : ESMF_VM, ESMF_VMAllreduce, ESMF_REDUCE_SUM + use ESMF , only : ESMF_VMGetCurrent, ESMF_VMBroadCast + use NUOPC , only : NUOPC_Advertise, NUOPC_IsConnected + use shr_kind_mod , only : r8=>shr_kind_r8, i8=>shr_kind_i8, cl=>shr_kind_cl, cs=>shr_kind_cs + use shr_sys_mod , only : shr_sys_abort + use shr_const_mod , only : SHR_CONST_RHOICE, SHR_CONST_RHOSW, SHR_CONST_REARTH, SHR_CONST_TKFRZ + use shr_const_mod , only : SHR_CONST_SPVAL + use shr_cal_mod , only : shr_cal_datetod2string + use dshr_methods_mod , only : dshr_state_getfldptr, dshr_fldbun_getfldptr, chkerr + use dshr_fldlist_mod , only : fldlist_type, dshr_fldlist_add + use dshr_strdata_mod , only : shr_strdata_type + use pio , only : file_desc_t, io_desc_t, var_desc_t, iosystem_desc_t + use pio , only : pio_openfile, pio_inq_varid, pio_inq_varndims, pio_inq_vardimid + use pio , only : pio_inq_dimlen, pio_initdecomp, pio_read_darray, pio_double + use pio , only : pio_closefile, pio_freedecomp, PIO_BCAST_ERROR, PIO_NOWRITE, PIO_CLOBBER + use pio , only : pio_createfile, pio_def_dim, pio_def_var, pio_put_att, pio_fill + use pio , only : pio_set_fill, pio_put_att, pio_enddef, pio_write_darray, PIO_GLOBAL + use pio , only : pio_seterrorhandling + + implicit none + private ! except + + public :: dglc_datamode_noevolve_advertise + public :: dglc_datamode_noevolve_init_pointers + public :: dglc_datamode_noevolve_advance + public :: dglc_datamode_noevolve_restart_write + public :: dglc_datamode_noevolve_restart_read + + logical :: initialized_noevolve = .false. + integer :: num_icesheets + real(r8) :: thk0 = 1._r8 + + ! Data structure to enable multiple ice sheets + type icesheet_ptr_t + real(r8), pointer :: ptr(:) => null() ! pointer to array + real(r8), pointer :: ptr2d(:,:) => null() ! pointer to 2d array + endtype icesheet_ptr_t + + ! Export fields + type(icesheet_ptr_t), allocatable :: Sg_area(:) + type(icesheet_ptr_t), allocatable :: Sg_topo(:) + type(icesheet_ptr_t), allocatable :: Sg_ice_covered(:) + type(icesheet_ptr_t), allocatable :: Sg_icemask(:) + type(icesheet_ptr_t), allocatable :: Sg_icemask_coupled_fluxes(:) + type(icesheet_ptr_t), allocatable :: Fgrg_rofi(:) + + ! Import fields + integer, parameter :: nlev_import = 30 + type(icesheet_ptr_t), allocatable :: Sl_tsrf(:) + type(icesheet_ptr_t), allocatable :: Flgl_qice(:) + ! type(icesheet_ptr_t), allocatable :: So_t(:) + ! type(icesheet_ptr_t), allocatable :: So_q(:) + + ! Export Field names + character(len=*), parameter :: field_out_area = 'Sg_area' + character(len=*), parameter :: field_out_topo = 'Sg_topo' + character(len=*), parameter :: field_out_ice_covered = 'Sg_ice_covered' + character(len=*), parameter :: field_out_icemask = 'Sg_icemask' + character(len=*), parameter :: field_out_icemask_coupled_fluxes = 'Sg_icemask_coupled_fluxes' + character(len=*), parameter :: field_out_rofi = 'Fgrg_rofi' + + ! Import Field names + character(len=*), parameter :: field_in_tsrf = 'Sl_tsrf' + character(len=*), parameter :: field_in_qice = 'Flgl_qice' + character(len=*), parameter :: field_in_so_t_depth = 'So_t_depth' + character(len=*), parameter :: field_in_so_s_depth = 'So_s_depth' + + character(*) , parameter :: nullstr = 'null' + character(*) , parameter :: rpfile = 'rpointer.glc' + character(*) , parameter :: u_FILE_u = & + __FILE__ + +!=============================================================================== +contains +!=============================================================================== + + subroutine dglc_datamode_noevolve_advertise(NStateExp, fldsexport, NStateImp, fldsimport, & + flds_scalar_name, rc) + + ! input/output variables + type(ESMF_State) , intent(inout) :: NStateExp(:) + type(fldlist_type), pointer :: fldsexport + type(ESMF_State) , intent(inout) :: NStateImp(:) + type(fldlist_type), pointer :: fldsimport + character(len=*) , intent(in) :: flds_scalar_name + integer , intent(out) :: rc + + ! local variables + integer :: ns + character(len=CS) :: cnum + type(fldlist_type), pointer :: fldList + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + !-------------------------------- + ! Create nested state for active ice sheets only + !-------------------------------- + + ! Set module variable for number of ice sheets + num_icesheets = size(NStateExp) + + ! Advertise export fields + call dshr_fldList_add(fldsExport, trim(flds_scalar_name)) + call dshr_fldList_add(fldsExport, field_out_area) + call dshr_fldList_add(fldsExport, field_out_ice_covered) + call dshr_fldList_add(fldsExport, field_out_topo) + call dshr_fldList_add(fldsExport, field_out_icemask) + call dshr_fldList_add(fldsExport, field_out_icemask_coupled_fluxes) + call dshr_fldList_add(fldsExport, field_out_rofi) + + do ns = 1,num_icesheets + write(cnum,'(i0)') ns + fldlist => fldsExport ! the head of the linked list + do while (associated(fldlist)) + call NUOPC_Advertise(NStateExp(ns), standardName=fldlist%stdname, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_LogWrite('(dglc_comp_advertise): Fr_glc'//trim(cnum)//"_"//trim(fldList%stdname), ESMF_LOGMSG_INFO) + fldList => fldList%next + end do + enddo + + ! Advertise import fields if appropriate + call dshr_fldList_add(fldsImport, trim(flds_scalar_name)) + call dshr_fldList_add(fldsImport, field_in_tsrf) + call dshr_fldList_add(fldsImport, field_in_qice) + ! TODO: Add namelist for this + call dshr_fldList_add(fldsImport, field_in_so_t_depth, ungridded_lbound=1, ungridded_ubound=nlev_import) + call dshr_fldList_add(fldsImport, field_in_so_s_depth, ungridded_lbound=1, ungridded_ubound=nlev_import) + + do ns = 1,num_icesheets + write(cnum,'(i0)') ns + fldlist => fldsImport ! the head of the linked list + do while (associated(fldlist)) + call NUOPC_Advertise(NStateImp(ns), standardName=fldlist%stdname, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_LogWrite('(dglc_comp_advertise): To_glc'//trim(cnum)//"_"//trim(fldList%stdname), ESMF_LOGMSG_INFO) + fldList => fldList%next + end do + enddo + + end subroutine dglc_datamode_noevolve_advertise + + !=============================================================================== + subroutine dglc_datamode_noevolve_init_pointers(NStateExp, NstateImp, rc) + + ! input/output variables + type(ESMF_State) , intent(inout) :: NStateExp(:) + type(ESMF_State) , intent(inout) :: NStateImp(:) + integer , intent(out) :: rc + + ! local variables + integer :: ns + character(len=*), parameter :: subname='(dglc_init_pointers): ' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! So this is tricky since you need pointers to fields in the nested state + ! So this will have to be done below in a loop + + ! initialize pointers to export fields + allocate(Sg_area(num_icesheets)) + allocate(Sg_topo(num_icesheets)) + allocate(Sg_ice_covered(num_icesheets)) + allocate(Sg_icemask(num_icesheets)) + allocate(Sg_icemask_coupled_fluxes(num_icesheets)) + allocate(Fgrg_rofi(num_icesheets)) + + do ns = 1,num_icesheets + call dshr_state_getfldptr(NStateExp(ns), field_out_area, & + fldptr1=Sg_area(ns)%ptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(NStateExp(ns), field_out_topo, & + fldptr1=Sg_topo(ns)%ptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(NStateExp(ns), field_out_ice_covered, & + fldptr1=Sg_ice_covered(ns)%ptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(NStateExp(ns), field_out_icemask, & + fldptr1=Sg_icemask(ns)%ptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(NStateExp(ns), field_out_icemask_coupled_fluxes, & + fldptr1=Sg_icemask_coupled_fluxes(ns)%ptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(NStateExp(ns), field_out_rofi, & + fldptr1=Fgrg_rofi(ns)%ptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + Fgrg_rofi(ns)%ptr(:) = 0._r8 + end do + + ! initialize pointers to import fields if appropriate + allocate(Sl_tsrf(num_icesheets)) + allocate(Flgl_qice(num_icesheets)) + + do ns = 1,num_icesheets + if (.not. NUOPC_IsConnected(NStateImp(ns), fieldName=field_in_tsrf)) then + ! NOTE: the field is connected ONLY if the MED->GLC entry is in the nuopc.runconfig file + ! This restriction occurs even if the field was advertised + call shr_sys_abort(trim(subname)//": MED->GLC must appear in run sequence") + end if + call dshr_state_getfldptr(NStateImp(ns), field_in_tsrf, fldptr1=Sl_tsrf(ns)%ptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(NStateImp(ns), field_in_qice, fldptr1=Flgl_qice(ns)%ptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + Sl_tsrf(ns)%ptr(:) = SHR_CONST_TKFRZ + Flgl_qice(ns)%ptr(:) = 0._r8 + end do + + end subroutine dglc_datamode_noevolve_init_pointers + + !=============================================================================== + subroutine dglc_datamode_noevolve_advance(gcomp, pio_subsystem, io_type, io_format, & + logunit, model_meshes, model_internal_gridsize, model_datafiles, rc) + + ! Assume that the model mesh is the same as the input data mesh + + ! input/output variables + type(ESMF_GridComp) :: gcomp + type(iosystem_desc_t) , pointer :: pio_subsystem ! pio info + integer , intent(in) :: io_type ! pio info + integer , intent(in) :: io_format ! pio info + integer , intent(in) :: logunit ! For writing logs + type(ESMF_Mesh) , intent(in) :: model_meshes(:) ! ice sheets meshes + real(r8) , intent(in) :: model_internal_gridsize(:) ! internal model gridsizes (m) + character(len=*) , intent(in) :: model_datafiles(:) ! input file names + integer , intent(out) :: rc + + ! local variables + type(ESMF_FieldBundle) :: fldbun_noevolve + type(ESMF_DistGrid) :: distgrid + type(ESMF_Field) :: field_noevolve + type(ESMF_VM) :: vm + type(file_desc_t) :: pioid + type(io_desc_t) :: pio_iodesc + integer :: ns ! ice sheet index + integer :: ng ! grid cell index + integer :: lsize ! local size + integer, pointer :: gindex(:) ! domain decomposition of data + integer :: ndims ! number of dims + integer, allocatable :: dimid(:) + type(var_desc_t) :: varid + integer :: rcode + integer :: nxg, nyg + real(r8), pointer :: topog(:) + real(r8), pointer :: thck(:) + logical :: exists + real(r8) :: rhoi ! density of ice ~ kg/m^3 + real(r8) :: rhoo ! density of sea water ~ kg/m^3 + real(r8) :: eus ! eustatic sea level + real(r8) :: lsrf ! lower surface elevation (m) on ice grid + real(r8) :: usrf ! upper surface elevation (m) on ice grid + real(r8) :: loc_pos_smb(1), Tot_pos_smb(1) ! Sum of positive smb values on each ice sheet for hole-filling + real(r8) :: loc_neg_smb(1), Tot_neg_smb(1) ! Sum of negative smb values on each ice sheet for hole-filling + real(r8) :: rat ! Ratio of hole-filling flux to apply + + character(len=*), parameter :: subname='(dglc_datamode_noevolve_advance): ' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + if (.not. initialized_noevolve) then + + ! Loop over ice sheets + do ns = 1,num_icesheets + + ! Determine lsize and gindex + call ESMF_MeshGet(model_meshes(ns), elementdistGrid=distGrid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_DistGridGet(distGrid, localDe=0, elementCount=lsize, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + allocate(gindex(lsize)) + call ESMF_DistGridGet(distGrid, localDe=0, seqIndexList=gindex, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Determine "glc_area" ; + ! Sg_areas is in radians + ! SHR_CONST_REARTH is the radius of earth in m + ! model_internal_gridsize is the internal model gridsize in m + do ng = 1,lsize + Sg_area(ns)%ptr(ng) = (model_internal_gridsize(ns)/SHR_CONST_REARTH)**2 + end do + + ! Create module level field bundle + fldbun_noevolve = ESMF_FieldBundleCreate(rc=rc) ! input field bundle + + ! "ice thickness" ; + field_noevolve = ESMF_FieldCreate(model_meshes(ns), ESMF_TYPEKIND_R8, & + name='thk', meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldBundleAdd(fldbun_noevolve, (/field_noevolve/), rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! "bed topography" ; + field_noevolve = ESMF_FieldCreate(model_meshes(ns), ESMF_TYPEKIND_R8, & + name='topg', meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldBundleAdd(fldbun_noevolve, (/field_noevolve/), rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Create pioid and pio_iodesc at the module level + inquire(file=trim(model_datafiles(ns)), exist=exists) + if (.not.exists) then + write(6,'(a)')' ERROR: model input file '//trim(model_datafiles(ns))//' does not exist' + call shr_sys_abort() + end if + rcode = pio_openfile(pio_subsystem, pioid, io_type, trim(model_datafiles(ns)), pio_nowrite) + call pio_seterrorhandling(pioid, PIO_BCAST_ERROR) + rcode = pio_inq_varid(pioid, 'thk', varid) + rcode = pio_inq_varndims(pioid, varid, ndims) + allocate(dimid(ndims)) + rcode = pio_inq_vardimid(pioid, varid, dimid(1:ndims)) + rcode = pio_inq_dimlen(pioid, dimid(1), nxg) + rcode = pio_inq_dimlen(pioid, dimid(2), nyg) + call pio_initdecomp(pio_subsystem, pio_double, (/nxg,nyg/), gindex, pio_iodesc) + deallocate(dimid) + + ! Read in the data into the appropriate field bundle pointers + ! Note that Sg_ice_covered(ns)%ptr points into the data for + ! the Sg_ice_covered field in NStateExp(ns) + ! Note that Sg_topo(ns)%ptr points into the data for + ! the Sg_topon NStateExp(ns) + ! Note that topog is bedrock topography + + call dshr_fldbun_getFldPtr(fldbun_noevolve, 'topg', topog, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + rcode = pio_inq_varid(pioid, 'topg', varid) + call pio_read_darray(pioid, varid, pio_iodesc, topog, rcode) + + call dshr_fldbun_getFldPtr(fldbun_noevolve, 'thk', thck, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + rcode = pio_inq_varid(pioid, 'thk', varid) + call pio_read_darray(pioid, varid, pio_iodesc, thck, rcode) + + rhoi = SHR_CONST_RHOICE ! 0.917e3 + rhoo = SHR_CONST_RHOSW ! 1.026e3 + eus = 0 + do ng = 1,lsize + if (topog(ng) - eus < (-rhoi/rhoo) * thck(ng)) then + lsrf = (-rhoi/rhoo) * thck(ng) + else + lsrf = topog(ng) + end if + usrf = max(0.d0, thck(ng) + lsrf) + + ! The export field 'ice_mask_coupled_fluxes' determines who is handling the + ! runoff associated with the surface mass balance + ! If its 0 -then ctsm needs to handle it. + ! Since we want dglc to handle it no evolve mode - then + ! ice_mask_coupled_fluxes to be identical to the mask + + if (is_in_active_grid(usrf)) then + Sg_icemask(ns)%ptr(ng) = 1.d0 + Sg_icemask_coupled_fluxes(ns)%ptr(ng) = 1.d0 + if (is_ice_covered(thck(ng))) then + Sg_ice_covered(ns)%ptr(ng) = 1.d0 + else + Sg_ice_covered(ns)%ptr(ng) = 0.d0 + end if + ! Note that we use the same method for computing topo whether this point is + ! ice-covered or ice-free. This is in contrast to the method for computing + ! ice-free topo in glint_upscaling_gcm. + Sg_topo(ns)%ptr(ng) = thk0 * usrf + else + ! Note that this logic implies that if (in theory) we had an ice-covered + ! point outside the "active grid", it will get classified as ice-free for + ! these purposes. This mimics the logic currently in glint_upscaling_gcm. + Sg_icemask(ns)%ptr(ng) = 0.d0 + Sg_icemask_coupled_fluxes(ns)%ptr(ng) = 0.d0 + Sg_ice_covered(ns)%ptr(ng) = 0.d0 + Sg_topo(ns)%ptr(ng) = 0.d0 + end if + end do + + call pio_closefile(pioid) + call pio_freedecomp(pio_subsystem, pio_iodesc) + + end do ! end loop over ice sheets + + end if + + if (initialized_noevolve) then + + ! Compute Fgrg_rofi + do ns = 1,num_icesheets + + ! Get number of grid cells per ice sheet + lsize = size(Fgrg_rofi(ns)%ptr) + + ! reset output variables to zero + Fgrg_rofi(ns)%ptr(:) = 0.d0 + loc_pos_smb(1) = 0.d0 + Tot_pos_smb(1) = 0.d0 + loc_neg_smb(1) = 0.d0 + Tot_neg_smb(1) = 0.d0 + rat = 0.d0 + + ! For No Evolve to reduce negative ice fluxes from DGLC, we will + ! Calculate the total positive and total negative fluxes on each + ! processor first (local totals). + do ng = 1,lsize + if (Sg_icemask_coupled_fluxes(ns)%ptr(ng) > 0.d0) then + if(Flgl_qice(ns)%ptr(ng) > 0.d0) then + loc_pos_smb(1) = loc_pos_smb(1)+Flgl_qice(ns)%ptr(ng)*Sg_area(ns)%ptr(ng) + end if + ! Ignore places that are exactly 0.d0 + if(Flgl_qice(ns)%ptr(ng) < 0.d0) then + loc_neg_smb(1) = loc_neg_smb(1)+Flgl_qice(ns)%ptr(ng)*Sg_area(ns)%ptr(ng) + end if + end if + end do + ! Now do two global sums to get the ice sheet total positive + ! and negative ice fluxes + call ESMF_GridCompGet(gcomp, vm=vm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllreduce(vm, senddata=loc_pos_smb, recvdata=Tot_pos_smb, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllreduce(vm, senddata=loc_neg_smb, recvdata=Tot_neg_smb, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! If there's more positive than negative, then set all + ! negative to zero and destribute the negative flux amount + ! across the positive values, scaled by the size of the + ! positive value. This section also applies to any chunks + ! where there is no negative smb. In that case the ice + ! runoff is exactly equal to the input smb. + if(abs(Tot_pos_smb(1)) >= abs(Tot_neg_smb(1))) then + do ng = 1,lsize + if (Sg_icemask_coupled_fluxes(ns)%ptr(ng) > 0.d0) then + if(Flgl_qice(ns)%ptr(ng) > 0.d0) then + rat = Flgl_qice(ns)%ptr(ng)/Tot_pos_smb(1) + Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_neg_smb(1) + else + Fgrg_rofi(ns)%ptr(ng) = 0.d0 + end if + else + Fgrg_rofi(ns)%ptr(ng) = 0.d0 + end if + end do + else + ! If there's more negative than positive, set the positive to zero + ! and distribute the positive amount to the negative spaces to + ! reduce their negativity a bit. This shouldn't happen often. + ! This section of code also applies if Tot_pos_smb is zero. + do ng = 1,lsize + if (Sg_icemask_coupled_fluxes(ns)%ptr(ng) > 0.d0) then + if(Flgl_qice(ns)%ptr(ng) < 0.d0) then + rat = Flgl_qice(ns)%ptr(ng)/Tot_neg_smb(1) + Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_pos_smb(1) + else + Fgrg_rofi(ns)%ptr(ng) = 0.d0 + end if + else + Fgrg_rofi(ns)%ptr(ng) = 0.d0 + end if + end do + + end if ! More neg or pos smb + + end do ! Each ice sheet + end if ! if initialized + + ! Set initialized flag + initialized_noevolve = .true. + + end subroutine dglc_datamode_noevolve_advance + + !=============================================================================== + logical function is_in_active_grid(usrf) + ! Return true if the given point is inside the "active grid". The active grid includes + ! any point that can receive a positive surface mass balance, which includes any + ! point classified as land or ice sheet. + + real(r8), intent(in) :: usrf ! surface elevation (m) + + if (thk0 * usrf > 0.d0) then + ! points not at sea level are assumed to be land or ice sheet + is_in_active_grid = .true. + else + is_in_active_grid = .false. + end if + end function is_in_active_grid + + !=============================================================================== + logical function is_ice_covered(thck) + ! Return true if the given point is ice-covered + + real(r8), intent(in) :: thck ! ice thickness (m) + real(r8), parameter :: min_thck = 0.d0 + + if (thk0 * thck > min_thck) then + is_ice_covered = .true. + else + is_ice_covered = .false. + end if + end function is_ice_covered + + !=============================================================================== + subroutine dglc_datamode_noevolve_restart_write(model_meshes, case_name, & + inst_suffix, ymd, tod, logunit, my_task, main_task, & + pio_subsystem, io_type, nx_global, ny_global, rc) + + ! input/output variables + type(ESMF_Mesh) , intent(in) :: model_meshes(:) ! ice sheets meshes + character(len=*) , intent(in) :: case_name + character(len=*) , intent(in) :: inst_suffix + integer , intent(in) :: ymd ! model date + integer , intent(in) :: tod ! model sec into model date + integer , intent(in) :: logunit + integer , intent(in) :: my_task + integer , intent(in) :: main_task + type(iosystem_desc_t) , pointer :: pio_subsystem ! pio info + integer , intent(in) :: io_type ! pio info + integer , intent(in) :: nx_global(:) + integer , intent(in) :: ny_global(:) + integer , intent(out) :: rc + + ! local variables + type(ESMF_DistGrid) :: distgrid + integer :: ns + character(len=CS) :: cnum + integer :: lsize + integer, pointer :: gindex(:) ! domain decomposition of data + integer :: nu + character(len=CL) :: rest_file_model + character(len=CS) :: date_str + type(file_desc_t) :: pioid + integer :: dimid2(2) + type(var_desc_t), allocatable :: varid(:) + type(io_desc_t), allocatable :: pio_iodesc(:) + integer :: oldmode + integer :: rcode + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + call shr_cal_datetod2string(date_str, ymd, tod) + write(rest_file_model ,"(7a)") trim(case_name),'.','dglc',trim(inst_suffix),'.r.',trim(date_str),'.nc' + ! write restart info to rpointer file + if (my_task == main_task) then + open(newunit=nu, file=trim(rpfile)//trim(inst_suffix), form='formatted') + write(nu,'(a)') rest_file_model + close(nu) + write(logunit,'(a,2x,i0,2x,i0)')' writing with no streams '//trim(rest_file_model), ymd, tod + endif + + ! write data model restart data + rcode = pio_createfile(pio_subsystem, pioid, io_type, trim(rest_file_model), pio_clobber) + allocate(varid(num_icesheets)) + do ns = 1,num_icesheets + ! Need to explicitly write restart since noevolve mode does not read a stream + write(cnum,'(i0)') ns + + rcode = pio_def_dim(pioid, '_nx'//trim(cnum), nx_global(ns), dimid2(1)) + rcode = pio_def_dim(pioid, '_ny'//trim(cnum), ny_global(ns), dimid2(2)) + rcode = pio_def_var(pioid, 'flgl_rofi'//cnum, PIO_DOUBLE, (/dimid2/), varid(ns)) + rcode = pio_put_att(pioid, varid(ns), "_FillValue", shr_const_spval) + rcode = pio_set_fill(pioid, PIO_FILL, oldmode) + rcode = pio_put_att(pioid, pio_global, "version", "nuopc_data_models_v0") + enddo + rcode = pio_enddef(pioid) + allocate(pio_iodesc(num_icesheets)) + do ns = 1,num_icesheets + + ! Determine gindex for this ice sheet + call ESMF_MeshGet(model_meshes(ns), elementdistGrid=distGrid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_DistGridGet(distGrid, localDe=0, elementCount=lsize, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + allocate(gindex(lsize)) + call ESMF_DistGridGet(distGrid, localDe=0, seqIndexList=gindex, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call pio_initdecomp(pio_subsystem, pio_double, (/nx_global(ns),ny_global(ns)/), gindex, pio_iodesc(ns)) + call pio_write_darray(pioid, varid(ns), pio_iodesc(ns), Fgrg_rofi(ns)%ptr, rcode, fillval=shr_const_spval) + + ! Deallocate gindex + deallocate (gindex) + end do + call pio_closefile(pioid) + do ns = 1,num_icesheets + call pio_freedecomp(pio_subsystem, pio_iodesc(ns)) + enddo + + end subroutine dglc_datamode_noevolve_restart_write + + !=============================================================================== + subroutine dglc_datamode_noevolve_restart_read(model_meshes, restfilem, & + inst_suffix, logunit, my_task, main_task, mpicom, & + pio_subsystem, io_type, nx_global, ny_global, rc) + + ! input/output arguments + type(ESMF_Mesh) , intent(in) :: model_meshes(:) ! ice sheets meshes + character(len=*) , intent(inout) :: restfilem + character(len=*) , intent(in) :: inst_suffix + integer , intent(in) :: logunit + integer , intent(in) :: my_task + integer , intent(in) :: main_task + integer , intent(in) :: mpicom + type(iosystem_desc_t) , pointer :: pio_subsystem ! pio info + integer , intent(in) :: io_type ! pio info + integer , intent(in) :: nx_global(:) + integer , intent(in) :: ny_global(:) + integer , intent(out) :: rc + + ! local variables + type(ESMF_DistGrid) :: distgrid + integer :: ns + character(len=CS) :: cnum + integer :: lsize + integer, pointer :: gindex(:) ! domain decomposition of data + type(ESMF_VM) :: vm + integer :: nu + logical :: exists ! file existance + type(file_desc_t) :: pioid + type(var_desc_t) :: varid + type(io_desc_t) :: pio_iodesc + integer :: rcode + integer :: tmp(1) + character(*), parameter :: subName = "(dglc_datamode_noevolve_restart_read) " + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + ! Determine restart file + + + if (trim(restfilem) == trim(nullstr)) then + exists = .false. + call ESMF_VMGetCurrent(vm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (my_task == main_task) then + write(logunit,'(a)') trim(subname)//' restart filename from rpointer' + inquire(file=trim(rpfile)//trim(inst_suffix), exist=exists) + if (.not.exists) then + write(logunit, '(a)') trim(subname)//' ERROR: rpointer file does not exist' + call shr_sys_abort(trim(subname)//' ERROR: rpointer file missing') + endif + open(newunit=nu, file=trim(rpfile)//trim(inst_suffix), form='formatted') + read(nu,'(a)') restfilem + close(nu) + inquire(file=trim(restfilem), exist=exists) + endif + call ESMF_VMBroadCast(vm, restfilem, CL, main_task, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + ! use namelist already read + if (my_task == main_task) then + write(logunit, '(a)') trim(subname)//' restart filenames from namelist ' + inquire(file=trim(restfilem), exist=exists) + endif + endif + tmp = 0 + if(exists) tmp=1 + exists = (tmp(1) == 1) + if (.not. exists .and. my_task == main_task) then + write(logunit, '(a)') trim(subname)//' file not found, skipping '//trim(restfilem) + return + end if + + ! Read restart file + if (my_task == main_task) then + write(logunit, '(a)') trim(subname)//' reading data model restart '//trim(restfilem) + end if + + rcode = pio_openfile(pio_subsystem, pioid, io_type, trim(restfilem), pio_nowrite) + do ns = 1,num_icesheets + + write(cnum,'(i0)') ns + + ! Determine gindex for this ice sheet + call ESMF_MeshGet(model_meshes(ns), elementdistGrid=distGrid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_DistGridGet(distGrid, localDe=0, elementCount=lsize, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + allocate(gindex(lsize)) + call ESMF_DistGridGet(distGrid, localDe=0, seqIndexList=gindex, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + + call pio_initdecomp(pio_subsystem, pio_double, (/nx_global(ns),ny_global(ns)/), gindex, pio_iodesc) + rcode = pio_inq_varid(pioid, 'flgl_rofi'//cnum, varid) + call pio_read_darray(pioid, varid, pio_iodesc, Fgrg_rofi(ns)%ptr, rcode) + call pio_freedecomp(pio_subsystem, pio_iodesc) + + ! Deallocate gindex + deallocate(gindex) + + end do ! loop over ice sheets + call pio_closefile(pioid) + + end subroutine dglc_datamode_noevolve_restart_read + +end module dglc_datamode_noevolve_mod diff --git a/dglc/glc_comp_nuopc.F90 b/dglc/glc_comp_nuopc.F90 new file mode 100644 index 000000000..68623456c --- /dev/null +++ b/dglc/glc_comp_nuopc.F90 @@ -0,0 +1,838 @@ +#ifdef CESMCOUPLED +module glc_comp_nuopc +#else +module cdeps_dglc_comp +#endif + + !---------------------------------------------------------------------------- + ! This is the NUOPC cap for DGLC + !---------------------------------------------------------------------------- + use ESMF , only : ESMF_VM, ESMF_VmGet, ESMF_VMBroadcast, ESMF_GridCompGet + use ESMF , only : ESMF_Mesh, ESMF_GridComp, ESMF_State, ESMF_Clock, ESMF_Time + use ESMF , only : ESMF_SUCCESS, ESMF_FAILURE, ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_METHOD_INITIALIZE + use ESMF , only : ESMF_TraceRegionEnter, ESMF_TraceRegionExit + use ESMF , only : ESMF_TimeGet, ESMF_TimeInterval, ESMF_TimeIntervalGet, ESMF_Field, ESMF_MAXSTR + use ESMF , only : ESMF_MethodRemove, ESMF_MethodAdd + use ESMF , only : ESMF_GridCompSetEntryPoint + use ESMF , only : ESMF_Alarm, ESMF_AlarmSet, ESMF_AlarmIsRinging, ESMF_ALARMLIST_ALL + use ESMF , only : ESMF_ClockGet, ESMF_ClockSet, ESMF_ClockAdvance, ESMF_ClockGetAlarm, ESMF_ClockGetAlarmList + use ESMF , only : ESMF_StateGet, operator(+), ESMF_AlarmRingerOff, ESMF_LogWrite + use ESMF , only : ESMF_Field, ESMF_FieldGet, ESMF_VmLogMemInfo + use ESMF , only : ESMF_MeshCreate, ESMF_FILEFORMAT_ESMFMESH + use NUOPC , only : NUOPC_CompDerive, NUOPC_CompSetEntryPoint, NUOPC_CompSpecialize + use NUOPC , only : NUOPC_Advertise, NUOPC_CompAttributeGet, NUOPC_CompAttributeSet + use NUOPC , only : NUOPC_AddNestedState, NUOPC_IsConnected + use NUOPC_Model , only : model_routine_SS => SetServices + use NUOPC_Model , only : model_label_Advance => label_Advance + use NUOPC_Model , only : model_label_SetRunClock => label_SetRunClock + use NUOPC_Model , only : model_label_Finalize => label_Finalize + use NUOPC_Model , only : NUOPC_ModelGet, SetVM + use shr_kind_mod , only : r8=>shr_kind_r8, i8=>shr_kind_i8, cl=>shr_kind_cl, cs=>shr_kind_cs + use shr_sys_mod , only : shr_sys_abort + use shr_cal_mod , only : shr_cal_ymd2date + use shr_log_mod , only : shr_log_setLogUnit + use shr_string_mod , only : shr_string_listGetNum, shr_string_listGetName +#ifdef CESMCOUPLED + use shr_pio_mod , only : shr_pio_getiosys, shr_pio_getiotype, shr_pio_getioformat +#endif + use dshr_methods_mod , only : dshr_state_diagnose, chkerr, memcheck + use dshr_strdata_mod , only : shr_strdata_type, shr_strdata_advance, shr_strdata_init_from_config + use dshr_mod , only : dshr_model_initphase, dshr_init, dshr_mesh_init, dshr_alarm_init + use dshr_mod , only : dshr_state_setscalar, dshr_set_runclock, dshr_check_restart_alarm + use dshr_dfield_mod , only : dfield_type, dshr_dfield_add, dshr_dfield_copy + use dshr_fldlist_mod , only : fldlist_type, dshr_fldlist_realize + + ! Datamode specialized modules + use dglc_datamode_noevolve_mod, only : dglc_datamode_noevolve_advertise + use dglc_datamode_noevolve_mod, only : dglc_datamode_noevolve_init_pointers + use dglc_datamode_noevolve_mod, only : dglc_datamode_noevolve_advance + use dglc_datamode_noevolve_mod, only : dglc_datamode_noevolve_restart_read + use dglc_datamode_noevolve_mod, only : dglc_datamode_noevolve_restart_write + + implicit none + private ! except + + public :: SetServices + public :: SetVM + private :: InitializeAdvertise + private :: InitializeRealize + private :: ModelSetRunClock + private :: ModelAdvance + private :: dglc_comp_run + private :: ModelFinalize + + !-------------------------------------------------------------------------- + ! Private module data + !-------------------------------------------------------------------------- + + character(*) , parameter :: nullstr = 'null' + integer , parameter :: max_icesheets = 10 ! maximum number of ice sheets for namelist input + integer :: num_icesheets ! actual number of ice sheets + + ! namelist input + character(CL) :: model_meshfiles(max_icesheets) = nullstr + character(CL) :: model_datafiles(max_icesheets) = nullstr + integer :: nx_global(max_icesheets) = 0 + integer :: ny_global(max_icesheets) = 0 + real(r8) :: model_internal_gridsize(max_icesheets) = 1.e36 + + ! module variables for multiple ice sheets + type(shr_strdata_type) , allocatable :: sdat(:) + type(ESMF_State) , allocatable :: NStateImp(:) + type(ESMF_State) , allocatable :: NStateExp(:) + type(ESMF_Mesh) , allocatable :: model_meshes(:) + + ! module variables common to all data models + character(CS) :: flds_scalar_name = '' + integer :: flds_scalar_num = 0 + integer :: flds_scalar_index_nx = 0 + integer :: flds_scalar_index_ny = 0 + integer :: mpicom ! mpi communicator + integer :: my_task ! my task in mpi communicator mpicom + logical :: mainproc ! true of my_task == main_task + character(len=16) :: inst_suffix = "" ! char string associated with instance (ie. "_0001" or "") + integer :: logunit ! logging unit number + logical :: restart_read ! start from restart + character(CL) :: case_name + + ! dglc_in namelist input + character(CL) :: streamfilename = nullstr ! filename to obtain stream info from + character(CL) :: nlfilename = nullstr ! filename to obtain namelist info from + character(CL) :: datamode = nullstr ! flags physics options wrt input data + character(CL) :: restfilm = nullstr ! model restart file namelist + logical :: skip_restart_read = .false. ! true => skip restart read in continuation run + logical :: export_all = .false. ! true => export all fields, do not check connected or not + + ! linked lists + type(fldList_type) , pointer :: fldsImport => null() + type(fldList_type) , pointer :: fldsExport => null() + + type dfields_icesheets_type + type(dfield_type), pointer :: dfields => null() + end type dfields_icesheets_type + type(dfields_icesheets_type), allocatable :: dfields_icesheets(:) + + ! constants + logical :: diagnose_data = .true. + integer , parameter :: main_task = 0 ! task number of main task +#ifdef CESMCOUPLED + character(*) , parameter :: module_name = "(glc_comp_nuopc)" +#else + character(*) , parameter :: module_name = "(cdeps_dglc_comp)" +#endif + character(*) , parameter :: modelname = 'dglc' + character(*) , parameter :: u_FILE_u = & + __FILE__ + +!=============================================================================== +contains +!=============================================================================== + + subroutine SetServices(gcomp, rc) + type(ESMF_GridComp) :: gcomp + integer, intent(out) :: rc + + ! Local varaibles + character(len=*),parameter :: subname=trim(module_name)//':(SetServices) ' + !-------------------------------- + + rc = ESMF_SUCCESS + call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO) + + ! the NUOPC gcomp component will register the generic methods + call NUOPC_CompDerive(gcomp, model_routine_SS, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! switching to IPD versions + call ESMF_GridCompSetEntryPoint(gcomp, ESMF_METHOD_INITIALIZE, & + userRoutine=dshr_model_initphase, phase=0, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! set entry point for methods that require specific implementation + call NUOPC_CompSetEntryPoint(gcomp, ESMF_METHOD_INITIALIZE, & + phaseLabelList=(/"IPDv01p1"/), userRoutine=InitializeAdvertise, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call NUOPC_CompSetEntryPoint(gcomp, ESMF_METHOD_INITIALIZE, & + phaseLabelList=(/"IPDv01p3"/), userRoutine=InitializeRealize, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! attach specializing method(s) + call NUOPC_CompSpecialize(gcomp, specLabel=model_label_Advance, specRoutine=ModelAdvance, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_MethodRemove(gcomp, label=model_label_SetRunClock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! The following specialization does not use dshr_set_runclock since we also need to add an alarm + ! to determine if the model inputs are valid + call NUOPC_CompSpecialize(gcomp, specLabel=model_label_SetRunClock, specRoutine=ModelSetRunClock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call NUOPC_CompSpecialize(gcomp, specLabel=model_label_Finalize, specRoutine=ModelFinalize, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_LogWrite(subname//' done', ESMF_LOGMSG_INFO) + + end subroutine SetServices + + !=============================================================================== + subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) + + use shr_nl_mod, only: shr_nl_find_group_name + + ! input/output variables + type(ESMF_GridComp) :: gcomp + type(ESMF_State) :: importState, exportState + type(ESMF_Clock) :: clock + integer, intent(out) :: rc + + ! local variables + type(ESMF_VM) :: vm + integer :: inst_index ! number of current instance (ie. 1) + integer :: nu ! unit number + integer :: ierr ! error code + integer :: bcasttmp(2) + integer :: ns + character(len=CS) :: cnum + character(len=ESMF_MAXSTR) :: model_datafiles_list ! colon separated string containing input datafiles + character(len=ESMF_MAXSTR) :: model_meshfiles_list ! colon separated string containing model meshfiles + character(len=*),parameter :: subname=trim(module_name)//':(InitializeAdvertise) ' + !------------------------------------------------------------------------------- + + ! Note that the suffix '-list' refers to a colon delimited string of names + namelist / dglc_nml / datamode, & + model_meshfiles_list, model_datafiles_list, model_internal_gridsize, nx_global, ny_global, & + restfilm, skip_restart_read + + rc = ESMF_SUCCESS + + call NUOPC_CompAttributeGet(gcomp, name='case_name', value=case_name, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Determine logical mainproc + mainproc = (my_task == main_task) + + ! Obtain flds_scalar values, mpi values, multi-instance values and + ! set logunit and set shr logging to my log file + call dshr_init(gcomp, 'GLC', mpicom, my_task, inst_index, inst_suffix, & + flds_scalar_name, flds_scalar_num, flds_scalar_index_nx, flds_scalar_index_ny, logunit, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Read dglc_nml from nlfilename + if (my_task == main_task) then + nlfilename = "dglc_in"//trim(inst_suffix) + open (newunit=nu,file=trim(nlfilename),status="old",action="read") + call shr_nl_find_group_name(nu, 'dglc_nml', status=ierr) + read (nu,nml=dglc_nml,iostat=ierr) + close(nu) + if (ierr > 0) then + write(logunit,'(a,i8)') 'ERROR: reading input namelist, '//trim(nlfilename)//' iostat=',ierr + call shr_sys_abort(subName//': namelist read error '//trim(nlfilename)) + end if + + ! Determine number of ice sheets + num_icesheets = shr_string_listGetNum(model_meshfiles_list) + + ! Determine array of model meshfile(s) and model input datafile(s) + do ns = 1,num_icesheets + ! determine mesh filename(s) + call shr_string_listGetName(model_meshfiles_list, ns, model_meshfiles(ns)) + ! determine input datafile name(s) + call shr_string_listGetName(model_datafiles_list, ns, model_datafiles(ns)) + end do + + ! Write diagnostics + write(logunit,'(a,a)')' case_name = ',trim(case_name) + write(logunit,'(a,a)')' datamode = ',trim(datamode) + do ns = 1,num_icesheets + write(logunit,'(a,i4 )') ' ice_sheet index = ',ns + write(logunit,'(a,a )') ' model_meshfile = ',trim(model_meshfiles(ns)) + write(logunit,'(a,a )') ' model_datafile = ',trim(model_datafiles(ns)) + write(logunit,'(a,d13.5)')' model_internal_gridsize = ',model_internal_gridsize(ns) + write(logunit,'(a,i10)') ' nx_global = ',nx_global(ns) + write(logunit,'(a,i10)') ' ny_global = ',ny_global(ns) + end do + write(logunit,'(a,a )')' restfilm = ',trim(restfilm) + write(logunit,'(a,l6)')' skip_restart_read = ',skip_restart_read + bcasttmp(1) = 0 + if(skip_restart_read) bcasttmp(1) = 1 + bcasttmp(2) = num_icesheets + endif + + ! Broadcast namelist input + call ESMF_GridCompGet(gcomp, vm=vm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMBroadcast(vm, datamode, CL, main_task, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMBroadcast(vm, restfilm, CL, main_task, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMBroadcast(vm, model_meshfiles, CL*max_icesheets, main_task, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMBroadcast(vm, model_datafiles, CL*max_icesheets, main_task, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMBroadcast(vm, model_internal_gridsize, max_icesheets, main_task, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMBroadcast(vm, nx_global, max_icesheets, main_task, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMBroadcast(vm, ny_global, max_icesheets, main_task, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMBroadcast(vm, bcasttmp, 3, main_task, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + skip_restart_read = (bcasttmp(1) == 1) + num_icesheets = bcasttmp(2) + + ! Validate datamode + if ( trim(datamode) == 'noevolve') then ! read stream, no import data + ! do nothing + else + call shr_sys_abort(' ERROR illegal dglc datamode = '//trim(datamode)) + endif + + ! Allocate module variables + allocate(sdat(num_icesheets)) + allocate(NStateImp(num_icesheets)) + allocate(NStateExp(num_icesheets)) + allocate(model_meshes(num_icesheets)) + + ! Create nested states + do ns = 1,num_icesheets + write(cnum,'(i0)') ns + call NUOPC_AddNestedState(importState, CplSet="GLC"//trim(cnum), nestedState=NStateImp(ns), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call NUOPC_AddNestedState(exportState, CplSet="GLC"//trim(cnum), nestedState=NStateExp(ns), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end do + + ! Advertise dglc fields + if (trim(datamode)=='noevolve') then + call dglc_datamode_noevolve_advertise(NStateExp, fldsexport, NStateImp, fldsimport, & + flds_scalar_name, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + + end subroutine InitializeAdvertise + + !=============================================================================== + subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) + + ! input/output variables + type(ESMF_GridComp) :: gcomp + type(ESMF_State) :: importState, exportState + type(ESMF_Clock) :: clock + integer, intent(out) :: rc + + ! local variables + type(ESMF_VM) :: vm + type(ESMF_Time) :: currTime + integer :: current_ymd ! model date + integer :: current_year ! model year + integer :: current_mon ! model month + integer :: current_day ! model day + integer :: current_tod ! model sec into model date + integer :: ns + character(CL) :: cvalue + character(CS) :: cns + logical :: ispresent, isset + logical :: exists + character(len=*), parameter :: subname=trim(module_name)//':(InitializeRealize) ' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! Initialize model mesh, restart flag, logunit, model_mask and model_frac + call ESMF_VMLogMemInfo("Entering "//trim(subname)) + call ESMF_TraceRegionEnter('dglc_strdata_init') + + ! Determine stream filename + streamfilename = trim(modelname)//'.streams'//trim(inst_suffix) +#ifndef DISABLE_FoX + streamfilename = trim(streamfilename)//'.xml' +#endif + + ! generate local mpi comm + call ESMF_GridCompGet(gcomp, vm=vm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMGet(vm, localPet=my_task, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + mainproc = (my_task == main_task) + call shr_log_setLogUnit(logunit) + + ! Set restart flag (sets module variable restart_read) + call NUOPC_CompAttributeGet(gcomp, name='read_restart', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + read(cvalue,*) restart_read + else + call shr_sys_abort(subname//' ERROR: read restart flag must be present') + end if + + ! Get the time to interpolate the stream data to + call ESMF_ClockGet(clock, currTime=currTime, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_TimeGet(currTime, yy=current_year, mm=current_mon, dd=current_day, s=current_tod, rc=rc ) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call shr_cal_ymd2date(current_year, current_mon, current_day, current_ymd) + + ! Loop over ice sheets + do ns = 1,num_icesheets + + write(cns,'(i0)') ns + + ! Initialize pio subsystem +#ifdef CESMCOUPLED + sdat(ns)%pio_subsystem => shr_pio_getiosys('GLC') + sdat(ns)%io_type = shr_pio_getiotype('GLC') + sdat(ns)%io_format = shr_pio_getioformat('GLC') +#else + call dshr_pio_init(gcomp, sdat(ns), logunit, rc) +#endif + + ! Check that model_meshfile exists + if (my_task == main_task) then + inquire(file=trim(model_meshfiles(ns)), exist=exists) + if (.not.exists) then + write(logunit,'(a)')' ERROR: model_meshfile '//trim(model_meshfiles(ns))//' does not exist' + call shr_sys_abort(trim(subname)//' ERROR: model_meshfile '//trim(model_meshfiles(ns))//' does not exist') + end if + endif + + ! Read in model mesh for given ice sheet + model_meshes(ns) = ESMF_MeshCreate(trim(model_meshfiles(ns)), fileformat=ESMF_FILEFORMAT_ESMFMESH, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Initialize stream data type + if (trim(datamode) /= 'noevolve') then + call shr_strdata_init_from_config(sdat(ns), streamfilename, model_meshes(ns), clock, 'GLC', logunit, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + + ! Realize the actively coupled fields, now that a mesh is established and + ! NUOPC_Realize "realizes" a previously advertised field in the importState and exportState + ! by replacing the advertised fields with the newly created fields of the same name. + + call ESMF_LogWrite(subname//' calling dshr_fldlist_realize export for ice sheet '//trim(cns), ESMF_LOGMSG_INFO) + call dshr_fldlist_realize( NStateExp(ns), fldsExport, flds_scalar_name, flds_scalar_num, model_meshes(ns), & + subname//trim(modelname)//':Export', export_all, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_LogWrite(subname//' calling dshr_fldlist_realize importfor ice sheet '//trim(cns), ESMF_LOGMSG_INFO) + call dshr_fldlist_realize( NStateImp(ns), fldsImport, flds_scalar_name, flds_scalar_num, model_meshes(ns), & + subname//trim(modelname)//':Import', .false., rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Add scalars to export state + call dshr_state_SetScalar(dble(nx_global(ns)),flds_scalar_index_nx, & + NStateExp(ns), flds_scalar_name, flds_scalar_num, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call dshr_state_SetScalar(dble(ny_global(ns)),flds_scalar_index_ny,& + NStateExp(ns), flds_scalar_name, flds_scalar_num, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + end do ! end loop over ice sheets + + ! Run dglc + call dglc_comp_run(gcomp, clock, current_ymd, current_tod, restart_write=.false., valid_inputs=.true., rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_TraceRegionExit('dglc_strdata_init') + call ESMF_VMLogMemInfo("Leaving "//trim(subname)) + + end subroutine InitializeRealize + + !=============================================================================== + subroutine ModelAdvance(gcomp, rc) + + ! input/output variables + type(ESMF_GridComp) :: gcomp + integer, intent(out) :: rc + + ! local variables + type(ESMF_Clock) :: clock + type(ESMF_TimeInterval) :: timeStep + type(ESMF_Time) :: currTime, nextTime + integer :: next_ymd ! model date + integer :: next_tod ! model sec into model date + integer :: yr ! year + integer :: mon ! month + integer :: day ! day in month + logical :: restart_write + type(ESMF_Alarm) :: valid_alarm + logical :: valid_inputs ! if true, inputs from mediator are valid + character(len=CS) :: timestring + character(len=*),parameter :: subname=trim(module_name)//':(ModelAdvance) ' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + call shr_log_setLogUnit(logunit) + + call ESMF_TraceRegionEnter(subname) + call memcheck(subname, 5, my_task == main_task) + + ! query the Component for its clock + call NUOPC_ModelGet(gcomp, modelClock=clock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Determine if inputs from mediator are valid + call ESMF_ClockGetAlarm(clock, alarmname='alarm_valid_inputs', alarm=valid_alarm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (ESMF_AlarmIsRinging(valid_alarm, rc=rc)) then + valid_inputs = .true. + call ESMF_AlarmRingerOff( valid_alarm, rc=rc ) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + valid_inputs = .false. + endif + + ! Need to advance nuopc one timestep ahead for shr_strdata time interpolation + call ESMF_ClockGet( clock, currTime=currTime, timeStep=timeStep, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + nextTime = currTime + timeStep + call ESMF_TimeGet( nextTime, yy=yr, mm=mon, dd=day, s=next_tod, rc=rc ) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call shr_cal_ymd2date(yr, mon, day, next_ymd) + if (my_task == main_task) then + call ESMF_TimeGet(currTime, timestring=timestring, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + write(logunit,'(a,l6)') trim(timestring)//': valid_input for dglc is ',valid_inputs + end if + + ! determine if will write restart + restart_write = dshr_check_restart_alarm(clock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (my_task == main_task) then + write(logunit,'(a,l6)') trim(timestring)//': restart write is ',restart_write + end if + + ! run dglc + call dglc_comp_run(gcomp, clock, next_ymd, next_tod, restart_write, valid_inputs, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_TraceRegionExit(subname) + + end subroutine ModelAdvance + + !=============================================================================== + subroutine dglc_comp_run(gcomp, clock, target_ymd, target_tod, restart_write, valid_inputs, rc) + + ! -------------------------- + ! advance dglc + ! -------------------------- + + ! input/output variables: + type(ESMF_GridComp) :: gcomp + type(ESMF_Clock) , intent(in) :: clock + integer , intent(in) :: target_ymd ! model date + integer , intent(in) :: target_tod ! model sec into model date + logical , intent(in) :: restart_write + logical , intent(in) :: valid_inputs + integer , intent(out) :: rc + + ! local variables + character(len=CS) :: cnum + integer :: ns ! ice sheet index + logical :: first_time = .true. + character(*), parameter :: subName = "(dglc_comp_run) " + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + call ESMF_TraceRegionEnter('DGLC_RUN') + + !-------------------- + ! First time initialization + !-------------------- + + if (first_time) then + ! Initialize dfields for all ice sheets + if (trim(datamode) /= 'noevolve') then + call dglc_init_dfields(rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + + ! Initialize datamode module ponters + select case (trim(datamode)) + case('noevolve') + call dglc_datamode_noevolve_init_pointers(NStateExp, NStateImp, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end select + + ! Read restart if needed + if (restart_read .and. .not. skip_restart_read) then + write(logunit,'(a)')' DEBUG: calling dglc_datamode_noevolve_restart_read' + call dglc_datamode_noevolve_restart_read(model_meshes, restfilm, & + inst_suffix, logunit, my_task, main_task, mpicom, & + sdat(1)%pio_subsystem, sdat(1)%io_type, nx_global, ny_global, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + + ! Reset first_time + first_time = .false. + end if + + !-------------------- + ! Update export (and possibly import data model states) + !-------------------- + + if (trim(datamode) /= 'noevolve') then + if (.not. allocated(dfields_icesheets)) then + allocate(dfields_icesheets(num_icesheets)) + end if + + ! Loop over ice sheets + do ns = 1,num_icesheets + ! Advance data model streams - time and spatially interpolate to model time and grid + ! Note that loop over ice sheets is done inside shr_strdata_advance + call ESMF_TraceRegionEnter('dglc_strdata_advance') + call shr_strdata_advance(sdat(ns), target_ymd, target_tod, logunit, 'dglc', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_TraceRegionExit('dglc_strdata_advance') + + ! Copy all fields from streams to export state as default + ! This automatically will update the fields in the export state + call ESMF_TraceRegionEnter('dglc_dfield_copy') + call dshr_dfield_copy(dfields_icesheets(ns)%dfields, sdat(ns), rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_TraceRegionExit('dglc_dfield_copy') + end do + end if + + ! Perform data mode specific calculations + select case (trim(datamode)) + case('noevolve') + if (valid_inputs) then + call dglc_datamode_noevolve_advance(gcomp, sdat(1)%pio_subsystem, sdat(1)%io_type, sdat(1)%io_format, & + logunit, model_meshes, model_internal_gridsize, model_datafiles, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + end select + + ! Write restarts if needed + + if (restart_write) then + if (trim(datamode) == 'noevolve') then + if (my_task == main_task) then + write(logunit,'(a)') 'calling dglc_datamode_noevolve_restart_write' + end if + call dglc_datamode_noevolve_restart_write(model_meshes, case_name, & + inst_suffix, target_ymd, target_tod, logunit, my_task, main_task, & + sdat(1)%pio_subsystem, sdat(1)%io_type, nx_global, ny_global, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + end if + + ! Write diagnostics + if (diagnose_data) then + do ns = 1,num_icesheets + write(cnum,'(i0)') ns + call dshr_state_diagnose(NStateExp(ns), flds_scalar_name, trim(subname)//':ES_'//trim(cnum), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end do + end if + + call ESMF_TraceRegionExit('DGLC_RUN') + + contains + + subroutine dglc_init_dfields(rc) + ! ----------------------------- + ! Initialize dfields arrays + ! ----------------------------- + + ! input/output variables + integer, intent(out) :: rc + + ! local variables + integer :: nf, ns + integer :: fieldcount + type(ESMF_Field) :: lfield + character(ESMF_MAXSTR) ,pointer :: lfieldnamelist(:) + character(*), parameter :: subName = "(dglc_init_dfields) " + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! Loop over ice sheets + ! Initialize dfields data type (to map streams to export state fields) + ! Create dfields linked list - used for copying stream fields to export state fields + do ns = 1,num_icesheets + call ESMF_StateGet(NStateExp(ns), itemCount=fieldCount, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + allocate(lfieldnamelist(fieldCount)) + call ESMF_StateGet(NStateExp(ns), itemNameList=lfieldnamelist, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + do nf = 1, fieldCount + call ESMF_StateGet(NStateExp(ns), itemName=trim(lfieldNameList(nf)), field=lfield, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + if (trim(lfieldnamelist(nf)) /= flds_scalar_name) then + call dshr_dfield_add( dfields_icesheets(ns)%dfields, sdat(ns), & + trim(lfieldnamelist(nf)), trim(lfieldnamelist(nf)), NStateExp(ns), logunit, mainproc, rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + end if + end do + deallocate(lfieldnamelist) + end do + end subroutine dglc_init_dfields + + end subroutine dglc_comp_run + + !=============================================================================== + + subroutine ModelSetRunClock(gcomp, rc) + + ! input/output variables + type(ESMF_GridComp) :: gcomp + integer, intent(out) :: rc + + ! local variables + type(ESMF_Clock) :: mclock, dclock + type(ESMF_Time) :: mcurrtime, dcurrtime + type(ESMF_Time) :: mstoptime + type(ESMF_TimeInterval) :: mtimestep, dtimestep + character(len=256) :: cvalue + character(len=256) :: restart_option ! Restart option units + integer :: restart_n ! Number until restart interval + integer :: restart_ymd ! Restart date (YYYYMMDD) + type(ESMF_ALARM) :: restart_alarm + character(len=256) :: stop_option ! Stop option units + integer :: stop_n ! Number until stop interval + integer :: stop_ymd ! Stop date (YYYYMMDD) + type(ESMF_ALARM) :: stop_alarm + integer :: alarmcount + character(len=CS) :: glc_avg_period ! averaging period in mediator + type(ESMF_ALARM) :: valid_alarm ! model alarm + integer :: dtime + character(len=*),parameter :: subname='glc_comp_nuopc:(dglc_set_runclock) ' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! query the Component for its clocks + call NUOPC_ModelGet(gcomp, driverClock=dclock, modelClock=mclock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_ClockGet(dclock, currTime=dcurrtime, timeStep=dtimestep, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_ClockGet(mclock, currTime=mcurrtime, timeStep=mtimestep, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! force model clock currtime and timestep to match driver and set stoptime + mstoptime = mcurrtime + dtimestep + call ESMF_ClockSet(mclock, currTime=dcurrtime, timeStep=dtimestep, stopTime=mstoptime, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! determine number of alarms + call ESMF_ClockGetAlarmList(mclock, alarmlistflag=ESMF_ALARMLIST_ALL, alarmCount=alarmCount, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + if (alarmCount == 0) then + + !---------------- + ! glc valid input alarm + !---------------- + call NUOPC_CompAttributeGet(gcomp, name="glc_avg_period", value=glc_avg_period, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + if (trim(glc_avg_period) == 'hour') then + call dshr_alarm_init(mclock, valid_alarm, 'nhours', opt_n=1, alarmname='alarm_valid_inputs', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else if (trim(glc_avg_period) == 'day') then + call dshr_alarm_init(mclock, valid_alarm, 'ndays' , opt_n=1, alarmname='alarm_valid_inputs', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else if (trim(glc_avg_period) == 'yearly') then + call dshr_alarm_init(mclock, valid_alarm, 'yearly', alarmname='alarm_valid_inputs', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else if (trim(glc_avg_period) == 'glc_coupling_period') then + call ESMF_TimeIntervalGet(mtimestep, s=dtime, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call dshr_alarm_init(mclock, valid_alarm, 'nseconds', opt_n=dtime, alarmname='alarm_valid_inputs', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + call ESMF_LogWrite(trim(subname)// ": ERROR glc_avg_period = "//trim(glc_avg_period)//" not supported", & + ESMF_LOGMSG_INFO, rc=rc) + rc = ESMF_FAILURE + RETURN + end if + + call ESMF_AlarmSet(valid_alarm, clock=mclock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + !---------------- + ! Restart alarm + !---------------- + call ESMF_LogWrite(subname//'setting restart alarm for dglc' , ESMF_LOGMSG_INFO) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call NUOPC_CompAttributeGet(gcomp, name="restart_option", value=restart_option, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call NUOPC_CompAttributeGet(gcomp, name="restart_n", value=cvalue, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) restart_n + + call NUOPC_CompAttributeGet(gcomp, name="restart_ymd", value=cvalue, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) restart_ymd + + call dshr_alarm_init(mclock, restart_alarm, restart_option, & + opt_n = restart_n, & + opt_ymd = restart_ymd, & + RefTime = mcurrTime, & + alarmname = 'alarm_restart', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_AlarmSet(restart_alarm, clock=mclock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + !---------------- + ! Stop alarm + !---------------- + call ESMF_LogWrite(subname//'setting stop alarm for dglc' , ESMF_LOGMSG_INFO) + call NUOPC_CompAttributeGet(gcomp, name="stop_option", value=stop_option, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call NUOPC_CompAttributeGet(gcomp, name="stop_n", value=cvalue, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) stop_n + + call NUOPC_CompAttributeGet(gcomp, name="stop_ymd", value=cvalue, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) stop_ymd + + call dshr_alarm_init(mclock, stop_alarm, stop_option, & + opt_n = stop_n, & + opt_ymd = stop_ymd, & + RefTime = mcurrTime, & + alarmname = 'alarm_stop', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_AlarmSet(stop_alarm, clock=mclock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + end if + + ! Advance model clock to trigger alarms then reset model clock back to currtime + call ESMF_ClockAdvance(mclock,rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_ClockSet(mclock, currTime=dcurrtime, timeStep=dtimestep, stopTime=mstoptime, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + end subroutine ModelSetRunClock + + !=============================================================================== + subroutine ModelFinalize(gcomp, rc) + type(ESMF_GridComp) :: gcomp + integer, intent(out) :: rc + !------------------------------------------------------------------------------- + rc = ESMF_SUCCESS + if (my_task == main_task) then + write(logunit,*) + write(logunit,*) 'dglc : end of main integration loop' + write(logunit,*) + end if + end subroutine ModelFinalize + +#ifdef CESMCOUPLED +end module glc_comp_nuopc +#else +end module cdeps_dglc_comp +#endif diff --git a/dice/cime_config/testdefs/testlist_dice.xml b/dice/cime_config/testdefs/testlist_dice.xml index eac359d4e..948dc092d 100644 --- a/dice/cime_config/testdefs/testlist_dice.xml +++ b/dice/cime_config/testdefs/testlist_dice.xml @@ -1,17 +1,19 @@ - + - + + - + - + + diff --git a/dice/ice_comp_nuopc.F90 b/dice/ice_comp_nuopc.F90 index 567a1943c..c347070af 100644 --- a/dice/ice_comp_nuopc.F90 +++ b/dice/ice_comp_nuopc.F90 @@ -202,7 +202,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ! Obtain flds_scalar values, mpi values, multi-instance values and ! set logunit and set shr logging to my log file - call dshr_init(gcomp, 'ICE', sdat, mpicom, my_task, inst_index, inst_suffix, & + call dshr_init(gcomp, 'ICE', mpicom, my_task, inst_index, inst_suffix, & flds_scalar_name, flds_scalar_num, flds_scalar_index_nx, flds_scalar_index_ny, & logunit, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return diff --git a/dlnd/cime_config/testdefs/testlist_dlnd.xml b/dlnd/cime_config/testdefs/testlist_dlnd.xml index de52bac4a..59026878c 100644 --- a/dlnd/cime_config/testdefs/testlist_dlnd.xml +++ b/dlnd/cime_config/testdefs/testlist_dlnd.xml @@ -1,9 +1,10 @@ - + - + + diff --git a/dlnd/lnd_comp_nuopc.F90 b/dlnd/lnd_comp_nuopc.F90 index 5fe855ddf..79bdf0621 100644 --- a/dlnd/lnd_comp_nuopc.F90 +++ b/dlnd/lnd_comp_nuopc.F90 @@ -187,7 +187,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ! Obtain flds_scalar values, mpi values, multi-instance values and ! set logunit and set shr logging to my log file - call dshr_init(gcomp, 'LND', sdat, mpicom, my_task, inst_index, inst_suffix, & + call dshr_init(gcomp, 'LND', mpicom, my_task, inst_index, inst_suffix, & flds_scalar_name, flds_scalar_num, flds_scalar_index_nx, flds_scalar_index_ny, & logunit, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return diff --git a/doc/source/datm.rst b/doc/source/datm.rst index ff067209c..d731bc689 100644 --- a/doc/source/datm.rst +++ b/doc/source/datm.rst @@ -60,7 +60,7 @@ ERA5 (``datm_datamode_era5_mod.F90``) .. note:: Due to the high temporal and spatial resoultion of ERA5 dataset, only 2019 - data is staged on NCAR's Cheyenne platform under + data is staged on NCAR's Derecho platform under `$CESMDATAROOT/inputdata/atm/datm7/ERA5` .. note:: @@ -176,6 +176,9 @@ DATM_CPLHIST_CASE DATM_YR_START - Starting year to loop data over +DATM_YR_START_FILENAME + - Start year listed in PLUMBER2 filename + DATM_YR_END - Ending year to loop data over diff --git a/doc/source/dglc.rst b/doc/source/dglc.rst new file mode 100644 index 000000000..51175681e --- /dev/null +++ b/doc/source/dglc.rst @@ -0,0 +1,117 @@ +.. _dglc: + +Data Land-Ice (DGLC) +================= + +DGLC is normally used as a substitute for CESM/CISM running in NOEVOLVE mode. +The various ways of running DGLC is referred to as its mode. + +NOEVOLVE mode is used in CESM as follows. +In typical runs, CISM is not evolving; CLM computes the surface mass +balance (SMB) and sends it to CISM, but CISM’s ice sheet geometry +remains fixed over the course of the run. In these runs, CISM serves +two roles in the system: + + - Over the CISM domain (typically Greenland in + CESM2), CISM dictates glacier areas and topographic elevations, + overriding the values on CLM’s surface dataset. CISM also dictates the + elevation of non-glacier land units in its domain, and only in this + domain are atmospheric fields downscaled to non-glacier land + units. (So if you run with a stub glacier model - SGLC - then glacier + areas and elevations will be taken entirely from CLM’s surface + dataset, and no downscaling will be done over non-glacier land units.) + + - CISM provides the grid onto which SMB is downscaled. (If you run with + SGLC then SMB will still be computed in CLM, but it won’t be + downscaled to a high-resolution ice sheet grid.) + +-------------------- +Supported Data Modes +-------------------- + +DGLC has the capability of supporting multiple ice sheets (as is the +case with CISM/CMEPS coupling). This is configured via the following +``dglc_in`` namelist settings: + + - ``model_meshfiles_list`` is a colon separated string containing model + meshfiles describing the different ice sheets. + + - ``model_datafiles_list`` is colon separated string containing + input datafiles that specify are used to obtain data for bedrock + topography and the ice thickness. + + - ``model_internal_gridsize`` is an array that is the size of the number of ice + sheets and that specifies the internal gridcell size that corresponds + what internal gridcell areas the prognostic land-ice component + uses internally (in this case CISM). From this value the internal grid areas in + radians squared are given by (model_internal_gridsize/radius_earth)**2. + Both model_internal_gridsize and radius_earth have units of meters. + + + - ``nx_global`` is an array that is the size of the number of ice + sheets and that specifies the global longitude dimension of the + each ice sheet. + + - ``ny_global`` is an array that is the size of the number of ice + sheets and that specifies the global latitude dimension of the + each ice sheet. + +.. note:: + Each element of ``model_data_filelist``, ``model_areas``, + ``nx_global`` and ``ny_global`` corresponds to a different ice + sheet mesh and should have the **same order** as those in the + ``model_meshfiles_list``. + +DGLC has its own set of supported ``datamode`` values that appears in +the ``dglc_in`` namelist input. The datamode specifies what additional +operations need to be done by DGLC on *ALL* of the streams in the +``dglc.streams.xml`` file. Each datamode value is also associated with +a DGLC source file that carries out these operations and these are +listed in parentheses next to the mode name. Currently, the only +supported ``datamode`` is ``noevolve``. + +noevolve (``dglc_datamode_noevolve_mod.F90``) + - This mode is an analytic mode that has no associated stream files. + This mode uses ``dglc_in`` namelist variables as follows: + +.. _dglc-cime-vars: + +--------------------------------------- +Configuring DGLC using CIME-CCS +--------------------------------------- + +If CDEPS is coupled to the CIME-CCS then the CIME ``$CASEROOT`` xml +variable ``DGLC_MODE`` will be generated based on the compset +specification ``DGLC%{DGLC_MODE}``. ``DGLC_MODE`` will in term be +used in the ``namelist_definition_dglc.xml`` file to determine the +collection of streams that are associated with DGLC and also sets the +dglc namelist variable ``datamode`` in the file ``dglc_in``. + +The following list describes the valid values of ``DGLC_MODE`` +(defined in the ``config_component.xml`` file for DGLC), and how they +relate to the associated input streams and the ``datamode`` namelist +variable. + +DGLC%NOEVOLVE + - fields sent to mediator are created analytically without stream + input + - dglc_mode: NOEVOLVE + - streams: none + - datamode: noevolve + +In addition, the following DGLC specific CIME-CCS xml variables will appear in ``$CASEROOT/env_run.xml``: + +DGLC_USE_GREENLAND + - Whether to include the Greenland Ice Sheet in this DGLC simulation + This should generally be set at create_newcase time (via the compset). In principle it + can be changed later, but great care is needed to change a number of other variables + to be consistent (GLC_GRID, GLC_DOMAIN_MESH and possibly others). + +DGLC_USE_ANTARCTICA + - Whether to include the Antarctic Ice Sheet in this DGLC simulation + This should generally be set at create_newcase time (via the compset). In principle it + can be changed later, but great care is needed to change a number of other variables + to be consistent (GLC_GRID, GLC_DOMAIN_MESH and possibly others). + +DGLC_SKIP_RESTART_READ + - If set to true, than dglc restarts will not be read on a continuation run. diff --git a/doc/source/index.rst b/doc/source/index.rst index 3990258c3..f1a826f07 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -29,4 +29,5 @@ Table of contents dlnd.rst drof.rst dwav.rst + dglc.rst extending.rst diff --git a/doc/source/introduction.rst b/doc/source/introduction.rst index b34a67419..dcdf2ef1b 100644 --- a/doc/source/introduction.rst +++ b/doc/source/introduction.rst @@ -59,6 +59,7 @@ cime_config CIME Case Control System cmake Build (can be used with or without CIME) datm Data atmosphere component dice Data sea-ice component +dglc Data land-ice component dlnd Data land component docn Data ocean component drof Data river component @@ -132,7 +133,7 @@ CDEPS and CIME Control System (CCS) If the CDEPS data models are used in conjunction with the CIME Case Control System (CCS) then the following will also hold: Each data model has an xml variable in ``env_run.xml`` that specifies the data model mode. -These are: ``DATM_MODE``, ``DICE_MODE``, ``DLND_MODE``, ``DOCN_MODE``, ``DROF_MODE``, ``DWAV_MODE``. +These are: ``DATM_MODE``, ``DICE_MODE``, ``DGLC_MODE``, ``DLND_MODE``, ``DOCN_MODE``, ``DROF_MODE``, ``DWAV_MODE``. Each data model mode specifies the streams that are associated with that data model. More details of the data model design are covered in :ref:`design details`. diff --git a/doc/source/streams.rst b/doc/source/streams.rst index 3ce4d8023..0b6fa7235 100644 --- a/doc/source/streams.rst +++ b/doc/source/streams.rst @@ -365,7 +365,7 @@ Data Model Namelist Input ------------------------- Each data model has an associated input namelist file, ``d{model_name}_in``, -where ``model_name=[datm,dlnd,dice,docn,drof,dwav]``. +where ``model_name=[datm,dlnd,dglc,dice,docn,drof,dwav]``. The following namelist variables appear in each data model namelist: @@ -381,6 +381,8 @@ The following namelist variables appear in each data model namelist: :ref:`Data Land ` + :ref:`Data Land-Ice ` + :ref:`Data Ocean ` :ref:`Data Runoff ` diff --git a/docn/CMakeLists.txt b/docn/CMakeLists.txt index 007d595ca..46e51007f 100644 --- a/docn/CMakeLists.txt +++ b/docn/CMakeLists.txt @@ -5,6 +5,8 @@ set(SRCFILES ocn_comp_nuopc.F90 docn_datamode_aquaplanet_mod.F90 docn_datamode_iaf_mod.F90 docn_datamode_cplhist_mod.F90 + docn_datamode_multilev_mod.F90 + docn_datamode_multilev_dom_mod.F90 docn_import_data_mod.F90) foreach(FILE ${SRCFILES}) diff --git a/docn/cime_config/config_component.xml b/docn/cime_config/config_component.xml index 3a98b7878..c5461d900 100644 --- a/docn/cime_config/config_component.xml +++ b/docn/cime_config/config_component.xml @@ -13,7 +13,7 @@ This file may have ocn desc entries. --> - DOCN + DOCN prescribed ocean mode slab ocean mode aquaplanet slab ocean mode @@ -32,6 +32,8 @@ file input aquaplanet sst globally constant SST for idealized experiments, such as RCE mediator history output for ocean fields imported to mediator + input stream files have multi level data + input stream files have multi level data and prescribed ocean SST @@ -45,7 +47,7 @@ char - prescribed,sst_aquap1,sst_aquap2,sst_aquap3,sst_aquap4,sst_aquap5,sst_aquap6,sst_aquap7,sst_aquap8,sst_aquap9,sst_aquap10,sst_aquapfile,som,som_aquap,sst_aquap_constant,interannual,cplhist + prescribed,sst_aquap1,sst_aquap2,sst_aquap3,sst_aquap4,sst_aquap5,sst_aquap6,sst_aquap7,sst_aquap8,sst_aquap9,sst_aquap10,sst_aquapfile,som,som_aquap,sst_aquap_constant,interannual,cplhist,multilev,multilev_dom prescribed prescribed @@ -65,6 +67,8 @@ sst_aquapfile sst_aquap_constant cplhist + multilev + multilev_dom run_component_docn env_run.xml diff --git a/docn/cime_config/namelist_definition_docn.xml b/docn/cime_config/namelist_definition_docn.xml index a44ae344d..3e6d12558 100644 --- a/docn/cime_config/namelist_definition_docn.xml +++ b/docn/cime_config/namelist_definition_docn.xml @@ -32,6 +32,8 @@ '' '' '' + sst_salinity_depth_blom,prescribed + sst_salinity_depth_blom @@ -39,7 +41,7 @@ char docn docn_nml - sstdata,sst_aquap1,sst_aquap2,sst_aquap3,sst_aquap4,sst_aquap5,sst_aquap6,sst_aquap7,sst_aquap8,sst_aquap9,sst_aquap10,sst_aquapfile,sst_aquap_constant,som,som_aquap,iaf,cplhist + sstdata,sst_aquap1,sst_aquap2,sst_aquap3,sst_aquap4,sst_aquap5,sst_aquap6,sst_aquap7,sst_aquap8,sst_aquap9,sst_aquap10,sst_aquapfile,sst_aquap_constant,som,som_aquap,iaf,cplhist,multilev,multilev_dom General method that operates on the data for a given docn_mode. ==> dataMode = "sstdata" @@ -107,6 +109,8 @@ sst_aquap_file sst_aquap_constant cplhist + multilev_dom + multilev diff --git a/docn/cime_config/stream_definition_docn.xml b/docn/cime_config/stream_definition_docn.xml index 0dd4baf54..4e4dd6655 100644 --- a/docn/cime_config/stream_definition_docn.xml +++ b/docn/cime_config/stream_definition_docn.xml @@ -206,6 +206,104 @@ 1.e30 - single + single + + + + + $DIN_LOC_ROOT/share/meshes/gx1v7_151008_ESMFmesh.nc + + + /glade/collections/cmip/CMIP6/OMIP/NCAR/CESM2/omip2/r1i1p1f1/Omon/thetao/gn/v20190802/thetao_Omon_CESM2_omip2_r1i1p1f1_gn_024501-030512.nc + + + thetao So_t_depth + + lev + + bilinear + + null + 1 + 245 + 245 + 0 + + linear + + + cycle + + + 1.5 + + single + + + + + $DIN_LOC_ROOT/share/meshes/gx1v7_151008_ESMFmesh.nc + + + /glade/collections/cmip/CMIP6/OMIP/NCAR/CESM2/omip2/r1i1p1f1/Omon/so/gn/v20190802/so_Omon_CESM2_omip2_r1i1p1f1_gn_024501-030512.nc + + + so So_s_depth + + lev + + bilinear + + null + 1 + 245 + 245 + 0 + + linear + + + cycle + + + 1.5 + + single + + + + + $DIN_LOC_ROOT/share/meshes/tnx1v4_20170601_cdf5_ESMFmesh.nc + + + $DIN_LOC_ROOT/ocn/docn7/MULTILEV/N1850frc2G_f09_tn14_gl4_SMB1_celo.blom.hm.2300.nc + $DIN_LOC_ROOT/ocn/docn7/MULTILEV/N1850frc2G_f09_tn14_gl4_SMB1_celo.blom.hm.2301.nc + $DIN_LOC_ROOT/ocn/docn7/MULTILEV/N1850frc2G_f09_tn14_gl4_SMB1_celo.blom.hm.2302.nc + $DIN_LOC_ROOT/ocn/docn7/MULTILEV/N1850frc2G_f09_tn14_gl4_SMB1_celo.blom.hm.2303.nc + + + templvl So_t_depth + salnlvl So_s_depth + + depth + + bilinear + + null + 1 + 2300 + 2303 + 0 + + linear + + + cycle + + + 1.5 + + single + diff --git a/docn/cime_config/testdefs/testlist_docn.xml b/docn/cime_config/testdefs/testlist_docn.xml index 5dea00d48..5c279cabf 100644 --- a/docn/cime_config/testdefs/testlist_docn.xml +++ b/docn/cime_config/testdefs/testlist_docn.xml @@ -1,29 +1,32 @@ - + - + + - + - + + - + - + + - + diff --git a/docn/docn_datamode_multilev_dom_mod.F90 b/docn/docn_datamode_multilev_dom_mod.F90 new file mode 100644 index 000000000..7cd7fd9eb --- /dev/null +++ b/docn/docn_datamode_multilev_dom_mod.F90 @@ -0,0 +1,264 @@ +module docn_datamode_multilev_dom_mod + + use ESMF , only : ESMF_State, ESMF_LOGMSG_INFO, ESMF_LogWrite, ESMF_SUCCESS + use NUOPC , only : NUOPC_Advertise + use shr_kind_mod , only : r8=>shr_kind_r8, i8=>shr_kind_i8, cl=>shr_kind_cl, cs=>shr_kind_cs + use shr_const_mod , only : shr_const_TkFrz, shr_const_pi, shr_const_ocn_ref_sal, shr_const_spval + use shr_sys_mod , only : shr_sys_abort + use dshr_methods_mod , only : dshr_state_getfldptr, dshr_fldbun_getfldptr, chkerr + use dshr_fldlist_mod , only : fldlist_type, dshr_fldlist_add + use dshr_mod , only : dshr_restart_read, dshr_restart_write + use dshr_strdata_mod , only : shr_strdata_get_stream_pointer, shr_strdata_type, shr_strdata_get_stream_count + + implicit none + private ! except + + public :: docn_datamode_multilev_dom_advertise + public :: docn_datamode_multilev_dom_init_pointers + public :: docn_datamode_multilev_dom_advance + public :: docn_datamode_multilev_dom_restart_read + public :: docn_datamode_multilev_dom_restart_write + + ! pointers to export fields + real(r8), pointer :: So_omask(:) => null() ! real ocean fraction sent to mediator + real(r8), pointer :: So_t_depth(:,:) => null() + real(r8), pointer :: So_s_depth(:,:) => null() + real(r8), pointer :: So_t(:) => null() + real(r8), pointer :: So_u(:) => null() + real(r8), pointer :: So_v(:) => null() + real(r8), pointer :: So_s(:) => null() + + ! pointers to stream fields + real(r8), pointer :: stream_So_t_depth(:,:) => null() + real(r8), pointer :: stream_So_s_depth(:,:) => null() + + integer, parameter :: nlev_export = 30 + real(r8) :: vertical_levels(nlev_export) = (/ & + 30., 90., 150., 210., 270., 330., 390., 450., 510., 570., & + 630., 690., 750., 810., 870., 930., 990., 1050., 1110., 1170., & + 1230., 1290., 1350., 1410., 1470., 1530., 1590., 1650., 1710., 1770. /) + + real(r8) , parameter :: tkfrz = shr_const_tkfrz ! freezing point, fresh water (kelvin) + real(r8) , parameter :: ocnsalt = shr_const_ocn_ref_sal ! ocean reference salinity + + ! constants + character(*) , parameter :: nullstr = 'null' + character(*) , parameter :: rpfile = 'rpointer.ocn' + character(*) , parameter :: u_FILE_u = & + __FILE__ + +!=============================================================================== +contains +!=============================================================================== + + subroutine docn_datamode_multilev_dom_advertise(exportState, fldsexport, flds_scalar_name, rc) + + ! input/output variables + type(esmf_State) , intent(inout) :: exportState + type(fldlist_type) , pointer :: fldsexport + character(len=*) , intent(in) :: flds_scalar_name + integer , intent(out) :: rc + + ! local variables + type(fldlist_type), pointer :: fldList + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! Advertise export fields + call dshr_fldList_add(fldsExport, trim(flds_scalar_name)) + call dshr_fldList_add(fldsExport, 'So_omask') + call dshr_fldList_add(fldsExport, 'So_t_depth', ungridded_lbound=1, ungridded_ubound=nlev_export) + call dshr_fldList_add(fldsExport, 'So_s_depth', ungridded_lbound=1, ungridded_ubound=nlev_export) + call dshr_fldList_add(fldsExport, 'So_t' ) + call dshr_fldList_add(fldsExport, 'So_s' ) + call dshr_fldList_add(fldsExport, 'So_u' ) + call dshr_fldList_add(fldsExport, 'So_v' ) + + fldlist => fldsExport ! the head of the linked list + do while (associated(fldlist)) + call NUOPC_Advertise(exportState, standardName=fldlist%stdname, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_LogWrite('(docn_comp_advertise): Fr_ocn'//trim(fldList%stdname), ESMF_LOGMSG_INFO) + fldList => fldList%next + enddo + + end subroutine docn_datamode_multilev_dom_advertise + + !=============================================================================== + subroutine docn_datamode_multilev_dom_init_pointers(exportState, sdat, ocn_fraction, rc) + + ! input/output variables + type(ESMF_State) , intent(inout) :: exportState + type(shr_strdata_type) , intent(in) :: sdat + real(r8) , intent(in) :: ocn_fraction(:) + integer , intent(out) :: rc + + ! local variables + character(len=*), parameter :: subname='(docn_init_pointers): ' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! initialize pointers to stream fields + ! this has the full set of leveles in the stream data + call shr_strdata_get_stream_pointer( sdat, 'So_t_depth', stream_So_t_depth, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call shr_strdata_get_stream_pointer( sdat, 'So_s_depth', stream_So_s_depth, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + call dshr_state_getfldptr(exportState, 'So_t', fldptr1=So_t, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'So_s', fldptr1=So_s, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'So_u', fldptr1=So_u, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'So_v', fldptr1=So_v, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! initialize pointers to export fields + ! the export state has only nlev_export levels + call dshr_state_getfldptr(exportState, 'So_omask' , fldptr1=So_omask , rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'So_t_depth' , fldptr2=So_t_depth , rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'So_s_depth' , fldptr2=So_s_depth , rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Initialize export state pointers to non-zero + So_t_depth(:,:) = shr_const_TkFrz + So_s_depth(:,:) = shr_const_ocn_ref_sal + + So_u(:) = 0.0_r8 + So_v(:) = 0.0_r8 + So_s(:) = ocnsalt + So_t(:) = TkFrz + + ! Set export state ocean fraction (So_omask) + So_omask(:) = ocn_fraction(:) + + end subroutine docn_datamode_multilev_dom_init_pointers + + !=============================================================================== + subroutine docn_datamode_multilev_dom_advance(sdat, logunit, mainproc, rc) + + ! input/output variables + type(shr_strdata_type) , intent(in) :: sdat + integer , intent(in) :: logunit + logical , intent(in) :: mainproc + integer , intent(out) :: rc + + ! local variables + integer :: i,ki,ko + integer :: nstreams + integer :: nlev_stream + integer :: stream_index + logical :: level_found + real(r8) :: factor + real(r8), allocatable :: stream_vlevs(:) + logical :: first_time = .true. + character(len=*), parameter :: subname='(docn_datamode_multilev_dom): ' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + So_t(:) = So_t(:) + TkFrz + + ! Determine number of vertical levels for multi level stream + nstreams = shr_strdata_get_stream_count(sdat) + nlev_stream = 0 + do stream_index = 1,nstreams + nlev_stream = sdat%pstrm(stream_index)%stream_nlev + if (nlev_stream > 1) exit + end do + if (nlev_stream == 0) then + call shr_sys_abort(trim(subname)//" could not find vertical levels greater than 0") + end if + allocate(stream_vlevs(nlev_stream)) + stream_vlevs(:) = sdat%pstrm(stream_index)%stream_vlevs(:) + + do ko = 1,nlev_export + level_found = .false. + do ki = 1,nlev_stream-1 + if (vertical_levels(ko) > stream_vlevs(ki) .and. vertical_levels(ko) <= stream_vlevs(ki+1)) then + if (mainproc .and. first_time) then + write(logunit,'(a,3(i5,2x),3(f13.5,2x))') & + 'vertical interpolation: ki,ko,ki+1,lev(ki),lev(ko),lev(ki+1) = ',& + ki,ko,ki+1,stream_vlevs(ki), vertical_levels(ko), stream_vlevs(ki+1) + end if + level_found = .true. + do i = 1,size(So_omask) + if (So_omask(i) == 0.) then + So_t_depth(ko,i) = shr_const_spval + So_s_depth(ko,i) = shr_const_spval + else + ! Assume input T forcing is in degrees C + if (stream_So_t_depth(ki+1,i) > 1.e10) then + if (stream_So_t_depth(ki,i) > 1.e10) then + So_t_depth(ko,i) = shr_const_spval + So_s_depth(ko,i) = shr_const_spval + else + So_t_depth(ko,i) = stream_So_t_depth(ki,i) + shr_const_tkfrz + So_s_depth(ko,i) = stream_So_s_depth(ki,i) + end if + else + factor = (stream_So_t_depth(ki+1,i)-stream_So_t_depth(ki,i))/(stream_vlevs(ki+1)-stream_vlevs(ki)) + So_t_depth(ko,i) = stream_So_t_depth(ki,i) + (vertical_levels(ko)-stream_vlevs(ki))*factor + So_t_depth(ko,i) = So_t_depth(ko,i) + shr_const_tkfrz + + factor = (stream_So_s_depth(ki+1,i)-stream_So_s_depth(ki,i))/(stream_vlevs(ki+1)-stream_vlevs(ki)) + So_s_depth(ko,i) = stream_So_s_depth(ki,i) + (vertical_levels(ko)-stream_vlevs(ki))*factor + end if + end if + end do + end if + end do + if (.not. level_found) then + call shr_sys_abort(trim(subname)//" could not find level bounds for vertical interpolation") + end if + end do + + first_time = .false. + + end subroutine docn_datamode_multilev_dom_advance + + !=============================================================================== + subroutine docn_datamode_multilev_dom_restart_write(case_name, inst_suffix, ymd, tod, & + logunit, my_task, sdat) + + ! write restart file + + ! input/output variables + character(len=*) , intent(in) :: case_name + character(len=*) , intent(in) :: inst_suffix + integer , intent(in) :: ymd ! model date + integer , intent(in) :: tod ! model sec into model date + integer , intent(in) :: logunit + integer , intent(in) :: my_task + type(shr_strdata_type) , intent(inout) :: sdat + !------------------------------------------------------------------------------- + + call dshr_restart_write(rpfile, case_name, 'docn', inst_suffix, ymd, tod, & + logunit, my_task, sdat) + + end subroutine docn_datamode_multilev_dom_restart_write + + !=============================================================================== + subroutine docn_datamode_multilev_dom_restart_read(rest_filem, inst_suffix, logunit, my_task, mpicom, sdat) + + ! read restart file + + ! input/output arguments + character(len=*) , intent(inout) :: rest_filem + character(len=*) , intent(in) :: inst_suffix + integer , intent(in) :: logunit + integer , intent(in) :: my_task + integer , intent(in) :: mpicom + type(shr_strdata_type) , intent(inout) :: sdat + !------------------------------------------------------------------------------- + + call dshr_restart_read(rest_filem, rpfile, inst_suffix, nullstr, logunit, my_task, mpicom, sdat) + + end subroutine docn_datamode_multilev_dom_restart_read + +end module docn_datamode_multilev_dom_mod diff --git a/docn/docn_datamode_multilev_mod.F90 b/docn/docn_datamode_multilev_mod.F90 new file mode 100644 index 000000000..897797f3f --- /dev/null +++ b/docn/docn_datamode_multilev_mod.F90 @@ -0,0 +1,229 @@ +module docn_datamode_multilev_mod + use ESMF , only : ESMF_State, ESMF_LOGMSG_INFO, ESMF_LogWrite, ESMF_SUCCESS + use NUOPC , only : NUOPC_Advertise + use shr_kind_mod , only : r8=>shr_kind_r8, i8=>shr_kind_i8, cl=>shr_kind_cl, cs=>shr_kind_cs + use shr_const_mod , only : shr_const_TkFrz, shr_const_pi, shr_const_ocn_ref_sal, shr_const_spval + use shr_sys_mod , only : shr_sys_abort + use dshr_methods_mod , only : dshr_state_getfldptr, dshr_fldbun_getfldptr, chkerr + use dshr_fldlist_mod , only : fldlist_type, dshr_fldlist_add + use dshr_mod , only : dshr_restart_read, dshr_restart_write + use dshr_strdata_mod , only : shr_strdata_get_stream_pointer, shr_strdata_type + + implicit none + private ! except + + public :: docn_datamode_multilev_advertise + public :: docn_datamode_multilev_init_pointers + public :: docn_datamode_multilev_advance + public :: docn_datamode_multilev_restart_read + public :: docn_datamode_multilev_restart_write + + ! pointers to export fields + real(r8), pointer :: So_omask(:) => null() ! real ocean fraction sent to mediator + real(r8), pointer :: So_t_depth(:,:) => null() + real(r8), pointer :: So_s_depth(:,:) => null() + + ! pointers to stream fields + real(r8), pointer :: stream_So_t_depth(:,:) => null() + real(r8), pointer :: stream_So_s_depth(:,:) => null() + + integer, parameter :: nlev_export = 30 + real(r8) :: vertical_levels(nlev_export) = (/ & + 30., 90., 150., 210., 270., 330., 390., 450., 510., 570., & + 630., 690., 750., 810., 870., 930., 990., 1050., 1110., 1170., & + 1230., 1290., 1350., 1410., 1470., 1530., 1590., 1650., 1710., 1770. /) + + ! constants + character(*) , parameter :: nullstr = 'null' + character(*) , parameter :: rpfile = 'rpointer.ocn' + character(*) , parameter :: u_FILE_u = & + __FILE__ + +!=============================================================================== +contains +!=============================================================================== + + subroutine docn_datamode_multilev_advertise(exportState, fldsexport, flds_scalar_name, rc) + + ! input/output variables + type(esmf_State) , intent(inout) :: exportState + type(fldlist_type) , pointer :: fldsexport + character(len=*) , intent(in) :: flds_scalar_name + integer , intent(out) :: rc + + ! local variables + type(fldlist_type), pointer :: fldList + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! Advertise export fields + call dshr_fldList_add(fldsExport, trim(flds_scalar_name)) + call dshr_fldList_add(fldsExport, 'So_omask') + call dshr_fldList_add(fldsExport, 'So_t_depth', ungridded_lbound=1, ungridded_ubound=nlev_export) + call dshr_fldList_add(fldsExport, 'So_s_depth', ungridded_lbound=1, ungridded_ubound=nlev_export) + + fldlist => fldsExport ! the head of the linked list + do while (associated(fldlist)) + call NUOPC_Advertise(exportState, standardName=fldlist%stdname, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_LogWrite('(docn_comp_advertise): Fr_ocn'//trim(fldList%stdname), ESMF_LOGMSG_INFO) + fldList => fldList%next + enddo + + end subroutine docn_datamode_multilev_advertise + + !=============================================================================== + subroutine docn_datamode_multilev_init_pointers(exportState, sdat, ocn_fraction, rc) + + ! input/output variables + type(ESMF_State) , intent(inout) :: exportState + type(shr_strdata_type) , intent(in) :: sdat + real(r8) , intent(in) :: ocn_fraction(:) + integer , intent(out) :: rc + + ! local variables + character(len=*), parameter :: subname='(docn_init_pointers): ' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! initialize pointers to stream fields + ! this has the full set of leveles in the stream data + call shr_strdata_get_stream_pointer( sdat, 'So_t_depth', stream_So_t_depth, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call shr_strdata_get_stream_pointer( sdat, 'So_s_depth', stream_So_s_depth, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! initialize pointers to export fields + ! the export state has only nlev_export levels + call dshr_state_getfldptr(exportState, 'So_omask' , fldptr1=So_omask , rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'So_t_depth' , fldptr2=So_t_depth , rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'So_s_depth' , fldptr2=So_s_depth , rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Initialize export state pointers to non-zero + So_t_depth(:,:) = shr_const_TkFrz + So_s_depth(:,:) = shr_const_ocn_ref_sal + + ! Set export state ocean fraction (So_omask) + So_omask(:) = ocn_fraction(:) + + end subroutine docn_datamode_multilev_init_pointers + + !=============================================================================== + subroutine docn_datamode_multilev_advance(sdat, logunit, mainproc, rc) + + ! input/output variables + type(shr_strdata_type) , intent(in) :: sdat + integer , intent(in) :: logunit + logical , intent(in) :: mainproc + integer , intent(out) :: rc + + ! local variables + integer :: i,ki,ko + integer :: nlev_stream + integer :: stream_index + logical :: level_found + real(r8) :: factor + real(r8), allocatable :: stream_vlevs(:) + logical :: first_time = .true. + character(len=*), parameter :: subname='(docn_datamode_multilev): ' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! For now assume that all the streams have the same vertical levels + stream_index = 1 + + nlev_stream = sdat%pstrm(stream_index)%stream_nlev + allocate(stream_vlevs(nlev_stream)) + stream_vlevs(:) = sdat%pstrm(stream_index)%stream_vlevs(:) + + do ko = 1,nlev_export + level_found = .false. + do ki = 1,nlev_stream-1 + if (vertical_levels(ko) > stream_vlevs(ki) .and. vertical_levels(ko) <= stream_vlevs(ki+1)) then + if (mainproc .and. first_time) then + write(logunit,'(a,3(i5,2x),3(f13.5,2x))') & + 'vertical interpolation: ki,ko,ki+1,lev(ki),lev(ko),lev(ki+1) = ',& + ki,ko,ki+1,stream_vlevs(ki), vertical_levels(ko), stream_vlevs(ki+1) + end if + level_found = .true. + do i = 1,size(So_omask) + if (So_omask(i) == 0.) then + So_t_depth(ko,i) = shr_const_spval + So_s_depth(ko,i) = shr_const_spval + else + ! Assume input T forcing is in degrees C + if (stream_So_t_depth(ki+1,i) > 1.e10) then + if (stream_So_t_depth(ki,i) > 1.e10) then + So_t_depth(ko,i) = shr_const_spval + So_s_depth(ko,i) = shr_const_spval + else + So_t_depth(ko,i) = stream_So_t_depth(ki,i) + shr_const_tkfrz + So_s_depth(ko,i) = stream_So_s_depth(ki,i) + end if + else + factor = (stream_So_t_depth(ki+1,i)-stream_So_t_depth(ki,i))/(stream_vlevs(ki+1)-stream_vlevs(ki)) + So_t_depth(ko,i) = stream_So_t_depth(ki,i) + (vertical_levels(ko)-stream_vlevs(ki))*factor + So_t_depth(ko,i) = So_t_depth(ko,i) + shr_const_tkfrz + + factor = (stream_So_s_depth(ki+1,i)-stream_So_s_depth(ki,i))/(stream_vlevs(ki+1)-stream_vlevs(ki)) + So_s_depth(ko,i) = stream_So_s_depth(ki,i) + (vertical_levels(ko)-stream_vlevs(ki))*factor + end if + end if + end do + end if + end do + if (.not. level_found) then + call shr_sys_abort("ERROR: could not find level bounds for vertical interpolation") + end if + end do + + first_time = .false. + + end subroutine docn_datamode_multilev_advance + + !=============================================================================== + subroutine docn_datamode_multilev_restart_write(case_name, inst_suffix, ymd, tod, & + logunit, my_task, sdat) + + ! write restart file + + ! input/output variables + character(len=*) , intent(in) :: case_name + character(len=*) , intent(in) :: inst_suffix + integer , intent(in) :: ymd ! model date + integer , intent(in) :: tod ! model sec into model date + integer , intent(in) :: logunit + integer , intent(in) :: my_task + type(shr_strdata_type) , intent(inout) :: sdat + !------------------------------------------------------------------------------- + + call dshr_restart_write(rpfile, case_name, 'docn', inst_suffix, ymd, tod, & + logunit, my_task, sdat) + + end subroutine docn_datamode_multilev_restart_write + + !=============================================================================== + subroutine docn_datamode_multilev_restart_read(rest_filem, inst_suffix, logunit, my_task, mpicom, sdat) + + ! read restart file + + ! input/output arguments + character(len=*) , intent(inout) :: rest_filem + character(len=*) , intent(in) :: inst_suffix + integer , intent(in) :: logunit + integer , intent(in) :: my_task + integer , intent(in) :: mpicom + type(shr_strdata_type) , intent(inout) :: sdat + !------------------------------------------------------------------------------- + + call dshr_restart_read(rest_filem, rpfile, inst_suffix, nullstr, logunit, my_task, mpicom, sdat) + + end subroutine docn_datamode_multilev_restart_read + +end module docn_datamode_multilev_mod diff --git a/docn/docn_datamode_som_mod.F90 b/docn/docn_datamode_som_mod.F90 index b91ce134d..672c58d4d 100644 --- a/docn/docn_datamode_som_mod.F90 +++ b/docn/docn_datamode_som_mod.F90 @@ -34,6 +34,7 @@ module docn_datamode_som_mod real(r8), pointer :: So_v(:) => null() real(r8), pointer :: So_dhdx(:) => null() real(r8), pointer :: So_dhdy(:) => null() + real(r8), pointer :: So_bldepth(:) => null() real(r8), pointer :: Fioo_q(:) => null() real(r8), pointer :: So_fswpen(:) => null() @@ -98,6 +99,7 @@ subroutine docn_datamode_som_advertise(importState, exportState, fldsimport, fld call dshr_fldList_add(fldsExport, 'So_v' ) call dshr_fldList_add(fldsExport, 'So_dhdx' ) call dshr_fldList_add(fldsExport, 'So_dhdy' ) + call dshr_fldList_add(fldsExport, 'So_bldepth' ) call dshr_fldList_add(fldsExport, 'Fioo_q' ) call dshr_fldList_add(fldsExport, 'So_fswpen' ) @@ -189,6 +191,8 @@ subroutine docn_datamode_som_init_pointers(importState, exportState, sdat, ocn_f if (chkerr(rc,__LINE__,u_FILE_u)) return call dshr_state_getfldptr(exportState, 'So_dhdy' , fldptr1=So_dhdy , rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'So_bldepth' , fldptr1=So_bldepth , rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return call dshr_state_getfldptr(exportState, 'Fioo_q' , fldptr1=Fioo_q , rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return @@ -289,6 +293,7 @@ subroutine docn_datamode_som_advance(importState, exportState, clock, restart_re ! save somtp to restart file somtp(n) = So_t(n) + So_bldepth(n) = strm_h(n) endif end do deallocate(tfreeze) diff --git a/docn/ocn_comp_nuopc.F90 b/docn/ocn_comp_nuopc.F90 index 8df6f9566..1fb8b7940 100644 --- a/docn/ocn_comp_nuopc.F90 +++ b/docn/ocn_comp_nuopc.F90 @@ -26,7 +26,7 @@ module cdeps_docn_comp use shr_kind_mod , only : r8=>shr_kind_r8, i8=>shr_kind_i8, cl=>shr_kind_cl, cs=>shr_kind_cs use shr_sys_mod , only : shr_sys_abort use shr_cal_mod , only : shr_cal_ymd2date - use shr_log_mod , only : shr_log_setLogUnit + use shr_log_mod , only : shr_log_setLogUnit use dshr_methods_mod , only : dshr_state_diagnose, chkerr, memcheck use dshr_strdata_mod , only : shr_strdata_type, shr_strdata_advance, shr_strdata_init_from_config use dshr_mod , only : dshr_model_initphase, dshr_init, dshr_mesh_init @@ -58,7 +58,17 @@ module cdeps_docn_comp use docn_datamode_cplhist_mod , only : docn_datamode_cplhist_advance use docn_datamode_cplhist_mod , only : docn_datamode_cplhist_restart_read use docn_datamode_cplhist_mod , only : docn_datamode_cplhist_restart_write - use docn_import_data_mod , only : docn_import_data_advertise + use docn_datamode_multilev_mod , only : docn_datamode_multilev_advertise + use docn_datamode_multilev_mod , only : docn_datamode_multilev_init_pointers + use docn_datamode_multilev_mod , only : docn_datamode_multilev_advance + use docn_datamode_multilev_mod , only : docn_datamode_multilev_restart_read + use docn_datamode_multilev_mod , only : docn_datamode_multilev_restart_write + use docn_datamode_multilev_dom_mod, only : docn_datamode_multilev_dom_advertise + use docn_datamode_multilev_dom_mod, only : docn_datamode_multilev_dom_init_pointers + use docn_datamode_multilev_dom_mod, only : docn_datamode_multilev_dom_advance + use docn_datamode_multilev_dom_mod, only : docn_datamode_multilev_dom_restart_read + use docn_datamode_multilev_dom_mod, only : docn_datamode_multilev_dom_restart_write + use docn_import_data_mod , only : docn_import_data_advertise implicit none private ! except @@ -179,6 +189,7 @@ end subroutine SetServices !=============================================================================== subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) use shr_nl_mod, only: shr_nl_find_group_name + ! input/output variables type(ESMF_GridComp) :: gcomp type(ESMF_State) :: importState, exportState @@ -212,7 +223,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ! Obtain flds_scalar values, mpi values, multi-instance values and ! set logunit and set shr logging to my log file - call dshr_init(gcomp, 'OCN', sdat, mpicom, my_task, inst_index, inst_suffix, & + call dshr_init(gcomp, 'OCN', mpicom, my_task, inst_index, inst_suffix, & flds_scalar_name, flds_scalar_num, flds_scalar_index_nx, flds_scalar_index_ny, logunit, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -300,7 +311,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) trim(datamode) == 'som_aquap' .or. & ! read stream, needs import data trim(datamode) == 'cplhist' .or. & ! read stream, needs import data trim(datamode) == 'sst_aquap_analytic' .or. & ! analytic, no streams, import or export data - trim(datamode) == 'sst_aquap_constant' ) then ! analytic, no streams, import or export data + trim(datamode) == 'sst_aquap_constant' .or. & ! analytic, no streams, import or export data + trim(datamode) == 'multilev' .or. & ! multilevel ocean input + trim(datamode) == 'multilev_dom') then ! multilevel ocean input and sst export ! success do nothing else call shr_sys_abort(' ERROR illegal docn datamode = '//trim(datamode)) @@ -323,6 +336,12 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) else if (trim(datamode) == 'cplhist') then call docn_datamode_cplhist_advertise(exportState, fldsExport, flds_scalar_name, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + else if (trim(datamode) == 'multilev') then + call docn_datamode_multilev_advertise(exportState, fldsExport, flds_scalar_name, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else if (trim(datamode) == 'multilev_dom') then + call docn_datamode_multilev_dom_advertise(exportState, fldsExport, flds_scalar_name, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return end if if (trim(import_data_fields) /= 'none') then @@ -550,6 +569,12 @@ subroutine docn_comp_run(importState, exportState, clock, target_ymd, target_tod case('cplhist') call docn_datamode_cplhist_init_pointers(exportState, model_frac, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + case('multilev') + call docn_datamode_multilev_init_pointers(exportState, sdat, model_frac, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + case('multilev_dom') + call docn_datamode_multilev_dom_init_pointers(exportState, sdat, model_frac, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return end select ! Read restart if needed @@ -607,6 +632,12 @@ subroutine docn_comp_run(importState, exportState, clock, target_ymd, target_tod case('cplhist') call docn_datamode_cplhist_advance(sst_constant_value=sst_constant_value, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + case('multilev') + call docn_datamode_multilev_advance(sdat, logunit, mainproc, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + case('multilev_dom') + call docn_datamode_multilev_dom_advance(sdat, logunit, mainproc, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return end select ! Write restarts if needed (no restarts for aquaplanet analytic or aquaplanet input file) @@ -650,8 +681,10 @@ subroutine docn_init_dfields(importState, exportState, rc) ! local variables integer :: n integer :: fieldcount + integer :: dimcount type(ESMF_Field) :: lfield character(ESMF_MAXSTR) ,pointer :: lfieldnamelist(:) + character(ESMF_MAXSTR) :: fieldname(1) character(*), parameter :: subName = "(docn_init_dfields) " !------------------------------------------------------------------------------- @@ -668,9 +701,18 @@ subroutine docn_init_dfields(importState, exportState, rc) call ESMF_StateGet(exportState, itemName=trim(lfieldNameList(n)), field=lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (trim(lfieldnamelist(n)) /= flds_scalar_name) then - call dshr_dfield_add( dfields, sdat, trim(lfieldnamelist(n)), trim(lfieldnamelist(n)), exportState, & - logunit, mainproc, rc) + call ESMF_FieldGet(lfield, dimcount=dimCount, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return + if (dimcount == 2) then + fieldname(1) = trim(lfieldnamelist(n)) + call dshr_dfield_add( dfields, sdat, trim(lfieldnamelist(n)), fieldname, exportState, & + logunit, mainproc, rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + else + call dshr_dfield_add( dfields, sdat, trim(lfieldnamelist(n)), trim(lfieldnamelist(n)), exportState, & + logunit, mainproc, rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + endif end if end do end subroutine docn_init_dfields diff --git a/drof/cime_config/testdefs/testlist_drof.xml b/drof/cime_config/testdefs/testlist_drof.xml index 0515f0fc6..6e423fabe 100644 --- a/drof/cime_config/testdefs/testlist_drof.xml +++ b/drof/cime_config/testdefs/testlist_drof.xml @@ -1,9 +1,10 @@ - + - + + diff --git a/drof/rof_comp_nuopc.F90 b/drof/rof_comp_nuopc.F90 index f6e9e7cd1..0808507a6 100644 --- a/drof/rof_comp_nuopc.F90 +++ b/drof/rof_comp_nuopc.F90 @@ -27,7 +27,7 @@ module cdeps_drof_comp use shr_const_mod , only : SHR_CONST_SPVAL use shr_sys_mod , only : shr_sys_abort use shr_cal_mod , only : shr_cal_ymd2date - use shr_log_mod , only : shr_log_setLogUnit + use shr_log_mod , only : shr_log_setLogUnit use dshr_methods_mod , only : dshr_state_getfldptr, dshr_state_diagnose, chkerr, memcheck use dshr_strdata_mod , only : shr_strdata_type, shr_strdata_advance, shr_strdata_get_stream_domain use dshr_strdata_mod , only : shr_strdata_init_from_config @@ -183,7 +183,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ! Obtain flds_scalar values, mpi values, multi-instance values and ! set logunit and set shr logging to my log file - call dshr_init(gcomp, 'ROF', sdat, mpicom, my_task, inst_index, inst_suffix, & + call dshr_init(gcomp, 'ROF', mpicom, my_task, inst_index, inst_suffix, & flds_scalar_name, flds_scalar_num, flds_scalar_index_nx, flds_scalar_index_ny, & logunit, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return diff --git a/dshr/dshr_dfield_mod.F90 b/dshr/dshr_dfield_mod.F90 index 2c2862042..2f3b33d69 100644 --- a/dshr/dshr_dfield_mod.F90 +++ b/dshr/dshr_dfield_mod.F90 @@ -1,7 +1,7 @@ module dshr_dfield_mod use ESMF , only : ESMF_State, ESMF_FieldBundle, ESMF_MAXSTR, ESMF_SUCCESS - use ESMF , only : ESMF_FieldBundleGet, ESMF_ITEMORDER_ADDORDER, ESMF_Field + use ESMF , only : ESMF_FieldBundleGet, ESMF_ITEMORDER_ADDORDER, ESMF_Field, ESMF_FieldGet use shr_kind_mod , only : r8=>shr_kind_r8, cs=>shr_kind_cs, cl=>shr_kind_cl, cxx=>shr_kind_cxx use shr_sys_mod , only : shr_sys_abort use dshr_strdata_mod , only : shr_strdata_type, shr_strdata_get_stream_count, shr_strdata_get_stream_fieldbundle @@ -438,9 +438,12 @@ subroutine dshr_dfield_copy(dfields, sdat, rc) type(ESMF_field) :: lfield type(dfield_type), pointer :: dfield real(r8), pointer :: data1d(:) + real(r8), pointer :: data2d(:,:) integer :: nf integer :: fldbun_index integer :: stream_index + integer :: ungriddedUBound_output(1) + integer :: ungriddedCount !------------------------------------------------------------------------------- rc = ESMF_SUCCESS @@ -464,13 +467,22 @@ subroutine dshr_dfield_copy(dfields, sdat, rc) do nf = 1,size(dfield%stream_indices) stream_index = dfield%stream_indices(nf) fldbun_index = dfield%fldbun_indices(nf) - if(stream_index > 0) then + if (stream_index > 0) then fldbun_model = shr_strdata_get_stream_fieldbundle(sdat, stream_index, 'model') call dshr_fldbun_getfieldn(fldbun_model, fldbun_index, lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - call dshr_field_getfldptr(lfield, fldptr1=data1d, rc=rc) + call ESMF_FieldGet(lfield, ungriddedUBound=ungriddedUBound_output, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - dfield%state_data2d(nf,:) = data1d(:) + ungriddedCount = ungriddedUBound_output(1) + if (ungriddedCount > 0) then + call dshr_field_getfldptr(lfield, fldptr2=data2d, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + dfield%state_data2d(:,:) = data2d(:,:) + else + call dshr_field_getfldptr(lfield, fldptr1=data1d, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + dfield%state_data2d(nf,:) = data1d(:) + end if endif end do end if diff --git a/dshr/dshr_mod.F90 b/dshr/dshr_mod.F90 index f253adc05..faa346b68 100644 --- a/dshr/dshr_mod.F90 +++ b/dshr/dshr_mod.F90 @@ -60,9 +60,9 @@ module dshr_mod public :: dshr_state_setscalar public :: dshr_orbital_update public :: dshr_orbital_init + public :: dshr_alarm_init ! initialize alarms private :: dshr_mesh_create_scol ! create mesh for single column mode - private :: dshr_alarm_init ! initialize alarms private :: dshr_time_init ! initialize time ! Note that gridTofieldMap = 2, therefore the ungridded dimension is innermost @@ -111,14 +111,13 @@ subroutine dshr_model_initphase(gcomp, importState, exportState, clock, rc) end subroutine dshr_model_initphase !=============================================================================== - subroutine dshr_init(gcomp, compname, sdat, mpicom, my_task, inst_index, inst_suffix, & + subroutine dshr_init(gcomp, compname, mpicom, my_task, inst_index, inst_suffix, & flds_scalar_name, flds_scalar_num, flds_scalar_index_nx, flds_scalar_index_ny, logunit, rc) #ifdef CESMCOUPLED use nuopc_shr_methods, only : set_component_logging #endif ! input/output variables type(ESMF_GridComp) :: gcomp - type(shr_strdata_type), intent(in) :: sdat ! No longer used character(len=*) , intent(in) :: compname !e.g. ATM, OCN, ... integer , intent(inout) :: mpicom integer , intent(out) :: my_task @@ -246,7 +245,7 @@ subroutine dshr_mesh_init(gcomp, sdat, nullstr, logunit, compname, model_nxg, mo ! input/output variables type(ESMF_GridComp) , intent(inout) :: gcomp - type(shr_strdata_type) , intent(inout) :: sdat + type(shr_strdata_type) , intent(inout) :: sdat integer , intent(in) :: logunit character(len=*) , intent(in) :: compname !e.g. ATM, OCN, ... character(len=*) , intent(in) :: nullstr diff --git a/dwav/cime_config/testdefs/testlist_dwav.xml b/dwav/cime_config/testdefs/testlist_dwav.xml index 0fed61eb2..1864398ad 100644 --- a/dwav/cime_config/testdefs/testlist_dwav.xml +++ b/dwav/cime_config/testdefs/testlist_dwav.xml @@ -1,9 +1,10 @@ - + - + + diff --git a/dwav/wav_comp_nuopc.F90 b/dwav/wav_comp_nuopc.F90 index 351e7c82e..ffabd3a49 100644 --- a/dwav/wav_comp_nuopc.F90 +++ b/dwav/wav_comp_nuopc.F90 @@ -181,7 +181,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ! Obtain flds_scalar values, mpi values, multi-instance values and ! set logunit and set shr logging to my log file - call dshr_init(gcomp, 'WAV', sdat, mpicom, my_task, inst_index, inst_suffix, & + call dshr_init(gcomp, 'WAV', mpicom, my_task, inst_index, inst_suffix, & flds_scalar_name, flds_scalar_num, flds_scalar_index_nx, flds_scalar_index_ny, & logunit, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return diff --git a/share/CMakeLists.txt b/share/CMakeLists.txt index 6d1c3607d..3a821304b 100644 --- a/share/CMakeLists.txt +++ b/share/CMakeLists.txt @@ -18,7 +18,8 @@ add_library(cdeps_share ${GenF90_SRCS} shr_log_mod.F90 shr_strconvert_mod.F90 shr_precip_mod.F90 - shr_string_mod.F90) + shr_string_mod.F90 + nuopc_shr_methods.F90) target_include_directories (cdeps_share PRIVATE ${CMAKE_CURRENT_SOURCE_DIR} ${ESMF_F90COMPILEPATHS} ${PIO_Fortran_INCLUDE_DIRS}) diff --git a/share/nuopc_shr_methods.F90 b/share/nuopc_shr_methods.F90 new file mode 100644 index 000000000..d65a9cc69 --- /dev/null +++ b/share/nuopc_shr_methods.F90 @@ -0,0 +1,778 @@ +module nuopc_shr_methods + + use ESMF , only : operator(<), operator(/=), operator(+) + use ESMF , only : operator(-), operator(*) , operator(>=) + use ESMF , only : operator(<=), operator(>), operator(==), MOD + use ESMF , only : ESMF_LOGERR_PASSTHRU, ESMF_LogFoundError, ESMF_LOGMSG_ERROR, ESMF_MAXSTR + use ESMF , only : ESMF_SUCCESS, ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_FAILURE + use ESMF , only : ESMF_State, ESMF_StateGet + use ESMF , only : ESMF_Field, ESMF_FieldGet + use ESMF , only : ESMF_GridComp, ESMF_GridCompGet, ESMF_GridCompSet + use ESMF , only : ESMF_GeomType_Flag, ESMF_FieldStatus_Flag + use ESMF , only : ESMF_Mesh, ESMF_MeshGet, ESMF_AlarmSet + use ESMF , only : ESMF_GEOMTYPE_MESH, ESMF_GEOMTYPE_GRID, ESMF_FIELDSTATUS_COMPLETE + use ESMF , only : ESMF_Clock, ESMF_ClockCreate, ESMF_ClockGet, ESMF_ClockSet + use ESMF , only : ESMF_ClockPrint, ESMF_ClockAdvance + use ESMF , only : ESMF_Alarm, ESMF_AlarmCreate, ESMF_AlarmGet, ESMF_AlarmSet + use ESMF , only : ESMF_Calendar, ESMF_CALKIND_NOLEAP, ESMF_CALKIND_GREGORIAN + use ESMF , only : ESMF_Time, ESMF_TimeGet, ESMF_TimeSet, ESMF_ClockGetAlarm + use ESMF , only : ESMF_TimeInterval, ESMF_TimeIntervalSet, ESMF_TimeIntervalGet + use ESMF , only : ESMF_VM, ESMF_VMGet, ESMF_VMBroadcast, ESMF_VMGetCurrent + use NUOPC , only : NUOPC_CompAttributeGet + use NUOPC_Model , only : NUOPC_ModelGet + use shr_kind_mod , only : r8 => shr_kind_r8, cl=>shr_kind_cl, cs=>shr_kind_cs + use shr_sys_mod , only : shr_sys_abort + use shr_log_mod , only : shr_log_setLogUnit + + implicit none + private + + public :: memcheck + public :: get_component_instance + public :: set_component_logging + public :: log_clock_advance + public :: state_getscalar + public :: state_setscalar + public :: state_diagnose + public :: alarmInit + public :: chkerr + + private :: timeInit + private :: field_getfldptr + + ! Module data + + ! Clock and alarm options shared with esm_time_mod along with dtime_driver which is initialized there. + ! Dtime_driver is needed here for setting alarm options which use the nstep option and is a module variable + ! to avoid requiring a change in all model caps. + character(len=*), public, parameter :: & + optNONE = "none" , & + optNever = "never" , & + optNSteps = "nstep" , & + optNSeconds = "nsecond" , & + optNMinutes = "nminute" , & + optNHours = "nhour" , & + optNDays = "nday" , & + optNMonths = "nmonth" , & + optNYears = "nyear" , & + optMonthly = "monthly" , & + optYearly = "yearly" , & + optDate = "date" , & + optEnd = "end" , & + optGLCCouplingPeriod = "glc_coupling_period" + + integer, public :: dtime_drv ! initialized in esm_time_mod.F90 + + integer, parameter :: SecPerDay = 86400 ! Seconds per day + integer, parameter :: memdebug_level=1 + character(len=1024) :: msgString + character(len=*), parameter :: u_FILE_u = & + __FILE__ + +!=============================================================================== +contains +!=============================================================================== + + subroutine memcheck(string, level, maintask) + + ! input/output variables + character(len=*) , intent(in) :: string + integer , intent(in) :: level + logical , intent(in) :: maintask + + ! local variables +#ifdef CESMCOUPLED + integer, external :: GPTLprint_memusage + integer :: ierr = 0 +#endif + !----------------------------------------------------------------------- + +#ifdef CESMCOUPLED + if ((maintask .and. memdebug_level > level) .or. memdebug_level > level+1) then + ierr = GPTLprint_memusage(string) + endif +#endif + + end subroutine memcheck + +!=============================================================================== + + subroutine get_component_instance(gcomp, inst_suffix, inst_index, rc) + + ! input/output variables + type(ESMF_GridComp) :: gcomp + character(len=*) , intent(out) :: inst_suffix + integer , intent(out) :: inst_index + integer , intent(out) :: rc + + ! local variables + logical :: isPresent + character(len=4) :: cvalue + !----------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + call NUOPC_CompAttributeGet(gcomp, name="inst_suffix", isPresent=isPresent, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + if (isPresent) then + call NUOPC_CompAttributeGet(gcomp, name="inst_suffix", value=inst_suffix, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + cvalue = inst_suffix(2:) + read(cvalue, *) inst_index + else + inst_suffix = "" + inst_index=1 + endif + + end subroutine get_component_instance + +!=============================================================================== + subroutine set_component_logging(gcomp, maintask, logunit, shrlogunit, rc) + use NUOPC, only: NUOPC_CompAttributeSet, NUOPC_CompAttributeAdd + ! input/output variables + type(ESMF_GridComp) :: gcomp + logical, intent(in) :: maintask + integer, intent(out) :: logunit + integer, intent(out) :: shrlogunit + integer, intent(out) :: rc + + ! local variables + character(len=CL) :: diro + character(len=CL) :: logfile + character(len=CL) :: inst_suffix + integer :: inst_index ! Not used here + integer :: n + character(len=CL) :: name + character(len=*), parameter :: subname = "("//__FILE__//": set_component_logging)" + !----------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + if (maintask) then + call NUOPC_CompAttributeGet(gcomp, name="diro", value=diro, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call NUOPC_CompAttributeGet(gcomp, name="logfile", value=logfile, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call get_component_instance(gcomp, inst_suffix, inst_index, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + ! Multiinstance logfile name needs a correction + if(len_trim(inst_suffix) > 0) then + n = index(logfile, '.') + logfile = logfile(1:n-1)//trim(inst_suffix)//logfile(n:) + endif + + open(newunit=logunit,file=trim(diro)//"/"//trim(logfile)) + + else + logUnit = 6 + endif + + call ESMF_GridCompGet(gcomp, name=name, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + call ESMF_LogWrite(trim(subname)//": setting logunit for component: "//trim(name), ESMF_LOGMSG_INFO) + call NUOPC_CompAttributeAdd(gcomp, (/"logunit"/), rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call NUOPC_CompAttributeSet(gcomp, "logunit", logunit, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call shr_log_setLogUnit (logunit) + ! Still need to set this return value + shrlogunit = logunit + call ESMF_LogWrite(trim(subname)//": done for component "//trim(name), ESMF_LOGMSG_INFO) + end subroutine set_component_logging + +!=============================================================================== + + subroutine log_clock_advance(clock, component, logunit, rc) + + ! input/output variables + type(ESMF_Clock) :: clock + character(len=*) , intent(in) :: component + integer , intent(in) :: logunit + integer , intent(out) :: rc + + ! local variables + character(len=CL) :: cvalue, prestring + !----------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + write(prestring, *) "------>Advancing ",trim(component)," from: " + call ESMF_ClockPrint(clock, options="currTime", unit=cvalue, preString=trim(prestring), rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + write(logunit, *) trim(cvalue) + + call ESMF_ClockPrint(clock, options="stopTime", unit=cvalue, & + preString="--------------------------------> to: ", rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + write(logunit, *) trim(cvalue) + + end subroutine log_clock_advance + +!=============================================================================== + + subroutine state_getscalar(state, scalar_id, scalar_value, flds_scalar_name, flds_scalar_num, rc) + + ! ---------------------------------------------- + ! Get scalar data from State for a particular name and broadcast it to all other pets + ! ---------------------------------------------- + + ! input/output variables + type(ESMF_State), intent(in) :: state + integer, intent(in) :: scalar_id + real(r8), intent(out) :: scalar_value + character(len=*), intent(in) :: flds_scalar_name + integer, intent(in) :: flds_scalar_num + integer, intent(inout) :: rc + + ! local variables + integer :: mytask + type(ESMF_VM) :: vm + type(ESMF_Field) :: field + real(r8), pointer :: farrayptr(:,:) + real(r8) :: tmp(1) + character(len=*), parameter :: subname='(state_getscalar)' + ! ---------------------------------------------- + + rc = ESMF_SUCCESS + + call ESMF_VMGetCurrent(vm, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + call ESMF_VMGet(vm, localPet=mytask, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + call ESMF_StateGet(State, itemName=trim(flds_scalar_name), field=field, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + if (mytask == 0) then + call ESMF_FieldGet(field, farrayPtr = farrayptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + if (scalar_id < 0 .or. scalar_id > flds_scalar_num) then + call ESMF_LogWrite(trim(subname)//": ERROR in scalar_id", ESMF_LOGMSG_INFO, line=__LINE__, file=u_FILE_u) + rc = ESMF_FAILURE + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=u_FILE_u)) return + endif + tmp(:) = farrayptr(scalar_id,:) + endif + call ESMF_VMBroadCast(vm, tmp, 1, 0, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + scalar_value = tmp(1) + + end subroutine state_getscalar + +!================================================================================ + + subroutine state_setscalar(scalar_value, scalar_id, State, flds_scalar_name, flds_scalar_num, rc) + + ! ---------------------------------------------- + ! Set scalar data from State for a particular name + ! ---------------------------------------------- + + ! input/output arguments + real(r8), intent(in) :: scalar_value + integer, intent(in) :: scalar_id + type(ESMF_State), intent(inout) :: State + character(len=*), intent(in) :: flds_scalar_name + integer, intent(in) :: flds_scalar_num + integer, intent(inout) :: rc + + ! local variables + integer :: mytask + type(ESMF_Field) :: lfield + type(ESMF_VM) :: vm + real(r8), pointer :: farrayptr(:,:) + character(len=*), parameter :: subname='(state_setscalar)' + ! ---------------------------------------------- + + rc = ESMF_SUCCESS + + call ESMF_VMGetCurrent(vm, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + call ESMF_VMGet(vm, localPet=mytask, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + call ESMF_StateGet(State, itemName=trim(flds_scalar_name), field=lfield, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + if (mytask == 0) then + call ESMF_FieldGet(lfield, farrayPtr = farrayptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + if (scalar_id < 0 .or. scalar_id > flds_scalar_num) then + call ESMF_LogWrite(trim(subname)//": ERROR in scalar_id", ESMF_LOGMSG_INFO) + rc = ESMF_FAILURE + return + endif + farrayptr(scalar_id,1) = scalar_value + endif + + end subroutine state_setscalar + +!=============================================================================== + + subroutine state_diagnose(State, string, rc) + + ! ---------------------------------------------- + ! Diagnose status of State + ! ---------------------------------------------- + + type(ESMF_State), intent(in) :: state + character(len=*), intent(in) :: string + integer , intent(out) :: rc + + ! local variables + integer :: n + type(ESMf_Field) :: lfield + integer :: fieldCount, lrank + character(ESMF_MAXSTR) ,pointer :: lfieldnamelist(:) + real(r8), pointer :: dataPtr1d(:) + real(r8), pointer :: dataPtr2d(:,:) + character(len=*),parameter :: subname='(state_diagnose)' + ! ---------------------------------------------- + + call ESMF_StateGet(state, itemCount=fieldCount, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + allocate(lfieldnamelist(fieldCount)) + + call ESMF_StateGet(state, itemNameList=lfieldnamelist, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + do n = 1, fieldCount + + call ESMF_StateGet(state, itemName=lfieldnamelist(n), field=lfield, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + call field_getfldptr(lfield, fldptr1=dataPtr1d, fldptr2=dataPtr2d, rank=lrank, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + if (lrank == 0) then + ! no local data + elseif (lrank == 1) then + if (size(dataPtr1d) > 0) then + write(msgString,'(A,a)') trim(string)//': for 1d field '//trim(lfieldnamelist(n)) + call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO) + write(msgString,'(A,3g14.7,i8)') trim(string)//': 1d field '//trim(lfieldnamelist(n)), & + minval(dataPtr1d), maxval(dataPtr1d), sum(dataPtr1d), size(dataPtr1d) + call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO) + else + write(msgString,'(A,a)') trim(string)//': '//trim(lfieldnamelist(n))," no data" + call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO) + endif + elseif (lrank == 2) then + if (size(dataPtr2d) > 0) then + write(msgString,'(A,a)') trim(string)//': for 2d field '//trim(lfieldnamelist(n)) + call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO) + write(msgString,'(A,3g14.7,i8)') trim(string)//': 2d field '//trim(lfieldnamelist(n)), & + minval(dataPtr2d), maxval(dataPtr2d), sum(dataPtr2d), size(dataPtr2d) + call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO) + else + write(msgString,'(A,a)') trim(string)//': '//trim(lfieldnamelist(n))," no data" + call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO) + endif + else + call ESMF_LogWrite(trim(subname)//": ERROR rank not supported ", ESMF_LOGMSG_ERROR) + rc = ESMF_FAILURE + return + endif + enddo + + deallocate(lfieldnamelist) + + end subroutine state_diagnose + +!=============================================================================== + + subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc) + + ! ---------------------------------------------- + ! for a field, determine rank and return fldptr1 or fldptr2 + ! abort is true by default and will abort if fldptr is not yet allocated in field + ! rank returns 0, 1, or 2. 0 means fldptr not allocated and abort=false + ! ---------------------------------------------- + + ! input/output variables + type(ESMF_Field) , intent(in) :: field + real(r8), pointer , intent(inout), optional :: fldptr1(:) + real(r8), pointer , intent(inout), optional :: fldptr2(:,:) + integer , intent(out) , optional :: rank + logical , intent(in) , optional :: abort + integer , intent(out) , optional :: rc + + ! local variables + type(ESMF_GeomType_Flag) :: geomtype + type(ESMF_FieldStatus_Flag) :: status + type(ESMF_Mesh) :: lmesh + integer :: lrank, nnodes, nelements + logical :: labort + character(len=*), parameter :: subname='(field_getfldptr)' + ! ---------------------------------------------- + + if (.not.present(rc)) then + call ESMF_LogWrite(trim(subname)//": ERROR rc not present ", & + ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u) + rc = ESMF_FAILURE + return + endif + + rc = ESMF_SUCCESS + + labort = .true. + if (present(abort)) then + labort = abort + endif + lrank = -99 + + call ESMF_FieldGet(field, status=status, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + if (status /= ESMF_FIELDSTATUS_COMPLETE) then + lrank = 0 + if (labort) then + call ESMF_LogWrite(trim(subname)//": ERROR data not allocated ", ESMF_LOGMSG_INFO, rc=rc) + rc = ESMF_FAILURE + return + else + call ESMF_LogWrite(trim(subname)//": WARNING data not allocated ", ESMF_LOGMSG_INFO, rc=rc) + endif + else + + call ESMF_FieldGet(field, geomtype=geomtype, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + if (geomtype == ESMF_GEOMTYPE_GRID) then + call ESMF_FieldGet(field, rank=lrank, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + elseif (geomtype == ESMF_GEOMTYPE_MESH) then + call ESMF_FieldGet(field, rank=lrank, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldGet(field, mesh=lmesh, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_MeshGet(lmesh, numOwnedNodes=nnodes, numOwnedElements=nelements, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + if (nnodes == 0 .and. nelements == 0) lrank = 0 + else + call ESMF_LogWrite(trim(subname)//": ERROR geomtype not supported ", & + ESMF_LOGMSG_INFO, rc=rc) + rc = ESMF_FAILURE + return + endif ! geomtype + + if (lrank == 0) then + call ESMF_LogWrite(trim(subname)//": no local nodes or elements ", & + ESMF_LOGMSG_INFO) + elseif (lrank == 1) then + if (.not.present(fldptr1)) then + call ESMF_LogWrite(trim(subname)//": ERROR missing rank=1 array ", & + ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u) + rc = ESMF_FAILURE + return + endif + call ESMF_FieldGet(field, farrayPtr=fldptr1, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + elseif (lrank == 2) then + if (.not.present(fldptr2)) then + call ESMF_LogWrite(trim(subname)//": ERROR missing rank=2 array ", & + ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u) + rc = ESMF_FAILURE + return + endif + call ESMF_FieldGet(field, farrayPtr=fldptr2, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + else + call ESMF_LogWrite(trim(subname)//": ERROR in rank ", & + ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u) + rc = ESMF_FAILURE + return + endif + + endif ! status + + if (present(rank)) then + rank = lrank + endif + + end subroutine field_getfldptr + +!=============================================================================== + + subroutine alarmInit( clock, alarm, option, & + opt_n, opt_ymd, opt_tod, RefTime, alarmname, advance_clock, rc) + use ESMF, only : ESMF_AlarmPrint + ! Setup an alarm in a clock + ! Notes: The ringtime sent to AlarmCreate MUST be the next alarm + ! time. If you send an arbitrary but proper ringtime from the + ! past and the ring interval, the alarm will always go off on the + ! next clock advance and this will cause serious problems. Even + ! if it makes sense to initialize an alarm with some reference + ! time and the alarm interval, that reference time has to be + ! advance forward to be >= the current time. In the logic below + ! we set an appropriate "NextAlarm" and then we make sure to + ! advance it properly based on the ring interval. + + ! input/output variables + type(ESMF_Clock) , intent(inout) :: clock ! clock + type(ESMF_Alarm) , intent(inout) :: alarm ! alarm + character(len=*) , intent(in) :: option ! alarm option + integer , optional , intent(in) :: opt_n ! alarm freq + integer , optional , intent(in) :: opt_ymd ! alarm ymd + integer , optional , intent(in) :: opt_tod ! alarm tod (sec) + type(ESMF_Time) , optional , intent(in) :: RefTime ! ref time + character(len=*) , optional , intent(in) :: alarmname ! alarm name + logical , optional , intent(in) :: advance_clock ! advance clock to trigger alarm + integer , intent(inout) :: rc ! Return code + + ! local variables + type(ESMF_Calendar) :: cal ! calendar + integer :: lymd ! local ymd + integer :: ltod ! local tod + integer :: cyy,cmm,cdd,csec ! time info + character(len=64) :: lalarmname ! local alarm name + logical :: update_nextalarm ! update next alarm + type(ESMF_Time) :: CurrTime ! Current Time + type(ESMF_Time) :: NextAlarm ! Next restart alarm time + type(ESMF_TimeInterval) :: AlarmInterval ! Alarm interval + type(ESMF_TimeInterval) :: TimeStepInterval ! Component timestep interval + character(len=*), parameter :: subname = '(alarmInit): ' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + lalarmname = 'alarm_unknown' + if (present(alarmname)) lalarmname = trim(alarmname) + ltod = 0 + if (present(opt_tod)) ltod = opt_tod + lymd = -1 + if (present(opt_ymd)) lymd = opt_ymd + + call ESMF_ClockGet(clock, CurrTime=CurrTime, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + call ESMF_TimeGet(CurrTime, yy=cyy, mm=cmm, dd=cdd, s=csec, rc=rc ) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! initial guess of next alarm, this will be updated below + if (present(RefTime)) then + NextAlarm = RefTime + else + NextAlarm = CurrTime + endif + + ! Determine calendar + call ESMF_ClockGet(clock, calendar=cal) + + ! Error checks + if (trim(option) == optdate) then + if (.not. present(opt_ymd)) then + call ESMF_LogWrite(trim(subname)//trim(option)//' requires opt_ymd', ESMF_LOGMSG_ERROR) + rc = ESMF_FAILURE + return + end if + if (lymd < 0 .or. ltod < 0) then + call ESMF_LogWrite(subname//trim(option)//'opt_ymd, opt_tod invalid', ESMF_LOGMSG_ERROR) + rc = ESMF_FAILURE + return + end if + else if (& + trim(option) == optNSteps .or. trim(option) == trim(optNSteps)//'s' .or. & + trim(option) == optNSeconds .or. trim(option) == trim(optNSeconds)//'s' .or. & + trim(option) == optNMinutes .or. trim(option) == trim(optNMinutes)//'s' .or. & + trim(option) == optNHours .or. trim(option) == trim(optNHours)//'s' .or. & + trim(option) == optNDays .or. trim(option) == trim(optNDays)//'s' .or. & + trim(option) == optNMonths .or. trim(option) == trim(optNMonths)//'s' .or. & + trim(option) == optNYears .or. trim(option) == trim(optNYears)//'s' ) then + if (.not.present(opt_n)) then + call ESMF_LogWrite(subname//trim(option)//' requires opt_n', ESMF_LOGMSG_ERROR) + rc = ESMF_FAILURE + return + end if + if (opt_n <= 0) then + call ESMF_LogWrite(subname//trim(option)//' invalid opt_n', ESMF_LOGMSG_ERROR) + rc = ESMF_FAILURE + return + end if + end if + + ! Determine inputs for call to create alarm + selectcase (trim(option)) + + case (optNONE) + call ESMF_TimeIntervalSet(AlarmInterval, yy=9999, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_TimeSet( NextAlarm, yy=9999, mm=12, dd=1, s=0, calendar=cal, rc=rc ) + if (chkerr(rc,__LINE__,u_FILE_u)) return + update_nextalarm = .false. + + case (optNever) + call ESMF_TimeIntervalSet(AlarmInterval, yy=9999, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_TimeSet( NextAlarm, yy=9999, mm=12, dd=1, s=0, calendar=cal, rc=rc ) + if (chkerr(rc,__LINE__,u_FILE_u)) return + update_nextalarm = .false. + + case (optDate) + call ESMF_TimeIntervalSet(AlarmInterval, yy=9999, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call timeInit(NextAlarm, lymd, cal, ltod, rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + update_nextalarm = .false. + + case (optNSteps,trim(optNSteps)//'s') + call ESMF_ClockGet(clock, TimeStep=TimestepInterval, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_TimeIntervalSet(AlarmInterval, s=dtime_drv, rc=rc ) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + AlarmInterval = AlarmInterval * opt_n + ! timestepinterval*0 is 0 of kind ESMF_TimeStepInterval + if (mod(AlarmInterval, TimestepInterval) /= (timestepinterval*0)) then + call ESMF_LogWrite(subname//'illegal Alarm setting for '//trim(alarmname), ESMF_LOGMSG_ERROR) + rc = ESMF_FAILURE + return + endif + update_nextalarm = .true. + + case (optNSeconds,trim(optNSeconds)//'s') + call ESMF_TimeIntervalSet(AlarmInterval, s=1, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + AlarmInterval = AlarmInterval * opt_n + update_nextalarm = .true. + + case (optNMinutes,trim(optNMinutes)//'s') + call ESMF_TimeIntervalSet(AlarmInterval, s=60, rc=rc) + AlarmInterval = AlarmInterval * opt_n + update_nextalarm = .true. + + case (optNHours,trim(optNHours)//'s') + call ESMF_TimeIntervalSet(AlarmInterval, s=3600, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + AlarmInterval = AlarmInterval * opt_n + update_nextalarm = .true. + + case (optNDays,trim(optNDays)//'s') + call ESMF_TimeIntervalSet(AlarmInterval, d=1, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + AlarmInterval = AlarmInterval * opt_n + update_nextalarm = .true. + + case (optNMonths,trim(optNMonths)//'s') + call ESMF_TimeIntervalSet(AlarmInterval, mm=1, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + AlarmInterval = AlarmInterval * opt_n + update_nextalarm = .true. + + case (optMonthly) + call ESMF_TimeIntervalSet(AlarmInterval, mm=1, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_TimeSet( NextAlarm, yy=cyy, mm=cmm, dd=1, s=0, calendar=cal, rc=rc ) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + update_nextalarm = .true. + + case (optNYears, trim(optNYears)//'s') + call ESMF_TimeIntervalSet(AlarmInterval, yy=1, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + AlarmInterval = AlarmInterval * opt_n + update_nextalarm = .true. + + case (optYearly) + call ESMF_TimeIntervalSet(AlarmInterval, yy=1, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_TimeSet( NextAlarm, yy=cyy, mm=1, dd=1, s=0, calendar=cal, rc=rc ) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + update_nextalarm = .true. + case default + call ESMF_LogWrite(subname//'unknown option '//trim(option), ESMF_LOGMSG_ERROR) + rc = ESMF_FAILURE + return + + end select + + ! -------------------------------------------------------------------------------- + ! --- AlarmInterval and NextAlarm should be set --- + ! -------------------------------------------------------------------------------- + + ! --- advance Next Alarm so it won't ring on first timestep for + ! --- most options above. go back one alarminterval just to be careful + + if (update_nextalarm) then + NextAlarm = NextAlarm - AlarmInterval + do while (NextAlarm <= CurrTime) + NextAlarm = NextAlarm + AlarmInterval + enddo + endif + alarm = ESMF_AlarmCreate( name=lalarmname, clock=clock, ringTime=NextAlarm, & + ringInterval=AlarmInterval, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Advance model clock to trigger alarm then reset model clock back to currtime + if (present(advance_clock)) then + if (advance_clock) then + call ESMF_AlarmSet(alarm, clock=clock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_ClockGet(clock, currTime=CurrTime, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_ClockAdvance(clock,rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_ClockSet(clock, currTime=currtime, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + end if + + end subroutine alarmInit + +!=============================================================================== + + subroutine timeInit( Time, ymd, cal, tod, rc) + + ! Create the ESMF_Time object corresponding to the given input time, + ! given in YMD (Year Month Day) and TOD (Time-of-day) format. + ! Set the time by an integer as YYYYMMDD and integer seconds in the day + + ! input/output parameters: + type(ESMF_Time) , intent(inout) :: Time ! ESMF time + integer , intent(in) :: ymd ! year, month, day YYYYMMDD + type(ESMF_Calendar) , intent(in) :: cal ! ESMF calendar + integer , intent(in) :: tod ! time of day in seconds + integer , intent(out) :: rc + + ! local variables + integer :: year, mon, day ! year, month, day as integers + integer :: tdate ! temporary date + character(len=*), parameter :: subname='(timeInit)' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + if ( (ymd < 0) .or. (tod < 0) .or. (tod > SecPerDay) )then + call shr_sys_abort( subname//'ERROR yymmdd is a negative number or time-of-day out of bounds' ) + end if + + tdate = abs(ymd) + year = int(tdate/10000) + if (ymd < 0) year = -year + mon = int( mod(tdate,10000)/ 100) + day = mod(tdate, 100) + + call ESMF_TimeSet( Time, yy=year, mm=mon, dd=day, s=tod, calendar=cal, rc=rc ) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + end subroutine timeInit + +!=============================================================================== + + logical function chkerr(rc, line, file) + + integer, intent(in) :: rc + integer, intent(in) :: line + character(len=*), intent(in) :: file + + integer :: lrc + + chkerr = .false. + lrc = rc + if (ESMF_LogFoundError(rcToCheck=lrc, msg=ESMF_LOGERR_PASSTHRU, line=line, file=file)) then + chkerr = .true. + endif + end function chkerr + +end module nuopc_shr_methods diff --git a/streams/dshr_methods_mod.F90 b/streams/dshr_methods_mod.F90 index b1972f7c8..0ce787596 100644 --- a/streams/dshr_methods_mod.F90 +++ b/streams/dshr_methods_mod.F90 @@ -541,6 +541,7 @@ subroutine dshr_field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc) integer :: ungriddedUBound(1) integer :: lrank logical :: labort + character(len=CS) :: name character(len=*), parameter :: subname='(field_getfldptr)' ! ---------------------------------------------- rc = ESMF_SUCCESS @@ -554,7 +555,9 @@ subroutine dshr_field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (status /= ESMF_FIELDSTATUS_COMPLETE) then if (labort) then - call ESMF_LogWrite(trim(subname)//": ERROR data not allocated ", ESMF_LOGMSG_ERROR, rc=rc) + call ESMF_FieldGet(field, name=name, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_LogWrite(trim(subname)//": field "//trim(name)//" has no data not allocated ", ESMF_LOGMSG_ERROR, rc=rc) rc = ESMF_FAILURE return else @@ -565,7 +568,7 @@ subroutine dshr_field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (ungriddedUBound(1) > 0) then if (.not.present(fldptr2)) then - call ESMF_LogWrite(trim(subname)//": ERROR missing rank=2 array ", & + call ESMF_LogWrite(trim(subname)//": ERROR missing rank=2 array for "//trim(name), & ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u) rc = ESMF_FAILURE return @@ -575,7 +578,7 @@ subroutine dshr_field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc) lrank = 2 else if (.not.present(fldptr1)) then - call ESMF_LogWrite(trim(subname)//": ERROR missing rank=1 array ", & + call ESMF_LogWrite(trim(subname)//": ERROR missing rank=1 array for "//trim(name), & ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u) rc = ESMF_FAILURE return diff --git a/streams/dshr_strdata_mod.F90 b/streams/dshr_strdata_mod.F90 index c056639d3..c9ba174a2 100644 --- a/streams/dshr_strdata_mod.F90 +++ b/streams/dshr_strdata_mod.F90 @@ -21,7 +21,7 @@ module dshr_strdata_mod use ESMF , only : ESMF_REGION_TOTAL, ESMF_FieldGet, ESMF_TraceRegionExit, ESMF_TraceRegionEnter use ESMF , only : ESMF_LOGMSG_INFO, ESMF_LogWrite use shr_kind_mod , only : r8=>shr_kind_r8, r4=>shr_kind_r4, i2=>shr_kind_I2 - use shr_kind_mod , only : cs=>shr_kind_cs, cl=>shr_kind_cl, cxx=>shr_kind_cxx + use shr_kind_mod , only : cs=>shr_kind_cs, cl=>shr_kind_cl, cxx=>shr_kind_cxx, cx=>shr_kind_cx use shr_sys_mod , only : shr_sys_abort use shr_const_mod , only : shr_const_pi, shr_const_cDay, shr_const_spval use shr_cal_mod , only : shr_cal_calendarname, shr_cal_timeSet @@ -51,7 +51,7 @@ module dshr_strdata_mod use pio , only : pio_inquire, pio_inq_varid, pio_inq_varndims, pio_inq_vardimid use pio , only : pio_inq_dimlen, pio_inq_vartype, pio_inq_dimname, pio_inq_dimid use pio , only : pio_double, pio_real, pio_int, pio_offset_kind, pio_get_var - use pio , only : pio_read_darray, pio_setframe, pio_fill_double, pio_get_att + use pio , only : pio_read_darray, pio_setframe, pio_fill_double, pio_get_att, pio_inq_att use pio , only : PIO_BCAST_ERROR, PIO_RETURN_ERROR, PIO_NOERR, PIO_INTERNAL_ERROR, PIO_SHORT implicit none @@ -94,6 +94,7 @@ module dshr_strdata_mod character(CL), allocatable :: fldlist_stream(:) ! names of stream file fields character(CL), allocatable :: fldlist_model(:) ! names of stream model fields integer :: stream_nlev ! number of vertical levels in stream + real(r8), allocatable :: stream_vlevs(:) ! values of vertical levels in stream integer :: stream_lb ! index of the Lowerbound (LB) in fldlist_stream integer :: stream_ub ! index of the Upperbound (UB) in fldlist_stream type(ESMF_Field) :: field_stream ! a field on the stream data domain @@ -389,7 +390,7 @@ subroutine shr_strdata_init(sdat, model_clock, stream_name, rc) character(CS) :: calendar ! calendar name integer :: ns ! stream index integer :: m ! generic index - character(CL) :: fileName ! generic file name + character(CX) :: fileName ! generic file name integer :: nfld ! loop stream field index type(ESMF_Field) :: lfield ! temporary type(ESMF_Field) :: lfield_dst ! temporary @@ -677,9 +678,12 @@ subroutine shr_strdata_get_stream_nlev(sdat, stream_index, rc) type(ESMF_VM) :: vm type(file_desc_t) :: pioid integer :: rcode - character(CL) :: filename + character(CX) :: filename integer :: dimid + type(var_desc_t) :: varid integer :: stream_nlev + integer :: old_handle ! previous setting of pio error handling + character(CS) :: units character(*), parameter :: subname = '(shr_strdata_set_stream_domain) ' ! ---------------------------------------------- @@ -694,14 +698,31 @@ subroutine shr_strdata_get_stream_nlev(sdat, stream_index, rc) if (sdat%mainproc) then call shr_stream_getData(sdat%stream(stream_index), 1, filename) end if - call ESMF_VMBroadCast(vm, filename, CL, 0, rc=rc) + call ESMF_VMBroadCast(vm, filename, CX, 0, rc=rc) rcode = pio_openfile(sdat%pio_subsystem, pioid, sdat%io_type, trim(filename), pio_nowrite) rcode = pio_inq_dimid(pioid, trim(sdat%stream(stream_index)%lev_dimname), dimid) rcode = pio_inq_dimlen(pioid, dimid, stream_nlev) + allocate(sdat%pstrm(stream_index)%stream_vlevs(stream_nlev)) + rcode = pio_inq_varid(pioid, trim(sdat%stream(stream_index)%lev_dimname), varid) + rcode = pio_get_var(pioid, varid, sdat%pstrm(stream_index)%stream_vlevs) + + ! Determine vertical coordinates units - assume that default is m + call pio_seterrorhandling(pioid, PIO_BCAST_ERROR, old_handle) + rcode = pio_inq_att(pioid, varid, 'units') + call pio_seterrorhandling(pioid, old_handle) + if (rcode == PIO_NOERR) then + rcode = pio_get_att(pioid, varid, 'units', units) + if (trim(units) == 'centimeters' .or. trim(units) == 'cm') then + sdat%pstrm(stream_index)%stream_vlevs(:) = sdat%pstrm(stream_index)%stream_vlevs(:) / 100. + end if + end if call pio_closefile(pioid) end if if (sdat%mainproc) then write(sdat%stream(1)%logunit,*) trim(subname)//' stream_nlev = ',stream_nlev + if (stream_nlev /= 1) then + write(sdat%stream(1)%logunit,*)' stream vertical levels = ',sdat%pstrm(stream_index)%stream_vlevs + end if end if ! Set stream_nlev in the per-stream sdat info @@ -726,7 +747,7 @@ subroutine shr_strdata_get_stream_domain(sdat, stream_index, fldname, flddata, r type(var_desc_t) :: varid type(file_desc_t) :: pioid integer :: rcode - character(CL) :: filename + character(CX) :: filename type(io_desc_t) :: pio_iodesc real(r4), allocatable :: data_real(:) real(r8), allocatable :: data_double(:) @@ -743,7 +764,7 @@ subroutine shr_strdata_get_stream_domain(sdat, stream_index, fldname, flddata, r if (sdat%mainproc) then call shr_stream_getData(sdat%stream(stream_index), 1, filename) end if - call ESMF_VMBroadCast(vm, filename, CL, 0, rc=rc) + call ESMF_VMBroadCast(vm, filename, CX, 0, rc=rc) ! Open the file rcode = pio_openfile(sdat%pio_subsystem, pioid, sdat%io_type, trim(filename), pio_nowrite) @@ -1305,10 +1326,10 @@ subroutine shr_strdata_readLBUB(sdat, ns, mDate, mSec, newData, istr, rc) real(r8) :: rDateM,rDateLB,rDateUB ! model,LB,UB dates with fractional days integer :: n_lb, n_ub integer :: i - character(CL) :: filename_lb - character(CL) :: filename_ub - character(CL) :: filename_next - character(CL) :: filename_prev + character(CX) :: filename_lb + character(CX) :: filename_ub + character(CX) :: filename_next + character(CX) :: filename_prev logical :: find_bounds character(*), parameter :: subname = '(shr_strdata_readLBUB) ' character(*), parameter :: F00 = "('(shr_strdata_readLBUB) ',8a)" @@ -1432,7 +1453,7 @@ subroutine shr_strdata_readstrm(sdat, per_stream, stream, fldbun_data, & ! local variables integer :: stream_nlev type(ESMF_Field) :: field_dst, field_vector_dst - character(CL) :: currfile + character(CX) :: currfile logical :: fileexists logical :: fileopen type(file_desc_t) :: pioid @@ -1933,6 +1954,7 @@ subroutine shr_strdata_set_stream_iodesc(sdat, per_stream, fldname, pioid, rc) character(*), parameter :: F00 = "('(shr_strdata_set_stream_iodesc) ',a,i8,2x,i8,2x,a)" character(*), parameter :: F01 = "('(shr_strdata_set_stream_iodesc) ',a,i8,2x,i8,2x,a)" character(*), parameter :: F02 = "('(shr_strdata_set_stream_iodesc) ',a,i8,2x,i8,2x,i8,2x,a)" + character(*), parameter :: F03 = "('(shr_strdata_set_stream_iodesc) ',a,i8,2x,a)" !------------------------------------------------------------------------------- rc = ESMF_SUCCESS @@ -1982,13 +2004,23 @@ subroutine shr_strdata_set_stream_iodesc(sdat, per_stream, fldname, pioid, rc) ! determine io descriptor if (ndims == 2) then - if (sdat%mainproc) then - write(sdat%stream(1)%logunit,F00) 'setting iodesc for : '//trim(fldname)// & - ' with dimlens(1), dimlens2 = ',dimlens(1),dimlens(2),& - ' variable has no time dimension ' + rcode = pio_inq_dimname(pioid, dimids(ndims), dimname) + if (trim(dimname) == 'time' .or. trim(dimname) == 'nt') then + if (sdat%mainproc) then + write(sdat%stream(1)%logunit,F03) 'setting iodesc for : '//trim(fldname)// & + ' with dimlens(1) = ',dimlens(1),' and the variable has a time dimension ' + end if + call pio_initdecomp(sdat%pio_subsystem, pio_iovartype, (/dimlens(1)/), compdof, & + per_stream%stream_pio_iodesc) + else + if (sdat%mainproc) then + write(sdat%stream(1)%logunit,F00) 'setting iodesc for : '//trim(fldname)// & + ' with dimlens(1), dimlens(2) = ',dimlens(1),dimlens(2),& + ' variable has no time dimension ' + end if + call pio_initdecomp(sdat%pio_subsystem, pio_iovartype, (/dimlens(1),dimlens(2)/), compdof, & + per_stream%stream_pio_iodesc) end if - call pio_initdecomp(sdat%pio_subsystem, pio_iovartype, (/dimlens(1),dimlens(2)/), compdof, & - per_stream%stream_pio_iodesc) else if (ndims == 3) then rcode = pio_inq_dimname(pioid, dimids(ndims), dimname) diff --git a/streams/dshr_stream_mod.F90 b/streams/dshr_stream_mod.F90 index 352bf6650..31ec21996 100644 --- a/streams/dshr_stream_mod.F90 +++ b/streams/dshr_stream_mod.F90 @@ -15,7 +15,7 @@ module dshr_stream_mod ! containing those dates. ! ------------------------------------------------------------------------------- - use shr_kind_mod , only : r8=>shr_kind_r8, cs=>shr_kind_cs, cl=>shr_kind_cl, cxx=>shr_kind_cxx + use shr_kind_mod , only : r8=>shr_kind_r8, cs=>shr_kind_cs, cl=>shr_kind_cl, cxx=>shr_kind_cxx, cx=>shr_kind_cx use shr_sys_mod , only : shr_sys_abort use shr_const_mod , only : shr_const_cday use shr_string_mod , only : shr_string_leftalign_and_convert_tabs, shr_string_parseCFtunit @@ -85,7 +85,7 @@ module dshr_stream_mod ! a useful derived type to use inside shr_streamType --- type shr_stream_file_type - character(CL) :: name = shr_stream_file_null ! the file name (full pathname) + character(CX) :: name = shr_stream_file_null ! the file name (full pathname) logical :: haveData = .false. ! has t-coord data been read in? integer :: nt = 0 ! size of time dimension integer ,allocatable :: date(:) ! t-coord date: yyyymmdd @@ -125,7 +125,7 @@ module dshr_stream_mod integer :: n_gvd = -1 ! file/sample of greatest valid date logical :: found_gvd = .false. ! T <=> k_gvd,n_gvd have been set logical :: fileopen = .false. ! is current file open - character(CL) :: currfile = ' ' ! current filename + character(CX) :: currfile = ' ' ! current filename integer :: nvars ! number of stream variables character(CL) :: stream_vectors = 'null' ! stream vectors names type(file_desc_t) :: currpioid ! current pio file desc @@ -379,7 +379,7 @@ subroutine shr_stream_init_from_xml(streamfilename, streamdat, isroot_task, logu allocate(streamdat(i)%varlist(streamdat(i)%nvars)) endif do n=1,streamdat(i)%nfiles - call ESMF_VMBroadCast(vm, streamdat(i)%file(n)%name, CL, 0, rc=rc) + call ESMF_VMBroadCast(vm, streamdat(i)%file(n)%name, CX, 0, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return enddo do n=1,streamdat(i)%nvars @@ -1217,7 +1217,7 @@ subroutine shr_stream_readTCoord(strm, k, isroot_task, rc) integer,optional ,intent(out) :: rc ! return code ! local variables - character(CL) :: fileName ! filename to read + character(CX) :: fileName ! filename to read integer :: nt integer :: num,n integer :: din,dout @@ -1511,7 +1511,7 @@ subroutine shr_stream_getCalendar(strm, k, calendar) type(ESMF_VM) :: vm integer :: myid integer :: vid, n - character(CL) :: fileName + character(CX) :: fileName character(CL) :: lcal integer(PIO_OFFSET_KIND) :: attlen integer :: old_handle @@ -1743,13 +1743,13 @@ subroutine shr_stream_restIO(pioid, streams, mode) integer :: maxnt = 0 integer, allocatable :: tmp(:) integer :: logunit - character(len=CL) :: fname, rfname, rsfname + character(len=CX) :: fname, rfname, rsfname !------------------------------------------------------------------------------- if (mode .eq. 'define') then - rcode = pio_def_dim(pioid, 'strlen', CL, dimid_str) + rcode = pio_def_dim(pioid, 'strlen', CX, dimid_str) do k=1,size(streams) ! maxnfiles is the maximum number of files across all streams logunit = streams(k)%logunit