#if defined(L08_1A) && !defined(RECON) ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* ! SUBROUTINE CALC_FSAT---------------------------------------------- ! Description: ! Calculates the surface saturation fraction. ! Subroutine Interface: SUBROUTINE calc_fsat(l_gamtot,soil_pts,soil_index, & npnts,ti_mean,ti_sig,wutot,top_crit,gamtot,fsat,fwetl) USE c_topog USE parkind1, ONLY: jprb, jpim USE yomhook, ONLY: lhook, dr_hook IMPLICIT NONE INTEGER & npnts & !IN No. of land points. ,soil_pts & ,soil_index(soil_pts) REAL & ti_mean(npnts) & !IN Gridbox mean topographic index. ,ti_sig(npnts) & !IN Standard deviation in topographic index. ,gamtot(npnts) & !IN Integrated complete Gamma function. ! ! Unless L_GAMTOT=TRUE, then INOUT ,top_crit(npnts) & !IN Critical topographic index required ! ! to calculate the surface saturation fraction. ,wutot(npnts) !IN UNFROZEN to TOTAL fraction at ZW. LOGICAL & l_gamtot !IN TRUE if need to calculate GAMTOT REAL & fsat(npnts) & !INOUT Surface saturation fraction. ,fwetl(npnts) !INOUT Wetland fraction. REAL & ti_max_use(npnts) & ! !WORK TI_MAX dependent on ti_sig ,dti_use(npnts) & !WORK increment which dependent on ti_sig ,alf & !WORK Parameter in incomplete Gamma function. ,alf_ksat & !WORK Parameter in incomplete Gamma function ! ! including horizontal ksat variability. ,cum & !WORK Integrated incomplete Gamma fn. ,ti_sc & !WORK Incremented Topographic index. ! ! Scaled by ti_mean ,dti_sc & !WORK Scale increment in TI ,ticr_sc & !WORK Critical topographic index. (Any value ! ! above this gives surface saturation). ! ! Scaled by TI_MEAN ,ticr_sc_w & !WORK As above, but for wetland fraction calc. ,fbox_s & !WORK Fraction of box for integration. ,fbox_w !WORK Fraction of box for integration. INTEGER & i,j & ! Loop counter for points. ,nti & ! Loop counter for topographic index. ,mti & ! Max loop count for topographic index. ,mti_w ! Max loop count for topographic index. INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 REAL(KIND=jprb) :: zhook_handle !---------------------------------------------------------------------- ! Set up maximum topographic index to integrate to and increment ! factor. These are topography dependent to maximise accuracy and ! minimise runtime: !---------------------------------------------------------------------- IF (lhook) CALL dr_hook('CALC_FSAT',zhook_in,zhook_handle) DO j=1,soil_pts i=soil_index(j) ti_max_use(i)=ti_mean(i)+5.*ti_sig(i) dti_use(i)=dti/(ti_sig(i)*ti_sig(i)) END DO !---------------------------------------------------------------------- ! Calculate the total integral under the Gamma function. Carried ! out via the reconfiguration: !---------------------------------------------------------------------- IF(l_gamtot)THEN DO j=1,soil_pts i=soil_index(j) mti=NINT(ti_max_use(i)/dti_use(i)) dti_sc=dti_use(i)/ti_mean(i) alf=(ti_mean(i)/ti_sig(i))**2 alf_ksat=1./(1./alf+(sigma_logk/ti_mean(i))**2) DO nti=1,mti ti_sc=(nti-0.5)*dti_sc gamtot(i)=gamtot(i) & +ti_sc**(alf_ksat-1.0)*EXP(-alf_ksat*ti_sc) END DO END DO ELSE !---------------------------------------------------------------------- ! Calculate the integral under the incomplete Gamma function for the ! saturation surface fraction: !---------------------------------------------------------------------- DO j=1,soil_pts i=soil_index(j) IF(top_crit(i) < ti_max_use(i))THEN ticr_sc=top_crit(i)/ti_mean(i)+1.0 ticr_sc_w=ticr_sc+ti_wetl/ti_mean(i) IF(ticr_sc*ti_mean(i) < ti_max_use(i))THEN alf=(ti_mean(i)/ti_sig(i))**2 alf_ksat=1./(1./alf+(sigma_logk/ti_mean(i))**2) dti_sc=dti_use(i)/ti_mean(i) mti=INT(ticr_sc/dti_sc) cum=0.0 DO nti=1,mti ti_sc=(nti-0.5)*dti_sc cum=cum+ti_sc**(alf_ksat-1.0) & *EXP(-alf_ksat*ti_sc) END DO ! Include the fractional increment: fbox_s=(ticr_sc-mti*dti_sc)/dti_sc ti_sc=(mti+0.5*fbox_s)*dti_sc cum=cum+ti_sc**(alf_ksat-1.0) & *EXP(-alf_ksat*ti_sc)*fbox_s fsat(i)=1.0-cum/gamtot(i) !---------------------------------------------------------------------- ! Calculate the integral under the incomplete Gamma function for the ! wetland fraction: !---------------------------------------------------------------------- IF(ticr_sc_w*ti_mean(i) < ti_max_use(i))THEN mti_w=INT(ticr_sc_w/dti_sc) ! Include the fractional increment: IF(mti_w == mti)THEN fbox_w=(ticr_sc_w-mti*dti_sc)/dti_sc ELSE fbox_w=1.0 END IF ti_sc=(mti+0.5*fbox_w)*dti_sc cum=cum+ti_sc**(alf_ksat-1.0) & *EXP(-alf_ksat*ti_sc)*(fbox_w-fbox_s) DO nti=mti+2,mti_w ti_sc=(nti-0.5)*dti_sc cum=cum+ti_sc**(alf_ksat-1.0) & *EXP(-alf_ksat*ti_sc) END DO ! Include the fractional increment: fbox_w=(ticr_sc_w-mti_w*dti_sc)/dti_sc ti_sc=(mti_w+0.5*fbox_w)*dti_sc cum=cum+ti_sc**(alf_ksat-1.0) & *EXP(-alf_ksat*ti_sc)*fbox_w fwetl(i)=-1.0+cum/gamtot(i)+fsat(i) END IF IF(fwetl(i) < 0.0)fwetl(i)=0.0 IF(fsat(i) < 0.0)fsat(i)=0.0 IF(fwetl(i) > fsat(i))fwetl(i)=fsat(i) !---------------------------------------------------------------------- ! Assume that in partially frozen water, flow is inhibited, therefore ! the critical flow for wetland area is no longer valid: !---------------------------------------------------------------------- fwetl(i)=fwetl(i)*wutot(i)+fsat(i)*(1.0-wutot(i)) END IF END IF END DO END IF IF (lhook) CALL dr_hook('CALC_FSAT',zhook_out,zhook_handle) RETURN END SUBROUTINE calc_fsat !----------------------------------------------------------------------- #endif