#include #include subroutine lsmdrv (londim ,latdim ,beglatspmd,endlatspmd, & pcnst ,pgcmxy ,thgcmxy ,ugcmxy , & vgcmxy ,tgcmxy ,qgcmxy ,flwdsxy , & zgcmxy ,prcxy ,prlxy ,solsxy , & sollxy ,solsdxy ,solldxy ,doalb , & nstep ,tsxy ,shxy ,lhxy , & cfxy ,tauxxy ,tauyxy ,trefxy , & asdirxy,asdifxy ,aldirxy ,aldifxy , & lwupxy ,snowxy ,eccen ,obliqr , & lambm0 ,mvelpp ) #include #include #include * ------------------------ code history --------------------------- * source file: lsmdrv.F * purpose: 2-d model driver using [lsmlon] by [lsmlat] grid * date last revised: March 1996 - lsm version 1 * author: Gordon Bonan * standardized: * reviewed: * * $Id: lsmdrv.F,v 1.1 1998/11/29 05:51:58 erik Exp $ * * ----------------------------------------------------------------- * ------------------------ notes ---------------------------------- * the land surface model works by gathering all the land points on a * [lsmlon] by [lsmlat] grid into a vector of [lpt] land points. this is * then expanded into a "big" vector of [kpt] subgrid points, allowing for * up to [msub] subgrid points per land point. the [kpt] subgrid points * are processed as [numlv] "little" vectors of [numkpt] points for [numlv] * calls to subroutine lsm. this subroutine: * o generates the calendar day (1.00 -> 365.99), month (1 -> 12), * and day (1 -> 31), which are used to calculate the surface * albedos and leaf and stem areas for the next time step * o maps atmospheric fields from the [lsmlon] by [lsmlat] grid to * subgrid vectors of length [kpt] * o calls the vector land surface model code [numlv] times in * strips of [numkpt] points * o maps fields from the subgrid vectors of length [kpt] to the * [lsmlon] by [lsmlat] grid * o determines if end of history interval and writes history * and restart files if appropriate * ----------------------------------------------------------------- * ------------------- input variables ----------------------------- integer londim !atm number of longitudes integer latdim !atm number of latitudes integer beglatspmd !beg lat on proc, degenerates to 1 for nonspmd integer endlatspmd !end lat on proc, degenerates to atmlat for nonspmd integer pcnst !atm number of constituents (water is first) integer nstep !atm time step index logical doalb !true if surface albedo calculation time step real eccen !Earth's orbital eccentricity real obliqr !Earth's obliquity in radians real lambm0 !Mean longitude of perihelion at the vernal !equinox (radians) real mvelpp !Earth's moving vernal equinox longitude of !perihelion plus pi (radians) real pgcmxy(londim,beglatspmd:endlatspmd) !atm bottom level pressure (pa) real thgcmxy(londim,beglatspmd:endlatspmd)!atm btm lvl pot temp (kelvin) real ugcmxy(londim,beglatspmd:endlatspmd) !atm btm level zonal wind (m/s) real vgcmxy(londim,beglatspmd:endlatspmd) !atm btm lev merid wind (m/s) real tgcmxy(londim,beglatspmd:endlatspmd) !atm bottom level temp (kelvin) real qgcmxy(londim,beglatspmd:endlatspmd) !atm btm lvl spc humidity(kg/kg) real zgcmxy(londim,beglatspmd:endlatspmd) !atm btm lev hght above srf (m) real prcxy(londim,latdim) !conv precip rate (mm h2o/s) real prlxy(londim,latdim) !large-scale prec rate(mm h2o/s) real flwdsxy(londim,latdim) !downward longwave rad onto surface (w/m**2) real solsxy (londim,latdim) !vis direct beam solar rad onto srf (w/m**2) real sollxy (londim,latdim) !nir direct beam solar rad onto srf (w/m**2) real solsdxy(londim,latdim) !vis diffuse solar rad onto srf (w/m**2) real solldxy(londim,latdim) !nir diffuse solar rad onto srf(w/m**2) * ----------------------------------------------------------------- * ------------------- output variables ---------------------------- real shxy(londim,beglatspmd:endlatspmd) !sens heat flux (w/m**2) C ![+ to atm] real lhxy(londim,beglatspmd:endlatspmd) !latent heat flux (w/m**2) C ![+ to atm] real cfxy(londim,pcnst,beglatspmd:endlatspmd)!constituent fluxes C !(kg/m**2/s) [+ to atm] real tauxxy(londim,beglatspmd:endlatspmd) !zonal surf stress (kg/m/s**2) real tauyxy(londim,beglatspmd:endlatspmd) !merid surf stress (kg/m/s**2) C real tsxy (londim,latdim) !surface radiative temperature (kelvin) real lwupxy (londim,latdim) !emitted longwave radiation (w/m**2) real trefxy(londim,beglatspmd:endlatspmd) !2 m height air temp (kelvin) C real asdirxy(londim,latdim) !albedo -- visible waveband, direct real asdifxy(londim,latdim) !albedo -- visible waveband, diffuse real aldirxy(londim,latdim) !albedo -- near infrared waveband, direct real aldifxy(londim,latdim) !albedo -- near infrared waveband, diffuse real snowxy (londim,latdim) !snow, liquid water equivalent (mm) * ----------------------------------------------------------------- * ------------------- common block variables ---------------------- #include #include #include #include #include #include * ----------------------------------------------------------------- * ---------------------- local variables -------------------------- integer i,j,k,l,m !loop/array indices real wt !subgrid weight * calendar info * ----------------------------------------------------------------- integer mdcur !current day (0, ...) integer mscur !current seconds of current day (0, ..., 86400) integer mcdate !current date (yymmdd format) (e.g., 030131) integer mcsec !current seconds of current date (0, ..., 86400) integer kyr !year (0, ...) integer kmo !month (1, ..., 12) integer kda !day of month (1, ..., 31) real calday !calendar day at greenwich (1.00, ..., 365.99) * constant atmospheric co2 and o2 * ----------------------------------------------------------------- real po2 !partial pressure o2 (mol/mol) real pco2 !partial pressure co2 (mol/mol) data po2,pco2 /0.209,355.e-06/ save po2,pco2 * ----------------------------------------------------------------- * begin executable code * determine if albedo calculation is done #ifndef COUP_CSM doalb = irad.eq.1 .or.(mod(nstep,irad).eq.0 .and. nstep+1.ne.1) #endif * ---------------------------------------------------------------------- * determine if end of history interval * ---------------------------------------------------------------------- call t_startf ('lsm_hist_single') call histend (nstep) * ---------------------------------------------------------------------- * calendar information for next time step * o calday = calendar day (1.00 -> 365.99) for cosine solar zenith angle * calday is based on greenwich time * o kmo = month (1 -> 12) for leaf area index and stem area index * o kda = day (1 -> 31) for leaf area index and stem area index * ---------------------------------------------------------------------- call calendr(nstep+1,dtlsm ,mdbase ,msbase ,mbdate , & mbsec ,mdcur ,mscur ,mcdate ,mcsec , & calday ) kyr = mcdate/10000 kmo = mod(mcdate,10000)/100 kda = mod(mcdate,100) call t_stopf ('lsm_hist_single') * ---------------------------------------------------------------------- * map atmospheric fields to force lsm: [lsmlon] x [lsmlat] grid -> * [lpt] vector of land points -> [kpt] vector of subgrid points * ---------------------------------------------------------------------- #if ( defined SPMD ) do k = begkptspmd(beglatspmd),endkptspmd(endlatspmd) #else C$DOACROSS LOCAL(L,I,J) do k = 1, kpt !lsm vector index #endif l = klnd(k) !land point index i = ixy(l) !longitude index j = jxy(l) !latitude index tgcm(k) = tgcmxy(i,j) hgcm(k) = zgcmxy(i,j) ugcm(k) = ugcmxy(i,j) vgcm(k) = vgcmxy(i,j) qgcm(k) = qgcmxy(i,j) pgcm(k) = pgcmxy(i,j) thgcm(k) = thgcmxy(i,j) qprecc(k) = prcxy(i,j) qprecl(k) = prlxy(i,j) solad(1,k) = solsxy(i,j) solad(2,k) = sollxy(i,j) solai(1,k) = solsdxy(i,j) solai(2,k) = solldxy(i,j) firgcm(k) = flwdsxy(i,j) * ---------------------------------------------------------------------- * potential temperature, vapor pressure, air density at atm reference * height. atmospheric co2 and o2 are currently constants * ---------------------------------------------------------------------- egcm(k) = qgcm(k)*pgcm(k) / (0.622+0.378*qgcm(k)) rhogcm(k) = (pgcm(k)-0.378*egcm(k)) / (rair*tgcm(k)) co2gcm(k) = pco2*pgcm(k) o2gcm(k) = po2*pgcm(k) end do * ---------------------------------------------------------------------- * process the "big" vector of [kpt] points as [numlv] "little" * vectors of [numkpt] points. [begkpt] is the starting location of the * [numkpt] points for the "little" vector in the "big" [kpt] vector. * multitask [numlv] processes * ---------------------------------------------------------------------- C Cray directive CMIC$ DO ALL SHARED(begkpt,nstep,numkpt,dtlsm,dtsoi,doalb,pgcm,tgcm, CMIC$C firgcm,lati,long,solad,solai,co2gcm,o2gcm,qprecc,qprecl,thgcm, CMIC$C rhogcm,egcm,igs,hsno,fsno,fwet,ugcm,vgcm,hgcm,ivt,ist,isc, CMIC$C dzsoi,zsoi,tsoi,foln,tg,tv,moz,h2osoi,h2ocan,h2osno,soot,eah, CMIC$C albd,albi,albgrd,albgri,fabd,fabi,ftdd,ftid,ftii,rootb,soilc, CMIC$C fsun,htop,elai,esai,tlai,fcev,fctr,fgev,qsoil,qvege,taux,tauy, CMIC$C trad,tsa,fpsn,frm,tksol,tkdry,fsh,qvegt,fire,tsai,stemb, CMIC$C fco2,fmicr,root,watsat,hksat,smpsat,bch,watdry,watopt,csol, CMIC$C frg,kmo,kda,calday,eccen,obliqr,mvelpp,lambm0,hydro, CMIC$C pergro,qgcm,conchk,qover,qdrai,numlv) CMIC$C PRIVATE(i,k) C C SGI directive C$DOACROSS LOCAL(I,K) C C HP directive C$DIR LOOP_PARALLEL,LOOP_PRIVATE(K) #if ( defined SPMD ) if (numks.gt.0) then #endif do i = 1, numlv k = begkpt(i) call lsm (nstep ,numkpt(i) ,dtlsm ,dtsoi , & doalb ,pgcm (k) ,tgcm(k) ,qgcm(k) ,firgcm(k) , & lati(k) ,long(k) ,solad(1,k) ,solai(1,k) ,co2gcm(k) , & o2gcm(k) ,qprecc(k) ,qprecl(k) ,thgcm(k) ,rhogcm(k) , & egcm(k) ,igs(k) ,hsno(k) ,fsno(k) ,fwet(k) , & ugcm(k) ,vgcm(k) ,hgcm(k) ,ivt(k) ,ist(k) , & isc(k) ,dzsoi(1,k) ,zsoi(1,k) ,tsoi(1,k) ,foln(k) , & tg(k) ,tv(k) ,moz(k) ,h2osoi(1,k),h2ocan(k) , & h2osno(k) ,soot(k) ,eah(k) ,fsh(k) ,frg(k) , & fire(k) ,albd(1,k) ,albi(1,k) ,albgrd(1,k),albgri(1,k), & fabd(1,k) ,fabi(1,k) ,ftdd(1,k) ,ftid(1,k) ,ftii(1,k) , & fsun(k) ,htop(k) ,elai(k) ,esai(k) ,tlai(k) , & fcev(k) ,fctr(k) ,fgev(k) ,qsoil(k) ,qvege(k) , & taux(k) ,tauy(k) ,trad(k) ,tsa(k) ,stemb(k) , & fpsn(k) ,frm(k) ,fco2(k) ,tsai(k) ,rootb(k) , & fmicr(k) ,root(1,k) ,watsat(k) ,hksat(k) ,soilc(k) , & smpsat(k) ,bch(k) ,watdry(k) ,watopt(k) ,csol(k) , & tksol(k) ,tkdry(k) ,qvegt(k) ,kmo ,kda , & calday ,eccen ,obliqr ,lambm0 ,mvelpp , & hydro ,pergro ,conchk ,qover(k) ,qdrai(k) , & i ,k ) end do #if ( defined SPMD ) endif #endif * ---------------------------------------------------------------------- * update basin runoff info * ---------------------------------------------------------------------- call t_startf ('lsm_single') call basindrv(beglatspmd, endlatspmd, qover, qdrai) call t_stopf ('lsm_single') * ---------------------------------------------------------------------- * accumulate counters. write lsm history and restart files * ---------------------------------------------------------------------- call t_startf ('lsm_hist_single') call histhan (nstep, nestep, beglatspmd, endlatspmd) call t_stopf ('lsm_hist_single') * ---------------------------------------------------------------------- * return required surface fields to atmospheric model: * [kpt] vector of subgrid points -> [lpt] vector of land points -> * [lsmlon] x [lsmlat] grid. this mapping is for land points only. * non-land points are undefined. * ---------------------------------------------------------------------- * need to initialize 2-d fields to zero for subgrid averaging, but * only for land points because other surface models may have already * calculated fields for non-land points. #if ( defined SPMD ) do j = beglatspmd, endlatspmd #else C$DOACROSS LOCAL(J,I) do j = 1, lsmlat #endif do i = 1, numlon(j) if (surf2d(i,j) .gt. 0) then shxy(i,j) = 0. lhxy(i,j) = 0. cfxy(i,1,j) = 0. tauxxy(i,j) = 0. tauyxy(i,j) = 0. tsxy(i,j) = 0. trefxy(i,j) = 0. lwupxy(i,j) = 0. asdirxy(i,j) = 0. asdifxy(i,j) = 0. aldirxy(i,j) = 0. aldifxy(i,j) = 0. snowxy(i,j) = 0. end if end do do m = 2, pcnst do i = 1, numlon(j) if (surf2d(i,j) .gt. 0) cfxy(i,m,j) = 0. end do end do end do * map fields from subgrid vector, with length [kpt], to * [lsmlon] x [lsmlat] surface, weighting by subgrid fraction. * for each land point on [lsmlon] x [lsmlat] surface, process * [msub] subgrid points. if the subgrid point is not valid, a * dummy index to the subgrid vector is used and the weight is zero. do m = 1, msub !subgrid points for each land point CDIR$ IVDEP #if ( defined SPMD ) do l = beglptspmd(beglatspmd),endlptspmd(endlatspmd) #else C$DOACROSS LOCAL(L,I,J,K,WT) do l = 1, lpt !land pnt index for lsmlon x lsmlat surface #endif i = ixy(l) !longitude index j = jxy(l) !latitude index #if ( defined HP ) k = kvec(m,l) !lsm subgrid vector index wt = wsg2g(m,l) !subgrid weights #else k = kvec(l,m) !lsm subgrid vector index wt = wsg2g(l,m) !subgrid weights #endif c c The following if test is so SPMD code will not use uninitialized stack c memory values for arrays like taux. c if (wt.ne.0.) then shxy(i,j) = shxy(i,j) + fsh(k)*wt lhxy(i,j) = lhxy(i,j) + (fgev(k)+fcev(k)+fctr(k))*wt cfxy(i,1,j)= cfxy(i,1,j) + (qsoil(k)+qvege(k)+qvegt(k))*wt tauxxy(i,j) = tauxxy(i,j) + taux(k)*wt tauyxy(i,j) = tauyxy(i,j) + tauy(k)*wt tsxy(i,j) = tsxy(i,j) + trad(k)*wt trefxy(i,j) = trefxy(i,j) + tsa(k)*wt lwupxy(i,j) = lwupxy(i,j) + fire(k)*wt asdirxy(i,j) = asdirxy(i,j)+ albd(1,k)*wt asdifxy(i,j) = asdifxy(i,j)+ albi(1,k)*wt aldirxy(i,j) = aldirxy(i,j)+ albd(2,k)*wt aldifxy(i,j) = aldifxy(i,j)+ albi(2,k)*wt snowxy(i,j) = snowxy(i,j) + h2osno(k)*wt * set cfxy = co2 fluxes for pcnst > 1, multiplying by 4.4e-08 * to convert from umol co2/m**2/s to kg co2/m**2/s c cfxy(i,2,j) = cfxy(i,2,j) + fco2(k)*wt * 4.4e-08 c cfxy(i,3,j) = cfxy(i,3,j) + fpsn(k)*wt * 4.4e-08 c cfxy(i,4,j) = cfxy(i,4,j) + (frm(k)+frg(k))*wt * 4.4e-08 c cfxy(i,5,j) = cfxy(i,5,j) + fmicr(k)*wt * 4.4e-08 end if end do end do return end