search

CABLE-LSM / CABLE-Trac-archive

Archive CABLE Trac contents as issues
Other
0 stars 0 forks source link

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

Open penguian opened 3 years ago

penguian commented 3 years ago

keyword_maygit owner:jxs599@nci.org.au type_model improvement | by srb001@csiro.au

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 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 commented 3 years ago

@jxs599@nci.org.au edited the issue description

penguian commented 3 years ago

@jxs599@nci.org.au edited the issue description

penguian commented 3 years ago

@jxs599@nci.org.au edited the issue description

penguian commented 3 years ago

@jxs599@nci.org.au 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 commented 3 years ago

@jxs599@nci.org.au changed _comment0 which not transferred by tractive

penguian commented 3 years ago

@jxs599@nci.org.au changed _comment1 which not transferred by tractive

penguian commented 3 years ago

@jxs599@nci.org.au 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 commented 3 years ago

@jxs599@nci.org.au edited the issue description

penguian commented 3 years ago

@jxs599@nci.org.au 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 commented 3 years ago

@jxs599@nci.org.au changed _comment0 which not transferred by tractive

penguian commented 3 years ago

@jxs599@nci.org.au edited the issue description

penguian commented 3 years ago

@jxs599@nci.org.au 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 commented 3 years ago

@jxs599@nci.org.au changed _comment0 which not transferred by tractive

penguian commented 3 years ago

@jxs599@nci.org.au changed _comment1 which not transferred by tractive

penguian commented 3 years ago

@jxs599@nci.org.au changed _comment2 which not transferred by tractive

penguian commented 3 years ago

@jxs599@nci.org.au edited the issue description

penguian commented 3 years ago

@jxs599@nci.org.au 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 commented 3 years ago

@jxs599@nci.org.au changed _comment0 which not transferred by tractive

penguian commented 3 years ago

@jxs599@nci.org.au changed _comment1 which not transferred by tractive

penguian commented 3 years ago

@jxs599@nci.org.au edited the issue description

penguian commented 3 years ago

@jxs599@nci.org.au 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 commented 3 years ago

@jxs599@nci.org.au changed _comment0 which not transferred by tractive

penguian commented 3 years ago

@jxs599@nci.org.au changed _comment1 which not transferred by tractive

penguian commented 3 years ago

@jxs599@nci.org.au changed _comment2 which not transferred by tractive

penguian commented 3 years ago

@jxs599@nci.org.au changed _comment3 which not transferred by tractive

penguian commented 3 years ago

@jxs599@nci.org.au 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 commented 1 year ago

@ccc561@nci.org.au 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 commented 1 year ago

@ccc561@nci.org.au set keywords to maygit