/wrfv2_fire/phys/module_bl_boulac.F
FORTRAN Legacy | 921 lines | 568 code | 152 blank | 201 comment | 14 complexity | 0590a77da7a2c9d54b43c8e91ec7c0bb MD5 | raw file
Possible License(s): AGPL-1.0
- MODULE module_bl_boulac
- !USE module_model_constants
-
- !------------------------------------------------------------------------
- ! Calculation of the tendency due to momentum, heat
- ! and moisture turbulent fluxes follwing the approach
- ! of Bougeault and Lacarrere, 1989 (MWR, 117, 1872-1890).
- ! The scheme computes a prognostic ecuation for TKE and derives
- ! dissipation and turbulent coefficients using length scales.
- !
- ! Subroutine written by Alberto Martilli, CIEMAT, Spain,
- ! e-mail:alberto_martilli@ciemat.es
- ! August 2006.
- !------------------------------------------------------------------------
- ! IN THIS VERSION TKE IS NOT ADVECTED!!!!
- ! TO BE CHANGED IN THE FUTURE
- !
- ! -----------------------------------------------------------------------
- ! Constant used in the module
- ! ck_b=constant used in the compuation of diffusion coefficients
- ! ceps_b=constant used inthe computation of dissipation
- ! temin= minimum value allowed for TKE
- ! vk=von karman constant
- ! -----------------------------------------------------------------------
- real ck_b,ceps_b,vk,temin ! constant for Bougeault and Lacarrere
- parameter(ceps_b=1/1.4,ck_b=0.4,temin=0.0001,vk=0.4) ! impose minimum values for tke similar to those of MYJ
- ! -----------------------------------------------------------------------
- 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
- !-----------------------------------------------------------------------
- ! Input
- !------------------------------------------------------------------------
- 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 !vertical resolution
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: qv_curr !moisture
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: RHO !air density
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: TH_PHY !potential temperature
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: U_PHY !x-component of wind
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: V_PHY !y-component of wind
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: ustar !friction velocity
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: hfx !sensible heat flux (W/m2) at surface
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: qfx !moisture flux at surface
- real, INTENT(IN ) :: g,cp !gravity and Cp
- REAL, INTENT(IN ):: DT ! Time step
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: FRC_URB2D !fraction cover urban
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: PBLH !PBL height
- !
- ! variable added for urban
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_u_bep ! Implicit component for the momemtum in X-direction
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_v_bep ! Implicit component for the momemtum in Y-direction
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_t_bep ! Implicit component for the Pot. Temp.
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_q_bep ! Implicit component for Moisture
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_e_bep ! Implicit component for the TKE
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_u_bep ! Explicit component for the momemtum in X-direction
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_v_bep ! Explicit component for the momemtum in Y-direction
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_t_bep ! Explicit component for the Pot. Temp.
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_q_bep ! Explicit component for Moisture
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_e_bep ! Explicit component for the TKE
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dlg_bep ! Height above ground (L_ground in formula (24) of the BLM paper).
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::dl_u_bep ! Length scale (lb in formula (22) ofthe BLM paper).
- ! urban surface and volumes
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::sf_bep ! surface of the urban grid cells
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::vl_bep ! volume of the urban grid cells
- LOGICAL, INTENT(IN) :: flag_bep !flag for BEP
-
- !
- !-----------------------------------------------------------------------
- ! Local, carried on from one timestep to the other
- !------------------------------------------------------------------------
- ! real, save, allocatable, dimension (:,:,:)::TKE ! Turbulent kinetic energy
- real, dimension (ims:ime, kms:kme, jms:jme) ::th_0 ! reference state for potential temperature
- !------------------------------------------------------------------------
- ! Output
- !------------------------------------------------------------------------
- real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: exch_h ! exchange coefficient for heat
- real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: exch_m ! exchange coefficient for momentum
- real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: tke ! Turbulence Kinetic Energy
- real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wu ! Turbulent flux of momentum (x)
- real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wv ! Turbulent flux of momentum (y)
- real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wt ! Turbulent flux of temperature
- real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wq ! Turbulent flux of water vapor
- real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: dlk ! Turbulent flux of water vapor
- ! only if idiff not equal 1:
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RUBLTEN !tendency for U_phy
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RVBLTEN !tendency for V_phy
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RTHBLTEN !tendency for TH_phy
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RQVBLTEN !tendency for QV_curr
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RQCBLTEN !tendency for QV_curr
- !--------------------------------------------------------------
- ! Local
- !--------------------------------------------------------------
- ! 1D array used for the input and output of the routine boulac1D
- real z1D(kms:kme) ! vertical coordinates (faces of the grid)
- real dz1D(kms:kme) ! vertical resolution
- real u1D(kms:kme) ! wind speed in the x directions
- real v1D(kms:kme) ! wind speed in the y directions
- real th1D(kms:kme) ! potential temperature
- real q1D(kms:kme) ! potential temperature
- real rho1D(kms:kme) ! air density
- real rhoz1D(kms:kme) ! air density at the faces
- real tke1D(kms:kme) ! air pressure
- real th01D(kms:kme) ! reference potential temperature
- real dlk1D(kms:kme) ! dlk
- real dls1D(kms:kme) ! dls
- real exch1D(kms:kme) ! exch
- real sf1D(kms:kme) ! surface of the grid cells
- real vl1D(kms:kme) ! volume of the grid cells
- real a_u1D(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
- real a_v1D(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
- real a_t1D(kms:kme) ! Implicit component of the heat sources or sinks
- real a_q1D(kms:kme) ! Implicit component of the moisture sources or sinks
- real a_e1D(kms:kme) ! Implicit component of the TKE sources or sinks
- real b_u1D(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
- real b_v1D(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
- real b_t1D(kms:kme) ! Explicit component of the heat sources or sinks
- real b_q1D(kms:kme) ! Explicit component of the moisture sources or sinks
- real b_e1D(kms:kme) ! Explicit component of the TKE sources or sinks
- real dlg1D(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
- real dl_u1D(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper)
- real sh1D(kms:kme) ! shear
- real bu1D(kms:kme) ! buoyancy
- real wu1D(kms:kme) ! turbulent flux of momentum (x component)
- real wv1D(kms:kme) ! turbulent flux of momentum (y component)
- real wt1D(kms:kme) ! turbulent flux of temperature
- real wq1D(kms:kme) ! turbulent flux of water vapor
- ! local added only for diagnostic output
- real a_e(ims:ime,kms:kme,jms:jme) ! implicit term in TKE
- real b_e(ims:ime,kms:kme,jms:jme) ! explicit term in TKE
- real bu(ims:ime,kms:kme,jms:jme) ! buoyancy term in TKE
- real sh(ims:ime,kms:kme,jms:jme) ! shear term in TKE
- real wrk(ims:ime) ! working array
- integer ix,iy,iz,id,iz_u,iw_u,ig,ir_u,ix1,iy1
- real ufrac_int ! urban fraction
- real vect,time_tke,hour,zzz
- real ustarf
- real summ1,summ2,summ3
- save time_tke,hour
- !
- ! store reference state for potential temperature
- ! and initialize tke
- !
- !here I fix the value of the reference state equal to the value of the potnetial temperature
- ! the only use of this variable in the code is to compute the paramter BETA = g/th0
- ! I fix it to 300K.
-
- do ix=its,ite
- do iy=jts,jte
- do iz=kts,kte
- ! th_0(ix,iz,iy)=th_phy(ix,iz,iy)
- 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. ! vertical coordinates (faces of the grid)
- dz1D=0. ! vertical resolution
- u1D =0. ! wind speed in the x directions
- v1D =0. ! wind speed in the y directions
- th1D=0. ! potential temperature
- q1D=0. ! potential temperature
- rho1D=0. ! air density
- rhoz1D=0. ! air density at the faces
- tke1D =0. ! air pressure
- th01D =0. ! reference potential temperature
- dlk1D =0. ! dlk
- dls1D =0. ! dls
- exch1D=0. ! exch
- sf1D =1. ! surface of the grid cells
- vl1D =1. ! volume of the grid cells
- a_u1D =0. ! Implicit component of the momentum sources or sinks in the X-direction
- a_v1D =0. ! Implicit component of the momentum sources or sinks in the Y-direction
- a_t1D =0. ! Implicit component of the heat sources or sinks
- a_q1D =0. ! Implicit component of the moisture sources or sinks
- a_e1D =0. ! Implicit component of the TKE sources or sinks
- b_u1D =0. ! Explicit component of the momentum sources or sinks in the X-direction
- b_v1D =0. ! Explicit component of the momentum sources or sinks in the Y-direction
- b_t1D =0. ! Explicit component of the heat sources or sinks
- b_q1D =0. ! Explicit component of the moisture sources or sinks
- b_e1D =0. ! Explicit component of the TKE sources or sinks
- dlg1D =0. ! Height above ground (L_ground in formula (24) of the BLM paper).
- dl_u1D=0. ! Length scale (lb in formula (22) ofthe BLM paper)
- sh1D =0. ! shear
- bu1D =0. ! buoyancy
- wu1D =0. ! turbulent flux of momentum (x component)
- wv1D =0. ! turbulent flux of momentum (y component)
- wt1D =0. ! turbulent flux of temperature
- wq1D =0. ! turbulent flux of water vapor
- ! local added only for diagnostic output
-
- ! loop over the columns.
- ! put variables in 1D temporary arrays
- !
- 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 the routine that will solve the turbulence in 1D for tke
- 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))
- ! store turbulent exchange coefficients, TKE, and other variables
-
- 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
- ! estimate the tendencies
- 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
- ! solve diffusion equation for momentum x component
- call diff(kms,kme,kts,kte,1,1,dt,u1D,rho1D,rhoz1D,exch1D,a_u1D,b_u1D,sf1D,vl1D,dz1D,wu1D)
-
- ! solve diffusion equation for momentum y component
- call diff(kms,kme,kts,kte,1,1,dt,v1D,rho1D,rhoz1D,exch1D,a_v1D,b_v1D,sf1D,vl1D,dz1D,wv1D)
- ! solve diffusion equation for potential temperature
- call diff(kms,kme,kts,kte,1,1,dt,th1D,rho1D,rhoz1D,exch1D,a_t1D,b_t1D,sf1D,vl1D,dz1D,wt1D)
- ! solve diffusion equation for water vapor mixing ratio
- call diff(kms,kme,kts,kte,1,1,dt,q1D,rho1D,rhoz1D,exch1D,a_q1D,b_q1D,sf1D,vl1D,dz1D,wq1D)
- ! compute the tendencies
-
- 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 ! iy
- enddo ! ix
-
- time_tke=time_tke+dt
-
- return
- end subroutine boulac
-
- ! ===6=8===============================================================72
- ! ===6=8===============================================================72
- 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)
- ! ----------------------------------------------------------------------
- ! 1D resolution of TKE following Bougeault and Lacarrere
- ! ----------------------------------------------------------------------
- implicit none
- integer iz,ix,iy
- ! ----------------------------------------------------------------------
- ! INPUT:
- ! ----------------------------------------------------------------------
- integer kms,kme,kts,kte
- real z(kms:kme) ! Altitude above the ground of the cell interfaces.
- real dz(kms:kme) ! vertical resolution
- real u(kms:kme) ! Wind speed in the x direction
- real v(kms:kme) ! Wind speed in the y direction
- real th(kms:kme) ! Potential temperature
- real rho(kms:kme) ! Air density
- real g ! gravity
- real cp !
- real te(kms:kme) ! turbulent kinetic energy
- real qa(kms:kme) ! air humidity
- real th0(kms:kme) ! Reference potential temperature
- real dt ! Time step
- real ustar ! ustar
- real hfx ! sensbile heat flux
- real qfx ! kinematic latent heat flux
- real sf(kms:kme) ! surface of the urban grid cells
- real vl(kms:kme) ! volume of the urban grid cells
- real a_e(kms:kme) ! Implicit component of the TKE sources or sinks
- real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
- real dlg(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
- real dl_u(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper)
- real ufrac_int ! urban fraction
- ! local variables not needed in principle, but that can be used to estimate the budget and turbulent fluxes
-
- real we(kms:kme),dwe(kms:kme)
-
- ! local variables
- real sh(kms:kme) ! shear term in TKE eqn.
- real bu(kms:kme) ! buoyancy term in TKE eqn.
- real td(kms:kme) ! dissipation term in TKE eqn.
- real exch(kms:kme) ! turbulent diffusion coefficients (defined at the faces)
- real dls(kms:kme) ! dissipation length scale
- real dlk(kms:kme) ! length scale used to estimate exch
- real dlu(kms:kme) ! l_up
- real dld(kms:kme) ! l_down
- real rhoz(kms:kme) !air density at the faces of the cell
- real tstar ! derived from hfx and ustar
- real beta
- real summ1,summ2,summ3,summ4
- ! interpolate air density at the faces
-
- ! estimation of tstar
- tstar=-hfx/rho(1)/cp/ustar
-
- ! first compute values of dlu and dld (length scales up and down).
-
- call dissip_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,te,dlu,dld,th,th0)
-
- !then average them to obtain dls and dlk (length scales for dissipation and eddy coefficients)
- call length_bougeault(ix,iy,kms,kme,kts,kte,dld,dlu,dlg,dl_u,dls,dlk)
-
- ! compute the turbulent diffusion coefficients exch
- call cdtur_bougeault(ix,iy,kms,kme,kts,kte,te,z,dz,exch,dlk)
- ! compute source and sink terms in the TKE equation (shear, buoyancy and dissipation)
-
- 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)
-
- ! solve for tke
-
- call diff(kms,kme,kts,kte,1,1,dt,te,rho,rhoz,exch,a_e,b_e,sf,vl,dz,we)
-
- ! avoid negative values for tke
-
- do iz=kts,kte
- if(te(iz).lt.temin) te(iz)=temin
- enddo
-
- return
- end subroutine boulac1d
- !
- ! ===6=8===============================================================72
- ! ===6=8===============================================================72
- subroutine dissip_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,te,dlu,dld,th,th0)
- ! compute the length scales up and down
- 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) !Buoyancy coefficient
- 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
- !
- ! ===6=8===============================================================72
- ! ===6=8===============================================================72
- subroutine length_bougeault(ix,iy,kms,kme,kts,kte,dld,dlu,dlg,dl_u,dls,dlk)
- ! compute the length scales for dissipation and turbulent coefficients
- 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
- !
- ! ===6=8===============================================================72
- ! ===6=8===============================================================72
- subroutine cdtur_bougeault(ix,iy,kms,kme,kts,kte,te,z,dz,exch,dlk)
- ! compute turbulent coefficients
- 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=2,nz-1
- 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.0001)
- exch(iz)=max(exch(iz),0.1)
- enddo
- exch(kte+1)=0.1
- return
- end subroutine cdtur_bougeault
- ! ===6=8===============================================================72
- ! ===6=8===============================================================72
- subroutine diff(kms,kme,kts,kte,iz1,izf,dt,co,rho,rhoz,cd,aa,bb,sf,vl,dz,fc)
- !------------------------------------------------------------------------
- ! Calculation of the diffusion in 1D
- !------------------------------------------------------------------------
- ! - Input:
- ! nz : number of points
- ! iz1 : first calculated point
- ! co : concentration of the variable of interest
- ! dz : vertical levels
- ! cd : diffusion coefficients
- ! dtext : external time step
- ! rho : density of the air at the center
- ! rhoz : density of the air at the face
- ! itest : if itest eq 1 then update co, else store in a flux array
- ! - Output:
- ! co :concentration of the variable of interest
- ! - Internal:
- ! cddz : constant terms in the equations
- ! dt : diffusion time step
- ! nt : number of the diffusion time steps
- ! cstab : ratio of the stability condition for the time step
- !---------------------------------------------------------------------
- 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)
-
- ! Compute cddz=2*cd/dz
-
- 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
-
- ! do iz=1,iz1
- ! df(iz)=0.
- ! enddo
- !
- ! do iz=iz1+1,nz-izf
- ! dzv=vl(iz)*dz(iz)
- ! df(iz)=+(co(iz-1)*cddz(iz)-co(iz)*(cddz(iz)+cddz(iz+1))+co(iz+1)*cddz(iz+1))/dzv/rho(iz)
- ! enddo
- !
- ! do iz=nz-izf,nz
- ! df(iz)=0.
- ! enddo
-
- return
- end subroutine diff
- ! ===6=8===============================================================72
- ! ===6=8===============================================================72
- subroutine buoy(ix,iy,g,kms,kme,kts,kte,th,th0,exch,dz,bu,ustar,tstar,ufrac_int)
- ! compute buoyancy term
- 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(1)=-ustar*tstar*g/th0(1)*(1.-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
- ! ===6=8===============================================================72
- ! ===6=8===============================================================72
- subroutine shear(ix,iy,g,kms,kme,kts,kte,u,v,cdua,dz,sh,ustar,tstar,th,ufrac_int)
- ! compute shear term
- 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
- ! al=vk*g*tstar/(th*(ustar**2.))
- ! if(al.ge.0.)phim=1.+4.7*dz(1)/2.*al
- ! if(al.lt.0.)phim=(1.-15*dz(1)/2.*al)**(-0.25)
- !
- ! sh(1)=(ustar**3.)/vk/(dz(1)/2.)*(1.-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
- ! ===6=8===============================================================72
- ! ===6=8===============================================================72
- subroutine invert(kms,kme,kts,kte,a,c,x)
-
- !ccccccccccccccccccccccccccccccc
- ! Aim: Inversion and resolution of a tridiagonal matrix
- ! A X = C
- ! Input:
- ! a(*,1) lower diagonal (Ai,i-1)
- ! a(*,2) principal diagonal (Ai,i)
- ! a(*,3) upper diagonal (Ai,i+1)
- ! c
- ! Output
- ! x results
- !ccccccccccccccccccccccccccccccc
- 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
-
- ! ===6=8===============================================================72
- ! ===6=8===============================================================72
-
- 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)
- ! in this routine the shear, buoyancy and part of the dissipation terms
- ! of the TKE equation are computed
- 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)
- ! this routine computes the PBL height
- ! with an approach similar to MYNN
- implicit none
- integer kms,kme,kts,kte,iz
- real z(kms:kme),dz(kms:kme),th(kms:kme),q(kms:kme)
- real pblh
- !Local
- real thv(kms:kme),zc(kms:kme)
- real thsfc
- ! compute the height of the center of the grid cells
- do iz=kts,kte
- zc(iz)=z(iz)+dz(iz)/2.
- enddo
- ! compute the virtual potential temperature
- do iz=kts,kte
- thv(iz)=th(iz)*(1.+0.61*q(iz))
- enddo
- ! now compute the PBL height
- 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))
- ! pblh=z(iz-1)+(thsfc-thv(iz-1))/(max(0.01,thv(iz)-thv(iz-1)))*(z(iz)-z(iz-1))
- endif
- enddo
- return
- end subroutine pbl_height
- ! ===6=8===============================================================72
- ! ===6=8===============================================================72
- 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