estm_module Module Reference

Functions/Subroutines
subroutine	suews_getestmdata (lunit)
subroutine	estm_initials
subroutine	load_gridlayout (gridiv, multiplelayoutfiles, diagnose)
subroutine	estm_ehc_initialise
subroutine	estm_ehc_finalise
subroutine	estm_translate (gridiv)
subroutine	estm (gridiv, tstep, avkdn, avu1, temp_c, zenith_deg, avrh, press_hpa, ldown, bldgh, ts5mindata_ir, tair_av, dataoutlineestm, qs)
elemental real(kind(1d0)) function	set_nan (x)

Function/Subroutine Documentation

estm()

subroutine estm_module::estm	(	integer, intent(in)	gridiv,
		integer, intent(in)	tstep,
		real(kind(1d0)), intent(in)	avkdn,
		real(kind(1d0)), intent(in)	avu1,
		real(kind(1d0)), intent(in)	temp_c,
		real(kind(1d0)), intent(in)	zenith_deg,
		real(kind(1d0)), intent(in)	avrh,
		real(kind(1d0)), intent(in)	press_hpa,
		real(kind(1d0)), intent(in)	ldown,
		real(kind(1d0)), intent(in)	bldgh,
		real(kind(1d0)), dimension(ncolsestmdata), intent(in)	ts5mindata_ir,
		real(kind(1d0)), intent(in)	tair_av,
		real(kind(1d0)), dimension(27), intent(out)	dataoutlineestm,
		real(kind(1d0)), intent(out)	qs )

Definition at line 1667 of file suews_phys_estm.f95.

      ! NB: HCW Questions:
      !                - should TFloor be set in namelist instead of hard-coded here?
      !                - zref used for radiation calculation and fair is set to 2*BldgH here. For compatibility with the rest of the
      !                  SUEWS model, should this be the (wind speed) measurement height z specified in RunControl.nml?
      !                - In SUEWS_translate, fwall=AreaWall/SurfaceArea. Is this correct?
      !                - If froof=1 (i.e. whole grid is building), is HW=0 correct?
      !                - Then is an IF(Fground ==0) option needed?
      !                - alb_wall=0.23 and em_wall=0.9 are set in LUMPS_module_constants. Shouldn't these be provided as input?
      !                - Do the LUP calculations here need to be compatible with those in LUMPS_NARP?
      !                - File opening rewritten using existing error handling in SUEWS - can delete mod_error module from SUEWS_ESTM_functions
      !                - In SUEWS_ESTM_v2016, the first row is set to -999. This may be acceptable at the
      !                  start of the run but should be handled properly between blocks of met data? - need to check what's actually happening here.
      !                - Many duplicate functions in SUEWS_ESTM_functions need changing to the existing SUEWS versions.
      !                - Are the following correctly initialised? T0_ibld, T0_ground, T0_wall, T0_roof, TN_wall, TN_roof, Tground, Twall, Troof, Tibld
      !                - What are Nalb, sumalb, Nemis and Sumemis for? Are they used correctly?
      !
 
      !SUEWS_ESTM_v2016
      ! revision:
      ! HCW 14 Jun 2016
      ! HCW 27 Jun 2016 Corrected iESTMcount bug - now increases for all grids together
      ! HCW 30 Jun 2016 Major changes to handle grids and met blocks
      ! TS  09 Oct 2017 Added explicit interface
 
      !Contains calculation for each time step
      !Calculate local scale heat storage from single building and surroundings observations
      !OFferle, May 2003
      !
      !MODIFICATION HISTORY
      !             15 DECEMBER 2003
      !             (1) CHANGED AIR EXCHANGE RATE TO ALSO BE DEPENDENT ON OUTSIDE AIR TEMPERATURE
      !             (2) ADDED SPECIFIC HEAT CALCULATION FOR AIR
      !             (3) RH ADDED AS VAR(14) IN INPUT FILE
      !  12 JANUARY 2004
      !             (1) ADDED ESTIMATED AVG NADIR LOOKING EXTERNAL SURFACE TEMPERATURE TO OUTPUT
      !                 WEIGHTED BY SURFACE FRACTION AND SVF. SVF_ROOF=1, SVF_WALLS = 1-SVF_ground
      !             (2) ADDED NET RADIATION (RN) CALCULATIONS FOR ALL SURFACES AND AVG RN TO OUTPUT BASED ON
      !                 RADIOMETER VIEW FROM ZREF
      !       14 JANUARY 2004
      !             (1) MOVED GRID SPECIFIC PARAMETERS OUT OF NAMELIST, INTO PARAMETER FILE E.G. ALB,EM,F
      !                 T0 MAKE CHANGES IN LOCATIONS EASIER.
      !
      !       23 JANUARY 2004 - Version 2
      !             (1) PUT ALL TEMPERATURES INTO K
      !             (2) added interpolation routine to run on shorter timesteps.
      !                 tested so that it doesn't change solution for forced temperatures.
      !                 however in version 2 the energy balance at the surface isn't correctly solved
      !
      !       25 JANUARY 2004 - Version 3
      !             (1) added solution to energy balance at surface
      !             (2) removed some extraneous code
      !             (3) added vegetation fractions for future development
      !             (4) some changes to input namelist.
      !             (5) need to add wind speed dependence for exchange coefficients
      !             (6) changed the way Rn_net is calculated. also radiometer view factor relationships
      !             (7) added a calculation for heat loss/gain to outside air (that going into building mass is storage)
      !                 this is labelled QFBLD which it is in a sense.
      !       6 FEBRUARY 2004
      !             (1) ADDED INTERNAL VIEW FACTOR FILE FOR INTERNAL GEOMETRY, INCLUDING FIXED FLOOR TEMPERATURE
      !             (2) added MeanU, MinWS to config, and U to ILOC.
      !
      !      11 FEBRUARY 2004
      !             (1) added site lat, long, elevation to inputs
      !             (2) need zenith angle for wall direct radiation interception
      !
      !      15 JUNE 2004
      !             (1) CORRECTED OUTPUT OF T0 FOR FORCED SURFACE TEMPERATURES
      !             (2) MADE SOME CHANGES TO RADIATION, COMPUTATION OF AVERAGE ALBEDO, ADDED AVERAGE EMISSIVITY
      !             (3) ADDED OUTPUT FILE FOR RADIATION COMPONENTS
      !             (4) CHANGED INPUT CONFIG SO CONVECTIVE EXCHANGE COEFFS ARE NOT USER SELECTABLE
      !             (5) MAY STILL BE PROBLEMS WITH WALL RADIATIVE EXCHANGE AND AVERAGE ALBEDO
      !             (6) CHANGED THE WAY HEAT STORAGE IS COMPUTED IN HEATCONDUCTION MODULE BUT THIS SHOULD NOT CHANGE RESULTS
      !
      !      16 NOVEMBER 2004
      !             (1) CHANGED AIR EXCHANGE RATE TO BE BASES ON DAILY TEMPERATURE CHANGES.
      !             (2) WRITES OUT FINAL LAYER TEMPERATURES AND INCLUDES OPTION TO READ AT BEGINNING.
      !
      !DESCRIPTION: uses explicit time differencing to compute heat conduction through roof, walls,
      !             internal mass, and grounds (elements). Air heat storage is computed from average air temperature
      !             Boundary conditions are determined by measured surface temperature(s) or computed
      !             from energy balance at the surface.
      !             Internal air temperature can either be fixed or allowed to evolve in response
      !             to air mass exchanges and convective heating from internal surfaces.
      !
      !INPUT:
      !FORCING DATA: DTIME,KDOWN,LDOWN,TSURF,TAIR_OUT,TAIR_IN,TROOF,TWALL_AVG,TWALL_N,TWALL_E,TWALL_S,TWALL_W,Tground,RH,U
      !NAMELIST    : HEATSTORAGE_vFO.NML
      !    &config
      !    ifile=FORCING DATA
      !    ofile=OUTPUT FILE
      !    pfile=HEAT STORAGE PARAMETER FILE
      !    Nibld = INTERNAL MASS LAYERS
      !    Nwall = EXTERNAL WALL LAYERS
      !    Nroof = ROOF LAYERS
      !    Nground = ground/SOIL LAYERS
      !    LBC_soil = LOWER BOUNDARY CONDITION FOR ground/SOIL
      !    iloc= INPUT COLUMNS IN DATA FILE
      !    evolveTibld= USE DIAGNOSTIC VALUE FOR INTERNAL BUILDING TEMPERATURE
      !                        0: don't use, use measured
      !                        1: TURN ON USE when temp goess ABOVE TINT_ON, off when temp is below TINT_OFF
      !                        2: always use diagnostic
      !    THEAT_ON= TEMPERATURE AT WHICH HEAT CONTROL IS TURNED ON
      !    THEAT_OFF= TEMPERATURE AT WHICH HEAT CONTROL IS TURNED OFF
      !       THEAT_FIX = Fixed internal temperature for climate control
      !    oneTsurf= USE SINGLE SURFACE TEMPERATURE TO DRIVE ALL LAYERS
      !    radforce= USE RADIATIVE ENERGY BALANCE TO DRIVE EXTERNAL TEMPERATURES
      !    maxtimestep=302, maximum time step in s for filling data gaps.
      !       Alt = STATION HEIGHT (m) FOR PRESSURE CALCULATION
      !       SPINUP = NUMBER OF LINES TO USE FOR SPINUP (REPEATS THESE LINES BUT ONLY OUTPUTS THE 2ND TIME)
      !    INITTEMP = if TRUE INITIALIZES TEMPERATURES TO THOSE IN FINALTEMP.TXT FILE
      !    CH_ibld = INTERNAL BUILDING CONVECTIVE EXCHANGE COEFFICIENT
      !       **** THESE SHOULD DEPEND ON WIND SPEED BUT CURRENTLY DO NOT ****
      !       CHAIR = CONVECTIVE EXCHANGE COEFFICIENT FOR ROOF
      !       chair_ground = ... FOR ground
      !       chair_wall = ... FOR WALL
      !    /
      ! ***************** PARAMETER FILE VARIABLES
      !               fveg = FRACTION OF ground SURFACE COVERED WITH VEG
      !               zveg = VEGETATION HEIGHT
      !               alb_veg = VEGETATION ALBEDO
      !               em_veg = VEGETATION EMISSIVITY
      !               ZREF = REFERENCE HEIGHT FOR FLUX CALCULATION
      !               BldgH    = mean building height
      !               HW    = CANYON ASPECT RATION
      !               f_X   = FRACTION OF X WHERE X IS INTERNAL, WALL, ROOF, ground
      !               Alb_x = ALBEDO OF X
      !               em_ibld = EMISSIVITY OF X
      !               TX    = INITIAL LAYER TEMPERATURES
      !               zX    = LAYER THICKNESS
      !               kX    = LAYER THERMAL CONDUCTIVITY
      !               ribld = LAYER VOLUMETRIC HEAT CAPACITY
      !
      !****************** INTERNAL VIEW FACTOR FILE
      !OUTPUT:      fixed format text with single header, heatstorage for all elements, and temperatures
      !             for each element-layer.
      !===============================================================================
 
      USE meteo, ONLY: pi, heatcapacity_air
      USE mod_solver, ONLY: newtonpolynomial
      ! USE modSolarCalc !!FO!! :modsolarcalc.f95
      ! USE MathConstants !!FO!! :MathConstants_module.f95
      USE physconstants, ONLY: c2k, sbconst
      USE heatflux, ONLY: heatcond1d
      USE estm_data
 
      IMPLICIT NONE
      INTEGER, PARAMETER :: ncolsESTMdata = 13
      ! INTEGER, PARAMETER:: ncolumnsDataOutESTM=32
      INTEGER, PARAMETER :: cTs_Tiair = 5
      INTEGER, PARAMETER :: cTs_Tsurf = 6
      INTEGER, PARAMETER :: cTs_Troof = 7
      INTEGER, PARAMETER :: cTs_Troad = 8
      INTEGER, PARAMETER :: cTs_Twall = 9
      INTEGER, PARAMETER :: cTs_Twall_n = 10
      INTEGER, PARAMETER :: cTs_Twall_e = 11
      INTEGER, PARAMETER :: cTs_Twall_s = 12
      INTEGER, PARAMETER :: cTs_Twall_w = 13
      REAL(KIND(1D0)), PARAMETER :: NAN = -999
 
      INTEGER, INTENT(in) :: Gridiv
      INTEGER, INTENT(in) :: tstep
      ! INTEGER,INTENT(in)::iy !Year
      ! INTEGER,INTENT(in)::id !Day of year
      ! INTEGER,INTENT(in)::it !Hour
      ! INTEGER,INTENT(in)::imin          !Minutes
 
      REAL(KIND(1D0)), INTENT(in) :: avkdn
      REAL(KIND(1D0)), INTENT(in) :: avu1
      REAL(KIND(1D0)), INTENT(in) :: temp_c
      REAL(KIND(1D0)), INTENT(in) :: zenith_deg
      REAL(KIND(1D0)), INTENT(in) :: avrh
      REAL(KIND(1D0)), INTENT(in) :: press_hpa
      REAL(KIND(1D0)), INTENT(in) :: ldown
      REAL(KIND(1D0)), INTENT(in) :: bldgh
      ! REAL(KIND(1d0)),INTENT(in):: dectime        !Decimal time
      REAL(KIND(1D0)), DIMENSION(ncolsESTMdata), INTENT(in) :: Ts5mindata_ir !surface temperature input data
      REAL(KIND(1D0)), INTENT(in) :: Tair_av ! mean air temperature of past 24hr
 
      REAL(KIND(1D0)), DIMENSION(27), INTENT(out) :: dataOutLineESTM
      !Output to SUEWS
      REAL(KIND(1D0)), INTENT(out) :: QS
      !Input from SUEWS, corrected as Gridiv by TS 09 Jun 2016
 
      !Use only in this subroutine
      INTEGER :: i, ii
      INTEGER :: Tair2Set = 0
      REAL(KIND(1D0)) :: AIREXHR, AIREXDT
      REAL(KIND(1D0)), DIMENSION(2) :: bc
      REAL(KIND(1D0)) :: chair_ground, chair_wall
      REAL(KIND(1D0)) :: EM_EQUIV
      REAL(KIND(1D0)) :: kdz
      REAL(KIND(1D0)) :: kup_estm, LUP_net, kdn_estm
      REAL(KIND(1D0)) :: QHestm
      REAL(KIND(1D0)) :: QFBld !Anthropogenic heat from HVAC
      REAL(KIND(1D0)) :: shc_airbld
      REAL(KIND(1D0)) :: sw_hor, sw_vert
      REAL(KIND(1D0)) :: T0
      REAL(KIND(1D0)) :: Tinternal, Tsurf_all, Troof_in, Troad, Twall_all, Tw_n, Tw_e, Tw_s, Tw_w
      REAL(KIND(1D0)) :: Twallout(5), Troofout(5), Tibldout(5), Tgroundout(5)
      REAL(KIND(1D0)) :: Tadd, Tveg
      REAL(KIND(1D0)) :: Tairmix
      REAL(KIND(1D0)) :: RN
      REAL(KIND(1D0)) :: Rs_roof, Rl_roof, RN_ROOF
      REAL(KIND(1D0)) :: Rs_wall, Rl_wall, RN_WALL
      REAL(KIND(1D0)) :: Rs_ground, Rl_ground, RN_ground
      REAL(KIND(1D0)) :: Rs_ibld, Rl_ibld
      REAL(KIND(1D0)) :: Rs_iroof, Rl_iroof
      REAL(KIND(1D0)) :: Rs_iwall, Rl_iwall
      REAL(KIND(1D0)) :: zenith_rad
      REAL(KIND(1D0)) :: dum(50)
      REAL(KIND(1D0)) :: bldgHX ! local bldgHX=max(bldgH,0.001)
      REAL(KIND(1D0)), PARAMETER :: WSmin = 0.1 ! Check why there is this condition. S.O.
      LOGICAL :: radforce, groundradforce
 
      radforce = .false.
      groundradforce = .false. !Close the radiation scheme in original ESTM S.O.O.
 
      bldghx = max(bldgh, 0.001) ! this is to prevent zero building height
 
      ! Set -999s for first row
      ! dum=(/-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,&
      !      -999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,&
      !      -999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,&
      !      -999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,&
      !      -999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999./)
      dum = [(-999, i=1, 50)]
 
      !External bulk exchange coefficients - set these somewhere more sensible***
      chr = 0.005
      chair = chr
      chair_ground = chair
      chair_wall = chair
 
      !Get met data for use in ESTM subroutine
      kdn_estm = avkdn
      ws = avu1
      IF (ws < wsmin) ws = wsmin
      tair1 = temp_c + c2k
      ! Set initial value of Tair2 to air temp
      IF (gridiv == 1) tair2set = tair2set + 1
      IF (tair2set == 1) THEN
         tair2 = temp_c + c2k
      ELSE
         tair2 = tair2_grids(gridiv)
         ! Also get other variables for this grid
         tievolve = tievolve_grids(gridiv)
         lup_ground = lup_ground_grids(gridiv)
         lup_wall = lup_wall_grids(gridiv)
         lup_roof = lup_roof_grids(gridiv)
         t0_ibld = t0_ibld_grids(gridiv)
         t0_ground = t0_ground_grids(gridiv)
         t0_wall = t0_wall_grids(gridiv)
         t0_roof = t0_roof_grids(gridiv)
         tn_wall = tn_wall_grids(gridiv)
         tn_roof = tn_roof_grids(gridiv)
         tground(:) = tground_grids(:, gridiv)
         twall(:) = twall_grids(:, gridiv)
         troof(:) = troof_grids(:, gridiv)
         tibld(:) = tibld_grids(:, gridiv)
         tw_4 = tw_4_grids(:, :, gridiv)
 
      END IF
 
      ! Get Ts from Ts5min data array
      ! Tinternal  = Ts5mindata(ir,cTs_Tiair)
      ! Tsurf_all  = Ts5mindata(ir,cTs_Tsurf)
      ! Troof_in   = Ts5mindata(ir,cTs_Troof)
      ! Troad      = Ts5mindata(ir,cTs_Troad)
      ! Twall_all  = Ts5mindata(ir,cTs_Twall)
      !
      ! Tw_n       = Ts5mindata(ir,cTs_Twall_n)
      ! Tw_e       = Ts5mindata(ir,cTs_Twall_e)
      ! Tw_s       = Ts5mindata(ir,cTs_Twall_s)
      ! Tw_w       = Ts5mindata(ir,cTs_Twall_w)
      tinternal = ts5mindata_ir(cts_tiair)
      tsurf_all = ts5mindata_ir(cts_tsurf)
      troof_in = ts5mindata_ir(cts_troof)
      troad = ts5mindata_ir(cts_troad)
      twall_all = ts5mindata_ir(cts_twall)
 
      tw_n = ts5mindata_ir(cts_twall_n)
      tw_e = ts5mindata_ir(cts_twall_e)
      tw_s = ts5mindata_ir(cts_twall_s)
      tw_w = ts5mindata_ir(cts_twall_w)
 
      !    if (any(isnan(Ts5mindata))) then                                    !!FO!! can't use data when time gap is too big (or neg.) or data is NaN
      !        if (spindone) then                                                  !!FO!! writes a line of NaNs
      !            write(20,'(1F8.4,I6,100f10.1)') dectime,it,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,&
      !                -0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,&
      !                -0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,&
      !                -0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.
      !        endif
      !        return ! changed from cycle            !!FO!! returns to beginning of do-loop, i.e. doesn't make any heat calculations
      !    endif
      !!FO!! this loop is used to run through the calculations (SPINUP number of times) to achieve better numerical stability
      !!FO!! when if condition is met the program starts over from the beginning of the input file and the calculations are performed
      !!FO!! ...once again, but this time saved to the output file
      !!FO!! it's because of this if statement arrangement that the output result starts at input time #2 (if not SPINUP=0 as originally in heatstorage_vFO.nml)
      !    IF (NLINESREAD>datalines.AND..NOT.SPINDONE) THEN                        !!FO!! at this stage spindone = .false.
      !        NLINESREAD=0
      !        SPINDONE=.TRUE.
      !PRINT*, "SPUNUP"
      !        return ! changed from cycle
      !    ENDIF
 
      !! Write first row of each met block as -999
      !IF (first) THEN  !Set to true in ESTM_initials
      !!   !Tair2=Temp_C+C2K !This is now set in SUEWS_translate for ir=0 only
      !!   ! first=.FALSE.
      !   IF(Gridiv == NumberOfGrids) first=.FALSE.  !Set to false only after all grids have run
      !   dataOutESTM(ir,1:32,Gridiv)=(/REAL(iy,KIND(1D0)),REAL(id,KIND(1D0)),&
      !        REAL(it,KIND(1D0)),REAL(imin,KIND(1D0)),dectime,(dum(ii),ii=1,27)/)
      !   RETURN
      !ENDIF
 
      ! What are these constants? - Need defining somewhere
      zenith_rad = zenith_deg/180*pi
      IF (zenith_rad > 0 .AND. zenith_rad < pi/2.-hw) THEN !ZENITH MUST BE HIGHER THAN BUILDINGS FOR DIRECT INTERCEPTION
         tanzenith = min(tan(zenith_rad), 5.67) !LIMITS TO ANGLES LESS THAN 80 EVEN FOR LOW HW
         tanzenith = tanzenith*kdn_estm/(1370*cos(zenith_rad)) !REDUCTION FACTOR FOR MAXIMUM
      ELSE
         tanzenith = 0.
      END IF
 
      shc_air = heatcapacity_air(tair1, avrh, press_hpa) ! Use SUEWS version
 
      !Evolution of building temperature from heat added by convection
      SELECT CASE (evolvetibld) !EvolveTiBld specifies which internal building temperature approach to use
      CASE (0)
         diagnoseti = .false.
         hvac = .false. !use data in file                                !!FO!! use measured indoor temperature (Tref in Lodz2002HS.txt)
      CASE (1) !!FO!! use of HVAC to counteract T changes
         diagnoseti = .true.
         IF (tievolve > theat_off) THEN !THEAT_OFF now converted to Kelvin in ESTM_initials - HCW 15 Jun 2016
            !IF (Tievolve>THEAT_OFF+C2K) THEN
            hvac = .false.
         ELSEIF (tievolve < theat_on) THEN !THEAT_OFF now converted to Kelvin in ESTM_initials - HCW 15 Jun 2016
            !ELSEIF (Tievolve<THEAT_ON+C2K) THEN
            hvac = .true.
         END IF
      CASE (2)
         diagnoseti = .true. !!FO!! convection between ibld and inside of external walls(?)
      END SELECT
 
      !ASSUME AIR MIXES IN PROPORTION TO # OF EXCHANGES
      IF (tair_av > 20.+c2k .AND. tievolve > 25.+c2k .AND. tair1 < tievolve .AND. .NOT. hvac) THEN
         airexhr = 2.0 !Windows or exterior doors on 3 sides (ASHRAE 1981 22.8)
      ELSEIF (tair_av < 17.+c2k .OR. hvac) THEN
         airexhr = 0.5 !No window or exterior doors, storm sash or weathertripped (ASHRAE 1981 22.8)
      ELSE
         airexhr = 1.0
      END IF
 
      airexdt = airexhr*(tstep/3600.0)
      shc_airbld = max(heatcapacity_air(tievolve, avrh, press_hpa), 0.00001) ! to avoid zero-division scenario TS 21 Oct 2017
      IF (shc_airbld < minshc_airbld) minshc_airbld = shc_airbld
 
      !internal convective exchange coefficients                         !!FO!! ibldCHmod = 0 originally
      !iBldCHmod specifies method for convective exchange coeffs
      IF (ibldchmod == 1) THEN !ASHRAE 2001
         ch_ibld = 1.31*(abs(t0_ibld - tievolve))**0.25/shc_airbld
         ch_iwall = 1.31*(abs(tn_wall - tievolve))**0.25/shc_airbld
         ch_iroof = 1.52*(abs(tn_roof - tievolve))**0.25/shc_airbld
         IF (abs(tn_roof - tievolve) > 0) ch_iroof = ch_iroof*0.39 !effect of convection is weaker downward
      ELSEIF (ibldchmod == 2) THEN !Awbi, H.B. 1998, Energy and Buildings 28: 219-227
         ch_ibld = 1.823*(abs(t0_ibld - tievolve))**0.293/shc_airbld
         ch_iwall = 1.823*(abs(tn_wall - tievolve))**0.293/shc_airbld
         ch_iroof = 2.175*(abs(tn_roof - tievolve))**0.308/shc_airbld
         IF (abs(tn_roof - tievolve) > 0) ch_iroof = 0.704*(abs(tn_roof - tievolve))**0.133/shc_airbld !effect of convection is weaker downward
      END IF
 
      !Evolving T = (Previous Temp + dT from Sensible heat flux) mixed with outside air
      !ASSUMES THE CH_BLD INCLUDES THE EFFECT OF VENTILATION RATE IN m/s (e.g. if a normal CH is .005 and
      !the value here is .003 the assumed ventilation is 0.6 m/s                                 !!FO!! CH_ibld=0.0015 from heatstorage_Barbican.nml => ventilation=0.3 m/s
      tairmix = (tievolve + tair1*airexdt)/(1.0 + airexdt)
      qfbld = froof*(tievolve - tairmix)*shc_airbld*bldghx/tstep !heat added or lost, requires cooling or heating if HVAC on
 
      !!FO!! CH_xxxx has unit [m/s]  !!**HCW what is going on with tstep here??
      tievolve = tairmix + tstep/bldghx/finternal & !!FO!! finternal(=froof+fibld+fwall) => normalisation of fractions
                 *(ch_ibld*fibld*(t0_ibld - tievolve) &
                   + ch_iroof*froof*(tn_roof - tievolve) &
                   + ch_iwall*fwall*(tn_wall - tievolve)) !!FO!! [K] = [K] + [s/m]*([m/s]*([K]))
 
      IF (.NOT. diagnoseti) tievolve = tinternal + c2k
      IF (hvac) THEN !Run up/down to set point +/- 1 degree with adjustment of 90% per hour
         tadd = (sign(-1.0d0, theat_fix - tievolve) + theat_fix - tievolve)*min(4.*tstep/3600.0, 0.9) !!**HCW check??
         tievolve = tievolve + tadd
      END IF
 
      !========>RADIATION<================
      IF (kdn_estm < 0) kdn_estm = 0. !set non-zero shortwave to zero  !Should this be moved up to line 183/4?
 
      !external components, no diffuse
      !for reflections complete absorption is assumed
      !for shortwave these are net values
      !for longwave these are incoming only
      !MUST DIVIDE SHORTWAVE INTO DIRECT AND DIFFUSE
      sw_hor = kdn_estm !incoming solar on horizontal surface
      sw_vert = kdn_estm*tanzenith !incoming solar on vertical surface = kdown(obs)*sin(zenith)/cos(zenith)
 
      rs_roof = svf_roof*(1.0 - alb_roof_estm)*sw_hor
      rl_roof = svf_roof*em_roof_estm*ldown
 
      rs_ground = svf_ground*(1.-alb_ground_estm)*sw_hor + &
                  zvf_ground*svf_wall*alb_wall_fix*sw_vert*(1 - alb_ground_estm) + &
                  zvf_ground*svf_ground*alb_ground_estm*sw_hor*xvf_wall*alb_wall_fix
 
      rl_ground = svf_ground*ldown*em_ground_estm + zvf_ground*(lup_wall + svf_wall*ldown*(1 - em_wall_fix))*em_ground_estm
 
      rs_wall = svf_wall*(1.-alb_wall_fix)*sw_vert + &
                zvf_wall*svf_wall*alb_wall_fix*sw_vert*(1.+zvf_wall*alb_wall_fix) + &
                xvf_wall*svf_ground*alb_ground_estm*sw_hor*(1 - alb_wall_fix) + &
                zvf_ground*xvf_wall*svf_ground*alb_ground_estm*sw_hor*alb_wall_fix
 
      !wall to wall exchange handled simultaneously with seb calc
      rl_wall = svf_wall*ldown*em_wall_fix + zvf_wall*svf_wall*ldown*(1 - em_wall_fix)*em_wall_fix + &
                xvf_wall*(lup_ground + svf_ground*ldown*(1 - em_ground_estm))*em_wall_fix
 
      !DIFFICULT TO DETERMINE WHAT THIS IS EXACTLY, DONT INCLUDE WALLS
      kup_estm = kdn_estm - rvf_roof*rs_roof - (rvf_ground + rvf_wall)*rs_ground/svf_ground - rvf_veg*alb_veg_estm*kdn_estm
      IF (kdn_estm > 10 .AND. kup_estm > 0) THEN
         alb_avg = kup_estm/kdn_estm
         sumalb = sumalb + alb_avg
         nalb = nalb + 1
      END IF
 
      !internal components
      rs_ibld = 0 ! This could change if there are windows (need solar angles or wall svf * fraction glazing * transmissivity)
      !internal incoming longwave terms do not include the view factors for its own surface e.g. for ibld and walls
      !added floor view factors
      rl_ibld = sbconst*(ivf_iw*em_w*tn_wall**4 + &
                         ivf_ir*em_r*tn_roof**4 + &
                         ivf_if*em_f*tfloor**4)
      rs_iwall = 0
      rl_iwall = sbconst*(ivf_wi*em_i*t0_ibld**4 + &
                          ivf_wr*em_r*tn_roof**4 + &
                          ivf_wf*em_f*tfloor**4)
      rs_iroof = 0
      rl_iroof = sbconst*(ivf_ri*em_i*t0_ibld**4 + &
                          ivf_rw*em_w*tn_wall**4 + &
                          ivf_rf*em_f*tfloor**4)
 
      !========>INTERNAL<================
      bctype = .false.
      kdz = 2*kibld(1)/zibld(1)
      pcoeff = (/em_ibld*sbconst*(1 - ivf_ii*em_ibld), 0.0d0, 0.0d0, kdz + shc_airbld*ch_ibld, &
                 -kdz*tibld(1) - shc_airbld*ch_ibld*tievolve - rs_ibld - rl_ibld/)
      t0_ibld = newtonpolynomial(t0_ibld, pcoeff, conv, maxiter)
 
      !!FO!! this leads to Tibld(1) = Tibld(3) , i.e. ...
      bc(1) = t0_ibld
 
      !!FO!! temperature equal on both sides of inside wall
      bc(2) = bc(1)
 
      ! print*,'ESTM Qsibld cal before',Qsibld,Tibld,zibld(1:Nibld)
      CALL heatcond1d(tibld, qsibld, zibld(1:ndepth_ibld), &
                      REAL(Tstep, KIND(1D0)), kibld(1:Ndepth_ibld), &
                      ribld(1:ndepth_ibld), bc, bctype)
 
      ! print*,'ESTM Qsibld cal after',Qsibld,Tibld
 
      !========>WALLS<================
      bctype = .false.
      kdz = 2*kwall(ndepth_wall)/zwall(ndepth_wall)
      pcoeff = (/em_ibld*sbconst*(1 - ivf_ww*em_ibld), 0.0d0, 0.0d0, kdz + shc_airbld*ch_iwall, &
                 -kdz*twall(ndepth_wall) - shc_airbld*ch_iwall*tievolve - rs_iwall - rl_iwall/)
      tn_wall = newtonpolynomial(tn_wall, pcoeff, conv, maxiter)
      bc(2) = tn_wall !!FO!! boundary condition #2 = inner surface Twall, originally from lodz_parms_ltm.txt or finaltemp.txt
 
      IF (tsurfchoice < 2 .OR. radforce) THEN
         IF (radforce) THEN !!FO!! 1st prio: radforce
            kdz = 2*kwall(1)/zwall(1)
            pcoeff = (/em_wall_fix*sbconst*(1 - zvf_wall*em_wall_fix), 0.0d0, 0.0d0, kdz + shc_air*chair_wall*ws, &
                       -kdz*twall(1) - shc_air*chair_wall*ws*tair1 - rs_wall - rl_wall/)
            t0_wall = newtonpolynomial(t0_wall, pcoeff, conv, maxiter)
            bc(1) = t0_wall !!FO!! boundary condition #1 = outer surface Twall, originally from lodz_parms_ltm.txt or finaltemp.txt
         ELSEIF (tsurfchoice == 0) THEN
            bc(1) = tsurf_all + c2k; t0_wall = bc(1)
         ELSEIF (tsurfchoice == 1) THEN
            bc(1) = twall_all + c2k; t0_wall = bc(1)
         END IF !!FO!! Tsoil in Lodz2002HS.txt NB => Lodz2002HS.txt doesn't work with onewall = TRUE
 
         CALL heatcond1d(twall, qswall, zwall(1:ndepth_wall), real(tstep, kind(1d0)), &
                         kwall(1:ndepth_wall), rwall(1:ndepth_wall), bc, bctype) !!FO!! new set of Twalls are calculated from heat conduction through wall
 
      ELSEIF (tsurfchoice == 2) THEN !SPECIAL FOR 4 WALLS
         t0_wall = 0.
         DO i = 1, 4 !do 4 walls
            bc(1) = tw_n + tw_e + tw_s + tw_w + c2k; t0_wall = t0_wall + bc(1)
            CALL heatcond1d(tw_4(:, i), qs_4(i), zwall(1:ndepth_wall), real(tstep, kind(1d0)), &
                            kwall(1:ndepth_wall), rwall(1:ndepth_wall), bc, bctype)
         END DO
         !Take average of 4 wall values
         t0_wall = t0_wall/4.
         qswall = sum(qs_4)/4.
         twall = sum(tw_4, 2)/4.
      END IF
 
      !========>ROOF<================
      bctype = .false.
      kdz = 2*kroof(ndepth_roof)/zroof(ndepth_roof)
      pcoeff = (/em_ibld*sbconst, 0.0d0, 0.0d0, kdz + shc_airbld*ch_iroof, &
                 -kdz*troof(ndepth_roof) - shc_airbld*ch_iroof*tievolve - rs_iroof - rl_iroof/)
      tn_roof = newtonpolynomial(tn_roof, pcoeff, conv, maxiter)
      bc(2) = tn_roof
 
      IF (radforce) THEN
         kdz = 2*kroof(1)/zroof(1)
         pcoeff = (/em_roof_estm*sbconst, 0.0d0, 0.0d0, kdz + shc_air*chair*ws, &
                    -kdz*troof(1) - shc_air*chair*ws*tair1 - rs_roof - rl_roof/)
         t0_roof = newtonpolynomial(t0_roof, pcoeff, conv, maxiter)
         bc(1) = t0_roof
      ELSEIF (tsurfchoice == 0) THEN
         bc(1) = tsurf_all + c2k; t0_roof = bc(1)
      ELSE
         bc(1) = troof_in + c2k; t0_roof = bc(1)
      END IF
 
      CALL heatcond1d(troof, qsroof, zroof(1:ndepth_roof), real(tstep, kind(1d0)), &
                      kroof(1:ndepth_roof), rroof(1:ndepth_roof), bc, bctype)
 
      !========>ground<================
      bctype = .false.
      kdz = 2*kground(1)/zground(1)
 
      IF (radforce .OR. groundradforce) THEN
         pcoeff = (/em_ground_estm*sbconst, 0.0d0, 0.0d0, kdz + shc_air*chair_ground*ws, &
                    -kdz*tground(1) - shc_air*chair_ground*ws*tair1 - rs_ground - rl_ground/)
         t0_ground = newtonpolynomial(t0_ground, pcoeff, conv, maxiter)
         bc(1) = t0_ground
      ELSEIF (tsurfchoice == 0) THEN
         bc(1) = tsurf_all + c2k; t0_ground = bc(1)
      ELSE
         bc(1) = troad + c2k; t0_ground = bc(1)
      END IF
 
      bc(2) = lbc_soil + c2k
      !     bc(2)=0.; bctype(2)=.t.
 
      IF (fground /= 0.) THEN ! check fground==0 scenario to avoid division-by-zero error, TS 21 Jul 2016
         CALL heatcond1d(tground, qsground, zground(1:ndepth_ground), &
                         REAL(Tstep, KIND(1D0)), kground(1:Ndepth_ground), rground(1:Ndepth_ground), bc, bctype)
      ELSE
         qsground = nan
      END IF
 
      qsair = fair*shc_air*(tair1 - tair2)/tstep
      qsibld = qsibld*fibld
      qswall = qswall*fwall
      qsroof = qsroof*froof
      qsground = qsground*fground
      qs = qsibld + qswall + qsroof + qsground !!FO!! QSair not included; called QS in output file (column #10)
 
      ! print*,'ESTM QS',qs,Qsibld,Qswall,Qsroof ,Qsground
      ! print*,'ESTM Qsibld',Qsibld,fibld
      !write(*,*) Qsair, QSibld, Qswall, Qsroof, Qsground, QS
 
      !========>Radiation<================
      !note that the LUP for individual components does not include reflected
      lup_ground = sbconst*em_ground_estm*t0_ground**4
      lup_wall = sbconst*em_wall_fix*t0_wall**4
      lup_roof = sbconst*em_roof_estm*t0_roof**4
      tveg = tair1
      lup_veg = sbconst*em_veg_estm*tveg**4
      t0 = rvf_ground*t0_ground + rvf_wall*t0_wall + rvf_roof*t0_roof + rvf_veg*tveg
      lup_net = rvf_ground*lup_ground + rvf_wall*lup_wall + rvf_roof*lup_roof + rvf_veg*lup_veg
      em_equiv = lup_net/(sbconst*t0**4) !!FO!! apparent emissivity of the atmosphere [cloudless sky: >� Ldown from gases in the lowest 100 m] calculated from surface at T0
      rn_ground = rs_ground + rl_ground - lup_ground
      rn_roof = rs_roof + rl_roof - lup_roof
      rn_wall = rs_wall + rl_wall - lup_wall*(1 - zvf_wall*em_wall_fix)
      rn = kdn_estm - kup_estm + ldown*em_equiv - lup_net !!FO!! average net radiation (at z > zref ????) = shortwave down - shortwave up + [longwave down * apparent emissivity] - longwave up
      qhestm = (t0 - tair1)*chair*shc_air*ws
      sumemis = sumemis + em_equiv
      nemis = nemis + 1
 
      ! IF (SPINDONE) THEN                                                      !!FO!! only the last set of values in the time interpolation loop is written to file
 
      IF (ndepth_wall < 5) THEN
         twallout = (/twall, (dum(ii), ii=1, (5 - ndepth_wall))/)
      ELSE
         twallout = twall
      END IF
 
      IF (ndepth_roof < 5) THEN
         troofout = (/troof, (dum(ii), ii=1, (5 - ndepth_roof))/)
      ELSE
         troofout = troof
      END IF
 
      IF (ndepth_ground < 5) THEN
         tgroundout = (/tground, (dum(ii), ii=1, (5 - ndepth_ground))/)
      ELSE
         tgroundout = tground
      END IF
 
      IF (ndepth_ibld < 5) THEN
         tibldout = (/tibld, (dum(ii), ii=1, (5 - ndepth_ibld))/)
      ELSE
         tibldout = tibld
      END IF
 
      ! dataOutESTM(ir,1:ncolumnsDataOutESTM,Gridiv)=[&
      !      REAL(iy,KIND(1D0)),REAL(id,KIND(1D0)),REAL(it,KIND(1D0)),REAL(imin,KIND(1D0)), dectime,&!5
      !      QS,Qsair,Qswall,Qsroof,Qsground,Qsibld,&!11
      !      Twallout,Troofout,Tgroundout,Tibldout,Tievolve]!32 !NB. These all have 5 elements except Tievolve (1).
      dataoutlineestm = [ &
                        qs, qsair, qswall, qsroof, qsground, qsibld, & !6
                        twallout, troofout, tgroundout, tibldout, tievolve] !27 !NB. These all have 5 elements except Tievolve (1).
      ! set invalid values to nan
      dataoutlineestm = set_nan(dataoutlineestm)
      ! call r8vec_print(ncolumnsDataOutESTM-5,dataOutESTM(ir,6:ncolumnsDataOutESTM,Gridiv),'dataOutESTM')
 
      tair2 = tair1
 
      ! Save variables for this grid
      tair2_grids(gridiv) = tair1
      lup_ground_grids(gridiv) = lup_ground
      lup_wall_grids(gridiv) = lup_wall
      lup_roof_grids(gridiv) = lup_roof
      tievolve_grids(gridiv) = tievolve
      t0_ibld_grids(gridiv) = t0_ibld
      t0_ground_grids(gridiv) = t0_ground
      t0_wall_grids(gridiv) = t0_wall
      t0_roof_grids(gridiv) = t0_roof
      tn_wall_grids(gridiv) = tn_wall
      tn_roof_grids(gridiv) = tn_roof
      tground_grids(:, gridiv) = tground(:)
      twall_grids(:, gridiv) = twall(:)
      troof_grids(:, gridiv) = troof(:)
      tibld_grids(:, gridiv) = tibld(:)
      tw_4_grids(:, :, gridiv) = tw_4(:, :)
 

References estm_data::alb_avg, estm_data::alb_ground_estm, estm_data::alb_roof_estm, estm_data::alb_veg_estm, estm_data::alb_wall_fix, estm_data::bctype, physconstants::c2k, estm_data::ch_ibld, estm_data::ch_iroof, estm_data::ch_iwall, estm_data::chair, estm_data::chr, estm_data::conv, estm_data::diagnoseti, estm_data::em_f, estm_data::em_ground_estm, estm_data::em_i, estm_data::em_ibld, estm_data::em_r, estm_data::em_roof_estm, estm_data::em_veg_estm, estm_data::em_w, estm_data::em_wall_fix, estm_data::evolvetibld, estm_data::fair, estm_data::fground, estm_data::fibld, estm_data::finternal, estm_data::froof, estm_data::fwall, meteo::heatcapacity_air(), heatflux::heatcond1d(), estm_data::hvac, estm_data::hw, estm_data::ibldchmod, estm_data::ivf_if, estm_data::ivf_ii, estm_data::ivf_ir, estm_data::ivf_iw, estm_data::ivf_rf, estm_data::ivf_ri, estm_data::ivf_rw, estm_data::ivf_wf, estm_data::ivf_wi, estm_data::ivf_wr, estm_data::ivf_ww, estm_data::kground, estm_data::kibld, estm_data::kroof, estm_data::kwall, estm_data::lbc_soil, estm_data::lup_ground, estm_data::lup_ground_grids, estm_data::lup_roof, estm_data::lup_roof_grids, estm_data::lup_veg, estm_data::lup_wall, estm_data::lup_wall_grids, estm_data::maxiter, estm_data::minshc_airbld, estm_data::nalb, estm_data::ndepth_ground, estm_data::ndepth_ibld, estm_data::ndepth_roof, estm_data::ndepth_wall, estm_data::nemis, mod_solver::newtonpolynomial(), estm_data::pcoeff, estm_data::qs_4, estm_data::qsair, estm_data::qsground, estm_data::qsibld, estm_data::qsroof, estm_data::qswall, estm_data::rroof, estm_data::rvf_ground, estm_data::rvf_roof, estm_data::rvf_veg, estm_data::rvf_wall, estm_data::rwall, physconstants::sbconst, set_nan(), estm_data::shc_air, estm_data::sumalb, estm_data::sumemis, estm_data::svf_ground, estm_data::svf_roof, estm_data::svf_wall, estm_data::t0_ground, estm_data::t0_ground_grids, estm_data::t0_ibld, estm_data::t0_ibld_grids, estm_data::t0_roof, estm_data::t0_roof_grids, estm_data::t0_wall, estm_data::t0_wall_grids, estm_data::tair1, estm_data::tair2, estm_data::tair2_grids, estm_data::tanzenith, estm_data::tfloor, estm_data::tground, estm_data::tground_grids, estm_data::theat_fix, estm_data::theat_off, estm_data::theat_on, estm_data::tibld, estm_data::tibld_grids, estm_data::tievolve, estm_data::tievolve_grids, estm_data::tn_roof, estm_data::tn_roof_grids, estm_data::tn_wall, estm_data::tn_wall_grids, estm_data::troof, estm_data::troof_grids, estm_data::tsurfchoice, estm_data::tw_4, estm_data::tw_4_grids, estm_data::twall, estm_data::twall_grids, estm_data::ws, estm_data::xvf_wall, estm_data::zground, estm_data::zibld, estm_data::zroof, estm_data::zvf_ground, estm_data::zvf_wall, and estm_data::zwall.

Referenced by suews_driver::suews_cal_qs(), and suews_driver::suews_cal_qs_dts().

Here is the call graph for this function:

Here is the caller graph for this function:

estm_ehc_finalise()

subroutine estm_module::estm_ehc_finalise

Definition at line 1279 of file suews_phys_estm.f95.

      USE allocatearray
      IMPLICIT NONE
 
      IF (ALLOCATED(nlayer_grids)) DEALLOCATE (nlayer_grids)
 
      IF (ALLOCATED(height_grids)) DEALLOCATE (height_grids)
      IF (ALLOCATED(building_frac_grids)) DEALLOCATE (building_frac_grids)
      IF (ALLOCATED(veg_frac_grids)) DEALLOCATE (veg_frac_grids)
      IF (ALLOCATED(building_scale_grids)) DEALLOCATE (building_scale_grids)
      IF (ALLOCATED(veg_scale_grids)) DEALLOCATE (veg_scale_grids)
 
      IF (ALLOCATED(sfr_roof_grids)) DEALLOCATE (sfr_roof_grids)
      IF (ALLOCATED(alb_roof_grids)) DEALLOCATE (alb_roof_grids)
      IF (ALLOCATED(emis_roof_grids)) DEALLOCATE (emis_roof_grids)
      IF (ALLOCATED(k_roof_grids)) DEALLOCATE (k_roof_grids)
      IF (ALLOCATED(cp_roof_grids)) DEALLOCATE (cp_roof_grids)
      IF (ALLOCATED(dz_roof_grids)) DEALLOCATE (dz_roof_grids)
      IF (ALLOCATED(tin_roof_grids)) DEALLOCATE (tin_roof_grids)
      IF (ALLOCATED(temp_roof_grids)) DEALLOCATE (temp_roof_grids)
      IF (ALLOCATED(state_roof_grids)) DEALLOCATE (state_roof_grids)
      IF (ALLOCATED(statelimit_roof_grids)) DEALLOCATE (statelimit_roof_grids)
      IF (ALLOCATED(wetthresh_roof_grids)) DEALLOCATE (wetthresh_roof_grids)
      IF (ALLOCATED(soilstore_roof_grids)) DEALLOCATE (soilstore_roof_grids)
      IF (ALLOCATED(soilstorecap_roof_grids)) DEALLOCATE (soilstorecap_roof_grids)
      IF (ALLOCATED(roof_albedo_dir_mult_fact_grids)) DEALLOCATE (roof_albedo_dir_mult_fact_grids)
 
      IF (ALLOCATED(sfr_wall_grids)) DEALLOCATE (sfr_wall_grids)
      IF (ALLOCATED(alb_wall_grids)) DEALLOCATE (alb_wall_grids)
      IF (ALLOCATED(emis_wall_grids)) DEALLOCATE (emis_wall_grids)
      IF (ALLOCATED(k_wall_grids)) DEALLOCATE (k_wall_grids)
      IF (ALLOCATED(cp_wall_grids)) DEALLOCATE (cp_wall_grids)
      IF (ALLOCATED(dz_wall_grids)) DEALLOCATE (dz_wall_grids)
      IF (ALLOCATED(tin_wall_grids)) DEALLOCATE (tin_wall_grids)
      IF (ALLOCATED(temp_wall_grids)) DEALLOCATE (temp_wall_grids)
      IF (ALLOCATED(state_wall_grids)) DEALLOCATE (state_wall_grids)
      IF (ALLOCATED(statelimit_wall_grids)) DEALLOCATE (statelimit_wall_grids)
      IF (ALLOCATED(wetthresh_wall_grids)) DEALLOCATE (wetthresh_wall_grids)
      IF (ALLOCATED(soilstore_wall_grids)) DEALLOCATE (soilstore_wall_grids)
      IF (ALLOCATED(soilstorecap_wall_grids)) DEALLOCATE (soilstorecap_wall_grids)
      IF (ALLOCATED(wall_specular_frac_grids)) DEALLOCATE (wall_specular_frac_grids)
 
      IF (ALLOCATED(k_surf_grids)) DEALLOCATE (k_surf_grids)
      IF (ALLOCATED(cp_surf_grids)) DEALLOCATE (cp_surf_grids)
      IF (ALLOCATED(dz_surf_grids)) DEALLOCATE (dz_surf_grids)
      ! IF (ALLOCATED(tin_surf_grids)) DEALLOCATE (tin_surf_grids)
      IF (ALLOCATED(temp_surf_grids)) DEALLOCATE (temp_surf_grids)
 

References allocatearray::alb_roof_grids, allocatearray::alb_wall_grids, allocatearray::building_frac_grids, allocatearray::building_scale_grids, allocatearray::cp_roof_grids, allocatearray::cp_surf_grids, allocatearray::cp_wall_grids, allocatearray::dz_roof_grids, allocatearray::dz_surf_grids, allocatearray::dz_wall_grids, allocatearray::emis_roof_grids, allocatearray::emis_wall_grids, allocatearray::height_grids, allocatearray::k_roof_grids, allocatearray::k_surf_grids, allocatearray::k_wall_grids, allocatearray::nlayer_grids, allocatearray::roof_albedo_dir_mult_fact_grids, allocatearray::sfr_roof_grids, allocatearray::sfr_wall_grids, allocatearray::soilstore_roof_grids, allocatearray::soilstore_wall_grids, allocatearray::soilstorecap_roof_grids, allocatearray::soilstorecap_wall_grids, allocatearray::state_roof_grids, allocatearray::state_wall_grids, allocatearray::statelimit_roof_grids, allocatearray::statelimit_wall_grids, allocatearray::temp_roof_grids, allocatearray::temp_surf_grids, allocatearray::temp_wall_grids, allocatearray::tin_roof_grids, allocatearray::tin_wall_grids, allocatearray::veg_frac_grids, allocatearray::veg_scale_grids, allocatearray::wall_specular_frac_grids, allocatearray::wetthresh_roof_grids, and allocatearray::wetthresh_wall_grids.

estm_ehc_initialise()

subroutine estm_module::estm_ehc_initialise

Definition at line 1199 of file suews_phys_estm.f95.

 
      ! Last modified HCW 30 Jun 2016 - reading in now done by SUEWS_GetESTMData subroutine.
      !                                 ESTM_initials now only runs once per run at the very start.
      ! Last modified HCW 15 Jun 2016 - code now reads in 5-min file (interpolation done beforehand, outside of SUEWS itself)
 
      ! USE defaultNotUsed
      USE physconstants, ONLY: c2k
      ! USE ESTM_data
      USE allocatearray
      USE data_in, ONLY: multiplelayoutfiles, diagnose
      USE initial, ONLY: numberofgrids
 
      IMPLICIT NONE
      LOGICAL :: flag_mutiple_layout_files
      INTEGER :: i_grid
 
      IF (multiplelayoutfiles == 0) THEN
         flag_mutiple_layout_files = .false.
      ELSE IF (multiplelayoutfiles == 1) THEN
         flag_mutiple_layout_files = .true.
      END IF
 
      ALLOCATE (nlayer_grids(numberofgrids))
 
      ALLOCATE (height_grids(numberofgrids, nlayer_max))
      ALLOCATE (building_frac_grids(numberofgrids, nlayer_max))
      ALLOCATE (veg_frac_grids(numberofgrids, nlayer_max))
      ALLOCATE (building_scale_grids(numberofgrids, nlayer_max))
      ALLOCATE (veg_scale_grids(numberofgrids, nlayer_max))
 
      ! allocate roof variables
      ALLOCATE (tsfc_roof_grids(numberofgrids, nlayer_max))
      ALLOCATE (sfr_roof_grids(numberofgrids, nlayer_max))
      ALLOCATE (alb_roof_grids(numberofgrids, nlayer_max))
      ALLOCATE (emis_roof_grids(numberofgrids, nlayer_max))
      ALLOCATE (state_roof_grids(numberofgrids, nlayer_max))
      ALLOCATE (statelimit_roof_grids(numberofgrids, nlayer_max))
      ALLOCATE (wetthresh_roof_grids(numberofgrids, nlayer_max))
      ALLOCATE (soilstore_roof_grids(numberofgrids, nlayer_max))
      ALLOCATE (soilstorecap_roof_grids(numberofgrids, nlayer_max))
      ALLOCATE (k_roof_grids(numberofgrids, nlayer_max, ndepth))
      ALLOCATE (cp_roof_grids(numberofgrids, nlayer_max, ndepth))
      ALLOCATE (dz_roof_grids(numberofgrids, nlayer_max, ndepth))
      ALLOCATE (tin_roof_grids(numberofgrids, nlayer_max))
      ALLOCATE (temp_roof_grids(numberofgrids, nlayer_max, ndepth))
      ALLOCATE (roof_albedo_dir_mult_fact_grids(numberofgrids, nspec, nlayer_max))
 
      ! allocate wall variables
      ALLOCATE (tsfc_wall_grids(numberofgrids, nlayer_max))
      ALLOCATE (sfr_wall_grids(numberofgrids, nlayer_max))
      ALLOCATE (alb_wall_grids(numberofgrids, nlayer_max))
      ALLOCATE (emis_wall_grids(numberofgrids, nlayer_max))
      ALLOCATE (state_wall_grids(numberofgrids, nlayer_max))
      ALLOCATE (statelimit_wall_grids(numberofgrids, nlayer_max))
      ALLOCATE (wetthresh_wall_grids(numberofgrids, nlayer_max))
      ALLOCATE (soilstore_wall_grids(numberofgrids, nlayer_max))
      ALLOCATE (soilstorecap_wall_grids(numberofgrids, nlayer_max))
      ALLOCATE (k_wall_grids(numberofgrids, nlayer_max, ndepth))
      ALLOCATE (cp_wall_grids(numberofgrids, nlayer_max, ndepth))
      ALLOCATE (dz_wall_grids(numberofgrids, nlayer_max, ndepth))
      ALLOCATE (tin_wall_grids(numberofgrids, nlayer_max))
      ALLOCATE (temp_wall_grids(numberofgrids, nlayer_max, ndepth))
      ALLOCATE (wall_specular_frac_grids(numberofgrids, nspec, nlayer_max))
 
      ! allocate surf variables
      ! ALLOCATE (tsfc_surf_grids(NumberOfGrids, nsurf))
      ALLOCATE (k_surf_grids(numberofgrids, nsurf, ndepth))
      ALLOCATE (cp_surf_grids(numberofgrids, nsurf, ndepth))
      ALLOCATE (dz_surf_grids(numberofgrids, nsurf, ndepth))
      ALLOCATE (tin_surf_grids(numberofgrids, nsurf))
      ALLOCATE (temp_surf_grids(numberofgrids, nsurf, ndepth))
 
      DO i_grid = 1, numberofgrids
         print *, 'Reading layout data for grid ', i_grid
         CALL load_gridlayout(i_grid, flag_mutiple_layout_files, diagnose)
      END DO
 

References allocatearray::alb_roof_grids, allocatearray::alb_wall_grids, allocatearray::building_frac_grids, allocatearray::building_scale_grids, physconstants::c2k, allocatearray::cp_roof_grids, allocatearray::cp_surf_grids, allocatearray::cp_wall_grids, data_in::diagnose, allocatearray::dz_roof_grids, allocatearray::dz_surf_grids, allocatearray::dz_wall_grids, allocatearray::emis_roof_grids, allocatearray::emis_wall_grids, allocatearray::height_grids, allocatearray::k_roof_grids, allocatearray::k_surf_grids, allocatearray::k_wall_grids, load_gridlayout(), data_in::multiplelayoutfiles, allocatearray::ndepth, allocatearray::nlayer_grids, allocatearray::nlayer_max, allocatearray::nspec, allocatearray::nsurf, initial::numberofgrids, allocatearray::roof_albedo_dir_mult_fact_grids, allocatearray::sfr_roof_grids, allocatearray::sfr_wall_grids, allocatearray::soilstore_roof_grids, allocatearray::soilstore_wall_grids, allocatearray::soilstorecap_roof_grids, allocatearray::soilstorecap_wall_grids, allocatearray::state_roof_grids, allocatearray::state_wall_grids, allocatearray::statelimit_roof_grids, allocatearray::statelimit_wall_grids, allocatearray::temp_roof_grids, allocatearray::temp_surf_grids, allocatearray::temp_wall_grids, allocatearray::tin_roof_grids, allocatearray::tin_surf_grids, allocatearray::tin_wall_grids, allocatearray::tsfc_roof_grids, allocatearray::tsfc_wall_grids, allocatearray::veg_frac_grids, allocatearray::veg_scale_grids, allocatearray::wall_specular_frac_grids, allocatearray::wetthresh_roof_grids, and allocatearray::wetthresh_wall_grids.

Here is the call graph for this function:

estm_initials()

subroutine estm_module::estm_initials

Definition at line 781 of file suews_phys_estm.f95.

 
      ! Last modified HCW 30 Jun 2016 - reading in now done by SUEWS_GetESTMData subroutine.
      !                                 ESTM_initials now only runs once per run at the very start.
      ! Last modified HCW 15 Jun 2016 - code now reads in 5-min file (interpolation done beforehand, outside of SUEWS itself)
 
      USE defaultnotused
      USE physconstants, ONLY: c2k
      USE estm_data
      USE allocatearray
      USE data_in, ONLY: fileinputpath
      USE initial, ONLY: numberofgrids
 
      IMPLICIT NONE
 
      !=====Read ESTMinput.nml================================
      namelist /estminput/ tsurfchoice, &
         evolvetibld, &
         ibldchmod, &
         lbc_soil, &
         theat_on, &
         theat_off, &
         theat_fix
 
      OPEN (511, file=trim(fileinputpath)//'ESTMinput.nml', status='old')
      READ (511, nml=estminput)
      CLOSE (511)
 
      !Convert specified temperatures to Kelvin
      theat_on = theat_on + c2k
      theat_off = theat_off + c2k
      theat_fix = theat_fix + c2k
 
      ALLOCATE (tair2_grids(numberofgrids))
      ALLOCATE (lup_ground_grids(numberofgrids))
      ALLOCATE (lup_wall_grids(numberofgrids))
      ALLOCATE (lup_roof_grids(numberofgrids))
      ALLOCATE (tievolve_grids(numberofgrids))
      ALLOCATE (t0_ibld_grids(numberofgrids))
      ALLOCATE (t0_ground_grids(numberofgrids))
      ALLOCATE (t0_wall_grids(numberofgrids))
      ALLOCATE (t0_roof_grids(numberofgrids))
      ALLOCATE (tn_wall_grids(numberofgrids))
      ALLOCATE (tn_roof_grids(numberofgrids))
 

References physconstants::c2k, estm_data::evolvetibld, data_in::fileinputpath, estm_data::ibldchmod, estm_data::lbc_soil, estm_data::lup_ground_grids, estm_data::lup_roof_grids, estm_data::lup_wall_grids, initial::numberofgrids, estm_data::t0_ground_grids, estm_data::t0_ibld_grids, estm_data::t0_roof_grids, estm_data::t0_wall_grids, estm_data::tair2_grids, estm_data::theat_fix, estm_data::theat_off, estm_data::theat_on, estm_data::tievolve_grids, estm_data::tn_roof_grids, estm_data::tn_wall_grids, and estm_data::tsurfchoice.

estm_translate()

subroutine estm_module::estm_translate

(

integer

gridiv

)

Definition at line 1330 of file suews_phys_estm.f95.

      ! HCW 30 Jun 2016
 
      USE defaultnotused, ONLY: nan
      USE physconstants, ONLY: c2k, sbconst
      USE estm_data
      USE allocatearray
      USE gis_data, ONLY: bldgh
      USE initial, ONLY: numberofgrids
 
      IMPLICIT NONE
      INTEGER :: i
      !REAL(KIND(1d0)) :: CFLval
      !REAL(KIND(1d0)) :: t5min
      REAL(KIND(1D0)) :: W, WB
      !CHARACTER (len=20)::FileCodeX
      !CHARACTER (len=150):: FileFinalTemp
      !LOGICAL:: inittemps=.FALSE.
      INTEGER :: ESTMStart = 0
      INTEGER :: Gridiv
 
      !Set initial values at the start of each run for each grid
      ! the initiliastion part is problematic:
      ! cannot be initialised under a multi-grid scenario
      IF (gridiv == 1) estmstart = estmstart + 1
      ! print *, 'gridiv in ESTM',gridiv,'ESTMStart',ESTMStart
      IF (estmstart == 1) THEN
 
         ! write(*,*) ' ESTMStart: ',ESTMStart, 'initialising ESTM for grid no. ', Gridiv
 
         tfloor = 20.0 ! This is used only when radforce =T  !TODO:  should be put in the namelist
         tfloor = tfloor + c2k
 
         ! Initial values
         tievolve = 20.0 + c2k
         shc_air = 1230.0
         minshc_airbld = 1300
 
         ! ---- Internal view factors ----
         !constant now but should be calculated in the future
         ivf_iw = 0.100000
         ivf_ir = 0.000000
         ivf_ii = 0.900000
         ivf_if = 0.000000
         ivf_ww = 0.050000
         ivf_wr = 0.000000
         ivf_wi = 0.950000
         ivf_wf = 0.000000
         ivf_rw = 0.050000
         ivf_ri = 0.950000
         ivf_rf = 0.000000
         ivf_fw = 0.050000
         ivf_fr = 0.000000
         ivf_fi = 0.950000
 
         tair24hr = c2k
 
         !Ts5mindata(1,ncolsESTMdata) = -999
         ! !Fill Ts5mindata for current grid and met block - this is done in SUEWS_translate
         ts5mindata(1, 1:ncolsestmdata) = estmforcingdata(1, 1:ncolsestmdata, gridiv)
 
         ! ---- Initialization of variables and parameters for first row of run for each grid ----
         ! N layers are calculated in SUEWS_translate
         IF (.NOT. ALLOCATED(tibld)) THEN
            ! print *, "Nibld", Nibld
            ! print *, "Nwall", Nwall
            ! print *, "Nroof", Nroof
            ! print *, "Nground", Nground
            ALLOCATE (tibld(ndepth_ibld), twall(ndepth_wall), troof(ndepth_roof), tground(ndepth_ground), tw_4(ndepth_wall, 4))
            ALLOCATE (tibld_grids(ndepth_ibld, numberofgrids), &
                      twall_grids(ndepth_wall, numberofgrids), &
                      troof_grids(ndepth_roof, numberofgrids), &
                      tground_grids(ndepth_ground, numberofgrids), &
                      tw_4_grids(ndepth_wall, 4, numberofgrids))
         END IF
 
         ! Transfer variables from Ts5mindata to variable names
         ! N.B. column numbers here for the following file format - need to change if input columns change!
         ! dectime iy id it imin Tiair Tsurf Troof Troad Twall Twall_n Twall_e Twall_s Twall_w
         !        1  2  3  4    5     6     7     8     9     10      11      12      13       !new
         !
         ! Calculate temperature of each layer in Kelvin
         !  QUESTION: what if (Nground/Nwall/Nroof-1)==0? TS 21 Oct 2017
         DO i = 1, ndepth_ground
            tground(i) = (lbc_soil - ts5mindata(1, cts_troad))*(i - 1)/(ndepth_ground - 1) + ts5mindata(1, cts_troad) + c2k
         END DO
         DO i = 1, ndepth_wall
            twall(i) = &
               (ts5mindata(1, cts_tiair) - ts5mindata(1, cts_twall))*(i - 1)/(ndepth_wall - 1) + ts5mindata(1, cts_twall) &
               + c2k
         END DO
         DO i = 1, ndepth_roof
            troof(i) = &
               (ts5mindata(1, cts_tiair) - ts5mindata(1, cts_troof))*(i - 1)/(ndepth_roof - 1) + ts5mindata(1, cts_troof) &
               + c2k
         END DO
         tibld(1:ndepth_ibld) = ts5mindata(1, cts_tiair) + c2k
 
      END IF !End of loop run only at start (for each grid)
 
      ! ---- Parameters related to land surface characteristics ----
      ! QUESTION: Would Zref=z be more appropriate?
      zref = 2.0*bldgh !!FO!! BldgH: mean bulding hight, zref: local scale reference height (local: ~ 10^2 x 10^2 -- 10^3 x 10^3 m^2)
 
      svf_ground = 1.0
      svf_roof = 1.0
 
      ! ==== roof (i.e. Bldgs)
      !froof=sfr_surf(BldgSurf)   ! Moved to SUEWS_translate HCW 16 Jun 2016
      alb_roof_estm = alb(bldgsurf)
      em_roof_estm = emis(bldgsurf)
 
      ! ==== vegetation (i.e. EveTr, DecTr, Grass)
      !fveg=sfr_surf(ConifSurf)+sfr_surf(DecidSurf)+sfr_surf(GrassSurf)  ! Moved to SUEWS_translate HCW 16 Jun 2016
      IF (fveg /= 0) THEN
         alb_veg_estm = dot_product(alb([conifsurf, decidsurf, grasssurf]), sfr_surf([conifsurf, decidsurf, grasssurf]))/fveg
         em_veg_estm = dot_product(emis([conifsurf, decidsurf, grasssurf]), sfr_surf([conifsurf, decidsurf, grasssurf]))/fveg
      ELSE ! check fveg==0 scenario to avoid division-by-zero error, TS 21 Oct 2017
         alb_veg_estm = nan
         em_veg_estm = nan
      END IF
 
      ! ==== ground (i.e. Paved, EveTr, DecTr, Grass, BSoil, Water - all except Bldgs)
      !fground=sfr_surf(ConifSurf)+sfr_surf(DecidSurf)+sfr_surf(GrassSurf)+sfr_surf(PavSurf)+sfr_surf(BsoilSurf)+sfr_surf(WaterSurf) ! Moved to SUEWS_translate HCW 16 Jun 2016
      IF (fground /= 0) THEN
         alb_ground_estm = (alb(conifsurf)*sfr_surf(conifsurf) + alb(decidsurf)*sfr_surf(decidsurf) &
                            + alb(grasssurf)*sfr_surf(grasssurf) + alb(pavsurf)*sfr_surf(pavsurf) &
                            + alb(bsoilsurf)*sfr_surf(bsoilsurf) + alb(watersurf)*sfr_surf(watersurf))/fground
         em_ground_estm = (emis(conifsurf)*sfr_surf(conifsurf) + emis(decidsurf)*sfr_surf(decidsurf) &
                           + emis(grasssurf)*sfr_surf(grasssurf) + emis(pavsurf)*sfr_surf(pavsurf) &
                           + emis(bsoilsurf)*sfr_surf(bsoilsurf) + emis(watersurf)*sfr_surf(watersurf))/fground
      ELSE ! check fground==0 scenario to avoid division-by-zero error, TS 21 Jul 2016
         alb_ground_estm = nan
         em_ground_estm = nan
      END IF
 
      IF (froof < 1.0) THEN
         hw = fwall/(2.0*(1.0 - froof))
      ELSE
         !  HW=0  !HCW if only roof, no ground
         hw = 0.00001 ! to avoid zero-HW scenario TS 21 Oct 2017
 
      END IF
      hw = max(0.00001, hw) ! to avoid zero-HW scenario TS 27 Oct 2017
 
      IF (fground == 1.0) THEN !!FO!! if only ground, i.e. no houses
         w = 1
         wb = 0
         svf_ground = 1.
         zvf_wall = 0.
         svf_wall = 0.
         svf_roof = 1.
         zvf_ground = 0.
         xvf_wall = 0.
         rvf_canyon = 1.
         rvf_ground = 1.-fveg
         rvf_roof = 0
         rvf_wall = 0
         rvf_veg = fveg
      ELSE IF (fground == 0.0) THEN !check fground==0 (or HW==0) scenario to avoid division-by-zero error, TS 21 Jul 2016
         ! the following values are calculated given HW=0
         w = 0
         wb = 1
         zvf_wall = 0 !COS(ATAN(2/HW))  when HW=0                                 !!FO!! wall view factor for wall
         hw = 0
         svf_ground = max(cos(atan(2*hw)), 0.00001) !!FO!! sky view factor for ground ! to avoid zero-division scenario TS 21 Oct 2017
         svf_wall = (1 - zvf_wall)/2 !!FO!! sky view factor for wall
         zvf_ground = 1 - svf_ground !!FO!! wall view factor for ground
         xvf_wall = svf_wall !!FO!! ground view factor
         !   RVF_CANYON=COS(ATAN(2*ZREF/W))
         !   RVF_ROOF=1-RVF_CANYON
         !   RVF_WALL=(COS(ATAN(2*(ZREF-BldgH)/W))-RVF_CANYON)*RVF_CANYON
         !   RVF_ground=RVF_CANYON-RVF_WALL
         rvf_ground = (fground - fveg)*svf_ground
         rvf_veg = fveg*svf_ground
         rvf_roof = froof
         rvf_wall = 1 - rvf_roof - rvf_ground - rvf_veg
      ELSE
         w = bldgh/hw !What about if HW = 0, need to add IF(Fground ==0) option ! fixed by setting a small number, TS 21 Oct 2017
         wb = w*sqrt(froof/fground)
         svf_ground = cos(atan(2*hw)) !!FO!! sky view factor for ground
         zvf_wall = cos(atan(2/hw)) !!FO!! wall view factor for wall
         svf_wall = (1 - zvf_wall)/2 !!FO!! sky view factor for wall
         zvf_ground = 1 - svf_ground !!FO!! wall view factor for ground
         xvf_wall = svf_wall !!FO!! ground view factor
         !   RVF_CANYON=COS(ATAN(2*ZREF/W))
         !   RVF_ROOF=1-RVF_CANYON
         !   RVF_WALL=(COS(ATAN(2*(ZREF-BldgH)/W))-RVF_CANYON)*RVF_CANYON
         !   RVF_ground=RVF_CANYON-RVF_WALL
         rvf_ground = (fground - fveg)*svf_ground
         rvf_veg = fveg*svf_ground
         rvf_roof = froof
         rvf_wall = 1 - rvf_roof - rvf_ground - rvf_veg
      END IF
 
      alb_avg = alb_ground_estm*rvf_ground + alb_wall_fix*rvf_wall + alb_roof_estm*rvf_roof + alb_veg_estm*rvf_veg
 
      sumalb = 0.; nalb = 0
      sumemis = 0.; nemis = 0
 
      !set emissivity for ceiling, wall and floor inside of buildings
      em_r = em_ibld; em_w = em_ibld; em_i = em_ibld; em_f = em_ibld
 
      !internal elements
      IF (nroom == 0) THEN
         fibld = (floor(bldgh/3.1 - 0.5) - 1)*froof
      ELSE
         fibld = (2.-2./nroom)*fwall + (floor(bldgh/3.1 - 0.5) - 1)*froof
      END IF
 
      IF (fibld == 0) fibld = 0.00001 !this just ensures a solution to radiation
      finternal = froof + fibld + fwall
      IF (finternal == 0) finternal = 0.00001 ! to avoid zero-devision error TS 21 Oct 2017
      fair = zref - bldgh*froof
      IF (fair == 0) fair = 0.00001 ! to avoid zero-devision error TS 21 Oct 2017
      !ivf_ii=1.-ivf_iw-ivf_ir-ivf_if    !S.O. I do not know these are should be calculated or read from input files
      !ivf_ww=1.-ivf_wi-ivf_wr-ivf_wf
      !ivf_rw=1.-ivf_ri-ivf_rf;
      !ivf_fr=ivf_rf;
 
      IF ((ivf_ii + ivf_iw + ivf_ir + ivf_if > 1.0001) .OR. &
          (ivf_wi + ivf_ww + ivf_wr + ivf_wf > 1.0001) .OR. &
          (ivf_ri + ivf_rw + ivf_rf > 1.0001) .OR. &
          (ivf_fi + ivf_fw + ivf_fr > 1.0001) .OR. &
          (ivf_ii + ivf_iw + ivf_ir + ivf_if < 0.9999) .OR. &
          (ivf_wi + ivf_ww + ivf_wr + ivf_wf < 0.9999) .OR. &
          (ivf_ri + ivf_rw + ivf_rf < 0.9999) .OR. &
          (ivf_fi + ivf_fw + ivf_fr < 0.9999)) THEN
         print *, "At least one internal view factor <> 1. Check ivf in ESTMinput.nml"
      END IF
 
      !=======Initial setting==============================================
      !! Rewritten by HCW 15 Jun 2016 to use existing SUEWS error handling
      !IF(inittemps) THEN
      !   write(*,*) 'inittemps:',inittemps
      !   FileFinalTemp=TRIM(FileOutputPath)//TRIM(FileCodeX)//'_ESTM_finaltemp.txt'
      !   OPEN(99,file=TRIM(FileFinalTemp),status='old',err=316)  ! Program stopped if error opening file
      !   READ(99,*) Twall,Troof,Tground,Tibld                    ! Twall, Troof, Tground & Tibld get new values
      !   CLOSE(99)
      !ENDIF
      !
      !!IF (inittemps) THEN                                                        !!FO!! inittemps=.true. set in nml file
      !!   OPEN(99,file='outputfiles/finaltemp.txt',status='old',iostat=ios)       !!FO!! has to exist
      !!
      !!   IF (ios/=0) CALL error('outputfiles/finaltemp.txt',ios,1)               !!FO!! calls mod_error.f95, writes that the opening failed and stops prg
      !!   IF (ios/=0) THEN
      !!      Twall   = (/273., 285., 291./)
      !!      Troof   = (/273., 285., 291./)
      !!      Tground = (/273., 275., 280., 290./)
      !!      Tibld   = (/293., 293., 293./)
      !!   ELSE
      !!      READ(99,*) Twall,Troof,Tground,Tibld                             !!FO!! if finaltemp.txt exists Twall[3], Troof[3], Tground[4] & Tibld[3] get new values
      !!      CLOSE(99)
      !!   ENDIF
      !!ENDIF
 
      !where (isnan(Twall))
      !    Twall = 273
      !endwhere
      !where (isnan(Troof))
      !    Troof = 273
      !endwhere
      !where (isnan(Tground))
      !    Tground = 281
      !endwhere
      !where (isnan(Tibld))
      !    Tibld = 293
      !endwhere
 
      IF (estmstart == 1) THEN
         DO i = 1, 4
            tw_4(:, i) = twall !!FO!! Tw_4 holds three differnet temp:s for each wall layer but the same set for all points of the compass
         END DO
 
         !initialize surface temperatures
         t0_ground = tground(1)
         t0_wall = twall(1)
         t0_roof = troof(1)
         t0_ibld = tibld(1)
         tn_roof = troof(ndepth_roof)
         tn_wall = twall(ndepth_wall)
 
         !initialize outgoing longwave   !QUESTION: HCW - Are these calculations compatible with those in LUMPS_NARP?
         lup_ground = sbconst*em_ground_estm*t0_ground**4
         lup_wall = sbconst*em_wall_fix*t0_wall**4
         lup_roof = sbconst*em_roof_estm*t0_roof**4
 
         !  PRINT*,"W,WB= ",W,WB
         !  PRINT*,'SVF_ground ','SVF_WALL ','zvf_WALL ','HW '
         !  PRINT*,SVF_ground,SVF_WALL,zvf_WALL,HW
         !  PRINT*,'RVF_ground ','RVF_WALL ','RVF_ROOF ','RVF_VEG'
         !  PRINT*,RVF_ground,RVF_WALL,RVF_ROOF,RVF_VEG
         !  print*,'Alb_avg (VF)=',alb_avg
         !  print*,'z0m, Zd', z0m, ZD
 
      END IF
 
      first = .true.
 
      !======Courant-Friedrichs-Lewy condition=================================
      !This is comment out by S.O. for now
      ! NB: should be recovered for calculation stability, FO and TS, 11 Oct 2017
      !   CFLval = minval(0.5*zibld*zibld*ribld/kibld)   !!FO!! z*z*r/k => unit [s]
      !   if (Tstep>CFLval) then !CFL condition   !!FO!! CFL condition:  Courant�Friedrichs�Lewy condition is a necessary condition for convergence while solving
      !      write(*,*) "IBLD: CFL condition: Tstep=",Tstep,">",CFLval !!FO!! certain partial differential equations numerically by the method of finite differences (like eq 5 in Offerle et al.,2005)
      !      CFLfail=.TRUE.
      !   endif
      !   CFLval = minval(0.5*zroof*zroof*rroof/kroof)
      !   if (Tstep>CFLval) then !CFL condition
      !      write(*,*) "ROOF: CFL condition: Tstep=",Tstep,">",CFLval
      !      CFLfail=.TRUE.
      !   endif
      !   CFLval = minval(0.5*zwall*zwall*rwall/kwall)
      !   if (Tstep>CFLval) then !CFL condition
      !      write(*,*) "WALL: CFL condition: Tstep=",Tstep,">",CFLval
      !      CFLfail=.TRUE.
      !   endif
      !   CFLval = minval(0.5*zground*zground*rground/kground)
      !   if (Tstep>CFLval) then !CFL condition
      !      write(*,*) "ground: CFL condition: Tstep=",Tstep,">",CFLval
      !      CFLfail=.TRUE.
      !   endif
      !   if (CFLfail) then
      !      write(*,*) "Increase dX or decrease maxtimestep. Hit any key to continue"
      !      read (*,*)
      !   endif
 
      ! Tiaircyc = (1+(LondonQSJune_Barbican.Tair-Tiair)./(5*Tiair)).*(Tiair + 0.4*sin(LondonQSJune_Barbican.HOUR*2*pi/24-10/24*2*pi))    !!FO!! outdoor temp affected
 
      RETURN
 
      !     315 CALL errorHint(11,TRIM(fileESTMTs),notUsed,notUsed,NotUsedI)
      ! 316 CALL errorHint(11,TRIM(fileFinalTemp),notUsed,notUsed,NotUsedI)
 

References allocatearray::alb, estm_data::alb_avg, estm_data::alb_ground_estm, estm_data::alb_roof_estm, estm_data::alb_veg_estm, estm_data::alb_wall_fix, gis_data::bldgh, allocatearray::bldgsurf, allocatearray::bsoilsurf, physconstants::c2k, allocatearray::conifsurf, allocatearray::cts_tiair, allocatearray::cts_troad, allocatearray::cts_troof, allocatearray::cts_twall, allocatearray::decidsurf, estm_data::em_f, estm_data::em_ground_estm, estm_data::em_i, estm_data::em_ibld, estm_data::em_r, estm_data::em_roof_estm, estm_data::em_veg_estm, estm_data::em_w, estm_data::em_wall_fix, allocatearray::emis, allocatearray::estmforcingdata, estm_data::fair, estm_data::fground, estm_data::fibld, estm_data::finternal, estm_data::first, estm_data::froof, estm_data::fveg, estm_data::fwall, allocatearray::grasssurf, estm_data::hw, estm_data::ivf_fi, estm_data::ivf_fr, estm_data::ivf_fw, estm_data::ivf_if, estm_data::ivf_ii, estm_data::ivf_ir, estm_data::ivf_iw, estm_data::ivf_rf, estm_data::ivf_ri, estm_data::ivf_rw, estm_data::ivf_wf, estm_data::ivf_wi, estm_data::ivf_wr, estm_data::ivf_ww, estm_data::lbc_soil, estm_data::lup_ground, estm_data::lup_roof, estm_data::lup_wall, estm_data::minshc_airbld, estm_data::nalb, defaultnotused::nan, allocatearray::ncolsestmdata, estm_data::ndepth_ground, estm_data::ndepth_ibld, estm_data::ndepth_roof, estm_data::ndepth_wall, estm_data::nemis, estm_data::nroom, initial::numberofgrids, allocatearray::pavsurf, estm_data::rvf_canyon, estm_data::rvf_ground, estm_data::rvf_roof, estm_data::rvf_veg, estm_data::rvf_wall, physconstants::sbconst, allocatearray::sfr_surf, estm_data::shc_air, estm_data::sumalb, estm_data::sumemis, estm_data::svf_ground, estm_data::svf_roof, estm_data::svf_wall, estm_data::t0_ground, estm_data::t0_ibld, estm_data::t0_roof, estm_data::t0_wall, allocatearray::tair24hr, estm_data::tfloor, estm_data::tground, estm_data::tground_grids, estm_data::tibld, estm_data::tibld_grids, estm_data::tievolve, estm_data::tn_roof, estm_data::tn_wall, estm_data::troof, estm_data::troof_grids, allocatearray::ts5mindata, estm_data::tw_4, estm_data::tw_4_grids, estm_data::twall, estm_data::twall_grids, allocatearray::watersurf, estm_data::xvf_wall, estm_data::zref, estm_data::zvf_ground, and estm_data::zvf_wall.

load_gridlayout()

subroutine estm_module::load_gridlayout	(	integer, intent(in)	gridiv,
		logical, intent(in)	multiplelayoutfiles,
		integer, intent(in)	diagnose )

Definition at line 828 of file suews_phys_estm.f95.

      USE allocatearray, ONLY: &
         ndepth, &
         nlayer, nlayer_grids, &
         height, height_grids, &
         building_frac, building_frac_grids, &
         veg_frac, veg_frac_grids, &
         building_scale, building_scale_grids, &
         veg_scale, veg_scale_grids, &
         roof_albedo_dir_mult_fact, roof_albedo_dir_mult_fact_grids, &
         wall_specular_frac, wall_specular_frac_grids, &
         sfr_roof, sfr_roof_grids, &
         tin_roof, tin_roof_grids, &
         alb_roof, alb_roof_grids, &
         emis_roof, emis_roof_grids, &
         state_roof, state_roof_grids, &
         statelimit_roof, statelimit_roof_grids, &
         wetthresh_roof, wetthresh_roof_grids, &
         soilstore_roof, soilstore_roof_grids, &
         soilstorecap_roof, soilstorecap_roof_grids, &
         k_roof, k_roof_grids, &
         cp_roof, cp_roof_grids, &
         dz_roof, dz_roof_grids, &
         temp_roof_grids, &
         sfr_wall, sfr_wall_grids, &
         tin_wall, tin_wall_grids, &
         alb_wall, alb_wall_grids, &
         emis_wall, emis_wall_grids, &
         state_wall, state_wall_grids, &
         statelimit_wall, statelimit_wall_grids, &
         wetthresh_wall, wetthresh_wall_grids, &
         soilstore_wall, soilstore_wall_grids, &
         soilstorecap_wall, soilstorecap_wall_grids, &
         k_wall, k_wall_grids, &
         cp_wall, cp_wall_grids, &
         dz_wall, dz_wall_grids, &
         temp_wall_grids, &
         nsurf, &
         dz_surf, dz_surf_grids, &
         k_surf, k_surf_grids, &
         cp_surf, cp_surf_grids, &
         tin_surf, tin_surf_grids, &
         nspec
      USE data_in, ONLY: fileinputpath, filecode
      USE strings, ONLY: writenum
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: gridIV
      INTEGER, INTENT(IN) :: diagnose
      LOGICAL, INTENT(IN) :: MultipleLayoutFiles
      CHARACTER(len=100) :: FileLayout, str_gridIV
      INTEGER :: istat, iunit, ERROR_UNIT, i
      INTEGER :: igroup, ilayer, idepth
      CHARACTER(len=1000) :: line
      REAL(KIND(1D0)) :: k, dz
 
      namelist &
         ! basic dimension info
         /dim/ &
         nlayer, &
         /geom/ &
         height, &
         building_frac, &
         veg_frac, &
         building_scale, &
         veg_scale, &
         ! read in roof part
         /roof/ &
         sfr_roof, &
         roof_albedo_dir_mult_fact, &
         tin_roof, &
         alb_roof, &
         emis_roof, &
         state_roof, &
         statelimit_roof, &
         wetthresh_roof, &
         soilstore_roof, &
         soilstorecap_roof, &
         dz_roof, &
         k_roof, &
         cp_roof &
         ! read in wall part
         /wall/ &
         sfr_wall, &
         wall_specular_frac, &
         tin_wall, &
         alb_wall, &
         emis_wall, &
         state_wall, &
         statelimit_wall, &
         wetthresh_wall, &
         soilstore_wall, &
         soilstorecap_wall, &
         dz_wall, &
         k_wall, &
         cp_wall, &
         /surf/ &
         tin_surf, &
         dz_surf, &
         k_surf, &
         cp_surf
 
      IF (multiplelayoutfiles) THEN
         CALL writenum(gridiv, str_gridiv, 'i4')
         filelayout = 'GridLayout'//trim(filecode)//trim(str_gridiv)//'.nml'
      ELSE
         filelayout = 'GridLayout'//trim(filecode)//'.nml'
      END IF
 
      IF (diagnose == 1) print *, 'Reading layout file: ', trim(fileinputpath)//filelayout
      iunit = 212
      OPEN (iunit, file=trim(fileinputpath)//trim(filelayout), status='old')
      IF (diagnose == 1) print *, 'Read dim info of GridLayout'
      READ (iunit, nml=dim, iostat=istat)
      IF (diagnose == 1) print *, 'Number of layers to read: ', nlayer
 
      IF (istat /= 0) THEN
         backspace(iunit)
         READ (iunit, fmt='(A)') line
         error_unit = 119
         WRITE (error_unit, '(A)') &
            'Invalid line in namelist: '//trim(line)
      END IF
      CLOSE (iunit)
 
      ALLOCATE (height(nlayer + 1))
      ALLOCATE (building_frac(nlayer))
      ALLOCATE (veg_frac(nlayer))
      ALLOCATE (building_scale(nlayer))
      ALLOCATE (veg_scale(nlayer))
 
      ALLOCATE (sfr_roof(nlayer))
      ALLOCATE (dz_roof(nlayer, ndepth))
      ALLOCATE (k_roof(nlayer, ndepth))
      ALLOCATE (cp_roof(nlayer, ndepth))
      ALLOCATE (tin_roof(nlayer))
      ALLOCATE (alb_roof(nlayer))
      ALLOCATE (emis_roof(nlayer))
      ALLOCATE (state_roof(nlayer))
      ALLOCATE (statelimit_roof(nlayer))
      ALLOCATE (wetthresh_roof(nlayer))
      ALLOCATE (soilstore_roof(nlayer))
      ALLOCATE (soilstorecap_roof(nlayer))
      ALLOCATE (roof_albedo_dir_mult_fact(nspec, nlayer))
 
      ALLOCATE (sfr_wall(nlayer))
      ALLOCATE (dz_wall(nlayer, ndepth))
      ALLOCATE (k_wall(nlayer, ndepth))
      ALLOCATE (cp_wall(nlayer, ndepth))
      ALLOCATE (tin_wall(nlayer))
      ALLOCATE (alb_wall(nlayer))
      ALLOCATE (emis_wall(nlayer))
      ALLOCATE (state_wall(nlayer))
      ALLOCATE (statelimit_wall(nlayer))
      ALLOCATE (wetthresh_wall(nlayer))
      ALLOCATE (soilstore_wall(nlayer))
      ALLOCATE (soilstorecap_wall(nlayer))
      ALLOCATE (wall_specular_frac(nspec, nlayer))
 
      ALLOCATE (dz_surf(nsurf, ndepth))
      ALLOCATE (k_surf(nsurf, ndepth))
      ALLOCATE (cp_surf(nsurf, ndepth))
      ALLOCATE (tin_surf(nsurf))
 
      OPEN (iunit, file=trim(fileinputpath)//trim(filelayout), status='old')
      IF (diagnose == 1) print *, 'Read geometry part of GridLayout'
      READ (iunit, nml=geom, iostat=istat)
      IF (diagnose == 1) print *, 'height', height
      IF (diagnose == 1) print *, 'building_frac', building_frac
      IF (diagnose == 1) print *, 'veg_frac', veg_frac
      IF (diagnose == 1) print *, 'building_scale', building_scale
      IF (diagnose == 1) print *, 'veg_scale', veg_scale
 
      IF (diagnose == 1) print *, 'Read roof part of GridLayout'
      READ (iunit, nml=roof, iostat=istat)
      IF (diagnose == 1) print *, 'sfr_roof', sfr_roof
      IF (diagnose == 1) print *, 'dz_roof', dz_roof
      IF (diagnose == 1) print *, 'k_roof', k_roof
      IF (diagnose == 1) print *, 'cp_roof', cp_roof
      IF (diagnose == 1) print *, 'tin_roof', tin_roof
      IF (diagnose == 1) print *, 'alb_roof', alb_roof
      IF (diagnose == 1) print *, 'emis_roof', emis_roof
      IF (diagnose == 1) print *, 'state_roof', state_roof
      IF (diagnose == 1) print *, 'statelimit_roof', statelimit_roof
      IF (diagnose == 1) print *, 'wetthresh_roof', wetthresh_roof
      IF (diagnose == 1) print *, 'soilstore_roof', soilstore_roof
      IF (diagnose == 1) print *, 'soilstorecap_roof', soilstorecap_roof
 
      IF (diagnose == 1) print *, 'Read wall part of GridLayout'
      READ (iunit, nml=wall, iostat=istat)
      IF (diagnose == 1) print *, 'sfr_wall', sfr_wall
      IF (diagnose == 1) print *, 'dz_wall', dz_wall
      IF (diagnose == 1) print *, 'k_wall', k_wall
      IF (diagnose == 1) print *, 'cp_wall', cp_wall
      IF (diagnose == 1) print *, 'tin_wall', tin_wall
      IF (diagnose == 1) print *, 'alb_wall', alb_wall
      IF (diagnose == 1) print *, 'emis_wall', emis_wall
      IF (diagnose == 1) print *, 'state_wall', state_wall
      IF (diagnose == 1) print *, 'statelimit_wall', statelimit_wall
      IF (diagnose == 1) print *, 'wetthresh_wall', wetthresh_wall
      IF (diagnose == 1) print *, 'soilstore_wall', soilstore_wall
      IF (diagnose == 1) print *, 'soilstorecap_wall', soilstorecap_wall
 
      ! PRINT *, 'Read surf part of GridLayout'
      READ (iunit, nml=surf, iostat=istat)
      IF (istat /= 0) THEN
         backspace(iunit)
         READ (iunit, fmt='(A)') line
         error_unit = 119
         WRITE (error_unit, '(A)') &
            'Invalid line in namelist: '//trim(line)
      END IF
      CLOSE (iunit)
 
      IF (diagnose == 1) print *, 'Read GridLayout'
      IF (diagnose == 1) print *, 'sfr_roof in load_GridLayout', sfr_roof
      IF (diagnose == 1) print *, 'sfr_wall in load_GridLayout', sfr_wall
      IF (diagnose == 1) print *, 'tin_surf in load_GridLayout', tin_surf
      IF (diagnose == 1) print *, 'dz_surf(1,:) in load_GridLayout', dz_surf(1, :)
      IF (diagnose == 1) print *, 'k_surf(1,:) in load_GridLayout', k_surf(1, :)
      IF (diagnose == 1) print *, 'dz_roof in load_GridLayout', dz_roof(1:nlayer, :ndepth)
      ! PRINT *, 'dz_roof in load_GridLayout', dz_roof(1, :ndepth)
 
      ! assign values to the grid container
      nlayer_grids(gridiv) = nlayer
      height_grids(gridiv, 1:nlayer + 1) = height(1:nlayer + 1)
      building_frac_grids(gridiv, 1:nlayer) = building_frac(1:nlayer)
      veg_frac_grids(gridiv, 1:nlayer) = veg_frac(1:nlayer)
      building_scale_grids(gridiv, 1:nlayer) = building_scale(1:nlayer)
      veg_scale_grids(gridiv, 1:nlayer) = veg_scale(1:nlayer)
 
      ! roof
      sfr_roof_grids(gridiv, 1:nlayer) = sfr_roof
      tin_roof_grids(gridiv, 1:nlayer) = tin_roof
      DO i = 1, nlayer
         temp_roof_grids(gridiv, i, :) = tin_roof(i)
      END DO
      dz_roof_grids(gridiv, 1:nlayer, 1:ndepth) = dz_roof(1:nlayer, 1:ndepth)
      k_roof_grids(gridiv, 1:nlayer, 1:ndepth) = k_roof(1:nlayer, 1:ndepth)
      cp_roof_grids(gridiv, 1:nlayer, 1:ndepth) = cp_roof(1:nlayer, 1:ndepth)
      alb_roof_grids(gridiv, 1:nlayer) = alb_roof(1:nlayer)
      emis_roof_grids(gridiv, 1:nlayer) = emis_roof(1:nlayer)
      state_roof_grids(gridiv, 1:nlayer) = state_roof(1:nlayer)
      statelimit_roof_grids(gridiv, 1:nlayer) = statelimit_roof(1:nlayer)
      wetthresh_roof_grids(gridiv, 1:nlayer) = wetthresh_roof(1:nlayer)
      soilstore_roof_grids(gridiv, 1:nlayer) = soilstore_roof(1:nlayer)
      soilstorecap_roof_grids(gridiv, 1:nlayer) = soilstorecap_roof(1:nlayer)
      roof_albedo_dir_mult_fact_grids(gridiv, 1:nspec, 1:nlayer) = roof_albedo_dir_mult_fact
      ! PRINT *, 'dz_roof_grids in loading', dz_roof_grids(Gridiv, 1:nroof, 1:ndepth)
      ! PRINT *, 'dz_roof_grids(1,1) in loading', dz_roof_grids(Gridiv, 1, 1:ndepth)
      ! tsfc_roof_grids(gridIV) = tsfc_roof
      ! tin_roof_grids(gridIV) = temp_in_roof
 
      ! wall
      sfr_wall_grids(gridiv, 1:nlayer) = sfr_wall
      tin_wall_grids(gridiv, 1:nlayer) = tin_wall
      DO i = 1, nlayer
         temp_wall_grids(gridiv, i, :) = tin_wall(i)
      END DO
      dz_wall_grids(gridiv, 1:nlayer, 1:ndepth) = dz_wall(1:nlayer, 1:ndepth)
      k_wall_grids(gridiv, 1:nlayer, 1:ndepth) = k_wall(1:nlayer, 1:ndepth)
      cp_wall_grids(gridiv, 1:nlayer, 1:ndepth) = cp_wall(1:nlayer, 1:ndepth)
      alb_wall_grids(gridiv, 1:nlayer) = alb_wall(1:nlayer)
      emis_wall_grids(gridiv, 1:nlayer) = emis_wall(1:nlayer)
      state_wall_grids(gridiv, 1:nlayer) = state_wall(1:nlayer)
      statelimit_wall_grids(gridiv, 1:nlayer) = statelimit_wall(1:nlayer)
      wetthresh_wall_grids(gridiv, 1:nlayer) = wetthresh_wall(1:nlayer)
      soilstore_wall_grids(gridiv, 1:nlayer) = soilstore_wall(1:nlayer)
      soilstorecap_wall_grids(gridiv, 1:nlayer) = soilstorecap_wall(1:nlayer)
      wall_specular_frac_grids(gridiv, 1:nspec, 1:nlayer) = wall_specular_frac
      ! PRINT *, 'dz_wall_grids in loading', dz_wall_grids(Gridiv, 1:nwall, 1:ndepth)
      ! PRINT *, 'dz_wall_grids(1,1) in loading', dz_wall_grids(Gridiv, 1, 1:ndepth)
      ! tsfc_wall_grids(gridIV) = tsfc_wall
      ! tin_wall_grids(gridIV) = temp_in_wall
 
      ! generic surfaces
      tin_surf_grids(gridiv, 1:nsurf) = tin_surf
      ! DO i = 1, nsurf
      !    tin_surf_grids(gridIV, i, :) = tin_surf(i)
      ! END DO
      dz_surf_grids(gridiv, 1:nsurf, 1:ndepth) = dz_surf(1:nsurf, 1:ndepth)
      k_surf_grids(gridiv, 1:nsurf, 1:ndepth) = k_surf(1:nsurf, 1:ndepth)
      cp_surf_grids(gridiv, 1:nsurf, 1:ndepth) = cp_surf(1:nsurf, 1:ndepth)
 
      ! ! check if CFL condition is met
      ! DO igroup = 1, 3
      !    DO ilayer = 1, merge(nlayer,7,igroup<3)
      !       idepth = 1
      !       ! DO idepth = 1, ndepth
      !       IF (igroup == 1) THEN
      !          k = k_roof_grids(gridIV, ilayer, idepth)
      !          dz = dz_roof_grids(gridIV, ilayer, idepth)
      !          str_igroup='roof'
      !       ELSE IF (igroup == 2) THEN
      !          k = k_wall_grids(gridIV, ilayer, idepth)
      !          dz = dz_wall_grids(gridIV, ilayer, idepth)
      !          str_igroup='wall'
      !       ELSE IF (igroup == 3) THEN
      !          k = k_surf_grids(gridIV, ilayer, idepth)
      !          dz = dz_surf_grids(gridIV, ilayer, idepth)
      !          str_igroup='surf'
      !       END IF
 
      !       IF (dz/k > 0.02) THEN
      !          PRINT *, 'CFL condition, dz/k < 0.02,  is not met for top layer of:'
      !           PRINT *, 'grid', gridIV
 
      !           PRINT *, 'group: ', str_igroup
      !           PRINT *, 'layer', ilayer
      !            PRINT *, 'dz/k=', dz/k
      !            PRINT *, 'dz=', dz
      !            PRINT *, 'k=', k
      !            PRINT *,  ''
      !          WRITE (str_gridIV, '(I5.5)') gridIV
      !          ! WRITE (str_igroup, '(I5.5)') igroup
      !          WRITE (str_idepth, '(I5.5)') idepth
      !          WRITE (str_ilayer, '(I5.5)') ilayer
 
      !          CALL ErrorHint(77, 'CFL condition is not met for '// &
      !                         'grid: '//TRIM(str_gridIV)// &
      !                         ', group: '//TRIM(str_igroup)// &
      !                         ', layer: '//TRIM(str_ilayer)// &
      !                         ', depth: '//TRIM(str_idepth) &
      !                         , dz/k, idepth*1.D0, ilayer)
      !       END IF
      !       ! END DO
      !    END DO
 
      ! END DO
 
      DEALLOCATE (height)
      DEALLOCATE (building_frac)
      DEALLOCATE (veg_frac)
      DEALLOCATE (building_scale)
      DEALLOCATE (veg_scale)
 
      DEALLOCATE (sfr_roof)
      DEALLOCATE (alb_roof)
      DEALLOCATE (emis_roof)
      DEALLOCATE (state_roof)
      DEALLOCATE (statelimit_roof)
      DEALLOCATE (wetthresh_roof)
      DEALLOCATE (soilstore_roof)
      DEALLOCATE (soilstorecap_roof)
      DEALLOCATE (dz_roof)
      DEALLOCATE (k_roof)
      DEALLOCATE (cp_roof)
      DEALLOCATE (tin_roof)
      DEALLOCATE (roof_albedo_dir_mult_fact)
 
      DEALLOCATE (sfr_wall)
      DEALLOCATE (alb_wall)
      DEALLOCATE (emis_wall)
      DEALLOCATE (state_wall)
      DEALLOCATE (statelimit_wall)
      DEALLOCATE (wetthresh_wall)
      DEALLOCATE (soilstore_wall)
      DEALLOCATE (soilstorecap_wall)
      DEALLOCATE (dz_wall)
      DEALLOCATE (k_wall)
      DEALLOCATE (cp_wall)
      DEALLOCATE (tin_wall)
      DEALLOCATE (wall_specular_frac)
 
      DEALLOCATE (tin_surf)
      DEALLOCATE (dz_surf)
      DEALLOCATE (k_surf)
      DEALLOCATE (cp_surf)
 

References allocatearray::alb_roof, allocatearray::alb_roof_grids, allocatearray::alb_wall, allocatearray::alb_wall_grids, allocatearray::building_frac, allocatearray::building_frac_grids, allocatearray::building_scale, allocatearray::building_scale_grids, allocatearray::cp_roof, allocatearray::cp_roof_grids, allocatearray::cp_surf, allocatearray::cp_surf_grids, allocatearray::cp_wall, allocatearray::cp_wall_grids, allocatearray::dz_roof, allocatearray::dz_roof_grids, allocatearray::dz_surf, allocatearray::dz_surf_grids, allocatearray::dz_wall, allocatearray::dz_wall_grids, allocatearray::emis_roof, allocatearray::emis_roof_grids, allocatearray::emis_wall, allocatearray::emis_wall_grids, data_in::filecode, data_in::fileinputpath, allocatearray::height, allocatearray::height_grids, allocatearray::k_roof, allocatearray::k_roof_grids, allocatearray::k_surf, allocatearray::k_surf_grids, allocatearray::k_wall, allocatearray::k_wall_grids, allocatearray::ndepth, allocatearray::nlayer, allocatearray::nlayer_grids, allocatearray::nspec, allocatearray::nsurf, allocatearray::roof_albedo_dir_mult_fact, allocatearray::roof_albedo_dir_mult_fact_grids, allocatearray::sfr_roof, allocatearray::sfr_roof_grids, allocatearray::sfr_wall, allocatearray::sfr_wall_grids, allocatearray::soilstore_roof, allocatearray::soilstore_roof_grids, allocatearray::soilstore_wall, allocatearray::soilstore_wall_grids, allocatearray::soilstorecap_roof, allocatearray::soilstorecap_roof_grids, allocatearray::soilstorecap_wall, allocatearray::soilstorecap_wall_grids, allocatearray::state_roof, allocatearray::state_roof_grids, allocatearray::state_wall, allocatearray::state_wall_grids, allocatearray::statelimit_roof, allocatearray::statelimit_roof_grids, allocatearray::statelimit_wall, allocatearray::statelimit_wall_grids, allocatearray::temp_roof_grids, allocatearray::temp_wall_grids, allocatearray::tin_roof, allocatearray::tin_roof_grids, allocatearray::tin_surf, allocatearray::tin_surf_grids, allocatearray::tin_wall, allocatearray::tin_wall_grids, allocatearray::veg_frac, allocatearray::veg_frac_grids, allocatearray::veg_scale, allocatearray::veg_scale_grids, allocatearray::wall_specular_frac, allocatearray::wall_specular_frac_grids, allocatearray::wetthresh_roof, allocatearray::wetthresh_roof_grids, allocatearray::wetthresh_wall, and allocatearray::wetthresh_wall_grids.

Referenced by estm_ehc_initialise().

Here is the caller graph for this function:

set_nan()

elemental real(kind(1d0)) function estm_module::set_nan

(

real(kind(1d0)), intent(in)

x

)

Definition at line 2663 of file suews_phys_estm.f95.

      IMPLICIT NONE
      REAL(KIND(1D0)), PARAMETER :: pNAN = 9999
      REAL(KIND(1D0)), PARAMETER :: NAN = -999
      REAL(KIND(1D0)), INTENT(in) :: x
      REAL(KIND(1D0)) :: xx
 
      IF (abs(x) > pnan) THEN
         xx = nan
      ELSE
         xx = x
      END IF
 

Referenced by estm().

Here is the caller graph for this function:

suews_getestmdata()

subroutine estm_module::suews_getestmdata

(

integer, intent(in)

lunit

)

Definition at line 719 of file suews_phys_estm.f95.

      USE allocatearray, ONLY: ncolsestmdata, estmforcingdata
      USE data_in, ONLY: fileestmts, skipheadermet
      USE sues_data, ONLY: tstep_real, tstep
      USE defaultnotused, ONLY: notused, ios_out
      USE initial, ONLY: skippedlines, readlinesmetdata, gridcounter
 
      IMPLICIT NONE
 
      INTEGER, INTENT(in) :: lunit
      INTEGER :: i, iyy !,RunNumber,NSHcounter
      INTEGER :: iostat_var
      REAL(KIND(1D0)), DIMENSION(ncolsESTMdata) :: ESTMArray
      REAL(KIND(1D0)) :: imin_prev, ih_prev, iday_prev, tstep_estm !For checks on temporal resolution of estm data
 
      ! initialise
      imin_prev = 0
      ih_prev = 0
      iday_prev = 0
 
      !---------------------------------------------------------------
 
      !Open the file for reading and read the actual data
      !write(*,*) FileESTMTs
      OPEN (lunit, file=trim(fileestmts), status='old', err=315)
      CALL skipheader(lunit, skipheadermet)
 
      ! Skip to the right place in the ESTM file, depending on how many chunks have been read already
      IF (skippedlines > 0) THEN
         DO iyy = 1, skippedlines
            READ (lunit, *)
         END DO
      END IF
 
      ! Read in next chunk of ESTM data and fill ESTMForcingData array with data for every timestep
      DO i = 1, readlinesmetdata
         READ (lunit, *, iostat=iostat_var) estmarray
         estmforcingdata(i, 1:ncolsestmdata, gridcounter) = estmarray
         ! Check timestamp of met data file matches TSTEP specified in RunControl
         IF (i == 1) THEN
            imin_prev = estmarray(4)
            ih_prev = estmarray(3)
            iday_prev = estmarray(2)
         ELSEIF (i == 2) THEN
            tstep_estm = ((estmarray(4) + 60*estmarray(3)) - (imin_prev + 60*ih_prev))*60 !tstep in seconds
            IF (tstep_estm /= tstep_real .AND. estmarray(2) == iday_prev) THEN
               CALL errorhint(39, 'TSTEP in RunControl does not match TSTEP of ESTM data (DOY).', &
                              REAL(tstep, KIND(1D0)), tstep_estm, INT(ESTMArray(2)))
            END IF
         END IF
      END DO
 
      CLOSE (lunit)
 
      RETURN
 
315   CALL errorhint(11, trim(fileestmts), notused, notused, ios_out)
 

References errorhint(), allocatearray::estmforcingdata, data_in::fileestmts, initial::gridcounter, defaultnotused::ios_out, allocatearray::ncolsestmdata, defaultnotused::notused, initial::readlinesmetdata, skipheader(), data_in::skipheadermet, initial::skippedlines, sues_data::tstep, and sues_data::tstep_real.

Here is the call graph for this function: