Correct regional boundary condition indices

[kaiyuan-cheng]

Description

Address the PR https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/pull/219.

Fixes # (issue)

1. ps_reg is initialized to signalling NaN
2. correct indices to avoid accessing uninitialized memory

How Has This Been Tested?

Tested on Gaea and Stellar. The fix does not change the answer.

Checklist:

Please check all whether they apply or not

• [X] My code follows the style guidelines of this project
• [X] I have performed a self-review of my own code
• [X] I have commented my code, particularly in hard-to-understand areas
• [X] I have made corresponding changes to the documentation
• [X] My changes generate no new warnings
• [X] Any dependent changes have been merged and published in downstream modules

[SamuelTrahanNOAA]

I'm testing this now on Jet and Hera to see if it fixes the boundary bugs.

[SamuelTrahanNOAA]

@kaiyuan-cheng These changes need to be on top of dev/emc for them to be in RRFS or UFS. Can you please open a PR to that branch?

[SamuelTrahanNOAA]

@kaiyuan-cheng The changes actually fail to compile in UFS, so I can't even test them. Also, this leaves me wondering if some of the bugs I'm seeing originate from differences between "main" and "dev/emc"

[SamuelTrahanNOAA]

I'm uncertain of what to do with this difference between the files (- lines are dev/emc, and + lines are your branch) but they do relate to indexing.

@@ -4281,9 +4307,9 @@
         enddo
       enddo
 !
-      ie=min(ubound(side_t0%w_BC,1),ubound(w,1))
-      je=min(ubound(side_t0%w_BC,2),ubound(w,2))
-      nz=ubound(w,3)
+      ie=min(ubound(side_t0%delp_BC,1),ubound(delp,1))
+      je=min(ubound(side_t0%delp_BC,2),ubound(delp,2))
+      nz=ubound(delp,3)
 !
       do nt=1,ntracers                                                                                                                                                                                                                      
         do k=1,nz                                                                                                                                                                                                                           
          do j=jstart,jend
          do i=i1,i2
            q(i,j,k,nt)=side_t0%q_BC(i,j,k,nt)                            &
                       +(side_t1%q_BC(i,j,k,nt)-side_t0%q_BC(i,j,k,nt))   &
                        *fraction_interval
          enddo
          enddo
        enddo
      enddo

[SamuelTrahanNOAA]

I also see this difference, and I don't know what to make of it. Again, - lines are dev/emc, and + lines are your branch:

-         call mappm(km, pe0, qp, npz, pe1,  qn1, is,ie, 0, 8, Atm%ptop)
+         call mappm(km, pe0, qp, npz, pe1,  qn1, is,ie, 0, 8)

This also has some indexing:

@@ -3765,11 +3787,9 @@
 ! If the source is from old GFS or operational GSM then the tracers will be fixed in the boundaries
 ! and may not provide a very good result
 !
+#ifndef SW_DYNAMICS
   if ( .not. data_source_fv3gfs ) then
-   if ( Atm%flagstruct%nwat .eq. 6 .or. Atm%flagstruct%nwat .eq. 7 ) then
-      if ( hailwat > 0 ) then
-        BC_side%q_BC(is:ie,j,1:npz,hailwat) = 0.
-      endif
- 
+   if ( Atm%flagstruct%nwat .eq. 6 ) then
       do k=1,npz
          do i=is,ie
             qn1(i,k) = BC_side%q_BC(i,j,k,liq_wat)

[SamuelTrahanNOAA]

I tried applying just the four lines you changed, to dev/emc, and the code fails here in fv_diagnostics.F90 with a "floating-point exception:"

Edit: It fails in debug mode (-DDEBUG=ON at ufs-weather-model level), not when optimized. The test used the gfortran compiler.

      do k=1,km
      do j=js,je
         do i=is,ie
!           qmin = min(qmin, q(i,j,k))                                                                                                                                                                                                       
!           qmax = max(qmax, q(i,j,k))                                                                                                                                                                                                       
            if( q(i,j,k) < qmin  ) then   ! <---- FLOATING-POINT EXCEPTION HERE
                qmin = q(i,j,k)
            elseif( q(i,j,k) > qmax ) then
                qmax = q(i,j,k)
            endif
          enddo
      enddo
      enddo

[bensonr]

@junwang-noaa @kaiyuan-cheng - is this PR still needed?

[bensonr]

@kaiyuan-cheng - it's been just about a year since you were last asked about the relevance of this PR. If it is no longer needed, please respond so I can close it or work to ensure it can be merged.

[kaiyuan-cheng]

@bensonr I don't think this PR is needed anymore. My proposed change was merged into the emc branch.