point source values replaced by FillValue

[jcwarner]

Another reason why the FillValues give me a headache: Most users place point sources at locations where the u- or v-mask is zero (i.e. land). During the simulation the model will actually compute a source value at that 'land' location and those transports are used in the simulation. However, when writing history or other output files, the u- or v-values at the ppoint source locations will be replaced with FillValues.

This does not allow the user to check if the correct transport is being introduced into the model. Please negate the FillValue at locations of the point sources. Suggested fix: make up a temporary u/vmask that does not mask the point sources in wrt_his.F send that temporary mask to nf_fwite3d instead of GRID(ng)%umask.

[arango]

Yes, thank you. This was a good one and it required complicate fixes and testing:

• Created several new masking arrays in mod_grids.F to process output history fields and time-averaged fields and time-averaged diagnostics:

  # ifdef AVERAGES
              GRID(ng) % pmask_avg(i,j)
              GRID(ng) % rmask_avg(i,j)
              GRID(ng) % umask_avg(i,j)
              GRID(ng) % vmask_avg(i,j)
  # endif
  # ifdef DIAGNOSTICS
              GRID(ng) % pmask_dia(i,j)
              GRID(ng) % rmask_dia(i,j)
              GRID(ng) % umask_dia(i,j)
              GRID(ng) % vmask_dia(i,j)
  # endif
              GRID(ng) % pmask_io(i,j)
              GRID(ng) % rmask_io(i,j)
              GRID(ng) % umask_io(i,j)
              GRID(ng) % vmask_io(i,j)
  

  These masking arrays include the regular grid masks, point sources locations, and wetting and drying time dependent masks.

• A new routine set_masks.F is introduced to initialize all the internal mask arrays. It is called during initialization after the point sources information is known.

        IF (iic(ng).eq.ntstart(ng)) THEN
          ...
  !
  !  Initialize other state variables.
  !
  !$OMP PARALLEL DO PRIVATE(thread,subs,tile) SHARED(ng,numthreads)
          DO thread=0,numthreads-1
            subs=NtileX(ng)*NtileE(ng)/numthreads
            DO tile=subs*(thread+1)-1,subs*thread,-1
              CALL ini_fields (ng, TILE, iNLM)
  # ifdef MASKING
              CALL set_masks (ng, TILE, iNLM)
  # endif
            END DO
          END DO
  !$OMP END PARALLEL DO
        END IF
  

• If wetting and drying (WET_DRY) is activated, the time-averaged fields are only accumulated and counted when the grid points are wet. During time averaging, the rmask_avg or rmask_dia and friends arrays contain the counter when each point is wet. It will be converted to time-averaged masks (dry=0, wet=1) at the end of the averaging period, so it can be used in output routines wrt_avg.F and wrt_diags.F as real FillValue masks. This is a correction to the code. We were not getting the correct time averages when WET_DRY is activated.

• The time-dependent rmask_io and friends are computed in routine wetdry.F when wetting and drying is activated.

• A parallel bug is corrected when computing umask_wet and vmask_wet. We needed to use the private mask array wetdry instead of rmask_wet in the IF statement. For example:

          DO j=JstrR,JendR
            DO i=Istr,IendR
              umask_wet(i,j)=1.0_r8
              IF (DU_avg1(i,j).eq.0.0_r8) THEN
                IF ((wetdry(i-1,j)+wetdry(i,j)).le.1.0_r8) THEN
                  umask_wet(i,j)=0.0_r8
                END IF
              END IF
            END DO
          END DO
  

• The accumulated fields are not converted to time averages in routine set_diags.F instead of wrt_diags.F. The following statement:

            scale=1.0_r8/(REAL(nDIA(ng),r8)*dt(ng))
  
  is changed to
  
            scale=1.0_r8/dt(ng)
  

  We not longer will get identical diagnostic fields because of this change of order in the computation due to round off. It gives an error of the order 1E-20.

• Warning: The time-averaged free-surface (zeta) is FillValue masked in wetting and drying, so it is processed as zero inside ROMS. If you use this time-averaged field as input initial conditions, it needs to be pre-processed to include a non-zero Dcrit(ng) value in dry points. See step2d for details. Time-averaged fields have a different meaning during wetting and drying. The instantaneous field in history and restart fieles, however, is processed correctly with the optional argument SetFillVal = .FALSE. argument to routine nf_fwrite2d.