module catch_diagn_routines ! subroutines for diagnostics of the Catchment model ! ! reichle, Feb 5, 2004 ! ! revised subroutine calc_soil_moist - koster+reichle, Mar 19, 2004 ! ! added calc_tsurf(), reichle, 31 Aug 2004 ! ! added calc_arX, reichle, 21 Jun 2005 ! ! reichle, 12 Sep 2007: ! deleted calc_soil_moist private because it is now provided with catchment.F90 ! added "use catchment_model" ! ! ------------------------------------------------------------------ use catch_constants, ONLY: & TF => MAPL_TICE, & N_snow => CATCH_N_SNOW use StieglitzSnow, ONLY: wemin use catch_types use catchment_model implicit none private public :: calc_tsurf, calc_arX, calc_asnow contains !!$#if 0 !!$ subroutine calc_soil_moist( & !!$ ncat,vegcls,vgwmax,cdcr1,cdcr2,wpwet,poros, & !!$ psis,bee,ars1,ars2,ars3,ara1,ara2, & !!$ ara3,ara4,arw1,arw2,arw3,arw4, & !!$ srfexc,rzexc,catdef, & !!$ sfmc, rzmc, prmc, & !!$ werror, sfmcun, rzmcun, prmcun ) !!$ !!$ ! Calculate diagnostic soil moisture content from prognostic !!$ ! excess/deficit variables. !!$ ! !!$ ! On input, also check validity of prognostic excess/deficit variables !!$ ! and modify if necessary. Perturbed or updated excess/deficit variables !!$ ! in data assimilation integrations may be unphysical. !!$ ! Optional output "werror" contains excess or missing water related !!$ ! to inconsistency. !!$ ! !!$ ! Optional outputs "smfcun", "rzmcun", "prmcun" are surface, !!$ ! root zone, and profile moisture content for unsaturated areas only, !!$ ! ie. excluding the saturated area of the catchment. !!$ ! !!$ ! NOTE: When calling with optional output arguments, use keywords !!$ ! unless arguments are in proper order! !!$ ! !!$ ! Example: !!$ ! (do not want "werror" as output, but want "*mcun" output) !!$ ! !!$ ! call calc_soil_moist( & !!$ ! ncat, ... & !!$ ! sfmc, rzmc, prmc, & !!$ ! sfmcun=sfmc_unsat, & !!$ ! rzmcun=rzmc_unsat, & !!$ ! prmcun=prmc_unsat ) !!$ ! !!$ ! replaces moisture_sep_22_2003.f (and older moisture.f) !!$ ! !!$ ! koster+reichle, Feb 5, 2004 !!$ ! !!$ ! revised - koster+reichle, Mar 19, 2004 !!$ ! !!$ ! added optional *un output - koster+reichle, Apr 6, 2004 !!$ ! !!$ ! ---------------------------------------------------------------- !!$ !!$ implicit none !!$ !!$ integer, intent(in) :: ncat !!$ integer, dimension(ncat), intent(in) :: vegcls !!$ !!$ real, dimension(ncat), intent(in) :: vgwmax,cdcr1,cdcr2 !!$ real, dimension(ncat), intent(in) :: wpwet,poros,psis !!$ real, dimension(ncat), intent(in) :: bee,ars1 !!$ real, dimension(ncat), intent(in) :: ars2,ars3,ara1,ara2,ara3 !!$ real, dimension(ncat), intent(in) :: ara4,arw1,arw2,arw3,arw4 !!$ !!$ real, dimension(ncat), intent(inout) :: srfexc, rzexc, catdef !!$ !!$ real, dimension(ncat), intent(out) :: sfmc, rzmc, prmc !!$ !!$ real, dimension(ncat), intent(out), optional :: werror !!$ !!$ real, dimension(ncat), intent(out), optional :: sfmcun !!$ real, dimension(ncat), intent(out), optional :: rzmcun !!$ real, dimension(ncat), intent(out), optional :: prmcun !!$ !!$ ! ---------------------------- !!$ ! !!$ ! local variables !!$ !!$ integer :: n !!$ !!$ real, parameter :: dtstep_dummy = -9999. !!$ !!$ real, dimension(ncat) :: rzeq, runsrf_dummy, catdef_dummy !!$ real, dimension(ncat) :: ar1, ar2, ar4, prmc_orig !!$ real, dimension(ncat) :: srfmn, srfmx, swsrf1, swsrf2, swsrf4 !!$ !!$ !!$ ! -------------------------------------------------------------------- !!$ ! !!$ ! compute soil water storage upon input [mm] !!$ !!$ do n=1,ncat !!$ prmc_orig(n) = & !!$ (cdcr2(n)/(1.-wpwet(n))-catdef(n)+rzexc(n)+srfexc(n)) !!$ enddo !!$ !!$ ! ----------------------------------- !!$ ! !!$ ! check limits of catchment deficit !!$ ! !!$ ! increased minimum catchment deficit from 0.01 to 1. to make the !!$ ! check work with perturbed parameters and initial condition !!$ ! reichle, 16 May 01 !!$ ! !!$ ! IT REALLY SHOULD WORK WITH catdef > 0 (rather than >1.) ???? !!$ ! reichle, 5 Feb 2004 !!$ ! !!$ ! back to catdef>=1. because of flaws in partition !!$ ! reichle, 26 March 2007 !!$ !!$ do n=1,ncat !!$ catdef(n)=max(1.,min(cdcr2(n),catdef(n))) !!$ end do !!$ !!$ ! ------------------------------------------------------------------ !!$ ! !!$ ! check limits of root zone excess !!$ ! !!$ ! calculate root zone equilibrium moisture for given catchment deficit !!$ !!$ call rzequil( & !!$ ncat, vegcls, catdef, vgwmax, & !!$ cdcr1, cdcr2, wpwet, & !!$ ars1, ars2, ars3, ara1, ara2, ara3, ara4, & !!$ arw1, arw2, arw3, arw4, & !!$ rzeq) !!$ !!$ ! assume srfexc=0 and constrain rzexc appropriately !!$ ! (iteration would be needed to contrain srfexc and rzexc simultaneously) !!$ !!$ do n=1,ncat !!$ rzexc(n)=max(wpwet(n)*vgwmax(n)-rzeq(n),min(vgwmax(n)-rzeq(n),rzexc(n))) !!$ end do !!$ !!$ ! this translates into: !!$ ! !!$ ! wilting level < rzmc < porosity !!$ ! !!$ ! or more precisely: wpwet*vgwmax < rzeq+rzexc < vgwmax !!$ ! !!$ ! NOTE: root zone moisture is not allowed to drop below wilting level !!$ !!$ ! ----------------------------------------------------------------- !!$ ! !!$ ! Call partition() for computation of surface moisture content. !!$ ! !!$ ! Call to partition() also checks limits of surface excess. !!$ ! !!$ ! Call partition with dtstep_dummy: !!$ ! In partition, dtstep is only used for a correction that !!$ ! puts water into runsrf (for which runsrf_dummy is used here). !!$ ! Also use catdef_dummy because partition() updates catdef !!$ ! whenever srfexc exceeds physical bounds, but this is not desired here. !!$ !!$ runsrf_dummy = 0. !!$ catdef_dummy = catdef !!$ !!$ call partition( & !!$ ncat,dtstep_dummy,vegcls,rzexc, & !!$ rzeq,vgwmax,cdcr1,cdcr2, & !!$ psis,bee,poros,wpwet, & !!$ ars1,ars2,ars3, & !!$ ara1,ara2,ara3,ara4, & !!$ arw1,arw2,arw3,arw4,.false., & !!$ srfexc,catdef_dummy,runsrf_dummy, & !!$ ar1, ar2, ar4,srfmx,srfmn, & !!$ swsrf1,swsrf2,swsrf4 & !!$ ) !!$ !!$ ! compute surface, root zone, and profile soil moisture !!$ !!$ do n=1,ncat !!$ !!$ sfmc(n) = poros(n) * & !!$ (swsrf1(n)*ar1(n) + swsrf2(n)*ar2(n) + swsrf4(n)*ar4(n)) !!$ !!$ rzmc(n) = (rzeq(n)+rzexc(n)+srfexc(n))*poros(n)/vgwmax(n) !!$ !!$ ! compute revised soil water storage [mm] !!$ !!$ prmc(n) = & !!$ (cdcr2(n)/(1.-wpwet(n))-catdef(n)+rzexc(n)+srfexc(n)) !!$ !!$ ! compute error in soil water storage [mm] (if argument is present) !!$ !!$ if (present(werror)) werror(n)=(prmc(n)-prmc_orig(n)) !!$ !!$ ! convert to volumetric soil moisture !!$ ! note: dzpr = (cdcr2/(1-wpwet)) / poros !!$ !!$ prmc(n) = prmc(n)*poros(n) / (cdcr2(n)/(1.-wpwet(n))) !!$ !!$ !!$ ! check for negative soil moisture !!$ !!$ if ( (sfmc(n)<.0) .or. (rzmc(n)<.0) .or. (prmc(n)<.0) ) then !!$ !!$ write (*,*) 'FOUND NEGATIVE SOIL MOISTURE CONTENT.... stopping' !!$ write (*,*) n, sfmc(n), rzmc(n), prmc(n) !!$ stop !!$ end if !!$ !!$ ! compute moisture content in unsaturated areas [m3/m3] (if arg present) !!$ !!$ if (ar1(n)<1.) then !!$ !!$ if (present(prmcun)) prmcun(n)=(prmc(n)-poros(n)*ar1(n))/(1.-ar1(n)) !!$ if (present(rzmcun)) rzmcun(n)=(rzmc(n)-poros(n)*ar1(n))/(1.-ar1(n)) !!$ if (present(sfmcun)) sfmcun(n)=(sfmc(n)-poros(n)*ar1(n))/(1.-ar1(n)) !!$ !!$ else !!$ !!$ if (present(prmcun)) prmcun(n)=poros(n) !!$ if (present(rzmcun)) rzmcun(n)=poros(n) !!$ if (present(sfmcun)) sfmcun(n)=poros(n) !!$ !!$ end if !!$ !!$ enddo !!$ !!$ return !!$ !!$ end subroutine calc_soil_moist !!$#endif ! ******************************************************************** subroutine calc_tsurf( N_cat, cat_param, cat_progn, tsurf ) ! Calculate diagnostic surface temperature "tsurf" from prognostics ! ! reichle, Aug 31, 2004 ! ! ---------------------------------------------------------------- implicit none integer, intent(in) :: N_cat type(cat_param_type), dimension(N_cat), intent(in) :: cat_param type(cat_progn_type), dimension(N_cat), intent(in) :: cat_progn real, dimension(N_cat), intent(out) :: tsurf ! ---------------------------- ! ! local variables integer :: n real, dimension(N_cat) :: ar1, ar2, ar4, prmc_orig real, dimension(N_cat) :: asnow real, dimension(1) :: tpsn1, real_dummy !!real, parameter :: tf = 273.16 ! ------------------------------------------------------------------ call calc_arX( N_cat, cat_param, cat_progn, ar1, ar2, ar4 ) ! Compute snow covered area call calc_asnow( N_cat, cat_progn, asnow ) ! Compute tsurf do n=1,N_cat tsurf(n) = ar1(n)*cat_progn(n)%tc1 + ar2(n)*cat_progn(n)%tc2 & + ar4(n)*cat_progn(n)%tc4 if (asnow(n)>0.) then ! get_tf_nd returns snow temperature in deg Celsius call get_tf_nd( 1, cat_progn(n)%htsn(1), cat_progn(n)%wesn(1), & tpsn1, real_dummy ) tsurf(n) = (1. - asnow(n))*tsurf(n) + asnow(n)*(tpsn1(1) + tf) end if enddo return end subroutine calc_tsurf ! ******************************************************************* subroutine calc_arX( N_cat, cat_param, cat_progn, ar1, ar2, ar4 ) ! reichle, 21 Jun 2005 ! reichle, 12 Sep 2007 - new interface for partition implicit none integer, intent(in) :: N_cat type(cat_param_type), dimension(N_cat), intent(in) :: cat_param type(cat_progn_type), dimension(N_cat), intent(in) :: cat_progn real, dimension(N_cat), intent(out) :: ar1, ar2, ar4 ! locals real, parameter :: dtstep_dummy = -9999. real, dimension(N_cat) :: rzeq, runsrf_dummy real, dimension(N_cat) :: catdef_dummy, rzexc_dummy, srfexc_dummy real, dimension(N_cat) :: srfmx, srfmn, swsrf1, swsrf2, swsrf4, rzi ! ------------------------------------------------------------------ ! ! Call partition to get saturated/unsaturated/wilting ! areas ar1, ar2, and ar4: ! ! Need to calculate root zone equilibrium moisture for given ! catchment deficit (needed for call to partition). call rzequil( & N_cat, cat_param%vegcls, cat_progn%catdef, cat_param%vgwmax, & cat_param%cdcr1, cat_param%cdcr2, cat_param%wpwet, & cat_param%ars1, cat_param%ars2, cat_param%ars3, & cat_param%ara1, cat_param%ara2, cat_param%ara3, cat_param%ara4, & cat_param%arw1, cat_param%arw2, cat_param%arw3, cat_param%arw4, & rzeq) ! Call partition with dtstep_dummy: ! In partition, dtstep is only used for a correction that ! puts water into runsrf (for which runsrf_dummy is used here). ! Also use catdef_dummy etc because partition() updates catdef ! whenever srfexc exceeds physical bounds, but this is not desired here. runsrf_dummy = 0. catdef_dummy = cat_progn%catdef rzexc_dummy = cat_progn%rzexc srfexc_dummy = cat_progn%srfexc call partition( & N_cat, dtstep_dummy, cat_param%vegcls, rzexc_dummy, & rzeq, cat_param%vgwmax, cat_param%cdcr1, cat_param%cdcr2, & cat_param%psis, cat_param%bee, cat_param%poros, cat_param%wpwet, & cat_param%ars1, cat_param%ars2, cat_param%ars3, & cat_param%ara1, cat_param%ara2, cat_param%ara3, cat_param%ara4, & cat_param%arw1, cat_param%arw2, cat_param%arw3, cat_param%arw4, & .false., & srfexc_dummy, catdef_dummy, runsrf_dummy, & ar1, ar2, ar4, srfmx, srfmn, swsrf1, swsrf2, swsrf4, rzi ) end subroutine calc_arX ! ******************************************************************* subroutine calc_asnow( N_cat, cat_progn, asnow ) ! Calculate diagnostic snow area from prognostic SWE ! ! reichle, Nov 3, 2004 ! ! ---------------------------------------------------------------- implicit none integer, intent(in) :: N_cat type(cat_progn_type), dimension(N_cat), intent(in) :: cat_progn real, dimension(N_cat), intent(out) :: asnow ! local variables integer :: n ! "wemin" MUST BE CONSISTENT WITH "wemin" IN SUBROUTINE SNOWRT()! !!real, parameter :: wemin = 13. ! [kg/m^2] ! ----------------------------------------------------------- do n=1,N_cat asnow(n) = min( sum(cat_progn(n)%wesn(1:N_snow))/wemin, 1.) end do end subroutine calc_asnow end module catch_diagn_routines ! *************************** EOF ****************************************