#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_BASEFLOW------------------------------------------ ! Description: ! Calculates methane emissions from wetland area. ! Subroutine Interface: SUBROUTINE calc_baseflow( & soil_pts,soil_index,npnts,nshyd & ,zdepth,ksz & ,b,fexp,ti_mean,zw,sthf,sthu & ,wutot,top_crit,qbase,qbase_l & ) USE c_topog USE soil_param, ONLY : dzsoil USE parkind1, ONLY: jprb, jpim USE yomhook, ONLY: lhook, dr_hook IMPLICIT NONE INTEGER & npnts & ! IN Number of gridpoints. ,nshyd & ! IN Number of soil moisture levels. ,soil_pts ! IN Number of soil points. ! Array arguments with intent(IN) : INTEGER & soil_index(npnts) ! IN Array of soil points. REAL & b(npnts) & ! IN Clapp-Hornberger exponent. ,fexp(npnts) & ! IN Decay factor in Sat. Conductivity ! ! in water table layer. ,ti_mean(npnts) & ! IN Mean topographic index. ,zw(npnts) & ! IN Water table depth (m). ,ksz(npnts,0:nshyd) & ! IN Saturated hydraulic condy ! ! for each layer (kg/m2/s). ,ksfz(npnts,nshyd+1) & ! WORK Function of sat. hydraulic ! ! conductivity, frozen soil ! ! and mean topographic index. ,sthf(npnts,nshyd) & ! INOUT Frozen soil moisture content of ! ! each layer as a fraction of ! ! saturation. ,sthu(npnts,nshyd) ! INOUT Unfrozen soil moisture content of ! ! each layer as a fraction of ! ! saturation. REAL & qbase(npnts) & ! OUT Total base flow (kg/m2/s). ,qbase_l(npnts,nshyd+1) & ! ! OUT Base flow from each layer (kg/m2/s). ,top_crit(npnts) & ! OUT Critical topographic index required ! ! to calc surface saturation frac. ,wutot(npnts) ! OUT UNFROZEN to TOTAL fraction at ZW. REAL & qbase_max(npnts) & ! WORK Max possible base flow (kg/m2/s). ,qbase_max_l(npnts,nshyd+1) & ! ! WORK Max possible base flow ! ! from each level (kg/m2/s). ,zdepth(0:nshyd) & ! WORK Lower soil layer boundary depth (m). ,qbase_min(npnts) ! WORK Residual base flow at zw_max ! Local scalars: INTEGER & i,j & ! WORK Loop counters. ,n & ! WORK Tile loop counter. ,errorstatus ! WORK Error status for ereport. INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 REAL(KIND=jprb) :: zhook_handle !----------------------------------------------------------------------- ! Initialise TOP_CRIT to maximum value. ! Initialise baseflow components to zero. !----------------------------------------------------------------------- IF (lhook) CALL dr_hook('CALC_BASEFLOW',zhook_in,zhook_handle) DO j=1,soil_pts i=soil_index(j) top_crit(i)=ti_max qbase_max(i)=0.0 DO n=1,nshyd qbase_max_l(i,n)=0.0 END DO END DO DO j=1,soil_pts i=soil_index(j) !----------------------------------------------------------------------- ! Calculate layer dependent variable with is dependent on ! effective saturated conductivity: !----------------------------------------------------------------------- DO n=1,nshyd IF(sthf(i,n) < 1.0)THEN ksfz(i,n)=0.5*(ksz(i,n-1)+ksz(i,n)) & *(1.0-sthf(i,n))**(2.0*b(i)+3.0)*EXP(-ti_mean(i)) ELSE ksfz(i,n)=0.0 END IF END DO IF(sthf(i,nshyd) < 1.0)THEN ksfz(i,nshyd+1)=ksz(i,nshyd)/fexp(i) & *(1.0-sthf(i,nshyd))**(2.0*b(i)+3.0) & *EXP(-ti_mean(i)) ELSE ksfz(i,nshyd+1)=0.0 END IF IF(EXP(-fexp(i)*(zw_max-zdepth(nshyd))) > 0.05)THEN WRITE(6,*)'CB_L1A: maximum water table depth is too low!!!' WRITE(6,*)'at soil point',i,fexp(i) WRITE(6,*)EXP(-fexp(i)*(zw_max-zdepth(nshyd))) WRITE(6,*)'try zw_max>',-LOG(0.05)/fexp(i)+zdepth(nshyd) errorstatus=1000 ! CALL EREPORT('CALC_BASEFLOW', ERRORSTATUS, & ! & 'ERROR ZW_MAX is TOO SMALL') END IF !----------------------------------------------------------------------- ! Calculate base flow between maximum allowed water table depth and ! "infinity": !----------------------------------------------------------------------- qbase_min(i)=ksfz(i,nshyd+1) * EXP(-fexp(i)*(zw_max-zdepth(nshyd))) !----------------------------------------------------------------------- ! Calculate maximum possible and actual base flow for each layer: !----------------------------------------------------------------------- DO n=1,nshyd qbase_max_l(i,n)=ksfz(i,n)*(zdepth(n)-zdepth(n-1)) IF(zw(i) <= zdepth(n-1)) qbase_l(i,n)=qbase_max_l(i,n) IF(zw(i) < zdepth(n).AND.zw(i) > zdepth(n-1))THEN qbase_l(i,n)=ksfz(i,n)*(zdepth(n)-zw(i)) IF(sthu(i,n)+sthf(i,n) > 0.0) & wutot(i)=sthu(i,n)/(sthu(i,n)+sthf(i,n)) END IF IF(n == 1.AND.zw(i) < zdepth(n))THEN qbase_l(i,n)=ksfz(i,n)*(zdepth(n)-zw(i)) IF(sthu(i,n)+sthf(i,n) > 0.0) & wutot(i)=sthu(i,n)/(sthu(i,n)+sthf(i,n)) END IF END DO qbase_max_l(i,nshyd+1)=ksfz(i,nshyd+1) - qbase_min(i) IF(zw(i) <= zdepth(nshyd))THEN qbase_l(i,nshyd+1)=qbase_max_l(i,nshyd+1) ELSE qbase_l(i,nshyd+1)=ksfz(i,nshyd+1) & *EXP(-fexp(i)*(zw(i)-zdepth(nshyd))) & -qbase_min(i) IF(sthu(i,nshyd)+sthf(i,nshyd) > 0.0) & wutot(i)=sthu(i,nshyd)/(sthu(i,nshyd)+sthf(i,nshyd)) END IF !----------------------------------------------------------------------- ! Calculate total possible and actual base flow: !----------------------------------------------------------------------- DO n=1,nshyd+1 qbase_l(i,n)=MAX(0.0,qbase_l(i,n)) qbase(i)=qbase(i)+qbase_l(i,n) qbase_max(i)=qbase_max(i)+qbase_max_l(i,n) END DO !----------------------------------------------------------------------- ! Calculate critical topographic index. !----------------------------------------------------------------------- IF(qbase(i) > qbase_max(i))qbase(i)=qbase_max(i) !Check that QBASE_MAX(I)/QBASE(I) will not underflow. IF(qbase_max(i) > EPSILON(qbase_max(i)) .AND. & qbase(i) > qbase_max(i)*(EPSILON(qbase(i)))) & top_crit(i) = LOG(qbase_max(i)/qbase(i)) END DO IF (lhook) CALL dr_hook('CALC_BASEFLOW',zhook_out,zhook_handle) RETURN END SUBROUTINE calc_baseflow #endif