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(:) :: CDS_TBL REAL, ALLOCATABLE, DIMENSION(:) :: AS_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 :: CAPR_DATA, CAPB_DATA, CAPG_DATA REAL :: AKSR_DATA, AKSB_DATA, AKSG_DATA REAL :: ALBR_DATA, ALBB_DATA, ALBG_DATA REAL :: EPSR_DATA, EPSB_DATA, EPSG_DATA REAL :: Z0R_DATA, Z0B_DATA, Z0G_DATA REAL :: Z0HR_DATA, Z0HB_DATA, Z0HG_DATA REAL :: TRLEND_DATA, TBLEND_DATA, TGLEND_DATA INTEGER :: BOUNDR_DATA,BOUNDB_DATA,BOUNDG_DATA INTEGER :: CH_SCHEME_DATA, TS_SCHEME_DATA 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 <-- urban_param.tbl ! |- urban_var_init ! ! Input Data from WRF [MKS unit]: ! ! UTYPE [-] : Urban type. 1=urban, 2=suburban, 3=rural ! 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: ! Following parameter are assigned in run/urban_param.tbl ! ! ZR [m] : roof level (building height) ! 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) ! SVF [-] : sky view factor. Calculated again in urban_param_init ! R [-] : building coverage ratio ! RW [-] : = 1 - R ! HGT [-] : normalized building height ! CDS [-] : drag coefficient by buildings ! AS [1/m] : buildings volumetric parameter ! AH [cal/cm/cm] : anthropogenic heat ! BETR [-] : minimum moisture availability of roof ! BETB [-] : minimum moisture availability of building wall ! BETG [-] : minimum moisture availability of road ! CAPR[cal/cm/cm/cm/degC]: heat capacity of roof ! CAPB[cal/cm/cm/cm/degC]: heat capacity of building wall ! CAPG[cal/cm/cm/cm/degC]: heat capacity of road ! AKSR [cal/cm/sec/degC] : thermal conductivity of roof ! AKSB [cal/cm/sec/degC] : thermal conductivity of building wall ! AKSG [cal/cm/sec/degC] : thermal conductivity of road ! 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 ! Z0R [m] : roughness length for momentum of roof ! 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) ! Z0HR [m] : roughness length for heat of roof ! Z0HB [m] : roughness length for heat of building wall (only for CH_SCHEME = 1) ! Z0HG [m] : roughness length for heat of roof ! 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 gound 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] ! ! ! 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 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 [urban=1, suburban=2, rural=3] 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] !------------------------------------------------------------------------------- ! 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, CDS, AS, AH REAL :: CAPR, CAPB, CAPG, AKSR, AKSB, AKSG, ALBR, ALBB, ALBG REAL :: EPSR, EPSB, EPSG, Z0R, Z0B, Z0G, Z0HR, Z0HB, Z0HG REAL :: TRLEND,TBLEND,TGLEND 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] !------------------------------------------------------------------------------- ! L: Local variables !------------------------------------------------------------------------------- REAL :: BETR, BETB, BETG REAL :: SX, SD, SQ, RX REAL :: UR, ZC, XLB, BB REAL :: Z, RIBR, 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 !------------------------------------------------------------------------------- ! Set parameters !------------------------------------------------------------------------------- CALL read_param(UTYPE,ZR,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,CDS,AS,AH, & CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG, & EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HR,Z0HB,Z0HG, & BETR,BETB,BETG,TRLEND,TBLEND,TGLEND, & BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME) IF( ZDC+Z0C+2. >= ZA) THEN PRINT *, 'ZDC + Z0C + 2m is larger than the 1st WRF level' PRINT *, 'Stop in the subroutine urban - change ZDC and Z0C' STOP 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 diffiusion 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 ! ZC=0.5*ZR XLB=0.4*(ZR-ZDC) BB=ZR*(CDS*AS/(2.*XLB**2.))**(1./3.) UC=UR*EXP(-BB*(1.-ZC/ZR)) ELSE print *,'ZR=',ZR, 'ZA=',ZA 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. 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) 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) 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 Fulxes from Urban Canopy !------------------------------------------------------------------------------- FLXUV = ( R*CDR + RW*CDC )*UA*UA FLXTH = ( R*FLXTHR + W*FLXTHB + RW*FLXTHG ) 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] RETURN 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, CAP, AKS, DELT,TSLEND INTEGER, INTENT(IN) :: KM, 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,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,CDS,AS,AH, & ! out CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG, & ! out EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HR,Z0HB,Z0HG, & ! out BETR,BETB,BETG,TRLEND,TBLEND,TGLEND, & ! out BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME) ! out INTEGER, INTENT(IN) :: UTYPE REAL, INTENT(OUT) :: ZR,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,CDS,AS,AH, & CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG, & EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HR,Z0HB,Z0HG, & BETR,BETB,BETG,TRLEND,TBLEND,TGLEND INTEGER, INTENT(OUT) :: BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME ZR = ZR_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) CDS= CDS_TBL(UTYPE) AS= AS_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_DATA CAPB= CAPB_DATA CAPG= CAPG_DATA AKSR= AKSR_DATA AKSB= AKSB_DATA AKSG= AKSG_DATA ALBR= ALBR_DATA ALBB= ALBB_DATA ALBG= ALBG_DATA EPSR= EPSR_DATA EPSB= EPSB_DATA EPSG= EPSG_DATA Z0R= Z0R_DATA Z0B= Z0B_DATA Z0G= Z0G_DATA Z0HR= Z0HR_DATA Z0HB= Z0HB_DATA Z0HG= Z0HG_DATA TRLEND= TRLEND_DATA TBLEND= TBLEND_DATA TGLEND= TGLEND_DATA BOUNDR= BOUNDR_DATA BOUNDB= BOUNDB_DATA BOUNDG= BOUNDG_DATA CH_SCHEME = CH_SCHEME_DATA TS_SCHEME = TS_SCHEME_DATA RETURN END SUBROUTINE read_param !=============================================================================== ! ! subroutine urban_param_init: Read parameters from urban_param.tbl ! !=============================================================================== SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers & ) ! 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 :: INDEX, 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 OPEN (UNIT=11, & FILE='urban_param.tbl', & ACCESS='SEQUENTIAL', & STATUS='OLD', & ACTION='READ', & POSITION='REWIND', & IOSTAT=IOSTATUS) IF (IOSTATUS > 0) STOP 'ERROR OPEN urban_param.tbl' READ(11,*) READ(11,'(A4)') LU_DATA_TYPE READ(11,*) ICATE ALLOCATE( ZR_TBL(ICATE), stat=allocate_status ) if(allocate_status == 0) THEN ALLOCATE( Z0C_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) stop 'error allocate Z0C_TBL in urban_param_init' IF( .NOT. ALLOCATED( Z0HC_TBL ) ) & ALLOCATE( Z0HC_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) stop 'error allocate Z0HC_TBL in urban_param_init' IF( .NOT. ALLOCATED( ZDC_TBL ) ) & ALLOCATE( ZDC_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) stop 'error allocate ZDC_TBL in urban_param_init' IF( .NOT. ALLOCATED( SVF_TBL ) ) & ALLOCATE( SVF_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) stop 'error allocate SVF_TBL in urban_param_init' IF( .NOT. ALLOCATED( R_TBL ) ) & ALLOCATE( R_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) stop 'error allocate R_TBL in urban_param_init' IF( .NOT. ALLOCATED( RW_TBL ) ) & ALLOCATE( RW_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) stop 'error allocate RW_TBL in urban_param_init' IF( .NOT. ALLOCATED( HGT_TBL ) ) & ALLOCATE( HGT_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) stop 'error allocate HGT_TBL in urban_param_init' IF( .NOT. ALLOCATED( CDS_TBL ) ) & ALLOCATE( CDS_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) stop 'error allocate CDS_TBL in urban_param_init' IF( .NOT. ALLOCATED( AS_TBL ) ) & ALLOCATE( AS_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) stop 'error allocate AS_TBL in urban_param_init' IF( .NOT. ALLOCATED( AH_TBL ) ) & ALLOCATE( AH_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) stop 'error allocate AH_TBL in urban_param_init' IF( .NOT. ALLOCATED( BETR_TBL ) ) & ALLOCATE( BETR_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) stop 'error allocate BETR_TBL in urban_param_init' IF( .NOT. ALLOCATED( BETB_TBL ) ) & ALLOCATE( BETB_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) stop 'error allocate BETB_TBL in urban_param_init' IF( .NOT. ALLOCATED( BETG_TBL ) ) & ALLOCATE( BETG_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) stop 'error allocate BETG_TBL in urban_param_init' ALLOCATE( FRC_URB_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) stop 'error allocate FRC_URB_TBL in urban_param_init' ENDIF DO LC = 1, ICATE READ(11,*) INDEX, & ZR_TBL(LC), & Z0C_TBL(LC), & Z0HC_TBL(LC), & ZDC_TBL(LC), & SVF_TBL(LC), & R_TBL(LC), & RW_TBL(LC), & HGT_TBL(LC), & CDS_TBL(LC), & AS_TBL(LC), & AH_TBL(LC), & BETR_TBL(LC), & BETB_TBL(LC), & BETG_TBL(LC), & FRC_URB_TBL(LC) END DO READ(11,*) READ(11,*) CAPR_DATA READ(11,*) READ(11,*) CAPB_DATA READ(11,*) READ(11,*) CAPG_DATA READ(11,*) READ(11,*) AKSR_DATA READ(11,*) READ(11,*) AKSB_DATA READ(11,*) READ(11,*) AKSG_DATA READ(11,*) READ(11,*) ALBR_DATA READ(11,*) READ(11,*) ALBB_DATA READ(11,*) READ(11,*) ALBG_DATA READ(11,*) READ(11,*) EPSR_DATA READ(11,*) READ(11,*) EPSB_DATA READ(11,*) READ(11,*) EPSG_DATA READ(11,*) READ(11,*) Z0R_DATA READ(11,*) READ(11,*) Z0B_DATA READ(11,*) READ(11,*) Z0G_DATA READ(11,*) READ(11,*) Z0HR_DATA READ(11,*) READ(11,*) Z0HB_DATA READ(11,*) READ(11,*) Z0HG_DATA READ(11,*) ! READ(11,*) num_roof_layers READ(11,*) dummy READ(11,*) ! READ(11,*) num_wall_layers READ(11,*) dummy READ(11,*) ! READ(11,*) num_road_layers READ(11,*) dummy num_roof_layers = num_soil_layers num_wall_layers = num_soil_layers num_road_layers = num_soil_layers DO K=1,num_roof_layers READ(11,*) READ(11,*) DZR(K) END DO DO K=1,num_wall_layers READ(11,*) READ(11,*) DZB(K) END DO DO K=1,num_road_layers READ(11,*) READ(11,*) DZG(K) END DO READ(11,*) READ(11,*) BOUNDR_DATA READ(11,*) READ(11,*) BOUNDB_DATA READ(11,*) READ(11,*) BOUNDG_DATA READ(11,*) READ(11,*) TRLEND_DATA READ(11,*) READ(11,*) TBLEND_DATA READ(11,*) READ(11,*) TGLEND_DATA READ(11,*) READ(11,*) CH_SCHEME_DATA READ(11,*) READ(11,*) TS_SCHEME_DATA CLOSE(11) ! Calculate Sky View Factor DO LC = 1, ICATE 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 END SUBROUTINE urban_param_init !=========================================================================== ! ! subroutine urban_var_init: initialization of urban state variables ! !=========================================================================== SUBROUTINE urban_var_init(TSURFACE0_URB,TLAYER0_URB,TDEEP0_URB,IVGTYP, & ! in ims,ime,jms,jme,num_soil_layers, & ! in ! num_roof_layers,num_wall_layers,num_road_layers, & ! 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, FRC_URB2D, UTYPE_URB2D) ! inout IMPLICIT NONE INTEGER, INTENT(IN) :: ims,ime,jms,jme,num_soil_layers ! 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 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 ! 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) == 1) THEN UTYPE_URB2D(I,J) = 2 ! for default. high-density UTYPE_URB = UTYPE_URB2D(I,J) ! for default. high-density FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB) ENDIF IF( IVGTYP(I,J) == 31) THEN UTYPE_URB2D(I,J) = 3 ! low-density residential UTYPE_URB = UTYPE_URB2D(I,J) ! low-density residential FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB) ENDIF IF( IVGTYP(I,J) == 32) THEN UTYPE_URB2D(I,J) = 2 ! high-density UTYPE_URB = UTYPE_URB2D(I,J) ! high-density 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 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 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 !=========================================================================== END MODULE module_sf_urban