WRF debug mode, 'Subscript upper bound' error in Thompson microphysics

[LluisFB]

Hi,

Running with a debug compiled WRF v4.5.1 version, I found an 'Subscript upper bound' error. The message I got is the following one:

(...)
WRF NUMBER OF TILES =   1
forrtl: severe (408): fort: (2): Subscript #3 of the array TMR_RACG has value 5 which is greater than the upper bound of 1

Image              PC                Routine            Line        Source
wrf.dbg.exe        000000000DB7F9CF  for_emit_diagnost     Unknown  Unknown
wrf.exe            0000000006261D66  module_mp_thompso        2474  module_mp_thompson.f90
wrf.exe            00000000061F82DE  module_mp_thompso        1257  module_mp_thompson.f90
wrf.exe            00000000050A91D2  module_microphysi        1005  module_microphysics_driver.f90

Looking into the module phys/module_mp_thompson.f90 at line #2474 (inside subroutine mp_thompson), the 3rd index of the array

tmr_racg[/ICODE] corresponds to [ICODE]idx_bg(k)[/ICODE]
[CODE]
           prg_rcg(k) = tmr_racg(idx_g1,idx_g,idx_bg(k),idx_r1,idx_r)  &
                         + tcr_gacr(idx_g1,idx_g,idx_bg(k),idx_r1,idx_r)

The tmr_racg is defined in line # 597 as:

      if (.NOT. ALLOCATED(tmr_racg)) ALLOCATE(tmr_racg(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r))

The vector idx_bg(k) is defined in line # 1870 as:

            idx_bg(k) = MAX(1,MIN(NINT(qg1d(k)/rb(k) *0.01)+1,NRHG))

At line #457 within thompson_init, one founds:

      if (PRESENT(ng)) then
         is_hail_aware = .TRUE.
         dimNRHG = NRHG
      else
         av_g(idx_bg1) = av_g_old
         bv_g(idx_bg1) = bv_g_old
         dimNRHG = NRHG1
      endif

Being at line #75:

     INTEGER, PARAMETER, PRIVATE:: NRHG = 9
     INTEGER, PARAMETER, PRIVATE:: NRHG1 = 1

The call to thompson_init is done from phys/physics_init.F, where the presence of ng, only happens with the microphysics option

THOMPSONGH[/ICODE]
[CODE]
     CASE (THOMPSONGH)
! Cycling the WRF forecast with moving nests will cause this initialization to be
! called for each nest move. This is potentially very computationally expensive.
         IF(start_of_simulation.or.restart.or.config_flags%cycling)     &
            CALL thompson_init(HGT=z_at_q,                              &
                          ORHO=inv_dens,                                &
                          NWFA2D=qnwfa2d,                               &
                          NWFA=scalar(ims,kms,jms,P_QNWFA),             &
                          NIFA=scalar(ims,kms,jms,P_QNIFA),             &
                          NBCA=scalar(ims,kms,jms,P_QNBCA),             &
                          wif_input_opt=config_flags%wif_input_opt,     &
                          FRC_URB2D=frc_urb2d,                          &
                          NG=scalar(ims,kms,jms,P_QNG),                 &
                          DX=DX, DY=DY,                                 &
                          is_start=start_of_simulation,                 &
                          IDS=ids, IDE=ide, JDS=jds, JDE=jde, KDS=kds, KDE=kde,   &
                          IMS=ims, IME=ime, JMS=jms, JME=jme, KMS=kms, KME=kme,   &
                          ITS=its, ITE=ite, JTS=jts, JTE=jte, KTS=kts, KTE=kte)

From a pure Fortran persepective, an easy fix, would be to look for the limit with dimNRHG, I think, at line # 1870:

            idx_bg(k) = MAX(1,MIN(NINT(qg1d(k)/rb(k) *0.01)+1,dimNRHG))

How ever, I am not an especialist in micro-physics neither in Thompson scheme, and I am not sure if the purposed solution has any physics sense

Thanks,

Lluís

P.S.: I am trying to do this execise to reach an error in a simulation with RRTMG. I am aware, that I might be filling the forum with too much technical issues which I think it is not the idea of the forum. If you would prefer that I post my findings directly in the comments of the WRFgit repository, please let me know and tell how I could do that.

 

[LluisFB]

Make sure, that you remove the files about density of water species that Thompson microphysics generates at the beginning of the simulations, otherwise, they will be configurated using the wrong dimension size.

$ rm qr_acr_q*

 

[LluisFB]

Sorry I rrealize there is a second place where idx_bg(k) is defined, some lines below at line #1944, #1950:

  dx_bg(k) = idx_bg1
(...)
if (.not. is_hail_aware) idx_bg(k) = idx_bg1

Where idx_bg1 is defined in line #79 as:

      INTEGER, PARAMETER :: idx_bg1 = 5 ! index for rhog when mp=8 or 28

So I would suggest, (again, from a pure Fortran perspective without knowing if it has any physical correct meaning):

  dx_bg(k) = MIN(idx_bg1, dimNRHG)
(...)
if (.not. is_hail_aware) idx_bg(k) = MIN(idx_bg1,dimNRHG)