-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathlightning.f90
2746 lines (2359 loc) · 106 KB
/
lightning.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
!//=========================================================================
!// Oslo CTM3
!//=========================================================================
!// Ole Amund Sovde, May 2015
!//=========================================================================
!// Lightning emissions.
!//=========================================================================
module lightning
!// ----------------------------------------------------------------------
!// MODULE: lightning
!// DESCRIPTION: Module conatining routines to calculate lightning NOx
!// emissions. Several options are available; from the first
!// CTM3 version (GMD2012) as described by Søvde et al. (2012,
!// GMD, doi:10.5194/gmd-7-175-2014) to the UCI2015 version
!// (not documented yet).
!//
!// I was not happy with the UCI2015 version, so I have combined
!// parts of it with the GMD2012 version; it produces better
!// annual cycle in total flash rates, and also acceptable
!// horizontal distribution.
!//
!// The lightning subroutine calculates the 3-D lightning source
!// for use in the SOURCE subroutine or as source term in Oslo
!// chemistry. It must be called for each met field after metdata
!// is updated.
!//
!// Contains:
!// - subroutine LIGHTNING_OAS2015
!// - subroutine LIGHTNING_GMD2012 (LIGHTNING_CDH)
!// - subroutine LIGHTNING_UCI2015
!// - subroutine LIGHTNING_ALLEN2002
!// - subroutine LIGHTDIST
!//
!// Ole Amund Sovde, March 2015, November 2014, August 2011
!// ----------------------------------------------------------------------
use cmn_precision, only: r8, r4
use cmn_parameters, only: secYear
use cmn_chem, only: INFILE_LIGHTNING !// Ott etal data
!// ----------------------------------------------------------------------
implicit none
!// ----------------------------------------------------------------------
integer, parameter :: NNLIGHT=3200
integer, parameter :: NLTYPE=4
real(r8), dimension(NNLIGHT,NLTYPE) :: LITPROFILE
!// Indices to indicate latitudes to regard as polar
integer :: JPS, JPN
!// Diagnose amount of lightning emissions
real(r8) :: totlitn, totlitfrq
!// Scaling factors.
!// The total accumulated modelled flash rates for ocean
!// and land, and the accumulated time, is needed to calculate
!// the scaling factors. As will be described in the subroutines,
!// the scaling factors convert from modeled flash rate to
!// climatologically normalised fraction of observed annual
!// flash rates.
real(r8) :: accFocean, accFland, accDT
!// Logical to flag if entrainment is included or not.
logical :: LentuExist
!// Seconds in a year
real(r8), parameter :: dtyear = secYear
!// Flag to turn on flashrate.dta production
logical, parameter :: getFlashFiles = .false.
!// ----------------------------------------------------------------------
character(len=*), parameter, private :: f90file = 'lightning.f90'
!// ----------------------------------------------------------------------
!// Save variables defined above
save
!// All is private
private
!// except the routines
public LIGHTNING_OAS2015, LIGHTNING_UCI2015, LIGHTNING_GMD2012, &
LIGHTNING_ALLEN2002
!// ----------------------------------------------------------------------
contains
!// ----------------------------------------------------------------------
subroutine getScaleFactors(scaleOcean, scaleLand)
!// --------------------------------------------------------------------
!// Scaling parameters depending on meteorological data.
!//
!// Price etal equations have units 1/minutes, while observations
!// are given in 1/s. This conversion is taken care of by the factors
!// scaleLand and scaleOcean. Thus, after multiplying with the
!// scaling factors, the units of flash rate is fraction of
!// all annual lightnings per second. Multiplying by time step
!// and annual amount of emission (e.g. 5Tg(N)), gives the LNOx
!// emissions for the time step.
!// --------------------------------------------------------------------
use cmn_size, only: IPARW, IDGRD
use cmn_met, only: metTYPE, metCYCLE, metREVNR
!// --------------------------------------------------------------------
implicit none
!// --------------------------------------------------------------------
real(r8), intent(out) :: scaleOcean, scaleLand
!// --------------------------------------------------------------------
!// ECMWF IFS cy36r1
!// T42L60: 1997 - 2010 cy36r1 average scaling factor (Ocean, Land)
real(r8), parameter :: Oecmwf_ifs_c36r1_T42 = 2.304304e-13_r8
real(r8), parameter :: Lecmwf_ifs_c36r1_T42 = 4.329157e-16_r8
!// T159L60: (T159 2007 / T42 2007 * clim.T42)
real(r8), parameter :: Oecmwf_ifs_c36r1_T159 = 3.594025e-14_r8
real(r8), parameter :: Lecmwf_ifs_c36r1_T159 = 7.061138e-17_r8
!// ECMWF oIFS cy38r1
!// T159L60: 1995 - 2005
real(r8), parameter :: Oecmwf_oifs_c38r1_T159 = 4.012695e-14_r8
real(r8), parameter :: Lecmwf_oifs_c38r1_T159 = 9.542578e-17_r8
!// --------------------------------------------------------------------
character(len=*), parameter :: subr = 'getScaleFactors'
!// --------------------------------------------------------------------
!// Initialise
scaleOcean = 0._r8
scaleLand = 0._r8
if (trim(metTYPE) .eq. 'ECMWF_IFS') then
!// ECMWF IFS
!// -----------------------------------------------------------------
if (metCYCLE .eq. 36 .and. metREVNR .eq. 1) then
!// Cycle 36r1
if (IPARW .eq. 128) then
if (IDGRD .eq. 1) then
scaleOcean = Oecmwf_ifs_c36r1_T42
scaleLand = Lecmwf_ifs_c36r1_T42
else if (IDGRD .eq. 2) then
!// Scaling factors based on cy36 2007 meteorology
scaleOcean = Oecmwf_ifs_c36r1_T42 * 1.642020_r8
scaleLand = Lecmwf_ifs_c36r1_T42 * 1.484240_r8
else if (IDGRD .eq. 4) then
!// Scaling factors based on cy36 2007 meteorology
scaleOcean = Oecmwf_ifs_c36r1_T42 * 2.309359_r8
scaleLand = Lecmwf_ifs_c36r1_T42 * 2.010008_r8
else
write(6,'(a,2i7)') f90file//':'//subr// &
': IPARW/IDGRD is unknown: ',iparw,idgrd
write(6,'(a)') ' metTYPE: '//trim(metTYPE)
write(6,'(a,i5)') ' metCYCLE:',metCYCLE
write(6,'(a,i5)') ' metREVNR:',metREVNR
stop 'STOP in '//subr
end if
else if (IPARW .eq. 320) then
if (IDGRD .eq. 1) then
scaleOcean = Oecmwf_ifs_c36r1_T159
scaleLand = Lecmwf_ifs_c36r1_T159
else if (IDGRD .eq. 2) then
!// Scaling factors based on cy36 2007 meteorology
scaleOcean = Oecmwf_ifs_c36r1_T159 * 2.210675_r8
scaleLand = Lecmwf_ifs_c36r1_T159 * 2.011487_r8
else if (IDGRD .eq. 4) then
!// Scaling factors based on cy36 2007 meteorology
scaleOcean = Oecmwf_ifs_c36r1_T159 * 3.350172_r8
scaleLand = Lecmwf_ifs_c36r1_T159 * 2.907207_r8
else
write(6,'(a,2i7)') f90file//':'//subr// &
': IPARW/IDGRD is unknown: ',iparw,idgrd
write(6,'(a)') ' metTYPE: '//trim(metTYPE)
write(6,'(a,i5)') ' metCYCLE:',metCYCLE
write(6,'(a,i5)') ' metREVNR:',metREVNR
stop 'STOP in '//subr
end if
else
write(6,'(a,2i7)') f90file//':'//subr// &
': scaleOcean and scaleLand not '// &
'defined for current horizontal resolution'
stop 'STOP in '//subr
end if
else !// if (metCYCLE .eq. 36 .and. metREVNR .eq. 1) then
write(6,'(a,2i7)') f90file//':'//subr// &
': metCYCLE and metREVNR not not defined'
write(6,'(a)') ' metTYPE: '//trim(metTYPE)
write(6,'(a,i5)') ' metCYCLE:',metCYCLE
write(6,'(a,i5)') ' metREVNR:',metREVNR
stop 'STOP in '//subr
end if
else if (trim(metTYPE) .eq. 'ECMWF_oIFS' .or. &
trim(metTYPE) .eq. 'ECMWF_oIFSnc4') then
!// ECMWF OpenIFS
if (metCYCLE .eq. 38 .and. metREVNR .eq. 1) then
!// Cycle 38r1
if (IPARW .eq. 320) then
if (IDGRD .eq. 1) then
scaleOcean = Oecmwf_oifs_c38r1_T159
scaleLand = Lecmwf_oifs_c38r1_T159
else if (IDGRD .eq. 2) then
!// Scaling factors based on cy38r1 2005 meteorology
scaleOcean = Oecmwf_oifs_c38r1_T159 * 2.319733_r8
scaleLand = Lecmwf_oifs_c38r1_T159 * 2.041230_r8
else if (IDGRD .eq. 4) then
!// Scaling factors based on cy38r1 2005 meteorology
scaleOcean = Oecmwf_oifs_c38r1_T159 * 5.221729_r8
scaleLand = Lecmwf_oifs_c38r1_T159 * 4.143434_r8
else
write(6,'(a,2i7)') f90file//':'//subr// &
': IPARW/IDGRD is unknown: ',iparw,idgrd
write(6,'(a)') ' metTYPE: '//trim(metTYPE)
write(6,'(a,i5)') ' metCYCLE:',metCYCLE
write(6,'(a,i5)') ' metREVNR:',metREVNR
stop 'STOP in '//subr
end if
else
write(6,'(a,2i7)') f90file//':'//subr// &
': scaleOcean and scaleLand not '// &
'defined for current horizontal resolution'
stop 'STOP in '//subr
end if
else if (metCYCLE .eq. 40 .and. metREVNR .eq. 1) then
!// Cycle 40r1 - assume same as cy38 for now
if (IPARW .eq. 320) then
if (IDGRD .eq. 1) then
scaleOcean = Oecmwf_oifs_c38r1_T159
scaleLand = Lecmwf_oifs_c38r1_T159
else if (IDGRD .eq. 2) then
!// Scaling factors based on cy38r1 2005 meteorology
scaleOcean = Oecmwf_oifs_c38r1_T159 * 2.319733_r8
scaleLand = Lecmwf_oifs_c38r1_T159 * 2.041230_r8
else if (IDGRD .eq. 4) then
!// Scaling factors based on cy38r1 2005 meteorology
scaleOcean = Oecmwf_oifs_c38r1_T159 * 5.221729_r8
scaleLand = Lecmwf_oifs_c38r1_T159 * 4.143434_r8
else
write(6,'(a,2i7)') f90file//':'//subr// &
': IPARW/IDGRD is unknown: ',iparw,idgrd
write(6,'(a)') ' metTYPE: '//trim(metTYPE)
write(6,'(a,i5)') ' metCYCLE:',metCYCLE
write(6,'(a,i5)') ' metREVNR:',metREVNR
stop 'STOP in '//subr
end if
else
write(6,'(a,2i7)') f90file//':'//subr// &
': scaleOcean and scaleLand not '// &
'defined for current horizontal resolution'
stop 'STOP in '//subr
end if
else !// if (metCYCLE .eq. 38 .and. metREVNR .eq. 1) then
write(6,'(a,2i7)') f90file//':'//subr// &
': metCYCLE and metREVNR not not defined'
write(6,'(a)') ' metTYPE: '//trim(metTYPE)
write(6,'(a,i5)') ' metCYCLE:',metCYCLE
write(6,'(a,i5)') ' metREVNR:',metREVNR
stop 'STOP in '//subr
end if
else
write(6,'(a,2i7)') f90file//':'//subr//': metTYPE not defined!'
write(6,'(a)') ' metTYPE: '//trim(metTYPE)
write(6,'(a,i5)') ' metCYCLE:',metCYCLE
write(6,'(a,i5)') ' metREVNR:',metREVNR
stop 'STOP in '//subr
end if !// if (trim(metType)...)
!// --------------------------------------------------------------------
end subroutine getScaleFactors
!// ----------------------------------------------------------------------
!// ----------------------------------------------------------------------
subroutine LIGHTNING_OAS2015(NDAY,NDAYI,dt_met,LNEW_MONTH)
!// --------------------------------------------------------------------
!// Lightning NOx
!//
!// Algorithm:
!// This routine combines parts of GMD2012 method with the
!// new UCI2015 treatment. It has better annual cycle in
!// flash rates.
!//
!// Ole Amund Sovde, March 2015
!// --------------------------------------------------------------------
use cmn_size, only: IPARW,IPAR,JPAR,LPAR,LPARW, LWEPAR, IDGRD
use cmn_ctm, only: ETAA, ETAB, XDGRD, YDGRD, XDEDG, YDEDG, &
AREAXY, PLAND
use cmn_chem, only: NEMLIT, LITSRC, LITFAC
use cmn_met, only: CWETE, CENTU, ZOFLE, T, PRECCNV, P
use utilities, only: get_free_fileid
!// --------------------------------------------------------------------
implicit none
!// --------------------------------------------------------------------
!// Input / Output
integer, intent(in) :: NDAY, NDAYI
logical, intent(in) :: LNEW_MONTH
real(r8), intent(in) :: dt_met
!// Locals
integer :: I, J, L
real(r8) :: CTH !// Cloud top height (top of convective mass flux)
real(r8) :: fPR !// Instant Price&Rind flash rate
!// Model flash rates over land and ocean. Their units do not
!// matter, since they are scaled to match observations, but
!// if you are curious the Price et al. equations are [fl/min].
real(r8) :: flashL(IPAR,JPAR), flashO(IPAR,JPAR)
!// We find the flash rates normalised to a climatological year.
!// The sum of normFrate*dt is 1 for a year when it represents the
!// climatological year. This is to distribute lightning emissions
!// (e.g. 5Tg(N)) using normFrate.
real(r8) :: normFrate(IPAR, JPAR)
!// normFrate is the sum of flashL and flashO scaled by their
!// respective factors, which also take care of the normalisation.
real(r8) :: scaleOcean, scaleLand
!// Vertical distribution of normFrate(i,j)
real(r8) :: vertprof(LWEPAR), finc, zfmass(LWEPAR)
!// Additional factors to calculate flash rates
real(r8) :: FACT1, FACT2
!// Meteorological variables
real(r8) :: LNCWET(LWEPAR), LNCENTU(LWEPAR), TEMP(LWEPAR)
real(r8) :: ZL(LWEPAR+1)
real(r8) :: PBOT !// Pressure at bottom of grid box
real(r8) :: MFLX !// Mass flux
real(r8) :: MENT !// Mass entrainment [kg/s]
real(r8) :: minT, maxT !// min/max temperatures in cloud
!// Diagnostics to get total emissions
real(r8) :: CNV_SUM, ntgyr
real(r8) :: NO_SUM(IPAR, JPAR)
real(r8) :: fLand, fOcean !// Total unscaled flash rates
!// Counters & indices
integer :: III, JJJ, IOS, ifnr, NLCL, NLCO !// Counters & indices
integer :: LFREE, CBL, CTL, LFLX !// Level indices
!// Pressures for defining shallow or deep convection.
real(r8), parameter :: pconvection = 500._r8
!// Observed lightning flash rated (1/s).
!// These are taken from LIS/OTD data LISOTD_HRAC_V2.3.2014.hdf,
!// as produced by Cecil et al. (2014, J. Atmos. Res.,
!// vol 135-136, 404-414, doi:10.1016/j.atmosres.2012.06.028).
!// With an annual global flash rate of 46fl/s, using land
!// fractions separates this into 9.1fl/s over ocean and
!// 36.9 over land. If we were to define all grid boxes having
!// land within 300km reach (roughly similar to Christian et al.
!// (2003, JGR, 108(D1),doi:10.1029/2002JD002347), these numbers
!// would be 4.25fl/s and 41.75fl/s, respectively.
!// Christian et al. (2003) suggested about 5fl/s over ocean.
!// However, we will use a stricter definition of land, so
!// we stick to the former.
real(r8), parameter :: obsFocean=9.1_r8, obsFland=36.9_r8
real(r8), parameter :: obsFall = obsFocean + obsFland
!// Scaling parameters depending on meteorological data.
!//
!// Price etal equations have units 1/minutes, while observations
!// are given in 1/s. This conversion is taken care of by the factors
!// scaleLand and scaleOcean. Thus, after multiplying with the
!// scaling factors, the units of flash rate is fraction of
!// all annual lightnings per second. Multiplying by time step
!// and annual amount of emission (e.g. 5Tg(N)), gives the LNOx
!// emissions for the time step.
! !//
! !// 1997 - 2010 cy36r1 average scaling factor
! real(r8), parameter :: scaleOceanT42 = 2.304304e-13_r8
! real(r8), parameter :: scaleLandT42 = 4.329157e-16_r8
!
! !//
! !// T159L60 cy36r1 (T159 2007 / T42 2007 * clim.T42)
! real(r8), parameter :: scaleOceanT159 = 3.594025e-14_r8
! real(r8), parameter :: scaleLandT159 = 7.061138e-17_r8
!// --------------------------------------------------------------------
character(len=*), parameter :: subr = 'LIGHTNING_OAS2015'
!// --------------------------------------------------------------------
!// If no lightning, return
if (NEMLIT .eq. 0) return
!// Initialize
CNV_SUM = 0._r8
NO_SUM(:,:) = 0._r8
call getScaleFactors(scaleOcean, scaleLand)
!// Get correct scaling factors
! if (IPARW .eq. 128) then
! if (IDGRD .eq. 1) then
! scaleOcean = scaleOceanT42
! scaleLand = scaleLandT42
! else if (IDGRD .eq. 2) then
! !// Scaling factors based on cy36 2007 meteorology
! scaleOcean = scaleOceanT42 * 1.642020_r8
! scaleLand = scaleLandT42 * 1.484240_r8
! else if (IDGRD .eq. 4) then
! !// Scaling factors based on cy36 2007 meteorology
! scaleOcean = scaleOceanT42 * 2.309359_r8
! scaleLand = scaleLandT42 * 2.010008_r8
! else
! print*,'lightning.f90: IDGRD is unknown:',idgrd
! stop 'in LIGHTNING_OAS2015'
! end if
! else if (IPARW .eq. 320) then
! if (IDGRD .eq. 1) then
! scaleOcean = scaleOceanT159
! scaleLand = scaleLandT159
! else if (IDGRD .eq. 2) then
! !// Scaling factors based on cy36 2007 meteorology
! scaleOcean = scaleOceanT159 * 2.210675_r8
! scaleLand = scaleLandT159 * 2.011487_r8
! else if (IDGRD .eq. 4) then
! !// Scaling factors based on cy36 2007 meteorology
! scaleOcean = scaleOceanT159 * 3.350172_r8
! scaleLand = scaleLandT159 * 2.907207_r8
! else
! print*,'lightning.f90: IDGRD is unknown:',idgrd
! stop 'in LIGHTNING_OAS2015'
! end if
! else
! print*,'* lightning.f90: scaleOcean and scaleLand not '// &
! 'defined for current horizontal resolution'
! stop 'in LIGHTNING_OAS2015'
! end if
!// Flash rates
normFrate(:,:) = 0._r8
FlashL(:,:) = 0._r8
FlashO(:,:) = 0._r8
!// Lightning source
LITSRC(:,:,:) = 0._r8
!=================================================================
! Read lightning CDF from Ott et al [JGR, 2010]. (ltm, 1/25/11)
!=================================================================
if ( LNEW_MONTH .and. NDAY.eq.NDAYI ) then
!// Initialize lightning profiles
LITPROFILE(:,:) = 0._r8
totlitfrq = 0._r8
totlitn = 0._r8
accFocean = 0._r8
accFland = 0._r8
accDT = 0._r8
!// Echo info
write(*,*) ' - INIT_LIGHTNING: Reading '//TRIM( INFILE_LIGHTNING )
!// Get unused file ID
IFNR = get_free_fileid()
!// Open file containing lightning PDF data
open( IFNR, file=trim( INFILE_LIGHTNING ), status='OLD', iostat=IOS)
if ( IOS /= 0 ) stop f90file//':'//subr//': lightdist not read:1'
!// Read 12 header lines
do III = 1, 12
read( IFNR, '(a)', IOSTAT=IOS )
if ( IOS /= 0 ) stop f90file//':'//subr//': lightdist not read:2'
end do
!// Read NNLIGHT types of lightning profiles
do III = 1, NNLIGHT
read( IFNR,*,IOSTAT=IOS) (LITPROFILE(III,JJJ),JJJ=1,NLTYPE)
end do
!// Close file
close( IFNR )
!// Generate flash files for calculating scaling factors
if (getFlashFiles) then
open(IFNR,file='flashrate.dta',form='unformatted')
write(IFNR) IPAR,JPAR
write(IFNR) XDGRD,YDGRD, AREAXY, PLAND,XDEDG,YDEDG
close(IFNR)
end if
!// Possible extended definition of land grid boxes when
!// considering lightning.
!call distland(landlit)
if (maxval(CENTU) .gt. 0._r8) then
LentuExist = .true.
else
LentuExist = .false.
end if
end if !// if ( LNEW_MONTH .and. NDAY.eq.NDAYI ) then
!// Locate layer of approximately 1500m above surface, which is
!// about layer 12 in native vertical L60. If collapsing this
!// will be layer 9 (as in UCI standard code).
LFREE = 12 - (LPARW-LPAR)
!=================================================================
! Flashrate: Horizontal and temporal distribution
!=================================================================
!$omp parallel private (I,J,L) &
!$omp private (CBL,CTL,LFLX) &
!$omp private (CTH,PBOT) &
!$omp private (fPR,finc,FACT1,FACT2,minT,maxT,MFLX,MENT) &
!$omp private (LNCWET,LNCENTU, ZL,TEMP) &
!$omp private (vertprof) &
!$omp shared (CWETE,CENTU,T,ZOFLE,PRECCNV) &
!$omp shared (P,ETAA,ETAB) &
!$omp shared (LITFAC,LITSRC,normFrate,NO_SUM) &
!$omp shared (LentuExist, LFREE, JPS, JPN, PLAND) &
!$omp shared (scaleLand, scaleOcean) &
!$omp shared (flashO, flashL) &
!$omp default(NONE)
!$omp do
do J = 1,JPAR
do I = 1,IPAR
!// Store updraft convection.
!// I do not know the reason behind this parameterisation yet.
!// Find convection at 850hPa. For all layers above this, we
!// do not allow convective mass flux to be larger.
LNCWET(1:(LFREE-1)) = 0._r8
LNCENTU(1:(LFREE-1)) = 0._r8
do L = LFREE, LWEPAR
LNCWET(L) = CWETE(I,J,L)
LNCENTU(L) = CENTU(I,J,L)
end do
!// Use rain as screening. There may be some convection
!// within large scale rain, without producing convective
!// rain falling to the ground. It is not clear to what extent
!// that can happen for vigorous convection. Therefore I
!// check only convective rain at the surface; it seems to
!// be a filter that works well to get temporal distribution
!// fairly correct.
if (PRECCNV(I,J,1) .eq. 0._r8) cycle
!// There may be some instances where convective mass flux
!// occur in single layers or separated into several intervals.
!// This is due to the nature of the IFS model physics.
!// If e.g. level 10 and 20 are the only layers with
!// convective mass flux, there should be no lightning.
call filterLNC(I, J, LFREE, LNCWET, LNCENTU, CBL, CTL)
!// Check if there are any entraining updrafts;
!// if none, then cycle to next box.
if (CTL .eq. 0) cycle
!// There are several scaling factors or additional filters
!// that should be included, to sort out the flash rates.
!//
!// 1. Temperature in cloud: Cloud should consist of both
!// liquid and ice. This can be checked with temperature.
!// 2. Cloud fraction scaling: The convective cell does not
!// cover the whole grid box.
!// 3. Updraft velocity: Lightning needs large vertical
!// velocities in the convective cell.
!// 4. Area scaling: A convective cloud with a certain cloud
!// top height may cover a larger area at low latitudes,
!// giving higher flash rate there.
!// 5. Flux scaling: This may take horizontal extent of the
!// convective cloud into account.
!// Final: scaleLand and scaleOcean: to scale to observed flash
!// rates (i.e. normalised to climatological year).
!// Retrieve temperature as 1D array
do L = 1, LWEPAR
TEMP(L) = T(I,J,L)
end do
!// 1. Temperature
!// ------------------------------------------------------------
!// Make sure the cloud is mixed phase, containing both
!// liquid and ice. We assume this occurs when minimum
!// temperature in the column is colder than -40C and
!// the maximum temperature is warmer than 0C.
!// To allow for some uncertainty, UCI scales the first
!// factor linearly from 0 at -30C to 1 at -40C, and the
!// latter from 1 at 0C to 0 at -20C.
!//
!// I think the first scaling is ok, but the latter scaling
!// in the UCI method, uses surface temperature. In any case
!// I think allowing lightning down to -20C sounds very weird.
!//
!// A better choice is to require the cloud to have max
!// temperatures above -2.5C before producing lightning;
!// taking into account some uncertainty in temperature: The
!// metdata temperature is gridbox average temperatures, not
!// necessarily the correct updraft temperatures.
!// Strict temperature control may not very physical, but could
!// be an alternative.
!//
!// Surface temperature scalings?
!// Using the surface temperature to scale lightning, will
!// probably give wrong diurnal cycle; there should be more
!// flashes at local night time. Instead of such a scaling,
!// I use a linear scale from 0 at -2.5C to 1 at 7.5C.
!// In this way we should increase the more vigorous lightning
!// and reduce the lightning from relatively cold clouds.
minT = minval(TEMP(CBL:CTL))
maxT = maxval(TEMP(CBL:CTL))
!// Linear scaling: T>=7.5: FACT1=1, T<-2.5: FACT1=0
FACT1 = max(0._r8, min(0.1_r8*(maxT - 270.5_r8), 1._r8) )
!// Check if criteria is met
if ( FACT1 .eq. 0._r8 ) cycle
!// Linear scaling for cloud minimum temperature:
!// T=<-40C: FACT2=1, T>-30C: FACT2=0
FACT2 = max(0._r8, min(-0.1_r8*(minT - 243._r8), 1._r8) )
!// Check if criteria is met
if ( FACT2 .eq. 0._r8 ) cycle
!// 2. Cloud fraction scaling?
!// ------------------------------------------------------------
!// The flux itself, in values of kg/s may be less meaningful;
!// when the same flux covers a smaller area, convection should
!// be more rigorous.
!// But convection covering a larger area is likely producing
!// more lightning. One way to reduce lightning frequency over
!// small areas is to multiply with cloud fraction.
!//
!// But, unfortunately, cloud fraction is an instant field,
!// whereas convective mass flux is based on accumulated values,
!// so there may be times where we have convective mass flux
!// but no clouds. Testing so far indicates that including
!// this factor makes the horizontal distribution worse.
!// 3. Updraft velocity
!// ------------------------------------------------------------
!// For a thunderstorm convective cell, updraft velocities
!// typically reach > 5m/s. However, since we do not know the
!// size of the convective cell (cloud fraction cannot be
!// used since it is instant and flux is accumulated), we
!// could in principle use the grid box mass. That would
!// not give a "real" updraft velocity, but how fast the
!// the mass in the grid box would move to match the flux.
!// -> This would be resolution dependent, so the velocity
!// would have different meaning in different resolutions.
!//
!// A possible work-around is to calculate the fraction of
!// moved mass to some fixed mass value that does not change
!// across resolutions. As a good proxy, we can use the total
!// mass per level. The total mass per level should
!// change minimally between resolutions.
!// But unfortunately, finding a certain limit that work well
!// in different resolutions is tricky. I tried it, but could
!// not get it working.
!// Restrict lightning to be produced only when convection
!// reaches a height where pressure is less than some
!// limit (pconvection).
!// Also require that entrainment occurs above this level (it
!// usually does).
LFLX = 0
do L = CBL, CTL
PBOT = ETAA(L) + P(I,J)*ETAB(L)
if (LentuExist) then
MENT = sum(LNCENTU(L:CTL))
if (PBOT .le. pconvection .and. MENT .gt. 0._r8) then
LFLX = L
exit !// Exit L-loop
end if
else
if (PBOT .le. pconvection) then
LFLX = L
exit !// Exit L-loop
end if
end if
end do
!// Skip lightning if convection does not reach the limit.
if (LFLX .eq. 0) cycle
!// 4. Area scaling?
!// ------------------------------------------------------------
!// We find frequency based on cloud top height, not considering
!// the horizontal extent of the convection. A cloud with a
!// certain size should have the same frequency at Equator as at
!// high latitudes,
!// However, at Equator the convective cloud may generally be
!// larger, or there may even be several convective clouds in
!// the same grid box. So it can be argued that the modelled
!// frequency should be weighted by the fractional area relative
!// to Equator.
!// Again, cloud fraction would solve this issue.
!// To take into account the horizontal extent of the convective
!// plume, we might scale with the convective mass flux at some
!// level. They should not be used together, because the flux
!// do carry some information about size.
!// 5. Flux scaling
!// ------------------------------------------------------------
!// If the flux is large, the convective cloud may be larger and
!// should produce more lightning. We can scale the flash rates
!// with the flux at a certain level, e.g. at 500hPa.
!// This produces more lightning at high latitudes than we want.
!//
!// I have also found that it could be smart to screen away
!// instances where the entrainment below this level
!// (i.e. at LFLX-1) is small, e.g. 10% of the flux at LFLX.
!// However, it seems that in cy36r1 had some changes to the
!// entrainment from 2009-2012. I remember the ECMWF corrected
!// a few things, but do not know if that should affect us.
!// So I will skip this screening for now, only require
!// entrainment > 0.
!// Flux into LFLX level
!MFLX = LNCWET(LFLX)
!// Entrainment into LFLX-1 level
!MENT = LNCENTU(LFLX-1)
!// Final: scaleLand and scaleOcean
!// ------------------------------------------------------------
!// Price etal equations have units 1/minutes. The final scaling
!// convert to 1/s, and also convert to fraction of climatological
!// mean, so that the sum for a year on average will be 1.
!// This is taken care of by the factors scaleLand and
!// scaleOcean.
!// Find layer bottom height above ground [km]
!// ------------------------------------------------------------
do L = 1, LWEPAR+1
ZL(L) = 1.e-3_r8 * (ZOFLE(L,I,J) - ZOFLE(1,I,J))
end do
!// Now use Price et al. equations
!// ------------------------------------------------------------
!// Cloud top height above ground is top of level CTL [km]
CTH = ZL(CTL+1)
if (PLAND(I,J) .gt. 0.25_r8) then
fPR = CTH**4.92_r8 !// 3.44d-5 included in scaleLand
finc = FACT1 * FACT2
normFrate(I,J) = finc * fPR * scaleLand
flashL(I,J) = finc * fPR !// Save unscaled land
else
fPR = CTH**1.73_r8 !// 6.4d-4 included in scaleOcean
finc = FACT1 * FACT2
normFrate(I,J) = finc * fPR * scaleOcean
flashO(I,J) = finc * fPR !// Save unscaled oceanic
end if
!---------------------------------------------------------
! NOx emission, vertical distribution
!---------------------------------------------------------
! If there's lightning in the column
if ( normFrate(I,J) .gt. 0._r8 ) then
!// Partition NOx emissions vertically using the PDFs from
!// Ott et al. (2010) and regional definitions from
!// Allen et al. (2010)
call lightdist(I, J, CTL, ZL, normFrate(I,J), vertprof)
!// Assign emissions into array for later use.
!// (units: fraction of annual emission per second)
do L=1,CTL
LITSRC(L,I,J) = vertprof(L)
!// Sum up NO to test 5 Tg(N)/year
NO_SUM(I,J) = NO_SUM(I,J) + (vertprof(L)*LITFAC(1))
end do
end if
end do
end do !// do J = 1,JPAR
!$omp end do
!$omp end parallel
!=================================================================
! Diagnostics for lightning production
!=================================================================
!// Total unscaled land and ocean flashes
fLand = sum(flashL)
fOcean = sum(flashO)
!// Flash files
if (getFlashFiles) then
IFNR = get_free_fileid()
open(IFNR,file='flashrate.dta',form='unformatted', position='append')
write(IFNR) real(flashL, r4), real(flashO, r4)
close(IFNR)
end if
!// How many grid boxes experienced lightning?
NLCL = 0
NLCO = 0
do J = 1, JPAR
do I = 1, IPAR
if (flashL(I,J) .gt. 0._r8) NLCL = NLCL + 1
if (flashO(I,J) .gt. 0._r8) NLCO = NLCO + 1
end do
end do
write(*,'(A,3i7)') 'Lightning: # grid boxes w/flash (L/O/L+O): ', &
NLCL, NLCO, NLCL+NLCO
!// Total global fraction of annual total flashrate. Must
!// multiply by sum up to make annual total 1
cnv_sum = sum(normFrate) * dt_met
!// Accumulate this fraction
totlitfrq = totlitfrq + cnv_sum
!// Sum up L-NOx until next update of LITSRC (duration dt_met)
totlitn = totlitn + sum(no_sum) * dt_met * 1.e-9_r8 * 14._r8/30._r8
!// Emissions this step, scaled to Tg(N)/yr
ntgyr = sum(NO_SUM) * 1.e-9_r8 * 14._r8/30._r8 * dtyear
!// Accumulate flashes during time step and total time.
accFland = accFland + fLand * dt_met
accFocean = accFocean + fOcean * dt_met
accDT = accDT + dt_met
!// Write total flash rate for ocean and land
! write(*,'(A,f12.3,1x,f12.3)')
! & 'Lightning: Flashrates (1/s) ocean/land:',
! & fOcean*scaleOcean, fLand*scaleLand
!// Fraction of annual total flashrates
write(*,'(A,es23.17)') 'Lightning: Fraction of annual flashes: ', cnv_sum
!// Emissions during this time step (as Tg(N)/yr) and
!// also the accumulated amount
write(*,'(A,2f12.5)') 'Lightning: Tg(N)/yr & accumulated Tg(N):', &
ntgyr, totlitn
!// Accumualted flash frequency fraction of annual total, both
!// new and old values. Their difference is cnv_sum.
write(*,'(A,es23.17,1x,es23.17)') 'Lightning: Accu.flashfreq n/o: ', &
totlitfrq,totlitfrq-cnv_sum
!// Scaling factors. These are the values you need to modify
!// scaleOean and scaleLand if you use other meteorological data.
!// IMPORTANT:
!// Remember that you need to use at least an annual mean, but
!// preferably a mean over several years of data.
!//
!// Note also that for a single year of meteorological data,
!// this printout may differ slightly from the applied scaling
!// factors. This does not mean that the scalings are wrong,
!// it may be that the applied scalings were calculated using
!// a different year or several years.
!//
!// Printout is done after calculating lightning for the
!// last meteorological time step of each day, giving the
!// average scalings for days from NDAYI through NDAY.
!//
!// Will find scaling so that
!// annual_modelflashes * scale = 1.
!// separated into land and ocean:
!// scaleO = obsFocean/obsFall / annual_model_Oflashes
!// scaleL = obsFland/obsFall / annual_model_Lflashes
if (mod(accDT,86400._r8) .eq. 0._r8) &
write(*,'(A,i4,es16.7,1x,es16.7)') &
'Lightning: Accu.scal. NDAY/ocean/land: ',NDAY, &
obsFocean/obsFall / (accFocean/accDT * dtyear), &
obsFland/obsFall / (accFland/accDT * dtyear)
!// Sensiblitity check
if (ntgyr .lt. 1._r8 .or. ntgyr .gt. 12._r8) then
print*,'*** WARNING: lightning.f90: Low/High LNOx: '// &
'Check scaleLand and scaleOcean!'
end if
!// --------------------------------------------------------------------
end subroutine LIGHTNING_OAS2015
!// ----------------------------------------------------------------------
!// ----------------------------------------------------------------------
subroutine LIGHTNING_UCI2015(NDAY,NDAYI,dt_met,LNEW_MONTH)
!// --------------------------------------------------------------------
!// Comment by Ole Amund Sovde:
!// This method is NOT to be used in CTM3. While the method
!// produces fairly ok annual horizontal distribution of flashes,
!// the temporal distribution is almost gone. In addition, the
!// number of lightning days are very high in many places, giving
!// small amounts of lightning emissions very often.
!// Better when filtering convective events using the filterLNC.
!//
!// About this routine and why to dismiss it:
!// This is from the qcode version 7.1. It uses the Price and Rind
!// (1992) equations slightly differently than before. Their method
!// was updated so that that the flash frequency in a 2x2 window
!// should be equal to the flash frequency summed up over its four
!// 1x1 boxes.
!// Clearly the max level of convection will be the same.
!// Say we have fluxes F1, F2, F3, F4 and cloud tops
!// H1, H2, H3, H4. Ideally, we could have:
!// FL = (F1*H1^m + F2*H2^m + F3*H3^m + F4*H4^m)/(F1+F2+F3+F4)
!// where m is the Price etal exponents.
!//
!// But how should we determine the individual F-values? When
!// running 2x2 we only have F=F1+F2+F3+F4 available.
!// The UCI method tries to look sum up several "clouds", where
!// it is assumed that each time the flux is reduced, a new
!// cloud top is reached. Thus, if the flux decreases in 10
!// consecutive layers, all 10 will count. Probably the largest
!// will give the most contribution. UCI also considers the
!// temperature of the layer above each cloud top, which needs
!// to be cold.
!// To me this seems to cover only a part of the picture. What they
!// do is a quasi-calculation of "detrainment", and seems not
!// completely correct.
!//
!// I have tried to calculate the true detrainment from LFLX
!// upwards, and then scale the normalized flashrate with the
!// sum of detrainment, but to no help. Dropping it seems better.
!//
!// For UCI method, using daily flashes f(i,j,d), scaled to match
!// 1.45d9 fl/yr, i.e. sum(f) = 1. We multiply f by 5d9*30/14 to
!// get total NO, and then plot e.g. Oslo (i=5,j=54 in T42).
!// This gives about 160 lightning days, where ~60 have 1d3kg(NO)
!// and ~5 have 1d4kg(NO).
!// 100kg(NO) per day occurs very frequently. Is that a problem?
!// Probably not very, but if the background is clean, it may very
!// well have an impact.
!//
!//
!// The global lightning flash rate is observed to be 46 flashes/s,
!// of which 76% (35 flashes/s) occurs over land. This is derived
!// from the space-based OTD and LIS sensors, which operated over
!// 1998-2005. We scale the flash rate derived from cloud top height
!// and other criteria to match this mean rate over land and ocean.
!// The scale factors need to be recalculated for each change of
!// meteorology version and resolution, or whenever other parts of
!// the algorithm change. CDH calculated factors for EC cycle 36r1,
!// T42L60, over 1999-2009.
!// --------------------------------------------------------------------
!// Horizontal distribution is based on:
!// Price & Rind (1992), doi:10.1029/92JD00719
!// Vertical distribution is based on:
!// Ott et al (2010), doi:10.1029/2009JD011880
!//
!// Ole Amund Sovde, March 2015
!// --------------------------------------------------------------------
use cmn_size, only: IPARW,IPAR,JPAR,LPAR,LPARW, LWEPAR, IDGRD
use cmn_ctm, only: ETAA, ETAB, XDGRD, YDGRD, XDEDG, YDEDG, &
AREAXY, PLAND
use cmn_chem, only: NEMLIT, LITSRC, LITFAC
use cmn_met, only: CWETE, CENTU, ZOFLE, T, PRECCNV, P
use utilities, only: get_free_fileid
!// --------------------------------------------------------------------
implicit none
!// --------------------------------------------------------------------
!// Input / Output
integer, intent(in) :: NDAY, NDAYI
logical, intent(in) :: LNEW_MONTH
real(r8), intent(in) :: dt_met
!// Locals
integer :: I, J, L
!// Model flash rates over land and ocean. Their units do not
!// matter, since they are scaled to match observations, but
!// if you are curious the Price et al. equations are [fl/min].
real(r8) :: flashL(IPAR,JPAR), flashO(IPAR,JPAR)
!// We find the flash rates normalised to a climatological year.
!// The sum of normFrate*dt is 1 for a year when it represents the
!// climatological year. This is to distribute lightning emissions
!// (e.g. 5Tg(N)) using normFrate.
real(r8) :: normFrate(IPAR, JPAR)
!// normFrate is the sum of flashL and flashO scaled by their
!// respective factors, which also take care of the normalisation.
real(r8) :: scaleOcean, scaleLand
!// Vertical distribution of normFrate(i,j)
real(r8) :: vertprof(LWEPAR)
!// Additional factors to calculate flash rates
real(r8) :: FACT1, FACT2, finc, fPR, MENT, PBOT
!// Meteorological variables
real(r8) :: LNCWET(LWEPAR), TEMP(LWEPAR)
real(r8) :: ZL(LWEPAR+1)
real(r8) :: LNCENTU(LWEPAR) !// Used in filterLNC
!// Diagnostics to get total emissions
real(r8) :: CNV_SUM, ntgyr
real(r8) :: NO_SUM(IPAR, JPAR)
real(r8) :: fLand, fOcean !// Total unscaled flash rates