collimated transmitted light problem

[xnwang]

In the abledo_surface.f90 in trunk, the subroutine multilvel_albedo calculates the proportions of lights distributed in several directions. Among them, the variable collim_trans_uncoll represents the collimated uncollided transmitted lights.

1) collim_trans_uncoll can be < 0
This is seen when cos_sun_angle (cosine of solar zenith angle) is too small, for example < 0.015. It can also happen when LAIeff = 0.

2) collim_trans_uncoll can have an important portion of data with very small positive values
Among the postive values, we can have about 15 % < 0.005, 7 % of postive data < 0.0005 for example. This variable collim_trans_uncoll is used in the computation of Collim_Abs_Tot (that is., the absorption of light by canopy). The latter can influence the vegetation growth in ORCHIDEE through the photosynthesis. The impact of small values of collim_trans_uncoll to the model may need some investigation, which help understand the sensitivity of other outputs to this variable.

[xnwang]

In the abledo_surface.f90 in trunk, the subroutine multilvel_albedo calculates the proportions of lights distributed in several directions. Among them, the variable collim_trans_uncoll represents the collimated uncollided transmitted lights, and collim_trans_coll represents the collimated collided transmitted lights.

1) collim_trans_coll can be < 0: often corresponds to collim_trans_uncoll = 1
This can be seen when LAIeff is < 1 sometimes, but with different cos_sun_angle.

2) Among the negative values of collim_trans_coll, a large proportion (> 99%) is within the range of -e-12 and 0. a small proportion is -0.999. This is with pft3 as an example. Need to understand why collim_trans_coll can be negative, because it is not physical.

3) when collim_trans_uncoll = 0.99 or 1 in the case of pft3, check the corresponding collim_trans_coll:
about 10% of collim_trans_coll < 0. about 15% of collim_trans_coll have very small values < e-5 but > 0.

3) when collim_trans_uncoll = 0.99 or 1:
This corresponds to 34% of data (pft3). It can be about 43% for grassland/cropland pfts.

Work is ongoing to understand why collim_trans_uncoll = 0.99 or 1 and collim_trans_coll < 0.

[xnwang]

To check why Collim_Trans_Coll can be -0.9, Mandresy printed some variables in subroutines such as propagate_fluxes. In propagate_fluxes subroutine:

Tran_UnColl_Tot(1:nlevels_tot)=Tran_UnColl_Tot_temp(1:nlevels_tot)
Tran_Coll_Tot(1:nlevels_tot)=Tran_Coll_Tot(1:nlevels_tot) - Tran_UnColl_Tot(1:nlevels_tot)

The negative values can come from the second line. This happens when laieff < 1 often, cos_sun_angle <= 0.001, and the vertical level is near the bottom of canopy.

It was found that Tran_UnColl_Tot_temp = 1, et Tran_Coll_Tot <= 0.1 in all the negative testing cases, so this is why Tran_Coll_Tot becomes negative after the substraction.

1)Before iterations:

Tran_Coll_Tot = 0, Tran_UnColl_Tot_temp = 1, and collim_down_cs(nlevels_tot) = 1, collim_down_cs(other leves)=0, iso_down_cs, isotrop_up_cs = 0 in the case of collimated lights.
Tran_Coll_Tot = 0, Tran_UnColl_Tot_temp = 1, and collim_down_cs = 0, isotrop_down_cs(top level)=1, isotrop_down_cs(other level) = 0 , isotropic_up_cs = 0 in the case of isotropic lights
collim_down_ns, iso_down_ns and iso_up_ns = 0 before iteration, for both cases

2) If the condition 'IF(collim_down_cs(ilevel) .GT. zero)' is satisfied at a quite late iteration step, that is, istep >= 11, or is never satisfied:

then the condition 'IF(ilevel == nlevels_tot-istep+1 .AND. .NOT. lisotrop)' will be never satifified for the bottom level of canopy. Then Tran_Uncoll_Tot_temp would be always 1 after all the iterations stop.

3) After each iteration, lconverged is set to false if the MAXVAL of collim_down_cs, isotrop_down_cs, and isotrop_up_cs are very small, ie., < 0.0000000001, then the iteration stops.

4) After each step of iteration, we update Tran_Coll_Tot(1:nlevels_tot), by adding collim_down_ns(0:nlevels_tot-1) + isotrop_down_ns(0:nlevels_tot-1). The index here is vertical levels.

Note that collim_down_ns(0) at an iteration step is related to its value at the previous iteration step and Collim_Tran_UnColl_Unscaled(1). isotrop_down_ns(0) are related to several variables, ie., Collim_Tran_Coll_Unscaled, Isotrop_Tran_UnColl_Unscaled, Isotrop_Tran_Coll_Unscaled, isotrop_up_cs, and Isotrop_Alb_Coll_Unscaled.

In case that isotrop_down_ns(0) and collim_down_ns(0) are very small during all the iterations, Tran_Coll_Tot(1) , that is Tran_Coll_Tot at the bottom canopy, becomes also small after the iterations stops.

5) If the two reasons in 2) and 4) can be confirmed, then this would explain why Collim_Trans_Coll can be -0.9. Mandresy is working on the confirmation by printing more variables.

6) why collim_down_ns can be 0 ? this is because Collim_Tran_UnColl_Unscaled = 0, which is due to too small cosine_sun_angle, according to the equation below:
Collim_Tran_UnColl_Unscaled(ilevel)=exp( - 0.5_r_std*laieff_collim_pft(ilevel)/cosine_sun_angle)

7) how to modify the code to deal with Collim_Trans_Coll -0.9, which is not a physical value ? if Collim_Tran_UnColl_Unscaled =0, we could set Collim_Trans_Coll=0 simply ?

[mrasolon]

For the point 2), we checked and saw that the 'IF(collim_down_cs(ilevel) .GT. zero)' is never satisfied.
For the point 4), isotrop_down_ns(0) is very small and collim_down_ns(0)=0 almost always. Therefore Tran_Coll_Tot(1) is very small. And we confirm the point 5), the reasons above explain why Collim_Trans_Coll can be -0.9.

[xnwang]

About the effects of small values on important variables such as gpp, albedo and lai, the sensivitiy test were made by Mandresy and were discussed and reported to the group already.

By setting collim_trans_xx < power(10,-5) to 0 while doing the spinup simulation, we reported very small differences in the gpp, lai and albedo. This point can be closed.

[mmcgrath]

I looked into propagate_fluxes in src_sechiba/albedo_surface.f90 of r7862.

For point (2), I am surprised. Which call of propagate_fluxes are you referring two? The subroutine is called twice: once for collimated light, once for isotropic light. For the first case (line 2176), there should be values in collim_down_cs. For the second case (line 2195), there should be none.

If you look at Figure 1 in ​https://gmd.copernicus.org/preprints/gmd-2016-280/gmd-2016-280.pdf, collim_down_cs represents the single arrow pointing downwards on the left hand graph. isotrop_down_cs are the groups of three arrows pointing down.

The amount of incoming light in the case of direct light (line 2176) and diffuse light (line 2195) is always equal to 1.0 (this is set in the lines just before the calls). Based on (6), if the sun angle is low, no direct light penetrates the layer. Physically this makes some sense (light coming in horizontally will not penetrate). Or, more accurately, light entering the canopy with a very shallow angle will only strike the ground kilometers away, after which it will have passed through kilometers of canopy and therefore will certainly collide with something before it arrives at the ground (becoming isotropic_down_cs).

Tran_Coll_Tot includes istrop_down_ns and collim_down_ns. Tran_Uncoll_tot_temp is initiallized to 1.0, which means that it's assumed that all incoming light makes it to the canopy floor without colliding with anything. I believe this is wrong, and the assumption should be that no light makes it to the floor without colliding, based on how Tran_Uncoll_Tot_temp is calculated (it is only calculated if collim/isotrop_down_cs is greater than zero, i.e., if there is actually downwelling light in this level). If there is no downwelling light, there can be no uncollimated transmitted light, either.

Consequently, I propose to change the line

Tran_Uncoll_Tot_temp(:)=un

to

Tran_Uncoll_Tot_temp(:)=zero

If you can put in an ipslerr to catch when values of Tran_Coll_Tot < 0.0, that could be useful to create a test case. I am a little concerned that these variables are used for something, given the comment in the code:

! now separate the collided from the uncollided light. Notice that
! this is not really needed for any purposes other than debugging, as the
! important quantity is the total amount of light striking the ground.

[xnwang]

About point 1) ‘before iteration’, missing information is added to the previous reply. The first calling of propagate_fluxes was examined more carefully, because it is used to compute the collimated lights.

1) About the proposition of Tran_Uncoll_Tot_temp(:)=zero :

Indeed, using Tran_Uncoll_Tot_temp(:)=zero would not change the total collimated transmitted lights. But does this variable has real physical meaning ? If yes, is the initialisation of Tran_Uncoll_Tot_temp to 0 physical ?

Note that Tran_Uncoll_Tot_temp is updated later by Tran_Uncoll_Tot_temp(ilevel)=Tran_Uncoll_Tot_temp(ilevel+1)*Collim_Tran_UnColl_Unscaled(ilevel).

Later another computation: Tran_Coll_Tot(1:nlevels_tot)=Tran_Coll_Tot(1:nlevels_tot)-Tran_UnColl_Tot_temp(1:nlevels_tot).
So does this equation has any real and useful meaning in physics ? If it has, then it seem better to update its value in a correct way.

2) I am wondering if it is a good idea to do in propagate_fluxes :

Tran_UnColl_Tot(1:nlevels_tot)=Tran_UnColl_Tot_temp(1:nlevels_tot)

adding one line: if at some levels where Tran_Coll_Tot <Tran_UnColl_Tot , we would set Tran_UnColl_Tot = 0 at these levels. Then keep the following line.
Tran_Coll_Tot(1:nlevels_tot)=Tran_Coll_Tot(1:nlevels_tot)-Tran_UnColl_Tot(1:nlevels_tot)

[mrasolon]

Few changes of code were proposed and tested.

1) In propagate_fluxes, what used to be defined as Tran_Coll_Tot and initialized to 0 in the beginning of the subroutine was changed to Tran_Tot. For more clarity. Tran_Coll_Tot is solely computed at the end of the subroutine in this line :
Tran_Coll_Tot(1:nlevels_tot)=Tran_Tot(1:nlevels_tot)-Tran_UnColl_Tot(1:nlevels_tot)

2) a) Right before this last computation, few lines have been added :

Tran_UnColl_Tot(1:nlevels_tot)=Tran_UnColl_Tot_temp(1:nlevels_tot)

Added lines below :

Tran_Coll_Tot_temp(:)=Tran_Tot(:)-Tran_UnColl_Tot(:)
WHERE((Tran_Coll_Tot_temp<0) .AND. (Tran_UnColl_Tot .eq. 1))

Tran_UnColl_Tot=0

ENDWHERE

And then the last computation :

Tran_Coll_Tot(1:nlevels_tot)=Tran_Tot(1:nlevels_tot)-Tran_UnColl_Tot(1:nlevels_tot)

2)b) Results : by adding these lines, we observed that there were no more Collim_Tran_Coll < -0.9. And it was checked for pft 3 but also for pft 2 (which had also this non-physical issue for Collim_Tran_Coll.