# Further Develop icon-liskov

###### tags: `functional cycle 13`

Developers: Samuel
Support: Matthias, (Christoph)
Appetite: Full Cycle

## Background

In cycle 12, icon-liskov, a DSL preprocessor was developed. It is currently albe to support ~70% of the icon-exclaim stencils. Hence, the goal of this project is to extend support of DSL preprocesor to the complete icon-exclaim codebase and make icon-liskov used in the integration.

This enables running CPU and plain OpenACC with the dsl branch.

## Known Steps

* Elicit requirements for the 30% missing stencils
* Come up with a plan to implement those requirements
* Make the icon-exclaim integration use icon-liskov
* Might be not so easy because icon-liskov generates new `*.f90` files. The icon build system itself needs to be instruted to disregard the old files but compile the new ones
* Maybe pp_ser can be used as a reference

## Non Goals

Do not try to solve stencil fusion

## Integration Code Targets

total: 111/111

### mo_nh_diffusion

17/17

**Supported**

- [x] IMPORTS statement
- [x] DECLARE statement
- [x] DATA CREATE statement
- [x] temporary_fields_for_turbulance_diagnostics
- [x] calculate_nabla4
- [x] apply_nabla2_and_nabla4_to_vn
- [x] apply_nabla2_to_vn_in_lateral_boundary
- [x] calculate_nabla2_for_w
- [x] calculate_horizontal_gradients_for_turbulance
- [x] apply_nabla2_to_w
- [x] apply_nabla2_to_w_in_upper_damping_layer
- [x] calculate_nabla2_for_z
- [x] calculate_nabla2_of_theta
- [x] update_theta_and_exner
- [x] calculate_nabla2_and_smag_coefficients_for_vn
- [x] calculate_diagnostics_for_turbulance
- [x] temporary_field_for_grid_point_cold_pools_enhancement
- [x] enhance_diffusion_coefficient_for_grid_point_cold_pools
- [x] mo_nh_diffusion_stencil_15
- [x] apply_nabla2_and_nabla4_global_to_vn

### Subroutine stencils

6/6

- [x] mo_math_divrot_rot_vertex_ri_dsl - src/shr_horizontal/mo_math_divrot.f90
- [x] mo_icon_interpolation_scalar_cells2verts_scalar_ri_dsl
- [x] mo_intp_rbf_rbf_vec_interpol_vertex - src/shr_horizontal/mo_intp_rbf.f90
- [x] mo_math_gradients_grad_green_gauss_cell_dsl - src/shr_horizontal/mo_math_gradients.f90
- [x] mo_advection_traj_btraj_compute_o1_dsl - src/advection/mo_advection_traj.f90

### mo_solve_nonhydro

70/70

**Supported**

- [x] IMPORTS statement
- [x] DATA CREATE statement
- [x] mo_solve_nonhydro_02
- [x] mo_solve_nonhydro_03
- [x] mo_solve_nonhydro_04
- [x] mo_solve_nonhydro_05
- [x] mo_solve_nonhydro_09
- [x] mo_solve_nonhydro_10
- [x] mo_solve_nonhydro_11_lower
- [x] mo_solve_nonhydro_11_upper
- [x] mo_solve_nonhydro_14
- [x] mo_solve_nonhydro_15
- [x] mo_solve_nonhydro_17
- [x] mo_solve_nonhydro_18
- [x] mo_solve_nonhydro_19
- [x] mo_solve_nonhydro_25
- [x] mo_solve_nonhydro_30
- [x] mo_solve_nonhydro_31
- [x] mo_solve_nonhydro_32
- [x] mo_solve_nonhydro_33
- [x] mo_solve_nonhydro_34
- [x] mo_solve_nonhydro_35
- [x] mo_solve_nonhydro_36
- [x] mo_solve_nonhydro_39
- [x] mo_solve_nonhydro_40
- [x] mo_solve_nonhydro_41
- [x] mo_solve_nonhydro_42
- [x] mo_solve_nonhydro_43
- [x] mo_solve_nonhydro_44
- [x] mo_solve_nonhydro_48
- [x] mo_solve_nonhydro_49
- [x] mo_solve_nonhydro_50
- [x] mo_solve_nonhydro_56_63
- [x] mo_solve_nonhydro_57
- [x] mo_solve_nonhydro_58
- [x] mo_solve_nonhydro_59
- [x] mo_solve_nonhydro_60
- [x] mo_solve_nonhydro_64
- [x] mo_solve_nonhydro_65
- [x] mo_solve_nonhydro_68
- [x] mo_solve_nonhydro_67
- [x] mo_solve_nonhydro_62
- [x] mo_solve_nonhydro_61
- [x] mo_solve_nonhydro_55
- [x] mo_solve_nonhydro_54
- [x] mo_solve_nonhydro_27
- [x] mo_solve_nonhydro_29
- [x] mo_solve_nonhydro_28
- [x] mo_solve_nonhydro_26
- [x] mo_solve_nonhydro_24
- [x] mo_solve_nonhydro_23
- [x] mo_solve_nonhydro_56_63
- [x] mo_solve_nonhydro_01
- [x] mo_solve_nonhydro_07
- [x] mo_solve_nonhydro_08
- [x] mo_solve_nonhydro_12
- [x] mo_solve_nonhydro_47
- [x] mo_solve_nonhydro_46
- [x] mo_solve_nonhydro_06
- [x] mo_solve_nonhydro_13
- [x] mo_solve_nonhydro_22
- [x] mo_solve_nonhydro_66
- [x] mo_solve_nonhydro_16
- [x] mo_solve_nonhydro_20
- [x] mo_solve_nonhydro_52
- [x] mo_solve_nonhydro_53
- [x] mo_solve_nonhydro_4th_order_divdamp
- [x] mo_solve_nonhydro_45
- [x] mo_solve_nonhydro_45_b
- [x] mo_solve_nonhydro_37
- [x] mo_solve_nonhydro_38
- [x] mo_solve_nonhydro_21

### mo_velocity_advection

20/20

**supported**

- [x] IMPORT statement
- [x] DATA CREATE statement
- [x] mo_velocity_advection_stencil_01
- [x] mo_velocity_advection_stencil_02
- [x] mo_velocity_advection_stencil_03
- [x] mo_velocity_advection_stencil_04
- [x] mo_velocity_advection_stencil_08
- [x] mo_velocity_advection_stencil_09
- [x] mo_velocity_advection_stencil_10
- [x] mo_velocity_advection_stencil_11
- [x] mo_velocity_advection_stencil_12
- [x] mo_velocity_advection_stencil_13
- [x] mo_velocity_advection_stencil_15
- [x] wrap_run_mo_velocity_advection_stencil_16
- [x] wrap_run_mo_velocity_advection_stencil_17
- [x] wrap_run_mo_velocity_advection_stencil_19
- [x] wrap_run_mo_velocity_advection_stencil_20
- [x] mo_velocity_advection_stencil_07
- [x] wrap_run_mo_velocity_advection_stencil_05
- [x] wrap_run_mo_velocity_advection_stencil_06
- [x] wrap_run_mo_velocity_advection_stencil_18
- [x] wrap_run_mo_velocity_advection_stencil_14

## Todo
- [x] Move setting `dsl_verify` to `START DATA CREATE`.
- [x] Research requirements for the implementation of the remaining stencil patterns.
- [x] Improve UnbalancedDirectiveError message.
- [x] Allow start and end stencil directives to be repeated.
- [x] Add `noendif` kwarg to `END STENCIL DATA`
- [x] Add `ENDIF` directive.
- [x] Rename stencils according to renaming changes in icon4py.
- [x] Ensure fields generated in `START DATA CREATE` are unique.
- [x] Ensure fields generated in `IMPORTS` are unique.
- [x] Add test for stencil with wrong kwargs for noendif, noprofile.
- [x] Improve UnknownStencilError message.
- [x] Use decorator for pipeline functions.
- [x] Update README.
- [x] Add generation metadata.
- [x] Build integration.
- [x] Add last remaining stencil.
- [x] Add metadata and profile flag to configure.ac

## Notes

- In Fortran variables must be declared before they are used. So setting `dsl_verify` should be done as part of `START CREATE()`.
- Should icon-liskov generate these includes as well?

USE cudafor
USE nvtx

## Currently unsupported patterns

### Call to func in DO loop *(fixed)*

CALL get_indices_e(p_patch, jb, i_startblk, i_endblk, i_startidx, i_endidx, rl_start, rl_end)
    call nvtxStartRange("enhance_diffusion_coefficient_for_grid_point_cold_pools")
    !$ACC PARALLEL LOOP DEFAULT(NONE) GANG VECTOR COLLAPSE(2) ASYNC(1) IF( i_am_accel_node .AND. acc_on )
       DO jk = nlev-1, nlev
         DO je = i_startidx, i_endidx
            kh_smag_e(je,jk,jb) = MAX(kh_smag_e(je,jk,jb), enh_diffu_3d(iecidx(je,jb,1),jk,iecblk(je,jb,1)), &
                 enh_diffu_3d(iecidx(je,jb,2),jk,iecblk(je,jb,2)) )
         ENDDO
       ENDDO

call nvtxEndRange()
     ENDDO
#endif
#endif

     rl_start = grf_bdywidth_e+1
     rl_end   = min_rledge_int

     i_startblk = p_patch%edges%start_block(rl_start)
     i_endblk   = p_patch%edges%end_block(rl_end)

     DO jb = i_startblk,i_endblk
       CALL get_indices_e(p_patch, jb, i_startblk, i_endblk, i_startidx, i_endidx, rl_start, rl_end)

       call wrap_run_enhance_diffusion_coefficient_for_grid_point_cold_pools(kh_smag_e=kh_smag_e(:,:,1), &
           enh_diffu_3d=enh_diffu_3d(:,:,1), kh_smag_e_before=kh_smag_e_before(:,:,1), &
           vertical_lower=nlev-1, vertical_upper=nlev, horizontal_lower=i_startidx, &
           horizontal_upper=i_endidx)
     ENDDO

- call wrap_run inside `DO` loop.

*Affected stencils*:
- enhance_diffusion_coefficient_for_grid_point_cold_pools

### Not yet ported to gt4py

A number of stencils are not yet ported to gt4py.

### Unterminated DO loops (in substitution mode)

call nvtxEndRange()
#endif
	ENDDO
	!$OMP END DO
	
! <wrap run func>

**Affected stencils**

- calculate_nabla2_and_smag_coefficients_for_vn
- wrap_run_calculate_diagnostics_for_turbulance
- wrap_run_temporary_field_for_grid_point_cold_pools_enhancement
- wrap_run_enhance_diffusion_coefficient_for_grid_point_cold_pools
- wrap_run_calculate_nabla2_for_z
- wrap_run_calculate_nabla2_of_theta
- wrap_run_update_theta_and_exner
- wrap_run_mo_velocity_advection_stencil_07

### Call to subroutine

CALL rbf_vec_interpol_vertex( p_nh_prog%vn, p_patch, p_int,             &
                                    u_vert, v_vert, opt_rlend=min_rlvert_int, &
                                    opt_acc_async=.TRUE. )

*Affected stencils:*

- wrap_run_mo_intp_rbf_rbf_vec_interpol_vertex `-->` rbf_vec_interpol_vertex
- wrap_run_mo_math_divrot_rot_vertex_ri_dsl `-->` rot_vertex_ri

### DECLARE DATA: Fields with suffix

REAL(wp), DIMENSION(nproma,p_patch%nlev,p_patch%nblks_e) :: z_hydro_corr_dsl 
REAL(wp), DIMENSION(nproma,p_patch%nblks_c) :: w_1
REAL(wp), DIMENSION(nproma, p_patch%nlev, p_patch%nblks_c) :: z_th_ddz_exner_c_before_new

- does not have `_before` suffix

*Affected modules:*

- mo_solve_nonhydro.f90
- mo_velocity_advection.f90

### DECLARE DATA: Fields with custom types

REAL(wp), DIMENSION(:,:,:,:), POINTER :: ikidx_dsl

- pattern of generated code is different uses pointer.

*Affected modules:*

- mo_solve_nonhydro.f90

### Stencil is used twice

- The exact same stencils can be found more than once in the source code. Currently this is not supported by icon-liskov triggering a `RepeatedDirectiveError`.

*Affected stencils*:
- wrap_run_mo_solve_nonhydro_stencil_56_63

### Invalid parsing of output field due to parentheses *(fixed)*

When parsed to generate the before field copies the following:

theta_v_new=p_nh%prog(nnew)%theta_v(:,:,1)

is turned into

theta_v_new_before(nno:)=p_nh%prog(nnow)%theta_v(:, :, 1(nno:)

when it should be

theta_v_new_before(:,:,:) = p_nh%prog(nnow)%theta_v(:, :, :)

The output fields in the wrap_fun call are also wrongly generated as:

theta_v_new_before=theta_v_new_before(nnew)

when it should be:

theta_v_new_before=theta_v_new_before(:, :, :)

### Custom output fields

Cases:

wrap_run_mo_solve_nonhydro_stencil_06

! wrap_run_field
z_dexner_dz_c_1=z_dexner_dz_c(:,:,1,1)

! before copy field:
z_dexner_dz_c_1_before(:,:,:) = z_dexner_dz_c(:,:,:,1)

wrap_run_mo_solve_nonhydro_stencil_01
wrap_run_mo_solve_nonhydro_stencil_07
wrap_run_mo_solve_nonhydro_stencil_08
wrap_run_mo_solve_nonhydro_stencil_13

! wrap_run_field
z_rth_pr_1_before(:,:,:) = z_rth_pr(:,:,:,1)
z_rth_pr_2_before(:,:,:) = z_rth_pr(:,:,:,2)

! before copy field:
z_rth_pr_1=z_rth_pr(:,:,1,1)
z_rth_pr_2=z_rth_pr(:,:,1,2)

wrap_run_mo_solve_nonhydro_stencil_12

! wrap_run_field
z_dexner_dz_c_2=z_dexner_dz_c(:,:,1,2)

! before copy field:
z_dexner_dz_c_2_before(:,:,:) = z_dexner_dz_c(:,:,:,2)

wrap_run_mo_solve_nonhydro_stencil_21

! wrap_run field
z_hydro_corr=z_hydro_corr_dsl(:,:,1), z_hydro_corr_before=z_hydro_corr_before(:,:,1)


! before copy field
z_hydro_corr_dsl(:,:,:) = 0._wp
z_hydro_corr_before(:,:,:) = 0._wp
z_hydro_corr_before(:,nlev,:) = z_hydro_corr(:,:)

wrap_run_mo_velocity_advection_stencil_16
wrap_run_mo_velocity_advection_stencil_17
wrap_run_mo_velocity_advection_stencil_18
wrap_run_mo_velocity_advection_stencil_19
wrap_run_mo_velocity_advection_stencil_20

! wrap_run field
ddt_w_adv = p_diag%ddt_w_adv_pc(:,:,1,ntnd)

! before copy field
ddt_w_adv_pc_before(:,:,:) = p_diag%ddt_w_adv_pc(:,:,:,ntnd)

- mo_solve_nonhydro_47 (int index)
- mo_solve_nonhydro_46 (int index)

### nvtx after get_indices

    call nvtxEndRange()
#endif
    rl_start_2 = grf_bdywidth_e+1
    rl_end_2   = min_rledge

    i_startblk_2 = p_patch%edges%start_block(rl_start_2)
    i_endblk_2   = p_patch%edges%end_block(rl_end_2)

    CALL get_indices_e(p_patch, 1, i_startblk_2, i_endblk_2, &
                       i_startidx_2, i_endidx_2, rl_start_2, rl_end_2)

    call wrap_run_mo_solve_nonhydro_stencil_22( &
      ipeidx_dsl=p_nh%metrics%pg_edgeidx_dsl(:,:,1), &
      pg_exdist=p_nh%metrics%pg_exdist_dsl(:,:,1), &
      z_hydro_corr=z_hydro_corr_dsl(:,nlev,1), &
      z_gradh_exner=z_gradh_exner(:,:,1), &
      z_gradh_exner_before=z_gradh_exner_before(:,:,1), &
      vertical_lower=1, &
      vertical_upper=nlev, &
      horizontal_lower=i_startidx_2, &
      horizontal_upper=i_endidx_2 &
    )

*Affected stencils:*

- mo_solve_nonhydro_22
- mo_solve_nonhydro_66

### Two consecutive stencils

call wrap_run_mo_solve_nonhydro_stencil_37(vn=p_nh%prog(nnew)%vn(:,:,1), vt=p_nh%diag%vt(:,:,1), &
          vn_ie=p_nh%diag%vn_ie(:,:,1), z_vt_ie=z_vt_ie(:,:,1), z_kin_hor_e=z_kin_hor_e(:,:,1), vn_ie_before=vn_ie_before(:,:,1), &
          z_vt_ie_before=z_vt_ie_before(:,:,1), z_kin_hor_e_before=z_kin_hor_e_before(:,:,1), vertical_lower=1, &
          vertical_upper=1, horizontal_lower=i_startidx, horizontal_upper=i_endidx)

call wrap_run_mo_solve_nonhydro_stencil_38(vn=p_nh%prog(nnew)%vn(:,:,1), wgtfacq_e=p_nh%metrics%wgtfacq_e_dsl(:,:,1), &
          vn_ie=p_nh%diag%vn_ie(:,:,1), vn_ie_before=vn_ie_before(:,:,1), vn_ie_rel_tol=1e-10_wp, vertical_lower=nlevp1, &
          vertical_upper=nlevp1, horizontal_lower=i_startidx, horizontal_upper=i_endidx)

*Affected stencils*

- mo_solve_nonhydro_37
- mo_solve_nonhydro_38
- mo_solve_nonhydro_45
- mo_solve_nonhydro_45_b
- wrap_run_mo_velocity_advection_stencil_05
- wrap_run_mo_velocity_advection_stencil_06

### Custom variables & pattern

wrap_run_mo_velocity_advection_stencil_14

! before fields
z_w_con_c_before(:,:) = z_w_con_c(:,:)
cfl_clipping(:,:) = .FALSE.
cfl_clipping_dsl(:,:) = 0._wp
cfl_clipping_before(:,:) = 0._wp
pre_levmask(:,:,:) = .FALSE.
pre_levelmask_dsl(:,:) = 0._wp
pre_levelmask_before(:,:) = 0._wp
vcfl_dsl(:,:) = 0._wp
vcfl_before(:,:) = 0._wp

- custom variables as before fields

call nvtxEndRange()
!$ACC PARALLEL IF( i_am_accel_node .AND. acc_on ) DEFAULT(NONE) ASYNC(1)
cfl_clipping_dsl(:,:) = real(cfl_clipping(:,:), wp)
!$ACC END PARALLEL
!DSL: Have to copy through explicit loop since the following:
!pre_levelmask_dsl(:,:) = real(pre_levmask(:,1,:), wp)
!results in -1.0 instead of 1.0 entries when the original logical is .TRUE.
!$ACC PARALLEL IF( i_am_accel_node .AND. acc_on )  DEFAULT(NONE) ASYNC(1)
      DO jk = MAX(3,nrdmax_jg-2), nlev-3
        DO jc = i_startidx, i_endidx
          IF (pre_levmask(jc,1,jk)) THEN
            pre_levelmask_dsl(jc,jk) = 1._wp
          ENDIF
        ENDDO
      ENDDO
!$ACC END PARALLEL
#endif
    call wrap_run_mo_velocity_advection_stencil_14(cfl_w_limit=cfl_w_limit, dtime=dtime, &
        ddqz_z_half=p_metrics%ddqz_z_half(:,:,1), z_w_con_c=z_w_con_c(:,:), &
        cfl_clipping=cfl_clipping_dsl(:,:), pre_levelmask=pre_levelmask_dsl(:,:), &
        vcfl=vcfl_dsl(:,:), z_w_con_c_before=z_w_con_c_before(:,:), &
        cfl_clipping_before=cfl_clipping_before(:,:), &
        pre_levelmask_before=pre_levelmask_before(:,:), vcfl_before=vcfl_before(:,:), &
        vertical_lower=MAX(3,nrdmax_jg-2), vertical_upper=nlev-3, horizontal_lower=i_startidx, &
        horizontal_upper=i_endidx)

        maxvcfl_dsl = 0
!$ACC PARALLEL IF( i_am_accel_node .AND. acc_on ) PRIVATE(vcfl) DEFAULT(NONE) ASYNC(1)
      !$ACC LOOP GANG VECTOR COLLAPSE(2) REDUCTION( max:maxvcfl_dsl )
      DO jk = MAX(3,nrdmax_jg-2), nlev-3
        DO jc = i_startidx, i_endidx
          IF (pre_levelmask_dsl(jc,jk)) THEN
            levelmask_dsl(jk) = 1.0_wp
          ENDIF

          IF ( cfl_clipping(jc,jk) ) THEN
            vcfl = vcfl_dsl(jc,jk)
            maxvcfl_dsl = MAX(maxvcfl_dsl, ABS(vcfl))
          ENDIF
        ENDDO
      ENDDO
!$ACC END PARALLEL

- custom pattern in wrap function call.

mo_solve_nonhydro_21:

- Custom before fields
- Custom Assign after Fortran stencil

#ifdef __DSL_VERIFY
  !$ACC PARALLEL IF( i_am_accel_node .AND. acc_on ) DEFAULT(NONE) ASYNC(1)
  z_hydro_corr_dsl(:,:,:) = 0._wp
  z_hydro_corr_before(:,:,:) = 0._wp
  z_hydro_corr_before(:,nlev,:) = z_hydro_corr(:,:)
  !$ACC END PARALLEL
  call nvtxStartRange("mo_solve_nonhydro_stencil_21")
!$ACC PARALLEL IF( i_am_accel_node .AND. acc_on )  DEFAULT(NONE) ASYNC(1)
            !$ACC LOOP GANG VECTOR PRIVATE(z_theta1,z_theta2)
            DO je = i_startidx, i_endidx

              z_theta1 = &
                 p_nh%prog(nnow)%theta_v(icidx(je,jb,1),ikidx(1,je,nlev,jb),icblk(je,jb,1)) +  &
                 p_nh%metrics%zdiff_gradp(1,je,nlev,jb)*                                       &
                (p_nh%diag%theta_v_ic(icidx(je,jb,1),ikidx(1,je,nlev,jb),  icblk(je,jb,1)) -   &
                 p_nh%diag%theta_v_ic(icidx(je,jb,1),ikidx(1,je,nlev,jb)+1,icblk(je,jb,1))) *  &
                 p_nh%metrics%inv_ddqz_z_full(icidx(je,jb,1),ikidx(1,je,nlev,jb),icblk(je,jb,1))

              z_theta2 = &
                 p_nh%prog(nnow)%theta_v(icidx(je,jb,2),ikidx(2,je,nlev,jb),icblk(je,jb,2)) +  &
                 p_nh%metrics%zdiff_gradp(2,je,nlev,jb)*                                       &
                (p_nh%diag%theta_v_ic(icidx(je,jb,2),ikidx(2,je,nlev,jb),  icblk(je,jb,2)) -   &
                 p_nh%diag%theta_v_ic(icidx(je,jb,2),ikidx(2,je,nlev,jb)+1,icblk(je,jb,2))) *  &
                 p_nh%metrics%inv_ddqz_z_full(icidx(je,jb,2),ikidx(2,je,nlev,jb),icblk(je,jb,2))

              z_hydro_corr(je,jb) = grav_o_cpd*p_patch%edges%inv_dual_edge_length(je,jb)*    &
                (z_theta2-z_theta1)*4._wp/(z_theta1+z_theta2)**2

            ENDDO
!$ACC END PARALLEL
    call nvtxEndRange()
    !$ACC PARALLEL IF( i_am_accel_node .AND. acc_on ) DEFAULT(NONE) ASYNC(1)
    z_hydro_corr_dsl(:,nlev,:) = z_hydro_corr(:,:)  !put fortran result in "blown up" dsl version of z_hydro_corr
    !$ACC END PARALLEL
#endif
    call wrap_run_mo_solve_nonhydro_stencil_21(grav_o_cpd=grav_o_cpd, theta_v=p_nh%prog(nnow)%theta_v(:,:,1), ikidx=ikidx_dsl(:,:,:,1), &
          zdiff_gradp=p_nh%metrics%zdiff_gradp_dsl(:,:,:,1), theta_v_ic=p_nh%diag%theta_v_ic(:,:,1), inv_ddqz_z_full=p_nh%metrics%inv_ddqz_z_full(:,:,1), &
          inv_dual_edge_length=p_patch%edges%inv_dual_edge_length(:,1), z_hydro_corr=z_hydro_corr_dsl(:,:,1), z_hydro_corr_before=z_hydro_corr_before(:,:,1))