SUEWS API Site
Documentation of SUEWS source code
suews_program.f95
Go to the documentation of this file.
1 !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
2 !Main program of SUEWS version 1.0 -------------------------------------------
3 ! Model calculations are made in two stages:
4 ! (1) initialise the run for each block of met data (iblock from 1 to ReadBlocksMetData)
5 ! (2) perform the actual model calculations (SUEWS_Calculations)
6 ! After reading in all input information from SiteInfo, the code loops
7 ! - over years
8 ! - then over blocks (met data read in and stored for each block)
9 ! - then over rows
10 ! - then over grids
11 !Last modified by TS 10 Apr 2017 - conditional compilation blocks added for netCDF adaptation
12 !Last modified by LJ 6 Apr 2017 - Snow initialisation, allocation and deallocation added
13 !Last modified by HCW 10 Feb 2017 - Disaggregation of met forcing data
14 !Last modified by HCW 12 Jan 2017 - Changes to InitialConditions
15 !Last modified by HCW 26 Aug 2016 - CO2 flux added
16 !Last modified by HCW 04 Jul 2016 - GridID can now be up to 10 digits long
17 !Last modified by HCW 29 Jun 2016 - Reversed over-ruling of ReadLinesMetdata so this is not restricted here to one day
18 !Last modified by HCW 27 Jun 2016 - Re-corrected grid number for output files. N.B. Gridiv seems to have been renamed iGrid
19 ! - Met file no longer has grid number attached if same met data used for all grids
20 !Last modified by HCW 24 May 2016 - InitialConditions file naming altered
21 ! Unused year_txt argument removed from InitialState
22 ! LJ 30 Mar 2016 - Grid run order changed from linear to non-linear
23 !Last modified by TS 14 Mar 2016 - Include AnOHM daily iteration
24 !Last modified by HCW 25 Jun 2015 - Fixed bug in LAI calculation at year change
25 !Last modified by HCW 12 Mar 2015
26 !Last modified by HCW 26 Feb 2015
27 !Last modified by HCW 03 Dec 2014
28 !
29 ! To do:
30 ! - Water movement between grids (GridConnections) not yet coded
31 !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
33 
34  USE allocatearray
36  USE data_in
37  USE defaultnotused
38  USE initial
39  USE metdisagg
40  USE sues_data
41  USE time
42  USE wherewhen
43  USE ctrl_output
46 
47  IMPLICIT NONE
48 
49  CHARACTER(len=4) :: year_txt !Year as a text string
50 
51  CHARACTER(len=20):: FileCodeX, & !Current file code
52  FileCodeXWY, & !File code without year
53  FileCodeXWG !File code without grid
54  CHARACTER(len=20):: grid_txt, & !Grid number as a text string (from FirstGrid to LastGrid)
55  tstep_txt, & !Model timestep (in minutes) as a text string (minutes)
56  ResIn_txt, ResInESTM_txt !Resolution of Original met/ESTM forcing file as a text string (minutes)
57 
58  INTEGER:: nlinesLimit, & !Max number of lines that can be read in one go for each grid
59  NumberOfYears !Number of years to be run
60 
61  INTEGER:: year_int, & ! Year as an integer (from SiteSelect rather than met forcing file)
62  igrid, & !Grid number (from 1 to NumberOfGrids)
63  iblock, & !Block number (from 1 to ReadBlocksMetData)
64  ir, irMax, & !Row number within each block (from 1 to irMax)
65  rr !Row of SiteSelect corresponding to current year and grid
66 
67  REAL:: timeStart, timeFinish ! profiling use, AnOHM TS
68 
69  !==========================================================================
70 
71  ! Start counting cpu time
72  CALL cpu_time(timestart)
73 
74  ! initialise simulation time
75  dt_since_start = 0
76 
77  WRITE (*, *) '========================================================'
78  WRITE (*, *) 'Running ', progname
79 
80  ! Initialise error file (0 -> problems.txt file will be newly created)
81  errorchoice = 0
82  ! Initialise error file (0 -> warnings.txt file will be newly created)
83  warningchoice = 0
84  ! Initialise OutputFormats to 1 so that output format is written out only once per run
85  outputformats = 1
86 
87  ! Initialise WhereWhen variables for error handling
88  gridid_text = '00000'
89  datetime = '00000'
90 
91  ! Read RunControl.nml and all .txt input files from SiteSelect spreadsheet
93 
94  WRITE (tstep_txt, '(I5)') tstep/60 !Get tstep (in minutes) as a text string
95  WRITE (resin_txt, '(I5)') resolutionfilesin/60 !Get ResolutionFilesIn (in minutes) as a text string
96  WRITE (resinestm_txt, '(I5)') resolutionfilesinestm/60
97 
98  ! Find first and last year of the current run
99  firstyear = minval(int(siteselect(:, c_year)))
100  lastyear = maxval(int(siteselect(:, c_year)))
101 
102  numberofyears = lastyear - firstyear + 1 !Find the number of years to run
103 
104  !Find the the number of grids within each year in SUEWS_SiteSelect.txt
105  ! N.B. need to have the same grids for each year
106  numberofgrids = int(nlinessiteselect/numberofyears)
107 
108  !! Find the first and last grid numbers (N.B. need to have the same grids for each year)
109  !FirstGrid = minval(int(SiteSelect(:,c_Grid)))
110  !LastGrid = maxval(int(SiteSelect(:,c_Grid)))
111  IF (numberofgrids > maxnumberofgrids) THEN
112  CALL errorhint(64, &
113  'No. of grids exceeds max. possible no. of grids.', &
114  REAL(MaxNumberOfGrids, KIND(1d0)), NotUsed, NumberOfGrids)
115  ENDIF
116 
117  ALLOCATE (grididmatrix(numberofgrids)) !Get the nGrid numbers correctly
118  ALLOCATE (grididmatrix0(numberofgrids)) !Get the nGrid numbers correctly
119 
120  DO igrid = 1, numberofgrids
121  grididmatrix(igrid) = int(siteselect(igrid, c_grid))
122  ENDDO
123 
124 ! #ifdef nc
125 ! ! sort grid matrix to conform the geospatial layout as in QGIS, TS 14 Dec 2016
126 ! IF (ncMode == 1) THEN
127 ! GridIDmatrix0 = GridIDmatrix
128 ! CALL sortGrid(GridIDmatrix0, GridIDmatrix, nRow, nCol)
129 ! ENDIF
130 ! ! GridIDmatrix0 stores the grid ID in the original order
131 ! #endif
132 
133  ! GridIDmatrix=GridIDmatrix0
134  WRITE (*, *) '--------------------------------------------'
135  WRITE (*, *) 'Years identified:', firstyear, 'to', lastyear
136  WRITE (*, *) 'No. grids identified:', numberofgrids, 'grids'
137  WRITE (*, *) 'Maximum No. grids allowed:', maxnumberofgrids, 'grids'
138 
139  ! Set limit on number of lines to read
140  nlineslimit = int(floor(maxlinesmet/REAL(numberofgrids, kind(1d0)))) !Uncommented HCW 29 Jun 2016
141  !nlinesLimit = 24*nsh !Commented out HCW 29 Jun 2016
142 
143  ! ---- Allocate arrays ----------------------------------------------------
144  ! Daily state needs to be outside year loop to transfer states between years
145  ALLOCATE (modeldailystate(numberofgrids, maxncols_cmds)) !DailyState
146  ALLOCATE (dailystatefirstopen(numberofgrids)) !Initialisation for header
147 
148  ! ---- Initialise arrays --------------------------------------------------
149  modeldailystate(:, :) = -999
150  dailystatefirstopen(:) = 1
151 
152  ! -------------------------------------------------------------------------
153  ! Initialise ESTM (reads ESTM nml, should only run once)
154  IF (storageheatmethod == 4 .OR. storageheatmethod == 14) THEN
155  IF (diagnose == 1) WRITE (*, *) 'Calling ESTM_initials...'
156  CALL estm_initials
157  ENDIF
158 
159  !==========================================================================
160  DO year_int = firstyear, lastyear !Loop through years
161 
162  WRITE (*, *) ' '
163  WRITE (year_txt, '(I4)') year_int !Get year as a text string
164 
165  ! Find number of days in the current year
166  CALL leapyearcalc(year_int, nofdaysthisyear)
167 
168  ! Prepare to disaggregate met data to model time-step (if required) ------
169  ! Find number of model time-steps per resolution of original met forcing file
170  nper_real = resolutionfilesin/REAL(tstep, kind(1d0))
171  nper = int(nper_real)
172 
173  IF (nper /= nper_real) THEN
174  CALL errorhint(2, 'Problem in SUEWS_Program: check resolution of met forcing data (ResolutionFilesIn)'// &
175  'and model time-step (Tstep).', &
176  REAL(Tstep, KIND(1d0)), NotUsed, ResolutionFilesIn)
177  ELSEIF (nper > 1) THEN
178  WRITE (*, *) 'Resolution of met forcing data: ', trim(adjustl(resin_txt)), ' min;', &
179  ' model time-step: ', trim(adjustl(tstep_txt)), ' min', ' -> SUEWS will perform disaggregation.'
180  IF (diagnose == 1) WRITE (*, *) 'Getting information for met disaggregation'
181  ! Get names of original met forcing file(s) to disaggregate (using first grid)
182  WRITE (grid_txt, '(I10)') grididmatrix(1) !Get grid as a text string
183 
184  ! Get met file name for this grid: SSss_YYYY_data_RR.txt
185  fileorigmet = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)//'_data_' &
186  //trim(adjustl(resin_txt))//'.txt'
187  ! But if each grid has the same met file, met file name does not include grid number
188  IF (multiplemetfiles /= 1) THEN
189  fileorigmet = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_data_' &
190  //trim(adjustl(resin_txt))//'.txt'
191  ENDIF
192 
193  nlinesorigmetdata = 0 !Initialise nlinesMetdata (total number of lines in met forcing file)
195  ! WRITE(*,*) 'nlinesOrigMetdata', nlinesOrigMetdata
196 
197  readlinesorigmetdata = nlinesorigmetdata !Initially set limit as the size of file
198  IF (nlinesorigmetdata*nper > nlineslimit) THEN !But restrict if this limit exceeds memory capacity
199  readlinesorigmetdata = int(nlineslimit/nper)
200  ENDIF
201  ! make sure the metblocks read in consists of complete diurnal cycles
202  nsdorig = nsd/nper
204  !WRITE(*,*) 'ReadlinesOrigMetdata', ReadlinesOrigMetdata
205  WRITE (*, *) 'Original met data will be read in chunks of ', readlinesorigmetdata, 'lines.'
206 
207  readblocksorigmetdata = int(ceiling(REAL(nlinesorigmetdata, kind(1d0))/REAL(readlinesorigmetdata, kind(1d0))))
208 
209  ! Set ReadLinesMetdata and ReadBlocksMetData
211  readblocksmetdata = int(ceiling(REAL(nlinesorigmetdata*nper, kind(1d0))/REAL(readlinesmetdata, kind(1d0))))
212  WRITE (*, *) 'Processing current year in ', readblocksmetdata, 'blocks.'
213 
215 
216  ELSEIF (nper == 1) THEN
217  WRITE (*, *) 'ResolutionFilesIn = Tstep: no disaggregation needed for met data.'
218 
219  !-----------------------------------------------------------------------
220  ! Find number of lines in met forcing file for current year (nlinesMetdata)
221  ! Need to know how many lines will be read each iteration
222  ! Use first grid as an example as the number of lines is the same for all grids
223  ! within one year
224  WRITE (grid_txt, '(I10)') grididmatrix(1) !Get grid as a text string
225 
226  ! Get met file name for this year for this grid
227  filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
228  filemet = trim(fileinputpath)//trim(filecodex)//'_data_'//trim(adjustl(tstep_txt))//'.txt'
229  !If each grid has the same met file, met file name does not include grid number (HCW 27 Jun 2016)
230  IF (multiplemetfiles /= 1) THEN
231  filecodexwg = trim(filecode)//'_'//trim(year_txt) !File code without grid
232  filemet = trim(fileinputpath)//trim(filecodexwg)//'_data_'//trim(adjustl(tstep_txt))//'.txt'
233  ENDIF
234 
235  nlinesmetdata = 0 !Initialise nlinesMetdata (total number of lines in met forcing file)
237  !-----------------------------------------------------------------------
238 
239  ! To conserve memory, read met data in blocks
240  ! Find number of lines that can be read in each block (i.e. read in at once)
241  readlinesmetdata = nlinesmetdata !Initially set limit as the size of the met file (N.B.solves problem with Intel fortran)
242  IF (nlinesmetdata > nlineslimit) THEN !But restrict if this limit exceeds memory capacity
243  readlinesmetdata = nlineslimit
244  ENDIF
245  ! make sure the metblocks read in consists of complete diurnal cycles, TS 08 Jul 2016
247 
248  WRITE (*, *) 'Met data will be read in blocks of ', readlinesmetdata, 'lines.'
249 
250  ! Find number of blocks of met data
251  readblocksmetdata = int(ceiling(REAL(nlinesmetdata, kind(1d0))/REAL(readlinesmetdata, kind(1d0))))
252  WRITE (*, *) 'Processing current year in ', readblocksmetdata, 'blocks.'
253 
254  ENDIF
255 
256  !write(*,*) ReadBlocksMetData, ReadBlocksOrigMetData
257 
258  ! ---- Allocate arrays--------------------------------------------------
259  IF (diagnose == 1) WRITE (*, *) 'Allocating arrays in SUEWS_Program.f95...'
260  ALLOCATE (surfacechar(numberofgrids, maxncols_c)) !Surface characteristics
261  ALLOCATE (metforcingdata(readlinesmetdata, ncolumnsmetforcingdata, numberofgrids)) !Met forcing data
262  ALLOCATE (metforcingdata_grid(readlinesmetdata, ncolumnsmetforcingdata)) !Met forcing data
263  ALLOCATE (modeloutputdata(0:readlinesmetdata, maxncols_cmod, numberofgrids)) !Data at model timestep
264  ALLOCATE (dataoutsuews(readlinesmetdata, ncolumnsdataoutsuews, numberofgrids)) !Main output array
265  dataoutsuews = nan ! initialise Main output array
266  ALLOCATE (dataoutrsl(readlinesmetdata, ncolumnsdataoutrsl, numberofgrids)) !RSL output array
267  dataoutrsl = nan ! initialise Main output array
268  ALLOCATE (dataoutdailystate(ndays, ncolumnsdataoutdailystate, numberofgrids)) !DailyState array
269  dataoutdailystate = nan ! initialise DailyState
270  ! IF (SOLWEIGuse == 1) ALLOCATE (dataOutSOLWEIG(ReadLinesMetdata, ncolumnsdataOutSOL, NumberOfGrids)) !SOLWEIG POI output
271  ALLOCATE (dataoutsolweig(readlinesmetdata, ncolumnsdataoutsol, numberofgrids)) !SOLWEIG POI output
273  IF (cbluse >= 1) ALLOCATE (dataoutbl(readlinesmetdata, ncolumnsdataoutbl, numberofgrids)) !CBL output
274  ! IF (SnowUse == 1) THEN
275  IF (.NOT. ALLOCATED(dataoutsnow)) ALLOCATE (dataoutsnow(readlinesmetdata, ncolumnsdataoutsnow, numberofgrids)) !Snow output
276  ! ALLOCATE(qn1_S_store(NSH,NumberOfGrids))
277  ! ALLOCATE(qn1_S_store_grid(NSH))
278  ! ALLOCATE(qn1_S_av_store(2*NSH+1,NumberOfGrids))
279  ! ALLOCATE(qn1_S_av_store_grid(2*NSH+1))
280  IF (.NOT. ALLOCATED(qn1_s_av_grids)) ALLOCATE (qn1_s_av_grids(numberofgrids))
281  IF (.NOT. ALLOCATED(dqnsdt_grids)) ALLOCATE (dqnsdt_grids(numberofgrids))
282  ! qn1_S_store(:,:) = NAN
283  ! qn1_S_av_store(:,:) = NaN
284  ! qn1_S_store_grid(:) = NAN
285  ! qn1_S_av_store_grid(:) = NaN
286  qn1_s_av_grids = 0 ! Initialise to 0
287  dqnsdt_grids = 0 ! Initialise to 0
288  ! ENDIF
289  ! IF (StorageHeatMethod==4 .OR. StorageHeatMethod==14) THEN
290  IF (.NOT. ALLOCATED(dataoutestm)) ALLOCATE (dataoutestm(readlinesmetdata, 32, numberofgrids)) !ESTM output
291  ! ENDIF
292  ! ALLOCATE(TstepProfiles(NumberOfGrids,14,24*NSH)) !Hourly profiles interpolated to model timestep
293  ! ALLOCATE(AHProf_tstep(24*NSH,2)) !Anthropogenic heat profiles at model timestep
294  ! ALLOCATE(WUProfM_tstep(24*NSH,2)) !Manual water use profiles at model timestep
295  ! ALLOCATE(WUProfA_tstep(24*NSH,2)) !Automatic water use profiles at model timestep
296  ! ALLOCATE(HumActivity_tstep(24*NSH,2)) !Human activity profiles at model timestep
297  ! ALLOCATE(TraffProf_tstep(24*NSH,2)) !Traffic profiles at model timestep
298  ! ALLOCATE(PopProf_tstep(24*NSH,2)) !Population profiles at model timestep
299  ! ALLOCATE(qn1_store(NSH,NumberOfGrids))
300  ! ALLOCATE(qn1_av_store(2*NSH+1,NumberOfGrids))
301  ! ALLOCATE(qn1_store_grid(NSH))
302  ! ALLOCATE(qn1_av_store_grid(2*NSH+1))
303  ALLOCATE (qhforcbl(numberofgrids))
304  ALLOCATE (qeforcbl(numberofgrids))
305 
306  ALLOCATE (tair_av_grids(numberofgrids))
307  ALLOCATE (qn1_av_grids(numberofgrids))
308  ALLOCATE (dqndt_grids(numberofgrids))
309  !! QUESTION: Add snow clearing (?)
310 
311  ! qn1_store(:,:) = NAN ! Initialise to -999
312  ! qn1_av_store(:,:) = NAN ! Initialise to -999
313  ! qn1_store_grid(:)=NAN
314  ! qn1_av_store_grid(:)=NAN
315  tair_av_grids = 273.15 ! Initialise to 273.15 K
316  qn1_av_grids = 0 ! Initialise to 0
317  dqndt_grids = 0 ! Initialise to 0
318 
319  qhforcbl(:) = nan
320  qeforcbl(:) = nan
321 
322  ! QUESTION: Initialise other arrays here?
323 
324  IF (storageheatmethod == 4 .OR. storageheatmethod == 14) THEN
325  ! Prepare to disaggregate ESTM data to model time-step (if required) ------
326  ! Find number of model time-steps per resolution of original met forcing file
327  nperestm_real = resolutionfilesinestm/REAL(tstep, kind(1d0))
328  nperestm = int(nperestm_real)
329  IF (nperestm /= nperestm_real) THEN
330  CALL errorhint(2, 'Problem in SUEWS_Program: check resolution of ESTM forcing data (ResolutionFilesInESTM)'// &
331  'and model time-step (Tstep).', &
332  REAL(Tstep, KIND(1d0)), NotUsed, ResolutionFilesInESTM)
333  ELSEIF (nperestm > 1) THEN
334  WRITE (*, *) 'Resolution of ESTM forcing data: ', trim(adjustl(resinestm_txt)), ' min;', &
335  ' model time-step: ', trim(adjustl(tstep_txt)), ' min', ' -> SUEWS will perform disaggregation.'
336  IF (diagnose == 1) WRITE (*, *) 'Getting information for ESTM disaggregation'
337  ! Get names of original met forcing file(s) to disaggregate (using first grid)
338  WRITE (grid_txt, '(I10)') grididmatrix(1) !Get grid as a text string
339 
340  ! Get met file name for this grid
341  fileestmts = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)//'_ESTM_Ts_data_' &
342  //trim(adjustl(resinestm_txt))//'.txt'
343  ! But if each grid has the same ESTM file, file name does not include grid number
344  IF (multipleestmfiles /= 1) THEN
345  fileestmts = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_ESTM_Ts_data_' &
346  //trim(adjustl(resinestm_txt))//'.txt'
347  ENDIF
348 
349  ! Find number of lines in orig ESTM file
350  nlinesorigestmdata = 0 !Initialise nlinesMetdata (total number of lines in met forcing file)
352 
353  ! Check ESTM data and met data will have the same length (so that ESTM file can be read in same blocks as met data)
355  CALL errorhint(66, &
356  'Downscaled ESTM and met input files will have different lengths', REAL(nlinesMetdata, KIND(1d0)), &
357  NotUsed, nlinesESTMdata*NperESTM)
358  ENDIF
359 
360  !write(*,*) 'nlinesOrigESTMdata', nlinesOrigESTMdata
361  ! Set number of lines to read from original ESTM file using met data blocks
363  !WRITE(*,*) 'ReadlinesOrigESTMdata', ReadlinesOrigESTMdata
364  WRITE (*, *) 'Original ESTM data will be read in chunks of ', readlinesorigestmdata, 'lines.'
365 
366  nlinesestmdata = nlinesorigestmdata*nperestm
367 
368  ELSEIF (nperestm == 1) THEN
369  WRITE (*, *) 'ResolutionFilesInESTM = Tstep: no disaggregation needed for met data.'
370 
371  !-----------------------------------------------------------------------
372  ! Find number of lines in ESTM forcing file for current year (nlinesESTMdata)
373  WRITE (grid_txt, '(I10)') grididmatrix(1) !Get grid as a text string (use first grid as example)
374  ! Get file name for this year for this grid
375  filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
376  fileestmts = trim(fileinputpath)//trim(filecodex)//'_ESTM_Ts_data_'//trim(adjustl(tstep_txt))//'.txt'
377  !If each grid has the same met file, met file name does not include grid number
378  IF (multipleestmfiles /= 1) THEN
379  filecodexwg = trim(filecode)//'_'//trim(year_txt) !File code without grid
380  fileestmts = trim(fileinputpath)//trim(filecodexwg)//'_ESTM_Ts_data_'//trim(adjustl(tstep_txt))//'.txt'
381  ENDIF
382 
383  ! Find number of lines in ESTM file
384  nlinesestmdata = 0 !Initialise nlinesMetdata (total number of lines in met forcing file)
385  nlinesestmdata = count_lines(trim(fileestmts))
386  !-----------------------------------------------------------------------
387 
388  ! Check ESTM data and met data are same length (so that ESTM file can be read in same blocks as met data)
389  IF (nlinesestmdata /= nlinesmetdata) THEN
390  CALL errorhint(66, 'ESTM input file different length to met forcing file', REAL(nlinesMetdata, KIND(1d0)), &
391  NotUsed, nlinesESTMdata)
392  ENDIF
393 
394  ENDIF
395 
396  ! Allocate arrays to receive ESTM forcing data
397  ! IF (.NOT. ALLOCATED(ESTMForcingData)) ALLOCATE(ESTMForcingData(1:ReadLinesMetdata,ncolsESTMdata,NumberOfGrids))
398  ! IF (.NOT. ALLOCATED(Ts5mindata)) ALLOCATE(Ts5mindata(1:ReadLinesMetdata,ncolsESTMdata))
399  ! IF (.NOT. ALLOCATED(Ts5mindata_ir)) ALLOCATE(Ts5mindata_ir(ncolsESTMdata))
400  !
401  ! IF (.NOT. ALLOCATED(Tair24HR)) ALLOCATE(Tair24HR(24*nsh))
402  ! if ( /= 0) print *, ": Deallocation request denied"
403  ! ALLOCATE(Tair24HR(24*nsh))
404 
405  ENDIF
406  IF (.NOT. ALLOCATED(estmforcingdata)) ALLOCATE (estmforcingdata(1:readlinesmetdata, ncolsestmdata, numberofgrids))
407  IF (.NOT. ALLOCATED(ts5mindata)) ALLOCATE (ts5mindata(1:readlinesmetdata, ncolsestmdata))
408  IF (.NOT. ALLOCATED(ts5mindata_ir)) ALLOCATE (ts5mindata_ir(ncolsestmdata))
409  IF (.NOT. ALLOCATED(tair24hr)) ALLOCATE (tair24hr(24*nsh))
410  ! ------------------------------------------------------------------------
411 
412  !-----------------------------------------------------------------------
413  skippedlines = 0 !Initialise lines to be skipped in met forcing file
414  skippedlinesorig = 0 !Initialise lines to be skipped in original met forcing file
415  skippedlinesorigestm = 0 !Initialise lines to be skipped in original met forcing file
416 
417  DO iblock = 1, readblocksmetdata !Loop through blocks of met data
418  ! WRITE(*,*) iblock,'/',ReadBlocksMetData
419 
420  ! Model calculations are made in two stages:
421  ! (1) initialise the run for each block of met data (iblock from 1 to ReadBlocksMetData)
422  ! (2) perform the actual model calculations (SUEWS_Calculations)
423 
424  gridcounter = 1 !Initialise counter for grids in each year
425  DO igrid = 1, numberofgrids !Loop through grids
426 
427  gridid = grididmatrix(igrid) !store grid here for referencing error codes
428  WRITE (grid_txt, '(I10)') grididmatrix(igrid) !Get grid ID as a text string
429 
430  ! (1) First stage: initialise run if this is the first iteration this year
431  ! (1a) Initialise surface characteristics
432  IF (iblock == 1) THEN
433  IF (diagnose == 1) WRITE (*, *) 'First block of data - doing initialisation'
434  ! (a) Transfer characteristics from SiteSelect to correct row of SurfaceChar
435  DO rr = 1, nlinessiteselect
436  !Find correct grid and year
437  IF (diagnose == 1) WRITE (*, *) 'grid found:', siteselect(rr, c_grid), 'grid needed:', grididmatrix(igrid)
438  IF (diagnose == 1) WRITE (*, *) 'year found:', siteselect(rr, c_year), 'year needed:', year_int
439  IF (siteselect(rr, c_grid) == grididmatrix(igrid) .AND. siteselect(rr, c_year) == year_int) THEN
440  IF (diagnose == 1) WRITE (*, *) 'Match found (grid and year) for rr = ', rr, 'of', nlinessiteselect
442  EXIT
443  ELSEIF (rr == nlinessiteselect) THEN
444  WRITE (*, *) 'Program stopped! Year', year_int, 'and/or grid', igrid, 'not found in SiteSelect.txt.'
445  CALL errorhint(59, 'Cannot find year and/or grid in SiteSelect.txt', REAL(igrid, KIND(1d0)), NotUsed, year_int)
446  ENDIF
447  ENDDO
448  ENDIF !end first block of met data
449 
450  ! (1b) Initialise met data
451  IF (nper > 1) THEN
452  ! Disaggregate met data ---------------------------------------------------
453 
454  ! Set maximum value for ReadLinesOrigMetData to handle end of file (i.e. small final block)
455  IF (iblock == readblocksmetdata) THEN !For last block of data in file
457  ELSE
459  ENDIF
460  !write(*,*) ReadLinesOrigMetDataMAX, ReadLinesOrigMetData
461  ! Get names of original met forcing file(s) to disaggregate
462  ! Get met file name for this grid: SSss_YYYY_data_RR.txt
463  IF (multiplemetfiles == 1) THEN !If each grid has its own met file
464  fileorigmet = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)//'_data_' &
465  //trim(adjustl(resin_txt))//'.txt'
466  ! Also set file name for downscaled file
467  filedscdmet = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)//'_data_' &
468  //trim(adjustl(tstep_txt))//'.txt'
469  ! Disaggregate met data
470  CALL disaggregatemet(iblock, igrid)
471  ELSE
472  ! If each grid has the same met file, met file name does not include grid number, and only need to disaggregate once
473  fileorigmet = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_data_' &
474  //trim(adjustl(resin_txt))//'.txt'
475  filedscdmet = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_data_' &
476  //trim(adjustl(tstep_txt))//'.txt'
477  IF (igrid == 1) THEN !Disaggregate for the first grid only
478  CALL disaggregatemet(iblock, igrid)
479  ELSE !Then for subsequent grids simply copy data
481  ENDIF
482  ENDIF
483 
484  ELSEIF (nper == 1) THEN
485  ! Get met forcing file name for this year for the first grid
486  ! Can be something else than 1
487  filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
488  filecodexwg = trim(filecode)//'_'//trim(year_txt) !File code without grid
489  ! IF(iblock==1) WRITE(*,*) 'Current FileCode: ', FileCodeX
490 
491  ! For every block of met data ------------------------------------
492  ! Initialise met forcing data into 3-dimensional matrix
493  !write(*,*) 'Initialising met data for block',iblock
494  IF (multiplemetfiles == 1) THEN !If each grid has its own met file
495  filemet = trim(fileinputpath)//trim(filecodex)//'_data_'//trim(adjustl(tstep_txt))//'.txt'
497  ELSE !If one met file used for all grids
498  !FileMet=TRIM(FileInputPath)//TRIM(FileCodeX)//'_data_'//TRIM(ADJUSTL(tstep_txt))//'.txt'
499  ! If one met file used for all grids, look for met file with no grid code (FileCodeXWG)
500  filemet = trim(fileinputpath)//trim(filecodexwg)//'_data_'//trim(adjustl(tstep_txt))//'.txt'
501  IF (igrid == 1) THEN !Read for the first grid only
503  ELSE !Then for subsequent grids simply copy data
505  ENDIF
506  ENDIF
507  ENDIF !end of nper statement
508 
509  ! Only for the first block of met data, read initial conditions (moved from above, HCW 12 Jan 2017)
510  IF (iblock == 1) THEN
511  filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
512  !write(*,*) ' Now calling InitialState'
513  CALL initialstate(filecodex, year_int, gridcounter, numberofgrids)
514  ENDIF
515 
516  ! Initialise ESTM if required, TS 05 Jun 2016; moved inside grid loop HCW 27 Jun 2016
517  IF (storageheatmethod == 4 .OR. storageheatmethod == 14) THEN
518  IF (nperestm > 1) THEN
519  ! Disaggregate ESTM data --------------------------------------------------
520  ! Set maximum value for ReadLinesOrigESTMData to handle end of file (i.e. small final block)
521  IF (iblock == readblocksmetdata) THEN !For last block of data in file
523  ELSE
525  ENDIF
526  !write(*,*) ReadLinesOrigESTMDataMAX, ReadLinesOrigESTMData
527  ! Get names of original ESTM forcing file(s) to disaggregate
528  ! Get ESTM file name for this grid: SSss_YYYY_ESTM_Ts_data_RR.txt
529  IF (multipleestmfiles == 1) THEN !If each grid has its own ESTM file
530  fileorigestm = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt) &
531  //'_ESTM_Ts_data_'//trim(adjustl(resinestm_txt))//'.txt'
532  ! Also set file name for downscaled file
533  filedscdestm = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt) &
534  //'_ESTM_Ts_data_'//trim(adjustl(tstep_txt))//'.txt'
535  ! Disaggregate ESTM data
536  CALL disaggregateestm(iblock)
537  ELSE
538  ! If each grid has the same ESTM file, ESTM file name does not include grid number, and only need to disaggregate once
539  fileorigestm = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_ESTM_Ts_data_' &
540  //trim(adjustl(resinestm_txt))//'.txt'
541  filedscdestm = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_ESTM_Ts_data_' &
542  //trim(adjustl(tstep_txt))//'.txt'
543  IF (igrid == 1) THEN !Disaggregate for the first grid only
544  CALL disaggregateestm(iblock)
545  ELSE !Then for subsequent grids simply copy data
547  1:ncolsestmdata, 1)
548  ENDIF
549  ENDIF
550 
551  ELSEIF (nperestm == 1) THEN
552  ! Get ESTM forcing file name for this year for the first grid
553  filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
554  filecodexwg = trim(filecode)//'_'//trim(year_txt) !File code without grid
555  ! For every block of ESTM data ------------------------------------
556  ! Initialise ESTM forcing data into 3-dimensional matrix
557  !write(*,*) 'Initialising ESTM data for block',iblock
558  IF (multipleestmfiles == 1) THEN !If each grid has its own met file
559  fileestmts = trim(fileinputpath)//trim(filecodex)//'_ESTM_Ts_data_'//trim(adjustl(tstep_txt))//'.txt'
560  !write(*,*) 'Calling GetESTMData...', FileCodeX, iblock, igrid
561  CALL suews_getestmdata(101)
562  ELSE !If one ESTM file used for all grids
563  fileestmts = trim(fileinputpath)//trim(filecodexwg)//'_ESTM_Ts_data_'//trim(adjustl(tstep_txt))//'.txt'
564  !write(*,*) 'Calling GetESTMData...', FileCodeX, iblock, igrid
565  IF (igrid == 1) THEN !Read for the first grid only
566  CALL suews_getestmdata(101)
567  ELSE !Then for subsequent grids simply copy data
569  1:ncolsestmdata, 1)
570  ENDIF
571  ENDIF
572  ENDIF !end of nperESTM statement
573  ENDIF
574 
575  gridcounter = gridcounter + 1 !Increase GridCounter by 1 for next grid
576 
577  ENDDO !end loop over grids
578  skippedlines = skippedlines + readlinesmetdata !Increase skippedLines ready for next block
579  skippedlinesorig = skippedlinesorig + readlinesorigmetdata !Increase skippedLinesOrig ready for next block
580  skippedlinesorigestm = skippedlinesorigestm + readlinesorigestmdata !Increase skippedLinesOrig ready for next block
581  !write(*,*) iblock
582  !write(*,*) ReadLinesMetdata, readlinesorigmetdata
583  !write(*,*) skippedLines, skippedLinesOrig, skippedLinesOrig*Nper
584 
585  ! Initialise the modules on the first day
586  ! Initialise CBL and SOLWEIG parts if required
587  IF (iblock == 1) THEN
588  IF ((cbluse == 1) .OR. (cbluse == 2)) CALL cbl_readinputdata(fileinputpath, qh_choice)
589  ENDIF
590 
591  ! NB: SOLWEIG is disabled in v2018a
592  ! QUESTION: not sure if solweig should be within iblock if statement too?
593  ! IF(SOLWEIGuse==1) CALL SOLWEIG_initial
594 
595  !write(*,*) 'Initialisation done'
596  ! First stage: initialisation done ----------------------------------
597 
598  ! (2) Second stage: do calculations at 5-min time-steps -------------
599  ! First set maximum value of ir
600  IF (iblock == readblocksmetdata) THEN !For last block of data in file
601  irmax = nlinesmetdata - (iblock - 1)*readlinesmetdata
602  ELSE
603  irmax = readlinesmetdata
604  ENDIF
605 
606  DO ir = 1, irmax !Loop through rows of current block of met data
607  gridcounter = 1 !Initialise counter for grids in each year
608 
609  DO igrid = 1, numberofgrids !Loop through grids
610  IF (diagnose == 1) WRITE (*, *) 'Row (ir):', ir, '/', irmax, 'of block (iblock):', iblock, '/', readblocksmetdata, &
611  'Grid:', grididmatrix(igrid)
612 
613  ! Call model calculation code
614  ! IF(ir==1) WRITE(*,*) 'Now running block ',iblock,'/',ReadBlocksMetData,' of year ',year_int,'...'
615  WRITE (grid_txt, '(I10)') grididmatrix(igrid) !Get grid ID as a text string
616  filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
617  IF (ir == 1 .AND. igrid == 1) THEN
618  WRITE (*, *) trim(adjustl(filecodex)), &
619  ': Now running block ', iblock, '/', readblocksmetdata, ' of ', trim(year_txt), '...'
620  ENDIF
621  IF (diagnose == 1) WRITE (*, *) 'Calling SUEWS_Calculations...'
622  ! print*, 'before cal:',sum(Tair24HR)
623  CALL suews_calculations(gridcounter, ir, iblock, irmax)
624  ! print*, 'after cal:',sum(Tair24HR)
625  IF (diagnose == 1) WRITE (*, *) 'SUEWS_Calculations finished...'
626 
627  ! Record iy and id for current time step to handle last row in yearly files (YYYY 1 0 0)
628  ! IF(GridCounter == NumberOfGrids) THEN !Adjust only when the final grid has been run for this time step
629  IF (igrid == numberofgrids) THEN !Adjust only when the final grid has been run for this time step
630  iy_prev_t = iy
631  id_prev_t = id
632  ENDIF
633 
634  ! Write state information to new InitialConditions files
635  IF (ir == irmax) THEN !If last row...
636  IF (iblock == readblocksmetdata) THEN !...of last block of met data
637  WRITE (grid_txt, '(I10)') grididmatrix(igrid)
638  ! WRITE(year_txtNext,'(I4)') year_int+1 !Get next year as a string format
639  ! FileCodeX = TRIM(FileCode)//TRIM(ADJUSTL(grid_txt))//'_'//TRIM(year_txt)
640  ! FileCodeXNext = TRIM(FileCode)//TRIM(ADJUSTL(grid_txt))//'_'//TRIM(year_txtNext)
641  ! CALL NextInitial(FileCodeXNext,year_int)
642  filecodexwy = trim(filecode)//trim(adjustl(grid_txt)) !File code without year (HCW 24 May 2016)
643  IF (diagnose == 1) WRITE (*, *) 'Calling NextInitial...'
644  CALL nextinitial(filecodexwy, year_int)
645 
646  ENDIF
647  ENDIF
648 
649  gridcounter = gridcounter + 1 !Increase GridCounter by 1 for next grid
650  ENDDO !end loop over grids
651 
652  ! update simulation time since start
654 
655  !! TODO: water movements between the grids needs to be taken into account here
656 
657  ENDDO !end loop over rows of met data
658 
659  ! Write output files in blocks --------------------------------
660 ! #ifdef nc
661 ! IF (ncMode == 0) THEN
662 ! #endif
663  DO igrid = 1, numberofgrids
664  IF (diagnose == 1) WRITE (*, *) 'Calling SUEWS_Output...'
665  CALL suews_output(irmax, iblock, igrid, year_int)
666  ENDDO
667 ! #ifdef nc
668 ! ENDIF
669 
670 ! IF (ncMode == 1) THEN
671 ! ! write resulst in netCDF
672 ! IF (Diagnose == 1) WRITE (*, *) 'Calling SUEWS_Output_nc...'
673 ! CALL SUEWS_Output(irMax)
674 ! ! write input information in netCDF as well for future development
675 ! ! IF ( iblock==1 ) THEN
676 ! ! CALL SiteSelect_txt2nc
677 ! ! ENDIF
678 ! ENDIF
679 ! #endif
680 
681  ENDDO !end loop over blocks of met data
682  !-----------------------------------------------------------------------
683 
684  ! ---- Decallocate arrays ----------------------------------------------
685  IF (diagnose == 1) WRITE (*, *) 'Deallocating arrays in SUEWS_Program.f95...'
686  DEALLOCATE (surfacechar)
687  DEALLOCATE (metforcingdata)
688  DEALLOCATE (metforcingdata_grid)
689  DEALLOCATE (modeloutputdata)
690  DEALLOCATE (dataoutsuews)
691  DEALLOCATE (dataoutrsl)
692  DEALLOCATE (dataoutsolweig)
693  DEALLOCATE (dataoutdailystate)
694  ! IF (SnowUse == 1) THEN
695  DEALLOCATE (dataoutsnow)
696  ! DEALLOCATE(qn1_S_store)
697  ! DEALLOCATE(qn1_S_av_store)
698  ! DEALLOCATE(qn1_S_store_grid)
699  ! DEALLOCATE(qn1_S_av_store_grid)
700  DEALLOCATE (qn1_s_av_grids)
701  DEALLOCATE (dqnsdt_grids)
702  ! ENDIF
703  ! IF (StorageHeatMethod==4 .OR. StorageHeatMethod==14) THEN
704  DEALLOCATE (dataoutestm) !ESTM output
705  DEALLOCATE (estmforcingdata)
706  DEALLOCATE (ts5mindata)
707  DEALLOCATE (ts5mindata_ir)
708  DEALLOCATE (tair24hr)
709  ! ENDIF
710  ! DEALLOCATE(TstepProfiles)
711  ! DEALLOCATE(AHProf_tstep)
712  ! DEALLOCATE(WUProfM_tstep)
713  ! DEALLOCATE(WUProfA_tstep)
714  ! DEALLOCATE(HumActivity_tstep)
715  ! DEALLOCATE(TraffProf_tstep)
716  ! DEALLOCATE(PopProf_tstep)
717  ! DEALLOCATE(qn1_store)
718  ! DEALLOCATE(qn1_av_store)
719  ! DEALLOCATE(qn1_store_grid)
720  ! DEALLOCATE(qn1_av_store_grid)
721  DEALLOCATE (qhforcbl)
722  DEALLOCATE (qeforcbl)
723  DEALLOCATE (tair_av_grids)
724  DEALLOCATE (qn1_av_grids)
725  DEALLOCATE (dqndt_grids)
726  IF (cbluse >= 1) THEN
727  DEALLOCATE (dataoutbl)
728  END IF
729  ! ----------------------------------------------------------------------
730 
731  ENDDO !end loop over years
732 
733  ! ---- Decallocate array --------------------------------------------------
734  ! Daily state needs to be outside year loop to transfer states between years
735  IF (ALLOCATED(modeldailystate)) DEALLOCATE (modeldailystate)
736  ! Also needs to happen at the end of the run
737  IF (ALLOCATED(usecolumnsdataout)) DEALLOCATE (usecolumnsdataout)
738  ! -------------------------------------------------------------------------
739 
740  ! get cpu time consumed
741  CALL cpu_time(timefinish)
742  WRITE (*, *) "Time = ", timefinish - timestart, " seconds."
743 
744  !Write to problems.txt that run has completed
745  IF (errorchoice == 0) THEN !if file has not been opened previously
746  OPEN (500, file='problems.txt')
747  errorchoice = 1
748  ELSE
749  OPEN (500, file='problems.txt', position="append")
750  ENDIF
751  !Writing of the problem file
752  WRITE (500, *) '--------------'
753  WRITE (500, *) 'Run completed.'
754  WRITE (500, *) '0' ! Write out error code 0 if run completed
755  CLOSE (500)
756 
757  ! Also print to screen
758  WRITE (*, *) "----- SUEWS run completed -----"
759 
760  stop 'finished'
761 
762  ! 313 CALL errorHint(11,TRIM(FileOrigMet),notUsed,notUsed,ios_out)
763  ! 314 CALL errorHint(11,TRIM(FileMet),notUsed,notUsed,ios_out)
764  ! 315 CALL errorHint(11,TRIM(fileESTMTs),notUsed,notUsed,NotUsedI)
765 
766 END PROGRAM suews_program
real(kind(1d0)) nper_real
integer readlinesorigmetdata
integer outputformats
real(kind(1d0)), dimension(:, :), allocatable metforcingdata_grid
integer iy_prev_t
character(len=90) progname
integer skippedlinesorig
integer numberofgrids
integer, parameter maxncols_cmds
character(len=15) datetime
real(kind(1d0)), dimension(:, :, :), allocatable dataoutestm
integer, parameter ncolumnsmetforcingdata
real(kind(1d0)), dimension(:), allocatable qn1_s_av_grids
integer, parameter ncolumnsdataoutdailystate
subroutine cbl_readinputdata(FileInputPath, qh_choice)
real(kind(1d0)), dimension(:), allocatable qn1_av_grids
real(kind(1d0)), dimension(:, :, :), allocatable dataoutbl
integer diagnose
real(kind(1d0)) nan
subroutine disaggregatemet(iBlock, igrid)
integer readlinesorigestmdatamax
subroutine suews_output(irMax, iv, Gridiv, iyr)
character(len=150) fileorigmet
real(kind(1d0)), dimension(:, :), allocatable ts5mindata
integer, dimension(:), allocatable usecolumnsdataout
subroutine leapyearcalc(year_int, nroDays)
integer nlinesorigmetdata
integer, parameter ncolumnsdataoutbl
integer multipleestmfiles
character(len=150) fileestmts
integer id_prev_t
subroutine suews_initializemetdata(lunit)
integer, dimension(:), allocatable grididmatrix
real(kind(1d0)), dimension(:, :), allocatable surfacechar
integer multiplemetfiles
subroutine initialstate(GridName, year_int, Gridiv, NumberOfGrids)
integer, parameter ncolumnsdataoutsol
real(kind(1d0)), dimension(:), allocatable dailystatefirstopen
integer resolutionfilesinestm
subroutine disaggregateestm(iBlock)
integer readlinesorigmetdatamax
integer id
character(len=150) filedscdestm
real(kind(1d0)), dimension(:), allocatable tair24hr
integer, parameter ndays
subroutine initializesurfacecharacteristics(Gridiv, rr)
integer, dimension(:), allocatable grididmatrix0
character(len=150) filedscdmet
subroutine nextinitial(GridName, year_int)
real(kind(1d0)) nperestm_real
character(len=150) fileinputpath
real(kind(1d0)), dimension(:, :, :), allocatable dataoutdailystate
integer, parameter ncolumnsdataoutrsl
integer, parameter ncolumnsdataoutsnow
real(kind(1d0)), dimension(:, :), allocatable siteselect
integer, parameter maxncols_cmod
integer nlinesorigestmdata
real(kind(1d0)), dimension(:, :), allocatable modeldailystate
subroutine estm_initials
integer iy
integer, parameter maxlinesmet
subroutine overallruncontrol
integer skippedlines
integer, parameter ncolsestmdata
character(len=150) filemet
integer dt_since_start
integer nofdaysthisyear
integer nlinessiteselect
character(len=20) filecode
real(kind(1d0)), dimension(:), allocatable qeforcbl
integer readblocksorigmetdata
character(len=10) gridid_text
integer nlinesmetdata
integer firstyear
real(kind(1d0)), dimension(:, :, :), allocatable dataoutsuews
integer function count_lines(filename)
integer cbluse
real(kind(1d0)), dimension(:, :, :), allocatable dataoutrsl
subroutine suews_getestmdata(lunit)
integer, parameter ncolumnsdataoutsuews
real(kind(1d0)), dimension(:), allocatable dqnsdt_grids
real(kind(1d0)), dimension(:), allocatable qhforcbl
integer readlinesorigestmdata
subroutine suews_calculations(Gridiv, ir, iMB, irMax)
integer readblocksmetdata
character(len=150) fileorigestm
integer resolutionfilesin
real(kind(1d0)), dimension(:, :, :), allocatable metforcingdata
real(kind(1d0)), dimension(:, :, :), allocatable dataoutsnow
integer, parameter maxncols_c
integer, parameter maxnumberofgrids
real(kind(1d0)), dimension(:, :, :), allocatable estmforcingdata
real(kind(1d0)), dimension(:), allocatable tair_av_grids
integer readlinesmetdata
integer lastyear
real(kind(1d0)), dimension(:), allocatable dqndt_grids
program suews_program
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
integer storageheatmethod
real(kind(1d0)), dimension(:, :, :), allocatable modeloutputdata
real(kind(1d0)), dimension(:), allocatable ts5mindata_ir
real(kind(1d0)), dimension(:, :, :), allocatable dataoutsolweig
integer gridcounter
integer skippedlinesorigestm