/wrfv2_fire/phys/module_fr_sfire_atm.F
FORTRAN Legacy | 2218 lines | 1267 code | 333 blank | 618 comment | 61 complexity | 01cc9230d1f3bc6663e5572fcf894d4f MD5 | raw file
Possible License(s): AGPL-1.0
Large files files are truncated, but you can click here to view the full file
- !WRF:MEDIATION_LAYER:FIRE_MODEL
- !*** Jan Mandel August 2007 - February 2010
- !*** email: Jan.Mandel@gmail.com
- ! Routines dealing with the atmosphere
- module module_fr_sfire_atm
- use module_model_constants, only: cp,xlv,g
- use module_fr_sfire_phys, only: fire_wind_height,have_fuel_cats,fcz0,fcwh,nfuelcats,no_fuel_cat,no_fuel_cat2,windrf
- use module_fr_sfire_util
- USE module_dm , only : wrf_dm_sum_reals
- use module_fr_sfire_phys, only: mfuelcats, nfuelcats ! for emissions
- USE module_state_description, only: num_tracer
- use module_fr_sfire_phys, only: fuel_name
- #ifdef WRF_CHEM
- USE module_state_description, only: num_chem
- USE module_configure, only: &
- p_co, &
- p_ch4, &
- p_h2, &
- p_no, &
- p_no2, &
- p_so2, &
- p_nh3, &
- p_p25, &
- p_p25i, &
- p_p25j, &
- p_oc1, &
- p_oc2, &
- p_bc1, &
- p_bc2, &
- p_ald, &
- p_csl, &
- p_eth, &
- p_hc3, &
- p_hc5, &
- p_hcho, &
- p_iso, &
- p_ket, &
- p_mgly, &
- p_ol2, &
- p_olt, &
- p_oli, &
- p_ora2,&
- p_tol, &
- p_xyl, &
- p_bigalk, &
- p_bigene, &
- p_c10h16, &
- p_c2h4, &
- p_c2h5oh, &
- p_c2h6, &
- p_c3h6, &
- p_c3h8, &
- p_ch3cooh, &
- p_ch3oh, &
- p_cres, &
- p_glyald, &
- ! p_hyac, &
- p_isopr, &
- p_macr, &
- p_mek, &
- p_mvk, &
- p_smoke ! tracer smoke exists only with CHEM
- #endif
- USE module_state_description, only: &
- p_tr17_1, &
- p_tr17_2, &
- p_tr17_3, &
- p_tr17_4, &
- p_tr17_5, &
- p_tr17_6, &
- p_tr17_7, &
- p_tr17_8
- implicit none
- #ifndef WRF_CHEM
- integer, parameter, private:: num_chem=0
- #endif
- logical, save :: have_wind_log_interpolation = .false. ! status
- ! emission tables
- REAL, dimension(mfuelcats), save:: &
- co=0., &
- ch4=0., &
- h2=0., &
- no=0., &
- no2=0., &
- so2=0., &
- nh3=0., &
- oc1=0., &
- oc2=0., &
- bc1=0., &
- bc2=0., &
- ald=0., &
- csl=0., &
- eth=0., &
- hc3=0., &
- p25=0., &
- p25i=0., &
- p25j=0., &
- hc5=0., &
- hcho=0., &
- iso=0., &
- ket=0., &
- mgly=0., &
- ol2=0., &
- olt=0., &
- oli=0., &
- ora2=0.,&
- tol=0., &
- xyl=0., &
- bigalk=0., &
- bigene=0., &
- c10h16=0., &
- c2h4=0., &
- c2h5oh=0., &
- c2h6=0., &
- c3h6=0., &
- c3h8=0., &
- ch3cooh=0., &
- ch3oh=0., &
- cres=0., &
- glyald=0., &
- ! hyac=0., &
- isopr=0., &
- macr=0., &
- mek=0., &
- mvk=0., &
- smoke=0., &
- tr17_1=0., &
- tr17_2=0., &
- tr17_3=0., &
- tr17_4=0., &
- tr17_5=0., &
- tr17_6=0., &
- tr17_7=0., &
- tr17_8=0.
- real, parameter:: & ! reciprocal molecular weights (mol/g)
- imw_co = 1./28.010,&
- imw_ch4 = 1./16.04,&
- imw_h2 = 1./2.016,&
- imw_no = 1./30.006,&
- imw_no2 = 1./46.006,&
- imw_so2 = 1./64.066,&
- imw_nh3 = 1./17.031
- real, parameter:: mw_air=28.97 ! molecular weight of air g/mol
- ! should be declared in the registry and stored in the state in future because of restarts
- real, pointer, save, dimension(:)::c_chem
- real, pointer, save, dimension(:)::c_fuel
- real, pointer, save, dimension(:)::c_tracer
- logical, save:: emis_read = .false.
- integer, save:: msglevel =1,printsums=0 ! when to print sums
- integer, parameter:: line=5 ! number of species, emissions, etc. per line
- contains
- ! chem arrays are chem tracer
- ! indices p_species are generated in inc/scalar_indices.inc and included in frame/module_configure.F
- subroutine read_emissions_table(chem_opt,tracer_opt)
- implicit none
- integer, intent(in)::chem_opt,tracer_opt
- logical, external:: wrf_dm_on_monitor
- external::wrf_dm_bcast_integer , wrf_dm_bcast_real
- integer, dimension(10)::compatible_chem_opt
- integer:: iounit,ierr,i
- character(len=128)::msg
- namelist/emissions/ compatible_chem_opt, printsums, &
- co, &
- ch4, &
- h2, &
- no, &
- no2, &
- so2, &
- nh3, &
- p25, &
- p25i, &
- p25j, &
- oc1, &
- oc2, &
- bc1, &
- bc2, &
- ald, &
- csl, &
- eth, &
- hc3, &
- hc5, &
- hcho, &
- iso, &
- ket, &
- mgly, &
- ol2, &
- olt, &
- oli, &
- ora2,&
- tol, &
- xyl, &
- bigalk, &
- bigene, &
- c10h16, &
- c2h4, &
- c2h5oh, &
- c2h6, &
- c3h6, &
- c3h8, &
- ch3cooh, &
- ch3oh, &
- cres, &
- glyald, &
- ! hyac, &
- isopr, &
- macr, &
- mek, &
- mvk, &
- smoke, &
- tr17_1, &
- tr17_2, &
- tr17_3, &
- tr17_4, &
- tr17_5, &
- tr17_6, &
- tr17_7, &
- tr17_8
- !$ if (OMP_GET_THREAD_NUM() .ne. 0)then
- !$ call crash('read_emissions_table: must be called from master thread')
- !$ endif
- IF ( wrf_dm_on_monitor() ) THEN
- ! we are the master task
- iounit=open_text_file('namelist.fire_emissions','read')
- compatible_chem_opt=0
- read(iounit,emissions,iostat=ierr)
- if(ierr.ne.0)call crash('read_emissions_table: error reading namelist emissions in file namelist.fire_emissions')
- CLOSE(iounit)
- write(msg,'(a,i3,a)')'reading emissions table for',nfuelcats,' fuel categories'
- call message(msg,level=0)
- if (.not.any(compatible_chem_opt.eq.chem_opt))then
- write(msg,'(a,i4,a)')'read_emissions_table: chem_opt=',chem_opt,' not between given compatible_chem_opt in namelist.fire_emissions'
- call message(msg,level=0)
- write(msg,'(a,10i4)')'compatible_chem_opt=', compatible_chem_opt
- call message(msg,level=0)
- call crash('chem_opt in namelist.input is not consistent with namelist.fire_emissions')
- endif
- ENDIF
- call wrf_dm_bcast_integer(printsums, 1)
- call wrf_dm_bcast_real(co, nfuelcats)
- call wrf_dm_bcast_real(ch4, nfuelcats)
- call wrf_dm_bcast_real(h2, nfuelcats)
- call wrf_dm_bcast_real(no, nfuelcats)
- call wrf_dm_bcast_real(no2, nfuelcats)
- call wrf_dm_bcast_real(so2, nfuelcats)
- call wrf_dm_bcast_real(nh3, nfuelcats)
- call wrf_dm_bcast_real(p25, nfuelcats)
- call wrf_dm_bcast_real(p25i, nfuelcats)
- call wrf_dm_bcast_real(p25j, nfuelcats)
- call wrf_dm_bcast_real(oc1, nfuelcats)
- call wrf_dm_bcast_real(oc2, nfuelcats)
- call wrf_dm_bcast_real(bc1, nfuelcats)
- call wrf_dm_bcast_real(bc2, nfuelcats)
- call wrf_dm_bcast_real(ald, nfuelcats)
- call wrf_dm_bcast_real(csl, nfuelcats)
- call wrf_dm_bcast_real(eth, nfuelcats)
- call wrf_dm_bcast_real(hc3, nfuelcats)
- call wrf_dm_bcast_real(hc5, nfuelcats)
- call wrf_dm_bcast_real(hcho, nfuelcats)
- call wrf_dm_bcast_real(iso, nfuelcats)
- call wrf_dm_bcast_real(ket, nfuelcats)
- call wrf_dm_bcast_real(mgly, nfuelcats)
- call wrf_dm_bcast_real(ol2, nfuelcats)
- call wrf_dm_bcast_real(olt, nfuelcats)
- call wrf_dm_bcast_real(oli, nfuelcats)
- call wrf_dm_bcast_real(ora2,nfuelcats)
- call wrf_dm_bcast_real(tol, nfuelcats)
- call wrf_dm_bcast_real(xyl, nfuelcats)
- call wrf_dm_bcast_real(bigalk, nfuelcats)
- call wrf_dm_bcast_real(bigene, nfuelcats)
- call wrf_dm_bcast_real(c10h16, nfuelcats)
- call wrf_dm_bcast_real(c2h4, nfuelcats)
- call wrf_dm_bcast_real(c2h5oh, nfuelcats)
- call wrf_dm_bcast_real(c2h6, nfuelcats)
- call wrf_dm_bcast_real(c3h6, nfuelcats)
- call wrf_dm_bcast_real(c3h8, nfuelcats)
- call wrf_dm_bcast_real(ch3cooh, nfuelcats)
- call wrf_dm_bcast_real(ch3oh, nfuelcats)
- call wrf_dm_bcast_real(cres, nfuelcats)
- call wrf_dm_bcast_real(glyald, nfuelcats)
- ! call wrf_dm_bcast_real(hyac, nfuelcats)
- call wrf_dm_bcast_real(isopr, nfuelcats)
- call wrf_dm_bcast_real(macr, nfuelcats)
- call wrf_dm_bcast_real(mek, nfuelcats)
- call wrf_dm_bcast_real(mvk, nfuelcats)
- call wrf_dm_bcast_real(smoke, nfuelcats)
- call wrf_dm_bcast_real(tr17_1, nfuelcats)
- call wrf_dm_bcast_real(tr17_2, nfuelcats)
- call wrf_dm_bcast_real(tr17_3, nfuelcats)
- call wrf_dm_bcast_real(tr17_4, nfuelcats)
- call wrf_dm_bcast_real(tr17_5, nfuelcats)
- call wrf_dm_bcast_real(tr17_6, nfuelcats)
- call wrf_dm_bcast_real(tr17_7, nfuelcats)
- call wrf_dm_bcast_real(tr17_8, nfuelcats)
- if(fire_print_msg .ge. msglevel .and.printsums .gt. 0)then
- ! should be stored in the registry future because of restarts
- write(msg,'(3(a,i3,1x))')'allocating c_chem size',num_chem,'c_tracer size',num_tracer,'c_fuel size',nfuelcats
- call message(msg,level=2)
- if(num_chem>0)then
- allocate(c_chem(num_chem))
- c_chem=0. ! cumulative burnt
- endif
- if(num_tracer>0)then
- allocate(c_tracer(num_tracer))
- c_tracer=0. ! cumulative burnt
- endif
- allocate(c_fuel(nfuelcats))
- c_fuel=0. ! total per timestep, rate burnt, cumulative burnt
- write(msg,'(a,i3,a,i3)')'allocated c_chem size',size(c_chem),' c_fuel size',size(c_fuel)
- call message(msg,level=2)
- endif
- emis_read=.true.
- end subroutine read_emissions_table
- subroutine add_fire_emissions(chem_opt,tracer_opt,dt,dx,dy, &
- ifms,ifme,jfms,jfme, &
- ifts,ifte,jtfs,jfte, &
- ids,ide,kds,kde,jds,jde, &
- ims,ime,kms,kme,jms,jme, &
- its,ite,kts,kte,jts,jte, &
- rho,dz8w, & ! input on atmosphere mesh
- fgip, fuel_frac_burnt, nfuel_cat, & ! input on fire mesh
- chem,tracer) ! output
- implicit none
- !*** purpose
- ! average fire emissions from fire mesh to coarser atmosphereic mesh and add to chemistry arrays
- !*** arguments
- ! the dimensions are in cells, not nodes!
- ! input
- integer, intent(in)::chem_opt,tracer_opt
- real, intent(in):: dt,dx,dy ! time step & mesh spacing
- integer, intent(in)::its,ite,kts,kte,jts,jte,ims,ime,kms,kme,jms,jme &! atm grid dims
- ,ids,ide,kds,kde,jds,jde
- integer, intent(in)::ifts,ifte,jtfs,jfte,ifms,ifme,jfms,jfme ! fire grid dims
- real, intent(in)::rho(ims:ime,kms:kme,jms:jme), & ! air density kg/m^3
- dz8w(ims:ime,kms:kme,jms:jme) ! layer height
- real, intent(in), dimension(ifms:ifme,jfms:jfme):: fgip, & ! initial fuel load kg/m^2
- fuel_frac_burnt, & ! fuel fraction burned this step kg/kg
- nfuel_cat ! fuel category (Anderson= 1 to 13)
- ! update
- real, intent(inout)::chem(ims:ime,kms:kme,jms:jme,num_chem),tracer(ims:ime,kms:kme,jms:jme,num_tracer)
- !*** local
- integer:: i,i_f,j,j_f,ir,jr,isz1,isz2,jsz1,jsz2,ioff,joff,ibase,jbase,cat,k1,areaw,m,k,k_p,errors
- real::fuel_burnt,vol,air,conv, avgw,emis
- character(len=128)msg
- real, dimension(mfuelcats)::s_fuel,t_fuel,r_fuel ! total per timestep, rate burnt
- #ifdef WRF_CHEM
- real, dimension(num_chem) ::s_chem,t_chem,r_chem ! total per timestep, rate burnt
- real, dimension(num_chem) ::a_chem,g_chem ! concentration in ground level 1
- integer, parameter:: chem_np=46
- integer:: chem_pointers(chem_np)
- character(len=8)::chem_names(chem_np)
- #endif
- real, dimension(num_tracer) ::s_tracer,t_tracer,r_tracer ! total per timestep, rate burnt
- integer, parameter:: tracer_np=8
- integer:: tracer_pointers(tracer_np)
- character(len=8)::tracer_names(tracer_np)
- !*** executable
- !check mesh dimensions and domain dimensions
- call check_mesh_2dim(its,ite,jts,jte,ims,ime,jms,jme)
- call check_mesh_2dim(ifts,ifte,jtfs,jfte,ifms,ifme,jfms,jfme)
- if(.not.emis_read)call crash('add_fire_emissions: read_emissions_table must be called first')
- write(msg,'(a,i3,a,i3,a,i3)')'add_fire_emissions: chem_opt=',chem_opt,' species ',num_chem,' tracers',num_tracer
- call message(msg)
- #ifdef WRF_CHEM
- chem_pointers= (/ &
- p_co, &
- p_ch4, &
- p_h2, &
- p_no, &
- p_no2, &
- p_so2, &
- p_nh3, &
- p_p25, &
- p_p25i, &
- p_p25j, &
- p_oc1, &
- p_oc2, &
- p_bc1, &
- p_bc2, &
- p_ald, &
- p_csl, &
- p_eth, &
- p_hc3, &
- p_hc5, &
- p_hcho, &
- p_iso, &
- p_ket, &
- p_mgly, &
- p_ol2, &
- p_olt, &
- p_oli, &
- p_ora2,&
- p_tol, &
- p_xyl, &
- p_bigalk, &
- p_bigene, &
- p_c10h16, &
- p_c2h4, &
- p_c2h5oh, &
- p_c2h6, &
- p_c3h6, &
- p_c3h8, &
- p_ch3cooh, &
- p_ch3oh, &
- p_cres, &
- p_glyald, &
- ! p_hyac, &
- p_isopr, &
- p_macr, &
- p_mek, &
- p_mvk, &
- p_smoke /)
- chem_names= (/ &
- 'co ', &
- 'ch4 ', &
- 'h2 ', &
- 'no ', &
- 'no2 ', &
- 'so2 ', &
- 'nh3 ', &
- 'p25 ', &
- 'p25i ', &
- 'p25j ', &
- 'oc1 ', &
- 'oc2 ', &
- 'bc1 ', &
- 'bc2 ', &
- 'ald ', &
- 'csl ', &
- 'eth ', &
- 'hc3 ', &
- 'hc5 ', &
- 'hcho ', &
- 'iso ', &
- 'ket ', &
- 'mgly ', &
- 'ol2 ', &
- 'olt ', &
- 'oli ', &
- 'ora2 ',&
- 'tol ', &
- 'xyl ', &
- 'bigalk ', &
- 'bigene ', &
- 'c10h16 ', &
- 'c2h4 ', &
- 'c2h5oh ', &
- 'c2h6 ', &
- 'c3h6 ', &
- 'c3h8 ', &
- 'ch3cooh ', &
- 'ch3oh ', &
- 'cres ', &
- 'glyald ', &
- ! 'hyac ', &
- 'isopr ', &
- 'macr ', &
- 'mek ', &
- 'mvk ', &
- 'smoke ' /)
- call check_pointers('chem',chem,chem_names,chem_pointers)
- #endif
- tracer_pointers= (/ &
- p_tr17_1, &
- p_tr17_2, &
- p_tr17_3, &
- p_tr17_4, &
- p_tr17_5, &
- p_tr17_6, &
- p_tr17_7, &
- p_tr17_8 /)
- tracer_names= (/ &
- 'tr17_1 ', &
- 'tr17_2 ', &
- 'tr17_3 ', &
- 'tr17_4 ', &
- 'tr17_5 ', &
- 'tr17_6 ', &
- 'tr17_7 ', &
- 'tr17_8 ' /)
- call check_pointers('tracer',tracer,tracer_names,tracer_pointers)
- ! compute mesh sizes
- isz1 = ite-its+1
- jsz1 = jte-jts+1
- isz2 = ifte-ifts+1
- jsz2 = jfte-jtfs+1
- ! check mesh sizes
- if(isz1.le.0.or.jsz1.le.0.or.isz2.le.0.or.jsz2.le.0)then
- call message('all mesh sizes must be positive',level=0)
- goto 9
- endif
- ! compute mesh ratios
- ir=isz2/isz1
- jr=jsz2/jsz1
- if(isz2.ne.isz1*ir .or. jsz2.ne.jsz1*jr)then
- call message('input mesh size must be multiple of output mesh size',level=0)
- goto 9
- endif
- avgw = 1.0/(ir*jr) ! averaging weight = 1/number of fire cells per atm cell
- ! initialize emissions statistics per timestep
- #ifdef WRF_CHEM
- t_chem = 0.
- #endif
- t_fuel = 0.
- do i=1,num_tracer
- t_tracer(i)=0.
- enddo
- ! conversion fuel_burnt kg/m^2 -> chem_X 1e6*mol/mol
- ! emis_X(g/kg)*fuel_burnt(kg/m^2)/mw_X(g/mol) = emissions (mol/m^2)
- ! rho(kg/m^3)*dz8w(m))/mw_air(28.97e-3 kg/mol) = dry air in the 1st layer (mol/m^2)
- k1 = kts
- write(msg,'(a,i3)')'Fire emissions inserted into atmosphere level',k1
- call message(msg,level=msglevel)
- #ifdef WRF_CHEM
- if(fire_print_msg .ge. msglevel .and.printsums .gt. 0)then
- ! sum ground concentrations and check for nans
- a_chem=0.0
- g_chem=0.0
- errors=0
- do j=jts,jte
- do k=1,chem_np
- k_p=chem_pointers(k)
- do i=its,ite
- if(chem(i,k1,j,k_p ).ne.chem(i,k1,j,k_p ))errors=errors+1
- a_chem(k_p)=a_chem(k_p)+chem(i,k1,j,k_p )
- enddo
- enddo
- enddo
- if(errors>0)call crash('NaN before chem update')
- call wrf_dm_sum_reals(a_chem,g_chem)
- call message('Layer1 raw sums before adding fire emissions',level=msglevel)
- do i=1,chem_np,line
- m=min(i+line-1,chem_np)
- write(msg,80)'Emissions ',(trim(chem_names(j)),j=i,m)
- call message(msg,level=msglevel)
- write(msg,81)'Layer1 beg',(g_chem(chem_pointers(j)),j=i,m)
- call message(msg,level=msglevel)
- call message(' ',level=msglevel)
- enddo
- endif
- #endif
- do j=max(jds+1,jts),min(jte,jde-1) ! safe distance from domain boundary
- jbase=jtfs+jr*(j-jts) ! indexing
- do i=max(ids+1,its),min(ite,ide-1)
- ibase=ifts+ir*(i-its) ! indexing
- !air = 1e6*mw_air/rho(i,kds,j) ! 1e6*mw_air/air density
- !vol = avgw/dz8w(i,kds,j) ! averaging volume factor / 1st layer depth
- do joff=0,jr-1
- j_f=joff+jbase
- do ioff=0,ir-1
- i_f=ioff+ibase
- !*** fire cell (i_f,j_f) contributes to atmosphere cell (i,j) at ground level
- fuel_burnt = fgip(i_f,j_f) * fuel_frac_burnt(i_f,j_f) ! kg/m^2
- cat = nfuel_cat(i_f, j_f) ! usually 1 to 13
- if(cat.lt.no_fuel_cat)t_fuel(cat)=t_fuel(cat) + fuel_burnt
- !*** chem compounds emissions given in g/kg
- ! fuel_burnt kg/m^2 * table g/kg -> ppmv = 1e6*mol/mol in 1st layer
- !AK rho is in kg/m3 so the conversion factor must be scaled by a 1000 to match
- !emissions in grams
- ! conv = avgw*1e6*mw_air/(rho(i,k1,j)*dz8w(i,k1,j))
- conv = avgw*1e3*mw_air/(rho(i,k1,j)*dz8w(i,k1,j))
- #ifdef WRF_CHEM
- emis=co (cat)*fuel_burnt ! emission from fire cell in g/m^2
- t_chem(p_co) = t_chem(p_co) + emis ! add to total
- ! if(isnan(chem(i,k1,j,p_co )))call crash('NaN before')
- chem(i,k1,j,p_co )=chem(i,k1,j,p_co ) + emis*conv*imw_co ! add to chem
- ! if(isnan(chem(i,k1,j,p_co )))call crash('NaN after')
- emis=ch4 (cat)*fuel_burnt ! emission from fire cell in g/m^2
- t_chem(p_ch4) = t_chem(p_ch4) + emis ! add to total
- chem(i,k1,j,p_ch4 )=chem(i,k1,j,p_ch4 ) + emis*conv*imw_ch4 ! add to chem
- emis=h2 (cat)*fuel_burnt ! emission from fire cell in g/m^2
- t_chem(p_h2) = t_chem(p_h2) + emis ! add to total
- chem(i,k1,j,p_h2 )=chem(i,k1,j,p_h2 ) + emis*conv*imw_h2 ! add to chem
- emis=no (cat)*fuel_burnt ! emission from fire cell in g/m^2
- t_chem(p_no) = t_chem(p_no) + emis ! add to total
- chem(i,k1,j,p_no )=chem(i,k1,j,p_no ) + emis*conv*imw_no ! add to chem
- emis=no2 (cat)*fuel_burnt ! emission from fire cell in g/m^2
- t_chem(p_no2) = t_chem(p_no2) + emis ! add to total
- chem(i,k1,j,p_no2 )=chem(i,k1,j,p_no2 ) + emis*conv*imw_no2 ! add to chem
- emis=so2 (cat)*fuel_burnt ! emission from fire cell in g/m^2
- t_chem(p_so2) = t_chem(p_so2) + emis ! ads to total
- chem(i,k1,j,p_so2 )=chem(i,k1,j,p_so2 ) + emis*conv*imw_so2 ! add to chem
- emis=nh3 (cat)*fuel_burnt ! emission from fire cell in g/m^2
- t_chem(p_nh3) = t_chem(p_nh3) + emis ! add to total
- chem(i,k1,j,p_nh3 )=chem(i,k1,j,p_nh3 ) + emis*conv*imw_nh3 ! add to chem
- !*** other emissions already given in mol/kg
- ! fuel_burnt kg/m^2 * table mol/kg -> ppmv = 1e6*mol/mol in 1st layer dry air in 1st layer
-
- ! same conversion factor but we will not divide by the molecular weight of the compound
- ! conv = avgw*1e6*mw_air/(rho(i,kds,j)*dz8w(i,kds,j))
- emis=ald (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_ald) = t_chem(p_ald) + emis ! add to total
- chem(i,k1,j,p_ald )=chem(i,k1,j,p_ald )+emis*conv
- emis=csl (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_csl) = t_chem(p_csl) + emis ! add to total
- chem(i,k1,j,p_csl )=chem(i,k1,j,p_csl )+emis*conv
- emis=eth (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_eth) = t_chem(p_eth) + emis ! add to total
- chem(i,k1,j,p_eth )=chem(i,k1,j,p_eth )+emis*conv
- emis=hc3 (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_hc3) = t_chem(p_hc3) + emis ! add to total
- chem(i,k1,j,p_hc3 )=chem(i,k1,j,p_hc3 )+emis*conv
- emis=hc5 (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_hc5) = t_chem(p_hc5) + emis ! add to total
- chem(i,k1,j,p_hc5 )=chem(i,k1,j,p_hc5 )+emis*conv
- emis=hcho (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_hcho) = t_chem(p_hcho) + emis ! add to total
- chem(i,k1,j,p_hcho)=chem(i,k1,j,p_hcho)+emis*conv
- emis=iso (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_iso) = t_chem(p_iso) + emis ! add to total
- chem(i,k1,j,p_iso )=chem(i,k1,j,p_iso )+emis*conv
- emis=ket (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_ket) = t_chem(p_ket) + emis ! add to total
- chem(i,k1,j,p_ket )=chem(i,k1,j,p_ket )+emis*conv
- emis=mgly (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_mgly) = t_chem(p_mgly) + emis ! add to total
- chem(i,k1,j,p_mgly)=chem(i,k1,j,p_mgly)+emis*conv
- emis=ol2 (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_ol2) = t_chem(p_ol2) + emis ! add to total
- chem(i,k1,j,p_ol2 )=chem(i,k1,j,p_ol2 )+emis*conv
- emis=olt (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_olt) = t_chem(p_olt) + emis ! add to total
- chem(i,k1,j,p_olt )=chem(i,k1,j,p_olt )+emis*conv
- emis=oli (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_oli) = t_chem(p_oli) + emis ! add to total
- chem(i,k1,j,p_oli )=chem(i,k1,j,p_oli )+emis*conv
- emis=ora2 (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_ora2) = t_chem(p_ora2) + emis ! add to total
- chem(i,k1,j,p_ora2)=chem(i,k1,j,p_ora2)+emis*conv
- emis=tol (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_tol) = t_chem(p_tol) + emis ! add to total
- chem(i,k1,j,p_tol )=chem(i,k1,j,p_tol )+emis*conv
- emis=xyl (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_xyl) = t_chem(p_xyl) + emis ! add to total
- chem(i,k1,j,p_xyl )=chem(i,k1,j,p_xyl )+emis*conv
- emis=bigalk (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_bigalk) = t_chem(p_bigalk) + emis ! add to total
- chem(i,k1,j,p_bigalk )=chem(i,k1,j,p_bigalk )+emis*conv
- emis=bigene (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_bigene) = t_chem(p_bigene) + emis ! add to total
- chem(i,k1,j,p_bigene )=chem(i,k1,j,p_bigene )+emis*conv
- emis=c10h16 (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_c10h16) = t_chem(p_c10h16) + emis ! add to total
- chem(i,k1,j,p_c10h16 )=chem(i,k1,j,p_c10h16 )+emis*conv
- emis=c2h4 (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_c2h4) = t_chem(p_c2h4) + emis ! add to total
- chem(i,k1,j,p_c2h4 )=chem(i,k1,j,p_c2h4 )+emis*conv
- emis=c2h5oh (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_c2h5oh) = t_chem(p_c2h5oh) + emis ! add to total
- chem(i,k1,j,p_c2h5oh )=chem(i,k1,j,p_c2h5oh )+emis*conv
- emis=c2h6 (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_c2h6) = t_chem(p_c2h6) + emis ! add to total
- chem(i,k1,j,p_c2h6 )=chem(i,k1,j,p_c2h6 )+emis*conv
- emis=c3h6 (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_c3h6) = t_chem(p_c3h6) + emis ! add to total
- chem(i,k1,j,p_c3h6 )=chem(i,k1,j,p_c3h6 )+emis*conv
- emis=c3h8 (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_c3h8) = t_chem(p_c3h8) + emis ! add to total
- chem(i,k1,j,p_c3h8 )=chem(i,k1,j,p_c3h8 )+emis*conv
- emis=ch3cooh (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_ch3cooh) = t_chem(p_ch3cooh) + emis ! add to total
- chem(i,k1,j,p_ch3cooh )=chem(i,k1,j,p_ch3cooh )+emis*conv
- emis=ch3oh (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_ch3oh) = t_chem(p_ch3oh) + emis ! add to total
- chem(i,k1,j,p_ch3oh )=chem(i,k1,j,p_ch3oh )+emis*conv
- emis=cres (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_cres) = t_chem(p_cres) + emis ! add to total
- chem(i,k1,j,p_cres )=chem(i,k1,j,p_cres )+emis*conv
- emis=glyald (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_glyald) = t_chem(p_glyald) + emis ! add to total
- chem(i,k1,j,p_glyald )=chem(i,k1,j,p_glyald )+emis*conv
- emis=isopr (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_isopr) = t_chem(p_isopr) + emis ! add to total
- chem(i,k1,j,p_isopr )=chem(i,k1,j,p_isopr )+emis*conv
- emis=macr (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_macr) = t_chem(p_macr) + emis ! add to total
- chem(i,k1,j,p_macr )=chem(i,k1,j,p_macr )+emis*conv
- emis=mek (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_mek) = t_chem(p_mek) + emis ! add to total
- chem(i,k1,j,p_mek )=chem(i,k1,j,p_mek )+emis*conv
- emis=mvk (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_mvk) = t_chem(p_mvk) + emis ! add to total
- chem(i,k1,j,p_mvk )=chem(i,k1,j,p_mvk )+emis*conv
- ! aerosols
- ! fuel_burnt kg/m^2 * table g/kg -> ug/kg dry air in 1st layer
- ! see also chem/emissions_driver.F line 515
- conv = avgw*1e6/(rho(i,k1,j)*dz8w(i,k1,j))
- emis=p25 (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_p25) = t_chem(p_p25) + emis ! add to total
- chem(i,k1,j,p_p25 )=chem(i,k1,j,p_p25 )+emis*conv
- emis=p25i (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_p25i) = t_chem(p_p25i) + emis ! add to total
- chem(i,k1,j,p_p25i )=chem(i,k1,j,p_p25i )+emis*conv
- emis=p25j (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_p25j) = t_chem(p_p25j) + emis ! add to total
- chem(i,k1,j,p_p25j )=chem(i,k1,j,p_p25j )+emis*conv
- emis=oc1 (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_oc1 ) = t_chem(p_oc1 ) + emis ! add to total
- chem(i,k1,j,p_oc1 )=chem(i,k1,j,p_oc1 )+emis*conv
- emis=oc2 (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_oc2 ) = t_chem(p_oc2 ) + emis ! add to total
- chem(i,k1,j,p_oc2 )=chem(i,k1,j,p_oc2 )+emis*conv
- emis=bc1 (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_bc1 ) = t_chem(p_bc1 ) + emis ! add to total
- chem(i,k1,j,p_bc1 )=chem(i,k1,j,p_bc1 )+emis*conv
- emis=bc2 (cat)*fuel_burnt ! emission from fire cell
- t_chem(p_bc2 ) = t_chem(p_bc2 ) + emis ! add to total
- chem(i,k1,j,p_bc2 )=chem(i,k1,j,p_bc2 )+emis*conv
- #endif
- if (num_tracer >0)then
- ! treat tracers exactly the same as aerosols, emissions g/kg fuel burned, tracer concentration ug/kg dry air
- conv = avgw*1e6/(rho(i,k1,j)*dz8w(i,k1,j))
- emis=tr17_1 (cat)*fuel_burnt ! emission from fire cell
- t_tracer(p_tr17_1) = t_tracer(p_tr17_1) + emis ! add to total
- tracer(i,k1,j,p_tr17_1 )=tracer(i,k1,j,p_tr17_1 )+emis*conv
- emis=tr17_2 (cat)*fuel_burnt ! emission from fire cell
- t_tracer(p_tr17_2) = t_tracer(p_tr17_2) + emis ! add to total
- tracer(i,k1,j,p_tr17_2 )=tracer(i,k1,j,p_tr17_2 )+emis*conv
- emis=tr17_3 (cat)*fuel_burnt ! emission from fire cell
- t_tracer(p_tr17_3) = t_tracer(p_tr17_3) + emis ! add to total
- tracer(i,k1,j,p_tr17_3 )=tracer(i,k1,j,p_tr17_3 )+emis*conv
- emis=tr17_4 (cat)*fuel_burnt ! emission from fire cell
- t_tracer(p_tr17_4) = t_tracer(p_tr17_4) + emis ! add to total
- tracer(i,k1,j,p_tr17_4 )=tracer(i,k1,j,p_tr17_4 )+emis*conv
- emis=tr17_5 (cat)*fuel_burnt ! emission from fire cell
- t_tracer(p_tr17_5) = t_tracer(p_tr17_5) + emis ! add to total
- tracer(i,k1,j,p_tr17_5 )=tracer(i,k1,j,p_tr17_5 )+emis*conv
- emis=tr17_6 (cat)*fuel_burnt ! emission from fire cell
- t_tracer(p_tr17_6) = t_tracer(p_tr17_6) + emis ! add to total
- tracer(i,k1,j,p_tr17_6 )=tracer(i,k1,j,p_tr17_6 )+emis*conv
- emis=tr17_7 (cat)*fuel_burnt ! emission from fire cell
- t_tracer(p_tr17_7) = t_tracer(p_tr17_7) + emis ! add to total
- tracer(i,k1,j,p_tr17_7 )=tracer(i,k1,j,p_tr17_7 )+emis*conv
- emis=tr17_8 (cat)*fuel_burnt ! emission from fire cell
- t_tracer(p_tr17_8) = t_tracer(p_tr17_8) + emis ! add to total
- tracer(i,k1,j,p_tr17_8 )=tracer(i,k1,j,p_tr17_8 )+emis*conv
- endif
- enddo
- enddo
- enddo
- enddo
- #ifdef WRF_CHEM
- if(fire_print_msg .ge. msglevel .and.printsums .gt. 0)then
- ! sum ground concentrations and check for nans
- a_chem=0.0
- g_chem=0.0
- errors=0
- do j=jts,jte
- do k=1,chem_np
- k_p=chem_pointers(k)
- do i=its,ite
- if(chem(i,k1,j,k_p ).ne.chem(i,k1,j,k_p ))errors=errors+1
- a_chem(k_p)=a_chem(k_p)+chem(i,k1,j,k_p )
- enddo
- enddo
- enddo
- if(errors>0)call crash('NaN after chem update')
- call wrf_dm_sum_reals(a_chem,g_chem)
- call message('Layer1 raw sums after adding fire emissions',level=msglevel)
- do i=1,chem_np,line
- m=min(i+line-1,chem_np)
- write(msg,80)'Emissions ',(trim(chem_names(j)),j=i,m)
- call message(msg,level=msglevel)
- write(msg,81)'Layer1 end',(g_chem(chem_pointers(j)),j=i,m)
- call message(msg,level=msglevel)
- call message(' ',level=msglevel)
- enddo
- endif
- #endif
- if(fire_print_msg .ge. msglevel .and.printsums .gt. 0)then
- ! sum over processes and add to cumulative sums
- ! fuel burned
- call wrf_dm_sum_reals(t_fuel,s_fuel)
- ! scale
- s_fuel = s_fuel*dx*dy
- ! get rates
- r_fuel = s_fuel/dt
- ! add to cumulative totals
- if(size(c_fuel).ne.nfuelcats)call crash('add_fire_emissions: bad size c_fuel')
- c_fuel = c_fuel + s_fuel
- #ifdef WRF_CHEM
- call wrf_dm_sum_reals(a_chem,g_chem)
- ! chem
- call wrf_dm_sum_reals(t_chem,s_chem)
- s_chem = s_chem*dx*dy
- r_chem = s_chem/dt
- if(size(c_chem).ne.num_chem)call crash('add_fire_emissions: bad size c_chem')
- c_chem = c_chem + s_chem
- call message('Total emissions in g or mol per the table',level=msglevel)
- do i=1,chem_np,line
- m=min(i+line-1,chem_np)
- write(msg,80)'Emissions ',(trim(chem_names(j)),j=i,m)
- call message(msg,level=msglevel)
- write(msg,81)'Cumulative',(c_chem(chem_pointers(j)),j=i,m)
- call message(msg,level=msglevel)
- write(msg,81)'Rate per s',(r_chem(chem_pointers(j)),j=i,m)
- call message(msg,level=msglevel)
- call message(' ',level=msglevel)
- enddo
- #endif
- if(num_tracer >0)then
- ! tracer
- call wrf_dm_sum_reals(t_tracer,s_tracer)
- s_tracer = s_tracer*dx*dy
- r_tracer = s_tracer/dt
- if(size(c_tracer).ne.num_tracer)call crash('add_fire_emissions: bad size c_tracer')
- c_tracer = c_tracer + s_tracer
-
- call message('Total emissions in g',level=msglevel)
- do i=1,tracer_np,line
- m=min(i+line-1,tracer_np)
- write(msg,80)'Emissions ',(trim(tracer_names(j)),j=i,m)
- call message(msg,level=msglevel)
- write(msg,81)'Cumulative',(c_tracer(tracer_pointers(j)),j=i,m)
- call message(msg,level=msglevel)
- write(msg,81)'Rate per s',(r_tracer(tracer_pointers(j)),j=i,m)
- call message(msg,level=msglevel)
- call message(' ',level=msglevel)
- enddo
- endif
- do cat=1,nfuelcats
- if(c_fuel(cat) > 0.)then
- write(msg,83)fuel_name(cat),' burned',c_fuel(cat),'kg, rate',r_fuel(cat),' kg/s'
- call message(msg,level=msglevel)
- endif
- enddo
- write(msg,83)'Total fuel',' burned',sum(c_fuel),'kg, rate',sum(r_fuel),' kg/s'
- call message(msg,level=msglevel)
- endif
- 80 format(a,8a11)
- 81 format(a,8e11.3)
- 83 format(a30,a,g14.4,a,g14.4,a,a)
- return
- 9 continue
- !$OMP CRITICAL(SFIRE_ATM_CRIT)
- write(msg,91)ifts,ifte,jtfs,jfte,ifms,ifme,jfms,jfme
- call message(msg,level=0)
- write(msg,91)its,ite,jts,jte,ims,ime,jms,jme
- call message(msg,level=0)
- write(msg,92)'input mesh size:',isz2,jsz2
- call message(msg,level=0)
- 91 format('dimensions: ',8i8)
- write(msg,92)'output mesh size:',isz1,jsz1
- call message(msg,level=0)
- 92 format(a,2i8)
- !$OMP END CRITICAL(SFIRE_ATM_CRIT)
- call crash('add_fire_emissions: bad mesh sizes')
- end subroutine add_fire_emissions
- !
- !***
- !
- subroutine check_pointers(array_name,array,pointer_names,pointers)
- implicit none
- !*** arguments
- character(len=*)::array_name
- real, dimension(:,:,:,:)::array
- character(len=*), dimension(:)::pointer_names
- integer, dimension(:)::pointers
- !*** local
- integer::np,i,m,j
- character(len=256)::msg
- !** executable
- np=ubound(pointers,1)
- 993 format(3a,4(1x,i3,':',i3))
- !$OMP CRITICAL(SFIRE_ATM_CRIT)
- write(msg,993)'array ',array_name,' has dimensions ',&
- lbound(array,1),ubound(array,1), &
- lbound(array,2),ubound(array,2), &
- lbound(array,3),ubound(array,3), &
- lbound(array,4),ubound(array,4)
- call message(msg)
- do i=1,np,line
- m=min(i+line-1,np)
- write(msg,'(a,8(1x,a8))')'Species',(trim(pointer_names(j)),j=i,m)
- call message(msg,level=msglevel)
- write(msg,'(a,8i9)') 'Pointer',(pointers(j),j=i,m)
- call message(msg,level=msglevel)
- call message(' ')
- enddo
- if(maxval(pointers) > ubound(array,4) .or. minval(pointers) < lbound(array,4))then
- write(msg,'(3a)')'add_fire_emissions: a ',array_name,' pointer is out of bounds'
- call crash(msg)
- endif
- !$OMP END CRITICAL(SFIRE_ATM_CRIT)
- end subroutine check_pointers
- !
- !***
- !
- SUBROUTINE fire_tendency( &
- ids,ide, kds,kde, jds,jde, & ! dimensions
- ims,ime, kms,kme, jms,jme, &
- its,ite, kts,kte, jts,jte, &
- grnhfx,grnqfx,canhfx,canqfx, & ! heat fluxes summed up to atm grid
- alfg,alfc,z1can, & ! coeffients, properties, geometry
- zs,z_at_w,dz8w,mu,rho, &
- rthfrten,rqvfrten) ! theta and Qv tendencies
- ! This routine is atmospheric physics
- ! it does NOT go into module_fr_sfire_phys because it is not fire physics
- ! taken from the code by Ned Patton, only order of arguments change to the convention here
- ! --- this routine takes fire generated heat and moisture fluxes and
- ! calculates their influence on the theta and water vapor
- ! --- note that these tendencies are valid at the Arakawa-A location
- IMPLICIT NONE
- ! --- incoming variables
- INTEGER , INTENT(in) :: ids,ide, kds,kde, jds,jde, &
- ims,ime, kms,kme, jms,jme, &
- its,ite, kts,kte, jts,jte
- REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: grnhfx,grnqfx ! W/m^2
- REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: canhfx,canqfx ! W/m^2
- REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: zs ! topography (m abv sealvl)
- REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: mu ! dry air mass (Pa)
- REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: z_at_w ! m abv sealvl
- REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: dz8w ! dz across w-lvl
- REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: rho ! density
- REAL, INTENT(in) :: alfg ! extinction depth ground fire heat (m)
- REAL, INTENT(in) :: alfc ! extinction depth crown fire heat (m)
- REAL, INTENT(in) :: z1can ! height of crown fire heat release (m)
- ! --- outgoing variables
- REAL, INTENT(out), DIMENSION( ims:ime,kms:kme,jms:jme ) :: &
- rthfrten, & ! theta tendency from fire (in mass units)
- rqvfrten ! Qv tendency from fire (in mass units)
- ! --- local variables
- INTEGER :: i,j,k
- INTEGER :: i_st,i_en, j_st,j_en, k_st,k_en
- REAL :: cp_i
- REAL :: rho_i
- REAL :: xlv_i
- REAL :: z_w
- REAL :: fact_g, fact_c
- REAL :: alfg_i, alfc_i
- REAL, DIMENSION( its:ite,kts:kte,jts:jte ) :: hfx,qfx
-
- !! character(len=128)::msg
- do j=jts,jte
- do k=kts,min(kte+1,kde)
- do i=its,ite
- rthfrten(i,k,j)=0.
- rqvfrten(i,k,j)=0.
- enddo
- enddo
- enddo
- ! --- set some local constants
-
- cp_i = 1./cp ! inverse of specific heat
- xlv_i = 1./xlv ! inverse of latent heat
- alfg_i = 1./alfg
- alfc_i = 1./alfc
- !!write(msg,'(8e11.3)')cp,cp_i,xlv,xlv_i,alfg,alfc,z1can
- !!call message(msg)
- call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnhfx,'fire_tendency:grnhfx')
- call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnqfx,'fire_tendency:grnqfx')
- ! --- set loop indicies : note that
- i_st = MAX(its,ids+1)
- i_en = MIN(ite,ide-1)
- k_st = kts
- k_en = MIN(kte,kde-1)
- j_st = MAX(jts,jds+1)
- j_en = MIN(jte,jde-1)
- ! --- distribute fluxes
- DO j = j_st,j_en
- DO k = k_st,k_en
- DO i = i_st,i_en
- ! --- set z (in meters above ground)
- z_w = z_at_w(i,k,j) - zs(i,j) ! should be zero when k=k_st
- ! --- heat flux
- fact_g = cp_i * EXP( - alfg_i * z_w )
- IF ( z_w < z1can ) THEN
- fact_c = cp_i
- ELSE
- fact_c = cp_i * EXP( - alfc_i * (z_w - z1can) )
- END IF
- hfx(i,k,j) = fact_g * grnhfx(i,j) + fact_c * canhfx(i,j)
- !! write(msg,2)i,k,j,z_w,grnhfx(i,j),hfx(i,k,j)
- !!2 format('hfx:',3i4,6e11.3)
- !! call message(msg)
- ! --- vapor flux
- fact_g = xlv_i * EXP( - alfg_i * z_w )
- IF (z_w < z1can) THEN
- fact_c = xlv_i
- ELSE
- fact_c = xlv_i * EXP( - alfc_i * (z_w - z1can) )
- END IF
- qfx(i,k,j) = fact_g * grnqfx(i,j) + fact_c * canqfx(i,j)
-
- !! if(hfx(i,k,j).ne.0. .or. qfx(i,k,j) .ne. 0.)then
- !! write(msg,1)i,k,j,hfx(i,k,j),qfx(i,k,j)
- !!1 format('tend:',3i6,2e11.3)
- !! call message(msg)
- ! endif
- END DO
- END DO
- END DO
- ! --- add flux divergence to tendencies
- !
- ! multiply by dry air mass (mu) to eliminate the need to
- ! call sr. calculate_phy_tend (in dyn_em/module_em.F)
- ! print *,'fire_tendency:',i_st,i_en,j_st,j_en
- DO j = j_st,j_en
- DO k = k_st,k_en-1
- DO i = i_st,i_en
- rho_i = 1./rho(i,k,j)
- rthfrten(i,k,j) = - mu(i,j) * rho_i * (hfx(i,k+1,j)-hfx(i,k,j)) / dz8w(i,k,j)
- rqvfrten(i,k,j) = - mu(i,j) * rho_i * (qfx(i,k+1,j)-qfx(i,k,j)) / dz8w(i,k,j)
-
- ! print *,i,j,k,rthfrten(i,k,j)
- END DO
- END DO
- END DO
- call print_3d_stats(its,ite,kts,kte,jts,jte,ims,ime,kms,kme,jms,jme,rthfrten,'fire_tendency:rthfrten')
- call print_3d_stats(its,ite,kts,kte,jts,jte,ims,ime,kms,kme,jms,jme,rqvfrten,'fire_tendency:rqvfrten')
- RETURN
- END SUBROUTINE fire_tendency
- !
- !***
- !
- subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no output
- ids,ide, kds,kde, jds,jde, & ! atm grid dimensions
- ims,ime, kms,kme, jms,jme, &
- ips,ipe,jps,jpe, &
- its,ite,jts,jte, &
- ifds, ifde, jfds, jfde, & ! fire grid dimensions
- ifms, ifme, jfms, jfme, &
- ifps, ifpe, jfps, jfpe, & ! fire patch bounds
- ifts,ifte,jfts,jfte, &
- ir,jr, & ! atm/fire grid ratio
- u_frame, v_frame, & ! velocity frame correction
- u,v, & ! atm grid arrays in
- ph,phb, &
- z0,zs, &
- uah,vah, &
- uf,vf) ! fire grid arrays out
-
- implicit none
- ! Jan Mandel, October 2010
- !*** purpose: interpolate winds and height
- !*** arguments
- integer, intent(in)::id
- integer, intent(in):: &
- ids,ide, kds,kde, jds,jde, & ! atm domain bounds
- ims,ime, kms,kme, jms,jme, & ! atm memory bounds
- ips,ipe,jps,jpe, &
- its,ite,jts,jte, & ! atm tile bounds
- ifds, ifde, jfds, jfde, & ! fire domain bounds
- ifms, ifme, jfms, jfme, & ! fire memory bounds
- ifps, ifpe, jfps, jfpe, & ! fire patch bounds
- ifts,ifte,jfts,jfte, & ! fire tile bounds
- ir,jr ! atm/fire grid refinement ratio
- real, intent(in):: u_frame, v_frame ! velocity frame correction
- real,intent(in),dimension(ims:ime,kms:kme,jms:jme)::&
- u,v, & ! atm wind velocity, staggered
- ph,phb ! geopotential
- real,intent(in),dimension(ims:ime,jms:jme)::&
- z0, & ! roughness height
- zs ! terrain height
- real,intent(out),dimension(ims:ime,jms:jme)::&
- uah, & ! atm wind at fire_wind_height, diagnostics
- vah !
- real,intent(out), dimension(ifms:ifme,jfms:jfme)::&
- uf,vf ! wind velocity fire grid nodes
-
-
- !*** local
- character(len=256)::msg
- #define TDIMS its-2,ite+2,jts-2,jte+2
- real, dimension(its-2:ite+2,jts-2:jte+2):: ua,va ! atm winds, interpolated over height, still staggered grid
- real, dimension(its-2:ite+2,kds:kde,jts-2:jte+2):: altw,altub,altvb,hgtu,hgtv ! altitudes
- integer:: i,j,k,ifts1,ifte1,jfts1,jfte1,ite1,jte1
- integer::itst,itet,jtst,jtet,itsu,iteu,jtsu,jteu,itsv,itev,jtsv,jtev
- integer::kdmax,its1,jts1,ips1,jps1
- integer::itsou,iteou,jtsou,jteou,itsov,iteov,jtsov,jteov
- real:: ground,loght,loglast,logz0,logfwh,ht,zr
- real::r_nan
- integer::i_nan
- equivalence (i_nan,r_nan)
- !*** executable
- ! debug init local arrays
- i_nan=2147483647
- ua=r_nan
- va=r_nan
- altw=r_nan
- altub=r_nan
- hgtu=r_nan
- hgtv=r_nan
- if(kds.ne.1)then
- !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
- write(msg,*)'WARNING: bottom index kds=',kds,' should be 1?'
- call message(msg)
- !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
- endif
- ! ^ j
- ! ------------ |
- ! | | ----> i
- ! u p |
- ! | | nodes in cell (i,j)
- ! ------v----- view from top
- !
- ! v(ide,jde+1)
- ! -------x------
- ! | |
- ! | |
- ! u(ide,jde) x x x u(ide+1,jde)
- ! | p(ide,hde) |
- ! | | p=ph,phb,z0,...
- ! -------x------
- ! v(ide,jde)
- !
- ! staggered values set u(ids:ide+1,jds:jde) v(ids:ide,jds:jde+1)
- ! p=ph+phb set at (ids:ide,jds:jde)
- ! location of u(i,j) needs p(i-1,j) and p(i,j)
- ! location of v(i,j) needs p(i,j-1) and p(i,j)
- ! *** NOTE need HALO in ph, phb
- ! so we can compute only u(ids+1:ide,jds:jde) v(ids:ide,jds+1,jde)
- ! unless we extend p at the boundary
- ! but because we care about the fire way in the inside it does not matter
- ! if the fire gets close to domain boundary the simulation is over anyway
- ite1=snode(ite,ide,1)
- jte1=snode(jte,jde,1)
- ! do this in any case to check for nans
- call print_3d_stats(its,ite1,kds,kde,jts,jte,ims,ime,kms,kme,jms,jme,u,'wind U in')
- call print_3d_stats(its,ite,kds,kde,jts,jte1,ims,ime,kms,kme,jms,jme,v,'wind V in')
- if(fire_print_msg.gt.0)then
- !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
- write(msg,'(a,f7.2,a)')'interpolate_atm2fire: log-interpolation of wind to',fire_wind_height,' m'
- call message(msg)
- !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
- endif
- ! indexing
-
- ! for w
- itst=ifval(ids.eq.its,its,its-1)
- itet=ifval(ide.eq.ite,ite,ite+1)
- jtst=ifval(jds.ge.jts,jts,jts-1)
- jtet=ifval(jde.eq.jte,jte,jte+1)
- ! for u
- itsu=ifval(ids.eq.its,its+1,its) ! staggered direction
- iteu=ifval(ide.eq.ite,ite,ite+1) ! staggered direction
- jtsu=ifval(jds.ge.jts,jts,jts-1)
- jteu=ifval(jde.eq.jte,jte,jte+1)
- ! for v
- jtsv=ifval(jds.eq.jts,jts+1,jts) ! staggered direction
- jtev=ifval(jde.eq.jte,jte,jte+1) ! staggered direction
- itsv=ifval(ids.ge.its,its,its-1)
- itev=ifval(ide.eq.ite,ite,ite+1)
- if(fire_print_msg.ge.1)then
- !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
- write(msg,7001)'atm input ','tile',its,ite,jts,jte
- call message(msg)
- write(msg,7001)'altw ','tile',itst,itet,jtst,jtet
- call message(msg)
- !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
- endif
- 7001 format(a,' dimensions ',a4,':',i6,' to ',i6,' by ',i6,' to ',i6)
- !**********************************************************
- !* *
- !* find the altitude of the w points *
- !* *
- !**********************************************************
- !! in future, replace by z8w & test if the same
- kdmax=kde-1 ! max layer to interpolate from, can be less
- do j = jtst,jtet
- do k=kds,kdmax+1
- do i = itst,itet
- altw(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g ! altitude of the bottom w-point
- enddo
- enddo
- enddo
- ! values at u points
- if(fire_print_msg.ge.1)then
- !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
- write(msg,7001)'u interp at','tile',itsu,iteu,jtsu,jteu
- call message(msg)
- !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
- endif
- !**********************************************************
- !* *
- !* interpolate the altitude from w points to u points *
- !* *
- !**********************************************************
- do j = jtsu,jteu
- do k=kds,kdmax+1
- do i = itsu,iteu
- altub(i,k,j)= 0.5*(altw(i-1,k,j)+altw(i,k,j)) ! altitude of the bottom point under u-point
- enddo
- enddo
- do k=kds,kdmax
- do i = itsu,iteu
- hgtu(i,k,j) = 0.5*(altub(i,k,j)+altub(i,k+1,j)) - altub(i,kds,j) ! height of the u-point above the ground
- enddo
- enddo
- enddo
- ! values at v points
- if(fire_print_msg.ge.1)then
- !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
- write(msg,7001)'v interp at','tile',itsv,itev,jtsv,jtev
- call message(msg)
- !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
- endif
- !**********************************************************
- !* *
- !* interpolate the altitude from w points to v points *
- !* *
- !**********************************************************
- do j = jtsv,jtev
- do k=kds,kdmax+1
- do i = itsv,itev
- altvb(i,k,j)= 0.5*(altw(i,k,j-1)+altw(i,k,j)) ! altitude of the bottom point under v-point
- enddo
- enddo
- do k=kds,kdmax
- do i = itsv,itev
- hgtv(i,k,j) = 0.5*(altvb(i,k,j)+altvb(i,k+1,j)) - altvb(i,kds,j) ! height of the v-point above the ground
- enddo
- enddo
- enddo
- #ifdef DEBUG_OUT
- call write_array_m3(itsu,iteu,kds,kdmax,jtsu,jteu,its-2,ite+2,kds,kde,jts-2,jte+2,altub,'altub',id)
- call write_array_m3(itsv,itev,kds,kdmax,jtsv,jtev,its-2,ite+2,kds,kde,jts-2,jte+2,altvb,'altvb',id)
- call write_array_m3(itsu,iteu,kds,kdmax,jtsu,jteu,its-2,…
Large files files are truncated, but you can click here to view the full file