-
Notifications
You must be signed in to change notification settings - Fork 4
/
uEMEP_crossreference_grids.f90
482 lines (411 loc) · 32.1 KB
/
uEMEP_crossreference_grids.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
!uEMEP_crossreference_grids
subroutine uEMEP_crossreference_grids
use uEMEP_definitions
implicit none
integer i,j,k
integer ii,jj
integer i_source
real x_temp,y_temp
!Cross referencing must be done for each new grid when using multiple grids
if (allocated(crossreference_target_to_emep_subgrid)) deallocate (crossreference_target_to_emep_subgrid)
if (allocated(crossreference_target_to_localfraction_subgrid)) deallocate (crossreference_target_to_localfraction_subgrid)
if (allocated(crossreference_integral_to_emep_subgrid)) deallocate (crossreference_integral_to_emep_subgrid)
if (allocated(crossreference_target_to_integral_subgrid)) deallocate (crossreference_target_to_integral_subgrid)
if (allocated(crossreference_target_to_emission_subgrid)) deallocate (crossreference_target_to_emission_subgrid)
if (allocated(crossreference_emission_to_EMEP_subgrid)) deallocate (crossreference_emission_to_EMEP_subgrid)
if (allocated(crossreference_integral_to_emission_subgrid)) deallocate (crossreference_integral_to_emission_subgrid)
if (allocated(crossreference_emission_to_integral_subgrid)) deallocate (crossreference_emission_to_integral_subgrid)
if (allocated(crossreference_target_to_population_subgrid)) deallocate (crossreference_target_to_population_subgrid)
if (use_alternative_meteorology_flag) then
if (allocated(crossreference_integral_to_meteo_nc_subgrid)) deallocate (crossreference_integral_to_meteo_nc_subgrid)
endif
if (calculate_deposition_flag) then
if (allocated(crossreference_emission_to_deposition_subgrid)) deallocate (crossreference_emission_to_deposition_subgrid)
if (allocated(crossreference_target_to_deposition_subgrid)) deallocate (crossreference_target_to_deposition_subgrid)
if (allocated(crossreference_deposition_to_emep_subgrid)) deallocate (crossreference_deposition_to_emep_subgrid)
endif
if (read_landuse_flag) then
if (allocated(crossreference_emission_to_landuse_subgrid)) deallocate (crossreference_emission_to_landuse_subgrid)
endif
allocate (crossreference_target_to_emep_subgrid(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),2))
allocate (crossreference_target_to_localfraction_subgrid(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),2,n_local_fraction_grids))
allocate (crossreference_integral_to_emep_subgrid(integral_subgrid_dim(x_dim_index),integral_subgrid_dim(y_dim_index),2))
allocate (crossreference_target_to_integral_subgrid(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),2))
allocate (crossreference_target_to_emission_subgrid(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),2,n_source_index))
allocate (crossreference_emission_to_EMEP_subgrid(emission_max_subgrid_dim(x_dim_index),emission_max_subgrid_dim(y_dim_index),2,n_source_index))
allocate (crossreference_integral_to_emission_subgrid(integral_subgrid_dim(x_dim_index),integral_subgrid_dim(y_dim_index),2,n_source_index))
allocate (crossreference_emission_to_integral_subgrid(emission_max_subgrid_dim(x_dim_index),emission_max_subgrid_dim(y_dim_index),2,n_source_index))
allocate (crossreference_target_to_population_subgrid(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),2))
if (use_alternative_meteorology_flag) then
allocate (crossreference_integral_to_meteo_nc_subgrid(integral_subgrid_dim(x_dim_index),integral_subgrid_dim(y_dim_index),2))
endif
if (calculate_deposition_flag) then
allocate (crossreference_emission_to_deposition_subgrid(emission_max_subgrid_dim(x_dim_index),emission_max_subgrid_dim(y_dim_index),2,n_source_index))
allocate (crossreference_target_to_deposition_subgrid(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),2))
allocate (crossreference_deposition_to_emep_subgrid(deposition_subgrid_dim(x_dim_index),deposition_subgrid_dim(y_dim_index),2))
endif
if (read_landuse_flag) then
allocate (crossreference_emission_to_landuse_subgrid(emission_max_subgrid_dim(x_dim_index),emission_max_subgrid_dim(y_dim_index),2,n_source_index))
endif
write(unit_logfile,'(A)')'Allocating EMEP grid index to subgrid index'
!Loop through subgrid and find those subgrids within EMEP grids and allocate concentrations directly from EMEP grids.
do j=1,subgrid_dim(y_dim_index)
do i=1,subgrid_dim(x_dim_index)
if (EMEP_projection_type.eq.LL_projection_index) then
ii=1+floor((lon_subgrid(i,j)-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((lat_subgrid(i,j)-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)+0.5)
elseif (EMEP_projection_type.eq.LCC_projection_index) then
!When EMEP is read as x,y projection then var1d_nc(:,lon/lat_nc_index) are the x, y projection indexes, actually
!if (use_alternative_LCC_projection_flag) then
call lb2lambert2_uEMEP(x_temp,y_temp,lon_subgrid(i,j),lat_subgrid(i,j),EMEP_projection_attributes)
!else
! call lb2lambert_uEMEP(x_temp,y_temp,lon_subgrid(i,j),lat_subgrid(i,j),real(EMEP_projection_attributes(3)),real(EMEP_projection_attributes(4)))
!endif
ii=1+floor((x_temp-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((y_temp-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)+0.5)
elseif (EMEP_projection_type.eq.PS_projection_index) then
call LL2PS_spherical(x_temp,y_temp,lon_subgrid(i,j),lat_subgrid(i,j),EMEP_projection_attributes)
ii=1+floor((x_temp-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((y_temp-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)+0.5)
else
write(unit_logfile,'(A)')'No valid projection in use. Stopping'
stop
endif
crossreference_target_to_emep_subgrid(i,j,x_dim_index)=ii
crossreference_target_to_emep_subgrid(i,j,y_dim_index)=jj
enddo
enddo
write(unit_logfile,'(A)')'Allocating EMEP local fraction grid index to subgrid index'
!Loop through subgrid and find those subgrids within EMEP grids and allocate concentrations directly from EMEP grids.
do k=1,n_local_fraction_grids
do j=1,subgrid_dim(y_dim_index)
do i=1,subgrid_dim(x_dim_index)
if (EMEP_projection_type.eq.LL_projection_index) then
ii=1+floor((lon_subgrid(i,j)-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)/local_fraction_grid_size(k)+0.5)
jj=1+floor((lat_subgrid(i,j)-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)/local_fraction_grid_size(k)+0.5)
elseif (EMEP_projection_type.eq.LCC_projection_index) then
!When EMEP is read as x,y projection then var1d_nc(:,lon/lat_nc_index) are the x, y projection indexes, actually
!if (use_alternative_LCC_projection_flag) then
call lb2lambert2_uEMEP(x_temp,y_temp,lon_subgrid(i,j),lat_subgrid(i,j),EMEP_projection_attributes)
!else
! call lb2lambert_uEMEP(x_temp,y_temp,lon_subgrid(i,j),lat_subgrid(i,j),real(EMEP_projection_attributes(3)),real(EMEP_projection_attributes(4)))
!endif
ii=1+floor((x_temp-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)/local_fraction_grid_size(k)+0.5)
jj=1+floor((y_temp-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)/local_fraction_grid_size(k)+0.5)
elseif (EMEP_projection_type.eq.PS_projection_index) then
call LL2PS_spherical(x_temp,y_temp,lon_subgrid(i,j),lat_subgrid(i,j),EMEP_projection_attributes)
ii=1+floor((x_temp-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)/local_fraction_grid_size(k)+0.5)
jj=1+floor((y_temp-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)/local_fraction_grid_size(k)+0.5)
else
write(unit_logfile,'(A)')'No valid projection in use. Stopping'
stop
endif
crossreference_target_to_localfraction_subgrid(i,j,x_dim_index,k)=ii
crossreference_target_to_localfraction_subgrid(i,j,y_dim_index,k)=jj
enddo
enddo
enddo
write(unit_logfile,'(A)')'Allocating integral grid index to subgrid index'
do j=1,subgrid_dim(y_dim_index)
do i=1,subgrid_dim(x_dim_index)
crossreference_target_to_integral_subgrid(i,j,x_dim_index)=1+floor((x_subgrid(i,j)-integral_subgrid_min(x_dim_index))/integral_subgrid_delta(x_dim_index))
crossreference_target_to_integral_subgrid(i,j,y_dim_index)=1+floor((y_subgrid(i,j)-integral_subgrid_min(y_dim_index))/integral_subgrid_delta(y_dim_index))
enddo
enddo
write(unit_logfile,'(A)')'Allocating population grid index to subgrid index'
do j=1,subgrid_dim(y_dim_index)
do i=1,subgrid_dim(x_dim_index)
crossreference_target_to_population_subgrid(i,j,x_dim_index)=1+floor((x_subgrid(i,j)-population_subgrid_min(x_dim_index))/population_subgrid_delta(x_dim_index))
crossreference_target_to_population_subgrid(i,j,y_dim_index)=1+floor((y_subgrid(i,j)-population_subgrid_min(y_dim_index))/population_subgrid_delta(y_dim_index))
enddo
enddo
write(unit_logfile,'(A)')'Allocating EMEP grid index to integral subgrid index'
do j=1,integral_subgrid_dim(y_dim_index)
do i=1,integral_subgrid_dim(x_dim_index)
if (EMEP_projection_type.eq.LL_projection_index) then
ii=1+floor((lon_integral_subgrid(i,j)-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((lat_integral_subgrid(i,j)-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)+0.5)
elseif (EMEP_projection_type.eq.LCC_projection_index) then
call lb2lambert2_uEMEP(x_temp,y_temp,lon_integral_subgrid(i,j),lat_integral_subgrid(i,j),EMEP_projection_attributes)
ii=1+floor((x_temp-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((y_temp-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)+0.5)
elseif (EMEP_projection_type.eq.PS_projection_index) then
call LL2PS_spherical(x_temp,y_temp,lon_integral_subgrid(i,j),lat_integral_subgrid(i,j),EMEP_projection_attributes)
ii=1+floor((x_temp-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((y_temp-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)+0.5)
else
write(unit_logfile,'(A)')'No valid projection in use. Stopping'
stop
endif
crossreference_integral_to_emep_subgrid(i,j,x_dim_index)=ii
crossreference_integral_to_emep_subgrid(i,j,y_dim_index)=jj
enddo
enddo
if (use_alternative_meteorology_flag) then
write(unit_logfile,'(A)')'Allocating alternative meteo nc grid index to integral subgrid index'
do j=1,integral_subgrid_dim(y_dim_index)
do i=1,integral_subgrid_dim(x_dim_index)
if (meteo_nc_projection_type.eq.LL_projection_index) then
ii=1+floor((lon_integral_subgrid(i,j)-meteo_var1d_nc(1,lon_nc_index))/meteo_dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((lat_integral_subgrid(i,j)-meteo_var1d_nc(1,lat_nc_index))/meteo_dgrid_nc(lat_nc_index)+0.5)
elseif (meteo_nc_projection_type.eq.LCC_projection_index) then
call lb2lambert2_uEMEP(x_temp,y_temp,lon_integral_subgrid(i,j),lat_integral_subgrid(i,j),meteo_nc_projection_attributes)
ii=1+floor((x_temp-meteo_var1d_nc(1,lon_nc_index))/meteo_dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((y_temp-meteo_var1d_nc(1,lat_nc_index))/meteo_dgrid_nc(lat_nc_index)+0.5)
elseif (meteo_nc_projection_type.eq.PS_projection_index) then
call LL2PS_spherical(x_temp,y_temp,lon_integral_subgrid(i,j),lat_integral_subgrid(i,j),meteo_nc_projection_attributes)
ii=1+floor((x_temp-meteo_var1d_nc(1,lon_nc_index))/meteo_dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((y_temp-meteo_var1d_nc(1,lat_nc_index))/meteo_dgrid_nc(lat_nc_index)+0.5)
else
write(unit_logfile,'(A)')'No valid projection in use. Stopping'
stop
endif
crossreference_integral_to_meteo_nc_subgrid(i,j,x_dim_index)=ii
crossreference_integral_to_meteo_nc_subgrid(i,j,y_dim_index)=jj
enddo
enddo
endif
do i_source=1,n_source_index
if (calculate_source(i_source)) then
!do j=1,emission_subgrid_dim(y_dim_index,i_source)
!do i=1,emission_subgrid_dim(x_dim_index,i_source)
! crossreference_emission_to_target_subgrid(i,j,x_dim_index,i_source)=1+floor((x_emission_subgrid(i,j,i_source)-subgrid_min(x_dim_index))/subgrid_delta(x_dim_index)+0.5)
! crossreference_emission_to_target_subgrid(i,j,y_dim_index,i_source)=1+floor((y_emission_subgrid(i,j,i_source)-subgrid_min(y_dim_index))/subgrid_delta(y_dim_index)+0.5)
!enddo
!enddo
do j=1,emission_subgrid_dim(y_dim_index,i_source)
do i=1,emission_subgrid_dim(x_dim_index,i_source)
if (EMEP_projection_type.eq.LL_projection_index) then
ii=1+floor((lon_emission_subgrid(i,j,i_source)-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((lat_emission_subgrid(i,j,i_source)-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)+0.5)
elseif (EMEP_projection_type.eq.LCC_projection_index) then
call lb2lambert2_uEMEP(x_temp,y_temp,lon_emission_subgrid(i,j,i_source),lat_emission_subgrid(i,j,i_source),EMEP_projection_attributes)
ii=1+floor((x_temp-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((y_temp-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)+0.5)
elseif (EMEP_projection_type.eq.PS_projection_index) then
call LL2PS_spherical(x_temp,y_temp,lon_emission_subgrid(i,j,i_source),lat_emission_subgrid(i,j,i_source),EMEP_projection_attributes)
ii=1+floor((x_temp-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((y_temp-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)+0.5)
else
write(unit_logfile,'(A)')'No valid projection in use. Stopping'
stop
endif
crossreference_emission_to_emep_subgrid(i,j,x_dim_index,i_source)=ii
crossreference_emission_to_emep_subgrid(i,j,y_dim_index,i_source)=jj
enddo
enddo
do j=1,subgrid_dim(y_dim_index)
do i=1,subgrid_dim(x_dim_index)
crossreference_target_to_emission_subgrid(i,j,x_dim_index,i_source)=1+floor((x_subgrid(i,j)-emission_subgrid_min(x_dim_index,i_source))/emission_subgrid_delta(x_dim_index,i_source))
crossreference_target_to_emission_subgrid(i,j,y_dim_index,i_source)=1+floor((y_subgrid(i,j)-emission_subgrid_min(y_dim_index,i_source))/emission_subgrid_delta(y_dim_index,i_source))
enddo
enddo
do j=1,emission_subgrid_dim(y_dim_index,i_source)
do i=1,emission_subgrid_dim(x_dim_index,i_source)
crossreference_emission_to_integral_subgrid(i,j,x_dim_index,i_source)=1+floor((x_emission_subgrid(i,j,i_source)-integral_subgrid_min(x_dim_index))/integral_subgrid_delta(x_dim_index))
crossreference_emission_to_integral_subgrid(i,j,y_dim_index,i_source)=1+floor((y_emission_subgrid(i,j,i_source)-integral_subgrid_min(y_dim_index))/integral_subgrid_delta(y_dim_index))
!At edge this can return negative distances due to the different sizes of emission and integral grids and buffer zones. Set the limits here. Should not be a problem
crossreference_emission_to_integral_subgrid(i,j,x_dim_index,i_source)=max(min(crossreference_emission_to_integral_subgrid(i,j,x_dim_index,i_source),integral_subgrid_dim(x_dim_index)),1)
crossreference_emission_to_integral_subgrid(i,j,y_dim_index,i_source)=max(min(crossreference_emission_to_integral_subgrid(i,j,y_dim_index,i_source),integral_subgrid_dim(y_dim_index)),1)
if (crossreference_emission_to_integral_subgrid(i,j,x_dim_index,i_source).lt.1.or.crossreference_emission_to_integral_subgrid(i,j,x_dim_index,i_source).gt.integral_subgrid_dim(x_dim_index) &
.or.crossreference_emission_to_integral_subgrid(i,j,y_dim_index,i_source).lt.1.or.crossreference_emission_to_integral_subgrid(i,j,y_dim_index,i_source).gt.integral_subgrid_dim(y_dim_index)) then
write(unit_logfile,'(A,4i,4f)') 'WARNING: crossreference_emission_to_integral_subgrid is out of bounds (i_emis,j_emis,i_integral,j_integral,x_emis,y_emis)',i,j, &
crossreference_emission_to_integral_subgrid(i,j,x_dim_index,i_source),crossreference_emission_to_integral_subgrid(i,j,y_dim_index,i_source) &
,x_emission_subgrid(i,j,i_source)/1000,y_emission_subgrid(i,j,i_source)/1000 &
,(x_emission_subgrid(i,j,i_source)-integral_subgrid_min(x_dim_index))/integral_subgrid_delta(x_dim_index)+0.5,(y_emission_subgrid(i,j,i_source)-integral_subgrid_min(y_dim_index))/integral_subgrid_delta(y_dim_index)+0.5
endif
enddo
enddo
do j=1,integral_subgrid_dim(y_dim_index)
do i=1,integral_subgrid_dim(x_dim_index)
crossreference_integral_to_emission_subgrid(i,j,x_dim_index,i_source)=1+floor((x_integral_subgrid(i,j)-emission_subgrid_min(x_dim_index,i_source))/emission_subgrid_delta(x_dim_index,i_source))
crossreference_integral_to_emission_subgrid(i,j,y_dim_index,i_source)=1+floor((y_integral_subgrid(i,j)-emission_subgrid_min(y_dim_index,i_source))/emission_subgrid_delta(y_dim_index,i_source))
enddo
enddo
if (calculate_deposition_flag) then
do j=1,emission_subgrid_dim(y_dim_index,i_source)
do i=1,emission_subgrid_dim(x_dim_index,i_source)
crossreference_emission_to_deposition_subgrid(i,j,x_dim_index,i_source)=1+floor((x_emission_subgrid(i,j,i_source)-deposition_subgrid_min(x_dim_index))/deposition_subgrid_delta(x_dim_index))
crossreference_emission_to_deposition_subgrid(i,j,y_dim_index,i_source)=1+floor((y_emission_subgrid(i,j,i_source)-deposition_subgrid_min(y_dim_index))/deposition_subgrid_delta(y_dim_index))
!At edge this can return negative distances due to the different sizes of emission and integral grids and buffer zones. Set the limits here. Should not be a problem
crossreference_emission_to_deposition_subgrid(i,j,x_dim_index,i_source)=max(min(crossreference_emission_to_deposition_subgrid(i,j,x_dim_index,i_source),deposition_subgrid_dim(x_dim_index)),1)
crossreference_emission_to_deposition_subgrid(i,j,y_dim_index,i_source)=max(min(crossreference_emission_to_deposition_subgrid(i,j,y_dim_index,i_source),deposition_subgrid_dim(y_dim_index)),1)
enddo
enddo
endif
if (read_landuse_flag) then
do j=1,emission_subgrid_dim(y_dim_index,i_source)
do i=1,emission_subgrid_dim(x_dim_index,i_source)
crossreference_emission_to_landuse_subgrid(i,j,x_dim_index,i_source)=1+floor((x_emission_subgrid(i,j,i_source)-landuse_subgrid_min(x_dim_index))/landuse_subgrid_delta(x_dim_index))
crossreference_emission_to_landuse_subgrid(i,j,y_dim_index,i_source)=1+floor((y_emission_subgrid(i,j,i_source)-landuse_subgrid_min(y_dim_index))/landuse_subgrid_delta(y_dim_index))
!At edge this can return negative distances due to the different sizes of emission and integral grids and buffer zones. Set the limits here. Should not be a problem
crossreference_emission_to_landuse_subgrid(i,j,x_dim_index,i_source)=max(min(crossreference_emission_to_landuse_subgrid(i,j,x_dim_index,i_source),landuse_subgrid_dim(x_dim_index)),1)
crossreference_emission_to_landuse_subgrid(i,j,y_dim_index,i_source)=max(min(crossreference_emission_to_landuse_subgrid(i,j,y_dim_index,i_source),landuse_subgrid_dim(y_dim_index)),1)
enddo
enddo
endif
endif
enddo
if (calculate_deposition_flag) then
do j=1,subgrid_dim(y_dim_index)
do i=1,subgrid_dim(x_dim_index)
crossreference_target_to_deposition_subgrid(i,j,x_dim_index)=1+floor((x_subgrid(i,j)-deposition_subgrid_min(x_dim_index))/deposition_subgrid_delta(x_dim_index))
crossreference_target_to_deposition_subgrid(i,j,y_dim_index)=1+floor((y_subgrid(i,j)-deposition_subgrid_min(y_dim_index))/deposition_subgrid_delta(y_dim_index))
!At edge this can return negative distances due to the different sizes of emission and integral grids and buffer zones. Set the limits here. Should not be a problem
crossreference_target_to_deposition_subgrid(i,j,x_dim_index)=max(min(crossreference_target_to_deposition_subgrid(i,j,x_dim_index),deposition_subgrid_dim(x_dim_index)),1)
crossreference_target_to_deposition_subgrid(i,j,y_dim_index)=max(min(crossreference_target_to_deposition_subgrid(i,j,y_dim_index),deposition_subgrid_dim(y_dim_index)),1)
enddo
enddo
write(unit_logfile,'(A)')'Allocating EMEP grid index to deposition subgrid index'
do j=1,deposition_subgrid_dim(y_dim_index)
do i=1,deposition_subgrid_dim(x_dim_index)
if (EMEP_projection_type.eq.LL_projection_index) then
ii=1+floor((lon_deposition_subgrid(i,j)-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((lat_deposition_subgrid(i,j)-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)+0.5)
elseif (EMEP_projection_type.eq.LCC_projection_index) then
call lb2lambert2_uEMEP(x_temp,y_temp,lon_deposition_subgrid(i,j),lat_deposition_subgrid(i,j),EMEP_projection_attributes)
ii=1+floor((x_temp-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((y_temp-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)+0.5)
elseif (EMEP_projection_type.eq.PS_projection_index) then
call LL2PS_spherical(x_temp,y_temp,lon_deposition_subgrid(i,j),lat_deposition_subgrid(i,j),EMEP_projection_attributes)
ii=1+floor((x_temp-var1d_nc(1,lon_nc_index))/dgrid_nc(lon_nc_index)+0.5)
jj=1+floor((y_temp-var1d_nc(1,lat_nc_index))/dgrid_nc(lat_nc_index)+0.5)
else
write(unit_logfile,'(A)')'No valid projection in use. Stopping'
stop
endif
crossreference_deposition_to_emep_subgrid(i,j,x_dim_index)=ii
crossreference_deposition_to_emep_subgrid(i,j,y_dim_index)=jj
enddo
enddo
endif
end subroutine uEMEP_crossreference_grids
subroutine uEMEP_assign_region_coverage_to_EMEP
use uEMEP_definitions
implicit none
integer i,j
integer ii,jj,iii,jjj,iiii,jjjj,ii_nc,jj_nc
real, allocatable :: count_EMEP_grid_fraction_in_region(:,:,:)
integer i_source, chosen_source
integer ii_start,jj_start
write(unit_logfile,'(A)') ''
write(unit_logfile,'(A)') '================================================================'
write(unit_logfile,'(A)') 'Assigning regional coverage to EMEP grids (uEMEP_assign_region_coverage_to_EMEP)'
write(unit_logfile,'(A)') '================================================================'
! xdist_centre_nc=1+dim_length_nc(xdist_dim_nc_index)/2
! ydist_centre_nc=1+dim_length_nc(ydist_dim_nc_index)/2
!allocate the fraction of each EMEP grid that is within a region, if required
if (trace_emissions_from_in_region) then
if (allocated(EMEP_grid_fraction_in_region)) deallocate (EMEP_grid_fraction_in_region)
if (.not.allocated(EMEP_grid_fraction_in_region)) allocate (EMEP_grid_fraction_in_region(dim_length_nc(x_dim_nc_index),dim_length_nc(y_dim_nc_index),n_source_nc_index,2))
EMEP_grid_fraction_in_region=0
if (allocated(lf_EMEP_grid_fraction_in_region)) deallocate (lf_EMEP_grid_fraction_in_region)
if (.not.allocated(lf_EMEP_grid_fraction_in_region)) allocate (lf_EMEP_grid_fraction_in_region(dim_length_nc(xdist_dim_nc_index),dim_length_nc(ydist_dim_nc_index),dim_length_nc(x_dim_nc_index),dim_length_nc(y_dim_nc_index),n_source_nc_index,2))
lf_EMEP_grid_fraction_in_region=0
if (allocated(count_EMEP_grid_fraction_in_region)) deallocate (count_EMEP_grid_fraction_in_region)
if (.not.allocated(count_EMEP_grid_fraction_in_region)) allocate (count_EMEP_grid_fraction_in_region(dim_length_nc(x_dim_nc_index),dim_length_nc(y_dim_nc_index),n_source_nc_index))
count_EMEP_grid_fraction_in_region=0
do i_source=1,n_source_index
if (calculate_source(i_source).or.calculate_EMEP_source(i_source)) then
!Loop through the subgrid and find the EMEP grid it is in
do j=1,emission_subgrid_dim(y_dim_index,i_source)
do i=1,emission_subgrid_dim(x_dim_index,i_source)
ii=crossreference_emission_to_emep_subgrid(i,j,x_dim_index,i_source)
jj=crossreference_emission_to_emep_subgrid(i,j,y_dim_index,i_source)
if (ii.ge.1.and.ii.le.dim_length_nc(x_dim_nc_index).and.jj.ge.1.and.jj.le.dim_length_nc(y_dim_nc_index)) then
if (use_subgrid_region(i,j,i_source)) then
EMEP_grid_fraction_in_region(ii,jj,i_source,1)=EMEP_grid_fraction_in_region(ii,jj,i_source,1)+1.
endif
!write(*,*) ii,jj,dim_length_nc(x_dim_nc_index),dim_length_nc(y_dim_nc_index),i_source
count_EMEP_grid_fraction_in_region(ii,jj,i_source)=count_EMEP_grid_fraction_in_region(ii,jj,i_source)+1.
!write(*,*) EMEP_grid_fraction_in_region(ii,jj,i_source,1),count_EMEP_grid_fraction_in_region(ii,jj,i_source)
endif
enddo
enddo
!write(*,*) i_source,sum(count_EMEP_grid_fraction_in_region(:,:,i_source)),sum(EMEP_grid_fraction_in_region(:,:,i_source,1))
where (count_EMEP_grid_fraction_in_region(:,:,i_source).gt.0) EMEP_grid_fraction_in_region(:,:,i_source,1)=EMEP_grid_fraction_in_region(:,:,i_source,1)/count_EMEP_grid_fraction_in_region(:,:,i_source)
!loop through all the EMEP grids and fill in the local fraction arrays
do ii=1,dim_length_nc(x_dim_nc_index)
do jj=1,dim_length_nc(y_dim_nc_index)
do iii=1,dim_length_nc(xdist_dim_nc_index)
do jjj=1,dim_length_nc(ydist_dim_nc_index)
iiii=iii-xdist_centre_nc
jjjj=jjj-ydist_centre_nc
ii_nc=ii+iiii
jj_nc=jj+jjjj
if (ii_nc.ge.1.and.ii_nc.le.dim_length_nc(x_dim_nc_index).and.jj_nc.ge.1.and.jj_nc.le.dim_length_nc(y_dim_nc_index)) then
lf_EMEP_grid_fraction_in_region(iii,jjj,ii,jj,:,1)=EMEP_grid_fraction_in_region(ii_nc,jj_nc,:,1)
endif
enddo
enddo
enddo
enddo
!Allocate the additional subgrid fraction based on the EMEP fraction
if (EMEP_additional_grid_interpolation_size.gt.0) then
!Use the starting position of the read in EMEP file to initialise the starting point
ii_start=mod(dim_start_EMEP_nc(x_dim_nc_index)-1,local_fraction_grid_size(2))
jj_start=mod(dim_start_EMEP_nc(y_dim_nc_index)-1,local_fraction_grid_size(2))
do ii=1,dim_length_nc(x_dim_nc_index)
do jj=1,dim_length_nc(y_dim_nc_index)
!Bottom left corner of the additional grid associatedd with thatEMEP grid
iii=int((ii-1+ii_start)/local_fraction_grid_size(2))*local_fraction_grid_size(2)+1-ii_start
jjj=int((jj-1+jj_start)/local_fraction_grid_size(2))*local_fraction_grid_size(2)+1-jj_start
if (iii.ge.1.and.iii.le.dim_length_nc(x_dim_nc_index).and.jjj.ge.1.and.jjj.le.dim_length_nc(y_dim_nc_index)) then
!Add the region fractions from EMEP to the additional
do iiii=iii,iii-1+local_fraction_grid_size(2)
do jjjj=jjj,jjj-1+local_fraction_grid_size(2)
if (iiii.ge.1.and.iiii.le.dim_length_nc(x_dim_nc_index).and.jjjj.ge.1.and.jjjj.le.dim_length_nc(y_dim_nc_index)) then
EMEP_grid_fraction_in_region(ii,jj,i_source,2)=EMEP_grid_fraction_in_region(ii,jj,i_source,2)+EMEP_grid_fraction_in_region(iiii,jjjj,i_source,1)
endif
enddo
enddo
endif
!Normalize. Does not account for edges but should not be a problem. For safety limit it to 1.
EMEP_grid_fraction_in_region(ii,jj,i_source,2)=min(1.,EMEP_grid_fraction_in_region(ii,jj,i_source,2)*local_fraction_grid_size(1)**2/local_fraction_grid_size(2)**2)
!write(*,'(8i,2f8.2)') ii,jj,iii,jjj,iii-1+local_fraction_grid_size(2),jjj-1+local_fraction_grid_size(2),ii_start,jj_start,EMEP_grid_fraction_in_region(ii,jj,i_source,1),EMEP_grid_fraction_in_region(ii,jj,i_source,2)
enddo
enddo
!When using additional grid
!loop through all the additional EMEP grids and fill in the local fraction arrays
do ii=1,dim_length_nc(x_dim_nc_index)
do jj=1,dim_length_nc(y_dim_nc_index)
do iii=1,dim_length_nc(xdist_dim_nc_index)
do jjj=1,dim_length_nc(ydist_dim_nc_index)
iiii=iii-xdist_centre_nc
jjjj=jjj-ydist_centre_nc
ii_nc=ii+iiii*local_fraction_grid_size(2)
jj_nc=jj+jjjj*local_fraction_grid_size(2)
if (ii_nc.ge.1.and.ii_nc.le.dim_length_nc(x_dim_nc_index).and.jj_nc.ge.1.and.jj_nc.le.dim_length_nc(y_dim_nc_index)) then
lf_EMEP_grid_fraction_in_region(iii,jjj,ii,jj,:,2)=EMEP_grid_fraction_in_region(ii_nc,jj_nc,:,2)
endif
enddo
enddo
enddo
enddo
endif
endif
enddo
if (allocated(count_EMEP_grid_fraction_in_region)) deallocate (count_EMEP_grid_fraction_in_region)
!Fill in missing sources for EMEP
!Choose a source with highest resolution, the last one in this case as they all should be the same
chosen_source=0
do i_source=1,n_source_index
if (calculate_source(i_source)) then
chosen_source=i_source
endif
enddo
!write(*,*) 'Chosen source: ',chosen_source
do i_source=1,n_source_nc_index
if (chosen_source.ne.i_source) then
EMEP_grid_fraction_in_region(:,:,i_source,:)=EMEP_grid_fraction_in_region(:,:,chosen_source,:)
lf_EMEP_grid_fraction_in_region(:,:,:,:,i_source,:)=lf_EMEP_grid_fraction_in_region(:,:,:,:,chosen_source,:)
endif
enddo
!EMEP_grid_fraction_in_region(:,:,allsource_index,:)=EMEP_grid_fraction_in_region(:,:,chosen_source,:)
!lf_EMEP_grid_fraction_in_region(:,:,:,:,allsource_index,:)=lf_EMEP_grid_fraction_in_region(:,:,:,:,chosen_source,:)
!write(*,*) sum(EMEP_grid_fraction_in_region(:,:,traffic_nonexhaust_index,1)),sum(EMEP_grid_fraction_in_region(:,:,traffic_index,1)),sum(EMEP_grid_fraction_in_region(:,:,traffic_exhaust_index,1))
!write(*,*) sum(EMEP_grid_fraction_in_region(:,:,traffic_nonexhaust_index,2)),sum(EMEP_grid_fraction_in_region(:,:,traffic_index,2)),sum(EMEP_grid_fraction_in_region(:,:,traffic_exhaust_index,2))
!write(*,*) sum(lf_EMEP_grid_fraction_in_region(:,:,:,:,traffic_nonexhaust_index,1)),sum(lf_EMEP_grid_fraction_in_region(:,:,:,:,traffic_index,1)),sum(lf_EMEP_grid_fraction_in_region(:,:,:,:,traffic_exhaust_index,1))
!write(*,*) sum(lf_EMEP_grid_fraction_in_region(:,:,:,:,traffic_nonexhaust_index,2)),sum(lf_EMEP_grid_fraction_in_region(:,:,:,:,traffic_index,2)),sum(lf_EMEP_grid_fraction_in_region(:,:,:,:,traffic_exhaust_index,2))
endif
end subroutine uEMEP_assign_region_coverage_to_EMEP