SUEWS API Site
Documentation of SUEWS source code
Functions/Subroutines
estm_module Module Reference

Functions/Subroutines

subroutine suews_getestmdata (lunit)
 
subroutine estm_initials
 
subroutine load_gridlayout (gridIV, MultipleLayoutFiles, diagnose)
 
subroutine estm_ext_initialise
 
subroutine estm_ext_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)
 
subroutine estm_ext (tstep, nlayer, QG_surf, qg_roof, qg_wall, tsfc_roof, tin_roof, temp_in_roof, k_roof, cp_roof, dz_roof, sfr_roof, tsfc_wall, tin_wall, temp_in_wall, k_wall, cp_wall, dz_wall, sfr_wall, tsfc_surf, tin_surf, temp_in_surf, k_surf, cp_surf, dz_surf, sfr_surf, temp_out_roof, QS_roof, temp_out_wall, QS_wall, temp_out_surf, QS_surf, 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 1582 of file suews_phys_estm.f95.

1589 ! NB: HCW Questions:
1590 ! - should TFloor be set in namelist instead of hard-coded here?
1591 ! - zref used for radiation calculation and fair is set to 2*BldgH here. For compatibility with the rest of the
1592 ! SUEWS model, should this be the (wind speed) measurement height z specified in RunControl.nml?
1593 ! - In SUEWS_translate, fwall=AreaWall/SurfaceArea. Is this correct?
1594 ! - If froof=1 (i.e. whole grid is building), is HW=0 correct?
1595 ! - Then is an IF(Fground ==0) option needed?
1596 ! - alb_wall=0.23 and em_wall=0.9 are set in LUMPS_module_constants. Shouldn't these be provided as input?
1597 ! - Do the LUP calculations here need to be compatible with those in LUMPS_NARP?
1598 ! - File opening rewritten using existing error handling in SUEWS - can delete mod_error module from SUEWS_ESTM_functions
1599 ! - In SUEWS_ESTM_v2016, the first row is set to -999. This may be acceptable at the
1600 ! start of the run but should be handled properly between blocks of met data? - need to check what's actually happening here.
1601 ! - Many duplicate functions in SUEWS_ESTM_functions need changing to the existing SUEWS versions.
1602 ! - Are the following correctly initialised? T0_ibld, T0_ground, T0_wall, T0_roof, TN_wall, TN_roof, Tground, Twall, Troof, Tibld
1603 ! - What are Nalb, sumalb, Nemis and Sumemis for? Are they used correctly?
1604 !
1605
1606 !SUEWS_ESTM_v2016
1607 ! revision:
1608 ! HCW 14 Jun 2016
1609 ! HCW 27 Jun 2016 Corrected iESTMcount bug - now increases for all grids together
1610 ! HCW 30 Jun 2016 Major changes to handle grids and met blocks
1611 ! TS 09 Oct 2017 Added explicit interface
1612
1613 !Contains calculation for each time step
1614 !Calculate local scale heat storage from single building and surroundings observations
1615 !OFferle, May 2003
1616 !
1617 !MODIFICATION HISTORY
1618 ! 15 DECEMBER 2003
1619 ! (1) CHANGED AIR EXCHANGE RATE TO ALSO BE DEPENDENT ON OUTSIDE AIR TEMPERATURE
1620 ! (2) ADDED SPECIFIC HEAT CALCULATION FOR AIR
1621 ! (3) RH ADDED AS VAR(14) IN INPUT FILE
1622 ! 12 JANUARY 2004
1623 ! (1) ADDED ESTIMATED AVG NADIR LOOKING EXTERNAL SURFACE TEMPERATURE TO OUTPUT
1624 ! WEIGHTED BY SURFACE FRACTION AND SVF. SVF_ROOF=1, SVF_WALLS = 1-SVF_ground
1625 ! (2) ADDED NET RADIATION (RN) CALCULATIONS FOR ALL SURFACES AND AVG RN TO OUTPUT BASED ON
1626 ! RADIOMETER VIEW FROM ZREF
1627 ! 14 JANUARY 2004
1628 ! (1) MOVED GRID SPECIFIC PARAMETERS OUT OF NAMELIST, INTO PARAMETER FILE E.G. ALB,EM,F
1629 ! T0 MAKE CHANGES IN LOCATIONS EASIER.
1630 !
1631 ! 23 JANUARY 2004 - Version 2
1632 ! (1) PUT ALL TEMPERATURES INTO K
1633 ! (2) added interpolation routine to run on shorter timesteps.
1634 ! tested so that it doesn't change solution for forced temperatures.
1635 ! however in version 2 the energy balance at the surface isn't correctly solved
1636 !
1637 ! 25 JANUARY 2004 - Version 3
1638 ! (1) added solution to energy balance at surface
1639 ! (2) removed some extraneous code
1640 ! (3) added vegetation fractions for future development
1641 ! (4) some changes to input namelist.
1642 ! (5) need to add wind speed dependence for exchange coefficients
1643 ! (6) changed the way Rn_net is calculated. also radiometer view factor relationships
1644 ! (7) added a calculation for heat loss/gain to outside air (that going into building mass is storage)
1645 ! this is labelled QFBLD which it is in a sense.
1646 ! 6 FEBRUARY 2004
1647 ! (1) ADDED INTERNAL VIEW FACTOR FILE FOR INTERNAL GEOMETRY, INCLUDING FIXED FLOOR TEMPERATURE
1648 ! (2) added MeanU, MinWS to config, and U to ILOC.
1649 !
1650 ! 11 FEBRUARY 2004
1651 ! (1) added site lat, long, elevation to inputs
1652 ! (2) need zenith angle for wall direct radiation interception
1653 !
1654 ! 15 JUNE 2004
1655 ! (1) CORRECTED OUTPUT OF T0 FOR FORCED SURFACE TEMPERATURES
1656 ! (2) MADE SOME CHANGES TO RADIATION, COMPUTATION OF AVERAGE ALBEDO, ADDED AVERAGE EMISSIVITY
1657 ! (3) ADDED OUTPUT FILE FOR RADIATION COMPONENTS
1658 ! (4) CHANGED INPUT CONFIG SO CONVECTIVE EXCHANGE COEFFS ARE NOT USER SELECTABLE
1659 ! (5) MAY STILL BE PROBLEMS WITH WALL RADIATIVE EXCHANGE AND AVERAGE ALBEDO
1660 ! (6) CHANGED THE WAY HEAT STORAGE IS COMPUTED IN HEATCONDUCTION MODULE BUT THIS SHOULD NOT CHANGE RESULTS
1661 !
1662 ! 16 NOVEMBER 2004
1663 ! (1) CHANGED AIR EXCHANGE RATE TO BE BASES ON DAILY TEMPERATURE CHANGES.
1664 ! (2) WRITES OUT FINAL LAYER TEMPERATURES AND INCLUDES OPTION TO READ AT BEGINNING.
1665 !
1666 !DESCRIPTION: uses explicit time differencing to compute heat conduction through roof, walls,
1667 ! internal mass, and grounds (elements). Air heat storage is computed from average air temperature
1668 ! Boundary conditions are determined by measured surface temperature(s) or computed
1669 ! from energy balance at the surface.
1670 ! Internal air temperature can either be fixed or allowed to evolve in response
1671 ! to air mass exchanges and convective heating from internal surfaces.
1672 !
1673 !INPUT:
1674 !FORCING DATA: DTIME,KDOWN,LDOWN,TSURF,TAIR_OUT,TAIR_IN,TROOF,TWALL_AVG,TWALL_N,TWALL_E,TWALL_S,TWALL_W,Tground,RH,U
1675 !NAMELIST : HEATSTORAGE_vFO.NML
1676 ! &config
1677 ! ifile=FORCING DATA
1678 ! ofile=OUTPUT FILE
1679 ! pfile=HEAT STORAGE PARAMETER FILE
1680 ! Nibld = INTERNAL MASS LAYERS
1681 ! Nwall = EXTERNAL WALL LAYERS
1682 ! Nroof = ROOF LAYERS
1683 ! Nground = ground/SOIL LAYERS
1684 ! LBC_soil = LOWER BOUNDARY CONDITION FOR ground/SOIL
1685 ! iloc= INPUT COLUMNS IN DATA FILE
1686 ! evolveTibld= USE DIAGNOSTIC VALUE FOR INTERNAL BUILDING TEMPERATURE
1687 ! 0: don't use, use measured
1688 ! 1: TURN ON USE when temp goess ABOVE TINT_ON, off when temp is below TINT_OFF
1689 ! 2: always use diagnostic
1690 ! THEAT_ON= TEMPERATURE AT WHICH HEAT CONTROL IS TURNED ON
1691 ! THEAT_OFF= TEMPERATURE AT WHICH HEAT CONTROL IS TURNED OFF
1692 ! THEAT_FIX = Fixed internal temperature for climate control
1693 ! oneTsurf= USE SINGLE SURFACE TEMPERATURE TO DRIVE ALL LAYERS
1694 ! radforce= USE RADIATIVE ENERGY BALANCE TO DRIVE EXTERNAL TEMPERATURES
1695 ! maxtimestep=302, maximum time step in s for filling data gaps.
1696 ! Alt = STATION HEIGHT (m) FOR PRESSURE CALCULATION
1697 ! SPINUP = NUMBER OF LINES TO USE FOR SPINUP (REPEATS THESE LINES BUT ONLY OUTPUTS THE 2ND TIME)
1698 ! INITTEMP = if TRUE INITIALIZES TEMPERATURES TO THOSE IN FINALTEMP.TXT FILE
1699 ! CH_ibld = INTERNAL BUILDING CONVECTIVE EXCHANGE COEFFICIENT
1700 ! **** THESE SHOULD DEPEND ON WIND SPEED BUT CURRENTLY DO NOT ****
1701 ! CHAIR = CONVECTIVE EXCHANGE COEFFICIENT FOR ROOF
1702 ! chair_ground = ... FOR ground
1703 ! chair_wall = ... FOR WALL
1704 ! /
1705 ! ***************** PARAMETER FILE VARIABLES
1706 ! fveg = FRACTION OF ground SURFACE COVERED WITH VEG
1707 ! zveg = VEGETATION HEIGHT
1708 ! alb_veg = VEGETATION ALBEDO
1709 ! em_veg = VEGETATION EMISSIVITY
1710 ! ZREF = REFERENCE HEIGHT FOR FLUX CALCULATION
1711 ! BldgH = mean building height
1712 ! HW = CANYON ASPECT RATION
1713 ! f_X = FRACTION OF X WHERE X IS INTERNAL, WALL, ROOF, ground
1714 ! Alb_x = ALBEDO OF X
1715 ! em_ibld = EMISSIVITY OF X
1716 ! TX = INITIAL LAYER TEMPERATURES
1717 ! zX = LAYER THICKNESS
1718 ! kX = LAYER THERMAL CONDUCTIVITY
1719 ! ribld = LAYER VOLUMETRIC HEAT CAPACITY
1720 !
1721 !****************** INTERNAL VIEW FACTOR FILE
1722 !OUTPUT: fixed format text with single header, heatstorage for all elements, and temperatures
1723 ! for each element-layer.
1724 !===============================================================================
1725
1726 USE meteo, ONLY: pi, heatcapacity_air
1727 USE mod_solver, ONLY: newtonpolynomial
1728 ! USE modSolarCalc !!FO!! :modsolarcalc.f95
1729 ! USE MathConstants !!FO!! :MathConstants_module.f95
1730 USE physconstants, ONLY: c2k, sbconst
1731 USE heatflux, ONLY: heatcond1d
1732 USE estm_data
1733
1734 IMPLICIT NONE
1735 INTEGER, PARAMETER :: ncolsESTMdata = 13
1736 ! INTEGER, PARAMETER:: ncolumnsDataOutESTM=32
1737 INTEGER, PARAMETER :: cTs_Tiair = 5
1738 INTEGER, PARAMETER :: cTs_Tsurf = 6
1739 INTEGER, PARAMETER :: cTs_Troof = 7
1740 INTEGER, PARAMETER :: cTs_Troad = 8
1741 INTEGER, PARAMETER :: cTs_Twall = 9
1742 INTEGER, PARAMETER :: cTs_Twall_n = 10
1743 INTEGER, PARAMETER :: cTs_Twall_e = 11
1744 INTEGER, PARAMETER :: cTs_Twall_s = 12
1745 INTEGER, PARAMETER :: cTs_Twall_w = 13
1746 REAL(KIND(1D0)), PARAMETER :: NAN = -999
1747
1748 INTEGER, INTENT(in) :: Gridiv
1749 INTEGER, INTENT(in) :: tstep
1750 ! INTEGER,INTENT(in)::iy !Year
1751 ! INTEGER,INTENT(in)::id !Day of year
1752 ! INTEGER,INTENT(in)::it !Hour
1753 ! INTEGER,INTENT(in)::imin !Minutes
1754
1755 REAL(KIND(1D0)), INTENT(in) :: avkdn
1756 REAL(KIND(1D0)), INTENT(in) :: avu1
1757 REAL(KIND(1D0)), INTENT(in) :: temp_c
1758 REAL(KIND(1D0)), INTENT(in) :: zenith_deg
1759 REAL(KIND(1D0)), INTENT(in) :: avrh
1760 REAL(KIND(1D0)), INTENT(in) :: press_hpa
1761 REAL(KIND(1D0)), INTENT(in) :: ldown
1762 REAL(KIND(1D0)), INTENT(in) :: bldgh
1763 ! REAL(KIND(1d0)),INTENT(in):: dectime !Decimal time
1764 REAL(KIND(1D0)), DIMENSION(ncolsESTMdata), INTENT(in) :: Ts5mindata_ir !surface temperature input data
1765 REAL(KIND(1D0)), INTENT(in) :: Tair_av ! mean air temperature of past 24hr
1766
1767 REAL(KIND(1D0)), DIMENSION(27), INTENT(out) :: dataOutLineESTM
1768 !Output to SUEWS
1769 REAL(KIND(1D0)), INTENT(out) :: QS
1770 !Input from SUEWS, corrected as Gridiv by TS 09 Jun 2016
1771
1772 !Use only in this subroutine
1773 INTEGER :: i, ii
1774 INTEGER :: Tair2Set = 0
1775 REAL(KIND(1D0)) :: AIREXHR, AIREXDT
1776 REAL(KIND(1D0)), DIMENSION(2) :: bc
1777 REAL(KIND(1D0)) :: chair_ground, chair_wall
1778 REAL(KIND(1D0)) :: EM_EQUIV
1779 REAL(KIND(1D0)) :: kdz
1780 REAL(KIND(1D0)) :: kup_estm, LUP_net, kdn_estm
1781 REAL(KIND(1D0)) :: QHestm
1782 REAL(KIND(1D0)) :: QFBld !Anthropogenic heat from HVAC
1783 REAL(KIND(1D0)) :: shc_airbld
1784 REAL(KIND(1D0)) :: sw_hor, sw_vert
1785 REAL(KIND(1D0)) :: T0
1786 REAL(KIND(1D0)) :: Tinternal, Tsurf_all, Troof_in, Troad, Twall_all, Tw_n, Tw_e, Tw_s, Tw_w
1787 REAL(KIND(1D0)) :: Twallout(5), Troofout(5), Tibldout(5), Tgroundout(5)
1788 REAL(KIND(1D0)) :: Tadd, Tveg
1789 REAL(KIND(1D0)) :: Tairmix
1790 REAL(KIND(1D0)) :: RN
1791 REAL(KIND(1D0)) :: Rs_roof, Rl_roof, RN_ROOF
1792 REAL(KIND(1D0)) :: Rs_wall, Rl_wall, RN_WALL
1793 REAL(KIND(1D0)) :: Rs_ground, Rl_ground, RN_ground
1794 REAL(KIND(1D0)) :: Rs_ibld, Rl_ibld
1795 REAL(KIND(1D0)) :: Rs_iroof, Rl_iroof
1796 REAL(KIND(1D0)) :: Rs_iwall, Rl_iwall
1797 REAL(KIND(1D0)) :: zenith_rad
1798 REAL(KIND(1D0)) :: dum(50)
1799 REAL(KIND(1D0)) :: bldgHX ! local bldgHX=max(bldgH,0.001)
1800 REAL(KIND(1D0)), PARAMETER :: WSmin = 0.1 ! Check why there is this condition. S.O.
1801 LOGICAL :: radforce, groundradforce
1802
1803 radforce = .false.
1804 groundradforce = .false. !Close the radiation scheme in original ESTM S.O.O.
1805
1806 bldghx = max(bldgh, 0.001) ! this is to prevent zero building height
1807
1808 ! Set -999s for first row
1809 ! dum=(/-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,&
1810 ! -999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,&
1811 ! -999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,&
1812 ! -999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,&
1813 ! -999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999./)
1814 dum = [(-999, i=1, 50)]
1815
1816 !External bulk exchange coefficients - set these somewhere more sensible***
1817 chr = 0.005
1818 chair = chr
1819 chair_ground = chair
1820 chair_wall = chair
1821
1822 !Get met data for use in ESTM subroutine
1823 kdn_estm = avkdn
1824 ws = avu1
1825 IF (ws < wsmin) ws = wsmin
1826 tair1 = temp_c + c2k
1827 ! Set initial value of Tair2 to air temp
1828 IF (gridiv == 1) tair2set = tair2set + 1
1829 IF (tair2set == 1) THEN
1830 tair2 = temp_c + c2k
1831 ELSE
1832 tair2 = tair2_grids(gridiv)
1833 ! Also get other variables for this grid
1834 tievolve = tievolve_grids(gridiv)
1836 lup_wall = lup_wall_grids(gridiv)
1837 lup_roof = lup_roof_grids(gridiv)
1838 t0_ibld = t0_ibld_grids(gridiv)
1839 t0_ground = t0_ground_grids(gridiv)
1840 t0_wall = t0_wall_grids(gridiv)
1841 t0_roof = t0_roof_grids(gridiv)
1842 tn_wall = tn_wall_grids(gridiv)
1843 tn_roof = tn_roof_grids(gridiv)
1844 tground(:) = tground_grids(:, gridiv)
1845 twall(:) = twall_grids(:, gridiv)
1846 troof(:) = troof_grids(:, gridiv)
1847 tibld(:) = tibld_grids(:, gridiv)
1848 tw_4 = tw_4_grids(:, :, gridiv)
1849
1850 END IF
1851
1852 ! Get Ts from Ts5min data array
1853 ! Tinternal = Ts5mindata(ir,cTs_Tiair)
1854 ! Tsurf_all = Ts5mindata(ir,cTs_Tsurf)
1855 ! Troof_in = Ts5mindata(ir,cTs_Troof)
1856 ! Troad = Ts5mindata(ir,cTs_Troad)
1857 ! Twall_all = Ts5mindata(ir,cTs_Twall)
1858 !
1859 ! Tw_n = Ts5mindata(ir,cTs_Twall_n)
1860 ! Tw_e = Ts5mindata(ir,cTs_Twall_e)
1861 ! Tw_s = Ts5mindata(ir,cTs_Twall_s)
1862 ! Tw_w = Ts5mindata(ir,cTs_Twall_w)
1863 tinternal = ts5mindata_ir(cts_tiair)
1864 tsurf_all = ts5mindata_ir(cts_tsurf)
1865 troof_in = ts5mindata_ir(cts_troof)
1866 troad = ts5mindata_ir(cts_troad)
1867 twall_all = ts5mindata_ir(cts_twall)
1868
1869 tw_n = ts5mindata_ir(cts_twall_n)
1870 tw_e = ts5mindata_ir(cts_twall_e)
1871 tw_s = ts5mindata_ir(cts_twall_s)
1872 tw_w = ts5mindata_ir(cts_twall_w)
1873
1874 ! if (any(isnan(Ts5mindata))) then !!FO!! can't use data when time gap is too big (or neg.) or data is NaN
1875 ! if (spindone) then !!FO!! writes a line of NaNs
1876 ! write(20,'(1F8.4,I6,100f10.1)') dectime,it,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,&
1877 ! -0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,&
1878 ! -0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,&
1879 ! -0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.
1880 ! endif
1881 ! return ! changed from cycle !!FO!! returns to beginning of do-loop, i.e. doesn't make any heat calculations
1882 ! endif
1883 !!FO!! this loop is used to run through the calculations (SPINUP number of times) to achieve better numerical stability
1884 !!FO!! when if condition is met the program starts over from the beginning of the input file and the calculations are performed
1885 !!FO!! ...once again, but this time saved to the output file
1886 !!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)
1887 ! IF (NLINESREAD>datalines.AND..NOT.SPINDONE) THEN !!FO!! at this stage spindone = .false.
1888 ! NLINESREAD=0
1889 ! SPINDONE=.TRUE.
1890 !PRINT*, "SPUNUP"
1891 ! return ! changed from cycle
1892 ! ENDIF
1893
1894 !! Write first row of each met block as -999
1895 !IF (first) THEN !Set to true in ESTM_initials
1896 !! !Tair2=Temp_C+C2K !This is now set in SUEWS_translate for ir=0 only
1897 !! ! first=.FALSE.
1898 ! IF(Gridiv == NumberOfGrids) first=.FALSE. !Set to false only after all grids have run
1899 ! dataOutESTM(ir,1:32,Gridiv)=(/REAL(iy,KIND(1D0)),REAL(id,KIND(1D0)),&
1900 ! REAL(it,KIND(1D0)),REAL(imin,KIND(1D0)),dectime,(dum(ii),ii=1,27)/)
1901 ! RETURN
1902 !ENDIF
1903
1904 ! What are these constants? - Need defining somewhere
1905 zenith_rad = zenith_deg/180*pi
1906 IF (zenith_rad > 0 .AND. zenith_rad < pi/2.-hw) THEN !ZENITH MUST BE HIGHER THAN BUILDINGS FOR DIRECT INTERCEPTION
1907 tanzenith = min(tan(zenith_rad), 5.67) !LIMITS TO ANGLES LESS THAN 80 EVEN FOR LOW HW
1908 tanzenith = tanzenith*kdn_estm/(1370*cos(zenith_rad)) !REDUCTION FACTOR FOR MAXIMUM
1909 ELSE
1910 tanzenith = 0.
1911 END IF
1912
1913 shc_air = heatcapacity_air(tair1, avrh, press_hpa) ! Use SUEWS version
1914
1915 !Evolution of building temperature from heat added by convection
1916 SELECT CASE (evolvetibld) !EvolveTiBld specifies which internal building temperature approach to use
1917 CASE (0)
1918 diagnoseti = .false.
1919 hvac = .false. !use data in file !!FO!! use measured indoor temperature (Tref in Lodz2002HS.txt)
1920 CASE (1) !!FO!! use of HVAC to counteract T changes
1921 diagnoseti = .true.
1922 IF (tievolve > theat_off) THEN !THEAT_OFF now converted to Kelvin in ESTM_initials - HCW 15 Jun 2016
1923 !IF (Tievolve>THEAT_OFF+C2K) THEN
1924 hvac = .false.
1925 ELSEIF (tievolve < theat_on) THEN !THEAT_OFF now converted to Kelvin in ESTM_initials - HCW 15 Jun 2016
1926 !ELSEIF (Tievolve<THEAT_ON+C2K) THEN
1927 hvac = .true.
1928 END IF
1929 CASE (2)
1930 diagnoseti = .true. !!FO!! convection between ibld and inside of external walls(?)
1931 END SELECT
1932
1933 !ASSUME AIR MIXES IN PROPORTION TO # OF EXCHANGES
1934 IF (tair_av > 20.+c2k .AND. tievolve > 25.+c2k .AND. tair1 < tievolve .AND. .NOT. hvac) THEN
1935 airexhr = 2.0 !Windows or exterior doors on 3 sides (ASHRAE 1981 22.8)
1936 ELSEIF (tair_av < 17.+c2k .OR. hvac) THEN
1937 airexhr = 0.5 !No window or exterior doors, storm sash or weathertripped (ASHRAE 1981 22.8)
1938 ELSE
1939 airexhr = 1.0
1940 END IF
1941
1942 airexdt = airexhr*(tstep/3600.0)
1943 shc_airbld = max(heatcapacity_air(tievolve, avrh, press_hpa), 0.00001) ! to avoid zero-division scenario TS 21 Oct 2017
1944 IF (shc_airbld < minshc_airbld) minshc_airbld = shc_airbld
1945
1946 !internal convective exchange coefficients !!FO!! ibldCHmod = 0 originally
1947 !iBldCHmod specifies method for convective exchange coeffs
1948 IF (ibldchmod == 1) THEN !ASHRAE 2001
1949 ch_ibld = 1.31*(abs(t0_ibld - tievolve))**0.25/shc_airbld
1950 ch_iwall = 1.31*(abs(tn_wall - tievolve))**0.25/shc_airbld
1951 ch_iroof = 1.52*(abs(tn_roof - tievolve))**0.25/shc_airbld
1952 IF (abs(tn_roof - tievolve) > 0) ch_iroof = ch_iroof*0.39 !effect of convection is weaker downward
1953 ELSEIF (ibldchmod == 2) THEN !Awbi, H.B. 1998, Energy and Buildings 28: 219-227
1954 ch_ibld = 1.823*(abs(t0_ibld - tievolve))**0.293/shc_airbld
1955 ch_iwall = 1.823*(abs(tn_wall - tievolve))**0.293/shc_airbld
1956 ch_iroof = 2.175*(abs(tn_roof - tievolve))**0.308/shc_airbld
1957 IF (abs(tn_roof - tievolve) > 0) ch_iroof = 0.704*(abs(tn_roof - tievolve))**0.133/shc_airbld !effect of convection is weaker downward
1958 END IF
1959
1960 !Evolving T = (Previous Temp + dT from Sensible heat flux) mixed with outside air
1961 !ASSUMES THE CH_BLD INCLUDES THE EFFECT OF VENTILATION RATE IN m/s (e.g. if a normal CH is .005 and
1962 !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
1963 tairmix = (tievolve + tair1*airexdt)/(1.0 + airexdt)
1964 qfbld = froof*(tievolve - tairmix)*shc_airbld*bldghx/tstep !heat added or lost, requires cooling or heating if HVAC on
1965
1966 !!FO!! CH_xxxx has unit [m/s] !!**HCW what is going on with tstep here??
1967 tievolve = tairmix + tstep/bldghx/finternal & !!FO!! finternal(=froof+fibld+fwall) => normalisation of fractions
1970 + ch_iwall*fwall*(tn_wall - tievolve)) !!FO!! [K] = [K] + [s/m]*([m/s]*([K]))
1971
1972 IF (.NOT. diagnoseti) tievolve = tinternal + c2k
1973 IF (hvac) THEN !Run up/down to set point +/- 1 degree with adjustment of 90% per hour
1974 tadd = (sign(-1.0d0, theat_fix - tievolve) + theat_fix - tievolve)*min(4.*tstep/3600.0, 0.9) !!**HCW check??
1975 tievolve = tievolve + tadd
1976 END IF
1977
1978 !========>RADIATION<================
1979 IF (kdn_estm < 0) kdn_estm = 0. !set non-zero shortwave to zero !Should this be moved up to line 183/4?
1980
1981 !external components, no diffuse
1982 !for reflections complete absorption is assumed
1983 !for shortwave these are net values
1984 !for longwave these are incoming only
1985 !MUST DIVIDE SHORTWAVE INTO DIRECT AND DIFFUSE
1986 sw_hor = kdn_estm !incoming solar on horizontal surface
1987 sw_vert = kdn_estm*tanzenith !incoming solar on vertical surface = kdown(obs)*sin(zenith)/cos(zenith)
1988
1989 rs_roof = svf_roof*(1.0 - alb_roof_estm)*sw_hor
1990 rl_roof = svf_roof*em_roof_estm*ldown
1991
1992 rs_ground = svf_ground*(1.-alb_ground_estm)*sw_hor + &
1995
1997
1998 rs_wall = svf_wall*(1.-alb_wall_fix)*sw_vert + &
2002
2003 !wall to wall exchange handled simultaneously with seb calc
2004 rl_wall = svf_wall*ldown*em_wall_fix + zvf_wall*svf_wall*ldown*(1 - em_wall_fix)*em_wall_fix + &
2006
2007 !DIFFICULT TO DETERMINE WHAT THIS IS EXACTLY, DONT INCLUDE WALLS
2008 kup_estm = kdn_estm - rvf_roof*rs_roof - (rvf_ground + rvf_wall)*rs_ground/svf_ground - rvf_veg*alb_veg_estm*kdn_estm
2009 IF (kdn_estm > 10 .AND. kup_estm > 0) THEN
2010 alb_avg = kup_estm/kdn_estm
2012 nalb = nalb + 1
2013 END IF
2014
2015 !internal components
2016 rs_ibld = 0 ! This could change if there are windows (need solar angles or wall svf * fraction glazing * transmissivity)
2017 !internal incoming longwave terms do not include the view factors for its own surface e.g. for ibld and walls
2018 !added floor view factors
2019 rl_ibld = sbconst*(ivf_iw*em_w*tn_wall**4 + &
2020 ivf_ir*em_r*tn_roof**4 + &
2021 ivf_if*em_f*tfloor**4)
2022 rs_iwall = 0
2023 rl_iwall = sbconst*(ivf_wi*em_i*t0_ibld**4 + &
2024 ivf_wr*em_r*tn_roof**4 + &
2025 ivf_wf*em_f*tfloor**4)
2026 rs_iroof = 0
2027 rl_iroof = sbconst*(ivf_ri*em_i*t0_ibld**4 + &
2028 ivf_rw*em_w*tn_wall**4 + &
2029 ivf_rf*em_f*tfloor**4)
2030
2031 !========>INTERNAL<================
2032 bctype = .false.
2033 kdz = 2*kibld(1)/zibld(1)
2034 pcoeff = (/em_ibld*sbconst*(1 - ivf_ii*em_ibld), 0.0d0, 0.0d0, kdz + shc_airbld*ch_ibld, &
2035 -kdz*tibld(1) - shc_airbld*ch_ibld*tievolve - rs_ibld - rl_ibld/)
2037
2038 !!FO!! this leads to Tibld(1) = Tibld(3) , i.e. ...
2039 bc(1) = t0_ibld
2040
2041 !!FO!! temperature equal on both sides of inside wall
2042 bc(2) = bc(1)
2043
2044 ! print*,'ESTM Qsibld cal before',Qsibld,Tibld,zibld(1:Nibld)
2046 REAL(Tstep, KIND(1D0)), kibld(1:Ndepth_ibld), &
2047 ribld(1:ndepth_ibld), bc, bctype)
2048
2049 ! print*,'ESTM Qsibld cal after',Qsibld,Tibld
2050
2051 !========>WALLS<================
2052 bctype = .false.
2054 pcoeff = (/em_ibld*sbconst*(1 - ivf_ww*em_ibld), 0.0d0, 0.0d0, kdz + shc_airbld*ch_iwall, &
2055 -kdz*twall(ndepth_wall) - shc_airbld*ch_iwall*tievolve - rs_iwall - rl_iwall/)
2057 bc(2) = tn_wall !!FO!! boundary condition #2 = inner surface Twall, originally from lodz_parms_ltm.txt or finaltemp.txt
2058
2059 IF (tsurfchoice < 2 .OR. radforce) THEN
2060 IF (radforce) THEN !!FO!! 1st prio: radforce
2061 kdz = 2*kwall(1)/zwall(1)
2062 pcoeff = (/em_wall_fix*sbconst*(1 - zvf_wall*em_wall_fix), 0.0d0, 0.0d0, kdz + shc_air*chair_wall*ws, &
2063 -kdz*twall(1) - shc_air*chair_wall*ws*tair1 - rs_wall - rl_wall/)
2065 bc(1) = t0_wall !!FO!! boundary condition #1 = outer surface Twall, originally from lodz_parms_ltm.txt or finaltemp.txt
2066 ELSEIF (tsurfchoice == 0) THEN
2067 bc(1) = tsurf_all + c2k; t0_wall = bc(1)
2068 ELSEIF (tsurfchoice == 1) THEN
2069 bc(1) = twall_all + c2k; t0_wall = bc(1)
2070 END IF !!FO!! Tsoil in Lodz2002HS.txt NB => Lodz2002HS.txt doesn't work with onewall = TRUE
2071
2072 CALL heatcond1d(twall, qswall, zwall(1:ndepth_wall), real(tstep, kind(1d0)), &
2073 kwall(1:ndepth_wall), rwall(1:ndepth_wall), bc, bctype) !!FO!! new set of Twalls are calculated from heat conduction through wall
2074
2075 ELSEIF (tsurfchoice == 2) THEN !SPECIAL FOR 4 WALLS
2076 t0_wall = 0.
2077 DO i = 1, 4 !do 4 walls
2078 bc(1) = tw_n + tw_e + tw_s + tw_w + c2k; t0_wall = t0_wall + bc(1)
2079 CALL heatcond1d(tw_4(:, i), qs_4(i), zwall(1:ndepth_wall), real(tstep, kind(1d0)), &
2081 END DO
2082 !Take average of 4 wall values
2083 t0_wall = t0_wall/4.
2084 qswall = sum(qs_4)/4.
2085 twall = sum(tw_4, 2)/4.
2086 END IF
2087
2088 !========>ROOF<================
2089 bctype = .false.
2091 pcoeff = (/em_ibld*sbconst, 0.0d0, 0.0d0, kdz + shc_airbld*ch_iroof, &
2092 -kdz*troof(ndepth_roof) - shc_airbld*ch_iroof*tievolve - rs_iroof - rl_iroof/)
2094 bc(2) = tn_roof
2095
2096 IF (radforce) THEN
2097 kdz = 2*kroof(1)/zroof(1)
2098 pcoeff = (/em_roof_estm*sbconst, 0.0d0, 0.0d0, kdz + shc_air*chair*ws, &
2099 -kdz*troof(1) - shc_air*chair*ws*tair1 - rs_roof - rl_roof/)
2101 bc(1) = t0_roof
2102 ELSEIF (tsurfchoice == 0) THEN
2103 bc(1) = tsurf_all + c2k; t0_roof = bc(1)
2104 ELSE
2105 bc(1) = troof_in + c2k; t0_roof = bc(1)
2106 END IF
2107
2108 CALL heatcond1d(troof, qsroof, zroof(1:ndepth_roof), real(tstep, kind(1d0)), &
2110
2111 !========>ground<================
2112 bctype = .false.
2113 kdz = 2*kground(1)/zground(1)
2114
2115 IF (radforce .OR. groundradforce) THEN
2116 pcoeff = (/em_ground_estm*sbconst, 0.0d0, 0.0d0, kdz + shc_air*chair_ground*ws, &
2117 -kdz*tground(1) - shc_air*chair_ground*ws*tair1 - rs_ground - rl_ground/)
2119 bc(1) = t0_ground
2120 ELSEIF (tsurfchoice == 0) THEN
2121 bc(1) = tsurf_all + c2k; t0_ground = bc(1)
2122 ELSE
2123 bc(1) = troad + c2k; t0_ground = bc(1)
2124 END IF
2125
2126 bc(2) = lbc_soil + c2k
2127 ! bc(2)=0.; bctype(2)=.t.
2128
2129 IF (fground /= 0.) THEN ! check fground==0 scenario to avoid division-by-zero error, TS 21 Jul 2016
2131 REAL(Tstep, KIND(1D0)), kground(1:Ndepth_ground), rground(1:Ndepth_ground), bc, bctype)
2132 ELSE
2133 qsground = nan
2134 END IF
2135
2136 qsair = fair*shc_air*(tair1 - tair2)/tstep
2141 qs = qsibld + qswall + qsroof + qsground !!FO!! QSair not included; called QS in output file (column #10)
2142
2143 ! print*,'ESTM QS',qs,Qsibld,Qswall,Qsroof ,Qsground
2144 ! print*,'ESTM Qsibld',Qsibld,fibld
2145 !write(*,*) Qsair, QSibld, Qswall, Qsroof, Qsground, QS
2146
2147 !========>Radiation<================
2148 !note that the LUP for individual components does not include reflected
2152 tveg = tair1
2153 lup_veg = sbconst*em_veg_estm*tveg**4
2156 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
2157 rn_ground = rs_ground + rl_ground - lup_ground
2158 rn_roof = rs_roof + rl_roof - lup_roof
2159 rn_wall = rs_wall + rl_wall - lup_wall*(1 - zvf_wall*em_wall_fix)
2160 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
2161 qhestm = (t0 - tair1)*chair*shc_air*ws
2162 sumemis = sumemis + em_equiv
2163 nemis = nemis + 1
2164
2165 ! IF (SPINDONE) THEN !!FO!! only the last set of values in the time interpolation loop is written to file
2166
2167 IF (ndepth_wall < 5) THEN
2168 twallout = (/twall, (dum(ii), ii=1, (5 - ndepth_wall))/)
2169 ELSE
2170 twallout = twall
2171 END IF
2172
2173 IF (ndepth_roof < 5) THEN
2174 troofout = (/troof, (dum(ii), ii=1, (5 - ndepth_roof))/)
2175 ELSE
2176 troofout = troof
2177 END IF
2178
2179 IF (ndepth_ground < 5) THEN
2180 tgroundout = (/tground, (dum(ii), ii=1, (5 - ndepth_ground))/)
2181 ELSE
2182 tgroundout = tground
2183 END IF
2184
2185 IF (ndepth_ibld < 5) THEN
2186 tibldout = (/tibld, (dum(ii), ii=1, (5 - ndepth_ibld))/)
2187 ELSE
2188 tibldout = tibld
2189 END IF
2190
2191 ! dataOutESTM(ir,1:ncolumnsDataOutESTM,Gridiv)=[&
2192 ! REAL(iy,KIND(1D0)),REAL(id,KIND(1D0)),REAL(it,KIND(1D0)),REAL(imin,KIND(1D0)), dectime,&!5
2193 ! QS,Qsair,Qswall,Qsroof,Qsground,Qsibld,&!11
2194 ! Twallout,Troofout,Tgroundout,Tibldout,Tievolve]!32 !NB. These all have 5 elements except Tievolve (1).
2195 dataoutlineestm = [ &
2196 qs, qsair, qswall, qsroof, qsground, qsibld, & !6
2197 twallout, troofout, tgroundout, tibldout, tievolve] !27 !NB. These all have 5 elements except Tievolve (1).
2198 ! set invalid values to nan
2199 dataoutlineestm = set_nan(dataoutlineestm)
2200 ! call r8vec_print(ncolumnsDataOutESTM-5,dataOutESTM(ir,6:ncolumnsDataOutESTM,Gridiv),'dataOutESTM')
2201
2202 tair2 = tair1
2203
2204 ! Save variables for this grid
2205 tair2_grids(gridiv) = tair1
2207 lup_wall_grids(gridiv) = lup_wall
2208 lup_roof_grids(gridiv) = lup_roof
2209 tievolve_grids(gridiv) = tievolve
2210 t0_ibld_grids(gridiv) = t0_ibld
2211 t0_ground_grids(gridiv) = t0_ground
2212 t0_wall_grids(gridiv) = t0_wall
2213 t0_roof_grids(gridiv) = t0_roof
2214 tn_wall_grids(gridiv) = tn_wall
2215 tn_roof_grids(gridiv) = tn_roof
2216 tground_grids(:, gridiv) = tground(:)
2217 twall_grids(:, gridiv) = twall(:)
2218 troof_grids(:, gridiv) = troof(:)
2219 tibld_grids(:, gridiv) = tibld(:)
2220 tw_4_grids(:, :, gridiv) = tw_4(:, :)
2221
integer tsurfchoice
real(kind(1d0)), dimension(:), allocatable lup_roof_grids
real(kind(1d0)) em_r
real(kind(1d0)) lbc_soil
real(kind(1d0)) theat_fix
real(kind(1d0)), dimension(:), allocatable t0_wall_grids
real(kind(1d0)), dimension(5) zground
real(kind(1d0)) em_roof_estm
real(kind(1d0)) ivf_ii
real(kind(1d0)) ivf_wf
integer ndepth_roof
real(kind(1d0)), dimension(:), allocatable t0_ground_grids
real(kind(1d0)) ivf_ri
real(kind(1d0)) froof
real(kind(1d0)) theat_off
real(kind(1d0)) lup_ground
real(kind(1d0)) t0_ground
real(kind(1d0)), dimension(:), allocatable tground
real(kind(1d0)) ivf_rf
real(kind(1d0)) rvf_roof
real(kind(1d0)) fair
real(kind(1d0)) zvf_ground
real(kind(1d0)) t0_ibld
logical diagnoseti
real(kind(1d0)) ch_iwall
integer ndepth_ibld
real(kind(1d0)) minshc_airbld
real(kind(1d0)), dimension(5) kwall
real(kind(1d0)) alb_ground_estm
real(kind(1d0)) tair1
real(kind(1d0)), dimension(:), allocatable twall
real(kind(1d0)) chr
real(kind(1d0)), dimension(5) kroof
real(kind(1d0)), dimension(5) kground
real(kind(1d0)) t0_roof
real(kind(1d0)) tn_roof
real(kind(1d0)) lup_roof
real(kind(1d0)) alb_veg_estm
integer ndepth_wall
real(kind(1d0)) tn_wall
real(kind(1d0)), dimension(:), allocatable tn_wall_grids
real(kind(1d0)) em_w
real(kind(1d0)), dimension(5) zibld
real(kind(1d0)), dimension(:), allocatable troof
real(kind(1d0)) tievolve
real(kind(1d0)) zvf_wall
real(kind(1d0)), dimension(4) qs_4
real(kind(1d0)) tanzenith
real(kind(1d0)), dimension(:), allocatable t0_ibld_grids
real(kind(1d0)), dimension(:), allocatable tair2_grids
real(kind(1d0)), dimension(5) rwall
real(kind(1d0)), dimension(:), allocatable t0_roof_grids
real(kind(1d0)) shc_air
real(kind(1d0)) fibld
real(kind(1d0)), dimension(:, :), allocatable tw_4
real(kind(1d0)), dimension(5) pcoeff
real(kind(1d0)) fwall
real(kind(1d0)) alb_roof_estm
real(kind(1d0)) em_i
real(kind(1d0)) qsibld
real(kind(1d0)) rvf_ground
real(kind(1d0)) lup_veg
real(kind(1d0)), dimension(:, :, :), allocatable tw_4_grids
real(kind(1d0)) lup_wall
real(kind(1d0)), dimension(:), allocatable lup_ground_grids
real(kind(1d0)) ch_ibld
real(kind(1d0)) fground
real(kind(1d0)) ivf_wi
real(kind(1d0)), dimension(:, :), allocatable twall_grids
real(kind(1d0)) qsroof
real(kind(1d0)) ivf_if
real(kind(1d0)) chair
real(kind(1d0)) tfloor
real(kind(1d0)) em_ibld
real(kind(1d0)) ch_iroof
real(kind(1d0)), dimension(5) rroof
real(kind(1d0)), dimension(:), allocatable lup_wall_grids
real(kind(1d0)) ivf_ww
real(kind(1d0)) qswall
real(kind(1d0)) qsground
real(kind(1d0)) finternal
real(kind(1d0)) ivf_ir
real(kind(1d0)) sumemis
real(kind(1d0)) rvf_veg
real(kind(1d0)) em_f
real(kind(1d0)) tair2
real(kind(1d0)) svf_wall
real(kind(1d0)), dimension(:), allocatable tibld
integer ibldchmod
real(kind(1d0)) ivf_iw
real(kind(1d0)), dimension(:), allocatable tn_roof_grids
integer evolvetibld
real(kind(1d0)) ivf_rw
real(kind(1d0)), parameter alb_wall_fix
real(kind(1d0)), dimension(5) kibld
real(kind(1d0)) em_veg_estm
integer ndepth_ground
real(kind(1d0)) em_ground_estm
real(kind(1d0)), parameter em_wall_fix
real(kind(1d0)), dimension(5) zwall
real(kind(1d0)), dimension(:, :), allocatable tibld_grids
real(kind(1d0)), parameter conv
real(kind(1d0)) hw
real(kind(1d0)) ivf_wr
real(kind(1d0)) alb_avg
real(kind(1d0)) svf_roof
real(kind(1d0)) qsair
real(kind(1d0)), dimension(:, :), allocatable troof_grids
real(kind(1d0)), dimension(:), allocatable tievolve_grids
real(kind(1d0)) xvf_wall
real(kind(1d0)) t0_wall
real(kind(1d0)), dimension(5) ribld
integer, parameter maxiter
real(kind(1d0)) sumalb
real(kind(1d0)), dimension(5) zroof
logical, dimension(2) bctype
real(kind(1d0)) ws
real(kind(1d0)) theat_on
real(kind(1d0)) svf_ground
real(kind(1d0)) rvf_wall
real(kind(1d0)), dimension(:, :), allocatable tground_grids
subroutine heatcond1d(T, Qs, dx, dt, k, rhocp, bc, bctype)
real(kind(1d0)) function heatcapacity_air(TK, RH, P)
real(kind(1d0)) function newtonpolynomial(x0, Pcoeff, conv, maxiter)
real(kind(1d0)), parameter c2k
real(kind(1d0)), parameter sbconst

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().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ estm_ext()

subroutine estm_module::estm_ext ( integer, intent(in)  tstep,
integer, intent(in)  nlayer,
real(kind(1d0)), dimension(nsurf), intent(in)  QG_surf,
real(kind(1d0)), dimension(nlayer), intent(in)  qg_roof,
real(kind(1d0)), dimension(nlayer), intent(in)  qg_wall,
real(kind(1d0)), dimension(nlayer), intent(in)  tsfc_roof,
real(kind(1d0)), dimension(nlayer), intent(in)  tin_roof,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in)  temp_in_roof,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in)  k_roof,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in)  cp_roof,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in)  dz_roof,
real(kind(1d0)), dimension(nlayer), intent(in)  sfr_roof,
real(kind(1d0)), dimension(nlayer), intent(in)  tsfc_wall,
real(kind(1d0)), dimension(nlayer), intent(in)  tin_wall,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in)  temp_in_wall,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in)  k_wall,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in)  cp_wall,
real(kind(1d0)), dimension(nlayer, ndepth), intent(in)  dz_wall,
real(kind(1d0)), dimension(nlayer), intent(in)  sfr_wall,
real(kind(1d0)), dimension(nsurf), intent(in)  tsfc_surf,
real(kind(1d0)), dimension(nsurf), intent(in)  tin_surf,
real(kind(1d0)), dimension(nsurf, ndepth), intent(in)  temp_in_surf,
real(kind(1d0)), dimension(nsurf, ndepth), intent(in)  k_surf,
real(kind(1d0)), dimension(nsurf, ndepth), intent(in)  cp_surf,
real(kind(1d0)), dimension(nsurf, ndepth), intent(in)  dz_surf,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr_surf,
real(kind(1d0)), dimension(nlayer, ndepth), intent(out)  temp_out_roof,
real(kind(1d0)), dimension(nlayer), intent(out)  QS_roof,
real(kind(1d0)), dimension(nlayer, ndepth), intent(out)  temp_out_wall,
real(kind(1d0)), dimension(nlayer), intent(out)  QS_wall,
real(kind(1d0)), dimension(nsurf, ndepth), intent(out)  temp_out_surf,
real(kind(1d0)), dimension(nsurf), intent(out)  QS_surf,
real(kind(1d0)), intent(out)  QS 
)

Definition at line 2229 of file suews_phys_estm.f95.

2240 USE allocatearray, ONLY: &
2241 nsurf, ndepth, &
2243 USE heatflux, ONLY: heatcond1d_ext
2244
2245 IMPLICIT NONE
2246 INTEGER, INTENT(in) :: tstep
2247 INTEGER, INTENT(in) :: nlayer ! number of vertical levels in urban canopy
2248
2249 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: QG_surf ! ground heat flux
2250 ! extended for ESTM_ext
2251
2252 ! keys:
2253 ! tsfc: surface temperature
2254 ! tin: indoor/deep bottom temperature
2255 ! temp_in: temperature at inner interfaces
2256 ! k: thermal conductivity
2257 ! cp: heat capacity
2258 ! dz: thickness of each layer
2259 ! roof/wall/surf: roof/wall/ground surface types
2260
2261 ! input arrays: roof facets
2262 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: qg_roof
2263 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: tsfc_roof
2264 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: tin_roof
2265 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: sfr_roof
2266 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: temp_in_roof
2267 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: k_roof
2268 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: cp_roof
2269 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: dz_roof
2270 ! input arrays: wall facets
2271 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: qg_wall
2272 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: tsfc_wall
2273 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: tin_wall
2274 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(in) :: sfr_wall
2275 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: temp_in_wall
2276 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: k_wall
2277 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: cp_wall
2278 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(in) :: dz_wall
2279 ! input arrays: standard suews surfaces
2280 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: tsfc_surf
2281 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: tin_surf
2282 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf
2283 REAL(KIND(1D0)), DIMENSION(nsurf, ndepth), INTENT(in) :: temp_in_surf
2284 REAL(KIND(1D0)), DIMENSION(nsurf, ndepth), INTENT(in) :: k_surf
2285 REAL(KIND(1D0)), DIMENSION(nsurf, ndepth), INTENT(in) :: cp_surf
2286 REAL(KIND(1D0)), DIMENSION(nsurf, ndepth), INTENT(in) :: dz_surf
2287
2288 ! output arrays
2289 ! roof facets
2290 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(out) :: QS_roof
2291 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(out) :: temp_out_roof !interface temperature between depth layers
2292 ! wall facets
2293 REAL(KIND(1D0)), DIMENSION(nlayer), INTENT(out) :: QS_wall
2294 REAL(KIND(1D0)), DIMENSION(nlayer, ndepth), INTENT(out) :: temp_out_wall !interface temperature between depth layers
2295 ! standard suews surfaces
2296 REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: QS_surf
2297 REAL(KIND(1D0)), DIMENSION(nsurf, ndepth), INTENT(out) :: temp_out_surf !interface temperature between depth layers
2298
2299 ! grid aggregated results
2300 REAL(KIND(1D0)), INTENT(out) :: QS
2301
2302 ! internal use
2303 ! temporary arrays in calculations
2304 ! note: nsurf is used as the maximum number of surfaces
2305 REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: tsfc_cal
2306 REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: tin_cal
2307 REAL(KIND(1D0)), DIMENSION(:, :), ALLOCATABLE :: temp_cal
2308 REAL(KIND(1D0)), DIMENSION(:, :), ALLOCATABLE :: k_cal
2309 REAL(KIND(1D0)), DIMENSION(:, :), ALLOCATABLE :: cp_cal
2310 REAL(KIND(1D0)), DIMENSION(:, :), ALLOCATABLE :: dz_cal
2311 REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: qs_cal
2312
2313 ! REAL(KIND(1D0)), DIMENSION(ndepth) :: temp_all_cal
2314 ! surface temperatures at innermost depth layer
2315 ! REAL(KIND(1D0)) :: temp_IBC
2316 INTEGER :: i_facet, i_group, nfacet, i_depth
2317
2318 ! settings for boundary conditions
2319 REAL(KIND(1D0)), DIMENSION(2) :: bc
2320
2321 ! use temperature as boundary conditions
2322 LOGICAL, DIMENSION(2) :: bctype
2323 LOGICAL :: debug
2324
2325 ! use finite depth heat conduction solver
2326 LOGICAL :: use_heatcond1d, use_heatcond1d_water
2327
2328 ! normalised surface fractions
2329 REAL(KIND(1D0)), DIMENSION(nlayer) :: sfr_roof_n
2330 REAL(KIND(1D0)), DIMENSION(nlayer) :: sfr_wall_n
2331
2332 ! initialise solver flags
2333 use_heatcond1d = .true.
2334 use_heatcond1d_water = .false.
2335 debug = .false.
2336
2337 ! normalised surface fractions
2338 ! NB: the sum of sfr_roof (sfr_wall) is NOT unity; so need to be normalised by the sum
2339 sfr_roof_n = sfr_roof/sum(sfr_roof)
2340 sfr_wall_n = sfr_wall/sum(sfr_wall)
2341 ! sub-facets of buildings: e.g. walls, roofs, etc.
2342 DO i_group = 1, 3
2343
2344 ! allocate arrays
2345 IF (i_group == 1) THEN
2346 ! PRINT *, 'group: ', 'roof'
2347 nfacet = nlayer
2348 ELSE IF (i_group == 2) THEN
2349 ! PRINT *, 'group: ', 'wall'
2350 nfacet = nlayer
2351 ELSE IF (i_group == 3) THEN
2352 ! PRINT *, 'group: ', 'surf'
2353 nfacet = nsurf
2354 END IF
2355 ! PRINT *, 'nfacet here: ', nfacet
2356 ALLOCATE (tsfc_cal(nfacet))
2357 ALLOCATE (tin_cal(nfacet))
2358 ALLOCATE (qs_cal(nfacet))
2359 ALLOCATE (temp_cal(nfacet, ndepth))
2360 ALLOCATE (k_cal(nfacet, ndepth))
2361 ALLOCATE (cp_cal(nfacet, ndepth))
2362 ALLOCATE (dz_cal(nfacet, ndepth))
2363 ! PRINT *, 'allocation done! '
2364
2365 ! translate input arrays of facet groups to internal use arrays
2366 IF (i_group == 1) THEN
2367 ! PRINT *, 'translation for roof! '
2368 ! TODO: to update with actual values from input files
2369 tsfc_cal(1:nfacet) = tsfc_roof(1:nfacet)
2370 ! PRINT *, 'tsfc_cal for roof! ', tsfc_cal
2371 tin_cal(1:nfacet) = tin_roof(1:nfacet)
2372 ! PRINT *, 'tin_cal for roof! ', tin_cal
2373 temp_cal(1:nfacet, 1:ndepth) = temp_in_roof(1:nfacet, 1:ndepth)
2374 ! PRINT *, 'temp_cal for roof! ',temp_cal
2375 k_cal(1:nfacet, 1:ndepth) = k_roof(1:nfacet, 1:ndepth)
2376 ! PRINT *, 'k_roof for roof! ',k_roof(1,:)
2377 cp_cal(1:nfacet, 1:ndepth) = cp_roof(1:nfacet, 1:ndepth)
2378 dz_cal(1:nfacet, 1:ndepth) = dz_roof(1:nfacet, 1:ndepth)
2379 ! qs_cal(1:nfacet) = qs_roof(1:nfacet)
2380 ELSE IF (i_group == 2) THEN
2381 ! PRINT *, 'translation for wall! '
2382 ! TODO: to update with actual values from input files
2383 tsfc_cal(1:nfacet) = tsfc_wall(1:nfacet)
2384 tin_cal(1:nfacet) = tin_wall(1:nfacet)
2385 temp_cal(1:nfacet, 1:ndepth) = temp_in_wall(1:nfacet, 1:ndepth)
2386 ! TODO: temporarily set this for testing
2387 k_cal(1:nfacet, 1:ndepth) = k_wall(1:nfacet, 1:ndepth)
2388 cp_cal(1:nfacet, 1:ndepth) = cp_wall(1:nfacet, 1:ndepth)
2389 dz_cal(1:nfacet, 1:ndepth) = dz_wall(1:nfacet, 1:ndepth)
2390 ! qs_cal(1:nfacet) = qs_wall(1:nfacet)
2391
2392 ELSE IF (i_group == 3) THEN
2393 ! PRINT *, 'translation for surf! '
2394 ! nfacet = nsurf
2395 tsfc_cal(1:nfacet) = tsfc_surf(1:nfacet)
2396 tin_cal(1:nfacet) = tin_surf(1:nfacet)
2397 temp_cal(1:nfacet, 1:ndepth) = temp_in_surf(1:nfacet, 1:ndepth)
2398 ! TODO: to update with actual values from input files
2399 k_cal(1:nfacet, 1:ndepth) = k_surf(1:nfacet, 1:ndepth)
2400 cp_cal(1:nfacet, 1:ndepth) = cp_surf(1:nfacet, 1:ndepth)
2401 dz_cal(1:nfacet, 1:ndepth) = dz_surf(1:nfacet, 1:ndepth)
2402 ! k_cal(1:nfacet, 1:ndepth) = 1.2
2403 ! cp_cal(1:nfacet, 1:ndepth) = 2E6
2404 ! dz_cal(1:nfacet, 1:ndepth) = 0.1
2405 ! qs_cal(1:nfacet) = QS_surf(1:nfacet)
2406 END IF
2407 ! PRINT *, 'translation done! '
2408 ! TODO: temporary setting
2409 ! k_cal(1:nfacet, 1:ndepth) = 1.2
2410 ! cp_cal(1:nfacet, 1:ndepth) = 2E6
2411 ! dz_cal(1:nfacet, 1:ndepth) = 0.2
2412
2413 ! PRINT *, 'nfacet: ', nfacet
2414 DO i_facet = 1, nfacet
2415 ! PRINT *, 'i_facet: ', i_facet
2416 ! ASSOCIATE (v => dz_cal(i_facet, 1:ndepth))
2417 ! PRINT *, 'dz_cal in estm_ext', v, SIZE(v)
2418 ! END ASSOCIATE
2419
2420 ! determine the calculation method
2421 IF (i_group == 3) THEN
2422 use_heatcond1d = .true.
2423 IF (i_facet == bldgsurf) THEN
2424 ! building surface needs a different treatment
2425 use_heatcond1d = .false.
2426 ELSE IF (i_facet == watersurf) THEN
2427 ! water surface needs a different treatment
2428 use_heatcond1d = .false.
2429 use_heatcond1d_water = .true.
2430 END IF
2431 ELSE
2432 use_heatcond1d = .true.
2433 END IF
2434
2435 ! actual heat conduction calculations
2436 IF (dz_cal(i_facet, 1) /= -999.0 .AND. use_heatcond1d) THEN
2437
2438 ! surface heat flux
2439 IF (i_group == 1) THEN
2440 bc(1) = qg_roof(i_facet)
2441 ELSE IF (i_group == 2) THEN
2442 bc(1) = qg_wall(i_facet)
2443 ELSE IF (i_group == 3) THEN
2444 bc(1) = qg_surf(i_facet)
2445 END IF
2446 ! bctype(1) = .TRUE.
2447 bc(1) = tsfc_cal(i_facet)
2448 bctype(1) = .false.
2449
2450 !TODO: should be a prescribed temperature of the innermost boundary
2451 bc(2) = tin_cal(i_facet)
2452 bctype(2) = .false.
2453
2454 ! IF (i_group == 3 .AND. i_facet == 3) THEN
2455 ! PRINT *, 'temp_cal before: ', temp_cal(i_facet, :)
2456 ! PRINT *, 'k_cal: ', k_cal(i_facet, 1:ndepth)
2457 ! PRINT *, 'cp_cal: ', cp_cal(i_facet, 1:ndepth)
2458 ! PRINT *, 'dz_cal: ', dz_cal(i_facet, 1:ndepth)
2459 ! PRINT *, 'bc: ', bc
2460
2461 ! END IF
2462
2463 ! 1D heat conduction for finite depth layers
2464
2465 ! IF ((i_group == 3) .AND. (i_facet == 1)) THEN
2466 ! debug = .FALSE.
2467 ! ELSE
2468 ! debug = .FALSE.
2469 ! END IF
2470 CALL heatcond1d_ext( &
2471 temp_cal(i_facet, :), &
2472 qs_cal(i_facet), &
2473 tsfc_cal(i_facet), &
2474 dz_cal(i_facet, 1:ndepth), &
2475 tstep*1.d0, &
2476 k_cal(i_facet, 1:ndepth), &
2477 cp_cal(i_facet, 1:ndepth), &
2478 bc, &
2479 bctype, debug)
2480
2481 ! update temperature at all inner interfaces
2482 ! tin_cal(i_facet, :) = temp_all_cal
2483 ! IF (i_group == 3 .AND. i_facet == 3) THEN
2484 ! PRINT *, 'temp_cal after: ', temp_cal(i_facet, :)
2485 ! PRINT *, 'QS_cal after: ', QS_cal(i_facet)
2486 ! PRINT *, 'k_cal: ', k_cal(i_facet, 1:ndepth)
2487 ! PRINT *, 'cp_cal: ', cp_cal(i_facet, 1:ndepth)
2488 ! PRINT *, 'dz_cal: ', dz_cal(i_facet, 1:ndepth)
2489 ! PRINT *, 'bc: ', bc
2490 ! PRINT *, ''
2491
2492 ! END IF
2493 END IF
2494
2495 IF (dz_cal(i_facet, 1) /= -999.0 .AND. use_heatcond1d_water) THEN
2496 ! temperatures at all interfaces, including the outmost surface
2497 ! temp_all_cal = temp_cal(i_facet, :)
2498
2499 ! outermost surface temperature
2500 ! bc(1) = tsfc_cal(i_facet)
2501 bc(1) = tsfc_cal(i_facet)
2502 bctype(1) = .false.
2503
2504 !TODO: should be a prescribed temperature of the innermost boundary
2505 bc(2) = tin_cal(i_facet)
2506 bctype(2) = .false.
2507
2508 ! 1D heat conduction for finite depth layers
2509 ! TODO: this should be a water specific heat conduction solver: to implement
2510 CALL heatcond1d_ext( &
2511 temp_cal(i_facet, :), &
2512 qs_cal(i_facet), &
2513 tsfc_cal(i_facet), &
2514 dz_cal(i_facet, 1:ndepth), &
2515 tstep*1.d0, &
2516 k_cal(i_facet, 1:ndepth), &
2517 cp_cal(i_facet, 1:ndepth), &
2518 bc, &
2519 bctype, debug)
2520
2521 ! ! update temperature at all inner interfaces
2522 ! temp_cal(i_facet, :) = temp_all_cal
2523 END IF
2524
2525 END DO ! end of i_facet loop
2526
2527 ! translate results to back to the output arrays of facet groups
2528 IF (i_group == 1) THEN
2529 qs_roof = qs_cal(1:nfacet)
2530 ! tsfc_roof = tsfc_cal(1:nfacet)
2531 temp_out_roof = temp_cal(:nfacet, :)
2532 ELSE IF (i_group == 2) THEN
2533 qs_wall = qs_cal(1:nfacet)
2534 ! tsfc_wall = tsfc_cal(1:nfacet)
2535 temp_out_wall = temp_cal(:nfacet, :)
2536 ELSE IF (i_group == 3) THEN
2537 qs_surf = qs_cal(1:nfacet)
2538 ! tsfc_surf = tsfc_cal(1:nfacet)
2539 temp_out_surf = temp_cal(:nfacet, :)
2540 END IF
2541
2542 ! deallocate memory
2543 DEALLOCATE (tsfc_cal)
2544 DEALLOCATE (tin_cal)
2545 DEALLOCATE (qs_cal)
2546 DEALLOCATE (temp_cal)
2547 DEALLOCATE (k_cal)
2548 DEALLOCATE (cp_cal)
2549 DEALLOCATE (dz_cal)
2550
2551 END DO ! end do i_group
2552
2553 ! aggregated results
2554 ! building surface
2555 ! PRINT *, 'QS_roof in estm_ext', DOT_PRODUCT(QS_roof, sfr_roof), 'for', sfr_roof
2556 ! PRINT *, 'QS_wall in estm_ext', DOT_PRODUCT(QS_wall, sfr_wall), 'for', sfr_wall
2557 qs_surf(bldgsurf) = (dot_product(qs_roof, sfr_roof) + dot_product(qs_wall, sfr_wall))/sfr_surf(bldgsurf)
2558 DO i_depth = 1, ndepth
2559 temp_out_surf(bldgsurf, i_depth) = &
2560 (dot_product(temp_out_roof(:, i_depth), sfr_roof) &
2561 + dot_product(temp_out_wall(:, i_depth), sfr_wall)) &
2562 /(sum(sfr_roof) + sum(sfr_wall))
2563 END DO
2564
2565 ! all standard suews surfaces
2566 qs = dot_product(qs_surf, sfr_surf)
2567
integer, parameter bldgsurf
real(kind(1d0)), dimension(:, :), allocatable dz_wall
integer, parameter conifsurf
real(kind(1d0)), dimension(:), allocatable tsfc_roof
real(kind(1d0)), dimension(:, :), allocatable dz_surf
real(kind(1d0)), dimension(nsurf) sfr_surf
real(kind(1d0)), dimension(:), allocatable sfr_roof
real(kind(1d0)), dimension(:, :), allocatable cp_wall
real(kind(1d0)), dimension(:), allocatable tsfc_surf
real(kind(1d0)), dimension(:, :), allocatable cp_surf
real(kind(1d0)), dimension(:, :), allocatable k_roof
real(kind(1d0)), dimension(:, :), allocatable k_surf
real(kind(1d0)), dimension(:), allocatable tin_surf
integer, parameter watersurf
real(kind(1d0)), dimension(:), allocatable tin_wall
integer, parameter grasssurf
real(kind(1d0)), dimension(:, :), allocatable dz_roof
integer, parameter bsoilsurf
integer, parameter nsurf
integer, parameter pavsurf
real(kind(1d0)), dimension(:, :), allocatable cp_roof
integer, parameter ndepth
real(kind(1d0)), dimension(:, :), allocatable k_wall
real(kind(1d0)), dimension(:), allocatable sfr_wall
real(kind(1d0)), dimension(:), allocatable tin_roof
integer, parameter decidsurf
real(kind(1d0)), dimension(:), allocatable tsfc_wall
recursive subroutine heatcond1d_ext(T, Qs, Tsfc, dx, dt, k, rhocp, bc, bctype, debug)

References allocatearray::bldgsurf, allocatearray::bsoilsurf, allocatearray::conifsurf, allocatearray::decidsurf, allocatearray::grasssurf, heatflux::heatcond1d_ext(), allocatearray::ndepth, allocatearray::nsurf, allocatearray::pavsurf, and allocatearray::watersurf.

Referenced by suews_driver::suews_cal_qs().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ estm_ext_finalise()

subroutine estm_module::estm_ext_finalise

Definition at line 1194 of file suews_phys_estm.f95.

1195 USE allocatearray
1196 IMPLICIT NONE
1197
1198 IF (ALLOCATED(nlayer_grids)) DEALLOCATE (nlayer_grids)
1199
1200 IF (ALLOCATED(height_grids)) DEALLOCATE (height_grids)
1201 IF (ALLOCATED(building_frac_grids)) DEALLOCATE (building_frac_grids)
1202 IF (ALLOCATED(veg_frac_grids)) DEALLOCATE (veg_frac_grids)
1203 IF (ALLOCATED(building_scale_grids)) DEALLOCATE (building_scale_grids)
1204 IF (ALLOCATED(veg_scale_grids)) DEALLOCATE (veg_scale_grids)
1205
1206 IF (ALLOCATED(sfr_roof_grids)) DEALLOCATE (sfr_roof_grids)
1207 IF (ALLOCATED(alb_roof_grids)) DEALLOCATE (alb_roof_grids)
1208 IF (ALLOCATED(emis_roof_grids)) DEALLOCATE (emis_roof_grids)
1209 IF (ALLOCATED(k_roof_grids)) DEALLOCATE (k_roof_grids)
1210 IF (ALLOCATED(cp_roof_grids)) DEALLOCATE (cp_roof_grids)
1211 IF (ALLOCATED(dz_roof_grids)) DEALLOCATE (dz_roof_grids)
1212 IF (ALLOCATED(tin_roof_grids)) DEALLOCATE (tin_roof_grids)
1213 IF (ALLOCATED(temp_roof_grids)) DEALLOCATE (temp_roof_grids)
1214 IF (ALLOCATED(state_roof_grids)) DEALLOCATE (state_roof_grids)
1215 IF (ALLOCATED(statelimit_roof_grids)) DEALLOCATE (statelimit_roof_grids)
1216 IF (ALLOCATED(wetthresh_roof_grids)) DEALLOCATE (wetthresh_roof_grids)
1217 IF (ALLOCATED(soilstore_roof_grids)) DEALLOCATE (soilstore_roof_grids)
1218 IF (ALLOCATED(soilstorecap_roof_grids)) DEALLOCATE (soilstorecap_roof_grids)
1220
1221 IF (ALLOCATED(sfr_wall_grids)) DEALLOCATE (sfr_wall_grids)
1222 IF (ALLOCATED(alb_wall_grids)) DEALLOCATE (alb_wall_grids)
1223 IF (ALLOCATED(emis_wall_grids)) DEALLOCATE (emis_wall_grids)
1224 IF (ALLOCATED(k_wall_grids)) DEALLOCATE (k_wall_grids)
1225 IF (ALLOCATED(cp_wall_grids)) DEALLOCATE (cp_wall_grids)
1226 IF (ALLOCATED(dz_wall_grids)) DEALLOCATE (dz_wall_grids)
1227 IF (ALLOCATED(tin_wall_grids)) DEALLOCATE (tin_wall_grids)
1228 IF (ALLOCATED(temp_wall_grids)) DEALLOCATE (temp_wall_grids)
1229 IF (ALLOCATED(state_wall_grids)) DEALLOCATE (state_wall_grids)
1230 IF (ALLOCATED(statelimit_wall_grids)) DEALLOCATE (statelimit_wall_grids)
1231 IF (ALLOCATED(wetthresh_wall_grids)) DEALLOCATE (wetthresh_wall_grids)
1232 IF (ALLOCATED(soilstore_wall_grids)) DEALLOCATE (soilstore_wall_grids)
1233 IF (ALLOCATED(soilstorecap_wall_grids)) DEALLOCATE (soilstorecap_wall_grids)
1234 IF (ALLOCATED(wall_specular_frac_grids)) DEALLOCATE (wall_specular_frac_grids)
1235
1236 IF (ALLOCATED(k_surf_grids)) DEALLOCATE (k_surf_grids)
1237 IF (ALLOCATED(cp_surf_grids)) DEALLOCATE (cp_surf_grids)
1238 IF (ALLOCATED(dz_surf_grids)) DEALLOCATE (dz_surf_grids)
1239 ! IF (ALLOCATED(tin_surf_grids)) DEALLOCATE (tin_surf_grids)
1240 IF (ALLOCATED(temp_surf_grids)) DEALLOCATE (temp_surf_grids)
1241
real(kind(1d0)), dimension(:, :), allocatable soilstorecap_wall_grids
real(kind(1d0)), dimension(:, :), allocatable soilstore_wall_grids
real(kind(1d0)), dimension(:, :, :), allocatable cp_wall_grids
real(kind(1d0)), dimension(:, :), allocatable emis_wall_grids
real(kind(1d0)), dimension(:, :, :), allocatable temp_surf_grids
real(kind(1d0)), dimension(:, :), allocatable tin_roof_grids
real(kind(1d0)), dimension(:, :), allocatable building_frac_grids
real(kind(1d0)), dimension(:, :, :), allocatable temp_wall_grids
real(kind(1d0)), dimension(:, :), allocatable tin_wall_grids
integer, dimension(:), allocatable nlayer_grids
real(kind(1d0)), dimension(:, :), allocatable alb_wall_grids
real(kind(1d0)), dimension(:, :, :), allocatable cp_surf_grids
real(kind(1d0)), dimension(:, :, :), allocatable wall_specular_frac_grids
real(kind(1d0)), dimension(:, :), allocatable alb_roof_grids
real(kind(1d0)), dimension(:, :, :), allocatable k_wall_grids
real(kind(1d0)), dimension(:, :), allocatable emis_roof_grids
real(kind(1d0)), dimension(:, :, :), allocatable temp_roof_grids
real(kind(1d0)), dimension(:, :, :), allocatable cp_roof_grids
real(kind(1d0)), dimension(:, :), allocatable sfr_wall_grids
real(kind(1d0)), dimension(:, :), allocatable state_wall_grids
real(kind(1d0)), dimension(:, :), allocatable building_scale_grids
real(kind(1d0)), dimension(:, :, :), allocatable dz_surf_grids
real(kind(1d0)), dimension(:, :), allocatable soilstore_roof_grids
real(kind(1d0)), dimension(:, :, :), allocatable k_roof_grids
real(kind(1d0)), dimension(:, :), allocatable statelimit_wall_grids
real(kind(1d0)), dimension(:, :), allocatable wetthresh_roof_grids
real(kind(1d0)), dimension(:, :), allocatable soilstorecap_roof_grids
real(kind(1d0)), dimension(:, :, :), allocatable dz_roof_grids
real(kind(1d0)), dimension(:, :, :), allocatable roof_albedo_dir_mult_fact_grids
real(kind(1d0)), dimension(:, :), allocatable height_grids
real(kind(1d0)), dimension(:, :), allocatable veg_scale_grids
real(kind(1d0)), dimension(:, :, :), allocatable k_surf_grids
real(kind(1d0)), dimension(:, :), allocatable veg_frac_grids
real(kind(1d0)), dimension(:, :), allocatable state_roof_grids
real(kind(1d0)), dimension(:, :), allocatable statelimit_roof_grids
real(kind(1d0)), dimension(:, :), allocatable sfr_roof_grids
real(kind(1d0)), dimension(:, :, :), allocatable dz_wall_grids
real(kind(1d0)), dimension(:, :), allocatable wetthresh_wall_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.

Referenced by suews_program().

Here is the caller graph for this function:

◆ estm_ext_initialise()

subroutine estm_module::estm_ext_initialise

Definition at line 1114 of file suews_phys_estm.f95.

1115
1116 ! Last modified HCW 30 Jun 2016 - reading in now done by SUEWS_GetESTMData subroutine.
1117 ! ESTM_initials now only runs once per run at the very start.
1118 ! Last modified HCW 15 Jun 2016 - code now reads in 5-min file (interpolation done beforehand, outside of SUEWS itself)
1119
1120 ! USE defaultNotUsed
1121 USE physconstants, ONLY: c2k
1122 ! USE ESTM_data
1123 USE allocatearray
1125 USE initial, ONLY: numberofgrids
1126
1127 IMPLICIT NONE
1128 LOGICAL :: flag_mutiple_layout_files
1129 INTEGER :: i_grid
1130
1131 IF (multiplelayoutfiles == 0) THEN
1132 flag_mutiple_layout_files = .false.
1133 ELSE IF (multiplelayoutfiles == 1) THEN
1134 flag_mutiple_layout_files = .true.
1135 END IF
1136
1137 ALLOCATE (nlayer_grids(numberofgrids))
1138
1144
1145 ! allocate roof variables
1161
1162 ! allocate wall variables
1178
1179 ! allocate surf variables
1180 ! ALLOCATE (tsfc_surf_grids(NumberOfGrids, nsurf))
1186
1187 DO i_grid = 1, numberofgrids
1188 print *, 'Reading layout data for grid ', i_grid
1189 CALL load_gridlayout(i_grid, flag_mutiple_layout_files, diagnose)
1190 END DO
1191
real(kind(1d0)), dimension(:, :), allocatable tsfc_wall_grids
real(kind(1d0)), dimension(:, :), allocatable tin_surf_grids
real(kind(1d0)), dimension(:, :), allocatable tsfc_roof_grids
integer, parameter nspec
integer, parameter nlayer_max
integer multiplelayoutfiles
integer diagnose
integer numberofgrids

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.

Referenced by suews_program().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ estm_initials()

subroutine estm_module::estm_initials

Definition at line 696 of file suews_phys_estm.f95.

697
698 ! Last modified HCW 30 Jun 2016 - reading in now done by SUEWS_GetESTMData subroutine.
699 ! ESTM_initials now only runs once per run at the very start.
700 ! Last modified HCW 15 Jun 2016 - code now reads in 5-min file (interpolation done beforehand, outside of SUEWS itself)
701
703 USE physconstants, ONLY: c2k
704 USE estm_data
705 USE allocatearray
706 USE data_in, ONLY: fileinputpath
707 USE initial, ONLY: numberofgrids
708
709 IMPLICIT NONE
710
711 !=====Read ESTMinput.nml================================
712 namelist /estminput/ tsurfchoice, &
713 evolvetibld, &
714 ibldchmod, &
715 lbc_soil, &
716 theat_on, &
717 theat_off, &
719
720 OPEN (511, file=trim(fileinputpath)//'ESTMinput.nml', status='old')
721 READ (511, nml=estminput)
722 CLOSE (511)
723
724 !Convert specified temperatures to Kelvin
728
729 ALLOCATE (tair2_grids(numberofgrids))
731 ALLOCATE (lup_wall_grids(numberofgrids))
732 ALLOCATE (lup_roof_grids(numberofgrids))
733 ALLOCATE (tievolve_grids(numberofgrids))
734 ALLOCATE (t0_ibld_grids(numberofgrids))
736 ALLOCATE (t0_wall_grids(numberofgrids))
737 ALLOCATE (t0_roof_grids(numberofgrids))
738 ALLOCATE (tn_wall_grids(numberofgrids))
739 ALLOCATE (tn_roof_grids(numberofgrids))
740
character(len=150) fileinputpath

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.

Referenced by suews_program().

Here is the caller graph for this function:

◆ estm_translate()

subroutine estm_module::estm_translate ( integer  Gridiv)

Definition at line 1245 of file suews_phys_estm.f95.

1246 ! HCW 30 Jun 2016
1247
1248 USE defaultnotused, ONLY: nan
1249 USE physconstants, ONLY: c2k, sbconst
1250 USE estm_data
1251 USE allocatearray
1252 USE gis_data, ONLY: bldgh
1253 USE initial, ONLY: numberofgrids
1254
1255 IMPLICIT NONE
1256 INTEGER :: i
1257 !REAL(KIND(1d0)) :: CFLval
1258 !REAL(KIND(1d0)) :: t5min
1259 REAL(KIND(1D0)) :: W, WB
1260 !CHARACTER (len=20)::FileCodeX
1261 !CHARACTER (len=150):: FileFinalTemp
1262 !LOGICAL:: inittemps=.FALSE.
1263 INTEGER :: ESTMStart = 0
1264 INTEGER :: Gridiv
1265
1266 !Set initial values at the start of each run for each grid
1267 ! the initiliastion part is problematic:
1268 ! cannot be initialised under a multi-grid scenario
1269 IF (gridiv == 1) estmstart = estmstart + 1
1270 ! print *, 'gridiv in ESTM',gridiv,'ESTMStart',ESTMStart
1271 IF (estmstart == 1) THEN
1272
1273 ! write(*,*) ' ESTMStart: ',ESTMStart, 'initialising ESTM for grid no. ', Gridiv
1274
1275 tfloor = 20.0 ! This is used only when radforce =T !TODO: should be put in the namelist
1276 tfloor = tfloor + c2k
1277
1278 ! Initial values
1279 tievolve = 20.0 + c2k
1280 shc_air = 1230.0
1281 minshc_airbld = 1300
1282
1283 ! ---- Internal view factors ----
1284 !constant now but should be calculated in the future
1285 ivf_iw = 0.100000
1286 ivf_ir = 0.000000
1287 ivf_ii = 0.900000
1288 ivf_if = 0.000000
1289 ivf_ww = 0.050000
1290 ivf_wr = 0.000000
1291 ivf_wi = 0.950000
1292 ivf_wf = 0.000000
1293 ivf_rw = 0.050000
1294 ivf_ri = 0.950000
1295 ivf_rf = 0.000000
1296 ivf_fw = 0.050000
1297 ivf_fr = 0.000000
1298 ivf_fi = 0.950000
1299
1300 tair24hr = c2k
1301
1302 !Ts5mindata(1,ncolsESTMdata) = -999
1303 ! !Fill Ts5mindata for current grid and met block - this is done in SUEWS_translate
1305
1306 ! ---- Initialization of variables and parameters for first row of run for each grid ----
1307 ! N layers are calculated in SUEWS_translate
1308 IF (.NOT. ALLOCATED(tibld)) THEN
1309 ! print *, "Nibld", Nibld
1310 ! print *, "Nwall", Nwall
1311 ! print *, "Nroof", Nroof
1312 ! print *, "Nground", Nground
1319 END IF
1320
1321 ! Transfer variables from Ts5mindata to variable names
1322 ! N.B. column numbers here for the following file format - need to change if input columns change!
1323 ! dectime iy id it imin Tiair Tsurf Troof Troad Twall Twall_n Twall_e Twall_s Twall_w
1324 ! 1 2 3 4 5 6 7 8 9 10 11 12 13 !new
1325 !
1326 ! Calculate temperature of each layer in Kelvin
1327 ! QUESTION: what if (Nground/Nwall/Nroof-1)==0? TS 21 Oct 2017
1328 DO i = 1, ndepth_ground
1329 tground(i) = (lbc_soil - ts5mindata(1, cts_troad))*(i - 1)/(ndepth_ground - 1) + ts5mindata(1, cts_troad) + c2k
1330 END DO
1331 DO i = 1, ndepth_wall
1332 twall(i) = &
1333 (ts5mindata(1, cts_tiair) - ts5mindata(1, cts_twall))*(i - 1)/(ndepth_wall - 1) + ts5mindata(1, cts_twall) &
1334 + c2k
1335 END DO
1336 DO i = 1, ndepth_roof
1337 troof(i) = &
1338 (ts5mindata(1, cts_tiair) - ts5mindata(1, cts_troof))*(i - 1)/(ndepth_roof - 1) + ts5mindata(1, cts_troof) &
1339 + c2k
1340 END DO
1342
1343 END IF !End of loop run only at start (for each grid)
1344
1345 ! ---- Parameters related to land surface characteristics ----
1346 ! QUESTION: Would Zref=z be more appropriate?
1347 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)
1348
1349 svf_ground = 1.0
1350 svf_roof = 1.0
1351
1352 ! ==== roof (i.e. Bldgs)
1353 !froof=sfr_surf(BldgSurf) ! Moved to SUEWS_translate HCW 16 Jun 2016
1356
1357 ! ==== vegetation (i.e. EveTr, DecTr, Grass)
1358 !fveg=sfr_surf(ConifSurf)+sfr_surf(DecidSurf)+sfr_surf(GrassSurf) ! Moved to SUEWS_translate HCW 16 Jun 2016
1359 IF (fveg /= 0) THEN
1362 ELSE ! check fveg==0 scenario to avoid division-by-zero error, TS 21 Oct 2017
1365 END IF
1366
1367 ! ==== ground (i.e. Paved, EveTr, DecTr, Grass, BSoil, Water - all except Bldgs)
1368 !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
1369 IF (fground /= 0) THEN
1376 ELSE ! check fground==0 scenario to avoid division-by-zero error, TS 21 Jul 2016
1379 END IF
1380
1381 IF (froof < 1.0) THEN
1382 hw = fwall/(2.0*(1.0 - froof))
1383 ELSE
1384 ! HW=0 !HCW if only roof, no ground
1385 hw = 0.00001 ! to avoid zero-HW scenario TS 21 Oct 2017
1386
1387 END IF
1388 hw = max(0.00001, hw) ! to avoid zero-HW scenario TS 27 Oct 2017
1389
1390 IF (fground == 1.0) THEN !!FO!! if only ground, i.e. no houses
1391 w = 1
1392 wb = 0
1393 svf_ground = 1.
1394 zvf_wall = 0.
1395 svf_wall = 0.
1396 svf_roof = 1.
1397 zvf_ground = 0.
1398 xvf_wall = 0.
1399 rvf_canyon = 1.
1400 rvf_ground = 1.-fveg
1401 rvf_roof = 0
1402 rvf_wall = 0
1403 rvf_veg = fveg
1404 ELSE IF (fground == 0.0) THEN !check fground==0 (or HW==0) scenario to avoid division-by-zero error, TS 21 Jul 2016
1405 ! the following values are calculated given HW=0
1406 w = 0
1407 wb = 1
1408 zvf_wall = 0 !COS(ATAN(2/HW)) when HW=0 !!FO!! wall view factor for wall
1409 hw = 0
1410 svf_ground = max(cos(atan(2*hw)), 0.00001) !!FO!! sky view factor for ground ! to avoid zero-division scenario TS 21 Oct 2017
1411 svf_wall = (1 - zvf_wall)/2 !!FO!! sky view factor for wall
1412 zvf_ground = 1 - svf_ground !!FO!! wall view factor for ground
1413 xvf_wall = svf_wall !!FO!! ground view factor
1414 ! RVF_CANYON=COS(ATAN(2*ZREF/W))
1415 ! RVF_ROOF=1-RVF_CANYON
1416 ! RVF_WALL=(COS(ATAN(2*(ZREF-BldgH)/W))-RVF_CANYON)*RVF_CANYON
1417 ! RVF_ground=RVF_CANYON-RVF_WALL
1420 rvf_roof = froof
1422 ELSE
1423 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
1424 wb = w*sqrt(froof/fground)
1425 svf_ground = cos(atan(2*hw)) !!FO!! sky view factor for ground
1426 zvf_wall = cos(atan(2/hw)) !!FO!! wall view factor for wall
1427 svf_wall = (1 - zvf_wall)/2 !!FO!! sky view factor for wall
1428 zvf_ground = 1 - svf_ground !!FO!! wall view factor for ground
1429 xvf_wall = svf_wall !!FO!! ground view factor
1430 ! RVF_CANYON=COS(ATAN(2*ZREF/W))
1431 ! RVF_ROOF=1-RVF_CANYON
1432 ! RVF_WALL=(COS(ATAN(2*(ZREF-BldgH)/W))-RVF_CANYON)*RVF_CANYON
1433 ! RVF_ground=RVF_CANYON-RVF_WALL
1436 rvf_roof = froof
1438 END IF
1439
1441
1442 sumalb = 0.; nalb = 0
1443 sumemis = 0.; nemis = 0
1444
1445 !set emissivity for ceiling, wall and floor inside of buildings
1447
1448 !internal elements
1449 IF (nroom == 0) THEN
1450 fibld = (floor(bldgh/3.1 - 0.5) - 1)*froof
1451 ELSE
1452 fibld = (2.-2./nroom)*fwall + (floor(bldgh/3.1 - 0.5) - 1)*froof
1453 END IF
1454
1455 IF (fibld == 0) fibld = 0.00001 !this just ensures a solution to radiation
1457 IF (finternal == 0) finternal = 0.00001 ! to avoid zero-devision error TS 21 Oct 2017
1458 fair = zref - bldgh*froof
1459 IF (fair == 0) fair = 0.00001 ! to avoid zero-devision error TS 21 Oct 2017
1460 !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
1461 !ivf_ww=1.-ivf_wi-ivf_wr-ivf_wf
1462 !ivf_rw=1.-ivf_ri-ivf_rf;
1463 !ivf_fr=ivf_rf;
1464
1465 IF ((ivf_ii + ivf_iw + ivf_ir + ivf_if > 1.0001) .OR. &
1466 (ivf_wi + ivf_ww + ivf_wr + ivf_wf > 1.0001) .OR. &
1467 (ivf_ri + ivf_rw + ivf_rf > 1.0001) .OR. &
1468 (ivf_fi + ivf_fw + ivf_fr > 1.0001) .OR. &
1469 (ivf_ii + ivf_iw + ivf_ir + ivf_if < 0.9999) .OR. &
1470 (ivf_wi + ivf_ww + ivf_wr + ivf_wf < 0.9999) .OR. &
1471 (ivf_ri + ivf_rw + ivf_rf < 0.9999) .OR. &
1472 (ivf_fi + ivf_fw + ivf_fr < 0.9999)) THEN
1473 print *, "At least one internal view factor <> 1. Check ivf in ESTMinput.nml"
1474 END IF
1475
1476 !=======Initial setting==============================================
1477 !! Rewritten by HCW 15 Jun 2016 to use existing SUEWS error handling
1478 !IF(inittemps) THEN
1479 ! write(*,*) 'inittemps:',inittemps
1480 ! FileFinalTemp=TRIM(FileOutputPath)//TRIM(FileCodeX)//'_ESTM_finaltemp.txt'
1481 ! OPEN(99,file=TRIM(FileFinalTemp),status='old',err=316) ! Program stopped if error opening file
1482 ! READ(99,*) Twall,Troof,Tground,Tibld ! Twall, Troof, Tground & Tibld get new values
1483 ! CLOSE(99)
1484 !ENDIF
1485 !
1486 !!IF (inittemps) THEN !!FO!! inittemps=.true. set in nml file
1487 !! OPEN(99,file='outputfiles/finaltemp.txt',status='old',iostat=ios) !!FO!! has to exist
1488 !!
1489 !! IF (ios/=0) CALL error('outputfiles/finaltemp.txt',ios,1) !!FO!! calls mod_error.f95, writes that the opening failed and stops prg
1490 !! IF (ios/=0) THEN
1491 !! Twall = (/273., 285., 291./)
1492 !! Troof = (/273., 285., 291./)
1493 !! Tground = (/273., 275., 280., 290./)
1494 !! Tibld = (/293., 293., 293./)
1495 !! ELSE
1496 !! READ(99,*) Twall,Troof,Tground,Tibld !!FO!! if finaltemp.txt exists Twall[3], Troof[3], Tground[4] & Tibld[3] get new values
1497 !! CLOSE(99)
1498 !! ENDIF
1499 !!ENDIF
1500
1501 !where (isnan(Twall))
1502 ! Twall = 273
1503 !endwhere
1504 !where (isnan(Troof))
1505 ! Troof = 273
1506 !endwhere
1507 !where (isnan(Tground))
1508 ! Tground = 281
1509 !endwhere
1510 !where (isnan(Tibld))
1511 ! Tibld = 293
1512 !endwhere
1513
1514 IF (estmstart == 1) THEN
1515 DO i = 1, 4
1516 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
1517 END DO
1518
1519 !initialize surface temperatures
1520 t0_ground = tground(1)
1521 t0_wall = twall(1)
1522 t0_roof = troof(1)
1523 t0_ibld = tibld(1)
1526
1527 !initialize outgoing longwave !QUESTION: HCW - Are these calculations compatible with those in LUMPS_NARP?
1531
1532 ! PRINT*,"W,WB= ",W,WB
1533 ! PRINT*,'SVF_ground ','SVF_WALL ','zvf_WALL ','HW '
1534 ! PRINT*,SVF_ground,SVF_WALL,zvf_WALL,HW
1535 ! PRINT*,'RVF_ground ','RVF_WALL ','RVF_ROOF ','RVF_VEG'
1536 ! PRINT*,RVF_ground,RVF_WALL,RVF_ROOF,RVF_VEG
1537 ! print*,'Alb_avg (VF)=',alb_avg
1538 ! print*,'z0m, Zd', z0m, ZD
1539
1540 END IF
1541
1542 first = .true.
1543
1544 !======Courant-Friedrichs-Lewy condition=================================
1545 !This is comment out by S.O. for now
1546 ! NB: should be recovered for calculation stability, FO and TS, 11 Oct 2017
1547 ! CFLval = minval(0.5*zibld*zibld*ribld/kibld) !!FO!! z*z*r/k => unit [s]
1548 ! if (Tstep>CFLval) then !CFL condition !!FO!! CFL condition: Courant�Friedrichs�Lewy condition is a necessary condition for convergence while solving
1549 ! 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)
1550 ! CFLfail=.TRUE.
1551 ! endif
1552 ! CFLval = minval(0.5*zroof*zroof*rroof/kroof)
1553 ! if (Tstep>CFLval) then !CFL condition
1554 ! write(*,*) "ROOF: CFL condition: Tstep=",Tstep,">",CFLval
1555 ! CFLfail=.TRUE.
1556 ! endif
1557 ! CFLval = minval(0.5*zwall*zwall*rwall/kwall)
1558 ! if (Tstep>CFLval) then !CFL condition
1559 ! write(*,*) "WALL: CFL condition: Tstep=",Tstep,">",CFLval
1560 ! CFLfail=.TRUE.
1561 ! endif
1562 ! CFLval = minval(0.5*zground*zground*rground/kground)
1563 ! if (Tstep>CFLval) then !CFL condition
1564 ! write(*,*) "ground: CFL condition: Tstep=",Tstep,">",CFLval
1565 ! CFLfail=.TRUE.
1566 ! endif
1567 ! if (CFLfail) then
1568 ! write(*,*) "Increase dX or decrease maxtimestep. Hit any key to continue"
1569 ! read (*,*)
1570 ! endif
1571
1572 ! 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
1573
1574 RETURN
1575
1576 ! 315 CALL errorHint(11,TRIM(fileESTMTs),notUsed,notUsed,NotUsedI)
1577 ! 316 CALL errorHint(11,TRIM(fileFinalTemp),notUsed,notUsed,NotUsedI)
1578
real(kind(1d0)), dimension(:, :), allocatable ts5mindata
integer, parameter ncolsestmdata
integer, parameter cts_twall
integer, parameter cts_troad
real(kind(1d0)), dimension(nsurf) emis
real(kind(1d0)), dimension(:), allocatable tair24hr
real(kind(1d0)), dimension(nsurf) alb
integer, parameter cts_tiair
real(kind(1d0)), dimension(:, :, :), allocatable estmforcingdata
integer, parameter cts_troof
real(kind(1d0)) nan
real(kind(1d0)) ivf_fr
real(kind(1d0)) nroom
real(kind(1d0)) zref
real(kind(1d0)) ivf_fw
real(kind(1d0)) ivf_fi
real(kind(1d0)) fveg
real(kind(1d0)) rvf_canyon
real(kind(1d0)) bldgh

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.

Referenced by suews_translate().

Here is the caller graph for this function:

◆ load_gridlayout()

subroutine estm_module::load_gridlayout ( integer, intent(in)  gridIV,
logical, intent(in)  MultipleLayoutFiles,
integer, intent(in)  diagnose 
)

Definition at line 743 of file suews_phys_estm.f95.

744 USE allocatearray, ONLY: &
745 ndepth, &
780 nsurf, &
785 nspec
786 USE data_in, ONLY: fileinputpath, filecode
787 USE strings, ONLY: writenum
788 IMPLICIT NONE
789 INTEGER, INTENT(IN) :: gridIV
790 INTEGER, INTENT(IN) :: diagnose
791 LOGICAL, INTENT(IN) :: MultipleLayoutFiles
792 CHARACTER(len=100) :: FileLayout, str_gridIV
793 INTEGER :: istat, iunit, ERROR_UNIT, i
794 INTEGER :: igroup, ilayer, idepth
795 CHARACTER(len=1000) :: line
796 REAL(KIND(1D0)) :: k, dz
797
798 namelist &
799 ! basic dimension info
800 /dim/ &
801 nlayer, &
802 /geom/ &
803 height, &
805 veg_frac, &
807 veg_scale, &
808 ! read in roof part
809 /roof/ &
810 sfr_roof, &
812 tin_roof, &
813 alb_roof, &
814 emis_roof, &
815 state_roof, &
820 dz_roof, &
821 k_roof, &
822 cp_roof &
823 ! read in wall part
824 /wall/ &
825 sfr_wall, &
827 tin_wall, &
828 alb_wall, &
829 emis_wall, &
830 state_wall, &
835 dz_wall, &
836 k_wall, &
837 cp_wall, &
838 /surf/ &
839 tin_surf, &
840 dz_surf, &
841 k_surf, &
842 cp_surf
843
844 IF (multiplelayoutfiles) THEN
845 CALL writenum(gridiv, str_gridiv, 'i4')
846 filelayout = 'GridLayout'//trim(filecode)//trim(str_gridiv)//'.nml'
847 ELSE
848 filelayout = 'GridLayout'//trim(filecode)//'.nml'
849 END IF
850
851 IF (diagnose == 1) print *, 'Reading layout file: ', trim(fileinputpath)//filelayout
852 iunit = 212
853 OPEN (iunit, file=trim(fileinputpath)//trim(filelayout), status='old')
854 IF (diagnose == 1) print *, 'Read dim info of GridLayout'
855 READ (iunit, nml=dim, iostat=istat)
856 IF (diagnose == 1) print *, 'Number of layers to read: ', nlayer
857
858 IF (istat /= 0) THEN
859 backspace(iunit)
860 READ (iunit, fmt='(A)') line
861 error_unit = 119
862 WRITE (error_unit, '(A)') &
863 'Invalid line in namelist: '//trim(line)
864 END IF
865 CLOSE (iunit)
866
867 ALLOCATE (height(nlayer + 1))
868 ALLOCATE (building_frac(nlayer))
869 ALLOCATE (veg_frac(nlayer))
870 ALLOCATE (building_scale(nlayer))
871 ALLOCATE (veg_scale(nlayer))
872
873 ALLOCATE (sfr_roof(nlayer))
874 ALLOCATE (dz_roof(nlayer, ndepth))
875 ALLOCATE (k_roof(nlayer, ndepth))
876 ALLOCATE (cp_roof(nlayer, ndepth))
877 ALLOCATE (tin_roof(nlayer))
878 ALLOCATE (alb_roof(nlayer))
879 ALLOCATE (emis_roof(nlayer))
880 ALLOCATE (state_roof(nlayer))
881 ALLOCATE (statelimit_roof(nlayer))
882 ALLOCATE (wetthresh_roof(nlayer))
883 ALLOCATE (soilstore_roof(nlayer))
884 ALLOCATE (soilstorecap_roof(nlayer))
886
887 ALLOCATE (sfr_wall(nlayer))
888 ALLOCATE (dz_wall(nlayer, ndepth))
889 ALLOCATE (k_wall(nlayer, ndepth))
890 ALLOCATE (cp_wall(nlayer, ndepth))
891 ALLOCATE (tin_wall(nlayer))
892 ALLOCATE (alb_wall(nlayer))
893 ALLOCATE (emis_wall(nlayer))
894 ALLOCATE (state_wall(nlayer))
895 ALLOCATE (statelimit_wall(nlayer))
896 ALLOCATE (wetthresh_wall(nlayer))
897 ALLOCATE (soilstore_wall(nlayer))
898 ALLOCATE (soilstorecap_wall(nlayer))
899 ALLOCATE (wall_specular_frac(nspec, nlayer))
900
901 ALLOCATE (dz_surf(nsurf, ndepth))
902 ALLOCATE (k_surf(nsurf, ndepth))
903 ALLOCATE (cp_surf(nsurf, ndepth))
904 ALLOCATE (tin_surf(nsurf))
905
906 OPEN (iunit, file=trim(fileinputpath)//trim(filelayout), status='old')
907 IF (diagnose == 1) print *, 'Read geometry part of GridLayout'
908 READ (iunit, nml=geom, iostat=istat)
909 IF (diagnose == 1) print *, 'height', height
910 IF (diagnose == 1) print *, 'building_frac', building_frac
911 IF (diagnose == 1) print *, 'veg_frac', veg_frac
912 IF (diagnose == 1) print *, 'building_scale', building_scale
913 IF (diagnose == 1) print *, 'veg_scale', veg_scale
914
915 IF (diagnose == 1) print *, 'Read roof part of GridLayout'
916 READ (iunit, nml=roof, iostat=istat)
917 IF (diagnose == 1) print *, 'sfr_roof', sfr_roof
918 IF (diagnose == 1) print *, 'dz_roof', dz_roof
919 IF (diagnose == 1) print *, 'k_roof', k_roof
920 IF (diagnose == 1) print *, 'cp_roof', cp_roof
921 IF (diagnose == 1) print *, 'tin_roof', tin_roof
922 IF (diagnose == 1) print *, 'alb_roof', alb_roof
923 IF (diagnose == 1) print *, 'emis_roof', emis_roof
924 IF (diagnose == 1) print *, 'state_roof', state_roof
925 IF (diagnose == 1) print *, 'statelimit_roof', statelimit_roof
926 IF (diagnose == 1) print *, 'wetthresh_roof', wetthresh_roof
927 IF (diagnose == 1) print *, 'soilstore_roof', soilstore_roof
928 IF (diagnose == 1) print *, 'soilstorecap_roof', soilstorecap_roof
929
930 IF (diagnose == 1) print *, 'Read wall part of GridLayout'
931 READ (iunit, nml=wall, iostat=istat)
932 IF (diagnose == 1) print *, 'sfr_wall', sfr_wall
933 IF (diagnose == 1) print *, 'dz_wall', dz_wall
934 IF (diagnose == 1) print *, 'k_wall', k_wall
935 IF (diagnose == 1) print *, 'cp_wall', cp_wall
936 IF (diagnose == 1) print *, 'tin_wall', tin_wall
937 IF (diagnose == 1) print *, 'alb_wall', alb_wall
938 IF (diagnose == 1) print *, 'emis_wall', emis_wall
939 IF (diagnose == 1) print *, 'state_wall', state_wall
940 IF (diagnose == 1) print *, 'statelimit_wall', statelimit_wall
941 IF (diagnose == 1) print *, 'wetthresh_wall', wetthresh_wall
942 IF (diagnose == 1) print *, 'soilstore_wall', soilstore_wall
943 IF (diagnose == 1) print *, 'soilstorecap_wall', soilstorecap_wall
944
945 ! PRINT *, 'Read surf part of GridLayout'
946 READ (iunit, nml=surf, iostat=istat)
947 IF (istat /= 0) THEN
948 backspace(iunit)
949 READ (iunit, fmt='(A)') line
950 error_unit = 119
951 WRITE (error_unit, '(A)') &
952 'Invalid line in namelist: '//trim(line)
953 END IF
954 CLOSE (iunit)
955
956 IF (diagnose == 1) print *, 'Read GridLayout'
957 IF (diagnose == 1) print *, 'sfr_roof in load_GridLayout', sfr_roof
958 IF (diagnose == 1) print *, 'sfr_wall in load_GridLayout', sfr_wall
959 IF (diagnose == 1) print *, 'tin_surf in load_GridLayout', tin_surf
960 IF (diagnose == 1) print *, 'dz_surf(1,:) in load_GridLayout', dz_surf(1, :)
961 IF (diagnose == 1) print *, 'k_surf(1,:) in load_GridLayout', k_surf(1, :)
962 IF (diagnose == 1) print *, 'dz_roof in load_GridLayout', dz_roof(1:nlayer, :ndepth)
963 ! PRINT *, 'dz_roof in load_GridLayout', dz_roof(1, :ndepth)
964
965 ! assign values to the grid container
966 nlayer_grids(gridiv) = nlayer
967 height_grids(gridiv, 1:nlayer + 1) = height(1:nlayer + 1)
969 veg_frac_grids(gridiv, 1:nlayer) = veg_frac(1:nlayer)
971 veg_scale_grids(gridiv, 1:nlayer) = veg_scale(1:nlayer)
972
973 ! roof
974 sfr_roof_grids(gridiv, 1:nlayer) = sfr_roof
975 tin_roof_grids(gridiv, 1:nlayer) = tin_roof
976 DO i = 1, nlayer
977 temp_roof_grids(gridiv, i, :) = tin_roof(i)
978 END DO
979 dz_roof_grids(gridiv, 1:nlayer, 1:ndepth) = dz_roof(1:nlayer, 1:ndepth)
980 k_roof_grids(gridiv, 1:nlayer, 1:ndepth) = k_roof(1:nlayer, 1:ndepth)
981 cp_roof_grids(gridiv, 1:nlayer, 1:ndepth) = cp_roof(1:nlayer, 1:ndepth)
982 alb_roof_grids(gridiv, 1:nlayer) = alb_roof(1:nlayer)
983 emis_roof_grids(gridiv, 1:nlayer) = emis_roof(1:nlayer)
990 ! PRINT *, 'dz_roof_grids in loading', dz_roof_grids(Gridiv, 1:nroof, 1:ndepth)
991 ! PRINT *, 'dz_roof_grids(1,1) in loading', dz_roof_grids(Gridiv, 1, 1:ndepth)
992 ! tsfc_roof_grids(gridIV) = tsfc_roof
993 ! tin_roof_grids(gridIV) = temp_in_roof
994
995 ! wall
996 sfr_wall_grids(gridiv, 1:nlayer) = sfr_wall
997 tin_wall_grids(gridiv, 1:nlayer) = tin_wall
998 DO i = 1, nlayer
999 temp_wall_grids(gridiv, i, :) = tin_wall(i)
1000 END DO
1001 dz_wall_grids(gridiv, 1:nlayer, 1:ndepth) = dz_wall(1:nlayer, 1:ndepth)
1002 k_wall_grids(gridiv, 1:nlayer, 1:ndepth) = k_wall(1:nlayer, 1:ndepth)
1003 cp_wall_grids(gridiv, 1:nlayer, 1:ndepth) = cp_wall(1:nlayer, 1:ndepth)
1004 alb_wall_grids(gridiv, 1:nlayer) = alb_wall(1:nlayer)
1005 emis_wall_grids(gridiv, 1:nlayer) = emis_wall(1:nlayer)
1012 ! PRINT *, 'dz_wall_grids in loading', dz_wall_grids(Gridiv, 1:nwall, 1:ndepth)
1013 ! PRINT *, 'dz_wall_grids(1,1) in loading', dz_wall_grids(Gridiv, 1, 1:ndepth)
1014 ! tsfc_wall_grids(gridIV) = tsfc_wall
1015 ! tin_wall_grids(gridIV) = temp_in_wall
1016
1017 ! generic surfaces
1018 tin_surf_grids(gridiv, 1:nsurf) = tin_surf
1019 ! DO i = 1, nsurf
1020 ! tin_surf_grids(gridIV, i, :) = tin_surf(i)
1021 ! END DO
1022 dz_surf_grids(gridiv, 1:nsurf, 1:ndepth) = dz_surf(1:nsurf, 1:ndepth)
1023 k_surf_grids(gridiv, 1:nsurf, 1:ndepth) = k_surf(1:nsurf, 1:ndepth)
1024 cp_surf_grids(gridiv, 1:nsurf, 1:ndepth) = cp_surf(1:nsurf, 1:ndepth)
1025
1026 ! ! check if CFL condition is met
1027 ! DO igroup = 1, 3
1028 ! DO ilayer = 1, merge(nlayer,7,igroup<3)
1029 ! idepth = 1
1030 ! ! DO idepth = 1, ndepth
1031 ! IF (igroup == 1) THEN
1032 ! k = k_roof_grids(gridIV, ilayer, idepth)
1033 ! dz = dz_roof_grids(gridIV, ilayer, idepth)
1034 ! str_igroup='roof'
1035 ! ELSE IF (igroup == 2) THEN
1036 ! k = k_wall_grids(gridIV, ilayer, idepth)
1037 ! dz = dz_wall_grids(gridIV, ilayer, idepth)
1038 ! str_igroup='wall'
1039 ! ELSE IF (igroup == 3) THEN
1040 ! k = k_surf_grids(gridIV, ilayer, idepth)
1041 ! dz = dz_surf_grids(gridIV, ilayer, idepth)
1042 ! str_igroup='surf'
1043 ! END IF
1044
1045 ! IF (dz/k > 0.02) THEN
1046 ! PRINT *, 'CFL condition, dz/k < 0.02, is not met for top layer of:'
1047 ! PRINT *, 'grid', gridIV
1048
1049 ! PRINT *, 'group: ', str_igroup
1050 ! PRINT *, 'layer', ilayer
1051 ! PRINT *, 'dz/k=', dz/k
1052 ! PRINT *, 'dz=', dz
1053 ! PRINT *, 'k=', k
1054 ! PRINT *, ''
1055 ! WRITE (str_gridIV, '(I5.5)') gridIV
1056 ! ! WRITE (str_igroup, '(I5.5)') igroup
1057 ! WRITE (str_idepth, '(I5.5)') idepth
1058 ! WRITE (str_ilayer, '(I5.5)') ilayer
1059
1060 ! CALL ErrorHint(77, 'CFL condition is not met for '// &
1061 ! 'grid: '//TRIM(str_gridIV)// &
1062 ! ', group: '//TRIM(str_igroup)// &
1063 ! ', layer: '//TRIM(str_ilayer)// &
1064 ! ', depth: '//TRIM(str_idepth) &
1065 ! , dz/k, idepth*1.D0, ilayer)
1066 ! END IF
1067 ! ! END DO
1068 ! END DO
1069
1070 ! END DO
1071
1072 DEALLOCATE (height)
1073 DEALLOCATE (building_frac)
1074 DEALLOCATE (veg_frac)
1075 DEALLOCATE (building_scale)
1076 DEALLOCATE (veg_scale)
1077
1078 DEALLOCATE (sfr_roof)
1079 DEALLOCATE (alb_roof)
1080 DEALLOCATE (emis_roof)
1081 DEALLOCATE (state_roof)
1082 DEALLOCATE (statelimit_roof)
1083 DEALLOCATE (wetthresh_roof)
1084 DEALLOCATE (soilstore_roof)
1085 DEALLOCATE (soilstorecap_roof)
1086 DEALLOCATE (dz_roof)
1087 DEALLOCATE (k_roof)
1088 DEALLOCATE (cp_roof)
1089 DEALLOCATE (tin_roof)
1090 DEALLOCATE (roof_albedo_dir_mult_fact)
1091
1092 DEALLOCATE (sfr_wall)
1093 DEALLOCATE (alb_wall)
1094 DEALLOCATE (emis_wall)
1095 DEALLOCATE (state_wall)
1096 DEALLOCATE (statelimit_wall)
1097 DEALLOCATE (wetthresh_wall)
1098 DEALLOCATE (soilstore_wall)
1099 DEALLOCATE (soilstorecap_wall)
1100 DEALLOCATE (dz_wall)
1101 DEALLOCATE (k_wall)
1102 DEALLOCATE (cp_wall)
1103 DEALLOCATE (tin_wall)
1104 DEALLOCATE (wall_specular_frac)
1105
1106 DEALLOCATE (tin_surf)
1107 DEALLOCATE (dz_surf)
1108 DEALLOCATE (k_surf)
1109 DEALLOCATE (cp_surf)
1110
real(kind(1d0)), dimension(:), allocatable statelimit_wall
real(kind(1d0)), dimension(:), allocatable soilstorecap_wall
real(kind(1d0)), dimension(:), allocatable alb_roof
real(kind(1d0)), dimension(:), allocatable state_wall
real(kind(1d0)), dimension(:), allocatable wetthresh_wall
real(kind(1d0)), dimension(:), allocatable state_roof
real(kind(1d0)), dimension(:), allocatable building_scale
real(kind(1d0)), dimension(:), allocatable soilstore_roof
real(kind(1d0)), dimension(:), allocatable height
real(kind(1d0)), dimension(:), allocatable emis_roof
real(kind(1d0)), dimension(:), allocatable wetthresh_roof
real(kind(1d0)), dimension(:), allocatable soilstorecap_roof
real(kind(1d0)), dimension(:), allocatable veg_frac
real(kind(1d0)), dimension(:, :), allocatable roof_albedo_dir_mult_fact
real(kind(1d0)), dimension(:), allocatable statelimit_roof
real(kind(1d0)), dimension(:), allocatable alb_wall
real(kind(1d0)), dimension(:, :), allocatable wall_specular_frac
real(kind(1d0)), dimension(:), allocatable emis_wall
real(kind(1d0)), dimension(:), allocatable building_frac
real(kind(1d0)), dimension(:), allocatable veg_scale
real(kind(1d0)), dimension(:), allocatable soilstore_wall
character(len=20) filecode

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_ext_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 2571 of file suews_phys_estm.f95.

2572 IMPLICIT NONE
2573 REAL(KIND(1D0)), PARAMETER :: pNAN = 9999
2574 REAL(KIND(1D0)), PARAMETER :: NAN = -999
2575 REAL(KIND(1D0)), INTENT(in) :: x
2576 REAL(KIND(1D0)) :: xx
2577
2578 IF (abs(x) > pnan) THEN
2579 xx = nan
2580 ELSE
2581 xx = x
2582 END IF
2583

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 634 of file suews_phys_estm.f95.

637 USE sues_data, ONLY: tstep_real, tstep
638 USE defaultnotused, ONLY: notused, ios_out
640
641 IMPLICIT NONE
642
643 INTEGER, INTENT(in) :: lunit
644 INTEGER :: i, iyy !,RunNumber,NSHcounter
645 INTEGER :: iostat_var
646 REAL(KIND(1D0)), DIMENSION(ncolsESTMdata) :: ESTMArray
647 REAL(KIND(1D0)) :: imin_prev, ih_prev, iday_prev, tstep_estm !For checks on temporal resolution of estm data
648
649 ! initialise
650 imin_prev = 0
651 ih_prev = 0
652 iday_prev = 0
653
654 !---------------------------------------------------------------
655
656 !Open the file for reading and read the actual data
657 !write(*,*) FileESTMTs
658 OPEN (lunit, file=trim(fileestmts), status='old', err=315)
659 CALL skipheader(lunit, skipheadermet)
660
661 ! Skip to the right place in the ESTM file, depending on how many chunks have been read already
662 IF (skippedlines > 0) THEN
663 DO iyy = 1, skippedlines
664 READ (lunit, *)
665 END DO
666 END IF
667
668 ! Read in next chunk of ESTM data and fill ESTMForcingData array with data for every timestep
669 DO i = 1, readlinesmetdata
670 READ (lunit, *, iostat=iostat_var) estmarray
671 estmforcingdata(i, 1:ncolsestmdata, gridcounter) = estmarray
672 ! Check timestamp of met data file matches TSTEP specified in RunControl
673 IF (i == 1) THEN
674 imin_prev = estmarray(4)
675 ih_prev = estmarray(3)
676 iday_prev = estmarray(2)
677 ELSEIF (i == 2) THEN
678 tstep_estm = ((estmarray(4) + 60*estmarray(3)) - (imin_prev + 60*ih_prev))*60 !tstep in seconds
679 IF (tstep_estm /= tstep_real .AND. estmarray(2) == iday_prev) THEN
680 CALL errorhint(39, 'TSTEP in RunControl does not match TSTEP of ESTM data (DOY).', &
681 REAL(tstep, KIND(1D0)), tstep_estm, INT(ESTMArray(2)))
682 END IF
683 END IF
684 END DO
685
686 CLOSE (lunit)
687
688 RETURN
689
690315 CALL errorhint(11, trim(fileestmts), notused, notused, ios_out)
691
character(len=150) fileestmts
integer skipheadermet
real(kind(1d0)) notused
integer readlinesmetdata
integer gridcounter
integer skippedlines
real(kind(1d0)) tstep_real
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
subroutine skipheader(lfn, skip)

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.

Referenced by suews_program().

Here is the call graph for this function:
Here is the caller graph for this function: