```phi :: (real ~ " m^2 s^-2" field (x(grid), y(grid), z(half)) on        (i = 1 .. nlon + 1 by         j = 1 .. nlat + 1 by         k = 1 .. nlev        )       )

Scalar Form

phi := ( phi_s       + sum(( R_d             * T             * (q * (1 / epsilon - 1) + 1)             * ( ((1 - ( (df(log(p), z ~ "1") * shift(p, z, - 1 / 2))                       / df(p, z ~ "1")                       ) if                   2 <= k and k <= nlev                  ) \\                  0.693147180559945                 )               + ((( df(log(p), z ~ "1")                   - ((1 - ( ( df(log(p), z ~ "1")                             * shift(p, z, - 1 / 2)                             )                           / df(p, z ~ "1")                           ) if                       2 <= k and k <= nlev                      ) \\                      0.693147180559945                     )                   ) if                   2 <= k and k <= nlev                  ) \\                  - 0.693147180559945                 )               )             ), k = k + 1 .. nlev)       + ( R_d         * T         * (q * (1 / epsilon - 1) + 1)         * ((1 - ( (df(log(p), z ~ "1") * shift(p, z, - 1 / 2))                 / df(p, z ~ "1")                 ) if             2 <= k and k <= nlev            ) \\            0.693147180559945           )         )       )
Dimensioned

phi := ( phi_s       + sum(( R_d             * T             * (q * (1 / epsilon - 1 ~ "1") + 1 ~ "1")             * ( ((( 1 ~ "1"                   - ( (df(log(p), z ~ "1") * shift(p, z, - 1 / 2))                     / df(p, z ~ "1")                     )                   ) if                   2 ~ "1" <= k and k <= nlev                  ) \\                  0.693147180559945 ~ "1"                 )               + ((( df(log(p), z ~ "1")                   - ((( 1 ~ "1"                       - ( ( df(log(p), z ~ "1")                           * shift(p, z, - 1 / 2)                           )                         / df(p, z ~ "1")                         )                       ) if                       2 ~ "1" <= k and k <= nlev                      ) \\                      0.693147180559945 ~ "1"                     )                   ) if                   2 ~ "1" <= k and k <= nlev                  ) \\                  - 0.693147180559945 ~ "1"                 )               )             ), k = k + 1 ~ "1" .. nlev)       + ( R_d         * T         * (q * (1 / epsilon - 1 ~ "1") + 1 ~ "1")         * ((( 1 ~ "1"             - ( (df(log(p), z ~ "1") * shift(p, z, - 1 / 2))               / df(p, z ~ "1")               )             ) if             2 ~ "1" <= k and k <= nlev            ) \\            0.693147180559945 ~ "1"           )         )       )
Discretized

phi := ( phi_s       + sum(( R_d             * T             * (q * (1 / epsilon - 1) + 1)             * ( ((1 - ( (df_g(log(p), z) * shift_g(p, z, - 1 / 2))                       / df_g(p, z)                       ) if                   2 <= k and k <= nlev                  ) \\                  0.693147180559945                 )               + ((( df_g(log(p), z)                   - ((1 - ( ( df_g(log(p), z)                             * shift_g(p, z, - 1 / 2)                             )                           / df_g(p, z)                           ) if                       2 <= k and k <= nlev                      ) \\                      0.693147180559945                     )                   ) if                   2 <= k and k <= nlev                  ) \\                  - 0.693147180559945                 )               )             ), k = k + 1 .. nlev)       + ( R_d         * T         * (q * (1 / epsilon - 1) + 1)         * ((1 - (df_g(log(p), z) * shift_g(p, z, - 1 / 2)) / df_g(p, z) if             2 <= k and k <= nlev            ) \\            0.693147180559945           )         )       )
CSE Eliminated

phi := ((( ( varidx(phi_s, [i, j])           + ( 198.960966707927             * varidx(T, [i, j, k])             * (0.607824693422519 * varidx(q, [i, j, 1]) + 1)             )           )         - 287.04 * varidx(t_20, [i, j, 1])         ) if         k <= 1        ) \\        ( varidx(phi_s, [i, j])        - ( 287.04          * ( varidx(t_20, [i, j, k])            - ( varidx(T, [i, j, k])              * varidx(t_26, [i, j, k])              * (0.607824693422519 * varidx(q, [i, j, k]) + 1)              )            )          )        )       )```