Skip to content

Commit 02ee145

Browse files
glemieuxglemieux
authored
Merge pull request #1277 from rgknox/two-stream-singularity-fix
comprehensive correction for singularity in two-stream
2 parents 0802e4b + e7b332f commit 02ee145

File tree

1 file changed

+59
-27
lines changed

1 file changed

+59
-27
lines changed

‎radiation/TwoStreamMLPEMod.F90‎

Lines changed: 59 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,8 @@ Module TwoStreamMLPEMod
3939
integer, parameter :: twostr_vis = 1 ! Named index of visible shortwave radiation
4040
integer, parameter :: twostr_nir = 2 ! Named index for near infrared shortwave radiation
4141

42-
42+
integer, parameter :: max_bands = 2 ! maximum number of bands (for scratch space)
43+
4344
! Allowable error, as a fraction of total incident for total canopy
4445
! radiation balance checks
4546

@@ -893,26 +894,31 @@ subroutine ZenithPrep(this,cosz_in)
893894
! notably the scattering coefficient "om".
894895

895896
class(twostream_type) :: this
896-
integer :: ib ! band index, matches indexing of rad_params
897+
897898
real(r8),intent(in) :: cosz_in ! Un-protected cosine of the zenith angle
898899

899-
real(r8) :: cosz ! the near-zero protected cosz
900-
integer :: ican ! scattering element canopy layer index (top down)
901-
integer :: icol ! scattering element column
902-
real(r8) :: asu ! single scattering albedo
900+
real(r8) :: cosz ! the near-zero protected cosz
901+
integer :: ican ! scattering element canopy layer index (top down)
902+
integer :: icol ! scattering element column
903+
integer :: ib ! band index, matches indexing of rad_params
904+
integer :: ib2 ! band inner loop index while testing for singularity
905+
real(r8) :: asu ! single scattering albedo
903906
real(r8) :: gdir
904907
real(r8) :: tmp0,tmp1,tmp2
905-
real(r8) :: betab_veg ! beam backscatter for vegetation (no snow)
906-
real(r8) :: betab_om ! multiplication of beam backscatter and reflectance
907-
real(r8) :: om_veg ! scattering coefficient for vegetation (no snow)
908-
real(r8) :: Kb_sing ! the KB_leaf that would generate a singularity
909-
! with the scelb%a parameter
910-
real(r8) :: Kb_stem ! actual optical depth of stem with not planar geometry effects
911-
! usually the base value
908+
real(r8) :: betab_veg ! beam backscatter for vegetation (no snow)
909+
real(r8) :: betab_om ! multiplication of beam backscatter and reflectance
910+
real(r8) :: om_veg ! scattering coefficient for vegetation (no snow)
911+
real(r8) :: Kb_sing(max_bands) ! the KB_leaf that would generate a singularity
912+
! with the scelb%a parameter
913+
real(r8) :: Kb_stem ! actual optical depth of stem with not planar geometry effects
914+
! usually the base value
912915
real(r8), parameter :: Kb_stem_base = 1.0_r8
913916
real(r8), parameter :: sing_tol = 0.01_r8 ! allowable difference between
914-
! the Kb_leaf that creates
915-
! a singularity and the actual
917+
! the Kb_leaf that creates
918+
! a singularity and the actual
919+
logical :: is_sing ! use this to control if we are actively trying to remove a singularity
920+
integer :: iter_sing ! iterator check to ensure we don't try to fix a singularity indefinitely
921+
real(r8) :: Kb_eff ! When testing for singularity, this is either the stem or stem and leaf optical depth
916922

917923
if( (cosz_in-1.0) > nearzero ) then
918924
write(log_unit,*)"The cosine of the zenith angle cannot exceed 1"
@@ -950,6 +956,7 @@ subroutine ZenithPrep(this,cosz_in)
950956
scelg%Kb_leaf = min(kb_max,rad_params%clumping_index(ft) * gdir / cosz)
951957

952958
! To avoid singularities, we need to make sure that Kb =/ a
959+
! for any of the bands...
953960
! If they are too similar, it will create a very large
954961
! term in the linear solution and generate solution errors
955962
! Lets identify the Kb_leaf that gives a singularity.
@@ -960,22 +967,47 @@ subroutine ZenithPrep(this,cosz_in)
960967
! (a*(lai+sai) - sai*kb_stem)/lai = Kb_sing
961968
! or.. adjust stem Kb?
962969
! (a*(lai+sai) - lai*kb_leaf)/sai = kb_stem_sing
970+
971+
if(scelg%lai>nearzero) then
972+
Kb_eff = scelg%Kb_leaf
973+
else
974+
Kb_eff = Kb_stem
975+
end if
976+
977+
! Assume there is a singularity so that we test for it
978+
is_sing = .true.
979+
iter_sing = 0
980+
981+
! Compute the singularity for all bands
982+
do ib = 1,this%n_bands
983+
Kb_sing(ib) = this%band(ib)%scelb(ican,icol)%a
984+
if (scelg%lai>nearzero) then
985+
Kb_sing(ib) = (Kb_sing(ib) * (scelg%lai+scelg%sai) - scelg%sai*Kb_stem)/scelg%lai
986+
end if
987+
end do
988+
989+
do_test_sing: do while(is_sing)
990+
! Now that we have commited to testing it, assume the solution works
991+
is_sing = .false.
992+
iter_sing = iter_sing
993+
if(iter_sing==10)then
994+
write(log_unit,*) 'error trying to remove singularity',iter_sing,scelg%Kb_leaf,Kb_stem,Kb_sing(:)
995+
call endrun(msg=errMsg(sourcefile, __LINE__))
996+
end if
997+
! Test to see if there is a singularity and make corrections if needed
998+
if (any((abs(Kb_sing(:) - Kb_eff)) < sing_tol)) then
999+
Kb_eff = Kb_eff + sing_tol
1000+
is_sing = .true.
1001+
end if
1002+
end do do_test_sing
1003+
9631004
if(scelg%lai>nearzero) then
964-
do ib = 1,this%n_bands
965-
Kb_sing = (this%band(ib)%scelb(ican,icol)%a*(scelg%lai+scelg%sai) - scelg%sai*Kb_stem)/scelg%lai
966-
if(abs(scelg%Kb_leaf - Kb_sing)<sing_tol)then
967-
scelg%Kb_leaf = Kb_sing + sing_tol
968-
end if
969-
end do
1005+
scelg%Kb_leaf = Kb_eff
9701006
else
971-
do ib = 1,this%n_bands
972-
Kb_sing = this%band(ib)%scelb(ican,icol)%a
973-
if(abs(Kb_stem - Kb_sing)<sing_tol)then
974-
Kb_stem = Kb_sing + sing_tol
975-
end if
976-
end do
1007+
Kb_stem = Kb_eff
9771008
end if
9781009

1010+
9791011
! RGK: My sense is that snow should be adding optical depth
9801012
! but we don't have any precedent for that in the FATES
9811013
! code or old CLM. Re-view this.

0 commit comments

Comments
 (0)