search

CABLE-LSM / CABLE-Trac-archive

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

ESM1.5 Reconciliation - esm15_canopy #348

Closed penguian closed 1 year ago

penguian commented 1 year ago

keyword_maygit owner:jxs599@nci.org.au resolution_fixed type_defect | by yh4968

Issue migrated from trac:348 at 2023-11-27 11:43:22 +1100

penguian commented 1 year ago

@yh4968@nci.org.au uploaded file switches-canopy.docx (49.2 KiB)

penguian commented 1 year ago

@yh4968@nci.org.au commented

https://trac.nci.org.au/trac/cable/browser/branches/Share/ESM15-CABLE3_draft1/science/canopy/cable_canopy.F90

penguian commented 1 year ago

@ccc561@nci.org.au commented

      ! Soil sensible heat:
      IF (cable_runtime%esm15_canopy) THEN

        canopy%fhs = air%rho*CCAPP*(ssnow%tss - met%tvair) /ssnow%rtsoil
        canopy%ga = canopy%fns-canopy%fhs-canopy%fes

      ELSE !(cable_runtime%esm15_canopy) THEN

        IF (cable_user%gw_model .OR. cable_user%or_evap) THEN
          canopy%fhs = air%rho*CCAPP*(ssnow%tss - met%tvair) / &
                       (ssnow%rtsoil + REAL(ssnow%rt_qh_sublayer))

        ELSEIF (cable_user%litter) THEN

          !! vh_js !! account for additional litter resistance to sensible heat transfer
          !! INH simplifying code using rhlitt
          canopy%fhs = air%rho*CCAPP*(ssnow%tss - met%tvair) / &
                       !(ssnow%rtsoil +  real((1-ssnow%isflag))*veg%clitt*0.003/canopy%kthLitt/(air%rho*CCAPP))
                       (ssnow%rtsoil + rhlitt)
        ELSE
          canopy%fhs = air%rho*CCAPP*(ssnow%tss - met%tvair) /ssnow%rtsoil
        ENDIF !IF (cable_user%gw_model .OR. cable_user%or_evap) THEN

        !! Ticket #90 ssnow%cls factor should be retained: required for energy balance
        !! INH: %cls factor included in %fes already - do not include here
        canopy%ga = canopy%fns-canopy%fhs-canopy%fes !*ssnow%cls

      ENDIF !(cable_runtime%esm15_canopy) THEN

Only added code for additional options, but no change to the case used by ESM1.5.

To accept.

penguian commented 1 year ago

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

penguian commented 1 year ago

@ccc561@nci.org.au commented

      IF (.NOT. cable_runtime%esm15_canopy) THEN
        write(6,*) "SLI is not an option right now"
          ! SLI SEB to get canopy%fhs, canopy%fess, canopy%ga
          ! (Based on old Tsoil, new canopy%tv, new canopy%fns)
          !H!CALL sli_main(1,dels,veg,soil,ssnow,met,canopy,air,rad,1)
      ENDIF !(cable_runtime%esm15_canopy) THEN

Only a WRITE statement added.

To accept.

penguian commented 1 year ago

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

penguian commented 1 year ago

@ccc561@nci.org.au commented

   met%dva = (qstvair - met%qvair) *  Crmair/Crmh2o * met%pmb * 100.0
   dsx = met%dva     ! init. leaf surface vpd
  IF( .NOT. cable_runtime%esm15_canopy ) THEN
    dsx= MAX(dsx,0.0)
  END IF   
      qstvair,       & ! sat spec humidity at leaf temperature

There is the same change in cbl_dryLeaf.F90 (#350) between ESM and trunk when dsx is updated:

            ! Update leaf surface vapour pressure deficit:
            dsx(i) = met%dva(i) + air%dsatdk(i) * (tlfx(i)-met%tvair(i))

            IF( .NOT. cable_runtime%esm15_dryLeaf ) THEN
              dsx(i)=  MAX(dsx(i),0.0)
            END IF
penguian commented 1 year ago

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

penguian commented 1 year ago

@ccc561@nci.org.au commented

  IF( .NOT. cable_runtime%esm15_canopy ) THEN
    !L_REV_CORR - initialise sensitivity/ACCESS correction terms
    !NB %fes_cor is NOT initialised to zero at this point
    canopy%fhs_cor = 0.0
    canopy%fns_cor = 0.0
    canopy%ga_cor = 0.0
    canopy%fes_cor = 0.0

    !L_REV_CORR - new working variables
    rttsoil = 0.
    rhlitt = 0.
    relitt = 0.
    alpm1  = 0.
    beta2  = 0.
  END IF

Should be harmless.

But it is a strange place to initialise the canopy% variables. They are not used in cable_canopy.F90 and seem to be given a value in cbl_soilsnow_main.F90 when running coupled to the UM.

penguian commented 1 year ago

@ccc561@nci.org.au commented

   canopy%zetar(:,1) = CZETA0 ! stability correction terms
   canopy%zetar(:,2) = CZETPOS + 1 
  IF( .NOT. cable_runtime%esm15_canopy ) THEN
    canopy%zetash(:,1) = CZETA0 ! stability correction terms
    canopy%zetash(:,2) = CZETPOS + 1
  END IF

No idea where it comes from

penguian commented 1 year ago

@ccc561@nci.org.au commented

            ! Forced convection boundary layer conductance                     
            ! (see Wang & Leuning 1998, AFM):
            IF( cable_runtime%esm15_canopy ) THEN
            gbhu(j,1) = gbvtop(j)*(1.0-EXP(-canopy%vlaiw(j)                    &
                        *(0.5*rough%coexp(j)+rad%extkb(j) ))) /                &
                        (rad%extkb(j)+0.5*rough%coexp(j))

            gbhu(j,2) = (2.0/rough%coexp(j))*gbvtop(j)*  &
                        (1.0-EXP(-0.5*rough%coexp(j)*canopy%vlaiw(j)))         &
                        - gbhu(j,1)
            ELSE
              !vh! inserted 'min' to avoid floating underflow
              gbhu(j,1) = gbvtop(j)*(1.0-EXP(-MIN(canopy%vlaiw(j)                    &
                              *(0.5*rough%coexp(j)+rad%extkb(j) ),20.0))) /            &
                              (rad%extkb(j)+0.5*rough%coexp(j))

              gbhu(j,2) = (2.0/rough%coexp(j))*gbvtop(j)*  &
                              (1.0-EXP(-MIN(0.5*rough%coexp(j)*canopy%vlaiw(j),20.0))) &
                              - gbhu(j,1)

            ENDIF

Same MIN() usage for the exponentials with the extinction coefficients as in #343 and some of #345. Should probably group with the other similar changes first.

penguian commented 1 year ago

@ccc561@nci.org.au commented

IF( cable_runtime%esm15_canopy ) THEN

  tv_denom_transd(j) = ( 1.0-rad%transd(j) )
  tv_denom(j) = 2.0*(tv_denom_transd(j)) * CSBOLTZ*CEMLEAF
  tv_frac(j) = rad%lwabv(j) / tv_denom(j)
  tv_tk(j) = met%tk(j) **4
  tv(j) = ( tv_frac(j) + tv_tk(j) )
  !effectivel MAX here does IF( |tv_frac| >  tv_tk ) tv_frac = -1.0 * tv_tk
  tv(j) = MAX( 0.0, tv(j) )
  tv(j) = tv(j) ** .25
  canopy%tv(j) = tv(j)

ELSE

             IF (  (rad%lwabv(j) / (2.0*(1.0-rad%transd(j))            &
                  * CSBOLTZ*CEMLEAF)+met%tvrad(j)**4) .GT. 0.0) THEN

                canopy%tv(j) = (rad%lwabv(j) / (2.0*(1.0-rad%transd(j))            &
                     * CSBOLTZ*CEMLEAF)+met%tvrad(j)**4)**0.25
             ELSE
                canopy%tv(j) = met%tvrad(j)
             ENDIF

END IF

It uses the same equation (several terms vs 1 line). But the behaviour if the equation is negative is different:

Note - potential bug in ESM1.5 if setting canopy%tv to absolute zero!! Perhaps condition isn't met very often.

penguian commented 1 year ago

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

penguian commented 1 year ago

@ccc561@nci.org.au commented

IF( cable_runtime%esm15_canopy ) THEN
           canopy%tv(j) = met%tk(j)
ELSE
  canopy%tv(j) = met%tvrad(j)
END IF

Same as in radiation (#346), need consistency in the choice for both.

penguian commented 1 year ago

@ccc561@nci.org.au commented

      ! Calculate soil sensible heat:
      IF (cable_runtime%esm15_canopy) THEN
        canopy%fhs = air%rho*CCAPP*(ssnow%tss - met%tk) /ssnow%rtsoil
      ELSE !(cable_runtime%esm15_canopy) THEN

        ! INH: I think this should be - met%tvair
        IF (cable_user%gw_model .OR. cable_user%or_evap) THEN
          canopy%fhs =  air%rho*CCAPP*(ssnow%tss - met%tk) / &
                     (ssnow%rtsoil + ssnow%rt_qh_sublayer)

          !note if or_evap and litter are true then litter resistance is
          !incluyded above in ssnow%rt_qh_sublayer
        ELSEIF (cable_user%litter) THEN
          !! vh_js !! account for additional litter resistance to sensible heat transfer
          !! INH simplifying code using rhlitt
           canopy%fhs = air%rho*CCAPP*(ssnow%tss - met%tk) / &
                        !(ssnow%rtsoil + real((1-ssnow%isflag))*veg%clitt*0.003/canopy%kthLitt/(air%rho*CCAPP))
                        (ssnow%rtsoil + rhlitt)
        ELSE
              canopy%fhs = air%rho*CCAPP*(ssnow%tss - met%tvair) /ssnow%rtsoil
        ENDIF !(cable_user%gw_model .OR. cable_user%or_evap) THEN

      ENDIF !(cable_runtime%esm15_canopy) THEN

Only additional code for new options. Accept trunk version.

penguian commented 1 year ago

@ccc561@nci.org.au commented

      IF (.NOT. cable_runtime%esm15_canopy) THEN

        write(6,*) "SLI is not an option right now"
        ! SLI SEB to get canopy%fhs, canopy%fess, canopy%ga
        ! (Based on old Tsoil, new canopy%tv, new canopy%fns)
            !H!CALL sli_main(1,dels,veg,soil,ssnow,met,canopy,air,rad,1)

      ENDIF !(cable_runtime%esm15_canopy) THEN

Only a WRITE statement added.

To accept.

penguian commented 1 year ago

@ccc561@nci.org.au commented

       !upwards flux density of water (kg/m2/s) - potential evapotranspiration 
       !radiation weighted soil and canopy contributions
       !Note: If PM routine corrected then match changes here
      IF ( cable_runtime%esm15_canopy) THEN
        canopy%epot = ((1.-rad%transd)*canopy%fevw_pot +                         &
                    rad%transd*ssnow%potev) * dels/air%rlam 
      ELSE 
       canopy%epot = ((1.-rad%transd)*canopy%fevw_pot +                         &
            rad%transd*ssnow%potev*ssnow%cls) * dels/air%rlam

       canopy%rniso = SUM(rad%rniso,2) + rad%qssabs + rad%transd*met%fld + &
            (1.0-rad%transd)*CEMLEAF* &
            CSBOLTZ*met%tvrad**4 - CEMSOIL*CSBOLTZ*met%tvrad**4

      END IF

For canopy%epot, this is likely linked to the %cls ticket (#137), should be in a flag with the #344 changes.

For canopy%rniso, this variable was added for cable_climate.F90, that's the only place it is used. To accept.

penguian commented 1 year ago

@ccc561@nci.org.au commented

          IF ( cable_runtime%esm15_canopy) THEN
                      r_sc(j) = term5(j) * LOG(zscl(j)/rough%z0soilsn(j)) *              &
                                ( EXP(2*CCSW*canopy%rghlai(j)) - term1(j) ) / term3(j)
          ELSE!( cable_runtime%esm15_canopy)
                       !Ticket #154
                       r_sc(j) = term5(j) * LOG(zscl(j)/rough%z0soilsn(j)) *              &
                            ( EXP(2*CCSW*canopy%rghlai(j)) - term2(j) ) / term3(j)
                       r_sc(j) = r_sc(j) + term5(j) * LOG(rough%disp(j)/rough%z0soilsn(j)) *  &
                            ( EXP(2*CCSW*canopy%rghlai(j)) - term1(j) ) / term3(j)
          END IF!( cable_runtime%esm15_canopy)

Bug fix from #154. Accept the trunk.

penguian commented 1 year ago

@ccc561@nci.org.au commented

            IF ( cable_runtime%esm15_canopy) THEN

              r_sc(j) = rough%rt0us(j) + rough%rt1usa(j) + rough%rt1usb(j) +     &
                      ( LOG( (zscl(j) - rough%disp(j)) /                       &
                      MAX( rough%zruffs(j)-rough%disp(j),                      &
                      rough%z0soilsn(j) ) ) - psis1( (zscl(j)-rough%disp(j))   &
                      / (rough%zref_tq(j)/canopy%zetar(j,iterplus) ) )         &
                      + psis1( (rough%zruffs(j) - rough%disp(j) )              &
                      / (rough%zref_tq(j)/canopy%zetar(j,iterplus ) ) ) )      &
                      / CVONK
          ELSE!( cable_runtime%esm15_canopy)

              !Ticket #67 - Modify order of operations to avoid potential error
             r_sc(j) = rough%rt0us(j) + rough%rt1usa(j) + rough%rt1usb(j) +     &
                  ( LOG( (zscl(j) - rough%disp(j)) /                       &
                  MAX( rough%zruffs(j)-rough%disp(j),                      &
                  rough%z0soilsn(j) ) ) - psis( (zscl(j)-rough%disp(j))    &
                                !Ticket #67 - change order of operations to avoid /0
                                !        / (rough%zref_tq(j)/canopy%zetar(j,iterplus) ) )        &
                  * canopy%zetar(j,iterplus)/rough%zref_tq(j) )            &
                  + psis( (rough%zruffs(j) - rough%disp(j) )               &
                                !        / (rough%zref_tq(j)/canopy%zetar(j,iterplus ) ) ) )     &
                  * canopy%zetar(j,iterplus)/rough%zref_tq(j) ) )          &
                  / CVONK

          ENDIF!( cable_runtime%esm15_canopy)

Considering the comments, we should probably accept the change in the order of operations.

2 functions psis in trunk and psis1 in ESM. Both are the same, psis1 added to cable_canopy.F90 but psis has been moved to cbl_friction_vel.F90.

penguian commented 1 year ago

@ccc561@nci.org.au commented

        IF ( cable_runtime%esm15_canopy) THEN
          canopy%tscrn(j) = ssnow%tss(j) + (met%tk(j) - ssnow%tss(j)) *          &
                          MIN(1.,r_sc(j) / MAX( 1.,                            &
                          rough%rt0us(j) + rough%rt1usa(j) + rough%rt1usb(j)   &
                          + rt1usc(j))) - Ctfrz 
        ELSE!( cable_runtime%esm15_canopy)

            !extensions for litter and Or evaporation model
            IF (cable_user%litter) THEN
               canopy%tscrn(j) = ssnow%tss(j) + (met%tk(j) - ssnow%tss(j)) *     &
                    MIN(1., ( (r_sc(j)+rhlitt(j)*canopy%us(j))  / MAX( 1.,          &
                    rough%rt0us(j) + rough%rt1usa(j) + rough%rt1usb(j)              &
                    + rt1usc(j) + rhlitt(j)*canopy%us(j) )) ) - Ctfrz
            ELSEIF (cable_user%or_evap .OR. cable_user%gw_model) THEN
               canopy%tscrn(j) = ssnow%tss(j) + (met%tk(j) - ssnow%tss(j)) *     &
                    MIN(1., ( (ssnow%rt_qh_sublayer(j)*canopy%us(j) + r_sc(j) ) /   &
                    MAX( 1., rough%rt0us(j) + rough%rt1usa(j) + rough%rt1usb(j)     &
                    + rt1usc(j) + ssnow%rt_qh_sublayer(j)*canopy%us(j) )) ) - Ctfrz
            ELSE
               canopy%tscrn(j) = ssnow%tss(j) + (met%tk(j) - ssnow%tss(j)) *      &
                    MIN(1., (r_sc(j) / MAX( 1.,                            &
                    rough%rt0us(j) + rough%rt1usa(j) + rough%rt1usb(j)   &
                    + rt1usc(j))) )  - Ctfrz
            ENDIF

          ENDIF!( cable_runtime%esm15_canopy)

Added code for additional options. Default code is the same.

To accept.

penguian commented 1 year ago

@ccc561@nci.org.au commented

        IF ( cable_runtime%esm15_canopy) THEN
          canopy%qscrn(j) = qsurf(j) + (met%qv(j) - qsurf(j)) * MIN( 1.,     &
                                      r_sc(j) / MAX( 1., rough%rt0us(j) +              &
                                      rough%rt1usa(j) + rough%rt1usb(j) + rt1usc(j) ) )
        ELSE!( cable_runtime%esm15_canopy)

          IF (cable_user%litter) THEN

            !extensions for litter and Or model
            canopy%qscrn(j) = qsurf(j) + (met%qv(j) - qsurf(j)) *             &
                     MIN(1., ( ( r_sc(j)+relitt(j)*canopy%us(j) ) / MAX( 1.,         &
                     rough%rt0us(j) + rough%rt1usa(j) + rough%rt1usb(j)            &
                     + rt1usc(j) + relitt(j)*canopy%us(j) )) )

          ELSEIF (cable_user%or_evap .OR. cable_user%gw_model) THEN
            !using alpm1 as a dumy variable
            alpm1(j) = REAL(&
                     ssnow%satfrac(j)/(REAL(ssnow%rtsoil(j),r_2)+&
                     ssnow%rtevap_sat(j)) &
                     + (1.0-ssnow%satfrac(j))/(REAL(ssnow%rtsoil(j),r_2)+ ssnow%rtevap_unsat(j)) &
                     )

            canopy%qscrn(j) = qsurf(j) + (met%qv(j) - qsurf(j)) *             &
                     MIN(1., ( (r_sc(j) + canopy%us(j)/alpm1(j) ) / MAX( 1.,         &
                     rough%rt0us(j) + rough%rt1usa(j) + rough%rt1usb(j)              &
                     + rt1usc(j) + canopy%us(j)/alpm1(j) )) )

          ELSE
            canopy%qscrn(j) = qsurf(j) + (met%qv(j) - qsurf(j)) *             &
                     MIN(1., (r_sc(j) / MAX( 1.,                                     &
                     rough%rt0us(j) + rough%rt1usa(j) + rough%rt1usb(j)              &
                     + rt1usc(j))) )
          ENDIF

        ENDIF!( cable_runtime%esm15_canopy)

Added code for additional options. Default code is the same.

To accept.

penguian commented 1 year ago

@ccc561@nci.org.au commented

   ! this change of units does not affect next timestep as canopy%through is
   ! re-calc in surf_wetness_fact routine
    IF ( cable_runtime%esm15_canopy) THEN
      canopy%through = canopy%through / dels   ! change units for stash output
    ENDIF!( cable_runtime%esm15_canopy)

canopy%through is used in GW and SLI, so the change of units in canopy could be catastrophic for these applications. Should we put the change of unit somewhere at the end of cbm instead?

REAL CHANGE IN MERGED CODE - Need to fix. Jhan to look at.

penguian commented 1 year ago

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

penguian commented 1 year ago

@ccc561@nci.org.au commented

IF ( cable_runtime%esm15_canopy) THEN
   ssnow%dfh_dtg = air%rho*CCAPP/ssnow%rtsoil      
   ssnow%dfe_ddq = ssnow%wetfac*air%rho*air%rlam*ssnow%cls/ssnow%rtsoil  
   ssnow%ddq_dtg = (Crmh2o/Crmair) /met%pmb * CTETENA*CTETENB * CTETENC   &
                   / ( ( CTETENC + ssnow%tss-Ctfrz )**2 )*EXP( CTETENB *       &
                   ( ssnow%tss-Ctfrz ) / ( CTETENC + ssnow%tss-Ctfrz ) )
   canopy%dgdtg = ssnow%dfn_dtg - ssnow%dfh_dtg - ssnow%dfe_ddq *    &
                  ssnow%ddq_dtg

   bal%drybal = REAL(ecy+hcy) - SUM(rad%rniso,2)                               &
                + CCAPP*Crmair*(tlfy-met%tk)*SUM(rad%gradis,2)  ! YP nov2009

bal%wetbal = canopy%fevw + canopy%fhvw - SUM(rad%rniso,2) * canopy%fwet      &
                + CCAPP*Crmair * (tlfy-met%tk) * SUM(rad%gradis,2) *          &
                canopy%fwet  ! YP nov2009

ENDIF!( cable_runtime%esm15_canopy)

IF ( .NOT. cable_runtime%esm15_canopy) THEN
       !INH: REV_CORR revised sensitivity terms working variable
       rttsoil = ssnow%rtsoil
       IF (cable_user%L_REV_CORR) THEN
          WHERE (canopy%vlaiw > CLAI_THRESH)
            !if %vlaiw<=%LAI_THRESH then %rt1 already added to %rtsoil
            rttsoil = rttsoil + rough%rt1
          ENDWHERE
       ENDIF

IF (cable_user%gw_model .or. cable_user%or_evap) THEN

  ssnow%dfh_dtg = air%rho*CCAPP/(ssnow%rtsoil+ real(ssnow%rt_qh_sublayer))

  !! INH simplifying code for legibility
  !ssnow%dfe_ddq = real(ssnow%satfrac)*air%rho*air%rlam*ssnow%cls/ &
  !     (ssnow%rtsoil+ real(ssnow%rtevap_sat))  +
  !     (1.0-real(ssnow%satfrac))*real(ssnow%rh_srf)*&
  !      air%rho*air%rlam*ssnow%cls/ (ssnow%rtsoil+
  !      real(ssnow%rtevap_unsat) )
   ssnow%dfe_ddq = real(ssnow%satfrac)/(ssnow%rtsoil+ real(ssnow%rtevap_sat))  &
              + (1.0-real(ssnow%satfrac))*real(ssnow%rh_srf)                   &
                   / (ssnow%rtsoil+ real(ssnow%rtevap_unsat) )

  !mrd561 fixes.  Do same thing as INH but has been tested.
  IF (cable_user%L_REV_CORR) THEN
    alpm1  = real(ssnow%satfrac/(real(ssnow%rtsoil,r_2)+ ssnow%rtevap_sat) +     &
                  (1.0-ssnow%satfrac) / (real(ssnow%rtsoil,r_2)+ ssnow%rtevap_unsat ) )
    beta2 = real(ssnow%satfrac/(real(ssnow%rtsoil,r_2)+ ssnow%rtevap_sat) +     &
                 (1.0-ssnow%satfrac) * ssnow%rh_srf                  &
                 / (real(ssnow%rtsoil,r_2)+ ssnow%rtevap_unsat ) )
    WHERE (canopy%vlaiw > CLAI_THRESH)
      alpm1 = alpm1 + 1._r_2/real(rough%rt1,r_2)
      beta_div_alpm  = beta2 / alpm1  !might need limit here
      rttsoil = ssnow%rtsoil + rough%rt1
    ELSEWHERE!if there is no canopy then qa should not change
      beta_div_alpm=0.0  !do not divide by aplm1 prevent issues
      rttsoil = ssnow%rtsoil 
    ENDWHERE
    ssnow%dfh_dtg = air%rho*CCAPP/(rttsoil +               & 
                              real(ssnow%rt_qh_sublayer))
    ssnow%dfe_ddq = real(ssnow%satfrac*(1.0-real(beta_div_alpm,r_2)) /        & 
                    (real(ssnow%rtsoil,r_2)+ ssnow%rtevap_sat) +           &
                    (1.0-ssnow%satfrac)* (ssnow%rh_srf - real(beta_div_alpm,r_2)) /    &
                    (real(ssnow%rtsoil,r_2)+ ssnow%rtevap_unsat ) )

  ELSE ! IF (cable_user%L_REV_CORR) THEN
    ssnow%dfh_dtg = air%rho*CCAPP/(ssnow%rtsoil+ real(ssnow%rt_qh_sublayer))
    ssnow%dfe_ddq = real(ssnow%satfrac)/(ssnow%rtsoil+ real(ssnow%rtevap_sat))  &
                         + (1.0-real(ssnow%satfrac))*real(ssnow%rh_srf)                   &
                              / (ssnow%rtsoil+ real(ssnow%rtevap_unsat) )
  ENDIF ! IF (cable_user%L_REV_CORR) THEN

  !cls applies for both REV_CORR false and true          
  ssnow%dfe_ddq = ssnow%dfe_ddq*air%rho*air%rlam*ssnow%cls

  !REV_CORR: factor %wetfac needed for potev>0. and gw_model &/or snow cover
  !NB %wetfac=1. if or_evap
  IF (cable_user%L_REV_CORR) THEN
    WHERE (ssnow%potev >= 0.)
      ssnow%dfe_ddq = ssnow%dfe_ddq*ssnow%wetfac
    ENDWHERE       
  ENDIF

ELSEIF (cable_user%litter) THEN ! IF (cable_user%gw_model .or. cable_user%or_evap) THEN
  !!vh_js!! INH simplifying code for legibility and REV_CORR
  !ssnow%dfh_dtg = air%rho*CCAPP/(ssnow%rtsoil+ &
  !     real((1-ssnow%isflag))*veg%clitt*0.003/canopy%kthLitt/(air%rho*CCAPP))
  !ssnow%dfe_ddq = ssnow%wetfac*air%rho*air%rlam*ssnow%cls/ &
  !     (ssnow%rtsoil+ real((1-ssnow%isflag))*veg%clitt*0.003/canopy%DvLitt)

  !recalculated - probably not needed 
  rhlitt = real((1-ssnow%isflag))*veg%clitt*0.003/canopy%kthLitt/(air%rho*CCAPP)
  relitt = real((1-ssnow%isflag))*veg%clitt*0.003/canopy%DvLitt

  !incorporates REV_CORR changes
  ssnow%dfh_dtg = air%rho*CCAPP/(rttsoil+rhlitt)
  ssnow%dfe_ddq = ssnow%wetfac*air%rho*air%rlam*ssnow%cls/(rttsoil+relitt)

  !REV_CORR: factor ssnow%wetfac is not applied if dew/frost i.e. potev<0
  IF (cable_user%L_REV_CORR) THEN
     WHERE (ssnow%potev < 0.)
         ssnow%dfe_ddq = air%rho*air%rlam*ssnow%cls/(rttsoil+relitt)
     ENDWHERE       
  ENDIF

ELSE ! i.e. NOT (%gw_model .or. %or_evap or SLI)
  !ssnow%dfh_dtg = air%rho*CCAPP/ssnow%rtsoil
  !ssnow%dfe_ddq = ssnow%wetfac*air%rho*air%rlam*ssnow%cls/ssnow%rtsoil

  !incorporates REV_CORR changes
  ssnow%dfh_dtg = air%rho*CCAPP/rttsoil
  ssnow%dfe_ddq = ssnow%wetfac*air%rho*air%rlam*ssnow%cls/rttsoil

  !REV_CORR: factor ssnow%wetfac is not applied if dew/frost i.e. potev<0
  IF (cable_user%L_REV_CORR) THEN
     WHERE (ssnow%potev < 0.)
         ssnow%dfe_ddq = air%rho*air%rlam*ssnow%cls/(rttsoil+relitt)
     ENDWHERE       
  ENDIF      

ENDIF ! IF (cable_user%gw_model .or. cable_user%or_evap) THEN

ssnow%ddq_dtg = (Crmh2o/Crmair) /met%pmb * CTETENA*CTETENB * CTETENC   &
     / ( ( CTETENC + ssnow%tss-Ctfrz )**2 )*EXP( CTETENB *       &
     ( ssnow%tss-Ctfrz ) / ( CTETENC + ssnow%tss-Ctfrz ) )

!canopy%dgdtg = ssnow%dfn_dtg - ssnow%dfh_dtg - ssnow%dfe_ddq *    &
!     ssnow%ddq_dtg

!INH: REV_CORR Rewritten for flexibility
ssnow%dfe_dtg = ssnow%dfe_ddq * ssnow%ddq_dtg
canopy%dgdtg = ssnow%dfn_dtg - ssnow%dfh_dtg - ssnow%dfe_dtg

bal%drybal = REAL(ecy+hcy) - SUM(rad%rniso,2)                               &
     + CCAPP*Crmair*(tlfy-met%tk)*SUM(rad%gradis,2)  ! YP nov2009

bal%wetbal = canopy%fevw + canopy%fhvw - SUM(rad%rniso,2) * canopy%fwet      &
     + CCAPP*Crmair * (tlfy-met%tk) * SUM(rad%gradis,2) *          &
     canopy%fwet  ! YP nov2009

ENDIF!( cable_runtime%esm15_canopy)

Calculations for bal% variables haven't changed, could these be taken out of the flag?

ssnow%dfh_dtg and ssnow%dfe_ddq: change from using ssnow%rtsoil to rttsoil. rttsoil and ssnow%rtsoil are the same if not using L_REV_CORR.

Additional condition on potev is only triggerd with L_REV_CORR. Safe to add.

canopy%dgdtg change from 1 liner to 2 lines.

penguian commented 1 year ago

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

penguian commented 1 year ago

@jxs599@nci.org.au changed status from new to closed

penguian commented 1 year ago

@jxs599@nci.org.au set resolution to fixed

penguian commented 1 year ago

@jxs599@nci.org.au set milestone to 1. Closed

penguian commented 1 year ago

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