-
Notifications
You must be signed in to change notification settings - Fork 312
/
SurfaceAlbedoMod.F90
1701 lines (1492 loc) · 90.1 KB
/
SurfaceAlbedoMod.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
module SurfaceAlbedoMod
#include "shr_assert.h"
!-----------------------------------------------------------------------
! !DESCRIPTION:
! Performs surface albedo calculations
!
! !PUBLIC TYPES:
use shr_kind_mod , only : r8 => shr_kind_r8
use shr_log_mod , only : errMsg => shr_log_errMsg
use decompMod , only : bounds_type, subgrid_level_patch
use abortutils , only : endrun
use landunit_varcon , only : istsoil, istcrop, istdlak
use clm_varcon , only : grlnd
use clm_varpar , only : numrad, nlevcan, nlevsno, nlevcan
use clm_varctl , only : fsurdat, iulog, use_snicar_frc, use_SSRE
use pftconMod , only : pftcon
use SnowSnicarMod , only : sno_nbr_aer, SNICAR_RT, DO_SNO_AER, DO_SNO_OC
use AerosolMod , only : aerosol_type
use CanopyStateType , only : canopystate_type
use LakeStateType , only : lakestate_type
use SurfaceAlbedoType , only : surfalb_type
use TemperatureType , only : temperature_type
use WaterStateBulkType , only : waterstatebulk_type
use WaterDiagnosticBulkType , only : waterdiagnosticbulk_type
use GridcellType , only : grc
use LandunitType , only : lun
use ColumnType , only : col
use PatchType , only : patch
!
implicit none
!
! !PUBLIC MEMBER FUNCTIONS:
public :: SurfaceAlbedo_readnl
public :: SurfaceAlbedoInitTimeConst
public :: SurfaceAlbedo ! Surface albedo and two-stream fluxes
!
! !PRIVATE MEMBER FUNCTIONS:
private :: SoilAlbedo ! Determine ground surface albedo
private :: TwoStream ! Two-stream fluxes for canopy radiative transfer
!
! !PUBLIC DATA MEMBERS:
! The CLM default albice values are too high.
! Full-spectral albedo for land ice is ~0.5 (Paterson, Physics of Glaciers, 1994, p. 59)
! This is the value used in CAM3 by Pritchard et al., GRL, 35, 2008.
! albedo land ice by waveband (1=vis, 2=nir)
real(r8), public :: albice(numrad) = (/ 0.80_r8, 0.55_r8 /)
! namelist default setting for inputting alblakwi
real(r8), public :: lake_melt_icealb(numrad) = (/ 0.10_r8, 0.10_r8/)
! albedo frozen lakes by waveband (1=vis, 2=nir)
! unclear what the reference is for this
real(r8), private :: alblak(numrad) = (/0.60_r8, 0.40_r8/)
! albedo of melting lakes due to puddling, open water, or white ice
! From D. Mironov (2010) Boreal Env. Research
! To revert albedo of melting lakes to the cold snow-free value, set
! lake_melt_icealb namelist to 0.60, 0.40 like alblak above.
real(r8), private :: alblakwi(numrad)
! Coefficient for calculating ice "fraction" for lake surface albedo
! From D. Mironov (2010) Boreal Env. Research
real(r8), parameter :: calb = 95.6_r8
!
! !PRIVATE DATA MEMBERS:
logical, private :: snowveg_affects_radiation = .true. ! Whether snow on the vegetation canopy affects the radiation/albedo calculations
!
! !PRIVATE DATA FUNCTIONS:
real(r8), allocatable, private :: albsat(:,:) ! wet soil albedo by color class and waveband (1=vis,2=nir)
real(r8), allocatable, private :: albdry(:,:) ! dry soil albedo by color class and waveband (1=vis,2=nir)
integer , allocatable, private :: isoicol(:) ! column soil color class
character(len=*), parameter, private :: sourcefile = &
__FILE__
!-----------------------------------------------------------------------
contains
!-----------------------------------------------------------------------
subroutine SurfaceAlbedo_readnl( NLFilename )
!
! !DESCRIPTION:
! Read the namelist for SurfaceAlbedo
!
! !USES:
use spmdMod , only : masterproc, mpicom
use fileutils , only : getavu, relavu, opnfil
use shr_nl_mod , only : shr_nl_find_group_name
use shr_mpi_mod , only : shr_mpi_bcast
!
! !ARGUMENTS:
character(len=*), intent(in) :: NLFilename ! Namelist filename
!
! !LOCAL VARIABLES:
integer :: ierr ! error code
integer :: unitn ! unit for namelist file
character(len=*), parameter :: nmlname = "surfacealbedo_inparm"
character(len=*), parameter :: subname = 'SurfaceAlbedo_readnl'
!-----------------------------------------------------------------------
namelist /surfacealbedo_inparm/ snowveg_affects_radiation
if (masterproc) then
unitn = getavu()
write(iulog,*) 'Read in '//nmlname//' namelist'
call opnfil (NLFilename, unitn, 'F')
call shr_nl_find_group_name(unitn, nmlname, status=ierr)
if (ierr == 0) then
read(unitn, nml=surfacealbedo_inparm, iostat=ierr)
if (ierr /= 0) then
call endrun(msg="ERROR reading "//nmlname//"namelist"//errmsg(sourcefile, __LINE__))
end if
else
call endrun(msg="ERROR could NOT find "//nmlname//"namelist"//errmsg(sourcefile, __LINE__))
end if
call relavu( unitn )
end if
call shr_mpi_bcast(snowveg_affects_radiation, mpicom)
if (masterproc) then
write(iulog,*)
write(iulog,*) nmlname, ' settings'
write(iulog,nml=surfacealbedo_inparm)
write(iulog,*)
end if
end subroutine SurfaceAlbedo_readnl
!-----------------------------------------------------------------------
subroutine SurfaceAlbedoInitTimeConst(bounds)
!
! !DESCRIPTION:
! Initialize module time constant variables
!
! !USES:
use shr_log_mod, only : errMsg => shr_log_errMsg
use fileutils , only : getfil
use abortutils , only : endrun
use ncdio_pio , only : file_desc_t, ncd_defvar, ncd_io, ncd_pio_openfile, ncd_pio_closefile
use spmdMod , only : masterproc
!
! !ARGUMENTS:
type(bounds_type), intent(in) :: bounds
!
! !LOCAL VARIABLES:
integer :: c,g ! indices
integer :: mxsoil_color ! maximum number of soil color classes
type(file_desc_t) :: ncid ! netcdf id
character(len=256) :: locfn ! local filename
integer :: ier ! error status
logical :: readvar
integer ,pointer :: soic2d (:) ! read in - soil color
!---------------------------------------------------------------------
! Allocate module variable for soil color
allocate(isoicol(bounds%begc:bounds%endc))
! Determine soil color and number of soil color classes
call getfil (fsurdat, locfn, 0)
call ncd_pio_openfile (ncid, locfn, 0)
call ncd_io(ncid=ncid, varname='mxsoil_color', flag='read', data=mxsoil_color, readvar=readvar)
if ( .not. readvar ) then
call endrun(msg=' ERROR: mxsoil_color NOT on surfdata file '//errMsg(sourcefile, __LINE__))
end if
allocate(soic2d(bounds%begg:bounds%endg))
call ncd_io(ncid=ncid, varname='SOIL_COLOR', flag='read', data=soic2d, dim1name=grlnd, readvar=readvar)
if (.not. readvar) then
call endrun(msg=' ERROR: SOIL_COLOR NOT on surfdata file'//errMsg(sourcefile, __LINE__))
end if
do c = bounds%begc, bounds%endc
g = col%gridcell(c)
isoicol(c) = soic2d(g)
end do
deallocate(soic2d)
call ncd_pio_closefile(ncid)
! Determine saturated and dry soil albedos for n color classes and
! numrad wavebands (1=vis, 2=nir)
allocate(albsat(mxsoil_color,numrad), albdry(mxsoil_color,numrad), stat=ier)
if (ier /= 0) then
write(iulog,*)'allocation error for albsat, albdry'
call endrun(msg=errMsg(sourcefile, __LINE__))
end if
if (masterproc) then
write(iulog,*) 'Attempting to read soil colo data .....'
end if
if (mxsoil_color == 8) then
albsat(1:8,1) = (/0.12_r8,0.11_r8,0.10_r8,0.09_r8,0.08_r8,0.07_r8,0.06_r8,0.05_r8/)
albsat(1:8,2) = (/0.24_r8,0.22_r8,0.20_r8,0.18_r8,0.16_r8,0.14_r8,0.12_r8,0.10_r8/)
albdry(1:8,1) = (/0.24_r8,0.22_r8,0.20_r8,0.18_r8,0.16_r8,0.14_r8,0.12_r8,0.10_r8/)
albdry(1:8,2) = (/0.48_r8,0.44_r8,0.40_r8,0.36_r8,0.32_r8,0.28_r8,0.24_r8,0.20_r8/)
else if (mxsoil_color == 20) then
albsat(1:20,1) = (/0.25_r8,0.23_r8,0.21_r8,0.20_r8,0.19_r8,0.18_r8,0.17_r8,0.16_r8,&
0.15_r8,0.14_r8,0.13_r8,0.12_r8,0.11_r8,0.10_r8,0.09_r8,0.08_r8,0.07_r8,0.06_r8,0.05_r8,0.04_r8/)
albsat(1:20,2) = (/0.50_r8,0.46_r8,0.42_r8,0.40_r8,0.38_r8,0.36_r8,0.34_r8,0.32_r8,&
0.30_r8,0.28_r8,0.26_r8,0.24_r8,0.22_r8,0.20_r8,0.18_r8,0.16_r8,0.14_r8,0.12_r8,0.10_r8,0.08_r8/)
albdry(1:20,1) = (/0.36_r8,0.34_r8,0.32_r8,0.31_r8,0.30_r8,0.29_r8,0.28_r8,0.27_r8,&
0.26_r8,0.25_r8,0.24_r8,0.23_r8,0.22_r8,0.20_r8,0.18_r8,0.16_r8,0.14_r8,0.12_r8,0.10_r8,0.08_r8/)
albdry(1:20,2) = (/0.61_r8,0.57_r8,0.53_r8,0.51_r8,0.49_r8,0.48_r8,0.45_r8,0.43_r8,&
0.41_r8,0.39_r8,0.37_r8,0.35_r8,0.33_r8,0.31_r8,0.29_r8,0.27_r8,0.25_r8,0.23_r8,0.21_r8,0.16_r8/)
else
write(iulog,*)'maximum color class = ',mxsoil_color,' is not supported'
call endrun(msg=errMsg(sourcefile, __LINE__))
end if
! Set alblakwi
alblakwi(:) = lake_melt_icealb(:)
end subroutine SurfaceAlbedoInitTimeConst
!-----------------------------------------------------------------------
subroutine SurfaceAlbedo(bounds,nc, &
num_nourbanc, filter_nourbanc, &
num_nourbanp, filter_nourbanp, &
num_urbanc , filter_urbanc, &
num_urbanp , filter_urbanp, &
nextsw_cday , declinp1, &
clm_fates, &
aerosol_inst, canopystate_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, &
lakestate_inst, temperature_inst, surfalb_inst)
!
! !DESCRIPTION:
! Surface albedo and two-stream fluxes
! Surface albedos. Also fluxes (per unit incoming direct and diffuse
! radiation) reflected, transmitted, and absorbed by vegetation.
! Calculate sunlit and shaded fluxes as described by
! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 and extended to
! a multi-layer canopy to calculate APAR profile
!
! The calling sequence is:
! -> SurfaceAlbedo: albedos for next time step
! -> SoilAlbedo: soil/lake/glacier/wetland albedos
! -> SNICAR_RT: snow albedos: direct beam (SNICAR)
! -> SNICAR_RT: snow albedos: diffuse (SNICAR)
! -> TwoStream: absorbed, reflected, transmitted solar fluxes (vis dir,vis dif, nir dir, nir dif)
!
! Note that this is called with the "inactive_and_active" version of the filters, because
! the variables computed here are needed over inactive points that might later become
! active (due to landuse change). Thus, this routine cannot depend on variables that are
! only computed over active points.
!
! !USES:
use shr_orb_mod
use clm_time_manager , only : get_nstep
use abortutils , only : endrun
use clm_varctl , only : use_subgrid_fluxes, use_snicar_frc, use_fates
use CLMFatesInterfaceMod, only : hlm_fates_interface_type
! !ARGUMENTS:
type(bounds_type) , intent(in) :: bounds ! bounds
integer , intent(in) :: nc ! clump index
integer , intent(in) :: num_nourbanc ! number of columns in non-urban filter
integer , intent(in) :: filter_nourbanc(:) ! column filter for non-urban points
integer , intent(in) :: num_nourbanp ! number of patches in non-urban filter
integer , intent(in) :: filter_nourbanp(:) ! patch filter for non-urban points
integer , intent(in) :: num_urbanc ! number of columns in urban filter
integer , intent(in) :: filter_urbanc(:) ! column filter for urban points
integer , intent(in) :: num_urbanp ! number of patches in urban filter
integer , intent(in) :: filter_urbanp(:) ! patch filter for rban points
real(r8) , intent(in) :: nextsw_cday ! calendar day at Greenwich (1.00, ..., days/year)
real(r8) , intent(in) :: declinp1 ! declination angle (radians) for next time step
type(hlm_fates_interface_type), intent(inout) :: clm_fates
type(aerosol_type) , intent(in) :: aerosol_inst
type(canopystate_type) , intent(in) :: canopystate_inst
type(waterstatebulk_type) , intent(in) :: waterstatebulk_inst
type(waterdiagnosticbulk_type) , intent(in) :: waterdiagnosticbulk_inst
type(lakestate_type) , intent(in) :: lakestate_inst
type(temperature_type) , intent(in) :: temperature_inst
type(surfalb_type) , intent(inout) :: surfalb_inst
!
! !LOCAL VARIABLES:
integer :: i ! index for layers [idx]
integer :: aer ! index for sno_nbr_aer
real(r8) :: extkn ! nitrogen allocation coefficient
integer :: fp,fc,g,c,p,iv ! indices
integer :: ib ! band index
integer :: ic ! 0=unit incoming direct; 1=unit incoming diffuse
real(r8) :: dinc ! lai+sai increment for canopy layer
real(r8) :: dincmax ! maximum lai+sai increment for canopy layer
real(r8) :: dincmax_sum ! cumulative sum of maximum lai+sai increment for canopy layer
real(r8) :: laisum ! sum of canopy layer lai for error check
real(r8) :: saisum ! sum of canopy layer sai for error check
integer :: flg_slr ! flag for SNICAR (=1 if direct, =2 if diffuse)
integer :: flg_snw_ice ! flag for SNICAR (=1 when called from CLM, =2 when called from sea-ice)
integer :: num_vegsol ! number of vegetated patches where coszen>0
integer :: num_novegsol ! number of vegetated patches where coszen>0
integer :: filter_vegsol (bounds%endp-bounds%begp+1) ! patch filter where vegetated and coszen>0
integer :: filter_novegsol (bounds%endp-bounds%begp+1) ! patch filter where vegetated and coszen>0
real(r8) :: wl (bounds%begp:bounds%endp) ! fraction of LAI+SAI that is LAI
real(r8) :: ws (bounds%begp:bounds%endp) ! fraction of LAI+SAI that is SAI
real(r8) :: blai(bounds%begp:bounds%endp) ! lai buried by snow: tlai - elai
real(r8) :: bsai(bounds%begp:bounds%endp) ! sai buried by snow: tsai - esai
real(r8) :: coszen_gcell (bounds%begg:bounds%endg) ! cosine solar zenith angle for next time step (grc)
real(r8) :: coszen_patch (bounds%begp:bounds%endp) ! cosine solar zenith angle for next time step (patch)
real(r8) :: rho(bounds%begp:bounds%endp,numrad) ! leaf/stem refl weighted by fraction LAI and SAI
real(r8) :: tau(bounds%begp:bounds%endp,numrad) ! leaf/stem tran weighted by fraction LAI and SAI
real(r8) :: h2osno_total (bounds%begc:bounds%endc) ! total snow water (mm H2O)
real(r8) :: albsfc (bounds%begc:bounds%endc,numrad) ! albedo of surface underneath snow (col,bnd)
real(r8) :: albsnd(bounds%begc:bounds%endc,numrad) ! snow albedo (direct)
real(r8) :: albsni(bounds%begc:bounds%endc,numrad) ! snow albedo (diffuse)
real(r8) :: albsnd_pur (bounds%begc:bounds%endc,numrad) ! direct pure snow albedo (radiative forcing)
real(r8) :: albsni_pur (bounds%begc:bounds%endc,numrad) ! diffuse pure snow albedo (radiative forcing)
real(r8) :: albsnd_bc (bounds%begc:bounds%endc,numrad) ! direct snow albedo without BC (radiative forcing)
real(r8) :: albsni_bc (bounds%begc:bounds%endc,numrad) ! diffuse snow albedo without BC (radiative forcing)
real(r8) :: albsnd_oc (bounds%begc:bounds%endc,numrad) ! direct snow albedo without OC (radiative forcing)
real(r8) :: albsni_oc (bounds%begc:bounds%endc,numrad) ! diffuse snow albedo without OC (radiative forcing)
real(r8) :: albsnd_dst (bounds%begc:bounds%endc,numrad) ! direct snow albedo without dust (radiative forcing)
real(r8) :: albsni_dst (bounds%begc:bounds%endc,numrad) ! diffuse snow albedo without dust (radiative forcing)
real(r8) :: flx_absd_snw (bounds%begc:bounds%endc,-nlevsno+1:1,numrad) ! flux absorption factor for just snow (direct) [frc]
real(r8) :: flx_absi_snw (bounds%begc:bounds%endc,-nlevsno+1:1,numrad) ! flux absorption factor for just snow (diffuse) [frc]
real(r8) :: foo_snw (bounds%begc:bounds%endc,-nlevsno+1:1,numrad) ! dummy array for forcing calls
real(r8) :: h2osno_liq (bounds%begc:bounds%endc,-nlevsno+1:0) ! liquid snow content (col,lyr) [kg m-2]
real(r8) :: h2osno_ice (bounds%begc:bounds%endc,-nlevsno+1:0) ! ice content in snow (col,lyr) [kg m-2]
integer :: snw_rds_in (bounds%begc:bounds%endc,-nlevsno+1:0) ! snow grain size sent to SNICAR (col,lyr) [microns]
real(r8) :: mss_cnc_aer_in_frc_pur (bounds%begc:bounds%endc,-nlevsno+1:0,sno_nbr_aer) ! mass concentration of aerosol species for forcing calculation (zero) (col,lyr,aer) [kg kg-1]
real(r8) :: mss_cnc_aer_in_frc_bc (bounds%begc:bounds%endc,-nlevsno+1:0,sno_nbr_aer) ! mass concentration of aerosol species for BC forcing (col,lyr,aer) [kg kg-1]
real(r8) :: mss_cnc_aer_in_frc_oc (bounds%begc:bounds%endc,-nlevsno+1:0,sno_nbr_aer) ! mass concentration of aerosol species for OC forcing (col,lyr,aer) [kg kg-1]
real(r8) :: mss_cnc_aer_in_frc_dst (bounds%begc:bounds%endc,-nlevsno+1:0,sno_nbr_aer) ! mass concentration of aerosol species for dust forcing (col,lyr,aer) [kg kg-1]
real(r8) :: mss_cnc_aer_in_fdb (bounds%begc:bounds%endc,-nlevsno+1:0,sno_nbr_aer) ! mass concentration of all aerosol species for feedback calculation (col,lyr,aer) [kg kg-1]
real(r8), parameter :: mpe = 1.e-06_r8 ! prevents overflow for division by zero
integer , parameter :: nband =numrad ! number of solar radiation waveband classes
!-----------------------------------------------------------------------
associate(&
rhol => pftcon%rhol , & ! Input: leaf reflectance: 1=vis, 2=nir
rhos => pftcon%rhos , & ! Input: stem reflectance: 1=vis, 2=nir
taul => pftcon%taul , & ! Input: leaf transmittance: 1=vis, 2=nir
taus => pftcon%taus , & ! Input: stem transmittance: 1=vis, 2=nir
tlai => canopystate_inst%tlai_patch , & ! Input: [real(r8) (:) ] one-sided leaf area index, no burying by snow
tsai => canopystate_inst%tsai_patch , & ! Input: [real(r8) (:) ] one-sided stem area index, no burying by snow
elai => canopystate_inst%elai_patch , & ! Input: [real(r8) (:) ] one-sided leaf area index with burying by snow
esai => canopystate_inst%esai_patch , & ! Input: [real(r8) (:) ] one-sided stem area index with burying by snow
frac_sno => waterdiagnosticbulk_inst%frac_sno_col , & ! Input: [real(r8) (:) ] fraction of ground covered by snow (0 to 1)
fcansno => waterdiagnosticbulk_inst%fcansno_patch , & ! Input: [real(r8) (:) ] fraction of canopy that is snow-covered (0 to 1)
h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water content (col,lyr) [kg/m2]
h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice lens content (col,lyr) [kg/m2]
snw_rds => waterdiagnosticbulk_inst%snw_rds_col , & ! Input: [real(r8) (:,:) ] snow grain radius (col,lyr) [microns]
mss_cnc_bcphi => aerosol_inst%mss_cnc_bcphi_col , & ! Input: [real(r8) (:,:) ] mass concentration of hydrophilic BC (col,lyr) [kg/kg]
mss_cnc_bcpho => aerosol_inst%mss_cnc_bcpho_col , & ! Input: [real(r8) (:,:) ] mass concentration of hydrophobic BC (col,lyr) [kg/kg]
mss_cnc_ocphi => aerosol_inst%mss_cnc_ocphi_col , & ! Input: [real(r8) (:,:) ] mass concentration of hydrophilic OC (col,lyr) [kg/kg]
mss_cnc_ocpho => aerosol_inst%mss_cnc_ocpho_col , & ! Input: [real(r8) (:,:) ] mass concentration of hydrophobic OC (col,lyr) [kg/kg]
mss_cnc_dst1 => aerosol_inst%mss_cnc_dst1_col , & ! Input: [real(r8) (:,:) ] mass concentration of dust aerosol species 1 (col,lyr) [kg/kg]
mss_cnc_dst2 => aerosol_inst%mss_cnc_dst2_col , & ! Input: [real(r8) (:,:) ] mass concentration of dust aerosol species 2 (col,lyr) [kg/kg]
mss_cnc_dst3 => aerosol_inst%mss_cnc_dst3_col , & ! Input: [real(r8) (:,:) ] mass concentration of dust aerosol species 3 (col,lyr) [kg/kg]
mss_cnc_dst4 => aerosol_inst%mss_cnc_dst4_col , & ! Input: [real(r8) (:,:) ] mass concentration of dust aerosol species 4 (col,lyr) [kg/kg]
fsun_z => surfalb_inst%fsun_z_patch , & ! Output: [real(r8) (:,:) ] sunlit fraction of canopy layer
tlai_z => surfalb_inst%tlai_z_patch , & ! Output: [real(r8) (:,:) ] tlai increment for canopy layer
tsai_z => surfalb_inst%tsai_z_patch , & ! Output: [real(r8) (:,:) ] tsai increment for canopy layer
vcmaxcintsun => surfalb_inst%vcmaxcintsun_patch , & ! Output: [real(r8) (:) ] leaf to canopy scaling coefficient, sunlit leaf vcmax
vcmaxcintsha => surfalb_inst%vcmaxcintsha_patch , & ! Output: [real(r8) (:) ] leaf to canopy scaling coefficient, shaded leaf vcmax
ncan => surfalb_inst%ncan_patch , & ! Output: [integer (:) ] number of canopy layers
nrad => surfalb_inst%nrad_patch , & ! Output: [integer (:) ] number of canopy layers, above snow for radiative transfer
coszen_col => surfalb_inst%coszen_col , & ! Output: [real(r8) (:) ] cosine of solar zenith angle
albgrd => surfalb_inst%albgrd_col , & ! Output: [real(r8) (:,:) ] ground albedo (direct)
albgri => surfalb_inst%albgri_col , & ! Output: [real(r8) (:,:) ] ground albedo (diffuse)
albsod => surfalb_inst%albsod_col , & ! Output: [real(r8) (:,:) ] direct-beam soil albedo (col,bnd) [frc]
albsoi => surfalb_inst%albsoi_col , & ! Output: [real(r8) (:,:) ] diffuse soil albedo (col,bnd) [frc]
albgrd_pur => surfalb_inst%albgrd_pur_col , & ! Output: [real(r8) (:,:) ] pure snow ground albedo (direct)
albgri_pur => surfalb_inst%albgri_pur_col , & ! Output: [real(r8) (:,:) ] pure snow ground albedo (diffuse)
albgrd_bc => surfalb_inst%albgrd_bc_col , & ! Output: [real(r8) (:,:) ] ground albedo without BC (direct)
albgri_bc => surfalb_inst%albgri_bc_col , & ! Output: [real(r8) (:,:) ] ground albedo without BC (diffuse)
albgrd_oc => surfalb_inst%albgrd_oc_col , & ! Output: [real(r8) (:,:) ] ground albedo without OC (direct)
albgri_oc => surfalb_inst%albgri_oc_col , & ! Output: [real(r8) (:,:) ] ground albedo without OC (diffuse)
albgrd_dst => surfalb_inst%albgrd_dst_col , & ! Output: [real(r8) (:,:) ] ground albedo without dust (direct)
albgri_dst => surfalb_inst%albgri_dst_col , & ! Output: [real(r8) (:,:) ] ground albedo without dust (diffuse)
albsnd_hst => surfalb_inst%albsnd_hst_col , & ! Output: [real(r8) (:,:) ] snow albedo, direct, for history files (col,bnd) [frc]
albsni_hst => surfalb_inst%albsni_hst_col , & ! Output: [real(r8) (:,:) ] snow ground albedo, diffuse, for history files (col,bnd) [frc]
albd => surfalb_inst%albd_patch , & ! Output: [real(r8) (:,:) ] surface albedo (direct)
albi => surfalb_inst%albi_patch , & ! Output: [real(r8) (:,:) ] surface albedo (diffuse)
albdSF => surfalb_inst%albdSF_patch , & ! Output: [real(r8) (:,:) ] diagnostic snow-free surface albedo (direct)
albiSF => surfalb_inst%albiSF_patch , & ! Output: [real(r8) (:,:) ] diagnostic snow-free surface albedo (diffuse)
fabd => surfalb_inst%fabd_patch , & ! Output: [real(r8) (:,:) ] flux absorbed by canopy per unit direct flux
fabd_sun => surfalb_inst%fabd_sun_patch , & ! Output: [real(r8) (:,:) ] flux absorbed by sunlit canopy per unit direct flux
fabd_sha => surfalb_inst%fabd_sha_patch , & ! Output: [real(r8) (:,:) ] flux absorbed by shaded canopy per unit direct flux
fabi => surfalb_inst%fabi_patch , & ! Output: [real(r8) (:,:) ] flux absorbed by canopy per unit diffuse flux
fabi_sun => surfalb_inst%fabi_sun_patch , & ! Output: [real(r8) (:,:) ] flux absorbed by sunlit canopy per unit diffuse flux
fabi_sha => surfalb_inst%fabi_sha_patch , & ! Output: [real(r8) (:,:) ] flux absorbed by shaded canopy per unit diffuse flux
ftdd => surfalb_inst%ftdd_patch , & ! Output: [real(r8) (:,:) ] down direct flux below canopy per unit direct flux
ftid => surfalb_inst%ftid_patch , & ! Output: [real(r8) (:,:) ] down diffuse flux below canopy per unit direct flux
ftii => surfalb_inst%ftii_patch , & ! Output: [real(r8) (:,:) ] down diffuse flux below canopy per unit diffuse flux
flx_absdv => surfalb_inst%flx_absdv_col , & ! Output: [real(r8) (:,:) ] direct flux absorption factor (col,lyr): VIS [frc]
flx_absdn => surfalb_inst%flx_absdn_col , & ! Output: [real(r8) (:,:) ] direct flux absorption factor (col,lyr): NIR [frc]
flx_absiv => surfalb_inst%flx_absiv_col , & ! Output: [real(r8) (:,:) ] diffuse flux absorption factor (col,lyr): VIS [frc]
flx_absin => surfalb_inst%flx_absin_col , & ! Output: [real(r8) (:,:) ] diffuse flux absorption factor (col,lyr): NIR [frc]
fabd_sun_z => surfalb_inst%fabd_sun_z_patch , & ! Output: [real(r8) (:,:) ] absorbed sunlit leaf direct PAR (per unit lai+sai) for each canopy layer
fabd_sha_z => surfalb_inst%fabd_sha_z_patch , & ! Output: [real(r8) (:,:) ] absorbed shaded leaf direct PAR (per unit lai+sai) for each canopy layer
fabi_sun_z => surfalb_inst%fabi_sun_z_patch , & ! Output: [real(r8) (:,:) ] absorbed sunlit leaf diffuse PAR (per unit lai+sai) for each canopy layer
fabi_sha_z => surfalb_inst%fabi_sha_z_patch & ! Output: [real(r8) (:,:) ] absorbed shaded leaf diffuse PAR (per unit lai+sai) for each canopy layer
)
! Cosine solar zenith angle for next time step
do g = bounds%begg,bounds%endg
coszen_gcell(g) = shr_orb_cosz (nextsw_cday, grc%lat(g), grc%lon(g), declinp1)
end do
do c = bounds%begc,bounds%endc
g = col%gridcell(c)
coszen_col(c) = coszen_gcell(g)
end do
do fp = 1,num_nourbanp
p = filter_nourbanp(fp)
g = patch%gridcell(p)
coszen_patch(p) = coszen_gcell(g)
end do
! Initialize output because solar radiation only done if coszen > 0
do ib = 1, numrad
do fc = 1,num_nourbanc
c = filter_nourbanc(fc)
albsod(c,ib) = 0._r8
albsoi(c,ib) = 0._r8
albgrd(c,ib) = 0._r8
albgri(c,ib) = 0._r8
albgrd_pur(c,ib) = 0._r8
albgri_pur(c,ib) = 0._r8
albgrd_bc(c,ib) = 0._r8
albgri_bc(c,ib) = 0._r8
albgrd_oc(c,ib) = 0._r8
albgri_oc(c,ib) = 0._r8
albgrd_dst(c,ib) = 0._r8
albgri_dst(c,ib) = 0._r8
do i=-nlevsno+1,1,1
flx_absdv(c,i) = 0._r8
flx_absdn(c,i) = 0._r8
flx_absiv(c,i) = 0._r8
flx_absin(c,i) = 0._r8
enddo
end do
do fp = 1,num_nourbanp
p = filter_nourbanp(fp)
albd(p,ib) = 1._r8
albi(p,ib) = 1._r8
if (use_SSRE) then
albdSF(p,ib) = 1._r8
albiSF(p,ib) = 1._r8
end if
fabd(p,ib) = 0._r8
fabd_sun(p,ib) = 0._r8
fabd_sha(p,ib) = 0._r8
fabi(p,ib) = 0._r8
fabi_sun(p,ib) = 0._r8
fabi_sha(p,ib) = 0._r8
ftdd(p,ib) = 0._r8
ftid(p,ib) = 0._r8
ftii(p,ib) = 0._r8
end do
end do ! end of numrad loop
! SoilAlbedo called before SNICAR_RT
! so that reflectance of soil beneath snow column is known
! ahead of time for snow RT calculation.
! Snow albedos
! Note that snow albedo routine will only compute nonzero snow albedos
! where h2osno> 0 and coszen > 0
! Ground surface albedos
! Note that ground albedo routine will only compute nonzero snow albedos
! where coszen > 0
call SoilAlbedo(bounds, &
num_nourbanc, filter_nourbanc, &
coszen_col(bounds%begc:bounds%endc), &
albsnd(bounds%begc:bounds%endc, :), &
albsni(bounds%begc:bounds%endc, :), &
lakestate_inst, temperature_inst, waterstatebulk_inst, surfalb_inst)
! set variables to pass to SNICAR.
flg_snw_ice = 1 ! calling from CLM, not CSIM
do c=bounds%begc,bounds%endc
albsfc(c,:) = albsoi(c,:)
h2osno_liq(c,:) = h2osoi_liq(c,-nlevsno+1:0)
h2osno_ice(c,:) = h2osoi_ice(c,-nlevsno+1:0)
snw_rds_in(c,:) = nint(snw_rds(c,:))
end do
! zero aerosol input arrays
do aer = 1, sno_nbr_aer
do i = -nlevsno+1, 0
do c = bounds%begc, bounds%endc
mss_cnc_aer_in_frc_pur(c,i,aer) = 0._r8
mss_cnc_aer_in_frc_bc(c,i,aer) = 0._r8
mss_cnc_aer_in_frc_oc(c,i,aer) = 0._r8
mss_cnc_aer_in_frc_dst(c,i,aer) = 0._r8
mss_cnc_aer_in_fdb(c,i,aer) = 0._r8
end do
end do
end do
! Set aerosol input arrays
! feedback input arrays have been zeroed
! set soot and dust aerosol concentrations:
if (DO_SNO_AER) then
mss_cnc_aer_in_fdb(bounds%begc:bounds%endc,:,1) = mss_cnc_bcphi(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_fdb(bounds%begc:bounds%endc,:,2) = mss_cnc_bcpho(bounds%begc:bounds%endc,:)
! DO_SNO_OC is set in SNICAR_varpar. Default case is to ignore OC concentrations because:
! 1) Knowledge of their optical properties is primitive
! 2) When 'water-soluble' OPAC optical properties are applied to OC in snow,
! it has a negligible darkening effect.
if (DO_SNO_OC) then
mss_cnc_aer_in_fdb(bounds%begc:bounds%endc,:,3) = mss_cnc_ocphi(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_fdb(bounds%begc:bounds%endc,:,4) = mss_cnc_ocpho(bounds%begc:bounds%endc,:)
endif
mss_cnc_aer_in_fdb(bounds%begc:bounds%endc,:,5) = mss_cnc_dst1(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_fdb(bounds%begc:bounds%endc,:,6) = mss_cnc_dst2(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_fdb(bounds%begc:bounds%endc,:,7) = mss_cnc_dst3(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_fdb(bounds%begc:bounds%endc,:,8) = mss_cnc_dst4(bounds%begc:bounds%endc,:)
endif
call waterstatebulk_inst%CalculateTotalH2osno(bounds, num_nourbanc, filter_nourbanc, &
caller = 'SurfaceAlbedo', &
h2osno_total = h2osno_total(bounds%begc:bounds%endc))
! If radiative forcing is being calculated, first estimate clean-snow albedo
if (use_snicar_frc) then
! 1. BC input array:
! set dust and (optionally) OC concentrations, so BC_FRC=[(BC+OC+dust)-(OC+dust)]
mss_cnc_aer_in_frc_bc(bounds%begc:bounds%endc,:,5) = mss_cnc_dst1(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_frc_bc(bounds%begc:bounds%endc,:,6) = mss_cnc_dst2(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_frc_bc(bounds%begc:bounds%endc,:,7) = mss_cnc_dst3(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_frc_bc(bounds%begc:bounds%endc,:,8) = mss_cnc_dst4(bounds%begc:bounds%endc,:)
if (DO_SNO_OC) then
mss_cnc_aer_in_frc_bc(bounds%begc:bounds%endc,:,3) = mss_cnc_ocphi(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_frc_bc(bounds%begc:bounds%endc,:,4) = mss_cnc_ocpho(bounds%begc:bounds%endc,:)
endif
! BC FORCING CALCULATIONS
flg_slr = 1; ! direct-beam
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
coszen_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
h2osno_total(bounds%begc:bounds%endc), &
snw_rds_in(bounds%begc:bounds%endc, :), &
mss_cnc_aer_in_frc_bc(bounds%begc:bounds%endc, :, :), &
albsfc(bounds%begc:bounds%endc, :), &
albsnd_bc(bounds%begc:bounds%endc, :), &
foo_snw(bounds%begc:bounds%endc, :, :), &
waterdiagnosticbulk_inst)
flg_slr = 2; ! diffuse
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
coszen_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
h2osno_total(bounds%begc:bounds%endc), &
snw_rds_in(bounds%begc:bounds%endc, :), &
mss_cnc_aer_in_frc_bc(bounds%begc:bounds%endc, :, :), &
albsfc(bounds%begc:bounds%endc, :), &
albsni_bc(bounds%begc:bounds%endc, :), &
foo_snw(bounds%begc:bounds%endc, :, :), &
waterdiagnosticbulk_inst)
! 2. OC input array:
! set BC and dust concentrations, so OC_FRC=[(BC+OC+dust)-(BC+dust)]
if (DO_SNO_OC) then
mss_cnc_aer_in_frc_oc(bounds%begc:bounds%endc,:,1) = mss_cnc_bcphi(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_frc_oc(bounds%begc:bounds%endc,:,2) = mss_cnc_bcpho(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_frc_oc(bounds%begc:bounds%endc,:,5) = mss_cnc_dst1(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_frc_oc(bounds%begc:bounds%endc,:,6) = mss_cnc_dst2(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_frc_oc(bounds%begc:bounds%endc,:,7) = mss_cnc_dst3(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_frc_oc(bounds%begc:bounds%endc,:,8) = mss_cnc_dst4(bounds%begc:bounds%endc,:)
! OC FORCING CALCULATIONS
flg_slr = 1; ! direct-beam
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
coszen_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
h2osno_total(bounds%begc:bounds%endc), &
snw_rds_in(bounds%begc:bounds%endc, :), &
mss_cnc_aer_in_frc_oc(bounds%begc:bounds%endc, :, :), &
albsfc(bounds%begc:bounds%endc, :), &
albsnd_oc(bounds%begc:bounds%endc, :), &
foo_snw(bounds%begc:bounds%endc, :, :), &
waterdiagnosticbulk_inst)
flg_slr = 2; ! diffuse
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
coszen_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
h2osno_total(bounds%begc:bounds%endc), &
snw_rds_in(bounds%begc:bounds%endc, :), &
mss_cnc_aer_in_frc_oc(bounds%begc:bounds%endc, :, :), &
albsfc(bounds%begc:bounds%endc, :), &
albsni_oc(bounds%begc:bounds%endc, :), &
foo_snw(bounds%begc:bounds%endc, :, :), &
waterdiagnosticbulk_inst)
endif
! 3. DUST input array:
! set BC and OC concentrations, so DST_FRC=[(BC+OC+dust)-(BC+OC)]
mss_cnc_aer_in_frc_dst(bounds%begc:bounds%endc,:,1) = mss_cnc_bcphi(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_frc_dst(bounds%begc:bounds%endc,:,2) = mss_cnc_bcpho(bounds%begc:bounds%endc,:)
if (DO_SNO_OC) then
mss_cnc_aer_in_frc_dst(bounds%begc:bounds%endc,:,3) = mss_cnc_ocphi(bounds%begc:bounds%endc,:)
mss_cnc_aer_in_frc_dst(bounds%begc:bounds%endc,:,4) = mss_cnc_ocpho(bounds%begc:bounds%endc,:)
endif
! DUST FORCING CALCULATIONS
flg_slr = 1; ! direct-beam
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
coszen_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
h2osno_total(bounds%begc:bounds%endc), &
snw_rds_in(bounds%begc:bounds%endc, :), &
mss_cnc_aer_in_frc_dst(bounds%begc:bounds%endc, :, :), &
albsfc(bounds%begc:bounds%endc, :), &
albsnd_dst(bounds%begc:bounds%endc, :), &
foo_snw(bounds%begc:bounds%endc, :, :), &
waterdiagnosticbulk_inst)
flg_slr = 2; ! diffuse
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
coszen_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
h2osno_total(bounds%begc:bounds%endc), &
snw_rds_in(bounds%begc:bounds%endc, :), &
mss_cnc_aer_in_frc_dst(bounds%begc:bounds%endc, :, :), &
albsfc(bounds%begc:bounds%endc, :), &
albsni_dst(bounds%begc:bounds%endc, :), &
foo_snw(bounds%begc:bounds%endc, :, :), &
waterdiagnosticbulk_inst)
! 4. ALL AEROSOL FORCING CALCULATION
! (pure snow albedo)
flg_slr = 1; ! direct-beam
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
coszen_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
h2osno_total(bounds%begc:bounds%endc), &
snw_rds_in(bounds%begc:bounds%endc, :), &
mss_cnc_aer_in_frc_pur(bounds%begc:bounds%endc, :, :), &
albsfc(bounds%begc:bounds%endc, :), &
albsnd_pur(bounds%begc:bounds%endc, :), &
foo_snw(bounds%begc:bounds%endc, :, :), &
waterdiagnosticbulk_inst)
flg_slr = 2; ! diffuse
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
coszen_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
h2osno_total(bounds%begc:bounds%endc), &
snw_rds_in(bounds%begc:bounds%endc, :), &
mss_cnc_aer_in_frc_pur(bounds%begc:bounds%endc, :, :), &
albsfc(bounds%begc:bounds%endc, :), &
albsni_pur(bounds%begc:bounds%endc, :), &
foo_snw(bounds%begc:bounds%endc, :, :), &
waterdiagnosticbulk_inst)
end if
! CLIMATE FEEDBACK CALCULATIONS, ALL AEROSOLS:
flg_slr = 1; ! direct-beam
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
coszen_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
h2osno_total(bounds%begc:bounds%endc), &
snw_rds_in(bounds%begc:bounds%endc, :), &
mss_cnc_aer_in_fdb(bounds%begc:bounds%endc, :, :), &
albsfc(bounds%begc:bounds%endc, :), &
albsnd(bounds%begc:bounds%endc, :), &
flx_absd_snw(bounds%begc:bounds%endc, :, :), &
waterdiagnosticbulk_inst)
flg_slr = 2; ! diffuse
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
coszen_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
h2osno_total(bounds%begc:bounds%endc), &
snw_rds_in(bounds%begc:bounds%endc, :), &
mss_cnc_aer_in_fdb(bounds%begc:bounds%endc, :, :), &
albsfc(bounds%begc:bounds%endc, :), &
albsni(bounds%begc:bounds%endc, :), &
flx_absi_snw(bounds%begc:bounds%endc, :, :), &
waterdiagnosticbulk_inst)
! ground albedos and snow-fraction weighting of snow absorption factors
do ib = 1, nband
do fc = 1,num_nourbanc
c = filter_nourbanc(fc)
if (coszen_col(c) > 0._r8) then
! ground albedo was originally computed in SoilAlbedo, but is now computed here
! because the order of SoilAlbedo and SNICAR_RT was switched for SNICAR.
albgrd(c,ib) = albsod(c,ib)*(1._r8-frac_sno(c)) + albsnd(c,ib)*frac_sno(c)
albgri(c,ib) = albsoi(c,ib)*(1._r8-frac_sno(c)) + albsni(c,ib)*frac_sno(c)
! albedos for radiative forcing calculations:
if (use_snicar_frc) then
! BC forcing albedo
albgrd_bc(c,ib) = albsod(c,ib)*(1.-frac_sno(c)) + albsnd_bc(c,ib)*frac_sno(c)
albgri_bc(c,ib) = albsoi(c,ib)*(1.-frac_sno(c)) + albsni_bc(c,ib)*frac_sno(c)
if (DO_SNO_OC) then
! OC forcing albedo
albgrd_oc(c,ib) = albsod(c,ib)*(1.-frac_sno(c)) + albsnd_oc(c,ib)*frac_sno(c)
albgri_oc(c,ib) = albsoi(c,ib)*(1.-frac_sno(c)) + albsni_oc(c,ib)*frac_sno(c)
endif
! dust forcing albedo
albgrd_dst(c,ib) = albsod(c,ib)*(1.-frac_sno(c)) + albsnd_dst(c,ib)*frac_sno(c)
albgri_dst(c,ib) = albsoi(c,ib)*(1.-frac_sno(c)) + albsni_dst(c,ib)*frac_sno(c)
! pure snow albedo for all-aerosol radiative forcing
albgrd_pur(c,ib) = albsod(c,ib)*(1.-frac_sno(c)) + albsnd_pur(c,ib)*frac_sno(c)
albgri_pur(c,ib) = albsoi(c,ib)*(1.-frac_sno(c)) + albsni_pur(c,ib)*frac_sno(c)
end if
! also in this loop (but optionally in a different loop for vectorized code)
! weight snow layer radiative absorption factors based on snow fraction and soil albedo
! (NEEDED FOR ENERGY CONSERVATION)
do i = -nlevsno+1,1,1
if (.not. use_subgrid_fluxes .or. lun%itype(col%landunit(c)) == istdlak) then
if (ib == 1) then
flx_absdv(c,i) = flx_absd_snw(c,i,ib)*frac_sno(c) + &
((1.-frac_sno(c))*(1-albsod(c,ib))*(flx_absd_snw(c,i,ib)/(1.-albsnd(c,ib))))
flx_absiv(c,i) = flx_absi_snw(c,i,ib)*frac_sno(c) + &
((1.-frac_sno(c))*(1-albsoi(c,ib))*(flx_absi_snw(c,i,ib)/(1.-albsni(c,ib))))
elseif (ib == 2) then
flx_absdn(c,i) = flx_absd_snw(c,i,ib)*frac_sno(c) + &
((1.-frac_sno(c))*(1-albsod(c,ib))*(flx_absd_snw(c,i,ib)/(1.-albsnd(c,ib))))
flx_absin(c,i) = flx_absi_snw(c,i,ib)*frac_sno(c) + &
((1.-frac_sno(c))*(1-albsoi(c,ib))*(flx_absi_snw(c,i,ib)/(1.-albsni(c,ib))))
endif
else
if (ib == 1) then
flx_absdv(c,i) = flx_absd_snw(c,i,ib)
flx_absiv(c,i) = flx_absi_snw(c,i,ib)
elseif (ib == 2) then
flx_absdn(c,i) = flx_absd_snw(c,i,ib)
flx_absin(c,i) = flx_absi_snw(c,i,ib)
endif
endif
enddo
endif
enddo
enddo
! For diagnostics, set snow albedo to spval over non-snow non-urban points
! so that it is not averaged in history buffer (OPTIONAL)
! TODO - this is set to 0 not spval - seems wrong since it will be averaged in
do ib = 1, nband
do fc = 1,num_nourbanc
c = filter_nourbanc(fc)
if ((coszen_col(c) > 0._r8) .and. (h2osno_total(c) > 0._r8)) then
albsnd_hst(c,ib) = albsnd(c,ib)
albsni_hst(c,ib) = albsni(c,ib)
else
albsnd_hst(c,ib) = 0._r8
albsni_hst(c,ib) = 0._r8
endif
enddo
enddo
! Create solar-vegetated filter for the following calculations
num_vegsol = 0
num_novegsol = 0
do fp = 1,num_nourbanp
p = filter_nourbanp(fp)
if (coszen_patch(p) > 0._r8) then
if ((lun%itype(patch%landunit(p)) == istsoil .or. &
lun%itype(patch%landunit(p)) == istcrop ) &
.and. (elai(p) + esai(p)) > 0._r8) then
num_vegsol = num_vegsol + 1
filter_vegsol(num_vegsol) = p
else
num_novegsol = num_novegsol + 1
filter_novegsol(num_novegsol) = p
end if
end if
end do
! Weight reflectance/transmittance by lai and sai
! Only perform on vegetated patches where coszen > 0
do fp = 1,num_vegsol
p = filter_vegsol(fp)
wl(p) = elai(p) / max( elai(p)+esai(p), mpe )
ws(p) = esai(p) / max( elai(p)+esai(p), mpe )
end do
do ib = 1, numrad
do fp = 1,num_vegsol
p = filter_vegsol(fp)
rho(p,ib) = max( rhol(patch%itype(p),ib)*wl(p) + rhos(patch%itype(p),ib)*ws(p), mpe )
tau(p,ib) = max( taul(patch%itype(p),ib)*wl(p) + taus(patch%itype(p),ib)*ws(p), mpe )
end do
end do
! Diagnose number of canopy layers for radiative transfer, in increments of dincmax.
! Add to number of layers so long as cumulative leaf+stem area does not exceed total
! leaf+stem area. Then add any remaining leaf+stem area to next layer and exit the loop.
! Do this first for elai and esai (not buried by snow) and then for the part of the
! canopy that is buried by snow.
! ------------------
! tlai_z = leaf area increment for a layer
! tsai_z = stem area increment for a layer
! nrad = number of canopy layers above snow
! ncan = total number of canopy layers
!
! tlai_z summed from 1 to nrad = elai
! tlai_z summed from 1 to ncan = tlai
! tsai_z summed from 1 to nrad = esai
! tsai_z summed from 1 to ncan = tsai
! ------------------
!
! Canopy layering needs to be done for all "num_nourbanp" not "num_vegsol"
! because layering is needed for all time steps regardless of radiation
!
! Sun/shade big leaf code uses only one layer (nrad = ncan = 1), triggered by
! nlevcan = 1
dincmax = 0.25_r8
do fp = 1,num_nourbanp
p = filter_nourbanp(fp)
if (nlevcan == 1) then
nrad(p) = 1
ncan(p) = 1
tlai_z(p,1) = elai(p)
tsai_z(p,1) = esai(p)
else if (nlevcan > 1) then
if (elai(p)+esai(p) == 0._r8) then
nrad(p) = 0
else
dincmax_sum = 0._r8
do iv = 1, nlevcan
dincmax_sum = dincmax_sum + dincmax
if (((elai(p)+esai(p))-dincmax_sum) > 1.e-06_r8) then
nrad(p) = iv
dinc = dincmax
tlai_z(p,iv) = dinc * elai(p) / max(elai(p)+esai(p), mpe)
tsai_z(p,iv) = dinc * esai(p) / max(elai(p)+esai(p), mpe)
else
nrad(p) = iv
dinc = dincmax - (dincmax_sum - (elai(p)+esai(p)))
tlai_z(p,iv) = dinc * elai(p) / max(elai(p)+esai(p), mpe)
tsai_z(p,iv) = dinc * esai(p) / max(elai(p)+esai(p), mpe)
exit
end if
end do
! Mimumum of 4 canopy layers
if (nrad(p) < 4) then
nrad(p) = 4
do iv = 1, nrad(p)
tlai_z(p,iv) = elai(p) / nrad(p)
tsai_z(p,iv) = esai(p) / nrad(p)
end do
end if
end if
end if
! Error check: make sure cumulative of increments does not exceed total
laisum = 0._r8
saisum = 0._r8
do iv = 1, nrad(p)
laisum = laisum + tlai_z(p,iv)
saisum = saisum + tsai_z(p,iv)
end do
if (abs(laisum-elai(p)) > 1.e-06_r8 .or. abs(saisum-esai(p)) > 1.e-06_r8) then
write (iulog,*) 'multi-layer canopy error 01 in SurfaceAlbedo: ',&
nrad(p),elai(p),laisum,esai(p),saisum
call endrun(subgrid_index=p, subgrid_level=subgrid_level_patch, msg=errmsg(sourcefile, __LINE__))
end if
! Repeat to find canopy layers buried by snow
if (nlevcan > 1) then
blai(p) = tlai(p) - elai(p)
bsai(p) = tsai(p) - esai(p)
if (blai(p)+bsai(p) == 0._r8) then
ncan(p) = nrad(p)
else
dincmax_sum = 0._r8
do iv = nrad(p)+1, nlevcan
dincmax_sum = dincmax_sum + dincmax
if (((blai(p)+bsai(p))-dincmax_sum) > 1.e-06_r8) then
ncan(p) = iv
dinc = dincmax
tlai_z(p,iv) = dinc * blai(p) / max(blai(p)+bsai(p), mpe)
tsai_z(p,iv) = dinc * bsai(p) / max(blai(p)+bsai(p), mpe)
else
ncan(p) = iv
dinc = dincmax - (dincmax_sum - (blai(p)+bsai(p)))
tlai_z(p,iv) = dinc * blai(p) / max(blai(p)+bsai(p), mpe)
tsai_z(p,iv) = dinc * bsai(p) / max(blai(p)+bsai(p), mpe)
exit
end if
end do
end if
! Error check: make sure cumulative of increments does not exceed total
laisum = 0._r8
saisum = 0._r8
do iv = 1, ncan(p)
laisum = laisum + tlai_z(p,iv)
saisum = saisum + tsai_z(p,iv)
end do
if (abs(laisum-tlai(p)) > 1.e-06_r8 .or. abs(saisum-tsai(p)) > 1.e-06_r8) then
write (iulog,*) 'multi-layer canopy error 02 in SurfaceAlbedo: ',nrad(p),ncan(p)
write (iulog,*) tlai(p),elai(p),blai(p),laisum,tsai(p),esai(p),bsai(p),saisum
call endrun(subgrid_index=p, subgrid_level=subgrid_level_patch, msg=errmsg(sourcefile, __LINE__))
end if
end if
end do
! Zero fluxes for active canopy layers
do fp = 1,num_nourbanp
p = filter_nourbanp(fp)
do iv = 1, nrad(p)
fabd_sun_z(p,iv) = 0._r8
fabd_sha_z(p,iv) = 0._r8
fabi_sun_z(p,iv) = 0._r8
fabi_sha_z(p,iv) = 0._r8
fsun_z(p,iv) = 0._r8
end do
end do
! Default leaf to canopy scaling coefficients, used when coszen <= 0.
! This is the leaf nitrogen profile integrated over the full canopy.
! Integrate exp(-kn*x) over x=0 to x=elai and assign to shaded canopy,
! because sunlit fraction is 0. Canopy scaling coefficients are set in
! TwoStream for coszen > 0. So kn must be set here and in TwoStream.
extkn = 0.30_r8
do fp = 1,num_nourbanp
p = filter_nourbanp(fp)
if (nlevcan == 1) then
vcmaxcintsun(p) = 0._r8
vcmaxcintsha(p) = (1._r8 - exp(-extkn*elai(p))) / extkn
if (elai(p) > 0._r8) then
vcmaxcintsha(p) = vcmaxcintsha(p) / elai(p)
else
vcmaxcintsha(p) = 0._r8
end if
else if (nlevcan > 1) then
vcmaxcintsun(p) = 0._r8
vcmaxcintsha(p) = 0._r8
end if
end do
! Calculate surface albedos and fluxes
! Only perform on vegetated pfts where coszen > 0
if (use_fates) then
call clm_fates%wrap_canopy_radiation(bounds, nc, &
num_vegsol, filter_vegsol, &