phys/module_fr_fire_atm.F

   1
   2 !WRF:MEDIATION_LAYER:FIRE_MODEL
   3 ! Routines dealing with the atmosphere
   4
   5 module module_fr_fire_atm
   6
   7 use module_model_constants, only: cp,xlv
   8 use module_fr_fire_util
   9 use module_state_description, only: num_tracer
  10 use module_state_description, only: p_fire_smoke
  11
  12 contains
  13
  14 ! DME subroutine for passive tracers
  15 subroutine add_fire_tracer_emissions(    &
  16        tracer_opt,dt,dx,dy,              &
  17        ifms,ifme,jfms,jfme,              &
  18        ifts,ifte,jtfs,jfte,              &
  19        ids,ide,kds,kde,jds,jde,          &
  20        ims,ime,kms,kme,jms,jme,          &
  21        its,ite,kts,kte,jts,jte,          &
  22        rho,dz8w,                         &
  23        burnt_area_dt,fgip,               &
  24        tracer,fire_tracer_smoke          &
  25 )
  26
  27 implicit none
  28 ! arguments
  29 integer,intent(in)::tracer_opt
  30 real,intent(in)::fire_tracer_smoke
  31 real,intent(in)::dt,dx,dy
  32 integer,intent(in)::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
  33 real,intent(in)::rho(ims:ime,kms:kme,jms:jme),dz8w(ims:ime,kms:kme,jms:jme)
  34 real,intent(in),dimension(ifms:ifme,jfms:jfme)::burnt_area_dt,fgip
  35 real,intent(inout)::tracer(ims:ime,kms:kme,jms:jme,num_tracer)
  36 ! local
  37 integer::isz1,jsz1,isz2,jsz2,ir,jr
  38 integer::i,j,ibase,jbase,i_f,ioff,j_f,joff
  39 real::avgw,emis,conv
  40
  41 isz1 = ite-its+1
  42 jsz1 = jte-jts+1
  43 isz2 = ifte-ifts+1
  44 jsz2 = jfte-jtfs+1
  45 ir=isz2/isz1
  46 jr=jsz2/jsz1
  47 avgw = 1.0/(ir*jr)
  48
  49 do j=max(jds+1,jts),min(jte,jde-2)
  50     jbase=jtfs+jr*(j-jts)
  51     do i=max(ids+1,its),min(ite,ide-2)
  52        ibase=ifts+ir*(i-its)
  53        do joff=0,jr-1
  54            j_f=joff+jbase
  55            do ioff=0,ir-1
  56                i_f=ioff+ibase
  57                if (num_tracer >0)then
  58                      emis=avgw*fire_tracer_smoke*burnt_area_dt(i_f,j_f)*fgip(i_f,j_f)*1000/(rho(i,kts,j)*dz8w(i,kts,j)) ! g_smoke/kg_air
  59                      tracer(i,kts,j,p_fire_smoke)=tracer(i,kts,j,p_fire_smoke)+emis
  60                 endif
  61            enddo
  62        enddo
  63     enddo
  64 enddo
  65
  66 end subroutine add_fire_tracer_emissions
  67
  68 !
  69 !***
  70 !
  71
  72 SUBROUTINE fire_tendency( &
  73     ids,ide, kds,kde, jds,jde,   & ! dimensions
  74     ims,ime, kms,kme, jms,jme,   &
  75     its,ite, kts,kte, jts,jte,   &
  76     grnhfx,grnqfx,canhfx,canqfx, & ! heat fluxes summed up to  atm grid
  77     alfg,alfc,z1can,             & ! coeffients, properties, geometry
  78     zs,z_at_w,dz8w,mu,c1h,c2h,rho, &
  79     rthfrten,rqvfrten)             ! theta and Qv tendencies
  80
  81 ! This routine is atmospheric physics
  82 ! it does NOT go into module_fr_fire_phys because it is not related to fire physical processes
  83
  84 ! --- this routine takes fire generated heat and moisture fluxes and
  85 !     calculates their influence on the theta and water vapor
  86 ! --- note that these tendencies are valid at the Arakawa-A location
  87
  88    IMPLICIT NONE
  89
  90 ! --- incoming variables
  91
  92    INTEGER , INTENT(in) :: ids,ide, kds,kde, jds,jde, &
  93                            ims,ime, kms,kme, jms,jme, &
  94                            its,ite, kts,kte, jts,jte
  95
  96    REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: grnhfx,grnqfx  ! W/m^2
  97    REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: canhfx,canqfx  ! W/m^2
  98    REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: zs  ! topography (m abv sealvl)
  99    REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: mu  ! dry air mass (Pa)
 100    REAL, INTENT(in), DIMENSION( kms:kme         ) :: c1h, c2h ! Hybrid coordinate weights
 101
 102    REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: z_at_w ! m abv sealvl
 103    REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: dz8w   ! dz across w-lvl
 104    REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: rho    ! density
 105
 106    REAL, INTENT(in) :: alfg ! extinction depth surface fire heat (m)
 107    REAL, INTENT(in) :: alfc ! extinction depth crown  fire heat (m)
 108    REAL, INTENT(in) :: z1can    ! height of crown fire heat release (m)
 109
 110 ! --- outgoing variables
 111
 112    REAL, INTENT(out), DIMENSION( ims:ime,kms:kme,jms:jme ) ::   &
 113        rthfrten, & ! theta tendency from fire (in mass units)
 114        rqvfrten    ! Qv tendency from fire (in mass units)
 115 ! --- local variables
 116
 117    INTEGER :: i,j,k
 118    INTEGER :: i_st,i_en, j_st,j_en, k_st,k_en
 119
 120    REAL :: cp_i
 121    REAL :: rho_i
 122    REAL :: xlv_i
 123    REAL :: z_w
 124    REAL :: fact_g, fact_c
 125    REAL :: alfg_i, alfc_i
 126
 127    REAL, DIMENSION( its:ite,kts:kte,jts:jte ) :: hfx,qfx
 128
 129 !!   character(len=128)::msg
 130
 131         do j=jts,jte
 132             do k=kts,min(kte+1,kde)
 133                do i=its,ite
 134                    rthfrten(i,k,j)=0.
 135                    rqvfrten(i,k,j)=0.
 136                enddo
 137             enddo
 138         enddo
 139
 140
 141 ! --- set some local constants
 142
 143
 144    cp_i = 1./cp     ! inverse of specific heat
 145    xlv_i = 1./xlv   ! inverse of latent heat
 146    alfg_i = 1./alfg
 147    alfc_i = 1./alfc
 148
 149 !!write(msg,'(8e11.3)')cp,cp_i,xlv,xlv_i,alfg,alfc,z1can
 150 !!call message(msg)
 151
 152    call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnhfx,'fire_tendency:grnhfx')
 153    call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnqfx,'fire_tendency:grnqfx')
 154
 155 ! --- set loop indicies : note that
 156
 157    i_st = MAX(its,ids+1)
 158    i_en = MIN(ite,ide-1)
 159    k_st = kts
 160    k_en = MIN(kte,kde-1)
 161    j_st = MAX(jts,jds+1)
 162    j_en = MIN(jte,jde-1)
 163
 164 ! --- distribute fluxes
 165
 166    DO j = j_st,j_en
 167       DO k = k_st,k_en
 168          DO i = i_st,i_en
 169
 170             ! --- set z (in meters above ground)
 171
 172             z_w = z_at_w(i,k,j) - zs(i,j) ! should be zero when k=k_st
 173
 174             ! --- heat flux
 175
 176             fact_g = cp_i * EXP( - alfg_i * z_w )
 177             IF ( z_w < z1can ) THEN
 178                fact_c = cp_i
 179             ELSE
 180                fact_c = cp_i * EXP( - alfc_i * (z_w - z1can) )
 181             END IF
 182             hfx(i,k,j) = fact_g * grnhfx(i,j) + fact_c * canhfx(i,j)
 183
 184 !!            write(msg,2)i,k,j,z_w,grnhfx(i,j),hfx(i,k,j)
 185 !!2           format('hfx:',3i4,6e11.3)
 186 !!            call message(msg)
 187
 188             ! --- vapor flux
 189
 190             fact_g = xlv_i * EXP( - alfg_i * z_w )
 191             IF (z_w < z1can) THEN
 192                fact_c = xlv_i
 193             ELSE
 194                fact_c = xlv_i * EXP( - alfc_i * (z_w - z1can) )
 195             END IF
 196             qfx(i,k,j) = fact_g * grnqfx(i,j) + fact_c * canqfx(i,j)
 197
 198 !!            if(hfx(i,k,j).ne.0. .or. qfx(i,k,j) .ne. 0.)then
 199 !!                write(msg,1)i,k,j,hfx(i,k,j),qfx(i,k,j)
 200 !!1               format('tend:',3i6,2e11.3)
 201 !!                call message(msg)
 202 !            endif
 203
 204          END DO
 205       END DO
 206    END DO
 207
 208 ! --- add flux divergence to tendencies
 209 !
 210 !   multiply by dry air mass (mu) to eliminate the need to
 211 !   call sr. calculate_phy_tend (in dyn_em/module_em.F)
 212
 213    DO j = j_st,j_en
 214       DO k = k_st,k_en-1
 215          DO i = i_st,i_en
 216
 217             rho_i = 1./rho(i,k,j)
 218
 219             rthfrten(i,k,j) = - (c1h(k)*mu(i,j)+c2h(k)) * rho_i * (hfx(i,k+1,j)-hfx(i,k,j)) / dz8w(i,k,j)
 220             rqvfrten(i,k,j) = - (c1h(k)*mu(i,j)+c2h(k)) * rho_i * (qfx(i,k+1,j)-qfx(i,k,j)) / dz8w(i,k,j)
 221
 222          END DO
 223       END DO
 224    END DO
 225
 226    call print_3d_stats(its,ite,kts,kte,jts,jte,ims,ime,kms,kme,jms,jme,rthfrten,'fire_tendency:rthfrten')
 227    call print_3d_stats(its,ite,kts,kte,jts,jte,ims,ime,kms,kme,jms,jme,rqvfrten,'fire_tendency:rqvfrten')
 228
 229    RETURN
 230
 231 END SUBROUTINE fire_tendency
 232
 233 !
 234 !***
 235 !
 236
 237 end module module_fr_fire_atm