BUGSrad
 All Classes Files Functions Variables
two_rt_sw_bs.f
Go to the documentation of this file.
1 
2 
3 ! CVS: $Id: two_rt_sw_bs.F,v 1.3 2003/11/11 21:55:13 norm Exp $
4 ! CVS: $Name: $
5 
6 !-----------------------------------------------------------------------
7 
8  subroutine two_rt_sw_bs
9  + ( ncol , nlm , mbs , ib
10  +, slr , amu0 , wc , wcclr
11  +, asym , asyclr , tau , tauclr
12  +, asdir , asdif , fudif , fddir
13  +, fddif ,sel_rules , b1 , b2
14  +, b3 , b4
15  + )
16 
17 
18  use kinds
19  use bandsolve
20 
21 
22 
23  implicit none
24 
25 !-----------------------------------------------------------------------
26 ! REFERENCES:
27 ! two_rt_sw replaces two_rt and add written by G. Stephens. two_rt_sw
28 ! computes the spectral fluxes using a two-stream approximation method.
29 ! Philip Partain, Philip Gabriel, and Laura D. Fowler/graben (09-08-99).
30 
31 ! MODIFICATIONS:
32 ! * changed declarations to adapt the code from BUGS4 to BUGS5.
33 ! Laura D. Fowler/slikrock (02-01-00).
34 
35 ! SUBROUTINES CALLED:
36 ! none.
37 
38 ! FUNCTIONS CALLED:
39 ! none.
40 
41 ! INCLUDED COMMONS:
42 ! none.
43 
44 ! ARGUMENT LIST VARIABLES:
45 ! All arrays indexed as nlm correspond to variables defined in the
46 ! middle of layers. All arrays indexed as nlm+1 correspond to variables
47 ! defined at levels at the top and bottom of layers.
48 
49 ! INPUT ARGUMENTS:
50 ! ----------------
51  logical (kind=log_kind), intent(in)::
52  & sel_rules
53 
54  integer (kind=int_kind), intent(in)::
55  & ncol !Length of sub-domain.
56  &, nlm !Number of layers.
57  &, mbs !Number of SW spectral intervals.
58  &, ib !Index of spectral interval.
59 
60  real (kind=dbl_kind), intent(in), dimension(ncol)::
61  & slr !Fraction of daylight (-).
62  &, amu0 !Cosine of the solar zenith angle (-).
63 
64  real (kind=dbl_kind), intent(in), dimension(ncol,mbs)::
65  & asdir !Spectral direct surface albedo (-).
66  &, asdif !Spectral diffuse surface albedo (-).
67 
68  real (kind=dbl_kind), intent(in), dimension(ncol,nlm)::
69  & wc !Single scattering albedo (-).
70  &, wcclr !Single scattering albedo (-).
71  &, asym !Asymmetry factor (-).
72  &, asyclr!Asymmetry factor (-).
73  &, tau !Optical depth (-).
74  &, tauclr!Optical depth (-).
75  &, b1 !Cloud overlap parameter (-).
76  &, b2 !Cloud overlap parameter (-).
77  &, b3 !Cloud overlap parameter (-).
78  &, b4 !Cloud overlap parameter (-).
79 
80 ! OUTPUT ARGUMENTS:
81 ! -----------------
82  real (kind=dbl_kind), intent(out), dimension(ncol,nlm+1)::
83  & fddir !Spectral direct downward flux (W/m^2).
84  &, fddif !Spectral diffuse downward flux (W/m^2).
85  &, fudif !Spectral diffuse upward flux (W/m^2).
86 
87 ! LOCAL VARIABLES:
88 ! ----------------
89  integer(kind=int_kind)
90  & i ! Horizontal index.
91  &, l ! Vertical index.
92 
93  integer (kind=int_kind), dimension(nlm*4+2)::
94  & indx !Vector used by bandec and banbks
95 
96  real (kind=dbl_kind), dimension(nlm)::
97  & rrcld !All sky global reflection
98  &, rrclr !Clear sky global reflection
99  &, trcld !All sky global transmission
100  &, trclr !Clear sky global transmission
101  &, sigucld !All sky upwelling source
102  &, siguclr !Clear sky upwelling source
103  &, sigdcld !All sky downwelling source
104  &, sigdclr !Clear sky downwelling source
105 
106  real (kind=dbl_kind), dimension(nlm*4+2,11)::
107  & a !Diagonal matrix for bandec and banbks
108  &, al !Matrix used by bandec and banbks
109 
110  real (kind=dbl_kind), dimension(nlm*4+2)::
111  & b !Vector of sources (for A*x = b)
112 
113  real(kind=dbl_kind)
114  & exptaucld
115  &, exptauclr
116  &, directcld(nlm+1)
117  &, directclr(nlm+1)
118  &, d
119 
120  real(kind=dbl_kind)
121  & aa , bb , cc , denom
122  &, eggtau , eps , g3 , g4
123  &, ggtau , kappa , r , rinf
124  &, t , oms , taus , fact
125  &, asy
126  data eps /1.e-02/
127 
128  real (kind=dbl_kind), dimension(nlm+1)::
129  & re , vd
130 
131 ! SELECTION RULE VARIABLES:
132 ! -------------------------
133  logical (kind=log_kind)::
134  & fail
135 
136  real (kind=dbl_kind)::
137  & tausthresh
138  &, wcthresh
139  &, tauscat
140 
141 !#ifdef usenewexp
142 ! real(kind=dbl_kind), external :: exp
143 !#endif
144 
145 ! data tausthresh / 0.001 /
146 ! data wcthresh / 0.975 /
147  data tausthresh / 0.01 /
148  data wcthresh / 0.98 /
149 
150 !----------------------------------------------------------------------
151 
152  do 1000 i = 1, ncol
153 
154 ! if (sel_rules) then
155 ! fail = .false.
156 ! tauscat = 0.0
157 ! do l = nlm, 1, -1
158 ! if (wc(i,l).gt.wcthresh) fail = .true.
159 ! tauscat = tauscat + wc(i,l)*tau(i,l)
160 ! enddo
161 ! if (fail.and.tauscat.ge.tausthresh) goto 2000
162 !
163 !!>> BEGIN SELECTION RULES <<
164 !! print *,'selection rules'
165 ! fddir(i,1) = amu0(i)*slr(i)
166 ! fddif(i,:) = 0.0
167 ! do l=1,nlm
168 ! fddir(i,l+1) = fddir(i,l) * exp(-1.*tau(i,l)/amu0(i))
169 ! enddo
170 !
171 ! fudif(i,nlm+1) = fddir(i,nlm+1) * asdir(i,ib)
172 !
173 ! do l=nlm,1,-1
174 ! fudif(i,l) = fudif(i,l+1) * exp(-2*tau(i,l))
175 ! enddo
176 !
177 ! cycle
178 ! endif
179 !!>> END SELECTION RULES <<
180 
182 2000 directcld(1) = 0.
183  directclr(1) = 1.
184  re(1) = 0.
185  vd(1) = 0.
186 ! print *,'full up calculation'
187 
188 !---- 1. DO SHORTWAVE:
189  !ALL SKY
190  do l = 1, nlm
191  fact = asym(i,l)*asym(i,l)
192  oms = ((1.-fact)*wc(i,l))/(1.-fact*wc(i,l))
193  taus = (1.-fact*wc(i,l))*tau(i,l)
194  asy = asym(i,l)/(1.+asym(i,l))
195 
196  exptaucld = exp(-taus/amu0(i))
197 
198 !--- local coefficients: delta-eddington
199  t = 0.25 * (7. - oms*(4.+3.*asy))
200  r = -0.25 * (1. - oms*(4.-3.*asy))
201  kappa = sqrt(t**2-r**2)
202  rinf = r/(kappa+t)
203  ggtau = kappa*taus
204 
205  eggtau = exp(-ggtau)
206  denom = (1.-rinf**2*eggtau**2)
207  trcld(l) = (1.-rinf**2)*eggtau/denom
208  rrcld(l) = rinf*(1.-eggtau**2)/denom
209 
210  if(abs(kappa**2-1./amu0(i)**2) .lt. eps) then
211  fact = 1./eps
212  else
213  fact = 1./(kappa**2-1./amu0(i)**2)
214  endif
215  !cc = oms*slr(i)*fact
216  cc = oms*fact
217  g3 = 0.5-0.75*asy*amu0(i)
218  g4 = 1.-g3
219  aa = g3*(t-1./amu0(i))+g4*r
220  bb = g4*(t+1./amu0(i))+g3*r
221 
222  sigucld(l) = cc*((aa-rrcld(l)*bb)-aa*trcld(l)*exptaucld) *
223  & (b3(i,l)*directcld(l)+(1.-b1(i,l))*directclr(l))
224  sigdcld(l) = cc*(-bb*trcld(l)+(bb-rrcld(l)*aa)*exptaucld) *
225  & (b3(i,l)*directcld(l)+(1.-b1(i,l))*directclr(l))
226 
227  !CLEAR SKY
228  fact = asyclr(i,l)*asyclr(i,l)
229  oms = ((1.-fact)*wcclr(i,l))/(1.-fact*wcclr(i,l))
230  taus = (1.-fact*wcclr(i,l))*tauclr(i,l)
231  asy = asyclr(i,l)/(1.+asyclr(i,l))
232 
233  exptauclr = exp(-taus/amu0(i))
234 
235 
236 !--- local coefficients: delta-eddington
237  t = 0.25 * (7. - oms*(4.+3.*asy))
238  r = -0.25 * (1. - oms*(4.-3.*asy))
239  kappa = sqrt(t**2-r**2)
240  rinf = r/(kappa+t)
241  ggtau = kappa*taus
242 
243  eggtau = exp(-ggtau)
244 
245  denom = (1.-rinf**2*eggtau**2)
246  trclr(l) = (1.-rinf**2)*eggtau/denom
247  rrclr(l) = rinf*(1.-eggtau**2)/denom
248 
249  if(abs(kappa**2-1./amu0(i)**2) .lt. eps) then
250  fact = 1./eps
251  else
252  fact = 1./(kappa**2-1./amu0(i)**2)
253  endif
254  !cc = oms*slr(i)*fact
255  cc = oms*fact
256  g3 = 0.5-0.75*asy*amu0(i)
257  g4 = 1.-g3
258  aa = g3*(t-1./amu0(i))+g4*r
259  bb = g4*(t+1./amu0(i))+g3*r
260 
261  siguclr(l) = cc*((aa-rrclr(l)*bb)-aa*trclr(l)*exptauclr) *
262  & (b1(i,l)*directclr(l)+(1.-b3(i,l))*directcld(l))
263  sigdclr(l) = cc*(-bb*trclr(l)+(bb-rrclr(l)*aa)*exptauclr) *
264  & (b1(i,l)*directclr(l)+(1.-b3(i,l))*directcld(l))
265 
266  directclr(l+1) = exptauclr *
267  & ((1.-b3(i,l))*directcld(l) + b1(i,l)*directclr(l))
268  directcld(l+1) = exptaucld *
269  & (b3(i,l)*directcld(l) + (1.-b1(i,l))*directclr(l))
270  enddo
271 
272 !---- 2. LOAD A MATRIX, B MATRIX
273  a(:,:) = 0.0
274 
275  a(1,6) = -1.0
276  a(1,10) = trcld(1) * b4(i,1)
277  a(1,11) = trcld(1) * (1.-b2(i,1))
278 
279  a(2,6) = -1.0
280  a(2,9) = trclr(1) * (1.-b4(i,1))
281  a(2,10) = trclr(1) * b2(i,1)
282 
283  a(3,6) = -1.0
284  a(3,8) = rrcld(1) * b4(i,1)
285  a(3,9) = rrcld(1) * (1.-b2(i,1))
286 
287  a(4,6) = -1.0
288  a(4,7) = rrclr(1) * (1.-b4(i,1))
289  a(4,8) = rrclr(1) * b2(i,1)
290 
291  do l = 2,nlm
292  a(l*4-3,4) = rrcld(l) * b3(i,l)
293  a(l*4-3,5) = rrcld(l) * (1.-b1(i,l))
294  a(l*4-3,6) = -1.0
295  a(l*4-3,10) = trcld(l) * b4(i,l)
296  a(l*4-3,11) = trcld(l) * (1.-b2(i,l))
297 
298  a(l*4-2,3) = rrclr(l) * (1.-b3(i,l))
299  a(l*4-2,4) = rrclr(l) * b1(i,l)
300  a(l*4-2,6) = -1.0
301  a(l*4-2,9) = trclr(l) * (1.-b4(i,l))
302  a(l*4-2,10) = trclr(l) * b2(i,l)
303 
304  a(l*4-1,2) = trcld(l) * b3(i,l)
305  a(l*4-1,3) = trcld(l) * (1.-b1(i,l))
306  a(l*4-1,6) = -1.0
307  a(l*4-1,8) = rrcld(l) * b4(i,l)
308  a(l*4-1,9) = rrcld(l) * (1.-b2(i,l))
309 
310  a(l*4,1) = trclr(l) * (1.-b3(i,l))
311  a(l*4,2) = trclr(l) * b1(i,l)
312  a(l*4,6) = -1.0
313  a(l*4,7) = rrclr(l) * (1.-b4(i,l))
314  a(l*4,8) = rrclr(l) * b2(i,l)
315  enddo
316 
317  a(nlm*4+1,4) = asdif(i,ib)
318  a(nlm*4+1,6) = -1.0
319  a(nlm*4+2,4) = asdif(i,ib)
320  a(nlm*4+2,6) = -1.0
321 
322  b(:) = 0.0
323  do l = 1,nlm
324 ! b(l*4-3) = -1.*(b3(i,l)*sigucld(l) + (1.-b1(i,l)) *siguclr(l))
325 ! b(l*4-2) = -1.*((1.-b3(i,l))*sigucld(l) + b1(i,l) *siguclr(l))
326 ! b(l*4-1) = -1.*(b3(i,l)*sigdcld(l) + (1.-b1(i,l)) *sigdclr(l))
327 ! b(l*4) = -1.*((1.-b3(i,l))*sigdcld(l) + b1(i,l) *sigdclr(l))
328  b(l*4-3) = -sigucld(l)
329  b(l*4-2) = -siguclr(l)
330  b(l*4-1) = -sigdcld(l)
331  b(l*4) = -sigdclr(l)
332  enddo
333  b(nlm*4+1) = -asdir(i,ib)*directcld(nlm+1)*amu0(i)
334  b(nlm*4+2) = -asdir(i,ib)*directclr(nlm+1)*amu0(i)
335 
336 ! do l = 1,nlm*4+2
337 ! print '(15E15.7)',(a(l,k),k=1,11)
338 ! enddo
339 
340 ! do l = 1,nlm*4+2
341 ! print *,l,b(l)
342 ! enddo
343 !
344  call bandec(a,nlm*4+2,5,5,nlm*4+2,11,al,11,indx,d)
345  call banbks(a,nlm*4+2,5,5,nlm*4+2,11,al,11,indx,b)
346 
347 !---- 3. SUM CLEAR AND CLOUDY FLUXES
348  do l = 1,nlm+1
349  fudif(i,l) = slr(i)*(b(l*4-3)+b(l*4-2))
350  enddo
351  do l = 1,nlm
352  fddif(i,l+1) = slr(i)*(b(l*4-1)+b(l*4))
353  enddo
354 
355  fddir(i,1) = amu0(i)*slr(i)
356  do l = 1, nlm
357  fddir(i,l+1) = amu0(i)*slr(i)*
358  & (directcld(l+1)+directclr(l+1))
359  enddo
360 ! do l = 1, nlm+1
361 ! print *,ib,l,fddir(l),fddif(l),fudif(l)
362 ! enddo
364 
365  1000 continue
366 
367  return
368  end subroutine two_rt_sw_bs
Generated on Tue Aug 13 2013 10:38:25 for BUGSrad by   doxygen 1.8.3.1