MODULE module_bl_boulac real ck_b,ceps_b,vk,temin parameter(ceps_b=1/1.4,ck_b=0.4,temin=0.0001,vk=0.4) CONTAINS subroutine boulac(frc_urb2d,idiff,flag_bep,dz8w,dt,u_phy,v_phy & ,th_phy,rho,qv_curr,hfx & ,qfx,ustar,cp,g & ,rublten,rvblten,rthblten & ,rqvblten,rqcblten & ,tke,dlk,wu,wv,wt,wq,exch_h,exch_m,pblh & ,a_u_bep,a_v_bep,a_t_bep,a_q_bep & ,a_e_bep,b_u_bep,b_v_bep & ,b_t_bep,b_q_bep,b_e_bep,dlg_bep & ,dl_u_bep,sf_bep,vl_bep & ,ids,ide, jds,jde, kds,kde & ,ims,ime, jms,jme, kms,kme & ,its,ite, jts,jte, kts,kte) implicit none INTEGER:: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte integer, INTENT(IN) :: idiff REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: DZ8W REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: qv_curr REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: RHO REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: TH_PHY REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: U_PHY REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: V_PHY REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: ustar REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: hfx REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: qfx real, INTENT(IN ) :: g,cp REAL, INTENT(IN ):: DT REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: FRC_URB2D REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: PBLH REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_u_bep REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_v_bep REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_t_bep REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_q_bep REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_e_bep REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_u_bep REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_v_bep REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_t_bep REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_q_bep REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_e_bep REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dlg_bep REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::dl_u_bep REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::sf_bep REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::vl_bep LOGICAL, INTENT(IN) :: flag_bep real, dimension (ims:ime, kms:kme, jms:jme) ::th_0 real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: exch_h real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: exch_m real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: tke real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wu real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wv real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wt real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wq real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: dlk REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RUBLTEN REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RVBLTEN REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RTHBLTEN REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RQVBLTEN REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RQCBLTEN real z1D(kms:kme) real dz1D(kms:kme) real u1D(kms:kme) real v1D(kms:kme) real th1D(kms:kme) real q1D(kms:kme) real rho1D(kms:kme) real rhoz1D(kms:kme) real tke1D(kms:kme) real th01D(kms:kme) real dlk1D(kms:kme) real dls1D(kms:kme) real exch1D(kms:kme) real sf1D(kms:kme) real vl1D(kms:kme) real a_u1D(kms:kme) real a_v1D(kms:kme) real a_t1D(kms:kme) real a_q1D(kms:kme) real a_e1D(kms:kme) real b_u1D(kms:kme) real b_v1D(kms:kme) real b_t1D(kms:kme) real b_q1D(kms:kme) real b_e1D(kms:kme) real dlg1D(kms:kme) real dl_u1D(kms:kme) real sh1D(kms:kme) real bu1D(kms:kme) real wu1D(kms:kme) real wv1D(kms:kme) real wt1D(kms:kme) real wq1D(kms:kme) real a_e(ims:ime,kms:kme,jms:jme) real b_e(ims:ime,kms:kme,jms:jme) real bu(ims:ime,kms:kme,jms:jme) real sh(ims:ime,kms:kme,jms:jme) real wrk(ims:ime) integer ix,iy,iz,id,iz_u,iw_u,ig,ir_u,ix1,iy1 real ufrac_int real vect,time_tke,hour,zzz real ustarf real summ1,summ2,summ3 save time_tke,hour do ix=its,ite do iy=jts,jte do iz=kts,kte th_0(ix,iz,iy)=300. enddo enddo enddo bu1d=0. sh1d=0. b_e1d=0. b_u1d=0. b_v1d=0. b_t1d=0. b_q1d=0. a_e1d=0. a_u1d=0. a_v1d=0. a_t1d=0. a_q1d=0. z1D=0. dz1D=0. u1D =0. v1D =0. th1D=0. q1D=0. rho1D=0. rhoz1D=0. tke1D =0. th01D =0. dlk1D =0. dls1D =0. exch1D=0. sf1D =1. vl1D =1. a_u1D =0. a_v1D =0. a_t1D =0. a_q1D =0. a_e1D =0. b_u1D =0. b_v1D =0. b_t1D =0. b_q1D =0. b_e1D =0. dlg1D =0. dl_u1D=0. sh1D =0. bu1D =0. wu1D =0. wv1D =0. wt1D =0. wq1D =0. do ix=its,ite do iy=jts,jte z1d(kts)=0. do iz= kts,kte u1D(iz)=u_phy(ix,iz,iy) v1D(iz)=v_phy(ix,iz,iy) th1D(iz)=th_phy(ix,iz,iy) q1D(iz)=qv_curr(ix,iz,iy) tke1D(iz)=tke(ix,iz,iy) rho1D(iz)=rho(ix,iz,iy) th01D(iz)=th_0(ix,iz,iy) dz1D(iz)=dz8w(ix,iz,iy) z1D(iz+1)=z1D(iz)+dz1D(iz) enddo rhoz1D(kts)=rho1D(kts) do iz=kts+1,kte rhoz1D(iz)=(rho1D(iz)*dz1D(iz-1)+rho1D(iz-1)*dz1D(iz))/(dz1D(iz-1)+dz1D(iz)) enddo rhoz1D(kte+1)=rho1D(kte) if(flag_bep)then do iz=kts,kte a_e1D(iz)=a_e_bep(ix,iz,iy) b_e1D(iz)=b_e_bep(ix,iz,iy) dlg1D(iz)=(z1D(iz)+z1D(iz+1))/2.*(1.-frc_urb2d(ix,iy))+dlg_bep(ix,iz,iy)*frc_urb2d(ix,iy) dl_u1D(iz)=dl_u_bep(ix,iz,iy) if((1.-frc_urb2d(ix,iy)).lt.1.)dl_u1D(iz)=dl_u1D(iz)/frc_urb2d(ix,iy) vl1D(iz)=vl_bep(ix,iz,iy) sf1D(iz)=sf_bep(ix,iz,iy) enddo ufrac_int=frc_urb2d(ix,iy) sf1D(kte+1)=sf_bep(ix,1,iy) else do iz=kts,kte a_e1D(iz)=0. b_e1D(iz)=0. dlg1D(iz)=(z1D(iz)+z1D(iz+1))/2. dl_u1D(iz)=0. vl1D(iz)=1. sf1D(iz)=1. enddo ufrac_int=0. sf1D(kte+1)=1. endif call boulac1D(ix,iy,ufrac_int,kms,kme,kts,kte,dz1d,z1D,dt,u1D,v1D,th1D,rho1D,rhoz1D,q1D,th01D,& tke1D,ustar(ix,iy),hfx(ix,iy),qfx(ix,iy),cp,g, & a_e1D,b_e1D, & dlg1D,dl_u1D,sf1D,vl1D,dlk1D,dls1D,exch1D,sh1D,bu1D) call pbl_height(kms,kme,kts,kte,dz1d,z1d,th1D,q1D,pblh(ix,iy)) do iz= kts,kte a_e(ix,iz,iy)=a_e1D(iz) b_e(ix,iz,iy)=b_e1D(iz) if(flag_bep)then dlg_bep(ix,iz,iy)=dlg1D(iz) endif tke(ix,iz,iy)=tke1D(iz) dlk(ix,iz,iy)=dlk1D(iz) sh(ix,iz,iy)=sh1D(iz) bu(ix,iz,iy)=bu1D(iz) exch_h(ix,iz,iy)=exch1D(iz) exch_m(ix,iz,iy)=exch1D(iz) enddo if(idiff.ne.1)then if(flag_bep)then do iz=kts,kte a_t1D(iz)=a_t_bep(ix,iz,iy) b_t1D(iz)=b_t_bep(ix,iz,iy) a_u1D(iz)=a_u_bep(ix,iz,iy) b_u1D(iz)=b_u_bep(ix,iz,iy) a_v1D(iz)=a_v_bep(ix,iz,iy) b_v1D(iz)=b_v_bep(ix,iz,iy) a_q1D(iz)=a_q_bep(ix,iz,iy) b_q1D(iz)=b_q_bep(ix,iz,iy) enddo else do iz=kts,kte a_t1D(iz)=0. b_t1D(iz)=0. a_u1D(iz)=0. b_u1D(iz)=0. a_v1D(iz)=0. b_v1D(iz)=0. a_q1D(iz)=0. b_q1D(iz)=0. enddo b_t1D(1)=hfx(ix,iy)/dz1D(1)/rho1D(1)/cp b_q1D(1)=qfx(ix,iy)/dz1D(1)/rho1D(1) a_u1D(1)=(-ustar(ix,iy)*ustar(ix,iy)/dz1D(1)/((u1D(1)**2.+v1D(1)**2.)**.5)) a_v1D(1)=(-ustar(ix,iy)*ustar(ix,iy)/dz1D(1)/((u1D(1)**2.+v1D(1)**2.)**.5)) endif call diff(kms,kme,kts,kte,1,1,dt,u1D,rho1D,rhoz1D,exch1D,a_u1D,b_u1D,sf1D,vl1D,dz1D,wu1D) call diff(kms,kme,kts,kte,1,1,dt,v1D,rho1D,rhoz1D,exch1D,a_v1D,b_v1D,sf1D,vl1D,dz1D,wv1D) call diff(kms,kme,kts,kte,1,1,dt,th1D,rho1D,rhoz1D,exch1D,a_t1D,b_t1D,sf1D,vl1D,dz1D,wt1D) call diff(kms,kme,kts,kte,1,1,dt,q1D,rho1D,rhoz1D,exch1D,a_q1D,b_q1D,sf1D,vl1D,dz1D,wq1D) do iz= kts,kte rthblten(ix,iz,iy)=rthblten(ix,iz,iy)+(th1D(iz)-th_phy(ix,iz,iy))/dt rqvblten(ix,iz,iy)=rqvblten(ix,iz,iy)+(q1D(iz)-qv_curr(ix,iz,iy))/dt rublten(ix,iz,iy)=rublten(ix,iz,iy)+(u1D(iz)-u_phy(ix,iz,iy))/dt rvblten(ix,iz,iy)=rvblten(ix,iz,iy)+(v1D(iz)-v_phy(ix,iz,iy))/dt wu(ix,iz,iy)=wu1D(iz) wv(ix,iz,iy)=wv1D(iz) wt(ix,iz,iy)=wt1D(iz) wq(ix,iz,iy)=wq1D(iz) enddo endif enddo enddo time_tke=time_tke+dt return end subroutine boulac subroutine boulac1D(ix,iy,ufrac_int,kms,kme,kts,kte,dz,z,dt,u,v,th,rho,rhoz,qa,th0,te, & ustar,hfx,qfx,cp,g, & a_e,b_e, & dlg,dl_u,sf,vl,dlk,dls,exch,sh,bu) implicit none integer iz,ix,iy integer kms,kme,kts,kte real z(kms:kme) real dz(kms:kme) real u(kms:kme) real v(kms:kme) real th(kms:kme) real rho(kms:kme) real g real cp real te(kms:kme) real qa(kms:kme) real th0(kms:kme) real dt real ustar real hfx real qfx real sf(kms:kme) real vl(kms:kme) real a_e(kms:kme) real b_e(kms:kme) real dlg(kms:kme) real dl_u(kms:kme) real ufrac_int real we(kms:kme),dwe(kms:kme) real sh(kms:kme) real bu(kms:kme) real td(kms:kme) real exch(kms:kme) real dls(kms:kme) real dlk(kms:kme) real dlu(kms:kme) real dld(kms:kme) real rhoz(kms:kme) real tstar real beta real summ1,summ2,summ3,summ4 tstar=-hfx/rho(1)/cp/ustar call dissip_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,te,dlu,dld,th,th0) call length_bougeault(ix,iy,kms,kme,kts,kte,dld,dlu,dlg,dl_u,dls,dlk) call cdtur_bougeault(ix,iy,kms,kme,kts,kte,te,z,dz,exch,dlk) call tke_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,vl,u,v,th,te,th0,ustar,tstar,exch,dls,td,sh,bu,b_e,a_e,sf,ufrac_int) call diff(kms,kme,kts,kte,1,1,dt,te,rho,rhoz,exch,a_e,b_e,sf,vl,dz,we) do iz=kts,kte if(te(iz).lt.temin) te(iz)=temin enddo return end subroutine boulac1d subroutine dissip_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,te,dlu,dld,th,th0) implicit none integer kms,kme,kts,kte,iz,izz,ix,iy real dzt,zup,beta,zup_inf,bbb,tl,zdo,zdo_sup,zzz,g real te(kms:kme),dlu(kms:kme),dld(kms:kme),dz(kms:kme) real z(kms:kme),th(kms:kme),th0(kms:kme) do iz=kts,kte zup=0. dlu(iz)=z(kte+1)-z(iz)-dz(iz)/2. zzz=0. zup_inf=0. beta=g/th0(iz) do izz=iz,kte-1 dzt=(dz(izz+1)+dz(izz))/2. zup=zup-beta*th(iz)*dzt zup=zup+beta*(th(izz+1)+th(izz))*dzt/2. zzz=zzz+dzt if(te(iz).lt.zup.and.te(iz).ge.zup_inf)then bbb=(th(izz+1)-th(izz))/dzt if(bbb.ne.0)then tl=(-beta*(th(izz)-th(iz))+sqrt( max(0.,(beta*(th(izz)-th(iz)))**2.+2.*bbb*beta*(te(iz)-zup_inf))))/bbb/beta else if(th(izz).ne.th(iz))then tl=(te(iz)-zup_inf)/(beta*(th(izz)-th(iz))) else tl=0. endif endif dlu(iz)=zzz-dzt+tl endif zup_inf=zup enddo zdo=0. zdo_sup=0. dld(iz)=z(iz)+dz(iz)/2. zzz=0. do izz=iz,kts+1,-1 dzt=(dz(izz-1)+dz(izz))/2. zdo=zdo+beta*th(iz)*dzt zdo=zdo-beta*(th(izz-1)+th(izz))*dzt/2. zzz=zzz+dzt if(te(iz).lt.zdo.and.te(iz).ge.zdo_sup)then bbb=(th(izz)-th(izz-1))/dzt if(bbb.ne.0.)then tl=(beta*(th(izz)-th(iz))+sqrt( max(0.,(beta*(th(izz)-th(iz)))**2.+2.*bbb*beta*(te(iz)-zdo_sup))))/bbb/beta else if(th(izz).ne.th(iz))then tl=(te(iz)-zdo_sup)/(beta*(th(izz)-th(iz))) else tl=0. endif endif dld(iz)=zzz-dzt+tl endif zdo_sup=zdo enddo enddo end subroutine dissip_bougeault subroutine length_bougeault(ix,iy,kms,kme,kts,kte,dld,dlu,dlg,dl_u,dls,dlk) implicit none integer kms,kme,kts,kte,iz,ix,iy real dlu(kms:kme),dld(kms:kme),dl_u(kms:kme) real dls(kms:kme),dlk(kms:kme),dlg(kms:kme) do iz=kts,kte dld(iz)=min(dld(iz),dlg(iz)) dls(iz)=sqrt(dlu(iz)*dld(iz)) dlk(iz)=min(dlu(iz),dld(iz)) if(dl_u(iz).gt.0.)then dls(iz)=1./(1./dls(iz)+1./dl_u(iz)) dlk(iz)=1./(1./dlk(iz)+1./dl_u(iz)) endif enddo return end subroutine length_bougeault subroutine cdtur_bougeault(ix,iy,kms,kme,kts,kte,te,z,dz,exch,dlk) implicit none integer iz,kms,kme,kts,kte,ix,iy real te_m,dlk_m real te(kms:kme),exch(kms:kme) real dz(kms:kme),z(kms:kme) real dlk(kms:kme) real fact exch(kts)=0. do iz=kts+1,kte te_m=(te(iz-1)*dz(iz)+te(iz)*dz(iz-1))/(dz(iz)+dz(iz-1)) dlk_m=(dlk(iz-1)*dz(iz)+dlk(iz)*dz(iz-1))/(dz(iz)+dz(iz-1)) exch(iz)=ck_b*dlk_m*sqrt(te_m) exch(iz)=max(exch(iz),0.1) enddo exch(kte+1)=0.1 return end subroutine cdtur_bougeault subroutine diff(kms,kme,kts,kte,iz1,izf,dt,co,rho,rhoz,cd,aa,bb,sf,vl,dz,fc) implicit none integer iz,iz1,izf integer kms,kme,kts,kte real dt,dzv real co(kms:kme),cd(kms:kme),dz(kms:kme) real rho(kms:kme),rhoz(kms:kme) real cddz(kms:kme+1),fc(kms:kme),df(kms:kme) real a(kms:kme,3),c(kms:kme) real sf(kms:kme),vl(kms:kme) real aa(kms:kme),bb(kms:kme) cddz(kts)=sf(kts)*rhoz(kts)*cd(kts)/dz(kts) do iz=kts+1,kte cddz(iz)=2.*sf(iz)*rhoz(iz)*cd(iz)/(dz(iz)+dz(iz-1)) enddo cddz(kte+1)=sf(kte+1)*rhoz(kte+1)*cd(kte+1)/dz(kte) do iz=kts,iz1-1 a(iz,1)=0. a(iz,2)=1. a(iz,3)=0. c(iz)=co(iz) enddo do iz=iz1,kte-izf dzv=vl(iz)*dz(iz) a(iz,1)=-cddz(iz)*dt/dzv/rho(iz) a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dzv/rho(iz)-aa(iz)*dt a(iz,3)=-cddz(iz+1)*dt/dzv/rho(iz) c(iz)=co(iz)+bb(iz)*dt enddo do iz=kte-(izf-1),kte a(iz,1)=0. a(iz,2)=1 a(iz,3)=0. c(iz)=co(iz) enddo call invert (kms,kme,kts,kte,a,c,co) do iz=kts,iz1 fc(iz)=0. enddo do iz=iz1+1,kte fc(iz)=-(cddz(iz)*(co(iz)-co(iz-1)))/rho(iz) enddo return end subroutine diff subroutine buoy(ix,iy,g,kms,kme,kts,kte,th,th0,exch,dz,bu,ustar,tstar,ufrac_int) implicit none integer kms,kme,kts,kte,iz,ix,iy real dtdz1,dtdz2,cdm,dtmdz,g real th(kms:kme),exch(kms:kme),dz(kms:kme),bu(kms:kme) real th0(kms:kme),ustar,tstar,ufrac_int bu(kts)=0. do iz=kts+1,kte-1 dtdz1=2.*(th(iz)-th(iz-1))/(dz(iz-1)+dz(iz)) dtdz2=2.*(th(iz+1)-th(iz))/(dz(iz+1)+dz(iz)) dtmdz=0.5*(dtdz1+dtdz2) cdm=0.5*(exch(iz+1)+exch(iz)) bu(iz)=-cdm*dtmdz*g/th0(iz) enddo bu(kte)=0. return end subroutine buoy subroutine shear(ix,iy,g,kms,kme,kts,kte,u,v,cdua,dz,sh,ustar,tstar,th,ufrac_int) implicit none integer kms,kme,kts,kte,iz,ix,iy real dudz1,dudz2,dvdz1,dvdz2,cdm,dumdz,ustar real tstar,th,al,phim,g real u(kms:kme),v(kms:kme),cdua(kms:kme),dz(kms:kme),sh(kms:kme) real u1,u2,v1,v2,ufrac_int sh(kts)=0. do iz=kts+1,kte-1 u2=(dz(iz+1)*u(iz)+dz(iz)*u(iz+1))/(dz(iz)+dz(iz+1)) u1=(dz(iz)*u(iz-1)+dz(iz-1)*u(iz))/(dz(iz-1)+dz(iz)) v2=(dz(iz+1)*v(iz)+dz(iz)*v(iz+1))/(dz(iz)+dz(iz+1)) v1=(dz(iz)*v(iz-1)+dz(iz-1)*v(iz))/(dz(iz-1)+dz(iz)) cdm=0.5*(cdua(iz)+cdua(iz+1)) dumdz=((u2-u1)/dz(iz))**2.+((v2-v1)/dz(iz))**2. sh(iz)=cdm*dumdz enddo sh(kte)=0. return end subroutine shear subroutine invert(kms,kme,kts,kte,a,c,x) implicit none integer in integer kts,kte,kms,kme real a(kms:kme,3),c(kms:kme),x(kms:kme) do in=kte-1,kts,-1 c(in)=c(in)-a(in,3)*c(in+1)/a(in+1,2) a(in,2)=a(in,2)-a(in,3)*a(in+1,1)/a(in+1,2) enddo do in=kts+1,kte c(in)=c(in)-a(in,1)*c(in-1)/a(in-1,2) enddo do in=kts,kte x(in)=c(in)/a(in,2) enddo return end subroutine invert subroutine tke_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,vl,u,v,th,te,th0,ustar,tstar,exch, & dls,td,sh,bu,b_e,a_e,sf,ufrac_int) implicit none integer kms,kme,kts,kte,iz,ix,iy real g,ustar,tstar,ufrac_int real z(kms:kme),dz(kms:kme),u(kms:kme),v(kms:kme),th(kms:kme),th0(kms:kme),te(kms:kme) real exch(kms:kme),dls(kms:kme),td(kms:kme),sh(kms:kme),bu(kms:kme) real a_e(kms:kme),b_e(kms:kme) real vl(kms:kme),sf(kms:kme) real te1,dl1 call shear(ix,iy,g,kms,kme,kts,kte,u,v,exch,dz,sh,ustar,tstar,th(kts),ufrac_int) call buoy(ix,iy,g,kms,kme,kts,kte,th,th0,exch,dz,bu,ustar,tstar,ufrac_int) do iz=kts,kte te1=max(te(iz),temin) dl1=max(dls(iz),0.1) td(iz)=-ceps_b*sqrt(te1)/dl1 sh(iz)=sh(iz)*sf(iz) bu(iz)=bu(iz)*sf(iz) a_e(iz)=a_e(iz)+td(iz) b_e(iz)=b_e(iz)+sh(iz)+bu(iz) enddo return end subroutine tke_bougeault subroutine pbl_height(kms,kme,kts,kte,dz,z,th,q,pblh) implicit none integer kms,kme,kts,kte,iz real z(kms:kme),dz(kms:kme),th(kms:kme),q(kms:kme) real pblh real thv(kms:kme),zc(kms:kme) real thsfc do iz=kts,kte zc(iz)=z(iz)+dz(iz)/2. enddo do iz=kts,kte thv(iz)=th(iz)*(1.+0.61*q(iz)) enddo pblh=0. thsfc=thv(kts)+0.5 do iz=kts+1,kte if(pblh.eq.0.and.thv(iz).gt.thsfc)then pblh=zc(iz-1)+(thsfc-thv(iz-1))/(max(0.01,thv(iz)-thv(iz-1)))*(zc(iz)-zc(iz-1)) endif enddo return end subroutine pbl_height SUBROUTINE BOULACINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, & & TKE_PBL,EXCH_H,RESTART,ALLOWED_TO_READ, & & IDS,IDE,JDS,JDE,KDS,KDE, & & IMS,IME,JMS,JME,KMS,KME, & & ITS,ITE,JTS,JTE,KTS,KTE ) IMPLICIT NONE LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, & & IMS,IME,JMS,JME,KMS,KME, & & ITS,ITE,JTS,JTE,KTS,KTE REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: EXCH_H, & & RUBLTEN, & & RVBLTEN, & & RTHBLTEN, & & RQVBLTEN, & & TKE_PBL INTEGER :: I,J,K,ITF,JTF,KTF JTF=MIN0(JTE,JDE-1) KTF=MIN0(KTE,KDE-1) ITF=MIN0(ITE,IDE-1) IF(.NOT.RESTART)THEN DO J=JTS,JTF DO K=KTS,KTF DO I=ITS,ITF TKE_PBL(I,K,J)=0.0001 RUBLTEN(I,K,J)=0. RVBLTEN(I,K,J)=0. RTHBLTEN(I,K,J)=0. RQVBLTEN(I,K,J)=0. EXCH_H(I,K,J)=0. ENDDO ENDDO ENDDO ENDIF END SUBROUTINE BOULACINIT END MODULE module_bl_boulac