bug in ana_psource.h with analytical initial conditions

Bug reports, work arounds and fixes

Moderators: arango, robertson

Post Reply
Message
Author
jpringle
Posts: 88
Joined: Sun Jul 27, 2003 6:49 pm
Location: UNH, USA

bug in ana_psource.h with analytical initial conditions

#1 Post by jpringle » Wed Jan 31, 2018 4:37 pm

I think I found a bug in the ana_psource.h when the model is initialized with ana_initial.h. In this case, the analytical forcing files are invoked in analytical.F before the first time step of integration as part of initializing the u, v etc fields. When this is done, the variable iic(ng) is 0.

In ana_psource.f, the location of the river fluxes is defined with the following code:

Code: Select all

!
!-----------------------------------------------------------------------
!  If initialization, set point Sources and/or Sinks locations.
!-----------------------------------------------------------------------
!
      IF (iic(ng).eq.ntstart(ng)) THEN
!
This test is not true when it is called from analytical.F, and the number of source locations Nsrc(ng) and the locations SOURCES(ng)%Isrc and Jsrc are set to however the compiler initializes them. As the code stands now, Nsrc(ng) is 200, and Jsrc starts out as 0. Further on in ana_psource.h these values are used. For example, in the RIVERPLUME2 case there is code like

Code: Select all

              SOURCES(ng)%Qshape(is,k)=cff*                             &
     &                                 (z_w(i,j-1,k    )-               &
     &                                  z_w(i,j-1,k-1  )+               &
...
j in this code is take from SOURCES(ng)%Jsrc, so z_w(i,j-1,k-1) is then an out of bounds reference.

To fix this, we can change the test for initialization in ana_psource.h to

Code: Select all

!
!-----------------------------------------------------------------------
!  If initialization, set point Sources and/or Sinks locations.
!-----------------------------------------------------------------------
!
      IF ((iic(ng).eq.ntstart(ng)).or.(iic(ng).eq.0)) THEN  
Or initialize Nsrc(ng) to 0.

Cheers,
Jamie

Post Reply