-
Notifications
You must be signed in to change notification settings - Fork 0
/
FastJ_mod.f90
3855 lines (3374 loc) · 152 KB
/
FastJ_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
! ---------- JX71_notes.f90
! same as JX70b except that test for solar flux = 0 so that can run 200 nm data set
! the 2200 nm data set cuts all J's at exactly 200 nm - used to combien with WACCM tables
! to allow J's above 60 km, including Lyman-alpha
!
! ----subroutines and calls:
! main standalone
! >>> call INIT_FJX (TITLJXX,NJX_,NJXX)
! call RD_JS_JX(9,'FJX_j2j.dat', TITLJXX,NJXX)
! call SOLAR_JX(GMTAU,IDAY,YGRD,XGRD, SZA,U0,SOLF)
! call ACLIM_FJX (YLAT, MONTH, PPP,TTT,ZZZ,DDD,OOO, L1_)
! call JP_ATM0(PPP,TTT,DDD,OOO,ZZZ, L_)
! >>> call PHOTO_JX(U0,SZA,REFLB,SOLF, LPRTJ, PPP,ZZZ,TTT,DDD,RRR,OOO,
! LWP,IWP,REFFL,REFFI, AERSP,NDXAER,L1U,ANU, VALJXX,NJXU)
!
! ----notes: >>> = only two essential fast-JX calls are denoted above with >>>
!
! >>> INIT_FJX (TITLJXX,NJX_,NJXX)
! called once to read in data files, returns thru calling sequence:
! TITLJXX(1:NJX_) the char*6 title for each fast-JX J-value
! NJXX the actual number of J-values that fast-JX will compute.
!
! RD_JS_JX(9,'FJX_j2j.dat', TITLJXX,NJXX)
! called once after INIT_FJX to map the CTM J's onto the fast-JX J's
! this is an example, see the data faile 'FJX_j2j.dat'
!
! >>> PHOTO_JX(U0,SZA,REFLB,SOLF, LPRTJ, PPP,ZZZ,TTT,DDD,RRR,OOO,
! CLDWP,AERSP,NDXCLD,NDXAER,L1_,AN_, VALJXX,NJX_)
! called every time for each Indep. Colm. Atmos. (ICA) to compute J's
! all the information is passed through the calling sequence:
! ----------------------------------------------------------------------------
! L1_ CTM has layers 1:L_, for JX must add top-of-atmos layer, L1_ = L_+1
! AN_ dimension of max number of aerosol types given to JX
! U0 cosine of the solar zenith angle
! SZA solar zenith angle (degrees)
! REFLB Lambertian reflectance of lower boundary
! (NB JX could be adapted to add lower layers for canopy/snow scatt/abs)
! SOLF Solar radiation factor, correcting for sun-earth distance
! (NB JX could be modified for solar cycle, but needs link to internal data)
! LPRTJ logical to produce internal printout (stdout) of fast-JX J's, fluxes,
! absorption, etc. This information is now only internal to PHOTO_JX
! PPP(1:L1_+1) edge press (hPa)
! ZZZ(1:L1_+1) edge altitude (cm)
! TTT(1:L1_) mid-layer temp (K)
! DDD(1:L1_) layer dens-path (# molec /cm2)
! RRR(1:L1_) mid-layer relative humidity (0.0 to 1.0)
! OOO(1:L1_) layer O3 path (# O3 /cm2)
! CLDWP(1:L1_) layer cloud water path (kg/m2), liquid and ice
! AERSP(1:L1_,1:AN_) aerosol path (g/m2)
! NDXCLD(1:L1_) layer cloud index (type)
! only a single cloud type is allowed for optical properties, pick
! the dominant one in terms of optical depth,
! see notes on cloud types allowed in fast-JX: 'FJX_scatt-cld.dat'
! NDXAER(1:L1_,1:AN_) aerosol index (type)
! sample aerosol types are in 'FJX_scat-aer.dat' and 'FJX_scat-UMa.dat'
! the UMa data allows for relative humidity to be included
! other aerosol types can be added.
! ----------------------------------------------------------------------------
! VALJXX(1:,NJX_,1:L) & NJX_ (first dimension of VALJXX) are returned
! VALJXX is the array of fast-JX J's, the second dimension is not given
! but is large enough to accommodate the CTM layers 1:L1_
! the main code must use the information calcualted by RD_JS_JX to
! re-map the VALJXX onto the CTM J's. A useful example is given.
! ----------------------------------------------------------------------------
!
! SOLAR_JX calculates solar zenith angle & distance correction (if needed)
!
! ACLIM_FJX fills in T & O3 from a climatology
! may be needed for the layer above the CTM to account for O3 & O2
!
! JP_ATM0 does a simple printout (stdout) of the atmosphere
!
!
! ---------- fjx70sub.f fast-JX core routines ver 7.0+ (10/2012, mjp)
!
! ----subroutines and calls: >>> only subroutines called from outside >>>
! one include 'cmn_FJX.f' is common to several and has parameters, etc.
! only other connection with CTM code is in call sequence and noted above.
!
! >>> subroutine INIT_FJX (TITLEJXX,NJXU,NJXX)
! call RD_XXX(JXUNIT,'FJX_spec.dat')
! call RD_MIE(JXUNIT,'FJX_scat.dat')
! call RD_UM (JXUNIT,'FJX_UMaer.dat')
! call RD_PROF(JXUNIT,'atmos_std.dat')
! subroutine RD_XXX(NJ1,NAMFIL)
! subroutine RD_MIE(NJ1,NAMFIL)
! subroutine RD_UM(NJ1,NAMFIL)
! subroutine RD_PROF(NJ2,NAMFIL)
! subroutine EXITC(T_EXIT)
! >>> subroutine SOLAR_JX(GMTIME,NDAY,YGRDJ,XGRDI, SZA,COSSZA,SOLFX)
! >>> subroutine ACLIM_FJX (YLATD, MONTH, PPP,TTT,ZZZ,DDD,OOO, L1U)
! >>> subroutine PHOTO_JX(U0,SZA,REFLB,SOLF,LPRTJ, PPP,ZZZ,TTT,DDD,RRR,OOO,
! CLDWP,AERSP,NDXCLD,NDXAER,L1U,ANU, VALJXX,NJXU)
! call SPHERE2 (U0,RAD,ZZJ,ZZHT,AMF2, L1U,JXL1_)
! call OPTICL (OPTX,SSAX,SLEGX, ODCLD,NDCLD)
! call OPTICA (OPTX,SSAX,SLEGX, PATH,RH, NAER)
! call OPTICM (OPTX,SSAX,SLEGX, PATH,RH,-NAER)
! call EXTRAL(OD600,L1U,L2U,N_,JTAUMX,ATAU,ATAU0, JXTRA)
! call X_interp (TTTX,XQO2, TQQ(1,1),QO2(K,1), TQQ(2,1),QO2(K,2),.,,)
! call X_interp (TTTX,XQO3, TQQ(1,2),QO3(K,1), TQQ(2,2),QO3(K,2),,,,)
! call OPMIE (DTAUX,POMEGAX,U0,RFL,AMF2,JXTRA,
! FJACT,FJTOP,FJBOT,FSBOT,FJFLX,FLXD,FLXD0, LU)
! call JRATET(PPJ,TTJ,FFF, VALJXX, LU,NJXU)
! call JP_ATM(PPJ,TTJ,DDJ,OOJ,ZZJ,DTAU600,POMG600,JXTRA, LU)
! subroutine OPMIE (DTAUX,POMEGAX,U0,RFL,AMF2,JXTRA,
! FJACT,FJTOP,FJBOT,FSBOT,FJFLX,FLXD,FLXD0, LU)
! call MIESCT(FJ,FJT,FJB,POMEGA,FZ,ZTAU,ZFLUX,RFL,U0,ND)
! subroutine MIESCT(FJ,FJT,FJB, POMEGA,FZ,ZTAU,ZFLUX,RFL,U0,ND)
! call LEGND0 (EMU(I),PM0,M2_)
! call LEGND0 (-U0,PM0,M2_)
! call BLKSLV(FJ,POMEGA,FZ,ZTAU,ZFLUX,RFL,PM,PM0,FJT,FJB, ND)
! subroutine LEGND0 (X,PL,N)
! subroutine BLKSLV(FJ,POMEGA,FZ,ZTAU,ZFLUX,RFL,PM,PM0,FJTOP,FJBOT,ND)
! call GEN_ID (POMEGA(1,1,K),FZ(1,K),ZTAU(1,K),ZFLUX(K),RFL(K),PM,PM0,
! B(1,1,1,K),CC(1,1,1,K),AA(1,1,1,K),A(1,1,K),H(1,1,K),C(1,1,K),ND)
! subroutine GEN_ID(POMEGA,FZ,ZTAU,ZFLUX,RFL,PM,PM0,B,CC,AA,A,H,C,ND)
! subroutine OPTICL (OPTD,SSALB,SLEG, ODCLD,NDCLD)
! subroutine OPTICA (OPTD,SSALB,SLEG, PATH,RELH,L)
! subroutine OPTICM (OPTD,SSALB,SLEG, PATH,RELH,L)
! subroutine JRATET(PPJ,TTJ,FFF, VALJL,LU,NJXU)
! call X_interp (TT,QO2TOT, TQQ(1,1),QO2(K,1),...)
! call X_interp (TT,QO3TOT, TQQ(1,2),QO3(K,1),...)
! call X_interp (TT,QO31DY, TQQ(1,3),Q1D(K,1),...)
! call X_interp (PP,QQQT, TQQ(1,J),QQQ(K,1,J),...)
! call X_interp (TT,QQQT, TQQ(1,J),QQQ(K,1,J),...)
! subroutine X_interp (TINT,XINT, T1,X1, T2,X2, T3,X3, L123)
! subroutine JP_ATM(PPJ,TTJ,DDJ,OOJ,ZZJ,DTAU6,POMEG6,JXTRA,LU)
! >>> subroutine JP_ATM0(PPJ,TTJ,DDJ,OOJ,ZZJ, LU)
! subroutine SPHERE2(U0,RAD,ZHL,ZZHT,AMF2, L1U,LJX1U)
! subroutine EXTRAL(DTAUX,L1X,L2X,NX,JTAUMX,ATAU,ATAU0, JXTRA)
! note - the calls to EXITC are not listed here
! e.g., call EXITC(' INIT_JX: invalid no. wavelengths')
!
!
!
! >>>>>>>>>>>>>>>>current code revised to JX ver 7.0+ (10/12)<<<<<<<<<<<<
!
! fastJX version 7.0+ (f90) - Prather notes (Jan 2013)
!
!---calculation of cloud optical depth in FJX-70b !!!
!--- assumes that clouds are 100% if in layer
!
! IWP = ice water path (in layer, in cloud) in g/m**2
! LWP = liquid water path (in layer, in cloud) in g/m**2
! REFFI = effective radius of ice cloud (microns)
! REFFL = effective radius of liquid cloud (microns)
!
!>>>>method for calculating cloud OD (in layer) used by FJX core or outside
!>>>>FJX core needs only the _WP and the REFF_
!>>>> note that FJX can use correct Q's and density updates if need be.
! ODL = LWP * 0.75 * 2.10 / REFFL
! ODI = IWP * 0.75 * 2.00 / (REFFI * 0.917)
!
!>>>R-effective determined by main code, not FJX
! REFF determined by user - some recommendations below (from Neu & Prather)
! REFFI is a simple function of ice water content IWC (g/m3, 0.0001 to 0.1)
! IWC = IWP / delta-Z (of layer in m, approx OK)
! REFFI = 50. * (1. + 8.333 * IWC)
! prefer Heymsfield++ 2003 JAM, log-log fit ext(/m) vs. IWC, Fig B1a, p.1389
! EXT (/m) = 1.7e-3 * (IWC/0.1)**0.77
! REFFI = 164. * IWC**0.23 (33 microns at 0.001 --- 164 at 1.0)
!
! REFFL is a simple function of pressure (PCLD):
! FACTOR = (PCLD - 610.) / 200.
! FACTOR = min(1.0, max(0.0, FACTOR))
! REFFL = 9.60*FACTOR + 12.68*(1.-FACTOR)
!
!>>>indices for cloud scattering determined by FJX core, not main code.
! NDX = cloud scattering index based on REFFL or TCLD(for ice cloud)
! NDXI = 6 ! ice hexag (cold)
! if (TCLD .ge. 233.15) then
! NDXI = 7 ! ice irreg
! end if
! NDXC = 1
! do I=2,4
! if (REFFL .gt. 0.5*(RCC(I-1)+RCC(I))) then
! NDXC = I
! end if
! end do
!
!
! version 7.0
! New modular structure:
! Designed for .f90 and CAM5
! All set up routines grouped into one super-call
! Separate interface with CTM and ICAs calls PHOT_J
! PHOT_JX only sees a single column atmosphere (or ICA) thru calling params
!
! version 6.8
! New layout and formatting of FJX_spec.dat (required)
! allows 1,2, or 3 T's (or P's) for all X-sects
! VOCs mostly use Pressure interp.
!
! version 6.7 (should output same J-values as ver 6.6 EXCEPT for VOCs) (3/12)
! New JPL-2010 update complete through all VOCs (see notes in FJX_spec.dat)
! Most important change is incorporation of Stern-Volmer pressure-dependent
! and wavelength-dependent quantum yields into the Temperature interpolation.
! Acetone now has only one pair of X sections for each pathway.
! Redo of mapping of (external) chemical model J's onto fastJX J's
! see examples of splitting a J with fixed branching q-yields
! see how chem model labels are not tied by the cross section labels in fastJX
! Changes in FX_spec.dat make it incompatible with earlier versions, although
! all of the Xsection data have the same format.
! Now the number of X sections read in equals the number of J's that fast-JX calculates.
! As before, all J's except for O2, O3, O3(1D) have only pairs of data at different T.
!
! version 6.6x
! N.B. SPECIAL FIX of tropospheric bin 11 used in FULL(W_=18) and TROP-ONLY(W_=12)
! Attenuation of this bin is too weak since it mixes 218 nm and 288 nm
! Thus it allows a low level of J-O2 (and related 218 nm Xsects) thru-out trop.
! The attenuation in strat (1e-4) is OK, but must zero out in trop (>100 hPa)
! The TROP-QUICK (W_=8) does not use bin#11 and is unaffected.
! Redo of the 4x4 matrix inversions from subroutine to inline code (faster)
! FJX_spec.dat: includes new JPL2010 and solar flux, but J-VOC unchanged (a to-do).
! The J-O2 and O3 steady-state are improved, but the NOy:N2O is too small.
! The J-NO appears too large based on both photocomp2008 & NOy values.
! Possibly lack of thermopsheric NO absorption, or correlation with S-R bands.
! ACTION: J-NO X-sections multiplied by 0.6 in new FJX_spec.dat
!
! version 6.5
! J-values are now averaged over CTM layer (and not just mid-layer point)
! This is important when cloud optical depth is large (~1).
!
! version 6.4
! allows for shortened, speeded up troposphere versions<<<<<<
! STD: W_=18
! identical results to v-6.2 if cloud OD is consistent
! TROP-ONLY: W_=12
! collapses the wavelength bins from 18 to 12 (5-6-7-8 & 11-18)
! drops many 'stratospheric' cross-sections (denoted by 'x' in 2nd title)
! allows use of single standard spectral data set: FJX_spec.dat
! results close to W_=18, largest difference is J-O2 (<1% in 13-18 km!!)
! This is recommended as accurate for troposphere only calculations.
! TROP-QUICK: W_=8
! reverts to original fast-J 7-bins (12-18) plus 1 scaled UV (5) for J-O2
! errors in 12-18 km range for J-O2, high sun are 10%, worse above.
! ***Photolysis of O2 in the upper tropical troposphere is an important
! source of O3. It needs to be included in tropospheric runs.
! TROP-ONLY is recommended, W_=8 is a quick fix if speed essential.
!
! Major rewrite of code to minimize calls and allow better vector-type ops.
! loop over wavelengths internal to Mie soln.
! Driven by profiling of CTM code, may still need optimization.
! Wavelengths can be contracted to W_=12 (trop only) and strat-only
! X-sections are dropped. With parm W_=18, the std fast-JX is retrieved.
! Many call eliminated and summations in BLKSLV and GEN_ID are explicit
! GEN_ID replaces GEN and calculates all matrix coeff's (1:L_) at once
! RD_XXX changed to collapse wavelengths & x-sections to Trop-only:
! WX_ = 18 (parm_CTM.f) should match the JX_spec.dat wavelengths
! W_ = 12 (Trop-only) or 18 (std) is set in (parm_MIE.f).
! if W_=12 then drop strat wavels, and drop x-sects (e.g. N2O, ...)
!
! version 6.3
! revise cloud/aerosol OD & wavelength properties for CTM link:
! OPTICL is new sub for cloud optical properties, but it
! now starts with cloud OD @ 600 nm and cloud NDX
! fast-JX now uses cloud NDX to scale OD to other wavelengths
! OPTICA & OPTICM are new subs to convert aerosol path (g/m2) to OD
! A is std UCI scat data
! M is U Michigan data tables for aerosols, includes Rel Hum effect
! drop sub GAUSSP and put into Parameter statement (parm_MIE.f)
!
! version 6.2
! corrects a long-standing problem at SZA > 89 degrees.
! In prior versions the ray-tracing of the path (and air-mass functions)
! back to the sun was done at the edges of the CTM layers (it was developed
! for the grid-point J-value code at Harvard/GISS/UCI). This left the
! interpolation to the mid-layer (needed for J's) open. The prior method
! gave irregular fluctuations in the direct solar beam at mid-layer for
! large SZA > 88. This is now corrected with exact ray-tracing from
! the mid-pt of each CTM layer. For small SZA, there is no effective
! difference, for large SZA, results could be erratic.
! v-6.2 fix should be easy if you have migrated to v6.1, else some minor
! caution may be needed:
! replace sub SPHERE with SPHERE2, AMF2 report factors for mid and egdes.
! replace sub OPMIE with new OPMIE, this uses the new AMF2 correctly.
! replace sub PHOTOJ with new PHOTOJ, this just hands off AMF2 from
! SPHERE2 to OPMIE.
!
! version 6.1 adds
! 6.1b simplifies calling sequences feeds solar factor, albedo, to PHOTOJ
! and read LAT, LNG directly. No substantive changes.
! new read-in of scat data for clouds/aerosols to allow for UMich data
! This has required substantial rewrite of some of the core subroutines:
! OPMIE is now called for each wavelength and without aersol/cloud data
! all subs below OPMIE are unchanged
! OPTICD & OPTICM are new subs to convert path (g/m2) to OD and phase fn
! D is std UCI scat data (re-ordered for clouds 1st)
! M is U Michigan data tables for aerosols, includes Rel Hum effect
! PHOTOJ now assembles the aerosol data (better for CTM implementation)
! This version can reproduce earlier versions exactly, but the test input
! is changed from OD and NDX to PATH (g/m2) and NDX.
!
! version 6.0 adds
! new 200-nm scattering data so that stratospheric aerosols can be done!
! version 5.7
!>>>>>>>>>>>>>>>>current code revised to JX ver 7.1 (10/13)<<<<<<<<<<<<
!------------------------------------------------------------------------------
! 'FJX_CMN_MOD.f90' for fast-JX code v 7.0+ (prather 9/12) !
!------------------------------------------------------------------------------
!
! NB - ALL of these common variables are set paramters,
! They are NOT to be used as variables for a local solution
! Thus this entire set is 'in' only after it is initialized
!-----------------------------------------------------------------------
!
! !DESCRIPTION: FJX_CMN contains fast-JX variables
!
!
! !INTERFACE:
!
MODULE FJX_CMN_MOD
IMPLICIT NONE
PUBLIC
!-----------------------------------------------------------------------
! JXL_: vertical(levels) dim for J-values computed within fast-JX
INTEGER, PARAMETER :: JXL_=100, JXL1_=JXL_+1
! JXL2_: 2*JXL_ + 2 = mx no. levels in the basic Fast-JX grid (mid-level)
INTEGER, PARAMETER :: JXL2_=2*JXL_+2
! WX_ = dim = no. of wavelengths in input file
INTEGER, PARAMETER :: WX_=18
! X_ = dim = max no. of X-section data sets (input data)
INTEGER, PARAMETER :: X_=72
! A_ = dim = no. of Aerosol/cloud Mie sets (input data)
INTEGER, PARAMETER :: A_=40
! C_ = dim = no. of cld-data sets (input data)
INTEGER, PARAMETER :: C_=16
! W_ = dim = no. of Wavelength bins: =18 std, =12 trop only
INTEGER, PARAMETER :: W_=18 ! W_=12 or 18
! N_ = no. of levels in Mie scattering arrays
! = 2*NC+1 = 4*(L_+1) + 1`+ 2*sum(JADDLV)
INTEGER, PARAMETER :: N_=601
! M_ = no. of Gauss points used, must = 4 in fast_JX (no option)
INTEGER, PARAMETER :: M_=4
! M2_ = 2*M_ = 8, replaces MFIT
INTEGER, PARAMETER :: M2_=2*M_
!-----------------------------------------------------------------------
! 4 Gauss pts = 8-stream
REAL*8, DIMENSION(M_), PARAMETER :: &
EMU = [.06943184420297d0, .33000947820757d0, &
.66999052179243d0, .93056815579703d0]
REAL*8, DIMENSION(M_), PARAMETER :: &
WT = [.17392742256873d0, .32607257743127d0, &
.32607257743127d0, .17392742256873d0]
!-----------------------------------------------------------------------
! ZZHT: scale height (cm)
REAL*8, PARAMETER :: ZZHT = 5.d5
! RAD: Radius of Earth (cm)
REAL*8, PARAMETER :: RAD = 6375.d5
! ATAU: heating rate (factor increase from one layer to the next)
REAL*8, PARAMETER :: ATAU = 1.120d0
! ATAU0: minimum heating rate
REAL*8, PARAMETER :: ATAU0 = 0.010d0
! JTAUMX = maximum number of divisions (i.e., may not get to ATAUMN)
INTEGER, PARAMETER :: JTAUMX = (N_ - 4*JXL_)/2
!---- Variables in file 'FJX_spec.dat' (RD_XXX)
! WBIN: Boundaries of wavelength bins
REAL*8 WBIN(WX_+1)
! WL: Centres of wavelength bins - 'effective wavelength'
REAL*8 WL(WX_)
! FL: Solar flux incident on top of atmosphere (cm-2.s-1)
REAL*8 FL(WX_)
REAL*8 QO2(WX_,3) ! QO2: O2 cross-sections
REAL*8 QO3(WX_,3) ! QO3: O3 cross-sections
REAL*8 Q1D(WX_,3) ! Q1D: O3 => O(1D) quantum yield
! QQQ: Supplied cross sections in each wavelength bin (cm2)
REAL*8 QQQ(WX_,3,X_)
! QRAYL: Rayleigh parameters (effective cross-section) (cm2)
REAL*8 QRAYL(WX_+1)
! TQQ: Temperature for supplied cross sections
REAL*8 TQQ(3,X_)
! LQQ = 1, 2, or 3 to determine interpolation with T or P
INTEGER LQQ(X_)
! TITLEJX: Title for supplied cross sections, from 'FJX_spec.dat'
CHARACTER*6 TITLEJX(X_)
! SQQ: Flag for supplied cross sections, from 'FJX_spec.dat'
CHARACTER*1 SQQ(X_)
!---- Variables in file 'FJX_scat-aer.dat' (RD_MIE)
! QAA: Aerosol scattering phase functions
REAL*8 QAA(5,A_)
! WAA: 5 Wavelengths for the supplied phase functions
REAL*8 WAA(5,A_)
! PAA: Phase function: first 8 terms of expansion
REAL*8 PAA(8,5,A_)
! RAA: Effective radius associated with aerosol type
REAL*8 RAA(A_)
! SAA: Single scattering albedo
REAL*8 SAA(5,A_)
! DAA: density (g/cm^3)
REAL*8 DAA(A_)
! NAA: Number of categories for scattering phase functions
INTEGER NAA
!---- Variables in file 'FJX_scat-cld.dat' (RD_CLD)
! QCC: Cloud scattering phase functions
REAL*8 QCC(5,C_)
! WCC: 5 Wavelengths for supplied phase functions
REAL*8 WCC(5,C_)
! PCC: Phase function: first 8 terms of expansion
REAL*8 PCC(8,5,C_)
! RCC: Effective radius associated with cloud type
REAL*8 RCC(C_)
! SCC: Single scattering albedo
REAL*8 SCC(5,C_)
! DCC: density (g/cm^3)
REAL*8 DCC(C_)
! NCC: Number of categories for cloud scattering phase functions
INTEGER NCC
!---- Variables in file 'FJX_scat-UMa.dat' (RD_CLD)
! WMM: U Michigan aerosol wavelengths
REAL*8 WMM(6)
! UMAER: U Michigan aerosol data sets
REAL*8 UMAER(3,6,21,33)
!---- Variables in file 'atmos_std.dat' (RD_PROF)
! T and O3 reference profiles
REAL*8, DIMENSION(51,18,12) :: TREF, OREF
INTEGER NJX,NW1,NW2
END MODULE FJX_CMN_MOD
!------------------------------------------------------------------------------
! 'cmn_CTM.f' for fast-JX code v7.1 !
!------------------------------------------------------------------------------
!
! !MODULE: FJX
!
! !DESCRIPTION: JX version 7.1 (10/13) works with FL=0 for <200 nm
!
!
! !INTERFACE:
!
MODULE FJX_SUB_MOD
!
! !USES:
!
USE FJX_CMN_MOD
IMPLICIT NONE
PUBLIC
!-----------------------------------------------------------------------
integer, parameter :: &
L_CTM=57, L1_CTM=L_CTM+1 & ! L_ = number of CTM layers
,L_=67, L1_=L_+1 & ! L_ = max number of layers EMEP+stratosphere
!The results should be independent of the values of L_
! ,L2_=2*L_+2 & ! no. levels in the Fast-JX grid that
! includes both layer edges and layer mid-points
,JVL_=57 & ! vertical(levels) dim for J-values sent to CTM
,JVN_=101 & ! max no. of J-values
,AN_=25 ! no of separate aerosols per layer (needs NDX for each)
!-----------------------------------------------------------------------
! variables used to map fast-JX J's onto CTM J's
!-----------------------------------------------------------------------
integer::debug
real*8 JFACTA(JVN_) ! multiplication factor for fast-JX calculated J
integer JIND(JVN_) ! index arrays that map Jvalue(j) onto rates
integer NRATJ ! number of Photolysis reactions in CTM chemistry,
! derived here NRATJ must be .le. JVN_
character*50 JLABEL(JVN_) ! label of J-value used in the main chem model
!
! !PUBLIC SUBROUTINES:
!
PUBLIC :: EXITC, SOLAR_JX, ACLIM_FJX, JP_ATM0, PHOTO_JX
CONTAINS
!-----------------------------------------------------------------------
subroutine SOLAR_JX(GMTIME,NDAY,YGRDJ,XGRDI, SZA,COSSZA,SOLFX)
!-----------------------------------------------------------------------
! GMTIME = UT for when J-values are wanted
! (for implicit solver this is at the end of the time step)
! NDAY = integer day of the year (used for solar lat and declin)
! YGRDJ = laitude (radians) for grid (I,J)
! XGDRI = longitude (radians) for grid (I,J)
!
! SZA = solar zenith angle in degrees
! COSSZA = U0 = cos(SZA)
!-----------------------------------------------------------------------
implicit none
real*8, intent(in) :: GMTIME,YGRDJ,XGRDI
integer, intent(in) :: NDAY
real*8, intent(out) :: SZA,COSSZA,SOLFX
!
real*8 PI, PI180, LOCT
real*8 SINDEC, SOLDEK, COSDEC, SINLAT, SOLLAT, COSLAT, COSZ
!
PI = 3.141592653589793d0
PI180 = PI/180.d0
SINDEC = 0.3978d0*sin(0.9863d0*(dble(NDAY)-80.d0)*PI180)
SOLDEK = asin(SINDEC)
COSDEC = cos(SOLDEK)
SINLAT = sin(YGRDJ)
SOLLAT = asin(SINLAT)
COSLAT = cos(SOLLAT)
!
LOCT = (((GMTIME)*15.d0)-180.d0)*PI180 + XGRDI
COSSZA = COSDEC*COSLAT*cos(LOCT) + SINDEC*SINLAT
SZA = acos(COSSZA)/PI180
!
! write(6,*) ' XGRDI,YGRDJ',XGRDI,YGRDJ
! write(6,*) ' LOCT (rad)',LOCT
! write(6,*) ' SINDEC,COSDEC', SINDEC,COSDEC
! write(6,*) ' SINLAT,COSLAT', SINLAT,COSLAT
! write(6,*) ' COS, SZA',COSSZA,SZA
!
SOLFX = 1.d0-(0.034d0*cos(dble(NDAY-186)*2.d0*PI/365.d0))
!
END SUBROUTINE SOLAR_JX
!-----------------------------------------------------------------------
subroutine ACLIM_FJX (YLATD, MONTH, PPP,TTT,ZZZ,DDD,OOO, L1U)
!-----------------------------------------------------------------------
! Load fast-JX climatology for latitude & month given pressure grid
!-----------------------------------------------------------------------
implicit none
real*8, intent(in) :: YLATD
integer, intent(in) :: MONTH, L1U
real*8, intent(in), dimension(L1U+1) :: PPP
real*8, intent(out), dimension(L1U+1) :: ZZZ
real*8, intent(out), dimension(L1U) :: TTT,DDD,OOO
real*8, dimension(51) :: OREF2,TREF2
real*8, dimension(52) :: PSTD
integer K, L, M, N
real*8 DLOGP,F0,T0,PB,PC,XC,MASFAC,SCALEH
! Select appropriate month
M = max(1,min(12,MONTH))
! Select appropriate latitudinal profiles
N = max(1, min(18, (int(YLATD+99)/10 )))
do K = 1,51
OREF2(K) = OREF(K,N,M)
TREF2(K) = TREF(K,N,M)
end do
! Apportion O3 and T on supplied climatology z levels onto CTM levels +1
! with mass (pressure) weighting, assuming constant mixing ratio and
! temperature half a layer on either side of the point supplied.
! PPP(L=1:L1_)=edge-pressure of CTM layer, PPP(L1_+1)=0 (top-of-atmos)
! Mass factor - delta-Pressure (mbars) to delta-Column (molecules.cm-2)
MASFAC = 100.d0*6.022d+23/(28.97d0*9.8d0*10.d0)
! Set up pressure levels for O3/T climatology - assume that value
! given for each 2 km z* level applies from 1 km below to 1 km above,
! so select pressures at these boundaries. Surface level values at
! 1000 mb are assumed to extend down to the actual PSURF (if > 1000)
PSTD(1) = max(PPP(1),1000.d0)
PSTD(2) = 1000.d0 * 10.d0**(-1.d0/16.d0)
DLOGP = 10.d0**(-2.d0/16.d0)
do K = 3,51
PSTD(K) = PSTD(K-1)*DLOGP
end do
PSTD(52) = 0.d0
do L = 1,L1U
F0 = 0.d0
T0 = 0.d0
do K = 1,51
PC = min(PPP(L),PSTD(K))
PB = max(PPP(L+1),PSTD(K+1))
if (PC .gt. PB) then
XC = (PC-PB)/(PPP(L)-PPP(L+1))
F0 = F0 + OREF2(K)*XC
T0 = T0 + TREF2(K)*XC
end if
end do
TTT(L) = T0
DDD(L) = (PPP(L)-PPP(L+1))*MASFAC
OOO(L) = F0*1.d-6*DDD(L)
end do
! Calculate effective altitudes using scale height at each level
ZZZ(1) = 0.d0
do L = 1,L1U-1
SCALEH = 1.3806d-19*MASFAC*TTT(L)
ZZZ(L+1) = ZZZ(L) -( LOG(PPP(L+1)/PPP(L)) * SCALEH )
end do
ZZZ(L1U+1) = ZZZ(L1U) + ZZHT
END SUBROUTINE ACLIM_FJX
!<<<<<<<<<<<<<<<<<<<end CTM-fastJX special call subroutines<<<<<<<<<<<<<
!<<<<<<<<<<<<<<<<<<<<<<<<begin fastJX core code<<<<<<<<<<<<<<<<<<<<<<<<<
!<<<<<<<<<<<<<<all outside calls go through PHOTO_JX<<<<<<<<<<<<<<<<<<<<
!-----------------------------------------------------------------------
! fastJX version 7.2 (f90) - Prather notes (Jan 2013)
!---calculation of cloud optical depth in FJX-72 !!!
!--- assumes that clouds are 100% if in layer
! IWP = ice water path (in layer, in cloud) in g/m**2
! LWP = liquid water path (in layer, in cloud) in g/m**2
! REFFI = effective radius of ice cloud (microns)
! REFFL = effective radius of liquid cloud (microns)
!>>>>method for calculating cloud OD (in layer) used by FJX core or outside
!>>>>FJX core needs only the _WP and the REFF_
!>>>> note that FJX can use correct Q's and density updates if need be.
! ODL = LWP * 0.75 * 2.10 / REFFL
! ODI = IWP * 0.75 * 2.00 / (REFFI * 0.917)
!>>>R-effective determined by main code, not FJX
! REFF determined by user - some recommendations below (from Neu & Prather)
! REFFI is a simplle function of ice water content IWC (g/m3, 0.0001 to 0.1)
! IWC = IWP / delta-Z (of layer in m, approx OK)
! REFFI = 50. * (1. + 8.333 * IWC)
! prefer Heymsfield++ 2003 JAM, log-log fit ext(/m) vs. IWC, Fig B1a, p.1389
! EXT (/m) = 1.7e-3 * (IWC/0.1)**0.77
! REFFI = 164. * IWC**0.23 (33 microns at 0.001 --- 164 at 1.0)
! REFFL is a simple function of pressure (PCLD):
! FACTOR = (PCLD - 610.) / 200.
! FACTOR = min(1.0, max(0.0, FACTOR))
! REFFL = 9.60*FACTOR + 12.68*(1.-FACTOR)
!>>>indices for cloud scattering determined by FJX core, not main code.
! NDX = cloud scattering index based on REFFL or TCLD(for ice cloud)
! NDXI = 6 ! ice hexag (cold)
! if (TCLD .ge. 233.15) then
! NDXI = 7 ! ice irreg
! end if
! NDXC = 1
! do I=2,4
! if (REFFL .gt. 0.5*(RCC(I-1)+RCC(I))) then
! NDXC = I
! end if
! end do
!>>>>>revise PHOTO_JX to drop NDXCLD, and to include
! LWP, REFFL, IWP, REFFI only, use T, etc to get NDX.
!-----------------------------------------------------------------------
subroutine PHOTO_JX &
(U0,SZA,REFLB,SOLF, LPRTJ,PPP,ZZZ,TTT,DDD,RRR,OOO, &
LWP,IWP,REFFL,REFFI, AERSP,NDXAER,L1U,L1_SIZE,ANU, VALJXX,NJXU, FJBOT,FSBOT,debug)
!-----------------------------------------------------------------------
!
! PHOTO_JX is the gateway to fast-JX calculations:
! calc J's for a single column atmosphere (aka Indep Colm Atmos or ICA)
! needs P, T, O3, clds, aersls; adds top-of-atmos layer from climatology
! needs day-fo-year for sun distance, SZA (not lat or long)
!-----------------------------------------------------------------------
implicit none
!---calling sequence variables
integer, intent(in) :: L1U,L1_SIZE,ANU,NJXU,debug
real*8, intent(in) :: U0,SZA,REFLB,SOLF
logical, intent(in) :: LPRTJ
real*8, intent(in),dimension(L1_SIZE+1) :: PPP,ZZZ
real*8, intent(in),dimension(L1_SIZE ) :: TTT,DDD,RRR,OOO
real*8, intent(in),dimension(L1_SIZE ) :: LWP,IWP,REFFL,REFFI
real*8, intent(in),dimension(L1_SIZE,ANU):: AERSP
integer, intent(in),dimension(L1_SIZE,ANU):: NDXAER
!---reports out the JX J-values, upper level program converts to CTM chemistry J's
real*8, intent(out), dimension(L1_SIZE-1,NJXU):: VALJXX
real*8, intent(out), dimension(W_) :: FJBOT,FSBOT
!-----------------------------------------------------------------------
!--------key LOCAL atmospheric data needed to solve plane-parallel J----
!-----these are dimensioned JXL_, and must have JXL_ .ge. L_
real*8, dimension(JXL1_+1) :: TTJ,DDJ,OOJ,ZZJ
real*8, dimension(JXL1_+1) :: PPJ,RELH
integer,dimension(JXL2_+1) :: JXTRA
real*8, dimension(W_) :: FJTOP,FLXD0,RFL
real*8, dimension(JXL_, W_) :: AVGF, FJFLX
real*8, dimension(JXL1_,W_) :: DTAUX, FLXD
real*8, dimension(8,JXL1_,W_) :: POMEGAX
real*8, dimension(JXL1_) :: DTAU600
real*8, dimension(8,JXL1_) :: POMG600
real*8, dimension(W_,JXL1_) :: FFX
real*8, dimension(W_,8) :: FFXNET
!---flux/heating arrays (along with FJFLX,FLXD,FLXD0)
real*8 FLXJ(JXL1_),FFX0,FXBOT,FABOT
real*8 :: ODABS,ODRAY,ODI,ODL
real*8 :: RFLECT,FREFS,FREFL,FREFI
real*8 :: AMF2(2*JXL1_+1,2*JXL1_+1)
!------------key SCATTERING arrays for clouds+aerosols------------------
real*8 :: OPTX(5),SSAX(5),SLEGX(8,5)
real*8 :: OD(5,JXL1_),SSA(5,JXL1_),SLEG(8,5,JXL1_)
real*8 :: OD600(JXL1_)
real*8 :: PATH,RH,XTINCT
!------------key arrays AFTER solving for J's---------------------------
real*8 FFF(W_,JXL_),VALJ(X_)
real*8 FLXUP(W_),FLXDN(W_),DIRUP(W_),DIRDN(W_)
real*8 VALJL(JXL_,X_) !2-D array of J_s returned by JRATET
integer LU,L2U, I,J,K,L,M,KMIE,KW,NAER,NDXI,NDXL, RATIO(W_)
real*8 XQO3,XQO2,DTAUC,WAVE, TTTX
!-----------------------------------------------------------------------
if (L1U .gt. JXL1_) then
call EXITC(' PHOTO_JX: not enough levels in JX')
end if
LU = L1U - 1
L2U = LU + LU + 2
FFF(:,:) = 0.d0
FREFI = 0.d0
FREFL = 0.d0
FREFS = 0.d0
!---check for dark conditions SZA > 98.0 deg => tan ht = 63 km
! or 99.0 80 km
if (SZA .gt. 98.d0)then
VALJXX=0.0
goto 99
end if
!---load the amtospheric column data
do L = 1,L1U
PPJ(L) = PPP(L)
TTJ(L) = TTT(L)
DDJ(L) = DDD(L)
OOJ(L) = OOO(L)
end do
PPJ(L1U+1) = 0.d0
!---calculate spherical weighting functions (AMF: Air Mass Factor)
do L = 1,L1U+1
ZZJ(L) = ZZZ(L)
end do
RFLECT = REFLB
!-----------------------------------------------------------------------
call SPHERE2 (U0,RAD,ZZJ,ZZHT,AMF2, L1U,JXL1_)
!-----------------------------------------------------------------------
!---calculate the optical properties (opt-depth, single-scat-alb, phase-fn(1:8))
!--- at the 5 std wavelengths 200-300-400-600-999 nm for cloud+aerosols
do L = 1,L1U
do K=1,5
OD(K,L) = 0.d0
SSA(K,L) = 0.d0
do I=1,8
SLEG(I,K,L) = 0.d0
end do
end do
end do
do L = 1,L1U
!---Liquid Water Cloud
if (LWP(L) .gt. 1.d-5 .and. REFFL(L) .gt. 0.1d0) then
!---extinction K(m2/g) = 3/4 * Q / [Reff(micron) * density(g/cm3)]
ODL = LWP(L) * 0.75d0 * 2.1d0 / REFFL(L)
NDXL = 1
do I=2,4
if (REFFL(L) .gt. 0.5*(RCC(I-1)+RCC(I))) then
NDXL = I
end if
end do
call OPTICL (OPTX,SSAX,SLEGX, ODL,NDXL)
do K=1,5
OD(K,L) = OD(K,L) + OPTX(K)
SSA(K,L) = SSA(K,L) + SSAX(K)*OPTX(K)
do I=1,8
SLEG(I,K,L)=SLEG(I,K,L) + SLEGX(I,K)*SSAX(K)*OPTX(K)
end do
end do
!>>>diagnostic print of cloud data:
! write(6,'(a,i3,2f8.2,f8.4,f8.2,f8.4,i4)') &
! 'Liq Cld',L,PPP(L),PPP(L+1),LWP(L),REFFL(L),ODL,NDXL
end if
!---Ice Water Cloud
if (IWP(L) .gt. 1.d-5 .and. REFFI(L) .gt. 0.1d0) then
ODI = IWP(L) * 0.75d0 * 2.0d0 / (REFFI(L) * 0.917d0)
if (TTT(L) .ge. 233.15d0) then
NDXI = 7 ! ice irreg
else
NDXI = 6 ! ice hexag (cold)
end if
call OPTICL (OPTX,SSAX,SLEGX, ODI,NDXI)
do K=1,5
OD(K,L) = OD(K,L) + OPTX(K)
SSA(K,L) = SSA(K,L) + SSAX(K)*OPTX(K)
do I=1,8
SLEG(I,K,L)=SLEG(I,K,L) + SLEGX(I,K)*SSAX(K)*OPTX(K)
end do
end do
!>>>diagnostic print of cloud data:
! write(6,'(a,i3,2f8.2,f8.4,f8.2,f8.4,i4)') &
! 'Ice Cld',L,PPP(L),PPP(L+1),IWP(L),REFFI(L),ODI,NDXI
end if
!---aerosols in layer: check aerosol index
!---this uses data from climatology OR from current CTM (STT of aerosols)
!---FIND useful way to sum over different aerosol types!
do M = 1,ANU
NAER = NDXAER(L,M)
PATH = AERSP(L,M)
RH = RRR(L)
!---subroutines OPTICA & OPTICM return the same information:
!--- optical depth (OPTX), single-scat albedo (SSAX) and phase fn (SLEGX(8))
!---subs have slightly different inputs:
!--- PATH is the g/m2 in the layer, NAER in the cloud/aerosol index
!--- UMich aerosols use relative humidity (RH)
if (PATH .gt. 0.d0) then
if (NAER .gt.0) then
call OPTICA (OPTX,SSAX,SLEGX, PATH,RH, NAER)
else
call OPTICM (OPTX,SSAX,SLEGX, PATH,RH,-NAER)
end if
do K=1,5
OD(K,L) = OD(K,L) + OPTX(K)
SSA(K,L) = SSA(K,L) + SSAX(K)*OPTX(K)
do I=1,8
SLEG(I,K,L)=SLEG(I,K,L) + SLEGX(I,K)*SSAX(K)*OPTX(K)
end do
end do
end if
end do
do K=1,5
if (OD(K,L) .gt. 0.d0) then
SSA(K,L) = SSA(K,L)/OD(K,L)
do I=1,8
SLEG(I,K,L) = SLEG(I,K,L)/OD(K,L)
end do
end if
end do
!---Include aerosol with cloud OD at 600 nm to determine added layers:
OD600(L) = OD(4,L)
end do
!---when combining with Rayleigh and O2-O3 abs, remember the SSA and
!--- phase fn SLEG are weighted by OD and OD*SSA, respectively.
!---Given the aerosol+cloud OD/layer in visible (600 nm) calculate how to add
! additonal levels at top of clouds (now uses log spacing)
!-----------------------------------------------------------------------
call EXTRAL(OD600,L1U,L2U,N_,JTAUMX,ATAU,ATAU0, JXTRA)
!-----------------------------------------------------------------------
!---set surface reflectance
RFL(:) = max(0.d0,min(1.d0,RFLECT))
!-----------------------------------------------------------------------
!---Loop over all wavelength bins to calc mean actinic flux AVGF(L)
!-----------------------------------------------------------------------
do K = 1,W_
if (FL(K) .gt. 1.d0) then
WAVE = WL(K)
!---Pick nearest Mie wavelength to get scattering properites------------
KMIE=1 ! use 200 nm prop for <255 nm
if( WAVE .gt. 255.d0 ) KMIE=2 ! use 300 nm prop for 255-355 nm
if( WAVE .gt. 355.d0 ) KMIE=3 ! use 400 nm prop for 355-500 nm
if( WAVE .gt. 500.d0 ) KMIE=4
if( WAVE .gt. 800.d0 ) KMIE=5
!---Combine: Rayleigh scatters & O2 & O3 absorbers to get optical properties
!---values at L1_=L_+1 are a pseudo/climatol layer above the top CTM layer (L_)
do L = 1,L1U
TTTX = TTJ(L)
call X_interp (TTTX,XQO2, TQQ(1,1),QO2(K,1), TQQ(2,1),QO2(K,2), &
TQQ(3,1),QO2(K,3), LQQ(1))
call X_interp (TTTX,XQO3, TQQ(1,2),QO3(K,1), TQQ(2,2),QO3(K,2), &
TQQ(3,2),QO3(K,3), LQQ(2))
ODABS = XQO3*OOJ(L) + XQO2*DDJ(L)*0.20948d0
ODRAY = DDJ(L)*QRAYL(K)
DTAUX(L,K) = OD(KMIE,L) + ODABS + ODRAY
do I=1,8
POMEGAX(I,L,K) = SLEG(I,KMIE,L)*OD(KMIE,L)
end do
POMEGAX(1,L,K) = POMEGAX(1,L,K) + 1.0d0*ODRAY
POMEGAX(3,L,K) = POMEGAX(3,L,K) + 0.5d0*ODRAY
do I=1,8
POMEGAX(I,L,K) = POMEGAX(I,L,K)/DTAUX(L,K)
end do
end do
end if
end do
!-----------------------------------------------------------------------
call OPMIE (DTAUX,POMEGAX,U0,RFL,AMF2,JXTRA, &
AVGF,FJTOP,FJBOT,FSBOT,FJFLX,FLXD,FLXD0, LU)
!-----------------------------------------------------------------------
FLXUP(:) = 0.d0
DIRUP(:) = 0.d0
FLXDN(:) = 0.d0
DIRDN(:) = 0.d0
FLXJ(:) = 0.d0
FFX(:,:) = 0.d0
FFXNET(:,:) = 0.d0
do K = 1,W_
if (FL(K) .gt. 1.d0) then
!----direct(DIR) and diffuse(FLX) fluxes at top(UP) (solar = negative by convention)
!---- also at bottom (DN), does not include diffuse reflected flux.
FLXUP(K) = FJTOP(K)
DIRUP(K) = -FLXD0(K)
FLXDN(K) = -FJBOT(K)
DIRDN(K) = -FSBOT(K)
do L = 1,LU
FFF(K,L) = FFF(K,L) + SOLF*FL(K)*AVGF(L,K)
end do
FREFI = FREFI + SOLF*FL(K)*FLXD0(K)/WL(K)
FREFL = FREFL + SOLF*FL(K)*FJTOP(K)/WL(K)
FREFS = FREFS + SOLF*FL(K)/WL(K)
!---for each wavelength calculate the flux budget/heating rates:
! FLXD(L) = direct flux deposited in layer L [approx = MU0*(F(L+1) -F(L))]
! but for spherical atmosphere!
! FJFLX(L) = diffuse flux across top of layer L
!---calculate divergence of diffuse flux in each CTM layer (& t-o-a)
!--- need special fix at top and bottom:
!---FABOT = total abs at L.B. & FXBOT = net diffusive flux at L.B.
FABOT = (1.d0-RFL(K))*(FJBOT(K)+FSBOT(K))
FXBOT = -FJBOT(K) + RFL(K)*(FJBOT(K)+FSBOT(K))
FLXJ(1) = FJFLX(1,K) - FXBOT
do L=2,LU
FLXJ(L) = FJFLX(L,K) - FJFLX(L-1,K)
end do
FLXJ(LU+1) = FJTOP(K) - FJFLX(LU,K)
!---calculate net flux deposited in each CTM layer (direct & diffuse):
FFX0 = 0.d0
do L=1,L1U
FFX(K,L) = FLXD(L,K) - FLXJ(L)
FFX0 = FFX0 + FFX(K,L)
end do
! NB: the radiation level ABOVE the top CTM level is included in these budgets
! these are the flux budget/heating terms for the column:
! FFXNET(K,1) = FLXD0 direct(solar) flux dep into atmos (spherical)
! FFXNET(K,2) = FSBOT direct(solar) flux dep onto LB (surface)
! FFXNET(K,3) = FLXD0+FSBOT TOTAL solar into atmopshere+surface
! FFXNET(K,4) = FJTOP diffuse flux leaving top-of-atmos
! FFXNET(K,5) = FFX0 diffuse flux absorbed in atmos
! FFXNET(K,6) = FABOT total (dir+dif) absorbed at LB (surface)
! these are surface fluxes to compare direct vs. diffuse:
! FFXNET(K,7) = FSBOT direct flux dep onto LB (surface) - for srf diags
! FFXNET(K,8) = FJBOT diffuse flux dep onto LB (surface)
FFXNET(K,1) = FLXD0(K)
FFXNET(K,2) = FSBOT(K)
FFXNET(K,3) = FLXD0(K) + FSBOT(K)
FFXNET(K,4) = FJTOP(K)
FFXNET(K,5) = FFX0
FFXNET(K,6) = FABOT
FFXNET(K,7) = FSBOT(K)
FFXNET(K,8) = FJBOT(K)
!-----------------------------------------------------------------------
end if
end do ! end loop over wavelength K
!-----------------------------------------------------------------------
FREFL = FREFL/FREFS !calculate reflected flux (energy weighted)
FREFI = FREFI/FREFS
!---NB UVB = 280-320 = bins 12:15, UVA = 320-400 = bins 16:17, VIS = bin 18 (++)
!-----------------------------------------------------------------------
call JRATET(PPJ,TTJ,FFF, VALJXX, LU,L1_SIZE-1,NJXU)
!-----------------------------------------------------------------------