BUGSrad
 All Classes Files Functions Variables
two_rt_sw_ocastrndm.f90
Go to the documentation of this file.
1 
2 
3 ! CVS: $Id: two_rt_sw_ocastrndm.F90,v 1.2 2006/11/16 21:11:23 norm Exp $
4 ! CVS: $Name: $
5 
6 !NBW This has been modified to apply descaling to the
7 !NBW delta-M scaled downwelling direct and diffuse fluxes
8 
9 subroutine two_rt_sw_gsolap &
10  ( ncol , nlm , mbs , ib , &
11  slr , amu0 , wc_clr , asym_clr , &
12  tau_clr , wc_cld , asym_cld , tau_cld , &
13  asdir , asdif , fddir_clr , fddif_clr , &
14  fudif_clr , fddir_all , fddif_all , fudif_all , &
15  cldamt , c_maximal , cf_max, cf_random )
16 
17  use kinds
18  implicit none
19 
20 !------------------------------------------------------------------------------
21 ! REFERENCES
22 ! Alternate cloud overlap calculation
23 !
24 ! SUBROUTINES CALLED:
25 !
26 ! FUNCTIONS CALLED:
27 !
28 ! INCLUDED COMMONS:
29 !
30 ! ARGUMENT LIST
31 ! Inputs
32  integer (kind=int_kind), intent(in):: &
33  ncol, & !Length of sub-domain
34  nlm, & !Number of layers
35  mbs, & !Number of SW spectral intervals
36  ib !Index of spectral interval
37 
38  real (kind=dbl_kind), intent(in), dimension(ncol):: &
39  slr, & !Fraction of daylight
40  amu0, & !Cosine of solar zenith angle
41  c_maximal !Maximally overlapped fraction
42 
43  real (kind=dbl_kind), intent(in), dimension(ncol,mbs):: &
44  asdir, & !Spectral direct surface albedo
45  asdif !Spectral diffuse serface albedo
46 
47  real (kind=dbl_kind), intent(in), dimension(ncol,nlm):: &
48  wc_clr, & !Clear-sky single scattering albedo
49  wc_cld, & !Cloudy-sky single scattering albedo
50  asym_clr, & !Clear-sky asymmetry parameter
51  asym_cld, & !Cloudy-sky asymmetry parameter
52  tau_clr, & !Clear-sky optical depth
53  tau_cld, & !Cloudy-sky optical depth
54  cldamt, & !Cloud fraction
55  cf_max, & !Maximally overlapped cloud fractions
56  cf_random !Randomly overlapped cloud fractions
57 
58 ! Outputs
59  real(kind=dbl_kind), intent(out), dimension(ncol, nlm+1):: &
60  fddir_clr, & !Clear-sky spectral direct downward flux (W/m^2)
61  fddir_all, & !All-sky spectral direct downward flux (W/m^2)
62  fddif_clr, & !Clear-sky spectral diffuse downward flux (W/m^2)
63  fddif_all, & !All-sky spectral diffuse downward flux (W/m^2)
64  fudif_clr, & !Clear-sky spectral diffuse upward flux (W/m^2)
65  fudif_all !All-sky spectral diffuse upward flux (W/m^2)
66 
67 
68 ! LOCAL VARIABLES
69  integer (kind=int_kind) :: &
70  i, & !Horizontal index
71  l !Vertical index
72 
73  real (kind=dbl_kind), parameter :: &
74  eps = 1.0e-2
75 
76  real (kind=dbl_kind) :: &
77  aa, & !
78  bb, & !
79  cc, & !
80  denom, & !
81  eggtau, & !
82  g3, & !
83  g4, & !
84  ggtau, & !
85  kappa, & !
86  r, & !
87  rinf, & !
88  t, & !
89  oms, & !
90  taus, & !
91  fact, & !
92  asy, & !
93  prop, & !
94  fd_total !
95 
96  real (kind=dbl_kind), dimension(nlm) :: &
97  exptau_s_clr, & !
98  exptau_s_cld, & !
99  exptau_us_clr,& !
100  exptau_us_cld,&
101  rr_clr, & !
102  rr_cld, & !
103  rr_tmp, & !
104  tr_clr, & !
105  tr_cld, & !
106  tr_tmp, & !
107  sigu, & !
108  sigd, & !
109  sigu_clr_fract,& !
110  sigu_cld_fract,& !
111  sigd_clr_fract,& !
112  sigd_cld_fract,& !
113  sigu_tmp_fract, & !
114  sigd_tmp_fract, & !
115  td, & !
116  vu !
117 
118  real (kind=dbl_kind), dimension(nlm+1) :: &
119  direct_s, & !
120  direct_us, &
121  re, & !
122  vd !
123 
124  real (kind=dbl_kind), dimension(nlm+1):: &
125  fddir_max, & !Maximum-overlap fraction spectral direct downward flux (W/m^2)
126  fddir_rndm, & !Random-overlap fraction spectral direct downward flux (W/m^2)
127  fddif_max, & !Maximum-overlap fraction spectral diffuse downward flux (W/m^2)
128  fddif_rndm, & !Random-overlap fraction spectral diffuse downward flux (W/m^2)
129  fudif_max, & !Maximum-overlap fraction spectral diffuse upward flux (W/m^2)
130  fudif_rndm !Random-overlap fraction spectral diffuse upward flux (W/m^2)
131 
132 
133  fddir_clr = 0.
134  fddir_all = 0.
135  fddif_clr = 0.
136  fddif_all = 0.
137  fudif_clr = 0.
138  fudif_all = 0.
139 
140  fddir_max = 0.
141  fddir_rndm = 0.
142  fddif_max = 0.
143  fddif_rndm = 0.
144  fudif_max = 0.
145  fudif_rndm = 0.
146 
147 ! Begin loop over all columns
148  do i = 1, ncol
149 
150  !Clear layer properties
151  direct_s(1) = 1.
152  direct_us(1) = 1.
153  do l = 1, nlm
154  !Delta-M on
155  fact = asym_clr(i,l)*asym_clr(i,l)
156  asy = asym_clr(i,l)/(1. + asym_clr(i,l))
157  !Delta-M off
158  !fact = 0.
159  !asy = asym_clr(i,l)
160  !End delta-M on/off
161  oms = ((1. - fact)*wc_clr(i,l))/(1. - fact*wc_clr(i,l))
162  taus = (1. - fact*wc_clr(i,l))*tau_clr(i,l)
163 
164 
165  exptau_s_clr(l) = exp(-taus/amu0(i))
166  direct_s(l+1) = exptau_s_clr(l)*direct_s(l)
167  exptau_us_clr(l) = exp(-tau_clr(i,l))
168  direct_us(l+1) = exptau_us_clr(l)*direct_us(l)
169  ! Local Eddington coefficients
170  t = 0.25*(7. - oms*(4. + 3.*asy))
171  r = -0.25*(1. - oms*(4. - 3.*asy))
172  kappa = sqrt(t*t - r*r)
173  rinf = r/(kappa+t)
174  ggtau = kappa*taus
175  eggtau = exp(-ggtau)
176  denom = (1. - rinf*rinf*eggtau*eggtau)
177  tr_clr(l) = (1. - rinf*rinf)*eggtau/denom
178  rr_clr(l) = rinf*(1. - eggtau*eggtau)/denom
179 
180  fact = kappa*kappa - 1./(amu0(i)*amu0(i))
181  if (abs(fact) .lt. eps) then
182  fact = 1./eps
183  else
184  fact = 1./fact
185  endif
186 
187  cc = oms*slr(i)*fact
188  g3 = 0.5 - 0.75*asy*amu0(i)
189  g4 = 1. - g3
190  aa = g3*(t - 1./amu0(i)) + g4*r
191  bb = g4*(t + 1./amu0(i)) + g3*r
192  sigu_clr_fract(l) = cc*((aa - rr_clr(l)*bb) - aa*tr_clr(l)*exptau_s_clr(l))
193  sigd_clr_fract(l) = cc*(-bb*tr_clr(l) + (bb-rr_clr(l)*aa)*exptau_s_clr(l))
194  enddo
195 
196  !CLEAR-SKY adding
197  re(1) = 0.
198  vd(1) = 0.
199  do l = 1, nlm
200  prop = 1./(1. - re(l)*rr_clr(l))
201  re(l+1) = rr_clr(l) + tr_clr(l)*tr_clr(l)*re(l)*prop
202  vd(l+1) = sigd_clr_fract(l)*direct_s(l) + (tr_clr(l)*vd(l) + tr_clr(l)*re(l)*sigu_clr_fract(l)*direct_s(l))*prop
203  vu(l) = (rr_clr(l)*vd(l) + sigu_clr_fract(l)*direct_s(l))*prop
204  td(l) = prop
205  enddo
206 
207  !CLEAR-SKY diffuse flux calculation
208  fddif_clr(i,1) = 0.
209  fudif_clr(i,nlm+1) = (asdif(i,ib)*vd(nlm+1) + &
210  asdir(i,ib)*slr(i)*amu0(i)*direct_s(nlm+1))/ &
211  (1. - asdif(i,ib)*re(nlm+1))
212  do l = nlm+1, 2, -1
213  fddif_clr(i,l) = re(l)*fudif_clr(i,l) + vd(l)
214  fudif_clr(i,l-1) = tr_clr(l-1)*fudif_clr(i,l)*td(l-1) + vu(l-1)
215  enddo
216 
217  !CLEAR-SKY direct flux calculation, with delta-M descaling
218  fddir_clr(i,1) = amu0(i)*slr(i)
219  do l = 1, nlm
220  fd_total = fddif_clr(i,1)+amu0(i)*slr(i)*direct_s(l+1)
221  fddir_clr(i,l+1) = amu0(i)*slr(i)*direct_us(l+1)
222  fddif_clr(i,l+1) = fd_total - fddir_clr(i,l+1)
223  enddo
224  !END CLEAR-SKY CALC
225 
226  !Cloudy layer properties
227  direct_s(1) = 1.
228  do l = 1, nlm
229  !Delta-M on
230  fact = asym_cld(i,l)*asym_cld(i,l)
231  asy = asym_cld(i,l)/(1. + asym_cld(i,l))
232  !Delta-M off
233  !fact = 0.
234  !asy = asym_cld(i,l)
235  !End Delta-M on/off
236  oms = ((1. - fact)*wc_cld(i,l))/(1. - fact*wc_cld(i,l))
237  taus = (1. - fact*wc_cld(i,l))*tau_cld(i,l)
238 
239 
240 
241  exptau_s_cld(l) = exp(-taus/amu0(i))
242  direct_s(l+1) = exptau_s_cld(l)*direct_s(l)
243  exptau_us_cld(l) = exp(-tau_cld(i,l)/amu0(i))
244  !direct_us_cld is not needed till later
245  ! Local eddington coefficients
246  t = 0.25*(7. - oms*(4. + 3.*asy))
247  r = -0.25*(1. - oms*(4. - 3.*asy))
248  kappa = sqrt(t*t - r*r)
249  rinf = r/(kappa+t)
250  ggtau = kappa*taus
251  eggtau = exp(-ggtau)
252  denom = (1. - rinf*rinf*eggtau*eggtau)
253  tr_cld(l) = (1. - rinf*rinf)*eggtau/denom
254  rr_cld(l) = rinf*(1. - eggtau*eggtau)/denom
255 
256  fact = kappa*kappa - 1./(amu0(i)*amu0(i))
257  if (abs(fact) .lt. eps) then
258  fact = 1./eps
259  else
260  fact = 1./fact
261  endif
262 
263  cc = oms*slr(i)*fact
264  g3 = 0.5 - 0.75*asy*amu0(i)
265  g4 = 1. - g3
266  aa = g3*(t - 1./amu0(i)) + g4*r
267  bb = g4*(t + 1./amu0(i)) + g3*r
268  sigu_cld_fract(l) = cc*((aa - rr_cld(l)*bb) - aa*tr_cld(l)*exptau_s_cld(l))
269  sigu(l) = sigu_cld_fract(l)*direct_s(l)
270  sigd_cld_fract(l) = cc*(-bb*tr_cld(l) + (bb-rr_cld(l)*aa)*exptau_s_cld(l))
271  sigd(l) = sigd_cld_fract(l)*direct_s(l)
272  enddo
273 
274  !BEGIN MAXIMALLY OVERLAPPED CALC
275  !For layers where 0. < cf_max < 1., treat as randomly combined
276  direct_s(1) = 1.
277  direct_us(1) = 1.
278  do l = 1,nlm
279  rr_tmp(l) = cf_max(i,l)*rr_cld(l) + (1. - cf_max(i,l))*rr_clr(l)
280  tr_tmp(l) = cf_max(i,l)*tr_cld(l) + (1. - cf_max(i,l))*tr_clr(l)
281  sigu_tmp_fract(l) = cf_max(i,l)*sigu_cld_fract(l) + (1. - cf_max(i,l))*sigu_clr_fract(l)
282  sigd_tmp_fract(l) = cf_max(i,l)*sigd_cld_fract(l) + (1. - cf_max(i,l))*sigd_clr_fract(l)
283  direct_s(l+1) = (cf_max(i,l)*exptau_s_cld(l) + (1. - cf_max(i,l))*exptau_s_clr(l))*direct_s(l)
284  direct_us(l+1) = (cf_max(i,l)*exptau_us_cld(l) + (1. - cf_max(i,l))*exptau_us_clr(l))*direct_us(l)
285  enddo
286  !MAXIMALLY OVERLAPPED adding
287  re(1) = 0.
288  vd(1) = 0.
289  do l = 1, nlm
290  prop = 1./(1. - re(l)*rr_tmp(l))
291  re(l+1) = rr_tmp(l) + tr_tmp(l)*tr_tmp(l)*re(l)*prop
292  vd(l+1) = sigd_tmp_fract(l)*direct_s(l) + (tr_tmp(l)*vd(l) + tr_tmp(l)*re(l)*sigu_tmp_fract(l)*direct_s(l))*prop
293  vu(l) = (rr_tmp(l)*vd(l) + sigu_tmp_fract(l)*direct_s(l))*prop
294  td(l) = prop
295  enddo
296 
297  !MAXIMALLY OVERLAPPED diffuse flux calculation
298  fddif_max(1) = 0.
299  fudif_max(nlm+1) = (asdif(i,ib)*vd(nlm+1) + &
300  asdir(i,ib)*slr(i)*amu0(i)*direct_s(nlm+1))/ &
301  (1. - asdif(i,ib)*re(nlm+1))
302  do l = nlm+1, 2, -1
303  fddif_max(l) = re(l)*fudif_max(l) + vd(l)
304  fudif_max(l-1) = tr_tmp(l-1)*fudif_max(l)*td(l-1) + vu(l-1)
305  enddo
306 
307  !MAXIMALLY OVERLAPPED direct flux calculation
308  fddir_max(1) = amu0(i)*slr(i)
309  do l = 1, nlm
310  fd_total = amu0(i)*slr(i)*direct_s(l+1) + fddif_max(l+1)
311  fddir_max(l+1) = amu0(i)*slr(i)*direct_us(l+1)
312  fddif_max(l+1) = fd_total - fddir_max(l+1)
313  enddo
314  !END MAXIMALLY OVERLAPPED CALC
315 
316  !BEGIN RANDOMLY-OVERLAPPED CALC
317  direct_s(1) = 1.
318  direct_us(1) = 1.
319  do l = 1, nlm
320  rr_tmp(l) = cf_random(i,l)*rr_cld(l) + (1. - cf_random(i,l))*rr_clr(l)
321  tr_tmp(l) = cf_random(i,l)*tr_cld(l) + (1. - cf_random(i,l))*tr_clr(l)
322  sigu_tmp_fract(l) = cf_random(i,l)*sigu_cld_fract(l) + (1. - cf_random(i,l))*sigu_clr_fract(l)
323  sigd_tmp_fract(l) = cf_random(i,l)*sigd_cld_fract(l) + (1. - cf_random(i,l))*sigd_clr_fract(l)
324  direct_s(l+1) = (cf_random(i,l)*exptau_s_cld(l) + (1. - cf_random(i,l))*exptau_s_clr(l))*direct_s(l)
325  direct_us(l+1) = (cf_random(i,l)*exptau_us_cld(l) + (1. - cf_random(i,l))*exptau_us_clr(l))*direct_us(l)
326  enddo
327  !RANDOMLY-OVERLAPPED adding
328  re(1) = 0.
329  vd(1) = 0.
330  do l = 1, nlm
331  prop = 1./(1. - re(l)*rr_tmp(l))
332  re(l+1) = rr_tmp(l) + tr_tmp(l)*tr_tmp(l)*re(l)*prop
333  vd(l+1) = sigd_tmp_fract(l)*direct_s(l) + (tr_tmp(l)*vd(l) + tr_tmp(l)*re(l)*sigu_tmp_fract(l)*direct_s(l))*prop
334  vu(l) = (rr_tmp(l)*vd(l) + sigu_tmp_fract(l)*direct_s(l))*prop
335  td(l) = prop
336  enddo
337 
338  !RANDOMLY-OVERLAPPED diffuse flux calculation
339  fddif_rndm(1) = 0.
340  fudif_rndm(nlm+1)= (asdif(i,ib)*vd(nlm+1) + &
341  asdir(i,ib)*slr(i)*amu0(i)*direct_s(nlm+1))/ &
342  (1. - asdif(i,ib)*re(nlm+1))
343  do l = nlm+1, 2, -1
344  fddif_rndm(l) = re(l)*fudif_rndm(l) + vd(l)
345  fudif_rndm(l-1) = tr_tmp(l-1)*fudif_rndm(l)*td(l-1) + vu(l-1)
346  enddo
347 
348  !RANDOMLY-OVERLAPPED direct flux calculation
349  fddir_rndm(1) = amu0(i)*slr(i)
350  do l = 1, nlm
351  fd_total = amu0(i)*slr(i)*direct_s(l+1) + fddif_rndm(l+1)
352  fddir_rndm(l+1) = amu0(i)*slr(i)*direct_us(l+1)
353  fddif_rndm(l+1) = fd_total - fddir_rndm(l+1)
354  enddo
355  !END RANDOMLY-OVERLAPPED CALC
356  !BEGIN ALL-SKY CALC
357  !This is the combined maximally- and randomly overlapped portions of the sky
358  fddif_all(i,1) = 0.
359  fddir_all(i,1) = amu0(i)*slr(i)
360  fudif_all(i,nlm+1) = (1. - c_maximal(i))*fudif_rndm(nlm+1) + c_maximal(i)*fudif_max(nlm+1)
361  do l = 1, nlm
362  fddif_all(i,l+1) = (1. - c_maximal(i))*fddif_rndm(l+1) + c_maximal(i)*fddif_max(l+1)
363  fddir_all(i,l+1) = (1. - c_maximal(i))*fddir_rndm(l+1) + c_maximal(i)*fddir_max(l+1)
364  fudif_all(i,l) = (1. - c_maximal(i))*fudif_rndm(l) + c_maximal(i)*fudif_max(l)
365  enddo
366  enddo ! Loop over columns
367  end subroutine
Generated on Tue Aug 13 2013 10:38:25 for BUGSrad by   doxygen 1.8.3.1