Scheduled Downtime
On Friday 21 April 2023 @ 5pm MT, this website will be down for maintenance and expected to return online the morning of 24 April 2023 at the latest

Possible bugs in ideal case profile initialisation

This post was from a previous version of the WRF&MPAS-A Support Forum. New replies have been disabled and if you have follow up questions related to this post, then please start a new thread from the forum home page.

Tim Raupach

New member
Hi WRF developers,

I was looking into how WRF deals with the input sounding in ideal cases, in module_initialize_ideal.F, and I found two possible bugs I'd like to query with you.

1. In subroutine get_sounding, the surface value for water vapour mixing ratio qv is read from the input profile into variable qv_surf, but never used, and instead the first-level value in qv_input(1) is repeated as the surface value (it is used three times in each profile initiation). The same thing is done in the calc_jet_sounding subroutine and in other initialization files. Is this a deliberate choice or a bug?

2. In each case, the input sounding is dealt with and grid%ph_1 is calculated. Then in at least some cases a perturbation bubble is added to kick off convection, which modifies grid%al, and grid%ph_1 is recalculated for hydrostatic balance. This recalculation differs from the first in which indices of c1h and c2h are used. In the first calculation of ph_1, c1h(kk-1) and c2h(kk-1) is used for calculation of ph_1(kk). In the second, c1h(k) and c2h(k) are used for calculation of ph_1(k). Is this correct, or should both calculations of ph_1 be the same?

Here are the relevant code snippets. Initial setting of ph_1:

    DO kk  = 2,kte
      k = kk-1
      grid%ph_1(i,kk,j) = grid%ph_1(i,kk-1,j) - (grid%dnw(kk-1))*(       &
                   ((grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))+(grid%c1h(k)*grid%mu_1(i,j)))*grid%al(i,kk-1,j)+ &
                    (grid%c1h(k)*grid%mu_1(i,j))*grid%alb(i,kk-1,j)  )

      grid%ph_2(i,kk,j) = grid%ph_1(i,kk,j)
      grid%ph0(i,kk,j) = grid%ph_1(i,kk,j) + grid%phb(i,kk,j)

Recalculation for hydrostatic balance in the quarter_ss case:

!  rebalance hydrostatically
      DO k  = 2,kte
        grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (grid%dnw(k-1))*(       &
                     ((grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))+(grid%c1h(k)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ &
                      (grid%c1h(k)*grid%mu_1(i,j))*grid%alb(i,k-1,j)  )

        grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
        grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)

Thanks very much for your help.

Good catch with the vertical indices!

The last letter of the c1h and c2h column arrays ("h"), refers to "half". We want the indices to be on the layer average coordinate locations, not the interfaces (the "f" for full levels c1f or c2f, where full levels are those that used for the geopotential, which we are vertically integrating in these loops).

Because of this half vs full designation, we want the vertical indexing for the c1h and c2h arrays to be identical to that of the other half-layer fields that are in the DO loop: dnw, al, alb. The indexing should indeed be k-1.

Would you be interested in putting a pull request together to fix this bug? I can help you through your rite of passage.
Hi Dave,

Thanks for the reply and info. I am happy to put together a pull request to fix the vertical indices.