Skip to content

SCOPS: top-level threshold uses ran(npoints) (one scalar for all gridpoints) instead of ran(1:npoints) #162

Description

@duncanwp

Location: subsample_and_optics_example/subcol/scops.F90, line 136 (master 5eb05e5), top-level threshold initialization for the non-maximum-overlap case:

DO ibox=1,ncol
   ran(1:npoints) = get_rng(RNGS)
   ! select random pixels from the non-convective
   ! part the gridbox ( some will be converted into
   ! convective pixels below )
   threshold(1:npoints,ibox) = conv(1:npoints,ilev)+(1-conv(1:npoints,ilev))*ran(npoints)
enddo

The right-hand side uses ran(npoints) — the last gridpoint's random number — for every gridpoint, rather than the elementwise ran(1:npoints) used everywhere else (e.g. the per-level reset at line 199). A whole vector of random numbers is drawn and all but one discarded.

Impact: at the top model level all gridpoints share a single random threshold per subcolumn, correlating the level-1 subcolumn assignment across the domain (and, through the overlap recursion, weakly imprinting on levels below the top). Statistically minor in practice because level 1 is rarely cloudy, but clearly not the intent, and it makes the top-level sampling depend on npoints (i.e. on chunking).

Suggested fix: ... * ran(1:npoints).

Found while building a JAX translation of COSP (jax-cosp), which draws an independent uniform per gridpoint.

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions