SUEWS API Site
Documentation of SUEWS source code
suews_ctrl_output.f95
Go to the documentation of this file.
1 MODULE ctrl_output
2  !===========================================================================================
3  ! generic output functions for SUEWS
4  ! authors: Ting Sun (ting.sun@reading.ac.uk)
5  !
6  ! disclamier:
7  ! This code employs the netCDF Fortran 90 API.
8  ! Full documentation of the netCDF Fortran 90 API can be found at:
9  ! https://www.unidata.ucar.edu/software/netcdf/netcdf-4/newdocs/netcdf-f90/
10  ! Part of the work is under the help of examples provided by the documentation.
11  !
12  ! purpose:
13  ! these subroutines write out the results of SUEWS in netCDF format.
14  !
15  !
16  ! history:
17  ! TS 20161209: initial version of netcdf function
18  ! TS 20161213: standalise the txt2nc procedure
19  ! TS 20170414: generic output procedures
20  ! TS 20171016: added support for DailyState
21  ! TS 20171017: combined txt and nc wrappers into one: reduced duplicate code at two places
22  !===========================================================================================
23 
24  USE allocatearray
25  USE cbl_module
26  USE data_in
27  ! USE defaultNotUsed
28  ! USE ESTM_data
29  USE gis_data
30  ! USE initial
31  USE sues_data
32  USE time
33  USE strings
34 
35  IMPLICIT NONE
36 
37  INTEGER :: n
38 
39  CHARACTER(len=10), PARAMETER:: & !Define useful formats here
40  fy = 'i0004,1X', & !4 digit integer for year
41  ft = 'i0004,1X', & !3 digit integer for id, it, imin
42  fd = 'f08.4,1X', & !3 digits + 4 dp for dectime
43  f94 = 'f09.4,1X', & !standard output format: 4 dp + 4 digits
44  f104 = 'f10.4,1X', & !standard output format: 4 dp + 5 digits
45  f106 = 'f10.6,1X', & !standard output format: 6 dp + 3 digits
46  f146 = 'f14.6,1X' !standard output format: 6 dp + 7 digits
47 
48  CHARACTER(len=1), PARAMETER:: & ! Define aggregation methods here
49  at = 'T', & !time columns
50  aa = 'A', & !average
51  as = 'S', & !sum
52  al = 'L' !last value
53 
54  CHARACTER(len=3):: itext
55 
56  ! define type: variable attributes
57  TYPE varattr
58  CHARACTER(len=15) :: header ! short name in headers
59  CHARACTER(len=12) :: unit ! unit
60  CHARACTER(len=10) :: fmt ! output format
61  CHARACTER(len=100) :: longnm ! long name for detailed description
62  CHARACTER(len=1) :: aggreg ! aggregation method
63  CHARACTER(len=10) :: group ! group: datetime, default, ESTM, Snow, etc.
64  INTEGER :: level ! output priority level: 0 for highest (defualt output)
65  END TYPE varattr
66 
67  ! initialise valist
68  TYPE(varattr) :: varlistall(500)
69 
70  ! datetime:
71  DATA(varlistall(n), n=1, 5)/ &
72  varattr('Year', 'YYYY', fy, 'Year', at, 'datetime', 0), &
73  varattr('DOY', 'DOY', ft, 'Day of Year', at, 'datetime', 0), &
74  varattr('Hour', 'HH', ft, 'Hour', at, 'datetime', 0), &
75  varattr('Min', 'MM', ft, 'Minute', at, 'datetime', 0), &
76  varattr('Dectime', '-', fd, 'Decimal time', at, 'datetime', 0) &
77  /
78 
79  ! defualt:
80  DATA(varlistall(n), &
81  n=5 + 1, &
83  varattr('Kdown', 'W m-2', f104, 'Incoming shortwave radiation', aa, 'SUEWS', 0), &
84  varattr('Kup', 'W m-2', f104, 'Outgoing shortwave radiation', aa, 'SUEWS', 0), &
85  varattr('Ldown', 'W m-2', f104, 'Incoming longwave radiation', aa, 'SUEWS', 0), &
86  varattr('Lup', 'W m-2', f104, 'Outgoing longwave radiation', aa, 'SUEWS', 0), &
87  varattr('Tsurf', 'degC', f104, 'Bulk surface temperature', aa, 'SUEWS', 0), &
88  varattr('QN', 'W m-2', f104, 'Net all-wave radiation', aa, 'SUEWS', 0), &
89  varattr('QF', 'W m-2', f104, 'Anthropogenic heat flux', aa, 'SUEWS', 0), &
90  varattr('QS', 'W m-2', f104, 'Net storage heat flux', aa, 'SUEWS', 0), &
91  varattr('QH', 'W m-2', f104, 'Sensible heat flux', aa, 'SUEWS', 0), &
92  varattr('QE', 'W m-2', f104, 'Latent heat flux', aa, 'SUEWS', 0), &
93  varattr('QHlumps', 'W m-2', f104, 'Sensible heat flux (using LUMPS)', aa, 'SUEWS', 1), &
94  varattr('QElumps', 'W m-2', f104, 'Latent heat flux (using LUMPS)', aa, 'SUEWS', 1), &
95  varattr('QHresis', 'W m-2', f104, 'Sensible heat flux (resistance method)', aa, 'SUEWS', 1), &
96  varattr('Rain', 'mm', f104, 'Rain', as, 'SUEWS', 0), &
97  varattr('Irr', 'mm', f104, 'Irrigation', as, 'SUEWS', 0), &
98  varattr('Evap', 'mm', f104, 'Evaporation', as, 'SUEWS', 0), &
99  varattr('RO', 'mm', f104, 'Runoff', as, 'SUEWS', 0), &
100  varattr('TotCh', 'mm', f146, 'Surface and soil moisture change', as, 'SUEWS', 0), &
101  varattr('SurfCh', 'mm', f146, 'Surface moisture change', as, 'SUEWS', 0), &
102  varattr('State', 'mm', f104, 'Surface Wetness State', al, 'SUEWS', 0), &
103  varattr('NWtrState', 'mm', f104, 'Surface wetness state (non-water surfaces)', al, 'SUEWS', 0), &
104  varattr('Drainage', 'mm', f104, 'Drainage', as, 'SUEWS', 0), &
105  varattr('SMD', 'mm', f94, 'Soil Moisture Deficit', al, 'SUEWS', 0), &
106  varattr('FlowCh', 'mm', f104, 'Additional flow into water body', as, 'SUEWS', 1), &
107  varattr('AddWater', 'mm', f104, 'Addtional water from other grids', as, 'SUEWS', 1), &
108  varattr('ROSoil', 'mm', f104, 'Runoff to soil', as, 'SUEWS', 1), &
109  varattr('ROPipe', 'mm', f104, 'Runoff to pipes', as, 'SUEWS', 1), &
110  varattr('ROImp', 'mm', f104, 'Runoff over impervious surfaces', as, 'SUEWS', 1), &
111  varattr('ROVeg', 'mm', f104, 'Runoff over vegetated surfaces', as, 'SUEWS', 1), &
112  varattr('ROWater', 'mm', f104, 'Runoff for water surface', as, 'SUEWS', 1), &
113  varattr('WUInt', 'mm', f94, 'InternalWaterUse', as, 'SUEWS', 1), &
114  varattr('WUEveTr', 'mm', f94, 'Water use for evergreen trees', as, 'SUEWS', 1), &
115  varattr('WUDecTr', 'mm', f94, 'Water use for deciduous trees', as, 'SUEWS', 1), &
116  varattr('WUGrass', 'mm', f94, 'Water use for grass', as, 'SUEWS', 1), &
117  varattr('SMDPaved', 'mm', f94, 'Soil moisture deficit for paved surface', al, 'SUEWS', 1), &
118  varattr('SMDBldgs', 'mm', f94, 'Soil moisture deficit for building surface', al, 'SUEWS', 1), &
119  varattr('SMDEveTr', 'mm', f94, 'Soil moisture deficit for evergreen tree surface', al, 'SUEWS', 1), &
120  varattr('SMDDecTr', 'mm', f94, 'Soil moisture deficit for deciduous tree surface', al, 'SUEWS', 1), &
121  varattr('SMDGrass', 'mm', f94, 'Soil moisture deficit for grass surface', al, 'SUEWS', 1), &
122  varattr('SMDBSoil', 'mm', f94, 'Soil moisture deficit for bare soil surface', al, 'SUEWS', 1), &
123  varattr('StPaved', 'mm', f94, 'Surface wetness state for paved surface', al, 'SUEWS', 1), &
124  varattr('StBldgs', 'mm', f94, 'Surface wetness state for building surface', al, 'SUEWS', 1), &
125  varattr('StEveTr', 'mm', f94, 'Surface wetness state for evergreen tree surface', al, 'SUEWS', 1), &
126  varattr('StDecTr', 'mm', f94, 'Surface wetness state for deciduous tree surface', al, 'SUEWS', 1), &
127  varattr('StGrass', 'mm', f94, 'Surface wetness state for grass surface', al, 'SUEWS', 1), &
128  varattr('StBSoil', 'mm', f94, 'Surface wetness state for bare soil surface', al, 'SUEWS', 1), &
129  varattr('StWater', 'mm', f104, 'Surface wetness state for water surface', al, 'SUEWS', 1), &
130  varattr('Zenith', 'degree', f104, 'Solar zenith angle', al, 'SUEWS', 0), &
131  varattr('Azimuth', 'degree', f94, 'Solar azimuth angle', al, 'SUEWS', 0), &
132  varattr('AlbBulk', '1', f94, 'Bulk albedo', aa, 'SUEWS', 0), &
133  varattr('Fcld', '1', f94, 'Cloud fraction', aa, 'SUEWS', 0), &
134  varattr('LAI', 'm2 m-2', f94, 'Leaf area index', aa, 'SUEWS', 0), &
135  varattr('z0m', 'm', f94, 'Roughness length for momentum', aa, 'SUEWS', 1), &
136  varattr('zdm', 'm', f94, 'Zero-plane displacement height', aa, 'SUEWS', 1), &
137  varattr('UStar', 'm s-1', f94, 'Friction velocity', aa, 'SUEWS', 0), &
138  varattr('Lob', 'm', f146, 'Obukhov length', aa, 'SUEWS', 0), &
139  varattr('RA', 's m-1', f104, 'Aerodynamic resistance', aa, 'SUEWS', 1), &
140  varattr('RS', 's m-1', f104, 'Surface resistance', aa, 'SUEWS', 1), &
141  varattr('Fc', 'umol m-2 s-1', f94, 'CO2 flux', aa, 'SUEWS', 0), &
142  varattr('FcPhoto', 'umol m-2 s-1', f94, 'CO2 flux from photosynthesis', aa, 'SUEWS', 1), &
143  varattr('FcRespi', 'umol m-2 s-1', f94, 'CO2 flux from respiration', aa, 'SUEWS', 1), &
144  varattr('FcMetab', 'umol m-2 s-1', f94, 'CO2 flux from metabolism', aa, 'SUEWS', 1), &
145  varattr('FcTraff', 'umol m-2 s-1', f94, 'CO2 flux from traffic', aa, 'SUEWS', 1), &
146  varattr('FcBuild', 'umol m-2 s-1', f94, 'CO2 flux from buildings', aa, 'SUEWS', 1), &
147  varattr('FcPoint', 'umol m-2 s-1', f94, 'CO2 flux from point source', aa, 'SUEWS', 1), &
148  varattr('QNSnowFr', 'W m-2', f94, 'Net all-wave radiation for non-snow area', aa, 'SUEWS', 2), &
149  varattr('QNSnow', 'W m-2', f94, 'Net all-wave radiation for snow area', aa, 'SUEWS', 2), &
150  varattr('AlbSnow', '-', f94, 'Snow albedo', aa, 'SUEWS', 2), &
151  varattr('QM', 'W m-2', f106, 'Snow-related heat exchange', aa, 'SUEWS', 2), &
152  varattr('QMFreeze', 'W m-2', f146, 'Internal energy change', aa, 'SUEWS', 2), &
153  varattr('QMRain', 'W m-2', f106, 'Heat released by rain on snow', aa, 'SUEWS', 2), &
154  varattr('SWE', 'mm', f104, 'Snow water equivalent', aa, 'SUEWS', 2), &
155  varattr('MeltWater', 'mm', f104, 'Meltwater', aa, 'SUEWS', 2), &
156  varattr('MeltWStore', 'mm', f104, 'Meltwater store', aa, 'SUEWS', 2), &
157  varattr('SnowCh', 'mm', f104, 'Change in snow pack', as, 'SUEWS', 2), &
158  varattr('SnowRPaved', 'mm', f94, 'Snow removed from paved surface', as, 'SUEWS', 2), &
159  varattr('SnowRBldgs', 'mm', f94, 'Snow removed from building surface', as, 'SUEWS', 2), &
160  varattr('Ts', 'degC', f94, 'Skin temperature', aa, 'SUEWS', 0), &
161  varattr('T2', 'degC', f94, 'Air temperature at 2 m', aa, 'SUEWS', 0), &
162  varattr('Q2', 'g kg-1', f94, 'Specific humidity at 2 m', aa, 'SUEWS', 0), &
163  varattr('U10', 'm s-1', f94, 'Wind speed at 10 m', aa, 'SUEWS', 0), &
164  varattr('RH2', '%', f94, 'Relative humidity at 2 m', aa, 'SUEWS', 0) &
165  /
166 
167  ! SOLWEIG:
168  DATA(varlistall(n), &
169  n=ncolumnsdataoutsuews + 1, &
171  varattr('azimuth', 'to_add', f106, 'azimuth', aa, 'SOLWEIG', 0), &
172  varattr('altitude', 'to_add', f106, 'altitude', aa, 'SOLWEIG', 0), &
173  varattr('GlobalRad', 'to_add', f106, 'GlobalRad', aa, 'SOLWEIG', 0), &
174  varattr('DiffuseRad', 'to_add', f106, 'DiffuseRad', aa, 'SOLWEIG', 0), &
175  varattr('DirectRad', 'to_add', f106, 'DirectRad', aa, 'SOLWEIG', 0), &
176  varattr('Kdown2d', 'to_add', f106, 'Kdown2d', aa, 'SOLWEIG', 0), &
177  varattr('Kup2d', 'to_add', f106, 'Kup2d', aa, 'SOLWEIG', 0), &
178  varattr('Ksouth', 'to_add', f106, 'Ksouth', aa, 'SOLWEIG', 0), &
179  varattr('Kwest', 'to_add', f106, 'Kwest', aa, 'SOLWEIG', 0), &
180  varattr('Knorth', 'to_add', f106, 'Knorth', aa, 'SOLWEIG', 0), &
181  varattr('Keast', 'to_add', f106, 'Keast', aa, 'SOLWEIG', 0), &
182  varattr('Ldown2d', 'to_add', f106, 'Ldown2d', aa, 'SOLWEIG', 0), &
183  varattr('Lup2d', 'to_add', f106, 'Lup2d', aa, 'SOLWEIG', 0), &
184  varattr('Lsouth', 'to_add', f106, 'Lsouth', aa, 'SOLWEIG', 0), &
185  varattr('Lwest', 'to_add', f106, 'Lwest', aa, 'SOLWEIG', 0), &
186  varattr('Lnorth', 'to_add', f106, 'Lnorth', aa, 'SOLWEIG', 0), &
187  varattr('Least', 'to_add', f106, 'Least', aa, 'SOLWEIG', 0), &
188  varattr('Tmrt', 'to_add', f106, 'Tmrt', aa, 'SOLWEIG', 0), &
189  varattr('I0', 'to_add', f106, 'I0', aa, 'SOLWEIG', 0), &
190  varattr('CI', 'to_add', f106, 'CI', aa, 'SOLWEIG', 0), &
191  varattr('gvf', 'to_add', f106, 'gvf', aa, 'SOLWEIG', 0), &
192  varattr('shadow', 'to_add', f106, 'shadow', aa, 'SOLWEIG', 0), &
193  varattr('svf', 'to_add', f106, 'svf', aa, 'SOLWEIG', 0), &
194  varattr('svfbuveg', 'to_add', f106, 'svfbuveg', aa, 'SOLWEIG', 0), &
195  varattr('Ta', 'to_add', f106, 'Ta', aa, 'SOLWEIG', 0), &
196  varattr('Tg', 'to_add', f106, 'Tg', aa, 'SOLWEIG', 0) &
197  /
198 
199  ! BL:
200  DATA(varlistall(n), &
203  varattr('z', 'to_add', f104, 'z', aa, 'BL', 0), &
204  varattr('theta', 'to_add', f104, 'theta', aa, 'BL', 0), &
205  varattr('q', 'to_add', f104, 'q', aa, 'BL', 0), &
206  varattr('theta+', 'to_add', f104, 'theta+', aa, 'BL', 0), &
207  varattr('q+', 'to_add', f104, 'q+', aa, 'BL', 0), &
208  varattr('Temp_C', 'to_add', f104, 'Temp_C', aa, 'BL', 0), &
209  varattr('rh', 'to_add', f104, 'rh', aa, 'BL', 0), &
210  varattr('QH_use', 'to_add', f104, 'QH_use', aa, 'BL', 0), &
211  varattr('QE_use', 'to_add', f104, 'QE_use', aa, 'BL', 0), &
212  varattr('Press_hPa', 'to_add', f104, 'Press_hPa', aa, 'BL', 0), &
213  varattr('avu1', 'to_add', f104, 'avu1', aa, 'BL', 0), &
214  varattr('UStar', 'to_add', f104, 'UStar', aa, 'BL', 0), &
215  varattr('avdens', 'to_add', f104, 'avdens', aa, 'BL', 0), &
216  varattr('lv_J_kg', 'to_add', f146, 'lv_J_kg', aa, 'BL', 0), &
217  varattr('avcp', 'to_add', f104, 'avcp', aa, 'BL', 0), &
218  varattr('gamt', 'to_add', f104, 'gamt', aa, 'BL', 0), &
219  varattr('gamq', 'to_add', f104, 'gamq', aa, 'BL', 0) &
220  /
221 
222  ! Snow:
223  DATA(varlistall(n), &
226  varattr('SWE_Paved', 'mm', f106, 'Snow water equivalent for paved surface', aa, 'snow', 0), &
227  varattr('SWE_Bldgs', 'mm', f106, 'Snow water equivalent for building surface', aa, 'snow', 0), &
228  varattr('SWE_EveTr', 'mm', f106, 'Snow water equivalent for evergreen tree surface', aa, 'snow', 0), &
229  varattr('SWE_DecTr', 'mm', f106, 'Snow water equivalent for deciduous tree surface', aa, 'snow', 0), &
230  varattr('SWE_Grass', 'mm', f106, 'Snow water equivalent for grass surface', aa, 'snow', 0), &
231  varattr('SWE_BSoil', 'mm', f106, 'Snow water equivalent for bare soil surface', aa, 'snow', 0), &
232  varattr('SWE_Water', 'mm', f106, 'Snow water equivalent for water surface', aa, 'snow', 0), &
233  varattr('Mw_Paved', 'mm', f106, 'Meltwater for paved surface', as, 'snow', 0), &
234  varattr('Mw_Bldgs', 'mm', f106, 'Meltwater for building surface', as, 'snow', 0), &
235  varattr('Mw_EveTr', 'mm', f106, 'Meltwater for evergreen tree surface', as, 'snow', 0), &
236  varattr('Mw_DecTr', 'mm', f106, 'Meltwater for deciduous tree surface', as, 'snow', 0), &
237  varattr('Mw_Grass', 'mm', f106, 'Meltwater for grass surface', as, 'snow', 0), &
238  varattr('Mw_BSoil', 'mm', f106, 'Meltwater for bare soil surface', as, 'snow', 0), &
239  varattr('Mw_Water', 'mm', f106, 'Meltwater for water surface', as, 'snow', 0), &
240  varattr('Qm_Paved', 'W m-2', f106, 'Snow-related heat exchange for paved surface', aa, 'snow', 0), &
241  varattr('Qm_Bldgs', 'W m-2', f106, 'Snow-related heat exchange for building surface', aa, 'snow', 0), &
242  varattr('Qm_EveTr', 'W m-2', f106, 'Snow-related heat exchange for evergreen tree surface', aa, 'snow', 0), &
243  varattr('Qm_DecTr', 'W m-2', f106, 'Snow-related heat exchange for deciduous tree surface', aa, 'snow', 0), &
244  varattr('Qm_Grass', 'W m-2', f106, 'Snow-related heat exchange for grass surface', aa, 'snow', 0), &
245  varattr('Qm_BSoil', 'W m-2', f106, 'Snow-related heat exchange for bare soil surface', aa, 'snow', 0), &
246  varattr('Qm_Water', 'W m-2', f106, 'Snow-related heat exchange for water surface', aa, 'snow', 0), &
247  varattr('Qa_Paved', 'W m-2', f106, 'Advective heat for paved surface', aa, 'snow', 0), &
248  varattr('Qa_Bldgs', 'W m-2', f106, 'Advective heat for building surface', aa, 'snow', 0), &
249  varattr('Qa_EveTr', 'W m-2', f106, 'Advective heat for evergreen tree surface', aa, 'snow', 0), &
250  varattr('Qa_DecTr', 'W m-2', f106, 'Advective heat for deciduous tree surface', aa, 'snow', 0), &
251  varattr('Qa_Grass', 'W m-2', f106, 'Advective heat for grass surface', aa, 'snow', 0), &
252  varattr('Qa_BSoil', 'W m-2', f106, 'Advective heat for bare soil surface', aa, 'snow', 0), &
253  varattr('Qa_Water', 'W m-2', f106, 'Advective heat for water surface', aa, 'snow', 0), &
254  varattr('QmFr_Paved', 'W m-2', f146, 'Heat related to freezing for paved surface', aa, 'snow', 0), &
255  varattr('QmFr_Bldgs', 'W m-2', f146, 'Heat related to freezing for building surface', aa, 'snow', 0), &
256  varattr('QmFr_EveTr', 'W m-2', f146, 'Heat related to freezing for evergreen tree surface', aa, 'snow', 0), &
257  varattr('QmFr_DecTr', 'W m-2', f146, 'Heat related to freezing for deciduous tree surface', aa, 'snow', 0), &
258  varattr('QmFr_Grass', 'W m-2', f146, 'Heat related to freezing for grass surface', aa, 'snow', 0), &
259  varattr('QmFr_BSoil', 'W m-2', f146, 'Heat related to freezing for bare soil surface', aa, 'snow', 0), &
260  varattr('QmFr_Water', 'W m-2', f146, 'Heat related to freezing for water surface', aa, 'snow', 0), &
261  varattr('fr_Paved', '1', f106, 'Fraction of snow for paved surface', aa, 'snow', 0), &
262  varattr('fr_Bldgs', '1', f106, 'Fraction of snow for building surface', aa, 'snow', 0), &
263  varattr('fr_EveTr', '1', f106, 'Fraction of snow for evergreen tree surface', aa, 'snow', 0), &
264  varattr('fr_DecTr', '1', f106, 'Fraction of snow for deciduous tree surface', aa, 'snow', 0), &
265  varattr('fr_Grass', '1', f106, 'Fraction of snow for grass surface', aa, 'snow', 0), &
266  varattr('fr_BSoil', '1', f106, 'Fraction of snow for bare soil surface', aa, 'snow', 0), &
267  varattr('RainSn_Paved', 'mm', f146, 'Rain on snow for paved surface', as, 'snow', 0), &
268  varattr('RainSn_Bldgs', 'mm', f146, 'Rain on snow for building surface', as, 'snow', 0), &
269  varattr('RainSn_EveTr', 'mm', f146, 'Rain on snow for evergreen tree surface', as, 'snow', 0), &
270  varattr('RainSn_DecTr', 'mm', f146, 'Rain on snow for deciduous tree surface', as, 'snow', 0), &
271  varattr('RainSn_Grass', 'mm', f146, 'Rain on snow for grass surface', as, 'snow', 0), &
272  varattr('RainSn_BSoil', 'mm', f146, 'Rain on snow for bare soil surface', as, 'snow', 0), &
273  varattr('RainSn_Water', 'mm', f146, 'Rain on snow for water surface', as, 'snow', 0), &
274  varattr('Qn_PavedSnow', 'W m-2', f146, 'Net all-wave radiation for snow paved surface', aa, 'snow', 0), &
275  varattr('Qn_BldgsSnow', 'W m-2', f146, 'Net all-wave radiation for snow building surface', aa, 'snow', 0), &
276  varattr('Qn_EveTrSnow', 'W m-2', f146, 'Net all-wave radiation for snow evergreen tree surface', aa, 'snow', 0), &
277  varattr('Qn_DecTrSnow', 'W m-2', f146, 'Net all-wave radiation for snow deciduous tree surface', aa, 'snow', 0), &
278  varattr('Qn_GrassSnow', 'W m-2', f146, 'Net all-wave radiation for snow grass surface', aa, 'snow', 0), &
279  varattr('Qn_BSoilSnow', 'W m-2', f146, 'Net all-wave radiation for snow bare soil surface', aa, 'snow', 0), &
280  varattr('Qn_WaterSnow', 'W m-2', f146, 'Net all-wave radiation for snow water surface', aa, 'snow', 0), &
281  varattr('kup_PavedSnow', 'W m-2', f146, 'Reflected shortwave radiation for snow paved surface', aa, 'snow', 0), &
282  varattr('kup_BldgsSnow', 'W m-2', f146, 'Reflected shortwave radiation for snow building surface', aa, 'snow', 0), &
283  varattr('kup_EveTrSnow', 'W m-2', f146, 'Reflected shortwave radiation for snow evergreen tree surface', aa, 'snow', 0), &
284  varattr('kup_DecTrSnow', 'W m-2', f146, 'Reflected shortwave radiation for snow deciduous tree surface', aa, 'snow', 0), &
285  varattr('kup_GrassSnow', 'W m-2', f146, 'Reflected shortwave radiation for snow grass surface', aa, 'snow', 0), &
286  varattr('kup_BSoilSnow', 'W m-2', f146, 'Reflected shortwave radiation for snow bare soil surface', aa, 'snow', 0), &
287  varattr('kup_WaterSnow', 'W m-2', f146, 'Reflected shortwave radiation for snow water surface', aa, 'snow', 0), &
288  varattr('frMelt_Paved', 'mm', f146, 'Amount of freezing melt water for paved surface', aa, 'snow', 0), &
289  varattr('frMelt_Bldgs', 'mm', f146, 'Amount of freezing melt water for building surface', aa, 'snow', 0), &
290  varattr('frMelt_EveTr', 'mm', f146, 'Amount of freezing melt water for evergreen tree surface', aa, 'snow', 0), &
291  varattr('frMelt_DecTr', 'mm', f146, 'Amount of freezing melt water for deciduous tree surface', aa, 'snow', 0), &
292  varattr('frMelt_Grass', 'mm', f146, 'Amount of freezing melt water for grass surface', aa, 'snow', 0), &
293  varattr('frMelt_BSoil', 'mm', f146, 'Amount of freezing melt water for bare soil surface', aa, 'snow', 0), &
294  varattr('frMelt_Water', 'mm', f146, 'Amount of freezing melt water for water surface', aa, 'snow', 0), &
295  varattr('MwStore_Paved', 'mm', f146, 'Meltwater store for paved surface', aa, 'snow', 0), &
296  varattr('MwStore_Bldgs', 'mm', f146, 'Meltwater store for building surface', aa, 'snow', 0), &
297  varattr('MwStore_EveTr', 'mm', f146, 'Meltwater store for evergreen tree surface', aa, 'snow', 0), &
298  varattr('MwStore_DecTr', 'mm', f146, 'Meltwater store for deciduous tree surface', aa, 'snow', 0), &
299  varattr('MwStore_Grass', 'mm', f146, 'Meltwater store for grass surface', aa, 'snow', 0), &
300  varattr('MwStore_BSoil', 'mm', f146, 'Meltwater store for bare soil surface', aa, 'snow', 0), &
301  varattr('MwStore_Water', 'mm', f146, 'Meltwater store for water surface', aa, 'snow', 0), &
302  varattr('DensSnow_Paved', 'kg m-3', f146, 'Snow density for paved surface', aa, 'snow', 0), &
303  varattr('DensSnow_Bldgs', 'kg m-3', f146, 'Snow density for building surface', aa, 'snow', 0), &
304  varattr('DensSnow_EveTr', 'kg m-3', f146, 'Snow density for evergreen tree surface', aa, 'snow', 0), &
305  varattr('DensSnow_DecTr', 'kg m-3', f146, 'Snow density for deciduous tree surface', aa, 'snow', 0), &
306  varattr('DensSnow_Grass', 'kg m-3', f146, 'Snow density for grass surface', aa, 'snow', 0), &
307  varattr('DensSnow_BSoil', 'kg m-3', f146, 'Snow density for bare soil surface', aa, 'snow', 0), &
308  varattr('DensSnow_Water', 'kg m-3', f146, 'Snow density for water surface', aa, 'snow', 0), &
309  varattr('Sd_Paved', 'mm', f106, 'Snow depth for paved surface', aa, 'snow', 0), &
310  varattr('Sd_Bldgs', 'mm', f106, 'Snow depth for building surface', aa, 'snow', 0), &
311  varattr('Sd_EveTr', 'mm', f106, 'Snow depth for evergreen tree surface', aa, 'snow', 0), &
312  varattr('Sd_DecTr', 'mm', f106, 'Snow depth for deciduous tree surface', aa, 'snow', 0), &
313  varattr('Sd_Grass', 'mm', f106, 'Snow depth for grass surface', aa, 'snow', 0), &
314  varattr('Sd_BSoil', 'mm', f106, 'Snow depth for bare soil surface', aa, 'snow', 0), &
315  varattr('Sd_Water', 'mm', f106, 'Snow depth for water surface', aa, 'snow', 0), &
316  varattr('Tsnow_Paved', 'degC', f146, 'Snow surface temperature for paved surface', aa, 'snow', 0), &
317  varattr('Tsnow_Bldgs', 'degC', f146, 'Snow surface temperature for building surface', aa, 'snow', 0), &
318  varattr('Tsnow_EveTr', 'degC', f146, 'Snow surface temperature for evergreen tree surface', aa, 'snow', 0), &
319  varattr('Tsnow_DecTr', 'degC', f146, 'Snow surface temperature for deciduous tree surface', aa, 'snow', 0), &
320  varattr('Tsnow_Grass', 'degC', f146, 'Snow surface temperature for grass surface', aa, 'snow', 0), &
321  varattr('Tsnow_BSoil', 'degC', f146, 'Snow surface temperature for bare soil surface', aa, 'snow', 0), &
322  varattr('Tsnow_Water', 'degC', f146, 'Snow surface temperature for water surface', aa, 'snow', 0), &
323  varattr('SnowAlb', '-', f146, 'Surface albedo for snow/ice', aa, 'snow', 0) &
324  /
325 
326  ! ESTM:
327  DATA(varlistall(n), &
330  varattr('QS', 'W m-2', f104, 'Total Storage', aa, 'ESTM', 0), &
331  varattr('QSAir', 'W m-2', f104, 'Storage air', aa, 'ESTM', 0), &
332  varattr('QSWall', 'W m-2', f104, 'Storage Wall', aa, 'ESTM', 0), &
333  varattr('QSRoof', 'W m-2', f104, 'Storage Roof', aa, 'ESTM', 0), &
334  varattr('QSGround', 'W m-2', f104, 'Storage Ground', aa, 'ESTM', 0), &
335  varattr('QSIBld', 'W m-2', f104, 'Storage Internal building', aa, 'ESTM', 0), &
336  varattr('TWALL1', 'degK', f104, 'Temperature in wall layer 1', aa, 'ESTM', 0), &
337  varattr('TWALL2', 'degK', f104, 'Temperature in wall layer 2', aa, 'ESTM', 0), &
338  varattr('TWALL3', 'degK', f104, 'Temperature in wall layer 3', aa, 'ESTM', 0), &
339  varattr('TWALL4', 'degK', f104, 'Temperature in wall layer 4', aa, 'ESTM', 0), &
340  varattr('TWALL5', 'degK', f104, 'Temperature in wall layer 5', aa, 'ESTM', 0), &
341  varattr('TROOF1', 'degK', f104, 'Temperature in roof layer 1', aa, 'ESTM', 0), &
342  varattr('TROOF2', 'degK', f104, 'Temperature in roof layer 2', aa, 'ESTM', 0), &
343  varattr('TROOF3', 'degK', f104, 'Temperature in roof layer 3', aa, 'ESTM', 0), &
344  varattr('TROOF4', 'degK', f104, 'Temperature in roof layer 4', aa, 'ESTM', 0), &
345  varattr('TROOF5', 'degK', f104, 'Temperature in roof layer 5', aa, 'ESTM', 0), &
346  varattr('TGROUND1', 'degK', f104, 'Temperature in ground layer 1', aa, 'ESTM', 0), &
347  varattr('TGROUND2', 'degK', f104, 'Temperature in ground layer 2', aa, 'ESTM', 0), &
348  varattr('TGROUND3', 'degK', f104, 'Temperature in ground layer 3', aa, 'ESTM', 0), &
349  varattr('TGROUND4', 'degK', f104, 'Temperature in ground layer 4', aa, 'ESTM', 0), &
350  varattr('TGROUND5', 'degK', f104, 'Temperature in ground layer 5', aa, 'ESTM', 0), &
351  varattr('TiBLD1', 'degK', f104, 'Temperature in internal building layer 1', aa, 'ESTM', 0), &
352  varattr('TiBLD2', 'degK', f104, 'Temperature in internal building layer 2', aa, 'ESTM', 0), &
353  varattr('TiBLD3', 'degK', f104, 'Temperature in internal building layer 3', aa, 'ESTM', 0), &
354  varattr('TiBLD4', 'degK', f104, 'Temperature in internal building layer 4', aa, 'ESTM', 0), &
355  varattr('TiBLD5', 'degK', f104, 'Temperature in internal building layer 5', aa, 'ESTM', 0), &
356  varattr('TaBLD', 'degK', f104, 'Indoor air temperature', aa, 'ESTM', 0) &
357  &/
358 
359  ! DailyState:
360  DATA(varlistall(n), &
365  + ncolumnsdataoutdailystate - 5)/ &
366  varattr('HDD1_h', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
367  varattr('HDD2_c', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
368  varattr('HDD3_Tmean', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
369  varattr('HDD4_T5d', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
370  varattr('P_day', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
371  varattr('DaysSR', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
372  varattr('GDD_EveTr', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
373  varattr('GDD_DecTr', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
374  varattr('GDD_Grass', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
375  varattr('SDD_EveTr', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
376  varattr('SDD_DecTr', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
377  varattr('SDD_Grass', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
378  varattr('Tmin', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
379  varattr('Tmax', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
380  varattr('DLHrs', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
381  varattr('LAI_EveTr', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
382  varattr('LAI_DecTr', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
383  varattr('LAI_Grass', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
384  varattr('DecidCap', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
385  varattr('Porosity', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
386  varattr('AlbEveTr', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
387  varattr('AlbDecTr', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
388  varattr('AlbGrass', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
389  varattr('WU_EveTr1', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
390  varattr('WU_EveTr2', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
391  varattr('WU_EveTr3', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
392  varattr('WU_DecTr1', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
393  varattr('WU_DecTr2', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
394  varattr('WU_DecTr3', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
395  varattr('WU_Grass1', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
396  varattr('WU_Grass2', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
397  varattr('WU_Grass3', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
398  varattr('deltaLAI', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
399  varattr('LAIlumps', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
400  varattr('AlbSnow', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
401  varattr('DensSnow_Paved', 'to be added', f146, 'to be added', al, 'DailyState', 0), &
402  varattr('DensSnow_Bldgs', 'to be added', f146, 'to be added', al, 'DailyState', 0), &
403  varattr('DensSnow_EveTr', 'to be added', f146, 'to be added', al, 'DailyState', 0), &
404  varattr('DensSnow_DecTr', 'to be added', f146, 'to be added', al, 'DailyState', 0), &
405  varattr('DensSnow_Grass', 'to be added', f146, 'to be added', al, 'DailyState', 0), &
406  varattr('DensSnow_BSoil', 'to be added', f146, 'to be added', al, 'DailyState', 0), &
407  varattr('DensSnow_Water', 'to be added', f146, 'to be added', al, 'DailyState', 0), &
408  varattr('a1', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
409  varattr('a2', 'to be added', f104, 'to be added', al, 'DailyState', 0), &
410  varattr('a3', 'to be added', f104, 'to be added', al, 'DailyState', 0) &
411  /
412 
413  ! RSL profiles
414  DATA(varlistall(n), &
418  + 1, &
422  + ncolumnsdataoutrsl - 5)/ &
423  varattr('z_1', 'm', f104, '0.1Zh', aa, 'RSL', 0), &
424  varattr('z_2', 'm', f104, '0.2Zh', aa, 'RSL', 0), &
425  varattr('z_3', 'm', f104, '0.3Zh', aa, 'RSL', 0), &
426  varattr('z_4', 'm', f104, '0.4Zh', aa, 'RSL', 0), &
427  varattr('z_5', 'm', f104, '0.5Zh', aa, 'RSL', 0), &
428  varattr('z_6', 'm', f104, '0.6Zh', aa, 'RSL', 0), &
429  varattr('z_7', 'm', f104, '0.7Zh', aa, 'RSL', 0), &
430  varattr('z_8', 'm', f104, '0.8Zh', aa, 'RSL', 0), &
431  varattr('z_9', 'm', f104, '0.9Zh', aa, 'RSL', 0), &
432  varattr('z_10', 'm', f104, 'Zh', aa, 'RSL', 0), &
433  varattr('z_11', 'm', f104, '1.1Zh', aa, 'RSL', 0), &
434  varattr('z_12', 'm', f104, '1.2Zh', aa, 'RSL', 0), &
435  varattr('z_13', 'm', f104, '1.3Zh', aa, 'RSL', 0), &
436  varattr('z_14', 'm', f146, '1.4Zh', aa, 'RSL', 0), &
437  varattr('z_15', 'm', f104, '1.5Zh', aa, 'RSL', 0), &
438  varattr('z_16', 'm', f104, '1.6Zh', aa, 'RSL', 0), &
439  varattr('z_17', 'm', f104, '1.7Zh', aa, 'RSL', 0), &
440  varattr('z_18', 'm', f104, '1.8Zh', aa, 'RSL', 0), &
441  varattr('z_19', 'm', f104, '1.9Zh', aa, 'RSL', 0), &
442  varattr('z_20', 'm', f104, '2.0Zh', aa, 'RSL', 0), &
443  varattr('z_21', 'm', f146, '2.1Zh', aa, 'RSL', 0), &
444  varattr('z_22', 'm', f104, '2.2Zh', aa, 'RSL', 0), &
445  varattr('z_23', 'm', f104, '2.3Zh', aa, 'RSL', 0), &
446  varattr('z_24', 'm', f104, '2.4Zh', aa, 'RSL', 0), &
447  varattr('z_25', 'm', f104, '2.5Zh', aa, 'RSL', 0), &
448  varattr('z_26', 'm', f104, '2.6Zh', aa, 'RSL', 0), &
449  varattr('z_27', 'm', f104, '2.7Zh', aa, 'RSL', 0), &
450  varattr('z_28', 'm', f104, '2.8Zh', aa, 'RSL', 0), &
451  varattr('z_29', 'm', f104, '2.9Zh', aa, 'RSL', 0), &
452  varattr('z_30', 'm', f104, '3.0Zh', aa, 'RSL', 0), &
453  varattr('U_1', 'm s-1', f104, 'U at 0.1Zh', aa, 'RSL', 0), &
454  varattr('U_2', 'm s-1', f104, 'U at 0.2Zh', aa, 'RSL', 0), &
455  varattr('U_3', 'm s-1', f104, 'U at 0.3Zh', aa, 'RSL', 0), &
456  varattr('U_4', 'm s-1', f104, 'U at 0.4Zh', aa, 'RSL', 0), &
457  varattr('U_5', 'm s-1', f104, 'U at 0.5Zh', aa, 'RSL', 0), &
458  varattr('U_6', 'm s-1', f104, 'U at 0.6Zh', aa, 'RSL', 0), &
459  varattr('U_7', 'm s-1', f104, 'U at 0.7Zh', aa, 'RSL', 0), &
460  varattr('U_8', 'm s-1', f104, 'U at 0.8Zh', aa, 'RSL', 0), &
461  varattr('U_9', 'm s-1', f104, 'U at 0.9Zh', aa, 'RSL', 0), &
462  varattr('U_10', 'm s-1', f104, 'U at Zh', aa, 'RSL', 0), &
463  varattr('U_11', 'm s-1', f104, 'U at 1.1Zh', aa, 'RSL', 0), &
464  varattr('U_12', 'm s-1', f104, 'U at 1.2Zh', aa, 'RSL', 0), &
465  varattr('U_13', 'm s-1', f104, 'U at 1.3Zh', aa, 'RSL', 0), &
466  varattr('U_14', 'm s-1', f146, 'U at 1.4Zh', aa, 'RSL', 0), &
467  varattr('U_15', 'm s-1', f104, 'U at 1.5Zh', aa, 'RSL', 0), &
468  varattr('U_16', 'm s-1', f104, 'U at 1.6Zh', aa, 'RSL', 0), &
469  varattr('U_17', 'm s-1', f104, 'U at 1.7Zh', aa, 'RSL', 0), &
470  varattr('U_18', 'm s-1', f104, 'U at 1.8Zh', aa, 'RSL', 0), &
471  varattr('U_19', 'm s-1', f104, 'U at 1.9Zh', aa, 'RSL', 0), &
472  varattr('U_20', 'm s-1', f104, 'U at 2.0Zh', aa, 'RSL', 0), &
473  varattr('U_21', 'm s-1', f146, 'U at 2.1Zh', aa, 'RSL', 0), &
474  varattr('U_22', 'm s-1', f104, 'U at 2.2Zh', aa, 'RSL', 0), &
475  varattr('U_23', 'm s-1', f104, 'U at 2.3Zh', aa, 'RSL', 0), &
476  varattr('U_24', 'm s-1', f104, 'U at 2.4Zh', aa, 'RSL', 0), &
477  varattr('U_25', 'm s-1', f104, 'U at 2.5Zh', aa, 'RSL', 0), &
478  varattr('U_26', 'm s-1', f104, 'U at 2.6Zh', aa, 'RSL', 0), &
479  varattr('U_27', 'm s-1', f104, 'U at 2.7Zh', aa, 'RSL', 0), &
480  varattr('U_28', 'm s-1', f104, 'U at 2.8Zh', aa, 'RSL', 0), &
481  varattr('U_29', 'm s-1', f104, 'U at 2.9Zh', aa, 'RSL', 0), &
482  varattr('U_30', 'm s-1', f104, 'U at 3.0Zh', aa, 'RSL', 0), &
483  varattr('T_1', 'degC', f104, 'T at 0.1Zh', aa, 'RSL', 0), &
484  varattr('T_2', 'degC', f104, 'T at 0.2Zh', aa, 'RSL', 0), &
485  varattr('T_3', 'degC', f104, 'T at 0.3Zh', aa, 'RSL', 0), &
486  varattr('T_4', 'degC', f104, 'T at 0.4Zh', aa, 'RSL', 0), &
487  varattr('T_5', 'degC', f104, 'T at 0.5Zh', aa, 'RSL', 0), &
488  varattr('T_6', 'degC', f104, 'T at 0.6Zh', aa, 'RSL', 0), &
489  varattr('T_7', 'degC', f104, 'T at 0.7Zh', aa, 'RSL', 0), &
490  varattr('T_8', 'degC', f104, 'T at 0.8Zh', aa, 'RSL', 0), &
491  varattr('T_9', 'degC', f104, 'T at 0.9Zh', aa, 'RSL', 0), &
492  varattr('T_10', 'degC', f104, 'T at Zh', aa, 'RSL', 0), &
493  varattr('T_11', 'degC', f104, 'T at 1.1Zh', aa, 'RSL', 0), &
494  varattr('T_12', 'degC', f104, 'T at 1.2Zh', aa, 'RSL', 0), &
495  varattr('T_13', 'degC', f104, 'T at 1.3Zh', aa, 'RSL', 0), &
496  varattr('T_14', 'degC', f146, 'T at 1.4Zh', aa, 'RSL', 0), &
497  varattr('T_15', 'degC', f104, 'T at 1.5Zh', aa, 'RSL', 0), &
498  varattr('T_16', 'degC', f104, 'T at 1.6Zh', aa, 'RSL', 0), &
499  varattr('T_17', 'degC', f104, 'T at 1.7Zh', aa, 'RSL', 0), &
500  varattr('T_18', 'degC', f104, 'T at 1.8Zh', aa, 'RSL', 0), &
501  varattr('T_19', 'degC', f104, 'T at 1.9Zh', aa, 'RSL', 0), &
502  varattr('T_20', 'degC', f104, 'T at 2.0Zh', aa, 'RSL', 0), &
503  varattr('T_21', 'degC', f146, 'T at 2.1Zh', aa, 'RSL', 0), &
504  varattr('T_22', 'degC', f104, 'T at 2.2Zh', aa, 'RSL', 0), &
505  varattr('T_23', 'degC', f104, 'T at 2.3Zh', aa, 'RSL', 0), &
506  varattr('T_24', 'degC', f104, 'T at 2.4Zh', aa, 'RSL', 0), &
507  varattr('T_25', 'degC', f104, 'T at 2.5Zh', aa, 'RSL', 0), &
508  varattr('T_26', 'degC', f104, 'T at 2.6Zh', aa, 'RSL', 0), &
509  varattr('T_27', 'degC', f104, 'T at 2.7Zh', aa, 'RSL', 0), &
510  varattr('T_28', 'degC', f104, 'T at 2.8Zh', aa, 'RSL', 0), &
511  varattr('T_29', 'degC', f104, 'T at 2.9Zh', aa, 'RSL', 0), &
512  varattr('T_30', 'degC', f104, 'T at 3.0Zh', aa, 'RSL', 0), &
513  varattr('q_1', 'g kg-1', f104, 'q at 0.1Zh', aa, 'RSL', 0), &
514  varattr('q_2', 'g kg-1', f104, 'q at 0.2Zh', aa, 'RSL', 0), &
515  varattr('q_3', 'g kg-1', f104, 'q at 0.3Zh', aa, 'RSL', 0), &
516  varattr('q_4', 'g kg-1', f104, 'q at 0.4Zh', aa, 'RSL', 0), &
517  varattr('q_5', 'g kg-1', f104, 'q at 0.5Zh', aa, 'RSL', 0), &
518  varattr('q_6', 'g kg-1', f104, 'q at 0.6Zh', aa, 'RSL', 0), &
519  varattr('q_7', 'g kg-1', f104, 'q at 0.7Zh', aa, 'RSL', 0), &
520  varattr('q_8', 'g kg-1', f104, 'q at 0.8Zh', aa, 'RSL', 0), &
521  varattr('q_9', 'g kg-1', f104, 'q at 0.9Zh', aa, 'RSL', 0), &
522  varattr('q_10', 'g kg-1', f104, 'q at Zh', aa, 'RSL', 0), &
523  varattr('q_11', 'g kg-1', f104, 'q at 1.1Zh', aa, 'RSL', 0), &
524  varattr('q_12', 'g kg-1', f104, 'q at 1.2Zh', aa, 'RSL', 0), &
525  varattr('q_13', 'g kg-1', f104, 'q at 1.3Zh', aa, 'RSL', 0), &
526  varattr('q_14', 'g kg-1', f146, 'q at 1.4Zh', aa, 'RSL', 0), &
527  varattr('q_15', 'g kg-1', f104, 'q at 1.5Zh', aa, 'RSL', 0), &
528  varattr('q_16', 'g kg-1', f104, 'q at 1.6Zh', aa, 'RSL', 0), &
529  varattr('q_17', 'g kg-1', f104, 'q at 1.7Zh', aa, 'RSL', 0), &
530  varattr('q_18', 'g kg-1', f104, 'q at 1.8Zh', aa, 'RSL', 0), &
531  varattr('q_19', 'g kg-1', f104, 'q at 1.9Zh', aa, 'RSL', 0), &
532  varattr('q_20', 'g kg-1', f104, 'q at 2.0Zh', aa, 'RSL', 0), &
533  varattr('q_21', 'g kg-1', f146, 'q at 2.1Zh', aa, 'RSL', 0), &
534  varattr('q_22', 'g kg-1', f104, 'q at 2.2Zh', aa, 'RSL', 0), &
535  varattr('q_23', 'g kg-1', f104, 'q at 2.3Zh', aa, 'RSL', 0), &
536  varattr('q_24', 'g kg-1', f104, 'q at 2.4Zh', aa, 'RSL', 0), &
537  varattr('q_25', 'g kg-1', f104, 'q at 2.5Zh', aa, 'RSL', 0), &
538  varattr('q_26', 'g kg-1', f104, 'q at 2.6Zh', aa, 'RSL', 0), &
539  varattr('q_27', 'g kg-1', f104, 'q at 2.7Zh', aa, 'RSL', 0), &
540  varattr('q_28', 'g kg-1', f104, 'q at 2.8Zh', aa, 'RSL', 0), &
541  varattr('q_29', 'g kg-1', f104, 'q at 2.9Zh', aa, 'RSL', 0), &
542  varattr('q_30', 'g kg-1', f104, 'q at 3.0Zh', aa, 'RSL', 0) &
543  /
544 
545 CONTAINS
546  ! main wrapper that handles both txt and nc files
547  SUBROUTINE suews_output(irMax, iv, Gridiv, iyr)
548  IMPLICIT NONE
549  INTEGER, INTENT(in) :: irMax
550 ! #ifdef nc
551 ! INTEGER, INTENT(in), OPTIONAL ::iv, Gridiv, iyr
552 ! #else
553  INTEGER, INTENT(in) ::iv, Gridiv, iyr
554 ! #endif
555 
556  INTEGER :: xx, err, outLevel, i
557  TYPE(varattr), DIMENSION(:), ALLOCATABLE::varListX
558  CHARACTER(len=10) :: grpList0(7)
559  CHARACTER(len=10), DIMENSION(:), ALLOCATABLE :: grpList
560  LOGICAL :: grpCond(7)
561 
562  ! determine outLevel
563  SELECT CASE (writeoutoption)
564  CASE (0) !all (not snow-related)
565  outlevel = 1
566  CASE (1) !all plus snow-related
567  outlevel = 2
568  CASE (2) !minimal output
569  outlevel = 0
570  END SELECT
571 
572  ! determine groups to output
573  ! TODO: needs to be smarter, automate this filtering
574  grplist0(1) = 'SUEWS'
575  grplist0(2) = 'SOLWEIG'
576  grplist0(3) = 'BL'
577  grplist0(4) = 'snow'
578  grplist0(5) = 'ESTM'
579  grplist0(6) = 'DailyState'
580  grplist0(7) = 'RSL'
581  grpcond = (/.true., &
582  solweigpoi_out == 1, &
583  cbluse >= 1, &
584  snowuse >= 1, &
585  storageheatmethod == 4 .OR. storageheatmethod == 14, &
586  .true., &
587  .true./)
588  xx = count(grpcond)
589 
590  ! PRINT*, grpList0,xx
591 
592  ALLOCATE (grplist(xx), stat=err)
593  IF (err /= 0) print *, "grpList: Allocation request denied"
594 
595  grplist = pack(grplist0, mask=grpcond)
596 
597  ! PRINT*, grpList,SIZE(grpList, dim=1)
598 
599  ! loop over all groups
600  DO i = 1, SIZE(grplist), 1
601  !PRINT*, 'i',i
602  xx = count(varlistall%group == trim(grplist(i)), dim=1)
603  ! PRINT*, 'number of variables:',xx, 'in group: ',grpList(i)
604  ! print*, 'all group names: ',varList%group
605  ALLOCATE (varlistx(5 + xx), stat=err)
606  IF (err /= 0) print *, "varListX: Allocation request denied"
607  ! datetime
608  varlistx(1:5) = varlistall(1:5)
609  ! variable
610  varlistx(6:5 + xx) = pack(varlistall, mask=(varlistall%group == trim(grplist(i))))
611 
612  IF (trim(varlistx(SIZE(varlistx))%group) /= 'DailyState') THEN
613  ! all output arrays but DailyState
614  ! all output frequency option:
615  ! as forcing:
616  IF (resolutionfilesout == tstep .OR. keeptstepfilesout == 1) THEN
617 ! #ifdef nc
618 ! IF (PRESENT(Gridiv)) THEN
619 ! #endif
620  CALL suews_output_txt_grp(iv, irmax, iyr, varlistx, gridiv, outlevel, tstep)
621 ! #ifdef nc
622 ! ELSE
623 ! CALL SUEWS_Output_nc_grp(irMax, varListX, outLevel, Tstep)
624 ! ENDIF
625 ! #endif
626 
627  ENDIF
628  ! as specified ResolutionFilesOut:
629  IF (resolutionfilesout /= tstep) THEN
630 ! #ifdef nc
631 ! IF (PRESENT(Gridiv)) THEN
632 ! #endif
633  CALL suews_output_txt_grp(iv, irmax, iyr, varlistx, gridiv, outlevel, resolutionfilesout)
634 ! #ifdef nc
635 ! ELSE
636 ! CALL SUEWS_Output_nc_grp(irMax, varListX, outLevel, ResolutionFilesOut)
637 ! ENDIF
638 ! #endif
639  ENDIF
640  ELSE
641  ! DailyState array, which does not need aggregation
642 ! #ifdef nc
643 ! IF (PRESENT(Gridiv)) THEN
644 ! #endif
645  CALL suews_output_txt_grp(iv, irmax, iyr, varlistx, gridiv, outlevel, tstep)
646 ! #ifdef nc
647 ! ELSE
648 ! CALL SUEWS_Output_nc_grp(irMax, varListX, outLevel, Tstep)
649 ! ENDIF
650 ! #endif
651  ENDIF
652 
653  IF (ALLOCATED(varlistx)) DEALLOCATE (varlistx, stat=err)
654  IF (err /= 0) print *, "varListX: Deallocation request denied"
655  ! PRINT*, 'i',i,'end'
656 
657  END DO
658  END SUBROUTINE suews_output
659 
660  ! output wrapper function for one group
661  SUBROUTINE suews_output_txt_grp(iv, irMax, iyr, varListX, Gridiv, outLevel, outFreq_s)
662  IMPLICIT NONE
663 
664  TYPE(varattr), DIMENSION(:), INTENT(in)::varListX
665  INTEGER, INTENT(in) :: iv, irMax, iyr, Gridiv, outLevel, outFreq_s
666 
667  INTEGER :: err
668 
669  INTEGER, DIMENSION(:), ALLOCATABLE ::id_seq ! id sequence as in the dataOutX/dataOutX_agg
670  REAL(KIND(1d0)), DIMENSION(:, :), ALLOCATABLE::dataOutX
671  REAL(KIND(1d0)), DIMENSION(:, :), ALLOCATABLE::dataOutX_agg
672 
673  IF (.NOT. ALLOCATED(dataoutx)) THEN
674  ALLOCATE (dataoutx(irmax, SIZE(varlistx)), stat=err)
675  IF (err /= 0) print *, "dataOutX: Allocation request denied"
676  ENDIF
677 
678  ! determine dataOutX array according to variable group
679  SELECT CASE (trim(varlistx(SIZE(varlistx))%group))
680  CASE ('SUEWS') !default
681  dataoutx = dataoutsuews(1:irmax, 1:SIZE(varlistx), gridiv)
682 
683  ! CASE ('SOLWEIG') !SOLWEIG
684  ! ! todo: inconsistent data structure
685  ! dataOutX = dataOutSOL(1:irMax, 1:SIZE(varListX), Gridiv)
686 
687  CASE ('BL') !BL
688  dataoutx = dataoutbl(1:irmax, 1:SIZE(varlistx), gridiv)
689 
690  CASE ('snow') !snow
691  dataoutx = dataoutsnow(1:irmax, 1:SIZE(varlistx), gridiv)
692 
693  CASE ('ESTM') !ESTM
694  dataoutx = dataoutestm(1:irmax, 1:SIZE(varlistx), gridiv)
695 
696  CASE ('RSL') !ESTM
697  dataoutx = dataoutrsl(1:irmax, 1:SIZE(varlistx), gridiv)
698 
699  CASE ('DailyState') !DailyState
700  ! get correct day index
701  CALL unique(int(pack(dataoutsuews(1:irmax, 2, gridiv), &
702  mask=(dataoutsuews(1:irmax, 3, gridiv) == 23 &
703  .AND. dataoutsuews(1:irmax, 4, gridiv) == (nsh - 1.)/nsh*60))), &
704  id_seq)
705 
706  IF (ALLOCATED(dataoutx)) THEN
707  DEALLOCATE (dataoutx)
708  IF (err /= 0) print *, "dataOutX: Deallocation request denied"
709  ENDIF
710 
711  IF (.NOT. ALLOCATED(dataoutx)) THEN
712  ALLOCATE (dataoutx(SIZE(id_seq), SIZE(varlistx)), stat=err)
713  IF (err /= 0) print *, "dataOutX: Allocation request denied"
714  ENDIF
715 
716  dataoutx = dataoutdailystate(id_seq, 1:SIZE(varlistx), gridiv)
717  ! print*, id_seq
718  ! print*, dataOutDailyState(id_seq,1:SIZE(varListX),Gridiv)
719  ! print*, 1/(nsh-nsh)
720  END SELECT
721 
722  ! aggregation:
723  ! aggregation is done for every group but 'DailyState'
724  IF (trim(varlistx(SIZE(varlistx))%group) /= 'DailyState') THEN
725 
726  CALL suews_output_agg(dataoutx_agg, dataoutx, varlistx, irmax, outfreq_s)
727  ELSE
728  IF (.NOT. ALLOCATED(dataoutx_agg)) THEN
729  ALLOCATE (dataoutx_agg(SIZE(dataoutx, dim=1), SIZE(varlistx)), stat=err)
730  IF (err /= 0) print *, ": Allocation request denied"
731  ENDIF
732  dataoutx_agg = dataoutx
733  ENDIF
734 
735  ! output:
736  ! initialise file when processing first metblock
737  IF (iv == 1) CALL suews_output_init(dataoutx_agg, varlistx, iyr, gridiv, outlevel)
738 
739  ! append the aggregated data to the specific txt file
740  CALL suews_write_txt(dataoutx_agg, varlistx, iyr, gridiv, outlevel)
741 
742  END SUBROUTINE suews_output_txt_grp
743 
744  ! initialise an output file with file name and headers
745  SUBROUTINE suews_output_init(dataOutX, varList, iyr, Gridiv, outLevel)
746  IMPLICIT NONE
747  REAL(KIND(1d0)), DIMENSION(:, :), INTENT(in)::dataOutX
748  TYPE(varattr), DIMENSION(:), INTENT(in)::varList
749  INTEGER, INTENT(in) :: iyr, Gridiv, outLevel
750 
751  TYPE(varattr), DIMENSION(:), ALLOCATABLE::varListSel
752  INTEGER :: xx, err, fn, i, nargs
753  CHARACTER(len=365) :: FileOutX
754  CHARACTER(len=3) :: itextX
755  CHARACTER(len=6) :: args(5)
756  CHARACTER(len=16*SIZE(varList)) :: FormatOut
757  CHARACTER(len=16) :: formatX
758  CHARACTER(len=16), DIMENSION(:), ALLOCATABLE:: headerOut
759 
760  ! select variables to output
761  xx = count((varlist%level <= outlevel), dim=1)
762  WRITE (itextx, '(i3)') xx
763  ALLOCATE (varlistsel(xx), stat=err)
764  IF (err /= 0) print *, "varListSel: Allocation request denied"
765  varlistsel = pack(varlist, mask=(varlist%level <= outlevel))
766 
767  ! generate file name
768  CALL filename_gen(dataoutx, varlist, iyr, gridiv, fileoutx)
769 
770  ! store right-aligned headers
771  ALLOCATE (headerout(xx), stat=err)
772  IF (err /= 0) print *, "headerOut: Allocation request denied"
773 
774  ! create format string:
775  DO i = 1, SIZE(varlistsel)
776  CALL parse(varlistsel(i)%fmt, 'if.,', args, nargs)
777  formatx = adjustl('(a'//trim(args(2))//',1x)')
778  ! adjust headers to right-aligned
779  WRITE (headerout(i), formatx) adjustr(trim(adjustl(varlistsel(i)%header)))
780  IF (i == 1) THEN
781  formatout = adjustl(trim(formatx))
782  ELSE
783  formatout = trim(formatout)//' '//adjustl(trim(formatx))
784  END IF
785  END DO
786  formatout = '('//trim(adjustl(formatout))//')'
787 
788  ! create file
789  fn = 9
790  OPEN (fn, file=trim(adjustl(fileoutx)), status='unknown')
791  ! PRINT*, 'FileOutX in SUEWS_Output_Init: ',FileOutX
792 
793  ! write out headers
794  WRITE (fn, formatout) headerout
795  CLOSE (fn)
796 
797  ! write out format file
798  CALL formatfile_gen(dataoutx, varlist, iyr, gridiv, outlevel)
799 
800  ! clean up
801  IF (ALLOCATED(varlistsel)) DEALLOCATE (varlistsel, stat=err)
802  IF (err /= 0) print *, "varListSel: Deallocation request denied"
803  IF (ALLOCATED(headerout)) DEALLOCATE (headerout, stat=err)
804  IF (err /= 0) print *, "headerOut: Deallocation request denied"
805 
806  END SUBROUTINE suews_output_init
807 
808  ! generate output format file
809  SUBROUTINE formatfile_gen(dataOutX, varList, iyr, Gridiv, outLevel)
810  IMPLICIT NONE
811  REAL(KIND(1d0)), DIMENSION(:, :), INTENT(in)::dataOutX
812  TYPE(varattr), DIMENSION(:), INTENT(in)::varList
813  INTEGER, INTENT(in) :: iyr, Gridiv, outLevel
814 
815  TYPE(varattr), DIMENSION(:), ALLOCATABLE::varListSel
816  INTEGER :: xx, err, fn, i
817  CHARACTER(len=365) :: FileOutX
818  CHARACTER(len=100*300) :: str_cat
819  CHARACTER(len=100) :: str_x = ''
820  CHARACTER(len=3) :: itextX
821 
822  ! get filename
823  CALL filename_gen(dataoutx, varlist, iyr, gridiv, fileoutx, 1)
824 
825  !select variables to output
826  xx = count((varlist%level <= outlevel), dim=1)
827  ALLOCATE (varlistsel(xx), stat=err)
828  IF (err /= 0) print *, "varListSel: Allocation request denied"
829  varlistsel = pack(varlist, mask=(varlist%level <= outlevel))
830 
831  ! create file
832  fn = 9
833  OPEN (fn, file=trim(adjustl(fileoutx)), status='unknown')
834 
835  ! write out format strings
836  ! column number:
837  str_cat = ''
838  DO i = 1, SIZE(varlistsel)
839  WRITE (itextx, '(i3)') i
840  IF (i == 1) THEN
841  str_cat = trim(adjustl(itextx))
842  ELSE
843  str_cat = trim(str_cat)//';'//adjustl(itextx)
844  ENDIF
845  END DO
846  WRITE (fn, '(a)') trim(str_cat)
847 
848  ! header:
849  str_cat = ''
850  DO i = 1, SIZE(varlistsel)
851  str_x = varlistsel(i)%header
852  IF (i == 1) THEN
853  str_cat = trim(adjustl(str_x))
854  ELSE
855  str_cat = trim(str_cat)//';'//adjustl(str_x)
856  ENDIF
857  END DO
858  WRITE (fn, '(a)') trim(str_cat)
859 
860  ! long name:
861  str_cat = ''
862  DO i = 1, SIZE(varlistsel)
863  str_x = varlistsel(i)%longNm
864  IF (i == 1) THEN
865  str_cat = trim(adjustl(str_x))
866  ELSE
867  str_cat = trim(str_cat)//';'//adjustl(str_x)
868  ENDIF
869  END DO
870  WRITE (fn, '(a)') trim(str_cat)
871 
872  ! unit:
873  str_cat = ''
874  DO i = 1, SIZE(varlistsel)
875  str_x = varlistsel(i)%unit
876  IF (i == 1) THEN
877  str_cat = trim(adjustl(str_x))
878  ELSE
879  str_cat = trim(str_cat)//';'//adjustl(str_x)
880  ENDIF
881  END DO
882  WRITE (fn, '(a)') trim(str_cat)
883 
884  ! format:
885  str_cat = ''
886  DO i = 1, SIZE(varlistsel)
887  str_x = varlistsel(i)%fmt
888  IF (i == 1) THEN
889  str_cat = trim(adjustl(str_x))
890  ELSE
891  str_cat = trim(str_cat)//';'//adjustl(str_x)
892  ENDIF
893  END DO
894  WRITE (fn, '(a)') trim(str_cat)
895 
896  ! aggregation method:
897  str_cat = ''
898  DO i = 1, SIZE(varlistsel)
899  str_x = varlistsel(i)%aggreg
900  IF (i == 1) THEN
901  str_cat = trim(adjustl(str_x))
902  ELSE
903  str_cat = trim(str_cat)//';'//adjustl(str_x)
904  ENDIF
905  END DO
906  WRITE (fn, '(a)') trim(str_cat)
907 
908  ! close file
909  CLOSE (fn)
910 
911  ! clean up
912  IF (ALLOCATED(varlistsel)) DEALLOCATE (varlistsel, stat=err)
913  IF (err /= 0) print *, "varListSel: Deallocation request denied"
914 
915  END SUBROUTINE formatfile_gen
916 
917  ! aggregate data to specified resolution
918  SUBROUTINE suews_output_agg(dataOut_agg, dataOutX, varList, irMax, outFreq_s)
919  IMPLICIT NONE
920  REAL(KIND(1d0)), DIMENSION(:, :), INTENT(in)::dataOutX
921  TYPE(varattr), DIMENSION(:), INTENT(in)::varList
922  INTEGER, INTENT(in) :: irMax, outFreq_s
923  REAL(KIND(1d0)), DIMENSION(:, :), ALLOCATABLE, INTENT(out)::dataOut_agg
924 
925  INTEGER :: nlinesOut, i, j, x
926  REAL(KIND(1d0))::dataOut_aggX(1:size(varlist))
927  REAL(KIND(1d0)), DIMENSION(:, :), ALLOCATABLE::dataOut_agg0
928  nlinesout = int(nsh/(60.*60/outfreq_s))
929  ! nGrid=SIZE(dataOutX, dim=3)
930 
931  ALLOCATE (dataout_agg(int(irmax/nlinesout), SIZE(varlist)))
932  ALLOCATE (dataout_agg0(nlinesout, SIZE(varlist)))
933 
934  DO i = nlinesout, irmax, nlinesout
935  x = i/nlinesout
936  dataout_agg0 = dataoutx(i - nlinesout + 1:i, :)
937  DO j = 1, SIZE(varlist), 1
938  ! aggregating different variables
939  SELECT CASE (varlist(j)%aggreg)
940  CASE (at) !time columns, aT
941  dataout_aggx(j) = dataout_agg0(nlinesout, j)
942  CASE (aa) !average, aA
943  dataout_aggx(j) = sum(dataout_agg0(:, j))/nlinesout
944  CASE (as) !sum, aS
945  dataout_aggx(j) = sum(dataout_agg0(:, j))
946  CASE (al) !last value, aL
947  dataout_aggx(j) = dataout_agg0(nlinesout, j)
948  END SELECT
949 
950  IF (diagnose == 1 .AND. i == irmax) THEN
951  ! IF ( i==irMax ) THEN
952  print *, 'raw data of ', j, ':'
953  print *, dataout_agg0(:, j)
954  print *, 'aggregated with method: ', varlist(j)%aggreg
955  print *, dataout_aggx(j)
956  print *, ''
957  END IF
958  END DO
959  dataout_agg(x, :) = dataout_aggx
960  END DO
961 
962  END SUBROUTINE suews_output_agg
963 
964  ! append output data to the specific file at the specified outLevel
965  SUBROUTINE suews_write_txt(dataOutX, varList, iyr, Gridiv, outLevel)
966  IMPLICIT NONE
967  REAL(KIND(1d0)), DIMENSION(:, :), INTENT(in)::dataOutX
968  TYPE(varattr), DIMENSION(:), INTENT(in)::varList
969  INTEGER, INTENT(in) :: iyr, Gridiv, outLevel
970 
971  REAL(KIND(1d0)), DIMENSION(:, :), ALLOCATABLE::dataOutSel
972  TYPE(varattr), DIMENSION(:), ALLOCATABLE::varListSel
973  CHARACTER(len=365) :: FileOutX
974  INTEGER :: fn, i, xx, err
975  INTEGER :: sizeVarListSel, sizedataOutX
976  CHARACTER(len=12*SIZE(varList)) :: FormatOut
977  ! LOGICAL :: initQ_file
978  formatout = ''
979 
980  IF (diagnose == 1) WRITE (*, *) 'Writting data of group: ', varlist(SIZE(varlist))%group
981 
982  !select variables to output
983  sizevarlistsel = count((varlist%level <= outlevel), dim=1)
984  ALLOCATE (varlistsel(sizevarlistsel), stat=err)
985  IF (err /= 0) print *, "varListSel: Allocation request denied"
986  varlistsel = pack(varlist, mask=(varlist%level <= outlevel))
987 
988  ! copy data accordingly
989  sizedataoutx = SIZE(dataoutx, dim=1)
990  ALLOCATE (dataoutsel(sizedataoutx, sizevarlistsel), stat=err)
991  IF (err /= 0) print *, "dataOutSel: Allocation request denied"
992  ! print*, SIZE(varList%level),PACK((/(i,i=1,SIZE(varList%level))/), varList%level <= outLevel)
993  ! print*, irMax,shape(dataOutX)
994  dataoutsel = dataoutx(:, pack((/(i, i=1, SIZE(varlist%level))/), varlist%level <= outlevel))
995 
996  ! create format string:
997  DO i = 1, sizevarlistsel
998  ! PRINT*,''
999  ! PRINT*,i
1000  ! PRINT*, LEN_TRIM(FormatOut),TRIM(FormatOut)
1001  ! PRINT*, LEN_TRIM(TRIM(FormatOut)//','),TRIM(FormatOut)//','
1002  IF (i == 1) THEN
1003  ! FormatOut=ADJUSTL(varListSel(i)%fmt)
1004  formatout = varlistsel(i)%fmt
1005  ELSE
1006 
1007  ! FormatOut=TRIM(FormatOut)//','//ADJUSTL(varListSel(i)%fmt)
1008  formatout = trim(formatout)//','//trim(varlistsel(i)%fmt)
1009  END IF
1010  ! PRINT*,''
1011  ! PRINT*,i
1012  ! PRINT*, 'FormatOut',FormatOut
1013  END DO
1014  formatout = '('//trim(adjustl(formatout))//')'
1015 
1016  ! get filename
1017  CALL filename_gen(dataoutsel, varlistsel, iyr, gridiv, fileoutx)
1018  ! PRINT*, 'FileOutX in SUEWS_Write_txt: ',FileOutX
1019 
1020  ! test if FileOutX has been initialised
1021  ! IF ( .NOT. initQ_file(FileOutX) ) THEN
1022  ! CALL SUEWS_Output_Init(dataOutSel,varListSel,Gridiv,outLevel)
1023  ! END IF
1024 
1025  ! write out data
1026  fn = 50
1027  OPEN (fn, file=trim(fileoutx), position='append')!,err=112)
1028  DO i = 1, sizedataoutx
1029  ! PRINT*, 'Writting',i
1030  ! PRINT*, 'FormatOut',FormatOut
1031  ! PRINT*, dataOutSel(i,1:sizeVarListSel)
1032  WRITE (fn, formatout) &
1033  (int(dataoutsel(i, xx)), xx=1, 4), &
1034  (dataoutsel(i, xx), xx=5, sizevarlistsel)
1035  ENDDO
1036  CLOSE (fn)
1037 
1038  IF (ALLOCATED(varlistsel)) DEALLOCATE (varlistsel, stat=err)
1039  IF (err /= 0) print *, "varListSel: Deallocation request denied"
1040 
1041  IF (ALLOCATED(dataoutsel)) DEALLOCATE (dataoutsel, stat=err)
1042  IF (err /= 0) print *, "dataOutSel: Deallocation request denied"
1043 
1044  END SUBROUTINE suews_write_txt
1045 
1046  SUBROUTINE filename_gen(dataOutX, varList, iyr, Gridiv, FileOutX, opt_fmt)
1048 
1049  IMPLICIT NONE
1050  REAL(KIND(1d0)), DIMENSION(:, :), INTENT(in)::dataOutX ! to determine year & output frequency
1051  TYPE(varattr), DIMENSION(:), INTENT(in)::varList ! to determine output group
1052  INTEGER, INTENT(in) :: iyr ! to determine year
1053  INTEGER, INTENT(in) :: Gridiv ! to determine grid name as in SiteSelect
1054  INTEGER, INTENT(in), OPTIONAL :: opt_fmt ! to determine if a format file
1055  CHARACTER(len=365), INTENT(out) :: FileOutX ! the output file name
1056 
1057  CHARACTER(len=20):: str_out_min, str_grid, &
1058  str_date, str_year, str_DOY, str_grp, str_sfx
1059  INTEGER :: year_int, DOY_int, val_fmt, delta_t_min
1060  TYPE(datetime) :: dt1, dt2
1061  TYPE(timedelta) :: dt_x
1062 
1063  ! initialise with a default value
1064  val_fmt = -999
1065 
1066  IF (PRESENT(opt_fmt)) val_fmt = opt_fmt
1067 
1068  ! PRINT*, varList(:)%header
1069  ! PRINT*, 'dataOutX(1)',dataOutX(1,:)
1070 
1071  ! date:
1072  doy_int = int(dataoutx(1, 2))
1073  WRITE (str_doy, '(i3.3)') doy_int
1074 
1075 ! #ifdef nc
1076 ! ! year for nc use that in dataOutX
1077 ! year_int = INT(dataOutX(1, 1))
1078 ! WRITE (str_year, '(i4)') year_int
1079 ! str_date = '_'//TRIM(ADJUSTL(str_year))
1080 ! ! add DOY as a specifier
1081 ! IF (ncMode == 1) str_date = TRIM(ADJUSTL(str_date))//TRIM(ADJUSTL(str_DOY))
1082 ! #endif
1083 
1084  ! year for txt use specified value to avoid conflicts when crossing years
1085  year_int = iyr
1086  WRITE (str_year, '(i4)') year_int
1087  str_date = '_'//trim(adjustl(str_year))
1088 
1089  ! output frequency in minute:
1090  IF (varlist(6)%group == 'DailyState') THEN
1091  str_out_min = '' ! ignore this for DailyState
1092  ELSE
1093  ! derive output frequency from output arrays
1094  ! dt_x=
1095  dt1 = datetime(int(dataoutx(1, 1)), 1, 1) + &
1096  timedelta(days=int(dataoutx(1, 2) - 1), &
1097  hours=int(dataoutx(1, 3)), &
1098  minutes=int(dataoutx(1, 4)))
1099 
1100  dt2 = datetime(int(dataoutx(2, 1)), 1, 1) + &
1101  timedelta(days=int(dataoutx(2, 2) - 1), &
1102  hours=int(dataoutx(2, 3)), &
1103  minutes=int(dataoutx(2, 4)))
1104 
1105  dt_x = dt2 - dt1
1106  delta_t_min = int(dt_x%total_seconds()/60)
1107  WRITE (str_out_min, '(i4)') delta_t_min
1108  str_out_min = '_'//trim(adjustl(str_out_min))
1109  ENDIF
1110 
1111  ! group: output type
1112  str_grp = varlist(6)%group
1113  IF (len(trim(str_grp)) > 0) str_grp = '_'//trim(adjustl(str_grp))
1114 
1115  ! grid name:
1116  WRITE (str_grid, '(i10)') grididmatrix(gridiv)
1117 ! #ifdef nc
1118 ! IF (ncMode == 1) str_grid = '' ! grid name not needed by nc files
1119 ! #endif
1120 
1121  ! suffix:
1122  str_sfx = '.txt'
1123 ! #ifdef nc
1124 ! IF (ncMode == 1) str_sfx = '.nc'
1125 ! #endif
1126 
1127  ! filename: FileOutX
1128  fileoutx = trim(fileoutputpath)// &
1129  trim(filecode)// &
1130  trim(adjustl(str_grid))// &
1131  trim(adjustl(str_date))// &
1132  trim(adjustl(str_grp))// &
1133  trim(adjustl(str_out_min))// &
1134  trim(adjustl(str_sfx))
1135 
1136  ! filename: format
1137  IF (val_fmt == 1) THEN
1138  fileoutx = trim(fileoutputpath)// &
1139  trim(filecode)// &
1140  trim(adjustl(str_grp))// &
1141  '_OutputFormat.txt'
1142  END IF
1143 
1144  END SUBROUTINE filename_gen
1145 
1146  SUBROUTINE unique(vec, vec_unique)
1147  ! Return only the unique values from vec.
1148 
1149  IMPLICIT NONE
1150 
1151  INTEGER, DIMENSION(:), INTENT(in) :: vec
1152  INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(out) :: vec_unique
1153 
1154  INTEGER :: i, num
1155  LOGICAL, DIMENSION(SIZE(vec)) :: mask
1156 
1157  mask = .false.
1158 
1159  DO i = 1, SIZE(vec)
1160 
1161  !count the number of occurrences of this element:
1162  num = count(vec(i) == vec)
1163 
1164  IF (num == 1) THEN
1165  !there is only one, flag it:
1166  mask(i) = .true.
1167  ELSE
1168  !flag this value only if it hasn't already been flagged:
1169  IF (.NOT. any(vec(i) == vec .AND. mask)) mask(i) = .true.
1170  END IF
1171 
1172  END DO
1173 
1174  !return only flagged elements:
1175  ALLOCATE (vec_unique(count(mask)))
1176  vec_unique = pack(vec, mask)
1177 
1178  !if you also need it sorted, then do so.
1179  ! For example, with slatec routine:
1180  !call ISORT (vec_unique, [0], size(vec_unique), 1)
1181 
1182  END SUBROUTINE unique
1183 
1184  ! test if a txt file has been initialised
1185  LOGICAL FUNCTION initq_file(FileName)
1186  IMPLICIT NONE
1187  CHARACTER(len=365), INTENT(in) :: FileName ! the output file name
1188  LOGICAL :: existQ
1189  CHARACTER(len=1000) :: longstring
1190 
1191  INQUIRE (file=trim(filename), exist=existq)
1192  IF (existq) THEN
1193  OPEN (10, file=trim(filename))
1194  READ (10, '(a)') longstring
1195  ! print*, 'longstring: ',longstring
1196  IF (verify(longstring, 'Year') == 0) initq_file = .false.
1197  CLOSE (unit=10)
1198  ELSE
1199  initq_file = .false.
1200  END IF
1201 
1202  END FUNCTION initq_file
1203 
1204  !========================================================================================
1205  FUNCTION count_lines(filename) RESULT(nlines)
1206  ! count the number of valid lines in a file
1207  ! invalid line starting with -9
1208 
1209  !========================================================================================
1210  IMPLICIT NONE
1211  CHARACTER(len=*) :: filename
1212  INTEGER :: nlines
1213  INTEGER :: io, iv
1214 
1215  OPEN (10, file=filename, iostat=io, status='old')
1216 
1217  ! if io error found, report iostat and exit
1218  IF (io /= 0) THEN
1219  print *, 'io', io, 'for', filename
1220  stop 'Cannot open file! '
1221  ENDIF
1222 
1223  nlines = 0
1224  DO
1225  READ (10, *, iostat=io) iv
1226  IF (io < 0 .OR. iv == -9) EXIT
1227 
1228  nlines = nlines + 1
1229  END DO
1230  CLOSE (10)
1231  nlines = nlines - 1 ! skip header
1232  END FUNCTION count_lines
1233 
1234  !===========================================================================!
1235  ! write the output of final SUEWS results in netCDF
1236  ! with spatial layout of QGIS convention
1237  ! the spatial matrix arranges successive rows down the page (i.e., north to south)
1238  ! and succesive columns across (i.e., west to east)
1239  ! the output file frequency is the same as metblocks in the main SUEWS loop
1240  !===========================================================================!
1241 
1242 ! #ifdef nc
1243 
1244 ! SUBROUTINE SUEWS_Output_nc_grp(irMax, varList, outLevel, outFreq_s)
1245 ! IMPLICIT NONE
1246 
1247 ! TYPE(varAttr), DIMENSION(:), INTENT(in)::varList
1248 ! INTEGER, INTENT(in) :: irMax, outLevel, outFreq_s
1249 
1250 ! REAL(KIND(1d0)), ALLOCATABLE::dataOutX(:, :, :)
1251 ! REAL(KIND(1d0)), ALLOCATABLE::dataOutX_agg(:, :, :), dataOutX_agg0(:, :)
1252 ! INTEGER :: iGrid, err, idMin, idMax
1253 ! INTEGER, DIMENSION(:), ALLOCATABLE ::id_seq
1254 
1255 ! IF (.NOT. ALLOCATED(dataOutX)) THEN
1256 ! ALLOCATE (dataOutX(irMax, SIZE(varList), NumberOfGrids), stat=err)
1257 ! IF (err /= 0) PRINT *, "dataOutX: Allocation request denied"
1258 ! ENDIF
1259 
1260 ! ! determine dataOutX array according to variable group
1261 ! SELECT CASE (TRIM(varList(SIZE(varList))%group))
1262 ! CASE ('SUEWS') !default
1263 ! dataOutX = dataOutSUEWS(1:irMax, 1:SIZE(varList), :)
1264 
1265 ! CASE ('SOLWEIG') !SOLWEIG
1266 ! ! todo: inconsistent data structure
1267 ! dataOutX = dataOutSOL(1:irMax, 1:SIZE(varList), :)
1268 
1269 ! CASE ('BL') !BL
1270 ! dataOutX = dataOutBL(1:irMax, 1:SIZE(varList), :)
1271 
1272 ! CASE ('snow') !snow
1273 ! dataOutX = dataOutSnow(1:irMax, 1:SIZE(varList), :)
1274 
1275 ! CASE ('ESTM') !ESTM
1276 ! dataOutX = dataOutESTM(1:irMax, 1:SIZE(varList), :)
1277 
1278 ! CASE ('DailyState') !DailyState
1279 ! ! get correct day index
1280 ! CALL unique(INT(PACK(dataOutSUEWS(1:irMax, 2, 1), &
1281 ! mask=(dataOutSUEWS(1:irMax, 3, Gridiv) == 23 &
1282 ! .AND. dataOutSUEWS(1:irMax, 4, Gridiv) == (nsh - 1)/nsh*60))), &
1283 ! id_seq)
1284 ! IF (ALLOCATED(dataOutX)) THEN
1285 ! DEALLOCATE (dataOutX)
1286 ! IF (err /= 0) PRINT *, "dataOutX: Deallocation request denied"
1287 ! ENDIF
1288 
1289 ! IF (.NOT. ALLOCATED(dataOutX)) THEN
1290 ! ALLOCATE (dataOutX(SIZE(id_seq), SIZE(varList), NumberOfGrids), stat=err)
1291 ! IF (err /= 0) PRINT *, "dataOutX: Allocation request denied"
1292 ! ENDIF
1293 
1294 ! dataOutX = dataOutDailyState(id_seq, 1:SIZE(varList), :)
1295 ! ! print*, 'idMin line',dataOutX(idMin,1:4,1)
1296 ! ! print*, 'idMax line',dataOutX(idMax,1:4,1)
1297 
1298 ! END SELECT
1299 
1300 ! ! aggregation:
1301 ! IF (TRIM(varList(SIZE(varList))%group) /= 'DailyState') THEN
1302 ! DO iGrid = 1, NumberOfGrids
1303 ! CALL SUEWS_Output_Agg(dataOutX_agg0, dataOutX(:, :, iGrid), varList, irMax, outFreq_s)
1304 ! IF (.NOT. ALLOCATED(dataOutX_agg)) THEN
1305 ! ALLOCATE (dataOutX_agg(SIZE(dataOutX_agg0, dim=1), SIZE(varList), NumberOfGrids), stat=err)
1306 ! IF (err /= 0) PRINT *, ": Allocation request denied"
1307 ! ENDIF
1308 ! dataOutX_agg(:, :, iGrid) = dataOutX_agg0
1309 ! END DO
1310 ! ELSE
1311 ! IF (.NOT. ALLOCATED(dataOutX_agg)) THEN
1312 ! ALLOCATE (dataOutX_agg(SIZE(dataOutX, dim=1), SIZE(varList), NumberOfGrids), stat=err)
1313 ! IF (err /= 0) PRINT *, ": Allocation request denied"
1314 ! ENDIF
1315 ! dataOutX_agg = dataOutX
1316 ! ENDIF
1317 
1318 ! ! write out data
1319 ! CALL SUEWS_Write_nc(dataOutX_agg, varList, outLevel)
1320 ! IF (ALLOCATED(dataOutX_agg)) THEN
1321 ! DEALLOCATE (dataOutX_agg)
1322 ! IF (err /= 0) PRINT *, "dataOutX_agg: Deallocation request denied"
1323 ! ENDIF
1324 ! END SUBROUTINE SUEWS_Output_nc_grp
1325 
1326 ! ! SUBROUTINE SUEWS_Write_nc(dataOutX, varList, outLevel)
1327 ! ! ! generic subroutine to write out data in netCDF format
1328 ! ! USE netCDF
1329 
1330 ! ! IMPLICIT NONE
1331 ! ! REAL(KIND(1d0)), DIMENSION(:, :, :), INTENT(in)::dataOutX
1332 ! ! TYPE(varAttr), DIMENSION(:), INTENT(in)::varList
1333 ! ! INTEGER, INTENT(in) :: outLevel
1334 
1335 ! ! CHARACTER(len=365):: fileOut
1336 ! ! REAL(KIND(1d0)), DIMENSION(:, :, :), ALLOCATABLE::dataOutSel
1337 ! ! TYPE(varAttr), DIMENSION(:), ALLOCATABLE::varListSel
1338 
1339 ! ! ! We are writing 3D data, {time, y, x}
1340 ! ! INTEGER, PARAMETER :: NDIMS = 3, iVarStart = 6
1341 ! ! INTEGER :: NX, NY, nTime, nVar, err
1342 
1343 ! ! ! When we create netCDF files, variables and dimensions, we get back
1344 ! ! ! an ID for each one.
1345 ! ! INTEGER :: ncID, varID, dimids(NDIMS), varIDGrid
1346 ! ! INTEGER :: x_dimid, y_dimid, time_dimid, iVar, varIDx, varIDy, varIDt, varIDCRS
1347 ! ! REAL(KIND(1d0)), ALLOCATABLE :: varOut(:, :, :), &
1348 ! ! varX(:, :), varY(:, :), &
1349 ! ! lat(:, :), lon(:, :), &
1350 ! ! varSeq0(:), varSeq(:), &
1351 ! ! xTime(:), xGridID(:, :)
1352 
1353 ! ! INTEGER :: idVar(iVarStart:SIZE(varList))
1354 ! ! CHARACTER(len=50):: header_str, longNm_str, unit_str
1355 ! ! CHARACTER(len=4) :: yrStr2
1356 ! ! CHARACTER(len=40) :: startStr2
1357 ! ! REAL(KIND(1d0)) :: minLat, maxLat, dLat, minLon, maxLon, dLon
1358 ! ! REAL(KIND(1d0)), DIMENSION(1:6) :: geoTrans
1359 ! ! CHARACTER(len=80) :: strGeoTrans
1360 
1361 ! ! ! determine number of times
1362 ! ! nTime = SIZE(dataOutX, dim=1)
1363 
1364 ! ! !select variables to output
1365 ! ! nVar = COUNT((varList%level <= outLevel), dim=1)
1366 ! ! ALLOCATE (varListSel(nVar), stat=err)
1367 ! ! IF (err /= 0) PRINT *, "varListSel: Allocation request denied"
1368 ! ! varListSel = PACK(varList, mask=(varList%level <= outLevel))
1369 
1370 ! ! ! copy data accordingly
1371 ! ! ALLOCATE (dataOutSel(nTime, nVar, NumberOfGrids), stat=err)
1372 ! ! IF (err /= 0) PRINT *, "dataOutSel: Allocation request denied"
1373 ! ! dataOutSel = dataOutX(:, PACK((/(i, i=1, SIZE(varList))/), varList%level <= outLevel), :)
1374 
1375 ! ! ! determine filename
1376 ! ! CALL filename_gen(dataOutSel(:, :, 1), varListSel, 1, FileOut)
1377 ! ! ! PRINT*, 'writing file:',TRIM(fileOut)
1378 
1379 ! ! ! set year string
1380 ! ! WRITE (yrStr2, '(i4)') INT(dataOutX(1, 1, 1))
1381 ! ! ! get start for later time unit creation
1382 ! ! startStr2 = TRIM(yrStr2)//'-01-01 00:00:00'
1383 
1384 ! ! ! define the dimension of spatial array/frame in the output
1385 ! ! nX = nCol
1386 ! ! nY = nRow
1387 
1388 ! ! ALLOCATE (varSeq0(nX*nY))
1389 ! ! ALLOCATE (varSeq(nX*nY))
1390 ! ! ALLOCATE (xGridID(nX, nY))
1391 ! ! ALLOCATE (lon(nX, nY))
1392 ! ! ALLOCATE (lat(nX, nY))
1393 ! ! ALLOCATE (varY(nX, nY))
1394 ! ! ALLOCATE (varX(nX, nY))
1395 ! ! ALLOCATE (xTime(nTime))
1396 
1397 ! ! ! GridID:
1398 ! ! varSeq = SurfaceChar(1:nX*nY, 1)
1399 ! ! ! CALL sortSeqReal(varSeq0,varSeq,nY,nX)
1400 ! ! xGridID = RESHAPE(varSeq, (/nX, nY/), order=(/1, 2/))
1401 ! ! ! PRINT*, 'before flipping:',lat(1:2,1)
1402 ! ! xGridID = xGridID(:, nY:1:-1)
1403 
1404 ! ! ! latitude:
1405 ! ! varSeq = SurfaceChar(1:nX*nY, 5)
1406 ! ! ! CALL sortSeqReal(varSeq0,varSeq,nY,nX)
1407 ! ! lat = RESHAPE(varSeq, (/nX, nY/), order=(/1, 2/))
1408 ! ! ! PRINT*, 'before flipping:',lat(1:2,1)
1409 ! ! lat = lat(:, nY:1:-1)
1410 ! ! ! PRINT*, 'after flipping:',lat(1:2,1)
1411 
1412 ! ! ! longitude:
1413 ! ! varSeq = SurfaceChar(1:nX*nY, 6)
1414 ! ! ! CALL sortSeqReal(varSeq0,varSeq,nY,nX)
1415 ! ! lon = RESHAPE(varSeq, (/nX, nY/), order=(/1, 2/))
1416 ! ! lon = lon(:, nY:1:-1)
1417 
1418 ! ! ! pass values to coordinate variables
1419 ! ! varY = lat
1420 ! ! varX = lon
1421 
1422 ! ! ! calculate GeoTransform array as needed by GDAL
1423 ! ! ! ref: http://www.perrygeo.com/python-affine-transforms.html
1424 ! ! ! the values below are different from the above ref,
1425 ! ! ! as the layout of SUEWS output is different from the schematic shown there
1426 ! ! ! SUEWS output is arranged northward down the page
1427 ! ! ! if data are formatted as a normal matrix
1428 ! ! minLat = lat(1, 1) ! the lower-left pixel
1429 ! ! maxLat = lat(1, NY) ! the upper-left pixel
1430 ! ! IF (nY > 1) THEN
1431 ! ! dLat = (maxLat - minLat)/(nY - 1) ! height of a pixel
1432 ! ! ELSE
1433 ! ! dLat = 1
1434 ! ! END IF
1435 
1436 ! ! ! PRINT*, 'lat:',minLat,maxLat,dLat
1437 ! ! minLon = lon(1, 1) ! the lower-left pixel
1438 ! ! maxLon = lon(NX, 1) ! the lower-right pixel
1439 ! ! IF (nY > 1) THEN
1440 ! ! dLon = (maxLon - minLon)/(nX - 1) ! width of a pixel
1441 ! ! ELSE
1442 ! ! dLon = 1
1443 ! ! END IF
1444 
1445 ! ! ! PRINT*, 'lon:',minLon,maxLon,dLon
1446 ! ! geoTrans(1) = minLon - dLon/2 ! x-coordinate of the lower-left corner of the lower-left pixel
1447 ! ! geoTrans(2) = dLon ! width of a pixel
1448 ! ! geoTrans(3) = 0. ! row rotation (typically zero)
1449 ! ! geoTrans(4) = minLat - dLat/2 ! y-coordinate of the of the lower-left corner of the lower-left pixel
1450 ! ! geoTrans(5) = 0. ! column rotation (typically zero)
1451 ! ! geoTrans(6) = dLat ! height of a pixel (typically negative, but here positive)
1452 ! ! ! write GeoTransform to strGeoTrans
1453 ! ! WRITE (strGeoTrans, '(6(f12.8,1x))') geoTrans
1454 
1455 ! ! ! Create the netCDF file. The nf90_clobber parameter tells netCDF to
1456 ! ! ! overwrite this file, if it already exists.
1457 ! ! CALL check(nf90_create(TRIM(fileOut), NF90_CLOBBER, ncID))
1458 
1459 ! ! ! put global attributes
1460 ! ! CALL check(nf90_put_att(ncID, NF90_GLOBAL, 'Conventions', 'CF1.6'))
1461 ! ! CALL check(nf90_put_att(ncID, NF90_GLOBAL, 'title', 'SUEWS output'))
1462 ! ! CALL check(nf90_put_att(ncID, NF90_GLOBAL, 'source', 'Micromet Group, University of Reading'))
1463 ! ! CALL check(nf90_put_att(ncID, NF90_GLOBAL, 'references', 'http://urban-climate.net/umep/SUEWS'))
1464 
1465 ! ! ! Define the dimensions. NetCDF will hand back an ID for each.
1466 ! ! ! nY = ncolumnsDataOutSUEWS-4
1467 ! ! ! nx = NumberOfGrids
1468 ! ! CALL check(nf90_def_dim(ncID, "time", NF90_UNLIMITED, time_dimid))
1469 ! ! CALL check(nf90_def_dim(ncID, "west_east", NX, x_dimid))
1470 ! ! CALL check(nf90_def_dim(ncID, "south_north", NY, y_dimid))
1471 ! ! ! PRINT*, 'good define dim'
1472 
1473 ! ! ! The dimids array is used to pass the IDs of the dimensions of
1474 ! ! ! the variables. Note that in fortran arrays are stored in
1475 ! ! ! column-major format.
1476 ! ! dimids = (/x_dimid, y_dimid, time_dimid/)
1477 
1478 ! ! ! write out each variable
1479 ! ! ALLOCATE (varOut(nX, nY, nTime))
1480 
1481 ! ! ! define all variables
1482 ! ! ! define time variable:
1483 ! ! CALL check(nf90_def_var(ncID, 'time', NF90_REAL, time_dimid, varIDt))
1484 ! ! CALL check(nf90_put_att(ncID, varIDt, 'units', 'minutes since '//startStr2))
1485 ! ! CALL check(nf90_put_att(ncID, varIDt, 'long_name', 'time'))
1486 ! ! CALL check(nf90_put_att(ncID, varIDt, 'standard_name', 'time'))
1487 ! ! CALL check(nf90_put_att(ncID, varIDt, 'calendar', 'gregorian'))
1488 ! ! CALL check(nf90_put_att(ncID, varIDt, 'axis', 'T'))
1489 
1490 ! ! ! define coordinate variables:
1491 ! ! CALL check(nf90_def_var(ncID, 'lon', NF90_REAL, (/x_dimid, y_dimid/), varIDx))
1492 ! ! CALL check(nf90_put_att(ncID, varIDx, 'units', 'degree_east'))
1493 ! ! CALL check(nf90_put_att(ncID, varIDx, 'long_name', 'longitude'))
1494 ! ! CALL check(nf90_put_att(ncID, varIDx, 'standard_name', 'longitude'))
1495 ! ! CALL check(nf90_put_att(ncID, varIDx, 'axis', 'X'))
1496 
1497 ! ! CALL check(nf90_def_var(ncID, 'lat', NF90_REAL, (/x_dimid, y_dimid/), varIDy))
1498 ! ! CALL check(nf90_put_att(ncID, varIDy, 'units', 'degree_north'))
1499 ! ! CALL check(nf90_put_att(ncID, varIDy, 'long_name', 'latitude'))
1500 ! ! CALL check(nf90_put_att(ncID, varIDy, 'standard_name', 'latitude'))
1501 ! ! CALL check(nf90_put_att(ncID, varIDy, 'axis', 'Y'))
1502 
1503 ! ! ! define coordinate referencing system:
1504 ! ! CALL check(nf90_def_var(ncID, 'crsWGS84', NF90_INT, varIDCRS))
1505 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'grid_mapping_name', 'latitude_longitude'))
1506 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'long_name', 'CRS definition'))
1507 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'longitude_of_prime_meridian', '0.0'))
1508 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'semi_major_axis', '6378137.0'))
1509 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'inverse_flattening', '298.257223563'))
1510 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'epsg_code', 'EPSG:4326'))
1511 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'GeoTransform', TRIM(strGeoTrans)))
1512 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'spatial_ref',&
1513 ! ! &'GEOGCS["WGS 84",&
1514 ! ! & DATUM["WGS_1984",&
1515 ! ! & SPHEROID["WGS 84",6378137,298.257223563,&
1516 ! ! & AUTHORITY["EPSG","7030"]],&
1517 ! ! & AUTHORITY["EPSG","6326"]],&
1518 ! ! & PRIMEM["Greenwich",0],&
1519 ! ! & UNIT["degree",0.0174532925199433],&
1520 ! ! & AUTHORITY["EPSG","4326"]]' &
1521 ! ! ))
1522 
1523 ! ! ! define grid_ID:
1524 ! ! CALL check(nf90_def_var(ncID, 'grid_ID', NF90_INT, (/x_dimid, y_dimid/), varIDGrid))
1525 ! ! CALL check(nf90_put_att(ncID, varIDGrid, 'coordinates', 'lon lat'))
1526 ! ! CALL check(nf90_put_att(ncID, varIDGrid, 'long_name', 'Grid ID as in SiteSelect'))
1527 ! ! CALL check(nf90_put_att(ncID, varIDGrid, 'grid_mapping', 'crsWGS84'))
1528 ! ! ! varIDGrid=varID
1529 
1530 ! ! ! define other 3D variables:
1531 ! ! DO iVar = iVarStart, nVar
1532 ! ! ! define variable name
1533 ! ! header_str = varListSel(iVar)%header
1534 ! ! unit_str = varListSel(iVar)%unit
1535 ! ! longNm_str = varListSel(iVar)%longNm
1536 
1537 ! ! ! Define the variable. The type of the variable in this case is
1538 ! ! ! NF90_REAL.
1539 
1540 ! ! CALL check(nf90_def_var(ncID, TRIM(ADJUSTL(header_str)), NF90_REAL, dimids, varID))
1541 
1542 ! ! CALL check(nf90_put_att(ncID, varID, 'coordinates', 'lon lat'))
1543 
1544 ! ! CALL check(nf90_put_att(ncID, varID, 'units', TRIM(ADJUSTL(unit_str))))
1545 
1546 ! ! CALL check(nf90_put_att(ncID, varID, 'long_name', TRIM(ADJUSTL(longNm_str))))
1547 
1548 ! ! CALL check(nf90_put_att(ncID, varID, 'grid_mapping', 'crsWGS84'))
1549 
1550 ! ! idVar(iVar) = varID
1551 ! ! END DO
1552 ! ! CALL check(nf90_enddef(ncID))
1553 ! ! ! End define mode. This tells netCDF we are done defining metadata.
1554 
1555 ! ! ! put all variable values into netCDF datasets
1556 ! ! ! put time variable in minute:
1557 ! ! xTime = (dataOutSel(1:nTime, 2, 1) - 1)*24*60 + dataOutSel(1:nTime, 3, 1)*60 + dataOutSel(1:nTime, 4, 1)
1558 ! ! CALL check(nf90_put_var(ncID, varIDt, xTime))
1559 
1560 ! ! ! put coordinate variables:
1561 ! ! CALL check(nf90_put_var(ncID, varIDx, varX))
1562 ! ! CALL check(nf90_put_var(ncID, varIDy, varY))
1563 
1564 ! ! ! put CRS variable:
1565 ! ! CALL check(nf90_put_var(ncID, varIDCRS, 9999))
1566 
1567 ! ! CALL check(NF90_SYNC(ncID))
1568 ! ! ! PRINT*, 'good put var'
1569 
1570 ! ! ! put grid_ID:
1571 ! ! CALL check(nf90_put_var(ncID, varIDGrid, xGridID))
1572 ! ! ! PRINT*, 'good put varIDGrid',varIDGrid
1573 
1574 ! ! CALL check(NF90_SYNC(ncID))
1575 
1576 ! ! ! then other 3D variables
1577 ! ! DO iVar = iVarStart, nVar
1578 ! ! ! reshape dataOutX to be aligned in checker board form
1579 ! ! varOut = RESHAPE(dataOutSel(1:nTime, iVar, :), (/nX, nY, nTime/), order=(/3, 1, 2/))
1580 ! ! varOut = varOut(:, nY:1:-1, :)
1581 ! ! ! get the variable id
1582 ! ! varID = idVar(iVar)
1583 
1584 ! ! CALL check(nf90_put_var(ncID, varID, varOut))
1585 
1586 ! ! CALL check(NF90_SYNC(ncID))
1587 ! ! END DO
1588 
1589 ! ! IF (ALLOCATED(varOut)) DEALLOCATE (varOut)
1590 ! ! IF (ALLOCATED(varSeq0)) DEALLOCATE (varSeq0)
1591 ! ! IF (ALLOCATED(varSeq)) DEALLOCATE (varSeq)
1592 ! ! IF (ALLOCATED(xGridID)) DEALLOCATE (xGridID)
1593 ! ! IF (ALLOCATED(lon)) DEALLOCATE (lon)
1594 ! ! IF (ALLOCATED(lat)) DEALLOCATE (lat)
1595 ! ! IF (ALLOCATED(varY)) DEALLOCATE (varY)
1596 ! ! IF (ALLOCATED(varX)) DEALLOCATE (varX)
1597 ! ! IF (ALLOCATED(xTime)) DEALLOCATE (xTime)
1598 
1599 ! ! ! Close the file. This frees up any internal netCDF resources
1600 ! ! ! associated with the file, and flushes any buffers.
1601 ! ! CALL check(nf90_close(ncID))
1602 
1603 ! ! ! PRINT*, "*** SUCCESS writing netCDF file:"
1604 ! ! ! PRINT*, FileOut
1605 ! ! END SUBROUTINE SUEWS_Write_nc
1606 
1607 ! !===========================================================================!
1608 ! ! convert a vector of grids to a matrix
1609 ! ! the grid IDs in seqGrid2Sort follow the QGIS convention
1610 ! ! the spatial matrix arranges successive rows down the page (i.e., north to south)
1611 ! ! and succesive columns across (i.e., west to east)
1612 ! ! seqGridSorted stores the grid IDs as aligned in matGrid but squeezed into a vector
1613 ! !===========================================================================!
1614 ! SUBROUTINE grid2mat(seqGrid2Sort, seqGridSorted, matGrid, nRow, nCol)
1615 
1616 ! IMPLICIT NONE
1617 
1618 ! INTEGER, DIMENSION(nRow*nCol) :: seqGrid2Sort, seqGridSorted
1619 ! INTEGER, DIMENSION(nRow, nCol) :: matGrid
1620 ! INTEGER :: nRow, nCol, i, j, loc
1621 
1622 ! CALL sortGrid(seqGrid2Sort, seqGridSorted, nRow, nCol)
1623 ! PRINT *, 'old:'
1624 ! PRINT *, seqGrid2Sort(1:5)
1625 ! PRINT *, 'sorted:'
1626 ! PRINT *, seqGridSorted(1:5)
1627 ! PRINT *, ''
1628 ! DO i = 1, nRow
1629 ! DO j = 1, nCol
1630 ! loc = (i - 1)*nCol + j
1631 ! matGrid(i, j) = seqGridSorted(loc)
1632 ! END DO
1633 ! END DO
1634 ! END SUBROUTINE grid2mat
1635 
1636 ! !===========================================================================!
1637 ! ! convert sequence of REAL values to a matrix
1638 ! ! the grid IDs in seqGrid2Sort follow the QGIS convention
1639 ! ! the spatial matrix arranges successive rows down the page (i.e., north to south)
1640 ! ! and succesive columns across (i.e., west to east)
1641 ! ! seqGridSorted stores the grid IDs as aligned in matGrid but squeezed into a vector
1642 ! !===========================================================================!
1643 ! SUBROUTINE seq2mat(seq2Sort, seqSorted, matGrid, nRow, nCol)
1644 
1645 ! IMPLICIT NONE
1646 
1647 ! REAL(KIND(1d0)), DIMENSION(nRow*nCol) :: seq2Sort, seqSorted
1648 ! REAL(KIND(1d0)), DIMENSION(nRow, nCol) :: matGrid
1649 ! INTEGER :: nRow, nCol, i, j, loc
1650 
1651 ! CALL sortSeqReal(seq2Sort, seqSorted, nRow, nCol)
1652 ! PRINT *, 'old:'
1653 ! PRINT *, seq2Sort(1:5)
1654 ! PRINT *, 'sorted:'
1655 ! PRINT *, seqSorted(1:5)
1656 ! PRINT *, ''
1657 ! DO i = 1, nRow
1658 ! DO j = 1, nCol
1659 ! loc = (i - 1)*nCol + j
1660 ! matGrid(i, j) = seqSorted(loc)
1661 ! END DO
1662 ! END DO
1663 ! END SUBROUTINE seq2mat
1664 
1665 ! !===========================================================================!
1666 ! ! sort a sequence of LONG values into the specially aligned sequence per QGIS
1667 ! !===========================================================================!
1668 ! SUBROUTINE sortGrid(seqGrid2Sort0, seqGridSorted, nRow, nCol)
1669 ! USE qsort_c_module
1670 ! ! convert a vector of grids to a matrix
1671 ! ! the grid IDs in seqGrid2Sort follow the QGIS convention
1672 ! ! the spatial matrix arranges successive rows down the page (i.e., north to south)
1673 ! ! and succesive columns across (i.e., west to east)
1674 ! ! seqGridSorted stores the grid IDs as aligned in matGrid but squeezed into a vector
1675 
1676 ! IMPLICIT NONE
1677 ! INTEGER :: nRow, nCol, i = 1, j = 1, xInd, len
1678 
1679 ! INTEGER, DIMENSION(nRow*nCol), INTENT(in) :: seqGrid2Sort0
1680 ! INTEGER, DIMENSION(nRow*nCol), INTENT(out) :: seqGridSorted
1681 ! INTEGER, DIMENSION(nRow*nCol) :: seqGrid2Sort, locSorted
1682 ! INTEGER :: loc
1683 ! REAL:: ind(nRow*nCol, 2)
1684 ! REAL, DIMENSION(nRow*nCol) :: seqGrid2SortReal, seqGridSortedReal
1685 ! REAL :: val
1686 
1687 ! ! number of grids
1688 ! len = nRow*nCol
1689 
1690 ! !sort the input array to make sure the grid order is in QGIS convention
1691 ! ! i.e., diagonally ascending
1692 ! seqGrid2SortReal = seqGrid2Sort0*1.
1693 ! CALL QsortC(seqGrid2SortReal)
1694 ! seqGrid2Sort = INT(seqGrid2SortReal)
1695 
1696 ! ! fill in an nRow*nCol array with values to determine sequence
1697 ! xInd = 1
1698 ! DO i = 1, nRow
1699 ! DO j = 1, nCol
1700 ! ! {row, col, value for sorting, index in new sequence}
1701 ! ind(xInd, :) = (/i + j + i/(nRow + 1.), xInd*1./)
1702 ! xInd = xInd + 1
1703 ! END DO
1704 ! END DO
1705 
1706 ! ! then sorted ind(:,3) will have the same order as seqGrid2Sort
1707 ! ! sort ind(:,3)
1708 ! seqGridSortedReal = ind(:, 1)*1.
1709 ! CALL QsortC(seqGridSortedReal)
1710 ! ! print*, 'sorted real:'
1711 ! ! print*, seqGridSortedReal
1712 
1713 ! ! get index of each element of old sequence in the sorted sequence
1714 ! DO i = 1, len
1715 ! ! value in old sequence
1716 ! ! val=ind(i,3)*1.
1717 ! val = seqGridSortedReal(i)
1718 ! DO j = 1, len
1719 ! IF (val == ind(j, 1)*1.) THEN
1720 ! ! location in sorted sequence
1721 ! locSorted(i) = j
1722 ! END IF
1723 ! END DO
1724 ! END DO
1725 
1726 ! ! put elements of old sequence in the sorted order
1727 ! DO i = 1, len
1728 ! loc = locSorted(i)
1729 ! seqGridSorted(loc) = seqGrid2Sort(i)
1730 ! END DO
1731 ! seqGridSorted = seqGridSorted(len:1:-1)
1732 
1733 ! END SUBROUTINE sortGrid
1734 
1735 ! !===========================================================================!
1736 ! ! sort a sequence of REAL values into the specially aligned sequence per QGIS
1737 ! !===========================================================================!
1738 ! SUBROUTINE sortSeqReal(seqReal2Sort, seqRealSorted, nRow, nCol)
1739 ! USE qsort_c_module
1740 ! ! convert a vector of grids to a matrix
1741 ! ! the grid IDs in seqReal2Sort follow the QGIS convention
1742 ! ! the spatial matrix arranges successive rows down the page (i.e., north to south)
1743 ! ! and succesive columns across (i.e., west to east)
1744 ! ! seqRealSorted stores the grid IDs as aligned in matGrid but squeezed into a vector
1745 
1746 ! IMPLICIT NONE
1747 ! INTEGER :: nRow, nCol, i = 1, j = 1, xInd, len
1748 
1749 ! REAL(KIND(1d0)), DIMENSION(nRow*nCol), INTENT(in) :: seqReal2Sort
1750 ! REAL(KIND(1d0)), DIMENSION(nRow*nCol), INTENT(out) :: seqRealSorted
1751 ! INTEGER(KIND(1d0)), DIMENSION(nRow*nCol) :: locSorted
1752 ! INTEGER(KIND(1d0)) :: loc
1753 ! REAL:: ind(nRow*nCol, 2)
1754 ! REAL :: seqRealSortedReal(nRow*nCol), val
1755 
1756 ! ! number of grids
1757 ! len = nRow*nCol
1758 
1759 ! ! fill in an nRow*nCol array with values to determine sequence
1760 ! xInd = 1
1761 ! DO i = 1, nRow
1762 ! DO j = 1, nCol
1763 ! ! {row, col, value for sorting, index in new sequence}
1764 ! ind(xInd, :) = (/i + j + i/(nRow + 1.), xInd*1./)
1765 ! xInd = xInd + 1
1766 ! END DO
1767 ! END DO
1768 
1769 ! ! then sorted ind(:,3) will have the same order as seqReal2Sort
1770 ! ! sort ind(:,3)
1771 ! seqRealSortedReal = ind(:, 1)*1.
1772 ! CALL QsortC(seqRealSortedReal)
1773 ! ! print*, 'sorted real:'
1774 ! ! print*, seqRealSortedReal
1775 
1776 ! ! get index of each element of old sequence in the sorted sequence
1777 ! DO i = 1, len
1778 ! ! value in old sequence
1779 ! ! val=ind(i,3)*1.
1780 ! val = seqRealSortedReal(i)
1781 ! DO j = 1, len
1782 ! IF (val == ind(j, 1)*1.) THEN
1783 ! ! location in sorted sequence
1784 ! locSorted(i) = j
1785 ! END IF
1786 ! END DO
1787 ! END DO
1788 
1789 ! ! put elements of old sequence in the sorted order
1790 ! DO i = 1, len
1791 ! loc = locSorted(i)
1792 ! seqRealSorted(loc) = seqReal2Sort(i)
1793 ! END DO
1794 ! seqRealSorted = seqRealSorted(len:1:-1)
1795 
1796 ! END SUBROUTINE sortSeqReal
1797 
1798 ! !===========================================================================!
1799 ! ! a wrapper for checking netCDF status
1800 ! !===========================================================================!
1801 
1802 ! SUBROUTINE check(status)
1803 ! USE netcdf
1804 ! IMPLICIT NONE
1805 
1806 ! INTEGER, INTENT(in) :: status
1807 
1808 ! IF (status /= nf90_noerr) THEN
1809 ! PRINT *, TRIM(nf90_strerror(status))
1810 ! STOP "Stopped"
1811 ! END IF
1812 ! END SUBROUTINE check
1813 ! #endif
1814 
1815 END MODULE ctrl_output
character(len=1), parameter al
integer keeptstepfilesout
integer, parameter ncolumnsdataoutestm
real(kind(1d0)), dimension(:, :, :), allocatable dataoutestm
subroutine filename_gen(dataOutX, varList, iyr, Gridiv, FileOutX, opt_fmt)
character(len=10), parameter ft
character(len=10), parameter fd
integer, parameter ncolumnsdataoutdailystate
integer resolutionfilesout
real(kind(1d0)), dimension(:, :, :), allocatable dataoutbl
integer diagnose
subroutine suews_output(irMax, iv, Gridiv, iyr)
integer solweigpoi_out
subroutine suews_write_txt(dataOutX, varList, iyr, Gridiv, outLevel)
integer snowuse
character(len=10), parameter f106
integer, parameter ncolumnsdataoutbl
integer writeoutoption
integer, dimension(:), allocatable grididmatrix
character(len=3) itext
subroutine suews_output_init(dataOutX, varList, iyr, Gridiv, outLevel)
integer, parameter ncolumnsdataoutsol
real(kind(1d0)), dimension(:, :, :), allocatable dataoutdailystate
character(len=1), parameter as
integer, parameter ncolumnsdataoutrsl
integer, parameter ncolumnsdataoutsnow
character(len=10), parameter f104
subroutine suews_output_txt_grp(iv, irMax, iyr, varListX, Gridiv, outLevel, outFreq_s)
type(varattr), dimension(500) varlistall
character(len=10), parameter f94
character(len=20) filecode
subroutine formatfile_gen(dataOutX, varList, iyr, Gridiv, outLevel)
character(len=1), parameter at
real(kind(1d0)), dimension(:, :, :), allocatable dataoutsuews
integer function count_lines(filename)
integer cbluse
character(len=10), parameter f146
real(kind(1d0)), dimension(:, :, :), allocatable dataoutrsl
character(len=1), parameter aa
integer, parameter ncolumnsdataoutsuews
subroutine parse(str, delims, args, nargs)
subroutine unique(vec, vec_unique)
logical function initq_file(FileName)
real(kind(1d0)), dimension(:, :, :), allocatable dataoutsnow
character(len=10), parameter fy
integer storageheatmethod
character(len=150) fileoutputpath
subroutine suews_output_agg(dataOut_agg, dataOutX, varList, irMax, outFreq_s)