ROMS crashing when using LwSrc

[aleeson]

Hi ROMS Users,

I am working on a model that uses both LuvSrc (for rivers) and LwSrc (for wastewater treatment plants, WWTPs) in the same grid. So far, I have not gotten the model to run consistently, and I hope that someone may have ideas of what could be the issue.

The model can run successfully for over a month of simulation time, then suddenly blow up near a WWTP. The symptoms include large velocities near the WWTP. Significantly, there is always a large negative w-velocity near the source. An example of surface velocities prior to model blowup are shown in the attached image. The green diamond denotes the location of the WWTP. When I remove all WWTPs (so there are no LwSrc's), the model runs fine.

Looking forward to hearing your thoughts and suggestions. I have attached the log file for this particular model run. Please let me know if you would like to see any other files.

Thanks,
Aurora

log.txt

(211.24 KiB) Downloaded 486 times

[jcwarner]

here are a few thoughts / suggestions.
1 what are the magnitudes of mass transport for the LwSRC? i see:
GET_NGFLD_NF90 - river runoff mass transport, 2021-02-20 12:00:00.00
(Grid= 01, Rec=5, Index=2, File: rivers.nc)
(Tmin= 18674.5000 Tmax= 18681.5000) t = 18678.5000
(Min = -9.00360189E+03 Max = 1.32969952E+03)
but that is probably for all the sources. make sure LwSrc is >0.

2 set Ninfo to provide more details. the model might have actually blown up early, but the call to write out the run details does a check. even try ninfo = 1.

3 did you look at the rst file? set nrst = 60 (that is 20 minutes, just a suggestion, you have it at 24 hours now)
nrst = 20*60/20 = 60

4 why not try hsimt for all the advections.

[arango]

It seems to me that your surface flux transport is too strong and violates the CFL condition. You can compute the CFL, given your horizontal resolution and timestep size. Then, you need to lower your timestep or reduce your point source transport. It is a basic numerical modeling strategy. The LwSRC is not coded implicitly. We coded a corant number depending on vertical velocity with explicit and implicit components. However, that capability has not been released yet.

[aleeson]

Hi,

Thank you for all of the replies/suggestions! I haven't yet run a test with the additional outputs, but intend to do so soon and will attach the new files.

In the meantime, I was hoping to get your thoughts on the issue of transport. The max transport of any LwSrc is 7 m3/s, and the source that is causing issues has a max transport of 0.1 m3/s. I have set Vshape such that the vertical point sources discharge solely from the bottom sigma layer.

In contrast, the LuvSrc's generally have much larger transport values. The nearby Whidbey east river is relatively small, with a transport rate of -4 m3/s.

My follow-up question is: does it make sense that a LwSrc of 0.1 m3/s could generate the large velocities shown in the original post, despite being much smaller than the nearby LuvSrc?

Thank you again!
Aurora

Here are the surface velocities for a run without any vertical sources:

[wilkin]

Brainstorming here ...

If your Vshape causes the source to flow entirely into the bottom cell, then a transport Q is going to drive vertical flow up out of that cell at velocity w = Q/(dx*dy). If the thickness of that cell, dz, is small, such that w*dt/dz > 1, then you might have a CFL instability of the vertical flow.

That said, there is no dynamical reason that the divergence should go entirely into the vertical velocity without some also flowing out the sides horizontally (thereby reducing the vertical velocity adjustment), but I'd have to think very hard about what is happening numerically with the order of operations in which the source modifies continuity, the vertical velocity, and the changing s-coordinate because of zeta.

[aleeson]

Hi everyone,

Again, thank you for your thoughts on this matter. I just completed another round of testing and wanted to share some new results.

Per your recommendations, I checked the vertical Courant number. They were indeed very high (but strangely only high near the surface sigma layers rather than the bottom sigma layer to which the vertical source discharges). In the latest test, I equally distributed Vshape of vertical sources amongst all sigma layers rather than discharging solely to the bottom sigma layer. Despite this update, ROMS still crashes near the site of the vertical source. So far the only remedy to the issue has been to remove the problematic point source.

Here I will note that this problem has happened before at an entirely different vertical source. After removing that vertical source, ROMS ran for 1.5 months before crashing at this second vertical source. Thus, the problem seems to be with vertical sources in general, rather than with individual sources. I could simply remove the problematic sources as they appear, but would prefer to understand the nature of the problem (and fix the issue, if possible).

Some other symptoms/observations:

• High vertical Courant number

• Vertical homogenization of velocity, temperature, biogeochemistry a few hours before blowup

• Sudden jump in SSH

• Large surface velocities near point source

Another note: forcing for vertical sources varies on a monthly basis. Thus, the input transport, temperature, biogeochem, etc. on the day of blowup is the same as the 18 days prior to blowup.

Please let me know if you have any other thoughts on this matter. I've exhausted most of my ideas and thus greatly appreciate any feedback. I have attached the most recent log file (for the case with vertical sources equally distributed amongst sigma layers). I have also included depth vs. time profiles at the location of the problematic vertical source on the day of blowup.

Thanks,
Aurora

log.txt

(211.62 KiB) Downloaded 439 times

[wilkin]

I think the sudden jump in sea level and surface velocities are symptoms of the instability, and not its origin.

Looking at the density profile and AKt vertical mixing eddy diffusivity immediately before and after the vertical homogenization of the water column can you convince yourself this is mixing due to static instability? If so, that should be OK and no indicative of an immediate problem. The sudden homogenization could, however, introduce large horizontal gradients in tracers and I see you have no explicit horizontal mixing of tracers to smooth that. Are you relying on the implicit diffusion in MPDATA to control all the small horizontal scales in tracers that can feed into instabilities?

Is this strictly a problem with the physics? i.e. nothing changes whether Biology is defined, or not?

[jcwarner]

so i can offer a few suggestions again

1- set Ninof to be 1, so we can see if there is more information.
2- use HSIMT for all the tracers including salt + temp
3 you say "Another note: forcing for vertical sources varies on a monthly basis. Thus, the input transport, temperature, biogeochem, etc. on the day of blowup is the same as the 18 days prior to blowup."
That is not completely true. ROMS is interpolating between 2 input times. you need to look at the month before and the next month to see how the fields are changing in time.
I see
GET_NGFLD_NF90 - river runoff mass transport, 2021-02-19 12:00:00.00
(Grid= 01, Rec=4, Index=1, File: rivers.nc)
(Tmin= 18674.5000 Tmax= 18681.5000) t = 18677.5000
(Min = -8.97510948E+03 Max = 1.35995384E+03)
GET_NGFLD_NF90 - river runoff potential temperature, 2021-02-19 12:00:00.00
(Grid= 01, Rec=4, Index=1, File: rivers.nc)
(Tmin= 18674.5000 Tmax= 18681.5000) t = 18677.5000
(Min = 1.81208333E+00 Max = 1.00000000E+01)
GET_NGFLD_NF90 - river runoff salinity, 2021-02-19 12:00:00.00
(Grid= 01, Rec=4, Index=1, File: rivers.nc)
(Tmin= 18674.5000 Tmax= 18681.5000) t = 18677.5000
(Min = 0.00000000E+00 Max = 0.00000000E+00)
.................
GET_NGFLD_NF90 - river runoff mass transport, 2021-02-20 12:00:00.00
(Grid= 01, Rec=5, Index=2, File: rivers.nc)
(Tmin= 18674.5000 Tmax= 18681.5000) t = 18678.5000
(Min = -9.00360189E+03 Max = 1.32969952E+03)
GET_NGFLD_NF90 - river runoff potential temperature, 2021-02-20 12:00:00.00
(Grid= 01, Rec=5, Index=2, File: rivers.nc)
(Tmin= 18674.5000 Tmax= 18681.5000) t = 18678.5000
(Min = 1.57416667E+00 Max = 1.00000000E+01)
GET_NGFLD_NF90 - river runoff salinity, 2021-02-20 12:00:00.00
(Grid= 01, Rec=5, Index=2, File: rivers.nc)
(Tmin= 18674.5000 Tmax= 18681.5000) t = 18678.5000
(Min = 0.00000000E+00 Max = 0.00000000E+00)

If the sal is coming in at 0 for all time at all locations, there is not really a sign of that in the plots below (nice pcolors by the way).
what it the temp coming in for this point source?

4- what is the time series of Qin (river flow) at the problematic point? can you plot that?

5- use KANTHA_CLAYSON, not CANUTO. my experience with canuto is that it mixes a lot. it allows mixing at a higher Richardson number.

[pmaccc]

John,

Parker here - I am part of the project that Aurora is writing about. We have tried running versions of our model with HSIMT for salt + temp (instead of MPDATA) and it has given really poor results, essentially eroding most of the stratification in a model of Hood Canal. Hence I am not excited about this option.

Thanks for helping with this. It has been a real puzzle.

My own guess is that there is some bug in the zeta calculation in the presence of a vertical source that leads to an instability. I can't see how a baroclinic error would show up as a huge SSH jump.

Parker

[arango]

Hi Parker, we have a new 2D kernel for ROMS that is more accurate and efficient. It is based on Shchepetkin and McWilliams (2005, 2009) formulation for the Generalized Forward-Backward 3rd-order Adams-Bashword 4th-order Adams-Moulton (FB AB3-AM4) time-stepping algorithm. It does a barotropic pressure gradient correction due to the split scheme we didn't have in our original ROMS version. The code is still in research mode, and we will release it to the community in the future. If you are interested, we can give you access to that research repository and see what solution you will get.

[wilkin]

My own guess is that there is some bug in the zeta calculation in the presence of a vertical source that leads to an instability. I can't see how a baroclinic error would show up as a huge SSH jump.

Hi Parker,

Once you get an instability in the velocity the extreme values drive extreme convergence and there goes zeta.

We like Akima 4th order with some explicit geopotential diffusion for the dynamic tracers. I recommend you at least test the sensitivity of the set-up to the advection scheme and some horizontal smoothing. You are supplying fresh, warm water to the bottom of the water column so that it going to drive convective mixing.

[jcwarner]

ok but can you still try my other suggestions:

1- set Ninof to be 1, so we can see if there is more information.
3 what it the temp coming in for this point source?
4- what is the time series of Qin (river flow) at the problematic point? can you plot that?
5- use KANTHA_CLAYSON, not CANUTO. my experience with canuto is that it mixes a lot. it allows mixing at a higher Richardson number.

[aleeson]

Hi all,

Thank you for the suggestions!

I've attached a depth vs. time plot of in-situ density and vertical eddy viscosity/diffusivity at the vertical source. There doesn't appear to be a static instability.
I've also included a transport timeseries of the problematic source. Discharge temperature from this source is always 10 degC.
The run seems to stop halfway when ninfo is too small, but I'll try to make it as small as possible.

Right now, I am making my way through testing different mixing and advection schemes to see if they may help.
More updates to come.

Thank you again for your thoughts,
Aurora

[jcwarner]

"The run seems to stop halfway when ninfo is too small, but I'll try to make it as small as possible."
When you set Ninfo == 1, roms will make a call to Nonlinear/diag.F every time step

IF (MOD(iic(ng)-1,ninfo(ng)).eq.0) THEN

in here , roms computes max values for vel and dens, compares them to some high level, and makes a call to stop if the values are too high.
so if your run is stopping early with a small ninfo value, it might not be the point source that is causing the issue. it might be coming from somewhere else, and just showing up at the point souce.
please set ninfo = 1, and when it crashes, it should write the model state to nrst time step 3. take a look at that.

[pmaccc]

Dear ROMS folks, we are still having the problem Aurora describes above when using LwScr. She has tried many things to fix this:
-carefully checking the forcing files
-injecting the source at all levels, or at the top level, instead of the bottom level
-trying different advection schemes (A4 instead of MPDATA)
-turning off biogeochemistry
In all cases the model still crashes, sometimes running for a few months but when we look carefully there are large disturbances in velocity around the vertical source locations. The source transports are small so we would expect very little modification to the circulation.

So, I am wondering if there may still be bugs in this relatively new ROMS capability. Here are notes about places I am suspicious of, all from a recent version in ROMS/Nonlinear:

In omega.F there is this section:

Code: Select all

!  Apply mass point sources (volume vertical influx), if any.
!
!  Overwrite W(Isrc,Jsrc,k) with the same divergence of Huon,Hvom as
!  above but add in point source Qsrc(k) and reaccumulate the vertical
!  sum to obtain the correct net Qbar given in user input - J. Levin
!  (Jupiter Intelligence Inc.) and J. Wilkin
!
!    Dsrc(is) = 2,  flow across grid cell w-face (positive or negative)
!
        IF (LwSrc(ng)) THEN
          DO is=1,Nsrc(ng)
            IF (INT(SOURCES(ng)%Dsrc(is)).eq.2) THEN
              ii=SOURCES(ng)%Isrc(is)
              jj=SOURCES(ng)%Jsrc(is)
              IF (((IstrR.le.ii).and.(ii.le.IendR)).and.                &
     &            ((JstrR.le.jj).and.(jj.le.JendR)).and.                &
     &            (j.eq.jj)) THEN
                DO k=1,N(ng)
                  W(ii,jj,k)=W(ii,jj,k-1)-                              &
     &                       (Huon(ii+1,jj,k)-Huon(ii,jj,k)+            &
     &                        Hvom(ii,jj+1,k)-Hvom(ii,jj,k))+           &
     &                       SOURCES(ng)%Qsrc(is,k)
                END DO
              END IF
            END IF
          END DO
        END IF
!
        DO i=Istr,Iend
          wrk(i)=W(i,j,N(ng))/(z_w(i,j,N(ng))-z_w(i,j,0))
        END DO

In step2d_LF_AM3.h there is this section:

Code: Select all

!  Apply mass point sources (volume vertical influx), if any.
!
!    Dsrc(is) = 2,  flow across grid cell w-face (positive or negative)
!
      IF (LwSrc(ng)) THEN
        DO is=1,Nsrc(ng)
          IF (INT(SOURCES(ng)%Dsrc(is)).eq.2) THEN
            i=SOURCES(ng)%Isrc(is)
            j=SOURCES(ng)%Jsrc(is)
            IF (((IstrR.le.i).and.(i.le.IendR)).and.                    &
     &          ((JstrR.le.j).and.(j.le.JendR))) THEN
              zeta(i,j,knew)=zeta(i,j,knew)+                            &
     &                       SOURCES(ng)%Qbar(is)*                      &
     &                       pm(i,j)*pn(i,j)*dtfast(ng)
            END IF
          END IF
        END DO
      END IF

And in step3d_t.F there are these two sections:
(1) if we are using MPDATA

Code: Select all

!  If MPDATA and cell-centered point source (LwSrc), augment the
!  intermediate tracer, Ta, with the tracer divergence. For other
!  advection schemes LwSrc is applied after the vertical advection
!  is completed (J. Wilkin).
!
!    Dsrc(is) = 2,  flow across grid cell w-face (positive or negative)
!
          IF (LwSrc(ng)) THEN
            IF (Hadvection(itrc,ng)%MPDATA) THEN
              DO is=1,Nsrc(ng)
                IF (INT(SOURCES(ng)%Dsrc(is)).eq.2) THEN
                  Isrc=SOURCES(ng)%Isrc(is)
                  Jsrc=SOURCES(ng)%Jsrc(is)
                  IF (((Istr.le.Isrc).and.(Isrc.le.Iend+1)).and.        &
     &                 ((Jstr.le.Jsrc).and.(Jsrc.le.Jend+1)).and.       &
     &                 (j.eq.Jsrc)) THEN
                    DO k=1,N(ng)
                      cff=dt(ng)*pm(i,j)*pn(i,j)
                      IF (LtracerSrc(itrc,ng)) THEN
                        cff3=SOURCES(ng)%Tsrc(is,k,itrc)
                      ELSE
                        cff3=t(Isrc,Jsrc,k,3,itrc)
                      END IF
                      Ta(Isrc,Jsrc,k,itrc)=Ta(Isrc,Jsrc,k,itrc)+        &
     &                                     cff*SOURCES(ng)%Qsrc(is,k)*  &
     &                                     cff3
                    END DO
                  END IF
                END IF
              END DO
            END IF
          END IF

(2) or if we are not using MADATA

Code: Select all

!-----------------------------------------------------------------------
!  Add tracer divergence due to cell-centered (LwSrc) point sources.
!-----------------------------------------------------------------------
!
!  When LTracerSrc is .true. the inflowing concentration is Tsrc.
!  When LtracerSrc is .false. we add tracer mass to compensate for the
!  added volume to keep the receiving cell concentration unchanged.
!  J. Levin (Jupiter Intelligence Inc.) and J. Wilkin
!
!    Dsrc(is) = 2,  flow across grid cell w-face (positive or negative)
!
      IF (LwSrc(ng)) THEN
        DO itrc=1,NT(ng)
          IF (.not.((Hadvection(itrc,ng)%MPDATA).and.                   &
     &              (Vadvection(itrc,ng)%MPDATA))) THEN
            DO is=1,Nsrc(ng)
             IF (INT(SOURCES(ng)%Dsrc(is)).eq.2) THEN
                Isrc=SOURCES(ng)%Isrc(is)
                Jsrc=SOURCES(ng)%Jsrc(is)
                IF (((Istr.le.Isrc).and.(Isrc.le.Iend+1)).and.          &
     &              ((Jstr.le.Jsrc).and.(Jsrc.le.Jend+1))) THEN
                  DO k=1,N(ng)
                    cff=dt(ng)*pm(i,j)*pn(i,j)
# ifdef SPLINES_VDIFF
                    cff=cff*oHz(Isrc,Jsrc,k)
# endif
                    IF (LtracerSrc(itrc,ng)) THEN
                      cff3=SOURCES(ng)%Tsrc(is,k,itrc)
                    ELSE
                      cff3=t(Isrc,Jsrc,k,3,itrc)
                    END IF
                    t(Isrc,Jsrc,k,nnew,itrc)=t(Isrc,Jsrc,k,nnew,itrc)+  &
     &                                       cff*SOURCES(ng)%Qsrc(is,k)*&
     &                                       cff3
                  END DO
                END IF
              END IF
            END DO
          END IF
        END DO
      END IF

There are two things I am curious about in these:
(i) in some cases the if clause looks like:
IF (((IstrR.le.ii).and.(ii.le.IendR)).and. &
& ((JstrR.le.jj).and.(jj.le.JendR)).and. &
& (j.eq.jj)) THEN
(ii) in other cases it looks like:
IF (((Istr.le.Isrc).and.(Isrc.le.Iend+1)).and. &
& ((Jstr.le.Jsrc).and.(Jsrc.le.Jend+1)).and. &
& (j.eq.Jsrc)) THEN
(iii) and in other cases is looks like:
IF (((Istr.le.Isrc).and.(Isrc.le.Iend+1)).and. &
& ((Jstr.le.Jsrc).and.(Jsrc.le.Jend+1))) THEN

So, there are differences in whether or not +1 is added to the I,J limits, and there are differences in whether of not we are there is an additional j.eq.[] condition. If someone (John Wilkin) can confirm that these are all as intended then we will look elsewhere in our debugging journey.

Thanks,

Parker and Aurora

[wilkin]

If there was a bug here, and I haven't yet tunneled into the code to double check, it would affect whether the source was identified in the active i,j tile, or not, for all time. I don't see how it could be the direct cause of an instability that appears well into a run.

[arango]

I will look at it when I get the chance. The ROMS parallel tile indices are shown below and include the interior (red points) and nesting contact points (blue) in the purple area. Things get complicated when dealing with nested grids. The diagram shows the southwestern and northeastern parallel tiles at the opposite corners of the application grid with their respective indices names and values. Every horizontal staggered grid has its index. At rho-points, the whole model stencil ranges from 0:L and 0:M (including the boundaries). The thick red line demarks the physical domain. The interior computational rho points range from 1:Lm and 1:Mm and do not include the boundary points.

In parallel, the tile indices range is Istr:Iend and Jstr:Jend. The indices IstrR, IendR, JstrR, JendR are only relevant at the boundary edges of the domain since in interior partitions IstrR=Istr, IendR=Iend, JstrR=Jstr, and JendR=Jend. You must look at the above diagram carefully to assimilate all this information.

It seems to me that the code in omega.F and step2d_LF_AM3.h is fine.

The conditional for the MPDATADA section of step3d_t.F has an additional suspicious condition

Code: Select all

                  IF (((Istr.le.Isrc).and.(Isrc.le.Iend+1)).and.        &
     &                 ((Jstr.le.Jsrc).and.(Jsrc.le.Jend+1)).and.       &
     &                 (j.eq.Jsrc)) THEN

I will have to think why (j.eq.Jsrc) is added. I don't recall. The Iend+1 and Jend+1 is OK because it operates on the local variable Ta, which we need for the call to the mpdata_adiff_tile routine for the anti-diffusion correction.

However, in the non-MPDATA block,

Code: Select all

                IF (((Istr.le.Isrc).and.(Isrc.le.Iend+1)).and.          &
     &              ((Jstr.le.Jsrc).and.(Jsrc.le.Jend+1))) THEN

we are operating on the state variable t, and seems like a parallel bug to me. We will need the following:

Code: Select all

                IF (((Istr.le.Isrc).and.(Isrc.le.Iend)).and.            &
     &              ((Jstr.le.Jsrc).and.(Jsrc.le.Jend))) THEN

Thus, we can have a parallel bug if you have a point source at the arbitrary tile boundary. As I said, this must be tested carefully, and I don't have the time now. However, you can try it to see what kind of behavior you get in your application solution.

Good luck!

[wilkin]

Parker and Aurora,

I don't think this is your issue, but could you just check whether the point source giving you trouble happens to fall on a parallel tile boundary? If so, or even if not, can you test whether the blow-up you get also occurs at the same time in the run for a different choice of NtileI,NtileJ?

- John.

[aleeson]

Hi everyone,

I am reviving this thread because we have recently resolved this issue, and we'd like to recommend an improvement to step3d_uv.F.

In most files in ROMS/Nonlinear, such as pre_step3d.F, we see a pattern of "IF Dsrc == 0 THEN ....; ELSE IF Dsrc == 1 THEN ..."

Code: Select all

!  Apply tracers point sources to the horizontal advection terms,
!  if any.
!
!    Dsrc(is) = 0,  flow across grid cell u-face (positive or negative)
!    Dsrc(is) = 1,  flow across grid cell v-face (positive or negative)
!
          IF (LuvSrc(ng)) THEN
            DO is=1,Nsrc(ng)
              Isrc=SOURCES(ng)%Isrc(is)
              Jsrc=SOURCES(ng)%Jsrc(is)
              IF (((Istr.le.Isrc).and.(Isrc.le.Iend+1)).and.            &
     &            ((Jstr.le.Jsrc).and.(Jsrc.le.Jend+1))) THEN
                IF (INT(SOURCES(ng)%Dsrc(is)).eq.0) THEN
                  IF (LtracerSrc(itrc,ng)) THEN
                    FX(Isrc,Jsrc)=Huon(Isrc,Jsrc,k)*                    &
     &                            SOURCES(ng)%Tsrc(is,k,itrc)
                  ELSE
                    FX(Isrc,Jsrc)=0.0_r8
                  END IF
                ELSE IF (INT(SOURCES(ng)%Dsrc(is)).eq.1) THEN
                  IF (LtracerSrc(itrc,ng)) THEN
                    FE(Isrc,Jsrc)=Hvom(Isrc,Jsrc,k)*                    &
     &                            SOURCES(ng)%Tsrc(is,k,itrc)
                  ELSE
                     FE(Isrc,Jsrc)=0.0_r8
                  END IF
                END IF
              END IF
            END DO
          END IF

However, in step3d_uv.F we have only ""IF Dsrc == 0 THEN ....; ELSE ..."

Code: Select all

!-----------------------------------------------------------------------
!  Apply momentum transport point sources (like river runoff), if any.
!-----------------------------------------------------------------------
!
      IF (LuvSrc(ng)) THEN
        DO is=1,Nsrc(ng)
          i=SOURCES(ng)%Isrc(is)
          j=SOURCES(ng)%Jsrc(is)
          IF (((IstrR.le.i).and.(i.le.IendR)).and.                      &
     &        ((JstrR.le.j).and.(j.le.JendR))) THEN
            IF (INT(SOURCES(ng)%Dsrc(is)).eq.0) THEN
              DO k=1,N(ng)
                cff1=1.0_r8/(on_u(i,j)*                                 &
     &                       0.5_r8*(z_w(i-1,j,k)-z_w(i-1,j,k-1)+       &
     &                               z_w(i  ,j,k)-z_w(i  ,j,k-1)))
                u(i,j,k,nnew)=SOURCES(ng)%Qsrc(is,k)*cff1
              END DO
            ELSE
              DO k=1,N(ng)
                cff1=1.0_r8/(om_v(i,j)*                                 &
     &                       0.5_r8*(z_w(i,j-1,k)-z_w(i,j-1,k-1)+       &
     &                               z_w(i,j  ,k)-z_w(i,j  ,k-1)))
                v(i,j,k,nnew)=SOURCES(ng)%Qsrc(is,k)*cff1
              END DO
            END IF
          END IF
        END DO
      END IF

I believe this ELSE statement is un-intentionally assigning v-momentum to vertical sources. I've recently run a couple of tests in which I've changed the ELSE statement in step3d_uv.F to "ELSE IF (INT(SOURCES(ng)%Dsrc(is)).eq.1) THEN" and they seem to run without issue.

Let me know if you have any questions or thoughts. Thank you all again for your feedback and suggestions!

Best,
Aurora