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
36 USE data_in
38 USE initial
39 USE metdisagg
40 USE sues_data
41 USE time
42 USE wherewhen
43 USE ctrl_output
44 USE estm_module, ONLY: &
50
51 IMPLICIT NONE
52
53 CHARACTER(len=4) :: year_txt !Year as a text string
54
55 CHARACTER(len=20) :: filecodex, & !Current file code
56 filecodexwy, & !File code without year
57 filecodexwg !File code without grid
58 CHARACTER(len=20) :: grid_txt, & !Grid number as a text string (from FirstGrid to LastGrid)
59 tstep_txt, & !Model timestep (in minutes) as a text string (minutes)
60 resin_txt, resinestm_txt !Resolution of Original met/ESTM forcing file as a text string (minutes)
61
62 INTEGER :: readlinesmetdata_read !Max number of lines that can be read in one go for each grid
63 INTEGER :: nlineslimit, & !Max number of lines that can be read in one go for each grid
64 numberofyears !Number of years to be run
65
66 INTEGER :: year_int, & ! Year as an integer (from SiteSelect rather than met forcing file)
67 igrid, & !Grid number (from 1 to NumberOfGrids)
68 iblock, & !Block number (from 1 to ReadBlocksMetData)
69 ir, irmax, & !Row number within each block (from 1 to irMax)
70 rr !Row of SiteSelect corresponding to current year and grid
71
72 REAL :: timestart, timefinish ! profiling use, AnOHM TS
73
74 ! Start counting cpu time
75 CALL cpu_time(timestart)
76
77 ! initialise simulation time
79
80 WRITE (*, *) '========================================================'
81 WRITE (*, *) 'Running ', progname
82 WRITE (*, *) 'Version commit: ', trim(git_commit)
83 WRITE (*, *) 'Compiler: ', trim(compiler_ver)
84
85 ! Initialise error file (0 -> problems.txt file will be newly created)
86 errorchoice = 0
87 ! Initialise error file (0 -> warnings.txt file will be newly created)
89 ! Initialise OutputFormats to 1 so that output format is written out only once per run
91
92 ! Initialise WhereWhen variables for error handling
93 gridid_text = '00000'
94 datetime = '00000'
95
96 ! Read RunControl.nml and all .txt input files from SiteSelect spreadsheet
98
99 WRITE (tstep_txt, '(I5)') tstep/60 !Get tstep (in minutes) as a text string
100 WRITE (resin_txt, '(I5)') resolutionfilesin/60 !Get ResolutionFilesIn (in minutes) as a text string
101 WRITE (resinestm_txt, '(I5)') resolutionfilesinestm/60
102
103 ! Find first and last year of the current run
104 firstyear = minval(int(siteselect(:, c_year)))
105 lastyear = maxval(int(siteselect(:, c_year)))
106
107 numberofyears = lastyear - firstyear + 1 !Find the number of years to run
108
109 !Find the the number of grids within each year in SUEWS_SiteSelect.txt
110 ! N.B. need to have the same grids for each year
111 numberofgrids = int(nlinessiteselect/numberofyears)
112
113 !! Find the first and last grid numbers (N.B. need to have the same grids for each year)
114 !FirstGrid = minval(int(SiteSelect(:,c_Grid)))
115 !LastGrid = maxval(int(SiteSelect(:,c_Grid)))
117 CALL errorhint(64, &
118 'No. of grids exceeds max. possible no. of grids.', &
119 REAL(maxnumberofgrids, kind(1d0)), notused, numberofgrids)
120 END IF
121
122 ALLOCATE (grididmatrix(numberofgrids)) !Get the nGrid numbers correctly
123 ALLOCATE (grididmatrix0(numberofgrids)) !Get the nGrid numbers correctly
124
125 DO igrid = 1, numberofgrids
126 grididmatrix(igrid) = int(siteselect(igrid, c_grid))
127 END DO
128
129! #ifdef nc
130! ! sort grid matrix to conform the geospatial layout as in QGIS, TS 14 Dec 2016
131! IF (ncMode == 1) THEN
132! GridIDmatrix0 = GridIDmatrix
133! CALL sortGrid(GridIDmatrix0, GridIDmatrix, nRow, nCol)
134! ENDIF
135! ! GridIDmatrix0 stores the grid ID in the original order
136! #endif
137
138 ! GridIDmatrix=GridIDmatrix0
139 WRITE (*, *) '--------------------------------------------'
140 WRITE (*, *) 'Years identified:', firstyear, 'to', lastyear
141 WRITE (*, *) 'No. grids identified:', numberofgrids, 'grids'
142 WRITE (*, *) 'Maximum No. grids allowed:', maxnumberofgrids, 'grids'
143
144 ! Set limit on number of lines to read
145 nlineslimit = int(floor(maxlinesmet/real(numberofgrids, kind(1d0)))) !Uncommented HCW 29 Jun 2016
146 !nlinesLimit = 24*nsh !Commented out HCW 29 Jun 2016
147
148 ! ---- Allocate arrays ----------------------------------------------------
149 ! Daily state needs to be outside year loop to transfer states between years
150 ALLOCATE (modeldailystate(numberofgrids, maxncols_cmds)) !DailyState
151 ALLOCATE (dailystatefirstopen(numberofgrids)) !Initialisation for header
152
153 ! ---- Initialise arrays --------------------------------------------------
154 modeldailystate(:, :) = -999
156
157 ! -------------------------------------------------------------------------
158 ! Initialise ESTM (reads ESTM nml, should only run once)
159 IF (storageheatmethod == 4 .OR. storageheatmethod == 14) THEN
160 IF (diagnose == 1) WRITE (*, *) 'Calling ESTM_initials...'
161 CALL estm_initials
162 END IF
163
164 ! -------------------------------------------------------------------------
165 ! Initialise ESTM (reads ESTM nml, should only run once)
166 IF (storageheatmethod == 5 .OR. netradiationmethod > 1000) THEN
167 IF (diagnose == 1) WRITE (*, *) 'Calling ESTM_ext_initialise...'
169 WRITE (*, *) 'No. vertical layers identified:', nlayer, 'layers'
170 END IF
171
172 ! -------------------------------------------------------------------------
173 ! Initialise SPARTACUS (reads SPARTACUS nml, should only run once)
174 IF (netradiationmethod > 1000) THEN
175 IF (diagnose == 1) WRITE (*, *) 'Calling ESTM_initials...'
177 END IF
178
179 !==========================================================================
180 DO year_int = firstyear, lastyear !Loop through years
181
182 WRITE (*, *) ' '
183 WRITE (year_txt, '(I4)') year_int !Get year as a text string
184
185 ! Find number of days in the current year
186 CALL leapyearcalc(year_int, nofdaysthisyear)
187
188 ! Prepare to disaggregate met data to model time-step (if required) ------
189 ! Find number of model time-steps per resolution of original met forcing file
190 npertstepin_real = resolutionfilesin/real(tstep, kind(1d0))
192
193 IF (npertstepin /= npertstepin_real) THEN
194 CALL errorhint(2, 'Problem in SUEWS_Program: check resolution of met forcing data (ResolutionFilesIn)'// &
195 'and model time-step (Tstep).', &
196 REAL(tstep, kind(1d0)), notused, resolutionfilesin)
197 ELSEIF (npertstepin > 1) THEN
198 WRITE (*, *) 'Resolution of met forcing data: ', trim(adjustl(resin_txt)), ' min;', &
199 ' model time-step: ', trim(adjustl(tstep_txt)), ' min', ' -> SUEWS will perform disaggregation.'
200 IF (diagnose == 1) WRITE (*, *) 'Getting information for met disaggregation'
201 ! Get names of original met forcing file(s) to disaggregate (using first grid)
202 WRITE (grid_txt, '(I10)') grididmatrix(1) !Get grid as a text string
203
204 ! Get met file name for this grid: SSss_YYYY_data_RR.txt
205 fileorigmet = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)//'_data_' &
206 //trim(adjustl(resin_txt))//'.txt'
207 ! But if each grid has the same met file, met file name does not include grid number
208 IF (multiplemetfiles /= 1) THEN
209 fileorigmet = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_data_' &
210 //trim(adjustl(resin_txt))//'.txt'
211 END IF
212
213 nlinesorigmetdata = 0 !Initialise nlinesMetdata (total number of lines in met forcing file)
215 ! WRITE(*,*) 'nlinesOrigMetdata', nlinesOrigMetdata
216
217 readlinesorigmetdata = nlinesorigmetdata !Initially set limit as the size of file
218 IF (nlinesorigmetdata*npertstepin > nlineslimit) THEN !But restrict if this limit exceeds memory capacity
219 readlinesorigmetdata = int(nlineslimit/npertstepin)
220 END IF
221 ! make sure the metblocks read in consists of complete diurnal cycles
224 !WRITE(*,*) 'ReadlinesOrigMetdata', ReadlinesOrigMetdata
225 WRITE (*, *) 'Original met data will be read in chunks of ', readlinesorigmetdata, 'lines.'
226
227 readblocksorigmetdata = int(ceiling(real(nlinesorigmetdata, kind(1d0))/real(readlinesorigmetdata, kind(1d0))))
228
229 ! Set ReadLinesMetdata and ReadBlocksMetData
231 readblocksmetdata = int(ceiling(real(nlinesorigmetdata*npertstepin, kind(1d0))/real(readlinesmetdata, kind(1d0))))
232 WRITE (*, *) 'Processing current year in ', readblocksmetdata, 'blocks.'
233
235
236 ELSEIF (npertstepin == 1) THEN
237 WRITE (*, *) 'ResolutionFilesIn = Tstep: no disaggregation needed for met data.'
238
239 !-----------------------------------------------------------------------
240 ! Find number of lines in met forcing file for current year (nlinesMetdata)
241 ! Need to know how many lines will be read each iteration
242 ! Use first grid as an example as the number of lines is the same for all grids
243 ! within one year
244 WRITE (grid_txt, '(I10)') grididmatrix(1) !Get grid as a text string
245
246 ! Get met file name for this year for this grid
247 filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
248 filemet = trim(fileinputpath)//trim(filecodex)//'_data_'//trim(adjustl(tstep_txt))//'.txt'
249 !If each grid has the same met file, met file name does not include grid number (HCW 27 Jun 2016)
250 IF (multiplemetfiles /= 1) THEN
251 filecodexwg = trim(filecode)//'_'//trim(year_txt) !File code without grid
252 filemet = trim(fileinputpath)//trim(filecodexwg)//'_data_'//trim(adjustl(tstep_txt))//'.txt'
253 END IF
254
255 nlinesmetdata = 0 !Initialise nlinesMetdata (total number of lines in met forcing file)
257 !-----------------------------------------------------------------------
258
259 ! To conserve memory, read met data in blocks
260 ! Find number of lines that can be read in each block (i.e. read in at once)
261 readlinesmetdata = nlinesmetdata !Initially set limit as the size of the met file (N.B.solves problem with Intel fortran)
262 IF (nlinesmetdata > nlineslimit) THEN !But restrict if this limit exceeds memory capacity
263 readlinesmetdata = nlineslimit
264 END IF
265 ! make sure the metblocks read in consists of complete diurnal cycles, TS 08 Jul 2016
266 ! ReadLinesMetdata = INT(MAX(nsd*(ReadLinesMetdata/nsd), nsd))
267
268 WRITE (*, *) 'Met data will be read in blocks of ', readlinesmetdata, 'lines.'
269
270 ! Find number of blocks of met data
271 readblocksmetdata = int(ceiling(real(nlinesmetdata, kind(1d0))/real(readlinesmetdata, kind(1d0))))
272 WRITE (*, *) 'Processing current year in ', readblocksmetdata, 'blocks.'
273
274 END IF
275
276 ! ---- Allocate arrays--------------------------------------------------
277 IF (diagnose == 1) WRITE (*, *) 'Allocating arrays in SUEWS_Program.f95...'
278 IF (.NOT. ALLOCATED(surfacechar)) ALLOCATE (surfacechar(numberofgrids, maxncols_c)) !Surface characteristics
279 IF (.NOT. ALLOCATED(metforcingdata)) ALLOCATE (metforcingdata(readlinesmetdata, ncolumnsmetforcingdata, numberofgrids)) !Met forcing data
280 IF (.NOT. ALLOCATED(metforcingdata_grid)) ALLOCATE (metforcingdata_grid(readlinesmetdata, ncolumnsmetforcingdata)) !Met forcing data
281 IF (.NOT. ALLOCATED(modeloutputdata)) ALLOCATE (modeloutputdata(0:readlinesmetdata, maxncols_cmod, numberofgrids)) !Data at model timestep
282 IF (.NOT. ALLOCATED(dataoutsuews)) ALLOCATE (dataoutsuews(readlinesmetdata, ncolumnsdataoutsuews, numberofgrids)) !Main output array
283 dataoutsuews = nan ! initialise Main output array
284 IF (.NOT. ALLOCATED(dataoutrsl)) ALLOCATE (dataoutrsl(readlinesmetdata, ncolumnsdataoutrsl, numberofgrids)) !RSL output array
285 dataoutrsl = nan ! initialise RSL array
286 IF (.NOT. ALLOCATED(dataoutdebug)) ALLOCATE (dataoutdebug(readlinesmetdata, ncolumnsdataoutdebug, numberofgrids)) !RSL output array
287 dataoutdebug = nan ! initialise Debug array
288 IF (.NOT. ALLOCATED(dataoutspartacus)) ALLOCATE (dataoutspartacus(readlinesmetdata, ncolumnsdataoutspartacus, numberofgrids)) !SPARTACUS output array
289 dataoutspartacus = nan ! initialise SPARTACUS array
290 IF (.NOT. ALLOCATED(dataoutdailystate)) ALLOCATE (dataoutdailystate(ndays, ncolumnsdataoutdailystate, numberofgrids)) !DailyState array
291 dataoutdailystate = nan ! initialise DailyState
292 IF (.NOT. ALLOCATED(dataoutbeers)) ALLOCATE (dataoutbeers(readlinesmetdata, ncolumnsdataoutbeers, numberofgrids)) !SOLWEIG POI output
294 IF (cbluse >= 1) ALLOCATE (dataoutbl(readlinesmetdata, ncolumnsdataoutbl, numberofgrids)) !CBL output
295 IF (.NOT. ALLOCATED(dataoutsnow)) ALLOCATE (dataoutsnow(readlinesmetdata, ncolumnsdataoutsnow, numberofgrids)) !Snow output
296
297 IF (.NOT. ALLOCATED(tsfc_surf_grids)) ALLOCATE (tsfc_surf_grids(numberofgrids, nsurf))
298 IF (.NOT. ALLOCATED(tin_surf_grids)) ALLOCATE (tin_surf_grids(numberofgrids, nsurf))
299 IF (.NOT. ALLOCATED(qn_s_av_grids)) ALLOCATE (qn_s_av_grids(numberofgrids))
300 IF (.NOT. ALLOCATED(dqnsdt_grids)) ALLOCATE (dqnsdt_grids(numberofgrids))
301 qn_s_av_grids = 0 ! Initialise to 0
302 dqnsdt_grids = 0 ! Initialise to 0
303
304 ! IF (StorageHeatMethod==4 .OR. StorageHeatMethod==14) THEN
305 IF (.NOT. ALLOCATED(dataoutestm)) ALLOCATE (dataoutestm(readlinesmetdata, ncolumnsdataoutestm, numberofgrids)) !ESTM output
306 IF (.NOT. ALLOCATED(dataoutestmext)) ALLOCATE (dataoutestmext(readlinesmetdata, ncolumnsdataoutestmext, numberofgrids)) !ESTM output
307 ! ENDIF
308
309 IF (.NOT. ALLOCATED(tair_av_grids)) ALLOCATE (tair_av_grids(numberofgrids))
310 IF (.NOT. ALLOCATED(qn_av_grids)) ALLOCATE (qn_av_grids(numberofgrids))
311 IF (.NOT. ALLOCATED(dqndt_grids)) ALLOCATE (dqndt_grids(numberofgrids))
312 !! QUESTION: Add snow clearing (?)
313 tair_av_grids = 273.15 ! Initialise to 273.15 K
314 qn_av_grids = 0 ! Initialise to 0
315 dqndt_grids = 0 ! Initialise to 0
316
317 IF (.NOT. ALLOCATED(qhforcbl)) ALLOCATE (qhforcbl(numberofgrids))
318 IF (.NOT. ALLOCATED(qeforcbl)) ALLOCATE (qeforcbl(numberofgrids))
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))
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 END IF
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)), &
358 END IF
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
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 END IF
382
383 ! Find number of lines in ESTM file
384 nlinesestmdata = 0 !Initialise nlinesMetdata (total number of lines in met forcing file)
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)), &
392 END IF
393
394 END IF
395
396 END IF
398 IF (.NOT. ALLOCATED(ts5mindata)) ALLOCATE (ts5mindata(1:readlinesmetdata, ncolsestmdata))
399 IF (.NOT. ALLOCATED(ts5mindata_ir)) ALLOCATE (ts5mindata_ir(ncolsestmdata))
400 IF (.NOT. ALLOCATED(tair24hr)) ALLOCATE (tair24hr(24*nsh))
401 ! ------------------------------------------------------------------------
402
403 !-----------------------------------------------------------------------
404 skippedlines = 0 !Initialise lines to be skipped in met forcing file
405 skippedlinesorig = 0 !Initialise lines to be skipped in original met forcing file
406 skippedlinesorigestm = 0 !Initialise lines to be skipped in original met forcing file
407
408 DO iblock = 1, readblocksmetdata !Loop through blocks of met data
409 ! WRITE(*,*) iblock,'/',ReadBlocksMetData
410
411 ! Model calculations are made in two stages:
412 ! (1) initialise the run for each block of met data (iblock from 1 to ReadBlocksMetData)
413 ! (2) perform the actual model calculations (SUEWS_Calculations)
414
415 gridcounter = 1 !Initialise counter for grids in each year
416 DO igrid = 1, numberofgrids !Loop through grids
417
418 gridid = grididmatrix(igrid) !store grid here for referencing error codes
419 WRITE (grid_txt, '(I10)') grididmatrix(igrid) !Get grid ID as a text string
420
421 ! (1) First stage: initialise run if this is the first iteration this year
422 ! (1a) Initialise surface characteristics
423 IF (iblock == 1) THEN
424 IF (diagnose == 1) WRITE (*, *) 'First block of data - doing initialisation'
425 ! (a) Transfer characteristics from SiteSelect to correct row of SurfaceChar
426 DO rr = 1, nlinessiteselect
427 !Find correct grid and year
428 IF (diagnose == 1) WRITE (*, *) 'grid found:', siteselect(rr, c_grid), 'grid needed:', grididmatrix(igrid)
429 IF (diagnose == 1) WRITE (*, *) 'year found:', siteselect(rr, c_year), 'year needed:', year_int
430 IF (siteselect(rr, c_grid) == grididmatrix(igrid) .AND. siteselect(rr, c_year) == year_int) THEN
431 IF (diagnose == 1) WRITE (*, *) 'Match found (grid and year) for rr = ', rr, 'of', nlinessiteselect
433 EXIT
434 ELSEIF (rr == nlinessiteselect) THEN
435 WRITE (*, *) 'Program stopped! Year', year_int, 'and/or grid', igrid, 'not found in SiteSelect.txt.'
436 CALL errorhint(59, 'Cannot find year and/or grid in SiteSelect.txt', real(igrid, kind(1d0)), notused, year_int)
437 END IF
438 END DO
439 END IF !end first block of met data
440
441 ! adjust ReadLinesMetdata for the last block of met data
442 IF (iblock == readblocksmetdata) THEN !For last block of data in file
443 readlinesmetdata_read = nlinesmetdata - (iblock - 1)*readlinesmetdata
444 ELSE
445 readlinesmetdata_read = readlinesmetdata
446 END IF
447
448 IF (igrid == 1) THEN
449 print *, 'Read in', readlinesmetdata_read, 'lines of met data in block', iblock, '/', readblocksmetdata
450 END IF
451
452 ! (1b) Initialise met data
453 IF (npertstepin > 1) THEN
454 ! Disaggregate met data ---------------------------------------------------
455
456 ! Set maximum value for ReadLinesOrigMetData to handle end of file (i.e. small final block)
457 IF (iblock == readblocksmetdata) THEN !For last block of data in file
459 ELSE
461 END IF
462 !write(*,*) ReadLinesOrigMetDataMAX, ReadLinesOrigMetData
463 ! Get names of original met forcing file(s) to disaggregate
464 ! Get met file name for this grid: SSss_YYYY_data_RR.txt
465 IF (multiplemetfiles == 1) THEN !If each grid has its own met file
466 fileorigmet = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)//'_data_' &
467 //trim(adjustl(resin_txt))//'.txt'
468 ! Also set file name for downscaled file
469 filedscdmet = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)//'_data_' &
470 //trim(adjustl(tstep_txt))//'.txt'
471 ! Disaggregate met data
472 CALL disaggregatemet(iblock, igrid)
473 ELSE
474 ! If each grid has the same met file, met file name does not include grid number, and only need to disaggregate once
475 fileorigmet = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_data_' &
476 //trim(adjustl(resin_txt))//'.txt'
477 filedscdmet = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_data_' &
478 //trim(adjustl(tstep_txt))//'.txt'
479 IF (igrid == 1) THEN !Disaggregate for the first grid only
480 CALL disaggregatemet(iblock, igrid)
481 ELSE !Then for subsequent grids simply copy data
482 metforcingdata(1:readlinesmetdata_read, 1:24, gridcounter) = metforcingdata(1:readlinesmetdata_read, 1:24, 1)
483 END IF
484 END IF
485
486 ELSEIF (npertstepin == 1) THEN
487 ! Get met forcing file name for this year for the first grid
488 ! Can be something else than 1
489 filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
490 filecodexwg = trim(filecode)//'_'//trim(year_txt) !File code without grid
491 ! IF(iblock==1) WRITE(*,*) 'Current FileCode: ', FileCodeX
492
493 ! For every block of met data ------------------------------------
494 ! Initialise met forcing data into 3-dimensional matrix
495 !write(*,*) 'Initialising met data for block',iblock
496 IF (multiplemetfiles == 1) THEN !If each grid has its own met file
497 filemet = trim(fileinputpath)//trim(filecodex)//'_data_'//trim(adjustl(tstep_txt))//'.txt'
498 CALL suews_initializemetdata(1, readlinesmetdata_read)
499 ELSE !If one met file used for all grids
500 !FileMet=TRIM(FileInputPath)//TRIM(FileCodeX)//'_data_'//TRIM(ADJUSTL(tstep_txt))//'.txt'
501 ! If one met file used for all grids, look for met file with no grid code (FileCodeXWG)
502 filemet = trim(fileinputpath)//trim(filecodexwg)//'_data_'//trim(adjustl(tstep_txt))//'.txt'
503 IF (igrid == 1) THEN !Read for the first grid only
504 CALL suews_initializemetdata(1, readlinesmetdata_read)
505 ELSE !Then for subsequent grids simply copy data
506 metforcingdata(1:readlinesmetdata_read, 1:24, gridcounter) = metforcingdata(1:readlinesmetdata_read, 1:24, 1)
507 END IF
508 END IF
509 END IF !end of nper statement
510
511 ! Only for the first block of met data, read initial conditions (moved from above, HCW 12 Jan 2017)
512 IF (iblock == 1) THEN
513 filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
514 !write(*,*) ' Now calling InitialState'
515 CALL initialstate(filecodex, year_int, gridcounter, numberofgrids)
516 END IF
517
518 ! Initialise ESTM if required, TS 05 Jun 2016; moved inside grid loop HCW 27 Jun 2016
519 IF (storageheatmethod == 4 .OR. storageheatmethod == 14) THEN
520 IF (nperestm > 1) THEN
521 ! Disaggregate ESTM data --------------------------------------------------
522 ! Set maximum value for ReadLinesOrigESTMData to handle end of file (i.e. small final block)
523 IF (iblock == readblocksmetdata) THEN !For last block of data in file
525 ELSE
527 END IF
528 !write(*,*) ReadLinesOrigESTMDataMAX, ReadLinesOrigESTMData
529 ! Get names of original ESTM forcing file(s) to disaggregate
530 ! Get ESTM file name for this grid: SSss_YYYY_ESTM_Ts_data_RR.txt
531 IF (multipleestmfiles == 1) THEN !If each grid has its own ESTM file
532 fileorigestm = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt) &
533 //'_ESTM_Ts_data_'//trim(adjustl(resinestm_txt))//'.txt'
534 ! Also set file name for downscaled file
535 filedscdestm = trim(fileinputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt) &
536 //'_ESTM_Ts_data_'//trim(adjustl(tstep_txt))//'.txt'
537 ! Disaggregate ESTM data
538 CALL disaggregateestm(iblock)
539 ELSE
540 ! If each grid has the same ESTM file, ESTM file name does not include grid number, and only need to disaggregate once
541 fileorigestm = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_ESTM_Ts_data_' &
542 //trim(adjustl(resinestm_txt))//'.txt'
543 filedscdestm = trim(fileinputpath)//trim(filecode)//'_'//trim(year_txt)//'_ESTM_Ts_data_' &
544 //trim(adjustl(tstep_txt))//'.txt'
545 IF (igrid == 1) THEN !Disaggregate for the first grid only
546 CALL disaggregateestm(iblock)
547 ELSE !Then for subsequent grids simply copy data
549 1:ncolsestmdata, 1)
550 END IF
551 END IF
552
553 ELSEIF (nperestm == 1) THEN
554 ! Get ESTM forcing file name for this year for the first grid
555 filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
556 filecodexwg = trim(filecode)//'_'//trim(year_txt) !File code without grid
557 ! For every block of ESTM data ------------------------------------
558 ! Initialise ESTM forcing data into 3-dimensional matrix
559 !write(*,*) 'Initialising ESTM data for block',iblock
560 IF (multipleestmfiles == 1) THEN !If each grid has its own met file
561 fileestmts = trim(fileinputpath)//trim(filecodex)//'_ESTM_Ts_data_'//trim(adjustl(tstep_txt))//'.txt'
562 !write(*,*) 'Calling GetESTMData...', FileCodeX, iblock, igrid
563 CALL suews_getestmdata(101)
564 ELSE !If one ESTM file used for all grids
565 fileestmts = trim(fileinputpath)//trim(filecodexwg)//'_ESTM_Ts_data_'//trim(adjustl(tstep_txt))//'.txt'
566 !write(*,*) 'Calling GetESTMData...', FileCodeX, iblock, igrid
567 IF (igrid == 1) THEN !Read for the first grid only
568 CALL suews_getestmdata(101)
569 ELSE !Then for subsequent grids simply copy data
571 1:ncolsestmdata, 1)
572 END IF
573 END IF
574 END IF !end of nperESTM statement
575 END IF
576
577 gridcounter = gridcounter + 1 !Increase GridCounter by 1 for next grid
578
579 END DO !end loop over grids
580 skippedlines = skippedlines + readlinesmetdata !Increase skippedLines ready for next block
581 skippedlinesorig = skippedlinesorig + readlinesorigmetdata !Increase skippedLinesOrig ready for next block
582 skippedlinesorigestm = skippedlinesorigestm + readlinesorigestmdata !Increase skippedLinesOrig ready for next block
583 !write(*,*) iblock
584 !write(*,*) ReadLinesMetdata, readlinesorigmetdata
585 !write(*,*) skippedLines, skippedLinesOrig, skippedLinesOrig*Nper
586
587 ! Initialise the modules on the first day
588 ! Initialise CBL if required
589 IF (iblock == 1) THEN
590 IF ((cbluse == 1) .OR. (cbluse == 2)) CALL cbl_readinputdata(fileinputpath, qh_choice)
591 END IF
592
593 ! First stage: initialisation done ----------------------------------
594
595 ! (2) Second stage: do calculations at 5-min time-steps -------------
596 ! First set maximum value of ir
597 IF (iblock == readblocksmetdata) THEN !For last block of data in file
598 irmax = nlinesmetdata - (iblock - 1)*readlinesmetdata
599 ELSE
600 irmax = readlinesmetdata
601 END IF
602
603 DO ir = 1, irmax !Loop through rows of current block of met data
604 ! GridCounter = 1 !Initialise counter for grids in each year
605 ! PRINT *, '*****************************************'
606 ! WRITE (*, *) 'ir here', ir, 'of', irMax, 'for block', iblock, 'of', ReadBlocksMetData
607 ! PRINT *, ''
608
609 ! quick stop : for testing
610 ! if ( ir>10 ) then
611 ! STOP 'testing finished'
612 ! end if
613
614 DO igrid = 1, numberofgrids !Loop through grids
615 IF (diagnose == 1) WRITE (*, *) 'Row (ir):', ir, '/', irmax, 'of block (iblock):', iblock, '/', readblocksmetdata, &
616 'Grid:', grididmatrix(igrid)
617
618 ! Call model calculation code
619 WRITE (grid_txt, '(I10)') grididmatrix(igrid) !Get grid ID as a text string
620 filecodex = trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(year_txt)
621 IF (ir == 1 .AND. igrid == 1) THEN
622 WRITE (*, *) trim(adjustl(filecodex)), &
623 ': Now running block ', iblock, '/', readblocksmetdata, ' of ', trim(year_txt), '...'
624 END IF
625
626 IF (diagnose == 1) WRITE (*, *) 'Calling SUEWS_Calculations...'
627 CALL suews_calculations(igrid, ir, iblock, irmax)
628 IF (diagnose == 1) WRITE (*, *) 'SUEWS_Calculations finished...'
629
630 ! Record iy and id for current time step to handle last row in yearly files (YYYY 1 0 0)
631 IF (igrid == numberofgrids) THEN !Adjust only when the final grid has been run for this time step
632 iy_prev_t = iy
633 id_prev_t = id
634 END IF
635
636 ! Write state information to new InitialConditions files
637 IF (ir == irmax) THEN !If last row...
638 IF (iblock == readblocksmetdata) THEN !...of last block of met data
639 WRITE (grid_txt, '(I10)') grididmatrix(igrid)
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 END IF
645 END IF
646
647 ! GridCounter = GridCounter + 1 !Increase GridCounter by 1 for next grid
648 END DO !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 END DO !end loop over rows of met data
656
657 ! Write output files in blocks --------------------------------
658 DO igrid = 1, numberofgrids
659 IF (diagnose == 1) WRITE (*, *) 'Calling SUEWS_Output...'
660 CALL suews_output(irmax, iblock, igrid, year_int)
661 END DO
662
663 END DO !end loop over blocks of met data
664 !-----------------------------------------------------------------------
665
666 ! ---- Decallocate arrays ----------------------------------------------
667 IF (diagnose == 1) WRITE (*, *) 'Deallocating arrays in SUEWS_Program.f95...'
668 DEALLOCATE (surfacechar)
669 DEALLOCATE (metforcingdata)
670 DEALLOCATE (metforcingdata_grid)
671 DEALLOCATE (modeloutputdata)
672 DEALLOCATE (dataoutsuews)
673 DEALLOCATE (dataoutrsl)
674 DEALLOCATE (dataoutdebug)
675 DEALLOCATE (dataoutspartacus)
676 DEALLOCATE (dataoutestmext)
677 DEALLOCATE (dataoutbeers)
678 DEALLOCATE (dataoutdailystate)
679 ! IF (SnowUse == 1) THEN
680 DEALLOCATE (dataoutsnow)
681 DEALLOCATE (tsfc_surf_grids)
682 DEALLOCATE (tin_surf_grids)
683 DEALLOCATE (qn_s_av_grids)
684 DEALLOCATE (dqnsdt_grids)
685 ! ENDIF
686 ! IF (StorageHeatMethod==4 .OR. StorageHeatMethod==14) THEN
687 DEALLOCATE (dataoutestm) !ESTM output
688 DEALLOCATE (estmforcingdata)
689 DEALLOCATE (ts5mindata)
690 DEALLOCATE (ts5mindata_ir)
691 DEALLOCATE (tair24hr)
692 ! ENDIF
693 DEALLOCATE (qhforcbl)
694 DEALLOCATE (qeforcbl)
695 DEALLOCATE (tair_av_grids)
696 DEALLOCATE (qn_av_grids)
697 DEALLOCATE (dqndt_grids)
698 IF (cbluse >= 1) THEN
699 DEALLOCATE (dataoutbl)
700 END IF
701 ! ----------------------------------------------------------------------
702
703 END DO !end loop over years
704
705 ! ---- Decallocate array --------------------------------------------------
706 ! Daily state needs to be outside year loop to transfer states between years
707 IF (ALLOCATED(modeldailystate)) DEALLOCATE (modeldailystate)
708 ! Also needs to happen at the end of the run
709 IF (ALLOCATED(usecolumnsdataout)) DEALLOCATE (usecolumnsdataout)
710
712 ! CALL spartacus_finalise
713 ! -------------------------------------------------------------------------
714
715 ! get cpu time consumed
716 CALL cpu_time(timefinish)
717 WRITE (*, *) "Time = ", timefinish - timestart, " seconds."
718
719 !Write to problems.txt that run has completed
720 IF (errorchoice == 0) THEN !if file has not been opened previously
721 OPEN (500, file='problems.txt')
722 errorchoice = 1
723 ELSE
724 OPEN (500, file='problems.txt', position="append")
725 END IF
726 !Writing of the problem file
727 WRITE (500, *) '--------------'
728 WRITE (500, *) 'Run completed.'
729 WRITE (500, *) '0' ! Write out error code 0 if run completed
730 CLOSE (500)
731
732 ! Also print to screen
733 WRITE (*, *) "----- SUEWS run completed -----"
734
735 stop 'finished'
736
737 ! 313 CALL errorHint(11,TRIM(FileOrigMet),notUsed,notUsed,ios_out)
738 ! 314 CALL errorHint(11,TRIM(FileMet),notUsed,notUsed,ios_out)
739 ! 315 CALL errorHint(11,TRIM(fileESTMTs),notUsed,notUsed,NotUsedI)
740
741END PROGRAM suews_program
real(kind(1d0)), dimension(:), allocatable dqndt_grids
real(kind(1d0)), dimension(:, :), allocatable modeldailystate
real(kind(1d0)), dimension(:, :, :), allocatable dataoutspartacus
real(kind(1d0)), dimension(:, :), allocatable siteselect
real(kind(1d0)), dimension(:, :), allocatable ts5mindata
real(kind(1d0)), dimension(:), allocatable dailystatefirstopen
real(kind(1d0)), dimension(:, :), allocatable tin_surf_grids
integer, parameter ncolsestmdata
integer, parameter maxncols_cmod
real(kind(1d0)), dimension(:, :), allocatable metforcingdata_grid
real(kind(1d0)), dimension(:, :, :), allocatable dataoutestmext
integer, parameter maxnumberofgrids
integer, parameter ncolumnsdataoutbl
real(kind(1d0)), dimension(:), allocatable dqnsdt_grids
integer, parameter ncolumnsmetforcingdata
real(kind(1d0)), dimension(:), allocatable qn_s_av_grids
real(kind(1d0)), dimension(:, :, :), allocatable dataoutsnow
integer, parameter ncolumnsdataoutsuews
integer, dimension(:), allocatable grididmatrix0
integer, dimension(:), allocatable usecolumnsdataout
integer, parameter ncolumnsdataoutestm
real(kind(1d0)), dimension(:, :, :), allocatable dataoutdebug
real(kind(1d0)), dimension(:, :), allocatable tsfc_surf_grids
integer, parameter maxncols_c
real(kind(1d0)), dimension(:), allocatable tair24hr
integer, parameter ncolumnsdataoutbeers
real(kind(1d0)), dimension(:, :, :), allocatable metforcingdata
integer, parameter ncolumnsdataoutspartacus
integer, parameter ncolumnsdataoutrsl
real(kind(1d0)), dimension(:, :), allocatable surfacechar
integer, dimension(:), allocatable grididmatrix
integer, parameter ncolumnsdataoutestmext
integer, parameter ncolumnsdataoutsnow
real(kind(1d0)), dimension(:, :, :), allocatable dataoutrsl
real(kind(1d0)), dimension(:, :, :), allocatable dataoutsuews
integer, parameter ncolumnsdataoutdailystate
integer, parameter maxncols_cmds
real(kind(1d0)), dimension(:), allocatable ts5mindata_ir
integer, parameter nsurf
integer, parameter ndays
real(kind(1d0)), dimension(:, :, :), allocatable dataoutbl
real(kind(1d0)), dimension(:, :, :), allocatable modeloutputdata
real(kind(1d0)), dimension(:, :, :), allocatable dataoutdailystate
real(kind(1d0)), dimension(:, :, :), allocatable dataoutestm
real(kind(1d0)), dimension(:), allocatable tair_av_grids
real(kind(1d0)), dimension(:, :, :), allocatable estmforcingdata
real(kind(1d0)), dimension(:, :, :), allocatable dataoutbeers
integer, parameter maxlinesmet
real(kind(1d0)), dimension(:), allocatable qn_av_grids
integer, parameter ncolumnsdataoutdebug
subroutine cbl_readinputdata(FileInputPath, qh_choice)
subroutine suews_output(irMax, iv, Gridiv, iyr)
integer function count_lines(filename)
character(len=150) filedscdestm
integer netradiationmethod
integer resolutionfilesin
character(len=20) filecode
integer resolutionfilesinestm
integer multiplemetfiles
character(len=150) fileestmts
character(len=90) progname
integer diagnose
character(len=150) filemet
integer storageheatmethod
character(len=150) fileinputpath
character(len=150) filedscdmet
integer multipleestmfiles
character(len=150) fileorigestm
character(len=150) fileorigmet
integer outputformats
real(kind(1d0)) nan
real(kind(1d0)) notused
subroutine estm_initials
subroutine suews_getestmdata(lunit)
subroutine estm_ext_finalise
subroutine estm_ext_initialise
integer skippedlinesorig
integer readlinesmetdata
integer gridcounter
integer nlinessiteselect
integer readblocksmetdata
integer nlinesestmdata
integer nlinesorigmetdata
integer numberofgrids
integer nlinesorigestmdata
integer lastyear
integer firstyear
integer skippedlines
integer readlinesorigestmdatamax
integer readlinesorigmetdata
integer readlinesorigmetdatamax
integer readlinesorigestmdata
integer nlinesmetdata
integer skippedlinesorigestm
integer readblocksorigmetdata
subroutine disaggregateestm(iBlock)
subroutine disaggregatemet(iBlock, igrid)
real(kind(1d0)), dimension(:), allocatable qeforcbl
real(kind(1d0)), dimension(:), allocatable qhforcbl
real(kind(1d0)) nperestm_real
real(kind(1d0)) npertstepin_real
integer dt_since_start
integer id_prev_t
integer nofdaysthisyear
integer iy_prev_t
integer iy
integer id
character(len=90) git_commit
character(len=90) compiler_ver
character(len=10) gridid_text
character(len=15) datetime
subroutine suews_calculations(Gridiv, ir, iMB, irMax)
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
subroutine overallruncontrol
subroutine initializesurfacecharacteristics(Gridiv, rr)
subroutine initialstate(GridName, year_int, Gridiv, NumberOfGrids)
subroutine suews_initializemetdata(lunit, ReadlinesMetdata_read)
subroutine nextinitial(GridName, year_int)
program suews_program
subroutine leapyearcalc(year_int, nroDays)