Model blow-up after adding river discharge

[txu68]

Hello,

I am using ROMS to simulate tidal flow in a coastal region containing a long and curved river. The grid also includes wetlands, so I have turned on WET_DRY. This setting runs smoothly to the end.

Now I want to add a point source to the upstream to account for the river discharge. I have turned on ANA_PSOURCE instead of making a Netcdf file of my own.

The related setting in ocean.in is

LuvSrc = T
LtracerSrc = F T

Related setting in ana_psource is
--------------------------------------------------------------------
IF (Master.and.DOMAIN(ng)%SouthWest_Test(tile)) THEN
Nsrc(ng)=2
SOURCES(ng)%Dsrc(1)=1.0_r8
SOURCES(ng)%Isrc(1)=2
SOURCES(ng)%Jsrc(1)=368
SOURCES(ng)%Dsrc(Nsrc(ng))=1.0_r8
SOURCES(ng)%Isrc(Nsrc(ng))=3
SOURCES(ng)%Jsrc(Nsrc(ng))=368
END IF
------------------------------------------------------------------
DO k=1,N(ng)
DO is=1,Nsrc(ng)
i=SOURCES(ng)%Isrc(is)
j=SOURCES(ng)%Jsrc(is)
IF (((IstrT.le.i).and.(i.le.IendT)).and. &
& ((JstrT.le.j).and.(j.le.JendT))) THEN
IF (ubar(i,j,knew).ne.0.0_r8) THEN
cff=ABS(u(i,j,k,nnew)/ubar(i,j,knew))
ELSE
cff=1.0_r8
END IF
SOURCES(ng)%Qshape(is,k)=cff* &
& (z_w(i-1,j,k )- &
& z_w(i-1,j,k-1 )+ &
& z_w(i ,j,k )- &
& z_w(i ,j,k-1 ))/ &
& (z_w(i-1,j,N(ng))- &
& z_w(i-1,j,0 )+ &
& z_w(i ,j,N(ng))- &
& z_w(i ,j,0 ))
END IF
END DO
END DO
-----------------------------------------------------------------------------
IF ((tdays(ng)-dstart).lt.0.5_r8) THEN
fac=1.0_r8+TANH((time(ng)-43200.0_r8)/43200.0_r8)
ELSE
fac=1.0_r8
END IF
DO is=1,Nsrc(ng)
SOURCES(ng)%Qbar(is)=fac*100.0_r8
END DO

The model only runs for 2 minutes before blow-up. The error message said,
DIAG speed trouble 24.64562562679086
DIAG speed ijk 34.85276857974117 2 368 8

When I plot the result from restart file, I can see the water level rising up from upstream. So I think I got the locations right for the point sources. Can anyone tell me how to deal with this problem?

Thanks,
Amy

[kate]

When you specify a river inflow, you have to provide values for all of the tracers. You can't have:

Code: Select all

LtracerSrc = F T

without developing an instability. I'm not sure your Qshape is done in a stable fashion either.

[txu68]

Hello Kate,

Thanks for the quick response. I really appreciate it.

Do you suggest I should change to

Code: Select all

LtracerSrc = T T

or

Code: Select all

LtracerSrc = T T T T

For the Qshape, I am simply using the example from riverplume1 in the ana_psource. What do you mean that it may not be stable?

Amy

[txu68]

I have changed the Qshape to

Code: Select all

       DO k=1,N(ng)
          DO is=1,Nsrc(ng)
            i=SOURCES(ng)%Isrc(is)
            j=SOURCES(ng)%Jsrc(is)
            IF (((IstrT.le.i).and.(i.le.IendT)).and.                    &
     &           ((JstrT.le.j).and.(j.le.JendT))) THEN
              IF (ubar(i,j,knew).ne.0.0_r8) THEN
                cff=ABS(u(i,j,k,nnew)/ubar(i,j,knew))
              ELSE
                cff=1.0_r8
              END IF
              SOURCES(ng)%Qshape(is,k)=cff*                             &
     &                                 (z_w(i,j-1,k    )-               &
     &                                  z_w(i,j-1,k-1  )+               &
     &                                  z_w(i,j  ,k    )-               &
     &                                  z_w(i,j  ,k-1  ))/              &
     &                                 (z_w(i,j-1,N(ng))-               &
     &                                  z_w(i,j-1,0    )+               &
     &                                  z_w(i,j  ,N(ng))-               &
     &                                  z_w(i,j  ,0    ))
            END IF
          END DO
        END DO

With

Code: Select all

LtracerSrc = T T

and,

Code: Select all

        IF (DOMAIN(ng)%NorthEast_Test(tile)) THEN
          DO k=1,N(ng)
            DO is=1,Nsrc(ng)
              SOURCES(ng)%Tsrc(is,k,itemp)=T0(ng)
              SOURCES(ng)%Tsrc(is,k,isalt)=0.0_r8
            END DO
          END DO
        END IF

The model runs for 30 minutes before blow-up. Any suggestions?

Thanks,
Amy

[kate]

How many tracers do you have? They all need values on inflow. Have you looked at the blow-up?

I don't know about the Qshape. I just specify mine to be a fixed function of k.

[kate]

Sure, why not? The trick is in finding an appropriate value.

Ken Coyle and I have had to change the ROMS "speed limit" for some major rivers. If you work out the velocity of the inflow as specified, is it ever above 20 m/s? If so, you can change the speed limit in mod_scalars.F:

Code: Select all

#ifdef NWGOA
        real(r8) :: max_speed = 80.0_r8         ! m/s
#else
        real(r8) :: max_speed = 20.0_r8         ! m/s
#endif

[txu68]

In ocean.in,

Code: Select all

NAT = 2

I don't have any other tracers. You said they all needed to be specified in the inflow. I don't know how to set the values in ana_psource, except changing LtracerSrc = T T T T in ocean.in.

I have looked up the blow up, but I don't see anything wield in salinity or temperature or zeta. ubar can be 10^37 in the masked region, but I don't think this is relevant.

Thanks,
Amy

[kate]

LtracerSrc = T T T T

This is four tracers, so the two NAT tracers plus two passive tracers. ROMS will likely ignore the last two if there are no passive tracers.

If the velocity is going bad, is it largest right at the source or a grid or two away? I sometimes add stations right there to watch what's going on. For a handful of stations in a short job like this, you can write to the stations file every timestep.