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, facd, faca, 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/faca(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/facd(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