text:
       subroutine rtout

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine summarizes data for reaches

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name          |units      |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ammonian(:)   |mg N/L     |ammonia concentration in reach
!!    bury          |mg pst     |loss of pesticide from active sediment layer
!!                              |by burial
!!    ch_l2(:)      |km         |length of main channel
!!    ch_w(2,:)     |m          |average width of main channel
!!    chlora(:)     |mg chl-a/L |chlorophyll-a concentration in reach
!!    difus         |mg pst     |diffusion of pesticide from sediment to reach
!!    disolvp(:)    |mg P/L     |dissolved phosphorus concentration in reach
!!    hbactlp(:)    |# cfu/100mL|less persistent bacteria in reach/outflow
!!                              |during hour
!!    hbactp(:)     |# cfu/100mL|persistent bacteria in reach/outflow during
!!                              |hour
!!    hbod(:)       |mg O2/L    |carbonaceous biochemical oxygen demand in
!!                              |reach at end of hour
!!    hchla(:)      |mg chl-a/L |chlorophyll-a concentration in reach at end of
!!                              |hour
!!    hdisox(:)     |mg O2/L    |dissolved oxygen concentration in reach at
!!                              |end of hour
!!    hnh4(:)       |mg N/L     |ammonia concentration in reach at end of hour
!!    hno2(:)       |mg N/L     |nitrite concentration in reach at end of hour
!!    hno3(:)       |mg N/L     |nitrate concentration in reach at end of hour
!!    horgn(:)      |mg N/L     |organic nitrogen concentration in reach at
!!                              |end of hour
!!    horgp(:)      |mg P/L     |organic phosphorus concentration in reach at
!!                              |end of hour
!!    hsedyld(:)    |metric tons|sediment transported out of reach during hour
!!    hsolp(:)      |mg P/L     |dissolved phosphorus concentration in reach at
!!                              |end of hour
!!    hsolpst(:)    |mg pst/m^3 |soluble pesticide concentration in outflow
!!                              |on day
!!    hsorpst(:)    |mg pst/m^3 |sorbed pesticide concentration in outflow
!!                              |on day
!!    hrtwtr(:)     |m^3 H2O    |water leaving reach during hour
!!    ievent      |none          |rainfall/runoff code
!!                               |0 daily rainfall/curve number technique
!!                               |1 sub-daily rainfall/Green&Ampt/hourly
!!                               |  routing
!!                               |3 sub-daily rainfall/Green&Ampt/hourly routing
!!    ihout         |none       |outflow hydrograph location
!!    inum1         |none       |reach number
!!    inum2         |none       |inflow hydrograph location
!!    nitraten(:)   |mg N/L     |nitrate concentration in reach
!!    nitriten(:)   |mg N/L     |nitrite concentration in reach
!!    organicn(:)   |mg N/L     |organic nitrogen concentration in reach
!!    organicp(:)   |mg P/L     |organic phosphorus concentration in reach
!!    rch_bactlp(:) |# cfu/100ml|less persistent bacteria in reach/outflow
!!                              |at end of day
!!    rch_bactp(:)  |# cfu/100ml|persistent bacteria in reach/outflow at end
!!                              |of day
!!    rch_cbod(:)   |mg O2/L    |carbonaceous biochemical oxygen demand in
!!                              |reach
!!    rch_dox(:)    |mg O2/L    |dissolved oxygen concentration in reach
!!    reactb        |mg pst     |amount of pesticide in sediment that is lost
!!                              |through reactions
!!    reactw        |mg pst     |amount of pesticide in reach that is lost
!!                              |through reactions
!!    resuspst      |mg pst     |amount of pesticide moving from sediment to
!!                              |reach due to resuspension
!!    rnum1         |none       |fraction of inflow that is overland flow
!!    rtevp         |m^3 H2O    |evaporation from reach on day
!!    rttlc         |m^3 H2O    |transmission losses from reach on day
!!    rtwtr         |m^3 H2O    |water leaving reach on day
!!    sedpst_act(:) |m          |depth of active sediment layer in reach for
!!                              |pesticide
!!    sedpst_conc(:)|mg/(m**3)  |inital pesticide concentration in river bed
!!                              |sediment
!!    sedrch        |metric tons|sediment transported out of channel
!!                              |during time step
!!    setlpst       |mg pst     |amount of pesticide moving from water to
!!                              |sediment due to settling
!!    solpesto      |mg pst/m^3 |soluble pesticide concentration in outflow
!!                              |on day
!!    sorpesto      |mg pst/m^3 |sorbed pesticide concentration in outflow
!!                              |on day
!!    volatpst      |mg pst     |amount of pesticide in reach lost by
!!                              |volatilization
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name             |units      |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    hhvaroute(1,:,:) |deg C      |temperature
!!    hhvaroute(2,:,:) |m^3 H2O    |water
!!    hhvaroute(3,:,:) |metric tons|sediment or suspended solid load
!!    hhvaroute(4,:,:) |kg N       |organic nitrogen
!!    hhvaroute(5,:,:) |kg P       |organic phosphorus
!!    hhvaroute(6,:,:) |kg N       |nitrate
!!    hhvaroute(7,:,:) |kg P       |mineral phosphorus
!!    hhvaroute(11,:,:)|mg pst     |pesticide in solution
!!    hhvaroute(12,:,:)|mg pst     |pesticide sorbed to sediment
!!    hhvaroute(13,:,:)|kg         |chlorophyll-a
!!    hhvaroute(16,:,:)|kg         |carbonaceous biological oxygen demand
!!    hhvaroute(17,:,:)|kg         |dissolved oxygen
!!    hhvaroute(18,:,:)|# cfu/100ml|persistent bacteria
!!    hhvaroute(19,:,:)|# cfu/100ml|less persistent bacteria
!!    hhvaroute(20,:,:)|kg         |conservative metal #1
!!    hhvaroute(21,:,:)|kg         |conservative metal #2
!!    hhvaroute(22,:,:)|kg         |conservative metal #3
!!    rchdy(1,:)       |m^3/s      |flow into reach on day
!!    rchdy(2,:)       |m^3/s      |flow out of reach on day
!!    rchdy(3,:)       |m^3/s      |evaporation from reach on day
!!    rchdy(4,:)       |m^3/s      |transmission losses from reach on day
!!    rchdy(5,:)       |metric tons|sediment transported into reach on day
!!    rchdy(6,:)       |metric tons|sediment transported out of reach on day
!!    rchdy(7,:)       |mg/L       |sediment concentration in outflow
!!    rchdy(8,:)       |kg N       |organic N transported into reach on day
!!    rchdy(9,:)       |kg N       |organic N transported out of reach on day
!!    rchdy(10,:)      |kg P       |organic P transported into reach on day
!!    rchdy(11,:)      |kg P       |organic P transported out of reach on day
!!    rchdy(12,:)      |kg N       |nitrate transported into reach on day
!!    rchdy(13,:)      |kg N       |nitrate transported out of reach on day
!!    rchdy(14,:)      |kg N       |ammonia transported into reach on day
!!    rchdy(15,:)      |kg N       |ammonia transported out of reach on day
!!    rchdy(16,:)      |kg N       |nitrite transported into reach on day
!!    rchdy(17,:)      |kg N       |nitrite transported out of reach on day
!!    rchdy(18,:)      |kg P       |soluble P transported into reach on day
!!    rchdy(19,:)      |kg P       |soluble P transported out of reach on day
!!    rchdy(20,:)      |kg chla    |chlorophyll-a transported into reach on day
!!    rchdy(21,:)      |kg chla    |chlorophyll-a transported out of reach on
!!                                 |day
!!    rchdy(22,:)      |kg O2      |CBOD transported into reach on day
!!    rchdy(23,:)      |kg O2      |CBOD transported out of reach on day
!!    rchdy(24,:)      |kg O2      |dissolved oxygen transported into reach on
!!                                 |day
!!    rchdy(25,:)      |kg O2      |dissolved oxygen transported out of reach on
!!                                 |day
!!    rchdy(26,:)      |mg pst     |soluble pesticide transported into reach on
!!                                 |day
!!    rchdy(27,:)      |mg pst     |soluble pesticide transported out of reach
!!                                 |on day
!!    rchdy(28,:)      |mg pst     |sorbed pesticide transported into reach on
!!                                 |day
!!    rchdy(29,:)      |mg pst     |sorbed pesticide transported out of reach on
!!                                 |day
!!    rchdy(30,:)      |mg pst     |amount of pesticide lost through reactions
!!                                 |in reach on day
!!    rchdy(31,:)      |mg pst     |amount of pesticide lost through
!!                                 |volatilization from reach on day
!!    rchdy(32,:)      |mg pst     |amount of pesticide settling out of reach to
!!                                 |bed sediment on day
!!    rchdy(33,:)      |mg pst     |amount of pesticide resuspended from bed
!!                                 |sediment to reach on day
!!    rchdy(34,:)      |mg pst     |amount of pesticide diffusing from reach to
!!                                 |bed sediment on day
!!    rchdy(35,:)      |mg pst     |amount of pesticide in sediment layer lost
!!                                 |through reactions on day
!!    rchdy(36,:)      |mg pst     |amount of pesticide in sediment layer lost
!!                                 |through burial on day
!!    rchdy(37,:)      |mg pst     |amount of pesticide stored in river bed
!!                                 |sediments
!!    rchdy(38,:)      |# cfu/100mL|persistent bacteria transported out of reach
!!                                 |on day
!!    rchdy(39,:)      |# cfu/100mL|less persistent bacteria transported out of
!!                                 |reach on day
!!    rchdy(40,:)      |kg         |amount of conservative metal #1 transported
!!                                 |out of reach on day
!!    rchdy(41,:)      |kg         |amount of conservative metal #2 transported
!!                                 |out of reach on day
!!    rchdy(42,:)      |kg         |amount of conservative metal #3 transported
!!                                 |out of reach on day
!!    rchmono(1,:)     |m^3/s      |flow into reach during month
!!    rchmono(2,:)     |m^3/s      |flow out of reach during month
!!    rchmono(3,:)     |metric tons|sediment transported into reach during month
!!    rchmono(4,:)     |metric tons|sediment transported out of reach during
!!                                 |month
!!    rchmono(5,:)     |mg/L       |sediment concentration in outflow during
!!                                 |month
!!    rchmono(6,:)     |kg N       |organic N transported into reach during
!!                                 |month
!!    rchmono(7,:)     |kg N       |organic N transported out of reach during
!!                                 |month
!!    rchmono(8,:)     |kg P       |organic P transported into reach during
!!                                 |month
!!    rchmono(9,:)     |kg P       |organic P transported out of reach during
!!                                 |month
!!    rchmono(10,:)    |m^3/s      |evaporation from reach during month
!!    rchmono(11,:)    |m^3/s      |transmission losses from reach during month
!!    rchmono(12,:)    |kg         |conservative metal #1 transported out of
!!                                 |reach during month
!!    rchmono(13,:)    |kg         |conservative metal #2 transported out of
!!                                 |reach during month
!!    rchmono(14,:)    |kg         |conservative metal #3 transported out of
!!                                 |reach during month
!!    rchmono(15,:)    |kg N       |nitrate transported into reach during month
!!    rchmono(16,:)    |kg N       |nitrate transported out of reach during
!!                                 |month
!!    rchmono(17,:)    |kg P       |soluble P transported into reach during
!!                                 |month
!!    rchmono(18,:)    |kg P       |soluble P transported out of reach during
!!                                 |month
!!    rchmono(19,:)    |mg pst     |soluble pesticide transported into reach
!!                                 |during month
!!    rchmono(20,:)    |mg pst     |soluble pesticide transported out of reach
!!                                 |during month
!!    rchmono(21,:)    |mg pst     |sorbed pesticide transported into reach
!!                                 |during month
!!    rchmono(22,:)    |mg pst     |sorbed pesticide transported out of reach
!!                                 |during month
!!    rchmono(23,:)    |mg pst     |amount of pesticide lost through reactions
!!                                 |in reach during month
!!    rchmono(24,:)    |mg pst     |amount of pesticide lost through
!!                                 |volatilization from reach during month
!!    rchmono(25,:)    |mg pst     |amount of pesticide settling out of reach to
!!                                 |bed sediment during month
!!    rchmono(26,:)    |mg pst     |amount of pesticide resuspended from bed
!!                                 |sediment to reach during month
!!    rchmono(27,:)    |mg pst     |amount of pesticide diffusing from reach to
!!                                 |bed sediment during month
!!    rchmono(28,:)    |mg pst     |amount of pesticide in sediment layer lost
!!                                 |through reactions during month
!!    rchmono(29,:)    |mg pst     |amount of pesticide in sediment layer lost
!!                                 |through burial during month
!!    rchmono(30,:)    |kg chla    |chlorophyll-a transported into reach during
!!                                 |month
!!    rchmono(31,:)    |kg chla    |chlorophyll-a transported out of reach
!!                                 |during month
!!    rchmono(32,:)    |kg N       |ammonia transported into reach during month
!!    rchmono(33,:)    |kg N       |ammonia transported out of reach during
!!                                 |month
!!    rchmono(34,:)    |kg N       |nitrite transported into reach during month
!!    rchmono(35,:)    |kg N       |nitrite transported out of reach during
!!                                 |month
!!    rchmono(36,:)    |kg O2      |CBOD transported into reach during month
!!    rchmono(37,:)    |kg O2      |CBOD transported out of reach during month
!!    rchmono(38,:)    |kg O2      |dissolved oxygen transported into reach
!!                                 |during month
!!    rchmono(39,:)    |kg O2      |dissolved oxygen transported out of reach
!!                                 |during month
!!    rchmono(40,:)    |# cfu/100mL|persistent bacteria transported out of reach
!!                                 |during month
!!    rchmono(41,:)    |# cfu/100mL|less persistent bacteria transported out of
!!                                 |reach during month
!!    varoute(2,:)     |m^3 H2O    |water
!!    varoute(3,:)     |metric tons|sediment or suspended solid load
!!    varoute(4,:)     |kg N       |organic nitrogen
!!    varoute(5,:)     |kg P       |organic phosphorus
!!    varoute(6,:)     |kg N       |nitrate
!!    varoute(7,:)     |kg P       |soluble phosphorus
!!    varoute(11,:)    |mg pst     |pesticide in solution
!!    varoute(12,:)    |mg pst     |pesticide sorbed to sediment
!!    varoute(13,:)    |kg         |chlorophyll-a
!!    varoute(14,:)    |kg N       |ammonia
!!    varoute(15,:)    |kg N       |nitrite
!!    varoute(16,:)    |kg         |carbonaceous biological oxygen demand
!!    varoute(17,:)    |kg         |dissolved oxygen
!!    varoute(18,:)    |# cfu/100mL|persistent bacteria
!!    varoute(19,:)    |# cfu/100mL|less persistent bacteria
!!    varoute(20,:)    |kg         |conservative metal #1
!!    varoute(21,:)    |kg         |conservative metal #2
!!    varoute(22,:)    |kg         |conservative metal #3
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    bedvol      |m^3           |volume of river bed sediment
!!    ii          |none          |counter
!!    jrch        |none          |reach number
!!    sedcon      |mg/L          |sediment concentration in outflow
!!    sedpest     |mg pst        |pesticide in river bed sediment
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use parm
      implicit none

      integer :: jrch, ii
      real*8 :: sedcon, bedvol, sedpest, wtmp

      jrch = 0
      jrch = inum1

      wtmp = 0.
      wtmp = 5.0 + 0.75 * tmpav(jrch)
!! set values for routing variables
      varoute(1,ihout) = wtmp
      varoute(2,ihout) = rtwtr
      varoute(3,ihout) = sedrch
      varoute(8,ihout) = 0.
      varoute(9,ihout) = 0.
      varoute(10,ihout) = 0.
      varoute(18,ihout) = rch_bactp(jrch)
      varoute(19,ihout) = rch_bactlp(jrch)
      varoute(20,ihout) = varoute(20,inum2) * (1. - rnum1)
      varoute(21,ihout) = varoute(21,inum2) * (1. - rnum1)
      varoute(22,ihout) = varoute(22,inum2) * (1. - rnum1)
!!    sediment routing
      varoute(23,ihout) = rch_san
      varoute(24,ihout) = rch_sil
      varoute(25,ihout) = rch_cla
      varoute(26,ihout) = rch_sag
      varoute(27,ihout) = rch_lag
      varoute(28,ihout) = rch_gra
      if (ievent == 0) then
       varoute(4,ihout) = organicn(jrch) * rtwtr / 1000. + ch_orgn(jrch)
       varoute(5,ihout) = organicp(jrch) * rtwtr / 1000. + ch_orgp(jrch)
       varoute(6,ihout) = nitraten(jrch) * rtwtr / 1000.
       varoute(7,ihout) = disolvp(jrch) * rtwtr / 1000.
       varoute(11,ihout) = solpesto * rtwtr
       varoute(12,ihout) = sorpesto * rtwtr
       varoute(13,ihout) = chlora(jrch) * rtwtr / 1000.
       varoute(14,ihout) = ammonian(jrch) * rtwtr / 1000.
       varoute(15,ihout) = nitriten(jrch) * rtwtr / 1000.
       varoute(16,ihout) = rch_cbod(jrch) *  rtwtr/ 1000.
       varoute(17,ihout) = rch_dox(jrch) *  rtwtr/ 1000.
      else
       do ii = 1, nstep
          hhvaroute(2,ihout,ii) = hrtwtr(ii)     ! urban modeling by J.Jeong
          hhvaroute(3,ihout,ii) = hsedyld(ii)  ! urban modeling by J.Jeong

	! From this point, check each variables if it is simulated at subdaily interval before using the output - Jaehak 9/11/09
         hhvaroute(1,ihout,ii) = 0.
!          hhvaroute(2,ihout,ii) = hrtwtr(ii)
!          hhvaroute(3,ihout,ii) = hsedyld(ii)
          hhvaroute(4,ihout,ii) = horgn(ii) * hrtwtr(ii) / 1000.
          hhvaroute(5,ihout,ii) = horgp(ii) *  hrtwtr(ii) / 1000.
          hhvaroute(6,ihout,ii) = hno3(ii) * hrtwtr(ii) / 1000.
          hhvaroute(7,ihout,ii) = hsolp(ii) * hrtwtr(ii) / 1000.
          hhvaroute(8,ihout,ii) = 0.
          hhvaroute(9,ihout,ii) = 0.
          hhvaroute(10,ihout,ii) = 0.
          hhvaroute(11,ihout,ii) = hsolpst(ii) * hrtwtr(ii)
          hhvaroute(12,ihout,ii) = hsorpst(ii) * hrtwtr(ii)
          hhvaroute(13,ihout,ii) = hchla(ii) * hrtwtr(ii) / 1000.
          hhvaroute(14,ihout,ii) = hnh4(ii) * hrtwtr(ii) / 1000.
          hhvaroute(15,ihout,ii) = hno2(ii) * hrtwtr(ii) / 1000.
          hhvaroute(16,ihout,ii) = hbod(ii) *  hrtwtr(ii)/ 1000.
          hhvaroute(17,ihout,ii) = hdisox(ii) *  hrtwtr(ii)/ 1000.
          hhvaroute(18,ihout,ii) = hbactp(ii)
          hhvaroute(19,ihout,ii) = hbactlp(ii)
          hhvaroute(20,ihout,ii) = hhvaroute(20,inum2,ii) * (1. - rnum1)
          hhvaroute(21,ihout,ii) = hhvaroute(21,inum2,ii) * (1. - rnum1)
          hhvaroute(22,ihout,ii) = hhvaroute(22,inum2,ii) * (1. - rnum1)

          varoute(4,ihout) = varoute(4,ihout) + hhvaroute(4,ihout,ii)
          varoute(5,ihout) = varoute(5,ihout) + hhvaroute(5,ihout,ii)
          varoute(6,ihout) = varoute(6,ihout) + hhvaroute(6,ihout,ii)
          varoute(7,ihout) = varoute(7,ihout) + hhvaroute(7,ihout,ii)
          varoute(11,ihout) = varoute(11,ihout) + hhvaroute(11,ihout,ii)
          varoute(12,ihout) = varoute(12,ihout) + hhvaroute(12,ihout,ii)
          varoute(13,ihout) = varoute(13,ihout) + hhvaroute(13,ihout,ii)
          varoute(14,ihout) = varoute(14,ihout) + hhvaroute(14,ihout,ii)
          varoute(15,ihout) = varoute(15,ihout) + hhvaroute(15,ihout,ii)
          varoute(16,ihout) = varoute(16,ihout) + hhvaroute(16,ihout,ii)
          varoute(17,ihout) = varoute(17,ihout) + hhvaroute(17,ihout,ii)
        end do
      end if

!! set subdaily reach output    - by jaehak jeong for urban project, subdaily output in output.rch file
	if (ievent==1.and.iprint==3) then
	  do ii=1,nstep
!! determine sediment concentration in outflow
          sedcon = 0.
          if (hrtwtr(ii) > 0.01) then
            sedcon = hsedyld(ii) / hrtwtr(ii) * 1.e6
          else
            sedcon = 0.
          end if
          rchhr(1,jrch,ii) = hhvaroute(2,inum2,ii) * (1. - rnum1)!!flow in (m^3/s)
     &      / (idt * 60.)
          rchhr(2,jrch,ii) = hrtwtr(ii) / (idt * 60.)            !!flow out (m^3/s)
          rchhr(3,jrch,ii) = hrtevp(ii) / (idt * 60.)            !!evap (m^3/s)
          rchhr(4,jrch,ii) = hrttlc(ii) / (idt * 60.)            !!tloss (m^3/s)
          rchhr(5,jrch,ii) = hhvaroute(3,inum2,ii) * (1. - rnum1)   !!sed in (tons)
          rchhr(6,jrch,ii) = hsedyld(ii)                         !!sed out (tons)
          rchhr(7,jrch,ii) = sedcon						       !!sed conc (mg/L)
	  end do
	endif

!! determine sediment concentration in outflow
      sedcon = 0.
      if (rtwtr > 0.01) then
        sedcon = sedrch / rtwtr * 1.e6
      else
        sedcon = 0.
      end if
      if (sedcon > 200000.) sedcon = 200000.

!! determine amount of pesticide in river bed sediments
      bedvol = 0.
      sedpest = 0.
      bedvol = ch_w(2,jrch) * ch_l2(jrch) * 1000. * sedpst_act(jrch)
      sedpest = sedpst_conc(jrch) * bedvol

!! set daily reach output
      rchdy(1,jrch) = varoute(2,inum2) * (1. - rnum1) / 86400. !!flow in (m^3/s)
      rchdy(2,jrch) = rtwtr / 86400. !!flow out (m^3/s)

      rchdy(3,jrch) = rtevp / 86400.                     !!evap (m^3/s)
      rchdy(4,jrch) = rttlc / 86400.                     !!tloss (m^3/s)
      rchdy(5,jrch) = varoute(3,inum2) * (1. - rnum1)    !!sed in (tons)
      rchdy(6,jrch) = sedrch                             !!sed out (tons)
      rchdy(7,jrch) = sedcon                             !!sed conc (mg/L)
      rchdy(8,jrch) = varoute(4,inum2) * (1. - rnum1)    !!orgN in (kg N)
      rchdy(9,jrch) = varoute(4,ihout)                   !!orgN out (kg N)
      rchdy(10,jrch) = varoute(5,inum2) * (1. - rnum1)   !!orgP in (kg P)
      rchdy(11,jrch) = varoute(5,ihout)                  !!orgP out (kg P)
      rchdy(12,jrch) = varoute(6,inum2) * (1. - rnum1)   !!NO3 in (kg N)
      rchdy(13,jrch) = varoute(6,ihout)                  !!NO3 out (kg N)
      rchdy(14,jrch) = varoute(14,inum2) * (1. - rnum1)  !!NH4 in (kg)
      rchdy(15,jrch) = varoute(14,ihout)                 !!NH4 out (kg)
      rchdy(16,jrch) = varoute(15,inum2) * (1. - rnum1)  !!NO2 in (kg)
      rchdy(17,jrch) = varoute(15,ihout)                 !!NO2 out (kg)
      rchdy(18,jrch) = varoute(7,inum2) * (1. - rnum1)   !!solP in (kg P)
      rchdy(19,jrch) = varoute(7,ihout)                  !!solP out (kg P)
      rchdy(20,jrch) = varoute(13,inum2) * (1. - rnum1)  !!chl-a in (kg)
      rchdy(21,jrch) = varoute(13,ihout)                 !!chl-a out (kg)
      rchdy(22,jrch) = varoute(16,inum2) * (1. - rnum1)  !!CBOD in (kg)
      rchdy(23,jrch) = varoute(16,ihout)                 !!CBOD out (kg)
      rchdy(24,jrch) = varoute(17,inum2) * (1. - rnum1)  !!dis O2 in (kg)
      rchdy(25,jrch) = varoute(17,ihout)                 !!dis O2 out (kg)
      rchdy(26,jrch) = varoute(11,inum2) * (1. - rnum1)  !!solpst in (mg pst)
      rchdy(27,jrch) = varoute(11,ihout)                 !!solpst out (mg pst)
      rchdy(28,jrch) = varoute(12,inum2) * (1. - rnum1)  !!srbpst in (mg pst)
      rchdy(29,jrch) = varoute(12,ihout)                 !!srbpst out (mg pst)
      rchdy(30,jrch) = reactw                            !!reacted pst (mg pst)
      rchdy(31,jrch) = volatpst                          !!volatilized pst (mg)
      rchdy(32,jrch) = setlpst                           !!pst settling (mg pst)
      rchdy(33,jrch) = resuspst                          !!pst resuspension (mg)
      rchdy(34,jrch) = -difus                            !!pst diffuse to sed mg
      rchdy(35,jrch) = reactb                            !!react pst/sed (mg)
      rchdy(36,jrch) = bury                              !!pst bury (mg)
      rchdy(37,jrch) = sedpest                           !!pst in rivbed sed mg
      rchdy(38,jrch) = varoute(18,ihout)                 !!persistent bact out
      rchdy(39,jrch) = varoute(19,ihout)                 !!lpersistent bact out
      rchdy(40,jrch) = varoute(20,ihout)                 !!cmetal #1
      rchdy(41,jrch) = varoute(21,ihout)                 !!cmetal #2
      rchdy(42,jrch) = varoute(22,ihout)                 !!cmetal #3
      rchdy(60,jrch) = varoute(1,ihout)                  !!water temp deg c
!!    sediment routing
!!    Assumed all silt for default sediment routine
!!    For other sediment routing models particle size are tracked
       if (ch_eqn(jrch) .NE. 0) then
       rchdy(43,jrch) = varoute(23,inum2) * (1. - rnum1)  !!sand in
       rchdy(44,jrch) = varoute(23,ihout)                 !!sand out
       rchdy(45,jrch) = varoute(24,inum2) * (1. - rnum1)  !!silt in
       rchdy(46,jrch) = varoute(24,ihout)                 !!silt out
       rchdy(47,jrch) = varoute(25,inum2) * (1. - rnum1)  !!clay in
       rchdy(48,jrch) = varoute(25,ihout)                 !!clay out
       rchdy(49,jrch) = varoute(26,inum2) * (1. - rnum1)  !!sm ag in
       rchdy(50,jrch) = varoute(26,ihout)                 !!sm ag out
       rchdy(51,jrch) = varoute(27,inum2) * (1. - rnum1)  !!lg ag in
       rchdy(52,jrch) = varoute(27,ihout)                 !!lg ag out
       rchdy(53,jrch) = varoute(28,inum2) * (1. - rnum1)  !!gravel in
       rchdy(54,jrch) = varoute(28,ihout)                 !!gravel out
	 else
       rchdy(43,jrch) = 0.                 !!sand in
       rchdy(44,jrch) = 0.                 !!sand out
       rchdy(45,jrch) = varoute(3,inum2) * (1. - rnum1)   !!silt in
       rchdy(46,jrch) = varoute(3,ihout)                  !!silt out
       rchdy(47,jrch) = 0.                 !!clay in
       rchdy(48,jrch) = 0.                 !!clay out
       rchdy(49,jrch) = 0.                 !!sm ag in
       rchdy(50,jrch) = 0.                 !!sm ag out
       rchdy(51,jrch) = 0.                 !!lg ag in
       rchdy(52,jrch) = 0.                 !!lg ag out
       rchdy(53,jrch) = 0.                 !!gravel in
       rchdy(54,jrch) = 0.                 !!gravel out
       end if

!! summarize monthly reach output
      rchmono(1,jrch) = rchmono(1,jrch) + rchdy(1,jrch)
      rchmono(2,jrch) = rchmono(2,jrch) + rchdy(2,jrch)
      rchmono(3,jrch) = rchmono(3,jrch) + rchdy(5,jrch)
      rchmono(4,jrch) = rchmono(4,jrch) + rchdy(6,jrch)
      rchmono(5,jrch) = rchmono(5,jrch) + rchdy(7,jrch)
      rchmono(6,jrch) = rchmono(6,jrch) + rchdy(8,jrch)
      rchmono(7,jrch) = rchmono(7,jrch) + rchdy(9,jrch)
      rchmono(8,jrch) = rchmono(8,jrch) + rchdy(10,jrch)
      rchmono(9,jrch) = rchmono(9,jrch) + rchdy(11,jrch)
      rchmono(10,jrch) = rchmono(10,jrch) + rchdy(3,jrch)
      rchmono(11,jrch) = rchmono(11,jrch) + rchdy(4,jrch)
      rchmono(12,jrch) = rchmono(12,jrch) + rchdy(40,jrch)
      rchmono(13,jrch) = rchmono(13,jrch) + rchdy(41,jrch)
      rchmono(14,jrch) = rchmono(14,jrch) + rchdy(42,jrch)
      rchmono(15,jrch) = rchmono(15,jrch) + rchdy(12,jrch)
      rchmono(16,jrch) = rchmono(16,jrch) + rchdy(13,jrch)
      rchmono(17,jrch) = rchmono(17,jrch) + rchdy(18,jrch)
      rchmono(18,jrch) = rchmono(18,jrch) + rchdy(19,jrch)
      rchmono(19,jrch) = rchmono(19,jrch) + rchdy(26,jrch)
      rchmono(20,jrch) = rchmono(20,jrch) + rchdy(27,jrch)
      rchmono(21,jrch) = rchmono(21,jrch) + rchdy(28,jrch)
      rchmono(22,jrch) = rchmono(22,jrch) + rchdy(29,jrch)
      rchmono(23,jrch) = rchmono(23,jrch) + rchdy(30,jrch)
      rchmono(24,jrch) = rchmono(24,jrch) + rchdy(31,jrch)
      rchmono(25,jrch) = rchmono(25,jrch) + rchdy(32,jrch)
      rchmono(26,jrch) = rchmono(26,jrch) + rchdy(33,jrch)
      rchmono(27,jrch) = rchmono(27,jrch) + rchdy(34,jrch)
      rchmono(28,jrch) = rchmono(28,jrch) + rchdy(35,jrch)
      rchmono(29,jrch) = rchmono(29,jrch) + rchdy(36,jrch)
      rchmono(30,jrch) = rchmono(30,jrch) + rchdy(20,jrch)
      rchmono(31,jrch) = rchmono(31,jrch) + rchdy(21,jrch)
      rchmono(32,jrch) = rchmono(32,jrch) + rchdy(14,jrch)
      rchmono(33,jrch) = rchmono(33,jrch) + rchdy(15,jrch)
      rchmono(34,jrch) = rchmono(34,jrch) + rchdy(16,jrch)
      rchmono(35,jrch) = rchmono(35,jrch) + rchdy(17,jrch)
      rchmono(36,jrch) = rchmono(36,jrch) + rchdy(22,jrch)
      rchmono(37,jrch) = rchmono(37,jrch) + rchdy(23,jrch)
      rchmono(38,jrch) = rchmono(38,jrch) + rchdy(24,jrch)
      rchmono(39,jrch) = rchmono(39,jrch) + rchdy(25,jrch)
      rchmono(40,jrch) = rchmono(40,jrch) + rchdy(38,jrch)
      rchmono(41,jrch) = rchmono(41,jrch) + rchdy(39,jrch)
!!    sediment routing
       rchmono(42,jrch) = rchmono(42,jrch) + rchdy(43,jrch)
       rchmono(43,jrch) = rchmono(43,jrch) + rchdy(44,jrch)
       rchmono(44,jrch) = rchmono(44,jrch) + rchdy(45,jrch)
       rchmono(45,jrch) = rchmono(45,jrch) + rchdy(46,jrch)
       rchmono(46,jrch) = rchmono(46,jrch) + rchdy(47,jrch)
       rchmono(47,jrch) = rchmono(47,jrch) + rchdy(48,jrch)
       rchmono(48,jrch) = rchmono(48,jrch) + rchdy(49,jrch)
       rchmono(49,jrch) = rchmono(49,jrch) + rchdy(50,jrch)
       rchmono(50,jrch) = rchmono(50,jrch) + rchdy(51,jrch)
       rchmono(51,jrch) = rchmono(51,jrch) + rchdy(52,jrch)
       rchmono(52,jrch) = rchmono(52,jrch) + rchdy(53,jrch)
       rchmono(53,jrch) = rchmono(53,jrch) + rchdy(54,jrch)
       rchmono(54,jrch) = rchmono(54,jrch) + rchdy(55,jrch)
       rchmono(55,jrch) = rchmono(55,jrch) + rchdy(56,jrch)
       rchmono(56,jrch) = rchmono(56,jrch) + rchdy(57,jrch)
       rchmono(57,jrch) = rchmono(57,jrch) + rchdy(58,jrch)
	   rchmono(58,jrch) = rchmono(58,jrch) + rchdy(59,jrch)

      return
      end

max_stars_repo_name: GISWAT/erosion-sediment

max_issues_repo_name: GISWAT/erosion-sediment

max_forks_repo_name: GISWAT/erosion-sediment


MODULE module_sf_sfclayrev

 REAL , PARAMETER :: VCONVC=1.
 REAL , PARAMETER :: CZO=0.0185
 REAL , PARAMETER :: OZO=1.59E-5

 REAL, DIMENSION(0:1000 ),SAVE :: PSIMTB,PSIHTB

CONTAINS

   SUBROUTINE SFCLAYREV(U3D,V3D,T3D,QV3D,P3D,dz8w, &
                     CP,G,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM, &
                     ZNT,UST,PBLH,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
                     FM,FH, &
                     XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,QGH,QSFC,RMOL, &
                     U10,V10,TH2,T2,Q2, &
                     GZ1OZ0,WSPD,BR,ISFFLX,DX, &
                     SVP1,SVP2,SVP3,SVPT0,EP1,EP2, &
                     KARMAN,EOMEG,STBOLT, &
                     P1000mb, &
                     ids,ide, jds,jde, kds,kde, &
                     ims,ime, jms,jme, kms,kme, &
                     its,ite, jts,jte, kts,kte, &
                     ustm,ck,cka,cd,cda,isftcflx,iz0tlnd )

      IMPLICIT NONE
      INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
                                        ims,ime, jms,jme, kms,kme, &
                                        its,ite, jts,jte, kts,kte
      INTEGER, INTENT(IN ) :: ISFFLX
      REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
      REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT
      REAL, INTENT(IN ) :: P1000mb
      REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
                INTENT(IN ) :: dz8w
      REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
                INTENT(IN ) :: QV3D, &
                                                              P3D, &
                                                              T3D
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(IN ) :: MAVAIL, &
                                                             PBLH, &
                                                            XLAND, &
                                                              TSK
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(OUT ) :: U10, &
                                                              V10, &
                                                              TH2, &
                                                               T2, &
                                                               Q2, &
                                                             QSFC
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(INOUT) :: REGIME, &
                                                              HFX, &
                                                              QFX, &
                                                               LH, &
                                                          MOL,RMOL
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(INOUT) :: GZ1OZ0,WSPD,BR, &
                                                  PSIM,PSIH,FM,FH
      REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
                INTENT(IN ) :: U3D, &
                                                              V3D
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(IN ) :: PSFC
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(INOUT) :: ZNT, &
                                                              ZOL, &
                                                              UST, &
                                                              CPM, &
                                                             CHS2, &
                                                             CQS2, &
                                                              CHS
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(INOUT) :: FLHC,FLQC
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(INOUT) :: &
                                                              QGH
      REAL, INTENT(IN ) :: CP,G,ROVCP,R,XLV,DX
      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(OUT) :: ck,cka,cd,cda,ustm
      INTEGER, OPTIONAL, INTENT(IN ) :: ISFTCFLX, IZ0TLND
      REAL, DIMENSION( its:ite ) :: U1D, &
                                                              V1D, &
                                                             QV1D, &
                                                              P1D, &
                                                              T1D
      REAL, DIMENSION( its:ite ) :: dz8w1d
      INTEGER :: I,J
      DO J=jts,jte
        DO i=its,ite
          dz8w1d(I) = dz8w(i,1,j)
        ENDDO
        DO i=its,ite
           U1D(i) =U3D(i,1,j)
           V1D(i) =V3D(i,1,j)
           QV1D(i)=QV3D(i,1,j)
           P1D(i) =P3D(i,1,j)
           T1D(i) =T3D(i,1,j)
        ENDDO
        CALL SFCLAYREV1D(J,U1D,V1D,T1D,QV1D,P1D,dz8w1d, &
                CP,G,ROVCP,R,XLV,PSFC(ims,j),CHS(ims,j),CHS2(ims,j),&
                CQS2(ims,j),CPM(ims,j),PBLH(ims,j), RMOL(ims,j), &
                ZNT(ims,j),UST(ims,j),MAVAIL(ims,j),ZOL(ims,j), &
                MOL(ims,j),REGIME(ims,j),PSIM(ims,j),PSIH(ims,j), &
                FM(ims,j),FH(ims,j), &
                XLAND(ims,j),HFX(ims,j),QFX(ims,j),TSK(ims,j), &
                U10(ims,j),V10(ims,j),TH2(ims,j),T2(ims,j), &
                Q2(ims,j),FLHC(ims,j),FLQC(ims,j),QGH(ims,j), &
                QSFC(ims,j),LH(ims,j), &
                GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),ISFFLX,DX, &
                SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT, &
                P1000mb, &
                ids,ide, jds,jde, kds,kde, &
                ims,ime, jms,jme, kms,kme, &
                its,ite, jts,jte, kts,kte &
                ,isftcflx,iz0tlnd, &
                USTM(ims,j),CK(ims,j),CKA(ims,j), &
                CD(ims,j),CDA(ims,j) &
                                                                   )
      ENDDO
   END SUBROUTINE SFCLAYREV
   SUBROUTINE SFCLAYREV1D(J,UX,VX,T1D,QV1D,P1D,dz8w1d, &
                     CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM,PBLH,RMOL, &
                     ZNT,UST,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH,FM,FH,&
                     XLAND,HFX,QFX,TSK, &
                     U10,V10,TH2,T2,Q2,FLHC,FLQC,QGH, &
                     QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX, &
                     SVP1,SVP2,SVP3,SVPT0,EP1,EP2, &
                     KARMAN,EOMEG,STBOLT, &
                     P1000mb, &
                     ids,ide, jds,jde, kds,kde, &
                     ims,ime, jms,jme, kms,kme, &
                     its,ite, jts,jte, kts,kte, &
                     isftcflx, iz0tlnd, &
                     ustm,ck,cka,cd,cda )
      IMPLICIT NONE
      REAL, PARAMETER :: XKA=2.4E-5
      REAL, PARAMETER :: PRT=1.
      INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
                                        ims,ime, jms,jme, kms,kme, &
                                        its,ite, jts,jte, kts,kte, &
                                        J
      INTEGER, INTENT(IN ) :: ISFFLX
      REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
      REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT
      REAL, INTENT(IN ) :: P1000mb
      REAL, DIMENSION( ims:ime ) , &
                INTENT(IN ) :: MAVAIL, &
                                                             PBLH, &
                                                            XLAND, &
                                                              TSK
      REAL, DIMENSION( ims:ime ) , &
                INTENT(IN ) :: PSFCPA
      REAL, DIMENSION( ims:ime ) , &
                INTENT(INOUT) :: REGIME, &
                                                              HFX, &
                                                              QFX, &
                                                         MOL,RMOL
      REAL, DIMENSION( ims:ime ) , &
                INTENT(INOUT) :: GZ1OZ0,WSPD,BR, &
                                                  PSIM,PSIH,FM,FH
      REAL, DIMENSION( ims:ime ) , &
                INTENT(INOUT) :: ZNT, &
                                                              ZOL, &
                                                              UST, &
                                                              CPM, &
                                                             CHS2, &
                                                             CQS2, &
                                                              CHS
      REAL, DIMENSION( ims:ime ) , &
                INTENT(INOUT) :: FLHC,FLQC
      REAL, DIMENSION( ims:ime ) , &
                INTENT(INOUT) :: &
                                                              QGH
      REAL, DIMENSION( ims:ime ) , &
                INTENT(OUT) :: U10,V10, &
                                                TH2,T2,Q2,QSFC,LH
      REAL, INTENT(IN ) :: CP,G,ROVCP,R,XLV,DX
      REAL, DIMENSION( its:ite ), INTENT(IN ) :: dz8w1d
      REAL, DIMENSION( its:ite ), INTENT(IN ) :: UX, &
                                                               VX, &
                                                             QV1D, &
                                                              P1D, &
                                                              T1D
      REAL, OPTIONAL, DIMENSION( ims:ime ) , &
                INTENT(OUT) :: ck,cka,cd,cda,ustm
      INTEGER, OPTIONAL, INTENT(IN ) :: ISFTCFLX, IZ0TLND
      REAL, DIMENSION( its:ite ) :: ZA, &
                                                        THVX,ZQKL, &
                                                           ZQKLP1, &
                                                           THX,QX, &
                                                            PSIH2, &
                                                            PSIM2, &
                                                           PSIH10, &
                                                           PSIM10, &
                                                           DENOMQ, &
                                                          DENOMQ2, &
                                                          DENOMT2, &
                                                            WSPDI, &
                                                           GZ2OZ0, &
                                                           GZ10OZ0
      REAL, DIMENSION( its:ite ) :: &
                                                      RHOX,GOVRTH, &
                                                            TGDSA
      REAL, DIMENSION( its:ite) :: SCR3,SCR4
      REAL, DIMENSION( its:ite ) :: THGB, PSFC
      INTEGER :: KL
      INTEGER :: N,I,K,KK,L,NZOL,NK,NZOL2,NZOL10
      REAL :: PL,THCON,TVCON,E1
      REAL :: ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10
      REAL :: DTG,PSIX,DTTHX,PSIX10,PSIT,PSIT2,PSIQ,PSIQ2,PSIQ10
      REAL :: FLUXC,VSGD,Z0Q,VISC,RESTAR,CZIL,GZ0OZQ,GZ0OZT
      REAL :: ZW, ZN1, ZN2
      REAL :: zolri,zolri2,zolzz,zol0
      REAL :: psih_stable,psim_stable,psih_unstable,psim_unstable
      REAL :: zl2,zl10,z0t
      REAL, DIMENSION( its:ite ) :: pq,pq2,pq10
      KL=kte
      DO i=its,ite
         PSFC(I)=PSFCPA(I)/1000.
      ENDDO
      DO 5 I=its,ite
        TGDSA(I)=TSK(I)
        THGB(I)=TSK(I)*(P1000mb/PSFCPA(I))**ROVCP
    5 CONTINUE
   10 CONTINUE
   26 CONTINUE
      DO 30 I=its,ite
         PL=P1D(I)/1000.
         SCR3(I)=T1D(I)
         THCON=(P1000mb*0.001/PL)**ROVCP
         THX(I)=SCR3(I)*THCON
         SCR4(I)=SCR3(I)
         THVX(I)=THX(I)
         QX(I)=0.
   30 CONTINUE
      DO I=its,ite
         QGH(I)=0.
         FLHC(I)=0.
         FLQC(I)=0.
         CPM(I)=CP
      ENDDO
      DO 50 I=its,ite
         QX(I)=QV1D(I)
         TVCON=(1.+EP1*QX(I))
         THVX(I)=THX(I)*TVCON
         SCR4(I)=SCR3(I)*TVCON
   50 CONTINUE
      DO 60 I=its,ite
        E1=SVP1*EXP(SVP2*(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3))
        if(xland(i).gt.1.5.or.qsfc(i).le.0.0)QSFC(I)=EP2*E1/(PSFC(I)-E1)
        E1=SVP1*EXP(SVP2*(T1D(I)-SVPT0)/(T1D(I)-SVP3))
        PL=P1D(I)/1000.
        QGH(I)=EP2*E1/(PL-E1)
        CPM(I)=CP*(1.+0.8*QX(I))
   60 CONTINUE
   80 CONTINUE
      DO 90 I=its,ite
        ZQKLP1(I)=0.
        RHOX(I)=PSFC(I)*1000./(R*SCR4(I))
   90 CONTINUE
      DO 110 I=its,ite
           ZQKL(I)=dz8w1d(I)+ZQKLP1(I)
  110 CONTINUE
      DO 120 I=its,ite
         ZA(I)=0.5*(ZQKL(I)+ZQKLP1(I))
  120 CONTINUE
      DO 160 I=its,ite
        GOVRTH(I)=G/THX(I)
  160 CONTINUE
      DO 260 I=its,ite
        GZ1OZ0(I)=ALOG((ZA(I)+ZNT(I))/ZNT(I))
        GZ2OZ0(I)=ALOG((2.+ZNT(I))/ZNT(I))
        GZ10OZ0(I)=ALOG((10.+ZNT(I))/ZNT(I))
        IF((XLAND(I)-1.5).GE.0)THEN
          ZL=ZNT(I)
        ELSE
          ZL=0.01
        ENDIF
        WSPD(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I))
        TSKV=THGB(I)*(1.+EP1*QSFC(I))
        DTHVDZ=(THVX(I)-TSKV)
        if (xland(i).lt.1.5) then
        fluxc = max(hfx(i)/rhox(i)/cp &
              + ep1*tskv*qfx(i)/rhox(i),0.)
        VCONV = vconvc*(g/tgdsa(i)*pblh(i)*fluxc)**.33
        else
        IF(-DTHVDZ.GE.0)THEN
          DTHVM=-DTHVDZ
        ELSE
          DTHVM=0.
        ENDIF
        VCONV = 2.*SQRT(DTHVM)
        endif
        VSGD = 0.32 * (max(dx/5000.-1.,0.))**.33
        WSPD(I)=SQRT(WSPD(I)*WSPD(I)+VCONV*VCONV+vsgd*vsgd)
        WSPD(I)=AMAX1(WSPD(I),0.1)
        BR(I)=GOVRTH(I)*ZA(I)*DTHVDZ/(WSPD(I)*WSPD(I))
        IF(MOL(I).LT.0.)BR(I)=AMIN1(BR(I),0.0)
        RMOL(I)=-GOVRTH(I)*DTHVDZ*ZA(I)*KARMAN
  260 CONTINUE
      DO 320 I=its,ite
      if (br(I).gt.0) then
        if (br(I).gt.250.0) then
        zol(I)=zolri(250.0,ZA(I),ZNT(I))
        else
        zol(I)=zolri(br(I),ZA(I),ZNT(I))
        endif
      endif
      if (br(I).lt.0) then
       IF(UST(I).LT.0.001)THEN
          ZOL(I)=BR(I)*GZ1OZ0(I)
        ELSE
        if (br(I).lt.-250.0) then
        zol(I)=zolri(-250.0,ZA(I),ZNT(I))
        else
        zol(I)=zolri(br(I),ZA(I),ZNT(I))
        endif
       ENDIF
      endif
        zolzz=zol(I)*(za(I)+znt(I))/za(I)
        zol10=zol(I)*(10.+znt(I))/za(I)
        zol2=zol(I)*(2.+znt(I))/za(I)
        zol0=zol(I)*znt(I)/za(I)
        ZL2=(2.)/ZA(I)*ZOL(I)
        ZL10=(10.)/ZA(I)*ZOL(I)
        IF((XLAND(I)-1.5).LT.0.)THEN
        ZL=(0.01)/ZA(I)*ZOL(I)
        ELSE
        ZL=ZOL0
        ENDIF
        IF(BR(I).LT.0.)GOTO 310
        IF(BR(I).EQ.0.)GOTO 280
        REGIME(I)=1.
        psim(I)=psim_stable(zolzz)-psim_stable(zol0)
        psih(I)=psih_stable(zolzz)-psih_stable(zol0)
        psim10(I)=psim_stable(zol10)-psim_stable(zol0)
        psih10(I)=psih_stable(zol10)-psih_stable(zol0)
        psim2(I)=psim_stable(zol2)-psim_stable(zol0)
        psih2(I)=psih_stable(zol2)-psih_stable(zol0)
        pq(I)=psih_stable(zol(I))-psih_stable(zl)
        pq2(I)=psih_stable(zl2)-psih_stable(zl)
        pq10(I)=psih_stable(zl10)-psih_stable(zl)
        RMOL(I)=ZOL(I)/ZA(I)
        GOTO 320
  280 REGIME(I)=3.
        PSIM(I)=0.0
        PSIH(I)=PSIM(I)
        PSIM10(I)=0.
        PSIH10(I)=PSIM10(I)
        PSIM2(I)=0.
        PSIH2(I)=PSIM2(I)
        pq(I)=PSIH(I)
        pq2(I)=PSIH2(I)
        pq10(I)=0.
        ZOL(I)=0.
        RMOL(I) = ZOL(I)/ZA(I)
        GOTO 320
  310 CONTINUE
        REGIME(I)=4.
        psim(I)=psim_unstable(zolzz)-psim_unstable(zol0)
        psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
        psim10(I)=psim_unstable(zol10)-psim_unstable(zol0)
        psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
        psim2(I)=psim_unstable(zol2)-psim_unstable(zol0)
        psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
        pq(I)=psih_unstable(zol(I))-psih_unstable(zl)
        pq2(I)=psih_unstable(zl2)-psih_unstable(zl)
        pq10(I)=psih_unstable(zl10)-psih_unstable(zl)
        PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
        PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
        PSIH2(I)=AMIN1(PSIH2(I),0.9*GZ2OZ0(I))
        PSIM10(I)=AMIN1(PSIM10(I),0.9*GZ10OZ0(I))
        PSIH10(I)=AMIN1(PSIH10(I),0.9*GZ10OZ0(I))
        RMOL(I) = ZOL(I)/ZA(I)
  320 CONTINUE
      DO 330 I=its,ite
        DTG=THX(I)-THGB(I)
        PSIX=GZ1OZ0(I)-PSIM(I)
        PSIX10=GZ10OZ0(I)-PSIM10(I)
       PSIT=GZ1OZ0(I)-PSIH(I)
       PSIT2=GZ2OZ0(I)-PSIH2(I)
        IF((XLAND(I)-1.5).GE.0)THEN
          ZL=ZNT(I)
        ELSE
          ZL=0.01
        ENDIF
        PSIQ=ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-pq(I)
        PSIQ2=ALOG(KARMAN*UST(I)*2./XKA+2./ZL)-pq2(I)
        PSIQ10=ALOG(KARMAN*UST(I)*10./XKA+10./ZL)-pq10(I)
        IF ( PRESENT(ISFTCFLX) ) THEN
           IF ( ISFTCFLX.EQ.1 .AND. (XLAND(I)-1.5).GE.0. ) THEN
              Z0Q = 1.e-4
           zolzz=zol(I)*(za(I)+z0q)/za(I)
           zol10=zol(I)*(10.+z0q)/za(I)
           zol2=zol(I)*(2.+z0q)/za(I)
           zol0=zol(I)*z0q/za(I)
              if (zol(I).gt.0.) then
              psih(I)=psih_stable(zolzz)-psih_stable(zol0)
              psih10(I)=psih_stable(zol10)-psih_stable(zol0)
              psih2(I)=psih_stable(zol2)-psih_stable(zol0)
              else
                if (zol(I).eq.0) then
                psih(I)=0.
                psih10(I)=0.
                psih2(I)=0.
                else
                psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
                psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
                psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
                endif
              endif
              PSIQ=ALOG((ZA(I)+z0q)/Z0Q)-PSIH(I)
              PSIT=PSIQ
              PSIQ2=ALOG((2.+z0q)/Z0Q)-PSIH2(I)
              PSIQ10=ALOG((10.+z0q)/Z0Q)-PSIH10(I)
              PSIT2=PSIQ2
           ENDIF
          IF ( ISFTCFLX.EQ.2 .AND. (XLAND(I)-1.5).GE.0. ) THEN
              VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5
              RESTAR=UST(I)*ZNT(I)/VISC
              GZ0OZT=0.40*(7.3*SQRT(SQRT(RESTAR))*SQRT(0.71)-5.)
              z0t=znt(I)/exp(GZ0OZT)
           zolzz=zol(I)*(za(I)+z0t)/za(I)
           zol10=zol(I)*(10.+z0t)/za(I)
           zol2=zol(I)*(2.+z0t)/za(I)
           zol0=zol(I)*z0t/za(I)
              if (zol(I).gt.0.) then
              psih(I)=psih_stable(zolzz)-psih_stable(zol0)
              psih10(I)=psih_stable(zol10)-psih_stable(zol0)
              psih2(I)=psih_stable(zol2)-psih_stable(zol0)
              else
                if (zol(I).eq.0) then
                psih(I)=0.
                psih10(I)=0.
                psih2(I)=0.
                else
                psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
                psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
                psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
                endif
              endif
              PSIT=ALOG((ZA(I)+z0t)/Z0t)-PSIH(I)
              PSIT2=ALOG((2.+z0t)/Z0t)-PSIH2(I)
              GZ0OZQ=0.40*(7.3*SQRT(SQRT(RESTAR))*SQRT(0.60)-5.)
              z0q=znt(I)/exp(GZ0OZQ)
           zolzz=zol(I)*(za(I)+z0q)/za(I)
           zol10=zol(I)*(10.+z0q)/za(I)
           zol2=zol(I)*(2.+z0q)/za(I)
           zol0=zol(I)*z0q/za(I)
              if (zol(I).gt.0.) then
              psih(I)=psih_stable(zolzz)-psih_stable(zol0)
              psih10(I)=psih_stable(zol10)-psih_stable(zol0)
              psih2(I)=psih_stable(zol2)-psih_stable(zol0)
              else
                if (zol(I).eq.0) then
                psih(I)=0.
                psih10(I)=0.
                psih2(I)=0.
                else
                psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
                psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
                psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
                endif
              endif
              PSIQ=ALOG((ZA(I)+z0q)/Z0q)-PSIH(I)
              PSIQ2=ALOG((2.+z0q)/Z0q)-PSIH2(I)
              PSIQ10=ALOG((10.+z0q)/Z0q)-PSIH10(I)
           ENDIF
        ENDIF
        IF(PRESENT(ck) .and. PRESENT(cd) .and. PRESENT(cka) .and. PRESENT(cda)) THEN
           Ck(I)=(karman/psix10)*(karman/psiq10)
           Cd(I)=(karman/psix10)*(karman/psix10)
           Cka(I)=(karman/psix)*(karman/psiq)
           Cda(I)=(karman/psix)*(karman/psix)
        ENDIF
        IF ( PRESENT(IZ0TLND) ) THEN
           IF ( IZ0TLND.EQ.1 .AND. (XLAND(I)-1.5).LE.0. ) THEN
              ZL=ZNT(I)
              VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5
              RESTAR=UST(I)*ZL/VISC
              CZIL = 10.0 ** ( -0.40 * ( ZL / 0.07 ) )
              z0t=znt(I)/exp(CZIL*KARMAN*SQRT(RESTAR))
           zolzz=zol(I)*(za(I)+z0t)/za(I)
           zol10=zol(I)*(10.+z0t)/za(I)
           zol2=zol(I)*(2.+z0t)/za(I)
           zol0=zol(I)*z0t/za(I)
              if (zol(I).gt.0.) then
              psih(I)=psih_stable(zolzz)-psih_stable(zol0)
              psih10(I)=psih_stable(zol10)-psih_stable(zol0)
              psih2(I)=psih_stable(zol2)-psih_stable(zol0)
              else
                if (zol(I).eq.0) then
                psih(I)=0.
                psih10(I)=0.
                psih2(I)=0.
                else
                psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
                psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
                psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
                endif
              endif
              PSIQ=ALOG((ZA(I)+z0t)/Z0t)-PSIH(I)
              PSIQ2=ALOG((2.+z0t)/Z0t)-PSIH2(I)
              PSIT=PSIQ
              PSIT2=PSIQ2
           ENDIF
        ENDIF
        UST(I)=0.5*UST(I)+0.5*KARMAN*WSPD(I)/PSIX
        WSPDI(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I))
        IF ( PRESENT(USTM) ) THEN
        USTM(I)=0.5*USTM(I)+0.5*KARMAN*WSPDI(I)/PSIX
        ENDIF
        U10(I)=UX(I)*PSIX10/PSIX
        V10(I)=VX(I)*PSIX10/PSIX
        TH2(I)=THGB(I)+DTG*PSIT2/PSIT
        Q2(I)=QSFC(I)+(QX(I)-QSFC(I))*PSIQ2/PSIQ
        T2(I) = TH2(I)*(PSFCPA(I)/P1000mb)**ROVCP
        IF((XLAND(I)-1.5).LT.0.)THEN
          UST(I)=AMAX1(UST(I),0.001)
        ENDIF
        MOL(I)=KARMAN*DTG/PSIT/PRT
        DENOMQ(I)=PSIQ
        DENOMQ2(I)=PSIQ2
        DENOMT2(I)=PSIT2
        FM(I)=PSIX
        FH(I)=PSIT
  330 CONTINUE
  335 CONTINUE
      DO i=its,ite
        QFX(i)=0.
        HFX(i)=0.
      ENDDO
      IF (ISFFLX.EQ.0) GOTO 410
      DO 360 I=its,ite
        IF((XLAND(I)-1.5).GE.0)THEN
          ZNT(I)=CZO*UST(I)*UST(I)/G+OZO
          IF ( PRESENT(ISFTCFLX) ) THEN
             IF ( ISFTCFLX.NE.0 ) THEN
                ZW = MIN((UST(I)/1.06)**(0.3),1.0)
                ZN1 = 0.011*UST(I)*UST(I)/G + OZO
                ZN2 = 10.*exp(-9.5*UST(I)**(-.3333)) + &
                       0.11*1.5E-5/AMAX1(UST(I),0.01)
                ZNT(I)=(1.0-ZW) * ZN1 + ZW * ZN2
                ZNT(I)=MIN(ZNT(I),2.85e-3)
                ZNT(I)=MAX(ZNT(I),1.27e-7)
             ENDIF
          ENDIF
          ZL = ZNT(I)
        ELSE
          ZL = 0.01
        ENDIF
        FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/DENOMQ(I)
        DTTHX=ABS(THX(I)-THGB(I))
        IF(DTTHX.GT.1.E-5)THEN
          FLHC(I)=CPM(I)*RHOX(I)*UST(I)*MOL(I)/(THX(I)-THGB(I))
 1001 format(f8.5,2x,f12.7,2x,f12.10,2x,f12.10,2x,f13.10,2x,f12.8,f12.8,2x,i3)
        ELSE
          FLHC(I)=0.
        ENDIF
  360 CONTINUE
      DO 370 I=its,ite
        QFX(I)=FLQC(I)*(QSFC(I)-QX(I))
        QFX(I)=AMAX1(QFX(I),0.)
        LH(I)=XLV*QFX(I)
  370 CONTINUE
  390 CONTINUE
      DO 400 I=its,ite
        IF(XLAND(I)-1.5.GT.0.)THEN
          HFX(I)=FLHC(I)*(THGB(I)-THX(I))
          IF ( PRESENT(ISFTCFLX) ) THEN
             IF ( ISFTCFLX.NE.0 ) THEN
                HFX(I)=HFX(I)+RHOX(I)*USTM(I)*USTM(I)*WSPDI(I)
             ENDIF
          ENDIF
        ELSEIF(XLAND(I)-1.5.LT.0.)THEN
          HFX(I)=FLHC(I)*(THGB(I)-THX(I))
          HFX(I)=AMAX1(HFX(I),-250.)
        ENDIF
  400 CONTINUE
      DO I=its,ite
         IF((XLAND(I)-1.5).GE.0)THEN
           ZL=ZNT(I)
         ELSE
           ZL=0.01
         ENDIF
         CHS(I)=UST(I)*KARMAN/DENOMQ(I)
         CQS2(I)=UST(I)*KARMAN/DENOMQ2(I)
         CHS2(I)=UST(I)*KARMAN/DENOMT2(I)
      ENDDO
  410 CONTINUE
   END SUBROUTINE SFCLAYREV1D
   SUBROUTINE sfclayrevinit( allowed_to_read )
   END SUBROUTINE sfclayrevinit
END MODULE module_sf_sfclayrev
      function zolri(ri,z,z0)
      real right,left,midpoint
      if (ri.lt.0.)then
      left=-1200.
      right=0.
      else
      left=0.
      right=20000.
      endif
      Do While (abs(right - left) > 0.01)
      midpoint=(right+left)/2.
      a=zolri2(left,ri,z,z0)
      b=zolri2(midpoint,ri,z,z0)
      c=zolri2(right,ri,z,z0)
      if ((a * b) < 0) then
      right = midpoint
      else
            if ((c * b) < 0) then
            left = midpoint
            else
            goto 11
            endif
      endif
      enddo
 11 continue
      zolri=midpoint
      return
      end
      function zolri2(zol2,ri2,z,z0)
      zol20=zol2*z0/z
      zol3=zol2+zol20
      if (ri2.lt.0) then
      psix2=log((z+z0)/z0)-(psim_unstable(zol3)-psim_unstable(zol20))
      psih2=log((z+z0)/z0)-(psih_unstable(zol3)-psih_unstable(zol20))
      else
      psix2=log((z+z0)/z0)-(psim_stable(zol3)-psim_stable(zol20))
      psih2=log((z+z0)/z0)-(psih_stable(zol3)-psih_stable(zol20))
      endif
      zolri2=zol2*psih2/psix2**2-ri2
      return
      end
      function psim_stable(zolf)
        psim_stable=-6.1*log(zolf+(1+zolf**2.5)**(1./2.5))
      return
      end
      function psih_stable(zolf)
        psih_stable=-5.3*log(zolf+(1+zolf**1.1)**(1./1.1))
      return
      end
      function psim_unstable(zolf)
        x=(1.-16.*zolf)**.25
        psimk=2*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))-2.*ATAN(X)+2.*ATAN(1.)
        ym=(1.-10.*zolf)**0.33
        psimc=(3./2.)*log((ym**2.+ym+1.)/3.)-sqrt(3.)*ATAN((2.*ym+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.)
        psim_unstable=(psimk+zolf**2*(psimc))/(1+zolf**2.)
      return
      end
      function psih_unstable(zolf)
        y=(1.-16.*zolf)**.5
        psihk=2.*log((1+y)/2.)
        yh=(1.-34.*zolf)**0.33
        psihc=(3./2.)*log((yh**2.+yh+1.)/3.)-sqrt(3.)*ATAN((2.*yh+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.)
        psih_unstable=(psihk+zolf**2*(psihc))/(1+zolf**2.)
      return
      end