Issue with MSKF when used with aercu_opt = 1 or 2

[sael9740]

While working with the MSKF/Morrison AA combo, I noticed that identical runs (same build, same arch, same number of MPI ranks, etc.) were not producing bit for bit results like I typically expect. This seemed to be the case for aercu_opt = 1 and aercu_opt = 2 but not for aercu_opt = 0. When digging deeper I noticed that QSATu is never initialized at the lower levels even though it is used later on in calculations for other variables:

QSATu initialization

updraft:    DO NK=K,KL-1
...
        QSATu(NK) = QSu/(1.+QSu)   ! saturated specific hum
...
            END DO updraft

QSATu used to calculate su, QSATZM and gamhat

            DO KQ=KTS,KTE
...
                   su(1,JK) = oldTU(KQ)*(1.0+0.622*QSATu(KQ)) + (G*zf_wrf(KQ))/CP !TWG updraft temperature calulation
...
                   QSATZM(1,JK) = QSATu(KQ)
...
                  gamhat(1,JK) = QSATu(KQ)*(1.+QSATu(KQ)/eps1u)   &
                 *eps1u*alatent/(R*oldTU(KQ)**2)*alatent/CP

              END DO

I am thinking that we need to add a QSATu(KQ) = ?? in the following loop to ensure it is initialized at the lower levels below NK=K: https://github.com/wrf-model/WRF/blob/21c72141142fc6c8d203d2bf79f1990e45a0aef8/phys/module_cu_mskf.F#L4717-L4739

I confirmed that adding a QSATu(KQ) = 0.0 here produces B4B results between runs but I'm wondering if the NCAR physics folks can help with what exactly this should be initialized to at the lower levels?

Steps to reproduce the behavior:

1. Use compiler and version 'nvidia' (probably does the same thing with other compilers)
2. Use namelist options 'cu_physics = 11, mp_physics = 40, aercu_opt = 1 or 2'

@weiwangncar , @dudhia

[weiwangncar]

@sael9740 I will take a look next week. Thanks for finding this.

[weiwangncar]

@sael9740 I tend to agree with you on your proposed fix to initialize qsatu in the loop you identified. Did you look at the model output with and without this fix to see if there is much impact?

[sael9740]

@weiwangncar

I tend to agree with you on your proposed fix to initialize qsatu in the loop you identified.

Is setting it to zero the correct thing to do? Or should it be initialized some other way here?

Did you look at the model output with and without this fix to see if there is much impact?

Assuming that we initialize qsat to zero it looks like it does have a significant impact in some cases. FYI - I typically do a three-way comparison of:

1. WRF built and run with -O3 (with bug fix)
2. WRF built and run with -O0 (with bug fix)
3. WRF built and run with -O3

This allows me to understand how much of an impact the changes have relative to the differences that occur due to simple floating-point error propagation, which for our purposes with GPU-WRF are extremely important but I'm explaining this here so that you can understand the plot I'm attaching:

The top row plots the fields for each run and the bottom row plots the differences between each pair. I'm seeing that the differences due to the code modification (bottom center plot) are notably larger than the differences between the -O3 and -O0 compilation (bottom left plot).

[weiwangncar]

@sael9740 I contacted the code contributor, and he suggested to initialize qsatu after line 4365 in v4.5 as

      DO K=1,KX
! initialize Qsatu
      QSATu(K) = 0.0
!
!  Saturation vapor pressure (ES) is calculated following Buck (1981)

Can you see if this works for bit-for-bit results? I think it should. Thanks!

[sael9740]

@weiwangncar - I just checked and it looks like adding QSATu(K) = 0.0 in the loop at phys/module_cu_mskf.F line 4365 also produces bit-for-bit results and has very similar effects to what we initially proposed:

This does produce different results than what we initially proposed for initializing QSATu but that makes sense after realizing all of the code I mentioned earlier is wrapped in the usl loop https://github.com/wrf-model/WRF/blob/21c72141142fc6c8d203d2bf79f1990e45a0aef8/phys/module_cu_mskf.F#L4424-L5369 and the starting index K for the updraft loop https://github.com/wrf-model/WRF/blob/21c72141142fc6c8d203d2bf79f1990e45a0aef8/phys/module_cu_mskf.F#L4741-L5026 can be different for different iterations of the usl loop. Therefore by initializing QSATu prior to the ust loop the assignments in https://github.com/wrf-model/WRF/blob/21c72141142fc6c8d203d2bf79f1990e45a0aef8/phys/module_cu_mskf.F#L5035-L5067 may be using some QSATu values that were initialized in a previous iteration of the usl loop.

My question now though is that if this is the contributor's intention for QSATu then why are we treating variables such as Aqnewlq, Aqnewic, etc. differently by initializing them here: https://github.com/wrf-model/WRF/blob/21c72141142fc6c8d203d2bf79f1990e45a0aef8/phys/module_cu_mskf.F#L4717-L4739

[weiwangncar]

@sael9740 The idea here is that the author may want to do further development of the scheme using qsatu without aerosol, and the use of qsatu may go outside the cloud layer (as you've found out). The variables aqnewlq and aqnewic are strictly in-cloud variables.

[sael9740]

@weiwangncar - That makes sense. Thanks for the explanation!

Should I create a pull request?

[weiwangncar]

@sael9740 Yes, please!