From 15a6f113185fc27e965402cabbb599a87be14305 Mon Sep 17 00:00:00 2001 From: crazyzlj Date: Tue, 1 Aug 2023 19:23:50 +0800 Subject: [PATCH 1/2] Version 689 Summary: - (Arun) Tropical growth code added including: mgtop 18 added (mgt_sched.f) 18 == monsoon planting for tropical growth 19 added " 19 == reset phenology during monsoon season for tropical perennial growth - (Arun) Auto irrigation problem resolved main.f - revision and date change Tropical growth: allocate_parms.f modparm.f readfcst.f readsub.f readwgn.f sched_mgt.f subbasin.f auto irrigatioin: autoirr.f --- VERSIONS | 2 +- src/allocate_parms.f | 7 ++++ src/autoirr.f | 9 +---- src/header.f | 2 +- src/main.f | 2 +- src/modparm.f | 19 +++++++++- src/readfcst.f | 2 +- src/readsub.f | 14 +++++++ src/readwgn.f | 19 +++++++++- src/rteinit.f | 4 +- src/sched_mgt.f | 30 ++++++++++++++- src/stdaa.f | 3 +- src/subbasin.f | 88 ++++++++++++++++++++++++++++++++++++++++---- src/writed.f | 1 - 14 files changed, 176 insertions(+), 26 deletions(-) diff --git a/VERSIONS b/VERSIONS index 9f75c79..8222130 100644 --- a/VERSIONS +++ b/VERSIONS @@ -1,3 +1,3 @@ VERSION_MAJOR 2012 -VERSION_MINOR 688 +VERSION_MINOR 689 VERSION_PATCH \ No newline at end of file diff --git a/src/allocate_parms.f b/src/allocate_parms.f index ce4c56a..276db01 100644 --- a/src/allocate_parms.f +++ b/src/allocate_parms.f @@ -60,6 +60,13 @@ subroutine allocate_parms mstdo = 113 motot = 1200 !!! 100 year limit +!! ppet arrays for tropical growth + allocate (ppet(mhru)) + do ihru = 1, mhru + allocate (ppet(ihru)%precip(ppet(ihru)%ndays)) + allocate (ppet(ihru)%pet(ppet(ihru)%ndays)) + end do + !! Srini 11_1_22 allocate (tmp_win1(mhru)) allocate (tmp_win2(mhru)) diff --git a/src/autoirr.f b/src/autoirr.f index 880e870..26c05d7 100644 --- a/src/autoirr.f +++ b/src/autoirr.f @@ -83,11 +83,6 @@ subroutine autoirr j = 0 j = ihru -!!!! Srin's irrigation source by each application changes - irrsc(j) = irr_sca(j) - irrno(j) = irr_noa(j) -!!!! Srin's irrigation source by each application changes - if ((wstrs_id(j) == 3) .or. & (wstrs_id(j) == 1 .and. strsw(j) < auto_wstr(j)) .or. & (wstrs_id(j)==2.and.sol_sumfc(j)-sol_sw(j)>auto_wstr(j))) then @@ -106,7 +101,7 @@ subroutine autoirr select case (irrsc(j)) case (3) !! shallow aquifer source do k = 1, nhru - if (hru_sub(k) == irrno(j)) then + if (k == irrno(j)) then cnv = 0. cnv = hru_ha(k) * 10. vmma = vmma + shallst(k) * cnv @@ -120,7 +115,7 @@ subroutine autoirr case (4) !! deep aquifer source do k = 1, nhru - if (hru_sub(k) == irrno(j)) then + if (k == irrno(j)) then cnv = 0. cnv = hru_ha(k) * 10. vmma = vmma + deepst(k) * cnv diff --git a/src/header.f b/src/header.f index e0bb5b1..a0ea253 100644 --- a/src/header.f +++ b/src/header.f @@ -93,7 +93,7 @@ subroutine header & " WTMPdegc"," Salt1 "," Salt2 ", & " Salt3 "," Salt4 "," Salt5 ", & " Salt6 "," Salt7 "," Salt8 ", - & " Salt9 "," Salt10 "," SAR ", + & " Salt9 "," Salt10 "," SAR ", & " EC "/) !! numbers printed to VB interface reach output file diff --git a/src/main.f b/src/main.f index 3d117dd..4927861 100644 --- a/src/main.f +++ b/src/main.f @@ -48,7 +48,7 @@ program main use parm implicit none !prog = "SWAT Dec 1 VER 2022/Merge Rev 663/Rev 687" - prog = "SWAT Jan 20 VER 2023/Rev 688" + prog = "SWAT May 10 VER 2023/Rev 689" write (*,1000) 1000 format(1x," SWAT2022 ",/, & " Rev. 688 ",/, diff --git a/src/modparm.f b/src/modparm.f index 4c814b5..afdfadd 100644 --- a/src/modparm.f +++ b/src/modparm.f @@ -2,6 +2,23 @@ module parm integer icalen real*8 :: prf_bsn + type precip_pet_moving_average !!for tropical plant growth + integer :: trop = 0 !! |1=tropical growth ?moisture driven + integer :: peren = 0 !! |0=annual crop; 1=perennial + integer :: mce = 0 !! |my current element -day in the p/pet arrays + integer :: mon_seas = 0 !! |0=not monsoon season; 1=in monsoon + integer :: ndays_mon = 0 !! |number of days in the monsoon period + integer :: curday_mon = 0 !! |current day into the monsoon period + integer :: ndays = 30 !! |number of days for precip/pet moving average + real :: trig = 0.5 !!mm/mm |precip/pet ratio to trigger plant/restart + real :: precip_sum !!mm |sum of precip during moving average period + real :: pet_sum !!mm |sum of pet during moving average period + real, dimension (:), allocatable :: precip !!mm |precip dimensioned by ndays + real, dimension (:), allocatable :: pet !!mm |pet dimensioned by ndays + real :: rto = 0 !! |sum of precip/sum of pet + end type precip_pet_moving_average + type (precip_pet_moving_average), dimension(:),allocatable :: ppet !dimensioned by subbasin + !! srin - co2 (EPA) real*8 :: co2_x2, co2_x @@ -11,7 +28,7 @@ module parm & tmp_spr2, tmp_fal1, tmp_fal2 real*8 :: wtmp - + real*8, dimension (12) :: pcpmm real*8, dimension (:), allocatable :: alph_e real*8, dimension (:), allocatable :: co_p, surlag, cdn, nperco real*8, dimension (:), allocatable :: cmn, phoskd, psp, sdnco diff --git a/src/readfcst.f b/src/readfcst.f index 1171326..c883dca 100644 --- a/src/readfcst.f +++ b/src/readfcst.f @@ -57,7 +57,7 @@ subroutine readfcst use parm character (len=80) :: titldum - real*8, dimension (12) :: pcpmm, pcpd + real*8, dimension (12) :: pcpd integer :: mon, mdays, j, fcstregtot diff --git a/src/readsub.f b/src/readsub.f index 531f6df..25b015f 100644 --- a/src/readsub.f +++ b/src/readsub.f @@ -454,6 +454,20 @@ subroutine readsub !!read in weather generator parameter values call readwgn plaps(i) = plaps(i) / pcpdays(i) + + !!initialize ppet for tropical growth + !! use average monthly December precip and pet from wgn - assume Jan 1 start + mo_ppet = 12 + ppet(ihru)%precip_sum = 0. + ppet(ihru)%pet_sum = 0. + do inext = 1, ppet(ihru)%ndays + !! set previous precip and pet to December ave daily precip + ppet(ihru)%precip(inext) = pcpmm(12) / 30. + ppet(ihru)%pet(inext) = pcpmm(12) / 30. + ppet(ihru)%precip_sum = ppet(ihru)%precip_sum + ppet(ihru)%precip(inext) + ppet(ihru)%pet_sum = ppet(ihru)%pet_sum + ppet(ihru)%pet(inext) + end do + !!read in subbasin impoundment parameter values call readpnd !!read in subbasin water use parameter values diff --git a/src/readwgn.f b/src/readwgn.f index ff4ff9b..a2b6e10 100644 --- a/src/readwgn.f +++ b/src/readwgn.f @@ -133,7 +133,7 @@ subroutine readwgn character (len=80) :: titldum real*8 :: xx, lattan, x1, x2, x3, tav, tmin, tmax, rain_yrs real*8 :: summx_t, summn_t, summm_p, sum, rnm2, r6, xlv, pcp - real*8, dimension (12) :: rainhhmx, rain_hhsm, pcpmm, pcpd + real*8, dimension (12) :: rainhhmx, rain_hhsm, pcpd real*8 :: tmpsoil, sffc, rndm1, dl integer :: mon, mdays, j, m1, nda, xrnd @@ -346,6 +346,23 @@ subroutine readwgn end do ffc(ihru) = sffc dormhr(ihru) = dl + !! set initial precip/pet values for tropical plant growth + !! compute potential et with Preistley-Taylor Method + mdays = ndays(13) - ndays(12) + tav = (tmpmx(12,i) + tmpmn(12,i)) / 2. + tk = tav + 273. + alb = .15 !tropical rainforests (0.05-0.15) + d = EXP(21.255 - 5304. / tk) * 5304. / tk ** 2 + gma = d / (d +.68) + ho = 23.9 * solarav(12,i) * (1. - alb) / 58.3 + aph = 1.28 + pet_dec = aph * ho * gma * 30. + ppet(j)%precip = pcpmm(12) / mdays + ppet(j)%pet = pet_dec / mdays + ppet(j)%precip_sum = pcpmm(12) + ppet(j)%pet_sum = pet_dec + ppet(j)%rto = ppet(j)%precip_sum / ppet(j)%pet_sum + end do close (114) diff --git a/src/rteinit.f b/src/rteinit.f index bb8bdb4..0b76c82 100644 --- a/src/rteinit.f +++ b/src/rteinit.f @@ -124,8 +124,8 @@ subroutine rteinit else if (bio_e(idplt(idum)) > 1.e-6) then subfr_nowtr(isb) = subfr_nowtr(isb) + hru_fr(idum) - end if - end if + end if + end if end do !! read in areas associated with .fig record files diff --git a/src/sched_mgt.f b/src/sched_mgt.f index 6bfe09e..0b495d2 100644 --- a/src/sched_mgt.f +++ b/src/sched_mgt.f @@ -61,6 +61,34 @@ subroutine sched_mgt & sol_sumsolp(j) end if + case (18) !! monsoon planting for tropical growth + igro(j) = 1 + lai_init = mgt5op(nop(j),j) + bio_init = mgt6op(nop(j),j) + ppet(j)%trop = 1 + ppet(j)%trig = mgt7op(nop(j),j) + bio_targ(j) = mgt8op(nop(j),j) * 1000. + cnop = mgt9op(nop(j),j) + ppet(j)%ndays_mon = mgt2iop(nop(j),j) + ppet(j)%peren = mgt3iop(nop(j),j) + if (curyr_mat(j) == 0) igrotree(j) = 1 + + idplt(j) = mgt1iop(nop(j),j) + ppet(j)%mon_seas = 1 + ppet(j)%curday_mon = 0 + + if (mgt4op(nop(j),j) < 700.) mgt4op(nop(j),j) = 1700. +! if (mgt4op(nop(j),j) > 5000.) mgt4op(nop(j),j) = 5000. + phu_plt(j) = mgt4op(nop(j),j) + + case (19) !! reset phenology during monsoon season for tropical perennial growth + ppet(j)%mon_seas = 1 + ppet(j)%trop = 1 + ppet(j)%peren = mgt3iop(nop(j),j) + ppet(j)%curday_mon = 0 + ppet(j)%trig = mgt7op(nop(j),j) + ppet(j)%ndays_mon = mgt2iop(nop(j),j) + case (2) !! irrigation operation irr_sc(ihru) = mgt2iop(nop(j),j) !!NUBZ @@ -234,7 +262,7 @@ subroutine sched_mgt if (imgt ==1) then write (143, 1010) subnum(j), hruno(j), iyr, i_mo, * iida, hru_km(j), " ", - * "SCHED AUTORR", phubase(j), phuacc(j), sol_sw(j), bio_ms(j), + * "SCHED AUTOIRR", phubase(j), phuacc(j), sol_sw(j), bio_ms(j), * sol_rsd(1,j), sol_sumno3(j),sol_sumsolp(j) 1010 format (a5,1x,a4,3i6,1x,e10.5,1x,2a15,7f10.2) end if diff --git a/src/stdaa.f b/src/stdaa.f index e091a54..c67fb3b 100644 --- a/src/stdaa.f +++ b/src/stdaa.f @@ -533,8 +533,7 @@ subroutine stdaa & 'AUTOPkh ',t84,'MIXEF',t90,'PRECmm',t97,'SURQGENmm',t109, & 'GWQmm',t118,'ETmm',t125,'SEDth ',t132,'NO3kgh ',t140, & 'ORGNkgh ',t148,'BIOMth',t157,'YLDth',t164,'SURQmm') - !1900 format (i7,i4,3x,a16,3x,e8.3,17f10.2) - 1900 format (i7,i4,3x,a8,3x,e8.3,17f8.2) + 1900 format (i7,i4,3x,a8,3x,e8.3,17f8.2) 2000 format (///,t17,'AVE MONTHLY BASIN VALUES',/t20,'SNOW',t46, & 'WATER',t66,'SED',/t3,'MON',t11,'RAIN',t20,'FALL',t27,'SURF Q', & t37,'LAT Q',t46,'YIELD',t58,'ET',t64,'YIELD',t75,'PET',/t11, diff --git a/src/subbasin.f b/src/subbasin.f index 3c4165e..94a0626 100644 --- a/src/subbasin.f +++ b/src/subbasin.f @@ -197,6 +197,79 @@ subroutine subbasin !! calculate soil temperature for soil layers call solt + !! compute evapotranspiration + call etpot +! if (pot_vol(j) < 1.e-6) call etact + call etact + + !! Update CMI and Precip minus PET 30 day moving sum + ppet(j)%mce = ppet(j)%mce + 1 + if (ppet(j)%mce > ppet(j)%ndays) ppet(j)%mce = 1 + ppet_mce = ppet(j)%mce + !! calculate climatic moisture index - cumulative p/pet + !! subtract the 30 day previous and add the current day precip/pet + ppet(j)%precip_sum = ppet(j)%precip_sum + precipday - ppet(j)%precip(ppet_mce) + ppet(j)%pet_sum = ppet(j)%pet_sum + pet_day - ppet(j)%pet(ppet_mce) + ppet(j)%rto = ppet(j)%precip_sum / ppet(j)%pet_sum + ppet(j)%precip(ppet_mce) = precipday + ppet(j)%pet(ppet_mce) = pet_day + if (ppet(j)%trop > 0) then + if (ppet(j)%mon_seas > 0) then + ppet(j)%curday_mon = ppet(j)%curday_mon + 1 + !! plant if last day of monsoon and not planted + if (ppet(j)%curday_mon == ppet(j)%ndays_mon) then + ppet(j)%mon_seas = 0 + if (ppet(j)%peren == 0) then + !! annual planting + call plantop + if (imgt == 1) then + write (143, 1000) subnum(j), hruno(j), iyr, i_mo, iida, + & hru_km(j),cpnm(idplt(j))," PLANT", phubase(j), phuacc(j), + & sol_sw(j),bio_ms(j), sol_rsd(1,j),sol_sumno3(j), + & sol_sumsolp(j) + end if + else + !! perennial phenology reset + igro(j) = 1 + idorm(j) = 0 + phuacc(j) = 0. + if (imgt == 1) then + write (143, 1000) subnum(j), hruno(j), iyr, i_mo, iida, + & hru_km(j),cpnm(idplt(j))," PHENO-RESET", phubase(j), phuacc(j), + & sol_sw(j),bio_ms(j), sol_rsd(1,j),sol_sumno3(j), + & sol_sumsolp(j) + end if + end if + end if + !else + !! check if ppet ratio is exceeded + if (ppet(j)%rto > ppet(j)%trig) then + ppet(j)%mon_seas = 0 + if (ppet(j)%peren == 0) then + !! annual planting + call plantop + if (imgt == 1) then + write (143, 1000) subnum(j), hruno(j), iyr, i_mo, iida, + & hru_km(j),cpnm(idplt(j))," PLANT", phubase(j), phuacc(j), + & sol_sw(j),bio_ms(j), sol_rsd(1,j),sol_sumno3(j), + & sol_sumsolp(j) + end if + else + !! perennial phenology reset + igro(j) = 1 + idorm(j) = 0 + phuacc(j) = 0. + if (imgt == 1) then + write (143, 1000) subnum(j), hruno(j), iyr, i_mo, iida, + & hru_km(j),cpnm(idplt(j))," PHENO-RESET", phubase(j), phuacc(j), + & sol_sw(j),bio_ms(j), sol_rsd(1,j),sol_sumno3(j), + & sol_sumsolp(j) + end if + end if + end if + end if + end if + call surface !! compute effective rainfall (amount that percs into soil) @@ -207,16 +280,16 @@ subroutine subbasin !! perform management operations if (yr_skip(j) == 0) call operatn +!!!! Srin's irrigation source by each application changes + irrsc(j) = irr_sca(j) + irrno(j) = irr_noa(j) +!!!! Srin's irrigation source by each application changes + if (irrsc(j) > 2) call autoirr !! perform soil water routing call percmain - !! compute evapotranspiration - call etpot -! if (pot_vol(j) < 1.e-6) call etact - call etact - !! compute water table depth using climate drivers call wattable @@ -251,8 +324,9 @@ subroutine subbasin !! compute crop growth call plantmod - !! check for dormancy - if (igro(j) == 1) call dormant + !! check for dormancy - not for tropical plants + if (igro(j) == 1 .and. ppet(j)%trop == 0) call dormant + !! compute actual ET for day in HRU etday = ep_day + es_day + canev diff --git a/src/writed.f b/src/writed.f index 4126db2..9a14bfa 100644 --- a/src/writed.f +++ b/src/writed.f @@ -183,6 +183,5 @@ subroutine writed 5001 format(2i5,500f12.4) 5100 format(1x,a5,a4,1x,i4,1x,i3,1x,250(e16.4,1x)) 5200 format(i7,i9,i6,i5,1x,e9.4,f12.3,f7.1,f14.3) -!!6200 format(i5,13f7.2,2f5.2,1x,5f8.2) 6200 format(i5,15f8.2,1x,4f8.2) end \ No newline at end of file From 41ba3dd60ed2c6175f0c0af5da618ccfdf708926 Mon Sep 17 00:00:00 2001 From: crazyzlj Date: Tue, 1 Aug 2023 19:25:12 +0800 Subject: [PATCH 2/2] Specify free line length to subbasin.f --- src/CMakeLists.txt | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c314e0e..c4e4c19 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -31,11 +31,16 @@ endif() ## Special sources that have fixed length of 72 set(LEN72_SRCS "grow" "tran") +## Special sources that have long length of 132 but retained the suffix of "*.f" +set(LENLONG_SRCS "subbasin") foreach(F77FILE ${F77SRCS}) get_filename_component(CORENAME ${F77FILE} NAME_WE) list(FIND LEN72_SRCS ${CORENAME} _FOUND_LEN72) + list(FIND LENLONG_SRCS ${CORENAME} _FOUND_LENLONG) if(${_FOUND_LEN72} GREATER -1) set_source_files_properties(${F77FILE} PROPERTIES COMPILE_FLAGS ${FLAG_L72}) + elseif(${_FOUND_LENLONG} GREATER -1) + set_source_files_properties(${F77FILE} PROPERTIES COMPILE_FLAGS ${FLAG_LLONG}) else() set_source_files_properties(${F77FILE} PROPERTIES COMPILE_FLAGS ${FLAG_L80}) endif()