# Port of interpolation and metric fields to python

- Shaped by: Magdalena
- Appetite (FTEs, weeks): full cycle
- Developers: <!-- Filled in at the betting table unless someone is specifically required here -->
- cycle 18: Andreas Jocksch, Jonas
- cycle 19: Jonas

## Problem
<!-- The raw idea, a use case, or something we’ve seen that motivates us to work on this -->

ICON precomputes a lot fields a model setup: those are geometric quantities describing the grid, horizontal and vertical interpolation coeffients and weights used later during the time loop.

Icon4py stencils use these fields as input argument to `gt4py` `field_operator`s and`program`s and for verification of ICON granules in the greenline, they are read from serialized data out of an ICON run.

The goal of this project is to port their calculation into python such that a python driven model can eventually be run independently of the a Fortran model run.

## Appetite
<!-- Explain how much time we want to spend and how that constrains the solution -->
Full cycle.
Fields can be added gradually. It is not strictly necessary to tackle the entire list below in one cycle/project.

## Solution
<!-- The core elements we came up with, presented in a form that’s easy for people to immediately understand -->
This project aims at porting the so called *interpolation fields* used in dycore and diffusion. Calculation should be done with `gt4py` where possible and can be verified against serialized data.

### **t_patch** fields
See fields marked with `y`in the `icon4py` column in [this document] (https://hackmd.io/zM0GivwaTx6iXvQB73Po3Q)

### **interpolation state** (`mo_intp_data_strc.f90`>`t_int_state`)
calculated in `shr_horizontal/mo_intp_coeffs.f90`

| Field | state | dimensions | dependencies |
|:------------------------------------------------------------- |:------------------------------------------------------------------------------------------ |:------------------------------ |:------------------------------------- |
| ~~p_int_state(jg)%rbf_vec_idx_v_1, v_2~~ | *done* | [VertexDim, V2EDim] | |
| p_int_state(jg)%rbf_vec_coeff_v_1, v_2 | *started* needs linear system solver | [VertexDim, V2EDim] | edges%cartesian_center,edges%center |
| p_int_state(jg)%rbf_vec_coeff_e_dsl | | [EdgeDim, E2C2EDim] | |
| ~~p_int_state(jg)%geofac_div~~ | *done* | [CEDim] = [CellDim, C2EDim] | |
| ~~p_int_state(jg)%geofac_n2s~~, separate into centers and nbh | *done* | [CellDim], [CellDim, C2E2CDim] | |
| ~~p_int_state(jg)%geofac_grg~~ | *done* | 2x [CellDim, C2E2CODim] | |
| ~~p_int_state(jg)%geofac_rot~~ | *done* | [VertexDim, V2EDim] | |
| ~~p_int_state(jg)%geofac_grdiv~~ | not clear, *done* , | [EdgeDim, E2C2EODim] | |
| ~~p_int_state(jg)%nudgecoeff_e~~ | *done* | [EdgeDim] | |
| ~~p_int_state(jg)%c_lin_e~~ | *done* (numpy) | [EdgeDim, E2CDim] | |
| ~~c_intp=p_int_state(jg)%cells_aw_verts~~ | *done* | [VertexDim, V2CDim] | |
| ~~p_int_state(jg)%pos_on_tplane_e_x, y~~ | *done* uses projections, can those be taken from a library? | [ECDim] = [EdgeDim, E2CDim] | |
| ~~p_int_state(jg)%e_flx_avg~~ | *done* | [EdgeDim, E2C2EODim] | ptr_int%c_bln_avg, ptr_int%geofac_div |
| ~~ptr_int%c_bln_avg~~ | *done* uses rotated lat, lon from `subroutine rotate_latlon` in `mo_math_utilities.f90` | | |
| ~~p_int_state(jg)%e_bln_c_s~~ | *done* | [CellDim, C2EDim] | |

* **metric state** `mo_nonhydro_types.f90`>`t_nh_metrics`
cacluated in `atm_dyn_iconam/mo_vertical_grid.f90`

| Field | state | Dimensions | dependencies | used in |
|:----------------------- |:--------------------------------------- |:---------------------- |:--------------------------------------------------------------------------------------------------------------------------- |:---------- |
| ~~wgtfac_c~~ done | done | [CellDim, KDim] nlevp1 | metrics%z_ifc | diffusion, |
| wgtfacq_c,wgtfacq_c_dsl | started Jonas | [CellDim, KDim] | metrics%z_ifc | dycore |
| wgtfacq_e_dsl | | [EdgeDim, KDim] | | dycore |
| wgtfac_e | | [EdgeDim, KDim] | exner_expol (config param), z_maxslp (ddxn_z_full), z_maxhgtd (ddxn_z_full, dual_edge_length): calculation straight forward | dycore |
| exner_exfac | | [CellDim, KDim] | metrics%ddxn_z_full, metrics%z_ifc | dycore |
| z_mc | done | [KDim] | metrics%z_ifc | |
| exner_ref_mc | done | [CellDim, KDim] | | |
| rho_ref_mc | done | [CellDim, KDim] | | |
| rho_ref_me | done (magdalena) | [CellDim, KDim] | metrics%z_me | |
| theta_ref_mc | done | [CellDim, KDim] | | |
| theta_ref_me | done (magdalena) | [EdgeDim, KDim] | | |
| theta_ref_ic | done (same computation as theta_ref_mc) | [CellDim, KDim] | | |
| d_exner_dz_ref_ic | done (magdalena) | [CellDim, KDim] | | |
| ddqz_z_half | done (magdalena) | [CellDim, KDim] | | |
| d2dexdz2_fac1_mc | | [CellDim, KDim] | | |
| ddxn_z_full | | [EdgeDim, KDim] | | |
| ddqz_z_full_e | | [EdgeDim, KDim] | | |
| ddxt_z_full | | [EdgeDim, KDim] | | |
| inv_ddqz_z_full | done | [CellDim, KDim] | | |
| coeff1_dwdz | *started* AJ | [CellDim, KDim] | | |
| scalfac_dd3d | | [KDim] | | |
| hmask_dd3d | | [EdgeDim] | | |
| vwind_impl_wgt | | [CellDim] | | |
| vwind_expl_wgt | | [CellDim] | | |
| pg_exdist | | [EdgeDim, KDim] | | |
| pg_edgeidx_dsl | | [EdgeDim, KDim] | | |
| mask_hdiff | | [CellDim, KDim] | | diffusion |
| zd_diffcoef | | [CellDim, KDim] | | diffusion |
| bdy_halo_c | | [CellDim] | | dycore |
| mask_prog_halo_c | | [CellDim, KDim] | | |
| coeff2_dwdz | | [CellDim, KDim] | | |
| rayleigh_w | | [KDim] | | |
| zd_intcoef | | [CECDim, KDim] | | |
| zd_vertoffset | | [CECDim, KDim] | | diffusion |
| coeff_gradekin | | [ECDim] | edges%edge_cell_length, edges%inv_dual_edge_length(je,jb) | |
| vertoffset_gradp | | [ECDim, KDim] | | |
| zdiff_gradp | | [ECDim, KDim] | | |
| | | | | |
| | | | | |

It is not currently possible to return a sparse field (Field with local dimension) from a field_operator in gt4py, however it might be possible to do so with embedded gt4py. The recommendation is look at the fields with non local dimensions (most of them in metric state) first.

## Rabbit holes
<!-- Details about the solution worth calling out to avoid problems -->

* for rbf interpolation weights verification might be difficult, maybe we want to keep them serialized for serialization

* it is not currently possible to return a sparse field (Field with local dimension) from a field_operator in gt4py, however it might be possible to do so with embedded gt4py. The recommendation is look at the fields with non local dimensions (most of them in metric state) first. Discuss with gt4py team a solution to overcome this issue (embedded should be come fully available in November with the workshop in Karlsruhe)
* wait for embedded
* use numpy until something else is available.
* ...

We need to discuss this with gt4py team to get a clearer picture on how and when it will be available.

## No-gos
<!-- Anything specifically excluded from the concept: functionality or use cases we intentionally aren’t covering to fit the ## appetite or make the problem tractable -->
We don't attempt to port all fields given in a certain Fortran derived type, rather we only port fields used in `solve_nonhydro` or `diffusion`, see the corresponding python modules or `serialbox_utils.py` for which those are.