Skip to content

Home

Jump to bottom
Jianjun Jin, NOAA EMC edited this page Jun 10, 2026 · 5 revisions

Lifecyle of qr in global hybrid-EnVar:

For global hybrid-EnVar with cwoption=3, the lifecycle of qr is:

  1. qr becomes a 3-D control variable only if it is present in cvars3d.
  2. Global background-error setup gives qr a neutral amplitude corz=1 and q-like length scales (hwll, vz) scaled by 0.5 when individual cloud species are used under cwoption=3.
  3. prewgt.f90 builds the initial dssv for qr from dsv * corz * as3d, which for this case is effectively dsv * 1 * as3d.
  4. compute_qvar3d.f90 then makes dssv(qr) state-dependent by multiplying it by 0.05 * max(qr,1e-10) above 15 hPa, and zeroing it below that pressure threshold.
  5. bkgvar.f90 later applies that final dssv(qr) directly to the qr control field during the background-error transform.
  6. Hybrid-EnVar does not appear to give qr a special ensemble/control-space treatment here; it is carried like any other active 3-D control variable, provided it exists in cvars3d.

Below is the full trace.

Note: GitHub code-search results are capped and may be incomplete. To inspect more matches in GitHub UI, use: https://github.com/search?q=repo%3ANOAA-EMC%2FGSI+qr+cwoption%3D%3D3&type=code

1) Where qr is activated

qr is not hardwired as always active. It is active if it appears in the configured 3-D control-variable list cvars3d. In several places, GSI checks this with getindex(cvars3d,'qr'). For example, the regional code explicitly does this:

  nrf3_ql  =getindex(cvars3d,'ql')
  nrf3_qi  =getindex(cvars3d,'qi')
  nrf3_qr  =getindex(cvars3d,'qr')
  nrf3_qs  =getindex(cvars3d,'qs')
  nrf3_qg  =getindex(cvars3d,'qg')

For global all-sky handling, compute_qvar3d does not hardcode qr; instead it loops through cvars3d and matches names against cloud_names_fwd. So qr participates when:

  • icloud_cv is true,
  • cwoption==3,
  • cw is not the active cloud control variable (nrf3_cw<=0 path),
  • and qr is one of the forward cloud names in cloud_names_fwd.

2) Global background-error statistics for qr under cwoption=3

In the global stats reader m_berror_stats.f90, when cwoption==1 .or. cwoption==3, and there is no total cloud-water control variable cw but there are individual cloud variables (n_clouds_fwd>0), each matching cloud variable such as qr is assigned:

  • corz(:,:,ivar)=1
  • hwll(:,:,ivar)=0.5*hwll(:,:,iq)
  • vz(:,:,ivar)=0.5*vz(:,:,iq)

where iq is the index of q. That is the key global setup for qr.

       if (n_clouds_fwd>0 .and. icw<=0) then
         do n=1,size(cvars3d)
            do ic=1,n_clouds_fwd
               if(trim(cvars3d(n))==trim(cloud_names_fwd(ic))) then
                  ivar=n
                  do k=1,nsig
                     do i=1,nlat
                        corz(i,k,ivar)=one
                     end do
                  end do
                  hwll(:,:,ivar)=0.5_r_kind*hwll(:,:,iq)
                  vz  (:,:,ivar)=0.5_r_kind*vz  (:,:,iq)
                  exit
               end if
            end do
         end do
       end if

So for global qr with cwoption=3, corz is deliberately neutralized to 1. The variance amplitude is not coming from a qr-specific climatological corz pattern; instead, the later state-dependent scaling in compute_qvar3d carries the important amplitude dependence.

3) How prewgt.f90 constructs initial dssv(qr) in the global path

In global prewgt.f90, the generic 3-D background-error amplitude is built as:

dssv(j,i,k,n)=dsv(i,k)*corz(jx,k,n)*as3d(n)

for all ordinary 3-D variables not handled specially as SF/VP/T/O3.

        else
           do k=1,nsig
              do i=1,lon2
                 dssv(j,i,k,n)=dsv(i,k)*corz(jx,k,n)*as3d(n)
              end do
           end do
        end if

Since global m_berror_stats.f90 already set corz(qr)=1 for this case, the initial qr amplitude is effectively:

dssv_init(qr) = dsv * as3d(qr)

with the vertical structure dsv itself determined from the cloud-adjusted vz field. Because vz(qr)=0.5*vz(q), qr inherits q-like but shorter vertical correlation scales. Likewise hwll(qr)=0.5*hwll(q) sets shorter horizontal scales for the later filter setup.

4) How compute_qvar3d.f90 modifies dssv(qr)

This is the step tied to your original line. compute_qvar3d is called from prewgt.f90 after the base dssv is built. The call site comment says it is a “Special case of dssv for qoption=2 and cw”, but the code path also handles individual hydrometeors when cwoption==3.

For the individual hydrometeor path:

  • if pressure is below 15 hPa: dssv(...)=0
  • else:
    • cwtmp = ges_var(i,j,k) i.e. for qr, cwtmp = ges_qr(i,j,k)
    • floor it at 1e-10
    • dn1 = 0.05 * cwtmp
    • dssv = dn1 * dssv

That is exactly the qr rule once ges_var is associated with qr.

  else
     if (cwoption/=3) return
     do n=1,size(cvars3d)
        do ic=1,n_clouds_fwd
           if(trim(cvars3d(n))==trim(cloud_names_fwd(ic))) then
              nrf3_var=n
              call gsi_bundlegetpointer (gsi_metguess_bundle(ntguessig), &
                             trim(cloud_names_fwd(ic)),ges_var,istatus)
              if (istatus/=0) return

              do k = 1,nsig
                 do j = 1,lon2
                    do i = 1,lat2
                       if (ges_prsl(i,j,k,ntguessig)<15.0_r_kind) then
                          dssv(i,j,k,nrf3_var)=zero
                       else
                          cwtmp=ges_var(i,j,k)
                          if (ges_var(i,j,k)<1.0e-10_r_kind) cwtmp=1.0e-10_r_kind
                          dn1=0.05_r_kind*cwtmp
                          dssv(i,j,k,nrf3_var)=dn1*dssv(i,j,k,nrf3_var)

So for qr specifically:

dssv_final(qr) =
    0,                              if p < 15 hPa
    (0.05 * max(qr_guess,1e-10)) *
    dsv * corz * as3d,              otherwise

and because global setup made corz(qr)=1:

dssv_final(qr) =
    (0.05 * max(qr_guess,1e-10)) * dsv * as3d(qr)

above 15 hPa.

Interpretation of dn1=0.05*cwtmp:

  • it makes the background-error standard deviation for qr proportional to the current guessed rainwater mixing ratio;
  • when qr is tiny, the factor is tiny but bounded below by the 1e-10 floor before scaling;
  • thus qr increments are strongly suppressed where no rainwater exists in the guess.

5) Where the guessed qr field comes from

compute_qvar3d gets ges_var from gsi_metguess_bundle(ntguessig) using the cloud variable name, so for qr it reads the current guess rainwater field from the metguess bundle.

More broadly, the codebase has multiple places where qr is read into bundles. For example, the FV3/GFS-side I/O accesses qr bundle pointers, which is consistent with qr being available as a metguess field when configured.

6) Where dssv(qr) is actually used

Once dssv is finalized, bkgvar.f90 applies it by multiplying the corresponding 3-D control field by dssv(j,i,k,n). This is the actual variance/preconditioning transform step.

  do n=1,cvec%n3d
     call gsi_bundlegetpointer ( cvec,cvec%r3(n)%shortname,ptr3d,istatus )
     do k=1,nsig
        do i=1,lon2
           do j=1,lat2
              ptr3d(j,i,k)  =ptr3d(j,i,k)*dssv(j,i,k,n)
           end do
        enddo
     enddo
  end do

So the lifecycle is:

  • stats setup defines corz, hwll, vz
  • prewgt turns them into dssv_init
  • compute_qvar3d turns that into state-dependent dssv_final(qr)
  • bkgvar uses dssv_final(qr) to scale the qr control variable in the transform

7) Global vs regional difference relevant to qr

For your question, the crucial distinction is:

  • Global: qr under cwoption=3 gets neutral corz=1 and q-derived scales in m_berror_stats.f90, then compute_qvar3d applies the state-dependent 0.05*qr factor.
  • Regional: there is additional special logic in prewgt_reg.f90 and m_berror_stats_reg.f90, including radar/cloud-temperature-dependent variance paths and explicit cloud variable indices. Regional can also use corz_cld and dsv_cld branches that do not exist in the same way in the global path.

So the global hybrid-EnVar qr path is simpler:

  • no special temperature-dependent corz_cld machinery;
  • baseline amplitude neutralized via corz=1;
  • state dependence comes primarily from compute_qvar3d.

8) Does hybrid-EnVar treat qr specially?

From the code I checked, not in the core variance-transform path. Hybrid-EnVar consistency checks exist generally, e.g. l_hyb_ens with TLNMC options in gsimod.F90, but I do not see a separate qr-specific hybrid weighting rule in the path we traced. qr appears to participate like any active 3-D control variable once it is in cvars3d.

There are ensemble-I/O paths that know about qr, especially for radar or cloud-aware applications, which means the ensemble side can carry qr when configured. But that is different from the static B-error / dssv lifecycle, where qr is not singled out by hybrid logic in the files examined.

Bottom-line formula for global qr, cwoption=3

For an active qr control variable in the global path in "prewgt_reg.f90" (https://github.com/NOAA-EMC/GSI/blob/develop/src/gsi/prewgt_reg.f90#L486-L558):

corz(qr) = 1
hwll(qr) = 0.5 * hwll(q)
vz(qr)   = 0.5 * vz(q)

dssv_init(qr)  = dsv(vz(qr)) * as3d(qr)
dssv_final(qr) = 0                                  if p < 15 hPa
               = 0.05 * max(qr_guess,1e-10) *
                 dsv(vz(qr)) * as3d(qr)            otherwise

and then bkgvar multiplies the qr control field by this dssv_final(qr).

If you want, I can next trace the exact call chain that makes qr appear in cloud_names_fwd and cvars3d for the global all-sky setup.

Clone this wiki locally