Updated public repository release, so it is synchronized with changes in several of my research branches:

• Updated mod_netcdf.F to include more information if a reading or writing error occurs. It now has the specific NetCDF library error by using the string from nf90_strerror(status).

• Changed some of the formatted integer fields to I0. The I0 edit descriptor was introduced in the 2003 Fortran Standard and produce a result with many digits as are necessary to represent the integer value. I love this capability. Well done Fortran Gods!

• Corrected a small bug in timers.F routine wclocks_off. We need have the following USE statement:

  #ifdef DISTRIBUTE
  !
        USE distribute_mod, ONLY : mp_barrier, mp_reduce
        USE distribute_mod, ONLY : mp_collect
        USE strings_mod,    ONLY : uppercase
  #endif
  

  The mp_collect was restricted to CPP ROMS_STDOUT in a recent update.

Many thanks to Frans van Eeden ​reporting this issue.

In ROMS, the river runoff forcing is modeled with point sources. There are two ways how is implemented as horizontal advection transport (LuvSrc) or volume vertical influx (LwSrc). These options are specified in roms.in:

! Logical switches (TRUE/FALSE) to activate horizontal momentum transport
! point Sources/Sinks (like river runoff transport) and mass point
! Sources/Sinks (like volume vertical influx), [1:Ngrids].

      LuvSrc == T                          ! horizontal momentum transport
       LwSrc == F                          !  volume vertical influx

The volume vertical influx was not working correctly, and it was reformulated. In omega.F now we have:

!
!  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
!
        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-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 DO
        END IF

Also, in step3d_t.F we have:

!
!-----------------------------------------------------------------------
!  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
!
      IF (LwSrc(ng)) THEN
        IF (.not.((Hadvection(itrc,ng)%MPDATA).and.                     &
     &            (Vadvection(itrc,ng)%MPDATA))) THEN
          DO itrc=1,NT(ng)
            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
                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 DO
          END DO
        END IF
      END IF

Many thanks to John Wilkin and Julia Levin for redesigning and testing the LwSrc scheme. John tested the scheme extensively with all the tracer advection options. Now, both LuvSrc and LwSrc have exact volume conservation.

Both methods produce similar results. In the RIVERPLUME1 test case, we get the following results with the HSIMT scheme for temperature and salinity:

The tangent linear (TLM), representer (RPM), and adjoint versions (ADM) of the scheme were updated accordingly. Extra terms were added to the TLM, RPM, and ADM formulation to allow point-source sensitivity studies and to be part of the 4D-Var control vector for analysis and forecast impacts.

Also, the tl_step3d_t.F, rp_step3d_t.F, and ad_step3d_t.F were cleaned by removing the unsupported code for MPDATA advection of tracer.

John Warner brought to my attention that the nesting algorithms were computing the vertical interpolation weights (Vweight) at every timestep for each nested grid.

Currently, the vertical interpolation weights are used in composite grids because their grids are usually not coincident. They are not needed in refinement grids because the donor and receiver grids have the same number of vertical levels and matching bathymetry. However, in the future, it is possible to have configurations that require vertical interpolation weights in refinement. The switch get_Vweights is introduced to control if such weights are computed or not. If false, it will accelerate computations because of less distributed-memory communications.

Therefore, in nesting.F we now have:

     IF ((isection.eq.nzwgt).and.get_Vweights) THEN
       DO tile=last_tile(ng),first_tile(ng),-1
         CALL z_weights (ng, model, tile)
       END DO
       RETURN
     END IF

Here, get_Vweiths is a new logical switch added in mod_nesting.F. It is initialized in set_contact.F as:

!
!  Set the switch to compute vertical interpolation weights. Currently,
!  they are only needed in non-coincident composite grids.
!
!
      IF (.not.ANY(Lcoincident).and.ANY(Lcomposite)) THEN
        get_Vweights=.TRUE.
      ELSE
        get_Vweights=.FALSE.
      END IF

where the coincident and composite variables are the set in the nested grid configuration NetCDF file and computed in Matlab script contact.m.

Profiling:

Several runs were made to measure the improvement in nested applications with refinement grids (no composite grids):

• The LAKE_JERSEY test case (grids a and d) with one refinement grid runs 21.11% faster on 12 processes.

• The LAKE_JERSEY test case (grids a, c, d, and e) with tree refinement grids runs 36.83% faster on 12 processes.

• Our operational US East Coast application (grids DOPPIO, PIONEER, and ARRAY) with 2 telescoping refinement grid runs 24.77% faster on 12 processes.

That's quite an improvement with such a simple change in the code. Obviously, as the number of nested grids are increased, the improvement is greater.