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', 'W m-2', f106, 'Global Irradiance', aa, 'SOLWEIG', 0), &
174  varattr('DiffuseRad', 'W m-2', f106, 'Diffuse Radiation', aa, 'SOLWEIG', 0), &
175  varattr('DirectRad', 'W m-2', f106, 'Direct Radiation', aa, 'SOLWEIG', 0), &
176  varattr('Kdown2d', 'W m-2', f106, 'Incoming shortwave radiation at POI', aa, 'SOLWEIG', 0), &
177  varattr('Kup2d', 'W m-2', f106, 'Outgoing shortwave radiation at POI', aa, 'SOLWEIG', 0), &
178  varattr('Ksouth', 'W m-2', f106, 'Shortwave radiation from south at POI', aa, 'SOLWEIG', 0), &
179  varattr('Kwest', 'W m-2', f106, 'Shortwave radiation from west at POI', aa, 'SOLWEIG', 0), &
180  varattr('Knorth', 'W m-2', f106, 'Shortwave radiation from north at POI', aa, 'SOLWEIG', 0), &
181  varattr('Keast', 'W m-2', f106, 'Shortwave radiation from east at POI', aa, 'SOLWEIG', 0), &
182  varattr('Ldown2d', 'W m-2', f106, 'Incoming longwave radiation at POI', aa, 'SOLWEIG', 0), &
183  varattr('Lup2d', 'W m-2', f106, 'Outgoing longwave radiation at POI', aa, 'SOLWEIG', 0), &
184  varattr('Lsouth', 'W m-2', f106, 'Longwave radiation from west at POI', aa, 'SOLWEIG', 0), &
185  varattr('Lwest', 'W m-2', f106, 'Longwave radiation from south at POI', aa, 'SOLWEIG', 0), &
186  varattr('Lnorth', 'W m-2', f106, 'Longwave radiation from north at POI', aa, 'SOLWEIG', 0), &
187  varattr('Least', 'W m-2', f106, 'Longwave radiation from east at POI', aa, 'SOLWEIG', 0), &
188  varattr('Tmrt', 'degC', f106, 'Mean Radiant Temperature', aa, 'SOLWEIG', 0), &
189  varattr('I0', 'W m-2', f106, 'theoretical value of maximum incoming solar radiation', aa, 'SOLWEIG', 0), &
190  varattr('CI', '', f106, 'clearness index for Ldown', aa, 'SOLWEIG', 0), &
191  varattr('gvf', '', f106, 'Ground view factor', aa, 'SOLWEIG', 0), &
192  varattr('shadow', '', f106, 'Shadow value (0= shadow, 1 = sun)', aa, 'SOLWEIG', 0), &
193  varattr('svf', '', f106, 'Sky View Factor from ground and buildings', aa, 'SOLWEIG', 0), &
194  varattr('svfbuveg', '', f106, 'Sky View Factor from ground, buildings and vegetation', aa, 'SOLWEIG', 0), &
195  varattr('Ta', 'degC', f106, 'Air temperature', aa, 'SOLWEIG', 0), &
196  varattr('Tg', 'degC', f106, 'Surface temperature', 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 = [ &
582  .true., &
583  .true., &
584  cbluse >= 1, &
585  snowuse >= 1, &
586  storageheatmethod == 4 .OR. storageheatmethod == 14, &
587  .true., &
588  .true.]
589  xx = count(grpcond)
590 
591  ! PRINT*, grpList0,xx
592 
593  ALLOCATE (grplist(xx), stat=err)
594  IF (err /= 0) print *, "grpList: Allocation request denied"
595 
596  grplist = pack(grplist0, mask=grpcond)
597 
598  ! PRINT*, grpList,SIZE(grpList, dim=1)
599 
600  ! loop over all groups
601  DO i = 1, SIZE(grplist), 1
602  !PRINT*, 'i',i
603  xx = count(varlistall%group == trim(grplist(i)), dim=1)
604  ! PRINT*, 'number of variables:',xx, 'in group: ',grpList(i)
605  ! print*, 'all group names: ',varList%group
606  ALLOCATE (varlistx(5 + xx), stat=err)
607  IF (err /= 0) print *, "varListX: Allocation request denied"
608  ! datetime
609  varlistx(1:5) = varlistall(1:5)
610  ! variable
611  varlistx(6:5 + xx) = pack(varlistall, mask=(varlistall%group == trim(grplist(i))))
612 
613  IF (trim(varlistx(SIZE(varlistx))%group) /= 'DailyState') THEN
614  ! all output arrays but DailyState
615  ! all output frequency option:
616  ! as forcing:
617  IF (resolutionfilesout == tstep .OR. keeptstepfilesout == 1) THEN
618  CALL suews_output_txt_grp(iv, irmax, iyr, varlistx, gridiv, outlevel, tstep)
619  ENDIF
620  ! as specified ResolutionFilesOut:
621  IF (resolutionfilesout /= tstep) THEN
622  CALL suews_output_txt_grp(iv, irmax, iyr, varlistx, gridiv, outlevel, resolutionfilesout)
623  ENDIF
624  ELSE
625  ! DailyState array, which does not need aggregation
626 
627  CALL suews_output_txt_grp(iv, irmax, iyr, varlistx, gridiv, outlevel, tstep)
628 
629  ENDIF
630 
631  IF (ALLOCATED(varlistx)) DEALLOCATE (varlistx, stat=err)
632  IF (err /= 0) print *, "varListX: Deallocation request denied"
633  ! PRINT*, 'i',i,'end'
634 
635  END DO
636  END SUBROUTINE suews_output
637 
638  ! output wrapper function for one group
639  SUBROUTINE suews_output_txt_grp(iv, irMax, iyr, varListX, Gridiv, outLevel, outFreq_s)
640  IMPLICIT NONE
641 
642  TYPE(varattr), DIMENSION(:), INTENT(in)::varListX
643  INTEGER, INTENT(in) :: iv, irMax, iyr, Gridiv, outLevel, outFreq_s
644 
645  INTEGER :: err
646 
647  INTEGER, DIMENSION(:), ALLOCATABLE ::id_seq ! id sequence as in the dataOutX/dataOutX_agg
648  REAL(KIND(1d0)), DIMENSION(:, :), ALLOCATABLE::dataOutX
649  REAL(KIND(1d0)), DIMENSION(:, :), ALLOCATABLE::dataOutX_agg
650 
651  IF (.NOT. ALLOCATED(dataoutx)) THEN
652  ALLOCATE (dataoutx(irmax, SIZE(varlistx)), stat=err)
653  IF (err /= 0) print *, "dataOutX: Allocation request denied"
654  ENDIF
655 
656  ! determine dataOutX array according to variable group
657  SELECT CASE (trim(varlistx(SIZE(varlistx))%group))
658  CASE ('SUEWS') !default
659  dataoutx = dataoutsuews(1:irmax, 1:SIZE(varlistx), gridiv)
660 
661  CASE ('SOLWEIG') !SOLWEIG
662  dataoutx = dataoutsolweig(1:irmax, 1:SIZE(varlistx), gridiv)
663 
664  CASE ('BL') !BL
665  dataoutx = dataoutbl(1:irmax, 1:SIZE(varlistx), gridiv)
666 
667  CASE ('snow') !snow
668  dataoutx = dataoutsnow(1:irmax, 1:SIZE(varlistx), gridiv)
669 
670  CASE ('ESTM') !ESTM
671  dataoutx = dataoutestm(1:irmax, 1:SIZE(varlistx), gridiv)
672 
673  CASE ('RSL') !ESTM
674  dataoutx = dataoutrsl(1:irmax, 1:SIZE(varlistx), gridiv)
675 
676  CASE ('DailyState') !DailyState
677  ! get correct day index
678  CALL unique(int(pack(dataoutsuews(1:irmax, 2, gridiv), &
679  mask=(dataoutsuews(1:irmax, 3, gridiv) == 23 &
680  .AND. dataoutsuews(1:irmax, 4, gridiv) == (nsh - 1.)/nsh*60))), &
681  id_seq)
682 
683  IF (ALLOCATED(dataoutx)) THEN
684  DEALLOCATE (dataoutx)
685  IF (err /= 0) print *, "dataOutX: Deallocation request denied"
686  ENDIF
687 
688  IF (.NOT. ALLOCATED(dataoutx)) THEN
689  ALLOCATE (dataoutx(SIZE(id_seq), SIZE(varlistx)), stat=err)
690  IF (err /= 0) print *, "dataOutX: Allocation request denied"
691  ENDIF
692 
693  dataoutx = dataoutdailystate(id_seq, 1:SIZE(varlistx), gridiv)
694  ! print*, id_seq
695  ! print*, dataOutDailyState(id_seq,1:SIZE(varListX),Gridiv)
696  ! print*, 1/(nsh-nsh)
697  END SELECT
698 
699  ! aggregation:
700  ! aggregation is done for every group but 'DailyState'
701  IF (trim(varlistx(SIZE(varlistx))%group) /= 'DailyState') THEN
702 
703  CALL suews_output_agg(dataoutx_agg, dataoutx, varlistx, irmax, outfreq_s)
704  ELSE
705  IF (.NOT. ALLOCATED(dataoutx_agg)) THEN
706  ALLOCATE (dataoutx_agg(SIZE(dataoutx, dim=1), SIZE(varlistx)), stat=err)
707  IF (err /= 0) print *, ": Allocation request denied"
708  ENDIF
709  dataoutx_agg = dataoutx
710  ENDIF
711 
712  ! output:
713  ! initialise file when processing first metblock
714  IF (iv == 1) CALL suews_output_init(dataoutx_agg, varlistx, iyr, gridiv, outlevel)
715 
716  ! append the aggregated data to the specific txt file
717  CALL suews_write_txt(dataoutx_agg, varlistx, iyr, gridiv, outlevel)
718 
719  END SUBROUTINE suews_output_txt_grp
720 
721  ! initialise an output file with file name and headers
722  SUBROUTINE suews_output_init(dataOutX, varList, iyr, Gridiv, outLevel)
723  IMPLICIT NONE
724  REAL(KIND(1d0)), DIMENSION(:, :), INTENT(in)::dataOutX
725  TYPE(varattr), DIMENSION(:), INTENT(in)::varList
726  INTEGER, INTENT(in) :: iyr, Gridiv, outLevel
727 
728  TYPE(varattr), DIMENSION(:), ALLOCATABLE::varListSel
729  INTEGER :: xx, err, fn, i, nargs
730  CHARACTER(len=365) :: FileOutX
731  CHARACTER(len=3) :: itextX
732  CHARACTER(len=6) :: args(5)
733  CHARACTER(len=16*SIZE(varList)) :: FormatOut
734  CHARACTER(len=16) :: formatX
735  CHARACTER(len=16), DIMENSION(:), ALLOCATABLE:: headerOut
736 
737  ! select variables to output
738  xx = count((varlist%level <= outlevel), dim=1)
739  WRITE (itextx, '(i3)') xx
740  ALLOCATE (varlistsel(xx), stat=err)
741  IF (err /= 0) print *, "varListSel: Allocation request denied"
742  varlistsel = pack(varlist, mask=(varlist%level <= outlevel))
743 
744  ! generate file name
745  CALL filename_gen(dataoutx, varlist, iyr, gridiv, fileoutx)
746 
747  ! store right-aligned headers
748  ALLOCATE (headerout(xx), stat=err)
749  IF (err /= 0) print *, "headerOut: Allocation request denied"
750 
751  ! create format string:
752  DO i = 1, SIZE(varlistsel)
753  CALL parse(varlistsel(i)%fmt, 'if.,', args, nargs)
754  formatx = adjustl('(a'//trim(args(2))//',1x)')
755  ! adjust headers to right-aligned
756  WRITE (headerout(i), formatx) adjustr(trim(adjustl(varlistsel(i)%header)))
757  IF (i == 1) THEN
758  formatout = adjustl(trim(formatx))
759  ELSE
760  formatout = trim(formatout)//' '//adjustl(trim(formatx))
761  END IF
762  END DO
763  formatout = '('//trim(adjustl(formatout))//')'
764 
765  ! create file
766  fn = 9
767  OPEN (fn, file=trim(adjustl(fileoutx)), status='unknown')
768  ! PRINT*, 'FileOutX in SUEWS_Output_Init: ',FileOutX
769 
770  ! write out headers
771  WRITE (fn, formatout) headerout
772  CLOSE (fn)
773 
774  ! write out format file
775  CALL formatfile_gen(dataoutx, varlist, iyr, gridiv, outlevel)
776 
777  ! clean up
778  IF (ALLOCATED(varlistsel)) DEALLOCATE (varlistsel, stat=err)
779  IF (err /= 0) print *, "varListSel: Deallocation request denied"
780  IF (ALLOCATED(headerout)) DEALLOCATE (headerout, stat=err)
781  IF (err /= 0) print *, "headerOut: Deallocation request denied"
782 
783  END SUBROUTINE suews_output_init
784 
785  ! generate output format file
786  SUBROUTINE formatfile_gen(dataOutX, varList, iyr, Gridiv, outLevel)
787  IMPLICIT NONE
788  REAL(KIND(1d0)), DIMENSION(:, :), INTENT(in)::dataOutX
789  TYPE(varattr), DIMENSION(:), INTENT(in)::varList
790  INTEGER, INTENT(in) :: iyr, Gridiv, outLevel
791 
792  TYPE(varattr), DIMENSION(:), ALLOCATABLE::varListSel
793  INTEGER :: xx, err, fn, i
794  CHARACTER(len=365) :: FileOutX
795  CHARACTER(len=100*300) :: str_cat
796  CHARACTER(len=100) :: str_x = ''
797  CHARACTER(len=3) :: itextX
798 
799  ! get filename
800  CALL filename_gen(dataoutx, varlist, iyr, gridiv, fileoutx, 1)
801 
802  !select variables to output
803  xx = count((varlist%level <= outlevel), dim=1)
804  ALLOCATE (varlistsel(xx), stat=err)
805  IF (err /= 0) print *, "varListSel: Allocation request denied"
806  varlistsel = pack(varlist, mask=(varlist%level <= outlevel))
807 
808  ! create file
809  fn = 9
810  OPEN (fn, file=trim(adjustl(fileoutx)), status='unknown')
811 
812  ! write out format strings
813  ! column number:
814  str_cat = ''
815  DO i = 1, SIZE(varlistsel)
816  WRITE (itextx, '(i3)') i
817  IF (i == 1) THEN
818  str_cat = trim(adjustl(itextx))
819  ELSE
820  str_cat = trim(str_cat)//';'//adjustl(itextx)
821  ENDIF
822  END DO
823  WRITE (fn, '(a)') trim(str_cat)
824 
825  ! header:
826  str_cat = ''
827  DO i = 1, SIZE(varlistsel)
828  str_x = varlistsel(i)%header
829  IF (i == 1) THEN
830  str_cat = trim(adjustl(str_x))
831  ELSE
832  str_cat = trim(str_cat)//';'//adjustl(str_x)
833  ENDIF
834  END DO
835  WRITE (fn, '(a)') trim(str_cat)
836 
837  ! long name:
838  str_cat = ''
839  DO i = 1, SIZE(varlistsel)
840  str_x = varlistsel(i)%longNm
841  IF (i == 1) THEN
842  str_cat = trim(adjustl(str_x))
843  ELSE
844  str_cat = trim(str_cat)//';'//adjustl(str_x)
845  ENDIF
846  END DO
847  WRITE (fn, '(a)') trim(str_cat)
848 
849  ! unit:
850  str_cat = ''
851  DO i = 1, SIZE(varlistsel)
852  str_x = varlistsel(i)%unit
853  IF (i == 1) THEN
854  str_cat = trim(adjustl(str_x))
855  ELSE
856  str_cat = trim(str_cat)//';'//adjustl(str_x)
857  ENDIF
858  END DO
859  WRITE (fn, '(a)') trim(str_cat)
860 
861  ! format:
862  str_cat = ''
863  DO i = 1, SIZE(varlistsel)
864  str_x = varlistsel(i)%fmt
865  IF (i == 1) THEN
866  str_cat = trim(adjustl(str_x))
867  ELSE
868  str_cat = trim(str_cat)//';'//adjustl(str_x)
869  ENDIF
870  END DO
871  WRITE (fn, '(a)') trim(str_cat)
872 
873  ! aggregation method:
874  str_cat = ''
875  DO i = 1, SIZE(varlistsel)
876  str_x = varlistsel(i)%aggreg
877  IF (i == 1) THEN
878  str_cat = trim(adjustl(str_x))
879  ELSE
880  str_cat = trim(str_cat)//';'//adjustl(str_x)
881  ENDIF
882  END DO
883  WRITE (fn, '(a)') trim(str_cat)
884 
885  ! close file
886  CLOSE (fn)
887 
888  ! clean up
889  IF (ALLOCATED(varlistsel)) DEALLOCATE (varlistsel, stat=err)
890  IF (err /= 0) print *, "varListSel: Deallocation request denied"
891 
892  END SUBROUTINE formatfile_gen
893 
894  ! aggregate data to specified resolution
895  SUBROUTINE suews_output_agg(dataOut_agg, dataOutX, varList, irMax, outFreq_s)
896  IMPLICIT NONE
897  REAL(KIND(1d0)), DIMENSION(:, :), INTENT(in)::dataOutX
898  TYPE(varattr), DIMENSION(:), INTENT(in)::varList
899  INTEGER, INTENT(in) :: irMax, outFreq_s
900  REAL(KIND(1d0)), DIMENSION(:, :), ALLOCATABLE, INTENT(out)::dataOut_agg
901 
902  INTEGER :: nlinesOut, i, j, x
903  REAL(KIND(1d0))::dataOut_aggX(1:size(varlist))
904  REAL(KIND(1d0)), DIMENSION(:, :), ALLOCATABLE::dataOut_agg0
905  nlinesout = int(nsh/(60.*60/outfreq_s))
906  ! nGrid=SIZE(dataOutX, dim=3)
907 
908  ALLOCATE (dataout_agg(int(irmax/nlinesout), SIZE(varlist)))
909  ALLOCATE (dataout_agg0(nlinesout, SIZE(varlist)))
910 
911  DO i = nlinesout, irmax, nlinesout
912  x = i/nlinesout
913  dataout_agg0 = dataoutx(i - nlinesout + 1:i, :)
914  DO j = 1, SIZE(varlist), 1
915  ! aggregating different variables
916  SELECT CASE (varlist(j)%aggreg)
917  CASE (at) !time columns, aT
918  dataout_aggx(j) = dataout_agg0(nlinesout, j)
919  CASE (aa) !average, aA
920  dataout_aggx(j) = sum(dataout_agg0(:, j))/nlinesout
921  CASE (as) !sum, aS
922  dataout_aggx(j) = sum(dataout_agg0(:, j))
923  CASE (al) !last value, aL
924  dataout_aggx(j) = dataout_agg0(nlinesout, j)
925  END SELECT
926 
927  IF (diagnose == 1 .AND. i == irmax) THEN
928  ! IF ( i==irMax ) THEN
929  print *, 'raw data of ', j, ':'
930  print *, dataout_agg0(:, j)
931  print *, 'aggregated with method: ', varlist(j)%aggreg
932  print *, dataout_aggx(j)
933  print *, ''
934  END IF
935  END DO
936  dataout_agg(x, :) = dataout_aggx
937  END DO
938 
939  END SUBROUTINE suews_output_agg
940 
941  ! append output data to the specific file at the specified outLevel
942  SUBROUTINE suews_write_txt(dataOutX, varList, iyr, Gridiv, outLevel)
943  IMPLICIT NONE
944  REAL(KIND(1d0)), DIMENSION(:, :), INTENT(in)::dataOutX
945  TYPE(varattr), DIMENSION(:), INTENT(in)::varList
946  INTEGER, INTENT(in) :: iyr, Gridiv, outLevel
947 
948  REAL(KIND(1d0)), DIMENSION(:, :), ALLOCATABLE::dataOutSel
949  TYPE(varattr), DIMENSION(:), ALLOCATABLE::varListSel
950  CHARACTER(len=365) :: FileOutX
951  INTEGER :: fn, i, xx, err
952  INTEGER :: sizeVarListSel, sizedataOutX
953  CHARACTER(len=12*SIZE(varList)) :: FormatOut
954  ! LOGICAL :: initQ_file
955  formatout = ''
956 
957  IF (diagnose == 1) WRITE (*, *) 'Writting data of group: ', varlist(SIZE(varlist))%group
958 
959  !select variables to output
960  sizevarlistsel = count((varlist%level <= outlevel), dim=1)
961  ALLOCATE (varlistsel(sizevarlistsel), stat=err)
962  IF (err /= 0) print *, "varListSel: Allocation request denied"
963  varlistsel = pack(varlist, mask=(varlist%level <= outlevel))
964 
965  ! copy data accordingly
966  sizedataoutx = SIZE(dataoutx, dim=1)
967  ALLOCATE (dataoutsel(sizedataoutx, sizevarlistsel), stat=err)
968  IF (err /= 0) print *, "dataOutSel: Allocation request denied"
969  ! print*, SIZE(varList%level),PACK((/(i,i=1,SIZE(varList%level))/), varList%level <= outLevel)
970  ! print*, irMax,shape(dataOutX)
971  dataoutsel = dataoutx(:, pack((/(i, i=1, SIZE(varlist%level))/), varlist%level <= outlevel))
972 
973  ! create format string:
974  DO i = 1, sizevarlistsel
975  ! PRINT*,''
976  ! PRINT*,i
977  ! PRINT*, LEN_TRIM(FormatOut),TRIM(FormatOut)
978  ! PRINT*, LEN_TRIM(TRIM(FormatOut)//','),TRIM(FormatOut)//','
979  IF (i == 1) THEN
980  ! FormatOut=ADJUSTL(varListSel(i)%fmt)
981  formatout = varlistsel(i)%fmt
982  ELSE
983 
984  ! FormatOut=TRIM(FormatOut)//','//ADJUSTL(varListSel(i)%fmt)
985  formatout = trim(formatout)//','//trim(varlistsel(i)%fmt)
986  END IF
987  ! PRINT*,''
988  ! PRINT*,i
989  ! PRINT*, 'FormatOut',FormatOut
990  END DO
991  formatout = '('//trim(adjustl(formatout))//')'
992 
993  ! get filename
994  CALL filename_gen(dataoutsel, varlistsel, iyr, gridiv, fileoutx)
995  ! PRINT*, 'FileOutX in SUEWS_Write_txt: ',FileOutX
996 
997  ! test if FileOutX has been initialised
998  ! IF ( .NOT. initQ_file(FileOutX) ) THEN
999  ! CALL SUEWS_Output_Init(dataOutSel,varListSel,Gridiv,outLevel)
1000  ! END IF
1001 
1002  ! write out data
1003  fn = 50
1004  OPEN (fn, file=trim(fileoutx), position='append')!,err=112)
1005  DO i = 1, sizedataoutx
1006  ! PRINT*, 'Writting',i
1007  ! PRINT*, 'FormatOut',FormatOut
1008  ! PRINT*, dataOutSel(i,1:sizeVarListSel)
1009  WRITE (fn, formatout) &
1010  (int(dataoutsel(i, xx)), xx=1, 4), &
1011  (dataoutsel(i, xx), xx=5, sizevarlistsel)
1012  ENDDO
1013  CLOSE (fn)
1014 
1015  IF (ALLOCATED(varlistsel)) DEALLOCATE (varlistsel, stat=err)
1016  IF (err /= 0) print *, "varListSel: Deallocation request denied"
1017 
1018  IF (ALLOCATED(dataoutsel)) DEALLOCATE (dataoutsel, stat=err)
1019  IF (err /= 0) print *, "dataOutSel: Deallocation request denied"
1020 
1021  END SUBROUTINE suews_write_txt
1022 
1023  SUBROUTINE filename_gen(dataOutX, varList, iyr, Gridiv, FileOutX, opt_fmt)
1025 
1026  IMPLICIT NONE
1027  REAL(KIND(1d0)), DIMENSION(:, :), INTENT(in)::dataOutX ! to determine year & output frequency
1028  TYPE(varattr), DIMENSION(:), INTENT(in)::varList ! to determine output group
1029  INTEGER, INTENT(in) :: iyr ! to determine year
1030  INTEGER, INTENT(in) :: Gridiv ! to determine grid name as in SiteSelect
1031  INTEGER, INTENT(in), OPTIONAL :: opt_fmt ! to determine if a format file
1032  CHARACTER(len=365), INTENT(out) :: FileOutX ! the output file name
1033 
1034  CHARACTER(len=20):: str_out_min, str_grid, &
1035  str_date, str_year, str_DOY, str_grp, str_sfx
1036  INTEGER :: year_int, DOY_int, val_fmt, delta_t_min
1037  TYPE(datetime) :: dt1, dt2
1038  TYPE(timedelta) :: dt_x
1039 
1040  ! initialise with a default value
1041  val_fmt = -999
1042 
1043  IF (PRESENT(opt_fmt)) val_fmt = opt_fmt
1044 
1045  ! PRINT*, varList(:)%header
1046  ! PRINT*, 'dataOutX(1)',dataOutX(1,:)
1047 
1048  ! date:
1049  doy_int = int(dataoutx(1, 2))
1050  WRITE (str_doy, '(i3.3)') doy_int
1051 
1052 ! #ifdef nc
1053 ! ! year for nc use that in dataOutX
1054 ! year_int = INT(dataOutX(1, 1))
1055 ! WRITE (str_year, '(i4)') year_int
1056 ! str_date = '_'//TRIM(ADJUSTL(str_year))
1057 ! ! add DOY as a specifier
1058 ! IF (ncMode == 1) str_date = TRIM(ADJUSTL(str_date))//TRIM(ADJUSTL(str_DOY))
1059 ! #endif
1060 
1061  ! year for txt use specified value to avoid conflicts when crossing years
1062  year_int = iyr
1063  WRITE (str_year, '(i4)') year_int
1064  str_date = '_'//trim(adjustl(str_year))
1065 
1066  ! output frequency in minute:
1067  IF (varlist(6)%group == 'DailyState') THEN
1068  str_out_min = '' ! ignore this for DailyState
1069  ELSE
1070  ! derive output frequency from output arrays
1071  ! dt_x=
1072  dt1 = datetime(int(dataoutx(1, 1)), 1, 1) + &
1073  timedelta(days=int(dataoutx(1, 2) - 1), &
1074  hours=int(dataoutx(1, 3)), &
1075  minutes=int(dataoutx(1, 4)))
1076 
1077  dt2 = datetime(int(dataoutx(2, 1)), 1, 1) + &
1078  timedelta(days=int(dataoutx(2, 2) - 1), &
1079  hours=int(dataoutx(2, 3)), &
1080  minutes=int(dataoutx(2, 4)))
1081 
1082  dt_x = dt2 - dt1
1083  delta_t_min = int(dt_x%total_seconds()/60)
1084  WRITE (str_out_min, '(i4)') delta_t_min
1085  str_out_min = '_'//trim(adjustl(str_out_min))
1086  ENDIF
1087 
1088  ! group: output type
1089  str_grp = varlist(6)%group
1090  IF (len(trim(str_grp)) > 0) str_grp = '_'//trim(adjustl(str_grp))
1091 
1092  ! grid name:
1093  WRITE (str_grid, '(i10)') grididmatrix(gridiv)
1094 ! #ifdef nc
1095 ! IF (ncMode == 1) str_grid = '' ! grid name not needed by nc files
1096 ! #endif
1097 
1098  ! suffix:
1099  str_sfx = '.txt'
1100 ! #ifdef nc
1101 ! IF (ncMode == 1) str_sfx = '.nc'
1102 ! #endif
1103 
1104  ! filename: FileOutX
1105  fileoutx = trim(fileoutputpath)// &
1106  trim(filecode)// &
1107  trim(adjustl(str_grid))// &
1108  trim(adjustl(str_date))// &
1109  trim(adjustl(str_grp))// &
1110  trim(adjustl(str_out_min))// &
1111  trim(adjustl(str_sfx))
1112 
1113  ! filename: format
1114  IF (val_fmt == 1) THEN
1115  fileoutx = trim(fileoutputpath)// &
1116  trim(filecode)// &
1117  trim(adjustl(str_grp))// &
1118  '_OutputFormat.txt'
1119  END IF
1120 
1121  END SUBROUTINE filename_gen
1122 
1123  SUBROUTINE unique(vec, vec_unique)
1124  ! Return only the unique values from vec.
1125 
1126  IMPLICIT NONE
1127 
1128  INTEGER, DIMENSION(:), INTENT(in) :: vec
1129  INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(out) :: vec_unique
1130 
1131  INTEGER :: i, num
1132  LOGICAL, DIMENSION(SIZE(vec)) :: mask
1133 
1134  mask = .false.
1135 
1136  DO i = 1, SIZE(vec)
1137 
1138  !count the number of occurrences of this element:
1139  num = count(vec(i) == vec)
1140 
1141  IF (num == 1) THEN
1142  !there is only one, flag it:
1143  mask(i) = .true.
1144  ELSE
1145  !flag this value only if it hasn't already been flagged:
1146  IF (.NOT. any(vec(i) == vec .AND. mask)) mask(i) = .true.
1147  END IF
1148 
1149  END DO
1150 
1151  !return only flagged elements:
1152  ALLOCATE (vec_unique(count(mask)))
1153  vec_unique = pack(vec, mask)
1154 
1155  !if you also need it sorted, then do so.
1156  ! For example, with slatec routine:
1157  !call ISORT (vec_unique, [0], size(vec_unique), 1)
1158 
1159  END SUBROUTINE unique
1160 
1161  ! test if a txt file has been initialised
1162  LOGICAL FUNCTION initq_file(FileName)
1163  IMPLICIT NONE
1164  CHARACTER(len=365), INTENT(in) :: FileName ! the output file name
1165  LOGICAL :: existQ
1166  CHARACTER(len=1000) :: longstring
1167 
1168  INQUIRE (file=trim(filename), exist=existq)
1169  IF (existq) THEN
1170  OPEN (10, file=trim(filename))
1171  READ (10, '(a)') longstring
1172  ! print*, 'longstring: ',longstring
1173  IF (verify(longstring, 'Year') == 0) initq_file = .false.
1174  CLOSE (unit=10)
1175  ELSE
1176  initq_file = .false.
1177  END IF
1178 
1179  END FUNCTION initq_file
1180 
1181  !========================================================================================
1182  FUNCTION count_lines(filename) RESULT(nlines)
1183  ! count the number of valid lines in a file
1184  ! invalid line starting with -9
1185 
1186  !========================================================================================
1187  IMPLICIT NONE
1188  CHARACTER(len=*) :: filename
1189  INTEGER :: nlines
1190  INTEGER :: io, iv
1191 
1192  OPEN (10, file=filename, iostat=io, status='old')
1193 
1194  ! if io error found, report iostat and exit
1195  IF (io /= 0) THEN
1196  print *, 'io', io, 'for', filename
1197  stop 'Cannot open file! '
1198  ENDIF
1199 
1200  nlines = 0
1201  DO
1202  READ (10, *, iostat=io) iv
1203  IF (io < 0 .OR. iv == -9) EXIT
1204 
1205  nlines = nlines + 1
1206  END DO
1207  CLOSE (10)
1208  nlines = nlines - 1 ! skip header
1209  END FUNCTION count_lines
1210 
1211  !===========================================================================!
1212  ! write the output of final SUEWS results in netCDF
1213  ! with spatial layout of QGIS convention
1214  ! the spatial matrix arranges successive rows down the page (i.e., north to south)
1215  ! and succesive columns across (i.e., west to east)
1216  ! the output file frequency is the same as metblocks in the main SUEWS loop
1217  !===========================================================================!
1218 
1219 ! #ifdef nc
1220 
1221 ! SUBROUTINE SUEWS_Output_nc_grp(irMax, varList, outLevel, outFreq_s)
1222 ! IMPLICIT NONE
1223 
1224 ! TYPE(varAttr), DIMENSION(:), INTENT(in)::varList
1225 ! INTEGER, INTENT(in) :: irMax, outLevel, outFreq_s
1226 
1227 ! REAL(KIND(1d0)), ALLOCATABLE::dataOutX(:, :, :)
1228 ! REAL(KIND(1d0)), ALLOCATABLE::dataOutX_agg(:, :, :), dataOutX_agg0(:, :)
1229 ! INTEGER :: iGrid, err, idMin, idMax
1230 ! INTEGER, DIMENSION(:), ALLOCATABLE ::id_seq
1231 
1232 ! IF (.NOT. ALLOCATED(dataOutX)) THEN
1233 ! ALLOCATE (dataOutX(irMax, SIZE(varList), NumberOfGrids), stat=err)
1234 ! IF (err /= 0) PRINT *, "dataOutX: Allocation request denied"
1235 ! ENDIF
1236 
1237 ! ! determine dataOutX array according to variable group
1238 ! SELECT CASE (TRIM(varList(SIZE(varList))%group))
1239 ! CASE ('SUEWS') !default
1240 ! dataOutX = dataOutSUEWS(1:irMax, 1:SIZE(varList), :)
1241 
1242 ! CASE ('SOLWEIG') !SOLWEIG
1243 ! ! todo: inconsistent data structure
1244 ! dataOutX = dataOutSOLWEIG(1:irMax, 1:SIZE(varList), :)
1245 
1246 ! CASE ('BL') !BL
1247 ! dataOutX = dataOutBL(1:irMax, 1:SIZE(varList), :)
1248 
1249 ! CASE ('snow') !snow
1250 ! dataOutX = dataOutSnow(1:irMax, 1:SIZE(varList), :)
1251 
1252 ! CASE ('ESTM') !ESTM
1253 ! dataOutX = dataOutESTM(1:irMax, 1:SIZE(varList), :)
1254 
1255 ! CASE ('DailyState') !DailyState
1256 ! ! get correct day index
1257 ! CALL unique(INT(PACK(dataOutSUEWS(1:irMax, 2, 1), &
1258 ! mask=(dataOutSUEWS(1:irMax, 3, Gridiv) == 23 &
1259 ! .AND. dataOutSUEWS(1:irMax, 4, Gridiv) == (nsh - 1)/nsh*60))), &
1260 ! id_seq)
1261 ! IF (ALLOCATED(dataOutX)) THEN
1262 ! DEALLOCATE (dataOutX)
1263 ! IF (err /= 0) PRINT *, "dataOutX: Deallocation request denied"
1264 ! ENDIF
1265 
1266 ! IF (.NOT. ALLOCATED(dataOutX)) THEN
1267 ! ALLOCATE (dataOutX(SIZE(id_seq), SIZE(varList), NumberOfGrids), stat=err)
1268 ! IF (err /= 0) PRINT *, "dataOutX: Allocation request denied"
1269 ! ENDIF
1270 
1271 ! dataOutX = dataOutDailyState(id_seq, 1:SIZE(varList), :)
1272 ! ! print*, 'idMin line',dataOutX(idMin,1:4,1)
1273 ! ! print*, 'idMax line',dataOutX(idMax,1:4,1)
1274 
1275 ! END SELECT
1276 
1277 ! ! aggregation:
1278 ! IF (TRIM(varList(SIZE(varList))%group) /= 'DailyState') THEN
1279 ! DO iGrid = 1, NumberOfGrids
1280 ! CALL SUEWS_Output_Agg(dataOutX_agg0, dataOutX(:, :, iGrid), varList, irMax, outFreq_s)
1281 ! IF (.NOT. ALLOCATED(dataOutX_agg)) THEN
1282 ! ALLOCATE (dataOutX_agg(SIZE(dataOutX_agg0, dim=1), SIZE(varList), NumberOfGrids), stat=err)
1283 ! IF (err /= 0) PRINT *, ": Allocation request denied"
1284 ! ENDIF
1285 ! dataOutX_agg(:, :, iGrid) = dataOutX_agg0
1286 ! END DO
1287 ! ELSE
1288 ! IF (.NOT. ALLOCATED(dataOutX_agg)) THEN
1289 ! ALLOCATE (dataOutX_agg(SIZE(dataOutX, dim=1), SIZE(varList), NumberOfGrids), stat=err)
1290 ! IF (err /= 0) PRINT *, ": Allocation request denied"
1291 ! ENDIF
1292 ! dataOutX_agg = dataOutX
1293 ! ENDIF
1294 
1295 ! ! write out data
1296 ! CALL SUEWS_Write_nc(dataOutX_agg, varList, outLevel)
1297 ! IF (ALLOCATED(dataOutX_agg)) THEN
1298 ! DEALLOCATE (dataOutX_agg)
1299 ! IF (err /= 0) PRINT *, "dataOutX_agg: Deallocation request denied"
1300 ! ENDIF
1301 ! END SUBROUTINE SUEWS_Output_nc_grp
1302 
1303 ! ! SUBROUTINE SUEWS_Write_nc(dataOutX, varList, outLevel)
1304 ! ! ! generic subroutine to write out data in netCDF format
1305 ! ! USE netCDF
1306 
1307 ! ! IMPLICIT NONE
1308 ! ! REAL(KIND(1d0)), DIMENSION(:, :, :), INTENT(in)::dataOutX
1309 ! ! TYPE(varAttr), DIMENSION(:), INTENT(in)::varList
1310 ! ! INTEGER, INTENT(in) :: outLevel
1311 
1312 ! ! CHARACTER(len=365):: fileOut
1313 ! ! REAL(KIND(1d0)), DIMENSION(:, :, :), ALLOCATABLE::dataOutSel
1314 ! ! TYPE(varAttr), DIMENSION(:), ALLOCATABLE::varListSel
1315 
1316 ! ! ! We are writing 3D data, {time, y, x}
1317 ! ! INTEGER, PARAMETER :: NDIMS = 3, iVarStart = 6
1318 ! ! INTEGER :: NX, NY, nTime, nVar, err
1319 
1320 ! ! ! When we create netCDF files, variables and dimensions, we get back
1321 ! ! ! an ID for each one.
1322 ! ! INTEGER :: ncID, varID, dimids(NDIMS), varIDGrid
1323 ! ! INTEGER :: x_dimid, y_dimid, time_dimid, iVar, varIDx, varIDy, varIDt, varIDCRS
1324 ! ! REAL(KIND(1d0)), ALLOCATABLE :: varOut(:, :, :), &
1325 ! ! varX(:, :), varY(:, :), &
1326 ! ! lat(:, :), lon(:, :), &
1327 ! ! varSeq0(:), varSeq(:), &
1328 ! ! xTime(:), xGridID(:, :)
1329 
1330 ! ! INTEGER :: idVar(iVarStart:SIZE(varList))
1331 ! ! CHARACTER(len=50):: header_str, longNm_str, unit_str
1332 ! ! CHARACTER(len=4) :: yrStr2
1333 ! ! CHARACTER(len=40) :: startStr2
1334 ! ! REAL(KIND(1d0)) :: minLat, maxLat, dLat, minLon, maxLon, dLon
1335 ! ! REAL(KIND(1d0)), DIMENSION(1:6) :: geoTrans
1336 ! ! CHARACTER(len=80) :: strGeoTrans
1337 
1338 ! ! ! determine number of times
1339 ! ! nTime = SIZE(dataOutX, dim=1)
1340 
1341 ! ! !select variables to output
1342 ! ! nVar = COUNT((varList%level <= outLevel), dim=1)
1343 ! ! ALLOCATE (varListSel(nVar), stat=err)
1344 ! ! IF (err /= 0) PRINT *, "varListSel: Allocation request denied"
1345 ! ! varListSel = PACK(varList, mask=(varList%level <= outLevel))
1346 
1347 ! ! ! copy data accordingly
1348 ! ! ALLOCATE (dataOutSel(nTime, nVar, NumberOfGrids), stat=err)
1349 ! ! IF (err /= 0) PRINT *, "dataOutSel: Allocation request denied"
1350 ! ! dataOutSel = dataOutX(:, PACK((/(i, i=1, SIZE(varList))/), varList%level <= outLevel), :)
1351 
1352 ! ! ! determine filename
1353 ! ! CALL filename_gen(dataOutSel(:, :, 1), varListSel, 1, FileOut)
1354 ! ! ! PRINT*, 'writing file:',TRIM(fileOut)
1355 
1356 ! ! ! set year string
1357 ! ! WRITE (yrStr2, '(i4)') INT(dataOutX(1, 1, 1))
1358 ! ! ! get start for later time unit creation
1359 ! ! startStr2 = TRIM(yrStr2)//'-01-01 00:00:00'
1360 
1361 ! ! ! define the dimension of spatial array/frame in the output
1362 ! ! nX = nCol
1363 ! ! nY = nRow
1364 
1365 ! ! ALLOCATE (varSeq0(nX*nY))
1366 ! ! ALLOCATE (varSeq(nX*nY))
1367 ! ! ALLOCATE (xGridID(nX, nY))
1368 ! ! ALLOCATE (lon(nX, nY))
1369 ! ! ALLOCATE (lat(nX, nY))
1370 ! ! ALLOCATE (varY(nX, nY))
1371 ! ! ALLOCATE (varX(nX, nY))
1372 ! ! ALLOCATE (xTime(nTime))
1373 
1374 ! ! ! GridID:
1375 ! ! varSeq = SurfaceChar(1:nX*nY, 1)
1376 ! ! ! CALL sortSeqReal(varSeq0,varSeq,nY,nX)
1377 ! ! xGridID = RESHAPE(varSeq, (/nX, nY/), order=(/1, 2/))
1378 ! ! ! PRINT*, 'before flipping:',lat(1:2,1)
1379 ! ! xGridID = xGridID(:, nY:1:-1)
1380 
1381 ! ! ! latitude:
1382 ! ! varSeq = SurfaceChar(1:nX*nY, 5)
1383 ! ! ! CALL sortSeqReal(varSeq0,varSeq,nY,nX)
1384 ! ! lat = RESHAPE(varSeq, (/nX, nY/), order=(/1, 2/))
1385 ! ! ! PRINT*, 'before flipping:',lat(1:2,1)
1386 ! ! lat = lat(:, nY:1:-1)
1387 ! ! ! PRINT*, 'after flipping:',lat(1:2,1)
1388 
1389 ! ! ! longitude:
1390 ! ! varSeq = SurfaceChar(1:nX*nY, 6)
1391 ! ! ! CALL sortSeqReal(varSeq0,varSeq,nY,nX)
1392 ! ! lon = RESHAPE(varSeq, (/nX, nY/), order=(/1, 2/))
1393 ! ! lon = lon(:, nY:1:-1)
1394 
1395 ! ! ! pass values to coordinate variables
1396 ! ! varY = lat
1397 ! ! varX = lon
1398 
1399 ! ! ! calculate GeoTransform array as needed by GDAL
1400 ! ! ! ref: http://www.perrygeo.com/python-affine-transforms.html
1401 ! ! ! the values below are different from the above ref,
1402 ! ! ! as the layout of SUEWS output is different from the schematic shown there
1403 ! ! ! SUEWS output is arranged northward down the page
1404 ! ! ! if data are formatted as a normal matrix
1405 ! ! minLat = lat(1, 1) ! the lower-left pixel
1406 ! ! maxLat = lat(1, NY) ! the upper-left pixel
1407 ! ! IF (nY > 1) THEN
1408 ! ! dLat = (maxLat - minLat)/(nY - 1) ! height of a pixel
1409 ! ! ELSE
1410 ! ! dLat = 1
1411 ! ! END IF
1412 
1413 ! ! ! PRINT*, 'lat:',minLat,maxLat,dLat
1414 ! ! minLon = lon(1, 1) ! the lower-left pixel
1415 ! ! maxLon = lon(NX, 1) ! the lower-right pixel
1416 ! ! IF (nY > 1) THEN
1417 ! ! dLon = (maxLon - minLon)/(nX - 1) ! width of a pixel
1418 ! ! ELSE
1419 ! ! dLon = 1
1420 ! ! END IF
1421 
1422 ! ! ! PRINT*, 'lon:',minLon,maxLon,dLon
1423 ! ! geoTrans(1) = minLon - dLon/2 ! x-coordinate of the lower-left corner of the lower-left pixel
1424 ! ! geoTrans(2) = dLon ! width of a pixel
1425 ! ! geoTrans(3) = 0. ! row rotation (typically zero)
1426 ! ! geoTrans(4) = minLat - dLat/2 ! y-coordinate of the of the lower-left corner of the lower-left pixel
1427 ! ! geoTrans(5) = 0. ! column rotation (typically zero)
1428 ! ! geoTrans(6) = dLat ! height of a pixel (typically negative, but here positive)
1429 ! ! ! write GeoTransform to strGeoTrans
1430 ! ! WRITE (strGeoTrans, '(6(f12.8,1x))') geoTrans
1431 
1432 ! ! ! Create the netCDF file. The nf90_clobber parameter tells netCDF to
1433 ! ! ! overwrite this file, if it already exists.
1434 ! ! CALL check(nf90_create(TRIM(fileOut), NF90_CLOBBER, ncID))
1435 
1436 ! ! ! put global attributes
1437 ! ! CALL check(nf90_put_att(ncID, NF90_GLOBAL, 'Conventions', 'CF1.6'))
1438 ! ! CALL check(nf90_put_att(ncID, NF90_GLOBAL, 'title', 'SUEWS output'))
1439 ! ! CALL check(nf90_put_att(ncID, NF90_GLOBAL, 'source', 'Micromet Group, University of Reading'))
1440 ! ! CALL check(nf90_put_att(ncID, NF90_GLOBAL, 'references', 'http://urban-climate.net/umep/SUEWS'))
1441 
1442 ! ! ! Define the dimensions. NetCDF will hand back an ID for each.
1443 ! ! ! nY = ncolumnsDataOutSUEWS-4
1444 ! ! ! nx = NumberOfGrids
1445 ! ! CALL check(nf90_def_dim(ncID, "time", NF90_UNLIMITED, time_dimid))
1446 ! ! CALL check(nf90_def_dim(ncID, "west_east", NX, x_dimid))
1447 ! ! CALL check(nf90_def_dim(ncID, "south_north", NY, y_dimid))
1448 ! ! ! PRINT*, 'good define dim'
1449 
1450 ! ! ! The dimids array is used to pass the IDs of the dimensions of
1451 ! ! ! the variables. Note that in fortran arrays are stored in
1452 ! ! ! column-major format.
1453 ! ! dimids = (/x_dimid, y_dimid, time_dimid/)
1454 
1455 ! ! ! write out each variable
1456 ! ! ALLOCATE (varOut(nX, nY, nTime))
1457 
1458 ! ! ! define all variables
1459 ! ! ! define time variable:
1460 ! ! CALL check(nf90_def_var(ncID, 'time', NF90_REAL, time_dimid, varIDt))
1461 ! ! CALL check(nf90_put_att(ncID, varIDt, 'units', 'minutes since '//startStr2))
1462 ! ! CALL check(nf90_put_att(ncID, varIDt, 'long_name', 'time'))
1463 ! ! CALL check(nf90_put_att(ncID, varIDt, 'standard_name', 'time'))
1464 ! ! CALL check(nf90_put_att(ncID, varIDt, 'calendar', 'gregorian'))
1465 ! ! CALL check(nf90_put_att(ncID, varIDt, 'axis', 'T'))
1466 
1467 ! ! ! define coordinate variables:
1468 ! ! CALL check(nf90_def_var(ncID, 'lon', NF90_REAL, (/x_dimid, y_dimid/), varIDx))
1469 ! ! CALL check(nf90_put_att(ncID, varIDx, 'units', 'degree_east'))
1470 ! ! CALL check(nf90_put_att(ncID, varIDx, 'long_name', 'longitude'))
1471 ! ! CALL check(nf90_put_att(ncID, varIDx, 'standard_name', 'longitude'))
1472 ! ! CALL check(nf90_put_att(ncID, varIDx, 'axis', 'X'))
1473 
1474 ! ! CALL check(nf90_def_var(ncID, 'lat', NF90_REAL, (/x_dimid, y_dimid/), varIDy))
1475 ! ! CALL check(nf90_put_att(ncID, varIDy, 'units', 'degree_north'))
1476 ! ! CALL check(nf90_put_att(ncID, varIDy, 'long_name', 'latitude'))
1477 ! ! CALL check(nf90_put_att(ncID, varIDy, 'standard_name', 'latitude'))
1478 ! ! CALL check(nf90_put_att(ncID, varIDy, 'axis', 'Y'))
1479 
1480 ! ! ! define coordinate referencing system:
1481 ! ! CALL check(nf90_def_var(ncID, 'crsWGS84', NF90_INT, varIDCRS))
1482 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'grid_mapping_name', 'latitude_longitude'))
1483 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'long_name', 'CRS definition'))
1484 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'longitude_of_prime_meridian', '0.0'))
1485 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'semi_major_axis', '6378137.0'))
1486 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'inverse_flattening', '298.257223563'))
1487 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'epsg_code', 'EPSG:4326'))
1488 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'GeoTransform', TRIM(strGeoTrans)))
1489 ! ! CALL check(nf90_put_att(ncID, varIDCRS, 'spatial_ref',&
1490 ! ! &'GEOGCS["WGS 84",&
1491 ! ! & DATUM["WGS_1984",&
1492 ! ! & SPHEROID["WGS 84",6378137,298.257223563,&
1493 ! ! & AUTHORITY["EPSG","7030"]],&
1494 ! ! & AUTHORITY["EPSG","6326"]],&
1495 ! ! & PRIMEM["Greenwich",0],&
1496 ! ! & UNIT["degree",0.0174532925199433],&
1497 ! ! & AUTHORITY["EPSG","4326"]]' &
1498 ! ! ))
1499 
1500 ! ! ! define grid_ID:
1501 ! ! CALL check(nf90_def_var(ncID, 'grid_ID', NF90_INT, (/x_dimid, y_dimid/), varIDGrid))
1502 ! ! CALL check(nf90_put_att(ncID, varIDGrid, 'coordinates', 'lon lat'))
1503 ! ! CALL check(nf90_put_att(ncID, varIDGrid, 'long_name', 'Grid ID as in SiteSelect'))
1504 ! ! CALL check(nf90_put_att(ncID, varIDGrid, 'grid_mapping', 'crsWGS84'))
1505 ! ! ! varIDGrid=varID
1506 
1507 ! ! ! define other 3D variables:
1508 ! ! DO iVar = iVarStart, nVar
1509 ! ! ! define variable name
1510 ! ! header_str = varListSel(iVar)%header
1511 ! ! unit_str = varListSel(iVar)%unit
1512 ! ! longNm_str = varListSel(iVar)%longNm
1513 
1514 ! ! ! Define the variable. The type of the variable in this case is
1515 ! ! ! NF90_REAL.
1516 
1517 ! ! CALL check(nf90_def_var(ncID, TRIM(ADJUSTL(header_str)), NF90_REAL, dimids, varID))
1518 
1519 ! ! CALL check(nf90_put_att(ncID, varID, 'coordinates', 'lon lat'))
1520 
1521 ! ! CALL check(nf90_put_att(ncID, varID, 'units', TRIM(ADJUSTL(unit_str))))
1522 
1523 ! ! CALL check(nf90_put_att(ncID, varID, 'long_name', TRIM(ADJUSTL(longNm_str))))
1524 
1525 ! ! CALL check(nf90_put_att(ncID, varID, 'grid_mapping', 'crsWGS84'))
1526 
1527 ! ! idVar(iVar) = varID
1528 ! ! END DO
1529 ! ! CALL check(nf90_enddef(ncID))
1530 ! ! ! End define mode. This tells netCDF we are done defining metadata.
1531 
1532 ! ! ! put all variable values into netCDF datasets
1533 ! ! ! put time variable in minute:
1534 ! ! xTime = (dataOutSel(1:nTime, 2, 1) - 1)*24*60 + dataOutSel(1:nTime, 3, 1)*60 + dataOutSel(1:nTime, 4, 1)
1535 ! ! CALL check(nf90_put_var(ncID, varIDt, xTime))
1536 
1537 ! ! ! put coordinate variables:
1538 ! ! CALL check(nf90_put_var(ncID, varIDx, varX))
1539 ! ! CALL check(nf90_put_var(ncID, varIDy, varY))
1540 
1541 ! ! ! put CRS variable:
1542 ! ! CALL check(nf90_put_var(ncID, varIDCRS, 9999))
1543 
1544 ! ! CALL check(NF90_SYNC(ncID))
1545 ! ! ! PRINT*, 'good put var'
1546 
1547 ! ! ! put grid_ID:
1548 ! ! CALL check(nf90_put_var(ncID, varIDGrid, xGridID))
1549 ! ! ! PRINT*, 'good put varIDGrid',varIDGrid
1550 
1551 ! ! CALL check(NF90_SYNC(ncID))
1552 
1553 ! ! ! then other 3D variables
1554 ! ! DO iVar = iVarStart, nVar
1555 ! ! ! reshape dataOutX to be aligned in checker board form
1556 ! ! varOut = RESHAPE(dataOutSel(1:nTime, iVar, :), (/nX, nY, nTime/), order=(/3, 1, 2/))
1557 ! ! varOut = varOut(:, nY:1:-1, :)
1558 ! ! ! get the variable id
1559 ! ! varID = idVar(iVar)
1560 
1561 ! ! CALL check(nf90_put_var(ncID, varID, varOut))
1562 
1563 ! ! CALL check(NF90_SYNC(ncID))
1564 ! ! END DO
1565 
1566 ! ! IF (ALLOCATED(varOut)) DEALLOCATE (varOut)
1567 ! ! IF (ALLOCATED(varSeq0)) DEALLOCATE (varSeq0)
1568 ! ! IF (ALLOCATED(varSeq)) DEALLOCATE (varSeq)
1569 ! ! IF (ALLOCATED(xGridID)) DEALLOCATE (xGridID)
1570 ! ! IF (ALLOCATED(lon)) DEALLOCATE (lon)
1571 ! ! IF (ALLOCATED(lat)) DEALLOCATE (lat)
1572 ! ! IF (ALLOCATED(varY)) DEALLOCATE (varY)
1573 ! ! IF (ALLOCATED(varX)) DEALLOCATE (varX)
1574 ! ! IF (ALLOCATED(xTime)) DEALLOCATE (xTime)
1575 
1576 ! ! ! Close the file. This frees up any internal netCDF resources
1577 ! ! ! associated with the file, and flushes any buffers.
1578 ! ! CALL check(nf90_close(ncID))
1579 
1580 ! ! ! PRINT*, "*** SUCCESS writing netCDF file:"
1581 ! ! ! PRINT*, FileOut
1582 ! ! END SUBROUTINE SUEWS_Write_nc
1583 
1584 ! !===========================================================================!
1585 ! ! convert a vector of grids to a matrix
1586 ! ! the grid IDs in seqGrid2Sort follow the QGIS convention
1587 ! ! the spatial matrix arranges successive rows down the page (i.e., north to south)
1588 ! ! and succesive columns across (i.e., west to east)
1589 ! ! seqGridSorted stores the grid IDs as aligned in matGrid but squeezed into a vector
1590 ! !===========================================================================!
1591 ! SUBROUTINE grid2mat(seqGrid2Sort, seqGridSorted, matGrid, nRow, nCol)
1592 
1593 ! IMPLICIT NONE
1594 
1595 ! INTEGER, DIMENSION(nRow*nCol) :: seqGrid2Sort, seqGridSorted
1596 ! INTEGER, DIMENSION(nRow, nCol) :: matGrid
1597 ! INTEGER :: nRow, nCol, i, j, loc
1598 
1599 ! CALL sortGrid(seqGrid2Sort, seqGridSorted, nRow, nCol)
1600 ! PRINT *, 'old:'
1601 ! PRINT *, seqGrid2Sort(1:5)
1602 ! PRINT *, 'sorted:'
1603 ! PRINT *, seqGridSorted(1:5)
1604 ! PRINT *, ''
1605 ! DO i = 1, nRow
1606 ! DO j = 1, nCol
1607 ! loc = (i - 1)*nCol + j
1608 ! matGrid(i, j) = seqGridSorted(loc)
1609 ! END DO
1610 ! END DO
1611 ! END SUBROUTINE grid2mat
1612 
1613 ! !===========================================================================!
1614 ! ! convert sequence of REAL values to a matrix
1615 ! ! the grid IDs in seqGrid2Sort follow the QGIS convention
1616 ! ! the spatial matrix arranges successive rows down the page (i.e., north to south)
1617 ! ! and succesive columns across (i.e., west to east)
1618 ! ! seqGridSorted stores the grid IDs as aligned in matGrid but squeezed into a vector
1619 ! !===========================================================================!
1620 ! SUBROUTINE seq2mat(seq2Sort, seqSorted, matGrid, nRow, nCol)
1621 
1622 ! IMPLICIT NONE
1623 
1624 ! REAL(KIND(1d0)), DIMENSION(nRow*nCol) :: seq2Sort, seqSorted
1625 ! REAL(KIND(1d0)), DIMENSION(nRow, nCol) :: matGrid
1626 ! INTEGER :: nRow, nCol, i, j, loc
1627 
1628 ! CALL sortSeqReal(seq2Sort, seqSorted, nRow, nCol)
1629 ! PRINT *, 'old:'
1630 ! PRINT *, seq2Sort(1:5)
1631 ! PRINT *, 'sorted:'
1632 ! PRINT *, seqSorted(1:5)
1633 ! PRINT *, ''
1634 ! DO i = 1, nRow
1635 ! DO j = 1, nCol
1636 ! loc = (i - 1)*nCol + j
1637 ! matGrid(i, j) = seqSorted(loc)
1638 ! END DO
1639 ! END DO
1640 ! END SUBROUTINE seq2mat
1641 
1642 ! !===========================================================================!
1643 ! ! sort a sequence of LONG values into the specially aligned sequence per QGIS
1644 ! !===========================================================================!
1645 ! SUBROUTINE sortGrid(seqGrid2Sort0, seqGridSorted, nRow, nCol)
1646 ! USE qsort_c_module
1647 ! ! convert a vector of grids to a matrix
1648 ! ! the grid IDs in seqGrid2Sort follow the QGIS convention
1649 ! ! the spatial matrix arranges successive rows down the page (i.e., north to south)
1650 ! ! and succesive columns across (i.e., west to east)
1651 ! ! seqGridSorted stores the grid IDs as aligned in matGrid but squeezed into a vector
1652 
1653 ! IMPLICIT NONE
1654 ! INTEGER :: nRow, nCol, i = 1, j = 1, xInd, len
1655 
1656 ! INTEGER, DIMENSION(nRow*nCol), INTENT(in) :: seqGrid2Sort0
1657 ! INTEGER, DIMENSION(nRow*nCol), INTENT(out) :: seqGridSorted
1658 ! INTEGER, DIMENSION(nRow*nCol) :: seqGrid2Sort, locSorted
1659 ! INTEGER :: loc
1660 ! REAL:: ind(nRow*nCol, 2)
1661 ! REAL, DIMENSION(nRow*nCol) :: seqGrid2SortReal, seqGridSortedReal
1662 ! REAL :: val
1663 
1664 ! ! number of grids
1665 ! len = nRow*nCol
1666 
1667 ! !sort the input array to make sure the grid order is in QGIS convention
1668 ! ! i.e., diagonally ascending
1669 ! seqGrid2SortReal = seqGrid2Sort0*1.
1670 ! CALL QsortC(seqGrid2SortReal)
1671 ! seqGrid2Sort = INT(seqGrid2SortReal)
1672 
1673 ! ! fill in an nRow*nCol array with values to determine sequence
1674 ! xInd = 1
1675 ! DO i = 1, nRow
1676 ! DO j = 1, nCol
1677 ! ! {row, col, value for sorting, index in new sequence}
1678 ! ind(xInd, :) = (/i + j + i/(nRow + 1.), xInd*1./)
1679 ! xInd = xInd + 1
1680 ! END DO
1681 ! END DO
1682 
1683 ! ! then sorted ind(:,3) will have the same order as seqGrid2Sort
1684 ! ! sort ind(:,3)
1685 ! seqGridSortedReal = ind(:, 1)*1.
1686 ! CALL QsortC(seqGridSortedReal)
1687 ! ! print*, 'sorted real:'
1688 ! ! print*, seqGridSortedReal
1689 
1690 ! ! get index of each element of old sequence in the sorted sequence
1691 ! DO i = 1, len
1692 ! ! value in old sequence
1693 ! ! val=ind(i,3)*1.
1694 ! val = seqGridSortedReal(i)
1695 ! DO j = 1, len
1696 ! IF (val == ind(j, 1)*1.) THEN
1697 ! ! location in sorted sequence
1698 ! locSorted(i) = j
1699 ! END IF
1700 ! END DO
1701 ! END DO
1702 
1703 ! ! put elements of old sequence in the sorted order
1704 ! DO i = 1, len
1705 ! loc = locSorted(i)
1706 ! seqGridSorted(loc) = seqGrid2Sort(i)
1707 ! END DO
1708 ! seqGridSorted = seqGridSorted(len:1:-1)
1709 
1710 ! END SUBROUTINE sortGrid
1711 
1712 ! !===========================================================================!
1713 ! ! sort a sequence of REAL values into the specially aligned sequence per QGIS
1714 ! !===========================================================================!
1715 ! SUBROUTINE sortSeqReal(seqReal2Sort, seqRealSorted, nRow, nCol)
1716 ! USE qsort_c_module
1717 ! ! convert a vector of grids to a matrix
1718 ! ! the grid IDs in seqReal2Sort follow the QGIS convention
1719 ! ! the spatial matrix arranges successive rows down the page (i.e., north to south)
1720 ! ! and succesive columns across (i.e., west to east)
1721 ! ! seqRealSorted stores the grid IDs as aligned in matGrid but squeezed into a vector
1722 
1723 ! IMPLICIT NONE
1724 ! INTEGER :: nRow, nCol, i = 1, j = 1, xInd, len
1725 
1726 ! REAL(KIND(1d0)), DIMENSION(nRow*nCol), INTENT(in) :: seqReal2Sort
1727 ! REAL(KIND(1d0)), DIMENSION(nRow*nCol), INTENT(out) :: seqRealSorted
1728 ! INTEGER(KIND(1d0)), DIMENSION(nRow*nCol) :: locSorted
1729 ! INTEGER(KIND(1d0)) :: loc
1730 ! REAL:: ind(nRow*nCol, 2)
1731 ! REAL :: seqRealSortedReal(nRow*nCol), val
1732 
1733 ! ! number of grids
1734 ! len = nRow*nCol
1735 
1736 ! ! fill in an nRow*nCol array with values to determine sequence
1737 ! xInd = 1
1738 ! DO i = 1, nRow
1739 ! DO j = 1, nCol
1740 ! ! {row, col, value for sorting, index in new sequence}
1741 ! ind(xInd, :) = (/i + j + i/(nRow + 1.), xInd*1./)
1742 ! xInd = xInd + 1
1743 ! END DO
1744 ! END DO
1745 
1746 ! ! then sorted ind(:,3) will have the same order as seqReal2Sort
1747 ! ! sort ind(:,3)
1748 ! seqRealSortedReal = ind(:, 1)*1.
1749 ! CALL QsortC(seqRealSortedReal)
1750 ! ! print*, 'sorted real:'
1751 ! ! print*, seqRealSortedReal
1752 
1753 ! ! get index of each element of old sequence in the sorted sequence
1754 ! DO i = 1, len
1755 ! ! value in old sequence
1756 ! ! val=ind(i,3)*1.
1757 ! val = seqRealSortedReal(i)
1758 ! DO j = 1, len
1759 ! IF (val == ind(j, 1)*1.) THEN
1760 ! ! location in sorted sequence
1761 ! locSorted(i) = j
1762 ! END IF
1763 ! END DO
1764 ! END DO
1765 
1766 ! ! put elements of old sequence in the sorted order
1767 ! DO i = 1, len
1768 ! loc = locSorted(i)
1769 ! seqRealSorted(loc) = seqReal2Sort(i)
1770 ! END DO
1771 ! seqRealSorted = seqRealSorted(len:1:-1)
1772 
1773 ! END SUBROUTINE sortSeqReal
1774 
1775 ! !===========================================================================!
1776 ! ! a wrapper for checking netCDF status
1777 ! !===========================================================================!
1778 
1779 ! SUBROUTINE check(status)
1780 ! USE netcdf
1781 ! IMPLICIT NONE
1782 
1783 ! INTEGER, INTENT(in) :: status
1784 
1785 ! IF (status /= nf90_noerr) THEN
1786 ! PRINT *, TRIM(nf90_strerror(status))
1787 ! STOP "Stopped"
1788 ! END IF
1789 ! END SUBROUTINE check
1790 ! #endif
1791 
1792 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)
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
real(kind(1d0)), dimension(:, :, :), allocatable dataoutsolweig
subroutine suews_output_agg(dataOut_agg, dataOutX, varList, irMax, outFreq_s)