100% allocation to reproduction in shrubs causes crash in FATES-CNP

[JessicaNeedham]

I’m getting a fail in global CNP simulations when DEBUG=TRUE due to a divide by zero error caused by 100% allocation to reproduction in shrubs.

https://github.com/NGEET/fates/blob/1be3962933351b94fef134d607ac7d821d884f96/parteh/PRTAllometricCNPMod.F90#L2544

The line above causes a ‘forrtl: error (65): floating invalid’ error which I guess is because repro_c_frac is 1.

If we aren’t using the TRS scheme, and the cohort has a dbh above their dbh threshold for increased allocation to reproduction, then repro_c_frac is defined here:

https://github.com/NGEET/fates/blob/1be3962933351b94fef134d607ac7d821d884f96/parteh/PRTAllometricCNPMod.F90#L2477

Looking at the default parameter file, shrubs would have 0.1 and 0.9 which triggers the error above.

https://github.com/NGEET/fates/blob/1be3962933351b94fef134d607ac7d821d884f96/parameter_files/fates_params_default.cdl#L1531

https://github.com/NGEET/fates/blob/1be3962933351b94fef134d607ac7d821d884f96/parameter_files/fates_params_default.cdl#L1534

This is only triggered when DEBUG=TRUE. When debug is false it keeps running, but I’m not sure what the behaviour is. This was with the intel compiler.

[JessicaNeedham]

In the case where repro_c_frac is 1 we can wrap it in an if statement and make repro_w = repro_c_frac. I'm not sure about the alternative though, when repro_c_frac is between 0 and 1. Currently, repro_w = 0 because total_w is 0 which doesn't seem right. Should repro_w = repro_c_frac in all cases?

[mpaiao]

I think the code is making repro_w = repro_c_frac, but not in the most straightforward way because it is increasing it by (1-repro_c_frac then dividing the ratios at the end. Perhaps for readability we could do something like this, although maybe more verbose and less efficient...

     ! Estimate the weights of each organ
     if (state_mask(leaf_id)) then
        leaf_w  = target_dcdd(leaf_organ)
     else
        leaf_w  = 0._r8
     end if
     if (state_mask(fnrt_id)) then
        fnrt_w  = target_dcdd(fnrt_organ)
     else
        fnrt_w  = 0._r8
     end if
     if (state_mask(sapw_id)) then
        sapw_w  = target_dcdd(sapw_organ)
     else
        sapw_w  = 0._r8
     end if
     if(state_mask(struct_id)) then
        struct_w = target_dcdd(struct_organ)
     else
        struct_w = 0._r8
     end if
     if(state_mask(store_id)) then
        store_w  = target_dcdd(store_organ)
        store_nc = this%GetNutrientTarget(nitrogen_element,store_organ,stoich_growth_min) / target_c(store_organ)
        store_pc = this%GetNutrientTarget(phosphorus_element,store_organ,stoich_growth_min) / target_c(store_organ)
     else
        store_w  = 0._r8
        store_nc = 0._r8
        store_pc = 0._r8
     end if

     ! Find the weighted averages
     total_w = leaf_w + fnrt_w + sapw_w + struct_w + store_w
     if (total_w < nearzero) then
        ! Weight for all other tissues is zero, assume everything goes to reproduction.
        ! Or crash if this should never happen.
        avg_nc = repro_nc
        avg_pc = repro_pc
     else
        ! Scale the weights for all the tissues and reproduction, so they always add up to 1.
        if (state_mask(repro_id)) then
           leaf_w   = leaf_w   / total_w * (1._r8 - repro_c_frac)
           fnrt_w   = fnrt_w   / total_w * (1._r8 - repro_c_frac)
           sapw_w   = sapw_w   / total_w * (1._r8 - repro_c_frac)
           struct_w = struct_w / total_w * (1._r8 - repro_c_frac)
           repro_w  = repro_c_frac
        else
           leaf_w   = leaf_w   / total_w
           fnrt_w   = fnrt_w   / total_w
           sapw_w   = sapw_w   / total_w
           struct_w = struct_w / total_w
           repro_w  = 0._r8
        end if

        ! Find averages N:C and P:C ratios
        avg_nc =                                                                                &
             leaf_w   * prt_params%nitr_stoich_p1(ipft,prt_params%organ_param_id(leaf_organ  )) &
           + fnrt_w   * prt_params%nitr_stoich_p1(ipft,prt_params%organ_param_id(fnrt_organ  )) &
           + sapw_w   * prt_params%nitr_stoich_p1(ipft,prt_params%organ_param_id(sapw_organ  )) &
           + struct_w * prt_params%nitr_stoich_p1(ipft,prt_params%organ_param_id(struct_organ)) &
           + store_w  * store_nc + repro_w * repro_nc
        avg_pc =                                                                                &
             leaf_w   * prt_params%phos_stoich_p1(ipft,prt_params%organ_param_id(leaf_organ  )) &
           + fnrt_w   * prt_params%phos_stoich_p1(ipft,prt_params%organ_param_id(fnrt_organ  )) &
           + sapw_w   * prt_params%phos_stoich_p1(ipft,prt_params%organ_param_id(sapw_organ  )) &
           + struct_w * prt_params%phos_stoich_p1(ipft,prt_params%organ_param_id(struct_organ)) &
           + store_w  * store_pc + repro_w * repro_pc
     end if

[JessicaNeedham]

Thanks Marcos, that looks like it might be worth doing at some point. For now I'm I've just wrapped the case where repro_c_frac is 1 in an if statement. I'd been misunderstanding state_mask, but it is a more trivial fix than I thought.