ESCOMP / CTSM

Community Terrestrial Systems Model (includes the Community Land Model of CESM)
http://www.cesm.ucar.edu/models/cesm2.0/land/
Other
305 stars 309 forks source link

problem with the bounds of some history fields #16

Open billsacks opened 6 years ago

billsacks commented 6 years ago

Bill Sacks < sacks@ucar.edu > - 2013-08-17 07:46:54 -0600 Bugzilla Id: 1786 Bugzilla CC: bandre@lbl.gov, muszala@ucar.edu, mvertens@ucar.edu, rfisher@ucar.edu,

There is a potential problem with the bounds of some history fields. My guess is that this doesn't cause any problems now, but could cause problems in the future, if either (1) hist_update_hbuf was called within a threaded region (right now it's not), or (2) assumptions were made about the lower bound of arrays in hist_update_hbuf.

The problem arises from associating a pointer with an array slice, as in:

ptr => target(:, 1:n)

When you do this, the lower bound of ptr is reset to 1. Contrast this to:

ptr => target

in which case ptr retains the lower bounds of target.

Specifically, this occurs in:

(1) hist_update_hbuf_field_2d

field          => clmptr_ra(hpindex)%ptr(:,1:num2d)

(2) histFldsMod; e.g.:

      data1dptr => ccs%decomp_cpools(:,l)

(and maybe elsewhere - I haven't done an extensive search)

I believe this can be solved with the following syntax:

! In the following, we need to explicitly set the lower bound of 'field', otherwise
! it gets set to 1 when it's associated with the slice of 'ptr'
arr_lbound = lbound(clmptr_ra(hpindex)%ptr, 1)
field(arr_lbound: , 1:) => clmptr_ra(hpindex)%ptr(:,1:num2d)

but I haven't tested this.

For now, in the interest of time, I am working around this problem simply by NOT explicitly specifying the bounds of the history fields in calls to p2g/c2g/l2g in hist_update_hbuf; e.g., I am calling these routines like:

      call p2g(bounds, &
           field, &
           field_gcell(bounds%begg:bounds%endg), &
           p2c_scale_type, c2l_scale_type, l2g_scale_type)

rather than like:

      call p2g(bounds, &
           field(bounds%begp:bounds%endp), &
           field_gcell(bounds%begg:bounds%endg), &
           p2c_scale_type, c2l_scale_type, l2g_scale_type)

I think this should be okay for now, but wouldn't work if this was called from within a threaded region.

billsacks commented 6 years ago

Bill Sacks < sacks@ucar.edu > - 2013-08-26 06:06:31 -0600

I think that the trickier part of this problem - the calls from histFldsMod (#2 in my original report) - can be solved by changing the dummy argument declarations in hist_addfld1d and hist_addfld2d, by adding lower bounds to these arguments:

real(r8)        , optional, pointer    :: ptr_gcell(:)   ! pointer to gridcell array
real(r8)        , optional, pointer    :: ptr_lunit(:)   ! pointer to landunit array
real(r8)        , optional, pointer    :: ptr_col(:)     ! pointer to column array
real(r8)        , optional, pointer    :: ptr_pft(:)     ! pointer to pft array

(I missed these in my initial rework of subroutine arguments.)

To do this, we would need to pass bounds to hist_addfld1d and hist_addfld2d - ensuring that these are the proc bounds rather than the clump bounds.

At that point, I believe that the only remaining problem would be #1 in my initial report, which could be solved via my original suggestion.

I'm not sure whether this refactoring is worth the time right now, though.

billsacks commented 6 years ago

Bill Sacks < sacks@ucar.edu > - 2013-09-19 14:31:41 -0600

I am adding more code in histFileMod where this issue arises - specifically, code to deal with multi-layer snow history fields.

Specifically, I have this code on my branch, in hist_update_hbuf_field_2d:

if (is_snow_layer_field) then
   ! For multi-layer snow fields, build a special output variable that handles
   ! missing snow layers appropriately

   ! Note that the following allocation is not what we would want if this routine
   ! were operating in a threaded region (or, more generally, within a loop over
   ! nclumps) - in that case we would want to use the bounds information for this
   ! clump. But currently that's not possible because the bounds of some fields have
   ! been reset to 1 - see also bug 1786.
   allocate(field(lbound(clmptr_ra(hpindex)%ptr, 1) : ubound(clmptr_ra(hpindex)%ptr, 1), 1:num2d))

Ideally, we would want field to be allocated to be just big enough for the given clump bounds (if we were within a loop over clumps). But that's not possible if the field's lower bound has been reset to 1.

So, as with my above notes, this code will only work as intended if the bounds passed to hist_update_hbuf_field_2d are proc bounds.

Once this bug is addressed, we should search for '1786' in histFileMod to find all notes about this problem.

billsacks commented 6 years ago

Bill Sacks < sacks@ucar.edu > - 2013-09-19 14:32:41 -0600

It looks like the syntax I suggested:

field(arr_lbound: , 1:) => clmptr_ra(hpindex)%ptr(:,1:num2d)

is a fortran2003 feature, which should now be supported by the major compilers.

billsacks commented 6 years ago

Bill Sacks < sacks@ucar.edu > - 2016-05-16 14:26:27 -0600

Probably deferring this until threading is reworked in CLM to be done at a higher level - at which point this bug will either go away or be changed considerably