Module module_chem_plumerise_scalar
! USE module_plumerise1
use module_model_constants
use module_zero_plumegen_coms
real,parameter :: rgas=r_d
real,parameter :: cpor=1./rcp
real,parameter :: p00=p1000mb
! real, external:: esat_pr
CONTAINS
subroutine plumerise(m1,m2,m3,ia,iz,ja,jz,firesize,mean_fct &
,nspecies,eburn_in,eburn_out &
,up,vp,wp,theta,pp,dn0,rv,zt_rams,zm_rams)
! use module_zero_plumegen_coms, only : ucon,vcon,wcon,thtcon,rvcon,picon,tmpcon&
! ,dncon,prcon,zcon,zzcon,scon
implicit none
integer :: ng,m1,m2,m3,ia,iz,ja,jz,ibcon,mynum,i,j,k,iveg_ag,&
imm,k1,k2,ixx,ispc,nspecies
integer :: ncall = 0
integer :: kmt
real,dimension(m1,nspecies), intent(out) :: eburn_out
real,dimension(nspecies), intent(in) :: eburn_in
real, dimension(m1,m2,m3) :: up, vp, wp,theta,pp,dn0,rv
real, dimension(m1) :: zt_rams,zm_rams
real :: burnt_area,STD_burnt_area,dz_flam,rhodzi,dzi
real, dimension(2) :: ztopmax
real :: q_smold_kgm2
! From plumerise1.F routine
integer, parameter :: nveg_agreg = 4
integer, parameter :: tropical_forest = 1
integer, parameter :: boreal_forest = 2
integer, parameter :: savannah = 3
integer, parameter :: grassland = 4
real, dimension(nveg_agreg) :: firesize,mean_fct
INTEGER, PARAMETER :: wind_eff = 1
!Fator de conversao de unidades
!!fcu=1. !=> kg [gas/part] /kg [ar]
!!fcu =1.e+12 !=> ng [gas/part] /kg [ar]
!!real,parameter :: fcu =1.e+6 !=> mg [gas/part] /kg [ar]
!----------------------------------------------------------------------
! indexacao para o array "plume(k,i,j)"
! k
! 1 => area media (m^2) dos focos em biomas floresta dentro do gribox i,j
! 2 => area media (m^2) dos focos em biomas savana dentro do gribox i,j
! 3 => area media (m^2) dos focos em biomas pastagem dentro do gribox i,j
! 4 => desvio padrao da area media (m^2) dos focos : floresta
! 5 => desvio padrao da area media (m^2) dos focos : savana
! 6 => desvio padrao da area media (m^2) dos focos : pastagem
! 7 a 9 => sem uso
!10(=k_CO_smold) => parte da emissao total de CO correspondente a fase smoldering
!11, 12 e 13 => este array guarda a relacao entre
! qCO( flaming, floresta) e a quantidade total emitida
! na fase smoldering, isto e;
! qCO( flaming, floresta) = plume(11,i,j)*plume(10,i,j)
! qCO( flaming, savana ) = plume(12,i,j)*plume(10,i,j)
! qCO( flaming, pastagem) = plume(13,i,j)*plume(10,i,j)
!20(=k_PM25_smold),21,22 e 23 o mesmo para PM25
!
!24-n1 => sem uso
!----------------------------------------------------------------------
! print *,' Plumerise_scalar 1',ncall
if (ncall == 0) then
ncall = 1
call zero_plumegen_coms
endif
! print *,' Plumerise_scalar 1',ncall
! print *,' Plumerise_scalar 2',m1
j=1
i=1
! do j = ja,jz ! loop em j
! do i = ia,iz ! loop em i
!- if the max value of flaming is close to zero => there is not emission with
!- plume rise => cycle
do k = 1,m1
ucon (k)=up(k,i,j) ! u wind
vcon (k)=vp(k,i,j) ! v wind
!wcon (k)=wp(k,i,j) ! w wind
thtcon(k)=theta(k,i,j) ! pot temperature
picon (k)=pp(k,i,j) ! exner function
!tmpcon(k)=thtcon(k)*picon(k)/cp ! temperature (K)
!dncon (k)=dn0(k,i,j) ! dry air density (basic state)
!prcon (k)=(picon(k)/cp)**cpor*p00 ! pressure (Pa)
rvcon (k)=rv(k,i,j) ! water vapor mixing ratio
zcon (k)=zt_rams(k) ! termod-point height
zzcon (k)=zm_rams(k) ! W-point height
! print*,'PL:',k,zcon(k),picon (k),thtcon(k),1000.*rvcon (k)
enddo
do ispc=1,nspecies
eburn_out(1,ispc) = eburn_in(ispc)
enddo
!- get envinronmental state (temp, water vapor mix ratio, ...)
call get_env_condition(1,m1,kmt,wind_eff)
!- loop nos 4 biomas agregados com possivel queimada
do iveg_ag=1,nveg_agreg
! print *,'iveg_ag = ',iveg_ag,mean_fct(iveg_ag)
!- verifica a existencia de emissao flaming para um bioma especifico
!orig: if( plume( k_CO_smold + iveg_ag ,i,j) < 1.e-6 ) cycle
if(mean_fct(iveg_ag) < 1.e-6 ) cycle
! burnt area and standard deviation
burnt_area = firesize(iveg_ag)
!not em use
!STD_burnt_area= plume(3+iveg_ag,i,j)
STD_burnt_area= 0.
!- loop nos valores minimo e maximo da taxa de calor
do imm=1,2
!--------------------
!ixx=iveg_ag*10 + imm
! print*,'i j veg=',i, j, iveg_ag,imm
!--------------------
!- get fire properties (burned area, plume radius, heating rates ...)
call get_fire_properties(imm,iveg_ag,burnt_area,STD_burnt_area)
!------ generates the plume rise ------
!-- only one value for eflux of GRASSLAND
! if(iveg_ag == GRASSLAND .and. imm == 2) then
if(iveg_ag == 4 .and. imm == 2) then
ztopmax(2)=ztopmax(1)
ztopmax(1)=zzcon(1)
! print *,'cycle',ztopmax(1),ztopmax(2)
cycle
endif
call makeplume (kmt,ztopmax(imm),ixx,imm)
enddo ! enddo do loop em imm
!- define o dominio vertical onde a emissao flaming ira ser colocada
call set_flam_vert(ztopmax,k1,k2,nkp,zzcon,W_VMD,VMD)
!- espessura da camada vertical
dz_flam=zzcon(k2)-zzcon(k1-1)
!- distribui a emissao flaming entre os niveis k1 e k2
! print *,'distribui, k1,k2,dz_flam',k1,k2,dz_flam
do k=k1,k2
!use this in case the emission src is already in mixing ratio
!rhodzi= 1./(dn0(k,i,j) * dz_flam)
!use this in case the emission src is tracer density
dzi= 1./( dz_flam)
do ispc = 1, nspecies
!- get back the smoldering emission in kg/m2 (actually in 1e-9 kg/m2)
!use this in case the emission src is already in mixing ratio
!q_smold_kgm2 = (1/dzt(2) * dn0(2,i,j) )* &
! chem1_src_g(bburn,ispc,ng)%sc_src(2,i,j)
!use this in case the emission src is tracer density
! q_smold_kgm2 = ((zt_rams(2)-zt_rams(1)) )* &
! eburn_in(ispc)
q_smold_kgm2 = eburn_in(ispc)
! units = already in ppbm, don't need "fcu" factor
eburn_out(k,ispc) = eburn_out(k,ispc) +&
mean_fct(iveg_ag) *&
q_smold_kgm2 * &
dzi !use this in case the emission src is tracer density
!rhodzi !use this in case the emission src is already in mixing ratio
! if(ispc.eq. 1) print *,' Plume: ',k,eburn_out(k,ispc),mean_fct(iveg_ag),q_smold_kgm2,dzi
!srcCO(k,i,j)= srcCO(k,i,j) + plume(k_CO_smold+iveg_ag,i,j)*&
! plume(k_CO_smold ,i,j)*&
! rhodzi*fcu
enddo
enddo
enddo ! enddo do loop em iveg_ag
!-----
!stop 333
!endif
!-----
! enddo ! loop em i
! enddo ! loop em j
!stop 4544
end subroutine plumerise
!-------------------------------------------------------------------------
subroutine get_env_condition(k1,k2,kmt,wind_eff)
!se module_zero_plumegen_coms
!use rconstants
implicit none
integer :: k1,k2,k,kcon,klcl,kmt,nk,nkmid,i
real :: znz,themax,tlll,plll,rlll,zlll,dzdd,dzlll,tlcl,plcl,dzlcl,dummy
integer :: n_setgrid = 0
integer :: wind_eff
if( n_setgrid == 0) then
n_setgrid = 1
call set_grid ! define vertical grid of plume model
! zt(k) = thermo and water levels
! zm(k) = dynamical levels
endif
znz=zcon(k2)
do k=nkp,1,-1
if(zt(k).lt.znz)go to 13
enddo
stop ' envir stop 12'
13 continue
!-srf-mb
kmt=min(k,nkp-1)
nk=k2-k1+1
!call htint(nk, wcon,zzcon,kmt,wpe,zt)
call htint(nk, ucon,zcon,kmt,upe,zt)
call htint(nk, vcon,zcon,kmt,vpe,zt)
call htint(nk,thtcon,zcon,kmt,the ,zt)
call htint(nk, rvcon,zcon,kmt,qvenv,zt)
do k=1,kmt
qvenv(k)=max(qvenv(k),1e-8)
enddo
pke(1)=picon(1)
do k=1,kmt
thve(k)=the(k)*(1.+.61*qvenv(k)) ! virtual pot temperature
enddo
do k=2,kmt
pke(k)=pke(k-1)-g*2.*(zt(k)-zt(k-1)) & ! exner function
/(thve(k)+thve(k-1))
enddo
do k=1,kmt
te(k) = the(k)*pke(k)/cp ! temperature (K)
pe(k) = (pke(k)/cp)**cpor*p00 ! pressure (Pa)
dne(k)= pe(k)/(rgas*te(k)*(1.+.61*qvenv(k))) ! dry air density (kg/m3)
! print*,'ENV=',qvenv(k)*1000., te(k)-273.15,zt(k)
!-srf-mb
vel_e(k) = sqrt(upe(k)**2+vpe(k)**2) !-env wind (m/s)
!print*,'k,vel_e(k),te(k)=',vel_e(k),te(k)
enddo
!-ewe - env wind effect
if(wind_eff < 1) vel_e(1:kmt) = 0.
!-use este para gerar o RAMS.out
! ------- print environment state
!print*,'k,zt(k),pe(k),te(k)-273.15,qvenv(k)*1000'
!do k=1,kmt
! write(*,100) k,zt(k),pe(k),te(k)-273.15,qvenv(k)*1000.
! 100 format(1x,I5,4f20.12)
!enddo
!stop 333
!--------- nao eh necessario este calculo
!do k=1,kmt
! call thetae(pe(k),te(k),qvenv(k),thee(k))
!enddo
!--------- converte press de Pa para kPa para uso modelo de plumerise
do k=1,kmt
pe(k) = pe(k)*1.e-3
enddo
return
end subroutine get_env_condition
!-------------------------------------------------------------------------
subroutine set_grid()
!use module_zero_plumegen_coms
implicit none
integer :: k,mzp
dz=100. ! set constant grid spacing of plume grid model(meters)
mzp=nkp
zt(1) = zsurf
zm(1) = zsurf
zt(2) = zt(1) + 0.5*dz
zm(2) = zm(1) + dz
do k=3,mzp
zt(k) = zt(k-1) + dz ! thermo and water levels
zm(k) = zm(k-1) + dz ! dynamical levels
enddo
!print*,zsurf
!Print*,zt(:)
do k = 1,mzp-1
dzm(k) = 1. / (zt(k+1) - zt(k))
enddo
dzm(mzp)=dzm(mzp-1)
do k = 2,mzp
dzt(k) = 1. / (zm(k) - zm(k-1))
enddo
dzt(1) = dzt(2) * dzt(2) / dzt(3)
! dzm(1) = 0.5/dz
! dzm(2:mzp) = 1./dz
return
end subroutine set_grid
!-------------------------------------------------------------------------
SUBROUTINE set_flam_vert(ztopmax,k1,k2,nkp,zzcon,W_VMD,VMD)
REAL , INTENT(IN) :: ztopmax(2)
INTEGER , INTENT(OUT) :: k1
INTEGER , INTENT(OUT) :: k2
! plumegen_coms
INTEGER , INTENT(IN) :: nkp
REAL , INTENT(IN) :: zzcon(nkp)
INTEGER imm,k
INTEGER, DIMENSION(2) :: k_lim
!- version 2
REAL , INTENT(IN) :: W_VMD(nkp,2)
REAL , INTENT(OUT) :: VMD(nkp,2)
real w_thresold,xxx
integer k_initial,k_final,ko,kk4,kl
!- version 1
DO imm=1,2
! checar
! do k=1,m1-1
DO k=1,nkp-1
IF(zzcon(k) > ztopmax(imm) ) EXIT
ENDDO
k_lim(imm) = k
ENDDO
k1=MAX(3,k_lim(1))
k2=MAX(3,k_lim(2))
IF(k2 < k1) THEN
!print*,'1: ztopmax k=',ztopmax(1), k1
!print*,'2: ztopmax k=',ztopmax(2), k2
k2=k1
!stop 1234
ENDIF
!- version 2
!- vertical mass distribution
!-
w_thresold = 1.
DO imm=1,2
VMD(1:nkp,imm)= 0.
xxx=0.
k_initial= 0
k_final = 0
!- define range of the upper detrainemnt layer
do ko=nkp-10,2,-1
if(w_vmd(ko,imm) < w_thresold) cycle
if(k_final==0) k_final=ko
if(w_vmd(ko,imm)-1. > w_vmd(ko-1,imm)) then
k_initial=ko
exit
endif
enddo
!- if there is a non zero depth layer, make the mass vertical distribution
if(k_final > 0 .and. k_initial > 0) then
k_initial=int((k_final+k_initial)*0.5)
!- parabolic vertical distribution between k_initial and k_final
kk4 = k_final-k_initial+2
do ko=1,kk4-1
kl=ko+k_initial-1
VMD(kl,imm) = 6.* float(ko)/float(kk4)**2 * (1. - float(ko)/float(kk4))
enddo
if(sum(VMD(1:NKP,imm)) .ne. 1.) then
xxx= ( 1.- sum(VMD(1:NKP,imm)) )/float(k_final-k_initial+1)
do ko=k_initial,k_final
VMD(ko,imm) = VMD(ko,imm)+ xxx !- values between 0 and 1.
enddo
! print*,'new mass=',sum(mass)*100.,xxx
!pause
endif
endif !k_final > 0 .and. k_initial >
ENDDO
END SUBROUTINE set_flam_vert
!-------------------------------------------------------------------------
subroutine get_fire_properties(imm,iveg_ag,burnt_area,STD_burnt_area)
!use module_zero_plumegen_coms
implicit none
integer :: moist, i, icount,imm,iveg_ag
real:: bfract, effload, heat, hinc ,burnt_area,STD_burnt_area,heat_fluxW
real, dimension(2,4) :: heat_flux
INTEGER, parameter :: use_last = 0
data heat_flux/ &
!---------------------------------------------------------------------
! heat flux !IGBP Land Cover !
! min ! max !Legend and ! reference
! kW/m^2 !description !
!--------------------------------------------------------------------
30.0, 80.0, &! Tropical Forest ! igbp 2 & 4
30.0, 80.0, &! Boreal forest ! igbp 1 & 3
4.4, 23.0, &! cerrado/woody savanna | igbp 5 thru 9
3.3, 3.3 /! Grassland/cropland ! igbp 10 thru 17
!--------------------------------------------------------------------
!-- fire at the surface
!
!area = 20.e+4 ! area of burn, m^2
area = burnt_area! area of burn, m^2
!fluxo de calor para o bioma
heat_fluxW = heat_flux(imm,iveg_ag) * 1000. ! converte para W/m^2
mdur = 53 ! duration of burn, minutes
bload = 10. ! total loading, kg/m**2
moist = 10 ! fuel moisture, %. average fuel moisture,percent dry
maxtime =mdur+2 ! model time, min
!heat = 21.e6 !- joules per kg of fuel consumed
!heat = 15.5e6 !joules/kg - cerrado
heat = 19.3e6 !joules/kg - floresta em alta floresta (mt)
!alpha = 0.1 !- entrainment constant
alpha = 0.05 !- entrainment constant
!-------------------- printout ----------------------------------------
!!WRITE ( * , * ) ' SURFACE =', ZSURF, 'M', ' LCL =', ZBASE, 'M'
!
!PRINT*,'======================================================='
!print * , ' FIRE BOUNDARY CONDITION :'
!print * , ' DURATION OF BURN, MINUTES =',MDUR
!print * , ' AREA OF BURN, HA =',AREA*1.e-4
!print * , ' HEAT FLUX, kW/m^2 =',heat_fluxW*1.e-3
!print * , ' TOTAL LOADING, KG/M**2 =',BLOAD
!print * , ' FUEL MOISTURE, % =',MOIST !average fuel moisture,percent dry
!print * , ' MODEL TIME, MIN. =',MAXTIME
!
!
!
! ******************** fix up inputs *********************************
!
!IF (MOD (MAXTIME, 2) .NE.0) MAXTIME = MAXTIME+1 !make maxtime even
MAXTIME = MAXTIME * 60 ! and put in seconds
!
RSURF = SQRT (AREA / 3.14159) !- entrainment surface radius (m)
FMOIST = MOIST / 100. !- fuel moisture fraction
!
!
! calculate the energy flux and water content at lboundary.
! fills heating() on a minute basis. could ask for a file at this po
! in the program. whatever is input has to be adjusted to a one
! minute timescale.
!
DO I = 1, ntime !- make sure of energy release
HEATING (I) = 0.0001 !- avoid possible divide by 0
enddo
!
TDUR = MDUR * 60. !- number of seconds in the burn
bfract = 1. !- combustion factor
EFFLOAD = BLOAD * BFRACT !- patchy burning
! spread the burning evenly over the interval
! except for the first few minutes for stability
ICOUNT = 1
!
if(MDUR > NTIME) STOP 'Increase time duration (ntime) in min - see file "plumerise_mod.f90"'
DO WHILE (ICOUNT.LE.MDUR)
! HEATING (ICOUNT) = HEAT * EFFLOAD / TDUR ! W/m**2
! HEATING (ICOUNT) = 80000. * 0.55 ! W/m**2
HEATING (ICOUNT) = heat_fluxW * 0.55 ! W/m**2 (0.55 converte para energia convectiva)
ICOUNT = ICOUNT + 1
ENDDO
! ramp for 5 minutes
IF(use_last /= 1) THEN
HINC = HEATING (1) / 4.
HEATING (1) = 0.1
HEATING (2) = HINC
HEATING (3) = 2. * HINC
HEATING (4) = 3. * HINC
ELSE
IF(imm==1) THEN
HINC = HEATING (1) / 4.
HEATING (1) = 0.1
HEATING (2) = HINC
HEATING (3) = 2. * HINC
HEATING (4) = 3. * HINC
ELSE
HINC = (HEATING (1) - heat_flux(imm-1,iveg_ag) * 1000. *0.55)/ 4.
HEATING (1) = heat_flux(imm-1,iveg_ag) * 1000. *0.55 + 0.1
HEATING (2) = HEATING (1)+ HINC
HEATING (3) = HEATING (2)+ HINC
HEATING (4) = HEATING (3)+ HINC
ENDIF
ENDIF
return
end subroutine get_fire_properties
!-------------------------------------------------------------------------------
!
SUBROUTINE MAKEPLUME ( kmt,ztopmax,ixx,imm)
!
! *********************************************************************
!
! EQUATION SOURCE--Kessler Met.Monograph No. 32 V.10 (K)
! Alan Weinstein, JAS V.27 pp 246-255. (W),
! Ogura and Takahashi, Monthly Weather Review V.99,pp895-911 (OT)
! Roger Pielke,Mesoscale Meteorological Modeling,Academic Press,1984
! Originally developed by: Don Latham (USFS)
!
!
! ************************ VARIABLE ID ********************************
!
! DT=COMPUTING TIME INCREMENT (SEC)
! DZ=VERTICAL INCREMENT (M)
! LBASE=LEVEL ,CLOUD BASE
!
! CONSTANTS:
! G = GRAVITATIONAL ACCELERATION 9.80796 (M/SEC/SEC).
! R = DRY AIR GAS CONSTANT (287.04E6 JOULE/KG/DEG K)
! CP = SPECIFIC HT. (1004 JOULE/KG/DEG K)
! HEATCOND = HEAT OF CONDENSATION (2.5E6 JOULE/KG)
! HEATFUS = HEAT OF FUSION (3.336E5 JOULE/KG)
! HEATSUBL = HEAT OF SUBLIMATION (2.83396E6 JOULE/KG)
! EPS = RATIO OF MOL.WT. OF WATER VAPOR TO THAT OF DRY AIR (0.622)
! DES = DIFFERENCE BETWEEN VAPOR PRESSURE OVER WATER AND ICE (MB)
! TFREEZE = FREEZING TEMPERATURE (K)
!
!
! PARCEL VALUES:
! T = TEMPERATURE (K)
! TXS = TEMPERATURE EXCESS (K)
! QH = HYDROMETEOR WATER CONTENT (G/G DRY AIR)
! QHI = HYDROMETEOR ICE CONTENT (G/G DRY AIR)
! QC = WATER CONTENT (G/G DRY AIR)
! QVAP = WATER VAPOR MIXING RATIO (G/G DRY AIR)
! QSAT = SATURATION MIXING RATIO (G/G DRY AIR)
! RHO = DRY AIR DENSITY (G/M**3) MASSES = RHO*Q'S IN G/M**3
! ES = SATURATION VAPOR PRESSURE (kPa)
!
! ENVIRONMENT VALUES:
! TE = TEMPERATURE (K)
! PE = PRESSURE (kPa)
! QVENV = WATER VAPOR (G/G)
! RHE = RELATIVE HUMIDITY FRACTION (e/esat)
! DNE = dry air density (kg/m^3)
!
! HEAT VALUES:
! HEATING = HEAT OUTPUT OF FIRE (WATTS/M**2)
! MDUR = DURATION OF BURN, MINUTES
!
! W = VERTICAL VELOCITY (M/S)
! RADIUS=ENTRAINMENT RADIUS (FCN OF Z)
! RSURF = ENTRAINMENT RADIUS AT GROUND (SIMPLE PLUME, TURNER)
! ALPHA = ENTRAINMENT CONSTANT
! MAXTIME = TERMINATION TIME (MIN)
!
!
!**********************************************************************
!**********************************************************************
!use module_zero_plumegen_coms
implicit none
!logical :: endspace
character (len=10) :: varn
integer :: izprint, iconv, itime, k, kk, kkmax, deltak,ilastprint,kmt &
,ixx,nrectotal,i_micro,n_sub_step
real :: vc, g, r, cp, eps, &
tmelt, heatsubl, heatfus, heatcond, tfreeze, &
ztopmax, wmax, rmaxtime, es, esat, heat,dt_save !ESAT_PR,
character (len=2) :: cixx
! Set threshold to be the same as dz=100., the constant grid spacing of plume grid model(meters) found in set_grid()
REAL :: DELZ_THRESOLD = 100.
INTEGER :: imm
! real, external:: esat_pr!
!
! ******************* SOME CONSTANTS **********************************
!
! XNO=10.0E06 median volume diameter raindrop (K table 4)
! VC = 38.3/(XNO**.125) mean volume fallspeed eqn. (K)
!
parameter (vc = 5.107387)
parameter (g = 9.80796, r = 287.04, cp = 1004., eps = 0.622, tmelt = 273.3)
parameter (heatsubl = 2.834e6, heatfus = 3.34e5, heatcond = 2.501e6)
parameter (tfreeze = 269.3)
!
tstpf = 2.0 !- timestep factor
viscosity = 500.!- viscosity constant (original value: 0.001)
nrectotal=150
!
!*************** PROBLEM SETUP AND INITIAL CONDITIONS *****************
mintime = 1
ztopmax = 0.
ztop = 0.
time = 0.
dt = 1.
wmax = 1.
kkmax = 10
deltaK = 20
ilastprint=0
L = 1 ! L initialization
!--- initialization
CALL INITIAL(kmt)
!--- initial print fields:
izprint = 0 ! if = 0 => no printout
if (izprint.ne.0) then
write(cixx(1:2),'(i2.2)') ixx
open(2, file = 'debug.'//cixx//'.dat')
open(19,file='plumegen9.'//cixx//'.gra', &
form='unformatted',access='direct',status='unknown', &
recl=4*nrectotal) !PC
! recl=1*nrectotal) !sx6 e tupay
call printout (izprint,nrectotal)
ilastprint=2
endif
! ******************* model evolution ******************************
rmaxtime = float(maxtime)
!
!print * ,' TIME=',time,' RMAXTIME=',rmaxtime
!print*,'======================================================='
DO WHILE (TIME.LE.RMAXTIME) !beginning of time loop
! do itime=1,120
!-- set model top integration
nm1 = min(kmt, kkmax + deltak)
!-- set timestep
!dt = (zm(2)-zm(1)) / (tstpf * wmax)
dt = min(5.,(zm(2)-zm(1)) / (tstpf * wmax))
!-- elapsed time, sec
time = time+dt
!-- elapsed time, minutes
mintime = 1 + int (time) / 60
wmax = 1. !no zeroes allowed.
!************************** BEGIN SPACE LOOP **************************
!-- zerout all model tendencies
call tend0_plumerise
!-- bounday conditions (k=1)
L=1
call lbound()
!-- dynamics for the level k>1
!-- W advection
! call vel_advectc_plumerise(NM1,WC,WT,DNE,DZM)
call vel_advectc_plumerise(NM1,WC,WT,RHO,DZM)
!-- scalars advection 1
call scl_advectc_plumerise('SC',NM1)
!-- scalars advection 2
!call scl_advectc_plumerise2('SC',NM1)
!-- scalars entrainment, adiabatic
call scl_misc(NM1)
!-- scalars dinamic entrainment
call scl_dyn_entrain(NM1,nkp,wbar,w,adiabat,alpha,radius,tt,t,te,qvt,qv,qvenv,qct,qc,qht,qh,qit,qi,&
vel_e,vel_p,vel_t,rad_p,rad_t)
!-- gravity wave damping using Rayleigh friction layer fot T
call damp_grav_wave(1,nm1,deltak,dt,zt,zm,w,t,tt,qv,qh,qi,qc,te,pe,qvenv)
!-- microphysics
! goto 101 ! bypass microphysics
dt_save=dt
n_sub_step=3
dt=dt/float(n_sub_step)
do i_micro=1,n_sub_step
!-- sedim ?
call fallpart(NM1)
!-- microphysics
do L=2,nm1-1
WBAR = 0.5*(W(L)+W(L-1))
ES = ESAT_PR (T(L)) !BLOB SATURATION VAPOR PRESSURE, EM KPA
QSAT(L) = (EPS * ES) / (PE(L) - ES) !BLOB SATURATION LWC G/G DRY AIR
EST (L) = ES
RHO (L) = 3483.8 * PE (L) / T (L) ! AIR PARCEL DENSITY , G/M**3
!srf18jun2005
! IF (W(L) .ge. 0.) DQSDZ = (QSAT(L ) - QSAT(L-1)) / (ZT(L ) -ZT(L-1))
! IF (W(L) .lt. 0.) DQSDZ = (QSAT(L+1) - QSAT(L )) / (ZT(L+1) -ZT(L ))
IF (W(L) .ge. 0.) then
DQSDZ = (QSAT(L+1) - QSAT(L-1)) / (ZT(L+1 )-ZT(L-1))
ELSE
DQSDZ = (QSAT(L+1) - QSAT(L-1)) / (ZT(L+1) -ZT(L-1))
ENDIF
call waterbal
enddo
enddo
dt=dt_save
!
101 continue
!
!-- W-viscosity for stability
call visc_W(nm1,deltak,kmt)
!-- update scalars
call update_plumerise(nm1,'S')
call hadvance_plumerise(1,nm1,dt,WC,WT,W,mintime)
!-- Buoyancy
call buoyancy_plumerise(NM1, T, TE, QV, QVENV, QH, QI, QC, WT, SCR1)
!-- Entrainment
call entrainment(NM1,W,WT,RADIUS,ALPHA)
!-- update W
call update_plumerise(nm1,'W')
call hadvance_plumerise(2,nm1,dt,WC,WT,W,mintime)
!-- misc
do k=2,nm1
! pe esta em kpa - esat do rams esta em mbar = 100 Pa = 0.1 kpa
! es = 0.1*esat (t(k)) !blob saturation vapor pressure, em kPa
! rotina do plumegen calcula em kPa
es = esat_pr (t(k)) !blob saturation vapor pressure, em kPa
qsat(k) = (eps * es) / (pe(k) - es) !blob saturation lwc g/g dry air
est (k) = es
txs (k) = t(k) - te(k)
rho (k) = 3483.8 * pe (k) / t (k) ! air parcel density , g/m**3
! no pressure diff with radius
if((abs(wc(k))).gt.wmax) wmax = abs(wc(k)) ! keep wmax largest w
enddo
! Gravity wave damping using Rayleigh friction layer for W
call damp_grav_wave(2,nm1,deltak,dt,zt,zm,w,t,tt,qv,qh,qi,qc,te,pe,qvenv)
!---
!- update radius
do k=2,nm1
radius(k) = rad_p(k)
enddo
!-- try to find the plume top (above surface height)
kk = 1
DO WHILE (w (kk) .GT. 1.)
kk = kk + 1
ztop = zm(kk)
!print*,'W=',w (kk)
ENDDO
!
ztop_(mintime) = ztop
ztopmax = MAX (ztop, ztopmax)
kkmax = MAX (kk , kkmax )
!print * ,'ztopmax=', mintime,'mn ',ztop_(mintime), ztopmax
!
! if the solution is going to a stationary phase, exit
IF(mintime > 10) THEN
! if(mintime > 20) then
! if( abs(ztop_(mintime)-ztop_(mintime-10)) < DZ ) exit
IF( ABS(ztop_(mintime)-ztop_(mintime-10)) < DELZ_THRESOLD) then
!- determine W parameter to determine the VMD
do k=2,nm1
W_VMD(k,imm) = w(k)
enddo
EXIT ! finish the integration
ENDIF
ENDIF
if(ilastprint == mintime) then
call printout (izprint,nrectotal)
ilastprint = mintime+1
endif
ENDDO !do next timestep
!print * ,' ztopmax=',ztopmax,'m',mintime,'mn '
!print*,'======================================================='
!
!the last printout
if (izprint.ne.0) then
call printout (izprint,nrectotal)
close (2)
close (19)
endif
RETURN
END SUBROUTINE MAKEPLUME
!-------------------------------------------------------------------------------
!
SUBROUTINE BURN(EFLUX, WATER)
!
!- calculates the energy flux and water content at lboundary
!use module_zero_plumegen_coms
!real, parameter :: HEAT = 21.E6 !Joules/kg
!real, parameter :: HEAT = 15.5E6 !Joules/kg - cerrado
real, parameter :: HEAT = 19.3E6 !Joules/kg - floresta em Alta Floresta (MT)
real :: eflux,water
!
! The emission factor for water is 0.5. The water produced, in kg,
! is then fuel mass*0.5 + (moist/100)*mass per square meter.
! The fire burns for DT out of TDUR seconds, the total amount of
! fuel burned is AREA*BLOAD*(DT/TDUR) kg. this amount of fuel is
! considered to be spread over area AREA and so the mass burned per
! unit area is BLOAD*(DT/TDUR), and the rate is BLOAD/TDUR.
!
IF (TIME.GT.TDUR) THEN !is the burn over?
EFLUX = 0.000001 !prevent a potential divide by zero
WATER = 0.
RETURN
ELSE
!
EFLUX = HEATING (MINTIME) ! Watts/m**2
! WATER = EFLUX * (DT / HEAT) * (0.5 + FMOIST) ! kg/m**2
WATER = EFLUX * (DT / HEAT) * (0.5 + FMOIST) /0.55 ! kg/m**2
WATER = WATER * 1000. ! g/m**2
!
! print*,'BURN:',time,EFLUX/1.e+9
ENDIF
!
RETURN
END SUBROUTINE BURN
!-------------------------------------------------------------------------------
!
SUBROUTINE LBOUND ()
!
! ********** BOUNDARY CONDITIONS AT ZSURF FOR PLUME AND CLOUD ********
!
! source of equations: J.S. Turner Buoyancy Effects in Fluids
! Cambridge U.P. 1973 p.172,
! G.A. Briggs Plume Rise, USAtomic Energy Commissio
! TID-25075, 1969, P.28
!
! fundamentally a point source below ground. at surface, this produces
! a velocity w(1) and temperature T(1) which vary with time. There is
! also a water load which will first saturate, then remainder go into
! QC(1).
! EFLUX = energy flux at ground,watt/m**2 for the last DT
!
!use module_zero_plumegen_coms
implicit none
real, parameter :: g = 9.80796, r = 287.04, cp = 1004.6, eps = 0.622,tmelt = 273.3
real, parameter :: tfreeze = 269.3, pi = 3.14159, e1 = 1./3., e2 = 5./3.
real :: es, esat, eflux, water, pres, c1, c2, f, zv, denscor, xwater !,ESAT_PR
! real, external:: esat_pr!
!
QH (1) = QH (2) !soak up hydrometeors
QI (1) = QI (2)
QC (1) = 0. !no cloud here
!
!
CALL BURN (EFLUX, WATER)
!
! calculate parameters at boundary from a virtual buoyancy point source
!
PRES = PE (1) * 1000. !need pressure in N/m**2
C1 = 5. / (6. * ALPHA) !alpha is entrainment constant
C2 = 0.9 * ALPHA
F = EFLUX / (PRES * CP * PI)
F = G * R * F * AREA !buoyancy flux
ZV = C1 * RSURF !virtual boundary height
W (1) = C1 * ( (C2 * F) **E1) / ZV**E1 !boundary velocity
DENSCOR = C1 * F / G / (C2 * F) **E1 / ZV**E2 !density correction
T (1) = TE (1) / (1. - DENSCOR) !temperature of virtual plume at zsurf
!
WC(1) = W(1)
VEL_P(1) = 0.
rad_p(1) = rsurf
!SC(1) = SCE(1)+F/1000.*dt ! gas/particle (g/g)
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! match dw/dz,dt/dz at the boundary. F is conserved.
!
!WBAR = W (1) * (1. - 1. / (6. * ZV) )
!ADVW = WBAR * W (1) / (3. * ZV)
!ADVT = WBAR * (5. / (3. * ZV) ) * (DENSCOR / (1. - DENSCOR) )
!ADVC = 0.
!ADVH = 0.
!ADVI = 0.
!ADIABAT = - WBAR * G / CP
VTH (1) = - 4.
VTI (1) = - 3.
TXS (1) = T (1) - TE (1)
VISC (1) = VISCOSITY
RHO (1) = 3483.8 * PE (1) / T (1) !air density at level 1, g/m**3
XWATER = WATER / (W (1) * DT * RHO (1) ) !firewater mixing ratio
QV (1) = XWATER + QVENV (1) !plus what's already there
! PE esta em kPa - ESAT do RAMS esta em mbar = 100 Pa = 0.1 kPa
! ES = 0.1*ESAT (T(1)) !blob saturation vapor pressure, em kPa
! rotina do plumegen ja calcula em kPa
ES = ESAT_PR (T(1)) !blob saturation vapor pressure, em kPa
EST (1) = ES
QSAT (1) = (EPS * ES) / (PE (1) - ES) !blob saturation lwc g/g dry air
IF (QV (1) .gt. QSAT (1) ) THEN
QC (1) = QV (1) - QSAT (1) + QC (1) !remainder goes into cloud drops
QV (1) = QSAT (1)
ENDIF
!
CALL WATERBAL
!
RETURN
END SUBROUTINE LBOUND
!-------------------------------------------------------------------------------
!
SUBROUTINE INITIAL ( kmt)
!
! ************* SETS UP INITIAL CONDITIONS FOR THE PROBLEM ************
!use module_zero_plumegen_coms
implicit none
real, parameter :: tfreeze = 269.3
integer :: isub, k, n1, n2, n3, lbuoy, itmp, isubm1 ,kmt
real :: xn1, xi, es, esat!,ESAT_PR
!
N=kmt
! initialize temperature structure,to the end of equal spaced sounding,
do k = 1, N
TXS (k) = 0.0
W (k) = 0.0
T (k) = TE(k) !blob set to environment
WC(k) = 0.0
WT(k) = 0.0
QV(k) = QVENV (k) !blob set to environment
VTH(k) = 0. !initial rain velocity = 0
VTI(k) = 0. !initial ice velocity = 0
QH(k) = 0. !no rain
QI(k) = 0. !no ice
QC(k) = 0. !no cloud drops
! PE esta em kPa - ESAT do RAMS esta em mbar = 100 Pa = 0.1 kPa
! ES = 0.1*ESAT (T(k)) !blob saturation vapor pressure, em kPa
! rotina do plumegen calcula em kPa
ES = ESAT_PR (T(k)) !blob saturation vapor pressure, em kPa
EST (k) = ES
QSAT (k) = (.622 * ES) / (PE (k) - ES) !saturation lwc g/g
RHO (k) = 3483.8 * PE (k) / T (k) !dry air density g/m**3
VEL_P(k) = 0.
rad_p(k) = 0.
enddo
! Initialize the entrainment radius, Turner-style plume
radius(1) = rsurf
do k=2,N
radius(k) = radius(k-1)+(6./5.)*alpha*(zt(k)-zt(k-1))
enddo
! Initialize the entrainment radius, Turner-style plume
radius(1) = rsurf
rad_p(1) = rsurf
DO k=2,N
radius(k) = radius(k-1)+(6./5.)*alpha*(zt(k)-zt(k-1))
rad_p(k) = radius(k)
ENDDO
! Initialize the viscosity
VISC (1) = VISCOSITY
do k=2,N
!VISC (k) = VISCOSITY!max(1.e-3,visc(k-1) - 1.* VISCOSITY/float(nkp))
VISC (k) = max(1.e-3,visc(k-1) - 1.* VISCOSITY/float(nkp))
enddo
!-- Initialize gas/concentration
!DO k =10,20
! SC(k) = 20.
!ENDDO
!stop 333
CALL LBOUND()
RETURN
END SUBROUTINE INITIAL
!-------------------------------------------------------------------------------
!
subroutine damp_grav_wave(ifrom,nm1,deltak,dt,zt,zm,w,t,tt,qv,qh,qi,qc,te,pe,qvenv)
implicit none
integer nm1,ifrom,deltak
real dt
real, dimension(nm1) :: w,t,tt,qv,qh,qi,qc,te,pe,qvenv,dummy,zt,zm
if(ifrom==1) then
call friction(ifrom,nm1,deltak,dt,zt,zm,t,tt ,te)
!call friction(ifrom,nm1,dt,zt,zm,qv,qvt,qvenv)
return
endif
dummy(:) = 0.
if(ifrom==2) call friction(ifrom,nm1,deltak,dt,zt,zm,w,dummy ,dummy)
!call friction(ifrom,nm1,dt,zt,zm,qi,qit ,dummy)
!call friction(ifrom,nm1,dt,zt,zm,qh,qht ,dummy)
!call friction(ifrom,nm1,dt,zt,zm,qc,qct ,dummy)
return
end subroutine damp_grav_wave
!-------------------------------------------------------------------------------
!
subroutine friction(ifrom,nm1,deltak,dt,zt,zm,var1,vart,var2)
implicit none
real, dimension(nm1) :: var1,var2,vart,zt,zm
integer k,nfpt,kf,nm1,ifrom,deltak
real zmkf,ztop,distim,c1,c2,dt
!nfpt=50
!kf = nm1 - nfpt
!kf = nm1 - int(deltak/2)
kf = nm1 - int(deltak)
zmkf = zm(kf) !old: float(kf )*dz
ztop = zm(nm1)
!distim = min(4.*dt,200.)
!distim = 60.
distim = min(3.*dt,60.)
c1 = 1. / (distim * (ztop - zmkf))
c2 = dt * c1
if(ifrom == 1) then
do k = nm1,2,-1
if (zt(k) .le. zmkf) cycle
vart(k) = vart(k) + c1 * (zt(k) - zmkf)*(var2(k) - var1(k))
enddo
elseif(ifrom == 2) then
do k = nm1,2,-1
if (zt(k) .le. zmkf) cycle
var1(k) = var1(k) + c2 * (zt(k) - zmkf)*(var2(k) - var1(k))
enddo
endif
return
end subroutine friction
!-------------------------------------------------------------------------------
!
subroutine vel_advectc_plumerise(m1,wc,wt,rho,dzm)
implicit none
integer :: k,m1
real, dimension(m1) :: wc,wt,flxw,dzm,rho
real, dimension(m1) :: dn0 ! var local
real :: c1z
!dzm(:)= 1./dz
dn0(1:m1)=rho(1:m1)*1.e-3 ! converte de cgs para mks
flxw(1) = wc(1) * dn0(1)
do k = 2,m1-1
flxw(k) = wc(k) * .5 * (dn0(k) + dn0(k+1))
enddo
! Compute advection contribution to W tendency
c1z = .5
do k = 2,m1-2
wt(k) = wt(k) &
+ c1z * dzm(k) / (dn0(k) + dn0(k+1)) * ( &
(flxw(k) + flxw(k-1)) * (wc(k) + wc(k-1)) &
- (flxw(k) + flxw(k+1)) * (wc(k) + wc(k+1)) &
+ (flxw(k+1) - flxw(k-1)) * 2.* wc(k) )
enddo
return
end subroutine vel_advectc_plumerise
!-------------------------------------------------------------------------------
!
subroutine hadvance_plumerise(iac,m1,dt,wc,wt,wp,mintime)
implicit none
integer :: k,iac
integer :: m1,mintime
real, dimension(m1) :: dummy, wc,wt,wp
real eps,dt
! It is here that the Asselin filter is applied. For the velocities
! and pressure, this must be done in two stages, the first when
! IAC=1 and the second when IAC=2.
eps = .2
if(mintime == 1) eps=0.5
! For both IAC=1 and IAC=2, call PREDICT for U, V, W, and P.
!
call predict_plumerise(m1,wc,wp,wt,dummy,iac,2.*dt,eps)
!print*,'mintime',mintime,eps
!do k=1,m1
! print*,'W-HAD',k,wc(k),wp(k),wt(k)
!enddo
return
end subroutine hadvance_plumerise
!-------------------------------------------------------------------------------
!
subroutine predict_plumerise(npts,ac,ap,fa,af,iac,dtlp,epsu)
implicit none
integer :: npts,iac,m
real :: epsu,dtlp
real, dimension(*) :: ac,ap,fa,af
! For IAC=3, this routine moves the arrays AC and AP forward by
! 1 time level by adding in the prescribed tendency. It also
! applies the Asselin filter given by:
! {AC} = AC + EPS * (AP - 2 * AC + AF)
! where AP,AC,AF are the past, current and future time levels of A.
! All IAC=1 does is to perform the {AC} calculation without the AF
! term present. IAC=2 completes the calculation of {AC} by adding
! the AF term only, and advances AC by filling it with input AP
! values which were already updated in ACOUSTC.
!
if (iac .eq. 1) then
do m = 1,npts
ac(m) = ac(m) + epsu * (ap(m) - 2. * ac(m))
enddo
return
elseif (iac .eq. 2) then
do m = 1,npts
af(m) = ap(m)
ap(m) = ac(m) + epsu * af(m)
enddo
!elseif (iac .eq. 3) then
! do m = 1,npts
! af(m) = ap(m) + dtlp * fa(m)
! enddo
! if (ngrid .eq. 1 .and. ipara .eq. 0) call cyclic(nzp,nxp,nyp,af,'T')
! do m = 1,npts
! ap(m) = ac(m) + epsu * (ap(m) - 2. * ac(m) + af(m))
! enddo
endif
do m = 1,npts
ac(m) = af(m)
enddo
return
end subroutine predict_plumerise
!-------------------------------------------------------------------------------
!
subroutine buoyancy_plumerise(m1, T, TE, QV, QVENV, QH, QI, QC, WT, scr1)
implicit none
integer :: k,m1
real, parameter :: g = 9.8, eps = 0.622, gama = 0.5 ! mass virtual coeff.
real, dimension(m1) :: T, TE, QV, QVENV, QH, QI, QC, WT, scr1
real :: TV,TVE,QWTOTL,umgamai
real, parameter :: mu = 0.15
!- orig
umgamai = 1./(1.+gama) ! compensa a falta do termo de aceleracao associado `as
! das pertubacoes nao-hidrostaticas no campo de pressao
!- new ! Siesbema et al, 2004
!umgamai = 1./(1.-2.*mu)
do k = 2,m1-1
TV = T(k) * (1. + (QV(k) /EPS))/(1. + QV(k) ) !blob virtual temp.
TVE = TE(k) * (1. + (QVENV(k)/EPS))/(1. + QVENV(k)) !and environment
QWTOTL = QH(k) + QI(k) + QC(k) ! QWTOTL*G is drag
!- orig
!scr1(k)= G*( umgamai*( TV - TVE) / TVE - QWTOTL)
scr1(k)= G* umgamai*( (TV - TVE) / TVE - QWTOTL)
!if(k .lt. 10)print*,'BT',k,TV,TVE,TVE,QWTOTL
enddo
do k = 2,m1-2
wt(k) = wt(k)+0.5*(scr1(k)+scr1(k+1))
! print*,'W-BUO',k,wt(k),scr1(k),scr1(k+1)
enddo
end subroutine buoyancy_plumerise
!-------------------------------------------------------------------------------
!
subroutine ENTRAINMENT(m1,w,wt,radius,ALPHA)
implicit none
integer :: k,m1
real, dimension(m1) :: w,wt,radius
REAL DMDTM,WBAR,RADIUS_BAR,umgamai,DYN_ENTR,ALPHA
real, parameter :: mu = 0.15 ,gama = 0.5 ! mass virtual coeff.
!- new - Siesbema et al, 2004
!umgamai = 1./(1.-2.*mu)
!- orig
!umgamai = 1
umgamai = 1./(1.+gama) ! compensa a falta do termo de aceleracao associado `as
! das pertubacoes nao-hidrostaticas no campo de pressao
!
!-- ALPHA/RADIUS(L) = (1/M)DM/DZ (W 14a)
do k=2,m1-1
!-- for W: WBAR is only W(k)
! WBAR=0.5*(W(k)+W(k-1))
WBAR=W(k)
RADIUS_BAR = 0.5*(RADIUS(k) + RADIUS(k-1))
! orig
!DMDTM = 2. * ALPHA * ABS (WBAR) / RADIUS_BAR != (1/M)DM/DT
DMDTM = umgamai * 2. * ALPHA * ABS (WBAR) / RADIUS_BAR != (1/M)DM/DT
!-- DMDTM*W(L) entrainment,
wt(k) = wt(k) - DMDTM*ABS (WBAR)
!print*,'W-ENTR=',k,w(k),- DMDTM*ABS (WBAR)
!if(VEL_P (k) - VEL_E (k) > 0.) cycle
!- dynamic entrainment
DYN_ENTR = (2./3.1416)*0.5*ABS (VEL_P(k)-VEL_E(k)+VEL_P(k-1)-VEL_E(k-1)) /RADIUS_BAR
wt(k) = wt(k) - DYN_ENTR*ABS (WBAR)
!- entraiment acceleration for output only
!dwdt_entr(k) = - DMDTM*ABS (WBAR)- DYN_ENTR*ABS (WBAR)
enddo
end subroutine ENTRAINMENT
!-------------------------------------------------------------------------------
!
subroutine scl_advectc_plumerise(varn,mzp)
!use module_zero_plumegen_coms
implicit none
integer :: mzp
character(len=*) :: varn
real :: dtlto2
integer :: k
! wp => w
!- Advect scalars
dtlto2 = .5 * dt
! vt3dc(1) = (w(1) + wc(1)) * dtlto2 * dne(1)
vt3dc(1) = (w(1) + wc(1)) * dtlto2 * rho(1)*1.e-3!converte de CGS p/ MKS
vt3df(1) = .5 * (w(1) + wc(1)) * dtlto2 * dzm(1)
do k = 2,mzp
! vt3dc(k) = (w(k) + wc(k)) * dtlto2 *.5 * (dne(k) + dne(k+1))
vt3dc(k) = (w(k) + wc(k)) * dtlto2 *.5 * (rho(k) + rho(k+1))*1.e-3
vt3df(k) = (w(k) + wc(k)) * dtlto2 *.5 * dzm(k)
!print*,'vt3df-vt3dc',k,vt3dc(k),vt3df(k)
enddo
!-srf-24082005
! do k = 1,mzp-1
do k = 1,mzp
vctr1(k) = (zt(k+1) - zm(k)) * dzm(k)
vctr2(k) = (zm(k) - zt(k)) * dzm(k)
! vt3dk(k) = dzt(k) / dne(k)
vt3dk(k) = dzt(k) /(rho(k)*1.e-3)
!print*,'VT3dk',k,dzt(k) , dne(k)
enddo
! scalarp => scalar_tab(n,ngrid)%var_p
! scalart => scalar_tab(n,ngrid)%var_t
!- temp advection tendency (TT)
scr1=T
call fa_zc_plumerise(mzp &
,T ,scr1 (1) &
,vt3dc (1) ,vt3df (1) &
,vt3dg (1) ,vt3dk (1) &
,vctr1,vctr2 )
call advtndc_plumerise(mzp,T,scr1(1),TT,dt)
!- water vapor advection tendency (QVT)
scr1=QV
call fa_zc_plumerise(mzp &
,QV ,scr1 (1) &
,vt3dc (1) ,vt3df (1) &
,vt3dg (1) ,vt3dk (1) &
,vctr1,vctr2 )
call advtndc_plumerise(mzp,QV,scr1(1),QVT,dt)
!- liquid advection tendency (QCT)
scr1=QC
call fa_zc_plumerise(mzp &
,QC ,scr1 (1) &
,vt3dc (1) ,vt3df (1) &
,vt3dg (1) ,vt3dk (1) &
,vctr1,vctr2 )
call advtndc_plumerise(mzp,QC,scr1(1),QCT,dt)
!- ice advection tendency (QIT)
scr1=QI
call fa_zc_plumerise(mzp &
,QI ,scr1 (1) &
,vt3dc (1) ,vt3df (1) &
,vt3dg (1) ,vt3dk (1) &
,vctr1,vctr2 )
call advtndc_plumerise(mzp,QI,scr1(1),QIT,dt)
!- hail/rain advection tendency (QHT)
! if(ak1 > 0. .or. ak2 > 0.) then
scr1=QH
call fa_zc_plumerise(mzp &
,QH ,scr1 (1) &
,vt3dc (1) ,vt3df (1) &
,vt3dg (1) ,vt3dk (1) &
,vctr1,vctr2 )
call advtndc_plumerise(mzp,QH,scr1(1),QHT,dt)
! endif
!- horizontal wind advection tendency (VEL_T)
scr1=VEL_P
call fa_zc_plumerise(mzp &
,VEL_P ,scr1 (1) &
,vt3dc (1) ,vt3df (1) &
,vt3dg (1) ,vt3dk (1) &
,vctr1,vctr2 )
call advtndc_plumerise(mzp,VEL_P,scr1(1),VEL_T,dt)
!- vertical radius transport
scr1=rad_p
call fa_zc_plumerise(mzp &
,rad_p ,scr1 (1) &
,vt3dc (1) ,vt3df (1) &
,vt3dg (1) ,vt3dk (1) &
,vctr1,vctr2 )
call advtndc_plumerise(mzp,rad_p,scr1(1),rad_t,dt)
return
!
!- gas/particle advection tendency (SCT)
! if(varn == 'SC')return
scr1=SC
call fa_zc_plumerise(mzp &
,SC ,scr1 (1) &
,vt3dc (1) ,vt3df (1) &
,vt3dg (1) ,vt3dk (1) &
,vctr1,vctr2 )
call advtndc_plumerise(mzp,SC,scr1(1),SCT,dt)
return
end subroutine scl_advectc_plumerise
!-------------------------------------------------------------------------------
!
subroutine fa_zc_plumerise(m1,scp,scr1,vt3dc,vt3df,vt3dg,vt3dk,vctr1,vctr2)
implicit none
integer :: m1,k
real :: dfact
real, dimension(m1) :: scp,scr1,vt3dc,vt3df,vt3dg,vt3dk
real, dimension(m1) :: vctr1,vctr2
dfact = .5
! Compute scalar flux VT3DG
do k = 1,m1-1
vt3dg(k) = vt3dc(k) &
* (vctr1(k) * scr1(k) &
+ vctr2(k) * scr1(k+1) &
+ vt3df(k) * (scr1(k) - scr1(k+1)))
enddo
! Modify fluxes to retain positive-definiteness on scalar quantities.
! If a flux will remove 1/2 quantity during a timestep,
! reduce to first order flux. This will remain positive-definite
! under the assumption that ABS(CFL(i)) + ABS(CFL(i-1)) < 1.0 if
! both fluxes are evacuating the box.
do k = 1,m1-1
if (vt3dc(k) .gt. 0.) then
if (vt3dg(k) * vt3dk(k) .gt. dfact * scr1(k)) then
vt3dg(k) = vt3dc(k) * scr1(k)
endif
elseif (vt3dc(k) .lt. 0.) then
if (-vt3dg(k) * vt3dk(k+1) .gt. dfact * scr1(k+1)) then
vt3dg(k) = vt3dc(k) * scr1(k+1)
endif
endif
enddo
! Compute flux divergence
do k = 2,m1-1
scr1(k) = scr1(k) &
+ vt3dk(k) * ( vt3dg(k-1) - vt3dg(k) &
+ scp (k) * ( vt3dc(k) - vt3dc(k-1)))
enddo
return
end subroutine fa_zc_plumerise
!-------------------------------------------------------------------------------
!
subroutine advtndc_plumerise(m1,scp,sca,sct,dtl)
implicit none
integer :: m1,k
real :: dtl,dtli
real, dimension(m1) :: scp,sca,sct
dtli = 1. / dtl
do k = 2,m1-1
sct(k) = sct(k) + (sca(k)-scp(k)) * dtli
enddo
return
end subroutine advtndc_plumerise
!-------------------------------------------------------------------------------
!
subroutine tend0_plumerise
!use module_zero_plumegen_coms, only: nm1,wt,tt,qvt,qct,qht,qit,sct
wt(1:nm1) = 0.
tt(1:nm1) = 0.
qvt(1:nm1) = 0.
qct(1:nm1) = 0.
qht(1:nm1) = 0.
qit(1:nm1) = 0.
vel_t(1:nm1) = 0.
rad_t(1:nm1) = 0.
!sct(1:nm1) = 0.
end subroutine tend0_plumerise
! ****************************************************************
subroutine scl_misc(m1)
!use module_zero_plumegen_coms
implicit none
real, parameter :: g = 9.81, cp=1004.
integer m1,k
real dmdtm
do k=2,m1-1
WBAR = 0.5*(W(k)+W(k-1))
!-- dry adiabat
ADIABAT = - WBAR * G / CP
!
!-- entrainment
DMDTM = 2. * ALPHA * ABS (WBAR) / RADIUS (k) != (1/M)DM/DT
!-- tendency temperature = adv + adiab + entrainment
TT(k) = TT(K) + ADIABAT - DMDTM * ( T (k) - TE (k) )
!-- tendency water vapor = adv + entrainment
QVT(K) = QVT(K) - DMDTM * ( QV (k) - QVENV (k) )
QCT(K) = QCT(K) - DMDTM * ( QC (k) )
QHT(K) = QHT(K) - DMDTM * ( QH (k) )
QIT(K) = QIT(K) - DMDTM * ( QI (k) )
!-- tendency horizontal speed = adv + entrainment
VEL_T(K) = VEL_T(K) - DMDTM * ( VEL_P (k) - VEL_E (k) )
!-- tendency horizontal speed = adv + entrainment
rad_t(K) = rad_t(K) + 0.5*DMDTM*(6./5.)*RADIUS (k)
!-- tendency gas/particle = adv + entrainment
! SCT(K) = SCT(K) - DMDTM * ( SC (k) - SCE (k) )
enddo
end subroutine scl_misc
! ****************************************************************
SUBROUTINE scl_dyn_entrain(m1,nkp,wbar,w,adiabat,alpha,radius,tt,t,te,qvt,qv,qvenv,qct,qc,qht,qh,qit,qi,&
vel_e,vel_p,vel_t,rad_p,rad_t)
implicit none
INTEGER , INTENT(IN) :: m1
! plumegen_coms
INTEGER , INTENT(IN) :: nkp
REAL , INTENT(INOUT) :: wbar
REAL , INTENT(IN) :: w(nkp)
REAL , INTENT(INOUT) :: adiabat
REAL , INTENT(IN) :: alpha
REAL , INTENT(IN) :: radius(nkp)
REAL , INTENT(INOUT) :: tt(nkp)
REAL , INTENT(IN) :: t(nkp)
REAL , INTENT(IN) :: te(nkp)
REAL , INTENT(INOUT) :: qvt(nkp)
REAL , INTENT(IN) :: qv(nkp)
REAL , INTENT(IN) :: qvenv(nkp)
REAL , INTENT(INOUT) :: qct(nkp)
REAL , INTENT(IN) :: qc(nkp)
REAL , INTENT(INOUT) :: qht(nkp)
REAL , INTENT(IN) :: qh(nkp)
REAL , INTENT(INOUT) :: qit(nkp)
REAL , INTENT(IN) :: qi(nkp)
REAL , INTENT(IN) :: vel_e(nkp)
REAL , INTENT(IN) :: vel_p(nkp)
REAL , INTENT(INOUT) :: vel_t(nkp)
REAL , INTENT(INOUT) :: rad_T(nkp)
REAL , INTENT(IN) :: rad_p(nkp)
real, parameter :: g = 9.81, cp=1004., pi=3.1416
integer k
real dmdtm
DO k=2,m1-1
!
!-- tendency horizontal radius from dyn entrainment
!rad_t(K) = rad_t(K) + (vel_e(k)-vel_p(k)) /pi
rad_t(K) = rad_t(K) + ABS((vel_e(k)-vel_p(k)))/pi
!-- entrainment
!DMDTM = (2./3.1416) * (VEL_E (k) - VEL_P (k)) / RADIUS (k)
DMDTM = (2./3.1416) * ABS(VEL_E (k) - VEL_P (k)) / RADIUS (k)
!-- tendency horizontal speed from dyn entrainment
VEL_T(K) = VEL_T(K) - DMDTM * ( VEL_P (k) - VEL_E (k) )
! if(VEL_P (k) - VEL_E (k) > 0.) cycle
!-- tendency temperature from dyn entrainment
TT(k) = TT(K) - DMDTM * ( T (k) - TE (k) )
!-- tendency water vapor from dyn entrainment
QVT(K) = QVT(K) - DMDTM * ( QV (k) - QVENV (k) )
QCT(K) = QCT(K) - DMDTM * ( QC (k) )
QHT(K) = QHT(K) - DMDTM * ( QH (k) )
QIT(K) = QIT(K) - DMDTM * ( QI (k) )
!-- tendency gas/particle from dyn entrainment
! SCT(K) = SCT(K) - DMDTM * ( SC (k) - SCE (k) )
ENDDO
END SUBROUTINE scl_dyn_entrain
! ****************************************************************
subroutine visc_W(m1,deltak,kmt)
!use module_zero_plumegen_coms
implicit none
integer m1,k,deltak,kmt,m2
real dz1t,dz1m,dz2t,dz2m,d2wdz,d2tdz ,d2qvdz ,d2qhdz ,d2qcdz ,d2qidz ,d2scdz, &
d2vel_pdz,d2rad_dz
!srf--- 17/08/2005
!m2=min(m1+deltak,kmt)
m2=min(m1,kmt)
!do k=2,m1-1
do k=2,m2-1
DZ1T = 0.5*(ZT(K+1)-ZT(K-1))
DZ2T = VISC (k) / (DZ1T * DZ1T)
DZ1M = 0.5*(ZM(K+1)-ZM(K-1))
DZ2M = VISC (k) / (DZ1M * DZ1M)
D2WDZ = (W (k + 1) - 2 * W (k) + W (k - 1) ) * DZ2M
D2TDZ = (T (k + 1) - 2 * T (k) + T (k - 1) ) * DZ2T
D2QVDZ = (QV (k + 1) - 2 * QV (k) + QV (k - 1) ) * DZ2T
D2QHDZ = (QH (k + 1) - 2 * QH (k) + QH (k - 1) ) * DZ2T
D2QCDZ = (QC (k + 1) - 2 * QC (k) + QC (k - 1) ) * DZ2T
D2QIDZ = (QI (k + 1) - 2 * QI (k) + QI (k - 1) ) * DZ2T
!D2SCDZ = (SC (k + 1) - 2 * SC (k) + SC (k - 1) ) * DZ2T
d2vel_pdz=(vel_P (k + 1) - 2 * vel_P (k) + vel_P (k - 1) ) * DZ2T
d2rad_dz =(rad_p (k + 1) - 2 * rad_p (k) + rad_p (k - 1) ) * DZ2T
WT(k) = WT(k) + D2WDZ
TT(k) = TT(k) + D2TDZ
QVT(k) = QVT(k) + D2QVDZ
QCT(k) = QCT(k) + D2QCDZ
QHT(k) = QHT(k) + D2QHDZ
QIT(k) = QIT(k) + D2QIDZ
vel_t(k) = vel_t(k) + d2vel_pdz
rad_t(k) = rad_t(k) + d2rad_dz
!SCT(k) = SCT(k) + D2SCDZ
!print*,'W-VISC=',k,D2WDZ
enddo
end subroutine visc_W
! ****************************************************************
subroutine update_plumerise(m1,varn)
!use module_zero_plumegen_coms
integer m1,k
character(len=*) :: varn
if(varn == 'W') then
do k=2,m1-1
W(k) = W(k) + WT(k) * DT
enddo
return
else
do k=2,m1-1
T(k) = T(k) + TT(k) * DT
QV(k) = QV(k) + QVT(k) * DT
QC(k) = QC(k) + QCT(k) * DT !cloud drops travel with air
QH(k) = QH(k) + QHT(k) * DT
QI(k) = QI(k) + QIT(k) * DT
! SC(k) = SC(k) + SCT(k) * DT
!srf---18jun2005
QV(k) = max(0., QV(k))
QC(k) = max(0., QC(k))
QH(k) = max(0., QH(k))
QI(k) = max(0., QI(k))
VEL_P(k) = VEL_P(k) + VEL_T(k) * DT
rad_p(k) = rad_p(k) + rad_t(k) * DT
! SC(k) = max(0., SC(k))
enddo
endif
end subroutine update_plumerise
!-------------------------------------------------------------------------------
!
subroutine fallpart(m1)
!use module_zero_plumegen_coms
integer m1,k
real vtc, dfhz,dfiz,dz1
!srf==================================
! verificar se o gradiente esta correto
!
!srf==================================
!
! XNO=1.E7 [m**-4] median volume diameter raindrop,Kessler
! VC = 38.3/(XNO**.125), median volume fallspeed eqn., Kessler
! for ice, see (OT18), use F0=0.75 per argument there. rho*q
! values are in g/m**3, velocities in m/s
real, PARAMETER :: VCONST = 5.107387, EPS = 0.622, F0 = 0.75
real, PARAMETER :: G = 9.81, CP = 1004.
!
do k=2,m1-1
VTC = VCONST * RHO (k) **.125 ! median volume fallspeed (KTable4)
! hydrometeor assembly velocity calculations (K Table4)
! VTH(k)=-VTC*QH(k)**.125 !median volume fallspeed, water
VTH (k) = - 4. !small variation with qh
VHREL = W (k) + VTH (k) !relative to surrounding cloud
! rain ventilation coefficient for evaporation
CVH(k) = 1.6 + 0.57E-3 * (ABS (VHREL) ) **1.5
!
! VTI(k)=-VTC*F0*QI(k)**.125 !median volume fallspeed,ice
VTI (k) = - 3. !small variation with qi
VIREL = W (k) + VTI (k) !relative to surrounding cloud
!
! ice ventilation coefficient for sublimation
CVI(k) = 1.6 + 0.57E-3 * (ABS (VIREL) ) **1.5 / F0
!
!
IF (VHREL.GE.0.0) THEN
DFHZ=QH(k)*(RHO(k )*VTH(k )-RHO(k-1)*VTH(k-1))/RHO(k-1)
ELSE
DFHZ=QH(k)*(RHO(k+1)*VTH(k+1)-RHO(k )*VTH(k ))/RHO(k)
ENDIF
!
!
IF (VIREL.GE.0.0) THEN
DFIZ=QI(k)*(RHO(k )*VTI(k )-RHO(k-1)*VTI(k-1))/RHO(k-1)
ELSE
DFIZ=QI(k)*(RHO(k+1)*VTI(k+1)-RHO(k )*VTI(k ))/RHO(k)
ENDIF
DZ1=ZM(K)-ZM(K-1)
qht(k) = qht(k) - DFHZ / DZ1 !hydrometeors don't
qit(k) = qit(k) - DFIZ / DZ1 !nor does ice? hail, what about
enddo
end subroutine fallpart
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!
subroutine printout (izprint,nrectotal)
!use module_zero_plumegen_coms
real, parameter :: tmelt = 273.3
integer, save :: nrec
data nrec/0/
integer :: ko,izprint,interval,nrectotal
real :: pea, btmp,etmp,vap1,vap2,gpkc,gpkh,gpki,deficit
interval = 1 !debug time interval,min
!
IF (IZPRINT.EQ.0) RETURN
IF(MINTIME == 1) nrec = 0
!
WRITE (2, 430) MINTIME, DT, TIME
WRITE (2, 431) ZTOP
WRITE (2, 380)
!
! do the print
!
DO 390 KO = 1, nrectotal, interval
PEA = PE (KO) * 10. !pressure is stored in decibars(kPa),print in mb;
BTMP = T (KO) - TMELT !temps in Celsius
ETMP = T (KO) - TE (KO) !temperature excess
VAP1 = QV (KO) * 1000. !printout in g/kg for all water,
VAP2 = QSAT (KO) * 1000. !vapor (internal storage is in g/g)
GPKC = QC (KO) * 1000. !cloud water
GPKH = QH (KO) * 1000. !raindrops
GPKI = QI (KO) * 1000. !ice particles
DEFICIT = VAP2 - VAP1 !vapor deficit
!
WRITE (2, 400) zt(KO)/1000., PEA, W (KO), BTMP, ETMP, VAP1, &
VAP2, GPKC, GPKH, GPKI, VTH (KO), SC(KO)
!
!
! !end of printout
390 CONTINUE
nrec=nrec+1
write (19,rec=nrec) (W (KO), KO=1,nrectotal)
nrec=nrec+1
write (19,rec=nrec) (T (KO), KO=1,nrectotal)
nrec=nrec+1
write (19,rec=nrec) (TE(KO), KO=1,nrectotal)
nrec=nrec+1
write (19,rec=nrec) (QV(KO)*1000., KO=1,nrectotal)
nrec=nrec+1
write (19,rec=nrec) (QC(KO)*1000., KO=1,nrectotal)
nrec=nrec+1
write (19,rec=nrec) (QH(KO)*1000., KO=1,nrectotal)
nrec=nrec+1
write (19,rec=nrec) (QI(KO)*1000., KO=1,nrectotal)
nrec=nrec+1
! write (19,rec=nrec) (SC(KO), KO=1,nrectotal)
write (19,rec=nrec) (QSAT(KO)*1000., KO=1,nrectotal)
nrec=nrec+1
write (19,rec=nrec) (QVENV(KO)*1000., KO=1,nrectotal)
!print*,'ntimes=',nrec/(9)
!
RETURN
!
! ************** FORMATS *********************************************
!
380 FORMAT(/,' Z(KM) P(MB) W(MPS) T(C) T-TE VAP SAT QC QH' &
' QI VTH(MPS) SCAL'/)
!
400 FORMAT(1H , F4.1,F7.2,F7.2,F6.1,6F6.2,F7.2,1X,F6.2)
!
430 FORMAT(1H ,//I5,' MINUTES DT= ',F6.2,' SECONDS TIME= ' &
,F8.2,' SECONDS')
431 FORMAT(' ZTOP= ',F10.2)
!
end subroutine printout
!
! *********************************************************************
SUBROUTINE WATERBAL
!use module_zero_plumegen_coms
!
IF (QC (L) .LE.1.0E-10) QC (L) = 0. !DEFEAT UNDERFLOW PROBLEM
IF (QH (L) .LE.1.0E-10) QH (L) = 0.
IF (QI (L) .LE.1.0E-10) QI (L) = 0.
!
CALL EVAPORATE !vapor to cloud,cloud to vapor
!
CALL SUBLIMATE !vapor to ice
!
CALL GLACIATE !rain to ice
CALL MELT !ice to rain
!
!if(ak1 > 0. .or. ak2 > 0.) &
CALL CONVERT () !(auto)conversion and accretion
!CALL CONVERT2 () !(auto)conversion and accretion
!
RETURN
END SUBROUTINE WATERBAL
! *********************************************************************
SUBROUTINE EVAPORATE
!
!- evaporates cloud,rain and ice to saturation
!
!use module_zero_plumegen_coms
implicit none
!
! XNO=10.0E06
! HERC = 1.93*1.E-6*XN035 !evaporation constant
!
real, PARAMETER :: HERC = 5.44E-4, CP = 1.004, HEATCOND = 2.5E3
real, PARAMETER :: HEATSUBL = 2834., TMELT = 273., TFREEZE = 269.3
real, PARAMETER :: FRC = HEATCOND / CP, SRC = HEATSUBL / CP
real :: evhdt, evidt, evrate, evap, sd, quant, dividend, divisor, devidt
!
!
SD = QSAT (L) - QV (L) !vapor deficit
IF (SD.EQ.0.0) RETURN
!IF (abs(SD).lt.1.e-7) RETURN
EVHDT = 0.
EVIDT = 0.
!evrate =0.; evap=0.; sd=0.0; quant=0.0; dividend=0.0; divisor=0.0; devidt=0.0
EVRATE = ABS (WBAR * DQSDZ) !evaporation rate (Kessler 8.32)
EVAP = EVRATE * DT !what we can get in DT
IF (SD.LE.0.0) THEN ! condense. SD is negative
IF (EVAP.GE.ABS (SD) ) THEN !we get it all
QC (L) = QC (L) - SD !deficit,remember?
QV (L) = QSAT(L) !set the vapor to saturation
T (L) = T (L) - SD * FRC !heat gained through condensation
!per gram of dry air
RETURN
ELSE
QC (L) = QC (L) + EVAP !get what we can in DT
QV (L) = QV (L) - EVAP !remove it from the vapor
T (L) = T (L) + EVAP * FRC !get some heat
RETURN
ENDIF
!
ELSE !SD is positive, need some water
!
! not saturated. saturate if possible. use everything in order
! cloud, rain, ice. SD is positive
IF (EVAP.LE.QC (L) ) THEN !enough cloud to last DT
!
IF (SD.LE.EVAP) THEN !enough time to saturate
QC (L) = QC (L) - SD !remove cloud
QV (L) = QSAT (L) !saturate
T (L) = T (L) - SD * FRC !cool the parcel
RETURN !done
!
ELSE !not enough time
SD = SD-EVAP !use what there is
QV (L) = QV (L) + EVAP !add vapor
T (L) = T (L) - EVAP * FRC !lose heat
QC (L) = QC (L) - EVAP !lose cloud
!go on to rain.
ENDIF
!
ELSE !not enough cloud to last DT
!
IF (SD.LE.QC (L) ) THEN !but there is enough to sat
QV (L) = QSAT (L) !use it
QC (L) = QC (L) - SD
T (L) = T (L) - SD * FRC
RETURN
ELSE !not enough to sat
SD = SD-QC (L)
QV (L) = QV (L) + QC (L)
T (L) = T (L) - QC (L) * FRC
QC (L) = 0.0 !all gone
ENDIF !on to rain
ENDIF !finished with cloud
!
! but still not saturated, so try to use some rain
! this is tricky, because we only have time DT to evaporate. if there
! is enough rain, we can evaporate it for dt. ice can also sublimate
! at the same time. there is a compromise here.....use rain first, then
! ice. saturation may not be possible in one DT time.
! rain evaporation rate (W12),(OT25),(K Table 4). evaporate rain first
! sd is still positive or we wouldn't be here.
IF (QH (L) .LE.1.E-10) GOTO 33
!srf-25082005
! QUANT = ( QC (L) + QV (L) - QSAT (L) ) * RHO (L) !g/m**3
QUANT = ( QSAT (L)- QC (L) - QV (L) ) * RHO (L) !g/m**3
!
EVHDT = (DT * HERC * (QUANT) * (QH (L) * RHO (L) ) **.65) / RHO (L)
! rain evaporation in time DT
IF (EVHDT.LE.QH (L) ) THEN !enough rain to last DT
IF (SD.LE.EVHDT) THEN !enough time to saturate
QH (L) = QH (L) - SD !remove rain
QV (L) = QSAT (L) !saturate
T (L) = T (L) - SD * FRC !cool the parcel
RETURN !done
!
ELSE !not enough time
SD = SD-EVHDT !use what there is
QV (L) = QV (L) + EVHDT !add vapor
T (L) = T (L) - EVHDT * FRC !lose heat
QH (L) = QH (L) - EVHDT !lose rain
ENDIF !go on to ice.
!
ELSE !not enough rain to last DT
!
IF (SD.LE.QH (L) ) THEN !but there is enough to sat
QV (L) = QSAT (L) !use it
QH (L) = QH (L) - SD
T (L) = T (L) - SD * FRC
RETURN
!
ELSE !not enough to sat
SD = SD-QH (L)
QV (L) = QV (L) + QH (L)
T (L) = T (L) - QH (L) * FRC
QH (L) = 0.0 !all gone
ENDIF !on to ice
!
ENDIF !finished with rain
!
!
! now for ice
! equation from (OT); correction factors for units applied
!
33 continue
IF (QI (L) .LE.1.E-10) RETURN !no ice there
!
DIVIDEND = ( (1.E6 / RHO (L) ) **0.475) * (SD / QSAT (L) &
- 1) * (QI (L) **0.525) * 1.13
DIVISOR = 7.E5 + 4.1E6 / (10. * EST (L) )
DEVIDT = - CVI(L) * DIVIDEND / DIVISOR !rate of change
EVIDT = DEVIDT * DT !what we could get
!
! logic here is identical to rain. could get fancy and make subroutine
! but duplication of code is easier. God bless the screen editor.
!
IF (EVIDT.LE.QI (L) ) THEN !enough ice to last DT
!
IF (SD.LE.EVIDT) THEN !enough time to saturate
QI (L) = QI (L) - SD !remove ice
QV (L) = QSAT (L) !saturate
T (L) = T (L) - SD * SRC !cool the parcel
RETURN !done
!
ELSE !not enough time
SD = SD-EVIDT !use what there is
QV (L) = QV (L) + EVIDT !add vapor
T (L) = T (L) - EVIDT * SRC !lose heat
QI (L) = QI (L) - EVIDT !lose ice
ENDIF !go on,unsatisfied
!
ELSE !not enough ice to last DT
!
IF (SD.LE.QI (L) ) THEN !but there is enough to sat
QV (L) = QSAT (L) !use it
QI (L) = QI (L) - SD
T (L) = T (L) - SD * SRC
RETURN
!
ELSE !not enough to sat
SD = SD-QI (L)
QV (L) = QV (L) + QI (L)
T (L) = T (L) - QI (L) * SRC
QI (L) = 0.0 !all gone
ENDIF !on to better things
!finished with ice
ENDIF
!
ENDIF !finished with the SD decision
!
RETURN
!
END SUBROUTINE EVAPORATE
!
! *********************************************************************
SUBROUTINE CONVERT ()
!
!- ACCRETION AND AUTOCONVERSION
!
!use module_zero_plumegen_coms
!
real, PARAMETER :: AK1 = 0.001 !conversion rate constant
real, PARAMETER :: AK2 = 0.0052 !collection (accretion) rate
real, PARAMETER :: TH = 0.5 !Kessler threshold
integer, PARAMETER ::iconv = 1 !- Kessler conversion (=0)
!real, parameter :: ANBASE = 50.!*1.e+6 !Berry-number at cloud base #/m^3(maritime)
real, parameter :: ANBASE =100000.!*1.e+6 !Berry-number at cloud base #/m^3(continental)
!real, parameter :: BDISP = 0.366 !Berry--size dispersion (maritime)
real, parameter :: BDISP = 0.146 !Berry--size dispersion (continental)
real, parameter :: TFREEZE = 269.3 !ice formation temperature
!
real :: accrete, con, q, h, bc1, bc2, total
IF (T (L) .LE. TFREEZE) RETURN !process not allowed above ice
!
IF (QC (L) .EQ. 0. ) RETURN
ACCRETE = 0.
CON = 0.
Q = RHO (L) * QC (L)
H = RHO (L) * QH (L)
!
! selection rules
!
!
IF (QH (L) .GT. 0. ) ACCRETE = AK2 * Q * (H**.875) !accretion, Kessler
!
IF (ICONV.NE.0) THEN !select Berry or Kessler
!
!old BC1 = 120.
!old BC2 = .0266 * ANBASE * 60.
!old CON = BDISP * Q * Q * Q / (BC1 * Q * BDISP + BC2)
CON = Q*Q*Q*BDISP/(60.*(5.*Q*BDISP+0.0366*ANBASE))
!
ELSE
!
! CON = AK1 * (Q - TH) !Kessler autoconversion rate
!
! IF (CON.LT.0.0) CON = 0.0 !havent reached threshold
CON = max(0.,AK1 * (Q - TH)) ! versao otimizada
!
ENDIF
!
!
TOTAL = (CON + ACCRETE) * DT / RHO (L)
!
IF (TOTAL.LT.QC (L) ) THEN
!
QC (L) = QC (L) - TOTAL
QH (L) = QH (L) + TOTAL !no phase change involved
RETURN
!
ELSE
!
QH (L) = QH (L) + QC (L) !uses all there is
QC (L) = 0.0
!
ENDIF
!
RETURN
!
END SUBROUTINE CONVERT
!
!**********************************************************************
!
SUBROUTINE CONVERT2 ()
!use module_zero_plumegen_coms
implicit none
LOGICAL AEROSOL
parameter(AEROSOL=.true.)
!
real, parameter :: TNULL=273.16, LAT=2.5008E6 &
,EPSI=0.622 ,DB=1. ,NB=1500. !ALPHA=0.2
real :: KA,KEINS,KZWEI,KDREI,VT
real :: A,B,C,D, CON,ACCRETE,total
real Y(6),ROH
A=0.
B=0.
Y(1) = T(L)
Y(4) = W(L)
y(2) = QC(L)
y(3) = QH(L)
Y(5) = RADIUS(L)
ROH = RHO(L)*1.e-3 ! dens (MKS) ??
! autoconversion
KA = 0.0005
IF( Y(1) .LT. 258.15 )THEN
! KEINS=0.00075
KEINS=0.0009
KZWEI=0.0052
KDREI=15.39
ELSE
KEINS=0.0015
KZWEI=0.00696
KDREI=11.58
ENDIF
! ROH=PE/RD/TE
VT=-KDREI* (Y(3)/ROH)**0.125
IF (Y(4).GT.0.0 ) THEN
IF (AEROSOL) THEN
A = 1/y(4) * y(2)*y(2)*1000./( 60. *( 5. + 0.0366*NB/(y(2)*1000.*DB) ) )
ELSE
IF (y(2).GT.(KA*ROH)) THEN
!print*,'1',y(2),KA*ROH
A = KEINS/y(4) *(y(2) - KA*ROH )
ENDIF
ENDIF
ELSE
A = 0.0
ENDIF
! accretion
IF(y(4).GT.0.0) THEN
B = KZWEI/(y(4) - VT) * MAX(0.,y(2)) * &
MAX(0.001,ROH)**(-0.875)*(MAX(0.,y(3)))**(0.875)
ELSE
B = 0.0
ENDIF
!PSATW=610.7*EXP( 17.25 *( Y(1) - TNULL )/( Y(1)-36. ) )
!PSATE=610.7*EXP( 22.33 *( Y(1) - TNULL )/( Y(1)- 2. ) )
!QSATW=EPSI*PSATW/( PE-(1.-EPSI)*PSATW )
!QSATE=EPSI*PSATE/( PE-(1.-EPSI)*PSATE )
!MU=2.*ALPHA/Y(5)
!C = MU*( ROH*QSATW - ROH*QVE + y(2) )
!D = ROH*LAT*QSATW*EPSI/Y1/Y1/RD *DYDX1
!DYDX(2) = - A - B - C - D ! d rc/dz
!DYDX(3) = A + B ! d rh/dz
! rc=rc+dydx(2)*dz
! rh=rh+dydx(3)*dz
CON = A
ACCRETE = B
TOTAL = (CON + ACCRETE) *(1/DZM(L)) /ROH ! DT / RHO (L)
!print*,'L=',L,total,QC(L),dzm(l)
!
IF (TOTAL.LT.QC (L) ) THEN
!
QC (L) = QC (L) - TOTAL
QH (L) = QH (L) + TOTAL !no phase change involved
RETURN
!
ELSE
!
QH (L) = QH (L) + QC (L) !uses all there is
QC (L) = 0.0
!
ENDIF
!
RETURN
!
END SUBROUTINE CONVERT2
! ice - effect on temperature
! TTD = 0.0
! TTE = 0.0
! CALL ICE(QSATW,QSATE,Y(1),Y(2),Y(3), &
! TTA,TTB,TTC,DZ,ROH,D,C,TTD,TTE)
! DYDX(1) = DYDX(1) + TTD + TTE ! DT/DZ on Temp
!
!**********************************************************************
!
SUBROUTINE SUBLIMATE
!
! ********************* VAPOR TO ICE (USE EQUATION OT22)***************
!use module_zero_plumegen_coms
!
real, PARAMETER :: EPS = 0.622, HEATFUS = 334., HEATSUBL = 2834., CP = 1.004
real, PARAMETER :: SRC = HEATSUBL / CP, FRC = HEATFUS / CP, TMELT = 273.3
real, PARAMETER :: TFREEZE = 269.3
real ::dtsubh, dividend,divisor, subl
!
DTSUBH = 0.
!
!selection criteria for sublimation
IF (T (L) .GT. TFREEZE ) RETURN
IF (QV (L) .LE. QSAT (L) ) RETURN
!
! from (OT); correction factors for units applied
!
DIVIDEND = ( (1.E6 / RHO (L) ) **0.475) * (QV (L) / QSAT (L) &
- 1) * (QI (L) **0.525) * 1.13
DIVISOR = 7.E5 + 4.1E6 / (10. * EST (L) )
!
DTSUBH = ABS (DIVIDEND / DIVISOR) !sublimation rate
SUBL = DTSUBH * DT !and amount possible
!
! again check the possibilities
!
IF (SUBL.LT.QV (L) ) THEN
!
QV (L) = QV (L) - SUBL !lose vapor
QI (L) = QI (L) + SUBL !gain ice
T (L) = T (L) + SUBL * SRC !energy change, warms air
!print*,'5',l,qi(l),SUBL
RETURN
!
ELSE
!
QI (L) = QV (L) !use what there is
T (L) = T (L) + QV (L) * SRC !warm the air
QV (L) = 0.0
!print*,'6',l,qi(l)
!
ENDIF
!
RETURN
END SUBROUTINE SUBLIMATE
!
! *********************************************************************
!
SUBROUTINE GLACIATE
!
! *********************** CONVERSION OF RAIN TO ICE *******************
! uses equation OT 16, simplest. correction from W not applied, but
! vapor pressure differences are supplied.
!
!use module_zero_plumegen_coms
!
real, PARAMETER :: HEATFUS = 334., CP = 1.004, EPS = 0.622, HEATSUBL = 2834.
real, PARAMETER :: FRC = HEATFUS / CP, FRS = HEATSUBL / CP, TFREEZE = 269.3
real, PARAMETER :: GLCONST = 0.025 !glaciation time constant, 1/sec
real dfrzh
!
DFRZH = 0. !rate of mass gain in ice
!
!selection rules for glaciation
IF (QH (L) .LE. 0. ) RETURN
IF (QV (L) .LT. QSAT (L) ) RETURN
IF (T (L) .GT. TFREEZE ) RETURN
!
! NT=TMELT-T(L)
! IF (NT.GT.50) NT=50
!
DFRZH = DT * GLCONST * QH (L) ! from OT(16)
!
IF (DFRZH.LT.QH (L) ) THEN
!
QI (L) = QI (L) + DFRZH
QH (L) = QH (L) - DFRZH
T (L) = T (L) + FRC * DFRZH !warms air
!print*,'7',l,qi(l),DFRZH
RETURN
!
ELSE
!
QI (L) = QI (L) + QH (L)
T (L) = T (L) + FRC * QH (L)
QH (L) = 0.0
!print*,'8',l,qi(l), QH (L)
!
ENDIF
!
RETURN
!
END SUBROUTINE GLACIATE
!
!
! *********************************************************************
SUBROUTINE MELT
!
! ******************* MAKES WATER OUT OF ICE **************************
!use module_zero_plumegen_coms
!
real, PARAMETER :: FRC = 332.27, TMELT = 273., F0 = 0.75 !ice velocity factor
real DTMELT
!
DTMELT = 0. !conversion,ice to rain
!
!selection rules
IF (QI (L) .LE. 0.0 ) RETURN
IF (T (L) .LT. TMELT) RETURN
!
!OT(23,24)
DTMELT = DT * (2.27 / RHO (L) ) * CVI(L) * (T (L) - TMELT) * ( (RHO(L) &
* QI (L) * 1.E-6) **0.525) * (F0** ( - 0.42) )
!after Mason,1956
!
! check the possibilities
!
IF (DTMELT.LT.QI (L) ) THEN
!
QH (L) = QH (L) + DTMELT
QI (L) = QI (L) - DTMELT
T (L) = T (L) - FRC * DTMELT !cools air
!print*,'9',l,qi(l),DTMELT
RETURN
!
ELSE
!
QH (L) = QH (L) + QI (L) !get all there is to get
T (L) = T (L) - FRC * QI (L)
QI (L) = 0.0
!print*,'10',l,qi(l)
!
ENDIF
!
RETURN
!
END SUBROUTINE MELT
SUBROUTINE htint (nzz1, vctra, eleva, nzz2, vctrb, elevb)
IMPLICIT NONE
INTEGER, INTENT(IN ) :: nzz1
INTEGER, INTENT(IN ) :: nzz2
REAL, INTENT(IN ) :: vctra(nzz1)
REAL, INTENT(OUT) :: vctrb(nzz2)
REAL, INTENT(IN ) :: eleva(nzz1)
REAL, INTENT(IN ) :: elevb(nzz2)
INTEGER :: l
INTEGER :: k
INTEGER :: kk
REAL :: wt
l=1
DO k=1,nzz2
DO
IF ( (elevb(k) < eleva(1)) .OR. &
((elevb(k) >= eleva(l)) .AND. (elevb(k) <= eleva(l+1))) ) THEN
wt = (elevb(k)-eleva(l))/(eleva(l+1)-eleva(l))
vctrb(k) = vctra(l)+(vctra(l+1)-vctra(l))*wt
EXIT
ELSE IF ( elevb(k) > eleva(nzz1)) THEN
wt = (elevb(k)-eleva(nzz1))/(eleva(nzz1-1)-eleva(nzz1))
vctrb(k) = vctra(nzz1)+(vctra(nzz1-1)-vctra(nzz1))*wt
EXIT
END IF
l=l+1
IF(l == nzz1) THEN
PRINT *,'htint:nzz1',nzz1
DO kk=1,l
PRINT*,'kk,eleva(kk),elevb(kk)',eleva(kk),elevb(kk)
END DO
STOP 'htint'
END IF
END DO
END DO
END SUBROUTINE htint
!-----------------------------------------------------------------------------
FUNCTION ESAT_PR (TEM)
!
! ******* Vapor Pressure A.L. Buck JAM V.20 p.1527. (1981) ***********
!
real, PARAMETER :: CI1 = 6.1115, CI2 = 22.542, CI3 = 273.48
real, PARAMETER :: CW1 = 6.1121, CW2 = 18.729, CW3 = 257.87, CW4 = 227.3
real, PARAMETER :: TMELT = 273.3
real ESAT_PR
real temc , tem,esatm
!
! formulae from Buck, A.L., JAM 20,1527-1532
! custom takes esat wrt water always. formula for h2o only
! good to -40C so:
!
!
TEMC = TEM - TMELT
IF (TEMC.GT. - 40.0) GOTO 230
ESATM = CI1 * EXP (CI2 * TEMC / (TEMC + CI3) ) !ice, millibars
ESAT_PR = ESATM / 10. !kPa
RETURN
!
230 ESATM = CW1 * EXP ( ( (CW2 - (TEMC / CW4) ) * TEMC) / (TEMC + CW3))
ESAT_PR = ESATM / 10. !kPa
RETURN
END function ESAT_PR
! ******************************************************************
! ------------------------------------------------------------------------
END Module module_chem_plumerise_scalar