WRF debug error: cumulus kfeta

[LluisFB]

Dear all,

I am running on a debug compiled version of WRF v4.5.1. I found an issue related to the Kain-Fritch cumulus scheme (mp_cumulus = 1). when running in debug mode I got the following error message:

(...)
forrtl: severe (408): fort: (3): Subscript #1 of the array THE0K has value 0 which is less than the lower bound of 1

Image              PC                Routine            Line        Source
wrf.dbg.exe        000000000DB8136F  for_emit_diagnost     Unknown  Unknown
wrf.exe            000000000B3D05BD  module_cu_kfeta_m        2739  module_cu_kfeta.f90
wrf.exe            000000000B370B00  module_cu_kfeta_m        1131  module_cu_kfeta.f90
(...)

If I look into the code corresponds to phys/module_cu_kfeta.f90 inside subroutine TPMIX2:

      bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)

Which corresponds in phys/module_cu_kfeta.F to:

!c******** LOOKUP TABLE VARIABLES... ****************************
!      parameter(kfnt=250,kfnp=220)
!c
!      COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp),
!     *              alu(200),rdpr,rdthk,plutop

!C***************************************************************
!c
!c***********************************************************************
!c     scaling pressure and tt table index                   
!c***********************************************************************
!c
      tp=(p-plutop)*rdpr
      qq=tp-aint(tp)
      iptb=int(tp)+1

!
!***********************************************************************
!              base and scaling factor for the                     
!***********************************************************************
!
!  scaling the and tt table index                                  
      bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)

According to the error message, iptb=0 or iptb+1 = 0, meaning that iptb = -1 and therfore int(tp) + 1 = -1

The subroutine TPMIX2 is being used by

phys/module_cu_kfcup.F[/ICODE] and [CODE]phys/module_cu_kfeta.F[/ICODE] as it is called as (e.g. line #1131):
[CODE]
call tpmix2(p0(nk1),theteu(nk1),tu(nk1),qu(nk1),qliq(nk1),        &
                     qice(nk1),qnewlq,qnewic,XLV1,XLV0)

where p: 3D air pressure, PLUTOP=5000.0, top pressure (pascals) rdpr=1./dpr: over pressure increment, DPR=(PBOT-PLUTOP)/REAL(KFNP-1), therefore the error occurs, because:

p <= plutot = 5000.
(p-plutop)*rdpr < -1.

Which would correspond to a cumulus higher than 50 hPa.

I do not know if this has too much sense. From a Fortran perspective, it could be fixed by filtering these cases and imposing iptb = 0. An specialist in cumulus scheme, maybe should find a more physically sound solution to the issue. The so-called look-up tables seem to be defined in subroutine

kf_lutab: Given a series of series of saturation equivalent potential temperatures, the temperature is calculated.

I introduce the following modification into the code (line #2730):

      iptb=int(tp)+1

      ! Lluis Fita, CIMA April 2024
      IF (iptb <= 0) iptb = 1

 

[Ming Chen]

Thank you for the detailed description of this issue. If you compile WRF without debug mode and run the same case, can the model run sucessfully?
Also, would you please upload the code that includes your fix? Thanks!