BUGSrad
 All Classes Files Functions Variables
two_rt_lw_iter.f
Go to the documentation of this file.
1 
2 
3 ! CVS: $Id: two_rt_lw_iter.F,v 1.2 2003/11/11 21:55:13 norm Exp $
4 ! CVS: $Name: $
5 
6 !-----------------------------------------------------------------------
7 
8  subroutine two_rt_lw_iter
9  + (
10  + ncol , nlm , mbs , mbir
11  +, ib , cldamt , wc , wcclr
12  +, asym , asyclr , tau , tauclr
13  +, es , bf , fu , fd
14  +, sel_rules , b1 , b2 , b3
15  +, b4
16  + )
17 
18  use kinds
19 
20 
21 
22  implicit none
23 
24 !-----------------------------------------------------------------------
25 ! REFERENCES:
26 ! two_rt_lw replaces two_rt and add written by G. Stephens. two_rt_lw
27 ! computes the spectral fluxes using a two-stream approximation method.
28 ! Philip Partain, Philip Gabriel, and Laura D. Fowler/graben (09-08-99).
29 
30 ! MODIFICATIONS:
31 ! * changed declarations to adapt the code from BUGS4 to BUGS5.
32 ! Laura D. Fowler/slikrock (02-01-00).
33 
34 ! SUBROUTINES CALLED:
35 ! none.
36 
37 ! FUNCTIONS CALLED:
38 ! none.
39 
40 ! INCLUDED COMMONS:
41 ! none.
42 
43 ! ARGUMENT LIST VARIABLES:
44 ! All arrays indexed as nlm correspond to variables defined in the
45 ! middle of layers. All arrays indexed as nlm+1 correspond to variables
46 ! defined at levels at the top and bottom of layers.
47 
48 ! INPUT ARGUMENTS:
49 ! ----------------
50  logical (kind=log_kind), intent(in)::
51  & sel_rules
52 
53  integer (kind=int_kind), intent(in)::
54  & ncol !Length of sub-domain.
55  &, nlm !Number of layers.
56  &, mbs !Number of SW spectral intervals.
57  &, mbir !Number of IR spectral intervals.
58  &, ib !Index of spectral interval.
59 
60  real (kind=dbl_kind), intent(in), dimension(ncol,mbir)::
61  & es !Spectral surface emissivity (-).
62 
63  real (kind=dbl_kind), intent(in), dimension(ncol,nlm)::
64  & cldamt !Cloud fraction (-).
65  &, wc !All sky single scattering albedo (-).
66  &, wcclr !Clear sky single scattering albedo (-).
67  &, asym !All sky asymmetry factor (-).
68  &, asyclr !Clear sky asymmetry factor (-).
69  &, tau !All sky optical depth (-).
70  &, tauclr !Clear sky optical depth (-).
71  &, b1 !Cloud overlap parameter (-).
72  &, b2 !Cloud overlap parameter (-).
73  &, b3 !Cloud overlap parameter (-).
74  &, b4 !Cloud overlap parameter (-).
75 
76  real (kind=dbl_kind), intent(in), dimension(ncol,nlm+1)::
77  & bf !Planck function (W/m^2).
78 
79 ! OUTPUT ARGUMENTS:
80 ! -----------------
81  real (kind=dbl_kind), intent(out), dimension(ncol,nlm+1)::
82  & fd !Spectral downward flux (W/m^2).
83  &, fu !Spectral upward flux (W/m^2).
84 
85 ! LOCAL VARIABLES:
86 
87  integer(kind=int_kind)
88  & i !Horizontal index.
89  &, l !Vertical index.
90  &, ibms !Index of spectral interval.
91  &, j
92  &, nsr
93  &, nsmx
94  &, n
95  &, ii
96  &, jj
97  &, kk
98  &, ir
99  &, iter
100  integer (kind=int_kind), dimension(16*nlm-6)::
101  & idc
102 
103  integer (kind=int_kind), dimension(4*nlm+2)::
104  & nir
105 
106  real (kind=dbl_kind), dimension(nlm)::
107  & rrcld !All sky global reflection (-).
108  &, rrclr !Clear sky global reflection (-).
109  &, trcld !All sky global transmission (-).
110  &, trclr !Clear sky global transmission (-).
111  &, sigucld !All sky upwelling source (-).
112  &, siguclr !Clear sky upwelling source (-).
113  &, sigdcld !All sky downwelling source (-).
114  &, sigdclr !Clear sky downwelling source (-).
115  &, exptau !
116 
117  real (kind=dbl_kind), dimension(4*nlm+2)::
118  & b
119  &, fvc
120  &, error
121 ! &, old_fvc
122 
123  real (kind=dbl_kind), dimension(16*nlm-6)::
124  & smx
125 
126  real(kind=dbl_kind)
127  & aa , bb , beta0 , cc
128  &, diffac , denom , fact , eggtau
129  &, ggtau
130  &, kappa , oms , prop ,r, rinf
131  &, t , taus , omega
132  data diffac /2./
133 
134 ! SELECTION RULE VARIABLES
135 
136  logical (kind=log_kind)::
137  & fail
138 
139  real (kind=dbl_kind)::
140  & tausthresh
141  &, wcthresh
142  &, tauscat
143  data tausthresh / 0.001 /
144  data wcthresh / 0.975 /
145 
146 !----------------------------------------------------------------------
147  nsr = 4*nlm + 2
148  nsmx = 16*nlm - 6
149 
150  fd(:,1) = 0.
151 
152  ibms = ib - mbs
153 
154  do i = 1, ncol
155 
156 !---- 2. DO LONGWAVE:
157  do l = 1, nlm
158  !ALL SKY
159  fact = asym(i,l)*asym(i,l)
160  oms = ((1.-fact)*wc(i,l))/(1.-fact*wc(i,l))
161  taus = (1.-fact*wc(i,l))*tau(i,l)
162 
163  beta0 = (4.+asym(i,l))/(8.*(1.+asym(i,l)))
164  t = diffac*(1.-oms*(1.-beta0)) !-0.25
165  r = diffac*oms*beta0 !-0.25
166  kappa = sqrt(t**2-r**2)
167  rinf = r/(kappa+t)
168  ggtau = kappa*taus
169 
170  eggtau = exp(-ggtau)
171 
172  denom = (1.-rinf**2*eggtau**2)
173  trcld(l) = (1.-rinf**2)*eggtau/denom
174  rrcld(l) = rinf*(1.-eggtau**2)/denom
175 
176  if(taus .lt. .8e-2) then
177  sigucld(l) = cldamt(i,l)*0.5*diffac*(bf(i,l)+
178  * bf(i,l+1))*taus
179  sigdcld(l) = cldamt(i,l)*sigucld(l)
180  else
181  aa = (t+r)*(1.-rrcld(l))-(1.+rrcld(l)-trcld(l))/taus
182  bb = -(t+r)*trcld(l)+(1.+rrcld(l)-trcld(l))/taus
183  cc = diffac*(1.-oms)/kappa**2
184  sigucld(l) = cldamt(i,l)*cc*(aa*bf(i,l)+bb*bf(i,l+1))
185  sigdcld(l) = cldamt(i,l)*cc*(bb*bf(i,l)+aa*bf(i,l+1))
186  endif
187 
188  !CLEAR SKY
189  fact = asyclr(i,l)*asyclr(i,l)
190  oms = ((1.-fact)*wcclr(i,l))/(1.-fact*wcclr(i,l))
191  taus = (1.-fact*wcclr(i,l))*tauclr(i,l)
192 
193  beta0 = (4.+asyclr(i,l))/(8.*(1.+asyclr(i,l)))
194  t = diffac*(1.-oms*(1.-beta0)) !-0.25
195  r = diffac*oms*beta0 !-0.25
196  kappa = sqrt(t**2-r**2)
197  rinf = r/(kappa+t)
198  ggtau = kappa*taus
199 
200  eggtau = exp(-ggtau)
201 
202  denom = (1.-rinf**2*eggtau**2)
203  trclr(l) = (1.-rinf**2)*eggtau/denom
204  rrclr(l) = rinf*(1.-eggtau**2)/denom
205 
206  if(taus .lt. .8e-2) then
207  siguclr(l) = (1.0-cldamt(i,l))*0.5*diffac*(bf(i,l)+
208  * bf(i,l+1))*taus
209  sigdclr(l) = (1.0-cldamt(i,l))*siguclr(l)
210  else
211  aa = (t+r)*(1.-rrclr(l))-(1.+rrclr(l)-trclr(l))/taus
212  bb = -(t+r)*trclr(l)+(1.+rrclr(l)-trclr(l))/taus
213  cc = diffac*(1.-oms)/kappa**2
214  siguclr(l) = (1.0-cldamt(i,l))*cc*(aa*bf(i,l)+
215  * bb*bf(i,l+1))
216  sigdclr(l) = (1.0-cldamt(i,l))*cc*(bb*bf(i,l)+
217  * aa*bf(i,l+1))
218  endif
219 
220  enddo
221 
222 
223 !---- 1. LOAD SMX VECTOR
224  nir(:) = 4
225 
226  idc(1) = 5
227  idc(2) = 6
228  smx(1) = -trcld(1) * b4(i,1)
229  smx(2) = -trcld(1) * (1.-b2(i,1))
230  nir(1) = 2
231 
232  idc(3) = 5
233  idc(4) = 6
234  smx(3) = -trclr(1) * (1.-b4(i,1))
235  smx(4) = -trclr(1) * b2(i,1)
236  nir(2) = 2
237 
238  idc(5) = 5
239  idc(6) = 6
240  smx(5) = -rrcld(1) * b4(i,1)
241  smx(6) = -rrcld(1) * (1.-b2(i,1))
242  nir(3) = 2
243 
244  idc(7) = 5
245  idc(8) = 6
246  smx(7) = -rrclr(1) * (1.-b4(i,1))
247  smx(8) = -rrclr(1) * b2(i,1)
248  nir(4) = 2
249 
250  do l = 1,nlm-1
251  n = (l-1)*16 + 9
252  ir = 4*l
253 
254  idc(n) = ir-1
255  idc(n+1) = ir
256  idc(n+2) = ir+5
257  idc(n+3) = ir+6
258  smx(n) = -rrcld(l+1) * b3(i,l+1)
259  smx(n+1) = -rrcld(l+1) * (1.-b1(i,l+1))
260  smx(n+2) = -trcld(l+1) * b4(i,l+1)
261  smx(n+3) = -trcld(l+1) * (1.-b2(i,l+1))
262 
263  idc(n+4) = ir-1
264  idc(n+5) = ir
265  idc(n+6) = ir+5
266  idc(n+7) = ir+6
267  smx(n+4) = -rrclr(l+1) * (1.-b3(i,l+1))
268  smx(n+5) = -rrclr(l+1) * b1(i,l+1)
269  smx(n+6) = -trclr(l+1) * (1.-b4(i,l+1))
270  smx(n+7) = -trclr(l+1) * b2(i,l+1)
271 
272  idc(n+8) = ir-1
273  idc(n+9) = ir
274  idc(n+10) = ir+5
275  idc(n+11) = ir+6
276  smx(n+8) = -trcld(l+1) * b3(i,l+1)
277  smx(n+9) = -trcld(l+1) * (1.-b1(i,l+1))
278  smx(n+10) = -rrcld(l+1) * b4(i,l+1)
279  smx(n+11) = -rrcld(l+1) * (1.-b2(i,l+1))
280 
281  idc(n+12) = ir-1
282  idc(n+13) = ir
283  idc(n+14) = ir+5
284  idc(n+15) = ir+6
285  smx(n+12) = -trclr(l+1) * (1.-b3(i,l+1))
286  smx(n+13) = -trclr(l+1) * b1(i,l+1)
287  smx(n+14) = -rrclr(l+1) * (1.-b4(i,l+1))
288  smx(n+15) = -rrclr(l+1) * b2(i,l+1)
289  enddo
290 
291  ir = 4*nlm
292 
293  idc(16*nlm-7) = 4*nlm-1
294  idc(16*nlm-6) = 4*nlm
295  smx(16*nlm-7) = 0.0 !1.-es(i,ibms)
296  smx(16*nlm-6) = 0.0 !1.-es(i,ibms)
297  nir(ir+1) = 1
298  nir(ir+2) = 1
299 
300  b(:) = 0.0
301  do l = 1,nlm
302  b(l*4-3) = sigucld(l)
303  b(l*4-2) = siguclr(l)
304  b(l*4-1) = sigdcld(l)
305  b(l*4) = sigdclr(l)
306  enddo
307  b(nlm*4+1) = cldamt(i,nlm)*es(i,ibms)*bf(i,nlm+1)
308  b(nlm*4+2) = (1.-cldamt(i,nlm))*es(i,ibms)*bf(i,nlm+1)
309 
310 
311 
312 
313 !-------------- GAUSS SEIDEL W/ OVERRELAXATION ------------------
314  omega = 1.0
315 
316  fvc(:) = b(:)
317 
318  do iter=1,200
319  kk = 1
320  do ii=1,nsr
321  t = 0.0
322  do j=1,nir(ii)
323  jj = idc(kk)
324  t = t + smx(kk) * fvc(jj)
325  kk = kk + 1
326  enddo
327  t = t + fvc(ii)
328  fvc(ii) = fvc(ii) + omega * (b(ii)-t)
329  error(ii) = b(ii) - t
330  enddo
331 
332  if (maxval(abs(error)) .le. 0.05) then
333  !print *,omega,iter,' iterations'
334  exit
335  endif
336  enddo
337 
338 
339 !---- 3. SUM CLEAR AND CLOUDY FLUXES
340  do l = 1,nlm+1
341  fu(i,l) = fvc(l*4-3)+fvc(l*4-2)
342  enddo
343  do l = 1,nlm
344  fd(i,l+1) = fvc(l*4-1)+fvc(l*4)
345  enddo
346 
347  end do !i=1,ncol
348 
349  return
350  end
Generated on Tue Aug 13 2013 10:38:25 for BUGSrad by   doxygen 1.8.3.1