icepack_algae underflow

[apcraig]

icepack_algae is underflowing in the second part of this line with gnu compiler and debug+bgcISPOL.

           dmobile(k) = mobile(mm)*(initcons_mobile(k)*(exp(-dt*rtau_ret( mm))-c1) + &
                             initcons_stationary(k)*(c1-exp(-dt*rtau_rel(mm))))

If anyone has any ideas, I'm happy to implement and test.

[eclare108213]

This needs to be recoded. Either rtau_ret or rtau_rel will always be zero, and so one of the exp() will always underflow. It would be better to put the calculation of the exp in the conditional block above this, where rtau_ret and rtau_rel are computed, and then just use the values here.

[apcraig]

I understand what you're getting at. When those values are zero, I assume the exp should be treated as zero?

So if rtau_ret is zero,
dmobile(k) = mobile(mm)( &
initcons_mobile(k)(-c1) + &
initcons_stationary(k)(c1-exp(-dtrtau_rel(mm))))
and if rtau_rel is zero,
dmobile(k) = mobile(mm)( &
initcons_mobile(k)(exp(-dt*rtau_ret( mm))-c1) + &
initcons_stationary(k)
??
Does that look correct?

[eclare108213]

If rtau is zero, then exp(dt*rtau)=1, and the expressions that have 1-exp() should be zero. So if rtau_ret=0, then the mobile term is zero and we use the stationary term, and vice versa if rtau_rel=0.
@njeffery please confirm.

[apcraig]

OK. There is another issue. I have found another case where rtau_ret is 1.0 and ts = 3600. The exp is under-flowing in that condition and then we'd be left with a "0" with the "exp" value so exp-1 or 1-exp would just be 1 or -1. Mostly, rtau_ret and rtau_rel are about 1e-4 which matches up OK with ts=3600. I need a little more help here about what to do.

There seem to be cases where rtau_* is zero and other cases where exp(-rtauts) truly is underflowing on large negative rtauts. Maybe there are other cases too. Can rtau*ts ever be negative?

[eclare108213]

Since rtau_ret and rtau_rel are both frequencies in units of 1/s, I don't think they should ever be negative. They are both initialized to 1 in icedrv_init_column.F90, but they should not stay that way unless bgc_tracer_type < 0. The code comments indicate that this is actually the default condition, so this will always produce underflows in the exponential. Now I am beginning to understand this better -- I think the intent is that the exponential should go to zero, and my comments above are not really the problem.

Is the logic is incomplete?
I would recode this section, but I'm not sure if will make a difference.

  !-----------------------------------------------------------------
  !   time constants for mobile/stationary phase changes
  !-----------------------------------------------------------------
   
     if (m .ne. nlt_bgc_N(1)) then  
        if (hin_old  > hin) then  !melting
           rtau_rel(m) = c1/tau_rel(m)
           rtau_ret(m) = c0
        else                              !not melting
           rtau_ret(m) = c1/tau_ret(m)
           rtau_rel(m) = c0
        endif  
     elseif (tr_bgc_N .and. hin_old > hin + algal_vel*dt) then
           rtau_rel(m) = c1/tau_rel(m)
           rtau_ret(m) = c0
     elseif (tr_bgc_N) then
           rtau_ret(m) = c1/tau_ret(m)
           rtau_rel(m) = c0
     endif

simplified logic:

  !-----------------------------------------------------------------
  !   time constants for mobile/stationary phase changes
  !-----------------------------------------------------------------
  
     rtau_ret(m) = c1/tau_ret(m)
     rtau_rel(m) = c0
     if ((m /= nlt_bgc_N(1)) .and.  (hin_old  > hin))  &  ! melting
         .or. (tr_bgc_N .and. hin_old > hin + algal_vel*dt)  then
           rtau_rel(m) = c1/tau_rel(m)
           rtau_ret(m) = c0
     endif

But again, I think that the exponentials should be calculated rather than the timescales. Does this capture the various cases? tau_ret and tau_rel are set using tau_min/max from namelist.

  !-----------------------------------------------------------------
  !   time constants for mobile/stationary phase changes
  !-----------------------------------------------------------------
  
  expt = c0
  expl = c0
  if (tau_ret(m) > c1) expt = exp(-dt/tau_ret(m)) ! assume tau_ret is O(dt) or larger
  if (tau_rel(m) > c1) expl = exp(-dt/tau_rel(m)) ! assume tau_rel is O(dt) or larger
     if ((m /= nlt_bgc_N(1)) .and.  (hin_old  > hin))  &  ! melting
         .or. (tr_bgc_N .and. hin_old > hin + algal_vel*dt)  then
         expt = c1
     else
         expl = c1
     endif

           dmobile(k) = mobile(mm)*(initcons_mobile(k)*(expt-c1) + &
                             initcons_stationary(k)*(c1-expl))

[apcraig]

So, I sort of did the same thing yesterday, following the argmax of some other implementations, but I think what you propose is better. Is the tau_ret>c1 the right logic? was c1 chosen because it was the default value? the default should be a negative number, say -99 so it can be flagged more clearly. And then that logic should be >c0. but it looks like the default is set outside the columnphysics, so how can we be sure any default is set, and that it will be 1. What if it's set to 2 by default by the driver? There is a bunch of dicey stuff here.

Maybe we need something like

  argmax = 10.
  expt = c0
  expl = c0
  if (tau_ret(m) > c0) then
    exparg = min(dt/tau_ret(m),argmax)
    expt = exp(-exparg)
  endif
  if (tau_rel(m) > c0) then
    exparg = min(dt/tau_rel(m),argmax)
    expl = exp(-exparg)
  endif

     if ((m /= nlt_bgc_N(1)) .and.  (hin_old  > hin))  &  ! melting
         .or. (tr_bgc_N .and. hin_old > hin + algal_vel*dt)  then
         expt = c1
     else
         expl = c1
     endif

           dmobile(k) = mobile(mm)*(initcons_mobile(k)*(expt-c1) + &
                             initcons_stationary(k)*(c1-expl))

Also, one final thought. e-10 is only about 5e-5. Why is argmax 10? Why not 20 (e-20=2e-9) or even a more. 10 is smaller than it could be, you are much greater than needed to protect the underflow. You can probably go to a couple hundred with double precision and still not underflow. That's a separate question though.

[eclare108213]

@njeffery We've merged this change but I'd still appreciate it if you'd check the logic in ice_algae.F90 (as above) to make sure it's doing what you intend. Thx!

[eclare108213]

Regarding this question:

Also, one final thought. e-10 is only about 5e-5. Why is argmax 10? Why not 20 (e-20=2e-9) or even a more. 10 is smaller than it could be, you are much greater than needed to protect the underflow. You can probably go to a couple hundred with double precision and still not underflow. That's a separate question though.

e-10 is a value that was (I think) chosen by Bruce Briegleb when he was implementing the delta-Eddington radiation. I'm not sure why he chose that value, but my impression was that it was as big as he was comfortable with given the data that he and Bonnie had to build the radiation scheme. So it's specific to the radiation. In my opinion, we could increase the maximum (negative) argument by an order of magnitude (to 100). We also don't have to use the same value everywhere in the code, though I prefer consistency for stuff like that.