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

Functions/Subroutines

subroutine suews_getestmdata (lunit)
 
subroutine estm_initials
 
subroutine estm_translate (Gridiv)
 
subroutine estm (Gridiv, tstep, avkdn, avu1, temp_c, zenith_deg, avrh, press_hpa, ldown, bldgh, Ts5mindata_ir, Tair_av, dataOutLineESTM, QS)
 
elemental real(kind(1d0)) function set_nan (x)
 

Function/Subroutine Documentation

◆ estm()

subroutine estm_module::estm ( integer, intent(in)  Gridiv,
integer, intent(in)  tstep,
real(kind(1d0)), intent(in)  avkdn,
real(kind(1d0)), intent(in)  avu1,
real(kind(1d0)), intent(in)  temp_c,
real(kind(1d0)), intent(in)  zenith_deg,
real(kind(1d0)), intent(in)  avrh,
real(kind(1d0)), intent(in)  press_hpa,
real(kind(1d0)), intent(in)  ldown,
real(kind(1d0)), intent(in)  bldgh,
real(kind(1d0)), dimension(ncolsestmdata), intent(in)  Ts5mindata_ir,
real(kind(1d0)), intent(in)  Tair_av,
real(kind(1d0)), dimension(27), intent(out)  dataOutLineESTM,
real(kind(1d0)), intent(out)  QS 
)

Definition at line 977 of file suews_phys_estm.f95.

References estm_data::alb_avg, estm_data::alb_ground, estm_data::alb_roof, estm_data::alb_veg, estm_data::alb_wall, 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_data::em_i, estm_data::em_ibld, estm_data::em_r, estm_data::em_roof, estm_data::em_veg, estm_data::em_w, estm_data::em_wall, 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::nemis, mod_solver::newtonpolynomial(), estm_data::nground, estm_data::nibld, estm_data::nroof, estm_data::nwall, estm_data::pcoeff, mathconstants::pi, estm_data::qs_4, estm_data::qsair, estm_data::qsground, estm_data::qsibld, estm_data::qsroof, estm_data::qswall, estm_data::rvf_ground, estm_data::rvf_roof, estm_data::rvf_veg, estm_data::rvf_wall, 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().

977  ! NB: HCW Questions:
978  ! - should TFloor be set in namelist instead of hard-coded here?
979  ! - zref used for radiation calculation and fair is set to 2*BldgH here. For compatibility with the rest of the
980  ! SUEWS model, should this be the (wind speed) measurement height z specified in RunControl.nml?
981  ! - In SUEWS_translate, fwall=AreaWall/SurfaceArea. Is this correct?
982  ! - If froof=1 (i.e. whole grid is building), is HW=0 correct?
983  ! - Then is an IF(Fground ==0) option needed?
984  ! - alb_wall=0.23 and em_wall=0.9 are set in LUMPS_module_constants. Shouldn't these be provided as input?
985  ! - Do the LUP calculations here need to be compatible with those in LUMPS_NARP?
986  ! - File opening rewritten using existing error handling in SUEWS - can delete mod_error module from SUEWS_ESTM_functions
987  ! - In SUEWS_ESTM_v2016, the first row is set to -999. This may be acceptable at the
988  ! start of the run but should be handled properly between blocks of met data? - need to check what's actually happening here.
989  ! - Many duplicate functions in SUEWS_ESTM_functions need changing to the existing SUEWS versions.
990  ! - Are the following correctly initialised? T0_ibld, T0_ground, T0_wall, T0_roof, TN_wall, TN_roof, Tground, Twall, Troof, Tibld
991  ! - What are Nalb, sumalb, Nemis and Sumemis for? Are they used correctly?
992  !
993 
994  !SUEWS_ESTM_v2016
995  ! revision:
996  ! HCW 14 Jun 2016
997  ! HCW 27 Jun 2016 Corrected iESTMcount bug - now increases for all grids together
998  ! HCW 30 Jun 2016 Major changes to handle grids and met blocks
999  ! TS 09 Oct 2017 Added explicit interface
1000 
1001  !Contains calculation for each time step
1002  !Calculate local scale heat storage from single building and surroundings observations
1003  !OFferle, May 2003
1004  !
1005  !MODIFICATION HISTORY
1006  ! 15 DECEMBER 2003
1007  ! (1) CHANGED AIR EXCHANGE RATE TO ALSO BE DEPENDENT ON OUTSIDE AIR TEMPERATURE
1008  ! (2) ADDED SPECIFIC HEAT CALCULATION FOR AIR
1009  ! (3) RH ADDED AS VAR(14) IN INPUT FILE
1010  ! 12 JANUARY 2004
1011  ! (1) ADDED ESTIMATED AVG NADIR LOOKING EXTERNAL SURFACE TEMPERATURE TO OUTPUT
1012  ! WEIGHTED BY SURFACE FRACTION AND SVF. SVF_ROOF=1, SVF_WALLS = 1-SVF_ground
1013  ! (2) ADDED NET RADIATION (RN) CALCULATIONS FOR ALL SURFACES AND AVG RN TO OUTPUT BASED ON
1014  ! RADIOMETER VIEW FROM ZREF
1015  ! 14 JANUARY 2004
1016  ! (1) MOVED GRID SPECIFIC PARAMETERS OUT OF NAMELIST, INTO PARAMETER FILE E.G. ALB,EM,F
1017  ! T0 MAKE CHANGES IN LOCATIONS EASIER.
1018  !
1019  ! 23 JANUARY 2004 - Version 2
1020  ! (1) PUT ALL TEMPERATURES INTO K
1021  ! (2) added interpolation routine to run on shorter timesteps.
1022  ! tested so that it doesn't change solution for forced temperatures.
1023  ! however in version 2 the energy balance at the surface isn't correctly solved
1024  !
1025  ! 25 JANUARY 2004 - Version 3
1026  ! (1) added solution to energy balance at surface
1027  ! (2) removed some extraneous code
1028  ! (3) added vegetation fractions for future development
1029  ! (4) some changes to input namelist.
1030  ! (5) need to add wind speed dependence for exchange coefficients
1031  ! (6) changed the way Rn_net is calculated. also radiometer view factor relationships
1032  ! (7) added a calculation for heat loss/gain to outside air (that going into building mass is storage)
1033  ! this is labelled QFBLD which it is in a sense.
1034  ! 6 FEBRUARY 2004
1035  ! (1) ADDED INTERNAL VIEW FACTOR FILE FOR INTERNAL GEOMETRY, INCLUDING FIXED FLOOR TEMPERATURE
1036  ! (2) added MeanU, MinWS to config, and U to ILOC.
1037  !
1038  ! 11 FEBRUARY 2004
1039  ! (1) added site lat, long, elevation to inputs
1040  ! (2) need zenith angle for wall direct radiation interception
1041  !
1042  ! 15 JUNE 2004
1043  ! (1) CORRECTED OUTPUT OF T0 FOR FORCED SURFACE TEMPERATURES
1044  ! (2) MADE SOME CHANGES TO RADIATION, COMPUTATION OF AVERAGE ALBEDO, ADDED AVERAGE EMISSIVITY
1045  ! (3) ADDED OUTPUT FILE FOR RADIATION COMPONENTS
1046  ! (4) CHANGED INPUT CONFIG SO CONVECTIVE EXCHANGE COEFFS ARE NOT USER SELECTABLE
1047  ! (5) MAY STILL BE PROBLEMS WITH WALL RADIATIVE EXCHANGE AND AVERAGE ALBEDO
1048  ! (6) CHANGED THE WAY HEAT STORAGE IS COMPUTED IN HEATCONDUCTION MODULE BUT THIS SHOULD NOT CHANGE RESULTS
1049  !
1050  ! 16 NOVEMBER 2004
1051  ! (1) CHANGED AIR EXCHANGE RATE TO BE BASES ON DAILY TEMPERATURE CHANGES.
1052  ! (2) WRITES OUT FINAL LAYER TEMPERATURES AND INCLUDES OPTION TO READ AT BEGINNING.
1053  !
1054  !DESCRIPTION: uses explicit time differencing to compute heat conduction through roof, walls,
1055  ! internal mass, and grounds (elements). Air heat storage is computed from average air temperature
1056  ! Boundary conditions are determined by measured surface temperature(s) or computed
1057  ! from energy balance at the surface.
1058  ! Internal air temperature can either be fixed or allowed to evolve in response
1059  ! to air mass exchanges and convective heating from internal surfaces.
1060  !
1061  !INPUT:
1062  !FORCING DATA: DTIME,KDOWN,LDOWN,TSURF,TAIR_OUT,TAIR_IN,TROOF,TWALL_AVG,TWALL_N,TWALL_E,TWALL_S,TWALL_W,Tground,RH,U
1063  !NAMELIST : HEATSTORAGE_vFO.NML
1064  ! &config
1065  ! ifile=FORCING DATA
1066  ! ofile=OUTPUT FILE
1067  ! pfile=HEAT STORAGE PARAMETER FILE
1068  ! Nibld = INTERNAL MASS LAYERS
1069  ! Nwall = EXTERNAL WALL LAYERS
1070  ! Nroof = ROOF LAYERS
1071  ! Nground = ground/SOIL LAYERS
1072  ! LBC_soil = LOWER BOUNDARY CONDITION FOR ground/SOIL
1073  ! iloc= INPUT COLUMNS IN DATA FILE
1074  ! evolveTibld= USE DIAGNOSTIC VALUE FOR INTERNAL BUILDING TEMPERATURE
1075  ! 0: don't use, use measured
1076  ! 1: TURN ON USE when temp goess ABOVE TINT_ON, off when temp is below TINT_OFF
1077  ! 2: always use diagnostic
1078  ! THEAT_ON= TEMPERATURE AT WHICH HEAT CONTROL IS TURNED ON
1079  ! THEAT_OFF= TEMPERATURE AT WHICH HEAT CONTROL IS TURNED OFF
1080  ! THEAT_FIX = Fixed internal temperature for climate control
1081  ! oneTsurf= USE SINGLE SURFACE TEMPERATURE TO DRIVE ALL LAYERS
1082  ! radforce= USE RADIATIVE ENERGY BALANCE TO DRIVE EXTERNAL TEMPERATURES
1083  ! maxtimestep=302, maximum time step in s for filling data gaps.
1084  ! Alt = STATION HEIGHT (m) FOR PRESSURE CALCULATION
1085  ! SPINUP = NUMBER OF LINES TO USE FOR SPINUP (REPEATS THESE LINES BUT ONLY OUTPUTS THE 2ND TIME)
1086  ! INITTEMP = if TRUE INITIALIZES TEMPERATURES TO THOSE IN FINALTEMP.TXT FILE
1087  ! CH_ibld = INTERNAL BUILDING CONVECTIVE EXCHANGE COEFFICIENT
1088  ! **** THESE SHOULD DEPEND ON WIND SPEED BUT CURRENTLY DO NOT ****
1089  ! CHAIR = CONVECTIVE EXCHANGE COEFFICIENT FOR ROOF
1090  ! chair_ground = ... FOR ground
1091  ! chair_wall = ... FOR WALL
1092  ! /
1093  ! ***************** PARAMETER FILE VARIABLES
1094  ! fveg = FRACTION OF ground SURFACE COVERED WITH VEG
1095  ! zveg = VEGETATION HEIGHT
1096  ! alb_veg = VEGETATION ALBEDO
1097  ! em_veg = VEGETATION EMISSIVITY
1098  ! ZREF = REFERENCE HEIGHT FOR FLUX CALCULATION
1099  ! BldgH = mean building height
1100  ! HW = CANYON ASPECT RATION
1101  ! f_X = FRACTION OF X WHERE X IS INTERNAL, WALL, ROOF, ground
1102  ! Alb_x = ALBEDO OF X
1103  ! em_ibld = EMISSIVITY OF X
1104  ! TX = INITIAL LAYER TEMPERATURES
1105  ! zX = LAYER THICKNESS
1106  ! kX = LAYER THERMAL CONDUCTIVITY
1107  ! ribld = LAYER VOLUMETRIC HEAT CAPACITY
1108  !
1109  !****************** INTERNAL VIEW FACTOR FILE
1110  !OUTPUT: fixed format text with single header, heatstorage for all elements, and temperatures
1111  ! for each element-layer.
1112  !===============================================================================
1113 
1114  USE meteo, ONLY: pi, heatcapacity_air
1115  USE mod_solver
1116  USE modsolarcalc !!FO!! :modsolarcalc.f95
1117  USE mathconstants !!FO!! :MathConstants_module.f95
1118  USE physconstants
1119  USE heatflux
1120  USE estm_data
1121 
1122  IMPLICIT NONE
1123  INTEGER, PARAMETER:: ncolsestmdata = 13
1124  ! INTEGER, PARAMETER:: ncolumnsDataOutESTM=32
1125  INTEGER, PARAMETER:: cts_tiair = 5
1126  INTEGER, PARAMETER:: cts_tsurf = 6
1127  INTEGER, PARAMETER:: cts_troof = 7
1128  INTEGER, PARAMETER:: cts_troad = 8
1129  INTEGER, PARAMETER:: cts_twall = 9
1130  INTEGER, PARAMETER:: cts_twall_n = 10
1131  INTEGER, PARAMETER:: cts_twall_e = 11
1132  INTEGER, PARAMETER:: cts_twall_s = 12
1133  INTEGER, PARAMETER:: cts_twall_w = 13
1134  REAL(KIND(1d0)), PARAMETER::nan = -999
1135 
1136  INTEGER, INTENT(in)::gridiv
1137  INTEGER, INTENT(in)::tstep
1138  ! INTEGER,INTENT(in)::iy !Year
1139  ! INTEGER,INTENT(in)::id !Day of year
1140  ! INTEGER,INTENT(in)::it !Hour
1141  ! INTEGER,INTENT(in)::imin !Minutes
1142 
1143  REAL(KIND(1d0)), INTENT(in)::avkdn
1144  REAL(KIND(1d0)), INTENT(in)::avu1
1145  REAL(KIND(1d0)), INTENT(in)::temp_c
1146  REAL(KIND(1d0)), INTENT(in)::zenith_deg
1147  REAL(KIND(1d0)), INTENT(in)::avrh
1148  REAL(KIND(1d0)), INTENT(in)::press_hpa
1149  REAL(KIND(1d0)), INTENT(in)::ldown
1150  REAL(KIND(1d0)), INTENT(in)::bldgh
1151  ! REAL(KIND(1d0)),INTENT(in):: dectime !Decimal time
1152  REAL(KIND(1d0)), DIMENSION(ncolsESTMdata), INTENT(in):: ts5mindata_ir !surface temperature input data
1153  REAL(KIND(1d0)), INTENT(in) :: tair_av ! mean air temperature of past 24hr
1154 
1155  REAL(KIND(1d0)), DIMENSION(27), INTENT(out):: dataoutlineestm
1156  !Output to SUEWS
1157  REAL(KIND(1d0)), INTENT(out)::qs
1158  !Input from SUEWS, corrected as Gridiv by TS 09 Jun 2016
1159 
1160  !Use only in this subroutine
1161  INTEGER::i, ii
1162  INTEGER:: tair2set = 0
1163  REAL(KIND(1d0))::airexhr, airexdt
1164  REAL(KIND(1d0)), DIMENSION(2)::bc
1165  REAL(KIND(1d0))::chair_ground, chair_wall
1166  REAL(KIND(1d0))::em_equiv
1167  REAL(KIND(1d0))::kdz
1168  REAL(KIND(1d0))::kup_estm, lup_net, kdn_estm
1169  REAL(KIND(1d0))::qhestm
1170  REAL(KIND(1d0))::qfbld !Anthropogenic heat from HVAC
1171  REAL(KIND(1d0))::shc_airbld
1172  REAL(KIND(1d0))::sw_hor, sw_vert
1173  REAL(KIND(1d0))::t0
1174  REAL(KIND(1d0))::tinternal, tsurf_all, troof_in, troad, twall_all, tw_n, tw_e, tw_s, tw_w
1175  REAL(KIND(1d0))::twallout(5), troofout(5), tibldout(5), tgroundout(5)
1176  REAL(KIND(1d0))::tadd, tveg
1177  REAL(KIND(1d0))::tairmix
1178  REAL(KIND(1d0))::rn
1179  REAL(KIND(1d0))::rs_roof, rl_roof, rn_roof
1180  REAL(KIND(1d0))::rs_wall, rl_wall, rn_wall
1181  REAL(KIND(1d0))::rs_ground, rl_ground, rn_ground
1182  REAL(KIND(1d0))::rs_ibld, rl_ibld
1183  REAL(KIND(1d0))::rs_iroof, rl_iroof
1184  REAL(KIND(1d0))::rs_iwall, rl_iwall
1185  REAL(KIND(1d0))::zenith_rad
1186  REAL(KIND(1d0))::dum(50)
1187  REAL(KIND(1d0))::bldghx ! local bldgHX=max(bldgH,0.001)
1188  REAL(KIND(1d0)), PARAMETER::wsmin = 0.1 ! Check why there is this condition. S.O.
1189  LOGICAL::radforce, groundradforce
1190 
1191  radforce = .false.
1192  groundradforce = .false. !Close the radiation scheme in original ESTM S.O.O.
1193 
1194  bldghx = max(bldgh, 0.001) ! this is to prevent zero building height
1195 
1196  ! Set -999s for first row
1197  ! dum=(/-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,&
1198  ! -999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,&
1199  ! -999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,&
1200  ! -999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,&
1201  ! -999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999./)
1202  dum = [(-999, i=1, 50)]
1203 
1204  !External bulk exchange coefficients - set these somewhere more sensible***
1205  chr = 0.005
1206  chair = chr
1207  chair_ground = chair
1208  chair_wall = chair
1209 
1210  !Get met data for use in ESTM subroutine
1211  kdn_estm = avkdn
1212  ws = avu1
1213  IF (ws < wsmin) ws = wsmin
1214  tair1 = temp_c + c2k
1215  ! Set initial value of Tair2 to air temp
1216  IF (gridiv == 1) tair2set = tair2set + 1
1217  IF (tair2set == 1) THEN
1218  tair2 = temp_c + c2k
1219  ELSE
1220  tair2 = tair2_grids(gridiv)
1221  ! Also get other variables for this grid
1222  tievolve = tievolve_grids(gridiv)
1223  lup_ground = lup_ground_grids(gridiv)
1224  lup_wall = lup_wall_grids(gridiv)
1225  lup_roof = lup_roof_grids(gridiv)
1226  t0_ibld = t0_ibld_grids(gridiv)
1227  t0_ground = t0_ground_grids(gridiv)
1228  t0_wall = t0_wall_grids(gridiv)
1229  t0_roof = t0_roof_grids(gridiv)
1230  tn_wall = tn_wall_grids(gridiv)
1231  tn_roof = tn_roof_grids(gridiv)
1232  tground(:) = tground_grids(:, gridiv)
1233  twall(:) = twall_grids(:, gridiv)
1234  troof(:) = troof_grids(:, gridiv)
1235  tibld(:) = tibld_grids(:, gridiv)
1236  tw_4 = tw_4_grids(:, :, gridiv)
1237 
1238  ENDIF
1239 
1240  ! Get Ts from Ts5min data array
1241  ! Tinternal = Ts5mindata(ir,cTs_Tiair)
1242  ! Tsurf_all = Ts5mindata(ir,cTs_Tsurf)
1243  ! Troof_in = Ts5mindata(ir,cTs_Troof)
1244  ! Troad = Ts5mindata(ir,cTs_Troad)
1245  ! Twall_all = Ts5mindata(ir,cTs_Twall)
1246  !
1247  ! Tw_n = Ts5mindata(ir,cTs_Twall_n)
1248  ! Tw_e = Ts5mindata(ir,cTs_Twall_e)
1249  ! Tw_s = Ts5mindata(ir,cTs_Twall_s)
1250  ! Tw_w = Ts5mindata(ir,cTs_Twall_w)
1251  tinternal = ts5mindata_ir(cts_tiair)
1252  tsurf_all = ts5mindata_ir(cts_tsurf)
1253  troof_in = ts5mindata_ir(cts_troof)
1254  troad = ts5mindata_ir(cts_troad)
1255  twall_all = ts5mindata_ir(cts_twall)
1256 
1257  tw_n = ts5mindata_ir(cts_twall_n)
1258  tw_e = ts5mindata_ir(cts_twall_e)
1259  tw_s = ts5mindata_ir(cts_twall_s)
1260  tw_w = ts5mindata_ir(cts_twall_w)
1261 
1262  ! if (any(isnan(Ts5mindata))) then !!FO!! can't use data when time gap is too big (or neg.) or data is NaN
1263  ! if (spindone) then !!FO!! writes a line of NaNs
1264  ! write(20,'(1F8.4,I6,100f10.1)') dectime,it,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,&
1265  ! -0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,&
1266  ! -0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,&
1267  ! -0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.,-0./0.
1268  ! endif
1269  ! return ! changed from cycle !!FO!! returns to beginning of do-loop, i.e. doesn't make any heat calculations
1270  ! endif
1271  !!FO!! this loop is used to run through the calculations (SPINUP number of times) to achieve better numerical stability
1272  !!FO!! when if condition is met the program starts over from the beginning of the input file and the calculations are performed
1273  !!FO!! ...once again, but this time saved to the output file
1274  !!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)
1275  ! IF (NLINESREAD>datalines.AND..NOT.SPINDONE) THEN !!FO!! at this stage spindone = .false.
1276  ! NLINESREAD=0
1277  ! SPINDONE=.TRUE.
1278  !PRINT*, "SPUNUP"
1279  ! return ! changed from cycle
1280  ! ENDIF
1281 
1282  !! Write first row of each met block as -999
1283  !IF (first) THEN !Set to true in ESTM_initials
1284  !! !Tair2=Temp_C+C2K !This is now set in SUEWS_translate for ir=0 only
1285  !! ! first=.FALSE.
1286  ! IF(Gridiv == NumberOfGrids) first=.FALSE. !Set to false only after all grids have run
1287  ! dataOutESTM(ir,1:32,Gridiv)=(/REAL(iy,KIND(1D0)),REAL(id,KIND(1D0)),&
1288  ! REAL(it,KIND(1D0)),REAL(imin,KIND(1D0)),dectime,(dum(ii),ii=1,27)/)
1289  ! RETURN
1290  !ENDIF
1291 
1292  ! What are these constants? - Need defining somewhere
1293  zenith_rad = zenith_deg/180*pi
1294  IF (zenith_rad > 0 .AND. zenith_rad < pi/2.-hw) THEN !ZENITH MUST BE HIGHER THAN BUILDINGS FOR DIRECT INTERCEPTION
1295  tanzenith = min(tan(zenith_rad), 5.67) !LIMITS TO ANGLES LESS THAN 80 EVEN FOR LOW HW
1296  tanzenith = tanzenith*kdn_estm/(1370*cos(zenith_rad)) !REDUCTION FACTOR FOR MAXIMUM
1297  ELSE
1298  tanzenith = 0.
1299  ENDIF
1300 
1301  shc_air = heatcapacity_air(tair1, avrh, press_hpa) ! Use SUEWS version
1302 
1303  !Evolution of building temperature from heat added by convection
1304  SELECT CASE (evolvetibld) !EvolveTiBld specifies which internal building temperature approach to use
1305  CASE (0)
1306  diagnoseti = .false.
1307  hvac = .false. !use data in file !!FO!! use measured indoor temperature (Tref in Lodz2002HS.txt)
1308  CASE (1) !!FO!! use of HVAC to counteract T changes
1309  diagnoseti = .true.
1310  IF (tievolve > theat_off) THEN !THEAT_OFF now converted to Kelvin in ESTM_initials - HCW 15 Jun 2016
1311  !IF (Tievolve>THEAT_OFF+C2K) THEN
1312  hvac = .false.
1313  ELSEIF (tievolve < theat_on) THEN !THEAT_OFF now converted to Kelvin in ESTM_initials - HCW 15 Jun 2016
1314  !ELSEIF (Tievolve<THEAT_ON+C2K) THEN
1315  hvac = .true.
1316  ENDIF
1317  CASE (2)
1318  diagnoseti = .true. !!FO!! convection between ibld and inside of external walls(?)
1319  END SELECT
1320 
1321  !ASSUME AIR MIXES IN PROPORTION TO # OF EXCHANGES
1322  IF (tair_av > 20.+c2k .AND. tievolve > 25.+c2k .AND. tair1 < tievolve .AND. .NOT. hvac) THEN
1323  airexhr = 2.0 !Windows or exterior doors on 3 sides (ASHRAE 1981 22.8)
1324  ELSEIF (tair_av < 17.+c2k .OR. hvac) THEN
1325  airexhr = 0.5 !No window or exterior doors, storm sash or weathertripped (ASHRAE 1981 22.8)
1326  ELSE
1327  airexhr = 1.0
1328  ENDIF
1329 
1330  airexdt = airexhr*(tstep/3600.0)
1331  shc_airbld = max(heatcapacity_air(tievolve, avrh, press_hpa), 0.00001) ! to avoid zero-division scenario TS 21 Oct 2017
1332  IF (shc_airbld < minshc_airbld) minshc_airbld = shc_airbld
1333 
1334  !internal convective exchange coefficients !!FO!! ibldCHmod = 0 originally
1335  !iBldCHmod specifies method for convective exchange coeffs
1336  IF (ibldchmod == 1) THEN !ASHRAE 2001
1337  ch_ibld = 1.31*(abs(t0_ibld - tievolve))**0.25/shc_airbld
1338  ch_iwall = 1.31*(abs(tn_wall - tievolve))**0.25/shc_airbld
1339  ch_iroof = 1.52*(abs(tn_roof - tievolve))**0.25/shc_airbld
1340  IF (abs(tn_roof - tievolve) > 0) ch_iroof = ch_iroof*0.39 !effect of convection is weaker downward
1341  ELSEIF (ibldchmod == 2) THEN !Awbi, H.B. 1998, Energy and Buildings 28: 219-227
1342  ch_ibld = 1.823*(abs(t0_ibld - tievolve))**0.293/shc_airbld
1343  ch_iwall = 1.823*(abs(tn_wall - tievolve))**0.293/shc_airbld
1344  ch_iroof = 2.175*(abs(tn_roof - tievolve))**0.308/shc_airbld
1345  IF (abs(tn_roof - tievolve) > 0) ch_iroof = 0.704*(abs(tn_roof - tievolve))**0.133/shc_airbld !effect of convection is weaker downward
1346  ENDIF
1347 
1348  !Evolving T = (Previous Temp + dT from Sensible heat flux) mixed with outside air
1349  !ASSUMES THE CH_BLD INCLUDES THE EFFECT OF VENTILATION RATE IN m/s (e.g. if a normal CH is .005 and
1350  !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
1351  tairmix = (tievolve + tair1*airexdt)/(1.0 + airexdt)
1352  qfbld = froof*(tievolve - tairmix)*shc_airbld*bldghx/tstep !heat added or lost, requires cooling or heating if HVAC on
1353 
1354  !!FO!! CH_xxxx has unit [m/s] !!**HCW what is going on with tstep here??
1355  tievolve = tairmix + tstep/bldghx/finternal & !!FO!! finternal(=froof+fibld+fwall) => normalisation of fractions
1356  *(ch_ibld*fibld*(t0_ibld - tievolve) &
1357  + ch_iroof*froof*(tn_roof - tievolve) &
1358  + ch_iwall*fwall*(tn_wall - tievolve)) !!FO!! [K] = [K] + [s/m]*([m/s]*([K]))
1359 
1360  IF (.NOT. diagnoseti) tievolve = tinternal + c2k
1361  IF (hvac) THEN !Run up/down to set point +/- 1 degree with adjustment of 90% per hour
1362  tadd = (sign(-1.0d0, theat_fix - tievolve) + theat_fix - tievolve)*min(4.*tstep/3600.0, 0.9) !!**HCW check??
1363  tievolve = tievolve + tadd
1364  ENDIF
1365 
1366  !========>RADIATION<================
1367  IF (kdn_estm < 0) kdn_estm = 0. !set non-zero shortwave to zero !Should this be moved up to line 183/4?
1368 
1369  !external components, no diffuse
1370  !for reflections complete absorption is assumed
1371  !for shortwave these are net values
1372  !for longwave these are incoming only
1373  !MUST DIVIDE SHORTWAVE INTO DIRECT AND DIFFUSE
1374  sw_hor = kdn_estm !incoming solar on horizontal surface
1375  sw_vert = kdn_estm*tanzenith !incoming solar on vertical surface = kdown(obs)*sin(zenith)/cos(zenith)
1376 
1377  rs_roof = svf_roof*(1.0 - alb_roof)*sw_hor
1378  rl_roof = svf_roof*em_roof*ldown
1379 
1380  rs_ground = svf_ground*(1.-alb_ground)*sw_hor + &
1381  zvf_ground*svf_wall*alb_wall*sw_vert*(1 - alb_ground) + &
1383 
1385 
1386  rs_wall = svf_wall*(1.-alb_wall)*sw_vert + &
1387  zvf_wall*svf_wall*alb_wall*sw_vert*(1.+zvf_wall*alb_wall) + &
1388  xvf_wall*svf_ground*alb_ground*sw_hor*(1 - alb_wall) + &
1390 
1391  !wall to wall exchange handled simultaneously with seb calc
1392  rl_wall = svf_wall*ldown*em_wall + zvf_wall*svf_wall*ldown*(1 - em_wall)*em_wall + &
1394 
1395  !DIFFICULT TO DETERMINE WHAT THIS IS EXACTLY, DONT INCLUDE WALLS
1396  kup_estm = kdn_estm - rvf_roof*rs_roof - (rvf_ground + rvf_wall)*rs_ground/svf_ground - rvf_veg*alb_veg*kdn_estm
1397  IF (kdn_estm > 10 .AND. kup_estm > 0) THEN
1398  alb_avg = kup_estm/kdn_estm
1399  sumalb = sumalb + alb_avg
1400  nalb = nalb + 1
1401  ENDIF
1402 
1403  !internal components
1404  rs_ibld = 0 ! This could change if there are windows (need solar angles or wall svf * fraction glazing * transmissivity)
1405  !internal incoming longwave terms do not include the view factors for its own surface e.g. for ibld and walls
1406  !added floor view factors
1407  rl_ibld = sbconst*(ivf_iw*em_w*tn_wall**4 + &
1408  ivf_ir*em_r*tn_roof**4 + &
1409  ivf_if*em_f*tfloor**4)
1410  rs_iwall = 0
1411  rl_iwall = sbconst*(ivf_wi*em_i*t0_ibld**4 + &
1412  ivf_wr*em_r*tn_roof**4 + &
1413  ivf_wf*em_f*tfloor**4)
1414  rs_iroof = 0
1415  rl_iroof = sbconst*(ivf_ri*em_i*t0_ibld**4 + &
1416  ivf_rw*em_w*tn_wall**4 + &
1417  ivf_rf*em_f*tfloor**4)
1418 
1419  !========>INTERNAL<================
1420  bctype = .false.
1421  kdz = 2*kibld(1)/zibld(1)
1422  pcoeff = (/em_ibld*sbconst*(1 - ivf_ii*em_ibld), 0.0d0, 0.0d0, kdz + shc_airbld*ch_ibld, &
1423  -kdz*tibld(1) - shc_airbld*ch_ibld*tievolve - rs_ibld - rl_ibld/)
1425 
1426  !!FO!! this leads to Tibld(1) = Tibld(3) , i.e. ...
1427  bc(1) = t0_ibld
1428 
1429  !!FO!! temperature equal on both sides of inside wall
1430  bc(2) = bc(1)
1431 
1432  ! print*,'ESTM Qsibld cal before',Qsibld,Tibld,zibld(1:Nibld)
1433  CALL heatcond1d(tibld, qsibld, zibld(1:nibld), REAL(Tstep, KIND(1d0)), kibld(1:nibld), ribld(1:nibld), bc, bctype)
1434 
1435  ! print*,'ESTM Qsibld cal after',Qsibld,Tibld
1436 
1437  !========>WALLS<================
1438  bctype = .false.
1439  kdz = 2*kwall(nwall)/zwall(nwall)
1440  pcoeff = (/em_ibld*sbconst*(1 - ivf_ww*em_ibld), 0.0d0, 0.0d0, kdz + shc_airbld*ch_iwall, &
1441  -kdz*twall(nwall) - shc_airbld*ch_iwall*tievolve - rs_iwall - rl_iwall/)
1443  bc(2) = tn_wall !!FO!! boundary condition #2 = inner surface Twall, originally from lodz_parms_ltm.txt or finaltemp.txt
1444 
1445  IF (tsurfchoice < 2 .OR. radforce) THEN
1446  IF (radforce) THEN !!FO!! 1st prio: radforce
1447  kdz = 2*kwall(1)/zwall(1)
1448  pcoeff = (/em_wall*sbconst*(1 - zvf_wall*em_wall), 0.0d0, 0.0d0, kdz + shc_air*chair_wall*ws, &
1449  -kdz*twall(1) - shc_air*chair_wall*ws*tair1 - rs_wall - rl_wall/)
1451  bc(1) = t0_wall !!FO!! boundary condition #1 = outer surface Twall, originally from lodz_parms_ltm.txt or finaltemp.txt
1452  ELSEIF (tsurfchoice == 0) THEN
1453  bc(1) = tsurf_all + c2k; t0_wall = bc(1)
1454  ELSEIF (tsurfchoice == 1) THEN
1455  bc(1) = twall_all + c2k; t0_wall = bc(1)
1456  ENDIF !!FO!! Tsoil in Lodz2002HS.txt NB => Lodz2002HS.txt doesn't work with onewall = TRUE
1457 
1458  CALL heatcond1d(twall, qswall, zwall(1:nwall), REAL(Tstep, KIND(1d0)), kwall(1:nwall), rwall(1:nwall), bc, bctype) !!FO!! new set of Twalls are calculated from heat conduction through wall
1459 
1460  ELSEIF (tsurfchoice == 2) then!SPECIAL FOR 4 WALLS
1461  t0_wall = 0.
1462  DO i = 1, 4 !do 4 walls
1463  bc(1) = tw_n + tw_e + tw_s + tw_w + c2k; t0_wall = t0_wall + bc(1)
1464  CALL heatcond1d(tw_4(:, i), qs_4(i), zwall(1:nwall), REAL(Tstep, KIND(1d0)), kwall(1:nwall), rwall(1:nwall), bc, bctype)
1465  ENDDO
1466  !Take average of 4 wall values
1467  t0_wall = t0_wall/4.
1468  qswall = sum(qs_4)/4.
1469  twall = sum(tw_4, 2)/4.
1470  ENDIF
1471 
1472  !========>ROOF<================
1473  bctype = .false.
1474  kdz = 2*kroof(nroof)/zroof(nroof)
1475  pcoeff = (/em_ibld*sbconst, 0.0d0, 0.0d0, kdz + shc_airbld*ch_iroof, &
1476  -kdz*troof(nroof) - shc_airbld*ch_iroof*tievolve - rs_iroof - rl_iroof/)
1478  bc(2) = tn_roof
1479 
1480  IF (radforce) THEN
1481  kdz = 2*kroof(1)/zroof(1)
1482  pcoeff = (/em_roof*sbconst, 0.0d0, 0.0d0, kdz + shc_air*chair*ws, &
1483  -kdz*troof(1) - shc_air*chair*ws*tair1 - rs_roof - rl_roof/)
1485  bc(1) = t0_roof
1486  ELSEIF (tsurfchoice == 0) THEN
1487  bc(1) = tsurf_all + c2k; t0_roof = bc(1)
1488  ELSE
1489  bc(1) = troof_in + c2k; t0_roof = bc(1)
1490  ENDIF
1491 
1492  CALL heatcond1d(troof, qsroof, zroof(1:nroof), REAL(Tstep, KIND(1d0)), kroof(1:nroof), rroof(1:nroof), bc, bctype)
1493 
1494  !========>ground<================
1495  bctype = .false.
1496  kdz = 2*kground(1)/zground(1)
1497 
1498  IF (radforce .OR. groundradforce) THEN
1499  pcoeff = (/em_ground*sbconst, 0.0d0, 0.0d0, kdz + shc_air*chair_ground*ws, &
1500  -kdz*tground(1) - shc_air*chair_ground*ws*tair1 - rs_ground - rl_ground/)
1502  bc(1) = t0_ground
1503  ELSEIF (tsurfchoice == 0) THEN
1504  bc(1) = tsurf_all + c2k; t0_ground = bc(1)
1505  ELSE
1506  bc(1) = troad + c2k; t0_ground = bc(1)
1507  ENDIF
1508 
1509  bc(2) = lbc_soil + c2k
1510  ! bc(2)=0.; bctype(2)=.t.
1511 
1512  IF (fground /= 0.) THEN ! check fground==0 scenario to avoid division-by-zero error, TS 21 Jul 2016
1514  REAL(Tstep, KIND(1d0)), kground(1:nground), rground(1:nground), bc, bctype)
1515  ELSE
1516  qsground = nan
1517  END IF
1518 
1520  qsibld = qsibld*fibld
1521  qswall = qswall*fwall
1522  qsroof = qsroof*froof
1524  qs = qsibld + qswall + qsroof + qsground !!FO!! QSair not included; called QS in output file (column #10)
1525 
1526  ! print*,'ESTM QS',qs,Qsibld,Qswall,Qsroof ,Qsground
1527  ! print*,'ESTM Qsibld',Qsibld,fibld
1528  !write(*,*) Qsair, QSibld, Qswall, Qsroof, Qsground, QS
1529 
1530  !========>Radiation<================
1531  !note that the LUP for individual components does not include reflected
1535  tveg = tair1
1536  lup_veg = sbconst*em_veg*tveg**4
1539  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
1540  rn_ground = rs_ground + rl_ground - lup_ground
1541  rn_roof = rs_roof + rl_roof - lup_roof
1542  rn_wall = rs_wall + rl_wall - lup_wall*(1 - zvf_wall*em_wall)
1543  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
1544  qhestm = (t0 - tair1)*chair*shc_air*ws
1545  sumemis = sumemis + em_equiv
1546  nemis = nemis + 1
1547 
1548  ! IF (SPINDONE) THEN !!FO!! only the last set of values in the time interpolation loop is written to file
1549 
1550  IF (nwall < 5) THEN
1551  twallout = (/twall, (dum(ii), ii=1, (5 - nwall))/)
1552  ELSE
1553  twallout = twall
1554  ENDIF
1555 
1556  IF (nroof < 5) THEN
1557  troofout = (/troof, (dum(ii), ii=1, (5 - nroof))/);
1558  ELSE
1559  troofout = troof
1560  ENDIF
1561 
1562  IF (nground < 5) THEN
1563  tgroundout = (/tground, (dum(ii), ii=1, (5 - nground))/)
1564  ELSE
1565  tgroundout = tground
1566  ENDIF
1567 
1568  IF (nibld < 5) THEN
1569  tibldout = (/tibld, (dum(ii), ii=1, (5 - nibld))/)
1570  ELSE
1571  tibldout = tibld
1572  ENDIF
1573 
1574  ! dataOutESTM(ir,1:ncolumnsDataOutESTM,Gridiv)=[&
1575  ! REAL(iy,KIND(1D0)),REAL(id,KIND(1D0)),REAL(it,KIND(1D0)),REAL(imin,KIND(1D0)), dectime,&!5
1576  ! QS,Qsair,Qswall,Qsroof,Qsground,Qsibld,&!11
1577  ! Twallout,Troofout,Tgroundout,Tibldout,Tievolve]!32 !NB. These all have 5 elements except Tievolve (1).
1578  dataoutlineestm = [ &
1579  qs, qsair, qswall, qsroof, qsground, qsibld, &!6
1580  twallout, troofout, tgroundout, tibldout, tievolve]!27 !NB. These all have 5 elements except Tievolve (1).
1581  ! set invalid values to nan
1582  dataoutlineestm = set_nan(dataoutlineestm)
1583  ! call r8vec_print(ncolumnsDataOutESTM-5,dataOutESTM(ir,6:ncolumnsDataOutESTM,Gridiv),'dataOutESTM')
1584 
1585  tair2 = tair1
1586 
1587  ! Save variables for this grid
1588  tair2_grids(gridiv) = tair1
1589  lup_ground_grids(gridiv) = lup_ground
1590  lup_wall_grids(gridiv) = lup_wall
1591  lup_roof_grids(gridiv) = lup_roof
1592  tievolve_grids(gridiv) = tievolve
1593  t0_ibld_grids(gridiv) = t0_ibld
1594  t0_ground_grids(gridiv) = t0_ground
1595  t0_wall_grids(gridiv) = t0_wall
1596  t0_roof_grids(gridiv) = t0_roof
1597  tn_wall_grids(gridiv) = tn_wall
1598  tn_roof_grids(gridiv) = tn_roof
1599  tground_grids(:, gridiv) = tground(:)
1600  twall_grids(:, gridiv) = twall(:)
1601  troof_grids(:, gridiv) = troof(:)
1602  tibld_grids(:, gridiv) = tibld(:)
1603  tw_4_grids(:, :, gridiv) = tw_4(:, :)
1604 
integer, parameter cts_tsurf
real(kind(1d0)) alb_avg
real(kind(1d0)) tair1
real(kind(1d0)) ch_iwall
real(kind(1d0)) hw
integer, parameter cts_troad
real(kind(1d0)), parameter sbconst
real(kind(1d0)) qsibld
real(kind(1d0)) rvf_wall
real(kind(1d0)) fwall
real(kind(1d0)) ch_ibld
real(kind(1d0)), dimension(:), allocatable lup_wall_grids
real(kind(1d0)), dimension(:, :), allocatable troof_grids
integer, parameter maxiter
real(kind(1d0)) froof
real(kind(1d0)) alb_veg
real(kind(1d0)) ws
real(kind(1d0)), dimension(:), allocatable lup_ground_grids
real(kind(1d0)) fground
real(kind(1d0)) function heatcapacity_air(TK, RH, P)
real(kind(1d0)), dimension(:, :, :), allocatable tw_4_grids
real(kind(1d0)) ivf_wf
real(kind(1d0)), parameter alb_wall
real(kind(1d0)) em_r
real(kind(1d0)) zenith_deg
real(kind(1d0)) ivf_wr
real(kind(1d0)) tn_wall
integer, parameter cts_troof
real(kind(1d0)) ivf_ir
real(kind(1d0)) avrh
real(kind(1d0)) nan
integer, parameter cts_twall
real(kind(1d0)), dimension(:), allocatable t0_roof_grids
real(kind(1d0)) zvf_wall
real(kind(1d0)) tfloor
real(kind(1d0)) temp_c
real(kind(1d0)), dimension(5) rroof
real(kind(1d0)) lbc_soil
real(kind(1d0)) qsroof
real(kind(1d0)) avu1
real(kind(1d0)), dimension(5) kibld
real(kind(1d0)), dimension(:), allocatable twall
real(kind(1d0)), dimension(:), allocatable lup_roof_grids
real(kind(1d0)) lup_roof
real(kind(1d0)), dimension(5) rwall
real(kind(1d0)), dimension(5) kwall
real(kind(1d0)), dimension(:), allocatable tibld
real(kind(1d0)) em_roof
real(kind(1d0)) ivf_rw
real(kind(1d0)), dimension(5) kroof
logical, dimension(2) bctype
integer, parameter cts_tiair
real(kind(1d0)) fair
real(kind(1d0)) function newtonpolynomial(x0, Pcoeff, conv, maxiter)
real(kind(1d0)) qswall
real(kind(1d0)), dimension(:), allocatable t0_ibld_grids
real(kind(1d0)) ldown
real(kind(1d0)) t0_ibld
real(kind(1d0)), dimension(27) dataoutlineestm
real(kind(1d0)) ch_iroof
real(kind(1d0)), dimension(5) zwall
real(kind(1d0)) rvf_roof
real(kind(1d0)) tair_av
real(kind(1d0)) rvf_ground
integer, parameter cts_twall_n
real(kind(1d0)) t0_ground
real(kind(1d0)), parameter em_wall
integer tsurfchoice
real(kind(1d0)) finternal
real(kind(1d0)) chair
real(kind(1d0)) theat_fix
integer, parameter cts_twall_w
real(kind(1d0)), dimension(:), allocatable tn_wall_grids
real(kind(1d0)) minshc_airbld
integer, parameter cts_twall_e
real(kind(1d0)) qsground
real(kind(1d0)), dimension(:, :), allocatable tground_grids
real(kind(1d0)) t0_wall
real(kind(1d0)) zvf_ground
real(kind(1d0)), dimension(:), allocatable tground
real(kind(1d0)), dimension(:, :), allocatable tibld_grids
real(kind(1d0)) em_veg
real(kind(1d0)) avkdn
real(kind(1d0)), dimension(5) zibld
real(kind(1d0)) alb_ground
real(kind(1d0)) ivf_ri
integer ibldchmod
real(kind(1d0)), dimension(:), allocatable tair2_grids
real(kind(1d0)), dimension(:, :), allocatable tw_4
real(kind(1d0)) ivf_iw
real(kind(1d0)) theat_on
real(kind(1d0)) svf_roof
integer, parameter ncolsestmdata
real(kind(1d0)) chr
real(kind(1d0)), dimension(:), allocatable tn_roof_grids
real(kind(1d0)) ivf_wi
real(kind(1d0)), dimension(5) ribld
real(kind(1d0)) em_i
real(kind(1d0)), dimension(:), allocatable troof
real(kind(1d0)) ivf_ww
real(kind(1d0)) rvf_veg
real(kind(1d0)) ivf_if
integer, parameter cts_twall_s
real(kind(1d0)) lup_wall
subroutine heatcond1d(T, Qs, dx, dt, k, rhocp, bc, bctype)
real(kind(1d0)), dimension(5) zroof
real(kind(1d0)) em_w
real(kind(1d0)), dimension(:), allocatable tievolve_grids
real(kind(1d0)) xvf_wall
real(kind(1d0)), dimension(5) rground
real(kind(1d0)), dimension(4) qs_4
real(kind(1d0)) em_ground
real(kind(1d0)) t0_roof
real(kind(1d0)), dimension(5) kground
logical diagnoseti
real(kind(1d0)) fibld
real(kind(1d0)) sumemis
real(kind(1d0)) tn_roof
real(kind(1d0)) qsair
real(kind(1d0)), dimension(5) zground
real(kind(1d0)) tair2
real(kind(1d0)), parameter conv
real(kind(1d0)) bldgh
real(kind(1d0)) tievolve
real(kind(1d0)) svf_wall
real(kind(1d0)) lup_veg
real(kind(1d0)), parameter pi
real(kind(1d0)), dimension(:), allocatable t0_ground_grids
real(kind(1d0)), dimension(5) pcoeff
real(kind(1d0)) theat_off
real(kind(1d0)), dimension(:), allocatable t0_wall_grids
integer evolvetibld
real(kind(1d0)) sumalb
real(kind(1d0)) alb_roof
real(kind(1d0)) shc_air
real(kind(1d0)) em_ibld
real(kind(1d0)), dimension(:), allocatable ts5mindata_ir
real(kind(1d0)) svf_ground
real(kind(1d0)) lup_ground
real(kind(1d0)), parameter c2k
real(kind(1d0)) ivf_ii
real(kind(1d0)) press_hpa
real(kind(1d0)) tanzenith
real(kind(1d0)) ivf_rf
integer nground
real(kind(1d0)) em_f
real(kind(1d0)), dimension(:, :), allocatable twall_grids
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 590 of file suews_phys_estm.f95.

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

590 
591  ! Last modified HCW 30 Jun 2016 - reading in now done by SUEWS_GetESTMData subroutine.
592  ! ESTM_initials now only runs once per run at the very start.
593  ! Last modified HCW 15 Jun 2016 - code now reads in 5-min file (interpolation done beforehand, outside of SUEWS itself)
594 
595  USE defaultnotused
596  USE physconstants, ONLY: c2k
597  USE estm_data
598  USE allocatearray
599  USE data_in, ONLY: fileinputpath
600  USE initial, ONLY: numberofgrids
601 
602  IMPLICIT NONE
603 
604  !=====Read ESTMinput.nml================================
605  namelist /estminput/ tsurfchoice, &
606  evolvetibld, &
607  ibldchmod, &
608  lbc_soil, &
609  theat_on, &
610  theat_off, &
611  theat_fix
612 
613  OPEN (511, file=trim(fileinputpath)//'ESTMinput.nml', status='old')
614  READ (511, nml=estminput)
615  CLOSE (511)
616 
617  !Convert specified temperatures to Kelvin
618  theat_on = theat_on + c2k
621 
622  ALLOCATE (tair2_grids(numberofgrids))
623  ALLOCATE (lup_ground_grids(numberofgrids))
624  ALLOCATE (lup_wall_grids(numberofgrids))
625  ALLOCATE (lup_roof_grids(numberofgrids))
626  ALLOCATE (tievolve_grids(numberofgrids))
627  ALLOCATE (t0_ibld_grids(numberofgrids))
628  ALLOCATE (t0_ground_grids(numberofgrids))
629  ALLOCATE (t0_wall_grids(numberofgrids))
630  ALLOCATE (t0_roof_grids(numberofgrids))
631  ALLOCATE (tn_wall_grids(numberofgrids))
632  ALLOCATE (tn_roof_grids(numberofgrids))
633 
integer numberofgrids
real(kind(1d0)), dimension(:), allocatable lup_wall_grids
real(kind(1d0)), dimension(:), allocatable lup_ground_grids
real(kind(1d0)), dimension(:), allocatable t0_roof_grids
real(kind(1d0)) lbc_soil
real(kind(1d0)), dimension(:), allocatable lup_roof_grids
real(kind(1d0)), dimension(:), allocatable t0_ibld_grids
integer tsurfchoice
real(kind(1d0)) theat_fix
real(kind(1d0)), dimension(:), allocatable tn_wall_grids
character(len=150) fileinputpath
integer ibldchmod
real(kind(1d0)), dimension(:), allocatable tair2_grids
real(kind(1d0)) theat_on
real(kind(1d0)), dimension(:), allocatable tn_roof_grids
real(kind(1d0)), dimension(:), allocatable tievolve_grids
real(kind(1d0)), dimension(:), allocatable t0_ground_grids
real(kind(1d0)) theat_off
real(kind(1d0)), dimension(:), allocatable t0_wall_grids
integer evolvetibld
real(kind(1d0)), parameter c2k
Here is the caller graph for this function:

◆ estm_translate()

subroutine estm_module::estm_translate ( integer  Gridiv)

Definition at line 638 of file suews_phys_estm.f95.

References allocatearray::alb, estm_data::alb_avg, estm_data::alb_ground, estm_data::alb_roof, estm_data::alb_veg, estm_data::alb_wall, 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_data::em_i, estm_data::em_ibld, estm_data::em_r, estm_data::em_roof, estm_data::em_veg, estm_data::em_w, estm_data::em_wall, 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::nemis, estm_data::nground, estm_data::nibld, estm_data::nroof, estm_data::nroom, initial::numberofgrids, estm_data::nwall, 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, 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().

638  ! HCW 30 Jun 2016
639 
640  USE defaultnotused, ONLY: nan
641  USE physconstants, ONLY: c2k, sbconst
642  USE estm_data
643  USE allocatearray
644  USE gis_data, ONLY: bldgh
645  USE initial, ONLY: numberofgrids
646 
647  IMPLICIT NONE
648  INTEGER :: i
649  !REAL(KIND(1d0)) :: CFLval
650  !REAL(KIND(1d0)) :: t5min
651  REAL(KIND(1d0))::w, wb
652  !CHARACTER (len=20)::FileCodeX
653  !CHARACTER (len=150):: FileFinalTemp
654  !LOGICAL:: inittemps=.FALSE.
655  INTEGER:: estmstart = 0
656  INTEGER:: gridiv
657 
658  !Set initial values at the start of each run for each grid
659  ! the initiliastion part is problematic:
660  ! cannot be initialised under a multi-grid scenario
661  IF (gridiv == 1) estmstart = estmstart + 1
662  ! print *, 'gridiv in ESTM',gridiv,'ESTMStart',ESTMStart
663  IF (estmstart == 1) THEN
664 
665  ! write(*,*) ' ESTMStart: ',ESTMStart, 'initialising ESTM for grid no. ', Gridiv
666 
667  tfloor = 20.0 ! This is used only when radforce =T !TODO: should be put in the namelist
668  tfloor = tfloor + c2k
669 
670  ! Initial values
671  tievolve = 20.0 + c2k
672  shc_air = 1230.0
673  minshc_airbld = 1300
674 
675  ! ---- Internal view factors ----
676  !constant now but should be calculated in the future
677  ivf_iw = 0.100000
678  ivf_ir = 0.000000
679  ivf_ii = 0.900000
680  ivf_if = 0.000000
681  ivf_ww = 0.050000
682  ivf_wr = 0.000000
683  ivf_wi = 0.950000
684  ivf_wf = 0.000000
685  ivf_rw = 0.050000
686  ivf_ri = 0.950000
687  ivf_rf = 0.000000
688  ivf_fw = 0.050000
689  ivf_fr = 0.000000
690  ivf_fi = 0.950000
691 
692  tair24hr = c2k
693 
694  !Ts5mindata(1,ncolsESTMdata) = -999
695  ! !Fill Ts5mindata for current grid and met block - this is done in SUEWS_translate
697 
698  ! ---- Initialization of variables and parameters for first row of run for each grid ----
699  ! N layers are calculated in SUEWS_translate
700  IF (.NOT. ALLOCATED(tibld)) THEN
701  ! print *, "Nibld", Nibld
702  ! print *, "Nwall", Nwall
703  ! print *, "Nroof", Nroof
704  ! print *, "Nground", Nground
705  ALLOCATE (tibld(nibld), twall(nwall), troof(nroof), tground(nground), tw_4(nwall, 4))
706  ALLOCATE (tibld_grids(nibld, numberofgrids), &
711  ENDIF
712 
713  ! Transfer variables from Ts5mindata to variable names
714  ! N.B. column numbers here for the following file format - need to change if input columns change!
715  ! dectime iy id it imin Tiair Tsurf Troof Troad Twall Twall_n Twall_e Twall_s Twall_w
716  ! 1 2 3 4 5 6 7 8 9 10 11 12 13 !new
717  !
718  ! Calculate temperature of each layer in Kelvin
719  ! QUESTION: what if (Nground/Nwall/Nroof-1)==0? TS 21 Oct 2017
720  DO i = 1, nground
721  tground(i) = (lbc_soil - ts5mindata(1, cts_troad))*(i - 1)/(nground - 1) + ts5mindata(1, cts_troad) + c2k
722  ENDDO
723  DO i = 1, nwall
724  twall(i) = (ts5mindata(1, cts_tiair) - ts5mindata(1, cts_twall))*(i - 1)/(nwall - 1) + ts5mindata(1, cts_twall) + c2k
725  ENDDO
726  DO i = 1, nroof
727  troof(i) = (ts5mindata(1, cts_tiair) - ts5mindata(1, cts_troof))*(i - 1)/(nroof - 1) + ts5mindata(1, cts_troof) + c2k
728  ENDDO
729  tibld(1:nibld) = ts5mindata(1, cts_tiair) + c2k
730 
731  ENDIF !End of loop run only at start (for each grid)
732 
733  ! ---- Parameters related to land surface characteristics ----
734  ! QUESTION: Would Zref=z be more appropriate?
735  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)
736 
737  svf_ground = 1.0
738  svf_roof = 1.0
739 
740  ! ==== roof (i.e. Bldgs)
741  !froof=sfr(BldgSurf) ! Moved to SUEWS_translate HCW 16 Jun 2016
744 
745  ! ==== vegetation (i.e. EveTr, DecTr, Grass)
746  !fveg=sfr(ConifSurf)+sfr(DecidSurf)+sfr(GrassSurf) ! Moved to SUEWS_translate HCW 16 Jun 2016
747  IF (fveg /= 0) THEN
750  ELSE ! check fveg==0 scenario to avoid division-by-zero error, TS 21 Oct 2017
751  alb_veg = nan
752  em_veg = nan
753  ENDIF
754 
755  ! ==== ground (i.e. Paved, EveTr, DecTr, Grass, BSoil, Water - all except Bldgs)
756  !fground=sfr(ConifSurf)+sfr(DecidSurf)+sfr(GrassSurf)+sfr(PavSurf)+sfr(BsoilSurf)+sfr(WaterSurf) ! Moved to SUEWS_translate HCW 16 Jun 2016
757  IF (fground /= 0) THEN
764  ELSE ! check fground==0 scenario to avoid division-by-zero error, TS 21 Jul 2016
765  alb_ground = nan
766  em_ground = nan
767  ENDIF
768 
769  IF (froof < 1.0) THEN
770  hw = fwall/(2.0*(1.0 - froof))
771  ELSE
772  ! HW=0 !HCW if only roof, no ground
773  hw = 0.00001 ! to avoid zero-HW scenario TS 21 Oct 2017
774 
775  END IF
776  hw = max(0.00001, hw)! to avoid zero-HW scenario TS 27 Oct 2017
777 
778  IF (fground == 1.0) THEN !!FO!! if only ground, i.e. no houses
779  w = 1
780  wb = 0
781  svf_ground = 1.
782  zvf_wall = 0.
783  svf_wall = 0.
784  svf_roof = 1.
785  zvf_ground = 0.
786  xvf_wall = 0.
787  rvf_canyon = 1.
788  rvf_ground = 1.-fveg
789  rvf_roof = 0
790  rvf_wall = 0
791  rvf_veg = fveg
792  ELSE IF (fground == 0.0) THEN !check fground==0 (or HW==0) scenario to avoid division-by-zero error, TS 21 Jul 2016
793  ! the following values are calculated given HW=0
794  w = 0
795  wb = 1
796  zvf_wall = 0 !COS(ATAN(2/HW)) when HW=0 !!FO!! wall view factor for wall
797  hw = 0
798  svf_ground = max(cos(atan(2*hw)), 0.00001)!!FO!! sky view factor for ground ! to avoid zero-division scenario TS 21 Oct 2017
799  svf_wall = (1 - zvf_wall)/2 !!FO!! sky view factor for wall
800  zvf_ground = 1 - svf_ground !!FO!! wall view factor for ground
801  xvf_wall = svf_wall !!FO!! ground view factor
802  ! RVF_CANYON=COS(ATAN(2*ZREF/W))
803  ! RVF_ROOF=1-RVF_CANYON
804  ! RVF_WALL=(COS(ATAN(2*(ZREF-BldgH)/W))-RVF_CANYON)*RVF_CANYON
805  ! RVF_ground=RVF_CANYON-RVF_WALL
808  rvf_roof = froof
810  ELSE
811  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
812  wb = w*sqrt(froof/fground)
813  svf_ground = cos(atan(2*hw)) !!FO!! sky view factor for ground
814  zvf_wall = cos(atan(2/hw)) !!FO!! wall view factor for wall
815  svf_wall = (1 - zvf_wall)/2 !!FO!! sky view factor for wall
816  zvf_ground = 1 - svf_ground !!FO!! wall view factor for ground
817  xvf_wall = svf_wall !!FO!! ground view factor
818  ! RVF_CANYON=COS(ATAN(2*ZREF/W))
819  ! RVF_ROOF=1-RVF_CANYON
820  ! RVF_WALL=(COS(ATAN(2*(ZREF-BldgH)/W))-RVF_CANYON)*RVF_CANYON
821  ! RVF_ground=RVF_CANYON-RVF_WALL
824  rvf_roof = froof
826  ENDIF
827 
829 
830  sumalb = 0.; nalb = 0
831  sumemis = 0.; nemis = 0
832 
833  !set emissivity for ceiling, wall and floor inside of buildings
835 
836  !internal elements
837  IF (nroom == 0) THEN
838  fibld = (floor(bldgh/3.1 - 0.5) - 1)*froof
839  ELSE
840  fibld = (2.-2./nroom)*fwall + (floor(bldgh/3.1 - 0.5) - 1)*froof
841  ENDIF
842 
843  IF (fibld == 0) fibld = 0.00001 !this just ensures a solution to radiation
844  finternal = froof + fibld + fwall
845  IF (finternal == 0) finternal = 0.00001 ! to avoid zero-devision error TS 21 Oct 2017
846  fair = zref - bldgh*froof
847  IF (fair == 0) fair = 0.00001 ! to avoid zero-devision error TS 21 Oct 2017
848  !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
849  !ivf_ww=1.-ivf_wi-ivf_wr-ivf_wf
850  !ivf_rw=1.-ivf_ri-ivf_rf;
851  !ivf_fr=ivf_rf;
852 
853  IF ((ivf_ii + ivf_iw + ivf_ir + ivf_if > 1.0001) .OR. &
854  (ivf_wi + ivf_ww + ivf_wr + ivf_wf > 1.0001) .OR. &
855  (ivf_ri + ivf_rw + ivf_rf > 1.0001) .OR. &
856  (ivf_fi + ivf_fw + ivf_fr > 1.0001) .OR. &
857  (ivf_ii + ivf_iw + ivf_ir + ivf_if < 0.9999) .OR. &
858  (ivf_wi + ivf_ww + ivf_wr + ivf_wf < 0.9999) .OR. &
859  (ivf_ri + ivf_rw + ivf_rf < 0.9999) .OR. &
860  (ivf_fi + ivf_fw + ivf_fr < 0.9999)) THEN
861  print *, "At least one internal view factor <> 1. Check ivf in ESTMinput.nml"
862  ENDIF
863 
864  !=======Initial setting==============================================
865  !! Rewritten by HCW 15 Jun 2016 to use existing SUEWS error handling
866  !IF(inittemps) THEN
867  ! write(*,*) 'inittemps:',inittemps
868  ! FileFinalTemp=TRIM(FileOutputPath)//TRIM(FileCodeX)//'_ESTM_finaltemp.txt'
869  ! OPEN(99,file=TRIM(FileFinalTemp),status='old',err=316) ! Program stopped if error opening file
870  ! READ(99,*) Twall,Troof,Tground,Tibld ! Twall, Troof, Tground & Tibld get new values
871  ! CLOSE(99)
872  !ENDIF
873  !
874  !!IF (inittemps) THEN !!FO!! inittemps=.true. set in nml file
875  !! OPEN(99,file='outputfiles/finaltemp.txt',status='old',iostat=ios) !!FO!! has to exist
876  !!
877  !! IF (ios/=0) CALL error('outputfiles/finaltemp.txt',ios,1) !!FO!! calls mod_error.f95, writes that the opening failed and stops prg
878  !! IF (ios/=0) THEN
879  !! Twall = (/273., 285., 291./)
880  !! Troof = (/273., 285., 291./)
881  !! Tground = (/273., 275., 280., 290./)
882  !! Tibld = (/293., 293., 293./)
883  !! ELSE
884  !! READ(99,*) Twall,Troof,Tground,Tibld !!FO!! if finaltemp.txt exists Twall[3], Troof[3], Tground[4] & Tibld[3] get new values
885  !! CLOSE(99)
886  !! ENDIF
887  !!ENDIF
888 
889  !where (isnan(Twall))
890  ! Twall = 273
891  !endwhere
892  !where (isnan(Troof))
893  ! Troof = 273
894  !endwhere
895  !where (isnan(Tground))
896  ! Tground = 281
897  !endwhere
898  !where (isnan(Tibld))
899  ! Tibld = 293
900  !endwhere
901 
902  IF (estmstart == 1) THEN
903  DO i = 1, 4
904  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
905  ENDDO
906 
907  !initialize surface temperatures
908  t0_ground = tground(1)
909  t0_wall = twall(1)
910  t0_roof = troof(1)
911  t0_ibld = tibld(1)
912  tn_roof = troof(nroof)
913  tn_wall = twall(nwall)
914 
915  !initialize outgoing longwave !QUESTION: HCW - Are these calculations compatible with those in LUMPS_NARP?
919 
920  ! PRINT*,"W,WB= ",W,WB
921  ! PRINT*,'SVF_ground ','SVF_WALL ','zvf_WALL ','HW '
922  ! PRINT*,SVF_ground,SVF_WALL,zvf_WALL,HW
923  ! PRINT*,'RVF_ground ','RVF_WALL ','RVF_ROOF ','RVF_VEG'
924  ! PRINT*,RVF_ground,RVF_WALL,RVF_ROOF,RVF_VEG
925  ! print*,'Alb_avg (VF)=',alb_avg
926  ! print*,'z0m, Zd', z0m, ZD
927 
928  ENDIF
929 
930  first = .true.
931 
932  !======Courant-Friedrichs-Lewy condition=================================
933  !This is comment out by S.O. for now
934  ! NB: should be recovered for calculation stability, FO and TS, 11 Oct 2017
935  ! CFLval = minval(0.5*zibld*zibld*ribld/kibld) !!FO!! z*z*r/k => unit [s]
936  ! if (Tstep>CFLval) then !CFL condition !!FO!! CFL condition: Courant�Friedrichs�Lewy condition is a necessary condition for convergence while solving
937  ! 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)
938  ! CFLfail=.TRUE.
939  ! endif
940  ! CFLval = minval(0.5*zroof*zroof*rroof/kroof)
941  ! if (Tstep>CFLval) then !CFL condition
942  ! write(*,*) "ROOF: CFL condition: Tstep=",Tstep,">",CFLval
943  ! CFLfail=.TRUE.
944  ! endif
945  ! CFLval = minval(0.5*zwall*zwall*rwall/kwall)
946  ! if (Tstep>CFLval) then !CFL condition
947  ! write(*,*) "WALL: CFL condition: Tstep=",Tstep,">",CFLval
948  ! CFLfail=.TRUE.
949  ! endif
950  ! CFLval = minval(0.5*zground*zground*rground/kground)
951  ! if (Tstep>CFLval) then !CFL condition
952  ! write(*,*) "ground: CFL condition: Tstep=",Tstep,">",CFLval
953  ! CFLfail=.TRUE.
954  ! endif
955  ! if (CFLfail) then
956  ! write(*,*) "Increase dX or decrease maxtimestep. Hit any key to continue"
957  ! read (*,*)
958  ! endif
959 
960  ! 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
961 
962  RETURN
963 
964  ! 315 CALL errorHint(11,TRIM(fileESTMTs),notUsed,notUsed,NotUsedI)
965  ! 316 CALL errorHint(11,TRIM(fileFinalTemp),notUsed,notUsed,NotUsedI)
966 
real(kind(1d0)) alb_avg
real(kind(1d0)) hw
integer, parameter cts_troad
real(kind(1d0)), parameter sbconst
real(kind(1d0)) rvf_wall
real(kind(1d0)) fwall
integer numberofgrids
real(kind(1d0)) ivf_fw
real(kind(1d0)), dimension(:, :), allocatable troof_grids
real(kind(1d0)) froof
real(kind(1d0)) alb_veg
real(kind(1d0)) fground
real(kind(1d0)), dimension(:, :, :), allocatable tw_4_grids
real(kind(1d0)) ivf_wf
real(kind(1d0)), parameter alb_wall
real(kind(1d0)) nroom
real(kind(1d0)) em_r
real(kind(1d0)) ivf_wr
real(kind(1d0)) tn_wall
integer, parameter cts_troof
real(kind(1d0)) ivf_ir
real(kind(1d0)) nan
integer, parameter bsoilsurf
integer, parameter cts_twall
real(kind(1d0)) zvf_wall
real(kind(1d0)) tfloor
real(kind(1d0)) lbc_soil
real(kind(1d0)), dimension(:), allocatable twall
real(kind(1d0)) lup_roof
real(kind(1d0)), dimension(:, :), allocatable ts5mindata
real(kind(1d0)), dimension(:), allocatable tibld
real(kind(1d0)) em_roof
real(kind(1d0)) ivf_rw
integer, parameter cts_tiair
real(kind(1d0)) fair
real(kind(1d0)) t0_ibld
real(kind(1d0)) rvf_roof
real(kind(1d0)) rvf_ground
real(kind(1d0)) ivf_fr
integer, parameter conifsurf
real(kind(1d0)) t0_ground
real(kind(1d0)) fveg
real(kind(1d0)), parameter em_wall
real(kind(1d0)) finternal
real(kind(1d0)) minshc_airbld
real(kind(1d0)), dimension(:), allocatable tair24hr
real(kind(1d0)), dimension(:, :), allocatable tground_grids
real(kind(1d0)) t0_wall
real(kind(1d0)) zvf_ground
real(kind(1d0)), dimension(:), allocatable tground
real(kind(1d0)), dimension(:, :), allocatable tibld_grids
real(kind(1d0)) em_veg
real(kind(1d0)) alb_ground
real(kind(1d0)), dimension(nsurf) emis
real(kind(1d0)) rvf_canyon
real(kind(1d0)) ivf_ri
integer, parameter grasssurf
real(kind(1d0)), dimension(:, :), allocatable tw_4
real(kind(1d0)) ivf_iw
real(kind(1d0)) svf_roof
integer, parameter ncolsestmdata
real(kind(1d0)) ivf_wi
real(kind(1d0)) em_i
real(kind(1d0)), dimension(:), allocatable troof
real(kind(1d0)) ivf_ww
real(kind(1d0)) rvf_veg
real(kind(1d0)) ivf_if
real(kind(1d0)) lup_wall
real(kind(1d0)) em_w
real(kind(1d0)) xvf_wall
real(kind(1d0)) em_ground
real(kind(1d0)) t0_roof
real(kind(1d0)) fibld
real(kind(1d0)) sumemis
real(kind(1d0)) tn_roof
real(kind(1d0)), dimension(nsurf) sfr
real(kind(1d0)), dimension(:, :, :), allocatable estmforcingdata
real(kind(1d0)) bldgh
real(kind(1d0)) tievolve
real(kind(1d0)) svf_wall
integer, parameter decidsurf
integer, parameter pavsurf
real(kind(1d0)) ivf_fi
real(kind(1d0)), dimension(nsurf) alb
real(kind(1d0)) sumalb
integer, parameter bldgsurf
integer, parameter watersurf
real(kind(1d0)) alb_roof
real(kind(1d0)) shc_air
real(kind(1d0)) em_ibld
real(kind(1d0)) svf_ground
real(kind(1d0)) lup_ground
real(kind(1d0)), parameter c2k
real(kind(1d0)) zref
real(kind(1d0)) ivf_ii
real(kind(1d0)) ivf_rf
integer nground
real(kind(1d0)) em_f
real(kind(1d0)), dimension(:, :), allocatable twall_grids
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 1609 of file suews_phys_estm.f95.

Referenced by estm().

1609  IMPLICIT NONE
1610  REAL(KIND(1d0)), PARAMETER::pnan = 9999
1611  REAL(KIND(1d0)), PARAMETER::nan = -999
1612  REAL(KIND(1d0)), INTENT(in)::x
1613  REAL(KIND(1d0))::xx
1614 
1615  IF (abs(x) > pnan) THEN
1616  xx = nan
1617  ELSE
1618  xx = x
1619  ENDIF
1620 
real(kind(1d0)) nan
Here is the caller graph for this function:

◆ suews_getestmdata()

subroutine estm_module::suews_getestmdata ( integer, intent(in)  lunit)

Definition at line 528 of file suews_phys_estm.f95.

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

529  USE data_in, ONLY: fileestmts, skipheadermet
530  USE sues_data, ONLY: tstep_real, tstep
531  USE defaultnotused, ONLY: notused, ios_out
533 
534  IMPLICIT NONE
535 
536  INTEGER, INTENT(in)::lunit
537  INTEGER::i, iyy !,RunNumber,NSHcounter
538  INTEGER :: iostat_var
539  REAL(KIND(1d0)), DIMENSION(ncolsESTMdata):: estmarray
540  REAL(KIND(1d0)):: imin_prev, ih_prev, iday_prev, tstep_estm !For checks on temporal resolution of estm data
541 
542  ! initialise
543  imin_prev = 0
544  ih_prev = 0
545  iday_prev = 0
546 
547  !---------------------------------------------------------------
548 
549  !Open the file for reading and read the actual data
550  !write(*,*) FileESTMTs
551  OPEN (lunit, file=trim(fileestmts), status='old', err=315)
552  CALL skipheader(lunit, skipheadermet)
553 
554  ! Skip to the right place in the ESTM file, depending on how many chunks have been read already
555  IF (skippedlines > 0) THEN
556  DO iyy = 1, skippedlines
557  READ (lunit, *)
558  ENDDO
559  ENDIF
560 
561  ! Read in next chunk of ESTM data and fill ESTMForcingData array with data for every timestep
562  DO i = 1, readlinesmetdata
563  READ (lunit, *, iostat=iostat_var) estmarray
564  estmforcingdata(i, 1:ncolsestmdata, gridcounter) = estmarray
565  ! Check timestamp of met data file matches TSTEP specified in RunControl
566  IF (i == 1) THEN
567  imin_prev = estmarray(4)
568  ih_prev = estmarray(3)
569  iday_prev = estmarray(2)
570  ELSEIF (i == 2) THEN
571  tstep_estm = ((estmarray(4) + 60*estmarray(3)) - (imin_prev + 60*ih_prev))*60 !tstep in seconds
572  IF (tstep_estm /= tstep_real .AND. estmarray(2) == iday_prev) THEN
573  CALL errorhint(39, 'TSTEP in RunControl does not match TSTEP of ESTM data (DOY).', &
574  REAL(tstep, KIND(1d0)), tstep_estm, int(estmarray(2)))
575  ENDIF
576  ENDIF
577  ENDDO
578 
579  CLOSE (lunit)
580 
581  RETURN
582 
583 315 CALL errorhint(11, trim(fileestmts), notused, notused, ios_out)
584 
real(kind(1d0)) notused
subroutine skipheader(lfn, skip)
character(len=150) fileestmts
integer skipheadermet
integer skippedlines
real(kind(1d0)) tstep_real
integer, parameter ncolsestmdata
real(kind(1d0)), dimension(:, :, :), allocatable estmforcingdata
integer readlinesmetdata
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
integer gridcounter
Here is the call graph for this function:
Here is the caller graph for this function: