Another bug in /Nonlinear/step3d_t.F when using Q_PSOURCE

[lanerolle]

I found another bug in step3d_t.F when using the Q_PSOURCE CPP options. If you look at the code segment:

Code: Select all

  +902  # ifdef Q_PSOURCE
  +903  !
  +904  !  Apply mass point sources - Volume influx.
  +905  !
  +906            DO is=1,Nsrc
  +907              i=Isrc(is)
  +908              IF (LtracerSrc(itrc,ng)) THEN
  +909                IF (((IstrR.le.i).and.(i.le.IendR)).and.                  &
  +910       &            ((JstrR.le.j).and.(j.le.JendR))) THEN
  +911                  IF (j.eq.Jsrc(is)) THEN
  +912                    DO k=1,N(ng)
  +913                      FC(i,k)=FC(i,k)+0.5_r8*                             &
  +914       &                      (Qsrc(is,k  )*Tsrc(is,k  ,itrc)+            &
  +915       &                       Qsrc(is,k+1)*Tsrc(is,k+1,itrc))
  +916                    END DO
  +917                  END IF
  +918                END IF
  +919              END IF
  +920            END DO
  +921  # endif

You will see that Qsrc(:,:,N+1), Tsrc(:,:,N+1) elements are accessed by ROMS but these variables are defined only to have N elements along their last dimension (see Modules/mod_sources.F). Hence, if you run ROMS with the Q_PSOURCE CPP option, it blows-up because these last elements have (machine) infinite values such as 1.000e+131! I guess if we use the -g -C -u Fortran compiler options this bug could be trapped and the code will exit but I used the -O option for convenience.

Could someone please let me know how to remedy this bug - can we simply not carry-out the averaging on lines 914-915 in step3d_t.F?

[wilkin]

The FC are fluxes at the cell faces (w-points) so it
looks to me like it's simply a matter of changing line 912 to

DO k=1,N(ng)-1

The boundary conditions on the fluxes at k=0 and k=N were set a few
lines above. Try it and tell us if you then have the correct next influx.

[arango]

Yes, good catch The k-loop needs to be k=1,N(ng)-1 as John mentioned. Please check the following trac ticket. You may update your code. Thank you for reporting this bug.

[lanerolle]

I tried John's suggestion and it works well - Thank you. The only issue I have is that I did a run using UV_PSOURCE and another with Q_PSOURCE (also used TS_PSOURCE for both of them) for a passive tracer and the vertical distributions from these two runs are very different. In the former, the passive tracer is distributed through the water column with some stratification. However, in the latter, it is pretty much confined to the top of the water column in a thinnish-layer. I used the same river.nc file for both runs (i.e. same river transport, river_Vshape, river_dye_NN values, etc.).

Now, in ROMS/Nonlinear/omega.F, when using the Q_PSOURCE CPP option, we make a contribution to the w-velocity as:

w=w+V*river_Vshape/(dx*dy) where V=V(t) is the volume transport (m3/s).

For my application, the additional contribution is ~0.00001 m/s which is smallish but I guess if it persists over a long period, it will cause things to move up the water column and be confined to the top.

What do you think?