Question about HSIMT in step3d_t.F (trunk version)

[stlaur]

Good morning,

I apologize in advance if I'm wrong about what follows. I could be misinterpreting the code.

Early on in step3d_t.F (lines 334-350), oHz is computed as:

Code: Select all

      IF (Lmpdata) THEN
        DO k=1,N(ng)
          DO j=Jstrm2,Jendp2
            DO i=Istrm2,Iendp2
              oHz(i,j,k)=1.0_r8/Hz(i,j,k)
            END DO
          END DO
        END DO
      ELSE
        DO k=1,N(ng)
          DO j=Jstr,Jend
            DO i=Istr,Iend
              oHz(i,j,k)=1.0_r8/Hz(i,j,k)
            END DO
          END DO
        END DO
      END IF

In the hypothetical case where MPDATA wasn't selected for any of the tracers, Lmpdata is false, and oHz is computed over the "regular" interval (Istr to Iend, Jstr to Jend).

However, if I look further down in step3d_t.F, it looks like HSIMT does require oHz to be computed over the same "extended stencil" as MPDATA?
For example, on lines 433-437 I see:

Code: Select all

            DO j=Jstr,Jend
              DO i=IstrU-1,Iendp2
                cff=0.125_r8*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))*   &
     &              dt(ng)
                cff1=cff*(oHz(i-1,j,k)+oHz(i,j,k))

Judging from the last line, it looks like HSIMT is using oHz up to i = Iendp2.
Similarly, on line 512, it looks like HSIMT is using oHz up to j = Jendp2.

(Again, I could be misinterpreting the code.)

[arango]

Yes, good catch Thank you. I have been looking for a bug in the HSIMT scheme on and off. We fixed the issue that it needs three ghost points a few days ago, but forgot to make changes to the range of oHz. We need have an additional switch, Lhsimt:

Code: Select all

      logical :: LapplySrc, Lhsimt, Lmpdata
      ...
      
      Lhsimt =ANY(Hadvection(:,ng)%HSIMT).and.                          &
     &        ANY(Vadvection(:,ng)%HSIMT)
      Lmpdata=ANY(Hadvection(:,ng)%MPDATA).and.                         &
     &        ANY(Vadvection(:,ng)%MPDATA)
      ...
      
!
!  Compute reciprocal thickness, 1/Hz.
!
      IF (Lmpdata.or.Lhsimt) THEN
        DO k=1,N(ng)
          DO j=Jstrm2,Jendp2
            DO i=Istrm2,Iendp2
              oHz(i,j,k)=1.0_r8/Hz(i,j,k)
            END DO
          END DO
        END DO
      ELSE
        DO k=1,N(ng)
          DO j=Jstr,Jend
            DO i=Istr,Iend
              oHz(i,j,k)=1.0_r8/Hz(i,j,k)
            END DO
          END DO
        END DO
      END IF

As you can see, the change is simple. I will test it and load it to the repository later. I am very busy at the moment.