Enable surface moments to be prescribed by the surface model

[megandevlan]

This relates to the CLASP CPT that aims to communicate subgrid heterogeneity in the land model up to the atmosphere. The proposed approach to achieve that is by computing the surface higher order moments that CLUBB needs within the surface model itself (in a way that accounts for heterogeneous surfaces), and pass those into CLUBB where they can then be prescribed.

In order to pass the moments into CLUBB however, there will need to be modifications to the clubb_api that allow additional inputs to be drawn from within camsrfexch.F90. After discussing with others (especially @Katetc, @adamrher, and @gold2718), it doesn't seem that there's a straightforward way to achieve the end goal without touching CLUBB, so as a first step towards implementing this my aim is to develop a method that works within the CLUBB_CESM repo (we'll need this to work in CESM even if it doesn't ultimately wind up on the UWM repo).

[adamrher]

Hi Megan,

Despite my remarks yesterday, I may not completely understand what code mods are needed in clubb for this work. Are you seeking to only pass the variance's at the surface from CTSM->CLUBB? Specifically, all these surface variance's that are now being computed by clubb:

         ! Diagnose surface variances based on surface fluxes.
        call calc_surface_varnce( upwp_sfc, vpwp_sfc, wpthlp_sfc, wprtp_sfc, &      ! intent(in)
                             um(2), vm(2), Lscale_up(2), wpsclrp_sfc,        &      ! intent(in)
                             wp2_splat(1), tau_zm(1),                        &      ! intent(in)
                             wp2(1), up2(1), vp2(1),                         &      ! intent(out)
                             thlp2(1), rtp2(1), rtpthlp(1),                  &      ! intent(out)
                             sclrp2(1,1:sclr_dim),                           &      ! intent(out)
                             sclrprtp(1,1:sclr_dim),                         &      ! intent(out)
                             sclrpthlp(1,1:sclr_dim) )                              ! intent(out)

Or is there more? I seem to recall a discussion with you about passing fluxes of variances to clubb (e.g., w'rt'2), but that seems a lot less straight forward to implement.

[megandevlan]

Hi Adam,

One way to put it is that I do want to pass all the variances at the surface from CTSM->CLUBB, but I'd like that to include the third order moments (including fluxes of variances), which aren't contained in the calc_surface_varnce module. And those moments are indeed way less straight forward to implement. Meng Huang at PNNL, who's also part of the CLASP project, figured out a way to do it by modifying advance_xm_wpxp_module.F90 and advance_xp2_xpyp_module.F90 so that the surface momentum-level values are set to whatever comes in from CTSM. So for example in advance_xm_wpxp_module, I've passed in CTSM calculated values for w'2rt' as wp2rtp_sfc and use it as below:

        if ( l_upwind_wpxp_ta ) then
           ! Interpolate coef_wp2rtp_implicit and term_wp2rtp_explicit
           ! to momentum levels as coef_wp2rtp_implicit_zm and
           ! term_wp2rtp_explicit_zm, respectively.
           coef_wp2rtp_implicit_zm = zt2zm( pdf_implicit_coefs_terms%coef_wp2rtp_implicit )
            term_wp2rtp_explicit_zm = zt2zm( pdf_implicit_coefs_terms%coef_wp2thlp_implicit )

            ! clasp
 !+++ MDF 
            if (wp2rtp_sfc .ne. -9999.0_core_rknd) then 
               term_wp2rtp_explicit_zm(1) = wp2rtp_sfc
            end if
 !--- MDF

But one caveat to what I've stated above is that we've noticed that the most important surface moments to prescribe from the land model are thl', rt' and rt'thl'. Those surface variances alone seem to control most of the response we get from a heterogeneous surface, so in my current implementation everything else is now commented out and we're only prescribing those three moments at the surface. I'm beginning to think we could go the route of only including mods for those three moments, which would probably be a lot simpler to get into CLUBB since we just need to set a logical and make some slightly less intrusive changes to advance_clubb_core_module, clubb_api_module, and the clubb_intr module in CAM.

Does that help at all?

Meg

[adamrher]

It would be convenient if the bulk of the improvements come from only passing thl'2, rt'2, rt'thl'. I'm skeptical of any benefit of trying to pass a w'2x' term (flux of a flux) since it is approximated by CLUBB as:

And so CLUBB already defines these as a proportional to latent/sensible heat fluxes, which seems about as good as you can do. Is something similar done for the fluxes of variances, e.g., w'x'2? These are approximated by CLUBB as:

And so I think since you are passing the variances x'2 already, that this would impact these fluxes of variances in a way that is consistent. It seems to me that if you just pass those three moments thl'2, rt'2, rt'thl', that things sort of just work out in your favor as they got propagated through to these third order moments. It might be good to bring Meng into this discussion (git handle?).

[megandevlan]

Correct, so the flux of a flux and the fluxes of variances are handled the same way in this scenario, prescribing those surface moments at the momentum level (@meng630, that's correct, right?). It does make sense that just supplying the moments thl'2, rt'2, and rt'thl' would then mostly capture the picture.

@meng630, do you know of any other reason why we might want to prescribe more than just those three moments at the surface in CLUBB? I guess if we fix the issue of z_const = 1m vs. = ZBOT in the use_andre formulation, that would help with the consistency of level at which the moments are calculated, without needing to actually pass them from the land model too.

[adamrher]

Unless @meng630 can convince us that the third order moments should be passed from CTSM to CLUBB, then I think the next move is to ask Vince to add a logical input argument to advance_clubb_core_api and advance_clubb_core that controls whether to call calc_surface_varnce.

Can you think of any reason for this not to be in the official UWM code? Basically adding functionality for the host model to provide surface variances instead of CLUBB computing them? It seems like this is a fundamental feature required of the CLASP project, so I don't think there is any harm in reaching out to him soon. What do @swrneale and @Katetc think? My reason for advocating for this is that it would make updating the CLUBB externals in CAM much easier since we wouldn't have to re-introduce this functionality every time, which Cheryl would appreciate.

Then this whole business of what z_cnst is in CLUBB is moot, right? It wouldn't be used?

[megandevlan]

I think calc_surface_varnce still needs to be called regardless, at the very least because we're not computing these moments in the ocean and ice models right now. So what I've been doing is letting CLUBB compute the surface variances as it sees fit, then overwriting with the CTSM values only if the flag is on and if the landfrac is above some threshold (I still need to refine this so that it uses a fractional amount of the CLUBB-computed value if the whole grid cell isn't land). There's also the issue of additional moments that are computed within calc_surface_varnce that we're not passing - wp2_sfc, up2_sfc, and vp2_sfc are the main ones, as far as I know.

That said, I can't see a reason to not add this functionality to the UWM code. I can share my branch that has some kind of namelist functionality like this, though my guess is that it could probably be done more cleanly as well.

I don't know that we've really solved the z_const issue by prescribing thl', rt', and rt'thl' alone. The reason I say that is when you turn on the use_andre option in calc_surface_varnce (which is how I've been running at least), z_const is used to define zeta and still factors into w'2 when zeta<0:

      thlp2_sfc   = reduce_coef  &
                    * ( wpthlp_sfc**2 / ustar**2 ) &
                    * four * ( one - 8.3_core_rknd * zeta )**(-two_thirds)

      rtp2_sfc    = reduce_coef  &
                    * ( wprtp_sfc**2 / ustar**2 ) &
                    * four * ( one - 8.3_core_rknd * zeta )**(-two_thirds)

      rtpthlp_sfc = reduce_coef  &
                    * ( wprtp_sfc * wpthlp_sfc / ustar**2 ) &
                    * four * ( one - 8.3_core_rknd * zeta )**(-two_thirds)

      wp2_sfc     = ( ustar**2 ) &
                    * ( 1.75_core_rknd + two * (-zeta)**(two_thirds) )

But I haven't checked if this is still as big of a problem if those other three moments are at least prescribed with a consistent level being used to compute zeta. So that would be something I need to check into a bit more next week.

[adamrher]

Can't you compute wp2_sfc, up2_sfc, and vp2_sfc in CTSM? To use Erik's framework, CTSM owns the roughness length, which is used to compute zeta. And CTSM has to be computing a ustar to get the surface fluxes in the first place. Maybe this is complicated by them only being available at the patch level? But there has to be a way to make it work. I'm just not comfortable using a different zeta's for different surface properties (variance vs. sensible heat flux).

[meng630]

In terms of passing the surface third order moments to CLUBB, the initial motivation came from the paper by Machulskaya and Mironov (2018, BLM) as they stated a zero third order flux condition is inappropriate over horizontally heterogeneous surfaces, and proposed an approach to computing the third order moments in a way that explicitly takes the surface heterogeneity into account. But as Meg said, we found these third order moments actually are not impacting that much as surface variance does. I agree that just supplying thl'2, rt'2, and rt'thl' at this point

[adamrher]

Thanks for chiming in @meng630, I was curious about this. Based on the approximations to the third order moments shown above, I would expect that since w'x' and x'2 don't have zero flux b.c.'s, then neither would the third order moments, since they are constructed from w'x' and x'2. Having not looked into the code, my guess would be that w'x' and x'2 are being interpolated to thermo. levs, and used to compute the third order moments.

[meng630]

@adamrher you're right. The computed third order moments via pdf closure module are non-zero at surface, while they may be imposed to zero in the next modules aimed for advancing the prognostic equations. That's why I have to modify advance_xm_wpxp_module.F90 and advance_xp2_xpyp_module.F90 to introduce the third order moments as we want such that they would not be overwritten in the prognostic subroutines. I was attaching a little part from turbulent_adv_pdf.F90 below. Note it's setting the lower boundary array to zero and starting from the second level to compute.

`
pure subroutine xpyp_term_ta_pdf_lhs_all( coef_wpxpyp_implicit, rho_ds_zt, & ! Intent(in)
invrs_rho_ds_zm, invrs_dzm, & ! Intent(in)
l_upwind_xpyp_turbulent_adv, & ! Intent(in)
sgn_turbulent_vel, & ! Intent(in)
coef_wpxpyp_implicit_zm, & ! Intent(in)
rho_ds_zm, invrs_dzt, & ! Intent(in)
lhs_ta ) ! Intent(out)

    !---------------- Begin Code -------------------
    ! Set lower boundary array to 0

    lhs_ta(:,1) = 0.0_core_rknd

    if ( .not. l_upwind_xpyp_turbulent_adv ) then

       do k = 2, gr%nz-1

           ! Momentum superdiagonal: [ x xpyp(k+1,<t+1>) ]
           lhs_ta(1,k) = invrs_rho_ds_zm(k) * invrs_dzm(k)          &
                         * rho_ds_zt(k+1) * coef_wpxpyp_implicit(k+1) &
                         * gr%weights_zm2zt(1,k+1)

`

[adamrher]

That's interesting. Sorry if it's late on Friday, but I wanted to get this out while it's on my mind. I do recognize that this is more-or-less academic at this point, since you and Megan have already agreed that providing third order moments at the surface does not matter much. That being said, I think that the b.c. for w'x' = surface fluxes means you can't evaluate a prognostic w'x' at momentum level=1. That the turbulent advection term is zero'd out at the surface doesn't mean that w'2x' is zero at the surface. I think that for momentum level=2, the third order moment at the surface would still enter into the ta term, as it would be needed to compute the vertical gradient. For example, here's an ugly description of the terms in the w'x' eqn:

The big term in brackets is the ta term, and you can see that it is evaluating w'3/w'2 * w'x' at k-1.

Apologies if you already know this and I'm missing the bigger picture.

[meng630]

Thanks @adamrher for pointing that out though I barely understood this equations :( Actually, I'm always confused when people say the surface for a third order moment in CLUBB. If you mean the 'real' surface level (i.e., zm(1)), I agree that those third order moments are nonzero and play a role in the ta term. But if you look at the first thermo. level zt(1), they are zero at least in E3SM SCM output (no idea how come). Sorry I'm not a CLUBB expert. Hopefully that makes better sense.

[Katetc]

Apologies if I am late to the conversation. I believe @megandevlan and I decided a few months ago that the easiest thing is to prototype the changes using the CESM CLUBB repo. Once she has a good idea of exactly what changes are needed, then we should absolutely work with Vince to add these changes into CLUBB. This code is likely far more extensive than any changes we would make during an NCAR import. And it sounds like it would be scientifically useful to many CLUBB users. I would bet that the UWM guys would include the changes after a pull request and review in their code repo.

[adamrher]

@meng630 it looks to that when a derivative of a third order moment is computed, it is expressed in terms of a second order moment, e.g, w'2x' = w'3/w'2 w'x', where the second order moment w'x' is on zm levels ... and so has boundary conditions since w'2(1)=w'2_sfc and w'x'(1) = w'x'_sfc.

@megandevlan in terms of getting towards, as Kate says "a good idea of exactly what changes are needed", I think the most consistent way of providing the surface variances is also the least invasive to clubb. Basically, define all surface variances/intent out's in calc_surface_varnce in CLM so that you are using identical stability functions for all surface terms. My thinking is that way we can would only need to pass a logical through to clubb on whether to call calc_surface_varnce. And another ask would be to allow us to pass a z_cnst as well, so that over the ocean, where you aren't computing surface variances, that you at least have a z_cnst consistent with what's used in the CPL to compute ocean fluxes. What are your thought on that?

[megandevlan]

We do now have a prototype version of the CESM CLUBB repo that allows us to use the surface moments that are computed in CTSM rather than the ones computed by CLUBB by default. One note is that the 2*dt oscillations in the surface winds tend to be stronger when prescribing the momentum variances rather than sticking to only the temperature/moisture variances. That said, with the changes in ordering @adamrher has suggested in the past, this might be a moot point (it would be good to confirm that the reordering does fix the issue here as well). I'd opted not to set w'2_sfc, u'2_sfc, and v'2_sfc (and not overwrite u'w', v'w', w'thl' and w'rt') for that reason, but in theory you can/probably should prescribe them to maintain as much consistency as possible.

So if I understand your suggestion @adamrher, we compute all the surface (second order) moments in CLM using consistent definitions for stability functions, etc. The module in CLM is indeed set up to do that, so all it's a matter of is prescribing them in CLUBB; that's somewhat easy to do by overwriting whatever calc_surface_varnce computes on its own right now. So at the moment, I've modified clubb_intr.f90 as below:

      !  Surface fluxes provided by host model                                                                  
      wpthlp_sfc = cam_in%shf(i)/(cpair*rho_ds_zm(1))       ! Sensible heat flux
      wprtp_sfc  = cam_in%cflx(i,1)/rho_ds_zm(1)            ! Moisture flux  (check rho)
      upwp_sfc   = cam_in%wsx(i)/rho_ds_zm(1)               ! Surface meridional momentum flux
      vpwp_sfc   = cam_in%wsy(i)/rho_ds_zm(1)               ! Surface zonal momentum flux  
! +++ MDF
      wp2_sfc       = -9999.0_r8
      thlp2_sfc     = -9999.0_r8
      rtp2_sfc      = -9999.0_r8
      rtpthlp_sfc   = -9999.0_r8
      up2_sfc       = -9999.0_r8
     ...
      ! CLASP: use moments from CTSM 
      ! TODO: set *_sfc using weighting of moments from CTSM and CLUBB's default calculation based on landfrac
      if (clubb_ctsm_moments .and. cam_in%landfrac(i)>0.75_r8) then
          wp2_sfc       = cam_in%wp2_clubb_sfc(i)
          thlp2_sfc     = cam_in%thlp2_clubb_sfc(i)
          rtp2_sfc      = cam_in%rtp2_clubb_sfc(i)
          rtpthlp_sfc   = cam_in%thlprtp_clubb_sfc(i)
          up2_sfc       = cam_in%up2_clubb_sfc(i)
          ...
      endif
! --- MDF

It's then a matter of adding a check in advance_clubb_core_module.F90 :

        ! Diagnose surface variances based on surface fluxes.
        call calc_surface_varnce( upwp_sfc, vpwp_sfc, wpthlp_sfc, wprtp_sfc, &      ! intent(in)
                             um(2), vm(2), Lscale_up(2), wpsclrp_sfc,        &      ! intent(in)
                             wp2_splat(1), tau_zm(1),                        &      ! intent(in)
                             wp2(1), up2(1), vp2(1),                         &      ! intent(out)
                             thlp2(1), rtp2(1), rtpthlp(1),                  &      ! intent(out)
                             sclrp2(1,1:sclr_dim),                           &      ! intent(out)
                             sclrprtp(1,1:sclr_dim),                         &      ! intent(out)
                             sclrpthlp(1,1:sclr_dim) )                              ! intent(out)

        ! clasp
!+++ MDF
        if (thlp2_sfc .ne. -9999.0_core_rknd) then
           thlp2(1)   = thlp2_sfc
           rtp2(1)    = rtp2_sfc
           rtpthlp(1) = rtpthlp_sfc
           !up2(1)     = up2_sfc
           !wp2(1)     = wp2_sfc
         end if
!--- MDF

None of that's pretty, but seems to work. If I understand your suggestion right @adamrher, the alternative would be to still define the moments in say clubb_intr, but pass a logical into CLUBB so that instead of overwriting what CLUBB computes, you adjust whether or not calc_surface_varnce is even called in the first place? The only issue I see is that some of the intent(out)'s aren't computed in CLM right now(sclrp2, sclrprtp, and sclrpthlp), but I don't know quite what those do and I'd assume they can be computed in CLM anyway.

On the z_const front, it does seem reasonable enough to at least pass that from the ocean model so the height's more consistent. I need to refresh on what height they've been using to calculate the fluxes at, but I seem to recall it wasn't necessarily ZBOT (wasn't it something like 2m or 10m?, which may be less inconsistent. I'll check into that later this week, but adding z_const as an input to calc_surface_varnce shouldn't be too hard.

[adamrher]

Hi Meg,

Yeah, I like the idea of passing a logical like l_host_provides_sfc_varnce into advance_clubb_core_api. But I now see a wrinkle (and why you code it up the way you do); while there are arguments to take in the usual surface fluxes (wpxp_sfc) in advance_clubb_core_api, there are none for the surface variances (xp2_sfc). I'm wondering if initializing the relevant input arrays in clubb_intr would do the trick:

      !  Surface fluxes provided by host model                                                                  
      wpthlp_sfc = cam_in%shf(i)/(cpair*rho_ds_zm(1))       ! Sensible heat flux
      wprtp_sfc  = cam_in%cflx(i,1)/rho_ds_zm(1)            ! Moisture flux  (check rho)
      upwp_sfc   = cam_in%wsx(i)/rho_ds_zm(1)               ! Surface meridional momentum flux
      vpwp_sfc   = cam_in%wsy(i)/rho_ds_zm(1)               ! Surface zonal momentum flux  
! +++ MDF
      l_host_provides_sfc_varnce = .false.
      ! CLASP: use moments from CTSM 
      ! TODO: set *_sfc using weighting of moments from CTSM and CLUBB's default calculation based on landfrac
      if (clubb_ctsm_moments .and. cam_in%landfrac(i)>0.75_r8) then
          l_host_provides_sfc_varnce = .true.
          thlp2_in(1)     = cam_in%thlp2_clubb_sfc(i)
          rtp2_in(1)      = cam_in%rtp2_clubb_sfc(i)
          rtpthlp_in(1)   = cam_in%thlprtp_clubb_sfc(i)
          wp2_in(1)       = cam_in%wp2_clubb_sfc(i)
          up2_in(1)       = cam_in%up2_clubb_sfc(i)
          vp2_in(1)       = cam_in%vp2_clubb_sfc(i)
          sclrp2(1,:)     = ...
          sclrprtp(1,:)   = ...
          sclrpthlp(1,:)  = ...
      endif
! --- MDF

And so advance_clubb_core would just look like this:

  if (.not.l_host_provides_sfc_varnce) then

        ! Diagnose surface variances based on surface fluxes.
        call calc_surface_varnce( upwp_sfc, vpwp_sfc, wpthlp_sfc, wprtp_sfc, &      ! intent(in)
                             um(2), vm(2), Lscale_up(2), wpsclrp_sfc,        &      ! intent(in)
                             wp2_splat(1), tau_zm(1),                        &      ! intent(in)
                             wp2(1), up2(1), vp2(1),                         &      ! intent(out)
                             thlp2(1), rtp2(1), rtpthlp(1),                  &      ! intent(out)
                             sclrp2(1,1:sclr_dim),                           &      ! intent(out)
                             sclrprtp(1,1:sclr_dim),                         &      ! intent(out)
                             sclrpthlp(1,1:sclr_dim) )                              ! intent(out)

  end if

That's what I'm thinking anyway. The least certain aspect of this approach is what the heck are going on with these sclrpxp variables. These refer to things like the covariance between any non-water tracer and thlp, rtp, or with itself. I suspect they are zero, but it would be good to verify this.

[adamrher]

oh while it's on my mind. You could experiment with my new physics re-ordering if you want to. Based on our mtg today, we are still weeks away from completing this pull request. But if the cam tag you are using isn't that much older than cam6_3_020, then you should be able to just copy the relevant source mods:

https://github.com/adamrher/CAM/blob/cam6_3_020.tphysac/src/physics/cam/physpkg.F90
https://github.com/adamrher/CAM/blob/cam6_3_020.tphysac/src/physics/cam/zm_conv_intr.F90

Do you know what cam tag are you using (I usually just look at the last entry in doc/ChangeLog)?

[megandevlan]

Hi Adam,

Sorry this took a bit to get to. I agree, passing in a single logical like l_host_provides_sfc_varnce is a lot less intrusive, and this is a really helpful write up of a good way to do that - thank you!

You're also right that the sclrpxp variables are zero. As it turns out, the equations for them in surface_varnce_module.F90 depend on a variable set in clubb_intr.F90, wpsclrp_sfc, which is coded to be zero along with those sclrpxp values:

!  higher order scalar stuff, put to zero
         sclrm(k,:)          = 0._r8
         wpsclrp(k,:)        = 0._r8
         sclrp2(k,:)         = 0._r8
         sclrp3(k,:)         = 0._r8
         sclrprtp(k,:)       = 0._r8
         sclrpthlp(k,:)      = 0._r8
         wpsclrp_sfc(:)      = 0._r8
         hydromet(k,:)       = 0._r8
         wphydrometp(k,:)    = 0._r8
         wp2hmp(k,:)         = 0._r8
         rtphmp_zt(k,:)      = 0._r8
         thlphmp_zt(k,:)     = 0._r8 

So I think we can pretty safely set those values to zero when initializing them in clubb_intr. I'm just running one more test to make sure they really are zero and then I'll be able to implement the new logical you suggest above.

I like the idea of incorporating your changes to the physics ordering. I'm using cam6_3_017; do you think that's too old to just copy in the source mods? Alternatively, I doubt updating cam6_3_020 would break anything dramatically either...

[adamrher]

Ah, good. I checked too and agree all those scalar terms retain the zero's all the way down the module. So we don't have to worry about it all. I thought it was weird that sclrm wouldn't take on the tracers, since I know clubb handles all the non-water species tracers as well. Turns out there is a separate variable that takes on tracer quantities to be diffused, called edsclr_in. Those probably have zero boundary conditions because fluxes of non-water tracers at the surface are simply dumped into the lowest model layer in the (non-clubb) vertical_diffusion module. And then when CLUBB see's them they get diffused. That's a dumb way of doing it but whatever ...

You should be fine just copying my physpkg.F90 and zm_conv_intr.F90 into your source mods. I don't think there were any changes to these between cam6_3_017 and cam6_3_020.