/wrfv2_fire/phys/module_sf_urban.F
FORTRAN Legacy | 2421 lines | 1481 code | 375 blank | 565 comment | 237 complexity | da74e6c6f59556f308ac988182e383d8 MD5 | raw file
Possible License(s): AGPL-1.0
- MODULE module_sf_urban
- !===============================================================================
- ! Single-Layer Urban Canopy Model for WRF Noah-LSM
- ! Original Version: 2002/11/06 by Hiroyuki Kusaka
- ! Last Update: 2006/08/24 by Fei Chen and Mukul Tewari (NCAR/RAL)
- !===============================================================================
- CHARACTER(LEN=4) :: LU_DATA_TYPE
- INTEGER :: ICATE
- REAL, ALLOCATABLE, DIMENSION(:) :: ZR_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: Z0C_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: Z0HC_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: ZDC_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: SVF_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: R_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: RW_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: HGT_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: AH_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: BETR_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: BETB_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: BETG_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: FRC_URB_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: COP_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: PWIN_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: BETA_TBL
- INTEGER, ALLOCATABLE, DIMENSION(:) :: SW_COND_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: TIME_ON_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: TIME_OFF_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: TARGTEMP_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: GAPTEMP_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: TARGHUM_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: GAPHUM_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: PERFLO_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: HSESF_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: CAPR_TBL, CAPB_TBL, CAPG_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: AKSR_TBL, AKSB_TBL, AKSG_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: ALBR_TBL, ALBB_TBL, ALBG_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: EPSR_TBL, EPSB_TBL, EPSG_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: Z0R_TBL, Z0B_TBL, Z0G_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: SIGMA_ZED_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: Z0HB_TBL, Z0HG_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: TRLEND_TBL, TBLEND_TBL, TGLEND_TBL
- REAL, ALLOCATABLE, DIMENSION(:) :: AKANDA_URBAN_TBL
- !for BEP
- ! MAXDIRS :: The maximum number of street directions we're allowed to define
- INTEGER, PARAMETER :: MAXDIRS = 3
- ! MAXHGTS :: The maximum number of building height bins we're allowed to define
- INTEGER, PARAMETER :: MAXHGTS = 50
- INTEGER, ALLOCATABLE, DIMENSION(:) :: NUMDIR_TBL
- REAL, ALLOCATABLE, DIMENSION(:,:) :: STREET_DIRECTION_TBL
- REAL, ALLOCATABLE, DIMENSION(:,:) :: STREET_WIDTH_TBL
- REAL, ALLOCATABLE, DIMENSION(:,:) :: BUILDING_WIDTH_TBL
- INTEGER, ALLOCATABLE, DIMENSION(:) :: NUMHGT_TBL
- REAL, ALLOCATABLE, DIMENSION(:,:) :: HEIGHT_BIN_TBL
- REAL, ALLOCATABLE, DIMENSION(:,:) :: HPERCENT_BIN_TBL
- !end BEP
- INTEGER :: BOUNDR_DATA,BOUNDB_DATA,BOUNDG_DATA
- INTEGER :: CH_SCHEME_DATA, TS_SCHEME_DATA
- INTEGER :: ahoption ! Miao, 2007/01/17, cal. ah
- REAL, DIMENSION(1:24) :: ahdiuprf ! ah diurnal profile, tloc: 1-24
- REAL, DIMENSION(1:24) :: hsequip_tbl
- INTEGER :: allocate_status
- ! INTEGER :: num_roof_layers
- ! INTEGER :: num_wall_layers
- ! INTEGER :: num_road_layers
- CONTAINS
- !===============================================================================
- !
- ! Author:
- ! Hiroyuki KUSAKA, PhD
- ! University of Tsukuba, JAPAN
- ! (CRIEPI, NCAR/MMM visiting scientist, 2002-2004)
- ! kusaka@ccs.tsukuba.ac.jp
- !
- ! Co-Researchers:
- ! Fei CHEN, PhD
- ! NCAR/RAP feichen@ucar.edu
- ! Mukul TEWARI, PhD
- ! NCAR/RAP mukul@ucar.edu
- !
- ! Purpose:
- ! Calculate surface temeprature, fluxes, canopy air temperature, and canopy wind
- !
- ! Subroutines:
- ! module_sf_urban
- ! |- urban
- ! |- read_param
- ! |- mos or jurges
- ! |- multi_layer or force_restore
- ! |- urban_param_init <-- URBPARM.TBL
- ! |- urban_var_init
- !
- ! Input Data from WRF [MKS unit]:
- !
- ! UTYPE [-] : Urban type. 1=Commercial/Industrial; 2=High-intensity residential;
- ! : 3=low-intensity residential
- ! TA [K] : Potential temperature at 1st wrf level (absolute temp)
- ! QA [kg/kg] : Mixing ratio at 1st atmospheric level
- ! UA [m/s] : Wind speed at 1st atmospheric level
- ! SSG [W/m/m] : Short wave downward radiation at a flat surface
- ! Note this is the total of direct and diffusive solar
- ! downward radiation. If without two components, the
- ! single solar downward can be used instead.
- ! SSG = SSGD + SSGQ
- ! LSOLAR [-] : Indicating the input type of solar downward radiation
- ! True: both direct and diffusive solar radiation
- ! are available
- ! False: only total downward ridiation is available.
- ! SSGD [W/m/m] : Direct solar radiation at a flat surface
- ! if SSGD is not available, one can assume a ratio SRATIO
- ! (e.g., 0.7), so that SSGD = SRATIO*SSG
- ! SSGQ [W/m/m] : Diffuse solar radiation at a flat surface
- ! If SSGQ is not available, SSGQ = SSG - SSGD
- ! LLG [W/m/m] : Long wave downward radiation at a flat surface
- ! RAIN [mm/h] : Precipitation
- ! RHOO [kg/m/m/m] : Air density
- ! ZA [m] : First atmospheric level
- ! as a lowest boundary condition
- ! DECLIN [rad] : solar declination
- ! COSZ : = sin(fai)*sin(del)+cos(fai)*cos(del)*cos(omg)
- ! OMG [rad] : solar hour angle
- ! XLAT [deg] : latitude
- ! DELT [sec] : Time step
- ! ZNT [m] : Roughnes length
- !
- ! Output Data to WRF [MKS unit]:
- !
- ! TS [K] : Surface potential temperature (absolute temp)
- ! QS [-] : Surface humidity
- !
- ! SH [W/m/m/] : Sensible heat flux, = FLXTH*RHOO*CPP
- ! LH [W/m/m] : Latent heat flux, = FLXHUM*RHOO*ELL
- ! LH_INEMATIC [kg/m/m/sec]: Moisture Kinematic flux, = FLXHUM*RHOO
- ! SW [W/m/m] : Upward shortwave radiation flux,
- ! = SSG-SNET*697.7*60. (697.7*60.=100.*100.*4.186)
- ! ALB [-] : Time-varying albedo
- ! LW [W/m/m] : Upward longwave radiation flux,
- ! = LNET*697.7*60.-LLG
- ! G [W/m/m] : Heat Flux into the Ground
- ! RN [W/m/m] : Net radiation
- !
- ! PSIM [-] : Diagnostic similarity stability function for momentum
- ! PSIH [-] : Diagnostic similarity stability function for heat
- !
- ! TC [K] : Diagnostic canopy air temperature
- ! QC [-] : Diagnostic canopy humidity
- !
- ! TH2 [K] : Diagnostic potential temperature at 2 m
- ! Q2 [-] : Diagnostic humidity at 2 m
- ! U10 [m/s] : Diagnostic u wind component at 10 m
- ! V10 [m/s] : Diagnostic v wind component at 10 m
- !
- ! CHS, CHS2 [m/s] : CH*U at ZA, CH*U at 2 m (not used)
- !
- ! Important parameters:
- !
- ! Morphology of the urban canyon:
- ! These parameters assigned in the URBPARM.TBL
- !
- ! ZR [m] : roof level (building height)
- ! SIGMA_ZED [m] : Standard Deviation of roof height
- ! ROOF_WIDTH [m] : roof (i.e., building) width
- ! ROAD_WIDTH [m] : road width
- !
- ! Parameters derived from the morphological terms above.
- ! These parameters are computed in the code.
- !
- ! HGT [-] : normalized building height
- ! SVF [-] : sky view factor
- ! R [-] : Normalized roof width (a.k.a. "building coverage ratio")
- ! RW [-] : = 1 - R
- ! Z0C [m] : Roughness length above canyon for momentum (1/10 of ZR)
- ! Z0HC [m] : Roughness length above canyon for heat (1/10 of Z0C)
- ! ZDC [m] : Zero plane displacement height (1/5 of ZR)
- !
- ! Following parameter are assigned in run/URBPARM.TBL
- !
- ! AH [ W m{-2} ] : anthropogenic heat ( W m{-2} in the table, converted internally to cal cm{-2} )
- ! CAPR[ J m{-3} K{-1} ] : heat capacity of roof ( units converted in code to [ cal cm{-3} deg{-1} ] )
- ! CAPB[ J m{-3} K{-1} ] : heat capacity of building wall ( units converted in code to [ cal cm{-3} deg{-1} ] )
- ! CAPG[ J m{-3} K{-1} ] : heat capacity of road ( units converted in code to [ cal cm{-3} deg{-1} ] )
- ! AKSR [ J m{-1} s{-1} K{-1} ] : thermal conductivity of roof ( units converted in code to [ cal cm{-1} s{-1} deg{-1} ] )
- ! AKSB [ J m{-1} s{-1} K{-1} ] : thermal conductivity of building wall ( units converted in code to [ cal cm{-1} s{-1} deg{-1} ] )
- ! AKSG [ J m{-1} s{-1} K{-1} ] : thermal conductivity of road ( units converted in code to [ cal cm{-1} s{-1} deg{-1} ] )
- ! ALBR [-] : surface albedo of roof
- ! ALBB [-] : surface albedo of building wall
- ! ALBG [-] : surface albedo of road
- ! EPSR [-] : surface emissivity of roof
- ! EPSB [-] : surface emissivity of building wall
- ! EPSG [-] : surface emissivity of road
- ! Z0B [m] : roughness length for momentum of building wall (only for CH_SCHEME = 1)
- ! Z0G [m] : roughness length for momentum of road (only for CH_SCHEME = 1)
- ! Z0HB [m] : roughness length for heat of building wall (only for CH_SCHEME = 1)
- ! Z0HG [m] : roughness length for heat of road
- ! num_roof_layers : number of layers within roof
- ! num_wall_layers : number of layers within building walls
- ! num_road_layers : number of layers within below road surface
- ! NOTE: for now, these layers are defined as same as the number of soil layers in namelist.input
- ! DZR [cm] : thickness of each roof layer
- ! DZB [cm] : thickness of each building wall layer
- ! DZG [cm] : thickness of each ground layer
- ! BOUNDR [integer 1 or 2] : Boundary Condition for Roof Layer Temp [1: Zero-Flux, 2: T = Constant]
- ! BOUNDB [integer 1 or 2] : Boundary Condition for Building Wall Layer Temp [1: Zero-Flux, 2: T = Constant]
- ! BOUNDG [integer 1 or 2] : Boundary Condition for Road Layer Temp [1: Zero-Flux, 2: T = Constant]
- ! TRLEND [K] : lower boundary condition of roof temperature
- ! TBLEND [K] : lower boundary condition of building temperature
- ! TGLEND [K] : lower boundary condition of ground temperature
- ! CH_SCHEME [integer 1 or 2] : Sfc exchange scheme used for building wall and road
- ! [1: M-O Similarity Theory, 2: Empirical Form (recommend)]
- ! TS_SCHEME [integer 1 or 2] : Scheme for computing surface temperature (for roof, wall, and road)
- ! [1: 4-layer model, 2: Force-Restore method]
- !
- !for BEP
- ! numdir [ - ] : Number of street directions defined for a particular urban category
- ! street_direction [ deg ] : Direction of streets for a particular urban category and a particular street direction
- ! street_width [ m ] : Width of street for a particular urban category and a particular street direction
- ! building_width [ m ] : Width of buildings for a particular urban category and a particular street direction
- ! numhgt [ - ] : Number of building height levels defined for a particular urban category
- ! height_bin [ m ] : Building height bins defined for a particular urban category.
- ! hpercent_bin [ % ] : Percentage of a particular urban category populated by buildings of particular height bins
- !end BEP
- ! Moved from URBPARM.TBL
- !
- ! BETR [-] : minimum moisture availability of roof
- ! BETB [-] : minimum moisture availability of building wall
- ! BETG [-] : minimum moisture availability of road
- ! Z0R [m] : roughness length for momentum of roof
- ! Z0HB [m] : roughness length for heat of building wall (only for CH_SCHEME = 1)
- ! Z0HG [m] : roughness length for heat of road
- ! num_roof_layers : number of layers within roof
- ! num_wall_layers : number of layers within building walls
- ! num_road_layers : number of layers within below road surface
- ! NOTE: for now, these layers are defined as same as the number of soil layers in namelist.input
- !
- ! References:
- ! Kusaka and Kimura (2004) J.Appl.Meteor., vol.43, p1899-1910
- ! Kusaka and Kimura (2004) J.Meteor.Soc.Japan, vol.82, p45-65
- ! Kusaka et al. (2001) Bound.-Layer Meteor., vol.101, p329-358
- !
- ! History:
- ! 2006/06 modified by H. Kusaka (Univ. Tsukuba), M. Tewari
- ! 2005/10/26, modified by Fei Chen, Mukul Tewari
- ! 2003/07/21 WRF , modified by H. Kusaka of CRIEPI (NCAR/MMM)
- ! 2001/08/26 PhD , modified by H. Kusaka of CRIEPI (Univ.Tsukuba)
- ! 1999/08/25 LCM , developed by H. Kusaka of CRIEPI (Univ.Tsukuba)
- !
- !===============================================================================
- !
- ! subroutine urban:
- !
- !===============================================================================
- SUBROUTINE urban(LSOLAR, & ! L
- num_roof_layers,num_wall_layers,num_road_layers, & ! I
- DZR,DZB,DZG, & ! I
- UTYPE,TA,QA,UA,U1,V1,SSG,SSGD,SSGQ,LLG,RAIN,RHOO, & ! I
- ZA,DECLIN,COSZ,OMG,XLAT,DELT,ZNT, & ! I
- CHS, CHS2, & ! I
- TR, TB, TG, TC, QC, UC, & ! H
- TRL,TBL,TGL, & ! H
- XXXR, XXXB, XXXG, XXXC, & ! H
- TS,QS,SH,LH,LH_KINEMATIC, & ! O
- SW,ALB,LW,G,RN,PSIM,PSIH, & ! O
- GZ1OZ0, & ! O
- CMR_URB,CHR_URB,CMC_URB,CHC_URB, & ! I/O
- U10,V10,TH2,Q2,UST & ! O
- )
- IMPLICIT NONE
- REAL, PARAMETER :: CP=0.24 ! heat capacity of dry air [cgs unit]
- REAL, PARAMETER :: EL=583. ! latent heat of vaporation [cgs unit]
- REAL, PARAMETER :: SIG=8.17E-11 ! stefun bolzman constant [cgs unit]
- REAL, PARAMETER :: SIG_SI=5.67E-8 ! [MKS unit]
- REAL, PARAMETER :: AK=0.4 ! kalman const. [-]
- REAL, PARAMETER :: PI=3.14159 ! pi [-]
- REAL, PARAMETER :: TETENA=7.5 ! const. of Tetens Equation [-]
- REAL, PARAMETER :: TETENB=237.3 ! const. of Tetens Equation [-]
- REAL, PARAMETER :: SRATIO=0.75 ! ratio between direct/total solar [-]
- REAL, PARAMETER :: CPP=1004.5 ! heat capacity of dry air [J/K/kg]
- REAL, PARAMETER :: ELL=2.442E+06 ! latent heat of vaporization [J/kg]
- REAL, PARAMETER :: XKA=2.4E-5
- !-------------------------------------------------------------------------------
- ! C: configuration variables
- !-------------------------------------------------------------------------------
- LOGICAL, INTENT(IN) :: LSOLAR ! logical [true=both, false=SSG only]
- ! The following variables are also model configuration variables, but are
- ! defined in the URBAN.TBL and in the contains statement in the top of
- ! the module_urban_init, so we should not declare them here.
- INTEGER, INTENT(IN) :: num_roof_layers
- INTEGER, INTENT(IN) :: num_wall_layers
- INTEGER, INTENT(IN) :: num_road_layers
- REAL, INTENT(IN), DIMENSION(1:num_roof_layers) :: DZR ! grid interval of roof layers [cm]
- REAL, INTENT(IN), DIMENSION(1:num_wall_layers) :: DZB ! grid interval of wall layers [cm]
- REAL, INTENT(IN), DIMENSION(1:num_road_layers) :: DZG ! grid interval of road layers [cm]
- !-------------------------------------------------------------------------------
- ! I: input variables from LSM to Urban
- !-------------------------------------------------------------------------------
- INTEGER, INTENT(IN) :: UTYPE ! urban type [1=Commercial/Industrial, 2=High-intensity residential,
- ! 3=low-intensity residential]
- REAL, INTENT(IN) :: TA ! potential temp at 1st atmospheric level [K]
- REAL, INTENT(IN) :: QA ! mixing ratio at 1st atmospheric level [kg/kg]
- REAL, INTENT(IN) :: UA ! wind speed at 1st atmospheric level [m/s]
- REAL, INTENT(IN) :: U1 ! u at 1st atmospheric level [m/s]
- REAL, INTENT(IN) :: V1 ! v at 1st atmospheric level [m/s]
- REAL, INTENT(IN) :: SSG ! downward total short wave radiation [W/m/m]
- REAL, INTENT(IN) :: LLG ! downward long wave radiation [W/m/m]
- REAL, INTENT(IN) :: RAIN ! precipitation [mm/h]
- REAL, INTENT(IN) :: RHOO ! air density [kg/m^3]
- REAL, INTENT(IN) :: ZA ! first atmospheric level [m]
- REAL, INTENT(IN) :: DECLIN ! solar declination [rad]
- REAL, INTENT(IN) :: COSZ ! sin(fai)*sin(del)+cos(fai)*cos(del)*cos(omg)
- REAL, INTENT(IN) :: OMG ! solar hour angle [rad]
- REAL, INTENT(IN) :: XLAT ! latitude [deg]
- REAL, INTENT(IN) :: DELT ! time step [s]
- REAL, INTENT(IN) :: ZNT ! roughness length [m]
- REAL, INTENT(IN) :: CHS,CHS2 ! CH*U at za and 2 m [m/s]
- REAL, INTENT(INOUT) :: SSGD ! downward direct short wave radiation [W/m/m]
- REAL, INTENT(INOUT) :: SSGQ ! downward diffuse short wave radiation [W/m/m]
- REAL, INTENT(INOUT) :: CMR_URB
- REAL, INTENT(INOUT) :: CHR_URB
- REAL, INTENT(INOUT) :: CMC_URB
- REAL, INTENT(INOUT) :: CHC_URB
- !-------------------------------------------------------------------------------
- ! O: output variables from Urban to LSM
- !-------------------------------------------------------------------------------
- REAL, INTENT(OUT) :: TS ! surface potential temperature [K]
- REAL, INTENT(OUT) :: QS ! surface humidity [K]
- REAL, INTENT(OUT) :: SH ! sensible heat flux [W/m/m]
- REAL, INTENT(OUT) :: LH ! latent heat flux [W/m/m]
- REAL, INTENT(OUT) :: LH_KINEMATIC ! latent heat, kinetic [kg/m/m/s]
- REAL, INTENT(OUT) :: SW ! upward short wave radiation flux [W/m/m]
- REAL, INTENT(OUT) :: ALB ! time-varying albedo [fraction]
- REAL, INTENT(OUT) :: LW ! upward long wave radiation flux [W/m/m]
- REAL, INTENT(OUT) :: G ! heat flux into the ground [W/m/m]
- REAL, INTENT(OUT) :: RN ! net radition [W/m/m]
- REAL, INTENT(OUT) :: PSIM ! similality stability shear function for momentum
- REAL, INTENT(OUT) :: PSIH ! similality stability shear function for heat
- REAL, INTENT(OUT) :: GZ1OZ0
- REAL, INTENT(OUT) :: U10 ! u at 10m [m/s]
- REAL, INTENT(OUT) :: V10 ! u at 10m [m/s]
- REAL, INTENT(OUT) :: TH2 ! potential temperature at 2 m [K]
- REAL, INTENT(OUT) :: Q2 ! humidity at 2 m [-]
- !m REAL, INTENT(OUT) :: CHS,CHS2 ! CH*U at za and 2 m [m/s]
- REAL, INTENT(OUT) :: UST ! friction velocity [m/s]
- !-------------------------------------------------------------------------------
- ! H: Historical (state) variables of Urban : LSM <--> Urban
- !-------------------------------------------------------------------------------
- ! TR: roof temperature [K]; TRP: at previous time step [K]
- ! TB: building wall temperature [K]; TBP: at previous time step [K]
- ! TG: road temperature [K]; TGP: at previous time step [K]
- ! TC: urban-canopy air temperature [K]; TCP: at previous time step [K]
- ! (absolute temperature)
- ! QC: urban-canopy air mixing ratio [kg/kg]; QCP: at previous time step [kg/kg]
- !
- ! XXXR: Monin-Obkhov length for roof [dimensionless]
- ! XXXB: Monin-Obkhov length for building wall [dimensionless]
- ! XXXG: Monin-Obkhov length for road [dimensionless]
- ! XXXC: Monin-Obkhov length for urban-canopy [dimensionless]
- !
- ! TRL, TBL, TGL: layer temperature [K] (absolute temperature)
- REAL, INTENT(INOUT):: TR, TB, TG, TC, QC, UC
- REAL, INTENT(INOUT):: XXXR, XXXB, XXXG, XXXC
- REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: TRL
- REAL, DIMENSION(1:num_wall_layers), INTENT(INOUT) :: TBL
- REAL, DIMENSION(1:num_road_layers), INTENT(INOUT) :: TGL
- !-------------------------------------------------------------------------------
- ! L: Local variables from read_param
- !-------------------------------------------------------------------------------
- REAL :: ZR, Z0C, Z0HC, ZDC, SVF, R, RW, HGT, AH
- REAL :: SIGMA_ZED
- REAL :: CAPR, CAPB, CAPG, AKSR, AKSB, AKSG, ALBR, ALBB, ALBG
- REAL :: EPSR, EPSB, EPSG, Z0R, Z0B, Z0G, Z0HB, Z0HG
- REAL :: TRLEND,TBLEND,TGLEND
- REAL :: T1VR, T1VC,TH2V
- REAL :: RLMO_URB
- REAL :: AKANDA_URBAN
-
- REAL :: TH2X !m
-
- INTEGER :: BOUNDR, BOUNDB, BOUNDG
- INTEGER :: CH_SCHEME, TS_SCHEME
- LOGICAL :: SHADOW ! [true=consider svf and shadow effects, false=consider svf effect only]
- !for BEP
- INTEGER :: NUMDIR
- REAL, DIMENSION ( MAXDIRS ) :: STREET_DIRECTION
- REAL, DIMENSION ( MAXDIRS ) :: STREET_WIDTH
- REAL, DIMENSION ( MAXDIRS ) :: BUILDING_WIDTH
- INTEGER :: NUMHGT
- REAL, DIMENSION ( MAXHGTS ) :: HEIGHT_BIN
- REAL, DIMENSION ( MAXHGTS ) :: HPERCENT_BIN
- !end BEP
- !-------------------------------------------------------------------------------
- ! L: Local variables
- !-------------------------------------------------------------------------------
- REAL :: BETR, BETB, BETG
- REAL :: SX, SD, SQ, RX
- REAL :: UR, ZC, XLB, BB
- REAL :: Z, RIBB, RIBG, RIBC, BHR, BHB, BHG, BHC
- REAL :: TSC, LNET, SNET, FLXUV, THG, FLXTH, FLXHUM, FLXG
- REAL :: W, VFGS, VFGW, VFWG, VFWS, VFWW
- REAL :: HOUI1, HOUI2, HOUI3, HOUI4, HOUI5, HOUI6, HOUI7, HOUI8
- REAL :: SLX, SLX1, SLX2, SLX3, SLX4, SLX5, SLX6, SLX7, SLX8
- REAL :: FLXTHR, FLXTHB, FLXTHG, FLXHUMR, FLXHUMB, FLXHUMG
- REAL :: SR, SB, SG, RR, RB, RG
- REAL :: SR1, SR2, SB1, SB2, SG1, SG2, RR1, RR2, RB1, RB2, RG1, RG2
- REAL :: HR, HB, HG, ELER, ELEB, ELEG, G0R, G0B, G0G
- REAL :: ALPHAC, ALPHAR, ALPHAB, ALPHAG
- REAL :: CHC, CHR, CHB, CHG, CDC, CDR, CDB, CDG
- REAL :: C1R, C1B, C1G, TE, TC1, TC2, QC1, QC2, QS0R, QS0B, QS0G,RHO,ES
- REAL :: DESDT
- REAL :: F
- REAL :: DQS0RDTR
- REAL :: DRRDTR, DHRDTR, DELERDTR, DG0RDTR
- REAL :: DTR, DFDT
- REAL :: FX, FY, GF, GX, GY
- REAL :: DTCDTB, DTCDTG
- REAL :: DQCDTB, DQCDTG
- REAL :: DRBDTB1, DRBDTG1, DRBDTB2, DRBDTG2
- REAL :: DRGDTB1, DRGDTG1, DRGDTB2, DRGDTG2
- REAL :: DRBDTB, DRBDTG, DRGDTB, DRGDTG
- REAL :: DHBDTB, DHBDTG, DHGDTB, DHGDTG
- REAL :: DELEBDTB, DELEBDTG, DELEGDTG, DELEGDTB
- REAL :: DG0BDTB, DG0BDTG, DG0GDTG, DG0GDTB
- REAL :: DQS0BDTB, DQS0GDTG
- REAL :: DTB, DTG, DTC
- REAL :: THEATAZ ! Solar Zenith Angle [rad]
- REAL :: THEATAS ! = PI/2. - THETAZ
- REAL :: FAI ! Latitude [rad]
- REAL :: CNT,SNT
- REAL :: PS ! Surface Pressure [hPa]
- REAL :: TAV ! Vertial Temperature [K]
- REAL :: XXX, X, Z0, Z0H, CD, CH
- REAL :: XXX2, PSIM2, PSIH2, XXX10, PSIM10, PSIH10
- REAL :: PSIX, PSIT, PSIX2, PSIT2, PSIX10, PSIT10
- REAL :: TRP, TBP, TGP, TCP, QCP, TST, QST
- INTEGER :: iteration, K
- INTEGER :: tloc
- !-------------------------------------------------------------------------------
- ! Set parameters
- !-------------------------------------------------------------------------------
- ! Miao, 2007/01/17, cal. ah
- if(ahoption==1) then
- tloc=mod(int(OMG/PI*180./15.+12.+0.5 ),24)
- if(tloc==0) tloc=24
- endif
- CALL read_param(UTYPE,ZR,SIGMA_ZED,Z0C,Z0HC,ZDC,SVF,R,RW,HGT, &
- AH,CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB, &
- ALBG,EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HB,Z0HG, &
- BETR,BETB,BETG,TRLEND,TBLEND,TGLEND, &
- !for BEP
- NUMDIR, STREET_DIRECTION, STREET_WIDTH, &
- BUILDING_WIDTH, NUMHGT, HEIGHT_BIN, &
- HPERCENT_BIN, &
- !end BEP
- BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME, &
- AKANDA_URBAN)
- ! Miao, 2007/01/17, cal. ah
- if(ahoption==1) AH=AH*ahdiuprf(tloc)
- IF( ZDC+Z0C+2. >= ZA) THEN
- CALL wrf_error_fatal ("ZDC + Z0C + 2m is larger than the 1st WRF level "// &
- "Stop in subroutine urban - change ZDC and Z0C" )
- END IF
- IF(.NOT.LSOLAR) THEN
- SSGD = SRATIO*SSG
- SSGQ = SSG - SSGD
- ENDIF
- SSGD = SRATIO*SSG ! No radiation scheme has SSGD and SSGQ.
- SSGQ = SSG - SSGD
- W=2.*1.*HGT
- VFGS=SVF
- VFGW=1.-SVF
- VFWG=(1.-SVF)*(1.-R)/W
- VFWS=VFWG
- VFWW=1.-2.*VFWG
- !-------------------------------------------------------------------------------
- ! Convert unit from MKS to cgs
- ! Renew surface and layer temperatures
- !-------------------------------------------------------------------------------
- SX=(SSGD+SSGQ)/697.7/60. ! downward short wave radition [ly/min]
- SD=SSGD/697.7/60. ! downward direct short wave radiation
- SQ=SSGQ/697.7/60. ! downward diffuse short wave radiation
- RX=LLG/697.7/60. ! downward long wave radiation
- RHO=RHOO*0.001 ! air density at first atmospheric level
- TRP=TR
- TBP=TB
- TGP=TG
- TCP=TC
- QCP=QC
- TAV=TA*(1.+0.61*QA)
- PS=RHOO*287.*TAV/100. ![hPa]
- !-------------------------------------------------------------------------------
- ! Canopy wind
- !-------------------------------------------------------------------------------
- IF ( ZR + 2. < ZA ) THEN
- UR=UA*LOG((ZR-ZDC)/Z0C)/LOG((ZA-ZDC)/Z0C)
- ZC=0.7*ZR
- XLB=0.4*(ZR-ZDC)
- ! BB formulation from Inoue (1963)
- BB = 0.4 * ZR / ( XLB * alog( ( ZR - ZDC ) / Z0C ) )
- UC=UR*EXP(-BB*(1.-ZC/ZR))
- ELSE
- ! PRINT *, 'Warning ZR + 2m is larger than the 1st WRF level'
- ZC=ZA/2.
- UC=UA/2.
- END IF
- !-------------------------------------------------------------------------------
- ! Net Short Wave Radiation at roof, wall, and road
- !-------------------------------------------------------------------------------
- SHADOW = .false.
- ! SHADOW = .true.
- IF (SSG > 0.0) THEN
- IF(.NOT.SHADOW) THEN ! no shadow effects model
- SR1=SX*(1.-ALBR)
- SG1=SX*VFGS*(1.-ALBG)
- SB1=SX*VFWS*(1.-ALBB)
- SG2=SB1*ALBB/(1.-ALBB)*VFGW*(1.-ALBG)
- SB2=SG1*ALBG/(1.-ALBG)*VFWG*(1.-ALBB)
- ELSE ! shadow effects model
- FAI=XLAT*PI/180.
- THEATAS=ABS(ASIN(COSZ))
- THEATAZ=ABS(ACOS(COSZ))
- SNT=COS(DECLIN)*SIN(OMG)/COS(THEATAS)
- CNT=(COSZ*SIN(FAI)-SIN(DECLIN))/COS(THEATAS)/COS(FAI)
- HOUI1=(SNT*COS(PI/8.) -CNT*SIN(PI/8.))
- HOUI2=(SNT*COS(2.*PI/8.) -CNT*SIN(2.*PI/8.))
- HOUI3=(SNT*COS(3.*PI/8.) -CNT*SIN(3.*PI/8.))
- HOUI4=(SNT*COS(4.*PI/8.) -CNT*SIN(4.*PI/8.))
- HOUI5=(SNT*COS(5.*PI/8.) -CNT*SIN(5.*PI/8.))
- HOUI6=(SNT*COS(6.*PI/8.) -CNT*SIN(6.*PI/8.))
- HOUI7=(SNT*COS(7.*PI/8.) -CNT*SIN(7.*PI/8.))
- HOUI8=(SNT*COS(8.*PI/8.) -CNT*SIN(8.*PI/8.))
- SLX1=HGT*ABS(TAN(THEATAZ))*ABS(HOUI1)
- SLX2=HGT*ABS(TAN(THEATAZ))*ABS(HOUI2)
- SLX3=HGT*ABS(TAN(THEATAZ))*ABS(HOUI3)
- SLX4=HGT*ABS(TAN(THEATAZ))*ABS(HOUI4)
- SLX5=HGT*ABS(TAN(THEATAZ))*ABS(HOUI5)
- SLX6=HGT*ABS(TAN(THEATAZ))*ABS(HOUI6)
- SLX7=HGT*ABS(TAN(THEATAZ))*ABS(HOUI7)
- SLX8=HGT*ABS(TAN(THEATAZ))*ABS(HOUI8)
- IF(SLX1 > RW) SLX1=RW
- IF(SLX2 > RW) SLX2=RW
- IF(SLX3 > RW) SLX3=RW
- IF(SLX4 > RW) SLX4=RW
- IF(SLX5 > RW) SLX5=RW
- IF(SLX6 > RW) SLX6=RW
- IF(SLX7 > RW) SLX7=RW
- IF(SLX8 > RW) SLX8=RW
- SLX=(SLX1+SLX2+SLX3+SLX4+SLX5+SLX6+SLX7+SLX8)/8.
- SR1=SD*(1.-ALBR)+SQ*(1.-ALBR)
- SG1=SD*(RW-SLX)/RW*(1.-ALBG)+SQ*VFGS*(1.-ALBG)
- SB1=SD*SLX/W*(1.-ALBB)+SQ*VFWS*(1.-ALBB)
- SG2=SB1*ALBB/(1.-ALBB)*VFGW*(1.-ALBG)
- SB2=SG1*ALBG/(1.-ALBG)*VFWG*(1.-ALBB)
- END IF
- SR=SR1
- SG=SG1+SG2
- SB=SB1+SB2
- SNET=R*SR+W*SB+RW*SG
- ELSE
- SR=0.
- SG=0.
- SB=0.
- SNET=0.
- END IF
- !-------------------------------------------------------------------------------
- ! Roof
- !-------------------------------------------------------------------------------
- !-------------------------------------------------------------------------------
- ! CHR, CDR, BETR
- !-------------------------------------------------------------------------------
- ! Z=ZA-ZDC
- ! BHR=LOG(Z0R/Z0HR)/0.4
- ! RIBR=(9.8*2./(TA+TRP))*(TA-TRP)*(Z+Z0R)/(UA*UA)
- ! CALL mos(XXXR,ALPHAR,CDR,BHR,RIBR,Z,Z0R,UA,TA,TRP,RHO)
- ! Alternative option for MOST using SFCDIF routine from Noah
- ! Virtual temperatures needed by SFCDIF
- T1VR = TRP* (1.0+ 0.61 * QA)
- TH2V = (TA + ( 0.0098 * ZA)) * (1.0+ 0.61 * QA)
-
- ! note that CHR_URB contains UA (=CHR_MOS*UA)
- RLMO_URB=0.0
- CALL SFCDIF_URB (ZA,Z0R,T1VR,TH2V,UA,AKANDA_URBAN,CMR_URB,CHR_URB,RLMO_URB,CDR)
- ALPHAR = RHO*CP*CHR_URB
- CHR=ALPHAR/RHO/CP/UA
- IF(RAIN > 1.) BETR=0.7
- IF (TS_SCHEME == 1) THEN
- !-------------------------------------------------------------------------------
- ! TR Solving Non-Linear Equation by Newton-Rapson
- ! TRL Solving Heat Equation by Tri Diagonal Matrix Algorithm
- !-------------------------------------------------------------------------------
- ! TSC=TRP-273.15
- ! ES=EXP(19.482-4303.4/(TSC+243.5)) ! WMO
- ! ES=6.11*10.**(TETENA*TSC/(TETENB+TSC)) ! Tetens
- ! DESDT=( 6.1078*(2500.-2.4*TSC)/ & ! Tetens
- ! (0.46151*(TSC+273.15)**2.) )*10.**(7.5*TSC/(237.3+TSC))
- ! ES=6.11*EXP((2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) ) ! Clausius-Clapeyron
- ! DESDT=(2.5*10.**6./461.51)*ES/(TRP**2.) ! Clausius-Clapeyron
- ! QS0R=0.622*ES/(PS-0.378*ES)
- ! DQS0RDTR = DESDT*0.622*PS/((PS-0.378*ES)**2.)
- ! DQS0RDTR = 17.269*(273.15-35.86)/((TRP-35.86)**2.)*QS0R
- ! TRP=350.
- DO ITERATION=1,20
- ES=6.11*EXP( (2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) )
- DESDT=(2.5*10.**6./461.51)*ES/(TRP**2.)
- QS0R=0.622*ES/(PS-0.378*ES)
- DQS0RDTR = DESDT*0.622*PS/((PS-0.378*ES)**2.)
- RR=EPSR*(RX-SIG*(TRP**4.)/60.)
- HR=RHO*CP*CHR*UA*(TRP-TA)*100.
- ELER=RHO*EL*CHR*UA*BETR*(QS0R-QA)*100.
- G0R=AKSR*(TRP-TRL(1))/(DZR(1)/2.)
- F = SR + RR - HR - ELER - G0R
- DRRDTR = (-4.*EPSR*SIG*TRP**3.)/60.
- DHRDTR = RHO*CP*CHR*UA*100.
- DELERDTR = RHO*EL*CHR*UA*BETR*DQS0RDTR*100.
- DG0RDTR = 2.*AKSR/DZR(1)
- DFDT = DRRDTR - DHRDTR - DELERDTR - DG0RDTR
- DTR = F/DFDT
- TR = TRP - DTR
- TRP = TR
- IF( ABS(F) < 0.000001 .AND. ABS(DTR) < 0.000001 ) EXIT
- END DO
- ! multi-layer heat equation model
- CALL multi_layer(num_roof_layers,BOUNDR,G0R,CAPR,AKSR,TRL,DZR,DELT,TRLEND)
- ELSE
- ES=6.11*EXP( (2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) )
- QS0R=0.622*ES/(PS-0.378*ES)
- RR=EPSR*(RX-SIG*(TRP**4.)/60.)
- HR=RHO*CP*CHR*UA*(TRP-TA)*100.
- ELER=RHO*EL*CHR*UA*BETR*(QS0R-QA)*100.
- G0R=SR+RR-HR-ELER
- CALL force_restore(CAPR,AKSR,DELT,SR,RR,HR,ELER,TRLEND,TRP,TR)
- TRP=TR
- END IF
- FLXTHR=HR/RHO/CP/100.
- FLXHUMR=ELER/RHO/EL/100.
- !-------------------------------------------------------------------------------
- ! Wall and Road
- !-------------------------------------------------------------------------------
- !-------------------------------------------------------------------------------
- ! CHC, CHB, CDB, BETB, CHG, CDG, BETG
- !-------------------------------------------------------------------------------
- ! Z=ZA-ZDC
- ! BHC=LOG(Z0C/Z0HC)/0.4
- ! RIBC=(9.8*2./(TA+TCP))*(TA-TCP)*(Z+Z0C)/(UA*UA)
- !
- ! CALL mos(XXXC,ALPHAC,CDC,BHC,RIBC,Z,Z0C,UA,TA,TCP,RHO)
- ! Virtual temperatures needed by SFCDIF routine from Noah
- T1VC = TCP* (1.0+ 0.61 * QA)
- RLMO_URB=0.0
- CALL SFCDIF_URB(ZA,Z0C,T1VC,TH2V,UA,AKANDA_URBAN,CMC_URB,CHC_URB,RLMO_URB,CDC)
- ALPHAC = RHO*CP*CHC_URB
- IF (CH_SCHEME == 1) THEN
- Z=ZDC
- BHB=LOG(Z0B/Z0HB)/0.4
- BHG=LOG(Z0G/Z0HG)/0.4
- RIBB=(9.8*2./(TCP+TBP))*(TCP-TBP)*(Z+Z0B)/(UC*UC)
- RIBG=(9.8*2./(TCP+TGP))*(TCP-TGP)*(Z+Z0G)/(UC*UC)
- CALL mos(XXXB,ALPHAB,CDB,BHB,RIBB,Z,Z0B,UC,TCP,TBP,RHO)
- CALL mos(XXXG,ALPHAG,CDG,BHG,RIBG,Z,Z0G,UC,TCP,TGP,RHO)
- ELSE
- ALPHAB=RHO*CP*(6.15+4.18*UC)/1200.
- IF(UC > 5.) ALPHAB=RHO*CP*(7.51*UC**0.78)/1200.
- ALPHAG=RHO*CP*(6.15+4.18*UC)/1200.
- IF(UC > 5.) ALPHAG=RHO*CP*(7.51*UC**0.78)/1200.
- END IF
- CHC=ALPHAC/RHO/CP/UA
- CHB=ALPHAB/RHO/CP/UC
- CHG=ALPHAG/RHO/CP/UC
- BETB=0.0
- IF(RAIN > 1.) BETG=0.7
- IF (TS_SCHEME == 1) THEN
- !-------------------------------------------------------------------------------
- ! TB, TG Solving Non-Linear Simultaneous Equation by Newton-Rapson
- ! TBL,TGL Solving Heat Equation by Tri Diagonal Matrix Algorithm
- !-------------------------------------------------------------------------------
- ! TBP=350.
- ! TGP=350.
- DO ITERATION=1,20
- ES=6.11*EXP( (2.5*10.**6./461.51)*(TBP-273.15)/(273.15*TBP) )
- DESDT=(2.5*10.**6./461.51)*ES/(TBP**2.)
- QS0B=0.622*ES/(PS-0.378*ES)
- DQS0BDTB=DESDT*0.622*PS/((PS-0.378*ES)**2.)
- ES=6.11*EXP( (2.5*10.**6./461.51)*(TGP-273.15)/(273.15*TGP) )
- DESDT=(2.5*10.**6./461.51)*ES/(TGP**2.)
- QS0G=0.622*ES/(PS-0.378*ES)
- DQS0GDTG=DESDT*0.22*PS/((PS-0.378*ES)**2.)
- RG1=EPSG*( RX*VFGS &
- +EPSB*VFGW*SIG*TBP**4./60. &
- -SIG*TGP**4./60. )
- RB1=EPSB*( RX*VFWS &
- +EPSG*VFWG*SIG*TGP**4./60. &
- +EPSB*VFWW*SIG*TBP**4./60. &
- -SIG*TBP**4./60. )
- RG2=EPSG*( (1.-EPSB)*(1.-SVF)*VFWS*RX &
- +(1.-EPSB)*(1.-SVF)*VFWG*EPSG*SIG*TGP**4./60. &
- +EPSB*(1.-EPSB)*(1.-SVF)*(1.-2.*VFWS)*SIG*TBP**4./60. )
- RB2=EPSB*( (1.-EPSG)*VFWG*VFGS*RX &
- +(1.-EPSG)*EPSB*VFGW*VFWG*SIG*(TBP**4.)/60. &
- +(1.-EPSB)*VFWS*(1.-2.*VFWS)*RX &
- +(1.-EPSB)*VFWG*(1.-2.*VFWS)*EPSG*SIG*EPSG*TGP**4./60. &
- +EPSB*(1.-EPSB)*(1.-2.*VFWS)*(1.-2.*VFWS)*SIG*TBP**4./60. )
- RG=RG1+RG2
- RB=RB1+RB2
- DRBDTB1=EPSB*(4.*EPSB*SIG*TB**3.*VFWW-4.*SIG*TB**3.)/60.
- DRBDTG1=EPSB*(4.*EPSG*SIG*TG**3.*VFWG)/60.
- DRBDTB2=EPSB*(4.*(1.-EPSG)*EPSB*SIG*TB**3.*VFGW*VFWG &
- +4.*EPSB*(1.-EPSB)*SIG*TB**3.*VFWW*VFWW)/60.
- DRBDTG2=EPSB*(4.*(1.-EPSB)*EPSG*SIG*TG**3.*VFWG*VFWW)/60.
- DRGDTB1=EPSG*(4.*EPSB*SIG*TB**3.*VFGW)/60.
- DRGDTG1=EPSG*(-4.*SIG*TG**3.)/60.
- DRGDTB2=EPSG*(4.*EPSB*(1.-EPSB)*SIG*TB**3.*VFWW*VFGW)/60.
- DRGDTG2=EPSG*(4.*(1.-EPSB)*EPSG*SIG*TG**3.*VFWG*VFGW)/60.
- DRBDTB=DRBDTB1+DRBDTB2
- DRBDTG=DRBDTG1+DRBDTG2
- DRGDTB=DRGDTB1+DRGDTB2
- DRGDTG=DRGDTG1+DRGDTG2
- HB=RHO*CP*CHB*UC*(TBP-TCP)*100.
- HG=RHO*CP*CHG*UC*(TGP-TCP)*100.
- DTCDTB=W*ALPHAB/(RW*ALPHAC+RW*ALPHAG+W*ALPHAB)
- DTCDTG=RW*ALPHAG/(RW*ALPHAC+RW*ALPHAG+W*ALPHAB)
- DHBDTB=RHO*CP*CHB*UC*(1.-DTCDTB)*100.
- DHBDTG=RHO*CP*CHB*UC*(0.-DTCDTG)*100.
- DHGDTG=RHO*CP*CHG*UC*(1.-DTCDTG)*100.
- DHGDTB=RHO*CP*CHG*UC*(0.-DTCDTB)*100.
- ELEB=RHO*EL*CHB*UC*BETB*(QS0B-QCP)*100.
- ELEG=RHO*EL*CHG*UC*BETG*(QS0G-QCP)*100.
- DQCDTB=W*ALPHAB*BETB*DQS0BDTB/(RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB)
- DQCDTG=RW*ALPHAG*BETG*DQS0GDTG/(RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB)
- DELEBDTB=RHO*EL*CHB*UC*BETB*(DQS0BDTB-DQCDTB)*100.
- DELEBDTG=RHO*EL*CHB*UC*BETB*(0.-DQCDTG)*100.
- DELEGDTG=RHO*EL*CHG*UC*BETG*(DQS0GDTG-DQCDTG)*100.
- DELEGDTB=RHO*EL*CHG*UC*BETG*(0.-DQCDTB)*100.
- G0B=AKSB*(TBP-TBL(1))/(DZB(1)/2.)
- G0G=AKSG*(TGP-TGL(1))/(DZG(1)/2.)
- DG0BDTB=2.*AKSB/DZB(1)
- DG0BDTG=0.
- DG0GDTG=2.*AKSG/DZG(1)
- DG0GDTB=0.
- F = SB + RB - HB - ELEB - G0B
- FX = DRBDTB - DHBDTB - DELEBDTB - DG0BDTB
- FY = DRBDTG - DHBDTG - DELEBDTG - DG0BDTG
- GF = SG + RG - HG - ELEG - G0G
- GX = DRGDTB - DHGDTB - DELEGDTB - DG0GDTB
- GY = DRGDTG - DHGDTG - DELEGDTG - DG0GDTG
- DTB = (GF*FY-F*GY)/(FX*GY-GX*FY)
- DTG = -(GF+GX*DTB)/GY
- TB = TBP + DTB
- TG = TGP + DTG
- TBP = TB
- TGP = TG
- TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB
- TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP
- TC=TC2/TC1
- QC1=RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB
- QC2=RW*ALPHAC*QA+RW*ALPHAG*BETG*QS0G+W*ALPHAB*BETB*QS0B
- QC=QC2/QC1
- DTC=TCP - TC
- TCP=TC
- QCP=QC
- IF( ABS(F) < 0.000001 .AND. ABS(DTB) < 0.000001 &
- .AND. ABS(GF) < 0.000001 .AND. ABS(DTG) < 0.000001 &
- .AND. ABS(DTC) < 0.000001) EXIT
- END DO
- CALL multi_layer(num_wall_layers,BOUNDB,G0B,CAPB,AKSB,TBL,DZB,DELT,TBLEND)
- CALL multi_layer(num_road_layers,BOUNDG,G0G,CAPG,AKSG,TGL,DZG,DELT,TGLEND)
- ELSE
- !-------------------------------------------------------------------------------
- ! TB, TG by Force-Restore Method
- !-------------------------------------------------------------------------------
- ES=6.11*EXP((2.5*10.**6./461.51)*(TBP-273.15)/(273.15*TBP) )
- QS0B=0.622*ES/(PS-0.378*ES)
- ES=6.11*EXP((2.5*10.**6./461.51)*(TGP-273.15)/(273.15*TGP) )
- QS0G=0.622*ES/(PS-0.378*ES)
- RG1=EPSG*( RX*VFGS &
- +EPSB*VFGW*SIG*TBP**4./60. &
- -SIG*TGP**4./60. )
- RB1=EPSB*( RX*VFWS &
- +EPSG*VFWG*SIG*TGP**4./60. &
- +EPSB*VFWW*SIG*TBP**4./60. &
- -SIG*TBP**4./60. )
- RG2=EPSG*( (1.-EPSB)*(1.-SVF)*VFWS*RX &
- +(1.-EPSB)*(1.-SVF)*VFWG*EPSG*SIG*TGP**4./60. &
- +EPSB*(1.-EPSB)*(1.-SVF)*(1.-2.*VFWS)*SIG*TBP**4./60. )
- RB2=EPSB*( (1.-EPSG)*VFWG*VFGS*RX &
- +(1.-EPSG)*EPSB*VFGW*VFWG*SIG*(TBP**4.)/60. &
- +(1.-EPSB)*VFWS*(1.-2.*VFWS)*RX &
- +(1.-EPSB)*VFWG*(1.-2.*VFWS)*EPSG*SIG*EPSG*TGP**4./60. &
- +EPSB*(1.-EPSB)*(1.-2.*VFWS)*(1.-2.*VFWS)*SIG*TBP**4./60. )
- RG=RG1+RG2
- RB=RB1+RB2
- HB=RHO*CP*CHB*UC*(TBP-TCP)*100.
- ELEB=RHO*EL*CHB*UC*BETB*(QS0B-QCP)*100.
- G0B=SB+RB-HB-ELEB
- HG=RHO*CP*CHG*UC*(TGP-TCP)*100.
- ELEG=RHO*EL*CHG*UC*BETG*(QS0G-QCP)*100.
- G0G=SG+RG-HG-ELEG
- CALL force_restore(CAPB,AKSB,DELT,SB,RB,HB,ELEB,TBLEND,TBP,TB)
- CALL force_restore(CAPG,AKSG,DELT,SG,RG,HG,ELEG,TGLEND,TGP,TG)
- TBP=TB
- TGP=TG
- TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB
- TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP
- TC=TC2/TC1
- QC1=RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB
- QC2=RW*ALPHAC*QA+RW*ALPHAG*BETG*QS0G+W*ALPHAB*BETB*QS0B
- QC=QC2/QC1
- TCP=TC
- QCP=QC
- END IF
- FLXTHB=HB/RHO/CP/100.
- FLXHUMB=ELEB/RHO/EL/100.
- FLXTHG=HG/RHO/CP/100.
- FLXHUMG=ELEG/RHO/EL/100.
- !-------------------------------------------------------------------------------
- ! Total Fluxes from Urban Canopy
- !-------------------------------------------------------------------------------
- FLXUV = ( R*CDR + RW*CDC )*UA*UA
- ! Miao, 2007/01/17, cal. ah
- if(ahoption==1) then
- FLXTH = ( R*FLXTHR + W*FLXTHB + RW*FLXTHG ) + AH/RHOO/CPP
- else
- FLXTH = ( R*FLXTHR + W*FLXTHB + RW*FLXTHG )
- endif
- FLXHUM = ( R*FLXHUMR + W*FLXHUMB + RW*FLXHUMG )
- FLXG = ( R*G0R + W*G0B + RW*G0G )
- LNET = R*RR + W*RB + RW*RG
- !----------------------------------------------------------------------------
- ! Convert Unit: FLUXES and u* T* q* --> WRF
- !----------------------------------------------------------------------------
- SH = FLXTH * RHOO * CPP ! Sensible heat flux [W/m/m]
- LH = FLXHUM * RHOO * ELL ! Latent heat flux [W/m/m]
- LH_KINEMATIC = FLXHUM * RHOO ! Latent heat, Kinematic [kg/m/m/s]
- LW = LLG - (LNET*697.7*60.) ! Upward longwave radiation [W/m/m]
- SW = SSG - (SNET*697.7*60.) ! Upward shortwave radiation [W/m/m]
- ALB = 0.
- IF( ABS(SSG) > 0.0001) ALB = SW/SSG ! Effective albedo [-]
- G = -FLXG*697.7*60. ! [W/m/m]
- RN = (SNET+LNET)*697.7*60. ! Net radiation [W/m/m]
- UST = SQRT(FLXUV) ! u* [m/s]
- TST = -FLXTH/UST ! T* [K]
- QST = -FLXHUM/UST ! q* [-]
- !------------------------------------------------------
- ! diagnostic GRID AVERAGED PSIM PSIH TS QS --> WRF
- !------------------------------------------------------
- Z0 = Z0C
- Z0H = Z0HC
- Z = ZA - ZDC
- XXX = 0.4*9.81*Z*TST/TA/UST/UST
- IF ( XXX >= 1. ) XXX = 1.
- IF ( XXX <= -5. ) XXX = -5.
- IF ( XXX > 0 ) THEN
- PSIM = -5. * XXX
- PSIH = -5. * XXX
- ELSE
- X = (1.-16.*XXX)**0.25
- PSIM = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + PI/2.
- PSIH = 2.*ALOG((1.+X*X)/2.)
- END IF
- GZ1OZ0 = ALOG(Z/Z0)
- CD = 0.4**2./(ALOG(Z/Z0)-PSIM)**2.
- !
- !m CH = 0.4**2./(ALOG(Z/Z0)-PSIM)/(ALOG(Z/Z0H)-PSIH)
- !m CHS = 0.4*UST/(ALOG(Z/Z0H)-PSIH)
- !m TS = TA + FLXTH/CH/UA ! surface potential temp (flux temp)
- !m QS = QA + FLXHUM/CH/UA ! surface humidity
- !
- TS = TA + FLXTH/CHS ! surface potential temp (flux temp)
- QS = QA + FLXHUM/CHS ! surface humidity
- !-------------------------------------------------------
- ! diagnostic GRID AVERAGED U10 V10 TH2 Q2 --> WRF
- !-------------------------------------------------------
- XXX2 = (2./Z)*XXX
- IF ( XXX2 >= 1. ) XXX2 = 1.
- IF ( XXX2 <= -5. ) XXX2 = -5.
- IF ( XXX2 > 0 ) THEN
- PSIM2 = -5. * XXX2
- PSIH2 = -5. * XXX2
- ELSE
- X = (1.-16.*XXX2)**0.25
- PSIM2 = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + 2.*ATAN(1.)
- PSIH2 = 2.*ALOG((1.+X*X)/2.)
- END IF
- !
- !m CHS2 = 0.4*UST/(ALOG(2./Z0H)-PSIH2)
- !
- XXX10 = (10./Z)*XXX
- IF ( XXX10 >= 1. ) XXX10 = 1.
- IF ( XXX10 <= -5. ) XXX10 = -5.
- IF ( XXX10 > 0 ) THEN
- PSIM10 = -5. * XXX10
- PSIH10 = -5. * XXX10
- ELSE
- X = (1.-16.*XXX10)**0.25
- PSIM10 = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + 2.*ATAN(1.)
- PSIH10 = 2.*ALOG((1.+X*X)/2.)
- END IF
- PSIX = ALOG(Z/Z0) - PSIM
- PSIT = ALOG(Z/Z0H) - PSIH
- PSIX2 = ALOG(2./Z0) - PSIM2
- PSIT2 = ALOG(2./Z0H) - PSIH2
- PSIX10 = ALOG(10./Z0) - PSIM10
- PSIT10 = ALOG(10./Z0H) - PSIH10
- U10 = U1 * (PSIX10/PSIX) ! u at 10 m [m/s]
- V10 = V1 * (PSIX10/PSIX) ! v at 10 m [m/s]
- ! TH2 = TS + (TA-TS)*(PSIT2/PSIT) ! potential temp at 2 m [K]
- ! TH2 = TS + (TA-TS)*(PSIT2/PSIT) ! Fei: this seems to be temp (not potential) at 2 m [K]
- !Fei: consistant with M-O theory
- TH2 = TS + (TA-TS) *(CHS/CHS2)
- Q2 = QS + (QA-QS)*(PSIT2/PSIT) ! humidity at 2 m [-]
- ! TS = (LW/SIG_SI/0.88)**0.25 ! Radiative temperature [K]
- END SUBROUTINE urban
- !===============================================================================
- !
- ! mos
- !
- !===============================================================================
- SUBROUTINE mos(XXX,ALPHA,CD,B1,RIB,Z,Z0,UA,TA,TSF,RHO)
- ! XXX: z/L (requires iteration by Newton-Rapson method)
- ! B1: Stanton number
- ! PSIM: = PSIX of LSM
- ! PSIH: = PSIT of LSM
- IMPLICIT NONE
- REAL, PARAMETER :: CP=0.24
- REAL, INTENT(IN) :: B1, Z, Z0, UA, TA, TSF, RHO
- REAL, INTENT(OUT) :: ALPHA, CD
- REAL, INTENT(INOUT) :: XXX, RIB
- REAL :: XXX0, X, X0, FAIH, DPSIM, DPSIH
- REAL :: F, DF, XXXP, US, TS, AL, XKB, DD, PSIM, PSIH
- INTEGER :: NEWT
- INTEGER, PARAMETER :: NEWT_END=10
- IF(RIB <= -15.) RIB=-15.
- IF(RIB < 0.) THEN
- DO NEWT=1,NEWT_END
- IF(XXX >= 0.) XXX=-1.E-3
- XXX0=XXX*Z0/(Z+Z0)
- X=(1.-16.*XXX)**0.25
- X0=(1.-16.*XXX0)**0.25
- PSIM=ALOG((Z+Z0)/Z0) &
- -ALOG((X+1.)**2.*(X**2.+1.)) &
- +2.*ATAN(X) &
- +ALOG((X+1.)**2.*(X0**2.+1.)) &
- -2.*ATAN(X0)
- FAIH=1./SQRT(1.-16.*XXX)
- PSIH=ALOG((Z+Z0)/Z0)+0.4*B1 &
- -2.*ALOG(SQRT(1.-16.*XXX)+1.) &
- +2.*ALOG(SQRT(1.-16.*XXX0)+1.)
- DPSIM=(1.-16.*XXX)**(-0.25)/XXX &
- -(1.-16.*XXX0)**(-0.25)/XXX
- DPSIH=1./SQRT(1.-16.*XXX)/XXX &
- -1./SQRT(1.-16.*XXX0)/XXX
- F=RIB*PSIM**2./PSIH-XXX
- DF=RIB*(2.*DPSIM*PSIM*PSIH-DPSIH*PSIM**2.) &
- /PSIH**2.-1.
- XXXP=XXX
- XXX=XXXP-F/DF
- IF(XXX <= -10.) XXX=-10.
- END DO
- ELSE IF(RIB >= 0.142857) THEN
- XXX=0.714
- PSIM=ALOG((Z+Z0)/Z0)+7.*XXX
- PSIH=PSIM+0.4*B1
- ELSE
- AL=ALOG((Z+Z0)/Z0)
- XKB=0.4*B1
- DD=-4.*RIB*7.*XKB*AL+(AL+XKB)**2.
- IF(DD <= 0.) DD=0.
- XXX=(AL+XKB-2.*RIB*7.*AL-SQRT(DD))/(2.*(RIB*7.**2-7.))
- PSIM=ALOG((Z+Z0)/Z0)+7.*MIN(XXX,0.714)
- PSIH=PSIM+0.4*B1
- END IF
- US=0.4*UA/PSIM ! u*
- IF(US <= 0.01) US=0.01
- TS=0.4*(TA-TSF)/PSIH ! T*
- CD=US*US/UA**2. ! CD
- ALPHA=RHO*CP*0.4*US/PSIH ! RHO*CP*CH*U
- RETURN
- END SUBROUTINE mos
- !===============================================================================
- !
- ! louis79
- !
- !===============================================================================
- SUBROUTINE louis79(ALPHA,CD,RIB,Z,Z0,UA,RHO)
- IMPLICIT NONE
- REAL, PARAMETER :: CP=0.24
- REAL, INTENT(IN) :: Z, Z0, UA, RHO
- REAL, INTENT(OUT) :: ALPHA, CD
- REAL, INTENT(INOUT) :: RIB
- REAL :: A2, XX, CH, CMB, CHB
- A2=(0.4/ALOG(Z/Z0))**2.
- IF(RIB <= -15.) RIB=-15.
- IF(RIB >= 0.0) THEN
- IF(RIB >= 0.142857) THEN
- XX=0.714
- ELSE
- XX=RIB*LOG(Z/Z0)/(1.-7.*RIB)
- END IF
- CH=0.16/0.74/(LOG(Z/Z0)+7.*MIN(XX,0.714))**2.
- CD=0.16/(LOG(Z/Z0)+7.*MIN(XX,0.714))**2.
- ELSE
- CMB=7.4*A2*9.4*SQRT(Z/Z0)
- CHB=5.3*A2*9.4*SQRT(Z/Z0)
- CH=A2/0.74*(1.-9.4*RIB/(1.+CHB*SQRT(-RIB)))
- CD=A2*(1.-9.4*RIB/(1.+CHB*SQRT(-RIB)))
- END IF
- ALPHA=RHO*CP*CH*UA
- RETURN
- END SUBROUTINE louis79
- !===============================================================================
- !
- ! louis82
- !
- !===============================================================================
- SUBROUTINE louis82(ALPHA,CD,RIB,Z,Z0,UA,RHO)
- IMPLICIT NONE
- REAL, PARAMETER :: CP=0.24
- REAL, INTENT(IN) :: Z, Z0, UA, RHO
- REAL, INTENT(OUT) :: ALPHA, CD
- REAL, INTENT(INOUT) :: RIB
- REAL :: A2, FM, FH, CH, CHH
- A2=(0.4/ALOG(Z/Z0))**2.
- IF(RIB <= -15.) RIB=-15.
- IF(RIB >= 0.0) THEN
- FM=1./((1.+(2.*5.*RIB)/SQRT(1.+5.*RIB)))
- FH=1./(1.+(3.*5.*RIB)*SQRT(1.+5.*RIB))
- CH=A2*FH
- CD=A2*FM
- ELSE
- CHH=5.*3.*5.*A2*SQRT(Z/Z0)
- FM=1.-(2.*5.*RIB)/(1.+3.*5.*5.*A2*SQRT(Z/Z0+1.)*(-RIB))
- FH=1.-(3.*5.*RIB)/(1.+CHH*SQRT(-RIB))
- CH=A2*FH
- CD=A2*FM
- END IF
- ALPHA=RHO*CP*CH*UA
- RETURN
- END SUBROUTINE louis82
- !===============================================================================
- !
- ! multi_layer
- !
- !===============================================================================
- SUBROUTINE multi_layer(KM,BOUND,G0,CAP,AKS,TSL,DZ,DELT,TSLEND)
- IMPLICIT NONE
- REAL, INTENT(IN) :: G0
- REAL, INTENT(IN) :: CAP
- REAL, INTENT(IN) :: AKS
- REAL, INTENT(IN) :: DELT ! Time step [ s ]
- REAL, INTENT(IN) :: TSLEND
- INTEGER, INTENT(IN) :: KM
- INTEGER, INTENT(IN) :: BOUND
- REAL, DIMENSION(KM), INTENT(IN) :: DZ
- REAL, DIMENSION(KM), INTENT(INOUT) :: TSL
- REAL, DIMENSION(KM) :: A, B, C, D, X, P, Q
- REAL :: DZEND
- INTEGER :: K
- DZEND=DZ(KM)
- A(1) = 0.0
- B(1) = CAP*DZ(1)/DELT &
- +2.*AKS/(DZ(1)+DZ(2))
- C(1) = -2.*AKS/(DZ(1)+DZ(2))
- D(1) = CAP*DZ(1)/DELT*TSL(1) + G0
- DO K=2,KM-1
- A(K) = -2.*AKS/(DZ(K-1)+DZ(K))
- B(K) = CAP*DZ(K)/DELT + 2.*AKS/(DZ(K-1)+DZ(K)) + 2.*AKS/(DZ(K)+DZ(K+1))
- C(K) = -2.*AKS/(DZ(K)+DZ(K+1))
- D(K) = CAP*DZ(K)/DELT*TSL(K)
- END DO
- IF(BOUND == 1) THEN ! Flux=0
- A(KM) = -2.*AKS/(DZ(KM-1)+DZ(KM))
- B(KM) = CAP*DZ(KM)/DELT + 2.*AKS/(DZ(KM-1)+DZ(KM))
- C(KM) = 0.0
- D(KM) = CAP*DZ(KM)/DELT*TSL(KM)
- ELSE ! T=constant
- A(KM) = -2.*AKS/(DZ(KM-1)+DZ(KM))
- B(KM) = CAP*DZ(KM)/DELT + 2.*AKS/(DZ(KM-1)+DZ(KM)) + 2.*AKS/(DZ(KM)+DZEND)
- C(KM) = 0.0
- D(KM) = CAP*DZ(KM)/DELT*TSL(KM) + 2.*AKS*TSLEND/(DZ(KM)+DZEND)
- END IF
- P(1) = -C(1)/B(1)
- Q(1) = D(1)/B(1)
- DO K=2,KM
- P(K) = -C(K)/(A(K)*P(K-1)+B(K))
- Q(K) = (-A(K)*Q(K-1)+D(K))/(A(K)*P(K-1)+B(K))
- END DO
- X(KM) = Q(KM)
- DO K=KM-1,1,-1
- X(K) = P(K)*X(K+1)+Q(K)
- END DO
- DO K=1,KM
- TSL(K) = X(K)
- END DO
- RETURN
- END SUBROUTINE multi_layer
- !===============================================================================
- !
- ! subroutine read_param
- !
- !===============================================================================
- SUBROUTINE read_param(UTYPE, & ! in
- ZR,SIGMA_ZED,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,AH, & ! out
- CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG, & ! out
- EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HB,Z0HG, & ! out
- BETR,BETB,BETG,TRLEND,TBLEND,TGLEND, & ! out
- !for BEP
- NUMDIR, STREET_DIRECTION, STREET_WIDTH, & ! out
- BUILDING_WIDTH, NUMHGT, HEIGHT_BIN, & ! out
- HPERCENT_BIN, & ! out
- !end BEP
- BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME, & ! out
- AKANDA_URBAN) ! out
- INTEGER, INTENT(IN) :: UTYPE
- REAL, INTENT(OUT) :: ZR,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,AH, &
- CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG, &
- SIGMA_ZED, &
- EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HB,Z0HG, &
- BETR,BETB,BETG,TRLEND,TBLEND,TGLEND
- REAL, INTENT(OUT) :: AKANDA_URBAN
- !for BEP
- INTEGER, INTENT(OUT) :: NUMDIR
- REAL, DIMENSION(MAXDIRS), INTENT(OUT) :: STREET_DIRECTION
- REAL, DIMENSION(MAXDIRS), INTENT(OUT) :: STREET_WIDTH
- REAL, DIMENSION(MAXDIRS), INTENT(OUT) :: BUILDING_WIDTH
- INTEGER, INTENT(OUT) :: NUMHGT
- REAL, DIMENSION(MAXHGTS), INTENT(OUT) :: HEIGHT_BIN
- REAL, DIMENSION(MAXHGTS), INTENT(OUT) :: HPERCENT_BIN
-
- !end BEP
- INTEGER, INTENT(OUT) :: BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME
- ZR = ZR_TBL(UTYPE)
- SIGMA_ZED = SIGMA_ZED_TBL(UTYPE)
- Z0C= Z0C_TBL(UTYPE)
- Z0HC= Z0HC_TBL(UTYPE)
- ZDC= ZDC_TBL(UTYPE)
- SVF= SVF_TBL(UTYPE)
- R= R_TBL(UTYPE)
- RW= RW_TBL(UTYPE)
- HGT= HGT_TBL(UTYPE)
- AH= AH_TBL(UTYPE)
- BETR= BETR_TBL(UTYPE)
- BETB= BETB_TBL(UTYPE)
- BETG= BETG_TBL(UTYPE)
- !m FRC_URB= FRC_URB_TBL(UTYPE)
- CAPR= CAPR_TBL(UTYPE)
- CAPB= CAPB_TBL(UTYPE)
- CAPG= CAPG_TBL(UTYPE)
- AKSR= AKSR_TBL(UTYPE)
- AKSB= AKSB_TBL(UTYPE)
- AKSG= AKSG_TBL(UTYPE)
- ALBR= ALBR_TBL(UTYPE)
- ALBB= ALBB_TBL(UTYPE)
- ALBG= ALBG_TBL(UTYPE)
- EPSR= EPSR_TBL(UTYPE)
- EPSB= EPSB_TBL(UTYPE)
- EPSG= EPSG_TBL(UTYPE)
- Z0R= Z0R_TBL(UTYPE)
- Z0B= Z0B_TBL(UTYPE)
- Z0G= Z0G_TBL(UTYPE)
- Z0HB= Z0HB_TBL(UTYPE)
- Z0HG= Z0HG_TBL(UTYPE)
- TRLEND= TRLEND_TBL(UTYPE)
- TBLEND= TBLEND_TBL(UTYPE)
- TGLEND= TGLEND_TBL(UTYPE)
- BOUNDR= BOUNDR_DATA
- BOUNDB= BOUNDB_DATA
- BOUNDG= BOUNDG_DATA
- CH_SCHEME = CH_SCHEME_DATA
- TS_SCHEME = TS_SCHEME_DATA
- AKANDA_URBAN = AKANDA_URBAN_TBL(UTYPE)
- !for BEP
- STREET_DIRECTION = -1.E36
- STREET_WIDTH = -1.E36
- BUILDING_WIDTH = -1.E36
- HEIGHT_BIN = -1.E36
- HPERCENT_BIN = -1.E36
- NUMDIR = NUMDIR_TBL ( UTYPE )
- STREET_DIRECTION(1:NUMDIR) = STREET_DIRECTION_TBL( 1:NUMDIR, UTYPE )
- STREET_WIDTH (1:NUMDIR) = STREET_WIDTH_TBL ( 1:NUMDIR, UTYPE )
- BUILDING_WIDTH (1:NUMDIR) = BUILDING_WIDTH_TBL ( 1:NUMDIR, UTYPE )
- NUMHGT = NUMHGT_TBL ( UTYPE )
- HEIGHT_BIN (1:NUMHGT) = HEIGHT_BIN_TBL ( 1:NUMHGT , UTYPE )
- HPERCENT_BIN (1:NUMHGT) = HPERCENT_BIN_TBL ( 1:NUMHGT , UTYPE )
- !end BEP
- END SUBROUTINE read_param
- !===============================================================================
- !
- ! subroutine urban_param_init: Read parameters from URBPARM.TBL
- !
- !===============================================================================
- SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers, &
- sf_urban_physics)
- ! num_roof_layers,num_wall_layers,num_road_layers)
- IMPLICIT NONE
- INTEGER, INTENT(IN) :: num_soil_layers
- ! REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: DZR
- ! REAL, DIMENSION(1:num_wall_layers), INTENT(INOUT) :: DZB
- ! REAL, DIMENSION(1:num_road_layers), INTENT(INOUT) :: DZG
- REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZR
- REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZB
- REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZG
- INTEGER, INTENT(IN) :: SF_URBAN_PHYSICS
- INTEGER :: LC, K
- INTEGER :: IOSTATUS, ALLOCATE_STATUS
- INTEGER :: num_roof_layers
- INTEGER :: num_wall_layers
- INTEGER :: num_road_layers
- INTEGER :: dummy
- REAL :: DHGT, HGT, VFWS, VFGS
- REAL, allocatable, dimension(:) :: ROOF_WIDTH
- REAL, allocatable, dimension(:) :: ROAD_WIDTH
- character(len=512) :: string
- character(len=128) :: name
- integer :: indx
- real, parameter :: VonK = 0.4
- real :: lambda_p
- real :: lambda_f
- real :: Cd
- real :: alpha_macd
- real :: beta_macd
- real :: lambda_fr
- !for BEP
- real :: dummy_hgt
- real :: dummy_pct
- real :: pctsum
- !end BEP
- num_roof_layers = num_soil_layers
- num_wall_layers = num_soil_layers
- num_road_layers = num_soil_layers
- ICATE=0
- OPEN (UNIT=11, &
- FILE='URBPARM.TBL', &
- ACCESS='SEQUENTIAL', &
- STATUS='OLD', &
- ACTION='READ', &
- POSITION='REWIND', &
- IOSTAT=IOSTATUS)
- IF (IOSTATUS > 0) THEN
- CALL wrf_error_fatal('ERROR OPEN URBPARM.TBL')
- ENDIF
- READLOOP : do
- read(11,'(A512)', iostat=iostatus) string
- if (iostatus /= 0) exit READLOOP
- if (string(1:1) == "#") cycle READLOOP
- if (trim(string) == "") cycle READLOOP
- indx = index(string,":")
- if (indx <= 0) cycle READLOOP
- name = trim(adjustl(string(1:indx-1)))
-
- ! Here are the variables we expect to be defined in the URBPARM.TBL:
- if (name == "Number of urban categories") then
- read(string(indx+1:),*) icate
- IF (.not. ALLOCATED(ZR_TBL)) then
- ALLOCATE( ZR_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ZR_TBL in urban_param_init')
- ALLOCATE( SIGMA_ZED_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0)CALL wrf_error_fatal('Error allocating SIGMA_ZED_TBL in urban_param_init')
- ALLOCATE( Z0C_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating Z0C_TBL in urban_param_init')
- ALLOCATE( Z0HC_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating Z0HC_TBL in urban_param_init')
- ALLOCATE( ZDC_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ZDC_TBL in urban_param_init')
- ALLOCATE( SVF_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating SVF_TBL in urban_param_init')
- ALLOCATE( R_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating R_TBL in urban_param_init')
- ALLOCATE( RW_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating RW_TBL in urban_param_init')
- ALLOCATE( HGT_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating HGT_TBL in urban_param_init')
- ALLOCATE( AH_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating AH_TBL in urban_param_init')
- ALLOCATE( BETR_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating BETR_TBL in urban_param_init')
- ALLOCATE( BETB_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating BETB_TBL in urban_param_init')
- ALLOCATE( BETG_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating BETG_TBL in urban_param_init')
- ALLOCATE( CAPR_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating CAPR_TBL in urban_param_init')
- ALLOCATE( CAPB_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating CAPB_TBL in urban_param_init')
- ALLOCATE( CAPG_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating CAPG_TBL in urban_param_init')
- ALLOCATE( AKSR_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating AKSR_TBL in urban_param_init')
- ALLOCATE( AKSB_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating AKSB_TBL in urban_param_init')
- ALLOCATE( AKSG_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating AKSG_TBL in urban_param_init')
- ALLOCATE( ALBR_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ALBR_TBL in urban_param_init')
- ALLOCATE( ALBB_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ALBB_TBL in urban_param_init')
- ALLOCATE( ALBG_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ALBG_TBL in urban_param_init')
- ALLOCATE( EPSR_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating EPSR_TBL in urban_param_init')
- ALLOCATE( EPSB_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating EPSB_TBL in urban_param_init')
- ALLOCATE( EPSG_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating EPSG_TBL in urban_param_init')
- ALLOCATE( Z0R_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating Z0R_TBL in urban_param_init')
- ALLOCATE( Z0B_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating Z0B_TBL in urban_param_init')
- ALLOCATE( Z0G_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating Z0G_TBL in urban_param_init')
- ALLOCATE( AKANDA_URBAN_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating AKANDA_URBAN_TBL in urban_param_init')
- ALLOCATE( Z0HB_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating Z0HB_TBL in urban_param_init')
- ALLOCATE( Z0HG_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating Z0HG_TBL in urban_param_init')
- ALLOCATE( TRLEND_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TRLEND_TBL in urban_param_init')
- ALLOCATE( TBLEND_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TBLEND_TBL in urban_param_init')
- ALLOCATE( TGLEND_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TGLEND_TBL in urban_param_init')
- ALLOCATE( FRC_URB_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating FRC_URB_TBL in urban_param_init')
- ! ALLOCATE( ROOF_WIDTH(ICATE), stat=allocate_status )
- ! if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ROOF_WIDTH in urban_param_init')
- ! ALLOCATE( ROAD_WIDTH(ICATE), stat=allocate_status )
- ! if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ROAD_WIDTH in urban_param_init')
- !for BEP
- ALLOCATE( NUMDIR_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating NUMDIR_TBL in urban_param_init')
- ALLOCATE( STREET_DIRECTION_TBL(MAXDIRS , ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating STREET_DIRECTION_TBL in urban_param_init')
- ALLOCATE( STREET_WIDTH_TBL(MAXDIRS , ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating STREET_WIDTH_TBL in urban_param_init')
- ALLOCATE( BUILDING_WIDTH_TBL(MAXDIRS , ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating BUILDING_WIDTH_TBL in urban_param_init')
- ALLOCATE( NUMHGT_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating NUMHGT_TBL in urban_param_init')
- ALLOCATE( HEIGHT_BIN_TBL(MAXHGTS , ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating HEIGHT_BIN_TBL in urban_param_init')
- ALLOCATE( HPERCENT_BIN_TBL(MAXHGTS , ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating HPERCENT_BIN_TBL in urban_param_init')
- ALLOCATE( COP_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating COP_TBL in urban_param_init')
- ALLOCATE( PWIN_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating PWIN_TBL in urban_param_init')
- ALLOCATE( BETA_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating BETA_TBL in urban_param_init')
- ALLOCATE( SW_COND_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating SW_COND_TBL in urban_param_init')
- ALLOCATE( TIME_ON_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TIME_ON_TBL in urban_param_init')
- ALLOCATE( TIME_OFF_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TIME_OFF_TBL in urban_param_init')
- ALLOCATE( TARGTEMP_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TARGTEMP_TBL in urban_param_init')
- ALLOCATE( GAPTEMP_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating GAPTEMP_TBL in urban_param_init')
- ALLOCATE( TARGHUM_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TARGHUM_TBL in urban_param_init')
- ALLOCATE( GAPHUM_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating GAPHUM_TBL in urban_param_init')
- ALLOCATE( PERFLO_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating PERFLO_TBL in urban_param_init')
- ALLOCATE( HSESF_TBL(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating HSESF_TBL in urban_param_init')
- endif
- numdir_tbl = 0
- street_direction_tbl = -1.E36
- street_width_tbl = 0
- building_width_tbl = 0
- numhgt_tbl = 0
- height_bin_tbl = -1.E36
- hpercent_bin_tbl = -1.E36
- !end BEP
- else if (name == "ZR") then
- read(string(indx+1:),*) zr_tbl(1:icate)
- else if (name == "SIGMA_ZED") then
- read(string(indx+1:),*) sigma_zed_tbl(1:icate)
- else if (name == "ROOF_WIDTH") then
- ALLOCATE( ROOF_WIDTH(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ROOF_WIDTH in urban_param_init')
- read(string(indx+1:),*) roof_width(1:icate)
- else if (name == "ROAD_WIDTH") then
- ALLOCATE( ROAD_WIDTH(ICATE), stat=allocate_status )
- if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ROAD_WIDTH in urban_param_init')
- read(string(indx+1:),*) road_width(1:icate)
- else if (name == "AH") then
- read(string(indx+1:),*) ah_tbl(1:icate)
- else if (name == "FRC_URB") then
- read(string(indx+1:),*) frc_urb_tbl(1:icate)
- else if (name == "CAPR") then
- read(string(indx+1:),*) capr_tbl(1:icate)
- ! Convert CAPR_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
- capr_tbl = capr_tbl * ( 1.0 / 4.1868 ) * 1.E-6
- else if (name == "CAPB") then
- read(string(indx+1:),*) capb_tbl(1:icate)
- ! Convert CABR_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
- capb_tbl = capb_tbl * ( 1.0 / 4.1868 ) * 1.E-6
- else if (name == "CAPG") then
- read(string(indx+1:),*) capg_tbl(1:icate)
- ! Convert CABG_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
- capg_tbl = capg_tbl * ( 1.0 / 4.1868 ) * 1.E-6
- else if (name == "AKSR") then
- read(string(indx+1:),*) aksr_tbl(1:icate)
- ! Convert AKSR_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
- AKSR_TBL = AKSR_TBL * ( 1.0 / 4.1868 ) * 1.E-2
- else if (name == "AKSB") then
- read(string(indx+1:),*) aksb_tbl(1:icate)
- ! Convert AKSB_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
- AKSB_TBL = AKSB_TBL * ( 1.0 / 4.1868 ) * 1.E-2
- else if (name == "AKSG") then
- read(string(indx+1:),*) aksg_tbl(1:icate)
- ! Convert AKSG_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
- AKSG_TBL = AKSG_TBL * ( 1.0 / 4.1868 ) * 1.E-2
- else if (name == "ALBR") then
- read(string(indx+1:),*) albr_tbl(1:icate)
- else if (name == "ALBB") then
- read(string(indx+1:),*) albb_tbl(1:icate)
- else if (name == "ALBG") then
- read(string(indx+1:),*) albg_tbl(1:icate)
- else if (name == "EPSR") then
- read(string(indx+1:),*) epsr_tbl(1:icate)
- else if (name == "EPSB") then
- read(string(indx+1:),*) epsb_tbl(1:icate)
- else if (name == "EPSG") then
- read(string(indx+1:),*) epsg_tbl(1:icate)
- else if (name == "AKANDA_URBAN") then
- read(string(indx+1:),*) akanda_urban_tbl(1:icate)
- else if (name == "Z0B") then
- read(string(indx+1:),*) z0b_tbl(1:icate)
- else if (name == "Z0G") then
- read(string(indx+1:),*) z0g_tbl(1:icate)
- else if (name == "DDZR") then
- read(string(indx+1:),*) dzr(1:num_roof_layers)
- ! Convert thicknesses from m to cm
- dzr = dzr * 100.0
- else if (name == "DDZB") then
- read(string(indx+1:),*) dzb(1:num_wall_layers)
- ! Convert thicknesses from m to cm
- dzb = dzb * 100.0
- else if (name == "DDZG") then
- read(string(indx+1:),*) dzg(1:num_road_layers)
- ! Convert thicknesses from m to cm
- dzg = dzg * 100.0
- else if (name == "BOUNDR") then
- read(string(indx+1:),*) boundr_data
- else if (name == "BOUNDB") then
- read(string(indx+1:),*) boundb_data
- else if (name == "BOUNDG") then
- read(string(indx+1:),*) boundg_data
- else if (name == "TRLEND") then
- read(string(indx+1:),*) trlend_tbl(1:icate)
- else if (name == "TBLEND") then
- read(string(indx+1:),*) tblend_tbl(1:icate)
- else if (name == "TGLEND") then
- read(string(indx+1:),*) tglend_tbl(1:icate)
- else if (name == "CH_SCHEME") then
- read(string(indx+1:),*) ch_scheme_data
- else if (name == "TS_SCHEME") then
- read(string(indx+1:),*) ts_scheme_data
- else if (name == "AHOPTION") then
- read(string(indx+1:),*) ahoption
- else if (name == "AHDIUPRF") then
- read(string(indx+1:),*) ahdiuprf(1:24)
- !for BEP
- else if (name == "STREET PARAMETERS") then
- STREETLOOP : do
- read(11,'(A512)', iostat=iostatus) string
- if (string(1:1) == "#") cycle STREETLOOP
- if (trim(string) == "") cycle STREETLOOP
- if (string == "END STREET PARAMETERS") exit STREETLOOP
- read(string, *) k ! , dirst, ws, bs
- numdir_tbl(k) = numdir_tbl(k) + 1
- read(string, *) k, street_direction_tbl(numdir_tbl(k),k), &
- street_width_tbl(numdir_tbl(k),k), &
- building_width_tbl(numdir_tbl(k),k)
- enddo STREETLOOP
- else if (name == "BUILDING HEIGHTS") then
- read(string(indx+1:),*) k
- HEIGHTLOOP : do
- read(11,'(A512)', iostat=iostatus) string
- if (string(1:1) == "#") cycle HEIGHTLOOP
- if (trim(string) == "") cycle HEIGHTLOOP
- if (string == "END BUILDING HEIGHTS") exit HEIGHTLOOP
- read(string,*) dummy_hgt, dummy_pct
- numhgt_tbl(k) = numhgt_tbl(k) + 1
- height_bin_tbl(numhgt_tbl(k), k) = dummy_hgt
- hpercent_bin_tbl(numhgt_tbl(k),k) = dummy_pct
-
- enddo HEIGHTLOOP
- pctsum = sum ( hpercent_bin_tbl(:,k) , mask=(hpercent_bin_tbl(:,k)>-1.E25 ) )
- if ( pctsum /= 100.) then
- write (*,'(//,"Building height percentages for category ", I2, " must sum to 100.0")') k
- write (*,'("Currently, they sum to ", F6.2,/)') pctsum
- CALL wrf_error_fatal('pctsum is not equal to 100.')
- endif
- else if ( name == "Z0R") then
- read(string(indx+1:),*) Z0R_tbl(1:icate)
- else if ( name == "COP") then
- read(string(indx+1:),*) cop_tbl(1:icate)
- else if ( name == "PWIN") then
- read(string(indx+1:),*) pwin_tbl(1:icate)
- else if ( name == "BETA") then
- read(string(indx+1:),*) beta_tbl(1:icate)
- else if ( name == "SW_COND") then
- read(string(indx+1:),*) sw_cond_tbl(1:icate)
- else if ( name == "TIME_ON") then
- read(string(indx+1:),*) time_on_tbl(1:icate)
- else if ( name == "TIME_OFF") then
- read(string(indx+1:),*) time_off_tbl(1:icate)
- else if ( name == "TARGTEMP") then
- read(string(indx+1:),*) targtemp_tbl(1:icate)
- else if ( name == "GAPTEMP") then
- read(string(indx+1:),*) gaptemp_tbl(1:icate)
- else if ( name == "TARGHUM") then
- read(string(indx+1:),*) targhum_tbl(1:icate)
- else if ( name == "GAPHUM") then
- read(string(indx+1:),*) gaphum_tbl(1:icate)
- else if ( name == "PERFLO") then
- read(string(indx+1:),*) perflo_tbl(1:icate)
- else if (name == "HSEQUIP") then
- read(string(indx+1:),*) hsequip_tbl(1:24)
- else if (name == "HSEQUIP_SCALE_FACTOR") then
- read(string(indx+1:),*) hsesf_tbl(1:icate)
- !end BEP
- else
- CALL wrf_error_fatal('URBPARM.TBL: Unrecognized NAME = "'//trim(name)//'" in Subr URBAN_PARAM_INIT')
- endif
- enddo READLOOP
- CLOSE(11)
- ! Assign a few table values that do not need to come from URBPARM.TBL
- Z0HB_TBL = 0.1 * Z0B_TBL
- Z0HG_TBL = 0.1 * Z0G_TBL
- DO LC = 1, ICATE
- ! HGT: Normalized height
- HGT_TBL(LC) = ZR_TBL(LC) / ( ROAD_WIDTH(LC) + ROOF_WIDTH(LC) )
- ! R: Normalized Roof Width (a.k.a. "building coverage ratio")
- R_TBL(LC) = ROOF_WIDTH(LC) / ( ROAD_WIDTH(LC) + ROOF_WIDTH(LC) )
- RW_TBL(LC) = 1.0 - R_TBL(LC)
- BETR_TBL(LC) = 0.0
- BETB_TBL(LC) = 0.0
- BETG_TBL(LC) = 0.0
- ! The following urban canyon geometry parameters are following Macdonald's (1998) formulations
-
- ! Lambda_P :: Plan areal fraction, which corresponds to R for a 2-d canyon.
- ! Lambda_F :: Frontal area index, which corresponds to HGT for a 2-d canyon
- ! Cd :: Drag coefficient ( 1.2 from Grimmond and Oke, 1998 )
- ! Alpha_macd :: Emperical coefficient ( 4.43 from Macdonald et al., 1998 )
- ! Beta_macd :: Correction factor for the drag coefficient ( 1.0 from Macdonald et al., 1998 )
- Lambda_P = R_TBL(LC)
- Lambda_F = HGT_TBL(LC)
- Cd = 1.2
- alpha_macd = 4.43
- beta_macd = 1.0
- ZDC_TBL(LC) = ZR_TBL(LC) * ( 1.0 + ( alpha_macd ** ( -Lambda_P ) ) * ( Lambda_P - 1.0 ) )
- Z0C_TBL(LC) = ZR_TBL(LC) * ( 1.0 - ZDC_TBL(LC)/ZR_TBL(LC) ) * &
- exp (-(0.5 * beta_macd * Cd / (VonK**2) * ( 1.0-ZDC_TBL(LC)/ZR_TBL(LC) ) * Lambda_F )**(-0.5))
- IF (SF_URBAN_PHYSICS == 1) THEN
- ! Include roof height variability in Macdonald
- ! to parameterize Z0R as a function of ZR_SD (Standard Deviation)
- Lambda_FR = SIGMA_ZED_TBL(LC) / ( ROAD_WIDTH(LC) + ROOF_WIDTH(LC) )
- Z0R_TBL(LC) = ZR_TBL(LC) * ( 1.0 - ZDC_TBL(LC)/ZR_TBL(LC) ) &
- * exp ( -(0.5 * beta_macd * Cd / (VonK**2) &
- * ( 1.0-ZDC_TBL(LC)/ZR_TBL(LC) ) * Lambda_FR )**(-0.5))
- ENDIF
- !
- ! Z0HC still one-tenth of Z0C, as before ?
- !
- Z0HC_TBL(LC) = 0.1 * Z0C_TBL(LC)
- !
- ! Calculate Sky View Factor:
- !
- DHGT=HGT_TBL(LC)/100.
- HGT=0.
- VFWS=0.
- HGT=HGT_TBL(LC)-DHGT/2.
- do k=1,99
- HGT=HGT-DHGT
- VFWS=VFWS+0.25*(1.-HGT/SQRT(HGT**2.+RW_TBL(LC)**2.))
- end do
- VFWS=VFWS/99.
- VFWS=VFWS*2.
- VFGS=1.-2.*VFWS*HGT_TBL(LC)/RW_TBL(LC)
- SVF_TBL(LC)=VFGS
- END DO
- deallocate(roof_width)
- deallocate(road_width)
- END SUBROUTINE urban_param_init
- !===========================================================================
- !
- ! subroutine urban_var_init: initialization of urban state variables
- !
- !===========================================================================
- SUBROUTINE urban_var_init(ISURBAN, TSURFACE0_URB,TLAYER0_URB,TDEEP0_URB,IVGTYP, & ! in
- ims,ime,jms,jme,kms,kme,num_soil_layers, & ! in
- ! num_roof_layers,num_wall_layers,num_road_layers, & ! in
- restart,sf_urban_physics, & !in
- XXXR_URB2D,XXXB_URB2D,XXXG_URB2D,XXXC_URB2D, & ! inout
- TR_URB2D,TB_URB2D,TG_URB2D,TC_URB2D,QC_URB2D, & ! inout
- TRL_URB3D,TBL_URB3D,TGL_URB3D, & ! inout
- SH_URB2D,LH_URB2D,G_URB2D,RN_URB2D, & ! inout
- TS_URB2D, & ! inout
- num_urban_layers, & ! in
- TRB_URB4D,TW1_URB4D,TW2_URB4D,TGB_URB4D, & ! inout
- TLEV_URB3D,QLEV_URB3D, & ! inout
- TW1LEV_URB3D,TW2LEV_URB3D, & ! inout
- TGLEV_URB3D,TFLEV_URB3D, & ! inout
- SF_AC_URB3D,LF_AC_URB3D,CM_AC_URB3D, & ! inout
- SFVENT_URB3D,LFVENT_URB3D, & ! inout
- SFWIN1_URB3D,SFWIN2_URB3D, & ! inout
- SFW1_URB3D,SFW2_URB3D,SFR_URB3D,SFG_URB3D, & ! inout
- A_U_BEP,A_V_BEP,A_T_BEP,A_Q_BEP, & ! inout multi-layer urban
- A_E_BEP,B_U_BEP,B_V_BEP, & ! inout multi-layer urban
- B_T_BEP,B_Q_BEP,B_E_BEP,DLG_BEP, & ! inout multi-layer urban
- DL_U_BEP,SF_BEP,VL_BEP, & ! inout multi-layer urban
- FRC_URB2D, UTYPE_URB2D) ! inout
- IMPLICIT NONE
- INTEGER, INTENT(IN) :: ISURBAN, sf_urban_physics
- INTEGER, INTENT(IN) :: ims,ime,jms,jme,kms,kme,num_soil_layers
- INTEGER, INTENT(IN) :: num_urban_layers !multi-layer urban
- ! INTEGER, INTENT(IN) :: num_roof_layers, num_wall_layers, num_road_layers
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: TSURFACE0_URB
- REAL, DIMENSION( ims:ime, 1:num_soil_layers, jms:jme ), INTENT(IN) :: TLAYER0_URB
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: TDEEP0_URB
- INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: IVGTYP
- LOGICAL , INTENT(IN) :: restart
-
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TR_URB2D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TB_URB2D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TG_URB2D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TC_URB2D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QC_URB2D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXR_URB2D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXB_URB2D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXG_URB2D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXC_URB2D
- ! REAL, DIMENSION(ims:ime, 1:num_roof_layers, jms:jme), INTENT(INOUT) :: TRL_URB3D
- ! REAL, DIMENSION(ims:ime, 1:num_wall_layers, jms:jme), INTENT(INOUT) :: TBL_URB3D
- ! REAL, DIMENSION(ims:ime, 1:num_road_layers, jms:jme), INTENT(INOUT) :: TGL_URB3D
- REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TRL_URB3D
- REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TBL_URB3D
- REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TGL_URB3D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SH_URB2D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LH_URB2D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: G_URB2D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: RN_URB2D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TS_URB2D
- ! multi-layer UCM variables
- REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: TRB_URB4D
- REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: TW1_URB4D
- REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: TW2_URB4D
- REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: TGB_URB4D
- REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: TLEV_URB3D
- REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: QLEV_URB3D
- REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: TW1LEV_URB3D
- REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: TW2LEV_URB3D
- REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: TGLEV_URB3D
- REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: TFLEV_URB3D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LF_AC_URB3D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SF_AC_URB3D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CM_AC_URB3D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SFVENT_URB3D
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LFVENT_URB3D
- REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: SFWIN1_URB3D
- REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: SFWIN2_URB3D
- REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: SFW1_URB3D
- REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: SFW2_URB3D
- REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: SFR_URB3D
- REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: SFG_URB3D
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_U_BEP
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_V_BEP
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_T_BEP
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_Q_BEP
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_E_BEP
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_U_BEP
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_V_BEP
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_T_BEP
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_Q_BEP
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_E_BEP
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: VL_BEP
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: DLG_BEP
- REAL, DIMENSION(ims:ime, kms:kme,jms:jme),INTENT(INOUT) :: SF_BEP
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: DL_U_BEP
- !
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FRC_URB2D
- INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UTYPE_URB2D
- INTEGER :: UTYPE_URB
- INTEGER :: I,J,K
- DO I=ims,ime
- DO J=jms,jme
- ! XXXR_URB2D(I,J)=0.
- ! XXXB_URB2D(I,J)=0.
- ! XXXG_URB2D(I,J)=0.
- ! XXXC_URB2D(I,J)=0.
- SH_URB2D(I,J)=0.
- LH_URB2D(I,J)=0.
- G_URB2D(I,J)=0.
- RN_URB2D(I,J)=0.
-
- !m
- FRC_URB2D(I,J)=0.
- UTYPE_URB2D(I,J)=0
- IF( IVGTYP(I,J) == ISURBAN) THEN
- UTYPE_URB2D(I,J) = 2 ! for default. high-intensity
- UTYPE_URB = UTYPE_URB2D(I,J) ! for default. high-intensity
- FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB)
- ENDIF
- IF( IVGTYP(I,J) == 31) THEN
- UTYPE_URB2D(I,J) = 3 ! low-intensity residential
- UTYPE_URB = UTYPE_URB2D(I,J) ! low-intensity residential
- FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB)
- ENDIF
- IF( IVGTYP(I,J) == 32) THEN
- UTYPE_URB2D(I,J) = 2 ! high-intensity
- UTYPE_URB = UTYPE_URB2D(I,J) ! high-intensity
- FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB)
- ENDIF
- IF( IVGTYP(I,J) == 33) THEN
- UTYPE_URB2D(I,J) = 1 ! Commercial/Industrial/Transportation
- UTYPE_URB = UTYPE_URB2D(I,J) ! Commercial/Industrial/Transportation
- FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB)
- ENDIF
- QC_URB2D(I,J)=0.01
- IF (.not.restart) THEN
- XXXR_URB2D(I,J)=0.
- XXXB_URB2D(I,J)=0.
- XXXG_URB2D(I,J)=0.
- XXXC_URB2D(I,J)=0.
- TC_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
- TR_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
- TB_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
- TG_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
- !
- TS_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
- ! DO K=1,num_roof_layers
- ! DO K=1,num_soil_layers
- ! TRL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
- ! TRL_URB3D(I,2,J)=TLAYER0_URB(I,2,J)+0.
- ! TRL_URB3D(I,3,J)=TLAYER0_URB(I,3,J)+0.
- ! TRL_URB3D(I,4,J)=TLAYER0_URB(I,4,J)+0.
- TRL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
- TRL_URB3D(I,2,J)=0.5*(TLAYER0_URB(I,1,J)+TLAYER0_URB(I,2,J))
- TRL_URB3D(I,3,J)=TLAYER0_URB(I,2,J)+0.
- TRL_URB3D(I,4,J)=TLAYER0_URB(I,2,J)+(TLAYER0_URB(I,3,J)-TLAYER0_URB(I,2,J))*0.29
- ! END DO
- ! DO K=1,num_wall_layers
- ! DO K=1,num_soil_layers
- !m TBL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
- !m TBL_URB3D(I,2,J)=TLAYER0_URB(I,2,J)+0.
- !m TBL_URB3D(I,3,J)=TLAYER0_URB(I,3,J)+0.
- !m TBL_URB3D(I,4,J)=TLAYER0_URB(I,4,J)+0.
- TBL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
- TBL_URB3D(I,2,J)=0.5*(TLAYER0_URB(I,1,J)+TLAYER0_URB(I,2,J))
- TBL_URB3D(I,3,J)=TLAYER0_URB(I,2,J)+0.
- TBL_URB3D(I,4,J)=TLAYER0_URB(I,2,J)+(TLAYER0_URB(I,3,J)-TLAYER0_URB(I,2,J))*0.29
- ! END DO
- ! DO K=1,num_road_layers
- DO K=1,num_soil_layers
- TGL_URB3D(I,K,J)=TLAYER0_URB(I,K,J)+0.
- END DO
-
- ! multi-layer urban
- ! IF( sf_urban_physics .EQ. 2)THEN
- IF((SF_URBAN_PHYSICS.eq.2).OR.(SF_URBAN_PHYSICS.eq.3)) THEN
- DO k=1,num_urban_layers
- ! TRB_URB4D(I,k,J)=TSURFACE0_URB(I,J)
- ! TW1_URB4D(I,k,J)=TSURFACE0_URB(I,J)
- ! TW2_URB4D(I,k,J)=TSURFACE0_URB(I,J)
- ! TGB_URB4D(I,k,J)=TSURFACE0_URB(I,J)
- !MT TRB_URB4D(I,K,J)=tlayer0_urb(I,1,J)
- !MT TW1_URB4D(I,K,J)=tlayer0_urb(I,1,J)
- !MT TW2_URB4D(I,K,J)=tlayer0_urb(I,1,J)
- IF (UTYPE_URB2D(I,J) > 0) THEN
- TRB_URB4D(I,K,J)=TBLEND_TBL(UTYPE_URB2D(I,J))
- TW1_URB4D(I,K,J)=TBLEND_TBL(UTYPE_URB2D(I,J))
- TW2_URB4D(I,K,J)=TBLEND_TBL(UTYPE_URB2D(I,J))
- ELSE
- TRB_URB4D(I,K,J)=tlayer0_urb(I,1,J)
- TW1_URB4D(I,K,J)=tlayer0_urb(I,1,J)
- TW2_URB4D(I,K,J)=tlayer0_urb(I,1,J)
- ENDIF
- TGB_URB4D(I,K,J)=tlayer0_urb(I,1,J)
- SFW1_URB3D(I,K,J)=0.
- SFW2_URB3D(I,K,J)=0.
- SFR_URB3D(I,K,J)=0.
- SFG_URB3D(I,K,J)=0.
- ENDDO
-
- ENDIF
-
- if (SF_URBAN_PHYSICS.EQ.3) then
- LF_AC_URB3D(I,J)=0.
- SF_AC_URB3D(I,J)=0.
- CM_AC_URB3D(I,J)=0.
- SFVENT_URB3D(I,J)=0.
- LFVENT_URB3D(I,J)=0.
- DO K=1,num_urban_layers
- TLEV_URB3D(I,K,J)=tlayer0_urb(I,1,J)
- TW1LEV_URB3D(I,K,J)=tlayer0_urb(I,1,J)
- TW2LEV_URB3D(I,K,J)=tlayer0_urb(I,1,J)
- TGLEV_URB3D(I,K,J)=tlayer0_urb(I,1,J)
- TFLEV_URB3D(I,K,J)=tlayer0_urb(I,1,J)
- QLEV_URB3D(I,K,J)=0.01
- SFWIN1_URB3D(I,K,J)=0.
- SFWIN2_URB3D(I,K,J)=0.
- !rm LF_AC_URB3D(I,J)=0.
- !rm SF_AC_URB3D(I,J)=0.
- !rm CM_AC_URB3D(I,J)=0.
- !rm SFVENT_URB3D(I,J)=0.
- !rm LFVENT_URB3D(I,J)=0.
- ENDDO
- endif
- ! IF( sf_urban_physics .EQ. 2 )THEN
- IF((SF_URBAN_PHYSICS.eq.2).OR.(SF_URBAN_PHYSICS.eq.3)) THEN
- DO K= KMS,KME
- SF_BEP(I,K,J)=1.
- VL_BEP(I,K,J)=1.
- A_U_BEP(I,K,J)=0.
- A_V_BEP(I,K,J)=0.
- A_T_BEP(I,K,J)=0.
- A_E_BEP(I,K,J)=0.
- A_Q_BEP(I,K,J)=0.
- B_U_BEP(I,K,J)=0.
- B_V_BEP(I,K,J)=0.
- B_T_BEP(I,K,J)=0.
- B_E_BEP(I,K,J)=0.
- B_Q_BEP(I,K,J)=0.
- DLG_BEP(I,K,J)=0.
- DL_U_BEP(I,K,J)=0.
- END DO
- ENDIF !sf_urban_physics=2
- ENDIF !restart
- END DO
- END DO
- RETURN
- END SUBROUTINE urban_var_init
- !===========================================================================
- !
- ! force_restore
- !
- !===========================================================================
- SUBROUTINE force_restore(CAP,AKS,DELT,S,R,H,LE,TSLEND,TSP,TS)
- REAL, INTENT(IN) :: CAP,AKS,DELT,S,R,H,LE,TSLEND,TSP
- REAL, INTENT(OUT) :: TS
- REAL :: C1,C2
- C2=24.*3600./2./3.14159
- C1=SQRT(0.5*C2*CAP*AKS)
- TS = TSP + DELT*( (S+R-H-LE)/C1 -(TSP-TSLEND)/C2 )
- END SUBROUTINE force_restore
- !===========================================================================
- !
- ! bisection (not used)
- !
- !==============================================================================
- SUBROUTINE bisection(TSP,PS,S,EPS,RX,SIG,RHO,CP,CH,UA,QA,TA,EL,BET,AKS,TSL,DZ,TS)
- REAL, INTENT(IN) :: TSP,PS,S,EPS,RX,SIG,RHO,CP,CH,UA,QA,TA,EL,BET,AKS,TSL,DZ
- REAL, INTENT(OUT) :: TS
- REAL :: ES,QS0,R,H,ELE,G0,F1,F
- TS1 = TSP - 5.
- TS2 = TSP + 5.
- DO ITERATION = 1,22
- ES=6.11*EXP( (2.5*10.**6./461.51)*(TS1-273.15)/(273.15*TS1) )
- QS0=0.622*ES/(PS-0.378*ES)
- R=EPS*(RX-SIG*(TS1**4.)/60.)
- H=RHO*CP*CH*UA*(TS1-TA)*100.
- ELE=RHO*EL*CH*UA*BET*(QS0-QA)*100.
- G0=AKS*(TS1-TSL)/(DZ/2.)
- F1= S + R - H - ELE - G0
- TS=0.5*(TS1+TS2)
- ES=6.11*EXP( (2.5*10.**6./461.51)*(TS-273.15)/(273.15*TS) )
- QS0=0.622*ES/(PS-0.378*ES)
- R=EPS*(RX-SIG*(TS**4.)/60.)
- H=RHO*CP*CH*UA*(TS-TA)*100.
- ELE=RHO*EL*CH*UA*BET*(QS0-QA)*100.
- G0=AKS*(TS-TSL)/(DZ/2.)
- F = S + R - H - ELE - G0
- IF (F1*F > 0.0) THEN
- TS1=TS
- ELSE
- TS2=TS
- END IF
- END DO
- RETURN
- END SUBROUTINE bisection
- !===========================================================================
- SUBROUTINE SFCDIF_URB (ZLM,Z0,THZ0,THLM,SFCSPD,AKANDA,AKMS,AKHS,RLMO,CD)
- ! ----------------------------------------------------------------------
- ! SUBROUTINE SFCDIF_URB (Urban version of SFCDIF_off)
- ! ----------------------------------------------------------------------
- ! CALCULATE SURFACE LAYER EXCHANGE COEFFICIENTS VIA ITERATIVE PROCESS.
- ! SEE CHEN ET AL (1997, BLM)
- ! ----------------------------------------------------------------------
- IMPLICIT NONE
- REAL WWST, WWST2, G, VKRM, EXCM, BETA, BTG, ELFC, WOLD, WNEW
- REAL PIHF, EPSU2, EPSUST, EPSIT, EPSA, ZTMIN, ZTMAX, HPBL, &
- & SQVISC
- REAL RIC, RRIC, FHNEU, RFC,RLMO_THR, RFAC, ZZ, PSLMU, PSLMS, PSLHU, &
- & PSLHS
- REAL XX, PSPMU, YY, PSPMS, PSPHU, PSPHS, ZLM, Z0, THZ0, THLM
- REAL SFCSPD, AKANDA, AKMS, AKHS, ZU, ZT, RDZ, CXCH
- REAL DTHV, DU2, BTGH, WSTAR2, USTAR, ZSLU, ZSLT, RLOGU, RLOGT
- REAL RLMO, ZETALT, ZETALU, ZETAU, ZETAT, XLU4, XLT4, XU4, XT4
- !CC ......REAL ZTFC
- REAL XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN, &
- & RLMA
- INTEGER ITRMX, ILECH, ITR
- REAL, INTENT(OUT) :: CD
- PARAMETER &
- & (WWST = 1.2,WWST2 = WWST * WWST,G = 9.8,VKRM = 0.40, &
- & EXCM = 0.001 &
- & ,BETA = 1./270.,BTG = BETA * G,ELFC = VKRM * BTG &
- & ,WOLD =.15,WNEW = 1. - WOLD,ITRMX = 05, &
- & PIHF = 3.14159265/2.)
- PARAMETER &
- & (EPSU2 = 1.E-4,EPSUST = 0.07,EPSIT = 1.E-4,EPSA = 1.E-8 &
- & ,ZTMIN = -5.,ZTMAX = 1.,HPBL = 1000.0 &
- & ,SQVISC = 258.2)
- PARAMETER &
- & (RIC = 0.183,RRIC = 1.0/ RIC,FHNEU = 0.8,RFC = 0.191 &
- & ,RLMO_THR = 0.001,RFAC = RIC / (FHNEU * RFC * RFC))
- ! ----------------------------------------------------------------------
- ! NOTE: THE TWO CODE BLOCKS BELOW DEFINE FUNCTIONS
- ! ----------------------------------------------------------------------
- ! LECH'S SURFACE FUNCTIONS
- ! ----------------------------------------------------------------------
- PSLMU (ZZ)= -0.96* log (1.0-4.5* ZZ)
- PSLMS (ZZ)= ZZ * RRIC -2.076* (1. -1./ (ZZ +1.))
- PSLHU (ZZ)= -0.96* log (1.0-4.5* ZZ)
- ! ----------------------------------------------------------------------
- ! PAULSON'S SURFACE FUNCTIONS
- ! ----------------------------------------------------------------------
- PSLHS (ZZ)= ZZ * RFAC -2.076* (1. -1./ (ZZ +1.))
- PSPMU (XX)= -2.* log ( (XX +1.)*0.5) - log ( (XX * XX +1.)*0.5) &
- & +2.* ATAN (XX) &
- &- PIHF
- PSPMS (YY)= 5.* YY
- PSPHU (XX)= -2.* log ( (XX * XX +1.)*0.5)
- ! ----------------------------------------------------------------------
- ! THIS ROUTINE SFCDIF CAN HANDLE BOTH OVER OPEN WATER (SEA, OCEAN) AND
- ! OVER SOLID SURFACE (LAND, SEA-ICE).
- ! ----------------------------------------------------------------------
- PSPHS (YY)= 5.* YY
- ! ----------------------------------------------------------------------
- ! ZTFC: RATIO OF ZOH/ZOM LESS OR EQUAL THAN 1
- ! C......ZTFC=0.1
- ! CZIL: CONSTANT C IN Zilitinkevich, S. S.1995,:NOTE ABOUT ZT
- ! ----------------------------------------------------------------------
- ILECH = 0
- ! ----------------------------------------------------------------------
- ! ZILFC = - CZIL * VKRM * SQVISC
- ! C.......ZT=Z0*ZTFC
- ZU = Z0
- RDZ = 1./ ZLM
- CXCH = EXCM * RDZ
- DTHV = THLM - THZ0
- ! ----------------------------------------------------------------------
- ! BELJARS CORRECTION OF USTAR
- ! ----------------------------------------------------------------------
- DU2 = MAX (SFCSPD * SFCSPD,EPSU2)
- !cc If statements to avoid TANGENT LINEAR problems near zero
- BTGH = BTG * HPBL
- IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
- WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
- ELSE
- WSTAR2 = 0.0
- END IF
- ! ----------------------------------------------------------------------
- ! ZILITINKEVITCH APPROACH FOR ZT
- ! ----------------------------------------------------------------------
- USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
- ! ----------------------------------------------------------------------
- ! KCL/TL Try Kanda approach instead (Kanda et al. 2007, JAMC)
- ! ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
- ZT = EXP (2.0-AKANDA*(SQVISC**2 * USTAR * Z0)**0.25)* Z0
-
- ZSLU = ZLM + ZU
- ZSLT = ZLM + ZT
- RLOGU = log (ZSLU / ZU)
- RLOGT = log (ZSLT / ZT)
- RLMO = ELFC * AKHS * DTHV / USTAR **3
- ! ----------------------------------------------------------------------
- ! 1./MONIN-OBUKKHOV LENGTH-SCALE
- ! ----------------------------------------------------------------------
- DO ITR = 1,ITRMX
- ZETALT = MAX (ZSLT * RLMO,ZTMIN)
- RLMO = ZETALT / ZSLT
- ZETALU = ZSLU * RLMO
- ZETAU = ZU * RLMO
- ZETAT = ZT * RLMO
- IF (ILECH .eq. 0) THEN
- IF (RLMO .lt. 0.0)THEN
- XLU4 = 1. -16.* ZETALU
- XLT4 = 1. -16.* ZETALT
- XU4 = 1. -16.* ZETAU
- XT4 = 1. -16.* ZETAT
- XLU = SQRT (SQRT (XLU4))
- XLT = SQRT (SQRT (XLT4))
- XU = SQRT (SQRT (XU4))
- XT = SQRT (SQRT (XT4))
- PSMZ = PSPMU (XU)
- SIMM = PSPMU (XLU) - PSMZ + RLOGU
- PSHZ = PSPHU (XT)
- SIMH = PSPHU (XLT) - PSHZ + RLOGT
- ELSE
- ZETALU = MIN (ZETALU,ZTMAX)
- ZETALT = MIN (ZETALT,ZTMAX)
- PSMZ = PSPMS (ZETAU)
- SIMM = PSPMS (ZETALU) - PSMZ + RLOGU
- PSHZ = PSPHS (ZETAT)
- SIMH = PSPHS (ZETALT) - PSHZ + RLOGT
- END IF
- ! ----------------------------------------------------------------------
- ! LECH'S FUNCTIONS
- ! ----------------------------------------------------------------------
- ELSE
- IF (RLMO .lt. 0.)THEN
- PSMZ = PSLMU (ZETAU)
- SIMM = PSLMU (ZETALU) - PSMZ + RLOGU
- PSHZ = PSLHU (ZETAT)
- SIMH = PSLHU (ZETALT) - PSHZ + RLOGT
- ELSE
- ZETALU = MIN (ZETALU,ZTMAX)
- ZETALT = MIN (ZETALT,ZTMAX)
- PSMZ = PSLMS (ZETAU)
- SIMM = PSLMS (ZETALU) - PSMZ + RLOGU
- PSHZ = PSLHS (ZETAT)
- SIMH = PSLHS (ZETALT) - PSHZ + RLOGT
- END IF
- ! ----------------------------------------------------------------------
- ! BELJAARS CORRECTION FOR USTAR
- ! ----------------------------------------------------------------------
- END IF
- USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
- !KCL/TL
- !ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
- ZT = EXP (2.0-AKANDA*(SQVISC**2 * USTAR * Z0)**0.25)* Z0
- ZSLT = ZLM + ZT
- RLOGT = log (ZSLT / ZT)
- USTARK = USTAR * VKRM
- AKMS = MAX (USTARK / SIMM,CXCH)
- AKHS = MAX (USTARK / SIMH,CXCH)
- !
- IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
- WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
- ELSE
- WSTAR2 = 0.0
- END IF
- !-----------------------------------------------------------------------
- RLMN = ELFC * AKHS * DTHV / USTAR **3
- !-----------------------------------------------------------------------
- ! IF(ABS((RLMN-RLMO)/RLMA).LT.EPSIT) GO TO 110
- !-----------------------------------------------------------------------
- RLMA = RLMO * WOLD+ RLMN * WNEW
- !-----------------------------------------------------------------------
- RLMO = RLMA
- END DO
- CD = USTAR*USTAR/SFCSPD**2
- ! ----------------------------------------------------------------------
- END SUBROUTINE SFCDIF_URB
- ! ----------------------------------------------------------------------
- !===========================================================================
- END MODULE module_sf_urban