Release points on land don't show up in initial checks

[salwood82]

Before moving the particles cms first checks (in loop.f90) whether a) there are nest files to move the particle, b) the particle release location is inside the nest grid and c) the particle release position is not on land. Whilst a and b work fine, c does not pick up on release points which are on land. loop.f90 therefore calls subroutine move (move.f90) to try to advect the particle, at which point move realises the particle is on land and can't move it.

The result is that, in the output, particles released outside of the grid exit immediately (exitcode -1), with only their release position given in the output file (cms also now prints a warning that this has happened). However, for a particle released on land, cms tries to move the particle for 1 time step before failing (exitcode -2), with 2 identical lat/lon positions given in the traj output file.

From testing with a particle that is definitely on land (using 2D nest files, no ibm, avoidcoast = true) I have got so far in identifying the issue...

1. loop.f90 calls subroutine rungakutta (rungekutta.f90) to check the start position, which in turn calls subroutine fld3_interp (fldinterp.f90) to get the u and v velocities at the start position.

2. At line 119, fld3_interp assigns all values of fld(i,j) to 0.00, instead of the CMS fill value (1.2676506E30):

DO counter = 1, nweight
      fld2d(i,j) = fld2d(i,j) + weight(counter) * fld(new_x,y,z,counter)
ENDDO

Gives...
counter = 1:
fld2d(i,j) = 0.0E+0 + -2.2245625E-12 * 0.0E+0 = 0.0E+0
counter = 2:
fld2d(i,j) = 0.0E+0 + 1.5304982E-41 * 0.0E+0 = 0.0E+0

None of the land flags are therefore triggered in fld3_interp, which returns u and v velocities of 0.00 to rungakutta, which in turn returns no errors to loop.f90, which just thinks that the particle is free to move.

3. loop.f90 then calls subroutine move (move.f90) to try to advect the particle. Move.f90 calls subroutine rungakutta which again calls fld3_interp to get the u and v values.

This time, at line 119, fld(i,j) is set to the correct fill value:

counter = 1:
fld2d(i,j) = 0.0E+0 + 0.9583333 * 1.2676506E+30 = 1.2148318E+30
counter = 2:
fld2d(i,j) = 1.2148318E+30 + 4.1666668E-02 * 1.2676506E+30 = 1.2676506E+30

This then, eventually further down the script, triggers landFlag = true, which is fed back to rungakutta and then move.f90.

4. move.f90 uses the landFlag trigger to set flag1 and flag2 to true (lines 111-124):

landFlag = .false.
CALL rungakutta(ngrid,xold,yold,zold,run_time,diam,dens,1,u1,v1,w1,ztemp,zsaln,zdens,flag,landFlag,r,n,grid_i,grid_j,grid_k)
     IF(flag(1) .and. flag(2)) THEN        **! these flags are not triggered**
      goto 100
     ENDIF
     IF (landFlag) THEN                          **! this flag is triggered**
      flag(1) = .true.
      flag(2) = .true.
      goto 100
     ENDIF

Is anyone able to help me figure out why in the first case (when called from loop>rungakutta>fldinterp), fld(new_x,y,z,counter) is returned as 0.0, whereas in the 2nd case (loop>move>rungakutta>fldinterp) fld(new_x,y,z,counter) is returned as the correct fill value of 1.2676506E+30?

It must be the different ways loop.f90 and move.f90 call rungakutta, as nothing has changed in terms of the particle between the 2 calls? But so far the only differences I can see is that a) loop.f90 calls rungakutta with landFlag = true (line 229), whilst move.f90 calls it with landFlag = false (line 111), and b) the run time is 0 in the call from loop.f90 and 3600s (my time-step) in the call from move.f90.

I have tried calling RK from move.f90 (RK step 1) with landFlag = true, and from loop.f90 with startsec = 3600 but this doesn't seem to make any difference.

So... my flow of events is:

 loop.f90 calls rungakutta with landFlag = true to check start position
    RK calls fld3_interp for u/v_vel
    NO landFlag triggered in rungakutta
All flags returned false to loop.f90

loop.f90 calls move...
    move.f90 1st call to RK with LF = false
        RK calls fld3_interp for u/v_vel
        landFlag triggered in rungakutta
        RK calls fld3_interp for wvel (why?)
    landFlag triggered in move.f90
move.f90 returns particle(r)%flag(n,1) .and. particle(r)%flag(n,2) = true to loop.f90
loop.f90 writes -2 exitcode to output file and sets particle(r)%move(n) = .false.

[salwood82]

Identified issue: when the initial check for the start position is carried out, no datafiles have been read in, so rungakutta was feeding a matrix of zeros into fld3_interp as the uvel/vvel matrices, therefore fld3_interp returned a value of 0.00 instead of the fill value.