! ! !----------------------------------- subroutine noahmpdrv & !................................... ! --- inputs: & ( im, km,itime,ps, u1, v1, t1, q1, soiltyp, vegtype, sigmaf, & & sfcemis, dlwflx, dswsfc, snet, delt, tg3, cm, ch, & & prsl1, prslki, zf, dry, wind, slopetyp, & & shdmin, shdmax, snoalb, sfalb, flag_iter, flag_guess, & & idveg,iopt_crs, iopt_btr, iopt_run, iopt_sfc, iopt_frz, & & iopt_inf,iopt_rad, iopt_alb, iopt_snf,iopt_tbot,iopt_stc, & & xlatin,xcoszin, iyrlen, julian,imon, & & rainn_mp,rainc_mp,snow_mp,graupel_mp,ice_mp, & ! --- in/outs: & weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc, & & canopy, trans, tsurf,zorl, & ! --- Noah MP specific & snowxy, tvxy, tgxy, canicexy,canliqxy, eahxy,tahxy,cmxy, & & chxy, fwetxy, sneqvoxy, alboldxy, qsnowxy, wslakexy, & & zwtxy, waxy, wtxy, tsnoxy,zsnsoxy, snicexy, snliqxy, & & lfmassxy, rtmassxy,stmassxy, woodxy, stblcpxy, fastcpxy, & & xlaixy,xsaixy,taussxy,smoiseq,smcwtdxy,deeprechxy,rechxy, & ! --- outputs: & sncovr1, qsurf, gflux, drain, evap, hflx, ep, runoff, & & cmm, chh, evbs, evcw, sbsno, snowc, stm, snohf, & & smcwlt2, smcref2,wet1,t2mmp,q2mp) ! ! use machine , only : kind_phys ! use date_def, only : idate use funcphys, only : fpvs use physcons, only : con_g, con_hvap, con_cp, con_jcal, & & con_eps, con_epsm1, con_fvirt, con_rd,con_hfus use module_sf_noahmplsm use module_sf_noahmp_glacier use noahmp_tables, only : isice_table, co2_table, o2_table, & & isurban_table,smcref_table,smcdry_table, & & smcmax_table,co2_table,o2_table, & & saim_table,laim_table implicit none ! --- constant parameters: real(kind=kind_phys), parameter :: cpinv = 1.0/con_cp real(kind=kind_phys), parameter :: hvapi = 1.0/con_hvap real(kind=kind_phys), parameter :: elocp = con_hvap/con_cp real(kind=kind_phys), parameter :: rhoh2o = 1000.0 real(kind=kind_phys), parameter :: convrad = con_jcal*1.e4/60.0 real(kind=kind_phys), parameter :: a2 = 17.2693882 real(kind=kind_phys), parameter :: a3 = 273.16 real(kind=kind_phys), parameter :: a4 = 35.86 real(kind=kind_phys), parameter :: a23m4 = a2*(a3-a4) ! ! --- ! real, parameter :: undefined = -1.e36 real :: dz8w = undefined real :: dx = undefined real :: qc = undefined real :: foln = 1.0 ! foliage integer :: nsoil = 4 ! hardwired to Noah integer :: nsnow = 3 ! max. snow layers integer :: ist = 1 ! soil type, 1 soil; 2 lake; 14 is water integer :: isc = 4 ! middle day soil color: soil 1-9 lightest real(kind=kind_phys), save :: zsoil(4),sldpth(4) data zsoil / -0.1, -0.4, -1.0, -2.0 / data sldpth /0.1, 0.3, 0.6, 1.0 / ! data dzs /0.1, 0.3, 0.6, 1.0 / ! ! --- input: ! integer, intent(in) :: im, km, itime,imon integer, dimension(im), intent(in) :: soiltyp, vegtype, slopetyp real (kind=kind_phys), dimension(im), intent(in) :: ps, u1, v1, & & t1, q1, sigmaf, sfcemis, dlwflx, dswsfc, snet, tg3, cm, & & ch, prsl1, prslki, wind, shdmin, shdmax, & & snoalb, sfalb, zf, & & rainn_mp,rainc_mp,snow_mp,graupel_mp,ice_mp logical, dimension(im), intent(in) :: dry real (kind=kind_phys),dimension(im),intent(in) :: xlatin,xcoszin integer, intent(in) :: idveg, iopt_crs,iopt_btr,iopt_run, & & iopt_sfc,iopt_frz,iopt_inf,iopt_rad, & & iopt_alb,iopt_snf,iopt_tbot,iopt_stc real (kind=kind_phys), intent(in) :: julian integer, intent(in) :: iyrlen real (kind=kind_phys), intent(in) :: delt logical, dimension(im), intent(in) :: flag_iter, flag_guess ! --- in/out: real (kind=kind_phys), dimension(im), intent(inout) :: weasd, & & snwdph, tskin, tprcp, srflag, canopy, trans, tsurf,zorl real (kind=kind_phys), dimension(im,km), intent(inout) :: & & smc, stc, slc real (kind=kind_phys), dimension(im), intent(inout) :: snowxy, & & tvxy,tgxy,canicexy,canliqxy,eahxy,tahxy, & & cmxy,chxy,fwetxy,sneqvoxy,alboldxy,qsnowxy, & & wslakexy,zwtxy,waxy,wtxy,lfmassxy,rtmassxy, & & stmassxy,woodxy,stblcpxy,fastcpxy,xlaixy, & & xsaixy,taussxy,smcwtdxy,deeprechxy,rechxy real (kind=kind_phys),dimension(im,-2:0),intent(inout) :: tsnoxy real (kind=kind_phys),dimension(im,-2:0),intent(inout) :: snicexy real (kind=kind_phys),dimension(im,-2:0),intent(inout) :: snliqxy real (kind=kind_phys),dimension(im,1:4), intent(inout) :: smoiseq real (kind=kind_phys),dimension(im,-2:4),intent(inout) :: zsnsoxy integer, dimension(im) :: jsnowxy real (kind=kind_phys),dimension(im) :: snodep real (kind=kind_phys),dimension(im,-2:4) :: tsnsoxy ! --- output: real (kind=kind_phys), dimension(im), intent(out) :: sncovr1, & & qsurf, gflux, drain, evap, hflx, ep, runoff, cmm, chh, & & evbs, evcw, sbsno, snowc, stm, snohf, smcwlt2, smcref2,wet1, & & t2mmp,q2mp ! --- locals: real (kind=kind_phys), dimension(im) :: rch, rho, & & q0, qs1, theta1, tv1, weasd_old, snwdph_old, & & tprcp_old, srflag_old, tskin_old, canopy_old real (kind=kind_phys), dimension(km) :: et,stsoil,smsoil, slsoil real (kind=kind_phys),dimension(im,km) :: smc_old,stc_old,slc_old real (kind=kind_phys), dimension(im) :: snow_old, tv_old,tg_old, & & canice_old,canliq_old,eah_old,tah_old,fwet_old,sneqvo_old, & & albold_old,qsnow_old,wslake_old,zwt_old,wa_old,wt_old, & & lfmass_old,rtmass_old,stmass_old,wood_old,stblcp_old, & & fastcp_old,xlai_old,xsai_old,tauss_old,smcwtd_old, & & deeprech_old,rech_old real(kind=kind_phys),dimension(im,1:4) :: smoiseq_old real(kind=kind_phys),dimension(im,-2:0) :: tsno_old real(kind=kind_phys),dimension(im,-2:0) :: snice_old real(kind=kind_phys),dimension(im,-2:0) :: snliq_old real(kind=kind_phys),dimension(im,-2:4) :: zsnso_old real(kind=kind_phys),dimension(im,-2:4) :: tsnso_old real (kind=kind_phys) :: alb, albedo, beta, chx, cmx, cmc, & & dew, drip, dqsdt2, ec, edir, ett, eta, esnow, etp, & & flx1, flx2, flx3, ffrozp, lwdn, pc, prcp, ptu, q2, & & q2sat, solnet, rc, rcs, rct, rcq, rcsoil, rsmin, & & runoff1, runoff2, runoff3, sfcspd, sfcprs, sfctmp, & & sfcems, sheat, shdfac, shdmin1d, shdmax1d, smcwlt, & & smcdry, smcref, smcmax, sneqv, snoalb1d, snowh, & & snomlt, sncovr, soilw, soilm, ssoil, tsea, th2, & & xlai, zlvl, swdn, tem, psfc,fdown,t2v,tbot real (kind=kind_phys) :: pconv,pnonc,pshcv,psnow,pgrpl,phail real (kind=kind_phys) :: lat,cosz,uu,vv,swe integer :: isnowx real (kind=kind_phys) :: tvx,tgx,canicex,canliqx,eahx, & & tahx,fwetx,sneqvox,alboldx,qsnowx,wslakex,zwtx, & & wax,wtx,lfmassx, rtmassx,stmassx, woodx,stblcpx, & & fastcpx,xlaix,xsaix,taussx,smcwtdx,deeprechx,rechx, & & qsfc1d real (kind=kind_phys), dimension(-2:0) :: tsnox, snicex, snliqx real (kind=kind_phys), dimension(-2:0) :: ficeold real (kind=kind_phys), dimension( km ) :: smoiseqx real (kind=kind_phys), dimension(-2:4) :: zsnsox real (kind=kind_phys), dimension(-2:4) :: tsnsox real (kind=kind_phys) :: z0wrf,fsa,fsr,fira,fsh,fcev,fgev, & & fctr,ecan,etran,trad,tgb,tgv,t2mv, & & t2mb,q2v,q2b,runsrf,runsub,apar, & & psn,sav,sag,fsno,nee,gpp,npp,fveg, & & qsnbot,ponding,ponding1,ponding2, & & rssun,rssha,bgap,wgap,chv,chb,emissi, & & shg,shc,shb,evg,evb,ghv,ghb,irg,irc, & & irb,tr,evc,chleaf,chuc,chv2,chb2, & & fpice,pahv,pahg,pahb,pah,co2pp,o2pp,ch2b integer :: i, k, ice, stype, vtype ,slope,nroot,couple logical :: flag(im) logical :: snowng,frzgra type(noahmp_parameters) :: parameters ! !===> ... begin here ! ! --- ... set flag for land points do i = 1, im flag(i) = dry(i) enddo ! --- ... save land-related prognostic fields for guess run do i = 1, im if (flag(i) .and. flag_guess(i)) then weasd_old(i) = weasd(i) snwdph_old(i) = snwdph(i) tskin_old(i) = tskin(i) canopy_old(i) = canopy(i) tprcp_old(i) = tprcp(i) srflag_old(i) = srflag(i) ! ! snow_old(i) = snowxy(i) tv_old(i) = tvxy(i) tg_old(i) = tgxy(i) canice_old(i) = canicexy(i) canliq_old(i) = canliqxy(i) eah_old(i) = eahxy(i) tah_old(i) = tahxy(i) fwet_old(i) = fwetxy(i) sneqvo_old(i) = sneqvoxy(i) albold_old(i) = alboldxy(i) qsnow_old(i) = qsnowxy(i) wslake_old(i) = wslakexy(i) zwt_old(i) = zwtxy(i) wa_old(i) = waxy(i) wt_old(i) = wtxy(i) lfmass_old(i) = lfmassxy(i) rtmass_old(i) = rtmassxy(i) stmass_old(i) = stmassxy(i) wood_old(i) = woodxy(i) stblcp_old(i) = stblcpxy(i) fastcp_old(i) = fastcpxy(i) xlai_old(i) = xlaixy(i) xsai_old(i) = xsaixy(i) tauss_old(i) = taussxy(i) smcwtd_old(i) = smcwtdxy(i) rech_old(i) = rechxy(i) deeprech_old(i) = deeprechxy(i) ! do k = 1, km smc_old(i,k) = smc(i,k) stc_old(i,k) = stc(i,k) slc_old(i,k) = slc(i,k) enddo ! do k = 1, km smoiseq_old(i,k) = smoiseq(i,k) enddo do k = -2,0 tsno_old(i,k) = tsnoxy(i,k) snice_old(i,k) = snicexy(i,k) snliq_old(i,k) = snliqxy(i,k) enddo do k = -2,4 zsnso_old (i,k) = zsnsoxy(i,k) enddo endif enddo ! ! call to init MP options ! ! &_________________________________________________________________ & ! --- ... initialization block do i = 1, im if (flag_iter(i) .and. flag(i)) then ep(i) = 0.0 evap (i) = 0.0 hflx (i) = 0.0 gflux(i) = 0.0 drain(i) = 0.0 canopy(i) = max(canopy(i), 0.0) evbs (i) = 0.0 evcw (i) = 0.0 trans(i) = 0.0 sbsno(i) = 0.0 snowc(i) = 0.0 snohf(i) = 0.0 endif enddo ! --- ... initialize variables do i = 1, im if (flag_iter(i) .and. flag(i)) then q0(i) = max(q1(i), 1.e-8) !* q1=specific humidity at level 1 (kg/kg) theta1(i) = t1(i) * prslki(i) !* adiabatic temp at level 1 (k) tv1(i) = t1(i) * (1.0 + con_fvirt*q0(i)) rho(i) = prsl1(i) / (con_rd * tv1(i)) qs1(i) = fpvs( t1(i) ) !* qs1=sat. humidity at level 1 (kg/kg) qs1(i) = con_eps*qs1(i) / (prsl1(i) + con_epsm1*qs1(i)) qs1(i) = max(qs1(i), 1.e-8) q0 (i) = min(qs1(i), q0(i)) if (vegtype(i) == isice_table ) then if (weasd(i) < 0.1) then weasd(i) = 0.1 endif endif endif enddo ! --- ... noah: prepare variables to run noah lsm ! 1. configuration information (c): ! ------------------------------ ! couple - couple-uncouple flag (=1: coupled, =0: uncoupled) ! ffrozp - fraction for snow-rain (1.=snow, 0.=rain, 0-1 mixed)) ! ice - sea-ice flag (=1: sea-ice, =0: land) ! dt - timestep (sec) (dt should not exceed 3600 secs) = delt ! zlvl - height (m) above ground of atmospheric forcing variables ! nsoil - number of soil layers (at least 2) ! sldpth - the thickness of each soil layer (m) do i = 1, im if (flag_iter(i) .and. flag(i)) then couple = 1 ice = 0 nsoil = km snowng = .false. frzgra = .false. ! if (srflag(i) == 1.0) then ! snow phase ! ffrozp = 1.0 ! elseif (srflag(i) == 0.0) then ! rain phase ! ffrozp = 0.0 ! endif ! use srflag directly to allow fractional rain/snow ffrozp = srflag(i) zlvl = zf(i) ! 2. forcing data (f): ! ----------------- ! lwdn - lw dw radiation flux (w/m2) ! solnet - net sw radiation flux (dn-up) (w/m2) ! sfcprs - pressure at height zlvl above ground (pascals) ! prcp - precip rate (kg m-2 s-1) ! sfctmp - air temperature (k) at height zlvl above ground ! th2 - air potential temperature (k) at height zlvl above ground ! q2 - mixing ratio at height zlvl above ground (kg kg-1) lat = xlatin(i) ! in radian cosz = xcoszin(i) lwdn = dlwflx(i) !..downward lw flux at sfc in w/m2 swdn = dswsfc(i) !..downward sw flux at sfc in w/m2 solnet = snet(i) !..net sw rad flx (dn-up) at sfc in w/m2 sfcems = sfcemis(i) sfctmp = t1(i) sfcprs = prsl1(i) psfc = ps(i) prcp = rhoh2o * tprcp(i) / delt if (prcp > 0.0) then if (ffrozp > 0.0) then ! rain/snow flag, one condition is enough? snowng = .true. qsnowxy(i) = ffrozp * prcp/10.0 !still use rho water? else if (sfctmp <= 275.15) frzgra = .true. endif endif th2 = theta1(i) q2 = q0(i) ! 3. other forcing (input) data (i): ! ------------------------------ ! sfcspd - wind speed (m s-1) at height zlvl above ground ! q2sat - sat mixing ratio at height zlvl above ground (kg kg-1) ! dqsdt2 - slope of sat specific humidity curve at t=sfctmp (kg kg-1 k-1) uu = u1(i) vv = v1(i) sfcspd = wind(i) q2sat = qs1(i) dqsdt2 = q2sat * a23m4/(sfctmp-a4)**2 ! 4. canopy/soil characteristics (s): ! -------------------------------- ! vegtyp - vegetation type (integer index) -> vtype ! soiltyp - soil type (integer index) -> stype ! slopetyp- class of sfc slope (integer index) -> slope ! shdfac - areal fractional coverage of green vegetation (0.0-1.0) ! shdmin - minimum areal fractional coverage of green vegetation -> shdmin1d ! ptu - photo thermal unit (plant phenology for annuals/crops) ! alb - backround snow-free surface albedo (fraction) ! snoalb - upper bound on maximum albedo over deep snow -> snoalb1d ! tbot - bottom soil temperature (local yearly-mean sfc air temp) vtype = vegtype(i) stype = soiltyp(i) slope = slopetyp(i) shdfac= sigmaf(i) shdmin1d = shdmin(i) shdmax1d = shdmax(i) snoalb1d = snoalb(i) alb = sfalb(i) tbot = tg3(i) ptu = 0.0 cmc = canopy(i)/1000. ! convert from mm to m tsea = tsurf(i) ! clu_q2m_iter snowh = snwdph(i) * 0.001 ! convert from mm to m sneqv = weasd(i) * 0.001 ! convert from mm to m ! 5. history (state) variables (h): ! ------------------------------ ! cmc - canopy moisture content (m) ! t1 - ground/canopy/snowpack) effective skin temperature (k) -> tsea ! stc(nsoil) - soil temp (k) -> stsoil ! smc(nsoil) - total soil moisture content (volumetric fraction) -> smsoil ! sh2o(nsoil)- unfrozen soil moisture content (volumetric fraction) -> slsoil ! snowh - actual snow depth (m) ! sneqv - liquid water-equivalent snow depth (m) ! albedo - surface albedo including snow effect (unitless fraction) ! ch - surface exchange coefficient for heat and moisture (m s-1) -> chx ! cm - surface exchange coefficient for momentum (m s-1) -> cmx isnowx = nint(snowxy(i)) tvx = tvxy(i) tgx = tgxy(i) canliqx = canliqxy(i) !in mm canicex = canicexy(i) eahxy(i) = (ps(i)*q2)/(0.622+q2) ! use q0 to reinit; eahx = eahxy(i) tahx = tahxy(i) co2pp = co2_table * sfcprs o2pp = o2_table * sfcprs fwetx = fwetxy(i) sneqvox = sneqvoxy(i) alboldx = alboldxy(i) qsnowx = qsnowxy(i) wslakex = wslakexy(i) zwtx = zwtxy(i) wax = waxy(i) wtx = waxy(i) do k = -2,0 tsnsoxy(i,k) = tsnoxy(i,k) enddo do k = 1,4 tsnsoxy(i,k) = stc(i,k) enddo do k = -2,0 snicex(k) = snicexy(i,k) ! in k/m3; mm snliqx(k) = snliqxy(i,k) ! in k/m3; mm tsnox (k) = tsnoxy(i,k) ficeold(k) = 0.0 ! derived if (snicex(k) > 0.0 ) then ficeold(k) = snicex(k) /(snicex(k)+snliqx(k)) endif enddo do k = -2, km zsnsox(k) = zsnsoxy(i,k) tsnsox(k) = tsnsoxy(i,k) enddo lfmassx = lfmassxy(i) rtmassx = rtmassxy(i) stmassx = stmassxy(i) woodx = woodxy(i) stblcpx = stblcpxy(i) fastcpx = fastcpxy(i) xsaix = xsaixy(i) xlaix = xlaixy(i) taussx = taussxy(i) qsfc1d = undefined ! derive later, it is an in/out? swe = weasd(i) do k = 1, km smoiseqx(k) = smoiseq(i,k) enddo smcwtdx = smcwtdxy(i) rechx = rechxy(i) deeprechx = deeprechxy(i) !-- ! the optional details for precip !-- ! pconv = 0. ! convective - may introduce later ! pnonc = (1 - ffrozp) * prcp ! large scale total in mm/s; ! pshcv = 0. ! psnow = ffrozp * prcp /10.0 ! snow = qsnowx? ! pgrpl = 0. ! phail = 0. pnonc = rainn_mp(i) pconv = rainc_mp(i) pshcv = 0. psnow = snow_mp(i) pgrpl = graupel_mp(i) phail = ice_mp(i) ! !-- old ! do k = 1, km ! stsoil(k) = stc(i,k) smsoil(k) = smc(i,k) slsoil(k) = slc(i,k) enddo snowh = snwdph(i) * 0.001 ! convert from mm to m if (swe /= 0.0 .and. snowh == 0.0) then snowh = 10.0 * swe /1000.0 endif chx = chxy(i) ! maybe chxy cmx = cmxy(i) chh(i) = ch(i) * wind(i) * rho(i) cmm(i) = cm(i) * wind(i) call transfer_mp_parameters(vtype,stype,slope,isc,parameters) call noahmp_options(idveg ,iopt_crs,iopt_btr,iopt_run,iopt_sfc, & & iopt_frz,iopt_inf,iopt_rad,iopt_alb,iopt_snf,iopt_tbot,iopt_stc) if ( vtype == isice_table ) then ice = -1 tbot = min(tbot,263.15) call noahmp_options_glacier & & (idveg ,iopt_crs ,iopt_btr, iopt_run ,iopt_sfc ,iopt_frz, & & iopt_inf ,iopt_rad ,iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc ) call noahmp_glacier ( & & i ,1 ,cosz ,nsnow ,nsoil ,delt , & ! in : time/space/model-related & sfctmp ,sfcprs ,uu ,vv ,q2 ,swdn , & ! in : forcing & prcp ,lwdn ,tbot ,zlvl ,ficeold ,zsoil , & ! in : forcing & qsnowx ,sneqvox ,alboldx ,cmx ,chx ,isnowx, & ! in/out :sneqvox + alboldx -LST & swe ,smsoil ,zsnsox ,snowh ,snicex ,snliqx , & ! in/out : sneqvx + snowhx are avgd & tgx ,tsnsox ,slsoil ,taussx ,qsfc1d , & ! in/out : & fsa ,fsr ,fira ,fsh ,fgev ,ssoil , & ! out : & trad ,edir ,runsrf ,runsub ,sag ,albedo , & ! out : albedo is surface albedo & qsnbot ,ponding ,ponding1,ponding2,t2mb ,q2b , & ! out : & emissi ,fpice ,ch2b ,esnow ) ! ! in/out and outs ! fsno = 1.0 tvx = undefined canicex = undefined canliqx = undefined eahx = undefined tahx = undefined fwetx = undefined wslakex = undefined zwtx = undefined wax = undefined wtx = undefined lfmassx = undefined rtmassx = undefined stmassx = undefined woodx = undefined stblcpx = undefined fastcpx = undefined xlaix = undefined xsaix = undefined smcwtdx = 0.0 rechx = 0.0 deeprechx = 0.0 do k = 1,4 smoiseqx(k) = smsoil(k) enddo fctr = undefined fcev = undefined z0wrf = 0.002 eta = fgev t2mmp(i) = t2mb q2mp(i) = q2b ! ! Non-glacial case ! else ice = 0 ! write(*,*)'tsnsox(1)=',tsnsox,'tgx=',tgx call noahmp_sflx (parameters ,& & i , 1 , lat , iyrlen , julian , cosz ,& ! in : time/space-related & delt , dx , dz8w , nsoil , zsoil , nsnow ,& ! in : model configuration & shdfac , shdmax1d, vtype , ice , ist ,& ! in : vegetation/soil & smoiseqx ,& ! in & sfctmp , sfcprs , psfc , uu , vv , q2 ,& ! in : forcing & qc , swdn , lwdn ,& ! in : forcing & pconv , pnonc , pshcv , psnow , pgrpl , phail ,& ! in : forcing & tbot , co2pp , o2pp , foln , ficeold , zlvl ,& ! in : forcing & alboldx , sneqvox ,& ! in/out : & tsnsox , slsoil , smsoil , tahx , eahx , fwetx ,& ! in/out : & canliqx , canicex , tvx , tgx , qsfc1d , qsnowx ,& ! in/out : & isnowx , zsnsox , snowh , swe , snicex , snliqx ,& ! in/out : & zwtx , wax , wtx , wslakex , lfmassx , rtmassx,& ! in/out : & stmassx , woodx , stblcpx , fastcpx , xlaix ,xsaix ,& ! in/out : & cmx , chx , taussx ,& ! in/out : & smcwtdx ,deeprechx, rechx ,& ! in/out : & z0wrf ,& ! out & fsa , fsr , fira , fsh , ssoil , fcev ,& ! out : & fgev , fctr , ecan , etran , edir , trad ,& ! out : & tgb , tgv , t2mv , t2mb , q2v , q2b ,& ! out : & runsrf , runsub , apar , psn , sav , sag ,& ! out : & fsno , nee , gpp , npp , fveg , albedo ,& ! out : & qsnbot , ponding , ponding1, ponding2, rssun , rssha ,& ! out : & bgap , wgap , chv , chb , emissi ,& ! out : & shg , shc , shb , evg , evb , ghv ,&! out : & ghb , irg , irc , irb , tr , evc ,& ! out : & chleaf , chuc , chv2 , chb2 , fpice , pahv ,& ! out & pahg , pahb , pah , esnow ) eta = fcev + fgev + fctr ! the flux w/m2 t2mmp(i) = t2mv*fveg+t2mb*(1-fveg) q2mp(i) = q2v*fveg+q2b*(1-fveg) endif ! glacial split ends ! ! mp in/out ! snowxy (i) = float(isnowx) tvxy (i) = tvx tgxy (i) = tgx canliqxy (i) = canliqx canicexy (i) = canicex eahxy (i) = eahx tahxy (i) = tahx cmxy (i) = cmx chxy (i) = chx fwetxy (i) = fwetx sneqvoxy (i) = sneqvox alboldxy (i) = alboldx qsnowxy (i) = qsnowx wslakexy (i) = wslakex zwtxy (i) = zwtx waxy (i) = wax wtxy (i) = wtx do k = -2,0 tsnoxy (i,k) = tsnsox(k) snicexy (i,k) = snicex (k) snliqxy (i,k) = snliqx (k) enddo do k = -2,4 zsnsoxy (i,k) = zsnsox(k) enddo lfmassxy (i) = lfmassx rtmassxy (i) = rtmassx stmassxy (i) = stmassx woodxy (i) = woodx stblcpxy (i) = stblcpx fastcpxy (i) = fastcpx xlaixy (i) = xlaix xsaixy (i) = xsaix taussxy (i) = taussx rechxy (i) = rechx deeprechxy(i) = deeprechx smcwtdxy(i) = smcwtdx smoiseq(i,1:4) = smoiseqx(1:4) ! ! generic in/outs ! do k = 1, km stc(i,k) = tsnsox(k) smc(i,k) = smsoil(k) slc(i,k) = slsoil(k) enddo canopy(i) = canicex + canliqx weasd(i) = swe snwdph(i) = snowh * 1000.0 ! write(*,*) 'swe,snowh,can' ! write (*,*) swe,snowh*1000.0,canopy(i) ! smcmax = smcmax_table(stype) smcref = smcref_table(stype) smcwlt = smcdry_table(stype) ! ! outs ! wet1(i) = smsoil(1) / smcmax smcwlt2(i) = smcwlt smcref2(i) = smcref runoff(i) = runsrf drain(i) = runsub zorl(i) = z0wrf * 100.0 sncovr1(i) = fsno snowc (i) = fsno sbsno(i) = esnow gflux(i) = -1.0*ssoil hflx(i) = fsh evbs(i) = fgev evcw(i) = fcev trans(i) = fctr evap(i) = eta ! write(*,*) 'vtype, stype are',vtype,stype ! write(*,*) 'fsh,gflx,eta',fsh,ssoil,eta ! write(*,*) 'esnow,runsrf,runsub',esnow,runsrf,runsub ! write(*,*) 'evbs,evcw,trans',fgev,fcev,fctr ! write(*,*) 'snowc',fsno tsurf(i) = trad stm(i) = (0.1*smsoil(1)+0.3*smsoil(2)+0.6*smsoil(3)+ & & 1.0*smsoil(4))*1000.0 ! unit conversion from m to kg m-2 ! snohf (i) = qsnbot * con_hfus ! only part of it but is diagnostic ! write(*,*) 'snohf',snohf(i) fdown = fsa + lwdn t2v = sfctmp * (1.0 + 0.61*q2) ! ssoil = -1.0 *ssoil call penman (sfctmp,sfcprs,chx,t2v,th2,prcp,fdown,ssoil, & & q2,q2sat,etp,snowng,frzgra,ffrozp,dqsdt2,emissi,fsno) ep(i) = etp endif ! end if_flag_iter_and_flag_block enddo ! end do_i_loop ! --- ... compute qsurf (specific humidity at sfc) do i = 1, im if (flag_iter(i) .and. flag(i)) then rch(i) = rho(i) * con_cp * ch(i) * wind(i) qsurf(i) = q1(i) + evap(i) / (elocp * rch(i)) endif enddo do i = 1, im if (flag_iter(i) .and. flag(i)) then tem = 1.0 / rho(i) hflx(i) = hflx(i) * tem * cpinv evap(i) = evap(i) * tem * hvapi endif enddo ! --- ... restore land-related prognostic fields for guess run do i = 1, im if (flag(i)) then if (flag_guess(i)) then weasd(i) = weasd_old(i) snwdph(i) = snwdph_old(i) tskin(i) = tskin_old(i) canopy(i) = canopy_old(i) tprcp(i) = tprcp_old(i) srflag(i) = srflag_old(i) snowxy(i) = snow_old(i) tvxy(i) = tv_old(i) tgxy(i) = tg_old(i) canicexy(i) = canice_old(i) canliqxy(i) = canliq_old(i) eahxy(i) = eah_old(i) tahxy(i) = tah_old(i) fwetxy(i) = fwet_old(i) sneqvoxy(i) = sneqvo_old(i) alboldxy(i) = albold_old(i) qsnowxy(i) = qsnow_old(i) wslakexy(i) = wslake_old(i) zwtxy(i) = zwt_old(i) waxy(i) = wa_old(i) wtxy(i) = wt_old(i) lfmassxy(i) = lfmass_old(i) rtmassxy(i) = rtmass_old(i) stmassxy(i) = stmass_old(i) woodxy(i) = wood_old(i) stblcpxy(i) = stblcp_old(i) fastcpxy(i) = fastcp_old(i) xlaixy(i) = xlai_old(i) xsaixy(i) = xsai_old(i) taussxy(i) = tauss_old(i) smcwtdxy(i) = smcwtd_old(i) deeprechxy(i) = deeprech_old(i) rechxy(i) = rech_old(i) do k = 1, km smc(i,k) = smc_old(i,k) stc(i,k) = stc_old(i,k) slc(i,k) = slc_old(i,k) enddo ! do k = 1, km smoiseq(i,k) = smoiseq_old(i,k) enddo do k = -2,0 tsnoxy(i,k) = tsno_old(i,k) snicexy(i,k) = snice_old(i,k) snliqxy(i,k) = snliq_old(i,k) enddo do k = -2,4 zsnsoxy(i,k) = zsnso_old(i,k) enddo else tskin(i) = tsurf(i) endif endif enddo ! return !................................... end subroutine noahmpdrv !----------------------------------- subroutine transfer_mp_parameters (vegtype,soiltype,slopetype, & & soilcolor,parameters) use noahmp_tables use module_sf_noahmplsm implicit none integer, intent(in) :: vegtype integer, intent(in) :: soiltype integer, intent(in) :: slopetype integer, intent(in) :: soilcolor type (noahmp_parameters), intent(out) :: parameters real :: refdk real :: refkdt real :: frzk real :: frzfact parameters%iswater = iswater_table parameters%isbarren = isbarren_table parameters%isice = isice_table parameters%eblforest = eblforest_table !-----------------------------------------------------------------------& parameters%urban_flag = .false. if( vegtype == isurban_table .or. vegtype == 31 & & .or.vegtype == 32 .or. vegtype == 33) then parameters%urban_flag = .true. endif !------------------------------------------------------------------------------------------! ! transfer veg parameters !------------------------------------------------------------------------------------------! parameters%ch2op = ch2op_table(vegtype) !maximum intercepted h2o per unit lai+sai (mm) parameters%dleaf = dleaf_table(vegtype) !characteristic leaf dimension (m) parameters%z0mvt = z0mvt_table(vegtype) !momentum roughness length (m) parameters%hvt = hvt_table(vegtype) !top of canopy (m) parameters%hvb = hvb_table(vegtype) !bottom of canopy (m) parameters%den = den_table(vegtype) !tree density (no. of trunks per m2) parameters%rc = rc_table(vegtype) !tree crown radius (m) parameters%mfsno = mfsno_table(vegtype) !snowmelt m parameter () parameters%saim = saim_table(vegtype,:) !monthly stem area index, one-sided parameters%laim = laim_table(vegtype,:) !monthly leaf area index, one-sided parameters%sla = sla_table(vegtype) !single-side leaf area per kg [m2/kg] parameters%dilefc = dilefc_table(vegtype) !coeficient for leaf stress death [1/s] parameters%dilefw = dilefw_table(vegtype) !coeficient for leaf stress death [1/s] parameters%fragr = fragr_table(vegtype) !fraction of growth respiration !original was 0.3 parameters%ltovrc = ltovrc_table(vegtype) !leaf turnover [1/s] parameters%c3psn = c3psn_table(vegtype) !photosynthetic pathway: 0. = c4, 1. = c3 parameters%kc25 = kc25_table(vegtype) !co2 michaelis-menten constant at 25c (pa) parameters%akc = akc_table(vegtype) !q10 for kc25 parameters%ko25 = ko25_table(vegtype) !o2 michaelis-menten constant at 25c (pa) parameters%ako = ako_table(vegtype) !q10 for ko25 parameters%vcmx25 = vcmx25_table(vegtype) !maximum rate of carboxylation at 25c (umol co2/m**2/s) parameters%avcmx = avcmx_table(vegtype) !q10 for vcmx25 parameters%bp = bp_table(vegtype) !minimum leaf conductance (umol/m**2/s) parameters%mp = mp_table(vegtype) !slope of conductance-to-photosynthesis relationship parameters%qe25 = qe25_table(vegtype) !quantum efficiency at 25c (umol co2 / umol photon) parameters%aqe = aqe_table(vegtype) !q10 for qe25 parameters%rmf25 = rmf25_table(vegtype) !leaf maintenance respiration at 25c (umol co2/m**2/s) parameters%rms25 = rms25_table(vegtype) !stem maintenance respiration at 25c (umol co2/kg bio/s) parameters%rmr25 = rmr25_table(vegtype) !root maintenance respiration at 25c (umol co2/kg bio/s) parameters%arm = arm_table(vegtype) !q10 for maintenance respiration parameters%folnmx = folnmx_table(vegtype) !foliage nitrogen concentration when f(n)=1 (%) parameters%tmin = tmin_table(vegtype) !minimum temperature for photosynthesis (k) parameters%xl = xl_table(vegtype) !leaf/stem orientation index parameters%rhol = rhol_table(vegtype,:) !leaf reflectance: 1=vis, 2=nir parameters%rhos = rhos_table(vegtype,:) !stem reflectance: 1=vis, 2=nir parameters%taul = taul_table(vegtype,:) !leaf transmittance: 1=vis, 2=nir parameters%taus = taus_table(vegtype,:) !stem transmittance: 1=vis, 2=nir parameters%mrp = mrp_table(vegtype) !microbial respiration parameter (umol co2 /kg c/ s) parameters%cwpvt = cwpvt_table(vegtype) !empirical canopy wind parameter parameters%wrrat = wrrat_table(vegtype) !wood to non-wood ratio parameters%wdpool = wdpool_table(vegtype) !wood pool (switch 1 or 0) depending on woody or not [-] parameters%tdlef = tdlef_table(vegtype) !characteristic t for leaf freezing [k] parameters%nroot = nroot_table(vegtype) !number of soil layers with root present parameters%rgl = rgl_table(vegtype) !parameter used in radiation stress function parameters%rsmin = rs_table(vegtype) !minimum stomatal resistance [s m-1] parameters%hs = hs_table(vegtype) !parameter used in vapor pressure deficit function parameters%topt = topt_table(vegtype) !optimum transpiration air temperature [k] parameters%rsmax = rsmax_table(vegtype) !maximal stomatal resistance [s m-1] !------------------------------------------------------------------------------------------! ! transfer rad parameters !------------------------------------------------------------------------------------------! parameters%albsat = albsat_table(soilcolor,:) parameters%albdry = albdry_table(soilcolor,:) parameters%albice = albice_table parameters%alblak = alblak_table parameters%omegas = omegas_table parameters%betads = betads_table parameters%betais = betais_table parameters%eg = eg_table !------------------------------------------------------------------------------------------! ! transfer global parameters !------------------------------------------------------------------------------------------! parameters%co2 = co2_table parameters%o2 = o2_table parameters%timean = timean_table parameters%fsatmx = fsatmx_table parameters%z0sno = z0sno_table parameters%ssi = ssi_table parameters%swemx = swemx_table ! ---------------------------------------------------------------------- ! transfer soil parameters ! ---------------------------------------------------------------------- parameters%bexp = bexp_table (soiltype) parameters%dksat = dksat_table (soiltype) parameters%dwsat = dwsat_table (soiltype) parameters%f1 = f1_table (soiltype) parameters%psisat = psisat_table (soiltype) parameters%quartz = quartz_table (soiltype) parameters%smcdry = smcdry_table (soiltype) parameters%smcmax = smcmax_table (soiltype) parameters%smcref = smcref_table (soiltype) parameters%smcwlt = smcwlt_table (soiltype) ! ---------------------------------------------------------------------- ! transfer genparm parameters ! ---------------------------------------------------------------------- parameters%csoil = csoil_table parameters%zbot = zbot_table parameters%czil = czil_table frzk = frzk_table refdk = refdk_table refkdt = refkdt_table parameters%kdt = refkdt * parameters%dksat / refdk parameters%slope = slope_table(slopetype) if(parameters%urban_flag)then ! hardcoding some urban parameters for soil parameters%smcmax = 0.45 parameters%smcref = 0.42 parameters%smcwlt = 0.40 parameters%smcdry = 0.40 parameters%csoil = 3.e6 endif ! adjust frzk parameter to actual soil type: frzk * frzfact !-----------------------------------------------------------------------& if(soiltype /= 14) then frzfact = (parameters%smcmax / parameters%smcref) & & * (0.412 / 0.468) parameters%frzx = frzk * frzfact end if end subroutine transfer_mp_parameters !-----------------------------------------------------------------------& subroutine penman (sfctmp,sfcprs,ch,t2v,th2,prcp,fdown,ssoil, & & q2,q2sat,etp,snowng,frzgra,ffrozp, & & dqsdt2,emissi_in,sncovr) ! etp is calcuated right after ssoil ! ---------------------------------------------------------------------- ! subroutine penman ! ---------------------------------------------------------------------- ! calculate potential evaporation for the current point. various ! partial sums/products are also calculated and passed back to the ! calling routine for later use. ! ---------------------------------------------------------------------- implicit none logical, intent(in) :: snowng, frzgra real, intent(in) :: ch, dqsdt2,fdown,prcp,ffrozp, & & q2, q2sat, ssoil, sfcprs, sfctmp, & & t2v, th2,emissi_in,sncovr real, intent(out) :: etp real :: epsca,flx2,rch,rr,t24 real :: a, delta, fnet,rad,rho,emissi,elcp1,lvs real, parameter :: elcp = 2.4888e+3, lsubc = 2.501000e+6,cp = 1004.6 real, parameter :: lsubs = 2.83e+6, rd = 287.05, cph2o = 4.1855e+3 real, parameter :: cpice = 2.106e+3, lsubf = 3.335e5 real, parameter :: sigma = 5.6704e-8 ! ---------------------------------------------------------------------- ! executable code begins here: ! ---------------------------------------------------------------------- ! ---------------------------------------------------------------------- ! prepare partial quantities for penman equation. ! ---------------------------------------------------------------------- emissi=emissi_in ! elcp1 = (1.0-sncovr)*elcp + sncovr*elcp*lsubs/lsubc lvs = (1.0-sncovr)*lsubc + sncovr*lsubs flx2 = 0.0 delta = elcp * dqsdt2 ! delta = elcp1 * dqsdt2 t24 = sfctmp * sfctmp * sfctmp * sfctmp rr = t24 * 6.48e-8 / (sfcprs * ch) + 1.0 ! rr = emissi*t24 * 6.48e-8 / (sfcprs * ch) + 1.0 rho = sfcprs / (rd * t2v) ! ---------------------------------------------------------------------- ! adjust the partial sums / products with the latent heat ! effects caused by falling precipitation. ! ---------------------------------------------------------------------- rch = rho * cp * ch if (.not. snowng) then if (prcp > 0.0) rr = rr + cph2o * prcp / rch else ! ---- ... fractional snowfall/rainfall rr = rr + (cpice*ffrozp+cph2o*(1.-ffrozp)) & & *prcp/rch end if ! ---------------------------------------------------------------------- ! include the latent heat effects of frzng rain converting to ice on ! impact in the calculation of flx2 and fnet. ! ---------------------------------------------------------------------- ! fnet = fdown - sigma * t24- ssoil fnet = fdown - emissi*sigma * t24- ssoil if (frzgra) then flx2 = - lsubf * prcp fnet = fnet - flx2 ! ---------------------------------------------------------------------- ! finish penman equation calculations. ! ---------------------------------------------------------------------- end if rad = fnet / rch + th2- sfctmp a = elcp * (q2sat - q2) ! a = elcp1 * (q2sat - q2) epsca = (a * rr + rad * delta) / (delta + rr) etp = epsca * rch / lsubc ! etp = epsca * rch / lvs ! ---------------------------------------------------------------------- end subroutine penman