ENERGY, keep open: Convert qneg4 routine to C++

[bogensch]

The qneg4 routine in the v0 Fortran code (located at components/eam/src/physics/cam/qneg4.F90) needs to be converted to c++. This code block checks to see if moisture flux into the ground (negative value of latent heat flux) exceeds the total moisture content of the lowest model layer. Currently, v1 produces negative values of water vapor near the surface (primarily during model spin-up) and this will provide a patch for that in an energy conserving manner.

Lines 64 through 77 represent the essential code that needs to be converted. The rest of the code is for bookkeeping purposes to keep track of the number incidences of clipping (not sure if we want to convert that portion).

In v1 this block of code should be called in the macrophysics (SHOC) interface directly before shoc_main is called.

https://github.com/E3SM-Project/scream/blob/1084bdc66d2f60855de49420226580405a56673c/components/eam/src/physics/cam/qneg4.F90#L58-L77

[oksanaguba]

A couple of suggestions: IIRC the code in qneg4 creates a leak if there is not enough mass to compensate for negative CFLX in the very bottom level. There is no reason not to try to compensate for negative CFLX in the first few bottom levels, or even the whole column, to make differences smaller. Also, let's not change SHFLX in this code as it is done above, this just creates yet another energy leak.

I wrote a simple substitute for qneg4 as below for the whole column,just to have it. I did not try to experiment using only a few bottom levels, Peter B or Peter C or anyone else will know if we should use just a subset of levels and how many. My simple code is below

do i=1,ncol
  cc = cam_in%cflx(i,1)*ztodt*gravit
  if (cc < 0.0_r8) then
    cc=abs(cc)
    mm=0.0_r8
    do k=1,pver
      mm = mm + state%q(i,k,1)*state%pdel(i,k)
    enddo
    cam_in%cflx(i,1) = 0.0_r8
    if (mm < cc) then
      state%q(i,1:pver,1) = 0.0_r8 !allow leak
    else
      do k=1,pver
        qq = state%q(i,k,1)
        pp = state%pdel(i,k)
        adjust = cc *  qq*pp/mm
        state%q(i,k,1) = (qq*pp - adjust)/pp
      enddo
    endif
!print *, "OG cc, mm, new mm",cc,mm,sum(state%q(i,:,1)*state%pdel(i,:))
!print *, 'OG before, after',mm-cc, sum(state%q(i,:,1)*state%pdel(i,:))
  endif
enddo

[PeterCaldwell]

Thanks Oksana. Your code looks good to me. I think we should just borrow mass from the entire column because I don't think there's a logical cutoff for where to stop and there isn't much qv at high levels anyways, so functionally we'll just be borrowing from near the surface anyways.

One question: will setting state%q to exactly zero in the case that column moisture<cflx deficiency cause the model to crash? I could imagine 0.0_r8 getting rounded to a negative number here even though avoiding negatives is our whole goal. I'm not sure what to do instead though.

[oksanaguba]

Great question, Peter. EAM sets the lower bound for q1 to 1e-12,

    real(r8),intent(in)    :: qminc  ! minimum value of mass mixing ratio (kg/kg)
                                     ! normally 0., except water 1.E-12, for radiation.

so, I am not sure what happens to with the original qneg4 code, when it is violated (my understanding it is violated, since 0 is used instead of 1e-12). When I ran very short ne30 eam run with the new version of qneg, the column mass was much bigger than (negative) cflx, so I would assume it won't be an issue. If there is ever such big negative cflx in the code that the whole column needs to dry out, then i would say either cflx computations are wrong or the column is already unrealistically dry?

[PeterCaldwell]

Yeah, it's hard for me to believe that cflx would ever be bigger than the total mass of moisture in the atmosphere. I think I'd rather the model errored out gracefully if that happened rather than setting qv to a value that will probably cause uncaught exceptions anyways.

[tcclevenger]

Just so I am caught up..

• We are now no longer porting qneg4.F90, but instead @oksanaguba's code.
• Error out if qv<0

If that is all correct..

• What should this new function be called?
• Should it still be called directly before shoc_main()?

[whannah1]

@tcclevenger qneg4 deals with surface fluxes, so this check needs to come before whatever applies the surface fluxes to the atmosphere, which in our case is SHOC.

I think "flux" should be in the name to communicate what the routine is doing. Maybe something like check_flux_state_consistency()

[oksanaguba]

We have been talking with Conrad more about this. One question we have for @bogensch whether he would want to see the new qneg limiter change q1 at the bottom layer or not.

Example: Consider (ignore time and density factors) cflx = -10, q1(bottom)=6. Then mass=-4 needs to be redistributed above the bottom. Do we want the limiter to do this, (version 1) cflx=0, q1=0, mass=-4 redistributed above or (version 2) cflx=-6, q1=6, mass=-4 redistributed above ?

Note that the code most likely will have a temp variable for cflx as a shoc input var, so that the issues about overwriting cflx in a macmic iteration will be avoided. We just do not know how else shoc uses cflx and whether setting it to 0 will be good enough.

[bogensch]

@oksanaguba thanks for asking. I have a preference for version 2.