CLUBB supersaturation issue

[adamrher]

CLUBB is supposed to remove super-saturation from the solution, but it is not very good at it (see attached pdf of RH computed just after exiting the first macmic iteration of CLUBB. Green 'cntl' is out of the box cam6).

The cause of the supersaturation in CLUBB is primarily due to the use of the trapezoidal rule to produce a layer mean cloud liquid. That is, CLUBB does its pdf closure/saturation adjustment to provide a cloud liquid at the mid-points. It then calls the pdf closure a second time to get a cloud liquid at the interfaces. The midpoint and 2 surrounding interface ql's are used by the trapezoidal rule to construct a piecewise linear subgrid reconstruction, which is integrated over the layer to give a layer-mean cloud liquid. This more often than not reduces the cloud liquid from the value produced from the pdf closure at the mid-points. When converting back to T from thetal using this new ql systematically lowers T from the value it would have had, had we used the mid-point ql from the pdf closure. The cooler T's result in a supersaturated state.

Taking the mid-point ql from the pdf closure as the layer mean ql solves this problem (light blue lines in pdf "notrap"), which is invoked by turning off the trapezoidal rule via the namelists:

clubb_l_trapezoidal_rule_zm = .false.
clubb_l_trapezoidal_rule_zt = .false.

It increases low level cloud liquid just about everywhere though ... but that is in a direction more in line with observations.

AMWG diagnostics of run with the two namelists set to false vs. cntl.

[vlarson]

Nice work, Adam.

I recommend trying the following flag settings, which are intended to go together:

l_vert_avg_closure  = .false.
l_trapezoidal_rule_zt = .false.
l_trapezoidal_rule_zm = .false.
l_call_pdf_closure_twice = .false.

Settings those flags to true increases computational cost and code complexity. Years ago, we set them to true because doing so improved the marine Sc a tad. But it seems that any improvement has vanished with the latest version of the code. So setting the flags to true no longer has any benefit.

[adamrher]

Hi Vince, I'm all for not calling the pdf_closure a second time as those calls are expensive (about 10% of the total clubb cost). But there still might be value in calling the pdf_closure on the interface levels, no? Or is interpreting the needed pdf_closure variables to the interface levels good enough? I'm thinking maybe I should run another 5 year F2000climo with:

l_trapezoidal_rule_zt = .false.
l_trapezoidal_rule_zm = .false.
l_call_pdf_closure_twice = .false.

And compare it against the run with just the trapezoidal rule's off. Thoughts? As for l_vert_avg_closure, is this conditioned on l_call_pdf_closure_twice=.true.?

[vlarson]

Hi Vince, I'm all for not calling the pdf_closure a second time as those calls are expensive (about 10% of the total clubb cost). But there still might be value in calling the pdf_closure on the interface levels, no?

There could be. UWM hasn't experimented with l_call_pdf_closure_twice independently of l_trapezoidal_rule_zt,zm. They were intended to be used together.

Or is interpreting the needed pdf_closure variables to the interface levels good enough? I'm thinking maybe I should run another 5 year F2000climo with:

l_trapezoidal_rule_zt = .false.
l_trapezoidal_rule_zm = .false.
l_call_pdf_closure_twice = .false.

And compare it against the run with just the trapezoidal rule's off. Thoughts?

Sounds prudent.

As for l_vert_avg_closure, is this conditioned on l_call_pdf_closure_twice=.true.?

In CLUBB-SCM, l_vert_avg_closure is used to set l_trapezoidal_rule_zt, l_trapezoidal_rule_zm, and l_call_pdf_closure_twice all to false. In CAM, l_vert_avg_closure may not be hooked up to do anything.

[adamrher]

Running two cases, w/ and w/o l_vert_avg_closure, the results are bit for bit. After looking through the code it is clear this namelist doesn't do anything and is indeed broken.

[adamrher]

AMWG diagnostics comparing the impact of l_call_pdf_closure_twice is rather surprising. Both runs have l_trapezoidal_rule_zt/zm = .false.

They indicate that calling the pdf_closure on interface levels results in more Sc compared with just interpolating the pdf_closure at mid-point levels to the interface levels. I don't understand why this occurs, and so I'll have to take a closer look...

[adamrher]

I'm going to post random plots here to document the analysis. l_call_pdf_closure_twice doesn't have (much) of an effect on the supersaturation issue.

[adamrher]

Summary of the climatological cam6 cntl, no trap & no trap+no calling pdf twice. No trap increases Sc quite a bit, but then turning off l_call_pdf_closure_twice reduces Sc. Small changes in RCMTEND correspond to large change in CLDLIQ. Is this the whole story, or are these namelist changes impacting MG autoconversion?

[adamrher]

Attaching the time mean (5day mean) cldliq tendency after each clubb & mg2 macmic iteration, expressed as differences (l_call_pdf_closure_twice = .true. minus l_call_pdf_closure_twice = .false.). It indicates that the incr. in cloud liquid from l_call_pdf_closure_twice = .true. may be due to "juicing" up clubb turbulence, which facilitates greater Sc cloud formation. I say this b/c for macmic /= 1, there are positive Sc cldliq tendency anomalies, which reflect pdf_closure/SVP adjustment responding to how the state is modified by the prior clubb/mg2 subcycle.

[vlarson]

The important variables that must be interpolated if l_call_pdf_closure_twice=.false. (Adam, please register at https://carson.math.uwm.edu/larson-group/clubb_site/signup/ if you don't have access) amount to wp4, wpthvp, thlpthvp, and rtpthvp. I'm guessing the latter 3 are more sensitive because they are spikier near the top of Sc clouds. These xpthvp variables are on located on the zm (interface) levels and are closed based on spiky variances, which are on the same zm levels. So maybe the interpolation smooths these variables enough to diminish the cloud. Just a conjecture.

[adamrher]

This may not perfectly encapsulate the sensitivity in the 3d simulations, but I took a look at this in the BOMEX case. The cloud fraction has an unusual maxima for 32lev & l_call_pdf_closure_twice = .false., which is opposite the behavior of the 3d sims ... just above the maxima the the 2X pdf call does give slightly more cldliq, cloud_frac. The xpthvp fluxes (pink) do look to have tighter gradients in general. The 128lev sims don't exhibit the same level of sensitivity as the 32lev runs, which is a nice sanity check.

[vlarson]

This may not perfectly encapsulate the sensitivity in the 3d simulations, but I took a look at this in the BOMEX case. The cloud fraction has an unusual maxima for 32lev & l_call_pdf_closure_twice = .false., which is opposite the behavior of the 3d sims ... just above the maxima the the 2X pdf call does give slightly more cldliq, cloud_frac.

Yes, the BOMEX 1D behavior looks different than the VOCALS 3D behavior. I wonder if a 1D Sc case like DYCOMS would behave more like the VOCALS region in your 3D runs.

[adamrher]

Thanks for the suggestion. I'm still ramping up my understanding of the different iop cases. I ran DYCOMS and it seems to capture the behavior in the VOCALS regions, although hard to know for sure. I've ran 32lev and 256 lev in the cntl config (trap rule's on + 2x call pdf), notrap, and notrap+no2Xpdf. The cldliq/frac's seem to mimic the 3D model behavior, with larger values for notrap, and then smaller values for notrap+no2xpdf.

The xpthvp terms look pretty bad with notrap+no2Xpdf (dark blue), with large fluxes well above the inversion layer and descending rather slowly w/ height thereafter, compared with the 256lev solution. My initial takeaway from this plot is that calling the pdf twice does appear to improve the solution in the 32lev configuration. Although this interpretation may not be supported by other metrics, such as w'3 (In which I'm having a hard time interpreting what constitutes as "better" for this var).

There is also an interesting sensitivity at 256levs, which may indicate the 2Xpdf call is leading to odd behavior (w/ a massive spike in xpthvp at the inversions). Perhaps the 2Xpdf call is less justified for higher-res grids, since it may break any vertical coherence in the pdf across a reasonably small vertical length scale.

[vlarson]

Thanks for the suggestion. I'm still ramping up my understanding of the different iop cases. I ran DYCOMS and it seems to capture the behavior in the VOCALS regions, although hard to know for sure. I've ran 32lev and 256 lev in the cntl config (trap rule's on + 2x call pdf), notrap, and notrap+no2Xpdf. The cldliq/frac's seem to mimic the 3D model behavior, with larger values for notrap, and then smaller values for notrap+no2xpdf.

Yes, the ordering is similar, although the effect is weak in the 1D simulation. However, I don't expect perfect agreement.

The xpthvp terms look pretty bad with notrap+no2Xpdf (dark blue), with large fluxes well above the inversion layer and descending rather slowly w/ height thereafter, compared with the 256lev solution. My initial takeaway from this plot is that calling the pdf twice does appear to improve the solution in the 32lev configuration. Although this interpretation may not be supported by other metrics, such as w'3 (In which I'm having a hard time interpreting what constitutes as "better" for this var).

The differences are small, but perhaps the improvement in wpthvp is significant. With 2xpdf, wpthvp goes negative at cloud top. That is, wpthvp kills TKE and thereby inhibits cloud-top entrainment.

There is also an interesting sensitivity at 256levs, which may indicate the 2Xpdf call is leading to odd behavior (w/ a massive spike in xpthvp at the inversions). Perhaps the 2Xpdf call is less justified for higher-res grids, since it may break any vertical coherence in the pdf across a reasonably small vertical length scale.

The cloud-top spike may be realistic. Fine-scale LES often show such a spike for Sc cases. See, e.g., Fig. 6 of Larson et al. (2012).

[adamrher]

The differences are small, but perhaps the improvement in wpthvp is significant. With 2xpdf, wpthvp goes negative at cloud top. That is, wpthvp kills TKE and thereby inhibits cloud-top entrainment.

@vlarson Since the cloud top is at 1km for all 32lev cases, and only the blue curve goes negative at 1km, you are referring to the no2xpdf case? I can see how w'2 would get "killed" by a negative buoyancy production, since it's proportional to wpthvp

For the light blue curve (notrap, 2xpdf), wpthvp goes negative within the cloud layer, at about 900m ... And this doesn't inhibit cloud top entrainment? Apologies for the naive question, maybe @mikaelwitte you can jump in and set me straight here?

[vlarson]

The differences are small, but perhaps the improvement in wpthvp is significant. With 2xpdf, wpthvp goes negative at cloud top. That is, wpthvp kills TKE and thereby inhibits cloud-top entrainment.

@vlarson Since the cloud top is at 1km for all 32lev cases, and only the blue curve goes negative at 1km, you are referring to the no2xpdf case? I can see how w'2 would get "killed" by a negative buoyancy production, since it's proportional to wpthvp

For the light blue curve (notrap, 2xpdf), wpthvp goes negative within the cloud layer, at about 900m ... And this doesn't inhibit cloud top entrainment? Apologies for the naive question, maybe @mikaelwitte you can jump in and set me straight here?

I was saying that the light blue curve shows negative wpthvp near cloud top, which should inhibit entrainment, producing more LWP. This is similar to what you find in global simulations, right?

[adamrher]

Hi Vince. Thanks for the clarification. Yes, this is consistent with what I see in the global sims. The buoyancy flux goes negative for light blue (notrap), but not for green (default), indicating this may explain the incr. in Sc in the vocals regions in the 3d sims. However, I wanted to include an updated dycoms plot here, as we've discovered that this SCAM case is not reading in the constant prescribed surface fluxes for this IOP, and instead being computed from the surface models. By reading in the surface fluxes, we can more easily compare to the UCONN LES results (attaching the pdf b/c it's higher res)

temp_prof.pdf

The buoyancy flux goes negative for both light blue and green, more so for light blue. I'm fine with the notrap soln. I'm not OK with the dark blue, which does not evaluate the PDF twice. The xpthvp fields look terrible when interpolated (dark blue) compared to the hi-res reference soln. I'd conclude that the 2X pdf call gives us a soln that is more closely aligned with the hi-res soln. And so I'd propose ditching the trapezoidal rule while retaining the 2X pdf call.

[vlarson]

Hi Vince. Thanks for the clarification. Yes, this is consistent with what I see in the global sims. The buoyancy flux goes negative for light blue (notrap), but not for green (default), indicating this may explain the incr. in Sc in the vocals regions in the 3d sims. However, I wanted to include an updated dycoms plot here, as we've discovered that this SCAM case is not reading in the constant prescribed surface fluxes for this IOP, and instead being computed from the surface models. By reading in the surface fluxes, we can more easily compare to the UCONN LES results (attaching the pdf b/c it's higher res)

temp_prof.pdf

The buoyancy flux goes negative for both light blue and green, more so for light blue. I'm fine with the notrap soln. I'm not OK with the dark blue, which does not evaluate the PDF twice. The xpthvp fields look terrible when interpolated (dark blue) compared to the hi-res reference soln. I'd conclude that the 2X pdf call gives us a soln that is more closely aligned with the hi-res soln. And so I'd propose ditching the trapezoidal rule while retaining the 2X pdf call.

That proposal sounds reasonable to me.

However, you could consider doing one more check. The green line (trap on, 2xpdf on?) also looks OK. How do the green-line and light-blue-line configurations compare in a global run with your new call ordering (CLUBB in tphysac)? I.e., in your new call order, does shutting off trap while keeping on 2xpdf still help? Or does it only help when CLUBB is in tphysbc?

[adamrher]

The call ordering of clubb in tphysac (surprisingly) doesn't make a very big difference ... although to be clear, the only call ordering I'm changing is when to pass fields to the coupler. The relative ordering of all the physics routines is otherwise unchanged. We were supposed to meet internally about this yesterday, but had to reschedule to next week at last minute. My guess is that the group will want to vet these changes in a Bcompset to be certain we aren't overlooking some adverse side effect of all these changes. And so yes, I will repeat these notrap/trap, no2Xpdf/2Xpdf tests for the new call ordering.

The cause of super-saturation is still that we are computing a saturation adjustment using thetal, qt to give us a ql, but then changing that value of ql via the trapezoidal rule, while not changing the value of thetal, qt that was used to compute ql. So the light blue line (notrap, 2Xpdf) is a fair alternative to the green line (trap,2Xpdf), in that it alleviates the super-saturation issue.

[cacraigucar]

I can bring this into an upcoming misc CAM PR that I am making. @adamrher - Could you explicitly say exactly what namelist settings you want changed?

[adamrher]

Sure, clubb_l_trapezoidal_rule_zm = .false. and clubb_l_trapezoidal_rule_zt = .false. That's all!