diff --git a/coupling.txt b/coupling.txt new file mode 100644 index 000000000..62ac04fbf --- /dev/null +++ b/coupling.txt @@ -0,0 +1,23 @@ +(internalstate allowed coupling) Allowed coupling flags + |from to -> med atm lnd ocn ice rof wav glc1 glc2 + |med - - - - - - - - - + |atm - - T T T - T - - + |lnd - T - - - T - T T + |ocn - T - - T - T T T + |ice - T - T - - T - - + |rof - - T T T - - - - + |wav - T - T T - - - - + |glc1 - - T - T T - - - + |glc2 - - T - T T - - - + + (internalstate allowed coupling) Active coupling flags + |from to -> med atm lnd ocn ice rof wav glc1 glc2 + |med - - - - - - - - - + |atm - - T T T - - - - + |lnd - T - - - T - T T + |ocn - T - - T - - T T + |ice - T - T - - - - - + |rof - - T T T - - - - + |wav - - - - - - - - - + |glc1 - - T - T T - - - + |glc2 - - T - T T - - - diff --git a/doc/source/developer/add_field.rst b/doc/source/developer/add_field.rst new file mode 100644 index 000000000..e5d47d21e --- /dev/null +++ b/doc/source/developer/add_field.rst @@ -0,0 +1,94 @@ +.. _add-field: + +====================== +Adding a coupled field +====================== + +Adding a new field for the mediator to exchange between two components is one of +the most common CMEPS tasks, and it can be done on its own without touching the +rest of the mediator. This page is a practical recipe. It assumes the +:ref:`field-exchange` concepts (advertise/realize, ``addfld`` / ``addmap`` / +``addmrg``) and the :ref:`fields` naming convention. + +What you will touch +=================== + +A new coupled field involves a handful of well-defined edits: + +#. the **field dictionary** — declare the field's standard name and units; +#. the host **exchange module** — advertise, map and merge it in + ``esmFldsExchange__mod.F90``; +#. the **component caps** — make sure each component actually advertises the + field so NUOPC connects it. + +The first two are inside CMEPS; the third is on the component side (and outside +CMEPS), but the connection will not form without it. + +Step 1 — add the field to the dictionary +======================================== + +Add an entry for the field to the host field dictionary (``fd_cesm.yaml`` for +CESM/NorESM). Give it a standard name that follows the :ref:`naming convention +`, its canonical units, and optionally an alias and description: + +.. code-block:: yaml + + - standard_name: Sa_foo + canonical_units: K + description: my new near-surface atmosphere state + +Step 2 — advertise, map and merge it +==================================== + +All of this goes in the host exchange module, ``esmFldsExchange__mod.F90``. + +**Advertise** the field on both sides, in the ``phase == 'advertise'`` block — +from the source component and to the destination component: + +.. code-block:: fortran + + call addfld_from(compatm, 'Sa_foo') ! atmosphere can send it + call addfld_to(complnd, 'Sa_foo') ! land can receive it + +**Map** the field from the source grid to the destination grid, in the +``phase == 'initialize'`` block. Choose the map type and normalization +appropriate to the field (see :ref:`mapping`): + +.. code-block:: fortran + + call addmap_from(compatm, 'Sa_foo', complnd, mapbilnr, 'none', 'unset') + +**Merge** the mapped field into the destination's export, also in the +``initialize`` block. For a field that comes from a single source, a plain +``copy`` is enough; a field built from several surfaces uses one ``addmrg_to`` +call per source with a fraction weight (see :ref:`merging`): + +.. code-block:: fortran + + call addmrg_to(complnd, 'Sa_foo', mrg_from=compatm, mrg_fld='Sa_foo', & + mrg_type='copy') + +Step 3 — advertise the field in the component caps +================================================== + +The mediator advertises a superset of possible fields, but a connection is only +made where the component on the other side advertises the field too +(:ref:`field-exchange`). So the field must also be advertised by the component +caps — the atmosphere cap must advertise ``Sa_foo`` as an **export** and fill it, +and the land cap must advertise it as an **import** and use it. This is +component code, outside CMEPS, but the coupling will silently not form without +it. + +Checklist +========= + +* Field added to ``fd_.yaml`` with correct name and units. +* ``addfld_from`` / ``addfld_to`` in the ``advertise`` phase. +* ``addmap_from`` with an appropriate map type and normalization. +* ``addmrg_to`` with the appropriate merge type (and fraction, for weighted + merges). +* The source and destination **component caps** advertise the field. + +If the field does not appear in the destination component, the first things to +check are that both caps advertised it (so NUOPC connected it) and that the name +matches the dictionary exactly. diff --git a/doc/source/developer/aofluxes.rst b/doc/source/developer/aofluxes.rst new file mode 100644 index 000000000..427532076 --- /dev/null +++ b/doc/source/developer/aofluxes.rst @@ -0,0 +1,100 @@ +.. _aofluxes: + +======================= +Atmosphere/ocean fluxes +======================= + +The mediator computes the turbulent air-sea fluxes itself, from the atmosphere +and ocean states brought onto a common grid, using a bulk flux algorithm. The +**concept** — +why the mediator owns this calculation, which fields it produces, and how they +are delivered and averaged — is covered in :ref:`concepts`. This page describes +the **mechanism**: the interface in ``med_phases_aofluxes_mod.F90`` and the grid +on which the calculation is carried out. + +Interface versus algorithm +========================== + +``med_phases_aofluxes_mod.F90`` is an **interface**: it prepares the inputs, +selects the grid, invokes the flux algorithm, and distributes the outputs. The +bulk-flux **algorithm** itself lives outside the shared mediator, in +host-specific code — ``cesm/flux_atmocn`` for CESM/NorESM and ``ufs/`` for UFS +(see :ref:`concepts`). The mediator reaches it through +``med_aofluxes_init`` / ``med_aofluxes_update``. + +The flux outputs are held in two internal-state field bundles, +``FBMed_aoflux_a`` and ``FBMed_aoflux_o`` — the flux fields on the atmosphere and +ocean grids. + +The run phase +============= + +``med_phases_aofluxes_run`` is the phase the run sequence calls +(``MED med_phases_aofluxes_run``). On the first call it initializes the flux +input/output field bundles (``med_phases_aofluxes_init_fldbuns``) and the flux +scheme; on every call it computes the fluxes with ``med_aofluxes_update``. + +Which bulk-flux scheme is used, and its options, come from run-configuration +attributes: + +* ``coupling_mode`` — the target modeling system that will compute the + atmosphere-ocean fluxes. Valid values are ``cesm``, ``noresm`` and ``ufs``. +* ``aoflux_grid`` — the grid the atmosphere-ocean flux computation is done on: + ``agrid`` (atmosphere), ``ogrid`` (ocean), or ``xgrid`` (the ESMF exchange + grid). This is described in detail in `The flux calculation grid`_ below. +* ``ocn_surface_flux_scheme`` — selects the ocean surface flux (bulk) + algorithm: + + * ``0`` — Large and Pond + * ``1`` — COARE algorithm + * ``2`` — UA algorithm + + CESM and UFS use the Large and Pond scheme (``0``); NorESM uses the COARE + scheme (``1``). +* ``aofluxes_use_shr_wv_sat`` — only applicable to CESM and NorESM. If + ``true``, ``shr_wv_sat_mod`` is used to calculate the saturation specific + humidity (``qsat``) for the atmosphere-ocean flux calculations; if ``false``, + an older inline calculation of ``qsat`` is used, which uses a different + formulation. +* ``add_gusts`` — only applicable to CESM and NorESM. Adds a wind gustiness + factor. This should be ``false`` for ``ocn_surface_flux_scheme`` settings of + ``1`` or ``2``. +* ``flux_max_iteration`` — the maximum number of iterations of the (iterative) + flux solver. + +The flux calculation grid +========================= + +The fluxes may be computed on the atmosphere grid, the ocean grid, or a dedicated +**ESMF exchange grid**, selected by the ``aoflux_grid`` attribute; for cases with +**all prognostic components, the exchange grid is the default**. Because the +calculation needs both atmosphere and ocean state on a single grid, the interface +maps the inputs onto the chosen grid and maps the outputs back out. The three +cases: + +``aoflux_grid = ogrid`` (ocean grid) + Requires an atmosphere-to-ocean mapping: map the atmosphere inputs to the + ocean grid, compute the fluxes there, and map the flux outputs back to the + atmosphere grid. + +``aoflux_grid = agrid`` (atmosphere grid) + Requires an ocean-to-atmosphere mapping: map the ocean inputs to the + atmosphere grid, compute the fluxes there, and map the flux outputs back to + the ocean grid. + +``aoflux_grid = xgrid`` (ESMF exchange grid) + Map **both** the atmosphere and ocean inputs onto the exchange grid, compute + the fluxes there, and map the flux outputs from the exchange grid to **both** + the atmosphere and ocean grids. + +The **ESMF exchange grid** is an ESMF construct formed from the overlap of the +atmosphere and ocean grids; the interface holds dedicated routehandles for it — +``rh_agrid2xgrid``, ``rh_ogrid2xgrid``, ``rh_xgrid2agrid`` and ``rh_xgrid2ogrid`` +(with second-order-conservative and bilinear atmosphere-to-exchange variants). +Computing on the exchange grid avoids biasing the fluxes toward either +component's grid, at the cost of the extra mappings. + +Whichever grid is chosen, the resulting fluxes end up in ``FBMed_aoflux_a`` and +``FBMed_aoflux_o`` for the ``prep_atm`` and ``prep_ocn`` phases to merge into the +atmosphere and ocean exports — on each component's grid and coupling interval, as +described in :ref:`concepts`. diff --git a/doc/source/developer/code_organization.rst b/doc/source/developer/code_organization.rst new file mode 100644 index 000000000..291dec74b --- /dev/null +++ b/doc/source/developer/code_organization.rst @@ -0,0 +1,198 @@ +.. _code-organization: + +================= +Code organization +================= + +This page describes how the mediator code is laid out and introduces the +**internal state** — the data structure at the heart of the mediator. It builds +on the high-level grouping in :ref:`architecture`; read that first if you have +not. As noted in the :ref:`Developer Guide introduction `, this page +assumes familiarity with ESMF and the NUOPC Layer. + +The ``mediator/`` directory +=========================== + +All of the shared mediator code lives in ``mediator/``. It falls into three +groups (described in the :ref:`architecture overview `): + +* **Generic infrastructure** — the reusable machinery: + + * the mediator component (``med.F90``) + * the internal state (``med_internalstate_mod.F90``) + * mapping (``med_map_mod.F90``) + * merging (``med_merge_mod.F90``) + * fractions (``med_fraction_mod.F90``) + * I/O (``med_io_mod.F90``) + * history output (``med_phases_history_mod.F90``) + * restarts (``med_phases_restart_mod.F90``) + * time management (``med_time_mod.F90``) + * methods and helpers (``med_methods_mod.F90``) + * shared utilities (``med_utils_mod.F90``) + * constants (``med_constants_mod.F90``) + * kind definitions (``med_kind_mod.F90``) +* **Application-specific exchange code** — decides *which* fields are exchanged + and *how* they are mapped and merged: + + * the per-host exchange modules (``esmFldsExchange__mod.F90``) + * shared exchange definitions (``esmFlds.F90``) + * the per-host field dictionaries (``fd_.yaml``) + +* **Phase code** — the per-component and cross-cutting work run each coupling + step: + + * per-component prep phases (``med_phases_prep__mod.F90``) + * per-component post phases (``med_phases_post__mod.F90``) + * the atmosphere/ocean flux interface (``med_phases_aofluxes_mod.F90``) + * the ocean albedo phase (``med_phases_ocnalb_mod.F90``) + * the diagnostics phase (``med_diag_mod.F90``) + * the timing/profiling phase (``med_phases_profile_mod.F90``) + +The host-specific driver (under ``cesm/``) and the UFS-specific flux code (under +``ufs/``) live outside ``mediator/`` and are not part of the shared mediator; see +:ref:`concepts` for the driver/mediator distinction. + +The mediator component +====================== + +``med.F90`` is the mediator's NUOPC component — a ``NUOPC_Mediator`` (its +``SetServices`` specializes the generic ``NUOPC_Mediator``). It registers the +initialization and run phases with NUOPC (the phases the run sequence names), +and it owns the mediator's **internal state**. The registration and the +initialization sequence are covered on the phase-lifecycle page (next); this +page focuses on the data the phases operate on. + +The internal state +================== + +The mediator's working data is held in a single **internal state** derived type +(defined in ``med_internalstate_mod.F90`` and carried on the component as +``is_local%wrap``). Understanding its field bundles is the key to reading almost +any phase. + +For a run with ``ncomps`` components, the central members are: + +Per-component states + * ``NStateImp(n)`` — the ESMF **import state** from component ``n``. + * ``NStateExp(n)`` — the ESMF **export state** to component ``n``. + +The import/export field bundles + * ``FBImp(n,k)`` — a **matrix** of field bundles: the fields imported from + component ``n``, interpolated to the grid of component ``k``. The diagonal + ``FBImp(n,n)`` points into component ``n``'s fields on its own grid (the + fields in ``NStateImp(n)``); the off-diagonal ``FBImp(n,k)`` holds those + fields after mapping onto grid ``k``. This matrix is what lets a ``prep`` + phase pull any source component's fields already on the destination grid. + * ``FBExp(n)`` — the field bundle of fields to be exported to component + ``n``, on component ``n``'s grid; its fields point into ``NStateExp(n)``, + the export state the connector sends to component ``n``. + +Mediator-computed field bundles + * ``FBMed_aoflux_a`` / ``FBMed_aoflux_o`` — the atmosphere/ocean flux outputs + on the atmosphere and ocean grids. + * ``FBMed_ocnalb_a`` / ``FBMed_ocnalb_o`` — the ocean albedos on the + atmosphere and ocean grids. + +Fractions and areas + * ``FBfrac(n)`` — the surface fractions for component ``n`` on its grid. + * ``FBArea(n)`` — component areas (used by the mediator history writes). + +Mapping + * ``RH(:,:,:)`` — the **routehandles** for each pair of components and each + mapping type (conservative, bilinear, and so on). These are created once + during initialization and reused every coupling step. + +Accumulators + * ``FBExpAccumOcn`` / ``FBExpAccumWav`` (with their counters + ``ExpAccumOcnCnt`` / ``ExpAccumWavCnt``) — accumulators used to average an + export over several fast coupling intervals before it is sent to a + component coupled less frequently (the accumulate-then-average pattern seen + in the :ref:`run sequence ` for the ocean). + +The coupling matrix +------------------- + +Which components exchange fields — and therefore which off-diagonal +``FBImp(n,k)`` entries are populated and which ``prep`` / ``post`` work runs — is +captured by a **coupling matrix** that the mediator prints at start-up. Rows are +the *source* (``from``) component and columns are the *destination* (``to``) +component; a ``T`` marks an active coupling from the row component to the column +component. + +CMEPS prints two versions, which mirror the :ref:`advertise/realize +` superset-vs-intersection idea. + +The **allowed** coupling matrix is every coupling CMEPS supports (the superset): + +.. code-block:: none + + from to -> med atm lnd ocn ice rof wav glc1 glc2 + med - - - - - - - - - + atm - - T T T - T - - + lnd - T - - - T - T T + ocn - T - - T - T T T + ice - T - T - - T - - + rof - - T T T - - - - + wav - T - T T - - - - + glc1 - - T - T T - - - + glc2 - - T - T T - - - + +The **active** coupling matrix is the subset actually turned on for a given +configuration — mirroring how NUOPC realizes only the connected fields. Two +examples follow. + +A fully coupled system without an active wave component (its row and column drop +out): + +.. code-block:: none + + from to -> med atm lnd ocn ice rof wav glc1 glc2 + med - - - - - - - - - + atm - - T T T - - - - + lnd - T - - - T - T T + ocn - T - - T - - T T + ice - T - T - - - - - + rof - - T T T - - - - + wav - - - - - - - - - + glc1 - - T - T T - - - + glc2 - - T - T T - - - + +A configuration with active **CAM**, **CLM** and **CICE** (atm, lnd, ice) and +**data** ocean (DOCN) and land-ice (DGLC): + +.. code-block:: none + + from to -> med atm lnd ocn ice rof wav glc1 glc2 + med - - - - - - - - - + atm - - T - T - - - - + lnd - T - - - T - T T + ocn - T - - T - - - - + ice - T - - - - - - - + rof - - T - T - - - - + wav - - - - - - - - - + glc1 - - T - T T - - - + glc2 - - T - T T - - - + +Notice how the **data** components change the picture. Nothing needs to *force* +the data ocean, so ``atm -> ocn``, ``ice -> ocn`` and ``rof -> ocn`` all go +inactive — but DOCN still provides its prescribed state *to* the atmosphere and +ice (``ocn -> atm``, ``ocn -> ice`` stay active). A data component is coupled +where it is a *source*, and decoupled where it would be a *destination* of +fields it does not consume. + +Reading the matrix ties directly back to the internal state. Each active +coupling is stored as one off-diagonal ``FBImp`` entry: + +.. code-block:: none + + an active coupling atm -> ocn (a "T" in the matrix above) + is stored as FBImp(atm, ocn) = atm fields on the ocn grid + +and it ties back to the phases. For each active coupling ``X -> Y``: + +* ``post_X`` maps ``X``'s fields onto ``Y``'s grid (filling ``FBImp(X,Y)``), and +* ``prep_Y`` merges those mapped fields — with the other active sources into + ``Y`` — into ``FBExp(Y)``. + +So a **row** shows everything a component feeds (its ``post`` targets), and a +**column** shows everything that feeds a component (its ``prep`` sources). diff --git a/doc/source/developer/conservation.rst b/doc/source/developer/conservation.rst new file mode 100644 index 000000000..03e82a00a --- /dev/null +++ b/doc/source/developer/conservation.rst @@ -0,0 +1,179 @@ +.. _conservation: + +=============================================== +Conservation, area corrections and accumulation +=============================================== + +Conservation of mass and energy is a defining requirement of the mediator: a +flux leaving one component must deposit the same integrated quantity in the +component that receives it, even after being mapped and merged across mismatched +grids and averaged over mismatched coupling intervals. Achieving this ties +together several details that must be considered *together* — mapping weights, +surface fractions, normalization, model-versus-ESMF areas, and accumulation. + +CMEPS uses **unmasked, unnormalized conservative weights**, maps with a +**fraction-weighted, normalized** approach, and carries out merges with +**fraction weights**. This page develops why that combination conserves; it is +the principle that the :ref:`mapping`, :ref:`merging` and +:ref:`fractions ` pages that follow then implement, and it +expands the conservation summary in :ref:`concepts`. + +The example below uses a single atmosphere cell ``a`` that overlaps four +ocean/ice cells ``1..4``; ``A`` is an overlap area, ``w`` a mapping weight, ``f`` +a surface fraction, ``F`` a field, and subscripts ``l/o/i`` denote land, ocean +and ice. + +Conservative mapping weights +============================ + +CMEPS uses area-overlap conservative weights, + +.. code-block:: none + + w1 = A1/Aa, w2 = A2/Aa, w3 = A3/Aa, w4 = A4/Aa + +which, being fractions of the atmosphere cell area, always sum to one: + +.. code-block:: none + + w1 + w2 + w3 + w4 = 1 + +Applying the weights is a sparse matrix multiply, ``Xa = [W] Xo``. These same +weights map the ocean/ice mask to obtain the land fraction, and map the surface +fractions themselves (see :ref:`fractions `). + +Fraction-weighted normalized mapping +==================================== + +Mapping a *masked* or *fraction-weighted* field — for example an SST that only +exists where there is ocean — cannot be done with a plain weighted average. A +direct average pulls in undefined land values; a mask-weighted average no longer +has weights that sum to one, so it is biased. The correct approach uses the +dynamically varying fraction exposed to the atmosphere: + +.. code-block:: none + + fao = w1*fo1 + w2*fo2 + w3*fo3 + w4*fo4 + Fao = (fo1*w1*Fo1 + fo2*w2*Fo2 + fo3*w3*Fo3 + fo4*w4*Fo4) / fao + +that is, the field is fraction-weighted, mapped, and then normalized by the +mapped fraction. In practice this is the four-step process implemented in +:ref:`mapping`: + +.. code-block:: none + + Fo' = fo * Fo (a) fraction-weight the field + fao = w1*fo1 + w2*fo2 + w3*fo3 + w4*fo4 (b) map the fraction + Fao' = w1*Fo1' + w2*Fo2' + w3*Fo3' + w4*Fo4' (c) map the weighted field + Fao = Fao' / fao (d) normalize by the mapped fraction + +Steps (b) and (c) are the conservative sparse-matrix multiply; step (a) +fraction-weights the field and step (d) normalizes. + +**Why it conserves.** Think of ``F`` as a flux and form the associated quantity. +On the ocean grid, + +.. code-block:: none + + Qo = (fo1*Fo1*A1 + fo2*Fo2*A2 + fo3*Fo3*A3 + fo4*Fo4*A4) * dt + +and on the atmosphere grid it is the mapped flux times the mapped fraction times +the atmosphere area, + +.. code-block:: none + + Qa = fao * Fao * Aa * dt + +Substituting the definitions of ``fao`` and ``Fao`` above (and ``w_i = A_i/Aa``) +gives ``Qo = Qa``: the integrated quantity is unchanged by the mapping. + +Merging with fractions +====================== + +The field the atmosphere sees over the cell is the fraction-weighted merge of +the three surfaces, + +.. code-block:: none + + Fa = fal*Fal + fao*Fao + fai*Fai, with fal + fao + fai = 1 + +Substituting the fraction-weighted normalized maps for ``Fao`` and ``Fai`` and +expanding, the whole expression collapses to a plain conservative mapping of the +field already merged *on the ocean grid*: + +.. code-block:: none + + Fa = w1*(fl1*Fl1 + fo1*Fo1 + fi1*Fi1) + + w2*(fl2*Fl2 + fo2*Fo2 + fi2*Fi2) + + w3*(fl3*Fl3 + fo3*Fo3 + fi3*Fi3) + + w4*(fl4*Fl4 + fo4*Fo4 + fi4*Fi4) + +which is manifestly conservative. Notice that the normalization at the end of the +mapping (dividing by ``fao``) and the fraction weighting in the merge +(multiplying by ``fao``) **cancel**. CMEPS nonetheless keeps both, because the +normalized field ``Fao`` is a meaningful gridcell-average value — useful for +history and diagnostics and required whenever a field is mapped but *not* merged +(such as an interpolated SST) — whereas the un-normalized ``Fao'`` is a subarea +average that is far less intuitive. + +Area corrections +================ + +The conservative weights above assume the overlap areas are computed +consistently. But each component computes its gridcell areas as part of its own +discretization (great-circle versus lon/lat approximations, different Earth radii, +and so on), and these **model areas do not, in general, match the ESMF areas** +used to build the mapping weights. Two grids can share identical cell corners and +centers yet disagree on area. + +Because a flux carries a quantity proportional to the *model* area, the fluxes +must be corrected by the ratio of the model and ESMF areas. Outgoing fluxes are +multiplied by ``model_area / ESMF_area`` and incoming fluxes by +``ESMF_area / model_area``: + +.. code-block:: none + + F1' = (A1m/A1) * F1 ! correct the source flux out of component 1 + F2' = w1 * F1' ! map + F2 = F2' * (A2/A2m) ! correct into component 2 + + => Q2 = F2*A2m = A1m*F1 = Q1 + +Three properties matter in practice: + +* area corrections apply to **fluxes only** (not states); +* they are computed **once at initialization** and **do not vary in time**; and +* the model areas must be passed to the mediator at initialization so the + corrections can be formed. + +Accumulation and averaging +========================== + +Components coupled at different frequencies require accumulation: a field is +accumulated over the fast coupling intervals and averaged before it is handed to +a component running on a slower interval. Two rules keep this conservative. + +**Accumulate the fraction-weighted flux, not the pieces separately.** For a +fraction-weighted flux the accumulation must be of the product ``f*F``, because + +.. code-block:: none + + sum_n(f*F) != sum_n(f) * sum_n(F) + +Accumulating the fraction and the flux separately and multiplying afterwards does +**not** give the correct time average and breaks conservation. + +**Accumulation and mapping commute, so accumulate first and map once.** Both the +accumulation ``sum_n`` and the mapping ``map`` are linear operations, so + +.. code-block:: none + + Fo = 1/n * sum_n( mapa2o(fao_a * Fao_a) ) + = mapa2o( 1/n * sum_n(fao_a * Fao_a) ) + +The fraction-weighted flux ``fao_a*Fao_a`` can therefore be accumulated on the +source (atmosphere) grid over the fast intervals and mapped **only once**, on the +slow coupling interval — the pattern the :ref:`run sequence ` uses +for the ocean (accumulate every fast step, average and send on the ocean +interval). CMEPS maintains dedicated accumulators for this +(``FBExpAccumOcn`` / ``FBExpAccumWav``; see :ref:`code-organization`). diff --git a/doc/source/developer/diagnostics.rst b/doc/source/developer/diagnostics.rst new file mode 100644 index 000000000..d2197afe6 --- /dev/null +++ b/doc/source/developer/diagnostics.rst @@ -0,0 +1,95 @@ +.. _diagnostics: + +======================== +Water and energy budgets +======================== + +One of the mediator's roles is to check that the coupled system conserves water +and energy. ``med_diag_mod.F90`` computes spatial and time averages of the +fluxed quantities needed for these **budget diagnostics**. + +Sign convention +=============== + +Fluxes are accounted with a **positive-downward** sign convention and a fixed +component hierarchy — atm, glc, lnd, rof, ice, ocn. With this convention a flux +that one component gains is the same flux another loses, so the totals across the +coupled system should close. A budget that does not close points to a +conservation error somewhere in the coupling, which makes these diagnostics a +primary tool for validating conservation. + +The budget phases +================= + +The budget work is spread across several phases that the run sequence calls: + +* per-component accumulation phases — ``med_phases_diag_atm``, + ``med_phases_diag_lnd``, ``med_phases_diag_ocn``, and so on — accumulate each + component's fluxed quantities as they pass through the mediator; +* ``med_phases_diag_accum`` accumulates those contributions over the averaging + periods; and +* ``med_phases_diag_print`` prints the budget tables. + +``med_diag_init`` and ``med_diag_zero`` set up and reset the accumulators. + +The budget tables printed by ``med_phases_diag_print`` are written to files of +the form ``diags.log..``. + +Reading a budget table +====================== + +For each averaging period the mediator prints a set of budgets — a **net area** +budget (m²/m²), a **net heat** budget (W/m²) and a **net water** budget +(kg/m²/s, scaled by 1e6). In each table the **columns** are the per-component +contributions (``atm``, ``lnd``, ``rof``, ``ocn``, ``ice nh``, ``ice sh``, +``glc``) plus a ``*SUM*`` column, and the **rows** are the individual flux terms +plus a ``*SUM*`` row. + +Conservation shows up directly in the sums: because a flux one component gains is +another's loss, the ``*SUM*`` column (across components) for each term, and the +``*SUM*`` row, should be **essentially zero**. A term whose row sum is far from +zero points at a non-conserving exchange. For example, from a NorESM run +(``diags.log`` excerpt), the annual **heat** budget (W/m²): + +.. code-block:: none + + (med_diag_print_summary) NET HEAT BUDGET (W/m2): period = annual: date = 11060101 0 + atm lnd rof ocn ice nh ice sh glc *SUM* + hfreeze 0.00000000 0.00000000 0.00000000 0.03292080 -0.00913038 -0.02379042 0.00000000 0.00000000 + hmelt 0.00000000 0.00000000 0.00000000 -0.89090241 0.32600589 0.56489652 0.00000000 0.00000000 + hnetsw -164.81223433 42.73795062 0.00000000 120.82313101 0.62392301 0.63004087 0.00000000 0.00281117 + hlwdn -336.69371959 87.02240716 0.00000000 239.28275576 4.67536186 5.71298315 0.00000000 -0.00021166 + hlwup 394.51485762 -107.98891375 0.00000000 -274.37211846 -5.48118805 -6.67262512 0.00000000 0.00001223 + hlatvap 85.44674862 -12.15707541 0.00000000 -73.10714372 -0.05853463 -0.12425713 0.00000000 -0.00026227 + hlatfus 0.82811998 -0.28266769 0.00000000 -0.36550691 -0.04546068 -0.13447928 0.00000000 0.00000542 + hiroff 0.00000000 0.00024877 0.00000000 -0.05026752 0.00000000 0.00000000 0.05001961 0.00000087 + hsen 20.86523625 -9.39739411 0.00000000 -11.46215568 -0.03417290 0.02776099 0.00000000 -0.00072545 + *SUM* 0.14900854 -0.06544441 0.00000000 -0.10928714 -0.00319588 -0.01947042 0.05001961 0.00163030 + +and the annual **water** budget (kg/m²/s, scaled by 1e6): + +.. code-block:: none + + (med_diag_print_summary) NET WATER BUDGET (kg/m2s*1e6): period = annual: date = 11060101 0 + atm lnd rof ocn ice nh ice sh glc *SUM* + wfreeze 0.00000000 0.00000000 0.00000000 -0.08728167 0.02420702 0.06307465 0.00000000 0.00000000 + wmelt 0.00000000 0.00000000 0.00000000 0.55264305 -0.16932198 -0.38332107 0.00000000 -0.00000000 + wrain -31.68859545 7.10225360 0.00000000 24.51612745 0.04130430 0.02895740 0.00000000 0.00004730 + wsnow -2.48163015 0.84707129 0.00000000 1.09531588 0.13623219 0.40299455 0.00000000 -0.00001623 + wevap 34.15015152 -4.85457580 0.00000000 -29.23116502 -0.02073919 -0.04377629 0.00000000 -0.00010478 + weqsaltf 0.00000000 0.00000000 0.00000000 0.00072424 0.00102805 -0.00175229 0.00000000 0.00000000 + wrunoff 0.00000000 -2.84704267 0.00760572 2.83936982 0.00000000 0.00000000 0.00000000 -0.00006712 + wfrzrof 0.00000000 -0.00074550 0.00000000 0.15063685 0.00000000 0.00000000 -0.14989395 -0.00000260 + *SUM* -0.02007408 0.24696091 0.00760572 -0.16362939 0.01271039 0.06617694 -0.14989395 -0.00014344 + +The rows are the individual flux terms — net shortwave, up/down longwave, latent +and sensible heat, freeze/melt and runoff for heat; rain, snow, evaporation, +runoff, freeze/melt for water. The small ``*SUM*`` values confirm the coupled +budgets close to round-off. The area budget has the same layout. + +Enabling the budgets +==================== + +The budgets are enabled by the ``do_budgets`` attribute, and the ``budget_*`` +settings select which averaging periods (instantaneous, daily, monthly, annual +and long-term) are printed. These are configured in :ref:`running-a-case`. diff --git a/doc/source/developer/field_exchange.rst b/doc/source/developer/field_exchange.rst new file mode 100644 index 000000000..ea2888ee9 --- /dev/null +++ b/doc/source/developer/field_exchange.rst @@ -0,0 +1,218 @@ +.. _field-exchange: + +============== +Field exchange +============== + +The mediator sets up all of its field exchange during initialization: which +fields it exchanges with each component, how each field is mapped between grids, +and how mapped fields are merged. To follow that setup, start with the NUOPC +two-phase field mechanism it is built on. + +Advertise and realize +===================== + +NUOPC (the coupling layer built on ESMF) splits component initialization into a +two-phase field setup, so components do not need to know each other's internals +up front: + +* **Advertise** — each component lists the *names* of the fields it can + potentially send or receive (for example ``sea_surface_temperature``, + ``air_pressure``), without allocating any real memory or grid data behind + them. Think of it as everyone posting a menu of "here is what I can offer / + here is what I need," so the coupler (the mediator and driver) can figure out + who talks to whom. +* **Realize** — after the coupler matches producers to consumers based on those + advertised names, each component goes back and actually creates the real data + structures for the fields that got matched: allocating arrays, attaching them + to a grid or mesh, and connecting them to memory. Fields that were advertised + but never matched to a partner are simply dropped instead of wastefully + allocated. + +Splitting setup this way keeps components decoupled from each other's grids and +internals during advertise, and only pays the memory and setup cost in realize, +once a field is known to be needed. + +Within this scheme the mediator and the components advertise differently: the +mediator advertises **all possible fields** that could be exchanged with a +component — a set that is kept consistent with the ``fd_.yaml`` field +dictionary — whereas each component advertises only the fields it actually +expects to send and receive. + +.. important:: + + Because the mediator advertises this superset while each component advertises + only its own subset, NUOPC connects only the **intersection** — a field is + realized only where both the mediator and a component advertised it. This is + why the ``addfld`` calls below advertise so many fields: the mediator offers + the entire menu, and the actual component set and configuration determine + which offerings are taken up. + +The field-exchange module +========================= + +The mediator's side of all this lives in one place: the host's **field-exchange +module**, ``esmFldsExchange__mod.F90``. It declares *which* fields are +exchanged and *how* they are mapped and merged, using three routines provided by +the generic ``esmFlds.F90``: + +* ``addfld`` — **advertise** a field the mediator can receive from, or send to, a + component. +* ``addmap`` — declare **how** a source field is mapped from its grid to a + destination grid. +* ``addmrg`` — declare **how** one or more mapped source fields are merged into a + destination field. + +Understanding these three routines is the key to the phases that follow: the +phases simply execute, every coupling step, the maps and merges declared here. + +When ``esmFldsExchange`` runs +============================= + +``esmFldsExchange_`` is called twice during initialization, selected by a +``phase`` argument: + +* ``phase == 'advertise'`` — makes all the ``addfld`` calls. This is the + mediator's contribution to the NUOPC **advertise** phase above: it advertises + the full superset of exchangeable fields, and NUOPC then determines which are + connected and realized. +* ``phase == 'initialize'`` — runs after the connected fields have been realized, + and makes the ``addmap`` and ``addmrg`` calls that set up mapping and merging + for the fields that actually got connected. + +Components are referred to by a **component index** — one of ``compatm``, +``compice``, ``compglc``, ``complnd``, ``compocn``, ``comprof`` or ``compwav``. + +``addfld`` — advertise fields +============================= + +The mediator advertises every field it *might* receive from or send to a +component. There are two directions: + +.. code-block:: fortran + + call addfld_from(comp_index, 'field_name') ! a field the mediator receives FROM comp_index + call addfld_to(comp_index, 'field_name') ! a field the mediator sends TO comp_index + +``addfld_from`` populates the "from" field list for a component +(``fldListFr(comp_index)``) and ``addfld_to`` the "to" field list +(``fldListTo(comp_index)``). Advertising a field does not guarantee it is used — +only advertised fields that are connected on both sides are realized. + +``addmap`` — map a source field to a destination +================================================ + +For each source field that must reach another component's grid, ``addmap`` +declares the mapping: + +.. code-block:: fortran + + call addmap_from(comp_index_src, 'field_name', comp_index_dst, maptype, mapnorm, mapfile) + +A single source field may be mapped to several destinations, so it may appear in +more than one ``addmap`` call. The arguments are: + +**maptype** — the interpolation method (from ESMF). Common values include: + +* ``mapbilnr`` — bilinear +* ``mapconsf`` — first-order conservative, conservative-fraction normalization +* ``mapconsd`` — first-order conservative, destination-area normalization +* ``mappatch`` — patch recovery +* ``mapfcopy`` — redistribution (identical source and destination grids) +* ``mapnstod`` — nearest source-to-destination (and the + ``mapnstod_consd`` / ``mapnstod_consf`` combinations) + +**mapnorm** — the mapping normalization: + +* ``unset`` — no normalization (only valid with ``mapfcopy``) +* ``none`` — no normalization (when ``maptype`` is not ``mapfcopy``) +* ``one`` — normalize by one +* ``lfrin`` / ``ifrac`` / ``ofrac`` — normalize by the land, ice or ocean + fraction in ``FBFrac(...)`` (used when mapping land, ice or ocean fields to the + atmosphere) +* ``custom`` — normalization is done by hand in the ``prep`` phase (used, for + example, mapping glc to land) + +.. note:: + + When a fractional ``mapnorm`` is used, the field is **scaled by the relevant + fraction before mapping and unscaled by the mapped fraction after mapping**. + This fraction weighting is what keeps ocean/ice-to-atmosphere mappings + accurate where surface fractions vary — see :ref:`surface fractions + `. + +**mapfile** — whether to build the routehandle online or read weights from a file: + +* ``unset`` — generate the routehandle at run time +* a pathname (supplied as a driver attribute) — read the mapping weights from + that file + +Example:: + + call addmap_from(compice, 'Si_snowh', compatm, mapconsf, 'ifrac', 'unset') + +This maps the sea-ice field ``Si_snowh`` conservatively to the atmosphere, using +ice-fraction normalization, with the routehandle generated at run time. + +``addmrg`` — merge mapped fields into a destination +=================================================== + +Once source fields have been mapped onto a destination grid, ``addmrg`` declares +how they combine into the field the mediator exports to that component: + +.. code-block:: fortran + + call addmrg_to(comp_index_dst, 'dst_field_name', & + mrg_from=comp_index_src, mrg_fld='src_field_name', & + mrg_type='...', mrg_fracname='...') + +The arguments are: + +* ``mrg_from`` — the source component index. +* ``mrg_fld`` — the source field name in that component's mapped field bundle. +* ``mrg_type`` — how the source contributes: + + * ``copy`` — copy the mapped source field into the destination. + * ``copy_with_weights`` — copy, weighting by the source fraction on the + destination grid (requires ``mrg_fracname``). + * ``sum`` — cumulative sum of the mapped source fields. + * ``sum_with_weights`` — cumulative sum, each source weighted by its fraction + (requires ``mrg_fracname``). + * ``merge`` — fraction-weighted merge of several sources into one destination. + +* ``mrg_fracname`` — the fraction field in ``FBFrac(comp_index_dst)`` used as the + weight for the weighted/merge types. + +**Multi-source merges** are expressed as **one ``addmrg_to`` call per source**, +all targeting the same destination field. For example, the merged surface +temperature the mediator sends to the atmosphere is built from land, sea ice and +ocean:: + + call addmrg_to(compatm, 'Sx_t', mrg_from=complnd, mrg_fld='Sl_t', & + mrg_type='merge', mrg_fracname=mrg_fracname_lnd2atm_state) + call addmrg_to(compatm, 'Sx_t', mrg_from=compice, mrg_fld='Si_t', & + mrg_type='merge', mrg_fracname='ifrac') + call addmrg_to(compatm, 'Sx_t', mrg_from=compocn, mrg_fld='So_t', & + mrg_type='merge', mrg_fracname='ofrac') + +Each call adds one fraction-weighted contribution; together they form the +``Sx_t`` field over each atmosphere cell (the land/ice/ocean fractions sum to +one). + +How it fits together +==================== + +For a single coupled field, the three routines line up like this: + +#. ``addfld_from`` / ``addfld_to`` advertise the field on the source and + destination sides. +#. ``addmap_from`` declares how the source field reaches each destination grid — + these become the routehandles (``RH``) and populate the off-diagonal + ``FBImp(n,k)`` bundles. +#. ``addmrg_to`` declares how the mapped contributions combine into ``FBExp`` for + the destination component. + +The run phases (next) do no field-selection logic of their own — they walk these +declarations and execute the maps and merges. See :ref:`code-organization` for +the ``FBImp`` / ``FBExp`` / ``FBFrac`` structures these routines populate, and +:ref:`fields` for the field-name convention. diff --git a/doc/source/developer/fractions.rst b/doc/source/developer/fractions.rst new file mode 100644 index 000000000..16cf8542a --- /dev/null +++ b/doc/source/developer/fractions.rst @@ -0,0 +1,154 @@ +.. _surface-fractions: + +================= +Surface fractions +================= + +Because a single grid cell can be part land, part ocean and part sea ice, the +mediator tracks the **fraction** of each surface type on every component grid. +These fractions are the weights in essentially every :ref:`merge ` and +in the fractional normalization of many :ref:`mappings `, so getting +them right is central to conservation (:ref:`concepts`). The machinery lives in +``med_fraction_mod.F90``, and the fractions are stored per grid in the +internal-state fraction bundles ``FBfrac(comp)``. + +Key assumptions +=============== + +The fraction bookkeeping rests on a few assumptions that are important to keep in +mind: + +* **The ocean and sea ice must share the same grid** (and the same mask). Ice and + ocean fractions are therefore computed together and mapped together. +* **The ocean grid has no partial land.** Every ocean gridcell is either entirely + ocean or entirely land — the ocean mask is binary (0 or 1) per cell, never a + fractional value. +* **The land fraction is derived from the ocean mask.** The fraction of land on + the land grid is obtained by mapping the ocean mask onto the land grid (land is + the complement of the mapped ocean fraction), rather than being prescribed + independently. +* **The cross-grid maps have specific forms.** The ocean-to-atmosphere and + ice-to-atmosphere maps are *masked* maps, while the land-to-atmosphere map is a + *global* map. + +The fraction fields +=================== + +The fractions for each component are held in an ESMF ``FieldBundle``. These are +the internal-state array ``FBfrac(:)`` — one bundle per component, allocated to +``ncomps`` in ``med_internalstate_mod.F90`` ("fraction data for various +components, on their grid"). They are initialized and subsequently updated in +``med_fraction_mod.F90``. Each such fraction bundle can carry the following +fields (not all on every grid): + +* ``lfrac`` — fraction of land on the grid +* ``ofrac`` — fraction of ocean on the grid +* ``ifrac`` — fraction of sea ice on the grid +* ``ifrad`` / ``ofrad`` — the ice / ocean fractions as of the last radiation + timestep + +The set carried on each grid is fixed per component, for example: + +.. code-block:: none + + atm : ifrac ofrac lfrac aofrac + ocn : ifrac ofrac ifrad ofrad + ice : ifrac ofrac + lnd : lfrac + glc : gfrac lfrac + rof : rfrac lfrac + +On the atmosphere grid the surface fractions are self-consistent and sum to one: + +.. math:: + + f_{lfrac} + f_{ofrac} + f_{ifrac} \approx 1 + +How the fractions are set +========================= + +Relative versus total fractions +------------------------------- + +Component fractions delivered at run time are **relative** fractions — the +fraction of the *active region* a surface covers, not the fraction of the whole +gridcell. For example, if a cell can be at most 50% ice and 50% land, an ice +model reporting "half covered" sends an ice fraction of 0.5 (a relative value +between 0 and 1). The mediator corrects these to **total** fractions so that, on +every grid, + +.. math:: + + f_{ifrac} + f_{ofrac} + f_{lfrac} = 1 + +The fractions are seeded from a few fields the surface components provide (read +from the diagonal ``FBImp`` bundles): + +* the ocean mask ``So_omask`` (``FBImp(compocn,compocn)``) +* the land-grid land fraction ``Sl_lfrin`` (``FBImp(complnd,complnd)``), itself a + map of the ocean mask onto the land grid +* the ice fraction ``Si_ifrac`` and ice mask ``Si_imask`` + (``FBImp(compice,compice)``) + +In the sequences below, ``mapX2Y`` denotes a mapping from the ``X`` grid to the +``Y`` grid — for example ``mapo2a`` (ocean → atmosphere), ``mapl2a`` (land → +atmosphere), ``mapa2l`` (atmosphere → land), ``mapi2o`` (ice → ocean), ``mapi2a`` +(ice → atmosphere), ``mapl2r`` (land → river) and ``mapl2g`` (land → land-ice). +These fraction and mask maps are all **first-order conservative** maps. + +Initialization (``med_fraction_init``) +-------------------------------------- + +At initialization the ice fraction is taken to be zero everywhere; the ocean and +land-grid fractions are taken from the components, mapped to the atmosphere, and +the remaining fractions derived: + +.. code-block:: none + + fractions_*(ifrac) = 0 ! ice fraction initially zero + fractions_o(ofrac) = So_omask ! ocean mask from the ocean + fractions_l(lfrac) = Sl_lfrin ! land fraction (ocean-mask complement + ! on the land grid) + + fractions_a(ofrac) = mapo2a(fractions_o(ofrac)) ! ocean fraction -> atmosphere grid + fractions_a(lfrac) = 1 - fractions_a(ofrac) ! land is the complement + ! (truncated to 0 below 0.001) + + fractions_r(lfrac) = mapl2r(fractions_l(lfrac)) ! land fraction -> river grid + fractions_g(lfrac) = mapl2g(fractions_l(lfrac)) ! land fraction -> land-ice grid + +The land-fraction truncation (values below 0.001 set to zero) helps preserve +non-land gridcells as fully non-land. + +Run time (``med_fraction_set``) +------------------------------- + +Every coupling step the time-evolving fractions are recomputed from the current +ice fraction — this is what the ``post_ice`` phase triggers when it "updates the +surface fractions": + +.. code-block:: none + + fractions_i(ifrac) = Si_ifrac ! ice fraction from ice model + fractions_i(ofrac) = 1 - fractions_i(ifrac) ! ocean = complement of ice + ! (relative -> total corrected) + fractions_o(ifrac) = mapi2o(fractions_i(ifrac)) ! -> ocean grid + fractions_o(ofrac) = mapi2o(fractions_i(ofrac)) + fractions_a(ifrac) = mapi2a(fractions_i(ifrac)) ! -> atmosphere grid + fractions_a(ofrac) = mapi2a(fractions_i(ofrac)) + +Because the land fraction is static, only the ice and ocean fractions are +updated at run time. + +Where the fractions are used +============================ + +The fractions feed the two operations of the previous pages: + +* **Merges.** A merge *to* the atmosphere is weighted by ``fractions_a`` (``lfrac``, + ``ofrac``, ``ifrac``); a merge *to* the ocean is weighted by ``fractions_o`` + (``ofrac``, ``ifrac``) normalized to one. +* **Mapping normalization.** Fraction-weighted mappings scale by a fraction + before mapping and unscale by the mapped fraction after (:ref:`mapping`) — for + example ``mapo2a`` scales by ``fractions_o(ofrac)`` and divides by + ``fractions_a(ofrac)``; ``mapi2a`` uses ``ifrac``; ``mapl2a`` uses ``lfrac``. diff --git a/doc/source/developer/history.rst b/doc/source/developer/history.rst new file mode 100644 index 000000000..357d9d023 --- /dev/null +++ b/doc/source/developer/history.rst @@ -0,0 +1,164 @@ +.. _history: + +============== +History output +============== + +``med_phases_history_mod.F90`` writes the mediator's history files — snapshots of +the coupled fields, on their grids, at a requested frequency. It is the coupler's +own diagnostic output, separate from each component's history. + +Instantaneous and time-averaged streams +======================================== + +History output comes in two kinds, maintained as separate streams: + +* **instantaneous** — the field values at the output time; and +* **time-averaged** — the fields averaged over the output interval (the mediator + accumulates them between writes and divides at write time). + +Each stream has its own output frequency, set by a ``freq_option`` / ``freq_n`` +pair (for example "every ``nmonths``", "every ``nsteps``"). These come from the +``history_option_*`` / ``history_n_*`` settings (see :ref:`running-a-case`). + +The write phases +================ + +The run sequence calls the appropriate write phase: + +* ``med_phases_history_write`` — instantaneous output of all variables; +* ``med_phases_history_write_comp`` — instantaneous, averaged and auxiliary + output for a component (``..._comp_inst`` and ``..._comp_avg`` do the actual + writes); +* ``med_phases_history_write_med`` — the mediator-computed fields (the + atmosphere/ocean fluxes and ocean albedos); +* ``med_phases_history_write_data2glc`` — the (normally yearly) average of the + fields sent to the land-ice component. + +``med_phases_history_init`` sets up the streams and their alarms. + +File names +========== + +The mediator history files are named after the case with a coupler (``.cpl``) +tag and the output timestamp, ending in ``.nc`` — for example +``.cpl[.].h*..nc``. + +Auxiliary history (``histaux``) files +===================================== + +A ``histaux`` stream records, at the frequency they are exchanged, the fields +that a source component sends into the mediator. For example, +``histaux_atm2med_*`` saves the atmosphere-to-mediator fields each coupling +period. + +Because these files capture the exact forcing a component provided, they can be +**replayed later**: a CDEPS data model reads them back as an input stream and +uses them to drive a *target* component, without rerunning the original source +component. This is a cheap way to **spin up** a slow component. + +A common example is the CESM/CTSM land biogeochemistry, which needs many years of +atmospheric forcing to equilibrate. Rather than rerunning CAM each time, one run +with an active CAM saves the atmosphere-to-mediator fields as ``histaux`` files. +CDEPS/DATM (a data atmosphere) then repeatedly cycles over those saved years — for +instance a 20-year block — reading the ``histaux`` files as its forcing to drive +CTSM on its own. + +Auxiliary history files carry an ``.hx.`` tag in their name. For example, the +land-ice files written by ``med_phases_history_write_data2glc`` name the source +of the mapping: ``.cpl.hx.lnd2glc.