-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathEmissions_mod.f90
executable file
·2268 lines (2042 loc) · 105 KB
/
Emissions_mod.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
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
module Emissions_mod
! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Calls up emission read/set routines
! This routine interfaces the stand-alone emission-file reading routines
! with the 3D model.
!_____________________________________________________________________________
use AirEmis_mod, only : airn
use Biogenics_mod, only: SoilNOx, AnnualNdep
use CheckStop_mod, only: CheckStop,StopAll
use ChemDims_mod, only: NSPEC_SHL, NSPEC_TOT,&
NEMIS_File ! No. emission files
use ChemSpecs_mod, only: NO2, SO2,species,species_adv
use Chemfields_mod, only: xn_adv
use Config_module,only: &
KMAX_MID, KMAX_BND, PT ,dt_advec, &
KCHEMTOP, &
emis_inputlist, & !TESTC
EmisDir, & ! template for emission path
DataDir, & ! template for path
EMIS_OUT, & ! output emissions in ASCII or not
! MONTHLY_GRIDEMIS, & !NML
INERIS_SNAP2 , & ! INERIS/TFMM HDD20 method
MasterProc, USES, & !
SEAFIX_GEA_NEEDED, & ! see below
EURO_SOILNOX_DEPSCALE,&! one or the other
USE_OCEAN_NH3,USE_OCEAN_DMS,FOUND_OCEAN_DMS,&
NPROC, EmisSplit_OUT,USE_uEMEP,uEMEP,SECTORS_NAME,USE_SECTORS_NAME,&
SecEmisOutWanted,MaxNSECTORS,&
AircraftEmis_FLFile,nox_emission_1996_2005File,RoadMapFile,&
AVG_SMI_2005_2010File,NdepFile,&
startdate, Emis_sourceFiles, EmisMask
use Country_mod, only: MAXNLAND,NLAND,Country,IC_NAT,IC_FI,IC_NO,IC_SE
use Country_mod, only: EU28,EUMACC2,IC_DUMMY
use Debug_module, only: DEBUG, & !DEBUG => DEBUG_EMISSIONS, &
DEBUG_EMISTIMEFACS
use EmisDef_mod, only: &
EMIS_FILE & ! Names of species ("sox ",...)
,NCMAX & ! Max. No. countries per grid
,ISNAP_DOM & ! snap index for domestic/resid emis
,ISNAP_TRAF & ! snap index for road-traffic (SNAP7)
,NROAD_FILES & ! No. road dust emis potential files
,ROAD_FILE & ! Names of road dust emission files
,NROADDUST & ! No. road dust components
,QROADDUST_FI & ! fine road dust emissions (PM2.5)
,QROADDUST_CO & ! coarse road dust emis
,ROADDUST_FINE_FRAC & ! fine (PM2.5) fraction of road dust emis
,ROADDUST_CLIMATE_FILE &! TEMPORARY! file for road dust climate factors
,nGridEmisCodes,GridEmisCodes,GridEmis,cdfemis&
,secemis,roaddust_emis_pot,SplitEmisOut,EmisOut&
,SecEmisOut,NSecEmisOutWanted,isec2SecOutWanted&
,nlandcode,landcode&
,road_nlandcode,road_landcode&
,gridrcemis,gridrcroadd,gridrcroadd0&
,O_NH3, O_DMS&
,Emis_4D,N_Emis_4D,Found_Emis_4D & !used for EEMEP
,KEMISTOP&
,MAXFEMISLONLAT,N_femis_lonlat,loc_frac &
,NSECTORS, N_HFAC, N_TFAC, N_SPLIT & ! No. emis sectors, height, time and split classes
,sec2tfac_map, sec2hfac_map, sec2split_map& !generic mapping of indices
,Nneighbors & !used for uemep/loc_frac
,NSECTORS_SNAP, SNAP_sec2tfac_map, SNAP_sec2hfac_map, SNAP_sec2split_map&!SNAP specific mapping
,NSECTORS_GNFR, GNFR_sec2tfac_map, GNFR_sec2hfac_map, GNFR_sec2split_map&!GNFR specific mapping
,NSECTORS_TEST, TEST_sec2tfac_map, TEST_sec2hfac_map, TEST_sec2split_map&!TEST specific mapping
,gnfr2snap,snap2gnfr&
,foundYearlySectorEmissions, foundMonthlySectorEmissions&
,Emis_mask, Emis_mask_allocate, MASK_LIMIT & !old format
, NEmisMask, EmisMaskValues & !new format
,Emis_field, Emis_id, NEmis_id &
,NEmisFile_sources, EmisFiles,NEmis_sources, Emis_source&
, Emis_source_2D, Emis_source_3D,ix3Dmap, NEmis_3Dsources
use EmisGet_mod, only: &
EmisSplit &
,EmisGetCdf &
,EmisGetASCII &
,femis & ! Gets scaling factors -> e_fact
,e_fact & ! scaling factors
,e_fact_lonlat & ! scaling factors
,EmisHeights & ! Generates vertical distrib
,nrcemis, nrcsplit, emisfrac & ! speciation routines and array
,nemis_kprofile, emis_kprofile &! Vertical emissions profile
,iqrc2itot,itot2iqrc & ! maps from split index to total index
,iqrc2iem,iemsplit2itot & ! maps from split index to emission index
,emis_masscorr & ! 1/molwt for most species
,emis_nsplit & ! No. species per emis file
,RoadDustGet &
,roaddust_masscorr & ! 1/200. Hard coded at the moment, needs proper setting in EmisGet_mod...
,femis_latmin,femis_latmax,femis_lonmin,femis_lonmax,femis_lonlat_ic&
,Emis_init_GetCdf, Emis_GetCdf, make_iland_for_time
use GridValues_mod, only: GRIDWIDTH_M & ! size of grid (m)
,xm2 & ! map factor squared
,debug_proc,debug_li,debug_lj &
,xmd,dA,dB,i_fdom,j_fdom,glon,glat
use Io_Nums_mod, only: IO_LOG, IO_DMS, IO_EMIS, IO_TMP
use Io_Progs_mod, only: ios, open_file, datewrite, PrintLog
use MetFields_mod, only: u_xmj, v_xmi, roa, ps, z_bnd, surface_precip,EtaKz ! ps in Pa, roa in kg/m3
use MetFields_mod, only: t2_nwp ! DS_TEST SOILNO - was zero!
use MPI_Groups_mod , only : MPI_BYTE, MPI_DOUBLE_PRECISION, MPI_REAL8, MPI_INTEGER&
,MPI_LOGICAL, MPI_SUM,MPI_COMM_CALC, IERROR
use NetCDF_mod, only: ReadField_CDF,ReadField_CDF_FL,ReadTimeCDF,IsCDFfractionFormat,&
GetCDF_modelgrid,PrintCDF,ReadSectorName,check,&
create_country_emission_file, output_country_emissions
use netcdf
use OwnDataTypes_mod, only: TXTLEN_FILE, TXTLEN_NAME,Emis_id_type, EmisFile_id_type, Emis_sourceFile_id_type
use Par_mod, only: MAXLIMAX,MAXLJMAX, GIMAX,GJMAX, IRUNBEG,JRUNBEG,&
me,limax,ljmax, MSG_READ1,MSG_READ7&
,gi0,gj0,li0,li1,lj0,lj1
use PhysicalConstants_mod,only: GRAV, AVOG, ATWAIR
use PointSource_mod, only: readstacks !MKPS
use ZchemData_mod,only: rcemis ! ESX
use SmallUtils_mod, only: find_index, key2str, trims
use TimeDate_mod, only: nydays, nmdays, date, current_date,&! No. days per
tdif_secs,timestamp,make_timestamp,daynumber,day_of_week ! year, date-type, weekday
use TimeDate_ExtraUtil_mod,only :nctime2date, date2string, to_idate,date_is_reached
use Timefactors_mod, only: &
NewDayFactors & ! subroutines
,DegreeDayFactors & ! degree-days used for SNAP-2
,Gridded_SNAP2_Factors, gridfac_HDD &
,fac_min,timefactors & ! subroutine
,fac_ehh24x7 ,fac_emm, fac_cemm, fac_edd, timefac & ! time-factors
,Read_monthly_emis_grid_fac &
,GridTfac &!array with monthly gridded time factors
,yearly_normalize !renormalize timefactors after reset
implicit none
private
! subroutines:
public :: Init_masks
public :: Init_Emissions ! defines emission setup and formats
public :: EmisUpdate ! update emission
public :: Emissions ! Main emissions module
public :: newmonth
public :: EmisSet ! Sets emission rates every hour/time-step
public :: EmisWriteOut ! Outputs emissions in ascii
! The main code does not need to know about the following
private :: expandcclist ! expands e.g. EU28, EUMACC2
private :: consistency_check ! Safety-checks
logical, save, private :: first_dms_read
! and for budgets (not yet used - not changed dimension)
real, public, save, dimension(NSPEC_SHL+1:NSPEC_TOT) :: totemadd
integer, private, save :: iemCO ! index of CO emissions, for debug
real ::TimesInDays(120),mpi_out
integer ::NTime_Read=-1,found
logical :: fileExists ! to test emission files
integer :: nin, nex ! include, exclude numbers for emis_inputlist
character(len=*), parameter :: sub='Emissions:'
logical, save :: mappingGNFR2SNAP = .false.
contains
!***********************************************************************
!new (Nov 2018) emission setup and formats
subroutine Init_Emissions
!loop through all sources to set parameters for each source.
!one source is defined as one two dimensional map giving emissions.
!one file can have many sources.
!each source has parameters. Parameters are constant through the run.
!There are several ways to set parameters. In increasing priority order:
!1) default
!2) read from global attributes in file
!3) read from variable attributes in file
!4) file parameter set in namelist (config_emep.nml)
!5) source parameter set in namelist (config_emep.nml)
!multiplicative factor on top of each other, from the most general:
!1) from femis.dat (can be switched off for specific files, but not for specific sources)
!2) global file attribute (can be overwritten by config)
!3) individual sources (can be overwritten by config)
!e_fact (femis) is applied through the source%factor in Emis_init_GetCdf
!Note on units "as SO2":
!if the sector is defined, units are defined "as SO2" and SO4 must include the
!factor 0.6666 (=mw(SO2)/mw(SO4)) to be right.
!If sector is defined as zero, no additional factor is required.
!note on apply_femis and some other parameters of type "logical":
!these cannot be defined on file and then possibly overwritten by config,
!because it is not possible to distiguish between default and set values from config.
!apply_femis must be either default (true) or set to false by config
!include_in_local_fractions must be either default (true) or set to false by config
integer, parameter ::maxnames=100
character(len=TXTLEN_FILE) :: fname, filename, names_in(maxnames)
integer :: i, ii, n, nn, ix, nemis_old, isource
integer :: isec, iland, iem, iqrc, itot, f
integer :: startsource(size(Emis_sourceFiles)), endsource(size(Emis_sourceFiles))
type(Emis_id_type):: Emis_id_undefined !to get undefined source values
type(EmisFile_id_type):: Emisfile_undefined !to get undefined file values
!Note must distinguish between default and not undefined values, to know when config has defined a value, even if it is the default.
type(Emis_id_type):: Emis_sources_defaults !set values when not specified otherwise
type(EmisFile_id_type):: Emisfile_defaults !set values when not specified otherwise
integer :: EmisFilesMap(0:size(Emis_sourceFiles)) !index of EmisFile given index of EmisFile_sources
integer :: max_levels3D
character(len=*),parameter :: dtxt='Ini_Em:'
logical :: found
logical, save :: first_call=.true., dbg
if ( first_call ) then
dbg = DEBUG%EMISSIONS .and. MasterProc
first_call = .false.
end if
!1) define lowest level default values
Emis_sources_defaults%units = 'mg/m2/h'
Emis_sources_defaults%country_ISO = 'N/A'
Emis_sources_defaults%sector = 0
Emis_sources_defaults%factor = 1.0
Emis_sources_defaults%include_in_local_fractions = .true.
Emis_sources_defaults%apply_femis = .true.
Emis_sources_defaults%injection_k = KMAX_MID
Emis_sources_defaults%is3D = .false.
Emis_sources_defaults%istart = 1
Emis_sources_defaults%jstart = 1
Emis_sources_defaults%kstart = 1
Emis_sources_defaults%kend = 1
Emis_sources_defaults%reversek = .false.
Emis_sources_defaults%species_ix = -1
Emis_source = Emis_sources_defaults!set all initial values to default
Emisfile_defaults%periodicity = 'time'
Emisfile_defaults%projection = 'lon lat'
Emisfile_defaults%factor = 1.0
!2) read from global attributes in file
!3) read from variable attributes in file
!Emis_sourceFiles is read from config and then not modified
!EmisFiles collect all valid data and sources
NEmis_sources = 0 !total number of valid emis (also 3D) variables across all files
NEmis_3Dsources = 0 !total number of valid 3D emis variables across all files
max_levels3D = 1
NEmisFile_sources = 0 !number of valid emission files
EmisFilesMap = 0 !index of EmisFile given index of EmisFile_sources
do n = 1, size(Emis_sourceFiles)
nemis_old = NEmis_sources
if(Emis_sourceFiles(n)%filename/='NOTSET')then
!Read all variables and set parameters as needed
!find which variables names are defined in config
names_in='NOTSET'
i=0
do nn = 1, size(Emis_sourceFiles(n)%source)
if(Emis_sourceFiles(n)%source(nn)%varname/='NOTSET')then
i= i + 1
names_in(i)=trim(Emis_sourceFiles(n)%source(nn)%varname)
endif
enddo
if(MasterProc)write(*,*)dtxt//"Initializing Emissions from ",&
trim(Emis_sourceFiles(n)%filename)
call Emis_init_GetCdf(Emis_sourceFiles(n), EmisFiles(NEmisFile_sources+1), names_in, i)
endif
startsource(n) = nemis_old + 1
endsource(n) = NEmis_sources
if(NEmis_sources>nemis_old)then
EmisFilesMap(NEmisFile_sources) = n
EmisFiles(NEmisFile_sources)%Nsources = NEmis_sources - nemis_old
EmisFiles(NEmisFile_sources)%source_start = nemis_old + 1
EmisFiles(NEmisFile_sources)%source_end = NEmis_sources
endif
enddo
allocate(Emis_source_2D(LIMAX,LJMAX,NEmis_sources))
!4)overwrite parameters with settings from config_emep.nml if they are set
!first overwrite the global attributes: projection and periodicity
do i = 1, size(Emis_sourceFiles)
n = EmisFilesMap(i)
if(n>0)then
if(Emis_sourceFiles(n)%periodicity /= Emisfile_undefined%periodicity) EmisFiles(i)%periodicity = Emis_sourceFiles(n)%periodicity
if(Emis_sourceFiles(n)%projection /= Emisfile_undefined%projection) EmisFiles(i)%projection = Emis_sourceFiles(n)%projection
if(Emis_sourceFiles(n)%grid_resolution /= Emisfile_undefined%grid_resolution) EmisFiles(i)%grid_resolution = Emis_sourceFiles(n)%grid_resolution
if(Emis_sourceFiles(n)%projection /= 'native')then
call CheckStop(EmisFiles(i)%grid_resolution <=1.0E-5,'Grid_resolution must be defined for '//trim(Emis_sourceFiles(n)%filename))
endif
if(Emis_sourceFiles(n)%mask_ID /= Emisfile_undefined%mask_ID) EmisFiles(i)%mask_ID = Emis_sourceFiles(n)%mask_ID
if(Emis_sourceFiles(n)%mask_ID_reverse /= Emisfile_undefined%mask_ID_reverse) EmisFiles(i)%mask_ID_reverse = Emis_sourceFiles(n)%mask_ID_reverse
if(Emis_sourceFiles(n)%species /= Emisfile_undefined%species) EmisFiles(i)%species = Emis_sourceFiles(n)%species
if(Emis_sourceFiles(n)%units /= Emisfile_undefined%units) EmisFiles(i)%units = Emis_sourceFiles(n)%units
if(Emis_sourceFiles(n)%country_ISO /= Emisfile_undefined%country_ISO) EmisFiles(i)%country_ISO = Emis_sourceFiles(n)%country_ISO
if(Emis_sourceFiles(n)%sector /= Emisfile_undefined%sector) EmisFiles(i)%sector = Emis_sourceFiles(n)%sector
endif
enddo
!then overwrite the variable attributes
do i = 1, NEmisFile_sources !loop over files
n = EmisFilesMap(i)
found = .false.
do ii = EmisFiles(i)%source_start, EmisFiles(i)%source_end !loop over sources found in the netcdf file
!set source default = file parameter if they are set. Defines default for sources in this file, can be also be redefined for individual sources
if(Emis_sourceFiles(n)%species /= Emisfile_undefined%species) Emis_source(ii)%species = Emis_sourceFiles(n)%species
if(Emis_sourceFiles(n)%units /= Emisfile_undefined%units) Emis_source(ii)%units = Emis_sourceFiles(n)%units
if(Emis_sourceFiles(n)%country_ISO /= Emisfile_undefined%country_ISO) Emis_source(ii)%country_ISO = Emis_sourceFiles(n)%country_ISO
if(Emis_sourceFiles(n)%sector /= Emisfile_undefined%sector) Emis_source(ii)%sector = Emis_sourceFiles(n)%sector
Emis_source(ii)%apply_femis = Emis_sourceFiles(n)%apply_femis
Emis_source(ii)%include_in_local_fractions = Emis_sourceFiles(n)%include_in_local_fractions
Emis_source(ii)%periodicity = Emis_sourceFiles(n)%periodicity
Emis_source(ii)%mask_ID = Emis_sourceFiles(n)%mask_ID
Emis_source(ii)%mask_ID_reverse = Emis_sourceFiles(n)%mask_ID_reverse
isource = Emis_source(ii)%ix_in
if(dbg) write(*,'(a,4i4,1x,a20)') 'wrting config attribute on '//trim(Emis_source(ii)%varname),n, ii, isource, NEmis_sources
if(isource>0)then
!source defined in config file
if(trim(Emis_sourceFiles(n)%source(isource)%varname)/=trim(Emis_source(ii)%varname))write(*,*)isource,'ERROR',trim(Emis_sourceFiles(n)%source(isource)%varname),' ',trim(Emis_source(ii)%varname),ii
found = .true.
if(Emis_sourceFiles(n)%source(isource)%species /= Emis_id_undefined%species) Emis_source(ii)%species = Emis_sourceFiles(n)%source(isource)%species
if(Emis_sourceFiles(n)%source(isource)%units /= Emis_id_undefined%units) Emis_source(ii)%units = Emis_sourceFiles(n)%source(isource)%units
if(Emis_sourceFiles(n)%source(isource)%sector /= Emis_id_undefined%sector) Emis_source(ii)%sector = Emis_sourceFiles(n)%source(isource)%sector
if(Emis_sourceFiles(n)%source(isource)%factor /= Emis_id_undefined%factor) Emis_source(ii)%factor = Emis_sourceFiles(n)%source(isource)%factor
if(Emis_sourceFiles(n)%source(isource)%country_ISO /= Emis_id_undefined%country_ISO) Emis_source(ii)%country_ISO = Emis_sourceFiles(n)%source(isource)%country_ISO
if(.not. Emis_sourceFiles(n)%source(isource)%include_in_local_fractions) Emis_source(ii)%include_in_local_fractions = .false.
if(.not. Emis_sourceFiles(n)%source(isource)%apply_femis) Emis_source(ii)%apply_femis = Emis_sourceFiles(n)%source(isource)%apply_femis
if(Emis_sourceFiles(n)%source(isource)%mask_ID /= Emis_id_undefined%mask_ID) Emis_source(ii)%mask_ID = Emis_sourceFiles(n)%source(isource)%mask_ID
if(Emis_sourceFiles(n)%source(isource)%mask_ID_reverse /= Emis_id_undefined%mask_ID_reverse) Emis_source(ii)%mask_ID_reverse = Emis_sourceFiles(n)%source(isource)%mask_ID_reverse
if(Emis_sourceFiles(n)%source(isource)%is3D .neqv. Emis_id_undefined%is3D) Emis_source(ii)%is3D = Emis_sourceFiles(n)%source(isource)%is3D
if(Emis_sourceFiles(n)%source(isource)%istart /= Emis_id_undefined%istart) Emis_source(ii)%istart = Emis_sourceFiles(n)%source(isource)%istart
if(Emis_sourceFiles(n)%source(isource)%jstart /= Emis_id_undefined%jstart) Emis_source(ii)%jstart = Emis_sourceFiles(n)%source(isource)%jstart
if(Emis_sourceFiles(n)%source(isource)%kstart /= Emis_id_undefined%kstart) Emis_source(ii)%kstart = Emis_sourceFiles(n)%source(isource)%kstart
if(Emis_sourceFiles(n)%source(isource)%kend /= Emis_id_undefined%kend) Emis_source(ii)%kend = Emis_sourceFiles(n)%source(isource)%kend
if(Emis_sourceFiles(n)%source(isource)%reversek .neqv. Emis_id_undefined%reversek) Emis_source(ii)%reversek = Emis_sourceFiles(n)%source(isource)%reversek
if(Emis_sourceFiles(n)%source(isource)%injection_k /= Emis_id_undefined%injection_k) Emis_source(ii)%injection_k = Emis_sourceFiles(n)%source(isource)%injection_k
endif
ix = find_index(trim(Emis_source(ii)%country_ISO) ,Country(:)%code, first_only=.true.)
if(ix<0)then
if(me==0)write(*,*)dtxt//'WARNING: country '//trim(Emis_source(ii)%country_ISO)//' not defined. '
else
Emis_source(ii)%country_ix = ix
if(dbg)write(*,*)dtxt//'country found '//trim(Emis_source(ii)%country_ISO), ix, NEmis_sources
endif
!find if it is defined as an individual species
ix = find_index(Emis_source(ii)%species, species(:)%name )
if(ix>0)then
Emis_source(ii)%species_ix = ix
if(dbg)write(*,'(a,i4,a)')dtxt//' species found '// &
trim(Emis_source(ii)%country_ISO), ix, ' '//trim(species(ix)%name)
if(Emis_source(ii)%include_in_local_fractions .and. USE_uEMEP )then
if(me==0)write(*,*)"WARNING: local fractions will not include single species "//Emis_source(ii)%species
endif
else ! ix<=0
if(dbg)write(*,'(a,i4,a)')dtxt//' species not found'// &
trim(Emis_source(ii)%country_ISO),ix,trim(Emis_source(ii)%species)
endif
max_levels3D=max(max_levels3D, Emis_source(ii)%kend - Emis_source(ii)%kstart + 1)
if(MasterProc .and. dbg)write(*,*)dtxt//"Final emission source parameters ",Emis_source(ii)
enddo
! if(.not. found .and. me==0)write(*,*)dtxt//'WARNING: did not find some of the emission sources defined in config in '&
! //trim(Emis_sourceFiles(n)%filename)
enddo
!include reduction factors
do n = 1, NEmis_sources
if(Emis_source(n)%apply_femis)then
isec = Emis_source(n)%sector
if(Emis_source(n)%sector>0 .and. Emis_source(n)%sector<=NSECTORS)then
iland = Emis_source(n)%country_ix
if(iland<0)iland=IC_DUMMY
isec = Emis_source(n)%sector
iem = find_index(Emis_source(n)%species,EMIS_FILE(:))
if(iem >0 )then
!apply femis
Emis_source(n)%factor = Emis_source(n)%factor * e_fact(isec,iland,iem)
else
!see if the species belongs to any of the splitted species
iqrc = 0
do iem = 1,NEMIS_FILE
do f = 1,emis_nsplit(iem)
iqrc = iqrc + 1
itot = iqrc2itot(iqrc)
if(trim(species(itot)%name)==trim(Emis_source(n)%species))then
if(dbg) write(*,'(a,5i4,a,f12.3)')&
trim(Emis_source(n)%species)//' included in '//trim(EMIS_FILE(iem)), n, &
isec, iland, Emis_source(n)%country_ix, iem, &
trim( species(itot)%name ), e_fact(isec,iland,iem)
Emis_source(n)%factor = Emis_source(n)%factor * e_fact(isec,iland,iem)
go to 888
endif
enddo
enddo ! iem
888 continue
endif
endif
endif
!define mask indices
if(Emis_source(n)%mask_ID /= Emis_id_undefined%mask_ID)then
ix = find_index(Emis_source(n)%mask_ID,EmisMask(:)%ID)
if(ix > 0)then
Emis_source(n)%mask_ix = ix
else
Emis_source(n)%mask_ix = Emis_id_undefined%mask_ix !in case somebody tries to set in config
if(MasterProc)write(*,*)'WARNING mask not defined ',trim(Emis_source(n)%mask_ID)
endif
endif
if(Emis_source(n)%mask_ID_reverse /= Emis_id_undefined%mask_ID_reverse)then
ix = find_index(Emis_source(n)%mask_ID_reverse,EmisMask(:)%ID)
if(ix > 0)then
Emis_source(n)%mask_reverse_ix = ix
else
Emis_source(n)%mask_reverse_ix = Emis_id_undefined%mask_reverse_ix !in case somebody tries to set in config
if(MasterProc)write(*,*)'WARNING mask not defined ',trim(Emis_source(n)%mask_ID_reverse)
endif
endif
enddo ! n = 1, NEmis_sources
!find and define the 3D emissions
ix=0
do n=1, NEmis_sources
if(Emis_source(n)%is3D)Nemis_3Dsources = Nemis_3Dsources + 1
ix3Dmap(n)=ix
enddo
if(Nemis_3Dsources>0)then
if(me==0)write(*,*)'found ',Nemis_3Dsources,' 3D sources'
allocate(Emis_source_3D(LIMAX,LJMAX,max_levels3D,Nemis_3Dsources))
endif
end subroutine Init_Emissions
subroutine Init_masks()
!sets and define the masks
real :: mask_cdf(LIMAX,LJMAX), xsum
logical :: found
integer :: iEmisMask = 0
integer :: i,ii,jj,ic
iEmisMask = 0
!1) find number of valid masks defined
do i = 1, size(EmisMask)
if(EmisMask(i)%filename /= 'NOTSET' .and. EmisMask(i)%cdfname /= 'NOTSET'&
.and. EmisMask(i)%ID /= 'NOTSET' ) then
iEmisMask = iEmisMask+1 !assumes the fields are defined, without checking
endif
enddo
NEmisMask = iEmisMask
allocate(EmisMaskValues(LIMAX,LJMAX,NEmisMask))
!now set the values for the actual masks
iEmisMask = 0
do i = 1, size(EmisMask)
if(EmisMask(i)%filename /= 'NOTSET' .and. EmisMask(i)%cdfname /= 'NOTSET'&
.and. EmisMask(i)%ID /= 'NOTSET' ) then
iEmisMask = iEmisMask+1
call ReadField_CDF(trim(EmisMask(i)%filename),trim(EmisMask(i)%cdfname),mask_cdf,1,&
interpol='conservative', needed=.false.,found = found, UnDef=-1.0E10, debug_flag=.false.)
if(found)then
!set mask value
ic = 0
xsum = 0.0
do jj = 1, LJMAX
do ii = 1, LIMAX
xsum = xsum + mask_cdf(ii,jj)
if(mask_cdf(ii,jj)>EmisMask(i)%threshold)then
EmisMaskValues(ii,jj,iEmisMask) = 0.0 !remove everything
ic = ic + 1
else
EmisMaskValues(ii,jj,iEmisMask) = 1.0 !keep everything
endif
end do
end do
if(MasterProc)write(*,*)'defined mask '//trim(EmisMask(i)%ID)//' based on '//trim(EmisMask(i)%cdfname)
if(ic>0)write(*,*)me,' masked ',ic,' cells'
else
call StopAll("Mask variable not found: "//trim(EmisMask(i)%cdfname))
EmisMask(i)%ID = 'NOTSET'!cannot be used anymore
endif
!to keep some compatibility with old format we also set old format masks
if(any(emis_inputlist(:)%use_mask))then
if(.not.allocated(Emis_mask))then
allocate(Emis_mask(LIMAX,LJMAX))
Emis_mask = .false.
endif
if(MasterProc)write(*,*)'defining mask for old formats ',trim(EmisMask(i)%cdfname)
do jj = 1, LJMAX
do ii = 1, LIMAX
if(EmisMaskValues(ii,jj,iEmisMask)<0.5)Emis_mask(ii,jj) = .true.
enddo
enddo
endif
endif
enddo
end subroutine Init_masks
!***********************************************************************
subroutine EmisUpdate
!Update emission arrays, and read new sets as required
integer :: n, i, j, ix, is, date_limit(5)
type(date) :: coming_date
real :: fac
TYPE(timestamp) :: ts1,ts2
ts1=make_timestamp(current_date)
coming_date = current_date
coming_date%seconds = coming_date%seconds + 1800!NB: end_of_validity_date is at end of period, for example 1-1-2018 for December 2017
!loop over all sources and see which one need to be reread from files
do n = 1, NEmisFile_sources
if(date_is_reached(to_idate(EmisFiles(n)%end_of_validity_date,5 )))then
if(me==0)write(*,*)'Emis: update date is reached ',EmisFiles(n)%end_of_validity_date,EmisFiles(n)%periodicity
!values are no more valid, fetch new one
do is = EmisFiles(n)%source_start,EmisFiles(n)%source_end
if(Emis_source(is)%is3D)then
ix = ix3Dmap(is)
Emis_source_3D(1:,1:,1:,ix)=0.0
call Emis_GetCdf(EmisFiles(n),Emis_source(is),Emis_source_3D(1,1,1,ix),coming_date)
else
if(me==0)write(*,*)is,' getemis '//trim(Emis_source(is)%units)//' '//trim(Emis_source(is)%varname)
Emis_source_2D(1:,1:,is)=0.0
call Emis_GetCdf(EmisFiles(n),Emis_source(is),Emis_source_2D(1,1,is),coming_date)
endif
!reduction factors
fac = EmisFiles(n)%factor
fac = fac* Emis_source(is)%factor
!unit and factor conversions
!convert into kg/m2/s
if(Emis_source(is)%units == 'kg/s' .or. Emis_source(is)%units == 'kg/m2/s')then
fac = fac
else if(Emis_source(is)%units == 'g/s' .or. Emis_source(is)%units == 'g/m2/s')then
fac = fac /(1000.0)
else if(Emis_source(is)%units == 'mg/s' .or. Emis_source(is)%units == 'mg/m2/s')then
fac = fac /(1000.0)
else
!depends on periodicity
if(EmisFiles(n)%periodicity == 'yearly')then
if(Emis_source(is)%units == 'tonnes/m2' .or. Emis_source(is)%units == 'tonnes/m2/year'&
.or. Emis_source(is)%units == 'tonnes' .or. Emis_source(is)%units == 'tonnes/year')then
fac = fac * 1000/(3600*24*nydays)
else if(Emis_source(is)%units == 'kg/m2' .or. Emis_source(is)%units == 'kg/m2/year'&
.or. Emis_source(is)%units == 'kg' .or. Emis_source(is)%units == 'kg/year')then
fac = fac /(3600*24*nydays)
else if(Emis_source(is)%units == 'g/m2' .or. Emis_source(is)%units == 'g/m2/year'&
.or. Emis_source(is)%units == 'g' .or. Emis_source(is)%units == 'g/year')then
fac = fac /(1000.0*3600*24*nydays)
else if(Emis_source(is)%units == 'mg/m2' .or. Emis_source(is)%units == 'mg/m2/year'&
.or. Emis_source(is)%units == 'mg' .or. Emis_source(is)%units == 'mg/year')then
fac = fac /(1.0e6*3600*24*nydays)
else
call StopAll("B Unit for emissions not recognized: "//trim(Emis_source(is)%units))
endif
else if(EmisFiles(n)%periodicity == 'monthly')then
if(Emis_source(is)%units == 'tonnes/m2' .or. Emis_source(is)%units == 'tonnes/m2/month'&
.or. Emis_source(is)%units == 'tonnes' .or. Emis_source(is)%units == 'tonnes/month')then
fac = fac *1000/(3600*24*nmdays(coming_date%month))
else if(Emis_source(is)%units == 'kg/m2' .or. Emis_source(is)%units == 'kg/m2/month'&
.or. Emis_source(is)%units == 'kg' .or. Emis_source(is)%units == 'kg/month')then
fac = fac /(3600*24*nmdays(coming_date%month))
else if(Emis_source(is)%units == 'g/m2' .or. Emis_source(is)%units == 'g/m2/month'&
.or. Emis_source(is)%units == 'g' .or. Emis_source(is)%units == 'g/month')then
fac = fac /(1000*3600*24*nmdays(coming_date%month))
else if(Emis_source(is)%units == 'mg/m2' .or. Emis_source(is)%units == 'mg/m2/month'&
.or. Emis_source(is)%units == 'mg' .or. Emis_source(is)%units == 'mg/month')then
fac = fac /(1.0e6*3600*24*nmdays(coming_date%month))
else
call StopAll("C Unit for emissions not recognized: "//trim(Emis_source(is)%units))
endif
else
!assume hourly
if(Emis_source(is)%units == 'mg/m2' .or. Emis_source(is)%units == 'mg/m2/h')then
!convert into kg/m2/s
fac = fac /(1.0e6*3600.0)
else if(Emis_source(is)%units == 'g/m2' .or. Emis_source(is)%units == 'g/m2/h')then
fac = fac /(1000.0*3600.0)
if(EmisFiles(n)%periodicity /= 'hourly' .and. Emis_source(is)%units == 'g/m2')then
call StopAll("Emis_source unit g/m2 only implemented for hourly, monthly or yearly. Found "//trim(EmisFiles(n)%periodicity))
endif
else if(Emis_source(is)%units == 'g/s')then
fac = fac /(1000.0)
else if(Emis_source(is)%units == 'kg/s')then
fac = fac
else if(Emis_source(is)%units == 'tonnes/s')then
fac = fac * 1000.0
else
call StopAll("Emis_source unit not implemented. Found "//trim(Emis_source(is)%units)//' '//trim(EmisFiles(n)%periodicity))
endif
!Note: easy to implement more unit choices. Just add "if" cases here
endif
endif
if(Emis_source(is)%units == 'tonnes' .or. Emis_source(is)%units == 'tonnes/s' &
.or.Emis_source(is)%units == 'tonnes/month' .or. Emis_source(is)%units == 'tonnes/year' &
.or. Emis_source(is)%units == 'kg' .or. Emis_source(is)%units == 'kg/s' &
.or. Emis_source(is)%units == 'kg/month' .or. Emis_source(is)%units == 'kg/year' &
.or. Emis_source(is)%units == 'g' .or. Emis_source(is)%units == 'g/s' &
.or. Emis_source(is)%units == 'g/month' .or. Emis_source(is)%units == 'g/year' &
.or. Emis_source(is)%units == 'mg' .or. Emis_source(is)%units == 'mg/s' &
.or. Emis_source(is)%units == 'mg/month' .or. Emis_source(is)%units == 'mg/year' &
.or. Emis_source(is)%units == 'g/h' .or. Emis_source(is)%units == 'mg/h')then
!divide by gridarea
fac = fac / (GRIDWIDTH_M * GRIDWIDTH_M)
do j = 1,ljmax
do i = 1,limax
Emis_source_2D(i,j,is) = Emis_source_2D(i,j,is) * xm2(i,j)
enddo
enddo
endif
do j = 1,ljmax
do i = 1,limax
Emis_source_2D(i,j,is) = Emis_source_2D(i,j,is) * fac
enddo
enddo
!now Emis_source_2D should be in kg/m2/s
!apply masks
!could easily be better CPU optimized if necessary by putting all factors in same i,j loop
if(Emis_source(is)%mask_ix>0)then
do j = 1,ljmax
do i = 1,limax
Emis_source_2D(i,j,is) = Emis_source_2D(i,j,is) * EmisMaskValues(i,j,Emis_source(is)%mask_ix)
enddo
enddo
endif
if(Emis_source(is)%mask_reverse_ix>0)then
do j = 1,ljmax
do i = 1,limax
Emis_source_2D(i,j,is) = Emis_source_2D(i,j,is) * (1.0-EmisMaskValues(i,j,Emis_source(is)%mask_ix))!take complementary
enddo
enddo
endif
enddo
!update date of valitdity
if(EmisFiles(n)%periodicity == 'yearly')then
!assumes only one record to read
EmisFiles(n)%end_of_validity_date = date(current_date%year,1,1,0,0)
EmisFiles(n)%end_of_validity_date%year = EmisFiles(n)%end_of_validity_date%year + 1
else if(EmisFiles(n)%periodicity == 'monthly')then
!assumes 12 records, one for each month
EmisFiles(n)%end_of_validity_date = date(current_date%year,1,1,0,0)
EmisFiles(n)%end_of_validity_date%month = current_date%month + 1
if(EmisFiles(n)%end_of_validity_date%month>12)then
EmisFiles(n)%end_of_validity_date = date(current_date%year+1,1,1,0,0)
endif
else
!the correct times must be written in the file and updated in Emis_GetCdf
endif
endif
enddo
end subroutine EmisUpdate
!***********************************************************************
subroutine Emissions(year)
! Initialize emission variables, and read yearly emissions
integer, intent(in) :: year ! Year ( 4-digit)
!-- local variables
integer :: i, j ! Loop variables
real :: tonne_to_kgm2s ! Converts tonnes/grid to kg/m2/s
real :: ccsum ! Sum of emissions for one country
! arrays for whole EMEP area:
! additional arrays on host only for landcode, nlandcode
! BIG arrays ... will be needed only on me=0. Make allocatable
! to reduce static memory requirements.
real, allocatable, dimension(:,:,:) :: globroad_dust_pot ! Road dust emission potentials
integer, allocatable, dimension(:,:) :: road_globnland
integer, allocatable, dimension(:,:,:) :: road_globland
real, allocatable, dimension(:,:) :: RoadDustEmis_climate_factor ! Climatic factor for scaling road dust emissions (in TNO model based on yearly average soil water)
integer :: err1, err2, err3, err4, err7, err8, err9 ! Error messages
integer :: fic ,insec,inland,iemis,iemislist
integer :: iic,ic,n ! country codes
integer :: isec ! loop variables: emission sectors
integer :: iem ! loop variable over pollutants (1..NEMIS_FILE)
integer :: icc ! loop variables over sources
character(len=300) :: fname ! txt, File name
logical :: fractionformat
! Emission sums (after e_fact adjustments):
real, dimension(NEMIS_FILE) :: emsum ! Sum emis over all countries
real, dimension(NLAND,NEMIS_FILE) :: sumemis, sumemis_local ! Sum of emissions per country
real, dimension(NEMIS_FILE) :: sumEU ! Sum of emissions in EU
! Road dust emission potential sums (just for testing the code, the actual emissions are weather dependent!)
real, dimension(NLAND,NROAD_FILES) :: sumroaddust ! Sum of emission potentials per country
real, dimension(NLAND,NROAD_FILES) :: sumroaddust_local ! Sum of emission potentials per country in subdomain
real :: fractions(LIMAX,LJMAX,NCMAX),SMI(LIMAX,LJMAX),SMI_roadfactor
logical ::SMI_defined=.false.
logical,save :: my_first_call=.true. ! Used for femis call
character(len=TXTLEN_NAME) :: varname, fmt,cdf_sector_name, species_name
integer :: allocerr, i_Emis_4D, sec_ix, ih, id, k, ncFileID, varID, iland
character(len=TXTLEN_FILE) ::fileName
integer :: emis_inputlist_NEMIS_FILE!number of files for each emis_inputlist(i)
real :: buffer(LIMAX,LJMAX)
logical country_owner_map(NLAND,NPROC)
if (MasterProc) write(6,*) "Reading emissions for year", year
ios = 0
country_owner_map = .false.
! initialize emis_inputlist
!>=========================================================
do iemislist = 1, size( emis_inputlist(:)%name )
fname = emis_inputlist(iemislist)%name
if(fname=="NOTSET") cycle
if(MasterProc)&
write(*,*)"Emission source number ", iemislist,"from ",sub//trim(fname)
if(emis_inputlist(iemislist)%type == "sectors".or.&
emis_inputlist(iemislist)%type == "GNFRsectors".or.&
emis_inputlist(iemislist)%type == "SNAPsectors")then ! Expand groups, e.g. EUMACC2
call expandcclist( emis_inputlist(iemislist)%incl , n)
emis_inputlist(iemislist)%Nincl = n
if(MasterProc .and. n>0) then
write(*,*) sub//trim(fname)//" INPUTLIST-INCL", n
write(*,*)'including only countries: ', (trim(emis_inputlist(iemislist)%incl(i))//' ',i=1,n)
endif
call expandcclist( emis_inputlist(iemislist)%excl , n)
emis_inputlist(iemislist)%Nexcl = n
if(MasterProc .and. n>0) then
write(*,*)'excluding countries: ', (trim(emis_inputlist(iemislist)%excl(i))//' ',i=1,n)
endif
end if
if(emis_inputlist(iemislist)%pollName(1)/='NOTSET')then
do iem = 1, NEMIS_FILE
if(all(emis_inputlist(iemislist)%pollName(:)/=trim(EMIS_FILE(iem))))cycle
if(Masterproc)write(*,"(A)")'including pollutant '//trim(EMIS_FILE(iem))//' from '//trim(fname)
enddo
else
!include all pollutants
endif
!replace keywords
22 format(5A)
if(MasterProc)write(*,22)'original emission name ',trim(fname)
fname = key2str(fname,'EmisDir',EmisDir)
fname = key2str(fname,'DataDir',DataDir)
fname = key2str(fname,'YYYY',year)
emis_inputlist(iemislist)%name=trim(fname)
if(MasterProc)write(*,22)'filename redefined as: ',&
trim(emis_inputlist(iemislist)%name)
cdf_sector_name='NOTSET'
fname = key2str(fname,'POLL',EMIS_FILE(1))
call ReadSectorname(fname,cdf_sector_name)
if(trim(cdf_sector_name)/='NOTSET')then
SECTORS_NAME=trim(cdf_sector_name)
if(Masterproc .and. USE_SECTORS_NAME =='NOTSET')&
write(*,*)"Switching sector categories to ",trim(SECTORS_NAME)
if(Masterproc .and. USE_SECTORS_NAME =='NOTSET')&
write(IO_LOG,*)"Switching sector categories to ",trim(SECTORS_NAME)
if(cdf_sector_name == 'GNFR')emis_inputlist(iemislist)%type = "GNFRsectors"
if(cdf_sector_name == 'SNAP')emis_inputlist(iemislist)%type = "SNAPsectors"
end if
if(IsCDFfractionFormat(trim(fname))) emis_inputlist(iemislist)%format='fractions'
if(emis_inputlist(iemislist)%set_mask.or.emis_inputlist(iemislist)%use_mask)Emis_mask_allocate = .true.
if(emis_inputlist(iemislist)%periodicity == 'NOTSET')then
emis_inputlist(iemislist)%periodicity = 'once' !default
if(index(emis_inputlist(iemislist)%name,".nc")>1)then
NTime_Read=-1
call ReadTimeCDF(trim(fname),TimesInDays,NTime_Read)
if(NTime_Read == 12)then
emis_inputlist(iemislist)%periodicity = "monthly"
endif
endif
endif
if(emis_inputlist(iemislist)%type == "OceanNH3")then
if(MasterProc)write(*,*)' using OceanNH3'
USE_OCEAN_NH3=.true.
O_NH3%index=find_index("NH3",species(:)%name)
call CheckStop(O_NH3%index<0,'NH3 not found. Needed for OceanNH3 emissions')
allocate(O_NH3%emis(LIMAX,LJMAX))
allocate(O_NH3%map(LIMAX,LJMAX))
O_NH3%emis=0.0
O_NH3%map=0.0
O_NH3%sum_month=0.0
O_NH3%sum_year=0.0
end if
if (emis_inputlist(iemislist)%type == "DMS")then
if(MasterProc)write(*,*)'using DMS'
USE_OCEAN_DMS=.true.
O_DMS%index=find_index("SO2",species(:)%name)
call CheckStop(O_DMS%index<0,'SO2 not found. Needed for DMS emissions')
allocate(O_DMS%emis(LIMAX,LJMAX))
allocate(O_DMS%map(LIMAX,LJMAX))
O_DMS%emis=0.0
O_DMS%map=0.0
O_DMS%sum_month=0.0
O_DMS%sum_year=0.0
endif
if(emis_inputlist(iemislist)%type == "Special_ShipEmis")then
write(*,*)"ERROR: Special_ShipEmis no more supported. Use new format"
write(*,*)"emissions from "//trim(fname)//" will not be included"
endif
end do ! iemislist
if(Emis_mask_allocate)then
if(.not.allocated(Emis_mask))then
allocate(Emis_mask(LIMAX,LJMAX))
Emis_mask = .false.
else
!the masks has been set by init_mask
endif
endif
! init_sectors
if(USE_SECTORS_NAME /='NOTSET')then
SECTORS_NAME = trim(USE_SECTORS_NAME)
call CheckStop((SECTORS_NAME /= 'GNFR' .and. SECTORS_NAME /= 'SNAP' .and. SECTORS_NAME /= 'TEST'), &
'Only SNAP and GNFR (and TEST) can be defined as sector names, not '//trim(SECTORS_NAME))
if(Masterproc)write(*,*)"Forcing sector categories to ",trim(SECTORS_NAME)
if(Masterproc .and. SECTORS_NAME == 'TEST')write(*,*)"WARNING: TEST sectors, requires to define sectors consistently yourself"
endif
if(SECTORS_NAME=='SNAP')then
!11 sectors defined in emissions
NSECTORS = NSECTORS_SNAP
!map timefactors onto SNAP map
sec2tfac_map => SNAP_sec2tfac_map
sec2hfac_map => SNAP_sec2hfac_map
sec2split_map => SNAP_sec2split_map
else if(SECTORS_NAME=='GNFR')then
!13 sectors defined in emissions
NSECTORS = NSECTORS_GNFR
!map timefactors onto GNFR map
sec2tfac_map => GNFR_sec2tfac_map
sec2hfac_map => GNFR_sec2hfac_map
sec2split_map => GNFR_sec2split_map
else if(SECTORS_NAME=='TEST')then
!11 sectors defined in emissions
NSECTORS = NSECTORS_TEST
!map timefactors onto TEST map
sec2tfac_map => TEST_sec2tfac_map
sec2hfac_map => TEST_sec2hfac_map
sec2split_map => TEST_sec2split_map
else
call StopAll("Sectors not defined")
end if
call CheckStop(NSECTORS>MaxNSECTORS, "redefine larger MaxNSECTORS")
allocate(Emis_field(LIMAX,LJMAX,10))
NEmis_id = 0
allocate(cdfemis(LIMAX,LJMAX))
allocate(nGridEmisCodes(LIMAX,LJMAX))
allocate(GridEmisCodes(LIMAX,LJMAX,NCMAX))
allocate(GridEmis(NSECTORS,LIMAX,LJMAX,NCMAX,NEMIS_FILE),stat=allocerr)
call CheckStop(allocerr /= 0, &
"EmisGet:Allocation error for GridEmis")
GridEmisCodes = -1
nGridEmisCodes = 0
GridEmis = 0.0
allocate(nlandcode(LIMAX,LJMAX),landcode(LIMAX,LJMAX,NCMAX))
nlandcode=0
landcode=0
allocate(road_nlandcode(LIMAX,LJMAX),road_landcode(LIMAX,LJMAX,NCMAX))
road_nlandcode=0
road_landcode=0
allocate(secemis(NSECTORS,LIMAX,LJMAX,NCMAX,NEMIS_FILE))
secemis=0.0
allocate(roaddust_emis_pot(LIMAX,LJMAX,NCMAX,NROAD_FILES))
roaddust_emis_pot=0.0
allocate(EmisOut(LIMAX,LJMAX,NEMIS_FILE))
EmisOut=0.0
allocate(e_fact(NSECTORS,NLAND,NEMIS_FILE))!NLAND defined in Country_Init()
e_fact=1.0
allocate(e_fact_lonlat(NSECTORS,MAXFEMISLONLAT,NEMIS_FILE))
e_fact_lonlat=1.0
if(.not.allocated(timefac))allocate(timefac(NLAND,N_TFAC,NEMIS_FILE))
if(.not.allocated(fac_ehh24x7))allocate(fac_ehh24x7(N_TFAC,24,7,NLAND))
if(.not.allocated(fac_emm))allocate(fac_emm(NLAND,12,N_TFAC,NEMIS_FILE))
if(.not.allocated(fac_min))allocate(fac_min(NLAND,N_TFAC,NEMIS_FILE))
if(.not.allocated(fac_edd))allocate(fac_edd(NLAND, 7,N_TFAC,NEMIS_FILE))
allocate(isec2SecOutWanted(0:NSECTORS))
isec2SecOutWanted = 0 !should never be used
isec2SecOutWanted(0) = 0
NSecEmisOutWanted = 0
do isec = 1, NSECTORS
if(SecEmisOutWanted(isec))then
NSecEmisOutWanted = NSecEmisOutWanted + 1
isec2SecOutWanted(isec) = NSecEmisOutWanted
endif
enddo
allocate(SecEmisOut(LIMAX,LJMAX,NEMIS_FILE,0:NSecEmisOutWanted))
SecEmisOut=0.0
call femis() ! emission factors (femis.dat file)
if(ios/=0) return
my_first_call = .false.
!=========================
call consistency_check() ! Below
!=========================
ios = 0
if(USES%DEGREEDAY_FACTORS) call DegreeDayFactors(0)! See if we have gridded SNAP-2
call EmisHeights() ! vertical emissions profile
KEMISTOP = KMAX_MID - nemis_kprofile + 1
if( USES%EMISSTACKS ) call readstacks(IO_EMIS)
if(MasterProc) then !::::::: ALL READ-INS DONE IN HOST PROCESSOR ::::
write(*,*) "Reading monthly and daily timefactors"
if(USES%GRIDDED_EMIS_MONTHLY_FACTOR)then
write(*,*)"Emissions using gridded monhtly timefactors "
write(IO_LOG,*)"Emissions using gridded monhtly timefactors "
end if
!=========================
call timefactors(year) ! => fac_emm, fac_edd
!=========================
end if
!=========================
call EmisSplit() ! In EmisGet_mod, => emisfrac
!=========================
!Must first call EmisSplit, to get nrcemis defined
if(EmisSplit_OUT)then
allocate(SplitEmisOut(LIMAX,LJMAX,nrcemis))
SplitEmisOut=0.0
end if
!=========================
call CheckStop(ios, "ioserror: EmisSplit")
! ####################################
! Broadcast monthly and Daily factors (and hourly factors if needed/wanted)
CALL MPI_BCAST(fac_cemm,8*12,MPI_BYTE,0,MPI_COMM_CALC,IERROR)
CALL MPI_BCAST(fac_emm,8*NLAND*12*N_TFAC*NEMIS_FILE,MPI_BYTE,0,MPI_COMM_CALC,IERROR)
CALL MPI_BCAST(fac_edd,8*NLAND*7*N_TFAC*NEMIS_FILE,MPI_BYTE,0,MPI_COMM_CALC,IERROR)
CALL MPI_BCAST(fac_ehh24x7,8*N_TFAC*24*7*NLAND,MPI_BYTE,0,MPI_COMM_CALC,IERROR)
!define fac_min for all processors
forall(iemis=1:NEMIS_FILE,insec=1:N_TFAC,inland=1:NLAND) &
fac_min(inland,insec,iemis) = minval(fac_emm(inland,:,insec,iemis))
if(INERIS_SNAP2) & ! INERIS do not use any base-line for SNAP2
fac_min(:,ISNAP_DOM,:) = 0.
! 4) Read emission files
! allocate for MasterProc (me:=0) only:
err1 = 0
if(MasterProc) then
if(USES%ROADDUST)then
allocate(road_globnland(GIMAX,GJMAX),stat=err7)
allocate(road_globland(GIMAX,GJMAX,NCMAX),stat=err8)
allocate(globroad_dust_pot(GIMAX,GJMAX,NCMAX),stat=err9)
allocate(RoadDustEmis_climate_factor(GIMAX,GJMAX),stat=err1)
call CheckStop(err7, "Allocation error 7 - globroadland")
call CheckStop(err8, "Allocation error 8 - globroadland")
call CheckStop(err9, "Allocation error 9 - globroad_dust_pot")
call CheckStop(err1, "Allocation error 1 - RoadDustEmis_climate_factor")
end if ! road dust
if(USES%ROADDUST)then
road_globnland(:,:)=0
road_globland(:,:,:)=0
globroad_dust_pot(:,:,:)=0.
RoadDustEmis_climate_factor(:,:)=1.0 ! default, no scaling
end if ! road dust
else
! needed for DEBUG=yes compilation options
if(USES%ROADDUST)then
allocate(road_globnland(1,1),road_globland(1,1,1),&
globroad_dust_pot(1,1,1),stat=err9)
call CheckStop(err9, "Allocation error 9 - dummy roadglob")
end if ! road dust
end if
emsum=0.0
Found_Emis_4D=0
do iemislist = 1, size( emis_inputlist(:)%name )
fractionformat = ( emis_inputlist(iemislist)%format=='fractions' )
fname=emis_inputlist(iemislist)%name
if ( fname == "NOTSET" ) cycle
if (emis_inputlist(iemislist)%periodicity /= 'once' ) cycle
38 FORMAT(A,I4,A)
if(MasterProc)write(*,38)sub//' reading emis_inputlist ',iemislist,trim(fname)
sumemis=0.0
sumemis_local(:,:)=0.0
nin = emis_inputlist(iemislist)%Nincl