updateGR Subroutine

   subroutine updateGR
!updates soil and vegetation resistance to evaporation
!updates soil moisture
!
! based on ERA40 surface scheme
! van den Hurk 2000
! plants
! E = max(0,vegetation% * rhoa * (qa-qsat(TGR)) * 1/(rc+ra)) !no dew!!
! rc=rsmin/LAI*f1(K)*f2(WGS)*f3(D)*f4(T)
! ra,qa,qsat
! f3(D) is 1 for small plants
! bare soil
! E = max(0,(1-vegetation%) * rhoa * (qa-qsat(TGR)*hu) * (1/(rs+ra))

      use modglobal, only:nfcts, rlv, rlvi, rhoa, cp, wfc, wwilt, wsoil, rsmin, GRLAI, tEB, rsmax, lconstW
      use initfac, only:netSW, faccth, fachurel, faclGR, facwsoil, facf, facef, facT, facefi, facqsat, facdi, facain, qsat

      integer :: n
      real :: vfraction = 0.8 !fraction of GR covered in vegetation, should be made into a proper model parameter (-> modglobal)
      real :: dum
      do n = 1, nfcts

         if (faclGR(n)) then
             !facefi is actually the accumulated moisture flux, has to be converted to energy flux to calculate temperature
             !yet actually the moisture flux is needed for water budget, i.e. currently many operations cancel each other e.g. X*Lv/Lv  
             !facefi is the sum over all gridcells of a facet, thus has to be averaged by dividing by number of cells in that facet
             !units of facefi are kgW/kgA*m/s 
             facefi(n) = facefi(n)/tEB/facain(n)*rhoa*rlv !mean heat flux since last EB calculation (time average)

            if (.not. lconstW) then !remove water from soil
               facwsoil(n) = max(facwsoil(n) + facefi(n)*tEB*rlvi*facdi(n, 1), 0.) !ils13, careful this assumes water only being present in the first layer!!!
            end if
            
            !update canopy resistance used in wf_gr
            fachurel(n) = max(min(1.0, 0.5*(1.0 - cos(3.14159*facwsoil(n)/wfc))), 0.) !relative humidity above soil
            facf(n, 1) = 1./min(1.0, (0.004*netSW(n) + 0.05)/(0.81*(0.004*netSW(n) + 1))) !f1
            facf(n, 2) = 1./min(max(0.001, (facwsoil(n) - wwilt)/(wfc - wwilt)), 1.0) !f2
            !f3 drops out because it is for high vegetation only
            facf(n, 3) = 1./max((1 - 0.0016*(298-facT(n, 1))**2), 0.001) !f4
            !store resistance for plants
            facf(n, 4) = min(rsmin/GRLAI*facf(n, 1)*facf(n, 2)*facf(n, 3), rsmax)
            !store resistance for soil
            facf(n, 5) = min(rsmin*facf(n, 2), rsmax)
            dum = facT(n, 1)
            facqsat(n) = qsat(dum)
         end if
      end do

   end subroutine updateGR