-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmodule_ra_kdm.f90
6655 lines (5601 loc) · 333 KB
/
module_ra_kdm.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
!WRF:MODEL_LAYER:PHYSICS
!c-mm!XREF:MEDIATION_LAYER:PHYSICS
!
!-----------------------------------------------------------------------
! This module contains most parameters needed to operate the Hadley !
! radiation scheme. For each variable type, there are two blocks of !
! variables. The first block contains fixed values that should never !
! change. The second block contains values that can be modified !
! according to need. !
!-----------------------------------------------------------------------
MODULE module_kdm_parameters
USE, INTRINSIC :: IEEE_ARITHMETIC ! to catch NaN and Inf
#ifndef mpas
USE module_model_constants
USE module_diag_common, only: test_for_bad, test_for_bad_detailed, test_for_bad_detailed_2d
#endif
IMPLICIT NONE
! -------- hardwired case-independent constants --- DO NOT CHANGE THESE!!!
INTEGER, PARAMETER :: & ! THESE VARIBLES SHOULD REMAIN FIXED
ip_scf_solar_up = 1, & ! Index of source coefficient for upward solar beam
ip_scf_solar_down = 2, & ! Index of source coefficient for downward solar beam
ip_scf_ir_1d = 1, & ! Index of source coefficient for 1st difference of Planckian
ip_scf_ir_2d = 2, & ! Index of source coefficient for 2nd difference of Planckian
npd_source_coeff = 2, & ! Number of coefficients for two-stream sources
ip_solar = 1, & ! Index for solar loop = 1
ip_infra_red = 2, & ! Index for IR loop = 2
ip_clear = 1, & ! Index for loop with clear-sky conditions = 1
ip_dust = 2, & ! Index for loop with dusty conditions = 2
kdm_do_vis = 0, & ! flag name: if (kdm_mode_flag == kdm_do_vis) => we do vis
kdm_do_ir = 1, & ! flag name: if (kdm_mode_flag == kdm_do_ir) => we do ir
n_band_ir = 7, & !c-mm Number of bands in the IR. Hardwired in KDM. Originally in module_driver_constants as max_n_band_ir
n_band_solar = 7 !c-mm Number of bands in the solar. Hardwired in KDM. Originally in module_driver_constants as max_n_band_solar
REAL(KIND(0.d0)), PARAMETER :: & ! THESE VARIABLES SHOULD REMAIN FIXED
tol_machine = 2.220446D-16, & ! Machine tolerance, also equal to EPSILON(x) (intrinsic) where x is virtually any
tol_div = 3.2D+01*tol_machine, & ! conceivable number used in this code. Architecture dependent, but this
sqrt_tol_machine = 1.490116D-08, & ! is a pretty good number for most.
log_tol_machine = -36.0437, &
quad_correct_limit = 6.828D-06 ! Equal to EXP(3.3D-01*log_tol_machine)
! -------- case-dependent constants --- CHANGE THESE IF YOU WILL, BUT ONLY IF YOU KNOW WHAT YOU'RE DOING!!
INTEGER, PARAMETER :: &
npd_band = max(n_band_ir,n_band_solar), & ! Basically equals MAX(n_band_ir,n_band_solar). Used for setting up arrays.
n_tot_bands = n_band_ir+n_band_solar ! Number of bands in k-coefficient array
#if ( WRF_MARS == 1 )
INTEGER, PARAMETER :: & ! THESE VARIABLES CAN BE MODIFIED
npd_species = 2, & ! Number of "gases" in atmospere. One for each variable gas plus one total for fixed gases
npd_aerosol_species = 3 ! Max number of "aerosols" in the atmosphere. Presently set to three. 1=dust, 2=ice, 3=co2ice
#elif ( WRF_TITAN == 1 )
INTEGER, PARAMETER :: & ! THESE VARIABLES CAN BE MODIFIED
npd_species = 2, & ! Number of "gases" in atmospere. One for each variable gas plus one total for fixed gases
npd_aerosol_species = 1 ! Max number of "aerosols" in the atmosphere. Presently set to one. 1=tholin
#else
INTEGER, PARAMETER :: & ! THESE VARIABLES CAN BE MODIFIED
npd_species = 1, & ! Number of "gases" in atmospere. One for each variable gas plus one total for fixed gases
npd_aerosol_species = 1 ! Max number of "aerosols" in the atmosphere. Presently set to three. 1=dust, 2=ice, 3=co2ice
#endif
REAL(KIND(0.d0)), PARAMETER :: & ! THESE VARIABLES SHOULD REMAIN FIXED
#if ( WRF_TITAN == 1 )
p_spacing = 0.3333, &
mix_spacing = 1.0, &
#elif ( WRF_MARS == 1 )
p_spacing = 0.4, &
mix_spacing = 1.0, &
zs_lv = -5.5, & ! Altitude of the NLTE transition region in terms of ALOG(p), where p is in nbar
zw_lv = 0.5, & ! Width of the NLTE transition region. The above value of zs_lv corresponds to
! approximately 82 km for reference atmosphere 1 in Lopez-Valverde.
#else
p_spacing = 0.3333, &
mix_spacing = 1.0, &
#endif
pressure_TOA = 1.d-8, & ! top of atmosphere pressure (in Pa) (machine small but non-zero) 10^-8 is about 300km on Mars
#if ( WRF_MARS == 1 )
pressure_AWDC = 5.d-4, & ! pressure at which to stop adding fake layers (awdc - above which dont care)
fake_scale = 0.85, & ! 0.85 is exp(-dz/H) where (e.g., for Mars 0.85 maps to dz=1.5km for H=10km)
#elif ( WRF_TITAN == 1 )
pressure_AWDC = 2.d-3, & ! pressure at which to stop adding fake layers (awdc - above which dont care)
fake_scale = 0.5, & ! 0.5 is exp(-dz/H) where (e.g., for Titan 0.5 maps to dz=15km for H=30km)
#else
pressure_AWDC = 2.d-4, & ! pressure at which to stop adding fake layers (awdc - above which dont care)
fake_scale = 0.5, & ! 0.5 is exp(-dz/H) where (e.g., for Titan 0.5 maps to dz=15km for H=30km)
#endif
math_pi = pi_d
INTEGER, PARAMETER :: max_fake_layers_visir = 50
LOGICAL :: verbose_debug = .false.
! -------- case-dependent constants that are set in kdminit
REAL :: T_ref, & ! these values are now set in kdminit
T_iso, & ! and pulled in from 1d state variables
lapse_rate
REAL(KIND(0.d0)), DIMENSION(n_band_solar+1) :: band_wn_solar
REAL(KIND(0.d0)), DIMENSION(n_band_ir+1) :: band_wn_ir
INTEGER :: n_temps, & ! Number of temperatures in k-coefficient array
n_press, & ! Number of pressures in k-coefficient array
n_mix, & ! Number of mixing ratios in k-coefficient array
n_terms ! Number of gaussian quadrature terms in each wl band
REAL(KIND(0.d0)), ALLOCATABLE, DIMENSION(:) :: kval_temps, & ! temperature values in kvals array
kval_press, & ! pressure values in kvals array
kval_mix, & ! mixing ratio values in kvals array
weights ! gaussian quadrature weights
REAL(KIND(0.d0)) :: kval_press_min, & ! min pressure value in kvals array
kval_press_max, & ! max pressure value in kvals array
kval_mix_min, & ! min mixing ratio value in kvals array
kval_mix_max, & ! max mixing ratio value in kvals array
kval_temps_min, & ! min temperature in kvals array
kval_temps_max ! max temperature in kvals array
#if ( WRF_MARS == 1 )
integer, parameter :: ko_debug = 0, &
ko_multimie_dust = 1, &
ko_multimie_water = 2, &
ko_multimie_co2 = 3, &
ko_dummy04 = 4, &
ko_dummy05 = 5, &
ko_dummy06 = 6, &
ko_dummy07 = 7, &
ko_dummy08 = 8, &
ko_dummy09 = 9, &
ko_dummy10 = 10, &
ko_dummy11 = 11, &
ko_dummy12 = 12, &
ko_dummy13 = 13, &
ko_dummy14 = 14, &
ko_dummy15 = 15
#endif
LOGICAL, PARAMETER :: gas_mix_fatal = .false.
END MODULE module_kdm_parameters
!-----------------------------------------------------------------------
! This is essentially the Hadley/Unified Model radiation code modified !
! to use a correlated-k scheme. The subroutine planetary_kdm is in !
! many ways a "harness" which sets up the variables needed in the !
! radiation code. Subroutine flux_calc is the nominal start of the !
! Hadley scheme. Note that it is called twice per implementation, !
! once for the solar, and once for the IR. !
! There will be a loop over 'j', or latitude, doing latitude strips !
! one at a time. !
! In this code, n_layer is at the surface, and 0 is above the top !
! layer. !
!-----------------------------------------------------------------------
MODULE module_ra_kdm
USE module_wrf_error
USE module_kdm_parameters
USE module_planet_utilities, ONLY: rsat_h2o
#if ( WRF_MARS == 1 )
USE module_state_description, ONLY : TES_LIMB_AVG
USE module_ra_valverde
USE module_state_description, ONLY: VALVERDE_LOOKUP, &
VALVERDE_FIT
#endif
IMPLICIT NONE
! Make sure the declared header common variables don't interfere
! with similarly-named variables in routines that USE this module
! by making all variables PRIVATE and then making only the
! necessary subroutines PUBLICally available
PRIVATE
PUBLIC :: planetary_kdm, kdminit, &
#if ( WRF_MARS == 1 )
aerosol_obs_scaling, build_ko_options, &
aerosol_obs_scaling_1particle, &
reference_opacity, reference_extinction, &
#endif
solar_fraction
REAL(KIND(0.d0)), ALLOCATABLE, DIMENSION(:,:,:,:,:) :: &
! DIMENSION(n_temps,n_press,n_terms,n_tot_bands,n_mix) :: & ! Full array of k-coefficients
kval_array
REAL(KIND(0.d0)), DIMENSION(n_band_solar) :: solar_flux_band, &
rayleigh_coefficient
#if ( WRF_MARS == 1 )
REAL(KIND(0.d0)), DIMENSION(n_band_solar) :: kCext_dust_solar, &
omega_dust_solar, &
asymmetry_dust_solar, &
kCext_ice_solar, &
omega_ice_solar, &
asymmetry_ice_solar, &
kCext_co2ice_solar, &
omega_co2ice_solar, &
asymmetry_co2ice_solar
REAL(KIND(0.d0)), DIMENSION(n_band_ir) :: &
kCext_dust_IR, &
omega_dust_IR, &
asymmetry_dust_IR, &
kCext_ice_IR, &
omega_ice_IR, &
asymmetry_ice_IR, &
kCext_co2ice_IR, &
omega_co2ice_IR, &
asymmetry_co2ice_IR
#endif
#if ( WRF_MARS == 1 )
REAL(KIND(0.d0)) :: kCext_dust_ref, &
kCext_dust_ref_ir, & ! reference at dust IR diagnostic wavelength
kCext_ice_ref, &
kCext_co2ice_ref
#endif
INTEGER :: n_fake_layers_needed_visir
REAL(KIND(0.d0)), DIMENSION(max_fake_layers_visir) :: &
p_fake_lower_edge, &
p_fake_centre, &
t_fake_lower_edge, &
t_fake_centre
LOGICAL :: l_double
LOGICAL :: kval_debug = .false.
#if ( WRF_MARS == 1 )
! Chicago water ice cloud layers
REAL, SAVE :: chi_cloud_prs, chi_ls1, chi_ls2, chi_lat1, chi_lat2, chi_lon1, chi_lon2, chi_ice_tau
REAL, SAVE :: chi_tval1, chi_tval2
INTEGER, SAVE :: chi_ncloud, chi_dayornight
!multiple particle size mie tables
real(kind(0.d0)), save, dimension(:,:,:,:), allocatable :: aerosol_mie
real(kind(0.d0)), save, dimension(:,:,:), allocatable :: aerosol_Qref
real(kind(0.d0)), save, dimension(:), allocatable :: Qref
real(kind(0.d0)) :: reffstart, reffstop, reffstep
integer :: nbins_mie
integer :: kdm_options
logical :: is_two_moment_dust, is_two_moment_water_ice
integer :: n_dust_data, n_ice_data, n_co2ice_data ! Number of rows in the various ice and dust parameter data files
#endif
#if ( WRF_TITAN == 1 )
integer, parameter :: n_tholin_rad_bins = 20
real, dimension(n_tholin_rad_bins) :: tholin_rad_bins
real(kind(0.d0)), save, dimension(n_band_solar,n_tholin_rad_bins) :: &
tholin_vis_ext_xsect, &
tholin_vis_ssa, &
tholin_vis_asym
real(kind(0.d0)), save, dimension(n_band_ir, n_tholin_rad_bins) :: &
tholin_ir_ext_xsect, &
tholin_ir_ssa, &
tholin_ir_asym
#endif
CONTAINS
!=======================================================================
SUBROUTINE planetary_kdm(rthratensw, rthratenlw, &
gsw, glw, &
toa_sw_u, toa_sw_d, toa_lw_u, toa_lw_d, &
flux_sw_u,flux_sw_d,flux_lw_u,flux_lw_d, &
xlat, xlong, albedo, &
emiss, ha, coszen, angslope, azmslope, &
dhr, diurnavg, rho_phy, t3d, tf3d, tsk, qv3d, &
qc3d, qr3d, qi3d, qs3d, qg3d, p3d, pf3d, pi3d, &
dz8w, julday, &
#ifdef mpas
p2si, &
#endif
hr_g_vis, hr_a_vis, hr_vis, &
hr_g_ir, hr_a_ir, hr_ir, &
sunfrac, julian, declin, solcon, &
l_s, hrang, &
f_qv, f_qc, f_qr, f_qi, f_qs, f_qg, & ! flags for whether moisture variables are present in WRF
i_2stream_solar, i_2stream_ir, &
nlte_physics, &
sw_correct, ra_du_physics, &
do_fake_layers_in, &
kdm_mode_flag, &
t_sounding, polar, p_scale, &
verbose_debug_in,nan_detail_check_in, &
#if ( WRF_MARS == 1 )
dust_array, cloud_array,co2_cloud_array, &
dust_array_ir, &
do_chicago_ice_layers, &
f_qst_co2ice, &
include_dust, include_water_ice, &
include_water_vap, include_co2_ice, &
rhod_h2o,rhod_co2,h2oice_radius,co2ice_radius, &
qst_co2ice, &
reff_dust, ndust, &
reff_ice, nice, &
cdod_scale, &
#endif
#if ( WRF_TITAN == 1 )
mck_haze_numden, mck_haze_radius, &
mck_haze_tau64, &
include_tholin, mp_physics_tholin, &
q_ch4_3d, f_q_ch4, include_ch4_vap, &
#endif
ra_kdm_gray_k_adhoc_ir, &
ra_kdm_gray_k_adhoc_vis, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte, &
ra_kdm_taulayer_calculation, &
ra_kdm_clearsky_calculation, &
swdown, swdowndir, &
swdownunscat &
)
USE module_ra_utils
#ifdef mpas
USE mpas_atmphys_constants,only: pa2bar, p0, rcp
#else
USE module_model_constants
#endif
#if ( WRF_MARS == 1 )
USE module_ra_mars_common, only: cp_mars
#endif
USE module_planet_utilities
IMPLICIT NONE
INTEGER, INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte, &
i_2stream_solar, & ! Two-stream method for solar wavelengths
i_2stream_ir, & ! Two-stream method for IR wavelengths
julday, &
nlte_physics, &
ra_du_physics
LOGICAL, INTENT(IN ) :: sw_correct, &
#if ( WRF_MARS == 1 )
include_dust, &
include_water_ice, &
include_water_vap, &
include_co2_ice, &
#endif
#if ( WRF_TITAN == 1 )
include_tholin, &
include_ch4_vap, &
#endif
do_fake_layers_in
LOGICAL, OPTIONAL :: f_qv,f_qc,f_qr,f_qi,f_qs,f_qg
#if ( WRF_MARS == 1 )
LOGICAL, OPTIONAL :: f_qst_co2ice
#endif
#if ( WRF_TITAN == 1 )
LOGICAL, OPTIONAL :: f_q_ch4
#endif
LOGICAL, INTENT(IN ) :: diurnavg
logical, intent(in) :: ra_kdm_taulayer_calculation, &
ra_kdm_clearsky_calculation
INTEGER, INTENT(IN ) :: kdm_mode_flag ! 0=vis, 1=IR see kdm_do_vis and kdm_do_ir definitions
integer, intent(in ) :: t_sounding ! method for temperature sounding, default = 2 (constant lapse rate), also used by fake layers
logical, intent(in ) :: polar ! gcm (.true.) or les (.false.) modeling, also used by fake layers
real, intent(in ) :: p_scale ! pressure scalar, default = 1., also used by fake layers
#if ( WRF_MARS == 1 )
logical, intent(in ) :: do_chicago_ice_layers ! use water ice opacity information from presc_aero block of the namelist
#endif
#if ( WRF_TITAN == 1 )
integer, intent(in ) :: mp_physics_tholin
#endif
REAL(KIND(0D0)), INTENT(IN ) :: dhr
REAL, DIMENSION(n_band_solar), INTENT(IN ) :: &
ra_kdm_gray_k_adhoc_vis
REAL, DIMENSION(n_band_ir), INTENT(IN ) :: &
ra_kdm_gray_k_adhoc_ir
REAL, INTENT(IN ) :: julian, &
declin, &
solcon, &
l_s
#ifdef mpas
REAL, INTENT(IN ) :: p2si
#endif
#if ( WRF_MARS == 1 )
REAL, INTENT(IN) :: rhod_h2o,rhod_co2,h2oice_radius,co2ice_radius
#endif
LOGICAL, INTENT(IN) :: verbose_debug_in, nan_detail_check_in
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: &
xlat, &
xlong, &
albedo, & ! Surface albedo
emiss, & ! Surface emissivity--generally pegged at 1.0 for IR
angslope, &
azmslope, &
sunfrac, &
tsk, &
ha, &
coszen, &
hrang
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
p3d, &
pf3d, &
pi3d, &
t3d, &
tf3d, &
rho_phy, &
dz8w
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, &
INTENT(IN ) :: &
qv3d, &
qc3d, &
qr3d, &
qi3d, &
qs3d, &
qg3d
#if ( WRF_MARS == 1 )
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, &
INTENT(IN ) :: &
qst_co2ice
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: &
dust_array, & ! dust_array is definitionally VIS (0.67micron) for MarsWRF
cloud_array, &
co2_cloud_array
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( OUT) :: &
dust_array_ir ! dust_array_ir is non-trivial for 2-moment (diagnostic)
#endif
#if ( WRF_TITAN == 1 )
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
mck_haze_numden, &
mck_haze_radius
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: &
mck_haze_tau64
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, &
INTENT(IN ) :: &
q_ch4_3d
#endif
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: &
rthratensw, &
rthratenlw, &
flux_sw_u, &
flux_sw_d, &
flux_lw_u, &
flux_lw_d
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
gsw, &
glw, &
toa_sw_d, & ! Top of Atmosphere Short Wave (TOASW) flux (down)for heat balance analysis.
toa_sw_u, & ! Top of Atmosphere Short Wave (TOASW) flux (up ) for heat balance analysis.
toa_lw_d, & ! Top of Atmosphere Long Wave (TOALW) flux (down) for heat balance analysis.
toa_lw_u ! Top of Atmosphere Long Wave (TOALW) flux (up) for heat balance analysis.
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( OUT) :: &
hr_g_vis, & ! Heating rate due to gas in solar wavelengths
hr_a_vis, & ! Heating rate due to aerosol in solar wavelengths
hr_vis, & ! Total heating rate in solar wavelengths
hr_g_ir, & ! Heating rate due to gas in IR wavelengths
hr_a_ir, & ! Heating rate due to aerosol in IR wavelengths
hr_ir ! Total heating rate in IR wavelengths
REAL, DIMENSION ( ims:ime, jms:jme), INTENT(INOUT), OPTIONAL :: &
swdown, & ! sw down at the surface (i.e. not net)
swdowndir, & ! direct beam sw down at the surface
swdownunscat ! unscattered beam sw down at the surface
#if ( WRF_MARS == 1 )
REAL, DIMENSION(ims:ime,kms:kme,jms:jme), intent(in), optional &
:: ndust, reff_dust
REAL, DIMENSION(ims:ime,kms:kme,jms:jme), intent(in), optional &
:: nice, reff_ice
REAL, DIMENSION(ims:ime,jms:jme), optional, intent(in) &
:: cdod_scale
#endif
!c-mm Local variables
! These variables are going to take the values:
! n_layer_wrf = kte-kts+1
! n_layer_visir_m = n_layer_wrf + max_fake_layers_visir
INTEGER :: n_layer_wrf , &
n_layer_visir_m ! for memory, not for looping
INTEGER :: n_layer_visir_l ! looping variable, needs creating from needed fake layers from init
! ADT: Enforced consistent dimensioning of arrays to ensure
! correct behavior on single and multiple processors
REAL(KIND(0.d0)), DIMENSION( ite-its+1, kte-kts+1+max_fake_layers_visir ) :: &
qv2d, &
qc2d, &
qr2d, &
qi2d, &
qs2d, &
#if ( WRF_MARS == 1 )
ndust2d, &
reff_dust2d, &
nice2d, &
reff_ice2d, &
#endif
#if ( WRF_TITAN == 1 )
n_tholin, &
r_tholin, &
#endif
qg2d
REAL(KIND(0.d0)), &
DIMENSION( ite-its+1, kte-kts+1+max_fake_layers_visir, npd_species ) :: & ! Array of gas mixing ratios for all gases in atmosphere
gas_mix_ratio
REAL, DIMENSION( kte-kts+1 ) :: hr_co2_nlte
! REAL, DIMENSION( kte-kts+2 ) :: t1d, &
! pf1d
REAL(KIND(0.d0)), &
DIMENSION( ite-its+1, 0:kte-kts+1+max_fake_layers_visir, npd_aerosol_species ):: & ! Array of aerosol opacities for all aerosol species in atmosphere
aerosol_opacity
#if ( WRF_TITAN == 1 )
REAL, DIMENSION( ite-its+1, kte-kts+1+max_fake_layers_visir ):: tau_haze_64
#endif
REAL, DIMENSION( ite-its+1, kte-kts+1+max_fake_layers_visir, jte-jts+1 ) :: &
hr_tot_vis, & ! Total (gas+dust+ice) heating rate in atmosphere for solar wavelengths
hr_tot_ir ! Total (gas+dust+ice) heating rate in atmosphere for IR wavelengths
REAL, DIMENSION( ite-its+1 ) :: xlat1d, & !-------!
xlong1d, & !
gsw1d, & !
ha1d, & ! One-dimensional subarrays of above 2-D arrays
alphaslp1d, & !
gammaslp1d, & !
sunfrac1d !-------!
REAL(KIND(0.d0)), DIMENSION(ite-its+1) :: &
solar_toa, & ! Incoming solar flux at top of the atmosphere
sec_0, & ! Solar zenith angle
tsk1d, & ! One-dimensional surface temperature
alb1d, & ! One-dimensional surface albedo
em1d, & ! One-dimensional surface emissivity
cosfac !
REAL(KIND(0.d0)), DIMENSION(kte-kts+1) :: xltecorrection
REAL(KIND(0.d0)), DIMENSION(ite-its+1,kte-kts+1) :: pnorm ! Two-dimensional layer pressure in normal model coordinates
REAL(KIND(0.d0)), DIMENSION(ite-its+1,kte-kts+1+max_fake_layers_visir) :: p2d, & ! Two-dimensional layer pressure
t2d ! Two-dimensional layer temperature
#if ( WRF_MARS == 1 )
REAL(KIND(0.d0)), DIMENSION(kte-kts+1+max_fake_layers_visir) :: z_lv, &
sigma_lv
#endif
REAL(KIND(0.d0)), DIMENSION(ite-its+1,0:kte-kts+1+max_fake_layers_visir) ::pf2d, & ! Two-dimensional level pressure
tf2d ! Two-dimensional level temperature
REAL(KIND(0.d0)), DIMENSION(ite-its+1,0:kte-kts+1+max_fake_layers_visir,2) :: &
flux_direct, & ! Direct fluxes for aerosol/no aerosol
flux_down, & ! Downward fluxes for aerosol/no aerosol
flux_up, & ! Upward fluxes for aerosol/no aerosol
flux_unscattered ! Unscattered downward fluxes for aerosol/no aerosol
REAL(KIND(0.d0)), DIMENSION(ite-its+1,kte-kts+1+max_fake_layers_visir) :: d_mass ! Mass of an atmospheric layer
REAL(KIND(0.d0)), DIMENSION(ite-its+1,kte-kts+1+max_fake_layers_visir) :: d_z ! Vertical distance thickness of an atmospheric layer
REAL(KIND(0.d0)), DIMENSION(ite-its+1,kte-kts+1+max_fake_layers_visir,npd_band) :: &
trans_layer
REAL :: juldayfrac, &
tsunfrac, &
cosfac0d, &
cpm ! heat capacity of actual air parcel inc. moisture
REAL(KIND(0.d0)) :: ref_extinction, top_n_potl, t_frost, top_n_kin, avg_weight
INTEGER :: i, j, k, ii, jj, kk, l, &
n_profile, &
n_layer, &
pos, &
n_to_average, &
n, j_end
LOGICAL :: l_rayleigh, bool, &
include_aerosols, &
verbose_debug, &
nan_detail_check, &
do_warn
REAL(KIND(0.d0)) :: p_top_check
INTEGER, DIMENSION(2) :: locations
INTEGER :: i_at_min, i_at_max, k_at_min, k_at_max
LOGICAL, PARAMETER :: from_profile = .false., & ! take T(z) from profile or (if false, then) from fit to top 5 prognostic layers
use_kinetic = .true. ! use kinetic temps from top 5 layers or (if false, then) use potential temps
LOGICAL, PARAMETER :: l_rayleigh_ir = .false. ! to treat rayleigh scattering set equal to true
LOGICAL, PARAMETER :: l_rayleigh_vis = .true. ! to treat rayleigh scattering set equal to true
LOGICAL :: do_fake_layers
! write(wrf_err_message,*) "solcon=",solcon
! call wrf_message(trim(wrf_err_message))
do_warn = .false.
nan_detail_check = nan_detail_check_in
verbose_debug = verbose_debug_in
if((wrf_debug_level==50) .or. (wrf_debug_level>= 999)) verbose_debug=.true.
if((wrf_debug_level==60) .or. (wrf_debug_level>= 999)) kval_debug=.true.
if(verbose_debug) nan_detail_check = .true.
include_aerosols = .false. ! false unless we actually request aerosols via include_dust, etc.
#if ( WRF_TITAN == 1 )
include_aerosols = include_tholin
#elif ( WRF_MARS == 1 )
include_aerosols = (include_dust).OR.(include_water_ice).OR.(include_co2_ice)
#endif
do_fake_layers = do_fake_layers_in
n_layer_wrf = kte-kts+1
!c-mm write(0,*) 'zz',ra_kdm_gray_k_adhoc_ir
#if ( WRF_MARS == 1 )
if((.not.include_dust).or.is_two_moment_dust) then
dust_array(its:ite,kts:kte,jts:jte) = 0.
dust_array_ir(its:ite,kts:kte,jts:jte) = 0.
endif
if((.not.include_water_ice).or.is_two_moment_water_ice) then
cloud_array(its:ite,kts:kte,jts:jte) = 0.
endif
if(.not.include_co2_ice) then
co2_cloud_array(its:ite,kts:kte,jts:jte) = 0.
endif
#endif
IF(do_fake_layers) THEN
!c-yl: prepare fake layers for KDM: set number of fake layers, and corresponding pressure and temperature
! (in wrf, pf3d at k=kde (not kme!) is the model top - kme might be the top, or there might be pad)
call set_fake_layers(p_top=minval(pf3d(its:ite,kde,jts:jte)), t_sounding=t_sounding, polar=polar, p_scale=p_scale)
ELSE
n_fake_layers_needed_visir = 0
ENDIF
p_top_check=MINVAL(pf3d(its:ite,kde,jts:jte))
n_layer_visir_l = n_layer_wrf + n_fake_layers_needed_visir
n_layer_visir_m = n_layer_wrf + max_fake_layers_visir ! for memory, not for looping
IF(do_fake_layers) THEN
if(n_layer_visir_l .gt. n_layer_visir_m) then
WRITE ( wrf_err_message ,*) 'ERROR: KDM has not set up enough total (real+fake) layers in memory, ',&
n_layer_visir_l,' > ',n_layer_visir_m
CALL wrf_error_fatal ( wrf_err_message )
ENDIF
if(n_fake_layers_needed_visir .le. 0) then
WRITE ( wrf_err_message ,*) 'ERROR: KDM needs at least one ("fake") layer (for vis and ir) above the model ptop, ',&
n_fake_layers_needed_visir
CALL wrf_error_fatal ( wrf_err_message )
ENDIF
ENDIF ! -- end of fake layers bit --
qv2d(:,:)= 0.d0 ! default all local water types to zero before reading any in over them
qc2d(:,:)= 0.d0 ! implicit looping is threadsafe as these are local and only defined on tile
qr2d(:,:)= 0.d0
qi2d(:,:)= 0.d0
qs2d(:,:)= 0.d0
qg2d(:,:)= 0.d0
gas_mix_ratio(:,:,:) = 0.d0 ! best not to just pick up junk lying around in memory
!c-mm call wrf_debug(500,'---> Entered SUBROUTINE planetary_kdm')
! Typically, in a WRF run, on the processors occupying the edge
! of a particular axis:
! kte = kde = # of Z points including extra stagger
! ite = ide = # of X points including extra stagger
! Don't calculate radiation at the extra x-stagger point:
#ifdef mpas
n_profile=ite-its+1
#else
n_profile=MIN(ite,ide-1)-its+1
#endif
! On the off chance the model is being decomposed in the
! Z-direction, kte may not equal kde. If kte does equal kde,
! then we want one less than the number including the extra
! z-stagger.
IF((kde-1) /= kte) THEN
WRITE ( wrf_err_message ,*) 'ERROR: KDM has not been setup to vertically tile, kte= ',kte,' kde= ',kde
CALL wrf_error_fatal ( wrf_err_message )
ENDIF
#if ( WRF_MARS == 1 )
! local arrays defined to have full extent same as tile, so no need to worry about writing outside of tile:
ndust2d(:,:) = 0.d0
reff_dust2d(:,:) = 0.d0
nice2d(:,:) =0.d0
reff_ice2d(:,:) = 0.d0
#endif
#if ( WRF_TITAN == 1 )
n_tholin(:,:) = 0.d0
r_tholin(:,:) = 1.d-3 ! minimum so scattering doesn't go crazy (in microns)
#endif
!c-mm I've determined the best place to calculate the cloud_array values is here, before we
!c-mm start looping. By doing it here, we have all the variables we need within the
!c-mm module, and don't have to go passing it around. qi3d and/or qst_co2ice are passed into
!c-mm this subroutine, which makes it the easiest way to do things.
!c-mm The data in kCext_ice_ref comes from a source which actually provides Qext (instead of k*Cext), so the
!c-mm name is a misnomer (although it's correct for the CO2 ice below). Just need to pass this value
!c-mm into the subroutine directly.
#if ( WRF_MARS == 1 )
if((ra_du_physics /= TES_LIMB_AVG) .and. present(f_qi) .and. present(qi3d) .and.(.not.do_chicago_ice_layers)) then
if(include_water_ice .and. f_qi) then
!c-mm call wrf_debug(500,'---> doing calc_ice_tau for h2o')
CALL CALC_ICE_TAU(qi3d, cloud_array, pf3d, &
kCext_ice_ref, &
H2OICE_RADIUS, RHOD_H2O, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte)
endif
endif
if(do_chicago_ice_layers .and. include_water_ice) then
!c-mm call wrf_debug(500,'---> doing ice_tau_chicago')
CALL ice_tau_chicago(cloud_array, pf3d, &
l_s, xlat, xlong, hrang, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte)
endif
!c-mm The 1.e6 multiplier of part_radius is because lmie code (offline) calculates in microns, and CO2ICE_RADIUS is in m.
!c-mm The value of kCext_ref is essentially in microns. Since we ultimately want Qext, which is kCext_ref/pi*r^2
!c-mm we need to convert the units of part_radius to coincide with those of kCext_ref.
!c-mm Unlike for water ice above, the data in the CO2 ice data file is truly k*Cext, so we need to calculate Qext from
!c-mm it and the geometric cross section of the ice particles (pi*r^2)
if(present(f_qst_co2ice) .and. present(qst_co2ice)) then
if(include_co2_ice .and. f_qst_co2ice) then
!c-mm call wrf_debug(500,'---> doing calc_ice_tau for co2')
CALL CALC_ICE_TAU(qst_co2ice, co2_cloud_array, pf3d, &
kCext_co2ice_ref, &
CO2ICE_RADIUS, RHOD_CO2, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte)
endif
endif
#endif
#ifdef mpas
j_end=jte
#else
j_end=MIN(jte,jde-1)
#endif
!c-mm call wrf_debug(500,"--- about to start j loop")
j_loop: DO j=jts,j_end
! Convert to lower bound of 1 (local) from lower bound of
! jts (WRF distributed memory, or serial with jts /= 1)
jj = j-jts+1
IF (kval_debug) then
write(wrf_err_message,*) '--- ---> j = ',j,' out of ',MIN(jte,jde-1)-jts+1
call wrf_message(trim(wrf_err_message))
ENDIF
XLTECORRECTION = 1.0
!c-mm Within the KDM code, the index 'i' is the local loop counter
!c-mm (only for those variables that are local to KDM).
!c-mm The variable 'ii' is the global loop counter, used for those
!c-mm variables that are passed into or out of KDM.
!c-mm **NOTE**, however, that the nomenclature is opposite for the j
!c-mm values. The index 'j' is the global loop counter and 'jj' is
!c-mm the local loop counter. I have no idea why the hell I did
!c-mm it this way...
DO i=1,n_profile
! Convert from lower bound of 1 (local) to lower bound of
! its (WRF distributed memory, or serial with its /= 1)
ii = i+its-1
! Variables on the z-stagger
! Now setup the vis and ir (main kdm) layer structure including 'fake' layers to treat the atmosphere above the dynamical p_top
! The KDM model layer structure and how it relates to the wrf layer structure:
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
! For the main vis and ir kdm (i.e. not the LTE correction)
!
! ===== is a layer edge
! ----- is a layer centre
!
! KDM layer <==> WRF layer
!
! space aka p=0 k=0 =========
! layer centre k=1 - - - - -
! layer edge k=1 =========
! ....
! layer centre k=n_fake_layer_needed_visir-1 - - - - -
! layer edge k=n_fake_layer_needed_visir-1 =========
! layer centre k=n_fake_layer_needed_visir - - - - -
! layer edge k=n_fake_layer_needed_visir ========= 0 ========= k=kte model top AKA p_top
! layer centre k=n_fake_layer_needed_visir+1 - - - - - 1 - - - - - k=kte-1
! layer edge k=n_fake_layer_needed_visir+1 ========= 1 ========= k=kte-1
! layer centre k=n_fake_layer_needed_visir+2 - - - - - 2 - - - - - k=kte-2
! ....
! layer centre k=n_fake_layer_needed_visir-1 - - - - - - - - - - k=kts+1
! layer edge k=n_fake_layer_needed_visir-1 ========= ========= k=kts+1
! layer centre k=n_fake_layer_needed_visir - - - - - - - - - - k=kts first layer centre
! surface k=n_layer_visir_l ========= n_layer_wrf ========= k=kts model bottom AKA psurf
!
! (above ascii arrow, these are KDM "real" layers /\ )
! figure out the simple (naive) average potential temperature in top n layers - where n = min( ! 1+nint(real(n_fake_layers_needed_visir)/4) , 5 )
! -- this formula has no deep meaning, it's arbitrary and just defined to grab a "reasonable" number of layers. It's max'ed to 5
! so as not to grab a crazy number of layers below the dynamical model top.
IF(do_fake_layers) THEN
n_to_average = min(1+int(real(n_fake_layers_needed_visir)/4),5)
top_n_potl = 0.d0
top_n_kin = 0.d0
avg_weight = 1.d0/real(n_to_average,kind(0.d0))
do kk=kde-1,kde-n_to_average,-1 ! this should be a loop from top down to n below top in wrf-space (layer centres)
top_n_potl = top_n_potl + REAL( (avg_weight * tf3d(ii,kk,j) * ( p0 / pf3d(ii,kk,j))**rcp) ,KIND(0.d0))
top_n_kin = top_n_kin + avg_weight * tf3d(ii,kk,j)
enddo
if(verbose_debug .and. (ii==(its+ite)/2) .and. (j==jts)) then
write(0,*) i,j,' number of points in top average=',n_to_average
write(0,*) ' p_top=',minval(pf3d(its:ite,kde,jts:jte)),' pressure_TOA=',pressure_TOA
write(0,*) ' top n potl, kin',top_n_potl,top_n_kin
write(0,*) ' at p_top, this is',real(top_n_potl)*((pf3d(ii,kde,j)/p0)**rcp)
endif
! setup fake vis+ir layer bottom edges
do k=1,n_fake_layers_needed_visir
pf2d(i,k) = p_fake_lower_edge(1+n_fake_layers_needed_visir - k)
if(from_profile) then
tf2d(i,k) = t_fake_lower_edge(1+n_fake_layers_needed_visir - k)
else
if(use_kinetic) then
tf2d(i,k) = top_n_kin
else
tf2d(i,k) = top_n_potl*((pf2d(i,k)/real(p0,kind(0.d0)))**real(rcp,kind(0.d0)))
endif
endif
aerosol_opacity(i,k,:)=0.
#if ( WRF_MARS == 1 )
t_frost = real(tfrost_co2(real(pf2d(i,k))),kind(0.d0)) ! function from module_planet_utilities
if(tf2d(i,k) .lt. t_frost) tf2d(i,k) = t_frost
#endif
enddo
ENDIF ! -- fake layer bit --
! Real layer edge info for the vis and ir
DO k=n_fake_layers_needed_visir+1,n_layer_visir_l ! this is a loop from the top->down for KDM
kk=kde-(k-n_fake_layers_needed_visir)
tf2d(i,k) = REAL(tf3d(ii,kk,j),KIND(0.d0))
pf2d(i,k) = REAL(pf3d(ii,kk,j),KIND(0.d0)) ! pf is pressure at the full levels (p8w in the regular code)
aerosol_opacity(i,k,:)=0. ! Initialize aerosol_opacity to zero everywhere (i.e. 'clear-sky' case)
#if ( WRF_MARS == 1 )
!c-mm So, what's the intent here? If we are including dust, then we want the internal KDM routine
!c-mm to receive the opacity values (from dust_array, which is passed in from radiation_driver)
!c-mm and operate accordingly. If we are *not* including dust, then we want to overwrite the
!c-mm dust_array values with zeroes at this stage. Just above, we initialize aerosol_opacity to
!c-mm zero, so the KDM routine will run with aerosol_opacity set to zero anyway. But, by setting
!c-mm dust_array to zero here, it passes it back into radiation_driver. Then, when we get to the
!c-mm diagnose_tau routines, the values of dust_array are zero, and it spits out zeros in the wrfout
!c-mm files. Otherwise, the KDM code will assume tau=0, but the wrfout output will still reflect
!c-mm the values of tau that were calculated in the du_phys routines.
IF (include_dust) THEN
aerosol_opacity(i,k,1)= &
REAL(dust_array(ii,kk,j),KIND(0.d0)) ! Take only a 2-D slice of dust_array
! can diagnose the IR (9.3micron) optical depth from the VIS (0.67micron) as we know that the optical depth
! is DEFINED as VIS for non-two-moment cases:
dust_array_ir(ii,kk,j) = dust_array(ii,kk,j) * kCext_dust_ref_ir/kCext_dust_ref
ELSE
dust_array(ii,kk,j)=0.0
dust_array_ir(ii,kk,j)=0.0
ENDIF
!c-mm Same as above, but for cloud_array. The values in cloud_array are either determined from the tes_limb
!c-mm routine (if du_phys=49) or else they're calculated above in this routine (from qi3d). Either way
!c-mm if water ice is not being used, we want to zero this out in cloud_array so radiation_driver properly
!c-mm sees zeroes.
aerosol_opacity(i,k,2)=0.d0
IF (include_water_ice) THEN ! can't check f_qi here since TES_LIMB can also set water ice opac
aerosol_opacity(i,k,2)= &
REAL(cloud_array(ii,kk,j),KIND(0.d0)) ! Take only a 2-D slice of cloud_array
ELSE
cloud_array(ii,kk,j)=0.0
ENDIF
!c-mm Same as above, but for co2_cloud_array. The values in co2_cloud_array are always calculated above, in this
!c-mm routine (from qst_co2ice). Make sure they get returned as zeroes if we are not considering co2 ice.
aerosol_opacity(i,k,3)=0.d0
IF (include_co2_ice .and. present(qst_co2ice) .and. present(f_qst_co2ice)) THEN
if(f_qst_co2ice) then
aerosol_opacity(i,k,3)= &
REAL(co2_cloud_array(ii,kk,j),KIND(0.d0)) ! Take only a 2-D slice of cloud_array
else
co2_cloud_array(ii,kk,j)=0.0
endif
ELSE
co2_cloud_array(ii,kk,j)=0.0
ENDIF
#endif
ENDDO
!c-mm call wrf_debug(500,"-- just about to do fake layers")
IF(do_fake_layers) THEN
#if ( WRF_MARS == 1 )
! don't need checks of f_qv, f_qst_co2ice, etc., below, since if everything is not 100% kosher, the respective _array = 0.
IF (include_dust) THEN ! fill fake layers with top layer aerosol so that we don't get a big
aerosol_opacity(i,0:n_fake_layers_needed_visir,1)= & ! "transparency" jump at top of prognostic model
REAL(dust_array(ii,kde-1,j),KIND(0.d0))
ENDIF
IF (include_water_ice) THEN
aerosol_opacity(i,0:n_fake_layers_needed_visir,2)= &
REAL(cloud_array(ii,kde-1,j),KIND(0.d0))
ENDIF
IF (include_co2_ice) THEN
aerosol_opacity(i,0:n_fake_layers_needed_visir,3)= &
REAL(co2_cloud_array(ii,kde-1,j),KIND(0.d0))
ENDIF
#endif
ENDIF ! -- end fake layers bit --
! Set Top Of Atmosphere values i.e. at top of a fake layer for the vis and ir
k=0 ! k=0 is the top of the atmosphere - no layer centre information is stored here, but the edge info (pf and tf) is stored here
tf2d(i,k) = tf2d(i,1)
pf2d(i,k) = pressure_TOA ! set TOA at vey small number (approx 0 Pa)
aerosol_opacity(i,k,:)=0. ! sat TOA opacity to zero for all aerosols
!c-mm aerosol_opacity(i,:,1)=aerosol_opacity(i,:,1)*0.05 !c-mm Hardwire an 95% reduction in opacity for testing, 12-30-17
! Varaibles in the middle of a layer
DO k = 1,n_layer_wrf
pnorm(i,k) = REAL(p3d(ii,k,j),KIND(0.d0)) ! pnorm, used in the NLTE code, has the same orientation as WRF, opposite of KDM
ENDDO
DO k=1,n_layer_visir_l
kk= kts + (n_layer_visir_l - k)
if (k.le.n_fake_layers_needed_visir) then
p2d(i,k) = p_fake_centre(n_fake_layers_needed_visir - (k-1))
if(from_profile) then
t2d(i,k) = t_fake_centre(n_fake_layers_needed_visir - (k-1))
else
if(use_kinetic) then
t2d(i,k) = top_n_kin
else
t2d(i,k) = top_n_potl*((p2d(i,k)/real(p0,kind(0.d0)))**real(rcp,kind(0.d0)))
endif
endif
#if ( WRF_MARS == 1 )
t_frost = real(tfrost_co2(real(p2d(i,k))),kind(0.d0)) ! function from module_planet_utilities
if(tf2d(i,k) .lt. t_frost) t2d(i,k) = t_frost
#endif
qv2d(i,k) = 0.d0 ! default all local water types to zero before reading any in over them
qc2d(i,k) = 0.d0
qr2d(i,k) = 0.d0
qi2d(i,k) = 0.d0
qs2d(i,k) = 0.d0
qg2d(i,k) = 0.d0
else
p2d(i,k) = REAL(p3d(ii,kk,j),KIND(0.d0))
t2d(i,k) = REAL(t3d(ii,kk,j),KIND(0.d0))
if(qv3d(ii,kk,j) .gt. 0.1) write(0,*) "qv3d bad at ",ii,kk,k,qv3d(ii,kk,j)