diff --git a/.github/workflows/extbuild.yml b/.github/workflows/extbuild.yml index 41708115b..8d2283d54 100644 --- a/.github/workflows/extbuild.yml +++ b/.github/workflows/extbuild.yml @@ -20,11 +20,11 @@ jobs: CPPFLAGS: "-I/usr/include -I/usr/local/include" # Versions of all dependencies can be updated here - ESMF_VERSION: v8.5.0 + ESMF_VERSION: v8.6.0 PNETCDF_VERSION: checkpoint.1.12.3 - NETCDF_FORTRAN_VERSION: v4.6.0 + NETCDF_FORTRAN_VERSION: v4.6.1 PIO_VERSION: pio2_6_2 - CDEPS_VERSION: cdeps1.0.21 + CDEPS_VERSION: cdeps1.0.26 steps: - uses: actions/checkout@v3 # Build the ESMF library, if the cache contains a previous build @@ -88,7 +88,7 @@ jobs: ref: ${{ env.CDEPS_VERSION }} - name: Build CDEPS if: steps.cache-cdeps.outputs.cache-hit != 'true' - uses: ESCOMP/CDEPS/.github/actions/buildcdeps@cdeps1.0.21 + uses: ESCOMP/CDEPS/.github/actions/buildcdeps@cdeps1.0.26 with: esmfmkfile: $ESMFMKFILE pio_path: ${GITHUB_WORKSPACE}/pio diff --git a/.github/workflows/srt.yml b/.github/workflows/srt.yml index 34252cb63..1044661ba 100644 --- a/.github/workflows/srt.yml +++ b/.github/workflows/srt.yml @@ -26,8 +26,8 @@ jobs: CPPFLAGS: "-I/usr/include -I/usr/local/include " LDFLAGS: "-L/usr/lib/x86_64-linux-gnu -lnetcdf -lnetcdff -lpnetcdf" # Versions of all dependencies can be updated here - ESMF_VERSION: v8.5.0 - PARALLELIO_VERSION: pio2_6_0 + ESMF_VERSION: v8.6.0 + PARALLELIO_VERSION: pio2_6_2 CIME_MODEL: cesm CIME_DRIVER: nuopc GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} @@ -79,7 +79,21 @@ jobs: - name: checkout externals run: | pushd cesm - ./manage_externals/checkout_externals ccs_config cdeps cime share mct cpl7 parallelio + ./manage_externals/checkout_externals ccs_config cdeps share mct cpl7 parallelio + cd ccs_config + git checkout main + cd ../ + git clone https://github.com/ESMCI/cime + cd cime + if [[ ! -e "${PWD}/.gitmodules.bak" ]] + then + echo "Converting git@github.com to https://github.com urls in ${PWD}/.gitmodules" + + sed -i".bak" "s/git@github.com:/https:\/\/github.com\//g" "${PWD}/.gitmodules" + fi + git submodule update --init + cd ../components/cdeps + git checkout main - name: Cache ESMF id: cache-esmf diff --git a/cesm/driver/ensemble_driver.F90 b/cesm/driver/ensemble_driver.F90 index 2656f10fc..5ca17b1b4 100644 --- a/cesm/driver/ensemble_driver.F90 +++ b/cesm/driver/ensemble_driver.F90 @@ -145,7 +145,7 @@ subroutine SetModelServices(ensemble_driver, rc) integer :: pio_asyncio_stride integer :: pio_asyncio_rootpe integer :: Global_Comm - character(len=CL) :: start_type ! Type of startup + character(len=CL) :: start_type ! Type of startup character(len=7) :: drvrinst character(len=5) :: inst_suffix character(len=CX) :: msgstr @@ -377,10 +377,8 @@ subroutine SetModelServices(ensemble_driver, rc) endif call shr_log_setLogUnit (logunit) ! Create a clock for each driver instance - call esm_time_clockInit(ensemble_driver, driver, logunit, localpet==petList(1), rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - enddo inst = localPet/(ntasks_per_member+pio_asyncio_ntasks) + 1 diff --git a/cesm/driver/esm.F90 b/cesm/driver/esm.F90 index b5207955a..759a4e986 100644 --- a/cesm/driver/esm.F90 +++ b/cesm/driver/esm.F90 @@ -1517,6 +1517,7 @@ subroutine esm_finalize(driver, rc) call ESMF_LogWrite(trim(subname)//": called", ESMF_LOGMSG_INFO) rc = ESMF_SUCCESS + call shr_log_setLogunit(logunit) call ESMF_GridCompGet(driver, vm=vm, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return diff --git a/cesm/driver/esm_time_mod.F90 b/cesm/driver/esm_time_mod.F90 index fc57eaf11..3b56bb953 100644 --- a/cesm/driver/esm_time_mod.F90 +++ b/cesm/driver/esm_time_mod.F90 @@ -282,7 +282,6 @@ subroutine esm_time_clockInit(ensemble_driver, instance_driver, logunit, maintas call ESMF_TimeSet( RefTime, yy=yr, mm=mon, dd=day, s=ref_tod, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_TimeIntervalSet( TimeStep, s=dtime_drv, rc=rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return diff --git a/cesm/flux_atmocn/shr_flux_mod.F90 b/cesm/flux_atmocn/shr_flux_mod.F90 index 741447d93..58f7ae923 100644 --- a/cesm/flux_atmocn/shr_flux_mod.F90 +++ b/cesm/flux_atmocn/shr_flux_mod.F90 @@ -133,7 +133,7 @@ end subroutine shr_flux_adjust_constants ! Thomas Toniazzo (Bjerknes Centre, Bergen) ” !=============================================================================== SUBROUTINE flux_atmOcn(logunit, nMax ,zbot ,ubot ,vbot ,thbot , & - & qbot ,s16O ,sHDO ,s18O ,rbot, & + & qbot, rainc ,s16O ,sHDO ,s18O ,rbot, & & tbot ,us ,vs, pslv, & & ts ,mask , seq_flux_atmocn_minwind, & & sen ,lat ,lwup , & @@ -141,7 +141,10 @@ SUBROUTINE flux_atmOcn(logunit, nMax ,zbot ,ubot ,vbot ,thbot , & & evap ,evap_16O, evap_HDO, evap_18O, & & taux ,tauy ,tref ,qref , & & ocn_surface_flux_scheme, & - & duu10n, ustar_sv ,re_sv ,ssq_sv, & + & add_gusts, & + & duu10n, & + & ugust_out, & + & ustar_sv ,re_sv ,ssq_sv, & & missval) ! !USES: @@ -156,11 +159,13 @@ SUBROUTINE flux_atmOcn(logunit, nMax ,zbot ,ubot ,vbot ,thbot , & integer(IN),intent(in) :: nMax ! data vector length integer(IN),intent(in) :: mask (nMax) ! ocn domain mask 0 <=> out of domain integer(IN),intent(in) :: ocn_surface_flux_scheme + logical ,intent(in) :: add_gusts real(R8) ,intent(in) :: zbot (nMax) ! atm level height (m) real(R8) ,intent(in) :: ubot (nMax) ! atm u wind (m/s) real(R8) ,intent(in) :: vbot (nMax) ! atm v wind (m/s) real(R8) ,intent(in) :: thbot(nMax) ! atm potential T (K) real(R8) ,intent(in) :: qbot (nMax) ! atm specific humidity (kg/kg) + real(R8) ,intent(in) :: rainc(nMax) ! atm precip for convective gustiness (kg/m^3) - RBN 24Nov2008/MDF 31Jan2022 real(R8) ,intent(in) :: s16O (nMax) ! atm H216O tracer conc. (kg/kg) real(R8) ,intent(in) :: sHDO (nMax) ! atm HDO tracer conc. (kg/kg) real(R8) ,intent(in) :: s18O (nMax) ! atm H218O tracer conc. (kg/kg) @@ -188,6 +193,7 @@ SUBROUTINE flux_atmOcn(logunit, nMax ,zbot ,ubot ,vbot ,thbot , & real(R8),intent(out) :: tref (nMax) ! diag: 2m ref height T (K) real(R8),intent(out) :: qref (nMax) ! diag: 2m ref humidity (kg/kg) real(R8),intent(out) :: duu10n(nMax) ! diag: 10m wind speed squared (m/s)^2 + real(R8),intent(out) :: ugust_out(nMax) ! diag: gustiness addition to U10 (m/s) real(R8),intent(out),optional :: ustar_sv(nMax) ! diag: ustar real(R8),intent(out),optional :: re_sv (nMax) ! diag: sqrt of exchange coefficient (water) @@ -257,6 +263,8 @@ SUBROUTINE flux_atmOcn(logunit, nMax ,zbot ,ubot ,vbot ,thbot , & real(R8) :: tdiff(nMax) ! tbot - ts real(R8) :: vscl + real(R8) :: ugust ! function: gustiness as a function of convective rainfall. + real(R8) :: gprec ! convective rainfall argument for ugust qsat(Tk) = 640380.0_R8 / exp(5107.4_R8/Tk) @@ -264,7 +272,7 @@ SUBROUTINE flux_atmOcn(logunit, nMax ,zbot ,ubot ,vbot ,thbot , & cdn(Umps) = 0.0027_R8 / min(33.0000_R8,Umps) + 0.000142_R8 + & 0.0000764_R8 * min(33.0000_R8,Umps) - 3.14807e-13_r8 * min(33.0000_R8,Umps)**6 ! Capped Large and Pond by wind - ! cdn(Umps) = 0.0027_R8 / min(30.0_R8,Umps) + 0.000142_R8 + 0.0000764_R8 * min(30.0_R8,Umps) + ! cdn(Umps) = 0.0027_R8 / min(30.0_R8,Umps) + 0.000142_R8 + 0.0000764_R8 * min(30.0_R8,Umps) ! Capped Large and Pond by Cd ! cdn(Umps) = min(0.0025_R8, (0.0027_R8 / Umps + 0.000142_R8 + 0.0000764_R8 * Umps )) ! Large and Pond @@ -273,6 +281,13 @@ SUBROUTINE flux_atmOcn(logunit, nMax ,zbot ,ubot ,vbot ,thbot , & psimhu(xd) = log((1.0_R8+xd*(2.0_R8+xd))*(1.0_R8+xd*xd)/8.0_R8) - 2.0_R8*atan(xd) + 1.571_R8 psixhu(xd) = 2.0_R8 * log((1.0_R8 + xd*xd)/2.0_R8) + ! Convective gustiness appropriate for input precipitation. + ! Following Regelsperger et al. (2000, J. Clim) + ! Ug = log(1.0+6.69R-0.476R^2) + ! Coefficients X by 8640 for mm/s (from cam) -> cm/day (for above forumla) + ugust(gprec) = log(1._R8+57801.6_r8*gprec-3.55332096e7_r8*(gprec**2)) + + !--- formats ---------------------------------------- character(*),parameter :: subName = '(flux_atmOcn) ' character(*),parameter :: F00 = "('(flux_atmOcn) ',4a)" @@ -327,7 +342,14 @@ SUBROUTINE flux_atmOcn(logunit, nMax ,zbot ,ubot ,vbot ,thbot , & if (mask(n) /= 0) then !--- compute some needed quantities --- - vmag = max(seq_flux_atmocn_minwind, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2) ) + if (add_gusts) then + vmag = max(seq_flux_atmocn_minwind, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2) + ugust(min(rainc(n),6.94444e-4_r8)) ) + ugust_out(n) = ugust(min(rainc(n),6.94444e-4_r8)) + else + vmag = max(seq_flux_atmocn_minwind, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2) ) + ugust_out(n) = 0.0_r8 + end if + if (use_coldair_outbreak_mod) then ! Cold Air Outbreak Modification: ! Increase windspeed for negative tbot-ts @@ -462,6 +484,7 @@ SUBROUTINE flux_atmOcn(logunit, nMax ,zbot ,ubot ,vbot ,thbot , & tref (n) = spval ! 2m reference height temperature (K) qref (n) = spval ! 2m reference height humidity (kg/kg) duu10n(n) = spval ! 10m wind speed squared (m/s)^2 + ugust_out(n) = spval ! gustiness addition (m/s) if (present(ustar_sv)) ustar_sv(n) = spval if (present(re_sv )) re_sv (n) = spval diff --git a/cime_config/buildnml b/cime_config/buildnml index 32be8ead4..ff2553be7 100755 --- a/cime_config/buildnml +++ b/cime_config/buildnml @@ -15,6 +15,7 @@ from CIME.case import Case from CIME.nmlgen import NamelistGenerator from CIME.utils import expect from CIME.utils import get_model, get_time_in_seconds, get_timestamp +from CIME.namelist import literal_to_python_value from CIME.buildnml import create_namelist_infile, parse_input from CIME.XML.files import Files @@ -105,7 +106,8 @@ def _create_drv_namelists(case, infile, confdir, nmlgen, files): config["COMP_OCN"] = case.get_value("COMP_OCN") config["COMP_ROF"] = case.get_value("COMP_ROF") config["COMP_WAV"] = case.get_value("COMP_WAV") - + config["CAMDEV"] = "True" if "CAM%DEV" in case.get_value("COMPSET") else "False" + if ( ( case.get_value("COMP_ROF") == "mosart" @@ -144,6 +146,11 @@ def _create_drv_namelists(case, infile, confdir, nmlgen, files): if config["COMP_OCN"] == "docn" and "aqua" in case.get_value("DOCN_MODE"): nmlgen.set_value("aqua_planet", value=".true.") + # make sure that variable add_gusts is only set to true if compset includes cam_dev + add_gusts = literal_to_python_value(nmlgen.get_value("add_gusts"), type_="logical") + if add_gusts: + expect("CAM%DEV" in case.get_value("COMPSET"),"ERROR: add_gusts can only be set if CAM%DEV in compset {}".format(case.get_value("COMPSET"))) + # -------------------------------- # Overwrite: set component coupling frequencies # -------------------------------- @@ -658,6 +665,7 @@ def buildnml(case, caseroot, component): create_namelist_infile(case, user_nl_file, namelist_infile, infile_text) infile = [namelist_infile] + # create the files nuopc.runconfig, nuopc.runseq, drv_in and drv_flds_in _create_drv_namelists(case, infile, confdir, nmlgen, files) diff --git a/cime_config/config_component.xml b/cime_config/config_component.xml index d73964961..938e0e31c 100644 --- a/cime_config/config_component.xml +++ b/cime_config/config_component.xml @@ -821,15 +821,6 @@ different MPI ranks to different GPUs within the same compute node - - logical - TRUE,FALSE - FALSE - build_def - env_build.xml - TRUE implies that at least one of the components is built threaded (DO NOT EDIT) - - logical TRUE,FALSE @@ -1035,23 +1026,6 @@ this to work. - - char - ESMF_LOGKIND_SINGLE,ESMF_LOGKIND_MULTI,ESMF_LOGKIND_NONE - ESMF_LOGKIND_NONE - run_flags - env_run.xml - - Determines what ESMF log files (if any) are generated when - USE_ESMF_LIB is TRUE. - ESMF_LOGKIND_SINGLE: Use a single log file, combining messages from - all of the PETs. Not supported on some platforms. - ESMF_LOGKIND_MULTI: Use multiple log files -- one per PET. - ESMF_LOGKIND_NONE: Do not issue messages to a log file. - By default, no ESMF log files are generated. - - - char off,low,high,max diff --git a/cime_config/namelist_definition_drv.xml b/cime_config/namelist_definition_drv.xml index 472889f5b..a861586cc 100644 --- a/cime_config/namelist_definition_drv.xml +++ b/cime_config/namelist_definition_drv.xml @@ -18,24 +18,6 @@ - - char - cime_pes - PELAYOUT_attributes - - Determines what ESMF log files (if any) are generated when - USE_ESMF_LIB is TRUE. - ESMF_LOGKIND_SINGLE: Use a single log file, combining messages from - all of the PETs. Not supported on some platforms. - ESMF_LOGKIND_MULTI: Use multiple log files — one per PET. - ESMF_LOGKIND_NONE: Do not issue messages to a log file. - By default, no ESMF log files are generated. - - - $ESMF_LOGFILE_KIND - - - integer pio @@ -978,6 +960,19 @@ + + logical + control + MED_attributes + + add a wind gustiness factor + + + .true. + .false. + + + logical budget diff --git a/cime_config/runseq/runseq_TG.py b/cime_config/runseq/runseq_TG.py index c0bb4ab92..dea8aede5 100644 --- a/cime_config/runseq/runseq_TG.py +++ b/cime_config/runseq/runseq_TG.py @@ -34,7 +34,7 @@ def gen_runseq(case, coupling_times): runseq.add_action ("MED med_phases_post_lnd" , run_lnd) runseq.add_action ("MED med_phases_prep_glc" , med_to_glc) runseq.add_action ("MED -> GLC :remapMethod=redist" , med_to_glc) - runseq.add_action ("GLC" , run_glc) + runseq.add_action ("GLC" , run_glc and med_to_glc) runseq.add_action ("GLC -> MED :remapMethod=redist" , run_glc) runseq.add_action ("MED med_phases_history_write" , True) diff --git a/cime_config/testdefs/testlist_drv.xml b/cime_config/testdefs/testlist_drv.xml index 985bd6ce9..e17b2ffcf 100644 --- a/cime_config/testdefs/testlist_drv.xml +++ b/cime_config/testdefs/testlist_drv.xml @@ -5,36 +5,36 @@ - + - + - + - + - + - + - + - + @@ -46,18 +46,18 @@ - + - + - + - + @@ -69,18 +69,18 @@ - + - + - + - + @@ -92,27 +92,27 @@ - + - + - + - + - + - + @@ -124,27 +124,27 @@ - + - + - + - + - + - + @@ -156,9 +156,9 @@ - + - + @@ -170,24 +170,24 @@ - + - + - + - + - + @@ -200,36 +200,36 @@ - + - + - + - + - + - + - + - + @@ -241,18 +241,18 @@ - + - + - + - + @@ -263,18 +263,18 @@ - + - + - + - + @@ -282,9 +282,9 @@ - + - + diff --git a/doc/source/addendum/req_attributes.rst b/doc/source/addendum/req_attributes.rst index 410303632..ed2130b32 100644 --- a/doc/source/addendum/req_attributes.rst +++ b/doc/source/addendum/req_attributes.rst @@ -6,12 +6,12 @@ The following attributes are obtained from the respective driver and available to all components that the driver uses. In the case of -NEMS, the NEMS driver ingests these attributes from the -``nems.configure`` file. In the case of CESM, the CESM driver ingests +UFS, the UFS driver ingests these attributes from the +``ufs.configure`` file. In the case of CESM, the CESM driver ingests these attributes from the ``nuopc.runconfig`` file. The list of attributes below are separated into application independent attributes and at this time additional attributes required by CESM. There are no -NEMS-specific attributes required by the NEMS application. +UFS-specific attributes required by the UFS application. General @@ -24,7 +24,7 @@ General CMEPS and is also leveraged in some of the custom calculations in the ``prep`` modules. - The currently supported values for ``coupling_mode`` are ``cesm``, ``nems_orig``, ``nems_frac`` and ``hafs``. + The currently supported values for ``coupling_mode`` are ``cesm``, ``ufs.(frac,nfrac).(aoflux)``, and ``hafs``. Scalar attributes ----------------- diff --git a/doc/source/esmflds.rst b/doc/source/esmflds.rst index 960789491..3289ddeb5 100644 --- a/doc/source/esmflds.rst +++ b/doc/source/esmflds.rst @@ -14,8 +14,8 @@ For each supported application, CMEPS contains two specific files that determine Three application specific versions are currently contained within CMEPS: * for CESM: **esmFldsExchange_cesm_mod.F90** and **fd_cesm.yaml** -* for UFS-S2S: **esmFldsExchange_nems_mod.F90** and **fd_nems.yaml** -* for UFS-HAFS: **esmFldsExchange_hafs_mod.F90** and **fd_hafs.yaml** +* for UFS-S2S: **esmFldsExchange_ufs_mod.F90** and **fd_ufs.yaml** +* for UFS-HAFS: **esmFldsExchange_hafs_mod.F90** and **fd_ufs.yaml** CMEPS advertises **all possible fields** that can be imported to and exported by the mediator for the target coupled system. Not all of @@ -23,10 +23,10 @@ these fields will be connected to the various components. The connections will be determined by what the components advertise in their respective advertise phase. -Across applications, component-specific names for the same fields may vary. The field +Across applications, component-specific names for the same fields may vary. The field dictionary is used to define how the application or component-specific name relates -to the name that the CMEPS mediator uses for that field. The mediator variable -names and their application specific aliases are found in the YAML field dictionary. +to the name that the CMEPS mediator uses for that field. The mediator variable +names and their application specific aliases are found in the YAML field dictionary. Details of the naming conventions and API's of this file can be found in the description of the :ref:`exchange of fields in @@ -38,7 +38,7 @@ Field Naming Convention The CMEPS field name convention in the YAML files is independent of the model components. The convention differentiates between variables that are state fields versus flux fields. The naming convention assumes the following one letter designation for the various components as -well as the mediator. +well as the mediator. **import to mediator**:: @@ -58,30 +58,30 @@ well as the mediator. State variables have a 3 character prefix followed by the state name. The prefix has the form ``S[a,i,l,g,o,r,w,x]_`` and is followed by - the field name. - + the field name. + As an example, ``Sx_t`` is the merged surface temperature from land, ice and ocean sent to the atmosphere for CESM. **Flux variables**: - Flux variables specify both source and destination components and have a - 5 character prefix followed by an identifier name of the flux. The first 5 - characters of the flux prefix ``Flmn_`` indicate a flux between - components l and m, computed by component n. The flux-prefix is followed - by the relevant flux-name. - + Flux variables specify both source and destination components and have a + 5 character prefix followed by an identifier name of the flux. The first 5 + characters of the flux prefix ``Flmn_`` indicate a flux between + components l and m, computed by component n. The flux-prefix is followed + by the relevant flux-name. + **mediator import flux prefixes**:: - + Faxa_, atm flux computed by atm Fall_, lnd-atm flux computed by lnd Fioi_, ice-ocn flux computed by ice Faii_, ice_atm flux computed by ice Flrr_, lnd-rof flux computed by rof Firr_, rof-ice flux computed by rof - + **mediator export flux prefixes**:: - + Faxx_, mediator merged fluxes sent to the atm Foxx_, mediator merged fluxes sent to the ocn Fixx_, mediator merged fluxes sent to the ice @@ -122,7 +122,7 @@ The API for this call is: .. code-block:: Fortran call addfld(fldListFr(comp_index)%flds, 'field_name') - call addfld(fldListTo(comp_index)%flds, 'field_name') + call addfld(fldListTo(comp_index)%flds, 'field_name') where: @@ -206,7 +206,7 @@ fraction corrections are not required in other mappings to improve accuracy beca call addmap(fldListFr(compice)%flds, 'Si_snowh', compatm, mapconsf, 'ifrac', 'unset') This will create an entry in ``fldListFr(compatm)`` specifying that the ``Si_snowh`` field from the ice should be mapped conservatively to the atmosphere using -fractional normalization where the ice fraction is obtained from ``FBFrac(compice)[snowh]``. The route handle for this mapping will be created at run time. +fractional normalization where the ice fraction is obtained from ``FBFrac(compice)[snowh]``. The route handle for this mapping will be created at run time. .. _addmrg: diff --git a/doc/source/generic.rst b/doc/source/generic.rst index 62055af1c..1be409a9a 100644 --- a/doc/source/generic.rst +++ b/doc/source/generic.rst @@ -15,7 +15,7 @@ application specific and provide general functionality. component state in the mediator's InternalState * initializing the mediator component specific fields via a call to - ``esmFldsExchange_xxx_`` (where currently xxx can be ``cesm``, ``nems`` or ``hafs``). + ``esmFldsExchange_xxx_`` (where currently xxx can be ``cesm``, ``ufs`` or ``hafs``). * determining which components are present diff --git a/doc/source/introduction.rst b/doc/source/introduction.rst index 3b79e1ed0..54a16761b 100644 --- a/doc/source/introduction.rst +++ b/doc/source/introduction.rst @@ -31,7 +31,7 @@ matching of standard field names. These standard names are defined in a field dictionary. Since CMEPS is a community mediator, these standard names are specific to each application. - + Organization of the CMEPS mediator code ####################################### @@ -39,28 +39,28 @@ Organization of the CMEPS mediator code When you check out the code you will files, which can be organized into three groups: -* totally generic components that carry out the mediator functionality such as mapping, - merging, restarts and history writes. Included here is a a "fraction" module that - determines the fractions of different source model components on every source +* totally generic components that carry out the mediator functionality such as mapping, + merging, restarts and history writes. Included here is a a "fraction" module that + determines the fractions of different source model components on every source destination mesh. -* application specific code that determines what fields are exchanged between +* application specific code that determines what fields are exchanged between components and how they are merged and mapped. -* prep phase modules that carry out the mapping and merging from one or more +* prep phase modules that carry out the mapping and merging from one or more source components to the destination component. =========================== ============================ =========================== Generic Code Application Specific Code Prep Phase Code =========================== ============================ =========================== med.F90 esmFldsExchange_cesm_mod.F90 med_phases_prep_atm_mod.F90 -esmFlds.F90 esmFldsExchange_nems_mod.F90 med_phases_prep_ice_mod.F90 +esmFlds.F90 esmFldsExchange_ufs_mod.F90 med_phases_prep_ice_mod.F90 med_map_mod.F90 esmFldsExchange_hafs_mod.F90 med_phases_prep_ocn_mod.F90 med_merge_mod.F90 fd_cesm.yaml med_phases_prep_glc_mod.F90 -med_frac_mod.F90 fd_nems.yaml med_phases_prep_lnd_mod.F90 -med_internalstate_mod.F90 fd_hafs.yaml med_phases_prep_rof_mod.F90 -med_methods_mod.F90. -med_phases_aofluxes_mod.F90 +med_frac_mod.F90 fd_ufs.yaml med_phases_prep_lnd_mod.F90 +med_internalstate_mod.F90 fd_ufs.yaml med_phases_prep_rof_mod.F90 +med_methods_mod.F90. +med_phases_aofluxes_mod.F90 med_phases_ocnalb_mod.F90 med_phases_history_mod.F90 med_phases_restart_mod.F90 @@ -78,8 +78,8 @@ Mapping and Merging Primer ####################################### This section provides a primer on mapping (interpolation) and merging of gridded -coupled fields. Masks, support for partial fractions on grids, weights generation, -and fraction +coupled fields. Masks, support for partial fractions on grids, weights generation, +and fraction weighted mapping and merging all play roles in the conservation and quality of the coupled fields. @@ -89,11 +89,11 @@ A pair of atmosphere and ocean/ice grids can be used to highlight the analysis. :width: 400 :alt: Sample CMEPS grids -The most general CMEPS mediator assumes the ocean and sea ice surface grids are +The most general CMEPS mediator assumes the ocean and sea ice surface grids are identical while the atmosphere and land grids are also identical. The ocean/ice grid defines the mask which means each ocean/ice gridcell is either a fully -active ocean/ice gridcell or not (i.e. land). Other configurations have been -and can be implemented and analyzed as well. +active ocean/ice gridcell or not (i.e. land). Other configurations have been +and can be implemented and analyzed as well. The ocean/ice mask interpolated to the atmosphere/land grid determines the complementary ocean/ice and land masks on the atmosphere grid. @@ -112,12 +112,12 @@ The gridcells can be labeled as follows. :width: 300 :alt: Sample CMEPS gridcell naming convention -The atmosphere gridcell is labeled "a". On the atmosphere gridcell (the red box), +The atmosphere gridcell is labeled "a". On the atmosphere gridcell (the red box), in general, there is a land fraction (fal), an ocean fraction (fao), and a sea ice fraction (fai). The sum of the surface fractions should always be 1.0 in these -conventions. There is also a gridbox average field on the atmosphere grid (Fa). -This could be a flux or a state that is +conventions. There is also a gridbox average field on the atmosphere grid (Fa). +This could be a flux or a state that is derived from the equivalent land (Fal), ocean (Fao), and sea ice (Fai) fields. The gridbox average field is computed by merging the various surfaces:: @@ -130,9 +130,9 @@ This is a standard merge where:: and each surface field, Fal, Fao, and Fai are the values of the surface fields on the atmosphere grid. -The ocean gridcells (blue boxes) are labeled 1, 2, 3, and 4 in this example. -In general, -each ocean/ice gridcell partially overlaps multiple atmosphere gridcells. +The ocean gridcells (blue boxes) are labeled 1, 2, 3, and 4 in this example. +In general, +each ocean/ice gridcell partially overlaps multiple atmosphere gridcells. Each ocean/ice gridcell has an overlapping Area (A) and a Mask (M) associated with it. In this example, land is colored green, ocean blue, and sea ice white so just for the figure depicted:: @@ -159,7 +159,7 @@ gridcells. Nonlinear interpolation is not yet supported in most coupled systems Mapping weights can be defined in a number of ways even beyond conservative or bilinear. They can be masked or normalized using multiple approaches. The -weights generation is intricately tied to other aspects of the coupling method. +weights generation is intricately tied to other aspects of the coupling method. In CMEPS, area-overlap conservative weights are defined as follows:: w1 = A1/Aa @@ -167,7 +167,7 @@ In CMEPS, area-overlap conservative weights are defined as follows:: w3 = A3/Aa w4 = A4/Aa -This simple approach which does not include any masking or normalization provides a +This simple approach which does not include any masking or normalization provides a number of useful attributes. The weights always add up to 1.0:: w1 + w2 + w3 + w4 = 1.0 @@ -207,11 +207,11 @@ And the equation for f_land and fal above are consistent if fl1=1-M1:: fal = w1 + w2 + w3 + w4 - (w1*M1 + w2*M2 + w3*M3 + w4*M4) fal = 1 - (w1*M1 + w2*M2 + w3*M3 + w4*M4) -Clearly defined and consistent weights, areas, fractions, and masks is critical +Clearly defined and consistent weights, areas, fractions, and masks is critical to generating conservation in the system. When mapping masked or fraction weighted fields, these weights require that the -mapped field be normalized by the mapped fraction. Consider a case where sea +mapped field be normalized by the mapped fraction. Consider a case where sea surface temperature (SST) is to be mapped to the atmosphere grid with:: M1 = 0; M2 = M3 = M4 = 1 @@ -223,15 +223,15 @@ because w1 is non-zero and Fo1 is underfined since it's a land gridcell on the ocean grid. A masked weighted average, **Fa = M1*w1*Fo1 + M2*w2*Fo2 + M3*w3*Fo3 + M4*w4*Fo4 is also problematic** because M1 is zero, so the contribution of the first term is zero. But the sum -of the remaining weights (M2*w2 + M3*w3 + M4*w4) is now not identically 1 -which means the weighted average is incorrect. (To test this, assume all the +of the remaining weights (M2*w2 + M3*w3 + M4*w4) is now not identically 1 +which means the weighted average is incorrect. (To test this, assume all the weights are each 0.25 and all the Fo values are 10 degC, Fa would then be 7.5 degC). Next consider a masked weighted normalized average, **f_ocean = (w1*M1 + w2*M2 + w3*M3 + w4*M4) combined with Fa = (M1*w1*Fo1 + M2*w2*Fo2 + M3*w3*Fo3 + M4*w4*Fo4) / (f_ocean) which produces a reasonable but incorrect result** because the weighted average uses the mask instead of the fraction. The mask only produces a correct result -in cases where there is no sea ice because sea ice impacts the surface fractions. +in cases where there is no sea ice because sea ice impacts the surface fractions. Finally, consider a fraction weighted normalized average using the dynamically varying ocean fraction that is exposed to the atmosphere:: @@ -249,9 +249,9 @@ fao is the mapped ocean fraction on the atmosphere gridcell, and Fa is the mapped SST. The ocean fractions are only defined where the ocean mask is 1, otherwise the ocean and sea ice fractions are zero. Now, the SST in each ocean gridcell is weighted by the fraction of the ocean -box exposed to the atmosphere and that weighted average is normalized by +box exposed to the atmosphere and that weighted average is normalized by the mapped dynamically varying fraction. This produces a reasonable result -as well as a conservative result. +as well as a conservative result. The conservation check involves thinking of Fo and Fa as a flux. On the ocean grid, the quantity associated with the flux is:: @@ -268,7 +268,7 @@ Via some simple math, it can be shown that Qo = Qa if:: fao = w1*fo1 + w2*fo2 + w3*fo3 + w4*fo4 Fao = (fo1*w1*Fo1 + fo2*w2*Fo2 + fo3*w3*Fo3 + fo4*w4*Fo4) / (fao) -In practice, the fraction weighted normlized mapping field is computed +In practice, the fraction weighted normlized mapping field is computed by mapping the ocean fraction and the fraction weighted field from the ocean to the atmosphere grid separately and then using the mapped fraction to normalize the field as a four step process:: @@ -278,19 +278,19 @@ using the mapped fraction to normalize the field as a four step process:: Fao' = w1*Fo1' + w2*Fo2' + w3*Fo3' + w4*Fo4' (c) Fao = Fao'/fao (d) -Steps (b) and (c) above are the sparse matrix multiply by the standard +Steps (b) and (c) above are the sparse matrix multiply by the standard conservative weights. -Step (a) fraction weighs the field and step (d) normalizes the mapped field. +Step (a) fraction weighs the field and step (d) normalizes the mapped field. Another way to think of this is that the mapped flux (Fao') is normalized by the -same fraction (fao) that is used in the merge, so they actually cancel. -Both the normalization at the end of the mapping and the fraction weighting +same fraction (fao) that is used in the merge, so they actually cancel. +Both the normalization at the end of the mapping and the fraction weighting in the merge can be skipped and the results should be identical. But then the mediator will carry around Fao' instead of Fao and that field is far less intuitive as it no longer represents the gridcell average value, but some subarea average value. In addition, that approach is only valid when carrying out full surface merges. If, -for instance, the SST is to be interpolated and not merged with anything, the field +for instance, the SST is to be interpolated and not merged with anything, the field must be normalized after mapping to be useful. The same mapping and merging process is valid for the sea ice:: @@ -311,23 +311,23 @@ where now:: Fao = (fo1*w1*Fo1 + fo2*w2*Fo2 + fo3*w3*Fo3 + fo4*w4*Fo4) / (fao) Fai = (fi1*w1*Fi1 + fi2*w2*Fi2 + fi3*w3*Fi3 + fi4*w4*Fi4) / (fai) -will simplify to an equation that contains twelve distinct terms for each of the +will simplify to an equation that contains twelve distinct terms for each of the four ocean gridboxes and the three different surfaces:: - Fa = (w1*fl1*Fl1 + w2*fl2*Fl2 + w3*fl3*Fl3 + w4*fl4*Fl4) + - (w1*fo1*Fo1 + w2*fo2*Fo2 + w3*fo3*Fo3 + w4*fo4*Fo4) + - (w1*fi1*Fi1 + w2*fi2*Fi2 + w3*fi3*Fi3 + w4*fi4*Fi4) + Fa = (w1*fl1*Fl1 + w2*fl2*Fl2 + w3*fl3*Fl3 + w4*fl4*Fl4) + + (w1*fo1*Fo1 + w2*fo2*Fo2 + w3*fo3*Fo3 + w4*fo4*Fo4) + + (w1*fi1*Fi1 + w2*fi2*Fi2 + w3*fi3*Fi3 + w4*fi4*Fi4) and this further simplifies to something that looks like a mapping of the field merged on the ocean grid:: - Fa = w1*(fl1*Fl1+fo1*Fo1+fi1*Fi1) + + Fa = w1*(fl1*Fl1+fo1*Fo1+fi1*Fi1) + w2*(fl2*Fl2+fo2*Fo2+fi2*Fi2) + - w3*(fl3*Fl3+fo3*Fo3+fi3*Fi3) + + w3*(fl3*Fl3+fo3*Fo3+fi3*Fi3) + w4*(fl4*Fl4+fo4*Fo4+fi4*Fi4) Like the exercise with Fao above, these equations can be shown to be -fully conservative. +fully conservative. To summarize, multiple features such as area calculations, weights, masking, normalization, fraction weighting, and merging approaches @@ -335,9 +335,9 @@ have to be considered together to ensure conservation. The CMEPS mediator uses unmasked and unnormalized weights and then generally maps using the fraction weighted normalized approach. Merges are carried out with fraction weights. -This is applied to both state and flux fields, with conservative, bilinear, +This is applied to both state and flux fields, with conservative, bilinear, and other mapping approaches, and for both merged and unmerged fields. -This ensures that the fields are always useful gridcell average values +This ensures that the fields are always useful gridcell average values when being coupled or analyzed throughout the coupling implementation. @@ -353,7 +353,7 @@ model discretization, they are NOT ad-hoc. If the previous section, areas and weights were introduced. Those areas were assumed to consist of the area overlaps between gridcells and were computed -using a consistent approach such that the areas conserve. ESMF is able to compute +using a consistent approach such that the areas conserve. ESMF is able to compute these area overlaps and the corresponding mapping weights such that fluxes can be mapped and quantities are conserved. @@ -369,13 +369,13 @@ ESMF are identical, and all the weights are 1.0. So:: F2*A2 = F1*A1 (conservation) Now lets assume that the two models have fundamentally different discretizations, -different area algorithms (i.e. great circle vs simpler lon/lat approximations), +different area algorithms (i.e. great circle vs simpler lon/lat approximations), or even different assumptions about the size and shape of the earth. The grids can be identical in -terms of the longitude and latitude of the +terms of the longitude and latitude of the gridcell corners and centers, but the areas can also -be different because of the underlying model implementation. When a flux is passed -to or from each component, the quantity associated with that flux is proportional to +be different because of the underlying model implementation. When a flux is passed +to or from each component, the quantity associated with that flux is proportional to the model area, so:: A1 = A2 (ESMF areas) @@ -385,7 +385,7 @@ the model area, so:: A1m != A2m (model areas) F1*A1m != F2*A2m (loss of conservation) -This can be corrected by multiplying the fluxes +This can be corrected by multiplying the fluxes by an area correction. For each model, outgoing fluxes should be multiplied by the model area divided by the ESMF area. Incoming fluxes should be multiplied by the ESMF area divided by the model area. So:: @@ -411,14 +411,14 @@ can actually be applied a number of ways. * Models can pass the areas to the mediator and the mediator can multiple fluxes by the source model area before mapping and divide by the destination model area area after mapping. * Models can pass the areas to the mediator and implement an area correction term on the incoming and outgoing fluxes that is the ratio of the model and ESMF areas. This is the approach shown above and is how CMEPS traditionally implements this feature. -Model areas should be passed to the mediator at initialization so the area corrections +Model areas should be passed to the mediator at initialization so the area corrections can be computed and applied. These area corrections do not vary in time. Lags, Accumulation and Averaging ####################################### -In a coupled model, the component model sequencing and coupling frequency tend to introduce +In a coupled model, the component model sequencing and coupling frequency tend to introduce some lags as well as a requirement to accumulate and average. This occurs when component models are running sequentially or concurrently. In general, the component models advance in time separately and the "current time" in each model becomes out of @@ -430,7 +430,7 @@ multiple timesteps are taken between coupling periods in a component model, the states should be averaged over those timesteps before being passed back out to the coupler. In the same way, the fluxes and states passed into the coupler should be averaged over shorter coupling periods for models that are coupled at longer coupling -periods. +periods. For conservation of mass and energy, the field that is accumluated should be consistent with the field that would be passed if there were no averaging required. Take for @@ -447,13 +447,13 @@ where sum_n represents the sum over n time periods. This can also be written as Fo = 1/n * (sum_n(fao*Fao) + sum_n(fio*Fio)) -So multiple terms can be summed and accumulated or the individual terms fao*Fao +So multiple terms can be summed and accumulated or the individual terms fao*Fao and fio*Fio can be accumulated and later summed and averaged in either order. Both approaches produce identical results. Finally, **it's important to note that sum_n(fao)*sum_n(Fao) does not produce the same results as the sum_n(fao*Fao)**. In other words, the fraction weighted flux has to be accumulated and NOT the fraction and flux separately. This is important for conservation -in flux coupling. The same approach should be taken with merged states to compute the +in flux coupling. The same approach should be taken with merged states to compute the most accurate representation of the average state over the slow coupling period. An analysis and review of each coupling field should be carried out to determine the most conservative and accurate representation of averaged fields. This is particularly @@ -493,14 +493,14 @@ Simplifying the above equation:: Fo = 1/n * sum_n(mapa2o(fao_a*Fao_a) -Accumulation (sum_n) and mapping (mapa2o) are both linear operations so this can +Accumulation (sum_n) and mapping (mapa2o) are both linear operations so this can be written as:: Fo = 1/n * mapa2o(sum_n(fao_a*Fao_a)) Fo = mapa2o(1/n*sum_n(fao_a*Fao_a)) which suggests that the accumulation can be done on the source side (i.e. atmosphere) -and only mapped on the slow coupling period. But again, fao_a*Fao_a has to be +and only mapped on the slow coupling period. But again, fao_a*Fao_a has to be accumulated and then when mapped, NO fraction would be applied to the merge as this is already included in the mapped field. In equation form, the full merged ocean field would be implemented as:: @@ -520,7 +520,7 @@ two atmosphere fields are mapped every fast coupling period, the merge is now fraction weighted for all terms, and the mapped fields, fao_o and Fao_o, have physically meaningful values. Fao'_o above does not. This implementation has a parallel with the normalization step. As suggested above, there are two -implementations for conservative mapping and merging in general. The one outlined +implementations for conservative mapping and merging in general. The one outlined above with fraction weighted normalized mapping and fraction weighted merging:: @@ -536,7 +536,7 @@ fraction is NOT applied during the merge:: These will produce identical results in the same way that their accumulated averages do. - + Flux Calculation Grid @@ -564,7 +564,7 @@ equations:: Fa = fl_a*Fal_a + fo_a*Fao_a + fi_a*Fai_a Fo = fo_o*Fao_o + fi_o*Fio_o -The above equations indicate that the land fraction on the atmosphere grid is the +The above equations indicate that the land fraction on the atmosphere grid is the complement of the mapped ocean mask and is static. The ice and ocean fractions are determined from the ice model and are dynamic. Both can be mapped to the atmosphere grid. Finally, the atmosphere flux is a three-way merge of the land, ocean, and @@ -572,7 +572,7 @@ ice terms on the atmosphere grid while the ocean flux is a two-way merge of the atmosphere and ice terms on the ocean grid. When the atmosphere/ocean and atmosphere/ice fluxes are both computed on the same -grid, at the same frequency, and both are mapped to the atmosphere grid, conservative +grid, at the same frequency, and both are mapped to the atmosphere grid, conservative mapping and merging is relatively straight-forward:: fo_a = mapo2a(fo_o) @@ -588,15 +588,15 @@ and everything conserves relatively directly:: fi_a*Fai_a = fi_o*Fai_o When the atmosphere/ice fluxes are computed on the ocean grid while -the atmosphere/ocean fluxes are computed on the atmosphere grid, +the atmosphere/ocean fluxes are computed on the atmosphere grid, extra care is needed with regard to fractions and conservation. In this case:: fo_a = mapo2a(fo_o) Fao_o = mapa2o(fo_a*Fao_a)/mapa2o(fo_a) fi_a = mapo2a(fi_o) Fai_a = mapo2a(fi_o*Fai_o)/fi_a - -fo_o, fi_o, Fai_o, and Fao_a are specified and Fao_o has to be computed. The most + +fo_o, fi_o, Fai_o, and Fao_a are specified and Fao_o has to be computed. The most important point here is that during the ocean merge, the mapped ocean fraction on the atmosphere grid is used so:: diff --git a/doc/source/prep.rst b/doc/source/prep.rst index 07595cb45..e75bf33a7 100644 --- a/doc/source/prep.rst +++ b/doc/source/prep.rst @@ -6,20 +6,20 @@ The following modules comprise the "prep phase" CMEPS code: -**med_phases_prep_atm_mod.F90**: prepares the mediator export state to the atmosphere component +**med_phases_prep_atm_mod.F90**: prepares the mediator export state to the atmosphere component + +**med_phases_prep_ice_mod.F90**: prepares the mediator export state to the sea-ice component + +**med_phases_prep_glc_mod.F90**: prepares the mediator export state to the land-ice component -**med_phases_prep_ice_mod.F90**: prepares the mediator export state to the sea-ice component - -**med_phases_prep_glc_mod.F90**: prepares the mediator export state to the land-ice component - **med_phases_prep_lnd_mod.F90**: prepares the mediator export state to the land component - + **med_phases_prep_ocn_mod.F90**: prepares the mediator export state to the ocean component **med_phases_prep_rof_mod.F90**: prepares the mediator export state to the river component - + **med_phases_prep_wav_mod.F90**: prepares the mediator export state to the wave component - + Each prep phase module has several sections: @@ -71,8 +71,7 @@ Each prep phase module has several sections: * ``med_phases_prep_ocn``: * computation of net shortwave that is sent to the ocean. - * apply precipitation fractor to scale rain and snow sent to ocean (for CESM) - * carry out custom merges for NEMS coupling modes (for NEMS) + * apply precipitation fractor to scale rain and snow sent to ocean * ``med_phases_prep_rof``: diff --git a/mediator/CMakeLists.txt b/mediator/CMakeLists.txt index 1ba39999d..031a23308 100644 --- a/mediator/CMakeLists.txt +++ b/mediator/CMakeLists.txt @@ -5,7 +5,7 @@ set(SRCFILES esmFldsExchange_cesm_mod.F90 med_fraction_mod.F90 med_phases_restart_mod.F90 esmFldsExchange_hafs_mod.F90 med_internalstate_mod.F90 med_phases_aofluxes_mod.F90 med_phases_prep_lnd_mod.F90 med_time_mod.F90 - esmFldsExchange_nems_mod.F90 med_io_mod.F90 + esmFldsExchange_ufs_mod.F90 med_io_mod.F90 med_phases_history_mod.F90 med_phases_prep_ocn_mod.F90 med_utils_mod.F90 esmFlds.F90 med_kind_mod.F90 med_phases_prep_rof_mod.F90 diff --git a/mediator/ESMFConvenienceMacros.h b/mediator/ESMFConvenienceMacros.h index 092760585..660583eaa 100644 --- a/mediator/ESMFConvenienceMacros.h +++ b/mediator/ESMFConvenienceMacros.h @@ -2,6 +2,6 @@ // ----------- ERROR handling macros ------------------------------------------ #endif -#define ESMF_ERR_ABORT(rc) if (ESMF_LogFoundError(rc, msg="Aborting NEMS", line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT) +#define ESMF_ERR_ABORT(rc) if (ESMF_LogFoundError(rc, msg="Aborting UFS", line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT) #define ESMF_ERR_RETURN(rc,rcOut) if (ESMF_LogFoundError(rc, msg="Breaking out of subroutine", line=__LINE__, file=__FILE__, rcToReturn=rcOut)) return diff --git a/mediator/ESMFVersionDefine.h b/mediator/ESMFVersionDefine.h index 0038c9db1..b9a5054c8 100644 --- a/mediator/ESMFVersionDefine.h +++ b/mediator/ESMFVersionDefine.h @@ -1,9 +1,8 @@ #if 0 // // Make this header file available as ESMFVersionDefine.h in order to build -// NEMS against an ESMF installation that contains a reference level NUOPC Layer. +// UFS against an ESMF installation that contains a reference level NUOPC Layer. // #endif #include "./ESMFConvenienceMacros.h" - diff --git a/mediator/Makefile b/mediator/Makefile index 126d040bd..990fe58eb 100644 --- a/mediator/Makefile +++ b/mediator/Makefile @@ -34,12 +34,12 @@ med_kind_mod.o : med_constants_mod.o : med_kind_mod.o esmFlds.o : med_kind_mod.o esmFldsExchange_cesm_mod.o : med_kind_mod.o med_methods_mod.o esmFlds.o med_internalstate_mod.o med_utils_mod.o -esmFldsExchange_nems_mod.o : med_kind_mod.o med_methods_mod.o esmFlds.o med_internalstate_mod.o med_utils_mod.o +esmFldsExchange_ufs_mod.o : med_kind_mod.o med_methods_mod.o esmFlds.o med_internalstate_mod.o med_utils_mod.o esmFldsExchange_hafs_mod.o : med_kind_mod.o med_methods_mod.o esmFlds.o med_internalstate_mod.o med_utils_mod.o med.o : med_kind_mod.o med_phases_profile_mod.o med_utils_mod.o med_phases_prep_rof_mod.o med_phases_aofluxes_mod.o \ med_phases_prep_ice_mod.o med_fraction_mod.o med_map_mod.o med_constants_mod.o med_phases_prep_wav_mod.o \ med_phases_prep_lnd_mod.o med_phases_history_mod.o med_phases_ocnalb_mod.o med_phases_restart_mod.o \ - med_time_mod.o med_internalstate_mod.o med_phases_prep_atm_mod.o esmFldsExchange_cesm_mod.o esmFldsExchange_nems_mod.o \ + med_time_mod.o med_internalstate_mod.o med_phases_prep_atm_mod.o esmFldsExchange_cesm_mod.o esmFldsExchange_ufs_mod.o \ esmFldsExchange_hafs_mod.o med_phases_prep_glc_mod.o esmFlds.o med_io_mod.o med_methods_mod.o med_phases_prep_ocn_mod.o \ med_phases_post_atm_mod.o med_phases_post_ice_mod.o med_phases_post_lnd_mod.o med_phases_post_glc_mod.o med_phases_post_rof_mod.o \ med_phases_post_wav_mod.o diff --git a/mediator/esmFldsExchange_cesm_mod.F90 b/mediator/esmFldsExchange_cesm_mod.F90 index 6e049a0a2..b18833fad 100644 --- a/mediator/esmFldsExchange_cesm_mod.F90 +++ b/mediator/esmFldsExchange_cesm_mod.F90 @@ -276,6 +276,7 @@ subroutine esmFldsExchange_cesm(gcomp, phase, rc) call addfld_from(compatm, 'Sa_shum') call addfld_from(compatm, 'Sa_ptem') call addfld_from(compatm, 'Sa_dens') + call addfld_from(compatm, 'Faxa_rainc') if (flds_wiso) then call addfld_from(compatm, 'Sa_shum_wiso') end if @@ -288,6 +289,7 @@ subroutine esmFldsExchange_cesm(gcomp, phase, rc) call addmap_from(compatm, 'Sa_u' , compocn, mappatch, 'one', atm2ocn_map) call addmap_from(compatm, 'Sa_v' , compocn, mappatch, 'one', atm2ocn_map) end if + call addmap_from(compatm, 'Faxa_rainc', compocn, mapconsf, 'one', atm2ocn_map) call addmap_from(compatm, 'Sa_z' , compocn, mapbilnr, 'one', atm2ocn_map) call addmap_from(compatm, 'Sa_tbot', compocn, mapbilnr, 'one', atm2ocn_map) call addmap_from(compatm, 'Sa_pbot', compocn, mapbilnr, 'one', atm2ocn_map) @@ -1365,6 +1367,24 @@ subroutine esmFldsExchange_cesm(gcomp, phase, rc) end if end if + ! --------------------------------------------------------------------- + ! to atm: unmerged ugust_out from ocn + ! --------------------------------------------------------------------- + if (phase == 'advertise') then + call addfld_aoflux('So_ugustOut') + call addfld_to(compatm, 'So_ugustOut') + else + if ( fldchk(is_local%wrap%FBexp(compatm), 'So_ugustOut', rc=rc)) then + if (fldchk(is_local%wrap%FBMed_aoflux_o, 'So_ugustOut', rc=rc)) then + if (trim(is_local%wrap%aoflux_grid) == 'ogrid') then + call addmap_aoflux('So_ugustOut', compatm, mapconsf, 'ofrac', ocn2atm_map) + end if + call addmrg_to(compatm , 'So_ugustOut', & + mrg_from=compmed, mrg_fld='So_ugustOut', mrg_type='merge', mrg_fracname='ofrac') + end if + end if + end if + ! --------------------------------------------------------------------- ! to atm: surface snow depth from ice (needed for cam) ! to atm: mean ice volume per unit area from ice diff --git a/mediator/esmFldsExchange_hafs_mod.F90 b/mediator/esmFldsExchange_hafs_mod.F90 index 1f645524e..b545b9b1c 100644 --- a/mediator/esmFldsExchange_hafs_mod.F90 +++ b/mediator/esmFldsExchange_hafs_mod.F90 @@ -13,6 +13,7 @@ module esmFldsExchange_hafs_mod use med_internalstate_mod , only : compwav use med_internalstate_mod , only : ncomps use med_internalstate_mod , only : coupling_mode + use esmFlds , only : addfld_ocnalb => med_fldList_addfld_ocnalb !--------------------------------------------------------------------- ! This is a mediator specific routine that determines ALL possible @@ -133,7 +134,7 @@ subroutine esmFldsExchange_hafs_advt(gcomp, phase, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return !===================================================================== - ! FIELDS TO MEDIATOR component (for fractions and atm/ocn flux calculation) + ! Mediator fields !===================================================================== !---------------------------------------------------------- @@ -146,6 +147,18 @@ subroutine esmFldsExchange_hafs_advt(gcomp, phase, rc) !---------------------------------------------------------- call addfld_to(compatm, 'So_ofrac') + !---------------------------------------------------------- + ! from med: ocean albedos (not sent to the ATM in UFS). + !---------------------------------------------------------- + if (trim(coupling_mode) == 'hafs.mom6') then + if (phase == 'advertise') then + call addfld_ocnalb('So_avsdr') + call addfld_ocnalb('So_avsdf') + call addfld_ocnalb('So_anidr') + call addfld_ocnalb('So_anidf') + end if + end if + !===================================================================== ! FIELDS TO ATMOSPHERE !===================================================================== @@ -154,28 +167,41 @@ subroutine esmFldsExchange_hafs_advt(gcomp, phase, rc) ! to atm: surface temperatures from ocn ! --------------------------------------------------------------------- if (hafs_attr%atm_present .and. hafs_attr%ocn_present) then - allocate(S_flds(1)) - S_flds = (/'So_t'/) ! sea_surface_temperature - do n = 1,size(S_flds) - fldname = trim(S_flds(n)) - call addfld_from(compocn, trim(fldname)) - call addfld_to(compatm, trim(fldname)) - end do - deallocate(S_flds) + if (trim(coupling_mode) == 'hafs') then + allocate(S_flds(1)) + S_flds = (/'So_t'/) ! sea_surface_temperature + do n = 1,size(S_flds) + fldname = trim(S_flds(n)) + call addfld_from(compocn, trim(fldname)) + call addfld_to(compatm, trim(fldname)) + end do + deallocate(S_flds) + else + allocate(S_flds(3)) + S_flds = (/'So_t', & ! sea_surface_temperature + 'So_u', & ! surface zonal current + 'So_v'/) ! surface meridional current + do n = 1,size(S_flds) + fldname = trim(S_flds(n)) + call addfld_from(compocn, trim(fldname)) + call addfld_to(compatm, trim(fldname)) + end do + deallocate(S_flds) + end if end if ! --------------------------------------------------------------------- ! to atm: surface roughness length ! --------------------------------------------------------------------- if (hafs_attr%atm_present .and. hafs_attr%wav_present) then - allocate(S_flds(1)) - S_flds = (/'Sw_z0'/) ! wave_z0_roughness_length - do n = 1,size(S_flds) - fldname = trim(S_flds(n)) - call addfld_from(compwav, trim(fldname)) - call addfld_to(compatm, trim(fldname)) - end do - deallocate(S_flds) + allocate(S_flds(1)) + S_flds = (/'Sw_z0'/) ! wave_z0_roughness_length + do n = 1,size(S_flds) + fldname = trim(S_flds(n)) + call addfld_from(compwav, trim(fldname)) + call addfld_to(compatm, trim(fldname)) + end do + deallocate(S_flds) end if !===================================================================== @@ -186,40 +212,72 @@ subroutine esmFldsExchange_hafs_advt(gcomp, phase, rc) ! to ocn: state fields ! --------------------------------------------------------------------- if (hafs_attr%atm_present .and. hafs_attr%ocn_present) then - allocate(S_flds(6)) - S_flds = (/'Sa_u10m', & ! inst_zonal_wind_height10m - 'Sa_v10m', & ! inst_merid_wind_height10m - 'Sa_t2m ', & ! inst_temp_height2m - 'Sa_q2m ', & ! inst_spec_humid_height2m - 'Sa_pslv', & ! inst_pres_height_surface - 'Sa_tskn' /) ! inst_temp_height_surface - do n = 1,size(S_flds) - fldname = trim(S_flds(n)) - call addfld_from(compatm, trim(fldname)) - call addfld_to(compocn, trim(fldname)) - end do - deallocate(S_flds) + if (trim(coupling_mode) == 'hafs') then + allocate(S_flds(6)) + S_flds = (/'Sa_u10m', & ! inst_zonal_wind_height10m + 'Sa_v10m', & ! inst_merid_wind_height10m + 'Sa_t2m ', & ! inst_temp_height2m + 'Sa_q2m ', & ! inst_spec_humid_height2m + 'Sa_pslv', & ! inst_pres_height_surface + 'Sa_tskn' /) ! inst_temp_height_surface + do n = 1,size(S_flds) + fldname = trim(S_flds(n)) + call addfld_from(compatm, trim(fldname)) + call addfld_to(compocn, trim(fldname)) + end do + deallocate(S_flds) + else + allocate(S_flds(1)) + S_flds = (/'Sa_pslv'/) ! inst_pres_height_surface + do n = 1,size(S_flds) + fldname = trim(S_flds(n)) + call addfld_from(compatm, trim(fldname)) + call addfld_to(compocn, trim(fldname)) + end do + deallocate(S_flds) + end if end if ! --------------------------------------------------------------------- ! to ocn: flux fields ! --------------------------------------------------------------------- if (hafs_attr%atm_present .and. hafs_attr%ocn_present) then - allocate(F_flds(7,2)) - F_flds(1,:) = (/'Faxa_taux ','Faxa_taux '/) ! mean_zonal_moment_flx_atm - F_flds(2,:) = (/'Faxa_tauy ','Faxa_tauy '/) ! mean_merid_moment_flx_atm - F_flds(3,:) = (/'Faxa_rain ','Faxa_rain '/) ! mean_prec_rate - F_flds(4,:) = (/'Faxa_swnet','Faxa_swnet'/) ! mean_net_sw_flx - F_flds(5,:) = (/'Faxa_lwnet','Faxa_lwnet'/) ! mean_net_lw_flx - F_flds(6,:) = (/'Faxa_sen ','Faxa_sen '/) ! mean_sensi_heat_flx - F_flds(7,:) = (/'Faxa_lat ','Faxa_lat '/) ! mean_laten_heat_flx - do n = 1,size(F_flds,1) - fldname1 = trim(F_flds(n,1)) - fldname2 = trim(F_flds(n,2)) - call addfld_from(compatm, trim(fldname1)) - call addfld_to(compocn, trim(fldname2)) - end do - deallocate(F_flds) + if (trim(coupling_mode) == 'hafs') then + allocate(F_flds(7,2)) + F_flds(1,:) = (/'Faxa_taux ','Faxa_taux '/) ! inst_zonal_moment_flx_atm + F_flds(2,:) = (/'Faxa_tauy ','Faxa_tauy '/) ! inst_merid_moment_flx_atm + F_flds(3,:) = (/'Faxa_rain ','Faxa_rain '/) ! inst_prec_rate + F_flds(4,:) = (/'Faxa_swnet','Faxa_swnet'/) ! inst_net_sw_flx + F_flds(5,:) = (/'Faxa_lwnet','Faxa_lwnet'/) ! inst_net_lw_flx + F_flds(6,:) = (/'Faxa_sen ','Faxa_sen '/) ! inst_sensi_heat_flx + F_flds(7,:) = (/'Faxa_lat ','Faxa_lat '/) ! inst_laten_heat_flx + do n = 1,size(F_flds,1) + fldname1 = trim(F_flds(n,1)) + fldname2 = trim(F_flds(n,2)) + call addfld_from(compatm, trim(fldname1)) + call addfld_to(compocn, trim(fldname2)) + end do + deallocate(F_flds) + else + allocate(F_flds(10,2)) + F_flds(1 ,:) = (/'Faxa_taux ','Foxx_taux '/) ! inst_zonal_moment_flx_atm + F_flds(2 ,:) = (/'Faxa_tauy ','Foxx_tauy '/) ! inst_merid_moment_flx_atm + F_flds(3 ,:) = (/'Faxa_rain ','Faxa_rain '/) ! inst_prec_rate + F_flds(4 ,:) = (/'Faxa_lwnet ','Foxx_lwnet '/) ! inst_net_lw_flx + F_flds(5 ,:) = (/'Faxa_sen ','Foxx_sen '/) ! inst_sensi_heat_flx + F_flds(6 ,:) = (/'Faxa_evap ','Foxx_evap '/) ! inst_evap_rate + F_flds(7 ,:) = (/'Faxa_swndr ','Foxx_swnet_idr'/) ! inst_down_sw_ir_dir_flx + F_flds(8 ,:) = (/'Faxa_swndf ','Foxx_swnet_idf'/) ! inst_down_sw_ir_dif_flx + F_flds(9 ,:) = (/'Faxa_swvdr ','Foxx_swnet_vdr'/) ! inst_down_sw_vis_dir_flx + F_flds(10,:) = (/'Faxa_swvdf ','Foxx_swnet_vdf'/) ! inst_down_sw_vis_dif_flx + do n = 1,size(F_flds,1) + fldname1 = trim(F_flds(n,1)) + fldname2 = trim(F_flds(n,2)) + call addfld_from(compatm, trim(fldname1)) + call addfld_to(compocn, trim(fldname2)) + end do + deallocate(F_flds) + end if end if !===================================================================== @@ -230,14 +288,14 @@ subroutine esmFldsExchange_hafs_advt(gcomp, phase, rc) ! to wav: 10-m wind components ! --------------------------------------------------------------------- if (hafs_attr%atm_present .and. hafs_attr%wav_present) then - allocate(S_flds(2)) - S_flds = (/'Sa_u10m', 'Sa_v10m'/) - do n = 1,size(S_flds) - fldname = trim(S_flds(n)) - call addfld_from(compatm, trim(fldname)) - call addfld_to(compwav, trim(fldname)) - end do - deallocate(S_flds) + allocate(S_flds(2)) + S_flds = (/'Sa_u10m', 'Sa_v10m'/) + do n = 1,size(S_flds) + fldname = trim(S_flds(n)) + call addfld_from(compatm, trim(fldname)) + call addfld_to(compwav, trim(fldname)) + end do + deallocate(S_flds) end if call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO) @@ -356,40 +414,59 @@ subroutine esmFldsExchange_hafs_init(gcomp, phase, rc) ! to atm: sea surface temperature ! --------------------------------------------------------------------- if (hafs_attr%atm_present .and. hafs_attr%ocn_present) then - allocate(S_flds(1)) - S_flds = (/'So_t'/) ! sea_surface_temperature - do n = 1,size(S_flds) - fldname = trim(S_flds(n)) - if (fldchk(is_local%wrap%FBExp(compatm),trim(fldname),rc=rc) .and. & - fldchk(is_local%wrap%FBImp(compocn,compocn),trim(fldname),rc=rc) & - ) then - call addmap_from(compocn, trim(fldname), compatm, & - mapfillv_bilnr, hafs_attr%mapnorm, hafs_attr%ocn2atm_smap) - call addmrg_to(compatm, trim(fldname), & - mrg_from=compocn, mrg_fld=trim(fldname), mrg_type='copy') - end if - end do - deallocate(S_flds) + if (trim(coupling_mode) == 'hafs') then + allocate(S_flds(1)) + S_flds = (/'So_t'/) ! sea_surface_temperature + do n = 1,size(S_flds) + fldname = trim(S_flds(n)) + if (fldchk(is_local%wrap%FBExp(compatm),trim(fldname),rc=rc) .and. & + fldchk(is_local%wrap%FBImp(compocn,compocn),trim(fldname),rc=rc) & + ) then + call addmap_from(compocn, trim(fldname), compatm, & + mapfillv_bilnr, hafs_attr%mapnorm, hafs_attr%ocn2atm_smap) + call addmrg_to(compatm, trim(fldname), & + mrg_from=compocn, mrg_fld=trim(fldname), mrg_type='copy') + end if + end do + deallocate(S_flds) + else + allocate(S_flds(3)) + S_flds = (/'So_t', & ! sea_surface_temperature + 'So_u', & ! surface zonal current + 'So_v'/) ! surface meridional current + do n = 1,size(S_flds) + fldname = trim(S_flds(n)) + if (fldchk(is_local%wrap%FBExp(compatm),trim(fldname),rc=rc) .and. & + fldchk(is_local%wrap%FBImp(compocn,compocn),trim(fldname),rc=rc) & + ) then + call addmap_from(compocn, trim(fldname), compatm, & + mapfillv_bilnr, hafs_attr%mapnorm, hafs_attr%ocn2atm_smap) + call addmrg_to(compatm, trim(fldname), & + mrg_from=compocn, mrg_fld=trim(fldname), mrg_type='copy') + end if + end do + deallocate(S_flds) + end if end if ! --------------------------------------------------------------------- ! to atm: surface roughness length ! --------------------------------------------------------------------- if (hafs_attr%atm_present .and. hafs_attr%wav_present) then - allocate(S_flds(1)) - S_flds = (/'Sw_z0'/) ! wave_z0_roughness_length - do n = 1,size(S_flds) - fldname = trim(S_flds(n)) - if (fldchk(is_local%wrap%FBExp(compatm),trim(fldname),rc=rc) .and. & - fldchk(is_local%wrap%FBImp(compwav,compwav),trim(fldname),rc=rc) & - ) then - call addmap_from(compwav, trim(fldname), compatm, & - mapfillv_bilnr, hafs_attr%mapnorm, hafs_attr%wav2atm_smap) - call addmrg_to(compatm, trim(fldname), & - mrg_from=compwav, mrg_fld=trim(fldname), mrg_type='copy') - end if - end do - deallocate(S_flds) + allocate(S_flds(1)) + S_flds = (/'Sw_z0'/) ! wave_z0_roughness_length + do n = 1,size(S_flds) + fldname = trim(S_flds(n)) + if (fldchk(is_local%wrap%FBExp(compatm),trim(fldname),rc=rc) .and. & + fldchk(is_local%wrap%FBImp(compwav,compwav),trim(fldname),rc=rc) & + ) then + call addmap_from(compwav, trim(fldname), compatm, & + mapfillv_bilnr, hafs_attr%mapnorm, hafs_attr%wav2atm_smap) + call addmrg_to(compatm, trim(fldname), & + mrg_from=compwav, mrg_fld=trim(fldname), mrg_type='copy') + end if + end do + deallocate(S_flds) end if !===================================================================== @@ -400,52 +477,96 @@ subroutine esmFldsExchange_hafs_init(gcomp, phase, rc) ! to ocn: state fields ! --------------------------------------------------------------------- if (hafs_attr%atm_present .and. hafs_attr%ocn_present) then - allocate(S_flds(6)) - S_flds = (/'Sa_u10m', & ! inst_zonal_wind_height10m - 'Sa_v10m', & ! inst_merid_wind_height10m - 'Sa_t2m ', & ! inst_temp_height2m - 'Sa_q2m ', & ! inst_spec_humid_height2m - 'Sa_pslv', & ! inst_pres_height_surface - 'Sa_tskn' /) ! inst_temp_height_surface - do n = 1,size(S_flds) - fldname = trim(S_flds(n)) - if (fldchk(is_local%wrap%FBExp(compocn),trim(fldname),rc=rc) .and. & - fldchk(is_local%wrap%FBImp(compatm,compatm),trim(fldname),rc=rc) & - ) then - call addmap_from(compatm, trim(fldname), compocn, & - mapfillv_bilnr, hafs_attr%mapnorm, hafs_attr%atm2ocn_smap) - call addmrg_to(compocn, trim(fldname), & - mrg_from=compatm, mrg_fld=trim(fldname), mrg_type='copy') - end if - end do - deallocate(S_flds) + if (trim(coupling_mode) == 'hafs') then + allocate(S_flds(6)) + S_flds = (/'Sa_u10m', & ! inst_zonal_wind_height10m + 'Sa_v10m', & ! inst_merid_wind_height10m + 'Sa_t2m ', & ! inst_temp_height2m + 'Sa_q2m ', & ! inst_spec_humid_height2m + 'Sa_pslv', & ! inst_pres_height_surface + 'Sa_tskn' /) ! inst_temp_height_surface + do n = 1,size(S_flds) + fldname = trim(S_flds(n)) + if (fldchk(is_local%wrap%FBExp(compocn),trim(fldname),rc=rc) .and. & + fldchk(is_local%wrap%FBImp(compatm,compatm),trim(fldname),rc=rc) & + ) then + call addmap_from(compatm, trim(fldname), compocn, & + mapfillv_bilnr, hafs_attr%mapnorm, hafs_attr%atm2ocn_smap) + call addmrg_to(compocn, trim(fldname), & + mrg_from=compatm, mrg_fld=trim(fldname), mrg_type='copy') + end if + end do + deallocate(S_flds) + else + allocate(S_flds(1)) + S_flds = (/'Sa_pslv'/) ! inst_pres_height_surface + do n = 1,size(S_flds) + fldname = trim(S_flds(n)) + if (fldchk(is_local%wrap%FBExp(compocn),trim(fldname),rc=rc) .and. & + fldchk(is_local%wrap%FBImp(compatm,compatm),trim(fldname),rc=rc) & + ) then + call addmap_from(compatm, trim(fldname), compocn, & + mapfillv_bilnr, hafs_attr%mapnorm, hafs_attr%atm2ocn_smap) + call addmrg_to(compocn, trim(fldname), & + mrg_from=compatm, mrg_fld=trim(fldname), mrg_type='copy') + end if + end do + deallocate(S_flds) + end if end if ! --------------------------------------------------------------------- ! to ocn: flux fields ! --------------------------------------------------------------------- if (hafs_attr%atm_present .and. hafs_attr%ocn_present) then - allocate(F_flds(7,2)) - F_flds(1,:) = (/'Faxa_taux ','Faxa_taux '/) ! mean_zonal_moment_flx_atm - F_flds(2,:) = (/'Faxa_tauy ','Faxa_tauy '/) ! mean_merid_moment_flx_atm - F_flds(3,:) = (/'Faxa_rain ','Faxa_rain '/) ! mean_prec_rate - F_flds(4,:) = (/'Faxa_swnet','Faxa_swnet'/) ! mean_net_sw_flx - F_flds(5,:) = (/'Faxa_lwnet','Faxa_lwnet'/) ! mean_net_lw_flx - F_flds(6,:) = (/'Faxa_sen ','Faxa_sen '/) ! mean_sensi_heat_flx - F_flds(7,:) = (/'Faxa_lat ','Faxa_lat '/) ! mean_laten_heat_flx - do n = 1,size(F_flds,1) - fldname1 = trim(F_flds(n,1)) - fldname2 = trim(F_flds(n,2)) - if (fldchk(is_local%wrap%FBExp(compocn),trim(fldname2),rc=rc) .and. & - fldchk(is_local%wrap%FBImp(compatm,compatm),trim(fldname1),rc=rc) & - ) then - call addmap_from(compatm, trim(fldname1), compocn, & - mapfillv_bilnr, hafs_attr%mapnorm, hafs_attr%atm2ocn_smap) - call addmrg_to(compocn, trim(fldname2), & - mrg_from=compatm, mrg_fld=trim(fldname1), mrg_type='copy') - end if - end do - deallocate(F_flds) + if (trim(coupling_mode) == 'hafs') then + allocate(F_flds(7,2)) + F_flds(1,:) = (/'Faxa_taux ','Faxa_taux '/) ! inst_zonal_moment_flx_atm + F_flds(2,:) = (/'Faxa_tauy ','Faxa_tauy '/) ! inst_merid_moment_flx_atm + F_flds(3,:) = (/'Faxa_rain ','Faxa_rain '/) ! inst_prec_rate + F_flds(4,:) = (/'Faxa_swnet','Faxa_swnet'/) ! inst_net_sw_flx + F_flds(5,:) = (/'Faxa_lwnet','Faxa_lwnet'/) ! inst_net_lw_flx + F_flds(6,:) = (/'Faxa_sen ','Faxa_sen '/) ! inst_sensi_heat_flx + F_flds(7,:) = (/'Faxa_lat ','Faxa_lat '/) ! inst_laten_heat_flx + do n = 1,size(F_flds,1) + fldname1 = trim(F_flds(n,1)) + fldname2 = trim(F_flds(n,2)) + if (fldchk(is_local%wrap%FBExp(compocn),trim(fldname2),rc=rc) .and. & + fldchk(is_local%wrap%FBImp(compatm,compatm),trim(fldname1),rc=rc) & + ) then + call addmap_from(compatm, trim(fldname1), compocn, & + mapfillv_bilnr, hafs_attr%mapnorm, hafs_attr%atm2ocn_smap) + call addmrg_to(compocn, trim(fldname2), & + mrg_from=compatm, mrg_fld=trim(fldname1), mrg_type='copy') + end if + end do + deallocate(F_flds) + else + allocate(F_flds(10,2)) + F_flds(1 ,:) = (/'Faxa_taux ','Foxx_taux '/) ! inst_zonal_moment_flx_atm + F_flds(2 ,:) = (/'Faxa_tauy ','Foxx_tauy '/) ! inst_merid_moment_flx_atm + F_flds(3 ,:) = (/'Faxa_rain ','Faxa_rain '/) ! inst_prec_rate + F_flds(4 ,:) = (/'Faxa_lwnet ','Foxx_lwnet '/) ! inst_net_lw_flx + F_flds(5 ,:) = (/'Faxa_sen ','Foxx_sen '/) ! inst_sensi_heat_flx + F_flds(6 ,:) = (/'Faxa_evap ','Foxx_evap '/) ! inst_evap_rate + F_flds(7 ,:) = (/'Faxa_swndr ','Foxx_swnet_idr'/) ! inst_down_sw_ir_dir_flx + F_flds(8 ,:) = (/'Faxa_swndf ','Foxx_swnet_idf'/) ! inst_down_sw_ir_dif_flx + F_flds(9 ,:) = (/'Faxa_swvdr ','Foxx_swnet_vdr'/) ! inst_down_sw_vis_dir_flx + F_flds(10,:) = (/'Faxa_swvdf ','Foxx_swnet_vdf'/) ! inst_down_sw_vis_dif_flx + do n = 1,size(F_flds,1) + fldname1 = trim(F_flds(n,1)) + fldname2 = trim(F_flds(n,2)) + if (fldchk(is_local%wrap%FBExp(compocn),trim(fldname2),rc=rc) .and. & + fldchk(is_local%wrap%FBImp(compatm,compatm),trim(fldname1),rc=rc) & + ) then + call addmap_from(compatm, trim(fldname1), compocn, & + mapfillv_bilnr, hafs_attr%mapnorm, hafs_attr%atm2ocn_smap) + call addmrg_to(compocn, trim(fldname2), & + mrg_from=compatm, mrg_fld=trim(fldname1), mrg_type='copy') + end if + end do + deallocate(F_flds) + end if end if !===================================================================== diff --git a/mediator/esmFldsExchange_nems_mod.F90 b/mediator/esmFldsExchange_ufs_mod.F90 similarity index 92% rename from mediator/esmFldsExchange_nems_mod.F90 rename to mediator/esmFldsExchange_ufs_mod.F90 index b5d3fbf8f..aa8088306 100644 --- a/mediator/esmFldsExchange_nems_mod.F90 +++ b/mediator/esmFldsExchange_ufs_mod.F90 @@ -1,4 +1,4 @@ -module esmFldsExchange_nems_mod +module esmFldsExchange_ufs_mod !--------------------------------------------------------------------- ! This is a mediator specific routine that determines ALL possible @@ -9,7 +9,7 @@ module esmFldsExchange_nems_mod implicit none public - public :: esmFldsExchange_nems + public :: esmFldsExchange_ufs character(*), parameter :: u_FILE_u = & __FILE__ @@ -18,7 +18,7 @@ module esmFldsExchange_nems_mod contains !================================================================================ - subroutine esmFldsExchange_nems(gcomp, phase, rc) + subroutine esmFldsExchange_ufs(gcomp, phase, rc) use ESMF use NUOPC @@ -54,7 +54,7 @@ subroutine esmFldsExchange_nems(gcomp, phase, rc) character(len=CL) :: cvalue character(len=CS) :: fldname character(len=CS), allocatable :: flds(:), oflds(:), aflds(:), iflds(:) - character(len=*) , parameter :: subname='(esmFldsExchange_nems)' + character(len=*) , parameter :: subname='(esmFldsExchange_ufs)' !-------------------------------------- rc = ESMF_SUCCESS @@ -68,7 +68,7 @@ subroutine esmFldsExchange_nems(gcomp, phase, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return ! Set maptype according to coupling_mode - if (trim(coupling_mode) == 'nems_orig' .or. trim(coupling_mode) == 'nems_orig_data') then + if (trim(coupling_mode) == 'ufs.nfrac' .or. trim(coupling_mode) == 'ufs.nfrac.aoflux') then maptype = mapnstod_consf else maptype = mapconsf @@ -76,7 +76,7 @@ subroutine esmFldsExchange_nems(gcomp, phase, rc) write(msgString,'(A,i6,A)') trim(subname)//': maptype is ',maptype,', '//mapnames(maptype) call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO) - if (trim(coupling_mode) == 'nems_orig_data' .or. trim(coupling_mode) == 'nems_frac_aoflux') then + if (trim(coupling_mode) == 'ufs.nfrac.aoflux' .or. trim(coupling_mode) == 'ufs.frac.aoflux') then med_aoflux_to_ocn = .true. else med_aoflux_to_ocn = .false. @@ -151,14 +151,6 @@ subroutine esmFldsExchange_nems(gcomp, phase, rc) call addfld_ocnalb('So_anidf') end if - ! Advertise the ocean albedos. These are not sent to the ATM in UFS. - if (phase == 'advertise') then - call addfld_ocnalb('So_avsdr') - call addfld_ocnalb('So_avsdf') - call addfld_ocnalb('So_anidr') - call addfld_ocnalb('So_anidf') - end if - !===================================================================== ! FIELDS TO ATMOSPHERE (compatm) !===================================================================== @@ -242,18 +234,46 @@ subroutine esmFldsExchange_nems(gcomp, phase, rc) end if end if - ! to atm: unmerged surface temperatures from lnd - if (phase == 'advertise') then - if (is_local%wrap%comp_present(complnd) .and. is_local%wrap%comp_present(compatm)) then - call addfld_from(complnd , 'Sl_t') - call addfld_to(compatm , 'Sl_t') + ! to atm: unmerged flux components from lnd + if (is_local%wrap%comp_present(complnd) .and. is_local%wrap%comp_present(compatm)) then + allocate(flds(6)) + flds = (/ 'lat ', 'sen ', 'evap', 'gflx', 'roff', 'soff' /) + if (phase == 'advertise') then + do n = 1,size(flds) + call addfld_from(complnd, 'Fall_'//trim(flds(n))) + call addfld_to(compatm, 'Fall_'//trim(flds(n))) + end do + else + do n = 1,size(flds) + if ( fldchk(is_local%wrap%FBexp(compatm) , 'Fall_'//trim(flds(n)), rc=rc) .and. & + fldchk(is_local%wrap%FBImp(complnd,complnd), 'Fall_'//trim(flds(n)), rc=rc)) then + call addmap_from(complnd, 'Fall_'//trim(flds(n)), compatm, maptype, 'lfrac', 'unset') + call addmrg_to(compatm, 'Fall_'//trim(flds(n)), mrg_from=complnd, mrg_fld='Fall_'//trim(flds(n)), mrg_type='copy') + end if + end do end if - else - if ( fldchk(is_local%wrap%FBexp(compatm) , 'Sl_t', rc=rc) .and. & - fldchk(is_local%wrap%FBImp(complnd,complnd), 'Sl_t', rc=rc)) then - call addmap_from(complnd, 'Sl_t', compatm, maptype, 'lfrin', 'unset') - call addmrg_to(compatm, 'Sl_t', mrg_from=complnd, mrg_fld='Sl_t', mrg_type='copy') + deallocate(flds) + end if + + ! to atm: unmerged state variables from lnd + if (is_local%wrap%comp_present(complnd) .and. is_local%wrap%comp_present(compatm)) then + allocate(flds(7)) + flds = (/ 'sfrac', 'tref ', 'qref ', 'q ', 'cmm ', 'chh ', 'zvfun' /) + if (phase == 'advertise') then + do n = 1,size(flds) + call addfld_from(complnd, 'Sl_'//trim(flds(n))) + call addfld_to(compatm, 'Sl_'//trim(flds(n))) + end do + else + do n = 1,size(flds) + if ( fldchk(is_local%wrap%FBexp(compatm) , 'Sl_'//trim(flds(n)), rc=rc) .and. & + fldchk(is_local%wrap%FBImp(complnd,complnd), 'Sl_'//trim(flds(n)), rc=rc)) then + call addmap_from(complnd, 'Sl_'//trim(flds(n)), compatm, maptype, 'lfrac', 'unset') + call addmrg_to(compatm, 'Sl_'//trim(flds(n)), mrg_from=complnd, mrg_fld='Sl_'//trim(flds(n)), mrg_type='copy') + end if + end do end if + deallocate(flds) end if ! to atm: unmerged from mediator, merge will be done under FV3/CCPP composite step @@ -721,21 +741,19 @@ subroutine esmFldsExchange_nems(gcomp, phase, rc) !===================================================================== ! to lnd - states and fluxes from atm - if ( trim(coupling_mode) == 'nems_orig_data') then + if ( trim(coupling_mode) == 'ufs.nfrac.aoflux') then allocate(flds(21)) flds = (/'Sa_z ', 'Sa_topo ', 'Sa_tbot ', 'Sa_pbot ', & - 'Sa_shum ', 'Sa_u ', 'Sa_v ', 'Faxa_lwdn ', & - 'Sa_ptem ', 'Sa_dens ', 'Faxa_swdn ', 'Sa_pslv ', & - 'Faxa_snowc', 'Faxa_snowl', 'Faxa_rainc', 'Faxa_rainl', & - 'Faxa_swndr', 'Faxa_swndf', 'Faxa_swvdr', 'Faxa_swvdf', & - 'Faxa_swnet'/) + 'Sa_shum ', 'Sa_u ', 'Sa_v ', 'Sa_pslv ', & + 'Faxa_lwdn ', 'Faxa_swdn ', 'Faxa_snowc', 'Faxa_snowl', & + 'Faxa_rainc', 'Faxa_rainl', 'Faxa_rain ', 'Faxa_swnet'/) else allocate(flds(18)) flds = (/'Sa_z ', 'Sa_ta ', 'Sa_pslv ', 'Sa_qa ', & - 'Sa_ua ', 'Sa_va ', 'Faxa_swdn ', 'Faxa_lwdn ', & - 'Faxa_swnet', 'Faxa_rain ', 'Sa_prsl ', 'vfrac ', & + 'Sa_u ', 'Sa_v ', 'Faxa_swdn ', 'Faxa_lwdn ', & + 'Faxa_swnet', 'Faxa_rain ', 'Sa_prsl ', 'Sa_vfrac ', & 'Faxa_snow ', 'Faxa_rainc', 'Sa_tskn ', 'Sa_exner ', & - 'Sa_ustar ', 'zorl ' /) + 'Sa_ustar ', 'Sa_zorl ' /) end if do n = 1,size(flds) fldname = trim(flds(n)) @@ -754,6 +772,6 @@ subroutine esmFldsExchange_nems(gcomp, phase, rc) end do deallocate(flds) - end subroutine esmFldsExchange_nems + end subroutine esmFldsExchange_ufs -end module esmFldsExchange_nems_mod +end module esmFldsExchange_ufs_mod diff --git a/mediator/fd_cesm.yaml b/mediator/fd_cesm.yaml index 06af37ed4..748c37711 100644 --- a/mediator/fd_cesm.yaml +++ b/mediator/fd_cesm.yaml @@ -522,6 +522,10 @@ canonical_units: m description: atmosphere import # + - standard_name: So_ugustOut + canonical_units: m/s + description: atmosphere import + # #----------------------------------- # section: land-ice export # Note that the fields sent from glc->med do NOT have elevation classes, diff --git a/mediator/med.F90 b/mediator/med.F90 index 3cd958d5c..087db69bc 100644 --- a/mediator/med.F90 +++ b/mediator/med.F90 @@ -48,7 +48,7 @@ module MED use esmFlds , only : med_fldList_GetNumFlds, med_fldList_GetFldNames, med_fldList_GetFldInfo use esmFlds , only : med_fldList_Document_Mapping, med_fldList_Document_Merging use esmFlds , only : med_fldList_GetfldListFr, med_fldList_GetfldListTo, med_fldList_Realize - use esmFldsExchange_nems_mod , only : esmFldsExchange_nems + use esmFldsExchange_ufs_mod , only : esmFldsExchange_ufs use esmFldsExchange_cesm_mod , only : esmFldsExchange_cesm use esmFldsExchange_hafs_mod , only : esmFldsExchange_hafs use med_phases_profile_mod , only : med_phases_profile_finalize @@ -123,6 +123,9 @@ subroutine SetServices(gcomp, rc) use med_diag_mod , only: med_phases_diag_ice_ice2med, med_phases_diag_ice_med2ice use med_fraction_mod , only: med_fraction_init, med_fraction_set use med_phases_profile_mod , only: med_phases_profile +#ifdef CDEPS_INLINE + use med_phases_cdeps_mod , only: med_phases_cdeps_run +#endif ! input/output variables type(ESMF_GridComp) :: gcomp @@ -505,6 +508,19 @@ subroutine SetServices(gcomp, rc) specPhaselabel="med_phases_diag_print", specRoutine=NUOPC_NoOp, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return +#ifdef CDEPS_INLINE + !------------------ + ! phase routine for cdeps inline capabilty + !------------------ + + call NUOPC_CompSetEntryPoint(gcomp, ESMF_METHOD_RUN, & + phaseLabelList=(/"med_phases_cdeps_run"/), userRoutine=mediator_routine_Run, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call NUOPC_CompSpecialize(gcomp, specLabel=mediator_label_Advance, & + specPhaseLabel="med_phases_cdeps_run", specRoutine=med_phases_cdeps_run, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return +#endif + !------------------ ! attach specializing method(s) ! -> NUOPC specializes by default --->>> first need to remove the default @@ -817,8 +833,8 @@ subroutine AdvertiseFields(gcomp, importState, exportState, clock, rc) if (trim(coupling_mode) == 'cesm') then call esmFldsExchange_cesm(gcomp, phase='advertise', rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - else if (trim(coupling_mode(1:4)) == 'nems') then - call esmFldsExchange_nems(gcomp, phase='advertise', rc=rc) + else if (trim(coupling_mode(1:3)) == 'ufs') then + call esmFldsExchange_ufs(gcomp, phase='advertise', rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return else if (trim(coupling_mode(1:4)) == 'hafs') then call esmFldsExchange_hafs(gcomp, phase='advertise', rc=rc) @@ -935,6 +951,7 @@ subroutine AdvertiseFields(gcomp, importState, exportState, clock, rc) endif endif + ! Should enthalpy be calculated call NUOPC_CompAttributeGet(gcomp, name="compute_enthalpy", value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -951,6 +968,23 @@ subroutine AdvertiseFields(gcomp, importState, exportState, clock, rc) write(logunit,*) ' Enthalpy calculation is OFF' endif endif + + ! Should target component use all data for first time step? + do ncomp = 1,ncomps + if (ncomp /= compmed) then + call NUOPC_CompAttributeGet(gcomp, name=trim(compname(ncomp))//"_use_data_first_import", value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + read(cvalue, *) is_local%wrap%med_data_force_first(ncomp) + else + is_local%wrap%med_data_force_first(ncomp) = .false. + endif + if (maintask) then + write(logunit,*) trim(compname(ncomp))//'_use_data_first_import is ', is_local%wrap%med_data_force_first(ncomp) + endif + end if + end do + if (profile_memory) call ESMF_VMLogMemInfo("Leaving "//trim(subname)) call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO) @@ -1040,7 +1074,6 @@ subroutine ModifyDecompofMesh(gcomp, importState, exportState, clock, rc) type(InternalState) :: is_local integer :: n1 character(len=*), parameter :: subname = '('//__FILE__//':ModifyDecompofMesh)' - !----------------------------------------------------------- call ESMF_LogWrite(trim(subname)//": called", ESMF_LOGMSG_INFO) @@ -1369,7 +1402,6 @@ subroutine RealizeFieldsWithTransferAccept(gcomp, importState, exportState, cloc type(InternalState) :: is_local integer :: n1 character(len=*), parameter :: subname = '('//__FILE__//':RealizeFieldsWithTransferAccept)' - !----------------------------------------------------------- call ESMF_LogWrite(trim(subname)//": called", ESMF_LOGMSG_INFO) @@ -1819,10 +1851,10 @@ subroutine DataInitialize(gcomp, rc) if (trim(coupling_mode) == 'cesm') then call esmFldsExchange_cesm(gcomp, phase='initialize', rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - else if (trim(coupling_mode(1:4)) == 'nems') then - call esmFldsExchange_nems(gcomp, phase='initialize', rc=rc) + else if (trim(coupling_mode(1:3)) == 'ufs') then + call esmFldsExchange_ufs(gcomp, phase='initialize', rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - else if (trim(coupling_mode) == 'hafs') then + else if (coupling_mode(1:4) == 'hafs') then call esmFldsExchange_hafs(gcomp, phase='initialize', rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if @@ -2239,7 +2271,6 @@ subroutine SetRunClock(gcomp, rc) integer :: stop_n, stop_ymd logical, save :: stopalarmcreated=.false. character(len=*), parameter :: subname = '('//__FILE__//':SetRunClock)' - !----------------------------------------------------------- rc = ESMF_SUCCESS diff --git a/mediator/med_fraction_mod.F90 b/mediator/med_fraction_mod.F90 index 7fe0315b6..2f7d43041 100644 --- a/mediator/med_fraction_mod.F90 +++ b/mediator/med_fraction_mod.F90 @@ -293,7 +293,7 @@ subroutine med_fraction_init(gcomp, rc) ! If ice and atm are on the same mesh - a redist route handle has already been created maptype = mapfcopy else - if (trim(coupling_mode) == 'nems_orig' ) then + if (trim(coupling_mode(1:9)) == 'ufs.nfrac' ) then maptype = mapnstod_consd else maptype = mapconsd @@ -345,7 +345,7 @@ subroutine med_fraction_init(gcomp, rc) ! If ocn and atm are on the same mesh - a redist route handle has already been created maptype = mapfcopy else - if (trim(coupling_mode) == 'nems_orig' ) then + if (trim(coupling_mode(1:9)) == 'ufs.nfrac' ) then maptype = mapnstod_consd else maptype = mapconsd @@ -756,7 +756,7 @@ subroutine med_fraction_set(gcomp, rc) call t_startf('MED:'//trim(subname)//' fbfrac(compatm)') ! Determine maptype - if (trim(coupling_mode) == 'nems_orig' ) then + if (trim(coupling_mode(1:9)) == 'ufs.nfrac' ) then maptype = mapnstod_consd else if (med_map_RH_is_created(is_local%wrap%RH(compice,compatm,:),mapfcopy, rc=rc)) then diff --git a/mediator/med_internalstate_mod.F90 b/mediator/med_internalstate_mod.F90 index 9d6029b1f..c03073b65 100644 --- a/mediator/med_internalstate_mod.F90 +++ b/mediator/med_internalstate_mod.F90 @@ -19,7 +19,7 @@ module med_internalstate_mod integer, public :: logunit ! logunit for mediator log output integer, public :: diagunit ! diagunit for budget output (med main only) - logical, public :: maintask=.false. ! is this the maintask + logical, public :: maintask = .false. ! is this the maintask integer, public :: med_id ! needed currently in med_io_mod and set in esm.F90 ! Components @@ -47,7 +47,7 @@ module med_internalstate_mod character(len=CS), public :: glc_name = '' ! Coupling mode - character(len=CS), public :: coupling_mode ! valid values are [cesm,nems_orig,nems_frac,nems_orig_data,hafs,nems_frac_aoflux,nems_frac_aoflux_sbs] + character(len=CS), public :: coupling_mode ! valid values are [cesm,ufs.nfrac,ufs.frac,ufs.nfrac.aoflux,ufs.frac.aoflux,hafs,hafs.mom6] ! Atmosphere-ocean flux algorithm character(len=CS), public :: aoflux_code ! valid values are [cesm,ccpp] @@ -119,10 +119,12 @@ module med_internalstate_mod type InternalStateStruct ! Present/allowed coupling/active coupling logical flags - logical, pointer :: comp_present(:) ! comp present flag - logical, pointer :: med_coupling_active(:,:) ! computes the active coupling - integer :: num_icesheets ! obtained from attribute - logical :: ocn2glc_coupling = .false. ! obtained from attribute + logical, pointer :: comp_present(:) ! comp present flag + logical, pointer :: med_coupling_active(:,:) ! computes the active coupling + logical, pointer :: med_data_active(:,:) ! uses stream data to provide background fill + logical, pointer :: med_data_force_first(:) ! force to use stream data for first coupling timestep + integer :: num_icesheets ! obtained from attribute + logical :: ocn2glc_coupling = .false. ! obtained from attribute logical :: lnd2glc_coupling = .false. logical :: accum_lnd2glc = .false. logical :: docn_present ! aoflux calc requires med_coupling_active true even for docn @@ -148,10 +150,10 @@ module med_internalstate_mod ! FBImp(n,n) = NState_Imp(n), copied in connector post phase ! FBImp(n,k) is the FBImp(n,n) interpolated to grid k ! Import/export States and field bundles (the field bundles have the scalar fields removed) - type(ESMF_State) , pointer :: NStateImp(:) ! Import data from various component, on their grid - type(ESMF_State) , pointer :: NStateExp(:) ! Export data to various component, on their grid - type(ESMF_FieldBundle) , pointer :: FBImp(:,:) ! Import data from various components interpolated to various grids - type(ESMF_FieldBundle) , pointer :: FBExp(:) ! Export data for various components, on their grid + type(ESMF_State) , pointer :: NStateImp(:) ! Import data from various component, on their grid + type(ESMF_State) , pointer :: NStateExp(:) ! Export data to various component, on their grid + type(ESMF_FieldBundle) , pointer :: FBImp(:,:) ! Import data from various components interpolated to various grids + type(ESMF_FieldBundle) , pointer :: FBExp(:) ! Export data for various components, on their grid ! Mediator field bundles for ocean albedo type(ESMF_FieldBundle) :: FBMed_ocnalb_o ! Ocn albedo on ocn grid @@ -174,6 +176,9 @@ module med_internalstate_mod ! Fractions type(ESMF_FieldBundle), pointer :: FBfrac(:) ! Fraction data for various components, on their grid + ! Data + type(ESMF_FieldBundle) , pointer :: FBData(:) ! Background data for various components, on their grid, provided by CDEPS inline + ! Accumulators for export field bundles type(ESMF_FieldBundle) :: FBExpAccumOcn ! Accumulator for Ocn export on Ocn grid integer :: ExpAccumOcnCnt = 0 ! Accumulator counter for FBExpAccumOcn @@ -306,6 +311,8 @@ subroutine med_internalstate_init(gcomp, rc) ! Allocate memory now that ncomps is determined allocate(is_local%wrap%med_coupling_active(ncomps,ncomps)) + allocate(is_local%wrap%med_data_active(ncomps,ncomps)) + allocate(is_local%wrap%med_data_force_first(ncomps)) allocate(is_local%wrap%nx(ncomps)) allocate(is_local%wrap%ny(ncomps)) allocate(is_local%wrap%NStateImp(ncomps)) @@ -319,6 +326,7 @@ subroutine med_internalstate_init(gcomp, rc) allocate(is_local%wrap%packed_data(ncomps,ncomps,nmappers)) allocate(is_local%wrap%FBfrac(ncomps)) allocate(is_local%wrap%FBArea(ncomps)) + allocate(is_local%wrap%FBData(ncomps)) allocate(is_local%wrap%mesh_info(ncomps)) ! Determine component names @@ -367,6 +375,15 @@ subroutine med_internalstate_init(gcomp, rc) write(msgString,*) trim(subname)//': Mediator dststatus_print is ',dststatus_print call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO) + ! Initialize flag for background fill using data + is_local%wrap%med_data_active(:,:) = .false. + is_local%wrap%med_data_active(compocn,compatm) = .true. + is_local%wrap%med_data_active(compatm,compocn) = .true. + is_local%wrap%med_data_active(compatm,compwav) = .true. + + ! Initialize flag to force using data in first coupling time step + is_local%wrap%med_data_force_first(:) = .false. + end subroutine med_internalstate_init !===================================================================== @@ -587,7 +604,7 @@ subroutine med_internalstate_defaultmasks(gcomp, rc) if (is_local%wrap%comp_present(compocn)) defaultMasks(compocn,:) = 0 if (is_local%wrap%comp_present(compice)) defaultMasks(compice,:) = 0 if (is_local%wrap%comp_present(compwav)) defaultMasks(compwav,:) = 0 - if ( trim(coupling_mode(1:4)) == 'nems') then + if ( trim(coupling_mode(1:3)) == 'ufs') then if (is_local%wrap%comp_present(compatm)) defaultMasks(compatm,:) = 1 endif if ( trim(coupling_mode) == 'hafs') then diff --git a/mediator/med_map_mod.F90 b/mediator/med_map_mod.F90 index 54bcbb154..bcf178fbd 100644 --- a/mediator/med_map_mod.F90 +++ b/mediator/med_map_mod.F90 @@ -408,13 +408,13 @@ subroutine med_map_routehandles_initfrom_field(n1, n2, fldsrc, flddst, mapindex, dstMaskValue = ispval_mask endif end if - if (trim(coupling_mode(1:4)) == 'nems') then + if (trim(coupling_mode(1:3)) == 'ufs') then if (n1 == compatm .and. n2 == complnd) then srcMaskValue = ispval_mask dstMaskValue = ispval_mask end if end if - if (trim(coupling_mode) == 'hafs') then + if (coupling_mode(1:4) == 'hafs') then if (n1 == compatm .and. n2 == compwav) then srcMaskValue = ispval_mask end if @@ -424,7 +424,7 @@ subroutine med_map_routehandles_initfrom_field(n1, n2, fldsrc, flddst, mapindex, call ESMF_LogWrite(trim(string), ESMF_LOGMSG_INFO) polemethod=ESMF_POLEMETHOD_ALLAVG - if (trim(coupling_mode) == 'cesm' .or. trim(coupling_mode(1:4)) == 'nems') then + if (trim(coupling_mode) == 'cesm' .or. trim(coupling_mode(1:3)) == 'ufs') then if (n1 == compwav .or. n2 == compwav) then polemethod = ESMF_POLEMETHOD_NONE ! todo: remove this when ESMF tripolar mapping fix is in place. endif @@ -893,12 +893,14 @@ subroutine med_map_packed_field_create(destcomp, flds_scalar_name, & if (npacked(mapindex) > 0) then ! Create the packed source field bundle for mapindex allocate(ptrsrc_packed(npacked(mapindex), lsize_src)) + ptrsrc_packed(npacked(mapindex),:) = 0._R8 packed_data(mapindex)%field_src = ESMF_FieldCreate(lmesh_src, & ptrsrc_packed, gridToFieldMap=(/2/), meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return ! Create the packed destination field bundle for mapindex allocate(ptrdst_packed(npacked(mapindex), lsize_dst)) + ptrdst_packed(npacked(mapindex),:) = 0._R8 packed_data(mapindex)%field_dst = ESMF_FieldCreate(lmesh_dst, & ptrdst_packed, gridToFieldMap=(/2/), meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return @@ -919,7 +921,7 @@ subroutine med_map_packed_field_create(destcomp, flds_scalar_name, & end subroutine med_map_packed_field_create !================================================================================ - subroutine med_map_field_packed(FBSrc, FBDst, FBFracSrc, field_normOne, packed_data, routehandles, rc) + subroutine med_map_field_packed(FBSrc, FBDst, FBFracSrc, FBDat, use_data, field_normOne, packed_data, routehandles, rc) ! ----------------------------------------------- ! Do regridding via packed field bundles @@ -929,33 +931,46 @@ subroutine med_map_field_packed(FBSrc, FBDst, FBFracSrc, field_normOne, packed_d use ESMF , only : ESMF_FieldBundle, ESMF_FieldBundleGet use ESMF , only : ESMF_FieldBundleIsCreated use ESMF , only : ESMF_FieldRedist, ESMF_RouteHandle + use ESMF , only : ESMF_FieldFill + use ESMF , only : ESMF_KIND_R8 + use ESMF , only : ESMF_Region_Flag, ESMF_REGION_SELECT, ESMF_REGION_TOTAL use med_internalstate_mod , only : nmappers, mapfcopy use med_internalstate_mod , only : mappatch_uv3d, mappatch, mapbilnr_uv3d, mapbilnr use med_internalstate_mod , only : packed_data_type + use med_methods_mod , only : Field_diagnose => med_methods_Field_diagnose ! input/output variables - type(ESMF_FieldBundle) , intent(in) :: FBSrc - type(ESMF_FieldBundle) , intent(inout) :: FBDst - type(ESMF_Field) , intent(in) :: field_normOne(:) ! array over mapping types - type(ESMF_FieldBundle) , intent(in) :: FBFracSrc ! fraction field bundle for source - type(packed_data_type) , intent(inout) :: packed_data(:) ! array over mapping types - type(ESMF_RouteHandle) , intent(inout) :: routehandles(:) - integer , intent(out) :: rc + type(ESMF_FieldBundle) , intent(in) :: FBSrc + type(ESMF_FieldBundle) , intent(inout) :: FBDst + type(ESMF_Field) , intent(in) :: field_normOne(:) ! array over mapping types + type(ESMF_FieldBundle) , intent(in) :: FBFracSrc ! fraction field bundle for source + type(packed_data_type) , intent(inout) :: packed_data(:) ! array over mapping types + type(ESMF_RouteHandle) , intent(inout) :: routehandles(:) + type(ESMF_FieldBundle), optional, intent(in) :: FBDat ! data field bundle + logical, optional , intent(in) :: use_data ! skip mapping and use data instead + integer, optional , intent(out) :: rc ! local variables - integer :: nf, nu, np, n - integer :: fieldcount - integer :: mapindex - integer :: ungriddedUBound(1) - real(r8), pointer :: dataptr1d(:) - real(r8), pointer :: dataptr2d(:,:) - real(r8), pointer :: dataptr2d_packed(:,:) - type(ESMF_Field) :: field_fracsrc - type(ESMF_Field), pointer :: fieldlist_src(:) - type(ESMF_Field), pointer :: fieldlist_dst(:) - real(r8), pointer :: data_norm(:) - real(r8), pointer :: data_dst(:,:) - character(len=*), parameter :: subname=' (med_map_mod:med_map_field_packed) ' + integer :: nf, nu, np, n, nfd + integer :: fieldcount, fieldcount_dat + integer :: mapindex + integer :: ungriddedUBound(1) + real(r8), pointer :: dataptr(:) + real(r8), pointer :: dataptr1d(:) + real(r8), pointer :: dataptr2d(:,:) + real(r8), pointer :: dataptr2d_packed(:,:) + type(ESMF_Field) :: field_fracsrc + type(ESMF_Field), pointer :: fieldlist_src(:) + type(ESMF_Field), pointer :: fieldlist_dst(:) + type(ESMF_Field), pointer :: fieldlist_dat(:) + real(r8), pointer :: data_norm(:) + real(r8), pointer :: data_dst(:,:) + character(cl) :: field_name + character(cl), allocatable :: field_namelist_dat(:) + logical :: skip_mapping + type(ESMF_Region_Flag) :: zeroregion + real(ESMF_KIND_R8), parameter :: fillValue = 9.99e20_ESMF_KIND_R8 + character(len=*), parameter :: subname=' (med_map_mod:med_map_field_packed) ' !----------------------------------------------------------- call t_startf('MED:'//subname) @@ -977,6 +992,22 @@ subroutine med_map_field_packed(FBSrc, FBDst, FBFracSrc, field_normOne, packed_d fieldcount=0 endif + ! Get field count for FBDat if it is given and created + fieldcount_dat = 0 + skip_mapping = .false. + if (present(FBdat)) then + if (ESMF_FieldBundleIsCreated(FBdat)) then + call ESMF_FieldBundleGet(FBDat, fieldCount=fieldcount_dat, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + allocate(fieldlist_dat(fieldcount_dat)) + allocate(field_namelist_dat(fieldcount_dat)) + call ESMF_FieldBundleGet(FBDat, fieldlist=fieldlist_dat, fieldNameList=field_namelist_dat, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + if (present(use_data)) skip_mapping = use_data + end if + end if ! Loop over mapping types do mapindex = 1,nmappers @@ -1012,6 +1043,7 @@ subroutine med_map_field_packed(FBSrc, FBDst, FBFracSrc, field_normOne, packed_d ! Get the indices into the packed data structure np = packed_data(mapindex)%fldindex(nf) if (np > 0) then + ! Fill packed source field call ESMF_FieldGet(fieldlist_src(nf), ungriddedUBound=ungriddedUBound, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (ungriddedUBound(1) > 0) then @@ -1027,8 +1059,87 @@ subroutine med_map_field_packed(FBSrc, FBDst, FBFracSrc, field_normOne, packed_d end if end if end do + + ! Nullify pointers + nullify(dataptr2d_packed) + nullify(dataptr2d) + nullify(dataptr1d) + call t_stopf('MED:'//trim(subname)//' copy from src') + ! ----------------------------------- + ! Fill destination field with background data provided by CDEPS inline + ! ----------------------------------- + + if (fieldcount_dat > 0) then + ! First get the pointer for the packed destination data + call ESMF_FieldGet(packed_data(mapindex)%field_dst, farrayptr=dataptr2d_packed, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Loop over fields and fill it if there is a match + do nf = 1,fieldcount + ! Get the indices into the packed data structure + np = packed_data(mapindex)%fldindex(nf) + if (np > 0) then + ! Get size of ungridded dimension and name of the field + call ESMF_FieldGet(fieldlist_dst(nf), ungriddedUBound=ungriddedUBound, name=field_name, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + if (maintask) write(logunit,'(a)') trim(subname)//" search "//trim(field_name)//" field for background fill." + + ! Check if field has match in data fields + do nfd = 1, fieldcount_dat + ! Debug output for checked fields to find match + if (maintask .and. dbug_flag > 1) write(logunit,'(a)') trim(field_name)//" - "//trim(field_namelist_dat(nfd)) + + if (trim(field_name) == trim(field_namelist_dat(nfd))) then + ! Debug output about match + if (maintask) write(logunit,'(a)') trim(subname)//" field "//trim(field_namelist_dat(nfd))//" is found!" + + ! Get pointer from data field + call ESMF_FieldGet(fieldlist_dat(nfd), farrayptr=dataptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + if (dbug_flag > 1) then + call Field_diagnose(packed_data(mapindex)%field_dst, trim(field_name), " --> before background fill: ", rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + end if + + ! Get pointer from destination field and fill it with data + if (ungriddedUBound(1) > 0) then + ! TODO: Currently assumes same data along the ungridded dimension + do nu = 1,ungriddedUBound(1) + dataptr2d_packed(np+nu-1,:) = dataptr(:) + end do + else + dataptr2d_packed(np,:) = dataptr(:) + end if + + if (dbug_flag > 1) then + call Field_diagnose(packed_data(mapindex)%field_dst, trim(field_name), " --> after background fill: ", rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + end if + + ! Exit from loop since match is already found + exit + end if + end do + end if + end do + + ! Set zeroregion option to select since we are blending data + zeroregion = ESMF_REGION_SELECT + else + ! Fill packed destination field/s with large value if data is unavailable + ! The data needs to be merged in the component side + ! This is also required for mapfillv_bilnr interpolation type + call ESMF_FieldFill(packed_data(mapindex)%field_dst, dataFillScheme="const", const1=fillValue, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Set zeroregion option to total since we have no data to blend + zeroregion = ESMF_REGION_TOTAL + end if + ! ----------------------------------- ! Do the mapping ! ----------------------------------- @@ -1049,6 +1160,7 @@ subroutine med_map_field_packed(FBSrc, FBDst, FBFracSrc, field_normOne, packed_d trim(packed_data(mapindex)%mapnorm) /= 'none') then ! Normalized mapping - assume that each packed field has only one normalization type + call ESMF_LogWrite(trim(subname)//": FB get "//trim(packed_data(mapindex)%mapnorm), ESMF_LOGMSG_INFO) call ESMF_FieldBundleGet(FBFracSrc, packed_data(mapindex)%mapnorm, field=field_fracsrc, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call med_map_field_normalized(& @@ -1062,30 +1174,36 @@ subroutine med_map_field_packed(FBSrc, FBDst, FBFracSrc, field_normOne, packed_d else if ( trim(packed_data(mapindex)%mapnorm) == 'one' .or. trim(packed_data(mapindex)%mapnorm) == 'none') then - ! Mapping with no normalization that is not redistribution - call med_map_field (& - field_src=packed_data(mapindex)%field_src, & - field_dst=packed_data(mapindex)%field_dst, & - routehandles=routehandles, & - maptype=mapindex, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - - ! Obtain unity normalization factor and multiply - ! interpolated field by reciprocal of normalization factor - if (trim(packed_data(mapindex)%mapnorm) == 'one') then - call ESMF_FieldGet(field_normOne(mapindex), farrayPtr=data_norm, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - call ESMF_FieldGet(packed_data(mapindex)%field_dst, farrayPtr=data_dst, rc=rc) + ! Skip mapping if it is requested + if (skip_mapping) then + if (maintask) write(logunit,'(a)') trim(subname)//" skip mapping since use_data is set to .true." + else + ! Mapping with no normalization that is not redistribution + call med_map_field (& + field_src=packed_data(mapindex)%field_src, & + field_dst=packed_data(mapindex)%field_dst, & + routehandles=routehandles, & + maptype=mapindex, & + zeroregiontype=zeroregion, & + rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - do n = 1,size(data_dst,dim=2) - if (data_norm(n) == 0.0_r8) then - data_dst(:,n) = 0.0_r8 - else - data_dst(:,n) = data_dst(:,n)/data_norm(n) - end if - end do - end if + ! Obtain unity normalization factor and multiply + ! interpolated field by reciprocal of normalization factor + if (trim(packed_data(mapindex)%mapnorm) == 'one') then + call ESMF_FieldGet(field_normOne(mapindex), farrayPtr=data_norm, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldGet(packed_data(mapindex)%field_dst, farrayPtr=data_dst, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + do n = 1,size(data_dst,dim=2) + if (data_norm(n) == 0.0_r8) then + data_dst(:,n) = 0.0_r8 + else + data_dst(:,n) = data_dst(:,n)/data_norm(n) + end if + end do + end if + end if end if call t_stopf('MED:'//trim(subname)//' map') @@ -1126,8 +1244,12 @@ subroutine med_map_field_packed(FBSrc, FBDst, FBFracSrc, field_normOne, packed_d end do ! end of loop over mapindex if (ESMF_FieldBundleIsCreated(FBsrc)) then - deallocate(fieldlist_src) - deallocate(fieldlist_dst) + deallocate(fieldlist_src) + deallocate(fieldlist_dst) + end if + if (fieldcount_dat > 0) then + deallocate(fieldlist_dat) + deallocate(field_namelist_dat) end if call t_stopf('MED:'//subname) @@ -1249,7 +1371,7 @@ subroutine med_map_field_normalized(field_src, field_dst, routehandles, maptype, end subroutine med_map_field_normalized !================================================================================ - subroutine med_map_field(field_src, field_dst, routehandles, maptype, fldname, rc) + subroutine med_map_field(field_src, field_dst, routehandles, maptype, fldname, zeroregiontype, rc) !--------------------------------------------------- ! map the source field to the destination field @@ -1257,29 +1379,29 @@ subroutine med_map_field(field_src, field_dst, routehandles, maptype, fldname, r use ESMF , only : ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_SUCCESS use ESMF , only : ESMF_LOGMSG_ERROR, ESMF_FAILURE, ESMF_MAXSTR - use ESMF , only : ESMF_KIND_R8 use ESMF , only : ESMF_Field, ESMF_FieldRegrid - use ESMF , only : ESMF_FieldFill - use ESMF , only : ESMF_TERMORDER_SRCSEQ, ESMF_Region_Flag, ESMF_REGION_TOTAL - use ESMF , only : ESMF_REGION_SELECT + use ESMF , only : ESMF_TERMORDER_SRCSEQ, ESMF_Region_Flag + use ESMF , only : ESMF_REGION_TOTAL, ESMF_REGION_SELECT use ESMF , only : ESMF_RouteHandle + use ESMF , only : ESMF_FieldWriteVTK use med_internalstate_mod , only : mapnstod_consd, mapnstod_consf, mapnstod_consd, mapnstod use med_internalstate_mod , only : mapconsd, mapconsf use med_internalstate_mod , only : mapfillv_bilnr use med_methods_mod , only : Field_diagnose => med_methods_Field_diagnose ! input/output variables - type(ESMF_Field) , intent(in) :: field_src - type(ESMF_Field) , intent(inout) :: field_dst - type(ESMF_RouteHandle) , intent(inout) :: routehandles(:) - integer , intent(in) :: maptype - character(len=*) , intent(in), optional :: fldname - integer , intent(out) :: rc + type(ESMF_Field) , intent(in) :: field_src + type(ESMF_Field) , intent(inout) :: field_dst + type(ESMF_RouteHandle) , intent(inout) :: routehandles(:) + integer , intent(in) :: maptype + character(len=*), optional , intent(in) :: fldname + type(ESMF_Region_Flag), optional, intent(in) :: zeroregiontype + integer, optional , intent(out) :: rc ! local variables logical :: checkflag = .false. character(len=CS) :: lfldname - real(ESMF_KIND_R8), parameter :: fillValue = 9.99e20_ESMF_KIND_R8 + type(ESMF_Region_Flag) :: zeroregion character(len=*), parameter :: subname='(med_map_mod:med_map_field) ' !--------------------------------------------------- @@ -1291,9 +1413,12 @@ subroutine med_map_field(field_src, field_dst, routehandles, maptype, fldname, r lfldname = 'unknown' if (present(fldname)) lfldname = trim(fldname) + zeroregion = ESMF_REGION_TOTAL + if (present(zeroregiontype)) zeroregion = zeroregiontype + if (maptype == mapnstod_consd) then call ESMF_FieldRegrid(field_src, field_dst, routehandle=RouteHandles(mapnstod), & - termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_TOTAL, rc=rc) + termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=zeroregion, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (dbug_flag > 1) then call Field_diagnose(field_dst, lfldname, " --> after nstod: ", rc=rc) @@ -1308,7 +1433,7 @@ subroutine med_map_field(field_src, field_dst, routehandles, maptype, fldname, r end if else if (maptype == mapnstod_consf) then call ESMF_FieldRegrid(field_src, field_dst, routehandle=RouteHandles(mapnstod), & - termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_TOTAL, rc=rc) + termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=zeroregion, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (dbug_flag > 1) then call Field_diagnose(field_dst, lfldname, " --> after nstod: ", rc=rc) @@ -1322,12 +1447,6 @@ subroutine med_map_field(field_src, field_dst, routehandles, maptype, fldname, r if (chkerr(rc,__LINE__,u_FILE_u)) return end if else if (maptype == mapfillv_bilnr) then - call ESMF_FieldFill(field_dst, dataFillScheme="const", const1=fillValue, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - if (dbug_flag > 1) then - call Field_diagnose(field_dst, lfldname, " --> after fillv: ", rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - end if call ESMF_FieldRegrid(field_src, field_dst, routehandle=RouteHandles(mapfillv_bilnr), & termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_SELECT, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return diff --git a/mediator/med_methods_mod.F90 b/mediator/med_methods_mod.F90 index 649c9c511..d4bdab2a7 100644 --- a/mediator/med_methods_mod.F90 +++ b/mediator/med_methods_mod.F90 @@ -44,6 +44,7 @@ module med_methods_mod public med_methods_FB_init_pointer public med_methods_FB_reset public med_methods_FB_diagnose + public med_methods_FB_write public med_methods_FB_FldChk public med_methods_FB_GetFldPtr public med_methods_FB_getNameN @@ -999,6 +1000,72 @@ subroutine med_methods_FB_diagnose(FB, string, rc) end subroutine med_methods_FB_diagnose !----------------------------------------------------------------------------- + + subroutine med_methods_FB_write(FB, string, rc) + ! ---------------------------------------------- + ! Diagnose status of FB + ! ---------------------------------------------- + + use ESMF, only : ESMF_FieldBundle, ESMF_FieldBundleGet + use ESMF, only : ESMF_Field, ESMF_FieldGet + use ESMF, only : ESMF_FieldWriteVTK + + type(ESMF_FieldBundle) , intent(inout) :: FB + character(len=*) , intent(in), optional :: string + integer , intent(out) :: rc + + ! local variables + integer :: n + integer :: fieldCount, lrank + character(ESMF_MAXSTR), pointer :: lfieldnamelist(:) + character(len=CL) :: lstring + type(ESMF_Field) :: lfield + character(len=*), parameter :: subname='(med_methods_FB_write)' + ! ---------------------------------------------- + + call ESMF_LogWrite(trim(subname)//": called", ESMF_LOGMSG_INFO) + rc = ESMF_SUCCESS + + lstring = '' + if (present(string)) then + lstring = trim(string) + endif + + ! Determine number of fields in field bundle and allocate memory for lfieldnamelist + call ESMF_FieldBundleGet(FB, fieldCount=fieldCount, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + allocate(lfieldnamelist(fieldCount)) + + ! Get the fields in the field bundle + call ESMF_FieldBundleGet(FB, fieldNameList=lfieldnamelist, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! For each field in the bundle, get its memory location and print out the field + do n = 1, fieldCount + call ESMF_FieldBundleGet(FB, fieldName=trim(lfieldnamelist(n)), field=lfield, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + call ESMF_FieldGet(lfield, rank=lrank, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + if (lrank == 1) then + call ESMF_FieldWriteVTK(lfield, trim(lfieldnamelist(n))//'_'//trim(lstring), rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + else + call ESMF_LogWrite(trim(subname)//": ERROR rank not supported ", ESMF_LOGMSG_ERROR) + rc = ESMF_FAILURE + return + endif + end do + + ! Deallocate memory + deallocate(lfieldnamelist) + + call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO) + + end subroutine med_methods_FB_write + + !----------------------------------------------------------------------------- #ifdef DIAGNOSE subroutine med_methods_Array_diagnose(array, string, rc) diff --git a/mediator/med_phases_aofluxes_mod.F90 b/mediator/med_phases_aofluxes_mod.F90 index 48055e92e..5252e6edc 100644 --- a/mediator/med_phases_aofluxes_mod.F90 +++ b/mediator/med_phases_aofluxes_mod.F90 @@ -78,6 +78,7 @@ module med_phases_aofluxes_mod logical :: compute_atm_dens logical :: compute_atm_thbot integer :: ocn_surface_flux_scheme ! use case + logical :: add_gusts character(len=CS), pointer :: fldnames_ocn_in(:) character(len=CS), pointer :: fldnames_atm_in(:) @@ -125,6 +126,7 @@ module med_phases_aofluxes_mod real(R8) , pointer :: shum_HDO (:) => null() ! atm HDO tracer real(R8) , pointer :: shum_18O (:) => null() ! atm H218O tracer real(R8) , pointer :: lwdn (:) => null() ! atm downward longwave heat flux + real(R8) , pointer :: rainc (:) => null() ! convective rain flux ! local size and computational mask and area: on aoflux grid integer :: lsize ! local size integer , pointer :: mask (:) => null() ! integer ocn domain mask: 0 <=> inactive cell @@ -146,6 +148,7 @@ module med_phases_aofluxes_mod real(R8) , pointer :: qref (:) => null() ! diagnostic: 2m ref Q real(R8) , pointer :: u10 (:) => null() ! diagnostic: 10m wind speed real(R8) , pointer :: duu10n (:) => null() ! diagnostic: 10m wind speed squared + real(R8) , pointer :: ugust_out (:) => null() ! diagnostic: gust wind added real(R8) , pointer :: ustar (:) => null() ! saved ustar real(R8) , pointer :: re (:) => null() ! saved re real(R8) , pointer :: ssq (:) => null() ! saved sq @@ -402,6 +405,14 @@ subroutine med_aofluxes_init(gcomp, aoflux_in, aoflux_out, rc) end if #endif + call NUOPC_CompAttributeGet(gcomp, name='add_gusts', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + read(cvalue,*) add_gusts + else + add_gusts = .false. + end if + ! bottom level potential temperature and/or botom level density ! will need to be computed if not received from the atm if (FB_fldchk(is_local%Wrap%FBImp(Compatm,Compatm), 'Sa_ptem', rc=rc)) then @@ -1024,7 +1035,7 @@ subroutine med_aofluxes_update(gcomp, aoflux_in, aoflux_out, rc) end if if (compute_atm_dens) then if (trim(aoflux_code) == 'ccpp' .and. & - (trim(coupling_mode) == 'nems_frac_aoflux' .or. trim(coupling_mode) == 'nems_frac_aoflux_sbs')) then + (trim(coupling_mode) == 'ufs.frac.aoflux')) then ! Add limiting factor to humidity to be consistent with UFS aoflux calculation do n = 1,aoflux_in%lsize if (aoflux_in%mask(n) /= 0.0_r8) then @@ -1052,6 +1063,7 @@ subroutine med_aofluxes_update(gcomp, aoflux_in, aoflux_out, rc) call flux_atmocn (logunit=logunit, & nMax=aoflux_in%lsize, & zbot=aoflux_in%zbot, ubot=aoflux_in%ubot, vbot=aoflux_in%vbot, thbot=aoflux_in%thbot, qbot=aoflux_in%shum, & + rainc=aoflux_in%rainc, & s16O=aoflux_in%shum_16O, sHDO=aoflux_in%shum_HDO, s18O=aoflux_in%shum_18O, rbot=aoflux_in%dens, & tbot=aoflux_in%tbot, us=aoflux_in%uocn, vs=aoflux_in%vocn, pslv=aoflux_in%psfc, ts=aoflux_in%tocn, & mask=aoflux_in%mask, seq_flux_atmocn_minwind=0.5_r8, & @@ -1060,7 +1072,10 @@ subroutine med_aofluxes_update(gcomp, aoflux_in, aoflux_out, rc) evap=aoflux_out%evap, evap_16O=aoflux_out%evap_16O, evap_HDO=aoflux_out%evap_HDO, evap_18O=aoflux_out%evap_18O, & taux=aoflux_out%taux, tauy=aoflux_out%tauy, tref=aoflux_out%tref, qref=aoflux_out%qref, & ocn_surface_flux_scheme=ocn_surface_flux_scheme, & - duu10n=aoflux_out%duu10n, ustar_sv=aoflux_out%ustar, re_sv=aoflux_out%re, ssq_sv=aoflux_out%ssq, & + add_gusts=add_gusts, & + duu10n=aoflux_out%duu10n, & + ugust_out = aoflux_out%ugust_out, & + ustar_sv=aoflux_out%ustar, re_sv=aoflux_out%re, ssq_sv=aoflux_out%ssq, & missval=0.0_r8) #else @@ -1084,7 +1099,8 @@ subroutine med_aofluxes_update(gcomp, aoflux_in, aoflux_out, rc) ocn_surface_flux_scheme=ocn_surface_flux_scheme, & sen=aoflux_out%sen, lat=aoflux_out%lat, lwup=aoflux_out%lwup, evap=aoflux_out%evap, & taux=aoflux_out%taux, tauy=aoflux_out%tauy, tref=aoflux_out%tref, qref=aoflux_out%qref, & - duu10n=aoflux_out%duu10n, missval=0.0_r8) + duu10n=aoflux_out%duu10n, & + missval=0.0_r8) #ifdef UFS_AOFLUX end if #endif @@ -1562,8 +1578,8 @@ subroutine set_aoflux_in_pointers(fldbun_a, fldbun_o, aoflux_in, lsize, xgrid, r lsize = size(aoflux_in%zbot) aoflux_in%lsize = lsize - ! bulk formula quantities for nems_orig_data - if (trim(coupling_mode) == 'nems_orig_data' .and. ocn_surface_flux_scheme == -1) then + ! bulk formula quantities for ufs non-frac with med-aoflux + if (trim(coupling_mode) == 'ufs.nfrac.aoflux' .and. ocn_surface_flux_scheme == -1) then call fldbun_getfldptr(fldbun_a, 'Sa_u10m', aoflux_in%ubot, xgrid=xgrid, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call fldbun_getfldptr(fldbun_a, 'Sa_v10m', aoflux_in%vbot, xgrid=xgrid, rc=rc) @@ -1581,10 +1597,14 @@ subroutine set_aoflux_in_pointers(fldbun_a, fldbun_o, aoflux_in, lsize, xgrid, r if (chkerr(rc,__LINE__,u_FILE_u)) return call fldbun_getfldptr(fldbun_a, 'Sa_shum', aoflux_in%shum, xgrid=xgrid, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return + if (add_gusts) then + call fldbun_getfldptr(fldbun_a, 'Faxa_rainc', aoflux_in%rainc, xgrid=xgrid, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + end if end if - ! extra fields for nems_frac_aoflux - if (trim(coupling_mode) == 'nems_frac_aoflux' .or. trim(coupling_mode) == 'nems_frac_aoflux_sbs') then + ! extra fields for ufs.frac.aoflux + if (trim(coupling_mode) == 'ufs.frac.aoflux') then call fldbun_getfldptr(fldbun_a, 'Sa_u10m', aoflux_in%usfc, xgrid=xgrid, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call fldbun_getfldptr(fldbun_a, 'Sa_v10m', aoflux_in%vsfc, xgrid=xgrid, rc=rc) @@ -1618,7 +1638,7 @@ subroutine set_aoflux_in_pointers(fldbun_a, fldbun_o, aoflux_in, lsize, xgrid, r if (compute_atm_dens .or. compute_atm_thbot) then call fldbun_getfldptr(fldbun_a, 'Sa_pbot', aoflux_in%pbot, xgrid=xgrid, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - if (trim(coupling_mode) == 'nems_frac_aoflux' .or. trim(coupling_mode) == 'nems_frac_aoflux_sbs') then + if (trim(coupling_mode) == 'ufs.frac.aoflux') then call fldbun_getfldptr(fldbun_a, 'Sa_pslv', aoflux_in%psfc, xgrid=xgrid, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return end if @@ -1704,6 +1724,7 @@ subroutine set_aoflux_out_pointers(fldbun, lsize, aoflux_out, xgrid, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call fldbun_getfldptr(fldbun, 'Faox_lwup', aoflux_out%lwup, xgrid=xgrid, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return + if (flds_wiso) then call fldbun_getfldptr(fldbun, 'Faox_evap_16O', aoflux_out%evap_16O, xgrid=xgrid, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return @@ -1717,6 +1738,13 @@ subroutine set_aoflux_out_pointers(fldbun, lsize, aoflux_out, xgrid, rc) allocate(aoflux_out%evap_HDO(lsize)); aoflux_out%evap_HDO(:) = 0._R8 end if + if (add_gusts) then + call fldbun_getfldptr(fldbun, 'So_ugustOut', aoflux_out%ugust_out, xgrid=xgrid, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + else + allocate(aoflux_out%ugust_out(lsize)); aoflux_out%ugust_out(:) = 0._R8 + end if + end subroutine set_aoflux_out_pointers !================================================================================ diff --git a/mediator/med_phases_cdeps_mod.F90 b/mediator/med_phases_cdeps_mod.F90 new file mode 100644 index 000000000..72ac560cc --- /dev/null +++ b/mediator/med_phases_cdeps_mod.F90 @@ -0,0 +1,292 @@ +module med_phases_cdeps_mod + + use ESMF, only: ESMF_Clock, ESMF_ClockGet, ESMF_Time, ESMF_TimeGet + use ESMF, only: ESMF_Mesh + use ESMF, only: ESMF_GridComp, ESMF_GridCompGet + use ESMF, only: ESMF_LogWrite + use ESMF, only: ESMF_Field, ESMF_FieldGet + use ESMF, only: ESMF_FieldBundleGet, ESMF_FieldBundleIsCreated + use ESMF, only: ESMF_FieldBundleCreate + use ESMF, only: ESMF_GridCompGetInternalState + use ESMF, only: ESMF_SUCCESS, ESMF_LOGMSG_INFO + + use med_internalstate_mod, only: InternalState + use med_internalstate_mod, only: logunit, maintask + use med_internalstate_mod, only: ncomps, compname, compatm, compocn + use perf_mod , only: t_startf, t_stopf + use med_kind_mod , only: cl => shr_kind_cl + use med_kind_mod , only: r8 => shr_kind_r8 + use med_constants_mod , only: dbug_flag => med_constants_dbug_flag + use med_utils_mod , only: chkerr => med_utils_ChkErr + use med_methods_mod , only: FB_FldChk => med_methods_FB_FldChk + use med_methods_mod , only: FB_getFieldN => med_methods_FB_getFieldN + use med_methods_mod , only: FB_getNumflds => med_methods_FB_getNumflds + use med_methods_mod , only: FB_init => med_methods_FB_Init + use med_methods_mod , only: FB_diagnose => med_methods_FB_diagnose + use med_methods_mod , only: FB_write => med_methods_FB_write + use med_methods_mod , only: FB_GetFldPtr => med_methods_FB_GetFldPtr + + use dshr_mod , only: dshr_pio_init + use dshr_strdata_mod , only: shr_strdata_type + use dshr_strdata_mod , only: shr_strdata_init_from_inline + use dshr_strdata_mod , only: shr_strdata_advance + use dshr_stream_mod , only: shr_stream_init_from_esmfconfig + + implicit none + private + + !-------------------------------------------------------------------------- + ! Public interfaces + !-------------------------------------------------------------------------- + + public med_phases_cdeps_run + + !-------------------------------------------------------------------------- + ! Private data + !-------------------------------------------------------------------------- + + type config + integer :: year_first + integer :: year_last + integer :: year_align + integer :: offset + real(r8) :: dtlimit + character(len=cl) :: mesh_filename + character(len=cl), allocatable :: data_filename(:) + character(len=cl), allocatable :: fld_list(:) + character(len=cl), allocatable :: fld_list_model(:) + character(len=cl) :: mapalgo + character(len=cl) :: taxmode + character(len=cl) :: tintalgo + character(len=cl) :: name + end type config + + type(config) :: stream ! stream configuration + type(shr_strdata_type), allocatable :: sdat(:,:) ! input data stream + + character(*),parameter :: u_FILE_u = __FILE__ + +!============================================================================ +contains +!============================================================================ + + subroutine med_phases_cdeps_run(gcomp, rc) + + !------------------------------------------------------------------------ + ! Use CDEPS inline capability to read in data + !------------------------------------------------------------------------ + + use ESMF, only : ESMF_GridComp + + ! input/output variables + type(ESMF_GridComp) :: gcomp + integer, intent(out) :: rc + + ! local variables + type(InternalState) :: is_local + type(ESMF_Clock) :: clock + type(ESMF_Time) :: currTime + type(ESMF_Mesh) :: meshdst + type(ESMF_Field) :: flddst + integer :: i, j, k, l, nflds, streamid + integer :: n1, n2, item, nstreams, localPet + integer :: curr_ymd, sec + integer :: year, month, day, hour, minute, second + logical :: found + logical, save :: first_time = .true. + character(len=cl), allocatable :: fileList(:), varList(:,:) + character(len=cl) :: streamfilename, suffix, fldname + type(shr_strdata_type) :: sdat_config + character(len=*), parameter :: subname = '(med_phases_cdeps_run)' + !--------------------------------------- + + rc = ESMF_SUCCESS + + call t_startf('MED:'//subname) + call ESMF_LogWrite(trim(subname)//": called", ESMF_LOGMSG_INFO) + + ! Get the internal state from gcomp + nullify(is_local%wrap) + call ESMF_GridCompGetInternalState(gcomp, is_local, rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Query component + call ESMF_GridCompGet(gcomp, clock=clock, localPet=localPet, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Initialize sdat streams + if (.not. allocated(sdat)) allocate(sdat(ncomps,ncomps)) + sdat(:,:)%mainproc = (localPet == 0) + + ! Initialize cdeps inline + if (first_time) then + ! Init PIO + call dshr_pio_init(gcomp, sdat_config, logunit, rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Read stream configuration file + ! TODO: At this point it only suports ESMF config format (XML?) + streamfilename = 'stream.config' + call shr_stream_init_from_esmfconfig(streamfilename, sdat_config%stream, logunit, & + sdat_config%pio_subsystem, sdat_config%io_type, sdat_config%io_format, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Get number of streams + nstreams = size(sdat_config%stream) + + ! Loop over coupling directions and try to find field match in given streams + do n1 = 1, ncomps + do n2 = 1, ncomps + ! Check for coupling direction and background fill + if (n1 /= n2 .and. is_local%wrap%med_coupling_active(n1,n2) .and. is_local%wrap%med_data_active(n1,n2)) then + ! Get number of fields + call FB_getNumflds(is_local%wrap%FBImp(n1,n2), trim(subname), nflds, rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Loop over fields and try to find it in the given stream + found = .false. + do i = 1, nflds + ! Query destination field + call FB_getFieldN(is_local%wrap%FBImp(n1,n2), i, flddst, rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Query destination field name and its mesh + call ESMF_FieldGet(flddst, mesh=meshdst, name=fldname, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + if (maintask) write(logunit,'(a)') trim(subname)//": extracting destination mesh from "//trim(fldname) + + ! Check if any field in FB in the given stream + ! NOTE: Single stream could provide multiple fields !!! + streamid = 0 + do j = 1, nstreams + do k = 1, sdat_config%stream(j)%nvars + if (trim(sdat_config%stream(j)%varlist(k)%nameinmodel) == trim(fldname)) then + streamid = j + end if + end do + end do + + ! If match is found and previously not initialized, then initialize cdeps inline for the stream + if (size(sdat(n1,n2)%stream) == 0 .and. streamid /= 0) then + ! Debug print + if (maintask) then + write(logunit,'(a,i3)') trim(subname)//": initialize stream ", streamid + end if + + ! Allocate temporary variable to store file names in the stream + allocate(fileList(sdat_config%stream(streamid)%nfiles)) + allocate(varList(sdat_config%stream(streamid)%nvars,2)) + + ! Fill file abd variable lists with data + do l = 1, sdat_config%stream(streamid)%nfiles + fileList(l) = trim(sdat_config%stream(streamid)%file(l)%name) + if (maintask) write(logunit,'(a,i2,x,a)') trim(subname)//": file ", l, trim(fileList(l)) + end do + do l = 1, sdat_config%stream(streamid)%nvars + varList(l,1) = trim(sdat_config%stream(streamid)%varlist(l)%nameinfile) + varList(l,2) = trim(sdat_config%stream(streamid)%varlist(l)%nameinmodel) + if (maintask) write(logunit,'(a,i2,x,a)') trim(subname)//": variable ", l, trim(varList(l,1))//" -> "//trim(varList(l,2)) + end do + + ! Set PIO related variables + sdat(n1,n2)%pio_subsystem => sdat_config%pio_subsystem + sdat(n1,n2)%io_type = sdat_config%io_type + sdat(n1,n2)%io_format = sdat_config%io_format + + ! Init stream + call shr_strdata_init_from_inline(sdat(n1,n2), my_task=localPet, logunit=logunit, & + compname = 'cmeps', model_clock=clock, model_mesh=meshdst, & + stream_meshfile=trim(sdat_config%stream(streamid)%meshfile), & + stream_filenames=fileList, & + stream_yearFirst=sdat_config%stream(streamid)%yearFirst, & + stream_yearLast=sdat_config%stream(streamid)%yearLast, & + stream_yearAlign=sdat_config%stream(streamid)%yearAlign, & + stream_fldlistFile=varList(:,1), & + stream_fldListModel=varList(:,2), & + stream_lev_dimname=trim(sdat_config%stream(streamid)%lev_dimname), & + stream_mapalgo=trim(sdat_config%stream(streamid)%mapalgo), & + stream_offset=sdat_config%stream(streamid)%offset, & + stream_taxmode=trim(sdat_config%stream(streamid)%taxmode), & + stream_dtlimit=sdat_config%stream(streamid)%dtlimit, & + stream_tintalgo=trim(sdat_config%stream(streamid)%tInterpAlgo), & + stream_name=trim(compname(n1))//'_'//trim(compname(n2)), & + stream_src_mask=sdat_config%stream(streamid)%src_mask_val, & + stream_dst_mask=sdat_config%stream(streamid)%dst_mask_val, & + rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Print out source and destination mask used in the stream + if (maintask) write(logunit,'(a,2i2)') trim(subname)//": mask values src, dst ", & + sdat_config%stream(streamid)%src_mask_val, sdat_config%stream(streamid)%dst_mask_val + + ! Remove temporary variables + deallocate(fileList) + deallocate(varList) + + ! Set flag + found = .true. + end if + end do ! nflds + + ! Create empty FB + if (.not. ESMF_FieldBundleIsCreated(is_local%wrap%FBData(n2), rc=rc) .and. found) then + is_local%wrap%FBData(n2) = ESMF_FieldBundleCreate(name="inline_"//trim(compname(n2)), rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + end if + end if + end do ! n2 + end do ! n1 + + ! Set flag to false + first_time = .false. + end if + + ! Get current time + call ESMF_ClockGet(clock, currTime=currTime, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Query current time + call ESMF_TimeGet(currTime, yy=year, mm=month, dd=day, h=hour, m=minute, s=second, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + curr_ymd = abs(year)*10000+month*100+day + sec = hour*3600+minute*60+second + + ! Read data if stream initialized + do n1 = 1, ncomps + do n2 = 1, ncomps + if (size(sdat(n1,n2)%stream) > 0) then + ! Debug print + if (maintask) then + write(logunit,'(a)') trim(subname)//": read stream "//trim(compname(n1))//" -> "//trim(compname(n2)) + end if + + ! Read data + call shr_strdata_advance(sdat(n1,n2), ymd=curr_ymd, tod=sec, logunit=logunit, & + istr=trim(compname(n1))//'_'//trim(compname(n2)), rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Check FB + call FB_diagnose(sdat(n1,n2)%pstrm(1)%fldbun_model, & + trim(subname)//':'//trim(compname(n1))//'_'//trim(compname(n2)), rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Point FB from internal one + is_local%wrap%FBData(n2) = sdat(n1,n2)%pstrm(1)%fldbun_model + + ! Write FB for debugging + if (dbug_flag > 10) then + write(suffix, fmt='(i4,a1,i2.2,a1,i2.2,a1,i5.5)') year, '-', month, '-', day, '-', sec + call FB_write(is_local%wrap%FBData(n2), suffix, rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + end if + end if + end do + end do + + call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO) + call t_stopf('MED:'//subname) + + end subroutine med_phases_cdeps_run + +end module med_phases_cdeps_mod diff --git a/mediator/med_phases_post_atm_mod.F90 b/mediator/med_phases_post_atm_mod.F90 index 9ed1b78d4..333497a69 100644 --- a/mediator/med_phases_post_atm_mod.F90 +++ b/mediator/med_phases_post_atm_mod.F90 @@ -65,6 +65,7 @@ subroutine med_phases_post_atm(gcomp, rc) FBSrc=is_local%wrap%FBImp(compatm,compatm), & FBDst=is_local%wrap%FBImp(compatm,compocn), & FBFracSrc=is_local%wrap%FBFrac(compatm), & + FBDat=is_local%wrap%FBData(compocn), & field_normOne=is_local%wrap%field_normOne(compatm,compocn,:), & packed_data=is_local%wrap%packed_data(compatm,compocn,:), & routehandles=is_local%wrap%RH(compatm,compocn,:), rc=rc) @@ -104,6 +105,7 @@ subroutine med_phases_post_atm(gcomp, rc) FBSrc=is_local%wrap%FBImp(compatm,compatm), & FBDst=is_local%wrap%FBImp(compatm,compwav), & FBFracSrc=is_local%wrap%FBFrac(compatm), & + FBDat=is_local%wrap%FBData(compwav), & field_normOne=is_local%wrap%field_normOne(compatm,compwav,:), & packed_data=is_local%wrap%packed_data(compatm,compwav,:), & routehandles=is_local%wrap%RH(compatm,compwav,:), rc=rc) diff --git a/mediator/med_phases_prep_atm_mod.F90 b/mediator/med_phases_prep_atm_mod.F90 index d9c9dfc55..edb74f38b 100644 --- a/mediator/med_phases_prep_atm_mod.F90 +++ b/mediator/med_phases_prep_atm_mod.F90 @@ -79,6 +79,7 @@ subroutine med_phases_prep_atm(gcomp, rc) FBSrc=is_local%wrap%FBImp(compocn,compocn), & FBDst=is_local%wrap%FBImp(compocn,compatm), & FBFracSrc=is_local%wrap%FBFrac(compocn), & + FBDat=is_local%wrap%FBData(compatm), & field_NormOne=is_local%wrap%field_normOne(compocn,compatm,:), & packed_data=is_local%wrap%packed_data(compocn,compatm,:), & routehandles=is_local%wrap%RH(compocn,compatm,:), rc=rc) @@ -113,8 +114,7 @@ subroutine med_phases_prep_atm(gcomp, rc) !--- map atm/ocn fluxes from ocn to atm grid if appropriate !--------------------------------------- if (trim(coupling_mode) == 'cesm' .or. & - trim(coupling_mode) == 'nems_frac_aoflux' .or. & - trim(coupling_mode) == 'nems_frac_aoflux_sbs') then + trim(coupling_mode) == 'ufs.frac.aoflux') then if (is_local%wrap%aoflux_grid == 'ogrid') then call med_aofluxes_map_ogrid2agrid_output(gcomp, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return diff --git a/mediator/med_phases_prep_ocn_mod.F90 b/mediator/med_phases_prep_ocn_mod.F90 index 987fee0c6..456cbf361 100644 --- a/mediator/med_phases_prep_ocn_mod.F90 +++ b/mediator/med_phases_prep_ocn_mod.F90 @@ -96,12 +96,40 @@ subroutine med_phases_prep_ocn_accum(gcomp, rc) rc = ESMF_SUCCESS call memcheck(subname, 5, maintask) - ! Get the internal state + !--------------------------------------- + ! --- Get the internal state + !--------------------------------------- nullify(is_local%wrap) call ESMF_GridCompGetInternalState(gcomp, is_local, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return fldList => med_fldList_GetfldListTo(compocn) - ! auto merges to ocn + + !--------------------------------------- + ! --- map atm to ocn, only if data stream is available + !--------------------------------------- + if (is_local%wrap%med_coupling_active(compatm,compocn) .and. & + is_local%wrap%med_data_active(compatm,compocn) .and. & + is_local%wrap%med_data_force_first(compocn)) then + call t_startf('MED:'//trim(subname)//' map_atm2ocn') + call med_map_field_packed( & + FBSrc=is_local%wrap%FBImp(compatm,compatm), & + FBDst=is_local%wrap%FBImp(compatm,compocn), & + FBFracSrc=is_local%wrap%FBFrac(compocn), & + FBDat=is_local%wrap%FBData(compocn), & + use_data=is_local%wrap%med_data_force_first(compocn), & + field_normOne=is_local%wrap%field_normOne(compatm,compocn,:), & + packed_data=is_local%wrap%packed_data(compatm,compocn,:), & + routehandles=is_local%wrap%RH(compatm,compocn,:), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call t_stopf('MED:'//trim(subname)//' map_atm2ocn') + + ! Reset flag to use data + is_local%wrap%med_data_force_first(compocn) = .false. + end if + + !--------------------------------------- + !--- merge all fields to ocn + !--------------------------------------- call med_merge_auto(& is_local%wrap%med_coupling_active(:,compocn), & is_local%wrap%FBExp(compocn), & @@ -355,11 +383,18 @@ subroutine med_phases_prep_ocn_custom(gcomp, rc) end do ! Compute sw export to ocean bands if required if (export_swnet_by_bands) then - c1 = 0.285; c2 = 0.285; c3 = 0.215; c4 = 0.215 - Foxx_swnet_vdr(:) = c1 * Foxx_swnet(:) - Foxx_swnet_vdf(:) = c2 * Foxx_swnet(:) - Foxx_swnet_idr(:) = c3 * Foxx_swnet(:) - Foxx_swnet_idf(:) = c4 * Foxx_swnet(:) + if (trim(coupling_mode) == 'cesm') then + c1 = 0.285; c2 = 0.285; c3 = 0.215; c4 = 0.215 + Foxx_swnet_vdr(:) = c1 * Foxx_swnet(:) + Foxx_swnet_vdf(:) = c2 * Foxx_swnet(:) + Foxx_swnet_idr(:) = c3 * Foxx_swnet(:) + Foxx_swnet_idf(:) = c4 * Foxx_swnet(:) + else + Foxx_swnet_vdr(:) = Faxa_swvdr(:) * (1.0_R8 - avsdr(:)) + Foxx_swnet_vdf(:) = Faxa_swvdf(:) * (1.0_R8 - avsdf(:)) + Foxx_swnet_idr(:) = Faxa_swndr(:) * (1.0_R8 - anidr(:)) + Foxx_swnet_idf(:) = Faxa_swndf(:) * (1.0_R8 - anidf(:)) + end if end if end if diff --git a/mediator/med_phases_prep_wav_mod.F90 b/mediator/med_phases_prep_wav_mod.F90 index c690aa522..93755d59c 100644 --- a/mediator/med_phases_prep_wav_mod.F90 +++ b/mediator/med_phases_prep_wav_mod.F90 @@ -19,7 +19,7 @@ module med_phases_prep_wav_mod use med_methods_mod , only : FB_reset => med_methods_FB_reset use med_methods_mod , only : FB_check_for_nans => med_methods_FB_check_for_nans use esmFlds , only : med_fldList_GetfldListTo - use med_internalstate_mod , only : compwav + use med_internalstate_mod , only : compatm, compwav use perf_mod , only : t_startf, t_stopf implicit none @@ -92,12 +92,39 @@ subroutine med_phases_prep_wav_accum(gcomp, rc) rc = ESMF_SUCCESS call memcheck(subname, 5, maintask) - ! Get the internal state + !--------------------------------------- + ! --- Get the internal state + !--------------------------------------- nullify(is_local%wrap) call ESMF_GridCompGetInternalState(gcomp, is_local, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! auto merges to wav + !--------------------------------------- + ! --- map atm to wav, only if data stream is available + !--------------------------------------- + if (is_local%wrap%med_coupling_active(compatm,compwav) .and. & + is_local%wrap%med_data_active(compatm,compwav) .and. & + is_local%wrap%med_data_force_first(compwav)) then + call t_startf('MED:'//trim(subname)//' map_atm2wav') + call med_map_field_packed( & + FBSrc=is_local%wrap%FBImp(compatm,compatm), & + FBDst=is_local%wrap%FBImp(compatm,compwav), & + FBFracSrc=is_local%wrap%FBFrac(compatm), & + FBDat=is_local%wrap%FBData(compwav), & + use_data=is_local%wrap%med_data_force_first(compwav), & + field_normOne=is_local%wrap%field_normOne(compatm,compwav,:), & + packed_data=is_local%wrap%packed_data(compatm,compwav,:), & + routehandles=is_local%wrap%RH(compatm,compwav,:), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call t_stopf('MED:'//trim(subname)//' map_atm2wav') + + ! Reset flag to use data + is_local%wrap%med_data_force_first(compwav) = .false. + end if + + !--------------------------------------- + !--- merge all fields to wav + !--------------------------------------- call med_merge_auto(& is_local%wrap%med_coupling_active(:,compwav), & is_local%wrap%FBExp(compwav), & diff --git a/ufs/ccpp/config/ccpp_prebuild_config.py b/ufs/ccpp/config/ccpp_prebuild_config.py index d2872972e..8d8963bad 100755 --- a/ufs/ccpp/config/ccpp_prebuild_config.py +++ b/ufs/ccpp/config/ccpp_prebuild_config.py @@ -25,7 +25,7 @@ VARIABLE_DEFINITION_FILES = [ # actual variable definition files '{}/ccpp/framework/src/ccpp_types.F90'.format(fv3_path), - '{}/ccpp/physics/physics/machine.F'.format(fv3_path), + '{}/ccpp/physics/physics/hooks/machine.F'.format(fv3_path), 'CMEPS/ufs/ccpp/data/MED_typedefs.F90', 'CMEPS/ufs/ccpp/data/MED_data.F90' ] @@ -58,13 +58,13 @@ # Add all physics scheme files relative to basedir SCHEME_FILES = [ - '{}/ccpp/physics/physics/sfc_ocean.F'.format(fv3_path), - '{}/ccpp/physics/physics/sfc_diff.f'.format(fv3_path), - '{}/ccpp/physics/physics/GFS_surface_loop_control_part1.F90'.format(fv3_path), - '{}/ccpp/physics/physics/GFS_surface_loop_control_part2.F90'.format(fv3_path), - '{}/ccpp/physics/physics/GFS_surface_composites_pre.F90'.format(fv3_path), - '{}/ccpp/physics/physics/GFS_surface_composites_post.F90'.format(fv3_path), - '{}/ccpp/physics/physics/sfc_diag.f'.format(fv3_path) + '{}/ccpp/physics/physics/SFC_Models/Ocean/UFS/sfc_ocean.F'.format(fv3_path), + '{}/ccpp/physics/physics/SFC_Layer/UFS/sfc_diff.f'.format(fv3_path), + '{}/ccpp/physics/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_loop_control_part1.F90'.format(fv3_path), + '{}/ccpp/physics/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_loop_control_part2.F90'.format(fv3_path), + '{}/ccpp/physics/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_composites_pre.F90'.format(fv3_path), + '{}/ccpp/physics/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_composites_post.F90'.format(fv3_path), + '{}/ccpp/physics/physics/SFC_Layer/UFS/sfc_diag.f'.format(fv3_path) ] # Default build dir, relative to current working directory, diff --git a/ufs/ccpp/data/MED_typedefs.F90 b/ufs/ccpp/data/MED_typedefs.F90 index e5e1b494f..e7c84506e 100644 --- a/ufs/ccpp/data/MED_typedefs.F90 +++ b/ufs/ccpp/data/MED_typedefs.F90 @@ -73,6 +73,7 @@ module MED_typedefs real (kind=kind_phys),pointer :: lake_q2m (:) => null() !< 2 meter humidity from CLM Lake model real(kind=kind_phys), pointer :: wind(:) => null() !< wind speed at lowest model level (m/s) logical, pointer :: flag_iter(:) => null() !< flag for iteration + logical, pointer :: flag_lakefreeze(:) => null() !< flag for lake freeze real(kind=kind_phys), pointer :: qss_water(:) => null() !< surface air saturation specific humidity over water (kg/kg) real(kind=kind_phys), pointer :: cmm_water(:) => null() !< momentum exchange coefficient over water (m/s) real(kind=kind_phys), pointer :: chh_water(:) => null() !< thermal exchange coefficient over water (kg/m2s) @@ -368,6 +369,8 @@ subroutine interstitial_create(interstitial, im) interstitial%wind = huge allocate(interstitial%flag_iter(im)) interstitial%flag_iter = .true. + allocate(interstitial%flag_lakefreeze(im)) + interstitial%flag_lakefreeze = .false. allocate(interstitial%qss_water(im)) interstitial%qss_water = huge allocate(interstitial%cmm_ice(im)) @@ -571,6 +574,7 @@ subroutine interstitial_phys_reset(interstitial) interstitial%flag_cice = .false. interstitial%flag_guess = .false. interstitial%flag_iter = .true. + interstitial%flag_lakefreeze = .false. interstitial%fm10_ice = huge interstitial%fm10_land = huge interstitial%fm10_water = huge diff --git a/ufs/ccpp/data/MED_typedefs.meta b/ufs/ccpp/data/MED_typedefs.meta index 271110e9c..439a617a3 100644 --- a/ufs/ccpp/data/MED_typedefs.meta +++ b/ufs/ccpp/data/MED_typedefs.meta @@ -237,6 +237,12 @@ units = flag dimensions = (horizontal_loop_extent) type = logical +[flag_lakefreeze] + standard_name = flag_for_lake_water_freeze + long_name = flag for lake water freeze + units = flag + dimensions = (horizontal_loop_extent) + type = logical [qss_water] standard_name = surface_specific_humidity_over_water long_name = surface air saturation specific humidity over water @@ -1312,7 +1318,7 @@ [ccpp-table-properties] name = MED_typedefs type = module - relative_path = ../../../../../FV3/ccpp/physics/physics + relative_path = ../../../../../FV3/ccpp/physics/physics/hooks dependencies = machine.F,physcons.F90 [ccpp-arg-table] diff --git a/ufs/ufs_io_mod.F90 b/ufs/ufs_io_mod.F90 index 8564be8e5..d89a6f014 100644 --- a/ufs/ufs_io_mod.F90 +++ b/ufs/ufs_io_mod.F90 @@ -22,21 +22,14 @@ module ufs_io_mod use ESMF, only : ESMF_FieldWrite, ESMF_FieldBundleRead, ESMF_FieldBundleWrite use ESMF, only : ESMF_REGRIDMETHOD_CONSERVE_2ND, ESMF_MeshCreate use ESMF, only : ESMF_TERMORDER_SRCSEQ, ESMF_REGION_TOTAL + use ESMF, only : ESMF_RouteHandle, ESMF_FieldBundleRedistStore + use ESMF, only : ESMF_FieldBundleRedist, ESMF_RouteHandleDestroy + use ESMF, only : ESMF_TYPEKIND_I4, ESMF_TYPEKIND_R4 use NUOPC, only : NUOPC_CompAttributeGet use NUOPC_Mediator, only : NUOPC_MediatorGet - use fms_mod, only : fms_init - use fms2_io_mod, only : open_file, FmsNetcdfFile_t - use mosaic2_mod, only : get_mosaic_ntiles, get_mosaic_grid_sizes - use mosaic2_mod, only : get_mosaic_contact, get_mosaic_ncontacts - use mpp_mod, only : mpp_pe, mpp_root_pe, mpp_error, FATAL - use mpp_domains_mod, only : mpp_define_layout, mpp_get_compute_domain - use mpp_domains_mod, only : mpp_domains_init, mpp_define_mosaic, domain2d - use mpp_io_mod, only : MPP_RDONLY, MPP_NETCDF, MPP_SINGLE, MPP_MULTI - use mpp_io_mod, only : mpp_get_info, mpp_get_fields, mpp_get_atts - use mpp_io_mod, only : mpp_open, mpp_read, fieldtype - - use med_kind_mod, only : r8=>SHR_KIND_R8, cs=>SHR_KIND_CS, cl=>SHR_KIND_CL + use med_kind_mod, only : r4=>SHR_KIND_R4, r8=>SHR_KIND_R8 + use med_kind_mod, only : cs=>SHR_KIND_CS, cl=>SHR_KIND_CL use med_utils_mod, only : chkerr => med_utils_chkerr use med_constants_mod, only : dbug_flag => med_constants_dbug_flag use med_internalstate_mod, only : InternalState, maintask, logunit @@ -62,24 +55,25 @@ module ufs_io_mod type(ESMF_Grid) :: grid ! ESMF grid object from mosaic file type(ESMF_Mesh) :: mesh ! ESMF mesh object from CS grid type(ESMF_RouteHandle) :: rh ! ESMF routehandle object to redist data from CS grid to mesh - type(domain2d) :: mosaic_domain ! domain object created by FMS integer :: layout(2) ! layout for domain decomposition - integer, allocatable :: nit(:) ! size of tile in i direction - integer, allocatable :: njt(:) ! size of tile in j direction integer :: ntiles ! number of tiles in case of having CS grid - integer :: ncontacts ! number of contacts in case of having CS grid - integer, allocatable :: tile1(:) ! list of tile numbers in tile 1 of each contact - integer, allocatable :: tile2(:) ! list of tile numbers in tile 2 of each contact - integer, allocatable :: istart1(:) ! list of starting i-index in tile 1 of each contact - integer, allocatable :: iend1(:) ! list of ending i-index in tile 1 of each contact - integer, allocatable :: jstart1(:) ! list of starting j-index in tile 1 of each contact - integer, allocatable :: jend1(:) ! list of ending j-index in tile 1 of each contact - integer, allocatable :: istart2(:) ! list of starting i-index in tile 2 of each contact - integer, allocatable :: iend2(:) ! list of ending i-index in tile 2 of each contact - integer, allocatable :: jstart2(:) ! list of starting j-index in tile 2 of each contact - integer, allocatable :: jend2(:) ! list of ending j-index in tile 2 of each contact end type domain_type + type field_type + real(r4), pointer :: ptr1r4(:) ! data pointer for 1d r4 + real(r8), pointer :: ptr1r8(:) ! data pointer for 1d r8 + integer , pointer :: ptr1i4(:) ! data pointer for 1d i4 + real(r4), pointer :: ptr2r4(:,:) ! data pointer for 2d r4 + real(r8), pointer :: ptr2r8(:,:) ! data pointer for 2d r8 + integer , pointer :: ptr2i4(:,:) ! data pointer for 2d i4 + character(len=128) :: short_name = "" ! variable short name + character(len=128) :: units = "" ! variable unit + character(len=128) :: long_name = "" ! variable long name + character(len=128) :: zaxis = "" ! name of z-axis + integer :: nlev ! number of layers in z-axis + integer :: nrec ! number of record in file (time axis) + end type field_type + character(cl) :: case_name = 'unset' ! case name character(*), parameter :: modName = "(ufs_io_mod)" @@ -106,10 +100,11 @@ subroutine read_initial(gcomp, ini_file, mosaic_file, input_dir, layout, rc) type(domain_type) :: domain type(InternalState) :: is_local type(ESMF_RouteHandle) :: rh - type(ESMF_Field) :: lfield, field, field_dst - real(ESMF_KIND_R8), pointer :: ptr(:) + type(ESMF_Field) :: field_src, field_dst + real(ESMF_KIND_R8), pointer :: ptr_src(:), ptr_dst(:) integer :: n - character(len=cs), allocatable :: flds(:) + character(len=cs) :: flds_name(2) + type(field_type) :: flds(1) character(len=*), parameter :: subname = trim(modName)//': (read_initial) ' !------------------------------------------------------------------------------- @@ -121,10 +116,14 @@ subroutine read_initial(gcomp, ini_file, mosaic_file, input_dir, layout, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return ! --------------------- - ! Create domain + ! Set domain specific parameters + ! TODO: This assumes global domain with six tile and needs to be + ! revisited to support regional apps with one tile ! --------------------- - call create_fms_domain(gcomp, domain, mosaic_file, layout, rc) + domain%ntiles = 6 + domain%layout(1) = layout(1) + domain%layout(2) = layout(2) ! --------------------- ! Create grid @@ -132,59 +131,88 @@ subroutine read_initial(gcomp, ini_file, mosaic_file, input_dir, layout, rc) call create_grid(gcomp, domain, mosaic_file, input_dir, rc) + ! --------------------- + ! Create field in source mesh + ! --------------------- + + ! create field + field_src = ESMF_FieldCreate(domain%mesh, ESMF_TYPEKIND_R8, name='field_src', & + meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! get pointer and init it + call ESMF_FieldGet(field_src, localDe=0, farrayPtr=ptr_src, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ptr_src(:) = 0.0_r8 + + ! --------------------- + ! Create field in destination mesh + ! --------------------- + + ! create destination field + field_dst = ESMF_FieldCreate(is_local%wrap%aoflux_mesh, ESMF_TYPEKIND_R8, & + name='field_dst', meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! get pointer and init it + call ESMF_FieldGet(field_dst, localDe=0, farrayPtr=ptr_dst, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ptr_src(:) = 0.0_r8 + + ! --------------------- + ! Create routehandle + ! --------------------- + + call ESMF_FieldRegridStore(field_src, field_dst, routehandle=rh, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + !---------------------- ! Read data !---------------------- - allocate(flds(2)) - flds = (/ 'zorl ', & - 'uustar' /) - do n = 1,size(flds) - ! read from tiled file - call read_tiled_file(gcomp, ini_file, trim(flds(n)), domain, field, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - ! create destination field - field_dst = ESMF_FieldCreate(is_local%wrap%aoflux_mesh, ESMF_TYPEKIND_R8, & - name=trim(flds(n)), meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! list of fields that need to be read + flds_name(1) = 'zorl' + flds_name(2) = 'uustar' - ! create rh - call ESMF_FieldRegridStore(field, field_dst, routehandle=rh, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! loop over fields and read them + do n = 1,size(flds) + ! read data + flds(1)%short_name = trim(flds_name(n)) + flds(1)%ptr1r8 => ptr_src + call read_tiled_file(domain, ini_file, flds, rc=rc) ! map field if (is_local%wrap%aoflux_grid == 'agrid') then - ! do nothing, just redist in case of haning different decomp. in here and aoflux mesh - call ESMF_FieldRedist(field, field_dst, rh, rc=rc) + ! do nothing, just redist in case of having different decomp. in here and aoflux mesh + call ESMF_FieldRedist(field_src, field_dst, rh, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return else ! remap from atm to ocn or exchange grid - call ESMF_FieldRegrid(field, field_dst, rh, termorderflag=ESMF_TERMORDER_SRCSEQ, zeroregion=ESMF_REGION_TOTAL, rc=rc) + call ESMF_FieldRegrid(field_src, field_dst, rh, termorderflag=ESMF_TERMORDER_SRCSEQ, zeroregion=ESMF_REGION_TOTAL, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if ! debug if (dbug_flag > 5) then - call ESMF_FieldWriteVTK(field_dst, 'ini_'//trim(flds(n))//'_'//trim(is_local%wrap%aoflux_grid), rc=rc) + call ESMF_FieldWriteVTK(field_dst, 'ini_'//trim(flds_name(n))//'_'//trim(is_local%wrap%aoflux_grid), rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if - ! return pointer and fill variable - call ESMF_FieldGet(field_dst, localDe=0, farrayPtr=ptr, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (maintask) write(logunit,'(a)') 'Reading: '//trim(flds(n)) - if (trim(flds(n)) == 'zorl' ) physics%sfcprop%zorl(:) = ptr(:) - if (trim(flds(n)) == 'uustar') physics%sfcprop%uustar(:)= ptr(:) - nullify(ptr) - - ! free memory - call ESMF_FieldDestroy(field_dst, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! fill variables + if (maintask) write(logunit,'(a)') 'Reading: '//trim(flds_name(n)) + if (trim(flds_name(n)) == 'zorl' ) physics%sfcprop%zorl(:) = ptr_dst(:) + if (trim(flds_name(n)) == 'uustar') physics%sfcprop%uustar(:)= ptr_dst(:) end do - ! free memory - if (allocated(flds)) deallocate(flds) + !---------------------- + ! Free memory + !---------------------- + + call ESMF_FieldDestroy(field_src, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_FieldDestroy(field_dst, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_LogWrite(subname//' done', ESMF_LOGMSG_INFO) @@ -333,181 +361,6 @@ subroutine read_restart(gcomp, rst_file, rc) end subroutine read_restart - !=============================================================================== - subroutine create_fms_domain(gcomp, domain, mosaic_file, layout, rc) - implicit none - - ! input/output variables - type(ESMF_GridComp), intent(in) :: gcomp - type(domain_type), intent(inout) :: domain - character(len=cl), intent(in) :: mosaic_file - integer :: layout(2) - integer, intent(inout) :: rc - - ! local variables - type(ESMF_VM) :: vm - type(FmsNetcdfFile_t) :: mosaic_fileobj - integer :: mpicomm, npes_per_tile - integer :: n, ntiles, npet - integer :: halo = 0 - integer :: global_indices(4,6) - integer :: layout2d(2,6) - integer, allocatable :: pe_start(:), pe_end(:) - character(len=cl) :: msg - character(len=*), parameter :: subname = trim(modName)//': (create_fms_domain) ' - !------------------------------------------------------------------------------- - - rc = ESMF_SUCCESS - call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO) - - ! --------------------- - ! Initialize FMS - ! --------------------- - - call ESMF_GridCompGet(gcomp, vm=vm, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - call ESMF_VMGet(vm=vm, mpiCommunicator=mpicomm, petCount=npet, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - call fms_init(mpicomm) - - ! --------------------- - ! Open mosaic file and query some information - ! --------------------- - - if (.not. open_file(mosaic_fileobj, trim(mosaic_file), 'read')) then - call ESMF_LogWrite(trim(subname)//'error in opening file '//trim(mosaic_file), ESMF_LOGMSG_ERROR) - rc = ESMF_FAILURE - return - end if - - ! query number of tiles - domain%ntiles = get_mosaic_ntiles(mosaic_fileobj) - - ! query domain sizes for each tile - if (.not. allocated(domain%nit)) allocate(domain%nit(domain%ntiles)) - if (.not. allocated(domain%njt)) allocate(domain%njt(domain%ntiles)) - call get_mosaic_grid_sizes(mosaic_fileobj, domain%nit, domain%njt) - - ! query number of contacts - domain%ncontacts = get_mosaic_ncontacts(mosaic_fileobj) - - ! allocate required arrays to create FMS domain from mosaic file - if (.not. allocated(domain%tile1)) allocate(domain%tile1(domain%ncontacts)) - if (.not. allocated(domain%tile2)) allocate(domain%tile2(domain%ncontacts)) - if (.not. allocated(domain%istart1)) allocate(domain%istart1(domain%ncontacts)) - if (.not. allocated(domain%iend1)) allocate(domain%iend1(domain%ncontacts)) - if (.not. allocated(domain%jstart1)) allocate(domain%jstart1(domain%ncontacts)) - if (.not. allocated(domain%jend1)) allocate(domain%jend1(domain%ncontacts)) - if (.not. allocated(domain%istart2)) allocate(domain%istart2(domain%ncontacts)) - if (.not. allocated(domain%iend2)) allocate(domain%iend2(domain%ncontacts)) - if (.not. allocated(domain%jstart2)) allocate(domain%jstart2(domain%ncontacts)) - if (.not. allocated(domain%jend2)) allocate(domain%jend2(domain%ncontacts)) - - ! query information about contacts - call get_mosaic_contact(mosaic_fileobj, domain%tile1, domain%tile2, & - domain%istart1, domain%iend1, domain%jstart1, domain%jend1, & - domain%istart2, domain%iend2, domain%jstart2, domain%jend2) - - ! print out debug information - if (dbug_flag > 2) then - do n = 1, domain%ncontacts - write(msg, fmt='(A,I2,A,2I5)') trim(subname)//' : tile1, tile2 (', n ,') = ', domain%tile1(n), domain%tile2(n) - call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO) - write(msg, fmt='(A,I2,A,4I5)') trim(subname)//' : istart1, iend1, jstart1, jend1 (', n ,') = ', & - domain%istart1(n), domain%iend1(n), domain%jstart1(n), domain%jend1(n) - call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO) - write(msg, fmt='(A,I2,A,4I5)') trim(subname)//' : istart2, iend2, jstart2, jend2 (', n ,') = ', & - domain%istart2(n), domain%iend2(n), domain%jstart2(n), domain%jend2(n) - call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO) - end do - end if - - !---------------------- - ! Initialize domain - !---------------------- - - call mpp_domains_init() - - !---------------------- - ! Find out layout that will be used to read the data - !---------------------- - - ! setup global indices - do n = 1, domain%ntiles - global_indices(1,n) = 1 - global_indices(2,n) = domain%nit(n) - global_indices(3,n) = 1 - global_indices(4,n) = domain%njt(n) - end do - - ! check total number of PETs - if (mod(npet, domain%ntiles) /= 0) then - write(msg, fmt='(A,I5)') trim(subname)//' : nPet should be multiple of 6 to read initial conditions but it is ', npet - call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_ERROR) - rc = ESMF_FAILURE - return - end if - - ! calculate layout if it is not provided as configuration option - if (layout(1) < 0 .and. layout(2) < 0) then - npes_per_tile = npet/domain%ntiles - call mpp_define_layout(global_indices(:,1), npes_per_tile, domain%layout) - else - domain%layout(:) = layout(:) - end if - - ! set layout and print out debug information - do n = 1, domain%ntiles - layout2d(:,n) = domain%layout(:) - if (dbug_flag > 2) then - write(msg, fmt='(A,I2,A,2I5)') trim(subname)//' layout (', n ,') = ', layout2d(1,n), layout2d(2,n) - call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO) - write(msg, fmt='(A,I2,A,4I5)') trim(subname)//' global_indices (', n,') = ', & - global_indices(1,n), global_indices(2,n), global_indices(3,n), global_indices(4,n) - call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO) - end if - enddo - - !---------------------- - ! Set pe_start, pe_end - !---------------------- - - allocate(pe_start(domain%ntiles)) - allocate(pe_end(domain%ntiles)) - do n = 1, domain%ntiles - pe_start(n) = mpp_root_pe()+(n-1)*domain%layout(1)*domain%layout(2) - pe_end(n) = mpp_root_pe()+n*domain%layout(1)*domain%layout(2)-1 - if (dbug_flag > 2) then - write(msg, fmt='(A,I2,A,2I5)') trim(subname)//' pe_start, pe_end (', n ,') = ', pe_start(n), pe_end(n) - call ESMF_LogWrite(trim(msg), ESMF_LOGMSG_INFO) - end if - enddo - - !---------------------- - ! Create FMS domain object - !---------------------- - - call mpp_define_mosaic(global_indices, layout2d, domain%mosaic_domain, & - domain%ntiles, domain%ncontacts, domain%tile1, domain%tile2, & - domain%istart1, domain%iend1, domain%jstart1, domain%jend1, & - domain%istart2, domain%iend2, domain%jstart2, domain%jend2, & - pe_start, pe_end, symmetry=.true., & - whalo=halo, ehalo=halo, shalo=halo, nhalo=halo, & - name='atm domain') - - !---------------------- - ! Deallocate temporary arrays - !---------------------- - - deallocate(pe_start) - deallocate(pe_end) - - call ESMF_LogWrite(subname//' done', ESMF_LOGMSG_INFO) - - end subroutine create_fms_domain - !=============================================================================== subroutine create_grid(gcomp, domain, mosaic_file, input_dir, rc) implicit none @@ -553,148 +406,316 @@ subroutine create_grid(gcomp, domain, mosaic_file, input_dir, rc) end subroutine create_grid !=============================================================================== - subroutine read_tiled_file(gcomp, filename, varname, domain, field_dst, rc) - implicit none + subroutine read_tiled_file(domain, filename, flds, rh, rc) ! input/output variables - type(ESMF_GridComp), intent(in) :: gcomp - character(len=*), intent(in) :: filename - character(len=*), intent(in) :: varname - type(domain_type), intent(inout) :: domain - type(ESMF_Field), intent(inout) :: field_dst - integer, intent(inout), optional :: rc + type(domain_type), intent(inout) :: domain + character(len=*), intent(in) :: filename + type(field_type), intent(in) :: flds(:) + type(ESMF_RouteHandle), optional, intent(in) :: rh + integer, optional, intent(inout) :: rc ! local variables - type(ESMF_Field) :: field_src, field_tmp + integer :: i, j, k, rank, fieldCount + integer, pointer :: ptr_i4(:) + real(r4), pointer :: ptr_r4(:) + real(r8), pointer :: ptr_r8(:) + type(ESMF_RouteHandle) :: rh_local + type(ESMF_FieldBundle) :: FBgrid, FBmesh type(ESMF_ArraySpec) :: arraySpec - type(InternalState) :: is_local - type(fieldtype), allocatable:: vars(:) - integer :: funit, my_tile - integer :: i, j, n - integer :: isc, iec, jsc, jec - integer :: ndim, nvar, natt, ntime - logical :: not_found, is_root_pe - real(ESMF_KIND_R8), pointer :: ptr2d(:,:) - real(r8), allocatable :: rdata(:,:) - character(len=cl) :: cname - character(len=*), parameter :: subname=trim(modName)//': (read_tiled_file) ' + type(ESMF_Field) :: fgrid, fmesh, ftmp + character(len=cl) :: fname + character(len=cl), allocatable :: fieldNameList(:) + character(len=*), parameter :: subname = trim(modName)//': (read_tiled_file) ' !------------------------------------------------------------------------------- rc = ESMF_SUCCESS - call ESMF_LogWrite(subname//' reading '//trim(varname), ESMF_LOGMSG_INFO) + call ESMF_LogWrite(subname//' called for '//trim(filename), ESMF_LOGMSG_INFO) !---------------------- - ! Get the internal state from the mediator component + ! Create field bundles !---------------------- - nullify(is_local%wrap) - call ESMF_GridCompGetInternalState(gcomp, is_local, rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + ! create empty field bundle on grid + FBgrid = ESMF_FieldBundleCreate(name="fields_on_grid", rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create empty field bundle on mesh + FBmesh = ESMF_FieldBundleCreate(name="fields_on_mesh", rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return !---------------------- - ! Set tile + ! Loop over fields and add them to the field bundles !---------------------- - my_tile = int(mpp_pe()/(domain%layout(1)*domain%layout(2)))+1 - is_root_pe = .false. - if (mpp_pe() == (my_tile-1)*(domain%layout(1)*domain%layout(2))) is_root_pe = .true. + do i = 1, size(flds) + ! 2d/r8 field (x,y) + if (associated(flds(i)%ptr1r8)) then + ! set field type + call ESMF_ArraySpecSet(arraySpec, typekind=ESMF_TYPEKIND_R8, rank=2, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create field on grid + fgrid = ESMF_FieldCreate(domain%grid, arraySpec, staggerloc=ESMF_STAGGERLOC_CENTER, & + indexflag=ESMF_INDEX_GLOBAL, name=trim(flds(i)%short_name), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create field on mesh + fmesh = ESMF_FieldCreate(domain%mesh, flds(i)%ptr1r8, meshloc=ESMF_MESHLOC_ELEMENT, & + name=trim(flds(i)%short_name), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! 2d/r4 field (x,y) + else if (associated(flds(i)%ptr1r4)) then + ! set field type + call ESMF_ArraySpecSet(arraySpec, typekind=ESMF_TYPEKIND_R4, rank=2, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create field on grid + fgrid = ESMF_FieldCreate(domain%grid, arraySpec, staggerloc=ESMF_STAGGERLOC_CENTER, & + indexflag=ESMF_INDEX_GLOBAL, name=trim(flds(i)%short_name), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create field on mesh + fmesh = ESMF_FieldCreate(domain%mesh, flds(i)%ptr1r4, meshloc=ESMF_MESHLOC_ELEMENT, & + name=trim(flds(i)%short_name), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! 2d/i4 field (x,y) + else if (associated(flds(i)%ptr1i4)) then + ! set field type + call ESMF_ArraySpecSet(arraySpec, typekind=ESMF_TYPEKIND_I4, rank=2, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create field on grid + fgrid = ESMF_FieldCreate(domain%grid, arraySpec, staggerloc=ESMF_STAGGERLOC_CENTER, & + indexflag=ESMF_INDEX_GLOBAL, name=trim(flds(i)%short_name), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create field on mesh + fmesh = ESMF_FieldCreate(domain%mesh, flds(i)%ptr1i4, meshloc=ESMF_MESHLOC_ELEMENT, & + name=trim(flds(i)%short_name), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! 3d/r8 field (x,y,rec) + else if (associated(flds(i)%ptr2r8)) then + ! set field type + call ESMF_ArraySpecSet(arraySpec, typekind=ESMF_TYPEKIND_R8, rank=3, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create field on grid + fgrid = ESMF_FieldCreate(domain%grid, arraySpec, staggerloc=ESMF_STAGGERLOC_CENTER, & + indexflag=ESMF_INDEX_GLOBAL, name=trim(flds(i)%short_name), ungriddedLbound=(/1/), & + ungriddedUbound=(/flds(i)%nrec/), gridToFieldMap=(/1,2/), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create field on mesh + fmesh = ESMF_FieldCreate(domain%mesh, flds(i)%ptr2r8, meshloc=ESMF_MESHLOC_ELEMENT, & + name=trim(flds(i)%short_name), gridToFieldMap=(/1/), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! 3d/r4 field (x,y,rec) + else if (associated(flds(i)%ptr2r4)) then + ! set field type + call ESMF_ArraySpecSet(arraySpec, typekind=ESMF_TYPEKIND_R4, rank=3, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create field on grid + fgrid = ESMF_FieldCreate(domain%grid, arraySpec, staggerloc=ESMF_STAGGERLOC_CENTER, & + indexflag=ESMF_INDEX_GLOBAL, name=trim(flds(i)%short_name), ungriddedLbound=(/1/), & + ungriddedUbound=(/flds(i)%nrec/), gridToFieldMap=(/1,2/), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create field on mesh + fmesh = ESMF_FieldCreate(domain%mesh, flds(i)%ptr2r4, meshloc=ESMF_MESHLOC_ELEMENT, & + name=trim(flds(i)%short_name), gridToFieldMap=(/1/), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! 3d/i4 field (x,y,rec) + else if (associated(flds(i)%ptr2i4)) then + ! set field type + call ESMF_ArraySpecSet(arraySpec, typekind=ESMF_TYPEKIND_I4, rank=3, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create field on grid + fgrid = ESMF_FieldCreate(domain%grid, arraySpec, staggerloc=ESMF_STAGGERLOC_CENTER, & + indexflag=ESMF_INDEX_GLOBAL, name=trim(flds(i)%short_name), ungriddedLbound=(/1/), & + ungriddedUbound=(/flds(i)%nrec/), gridToFieldMap=(/1,2/), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create field on mesh + fmesh = ESMF_FieldCreate(domain%mesh, flds(i)%ptr2i4, meshloc=ESMF_MESHLOC_ELEMENT, & + name=trim(flds(i)%short_name), gridToFieldMap=(/1/), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + + ! debug print + call ESMF_LogWrite(trim(subname)//' adding '//trim(flds(i)%short_name)//' to FB', ESMF_LOGMSG_INFO) + + ! add it to the field bundle on grid + call ESMF_FieldBundleAdd(FBgrid, [fgrid], rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! add it to the field bundle on mesh + call ESMF_FieldBundleAdd(FBmesh, [fmesh], rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end do !---------------------- - ! Open file and query file attributes + ! Read data !---------------------- - - write(cname, fmt='(A,I1,A)') trim(filename), my_tile, '.nc' - call mpp_open(funit, trim(cname), action=MPP_RDONLY, form=MPP_NETCDF, threading=MPP_MULTI, fileset=MPP_SINGLE, is_root_pe=is_root_pe) - call mpp_get_info(funit, ndim, nvar, natt, ntime) - allocate(vars(nvar)) - call mpp_get_fields(funit, vars(:)) + + call ESMF_FieldBundleRead(FBgrid, fileName=trim(filename), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return !---------------------- - ! Find and read requested variable + ! Create routehandle if it is not provided to transfer data from grid to mesh !---------------------- - not_found = .true. - do n = 1, nvar - ! get variable name - call mpp_get_atts(vars(n), name=cname) + if (present(rh)) then + rh_local = rh + else + call ESMF_FieldBundleRedistStore(FBgrid, FBmesh, routehandle=rh_local, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if - ! check variable name - if (trim(cname) == trim(varname)) then - ! get array bounds or domain - call mpp_get_compute_domain(domain%mosaic_domain, isc, iec, jsc, jec) + !---------------------- + ! Move data from ESMF grid to mesh + !---------------------- - ! allocate data array and set initial value - allocate(rdata(isc:iec,jsc:jec)) - rdata(:,:) = 0.0_r8 + call ESMF_FieldBundleRedist(FBgrid, FBmesh, rh_local, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! read data - call mpp_read(funit, vars(n), domain%mosaic_domain, rdata, 1) + !---------------------- + ! Debug output + !---------------------- - ! set missing values to zero - where (rdata == 1.0e20) - rdata(:,:) = 0.0_r8 - end where - end if + if (dbug_flag > 5) then + do i = 1, size(flds) + ! get field from FB + call ESMF_FieldBundleGet(FBmesh, fieldName=trim(flds(i)%short_name), field=fmesh, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return - not_found = .false. - end do + ! check its rank + call ESMF_FieldGet(fmesh, rank=rank, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (not_found) then - call mpp_error(FATAL, 'File being read is not the expected one. '//trim(varname)//' is not found.') + ! TODO: ESMF_FieldWriteVTK() call does not support ungridded dimension + ! The workaround is implemented in here but it would be nice to extend + ! ESMF_FieldWriteVTK() call to handle it. + if (rank > 1) then + ! create temporary field + if (associated(flds(i)%ptr2r4)) then + ftmp = ESMF_FieldCreate(domain%mesh, typekind=ESMF_TYPEKIND_R4, & + name=trim(flds(i)%short_name), meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_FieldGet(ftmp, localDe=0, farrayPtr=ptr_r4, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else if (associated(flds(i)%ptr2r8)) then + ftmp = ESMF_FieldCreate(domain%mesh, typekind=ESMF_TYPEKIND_R8, & + name=trim(flds(i)%short_name), meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_FieldGet(ftmp, localDe=0, farrayPtr=ptr_r8, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else if (associated(flds(i)%ptr2i4)) then + ftmp = ESMF_FieldCreate(domain%mesh, typekind=ESMF_TYPEKIND_I4, & + name=trim(flds(i)%short_name), meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_FieldGet(ftmp, localDe=0, farrayPtr=ptr_i4, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + + ! write all record to seperate VTK file + do j = 1, flds(i)%nrec + if (associated(flds(i)%ptr2i4)) ptr_i4(:) = flds(i)%ptr2i4(:,j) + if (associated(flds(i)%ptr2r4)) ptr_r4(:) = flds(i)%ptr2r4(:,j) + if (associated(flds(i)%ptr2r8)) ptr_r8(:) = flds(i)%ptr2r8(:,j) + write(fname, fmt='(A,I2.2)') trim(flds(i)%short_name)//'_rec', j + call ESMF_FieldWriteVTK(ftmp, trim(fname), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end do + + ! delete temporary field + call ESMF_FieldDestroy(ftmp, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + ! write field to VTK file + call ESMF_FieldWriteVTK(fmesh, trim(flds(i)%short_name), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + end do end if !---------------------- - ! Move data from grid to mesh + ! Empty FBs and destroy them !---------------------- - ! set type and rank for ESMF arrayspec - call ESMF_ArraySpecSet(arraySpec, typekind=ESMF_TYPEKIND_R8, rank=2, rc=rc) + ! FB grid + call ESMF_FieldBundleGet(FBgrid, fieldCount=fieldCount, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! create source field - field_src = ESMF_FieldCreate(domain%grid, arraySpec, staggerloc=ESMF_STAGGERLOC_CENTER, & - indexflag=ESMF_INDEX_GLOBAL, name=trim(varname), rc=rc) + allocate(fieldNameList(fieldCount)) + call ESMF_FieldBundleGet(FBgrid, fieldNameList=fieldNameList, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! get pointer and fill it - call ESMF_FieldGet(field_src, localDe=0, farrayPtr=ptr2d, rc=rc) + do i = 1, fieldCount + ! pull field from FB + call ESMF_FieldBundleGet(FBgrid, fieldName=trim(fieldNameList(i)), field=ftmp, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! destroy field + call ESMF_FieldDestroy(ftmp, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! remove field from FB + call ESMF_FieldBundleRemove(FBgrid, fieldNameList=[trim(fieldNameList(i))], rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end do + deallocate(fieldNameList) + + ! destroy grid FB + call ESMF_FieldBundleDestroy(FBgrid, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ptr2d(:,:) = rdata(:,:) - nullify(ptr2d) - if (allocated(rdata)) deallocate(rdata) - ! create destination field - field_dst = ESMF_FieldCreate(domain%mesh, ESMF_TYPEKIND_R8, name=trim(varname), & - meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + ! FB mesh + call ESMF_FieldBundleGet(FBmesh, fieldCount=fieldCount, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + allocate(fieldNameList(fieldCount)) + call ESMF_FieldBundleGet(FBmesh, fieldNameList=fieldNameList, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! create routehandle from grid to mesh - if (.not. ESMF_RouteHandleIsCreated(domain%rh, rc=rc)) then - call ESMF_FieldRegridStore(field_src, field_dst, routehandle=domain%rh, rc=rc) + do i = 1, fieldCount + ! pull field from FB + call ESMF_FieldBundleGet(FBmesh, fieldName=trim(fieldNameList(i)), field=ftmp, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if - ! redist field from ESMF Grid to Mesh - call ESMF_FieldRedist(field_src, field_dst, domain%rh, rc=rc) + ! destroy field + call ESMF_FieldDestroy(ftmp, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! remove field from FB + call ESMF_FieldBundleRemove(FBmesh, fieldNameList=[trim(fieldNameList(i))], rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end do + deallocate(fieldNameList) + + ! destroy grid FB + call ESMF_FieldBundleDestroy(FBmesh, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return !---------------------- - ! Output result field for debugging purpose + ! Destroy route handle if it is created locally !---------------------- - if (dbug_flag > 2) then - call ESMF_FieldWrite(field_dst, trim(varname)//'_agrid.nc', variableName=trim(varname), overwrite=.true., rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if - - if (dbug_flag > 5) then - call ESMF_FieldWriteVTK(field_dst, trim(varname)//'_agrid', rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (.not. present(rh)) then + call ESMF_RouteHandleDestroy(rh_local, rc=rc) end if - ! clean memory - call ESMF_FieldDestroy(field_src, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_LogWrite(subname//' done', ESMF_LOGMSG_INFO) end subroutine read_tiled_file