BUGSrad
 All Classes Files Functions Variables
two_rt_sw.f
Go to the documentation of this file.
1 
2 
3 ! CVS: $Id: two_rt_sw.F,v 1.7 2006/11/16 21:19:38 norm Exp $
4 ! CVS: $Name: $
5 
6 !NBW This has been modified to apply delta-M descaling
7 !NBW to the direct and diffuse downwelling fluxes
8 
9 
10 !-----------------------------------------------------------------------
11 
12  subroutine two_rt_sw
13  + ( ncol , nlm , mbs , ib
14  +, slr , amu0 , wc , asym
15  +, tau , asdir , asdif , fudif
16  +, fddir , fddif ,sel_rules
17  + )
18 
19  use kinds
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  &, asym !Asymmetry factor (-).
71  &, tau !Optical depth (-).
72 
73 ! OUTPUT ARGUMENTS:
74 ! -----------------
75  real (kind=dbl_kind), intent(out), dimension(ncol,nlm+1)::
76  & fddir !Spectral direct downward flux (W/m^2).
77  &, fddif !Spectral diffuse downward flux (W/m^2).
78  &, fudif !Spectral diffuse upward flux (W/m^2).
79 
80 ! LOCAL VARIABLES:
81 
82  integer(kind=int_kind)
83  & i ! Horizontal index.
84  &, l ! Vertical index.
85 
86  real (kind=dbl_kind), dimension(nlm)::
87  & rr !
88  &, tr !
89  &, sigu !
90  &, sigd !
91 
92  real(kind=dbl_kind)
93  & exptau_s(nlm),direct_s(nlm+1),exptau_us(nlm),direct_us(nlm+1)
94 
95  real(kind=dbl_kind)
96  & aa , bb , cc , denom
97  &, eggtau , eps , g3 , g4
98  &, ggtau , kappa , r , rinf
99  &, t , oms , taus , fact
100  &, asy , prop, fd_total
101  data eps /1.e-02/
102 
103  real (kind=dbl_kind), dimension(nlm)::
104  & td , vu
105 
106  real (kind=dbl_kind), dimension(nlm+1)::
107  & re , vd
108 
109 ! SELECTION RULE VARIABLES
110 
111  logical (kind=log_kind)::
112  & fail
113 
114  real (kind=dbl_kind)::
115  & tausthresh
116  &, wcthresh
117  &, tauscat
118 
119 
120 ! data tausthresh / 0.001 /
121 ! data wcthresh / 0.975 /
122  data tausthresh / 0.01 /
123  data wcthresh / 0.98 /
124 
125 !----------------------------------------------------------------------
126 
127  do 1000 i = 1, ncol
128 
129  if(sel_rules) then
130 
131  fail = .false.
132  tauscat = 0.0
133  do l = nlm, 1, -1
134  if (wc(i,l).gt.wcthresh) fail = .true.
135  tauscat = tauscat + wc(i,l)*tau(i,l)
136  enddo
137 
138  if (fail.and.tauscat.ge.tausthresh) goto 2000
139 
141  print *,'selection rules'
142  fddir(i,1) = amu0(i)*slr(i)
143  fddif(i,:) = 0.0
144  do l = 1, nlm
145  fddir(i,l+1) = fddir(i,l) * exp(-1.*tau(i,l)/amu0(i))
146  enddo
147  fudif(i,nlm+1) = fddir(i,nlm+1) * asdir(i,ib)
148 
149  do l = nlm, 1, -1
150  fudif(i,l) = fudif(i,l+1) * exp(-2*tau(i,l))
151  enddo
152 
153  cycle
154 
155  endif
157 
159 2000 direct_s(1) = 1.
160  direct_us(1) = 1.
161  re(1) = 0.
162  vd(1) = 0.
163 
164  do l = 1, nlm
165  !delta-M on
166  fact = asym(i,l)*asym(i,l)
167  asy = asym(i,l)/(1.+asym(i,l)) !Orig w/ delta M
168  !This is the delta-M scaling approx
169  !delta-M off
170  !fact = 0.
171  !asy = asym(i,l)
172  !End delta-M on/off
173  oms = ((1.-fact)*wc(i,l))/(1.-fact*wc(i,l))
174  taus = (1.-fact*wc(i,l))*tau(i,l)
175 
176 
177  exptau_s(l) = exp(-taus/amu0(i)) !delta-M scaled
178  direct_s(l+1) = exptau_s(l)*direct_s(l) !delta-M scaled
179  exptau_us(l) = exp(-tau(i,l)/amu0(i)) !unscaled
180  direct_us(l+1) = exptau_us(l)*direct_us(l) !unscaled
181 
182 !---- local coefficients: delta-eddington
183  t = 0.25 * (7. - oms*(4.+3.*asy))
184  r = -0.25 * (1. - oms*(4.-3.*asy))
185  kappa = sqrt(t**2-r**2)
186  rinf = r/(kappa+t)
187  ggtau = kappa*taus
188  eggtau = exp(-ggtau)
189  denom = (1.-rinf**2*eggtau**2)
190  tr(l) = (1.-rinf**2)*eggtau/denom
191  rr(l) = rinf*(1.-eggtau**2)/denom
192 
193  if(abs(kappa**2-1./amu0(i)**2) .lt. eps) then
194  fact = 1./eps
195  else
196  fact = 1./(kappa**2-1./amu0(i)**2)
197  endif
198  cc = oms*slr(i)*fact
199  g3 = 0.5-0.75*asy*amu0(i)
200  g4 = 1.-g3
201  aa = g3*(t-1./amu0(i))+g4*r
202  bb = g4*(t+1./amu0(i))+g3*r
203  sigu(l) = cc*((aa-rr(l)*bb)-aa*tr(l)*exptau_s(l))
204  + * direct_s(l)
205  sigd(l) = cc*(-bb*tr(l)+(bb-rr(l)*aa)*exptau_s(l))
206  + * direct_s(l)
207  enddo
208 
209 !---- 1. do adding, going from top down:
210 
211  do l = 1, nlm
212  prop = 1. / (1. - re(l)*rr(l))
213  re(l+1) = rr(l) + tr(l)**2*re(l)*prop
214  vd(l+1) = sigd(l) + (tr(l)*vd(l)
215  + + tr(l)*re(l)*sigu(l))*prop
216  vu(l) = (rr(l)*vd(l) + sigu(l))*prop
217  td(l) = prop
218  enddo
219 
220 !---- 2. calculate diffuse fluxes:
221 
222  fddif(i,1) = 0.
223  fudif(i,nlm+1) = (asdif(i,ib)*vd(nlm+1)
224  + + asdir(i,ib)*slr(i)*amu0(i)*direct_s(nlm+1))
225  + / (1.-asdif(i,ib)*re(nlm+1))
226 
227  do l = nlm+1, 2, -1
228  fddif(i,l) = re(l)*fudif(i,l) + vd(l)
229  fudif(i,l-1) = tr(l-1)*fudif(i,l)*td(l-1) + vu(l-1)
230  enddo
231 
232 !---- 3. Compute direct flux and descale the direct and diffuse fluxes
233  fddir(i,1) = amu0(i)*slr(i)
234  do l = 1, nlm
235  fd_total = amu0(i)*slr(i)*direct_s(l+1) + fddif(i,l+1)
236  fddir(i,l+1) = amu0(i)*slr(i)*direct_us(l+1)
237  fddif(i,l+1) = fd_total - fddir(i,l+1)
238  enddo
239 
241 
242  1000 continue
243  return
244  end subroutine two_rt_sw
245 
246 !------------------------------------------------------------------------
Generated on Tue Aug 13 2013 10:38:25 for BUGSrad by   doxygen 1.8.3.1