Input river flow with LwSrc instead of LuvSrc

[julia]

Most of ROMS users use LuvSrc=T and LwSrc=F to input river flow, and we are advised to use LuvSrc=T. But when a setup has many rivers (ex. fine resolution estuary and river), it may become preferable to use LwSrc=T and LuvSrc=F: it is easier to set up rivers this way; since river flow goes directly into a water column on rho-grid, there is no need to bother with river direction (u or v) and mass flux has the same sign. LwSrc=T is also more stable, since no horizontal momentum is added to the system.

Unfortunately, this option is buggy, and in a case when river flow is strong, it can lead to temperatures quickly going out of range. Looking at the code together with John Wilkin and Chuning Wang, we discovered the source of the problem and came up with the modifications to the code. Below we discuss the modifications. We also attach the files with the modifications.

There are several routines that are responsible for implementation of LwSrc. First, volume flux is added directly to the water column by increasing surface height in Nonlinear/step2d_LF_AM3.h.

Then, since this operation influences flow divergence, vertical velocity needs to be adjusted. This is done in omega.F, where velocity divergence is integrated from the bottom up, but the mass source is not integrated from the bottom up, it just added on each sigma level, as follows:

Code: Select all

        IF (LwSrc(ng)) THEN
          DO is=1,Nsrc(ng)
            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)+SOURCES(ng)%Qsrc(is,k)
              END DO
            END IF
          END DO

In order to obtain correct vertical integral of W, the code above has to be changed to this:

Code: Select all

        IF (LwSrc(ng)) THEN
          DO is=1,Nsrc(ng)
            ii=SOURCES(ng)%Isrc(is)
            jj=SOURCES(ng)%Jsrc(is)
            IF (((IstrR.le.ii).and.(ii.le.IendR)).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)
 

In this case, at river locations we recompute horizontal flow divergence with additional river source from the bottom up.

Heat and fresh water fluxes from LwSrc are done in pre_step3d.F and stepd3d_t.F. They are currently done through vertical advection, like this:

Code: Select all

         IF (LwSrc(ng).and.ANY(LtracerSrc(:,ng))) THEN
            DO is=1,Nsrc(ng)
              Isrc=SOURCES(ng)%Isrc(is)
              Jsrc=SOURCES(ng)%Jsrc(is)
              IF (LtracerSrc(itrc,ng).and.                              &
     &            ((IstrR.le.Isrc).and.(Isrc.le.IendR)).and.            &
     &            (j.eq.Jsrc)) THEN
                DO k=1,N(ng)-1
                  FC(Isrc,k)=FC(Isrc,k)+0.5_r8*                         &
     &                       (SOURCES(ng)%Qsrc(is,k  )*                 &
     &                        SOURCES(ng)%Tsrc(is,k  ,itrc)+            &
     &                        SOURCES(ng)%Qsrc(is,k+1)*                 &
     &                        SOURCES(ng)%Tsrc(is,k+1,itrc))

It can be easily seen that the fluxes do no add up properly. For example, if a river contributes uniformly to a water column (i.e. Qsrc are all equal, and Tsrc is constant), then after taking the difference of FC terms, there are no contribution from the river to the T and S at all. Instead of implementing these river fluxes through vertical advection, we suggest adding them as a source term to the tracer equations. To do this, we need to remove the LwSrc portion of the code in pre_step3d.F and step3d_t.F . Instead in step3d_t.F we add the source term coded like this

Code: Select all

          IF (LwSrc(ng)) THEN
            DO is=1,Nsrc(ng)
              Isrc=SOURCES(ng)%Isrc(is)
              Jsrc=SOURCES(ng)%Jsrc(is)
              IF (                                                      &
     &            ((IstrR.le.Isrc).and.(Isrc.le.IendR)).and.            &
     &            ((JstrR.le.Jsrc).and.(Jsrc.le.JendR))) THEN
                  cff=dt(ng)*pm(Isrc,Jsrc)*pn(Isrc,Jsrc)
                  IF (LtracerSrc(itrc,ng)) THEN
                    t(Isrc,Jsrc,k,nnew,itrc)=t(Isrc,Jsrc,k,nnew,itrc)+  &
     &                           cff*SOURCES(ng)%Qsrc(is,k)*            &            
     &                           SOURCES(ng)%Tsrc(is,k,itrc)
                  ELSE
                    t(Isrc,Jsrc,k,nnew,itrc)=t(Isrc,Jsrc,k,nnew,itrc)+  &
     &                           cff*SOURCES(ng)%Qsrc(is,k)*            &            
     &                           t(Isrc,Jsrc,k,3,itrc)

and place this code in front of the section where time stepping of vertical advection is performed. This code also takes into account whether we use ambient T/S, or the ones from the river.

This seems to have made the difference in my runs with LwSrc=T. The temperature stopped going out of range, and all the fluxes add up correctly to agree with the LuvSrc=T results.

Acknowledgement: Jupiter Intelligence INC has paid for my work on discovering and correcting this issue. Jupiter agrees to release this code to the ROMS community.