diff --git a/wrfv2_fire/Registry/registry.fire_sfc_hypgrid b/wrfv2_fire/Registry/registry.fire_sfc_hypgrid new file mode 100644 index 00000000..ff303c64 --- /dev/null +++ b/wrfv2_fire/Registry/registry.fire_sfc_hypgrid @@ -0,0 +1,177 @@ +# +# ---------------------------------------- +# begin fire variables and configuration +# ---------------------------------------- +# +# declare fire package and choose which fire scheme +# +# +# name> namelist choice> state vars> +# +package fire_sfire ifire==2 - state:nfuel_cat,zsf,tign_g,rthfrten,rqvfrten,grnhfx,grnqfx,canhfx,canqfx,lfn,fuel_frac,fire_area,uf,vf,fgrnhfx,fgrnqfx,fcanhfx,fcanhfx,fcanqfx,fxlong,fxlat,fuel_time,bbb,betafl,phiwc,r_0,fgip,ischap + +# fire variables on fire grid +# +# +state real nfuel_cat *i*j fire 1 z i012hr "NFUEL_CAT" "fuel data" +state real zsf *i*j fire 1 z i012hr "ZSF" "height of surface above sea level" "m" +state real tign_g *i*j fire 1 z hr "TIGN_G" "ignition time on ground" "s" + +# fire variables on atm grid +# +state real rthfrten ikj fire 1 z hr "RTHFRTEN" "temperature tendency" "K/s" +state real rqvfrten ikj fire 1 z hr "RQVFRTEN" "humidity tendency" +state real grnhfx ij fire 1 z hr "GRNHFX" "heat flux from ground fire" "W/m^2" +state real grnqfx ij fire 1 z hr "GRNQFX" "moisture flux from ground fire" "W/m^2" +state real canhfx ij fire 1 z hr "CANHFX" "heat flux from crown fire" "W/m^2" +state real canqfx ij fire 1 z hr "CANQFX" "moisture flux from crown fire" "W/m^2" + +# sfire variables on fire grid +# (also using tign_g,zs,z_at_w,dz8w,nfuel_cat,fluxes,zsf) +# +state real lfn *i*j fire 1 z hr "LFN" "level function" "1" +state real fuel_frac *i*j fire 1 z hr "FUEL_FRAC" "fuel remaining" "1" +state real fire_area *i*j fire 1 z hr "FIRE_AREA" "fraction of cell area on fire" "1" +state real uf *i*j fire 1 z hr "UF" "fire wind" "m/s" +state real vf *i*j fire 1 z hr "VF" "fire wind" "m/s" +state real fgrnhfx *i*j fire 1 z hr "FGRNHFX" "heat flux from ground fire" "W/m^2" +state real fgrnqfx *i*j fire 1 z hr "FGRNQFX" "moisture flux from ground fire" "W/m^2" +state real fcanhfx *i*j fire 1 z hr "FCANHFX" "heat flux from crown fire" "W/m^2" +state real fcanqfx *i*j fire 1 z hr "FCANQFX" "moisture flux from crown fire" "W/m^2" +state real ros *i*j fire 1 z hr "ROS" "rate of spread" "m/s" + +# constant data arrays +state real fxlong *i*j fire 1 z hr "FXLONG" "longitude of midpoints of fire cells" "degrees" +state real fxlat *i*j fire 1 z hr "FXLAT" "latitude of midpoints of fire cells" "degrees" +state real fuel_time *i*j fire 1 z hr "FUEL_TIME" "fuel" +state real bbb *i*j fire 1 z hr "BBB" "fuel" +state real betafl *i*j fire 1 z hr "BETAFL" "fuel" +state real phiwc *i*j fire 1 z hr "PHIWC" "fuel" +state real r_0 *i*j fire 1 z hr "R_0" "fuel" +state real fgip *i*j fire 1 z hr "FGIP" "fuel" +state integer ischap *i*j fire 1 z hr "ISCHAP" "fuel" + +# +# fire configure namelist variables +# +#
+rconfig integer ifire namelist,fire max_domains 0 +rconfig integer fire_boundary_guard namelist,fire max_domains 2 - "fire_boundary_guard" "cells to stop when fire close to domain boundary" +# ignition for sfire +rconfig integer fire_num_ignitions namelist,fire max_domains 0. - "fire_num_ignitions" "number of ignition lines" +rconfig real fire_ignition_start_lon1 namelist,fire max_domains 0. - "fire_ignition_start_long1" "long coord of start of ignition line" "m" +rconfig real fire_ignition_start_lat1 namelist,fire max_domains 0. - "fire_ignition_start_lat1" "lat coord of start of ignition line" "m" +rconfig real fire_ignition_end_lon1 namelist,fire max_domains 0. - "fire_ignition_end_long1" "long coord of end of ignition line" "m" +rconfig real fire_ignition_end_lat1 namelist,fire max_domains 0. - "fire_ignition_end_lat1" "lat coord of end of ignition line" "m" +rconfig real fire_ignition_radius1 namelist,fire max_domains 0. - "fire_ignition_radius1" "ignite all within the radius" "m" +rconfig real fire_ignition_time1 namelist,fire max_domains 0. - "fire_ignition_time1" "ignition time" "s" +rconfig real fire_ignition_start_lon2 namelist,fire max_domains 0. - "fire_ignition_start_long2" "long coord of start of ignition line" "m" +rconfig real fire_ignition_start_lat2 namelist,fire max_domains 0. - "fire_ignition_start_lat2" "lat coord of start of ignition line" "m" +rconfig real fire_ignition_end_lon2 namelist,fire max_domains 0. - "fire_ignition_end_long2" "long coord of end of ignition line" "m" +rconfig real fire_ignition_end_lat2 namelist,fire max_domains 0. - "fire_ignition_end_lat2" "lat coord of end of ignition line" "m" +rconfig real fire_ignition_radius2 namelist,fire max_domains 0. - "fire_ignition_radius2" "ignite all within the radius" "m" +rconfig real fire_ignition_time2 namelist,fire max_domains 0. - "fire_ignition_time2" "ignition time" "s" +rconfig real fire_ignition_start_lon3 namelist,fire max_domains 0. - "fire_ignition_start_long3" "long coord of start of ignition line" "m" +rconfig real fire_ignition_start_lat3 namelist,fire max_domains 0. - "fire_ignition_start_lat3" "lat coord of start of ignition line" "m" +rconfig real fire_ignition_end_lon3 namelist,fire max_domains 0. - "fire_ignition_end_long3" "long coord of end of ignition line" "m" +rconfig real fire_ignition_end_lat3 namelist,fire max_domains 0. - "fire_ignition_end_lat3" "lat coord of end of ignition line" "m" +rconfig real fire_ignition_radius3 namelist,fire max_domains 0. - "fire_ignition_radius3" "ignite all within the radius" "m" +rconfig real fire_ignition_time3 namelist,fire max_domains 0. - "fire_ignition_time3" "ignition time" "s" +rconfig real fire_ignition_start_lon4 namelist,fire max_domains 0. - "fire_ignition_start_long4" "long coord of start of ignition line" "m" +rconfig real fire_ignition_start_lat4 namelist,fire max_domains 0. - "fire_ignition_start_lat4" "lat coord of start of ignition line" "m" +rconfig real fire_ignition_end_lon4 namelist,fire max_domains 0. - "fire_ignition_end_long4" "long coord of end of ignition line" "m" +rconfig real fire_ignition_end_lat4 namelist,fire max_domains 0. - "fire_ignition_end_lat4" "lat coord of end of ignition line" "m" +rconfig real fire_ignition_radius4 namelist,fire max_domains 0. - "fire_ignition_radius4" "ignite all within the radius" "m" +rconfig real fire_ignition_time4 namelist,fire max_domains 0. - "fire_ignition_time4" "ignition time" "s" +rconfig real fire_ignition_start_lon5 namelist,fire max_domains 0. - "fire_ignition_start_long5" "long coord of start of ignition line" "m" +rconfig real fire_ignition_start_lat5 namelist,fire max_domains 0. - "fire_ignition_start_lat5" "lat coord of start of ignition line" "m" +rconfig real fire_ignition_end_lon5 namelist,fire max_domains 0. - "fire_ignition_end_long5" "long coord of end of ignition line" "m" +rconfig real fire_ignition_end_lat5 namelist,fire max_domains 0. - "fire_ignition_end_lat5" "lat coord of end of ignition line" "m" +rconfig real fire_ignition_radius5 namelist,fire max_domains 0. - "fire_ignition_radius5" "ignite all within the radius" "m" +rconfig real fire_ignition_time5 namelist,fire max_domains 0. - "fire_ignition_time5" "ignition time" "s" +rconfig real fire_ignition_start_x1 namelist,fire max_domains 0. - "fire_ignition_start_x1" "x coord of start of ignition line" "m" +rconfig real fire_ignition_start_y1 namelist,fire max_domains 0. - "fire_ignition_start_y1" "y coord of start of ignition line" "m" +rconfig real fire_ignition_end_x1 namelist,fire max_domains 0. - "fire_ignition_end_x1" "x coord of end of ignition line" "m" +rconfig real fire_ignition_end_y1 namelist,fire max_domains 0. - "fire_ignition_end_y1" "y coord of end of ignition line" "m" +rconfig real fire_ignition_start_x2 namelist,fire max_domains 0. - "fire_ignition_start_x2" "x coord of start of ignition line" "m" +rconfig real fire_ignition_start_y2 namelist,fire max_domains 0. - "fire_ignition_start_y2" "y coord of start of ignition line" "m" +rconfig real fire_ignition_end_x2 namelist,fire max_domains 0. - "fire_ignition_end_x2" "x coord of end of ignition line" "m" +rconfig real fire_ignition_end_y2 namelist,fire max_domains 0. - "fire_ignition_end_y2" "y coord of end of ignition line" "m" +rconfig real fire_ignition_start_x3 namelist,fire max_domains 0. - "fire_ignition_start_x3" "x coord of start of ignition line" "m" +rconfig real fire_ignition_start_y3 namelist,fire max_domains 0. - "fire_ignition_start_y3" "y coord of start of ignition line" "m" +rconfig real fire_ignition_end_x3 namelist,fire max_domains 0. - "fire_ignition_end_x3" "x coord of end of ignition line" "m" +rconfig real fire_ignition_end_y3 namelist,fire max_domains 0. - "fire_ignition_end_y3" "y coord of end of ignition line" "m" +rconfig real fire_ignition_start_x4 namelist,fire max_domains 0. - "fire_ignition_start_x4" "x coord of start of ignition line" "m" +rconfig real fire_ignition_start_y4 namelist,fire max_domains 0. - "fire_ignition_start_y4" "y coord of start of ignition line" "m" +rconfig real fire_ignition_end_x4 namelist,fire max_domains 0. - "fire_ignition_end_x4" "x coord of end of ignition line" "m" +rconfig real fire_ignition_end_y4 namelist,fire max_domains 0. - "fire_ignition_end_y4" "y coord of end of ignition line" "m" +rconfig real fire_ignition_start_x5 namelist,fire max_domains 0. - "fire_ignition_start_x5" "x coord of start of ignition line" "m" +rconfig real fire_ignition_start_y5 namelist,fire max_domains 0. - "fire_ignition_start_y5" "y coord of start of ignition line" "m" +rconfig real fire_ignition_end_x5 namelist,fire max_domains 0. - "fire_ignition_end_x5" "x coord of end of ignition line" "m" +rconfig real fire_ignition_end_y5 namelist,fire max_domains 0. - "fire_ignition_end_y5" "y coord of end of ignition line" "m" +# variables from old cawfe code +rconfig real fire_lat_init namelist,fire max_domains 0. - "fire_lat_init" "latitude to start fire" "degrees" +rconfig real fire_lon_init namelist,fire max_domains 0. - "fire_lon_init" "longitude to start fire" "degrees" +rconfig real fire_ign_time namelist,fire max_domains 0. - "fire_ign_time" "time when fire should be ignited" "min" +rconfig integer fire_shape namelist,fire max_domains 0 - "fire_shape" "fire shape" "" +rconfig integer fire_sprd_mdl namelist,fire max_domains 1 - "fire_sprd_mdl" "which spread rate formula: if 0, Macarthur; if 1, BEHAVE" "" +rconfig real fire_crwn_hgt namelist,fire max_domains 15. - "fire_crwn_hgt" "height that heat from crown fire is released" "m" +rconfig real fire_ext_grnd namelist,fire max_domains 50. - "fire_ext_grnd" "extinction depth of sfc fire heat" "m" +rconfig real fire_ext_crwn namelist,fire max_domains 50. - "fire_ext_crwn" "extinction depth of crown fire heat" "m" +rconfig integer fire_fuel_read namelist,fire max_domains -1 - "fire_fuel_read" "fuel categories are set by: if 0, uniform; if 1, user-presc; if 2, read from file" "" +rconfig integer fire_fuel_cat namelist,fire max_domains 1 - "fire_fuel_cat" "fuel category if ifuelread=0" "" +# sfire switches +rconfig integer fire_print_msg namelist,fire max_domains 0 - "fire_write_msg" "write fire statistics, 0 no writes, 1+ for more" "" +rconfig integer fire_print_file namelist,fire max_domains 0 - "fire_write_file" "write fire output text files, 0 no writes, 1+ for more" "" +# method selection` +rconfig integer fire_fuel_left_method namelist,fire max_domains 1 - "fire_fuel_left_method" "1 or 2, compute fuel_left" "" +rconfig integer fire_fuel_left_irl namelist,fire max_domains 2 - "fire_fuel_left_irl" "submesh to compute fuel lwft, even, at least 2" "" +rconfig integer fire_fuel_left_jrl namelist,fire max_domains 2 - "fire_fuel_left_jrl" "submesh to compute fuel lwft, even, at least 2" "" +rconfig real fire_back_weight namelist,fire max_domains 0.5 - "fire_back_weight" "RK timestepping coefficient, 0=forward, 0.5=Heun" "1" +rconfig integer fire_grows_only namelist,fire max_domains 1 - "fire_grows_only" "if >0 level set function cannot increase = fire can only grow" "1" +rconfig integer fire_upwinding namelist,fire max_domains 3 - "fire_upwinding" "upwind normal spread: 1=standard, 2=godunov, 3=eno, 4=sethian" "1" +rconfig integer fire_upwind_split namelist,fire max_domains 0 - "fire_upwind_split" "1=upwind advection separately from normal direction spread" "1" +rconfig real fire_viscosity namelist,fire max_domains 0.4 - "fire_viscosity" "artificial viscosity in level set method" "1" +rconfig real fire_lfn_ext_up namelist,fire max_domains 1.0 - "fire_lfn_ext_up" "0.=extend level set function at boundary by reflection, 1.=always up" "1" +rconfig integer fire_topo_from_atm namelist,fire max_domains 1 - "fire_topo_from_atm" "0 = do nothing, 1 = populate ZSF by interpolating from atmosphere" "1" +rconfig integer fire_advection namelist,fire max_domains 1 - "fire_advection" "0 = fire spread computed from normal wind speed/slope, 1 = fireline particle speed projected on normal" "0" +# experiments +rconfig integer fire_test_steps namelist,fire max_domains 0 - "fire_test_steps" ">0 = on first call, do specified number of steps and terminate (testing only)" "1" +rconfig real fire_const_time namelist,fire max_domains -1. - "fire_const_time" "time from ignition to freeze fire, <0 never" "s" +rconfig real fire_const_grnhfx namelist,fire max_domains 0. - "fire_const_grnhfx" "if both >=0, the amount of constant heat flux" "1" +rconfig real fire_const_grnqfx namelist,fire max_domains 0. - "fire_const_grnqfx" "if both >=0, the amount of constant heat flux" "1" +rconfig real fire_atm_feedback namelist,fire max_domains 1. - "fire_atm_feedback" "the heat fluxes to the atmosphere are multiplied by this" "1" +rconfig integer fire_mountain_type namelist,fire max_domains 0 - "fire_mountain_type" "in ideal: 0=none, 1=COS hill, 2=EW ridge, 3=NS ridge" "1" +rconfig real fire_mountain_height namelist,fire max_domains 500. - "fire_mountain_height" "ideal mountain height" "m" +rconfig real fire_mountain_start_x namelist,fire max_domains 100. - "fire_mountain_start_x" "x coord of start of the mountain" "m" +rconfig real fire_mountain_start_y namelist,fire max_domains 100. - "fire_mountain_start_y" "y coord of start of the mountain" "m" +rconfig real fire_mountain_end_x namelist,fire max_domains 100. - "fire_mountain_end_x" "x coord of end of the mountain" "m" +rconfig real fire_mountain_end_y namelist,fire max_domains 100. - "fire_mountain_end_y" "y coord of end of the mountain" "m" +rconfig real delt_perturbation namelist,fire max_domains 3.0 - "delt_perturbation" "temperature perturbation for cold (-) /warm (+) bubble" "K" +rconfig real xrad_perturbation namelist,fire max_domains 10000.0 - "xrad_perturbation" "horizontal radius of the perturbation in E-W direction" "m" +rconfig real yrad_perturbation namelist,fire max_domains 10000.0 - "yrad_perturbation" "horizontal radius of the perturbation in N-S direction" "m" +rconfig real zrad_perturbation namelist,fire max_domains 1500.0 - "zrad_perturbation" "vertical radius of the perturbation (bubble) direction" "m" +rconfig real hght_perturbation namelist,fire max_domains 1500.0 - "hght_perturbation" "height at which the perturbation (bubble) will be suspended" "m" +rconfig real z_grd_scale namelist,fire max_domains 0.4 - "z_grd_scale" "z scale for grid stretching" "m" +rconfig logical stretch_grd namelist,fire max_domains .true. - "switch for turning on grid stretching" "" +rconfig logical stretch_hyp namelist,fire max_domains .false. - "switch for turning on hyperbolic grid stretching" "m" + + + +# +# Fire halo descriptions +# +halo HALO_FIRE_LFN dyn_em 24:lfn +halo HALO_FIRE_TIGN dyn_em 8:tign_g +halo HALO_FIRE_HT dyn_em 8:ht +halo HALO_FIRE_WIND_F dyn_em 12:uf,vf +halo HALO_FIRE_LONGLAT dyn_em 24:xlong,xlat +halo HALO_FIRE_WIND_A dyn_em 8:u_2,v_2 +halo HALO_FIRE_ZSF dyn_em 24:zsf +halo HALO_FIRE_FUEL dyn_em 8:fuel_frac,fuel_time,bbb,betafl,phiwc,r_0,fgip,ischap,nfuel_cat +# +# ---------------------------------------- +# end fire variables and configuration +# ---------------------------------------- + +## diff --git a/wrfv2_fire/dyn_em/module_initialize_fire_sfc_hypgrid.F b/wrfv2_fire/dyn_em/module_initialize_fire_sfc_hypgrid.F new file mode 100644 index 00000000..6276f64b --- /dev/null +++ b/wrfv2_fire/dyn_em/module_initialize_fire_sfc_hypgrid.F @@ -0,0 +1,983 @@ +!IDEAL:MODEL_LAYER:INITIALIZATION +! + +! This MODULE holds the routines which are used to perform various initializations +! for the individual domains. + +! This MODULE CONTAINS the following routines: + +! initialize_field_test - 1. Set different fields to different constant +! values. This is only a test. If the correct +! domain is not found (based upon the "id") +! then a fatal error is issued. + +!----------------------------------------------------------------------- + +MODULE module_initialize_ideal + + USE module_domain + USE module_io_domain + USE module_state_description + USE module_model_constants + USE module_bc + USE module_timing + USE module_configure + USE module_init_utilities + USE module_soil_pre !AK/ak +#ifdef DM_PARALLEL + USE module_dm +#endif + + +CONTAINS + + +!------------------------------------------------------------------- +! this is a wrapper for the solver-specific init_domain routines. +! Also dereferences the grid variables and passes them down as arguments. +! This is crucial, since the lower level routines may do message passing +! and this will get fouled up on machines that insist on passing down +! copies of assumed-shape arrays (by passing down as arguments, the +! data are treated as assumed-size -- ie. f77 -- arrays and the copying +! business is avoided). Fie on the F90 designers. Fie and a pox. + + SUBROUTINE init_domain ( grid ) + + IMPLICIT NONE + + ! Input data. + TYPE (domain), POINTER :: grid + ! Local data. + INTEGER :: idum1, idum2 + + CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 ) + + CALL init_domain_rk( grid & +! +#include +! + ) + + END SUBROUTINE init_domain + +!------------------------------------------------------------------- + + SUBROUTINE init_domain_rk ( grid & +! +# include +! +) + IMPLICIT NONE + + ! Input data. + TYPE (domain), POINTER :: grid + +# include + + TYPE (grid_config_rec_type) :: config_flags + + ! Local data + INTEGER :: & + ids, ide, jds, jde, kds, kde, & + ims, ime, jms, jme, kms, kme, & + its, ite, jts, jte, kts, kte, & + i, j, k + + ! Local data + + INTEGER, PARAMETER :: nl_max = 1000 + REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in + INTEGER :: nl_in + + + INTEGER :: icm,jcm, ii, im1, jj, jm1, loop, error, fid, nxc, nyc + REAL :: u_mean,v_mean, f0, p_surf, p_level, qvf, z_at_v, z_at_u + REAL :: z_scale, xrad, yrad, zrad, rad, delt, cof1, cof2 + REAL :: x_rad, y_rad, z_rad, hght_pert !Ak/ak + character (len=256) :: mminlu2 !AK/ak + CHARACTER (len=80) :: etafile !AK/ak +! REAL, EXTERNAL :: interp_0 + REAL :: hm + REAL :: pi + +! stuff from original initialization that has been dropped from the Registry + REAL :: vnu, xnu, xnus, dinit0, cbh, p0_temp, t0_temp, zd, zt + REAL :: qvf1, qvf2, pd_surf + INTEGER :: it + real :: thtmp, ptmp, temp(3) + + LOGICAL :: moisture_init + LOGICAL :: stretch_grid, dry_sounding + LOGICAL :: stretch_hyp + + INTEGER :: xs , xe , ys , ye + INTEGER :: mtn_type + INTEGER :: & ! fire mesh sizes + ifds,ifde, kfds,kfde, jfds,jfde, & + ifms,ifme, kfms,kfme, jfms,jfme, & + ifts,ifte, kfts,kfte, jfts,jfte + + REAL :: mtn_ht, mtn_max, mtn_fdx, mtn_fdy, mtn_x, mtn_y, mtn_z + REAL :: mtn_axs, mtn_ays, mtn_axe, mtn_aye + REAL :: mtn_fxs, mtn_fys, mtn_fxe, mtn_fye + REAL :: mtn_xs, mtn_ys, mtn_xe, mtn_ye + LOGICAL, EXTERNAL :: wrf_dm_on_monitor + + SELECT CASE ( model_data_order ) + CASE ( DATA_ORDER_ZXY ) + kds = grid%sd31 ; kde = grid%ed31 ; + ids = grid%sd32 ; ide = grid%ed32 ; + jds = grid%sd33 ; jde = grid%ed33 ; + + kms = grid%sm31 ; kme = grid%em31 ; + ims = grid%sm32 ; ime = grid%em32 ; + jms = grid%sm33 ; jme = grid%em33 ; + + kts = grid%sp31 ; kte = grid%ep31 ; ! note that tile is entire patch + its = grid%sp32 ; ite = grid%ep32 ; ! note that tile is entire patch + jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch + CASE ( DATA_ORDER_XYZ ) + ids = grid%sd31 ; ide = grid%ed31 ; + jds = grid%sd32 ; jde = grid%ed32 ; + kds = grid%sd33 ; kde = grid%ed33 ; + + ims = grid%sm31 ; ime = grid%em31 ; + jms = grid%sm32 ; jme = grid%em32 ; + kms = grid%sm33 ; kme = grid%em33 ; + + its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch + jts = grid%sp32 ; jte = grid%ep32 ; ! note that tile is entire patch + kts = grid%sp33 ; kte = grid%ep33 ; ! note that tile is entire patch + CASE ( DATA_ORDER_XZY ) + ids = grid%sd31 ; ide = grid%ed31 ; + kds = grid%sd32 ; kde = grid%ed32 ; + jds = grid%sd33 ; jde = grid%ed33 ; + + ims = grid%sm31 ; ime = grid%em31 ; + kms = grid%sm32 ; kme = grid%em32 ; + jms = grid%sm33 ; jme = grid%em33 ; + + its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch + kts = grid%sp32 ; kte = grid%ep32 ; ! note that tile is entire patch + jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch + + END SELECT + +! z_scale = .40 + pi = 2.*asin(1.0) + write(6,*) ' pi is ',pi + nxc = (ide-ids)/2 + nyc = (jde-jds)/2 + + CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags ) + +! here we check to see if the boundary conditions are set properly + + CALL boundary_condition_check( config_flags, bdyzone, error, grid%id ) + + delt = config_flags%delt_perturbation + x_rad = config_flags%xrad_perturbation + y_rad = config_flags%yrad_perturbation + z_rad = config_flags%zrad_perturbation + hght_pert = config_flags%hght_perturbation + + stretch_grid = config_flags%stretch_grd + stretch_hyp = config_flags%stretch_hyp + z_scale = config_flags%z_grd_scale +! delt = 3. +! x_rad = 10000. +! y_rad = 10000. +! z_rad = 1500. +! hght_pert = 1500. + + moisture_init = .true. + + grid%itimestep=0 + +#ifdef DM_PARALLEL + CALL wrf_dm_bcast_bytes( icm , IWORDSIZE ) + CALL wrf_dm_bcast_bytes( jcm , IWORDSIZE ) +#endif + +!AK/ak land use initialization + mminlu2 = ' ' !Ak/ak + mminlu2(1:4) = 'USGS' !Ak/ak + CALL nl_set_mminlu(1, mminlu2) !Ak/ak + CALL nl_set_iswater(1,0) + CALL nl_set_cen_lat(1,40.) + CALL nl_set_cen_lon(1,-105.) + CALL nl_set_truelat1(1,0.) + CALL nl_set_truelat2(1,0.) + CALL nl_set_moad_cen_lat (1,0.) + CALL nl_set_stand_lon (1,0.) +! CALL nl_set_pole_lon (1,0.) !Ak/ak +! CALL nl_set_pole_lat (1,90.) !Ak/ak + CALL nl_set_map_proj(1,0) + CALL nl_get_iswater(1,grid%iswater) ! Ak/sk + +! here we initialize data we currently is not initialized +! in the input data + + DO j = jts, jte + DO i = its, ite + grid%msftx(i,j) = 1. + grid%msfty(i,j) = 1. + grid%msfux(i,j) = 1. + grid%msfuy(i,j) = 1. + grid%msfvx(i,j) = 1. + grid%msfvx_inv(i,j)= 1. + grid%msfvy(i,j) = 1. + grid%sina(i,j) = 0. + grid%cosa(i,j) = 1. + grid%e(i,j) = 0. + grid%f(i,j) = 0. + + grid%xlat(i,j) = 35. !Ak/sk + grid%xlong(i,j) = -111. !Ak/ak + grid%xland(i,j) = 1. !Ak/ak + grid%lu_index(i,j) = 28 !AK/ak + grid%tsk(i,j) = 285.0 !AK/ak + grid%tmn(i,j) = 285.0 !AK/ak + END DO + END DO + + +! for Noah LSM, additional variables need to be initializedi !AK/ak | + + other_masked_fields : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) ) + + CASE (SLABSCHEME) + + CASE (LSMSCHEME) + + DO j = jts , MIN(jde-1,jte) + DO i = its , MIN(ide-1,ite) + grid%vegfra(i,j) = 0.5 + grid%canwat(i,j) = 0. + grid%ivgtyp(i,j) = 18 + grid%isltyp(i,j) = 7 + grid%xice(i,j) = 0. + grid%snow(i,j) = 0. + END DO + END DO + + CASE (RUCLSMSCHEME) + ! ^ + END SELECT other_masked_fields !AK/ak | + + + DO j = jts, jte + DO k = kts, kte + DO i = its, ite + grid%ww(i,k,j) = 0. + END DO + END DO + END DO + + grid%step_number = 0 + +! Process the soil; note that there are some things hard-wired into share/module_soil_pre.F + CALL process_soil_ideal(grid%xland,grid%xice,grid%vegfra,grid%snow,grid%canwat, & + grid%ivgtyp,grid%isltyp,grid%tslb,grid%smois, & + grid%tsk,grid%tmn,grid%zs,grid%dzs,model_config_rec%num_soil_layers, & + model_config_rec%sf_surface_physics(grid%id), & + ids,ide, jds,jde, kds,kde,& + ims,ime, jms,jme, kms,kme,& + its,ite, jts,jte, kts,kte ) + +! set up the grid + IF (stretch_grid) THEN ! exponential or hyperbolic tangential stretch for eta + write(6,*) ' stretched grid with z_scale =',z_scale + IF (stretch_hyp) THEN ! hyperbolic tangential stretch (more levels at the surface) + write(6,*) ' hyperbolic stretching activated' + DO k=1, kde + grid%znw(k) = -1.* (tanh(z_scale*(float(k-1) / float(kde-1) -1.)))/ & + (tanh(z_scale)) + + ENDDO + ELSE ! exponential stretch for eta (nearly constant dz) + + DO k=1, kde + grid%znw(k) = (exp(-(k-1)/float(kde-1)/z_scale) - exp(-1./z_scale))/ & + (1.-exp(-1./z_scale)) + ENDDO + ENDIF + ELSE + write(6,*) 'no grid stretching' + DO k=1, kde + grid%znw(k) = 1. - float(k-1)/float(kde-1) + ENDDO + ENDIF + + + DO k=1, kde-1 + grid%dnw(k) = grid%znw(k+1) - grid%znw(k) + grid%rdnw(k) = 1./grid%dnw(k) + grid%znu(k) = 0.5*(grid%znw(k+1)+grid%znw(k)) + ENDDO + DO k=2, kde-1 + grid%dn(k) = 0.5*(grid%dnw(k)+grid%dnw(k-1)) + grid%rdn(k) = 1./grid%dn(k) + grid%fnp(k) = .5* grid%dnw(k )/grid%dn(k) + grid%fnm(k) = .5* grid%dnw(k-1)/grid%dn(k) + ENDDO + + cof1 = (2.*grid%dn(2)+grid%dn(3))/(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(2) + cof2 = grid%dn(2) /(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(3) + grid%cf1 = grid%fnp(2) + cof1 + grid%cf2 = grid%fnm(2) - cof1 - cof2 + grid%cf3 = cof2 + + grid%cfn = (.5*grid%dnw(kde-1)+grid%dn(kde-1))/grid%dn(kde-1) + grid%cfn1 = -.5*grid%dnw(kde-1)/grid%dn(kde-1) + grid%rdx = 1./config_flags%dx + grid%rdy = 1./config_flags%dy + +! get the sounding from the ascii sounding file, first get dry sounding and +! calculate base state + + dry_sounding = .true. + IF ( wrf_dm_on_monitor() ) THEN + write(6,*) ' getting dry sounding for base state ' + + CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in ) + ENDIF + CALL wrf_dm_bcast_real( zk , nl_max ) + CALL wrf_dm_bcast_real( p_in , nl_max ) + CALL wrf_dm_bcast_real( pd_in , nl_max ) + CALL wrf_dm_bcast_real( theta , nl_max ) + CALL wrf_dm_bcast_real( rho , nl_max ) + CALL wrf_dm_bcast_real( u , nl_max ) + CALL wrf_dm_bcast_real( v , nl_max ) + CALL wrf_dm_bcast_real( qv , nl_max ) + CALL wrf_dm_bcast_integer ( nl_in , 1 ) + + write(6,*) ' returned from reading sounding, nl_in is ',nl_in + +! find ptop for the desired ztop (ztop is input from the namelist), +! and find surface pressure + + grid%p_top = interp_0( p_in, zk, config_flags%ztop, nl_in ) + +! get fire mesh dimensions + CALL get_ijk_from_subgrid ( grid , & + ifds,ifde, jfds,jfde,kfds,kfde, & + ifms,ifme, jfms,jfme,kfms,kfme, & + ifts,ifte, jfts,jfte,kfts,kfte) + + write (6,*)' ******** SFIRE ideal initialization ********' + + ! fire grid step size + mtn_fdx = grid%dx/grid%sr_x + mtn_fdy = grid%dy/grid%sr_y + + write (6,*)' atm mesh step ',grid%dx,grid%dy + write (6,*)' fire mesh step ',mtn_fdx,mtn_fdy + write (6,*)' refinement ratio ',grid%sr_x,grid%sr_y + write (6,*)' atm domain bounds ',ids,ide, jds,jde,kds,kde + write (6,*)' atm memory bounds ',ims,ime, jms,jme,kms,kme + write (6,*)' atm tile bounds ',its,ite, jts,jte,kts,kte + write (6,*)' fire domain bounds ',ifds,ifde, jfds,jfde,kfds,kfde + write (6,*)' fire memory bounds ',ifms,ifme, jfms,jfme,kfms,kfme + write (6,*)' fire tile bounds ',ifts,ifte, jfts,jfte,kfts,kfte + write (6,*)' Note that atm mesh and fire mesh are cell-centered' + + DO j=jts,jte + DO i=its,ite + grid%ht(i,j) = 0. + ENDDO + ENDDO + + !******* set terrain height + + ! copy params from the namelist + mtn_type = config_flags%fire_mountain_type + mtn_xs = config_flags%fire_mountain_start_x + mtn_ys = config_flags%fire_mountain_start_y + mtn_xe = config_flags%fire_mountain_end_x + mtn_ye = config_flags%fire_mountain_end_y + mtn_ht = config_flags%fire_mountain_height + + IF(mtn_type .ne. 0)THEN + + ! atmospheric grid coordinates of the mountain + mtn_axs = mtn_xs/grid%dx + ids - 0.5 + mtn_axe = mtn_xe/grid%dx + ids - 0.5 + mtn_ays = mtn_ys/grid%dy + jds - 0.5 + mtn_aye = mtn_ye/grid%dy + jds - 0.5 + + ! fire grid coordinates of the mountain + mtn_fxs = mtn_xs/mtn_fdx + ifds - 0.5 + mtn_fxe = mtn_xe/mtn_fdx + ifds - 0.5 + mtn_fys = mtn_ys/mtn_fdy + jfds - 0.5 + mtn_fye = mtn_ye/mtn_fdy + jfds - 0.5 + + write(6,*)' Mountain height ',mtn_ht,' type',mtn_type + write(6,*)' Mountain (m) LL=(0,0) ',mtn_xs,':',mtn_xe,' ',mtn_ys,':',mtn_ye + write(6,*)' Mountain on atm grid ',mtn_axs,':',mtn_axe,' ',mtn_ays,':',mtn_aye + write(6,*)' Mountain on fire grid ',mtn_fxs,':',mtn_fxe,' ',mtn_fys,':',mtn_fye + + mtn_max = 0. + DO j=jts,jte + DO i=its,ite + mtn_x = pi + 2*pi* max(0. , min( (i - mtn_axs)/(mtn_axe - mtn_axs), 1. )) + mtn_y = pi + 2*pi* max(0. , min( (j - mtn_ays)/(mtn_aye - mtn_ays), 1. )) + SELECT CASE(mtn_type) + CASE (1) ! circ/elliptic mountain + mtn_z = mtn_ht * 0.25 * (1. + COS(mtn_x))*(1. + COS(mtn_y)) + CASE (2) ! EW ridge + mtn_z = mtn_ht * 0.5 * (1. + COS(mtn_y)) + CASE (3) ! NS ridge + mtn_z = mtn_ht * 0.5 * (1. + COS(mtn_x)) + CASE DEFAULT + call wrf_error_fatal ( ' bad fire_mountain_type ' ) + END SELECT + mtn_max = max(mtn_max, mtn_z) + grid%ht(i,j) = mtn_z + ENDDO + ENDDO + + write(6, *)' Atm tile ',its,':',ite,' ',jts,':',jte,' max terrain height ',mtn_max + + IF (config_flags%fire_topo_from_atm == 0) THEN + + mtn_max = 0. + DO j=jfts,jfte + DO i=ifts,ifte + mtn_x = pi + 2*pi* max(0. , min( (i - mtn_fxs)/(mtn_fxe - mtn_fxs), 1. )) + mtn_y = pi + 2*pi* max(0. , min( (j - mtn_fys)/(mtn_fye - mtn_fys), 1. )) + SELECT CASE(mtn_type) + CASE (1) ! circ/elliptic mountain + mtn_z = mtn_ht * 0.25 * (1. + COS(mtn_x))*(1. + COS(mtn_y)) + CASE (2) ! EW ridge + mtn_z = mtn_ht * 0.5 * (1. + COS(mtn_y)) + CASE (3) ! NS ridge + mtn_z = mtn_ht * 0.5 * (1. + COS(mtn_x)) + CASE DEFAULT + call wrf_error_fatal ( ' bad fire_mountain_type ' ) + END SELECT + mtn_max = max(mtn_max, mtn_z) + grid%zsf(i,j) = mtn_z + ENDDO + ENDDO + + write(6, *)' Fire tile ',ifts,':',ifte,' ',jfts,':',jfte,' max terrain height ',mtn_max + write(6, *)' Fine-resolution terrain will be used.' + + ELSE + + write(6, *)' ******** WARNING ********' + write(6, *)' Terrain height on the fire grid will be interpolated from the atmosphere grid.' + write(6, *)' Set fire_topo_from_atm=0 in the namelist to use fine-resolution terrain.' + write(6, *)' ******** WARNING ********' + + ENDIF + ENDIF + + DO j=jts,jte + DO i=its,ite + grid%phb(i,1,j) = g * grid%ht(i,j) + grid%ph0(i,1,j) = g * grid%ht(i,j) + ENDDO + ENDDO + + DO J = jts, jte + DO I = its, ite + + p_surf = interp_0( p_in, zk, grid%phb(i,1,j)/g, nl_in ) + grid%mub(i,j) = p_surf-grid%p_top + +! this is dry hydrostatic sounding (base state), so given grid%p (coordinate), +! interp theta (from interp) and compute 1/rho from eqn. of state + + DO K = 1, kte-1 + p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top + grid%pb(i,k,j) = p_level + grid%t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0 + grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm + ENDDO + +! calc hydrostatic balance (alternatively we could interp the geopotential from the +! sounding, but this assures that the base state is in exact hydrostatic balance with +! respect to the model eqns. + + DO k = 2,kte + grid%phb(i,k,j) = grid%phb(i,k-1,j) - grid%dnw(k-1)*grid%mub(i,j)*grid%alb(i,k-1,j) + ENDDO + + ENDDO + ENDDO + + IF ( wrf_dm_on_monitor() ) THEN + write(6,*) ' ptop is ',grid%p_top + write(6,*) ' base state grid%mub(1,1), p_surf is ',grid%mub(1,1),grid%mub(1,1)+grid%p_top + ENDIF + +! calculate full state for each column - this includes moisture. + + write(6,*) ' getting moist sounding for full state ' + dry_sounding = .false. + CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in ) + + DO J = jts, min(jde-1,jte) + DO I = its, min(ide-1,ite) + +! At this point grid%p_top is already set. find the DRY mass in the column +! by interpolating the DRY pressure. + + pd_surf = interp_0( pd_in, zk, grid%phb(i,1,j)/g, nl_in ) + +! compute the perturbation mass and the full mass + + grid%mu_1(i,j) = pd_surf-grid%p_top - grid%mub(i,j) + grid%mu_2(i,j) = grid%mu_1(i,j) + grid%mu0(i,j) = grid%mu_1(i,j) + grid%mub(i,j) + +! given the dry pressure and coordinate system, interp the potential +! temperature and qv + + do k=1,kde-1 + + p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top + + moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in ) + grid%t_1(i,k,j) = interp_0( theta, pd_in, p_level, nl_in ) - t0 + grid%t_2(i,k,j) = grid%t_1(i,k,j) + + + enddo + +! integrate the hydrostatic equation (from the RHS of the bigstep +! vertical momentum equation) down from the top to get grid%p. +! first from the top of the model to the top pressure + + k = kte-1 ! top level + + qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV)) + qvf2 = 1./(1.+qvf1) + qvf1 = qvf1*qvf2 + +! grid%p(i,k,j) = - 0.5*grid%mu_1(i,j)/grid%rdnw(k) + grid%p(i,k,j) = - 0.5*(grid%mu_1(i,j)+qvf1*grid%mub(i,j))/grid%rdnw(k)/qvf2 + qvf = 1. + rvovrd*moist(i,k,j,P_QV) + grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* & + (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm) + grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j) + +! down the column + + do k=kte-2,1,-1 + qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV)) + qvf2 = 1./(1.+qvf1) + qvf1 = qvf1*qvf2 + grid%p(i,k,j) = grid%p(i,k+1,j) - (grid%mu_1(i,j) + qvf1*grid%mub(i,j))/qvf2/grid%rdn(k+1) + qvf = 1. + rvovrd*moist(i,k,j,P_QV) + grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* & + (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm) + grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j) + enddo + +! this is the hydrostatic equation used in the model after the +! small timesteps. In the model, grid%al (inverse density) +! is computed from the geopotential. + + + grid%ph_1(i,1,j) = 0. + DO k = 2,kte + grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*( & + (grid%mub(i,j)+grid%mu_1(i,j))*grid%al(i,k-1,j)+ & + grid%mu_1(i,j)*grid%alb(i,k-1,j) ) + + grid%ph_2(i,k,j) = grid%ph_1(i,k,j) + grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j) + ENDDO + + IF ( wrf_dm_on_monitor() ) THEN + if((i==2) .and. (j==2)) then + write(6,*) ' grid%ph_1 calc ',grid%ph_1(2,1,2),grid%ph_1(2,2,2),& + grid%mu_1(2,2)+grid%mub(2,2),grid%mu_1(2,2), & + grid%alb(2,1,2),grid%al(1,2,1),grid%rdnw(1) + endif + ENDIF + + ENDDO + ENDDO + +! checking if the perturbation (bubble) is to be applied + IF ((delt/=0.) .and. (x_rad > 0.) & + .and. (y_rad > 0.) & + .and. (z_rad > 0.)) THEN +! thermal perturbation to kick off convection + + write(6,*) ' nxc, nyc for perturbation ',nxc,nyc + write(6,('A23,f18.16')) ' delt for perturbation ',delt + write(6,('A30,f18.12')) ' x radius of the perturbation ' ,x_rad + write(6,('A30,f18.12')) ' y radius of the perturbation ' ,y_rad + write(6,('A30,f18.12')) ' z radius of the perturbation ' ,z_rad + write(6,('A30,f18.12')) ' height of the perturbation ' ,hght_pert + + DO J = jts, min(jde-1,jte) + yrad = config_flags%dy*float(j-nyc)/y_rad +! yrad = 0. + DO I = its, min(ide-1,ite) + xrad = config_flags%dx*float(i-nxc)/x_rad +! xrad = 0. + DO K = 1, kte-1 + +! put in preturbation theta (bubble) and recalc density. note, +! the mass in the column is not changing, so when theta changes, +! we recompute density and geopotential + + zrad = 0.5*(grid%ph_1(i,k,j)+grid%ph_1(i,k+1,j) & + +grid%phb(i,k,j)+grid%phb(i,k+1,j))/g + zrad = (zrad-hght_pert)/z_rad + RAD=SQRT(xrad*xrad+yrad*yrad+zrad*zrad) + IF(RAD <= 1.) THEN + grid%t_1(i,k,j)=grid%t_1(i,k,j)+delt*COS(.5*PI*RAD)**2 + grid%t_2(i,k,j)=grid%t_1(i,k,j) + qvf = 1. + rvovrd*moist(i,k,j,P_QV) + grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* & + (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm) + grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j) + ENDIF + ENDDO + +! rebalance hydrostatically + + DO k = 2,kte + grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*( & + (grid%mub(i,j)+grid%mu_1(i,j))*grid%al(i,k-1,j)+ & + grid%mu_1(i,j)*grid%alb(i,k-1,j) ) + + grid%ph_2(i,k,j) = grid%ph_1(i,k,j) + grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j) + ENDDO + + ENDDO + ENDDO + +!End of setting up the perturbation (bubble) +ENDIF + + IF ( wrf_dm_on_monitor() ) THEN + write(6,*) ' grid%mu_1 from comp ', grid%mu_1(1,1) + write(6,*) ' full state sounding from comp, ph, grid%p, grid%al, grid%t_1, qv ' + do k=1,kde-1 + write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1)+grid%phb(1,k,1), & + grid%p(1,k,1)+grid%pb(1,k,1), grid%alt(1,k,1), & + grid%t_1(1,k,1)+t0, moist(1,k,1,P_QV) + enddo + + write(6,*) ' pert state sounding from comp, grid%ph_1, pp, alp, grid%t_1, qv ' + do k=1,kde-1 + write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1), & + grid%p(1,k,1), grid%al(1,k,1), & + grid%t_1(1,k,1), moist(1,k,1,P_QV) + enddo + ENDIF + +! interp v + + DO J = jts, jte + DO I = its, min(ide-1,ite) + + IF (j == jds) THEN + z_at_v = grid%phb(i,1,j)/g + ELSE IF (j == jde) THEN + z_at_v = grid%phb(i,1,j-1)/g + ELSE + z_at_v = 0.5*(grid%phb(i,1,j)+grid%phb(i,1,j-1))/g + END IF + p_surf = interp_0( p_in, zk, z_at_v, nl_in ) + + DO K = 1, kte-1 + p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top + grid%v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in ) + grid%v_2(i,k,j) = grid%v_1(i,k,j) + ENDDO + + ENDDO + ENDDO + +! interp u + + DO J = jts, min(jde-1,jte) + DO I = its, ite + + IF (i == ids) THEN + z_at_u = grid%phb(i,1,j)/g + ELSE IF (i == ide) THEN + z_at_u = grid%phb(i-1,1,j)/g + ELSE + z_at_u = 0.5*(grid%phb(i,1,j)+grid%phb(i-1,1,j))/g + END IF + + p_surf = interp_0( p_in, zk, z_at_u, nl_in ) + + DO K = 1, kte-1 + p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top + grid%u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in ) + grid%u_2(i,k,j) = grid%u_1(i,k,j) + ENDDO + + ENDDO + ENDDO + +! set w + + DO J = jts, min(jde-1,jte) + DO K = kts, kte + DO I = its, min(ide-1,ite) + grid%w_1(i,k,j) = 0. + grid%w_2(i,k,j) = 0. + ENDDO + ENDDO + ENDDO + +! set a few more things + + DO J = jts, min(jde-1,jte) + DO K = kts, kte-1 + DO I = its, min(ide-1,ite) + grid%h_diabatic(i,k,j) = 0. + ENDDO + ENDDO + ENDDO + + IF ( wrf_dm_on_monitor() ) THEN + DO k=1,kte-1 + grid%t_base(k) = grid%t_1(1,k,1) + grid%qv_base(k) = moist(1,k,1,P_QV) + grid%u_base(k) = grid%u_1(1,k,1) + grid%v_base(k) = grid%v_1(1,k,1) + grid%z_base(k) = 0.5*(grid%phb(1,k,1)+grid%phb(1,k+1,1)+grid%ph_1(1,k,1)+grid%ph_1(1,k+1,1))/g + ENDDO + ENDIF + CALL wrf_dm_bcast_real( grid%t_base , kte ) + CALL wrf_dm_bcast_real( grid%qv_base , kte ) + CALL wrf_dm_bcast_real( grid%u_base , kte ) + CALL wrf_dm_bcast_real( grid%v_base , kte ) + CALL wrf_dm_bcast_real( grid%z_base , kte ) + + DO J = jts, min(jde-1,jte) + DO I = its, min(ide-1,ite) + thtmp = grid%t_2(i,1,j)+t0 + ptmp = grid%p(i,1,j)+grid%pb(i,1,j) + temp(1) = thtmp * (ptmp/p1000mb)**rcp + thtmp = grid%t_2(i,2,j)+t0 + ptmp = grid%p(i,2,j)+grid%pb(i,2,j) + temp(2) = thtmp * (ptmp/p1000mb)**rcp + thtmp = grid%t_2(i,3,j)+t0 + ptmp = grid%p(i,3,j)+grid%pb(i,3,j) + temp(3) = thtmp * (ptmp/p1000mb)**rcp + + grid%tsk(I,J)=grid%cf1*temp(1)+grid%cf2*temp(2)+grid%cf3*temp(3) + grid%tmn(I,J)=grid%tsk(I,J)-0.5 + ENDDO + ENDDO + + END SUBROUTINE init_domain_rk + + SUBROUTINE init_module_initialize + END SUBROUTINE init_module_initialize + +!--------------------------------------------------------------------- + +! test driver for get_sounding +! +! implicit none +! integer n +! parameter(n = 1000) +! real zk(n),p(n),theta(n),rho(n),u(n),v(n),qv(n),pd(n) +! logical dry +! integer nl,k +! +! dry = .false. +! dry = .true. +! call get_sounding( zk, p, pd, theta, rho, u, v, qv, dry, n, nl ) +! write(6,*) ' input levels ',nl +! write(6,*) ' sounding ' +! write(6,*) ' k height(m) press (Pa) pd(Pa) theta (K) den(kg/m^3) u(m/s) v(m/s) qv(g/g) ' +! do k=1,nl +! write(6,'(1x,i3,8(1x,1pe10.3))') k, zk(k), p(k), pd(k), theta(k), rho(k), u(k), v(k), qv(k) +! enddo +! end +! +!--------------------------------------------------------------------------- + + subroutine get_sounding( zk, p, p_dry, theta, rho, & + u, v, qv, dry, nl_max, nl_in ) + implicit none + + integer nl_max, nl_in + real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), & + u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max) + logical dry + + integer n + parameter(n=1000) + logical debug + parameter( debug = .true.) + +! input sounding data + + real p_surf, th_surf, qv_surf + real pi_surf, pi(n) + real h_input(n), th_input(n), qv_input(n), u_input(n), v_input(n) + +! diagnostics + + real rho_surf, p_input(n), rho_input(n) + real pm_input(n) ! this are for full moist sounding + +! local data + + real r + parameter (r = r_d) + integer k, it, nl + real qvf, qvf1, dz + +! first, read the sounding + + call read_sounding( p_surf, th_surf, qv_surf, & + h_input, th_input, qv_input, u_input, v_input,n, nl, debug ) + + if(dry) then + do k=1,nl + qv_input(k) = 0. + enddo + endif + + if(debug) write(6,*) ' number of input levels = ',nl + + nl_in = nl + if(nl_in .gt. nl_max ) then + write(6,*) ' too many levels for input arrays ',nl_in,nl_max + call wrf_error_fatal ( ' too many levels for input arrays ' ) + end if + +! compute diagnostics, +! first, convert qv(g/kg) to qv(g/g) + + do k=1,nl + qv_input(k) = 0.001*qv_input(k) + enddo + + p_surf = 100.*p_surf ! convert to pascals + qvf = 1. + rvovrd*qv_input(1) + rho_surf = 1./((r/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm)) + pi_surf = (p_surf/p1000mb)**(r/cp) + + if(debug) then + write(6,*) ' surface density is ',rho_surf + write(6,*) ' surface pi is ',pi_surf + end if + + +! integrate moist sounding hydrostatically, starting from the +! specified surface pressure +! -> first, integrate from surface to lowest level + + qvf = 1. + rvovrd*qv_input(1) + qvf1 = 1. + qv_input(1) + rho_input(1) = rho_surf + dz = h_input(1) + do it=1,10 + pm_input(1) = p_surf & + - 0.5*dz*(rho_surf+rho_input(1))*g*qvf1 + rho_input(1) = 1./((r/p1000mb)*th_input(1)*qvf*((pm_input(1)/p1000mb)**cvpm)) + enddo + +! integrate up the column + + do k=2,nl + rho_input(k) = rho_input(k-1) + dz = h_input(k)-h_input(k-1) + qvf1 = 0.5*(2.+(qv_input(k-1)+qv_input(k))) + qvf = 1. + rvovrd*qv_input(k) ! qv is in g/kg here + + do it=1,10 + pm_input(k) = pm_input(k-1) & + - 0.5*dz*(rho_input(k)+rho_input(k-1))*g*qvf1 + rho_input(k) = 1./((r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm)) + enddo + enddo + +! we have the moist sounding + +! next, compute the dry sounding using p at the highest level from the +! moist sounding and integrating down. + + p_input(nl) = pm_input(nl) + + do k=nl-1,1,-1 + dz = h_input(k+1)-h_input(k) + p_input(k) = p_input(k+1) + 0.5*dz*(rho_input(k)+rho_input(k+1))*g + enddo + + + do k=1,nl + + zk(k) = h_input(k) + p(k) = pm_input(k) + p_dry(k) = p_input(k) + theta(k) = th_input(k) + rho(k) = rho_input(k) + u(k) = u_input(k) + v(k) = v_input(k) + qv(k) = qv_input(k) + + enddo + + if(debug) then + write(6,*) ' sounding ' + write(6,*) ' k height(m) press (Pa) pd(Pa) theta (K) den(kg/m^3) u(m/s) v(m/s) qv(g/g) ' + do k=1,nl + write(6,'(1x,i3,8(1x,1pe10.3))') k, zk(k), p(k), p_dry(k), theta(k), rho(k), u(k), v(k), qv(k) + enddo + + end if + + end subroutine get_sounding + +!------------------------------------------------------- + + subroutine read_sounding( ps,ts,qvs,h,th,qv,u,v,n,nl,debug ) + implicit none + integer n,nl + real ps,ts,qvs,h(n),th(n),qv(n),u(n),v(n) + logical end_of_file + logical debug + + integer k + + open(unit=10,file='input_sounding',form='formatted',status='old') + rewind(10) + read(10,*) ps, ts, qvs + if(debug) then + write(6,*) ' input sounding surface parameters ' + write(6,*) ' surface pressure (mb) ',ps + write(6,*) ' surface pot. temp (K) ',ts + write(6,*) ' surface mixing ratio (g/kg) ',qvs + end if + + end_of_file = .false. + k = 0 + + do while (.not. end_of_file) + + read(10,*,end=100) h(k+1), th(k+1), qv(k+1), u(k+1), v(k+1) + k = k+1 + if(debug) write(6,'(1x,i3,5(1x,e10.3))') k, h(k), th(k), qv(k), u(k), v(k) + go to 110 + 100 end_of_file = .true. + 110 continue + enddo + + nl = k + + close(unit=10,status = 'keep') + + end subroutine read_sounding + +END MODULE module_initialize_ideal diff --git a/wrfv2_fire/test/em_fire/small/namelist.input_hypgrid b/wrfv2_fire/test/em_fire/small/namelist.input_hypgrid new file mode 100644 index 00000000..a1bd645e --- /dev/null +++ b/wrfv2_fire/test/em_fire/small/namelist.input_hypgrid @@ -0,0 +1,181 @@ + &time_control + run_days = 0, + run_hours = 0, + run_minutes = 5, + run_seconds = 0, + start_year = 0001, 0001, 0001, + start_month = 01, 01, 01, + start_day = 01, 01, 01, + start_hour = 00, 00, 00, + start_minute = 00, 00, 00, + start_second = 00, 00, 00, + end_year = 0001, 0001, 0001, + end_month = 01, 01, 01, + end_day = 01, 01, 01, + end_hour = 00, 00, 00, + end_minute = 600, 600, 600, + end_second = 00, 00, 00, + history_interval = 5, 30, 30, + frames_per_outfile = 1000, 1000, 1000, + restart = .false., + restart_interval = 1 + io_form_history = 2 + io_form_restart = 2 + io_form_input = 2 + io_form_boundary = 2 + debug_level = 101 + / + + &domains + time_step = 0, + !time_step = 5, + time_step_fract_num = 25, + time_step_fract_den = 100, + max_dom = 1, + s_we = 1, 1, 1, + e_we = 42, 43, 43, + s_sn = 1, 1, 1, + e_sn = 42, 43, 43, + s_vert = 1, 1, 1, + e_vert = 41, 41, 41, + dx = 60, 30, 10, + dy = 60, 30, 10, + ztop = 1500, 1500, 1500, + grid_id = 1, 2, 3, + parent_id = 0, 1, 2, + i_parent_start = 0, 1, 1, + j_parent_start = 0, 1, 1, + parent_grid_ratio = 1, 2, 3, + parent_time_step_ratio = 1, 2, 3, + feedback = 1, + smooth_option = 0 + sr_x = 10, 0, 0 + sr_y = 10, 0, 0 + / + + &physics + mp_physics = 2, 0, 0, + ra_lw_physics = 1, 0, 0, + ra_sw_physics = 1, 0, 0, + radt = 30, 30, 30, + sf_sfclay_physics = 1, 0, 0, + sf_surface_physics = 1, 0, 0, + bl_pbl_physics = 0, 0, 0, + bldt = 0, 0, 0, + cu_physics = 0, 0, 0, + cudt = 0, 0, 0, + isfflx = 1, + ifsnow = 0, + icloud = 0, + num_soil_layers = 5, + mp_zero_out = 0, + / + + &fdda + / + + &dynamics + rk_ord = 3, + diff_opt = 2, + km_opt = 2, + damp_opt = 0, + zdamp = 5000., 5000., 5000., + dampcoef = 0.2, 0.2, 0.2 + khdif = 0.05, 0.05, 0.05, + kvdif = 0.05, 0.05, 0.05, + smdiv = 0.1, 0.1, 0.1, + emdiv = 0.01, 0.01, 0.01, + epssm = 0.1, 0.1, 0.1 + mix_full_fields = .true., .true., .true., + non_hydrostatic = .true., .true., .true., + h_mom_adv_order = 5, 5, 5, + v_mom_adv_order = 3, 3, 3, + h_sca_adv_order = 5, 5, 5, + v_sca_adv_order = 3, 3, 3, + time_step_sound = 20, 20, 20, + moist_adv_opt = 1, 1, 1, + scalar_adv_opt = 1, 1, 1, + / + + &bdy_control + periodic_x = .false.,.false.,.false., + symmetric_xs = .false.,.false.,.false., + symmetric_xe = .false.,.false.,.false., + open_xs = .true., .false.,.false., + open_xe = .true., .false.,.false., + periodic_y = .false.,.false.,.false., + symmetric_ys = .false.,.false.,.false., + symmetric_ye = .false.,.false.,.false., + open_ys = .true., .false.,.false., + open_ye = .true., .false.,.false., + nested = .false., .true., .true., + / + + &grib2 + / + + &namelist_quilt + nio_tasks_per_group = 0, + nio_groups = 1, + / + + &fire ! be sure to set sr_x,sr_y in domains-namelist (to set refinement in x,y) + ifire = 2, ! integer, = 0: no fire, 2=turn on fire model + fire_fuel_read = 0, ! integer, -1: from WPS, 0= use fire_fuel_cat, 1= by altitude + fire_fuel_cat = 3, ! integer, if specified which fuel category? +! ignition + fire_num_ignitions = 3, ! integer, only the first fire_num_ignition used, up to 5 allowed + fire_ignition_start_x1 = 1000, ! start points of ignition lines, in m from lower left corner + fire_ignition_start_y1 = 500, ! start points of ignition lines, in m from lower left corner + fire_ignition_end_x1 = 1000, ! end points of ignition lines, in m from lower left corner + fire_ignition_end_y1 = 1900, ! end points of ignition lines, in m from lower left corner + fire_ignition_radius1 = 18, ! all within this radius will ignite, > fire mesh step + fire_ignition_time1 = 2, ! sec for ignition from the start + fire_ignition_start_x2 = 1500, ! start points of ignition lines, in m from lower left corner + fire_ignition_start_y2 = 500, ! start points of ignition lines, in m from lower left corner + fire_ignition_end_x2 = 1500, ! end points of ignition lines, in m from lower left corner + fire_ignition_end_y2 = 1900, ! end points of ignition lines, in m from lower left corner + fire_ignition_radius2 = 18, ! all within this radius will ignite, > fire mesh step + fire_ignition_time2 = 3, ! sec for ignition from the start! end ignition for sfire + fire_ignition_start_x3 = 1400, ! start points of ignition lines, in m from lower left corner + fire_ignition_start_y3 = 1400, ! start points of ignition lines, in m from lower left corner + fire_ignition_end_x3 = 1400, ! end points of ignition lines, in m from lower left corner + fire_ignition_end_y3 = 1400, ! end points of ignition lines, in m from lower left corner + fire_ignition_radius3 = 50, ! all within this radius will ignite, > fire mesh step + fire_ignition_time3 = 4, ! sec for ignition from the start! end ignition for sfire +! +! verbosity + fire_print_msg = 1, ! 1 print fire debugging messages + fire_print_file = 0, ! 1 write files for matlab +! +! experiments +! + stretch_grd = .true., + stretch_hyp = .false., + fire_const_time = -1., ! (s) if >0, time from start to stop fire evolution and keep heat output constant + fire_const_grnhfx = -1, ! (W/s) if both >=0, use this flux (meant to be used when fire_const_time=ignition time) + fire_const_grnqfx = -1, ! (W/s) if both >=0, use this flux (meant to be used when fire_const_time=ignition time) + fire_test_steps=0, ! >0 = on first call, do specified number of steps and terminate (testing only) + fire_mountain_type=0, ! in ideal: 0=none, 1= hill, 2=EW ridge, 3=NS ridge + fire_mountain_height=500., ! (m) ideal mountain height + fire_mountain_start_x=1000., ! (m) coord of start of the mountain from lower left corder (just like ignition) + fire_mountain_start_y=1100., ! (m) coord of start of the mountain from lower left corder (just like ignition) + fire_mountain_end_x=1500., ! (m) coord of end of the mountain from lower left corder (just like ignition) + fire_mountain_end_y=1400., ! (m) coord of end of the mountain from lower left corder (just like ignition) + fire_topo_from_atm=1, ! 0 = fire mesh topo set from fine-res data, 1 = populate by interpolating from atmosphere + +! +! method switches for developers only, do not change! +! + fire_boundary_guard = -1, ! integer, number of cells to stop when fire close to the domain boundary, -1 turn off + fire_fuel_left_irl=2, ! refinement to integrate fuel_left, must be even + fire_fuel_left_jrl=2, ! refinement to integrate fuel_left, must be even + fire_atm_feedback=1.0, ! real, multiplier for heat fluxes, 1.=normal, 0.=turn off two-way coupling + fire_back_weight=0.5, ! RK timestepping coefficient, 0=forward, 0.5=Heun + fire_grows_only=1, ! if >0 level set function cannot increase = fire can only grow + fire_viscosity=0.4, ! artificial viscosity in level set method (max 1, needed with fire_upwinding=0) + fire_upwinding=3, ! 0=none, 1=standard, 2=godunov, 3=eno, 4=sethian + fire_fuel_left_method=1, ! for now, use 1 only + fire_lfn_ext_up=1.0, ! 0.=extend level set function at boundary by reflection, 1.=always up + fire_advection=0, ! 0 = cawfe, 1 = use abs speed/slope in spread rate, then project on normal to fireline +/