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 (dataOutSOL(ReadLinesMetdata, ncolumnsdataOutSOL, NumberOfGrids)) !SOLWEIG POI output
271  IF (cbluse >= 1) ALLOCATE (dataoutbl(readlinesmetdata, ncolumnsdataoutbl, numberofgrids)) !CBL output
272  ! IF (SnowUse == 1) THEN
273  IF (.NOT. ALLOCATED(dataoutsnow)) ALLOCATE (dataoutsnow(readlinesmetdata, ncolumnsdataoutsnow, numberofgrids)) !Snow output
274  ! ALLOCATE(qn1_S_store(NSH,NumberOfGrids))
275  ! ALLOCATE(qn1_S_store_grid(NSH))
276  ! ALLOCATE(qn1_S_av_store(2*NSH+1,NumberOfGrids))
277  ! ALLOCATE(qn1_S_av_store_grid(2*NSH+1))
278  IF (.NOT. ALLOCATED(qn1_s_av_grids)) ALLOCATE (qn1_s_av_grids(numberofgrids))
279  IF (.NOT. ALLOCATED(dqnsdt_grids)) ALLOCATE (dqnsdt_grids(numberofgrids))
280  ! qn1_S_store(:,:) = NAN
281  ! qn1_S_av_store(:,:) = NaN
282  ! qn1_S_store_grid(:) = NAN
283  ! qn1_S_av_store_grid(:) = NaN
284  qn1_s_av_grids = 0 ! Initialise to 0
285  dqnsdt_grids = 0 ! Initialise to 0
286  ! ENDIF
287  ! IF (StorageHeatMethod==4 .OR. StorageHeatMethod==14) THEN
288  IF (.NOT. ALLOCATED(dataoutestm)) ALLOCATE (dataoutestm(readlinesmetdata, 32, numberofgrids)) !ESTM output
289  ! ENDIF
290  ! ALLOCATE(TstepProfiles(NumberOfGrids,14,24*NSH)) !Hourly profiles interpolated to model timestep
291  ! ALLOCATE(AHProf_tstep(24*NSH,2)) !Anthropogenic heat profiles at model timestep
292  ! ALLOCATE(WUProfM_tstep(24*NSH,2)) !Manual water use profiles at model timestep
293  ! ALLOCATE(WUProfA_tstep(24*NSH,2)) !Automatic water use profiles at model timestep
294  ! ALLOCATE(HumActivity_tstep(24*NSH,2)) !Human activity profiles at model timestep
295  ! ALLOCATE(TraffProf_tstep(24*NSH,2)) !Traffic profiles at model timestep
296  ! ALLOCATE(PopProf_tstep(24*NSH,2)) !Population profiles at model timestep
297  ! ALLOCATE(qn1_store(NSH,NumberOfGrids))
298  ! ALLOCATE(qn1_av_store(2*NSH+1,NumberOfGrids))
299  ! ALLOCATE(qn1_store_grid(NSH))
300  ! ALLOCATE(qn1_av_store_grid(2*NSH+1))
301  ALLOCATE (qhforcbl(numberofgrids))
302  ALLOCATE (qeforcbl(numberofgrids))
303 
304  ALLOCATE (tair_av_grids(numberofgrids))
305  ALLOCATE (qn1_av_grids(numberofgrids))
306  ALLOCATE (dqndt_grids(numberofgrids))
307  !! QUESTION: Add snow clearing (?)
308 
309  ! qn1_store(:,:) = NAN ! Initialise to -999
310  ! qn1_av_store(:,:) = NAN ! Initialise to -999
311  ! qn1_store_grid(:)=NAN
312  ! qn1_av_store_grid(:)=NAN
313  tair_av_grids = 273.15 ! Initialise to 273.15 K
314  qn1_av_grids = 0 ! Initialise to 0
315  dqndt_grids = 0 ! Initialise to 0
316 
317  qhforcbl(:) = nan
318  qeforcbl(:) = nan
319 
320  ! QUESTION: Initialise other arrays here?
321 
322  IF (storageheatmethod == 4 .OR. storageheatmethod == 14) THEN
323  ! Prepare to disaggregate ESTM data to model time-step (if required) ------
324  ! Find number of model time-steps per resolution of original met forcing file
325  nperestm_real = resolutionfilesinestm/REAL(tstep, kind(1d0))
326  nperestm = int(nperestm_real)
327  IF (nperestm /= nperestm_real) THEN
328  CALL errorhint(2, 'Problem in SUEWS_Program: check resolution of ESTM forcing data (ResolutionFilesInESTM)'// &
329  'and model time-step (Tstep).', &
330  REAL(Tstep, KIND(1d0)), NotUsed, ResolutionFilesInESTM)
331  ELSEIF (nperestm > 1) THEN
332  WRITE (*, *) 'Resolution of ESTM forcing data: ', trim(adjustl(resinestm_txt)), ' min;', &
333  ' model time-step: ', trim(adjustl(tstep_txt)), ' min', ' -> SUEWS will perform disaggregation.'
334  IF (diagnose == 1) WRITE (*, *) 'Getting information for ESTM disaggregation'
335  ! Get names of original met forcing file(s) to disaggregate (using first grid)
336  WRITE (grid_txt, '(I10)') grididmatrix(1) !Get grid as a text string
337 
338  ! Get met file name for this grid
339  fileestmts = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)//'_ESTM_Ts_data_' &
340  //trim(adjustl(resinestm_txt))//'.txt'
341  ! But if each grid has the same ESTM file, file name does not include grid number
342  IF (multipleestmfiles /= 1) THEN
343  fileestmts = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_ESTM_Ts_data_' &
344  //trim(adjustl(resinestm_txt))//'.txt'
345  ENDIF
346 
347  ! Find number of lines in orig ESTM file
348  nlinesorigestmdata = 0 !Initialise nlinesMetdata (total number of lines in met forcing file)
350 
351  ! Check ESTM data and met data will have the same length (so that ESTM file can be read in same blocks as met data)
353  CALL errorhint(66, &
354  'Downscaled ESTM and met input files will have different lengths', REAL(nlinesMetdata, KIND(1d0)), &
355  NotUsed, nlinesESTMdata*NperESTM)
356  ENDIF
357 
358  !write(*,*) 'nlinesOrigESTMdata', nlinesOrigESTMdata
359  ! Set number of lines to read from original ESTM file using met data blocks
361  !WRITE(*,*) 'ReadlinesOrigESTMdata', ReadlinesOrigESTMdata
362  WRITE (*, *) 'Original ESTM data will be read in chunks of ', readlinesorigestmdata, 'lines.'
363 
364  nlinesestmdata = nlinesorigestmdata*nperestm
365 
366  ELSEIF (nperestm == 1) THEN
367  WRITE (*, *) 'ResolutionFilesInESTM = Tstep: no disaggregation needed for met data.'
368 
369  !-----------------------------------------------------------------------
370  ! Find number of lines in ESTM forcing file for current year (nlinesESTMdata)
371  WRITE (grid_txt, '(I10)') grididmatrix(1) !Get grid as a text string (use first grid as example)
372  ! Get file name for this year for this grid
373  filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
374  fileestmts = trim(fileinputpath)//trim(filecodex)//'_ESTM_Ts_data_'//trim(adjustl(tstep_txt))//'.txt'
375  !If each grid has the same met file, met file name does not include grid number
376  IF (multipleestmfiles /= 1) THEN
377  filecodexwg = trim(filecode)//'_'//trim(year_txt) !File code without grid
378  fileestmts = trim(fileinputpath)//trim(filecodexwg)//'_ESTM_Ts_data_'//trim(adjustl(tstep_txt))//'.txt'
379  ENDIF
380 
381  ! Find number of lines in ESTM file
382  nlinesestmdata = 0 !Initialise nlinesMetdata (total number of lines in met forcing file)
383  nlinesestmdata = count_lines(trim(fileestmts))
384  !-----------------------------------------------------------------------
385 
386  ! Check ESTM data and met data are same length (so that ESTM file can be read in same blocks as met data)
387  IF (nlinesestmdata /= nlinesmetdata) THEN
388  CALL errorhint(66, 'ESTM input file different length to met forcing file', REAL(nlinesMetdata, KIND(1d0)), &
389  NotUsed, nlinesESTMdata)
390  ENDIF
391 
392  ENDIF
393 
394  ! Allocate arrays to receive ESTM forcing data
395  ! IF (.NOT. ALLOCATED(ESTMForcingData)) ALLOCATE(ESTMForcingData(1:ReadLinesMetdata,ncolsESTMdata,NumberOfGrids))
396  ! IF (.NOT. ALLOCATED(Ts5mindata)) ALLOCATE(Ts5mindata(1:ReadLinesMetdata,ncolsESTMdata))
397  ! IF (.NOT. ALLOCATED(Ts5mindata_ir)) ALLOCATE(Ts5mindata_ir(ncolsESTMdata))
398  !
399  ! IF (.NOT. ALLOCATED(Tair24HR)) ALLOCATE(Tair24HR(24*nsh))
400  ! if ( /= 0) print *, ": Deallocation request denied"
401  ! ALLOCATE(Tair24HR(24*nsh))
402 
403  ENDIF
404  IF (.NOT. ALLOCATED(estmforcingdata)) ALLOCATE (estmforcingdata(1:readlinesmetdata, ncolsestmdata, numberofgrids))
405  IF (.NOT. ALLOCATED(ts5mindata)) ALLOCATE (ts5mindata(1:readlinesmetdata, ncolsestmdata))
406  IF (.NOT. ALLOCATED(ts5mindata_ir)) ALLOCATE (ts5mindata_ir(ncolsestmdata))
407  IF (.NOT. ALLOCATED(tair24hr)) ALLOCATE (tair24hr(24*nsh))
408  ! ------------------------------------------------------------------------
409 
410  !-----------------------------------------------------------------------
411  skippedlines = 0 !Initialise lines to be skipped in met forcing file
412  skippedlinesorig = 0 !Initialise lines to be skipped in original met forcing file
413  skippedlinesorigestm = 0 !Initialise lines to be skipped in original met forcing file
414 
415  DO iblock = 1, readblocksmetdata !Loop through blocks of met data
416  ! WRITE(*,*) iblock,'/',ReadBlocksMetData
417 
418  ! Model calculations are made in two stages:
419  ! (1) initialise the run for each block of met data (iblock from 1 to ReadBlocksMetData)
420  ! (2) perform the actual model calculations (SUEWS_Calculations)
421 
422  gridcounter = 1 !Initialise counter for grids in each year
423  DO igrid = 1, numberofgrids !Loop through grids
424 
425  gridid = grididmatrix(igrid) !store grid here for referencing error codes
426  WRITE (grid_txt, '(I10)') grididmatrix(igrid) !Get grid ID as a text string
427 
428  ! (1) First stage: initialise run if this is the first iteration this year
429  ! (1a) Initialise surface characteristics
430  IF (iblock == 1) THEN
431  IF (diagnose == 1) WRITE (*, *) 'First block of data - doing initialisation'
432  ! (a) Transfer characteristics from SiteSelect to correct row of SurfaceChar
433  DO rr = 1, nlinessiteselect
434  !Find correct grid and year
435  IF (diagnose == 1) WRITE (*, *) 'grid found:', siteselect(rr, c_grid), 'grid needed:', grididmatrix(igrid)
436  IF (diagnose == 1) WRITE (*, *) 'year found:', siteselect(rr, c_year), 'year needed:', year_int
437  IF (siteselect(rr, c_grid) == grididmatrix(igrid) .AND. siteselect(rr, c_year) == year_int) THEN
438  IF (diagnose == 1) WRITE (*, *) 'Match found (grid and year) for rr = ', rr, 'of', nlinessiteselect
440  EXIT
441  ELSEIF (rr == nlinessiteselect) THEN
442  WRITE (*, *) 'Program stopped! Year', year_int, 'and/or grid', igrid, 'not found in SiteSelect.txt.'
443  CALL errorhint(59, 'Cannot find year and/or grid in SiteSelect.txt', REAL(igrid, KIND(1d0)), NotUsed, year_int)
444  ENDIF
445  ENDDO
446  ENDIF !end first block of met data
447 
448  ! (1b) Initialise met data
449  IF (nper > 1) THEN
450  ! Disaggregate met data ---------------------------------------------------
451 
452  ! Set maximum value for ReadLinesOrigMetData to handle end of file (i.e. small final block)
453  IF (iblock == readblocksmetdata) THEN !For last block of data in file
455  ELSE
457  ENDIF
458  !write(*,*) ReadLinesOrigMetDataMAX, ReadLinesOrigMetData
459  ! Get names of original met forcing file(s) to disaggregate
460  ! Get met file name for this grid: SSss_YYYY_data_RR.txt
461  IF (multiplemetfiles == 1) THEN !If each grid has its own met file
462  fileorigmet = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)//'_data_' &
463  //trim(adjustl(resin_txt))//'.txt'
464  ! Also set file name for downscaled file
465  filedscdmet = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)//'_data_' &
466  //trim(adjustl(tstep_txt))//'.txt'
467  ! Disaggregate met data
468  CALL disaggregatemet(iblock, igrid)
469  ELSE
470  ! If each grid has the same met file, met file name does not include grid number, and only need to disaggregate once
471  fileorigmet = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_data_' &
472  //trim(adjustl(resin_txt))//'.txt'
473  filedscdmet = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_data_' &
474  //trim(adjustl(tstep_txt))//'.txt'
475  IF (igrid == 1) THEN !Disaggregate for the first grid only
476  CALL disaggregatemet(iblock, igrid)
477  ELSE !Then for subsequent grids simply copy data
479  ENDIF
480  ENDIF
481 
482  ELSEIF (nper == 1) THEN
483  ! Get met forcing file name for this year for the first grid
484  ! Can be something else than 1
485  filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
486  filecodexwg = trim(filecode)//'_'//trim(year_txt) !File code without grid
487  ! IF(iblock==1) WRITE(*,*) 'Current FileCode: ', FileCodeX
488 
489  ! For every block of met data ------------------------------------
490  ! Initialise met forcing data into 3-dimensional matrix
491  !write(*,*) 'Initialising met data for block',iblock
492  IF (multiplemetfiles == 1) THEN !If each grid has its own met file
493  filemet = trim(fileinputpath)//trim(filecodex)//'_data_'//trim(adjustl(tstep_txt))//'.txt'
495  ELSE !If one met file used for all grids
496  !FileMet=TRIM(FileInputPath)//TRIM(FileCodeX)//'_data_'//TRIM(ADJUSTL(tstep_txt))//'.txt'
497  ! If one met file used for all grids, look for met file with no grid code (FileCodeXWG)
498  filemet = trim(fileinputpath)//trim(filecodexwg)//'_data_'//trim(adjustl(tstep_txt))//'.txt'
499  IF (igrid == 1) THEN !Read for the first grid only
501  ELSE !Then for subsequent grids simply copy data
503  ENDIF
504  ENDIF
505  ENDIF !end of nper statement
506 
507  ! Only for the first block of met data, read initial conditions (moved from above, HCW 12 Jan 2017)
508  IF (iblock == 1) THEN
509  filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
510  !write(*,*) ' Now calling InitialState'
511  CALL initialstate(filecodex, year_int, gridcounter, numberofgrids)
512  ENDIF
513 
514  ! Initialise ESTM if required, TS 05 Jun 2016; moved inside grid loop HCW 27 Jun 2016
515  IF (storageheatmethod == 4 .OR. storageheatmethod == 14) THEN
516  IF (nperestm > 1) THEN
517  ! Disaggregate ESTM data --------------------------------------------------
518  ! Set maximum value for ReadLinesOrigESTMData to handle end of file (i.e. small final block)
519  IF (iblock == readblocksmetdata) THEN !For last block of data in file
521  ELSE
523  ENDIF
524  !write(*,*) ReadLinesOrigESTMDataMAX, ReadLinesOrigESTMData
525  ! Get names of original ESTM forcing file(s) to disaggregate
526  ! Get ESTM file name for this grid: SSss_YYYY_ESTM_Ts_data_RR.txt
527  IF (multipleestmfiles == 1) THEN !If each grid has its own ESTM file
528  fileorigestm = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt) &
529  //'_ESTM_Ts_data_'//trim(adjustl(resinestm_txt))//'.txt'
530  ! Also set file name for downscaled file
531  filedscdestm = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt) &
532  //'_ESTM_Ts_data_'//trim(adjustl(tstep_txt))//'.txt'
533  ! Disaggregate ESTM data
534  CALL disaggregateestm(iblock)
535  ELSE
536  ! If each grid has the same ESTM file, ESTM file name does not include grid number, and only need to disaggregate once
537  fileorigestm = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_ESTM_Ts_data_' &
538  //trim(adjustl(resinestm_txt))//'.txt'
539  filedscdestm = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_ESTM_Ts_data_' &
540  //trim(adjustl(tstep_txt))//'.txt'
541  IF (igrid == 1) THEN !Disaggregate for the first grid only
542  CALL disaggregateestm(iblock)
543  ELSE !Then for subsequent grids simply copy data
545  1:ncolsestmdata, 1)
546  ENDIF
547  ENDIF
548 
549  ELSEIF (nperestm == 1) THEN
550  ! Get ESTM forcing file name for this year for the first grid
551  filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
552  filecodexwg = trim(filecode)//'_'//trim(year_txt) !File code without grid
553  ! For every block of ESTM data ------------------------------------
554  ! Initialise ESTM forcing data into 3-dimensional matrix
555  !write(*,*) 'Initialising ESTM data for block',iblock
556  IF (multipleestmfiles == 1) THEN !If each grid has its own met file
557  fileestmts = trim(fileinputpath)//trim(filecodex)//'_ESTM_Ts_data_'//trim(adjustl(tstep_txt))//'.txt'
558  !write(*,*) 'Calling GetESTMData...', FileCodeX, iblock, igrid
559  CALL suews_getestmdata(101)
560  ELSE !If one ESTM file used for all grids
561  fileestmts = trim(fileinputpath)//trim(filecodexwg)//'_ESTM_Ts_data_'//trim(adjustl(tstep_txt))//'.txt'
562  !write(*,*) 'Calling GetESTMData...', FileCodeX, iblock, igrid
563  IF (igrid == 1) THEN !Read for the first grid only
564  CALL suews_getestmdata(101)
565  ELSE !Then for subsequent grids simply copy data
567  1:ncolsestmdata, 1)
568  ENDIF
569  ENDIF
570  ENDIF !end of nperESTM statement
571  ENDIF
572 
573  gridcounter = gridcounter + 1 !Increase GridCounter by 1 for next grid
574 
575  ENDDO !end loop over grids
576  skippedlines = skippedlines + readlinesmetdata !Increase skippedLines ready for next block
577  skippedlinesorig = skippedlinesorig + readlinesorigmetdata !Increase skippedLinesOrig ready for next block
578  skippedlinesorigestm = skippedlinesorigestm + readlinesorigestmdata !Increase skippedLinesOrig ready for next block
579  !write(*,*) iblock
580  !write(*,*) ReadLinesMetdata, readlinesorigmetdata
581  !write(*,*) skippedLines, skippedLinesOrig, skippedLinesOrig*Nper
582 
583  ! Initialise the modules on the first day
584  ! Initialise CBL and SOLWEIG parts if required
585  IF (iblock == 1) THEN
586  IF ((cbluse == 1) .OR. (cbluse == 2)) CALL cbl_readinputdata(fileinputpath, qh_choice)
587  ENDIF
588 
589  ! NB: SOLWEIG is disabled in v2018a
590  ! QUESTION: not sure if solweig should be within iblock if statement too?
591  ! IF(SOLWEIGuse==1) CALL SOLWEIG_initial
592 
593  !write(*,*) 'Initialisation done'
594  ! First stage: initialisation done ----------------------------------
595 
596  ! (2) Second stage: do calculations at 5-min time-steps -------------
597  ! First set maximum value of ir
598  IF (iblock == readblocksmetdata) THEN !For last block of data in file
599  irmax = nlinesmetdata - (iblock - 1)*readlinesmetdata
600  ELSE
601  irmax = readlinesmetdata
602  ENDIF
603 
604  DO ir = 1, irmax !Loop through rows of current block of met data
605  gridcounter = 1 !Initialise counter for grids in each year
606 
607  DO igrid = 1, numberofgrids !Loop through grids
608  IF (diagnose == 1) WRITE (*, *) 'Row (ir):', ir, '/', irmax, 'of block (iblock):', iblock, '/', readblocksmetdata, &
609  'Grid:', grididmatrix(igrid)
610 
611  ! Call model calculation code
612  ! IF(ir==1) WRITE(*,*) 'Now running block ',iblock,'/',ReadBlocksMetData,' of year ',year_int,'...'
613  WRITE (grid_txt, '(I10)') grididmatrix(igrid) !Get grid ID as a text string
614  filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
615  IF (ir == 1 .AND. igrid == 1) THEN
616  WRITE (*, *) trim(adjustl(filecodex)), &
617  ': Now running block ', iblock, '/', readblocksmetdata, ' of ', trim(year_txt), '...'
618  ENDIF
619  IF (diagnose == 1) WRITE (*, *) 'Calling SUEWS_Calculations...'
620  ! print*, 'before cal:',sum(Tair24HR)
621  CALL suews_calculations(gridcounter, ir, iblock, irmax)
622  ! print*, 'after cal:',sum(Tair24HR)
623  IF (diagnose == 1) WRITE (*, *) 'SUEWS_Calculations finished...'
624 
625  ! Record iy and id for current time step to handle last row in yearly files (YYYY 1 0 0)
626  ! IF(GridCounter == NumberOfGrids) THEN !Adjust only when the final grid has been run for this time step
627  IF (igrid == numberofgrids) THEN !Adjust only when the final grid has been run for this time step
628  iy_prev_t = iy
629  id_prev_t = id
630  ENDIF
631 
632  ! Write state information to new InitialConditions files
633  IF (ir == irmax) THEN !If last row...
634  IF (iblock == readblocksmetdata) THEN !...of last block of met data
635  WRITE (grid_txt, '(I10)') grididmatrix(igrid)
636  ! WRITE(year_txtNext,'(I4)') year_int+1 !Get next year as a string format
637  ! FileCodeX = TRIM(FileCode)//TRIM(ADJUSTL(grid_txt))//'_'//TRIM(year_txt)
638  ! FileCodeXNext = TRIM(FileCode)//TRIM(ADJUSTL(grid_txt))//'_'//TRIM(year_txtNext)
639  ! CALL NextInitial(FileCodeXNext,year_int)
640  filecodexwy = trim(filecode)//trim(adjustl(grid_txt)) !File code without year (HCW 24 May 2016)
641  IF (diagnose == 1) WRITE (*, *) 'Calling NextInitial...'
642  CALL nextinitial(filecodexwy, year_int)
643 
644  ENDIF
645  ENDIF
646 
647  gridcounter = gridcounter + 1 !Increase GridCounter by 1 for next grid
648  ENDDO !end loop over grids
649 
650  ! update simulation time since start
652 
653  !! TODO: water movements between the grids needs to be taken into account here
654 
655  ENDDO !end loop over rows of met data
656 
657  ! Write output files in blocks --------------------------------
658 ! #ifdef nc
659 ! IF (ncMode == 0) THEN
660 ! #endif
661  DO igrid = 1, numberofgrids
662  IF (diagnose == 1) WRITE (*, *) 'Calling SUEWS_Output...'
663  CALL suews_output(irmax, iblock, igrid, year_int)
664  ENDDO
665 ! #ifdef nc
666 ! ENDIF
667 
668 ! IF (ncMode == 1) THEN
669 ! ! write resulst in netCDF
670 ! IF (Diagnose == 1) WRITE (*, *) 'Calling SUEWS_Output_nc...'
671 ! CALL SUEWS_Output(irMax)
672 ! ! write input information in netCDF as well for future development
673 ! ! IF ( iblock==1 ) THEN
674 ! ! CALL SiteSelect_txt2nc
675 ! ! ENDIF
676 ! ENDIF
677 ! #endif
678 
679  ENDDO !end loop over blocks of met data
680  !-----------------------------------------------------------------------
681 
682  ! ---- Decallocate arrays ----------------------------------------------
683  IF (diagnose == 1) WRITE (*, *) 'Deallocating arrays in SUEWS_Program.f95...'
684  DEALLOCATE (surfacechar)
685  DEALLOCATE (metforcingdata)
686  DEALLOCATE (metforcingdata_grid)
687  DEALLOCATE (modeloutputdata)
688  DEALLOCATE (dataoutsuews)
689  DEALLOCATE (dataoutrsl)
690  DEALLOCATE (dataoutdailystate)
691  ! IF (SnowUse == 1) THEN
692  DEALLOCATE (dataoutsnow)
693  ! DEALLOCATE(qn1_S_store)
694  ! DEALLOCATE(qn1_S_av_store)
695  ! DEALLOCATE(qn1_S_store_grid)
696  ! DEALLOCATE(qn1_S_av_store_grid)
697  DEALLOCATE (qn1_s_av_grids)
698  DEALLOCATE (dqnsdt_grids)
699  ! ENDIF
700  ! IF (StorageHeatMethod==4 .OR. StorageHeatMethod==14) THEN
701  DEALLOCATE (dataoutestm) !ESTM output
702  DEALLOCATE (estmforcingdata)
703  DEALLOCATE (ts5mindata)
704  DEALLOCATE (ts5mindata_ir)
705  DEALLOCATE (tair24hr)
706  ! ENDIF
707  ! DEALLOCATE(TstepProfiles)
708  ! DEALLOCATE(AHProf_tstep)
709  ! DEALLOCATE(WUProfM_tstep)
710  ! DEALLOCATE(WUProfA_tstep)
711  ! DEALLOCATE(HumActivity_tstep)
712  ! DEALLOCATE(TraffProf_tstep)
713  ! DEALLOCATE(PopProf_tstep)
714  ! DEALLOCATE(qn1_store)
715  ! DEALLOCATE(qn1_av_store)
716  ! DEALLOCATE(qn1_store_grid)
717  ! DEALLOCATE(qn1_av_store_grid)
718  DEALLOCATE (qhforcbl)
719  DEALLOCATE (qeforcbl)
720  DEALLOCATE (tair_av_grids)
721  DEALLOCATE (qn1_av_grids)
722  DEALLOCATE (dqndt_grids)
723  IF (cbluse >= 1) THEN
724  DEALLOCATE (dataoutbl)
725  END IF
726  ! ----------------------------------------------------------------------
727 
728  ENDDO !end loop over years
729 
730  ! ---- Decallocate array --------------------------------------------------
731  ! Daily state needs to be outside year loop to transfer states between years
732  IF (ALLOCATED(modeldailystate)) DEALLOCATE (modeldailystate)
733  ! Also needs to happen at the end of the run
734  IF (ALLOCATED(usecolumnsdataout)) DEALLOCATE (usecolumnsdataout)
735  ! -------------------------------------------------------------------------
736 
737  ! get cpu time consumed
738  CALL cpu_time(timefinish)
739  WRITE (*, *) "Time = ", timefinish - timestart, " seconds."
740 
741  !Write to problems.txt that run has completed
742  IF (errorchoice == 0) THEN !if file has not been opened previously
743  OPEN (500, file='problems.txt')
744  errorchoice = 1
745  ELSE
746  OPEN (500, file='problems.txt', position="append")
747  ENDIF
748  !Writing of the problem file
749  WRITE (500, *) '--------------'
750  WRITE (500, *) 'Run completed.'
751  WRITE (500, *) '0' ! Write out error code 0 if run completed
752  CLOSE (500)
753 
754  ! Also print to screen
755  WRITE (*, *) "----- SUEWS run completed -----"
756 
757  stop 'finished'
758 
759  ! 313 CALL errorHint(11,TRIM(FileOrigMet),notUsed,notUsed,ios_out)
760  ! 314 CALL errorHint(11,TRIM(FileMet),notUsed,notUsed,ios_out)
761  ! 315 CALL errorHint(11,TRIM(fileESTMTs),notUsed,notUsed,NotUsedI)
762 
763 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)
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
integer gridcounter
integer skippedlinesorigestm