Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Calculating the surface albedo; re-writing cable_albedo.F90 #266

Open
penguian opened this issue Sep 2, 2020 · 28 comments
Open

Calculating the surface albedo; re-writing cable_albedo.F90 #266

penguian opened this issue Sep 2, 2020 · 28 comments
Assignees
@JhanSrbinovsky
Labels
component: JAC readiness major new
Milestone
6. Report

Comments

@penguian
Copy link
Collaborator

penguian commented Sep 2, 2020

keyword_maygit owner:[email protected] type_model improvement | by [email protected]

In JAC, the surface albedo is required before CABLE has more generally been
initialized. Closer inspection shows that we can compute the surface albedo given
only a few veg and soil parameters - So long as we move to accessing
whether a cell is sunlit or not given the zenith angle rather than downward SW,
which in a circular fashion precluding calculation of the surface albedo has not
yet been determined, requiring surface albedo in its determination.
Closer inspection also reveals that this pathway to get the surface albedo is
very un-neccessarily convoluted and hard to follow with non-descriptive varible
names and a haphazzard order of operation.

(surface)Albedo is CALLed from cable_cbm, per timestep, and is in the file cable_albedo.F90 which contains 2 further subroutines that are only called from (surface)Albedo.

This file is split into 3 files.

    : cbl_albedo.F90

    : cbl_snow_albedo.F90

    : cbl_soilColour_albedo.F90

The "surface_albedo" routine in particular is re-written to carify what is actually
needed and when. See comment:4.

In comment:5 we discuss the beam fraction and beam fraction weighted albedo further.

In comment:6 we present the trunk@516aafe393e159a1fe2d36044f8f3d687d40687a version of the (surface-)albedo subroutine.

In comment:7 we begin our discussion of sequential elements in the trunk@516aafe393e159a1fe2d36044f8f3d687d40687a version of the (surface-)albedo subroutine.

This then dictated a re-write of cable_init_radiation following a
similar strategy as here and for analagous reasons. (see #267)

The below issues are addressed in Ticket #263

To calculate the surface albedo CABLE needds to compute the canopy
reflectivity for which it needs to know the effective LAI, i.e. the LAI that
may actually be seen by SW radiation, IF the canopy is covered in snow.
This calculation has been extracted from the model so that it can be
specifically called in isolation: cbl_lai_eff.F90

  • in contrast to the previous scenario where it was embedded in a subroutine
    which DOES require CABLE to be properly initialised - with forcing that does
    not even exist at the stage where surface albedo is required.

In order to calculate the effective LAI CABLE needs to know the height of the
canopy above snow. This calculation has been extracted from the model so that
it can be specifically called in isolation: cbl_hgtAbove_snow.F90

In order to calculate the effective LAI, clearly CABLE needs to know the height of the
canopy and the LAI of the PFT at that cell. This calculation has been extracted from
the model so that it can be specifically called in isolation: cbl_LAI_canopy_height.F90

In the JAC version of CABLE, we get the LAI and canopy height from UM I/O spatial maps. In the interface we limit these values. However it is not even clear any longer if this is actually needed any more.

Issue migrated from trac:266 at 2023-11-27 11:34:20 +1100
@penguian
Copy link
Collaborator Author

penguian commented Sep 2, 2020

@[email protected] edited the issue description

2 similar comments
@penguian
Copy link
Collaborator Author

@[email protected] edited the issue description

@penguian
Copy link
Collaborator Author

@[email protected] edited the issue description

@penguian
Copy link
Collaborator Author

@[email protected] commented

Rename variables to have readable, meaningful names

The extinction co-efficient to diffuse radiation, rad%extkd is renamed as ExtCoeff_dif.

The extinction co-efficient to direct radiation, rad%extkb is renamed as ExtCoeff_beam.

The effective extinction co-efficient to diffuse radiation, rad%extkdm is renamed as EffExtCoeff_dif.

The effective extinction co-efficient to direct radiation, rad%extkbm is renamed as EffExtCoeff_beam.

The canopy transmittance fraction to diffuse radiation, rad%cexpkdm is renamed as CanopyTransmit_dif.

The canopy transmittance fraction to beam radiation, rad%cexpkbm is renamed as CanopyTransmit_beam.

The canopy reflectance fraction to diffuse radiation, rad%rhocdf is renamed as CanopyTransmit_dif.

The canopy reflectance fraction to beam radiation, rad%rhocbm is renamed as CanopyTransmit_beam.

The albedo of the surface as seen by the atmosphere due to both the soil and canopy reflection to diffuse radiation, rad%reffdf is renamed as EffSurfRefl_dif.

The albedo of the surface as effectively seen by the atmosphere due to both the soil and canopy reflection to beam radiation, rad%reffbm is renamed as EffSurfRefl_beam.

The beam fraction of downward SW (in both bands [VIS/NIR]), rad%fbeam is renamed as RadFbeam.

The surface albedo (EffSurfRefl_*) weighted by beam fraction **rad%albedo **, is renamed as RadAlbedo.

The reduced LAI due to snow cover canopy%vlaiw, is renamed as reducedLAIdue2snow.

See comment:5 for further discussion of RadAlbedo and RadFbeam.

@penguian
Copy link
Collaborator Author

@[email protected] changed _comment0 which not transferred by tractive

@penguian
Copy link
Collaborator Author

@[email protected] changed _comment1 which not transferred by tractive

@penguian
Copy link
Collaborator Author

@[email protected] commented

Note: RadFbeam and RadAlbedo

The JAC implementation obviously is based on the same interface as ACCESS-CM2. We dodge calculating beam fraction on the Albedo (surf_couple_radiation) pathway for two reasons:

  1. The forced model does not use the albedo in calculation of downwards SW anyway.

  2. Calculation of beam fraction via the spitter function requires downwards SW which is not readily available at this point in the JULES framework (they have no need for it AND in a circular fashion the albedo is being calculated to determine the downwards SW)

The forced model only uses the calculated albedo as a diagnostic.
Finally, for comparisson purposes, from memory JULES outputs an albedo analgous to the un-weighted surface albedo (EffSurfRefl_*). CABLE offline outputs rad%albedo.

@penguian
Copy link
Collaborator Author

@[email protected] edited the issue description

@penguian
Copy link
Collaborator Author

@[email protected] commented

Trunk@516aafe393e159a1fe2d36044f8f3d687d40687a version of (surface_)albedo subroutine

This is a verbatim copy of the code existing in the trunk for this subroutine following the call to surface_albedosn which shall be discussed elsewhere.

CALL calc_rhoch( c1,rhoch, mp, nrb, VegTaul, VegRefl )

    ! Update extinction coefficients and fractional transmittance for
    ! leaf transmittance and reflection (ie. NOT black leaves):
    !---1 # visible, 2 nir radiaition
    DO b = 1, 2

       EffExtCoeff_dif(:,b) = ExtCoeff_dif * c1(:,b)

       !--Define canopy diffuse transmittance (fraction):
       CanopyTransmit_dif(:,b) = EXP(-EffExtCoeff_dif(:,b) * reducedLAIdue2snow)

       !---Calculate effective diffuse reflectance (fraction):
       WHERE( veg_mask )                                             &
            EffSurfRefl_dif(:,b) = CanopyRefl_dif(:,b) + (AlbSnow(:,b)             &
            - CanopyRefl_dif(:,b)) * CanopyTransmit_dif(:,b)**2

       !---where vegetated and sunlit
       WHERE (sunlit_veg_mask)

          EffExtCoeff_beam(:,b) = ExtCoeff_beam * c1(:,b)

          ! Canopy reflection (6.21) beam:
          CanopyRefl_beam(:,b) = 2. * ExtCoeff_beam / ( ExtCoeff_beam + ExtCoeff_dif )          &
               * rhoch(:,b)

          ! Canopy beam transmittance (fraction):
          dummy2 = MIN(EffExtCoeff_beam(:,b)*reducedLAIdue2snow, 20.)
          dummy  = EXP(-dummy2)
          CanopyTransmit_beam(:,b) = REAL(dummy)

          ! Calculate effective beam reflectance (fraction):
          EffSurfRefl_beam(:,b) = CanopyRefl_beam(:,b) + (AlbSnow(:,b)             &
               - CanopyRefl_beam(:,b))*CanopyTransmit_beam(:,b)**2

       END WHERE

       ! Define albedo:
       WHERE( veg_mask )                                      &
            RadAlbedo(:,b) = ( 1. - RadFbeam(:,b) )*EffSurfRefl_dif(:,b) +           &
            RadFbeam(:,b) * EffSurfRefl_beam(:,b)

    END DO

@penguian
Copy link
Collaborator Author

@[email protected] changed _comment0 which not transferred by tractive

@penguian
Copy link
Collaborator Author

@[email protected] edited the issue description

@penguian
Copy link
Collaborator Author

@[email protected] commented

Sequemntial discussion of elements of Trunk@516aafe393e159a1fe2d36044f8f3d687d40687a version of (surface_)albedo subroutine

CALL calc_rhoch( c1,rhoch, mp, nrb, VegTaul, VegRefl )

c1 and rhoch are additional co-efficients used in conjunction with extinction co-efficients that require definition of per PFT Refl and Tau in vegetation parameters (namelist or file etc).

The standout point to make here is that this is flagrantly wrong. This is the 2nd time this is called deeming it inefficient at best AND the same subroutine appears in two files (cable_radiation being the other) which is ....less than good.

The ONLY reason in support of this is that the infrastructure does not exist to call the one instance of calc_rhoch(), once only.

We rectify this by eliminating both instances of calc_rhoch() and creating a new file cbl_rhoch.F90, containing only calc_rhoch(). The variables required to be returned by calc_rhoch(), c1 and rhoch, ar the declared at a higher level and passed both to the init_radiation() subroutine where calc_rhoch() willl be CALLed and the to albedo(), where they will be used, again.

@penguian
Copy link
Collaborator Author

@[email protected] changed _comment0 which not transferred by tractive

@penguian
Copy link
Collaborator Author

@[email protected] changed _comment1 which not transferred by tractive

@penguian
Copy link
Collaborator Author

@[email protected] changed _comment2 which not transferred by tractive

@penguian
Copy link
Collaborator Author

@[email protected] edited the issue description

@penguian
Copy link
Collaborator Author

@[email protected] commented

The remainder of the calculations exist inside a do loop over the radiation bands (VIS/NIR)

    DO b = 1, 2
....
....
....
       ! Define albedo:
       WHERE( veg_mask )                                      &
            RadAlbedo(:,b) = ( 1. - RadFbeam(:,b) )*EffSurfRefl_dif(:,b) +           &
            RadFbeam(:,b) * EffSurfRefl_beam(:,b)

    END DO

Where here we have shown the very last thing computed within that loop, before leaving the subroutine.

As discussed above (comment:5), neither of these is computed in the radiation/albedo pathway in JAC and it seems that RadAlbedo is merely a diagnostic anyway, not feeding into the model at all.

@penguian
Copy link
Collaborator Author

@[email protected] changed _comment0 which not transferred by tractive

@penguian
Copy link
Collaborator Author

@[email protected] changed _comment1 which not transferred by tractive

@penguian
Copy link
Collaborator Author

@[email protected] edited the issue description

@penguian
Copy link
Collaborator Author

@[email protected] commented

Following comment:7, within the DO loop (and omitting calc in comment:8) and visually separating the diffuse light calcs from the direct beam calcs (and removing Fortran !comments):

       EffExtCoeff_dif(:,b) = ExtCoeff_dif * c1(:,b)

       CanopyTransmit_dif(:,b) = EXP(-EffExtCoeff_dif(:,b) * reducedLAIdue2snow)

       WHERE( veg_mask )                                             &
            EffSurfRefl_dif(:,b) = CanopyRefl_dif(:,b) + (AlbSnow(:,b)             &
            - CanopyRefl_dif(:,b)) * CanopyTransmit_dif(:,b)**2

       WHERE (sunlit_veg_mask)

          EffExtCoeff_beam(:,b) = ExtCoeff_beam * c1(:,b)

          CanopyRefl_beam(:,b) = 2. * ExtCoeff_beam / ( ExtCoeff_beam + ExtCoeff_dif )          &
               * rhoch(:,b)

          dummy2 = MIN(EffExtCoeff_beam(:,b)*reducedLAIdue2snow, 20.)
          dummy  = EXP(-dummy2)
          CanopyTransmit_beam(:,b) = REAL(dummy)

          EffSurfRefl_beam(:,b) = CanopyRefl_beam(:,b) + (AlbSnow(:,b)             &
               - CanopyRefl_beam(:,b))*CanopyTransmit_beam(:,b)**2

       END WHERE

In summary, using Pseudo code and assuming that we do this over radiation bands, and for diffuse and direct radiation:

Calc  Effective Exinction co-efficients

Calc  Canopy Transmitance

Calc  Canopy Reflectance

Calc  **Effective Surface Reflectance** - this being what we actually need to pass back to th atmosphere/radiation scheme

Note: Currently the canopy reflectance for the diffuse component is mysteriously calculated in init_radiation(). It is only used here and indeed belongs here. Thus it is moved here.

@penguian
Copy link
Collaborator Author

@[email protected] changed _comment0 which not transferred by tractive

@penguian
Copy link
Collaborator Author

@[email protected] changed _comment1 which not transferred by tractive

@penguian
Copy link
Collaborator Author

@[email protected] changed _comment2 which not transferred by tractive

@penguian
Copy link
Collaborator Author

@[email protected] changed _comment3 which not transferred by tractive

@penguian
Copy link
Collaborator Author

@[email protected] commented

Change occurred in the 5th decimal place of sumbal following eliminating strict definition of variable dummy as KIND=r_2, used in calc of canopy transmittance. This is because we wan't to get rid of r_2. No change detectable in fluxes - as discussed before sumbal is way way way more sensitive to any change than are plotted fluxes etc. @22e293bd67d34b51639d9645fa05023a0e4d55b5

@penguian
Copy link
Collaborator Author

penguian commented Feb 2, 2023

@[email protected] commented

Should we close this? We need to put the list of renamed variables in the new CABLE documentation instead of a ticket/issue.

@penguian
Copy link
Collaborator Author

penguian commented Feb 7, 2023

@[email protected] set keywords to maygit

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@JhanSrbinovsky JhanSrbinovsky

Labels
component: JAC readiness major new
Projects
None yet
Milestone
6. Report
Development

No branches or pull requests
2 participants
@penguian @JhanSrbinovsky