chem/rtrans.F

   1 !=============================================================================
   2 ! This file contains the following subroutines, related to the solution of
   3 ! the equation of radiative transfer in multiple homogeneous layers.
   4 !     rtlink
   5 !     ps2str
   6 !        tridag
   7 !=============================================================================
   8
   9       MODULE rad_trans
  10
  11       IMPLICIT NONE
  12
  13       private
  14       public :: rtlink
  15
  16       CONTAINS
  17
  18       SUBROUTINE rtlink( &
  19            nstr, nlev, nlyr, nwave, &
  20            iw, albedo, zen, &
  21            dsdh, nid, &
  22            dtrl,  &
  23            dto3,  &
  24            dto2, &
  25            dtso2, &
  26            dtno2,  &
  27            dtcld, omcld, gcld, &
  28            dtaer, omaer, gaer, &
  29            dtsnw, omsnw, gsnw, &
  30 #ifdef SW_DEBUG
  31            edir, edn, eup, fdir, fdn, fup, tuv_diags )
  32 #else
  33            edir, edn, eup, fdir, fdn, fup )
  34 #endif
  35
  36       use params_mod, only : largest, pi
  37
  38 !---------------------------------------------------------------------
  39 !     ... dummy arguments
  40 !---------------------------------------------------------------------
  41       INTEGER, intent(in) :: nstr
  42       INTEGER, intent(in) :: nlev, nlyr
  43       INTEGER, intent(in) :: nwave, iw
  44       REAL, intent(in)    :: albedo
  45       REAL, intent(in)    :: zen
  46       REAL, intent(in)    :: dsdh(0:nlyr,nlyr)
  47       INTEGER, intent(in) :: nid(0:nlyr)
  48
  49       REAL, intent(in)    :: dtrl(nlyr,nwave)
  50       REAL, intent(in)    :: dto3(nlyr,nwave), dto2(nlyr,nwave)
  51       REAL, intent(in)    :: dtso2(nlyr,nwave), dtno2(nlyr,nwave)
  52       REAL, intent(in)    :: dtcld(nlyr,nwave), omcld(nlyr,nwave), gcld(nlyr,nwave)
  53       REAL, intent(in)    :: dtaer(nlyr,nwave), omaer(nlyr,nwave), gaer(nlyr,nwave)
  54       REAL, intent(in)    :: dtsnw(nlyr,nwave), omsnw(nlyr,nwave), gsnw(nlyr,nwave)
  55
  56       REAL, intent(out)   :: edir(nlev), edn(nlev), eup(nlev)
  57       REAL, intent(out)   :: fdir(nlev), fdn(nlev), fup(nlev)
  58
  59 #ifdef SW_DEBUG
  60       LOGICAL, intent(in) :: tuv_diags
  61 #endif
  62
  63
  64 !---------------------------------------------------------------------
  65 !     ... local variables
  66 !---------------------------------------------------------------------
  67       REAL, PARAMETER :: dr = pi/180.
  68
  69       INTEGER :: k, kk
  70       REAL    :: dtabs,dtsct,dscld,dsaer,dssnw,dacld,daaer,dasnw
  71       REAL    :: dt(nlyr), om(nlyr), g(nlyr)
  72
  73
  74 !---------------------------------------------------------------------
  75 !     ... specific to ps2str
  76 !---------------------------------------------------------------------
  77       LOGICAL, parameter :: delta = .true.
  78       REAL ediri(nlev), edni(nlev), eupi(nlev)
  79       REAL fdiri(nlev), fdni(nlev), fupi(nlev)
  80
  81 !---------------------------------------------------------------------
  82 ! initialize:
  83 !---------------------------------------------------------------------
  84       fdir(1:nlev) = 0.
  85       fup(1:nlev)  = 0.
  86       fdn(1:nlev)  = 0.
  87       edir(1:nlev) = 0.
  88       eup(1:nlev)  = 0.
  89       edn(1:nlev)  = 0.
  90
  91       DO k = 1, nlyr
  92         dscld = dtcld(k,iw)*omcld(k,iw)
  93         dacld = dtcld(k,iw)*(1.-omcld(k,iw))
  94
  95         dsaer = dtaer(k,iw)*omaer(k,iw)
  96         daaer = dtaer(k,iw)*(1.-omaer(k,iw))
  97
  98         dssnw = dtsnw(k,iw)*omsnw(k,iw)
  99         dasnw = dtsnw(k,iw)*(1.-omsnw(k,iw))
 100
 101         dtsct = dtrl(k,iw) + dscld + dsaer + dssnw
 102         dtabs = dtso2(k,iw) + dto2(k,iw) + dto3(k,iw) &
 103               + dtno2(k,iw) + dacld + daaer + dasnw
 104
 105         dtabs = max( dtabs,1./largest )
 106         dtsct = max( dtsct,1./largest )
 107
 108 !---------------------------------------------------------------------
 109 ! from bottom-up to top-down
 110 !---------------------------------------------------------------------
 111         kk = nlyr - k + 1
 112         dt(kk) = dtsct + dtabs
 113         g(kk)  = (gcld(k,iw)*dscld + gsnw(k,iw)*dssnw + gaer(k,iw)*dsaer)/dtsct
 114         IF( dtsct /= 1./largest ) then
 115           om(kk) = dtsct/(dtsct + dtabs)
 116         ELSE
 117           om(kk) = 1./largest
 118         ENDIF
 119       END DO
 120
 121 #ifdef SW_DEBUG
 122       if( tuv_diags .and. iw == 100 ) then
 123         open(unit=33,file='WRF-TUV.dbg1.out')
 124         do k = 1,nlyr,5
 125           write(33,'(1p,5g15.7)') dt(k:min(k+4,nlyr))
 126         end do
 127         do k = 1,nlyr,5
 128           write(33,'(1p,5g15.7)') g(k:min(k+4,nlyr))
 129         end do
 130         do k = 1,nlyr,5
 131           write(33,'(1p,5g15.7)') om(k:min(k+4,nlyr))
 132         end do
 133         close(33)
 134       endif
 135 #endif
 136
 137       CALL ps2str( nlyr, nlev, zen, albedo, &
 138                    dt, om, g, &
 139                    dsdh, nid, delta, &
 140                    fdiri, fupi, fdni, ediri, eupi, edni)
 141
 142 !---------------------------------------------------------------------
 143 ! from top-down to bottom-up
 144 !---------------------------------------------------------------------
 145       fdir(1:nlev) = fdiri(nlev:1:-1)
 146       fup(1:nlev)  = fupi(nlev:1:-1)
 147       fdn(1:nlev)  = fdni(nlev:1:-1)
 148       edir(1:nlev) = ediri(nlev:1:-1)
 149       eup(1:nlev)  = eupi(nlev:1:-1)
 150       edn(1:nlev)  = edni(nlev:1:-1)
 151
 152       END SUBROUTINE rtlink
 153
 154       SUBROUTINE ps2str( &
 155            nlyr, nlev, zen, rsfc, &
 156            tauu, omu, gu, &
 157            dsdh, nid, delta, &
 158            fdr, fup, fdn, edr, eup, edn)
 159 !-----------------------------------------------------------------------------*
 160 !=  PURPOSE:                                                                 =*
 161 !=  Solve two-stream equations for multiple layers.  The subroutine is based =*
 162 !=  on equations from:  Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.=*
 163 !=  It contains 9 two-stream methods to choose from.  A pseudo-spherical     =*
 164 !=  correction has also been added.                                          =*
 165 !-----------------------------------------------------------------------------*
 166 !=  PARAMETERS:                                                              =*
 167 !=  NLEVEL  - INTEGER, number of specified altitude levels in the working (I)=*
 168 !=            grid                                                           =*
 169 !=  ZEN     - REAL, solar zenith angle (degrees)                          (I)=*
 170 !=  RSFC    - REAL, surface albedo at current wavelength                  (I)=*
 171 !=  TAUU    - REAL, unscaled optical depth of each layer                  (I)=*
 172 !=  OMU     - REAL, unscaled single scattering albedo of each layer       (I)=*
 173 !=  GU      - REAL, unscaled asymmetry parameter of each layer            (I)=*
 174 !=  DSDH    - REAL, slant path of direct beam through each layer crossed  (I)=*
 175 !=            when travelling from the top of the atmosphere to layer i;     =*
 176 !=            DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1                            =*
 177 !=  NID     - INTEGER, number of layers crossed by the direct beam when   (I)=*
 178 !=            travelling from the top of the atmosphere to layer i;          =*
 179 !=            NID(i), i = 0..NZ-1                                            =*
 180 !=  DELTA   - LOGICAL, switch to use delta-scaling                        (I)=*
 181 !=            .TRUE. -> apply delta-scaling                                  =*
 182 !=            .FALSE.-> do not apply delta-scaling                           =*
 183 !=  FDR     - REAL, contribution of the direct component to the total     (O)=*
 184 !=            actinic flux at each altitude level                            =*
 185 !=  FUP     - REAL, contribution of the diffuse upwelling component to    (O)=*
 186 !=            the total actinic flux at each altitude level                  =*
 187 !=  FDN     - REAL, contribution of the diffuse downwelling component to  (O)=*
 188 !=            the total actinic flux at each altitude level                  =*
 189 !=  EDR     - REAL, contribution of the direct component to the total     (O)=*
 190 !=            spectral irradiance at each altitude level                     =*
 191 !=  EUP     - REAL, contribution of the diffuse upwelling component to    (O)=*
 192 !=            the total spectral irradiance at each altitude level           =*
 193 !=  EDN     - REAL, contribution of the diffuse downwelling component to  (O)=*
 194 !=            the total spectral irradiance at each altitude level           =*
 195 !-----------------------------------------------------------------------------*
 196
 197       use params_mod, only : pi, precis, largest
 198
 199 !-----------------------------------------------------------------------------
 200 !     ... dummy arguments
 201 !-----------------------------------------------------------------------------
 202       INTEGER, intent(in) :: nlyr, nlev
 203       REAL, intent(in)    :: zen, rsfc
 204       REAL, intent(in)    :: tauu(nlyr), omu(nlyr), gu(nlyr)
 205       REAL, intent(in)    :: dsdh(0:nlyr,nlyr)
 206       INTEGER, intent(in) :: nid(0:nlyr)
 207       LOGICAL, intent(in) :: delta
 208
 209       REAL, intent(out) :: fup(nlev), fdn(nlev), fdr(nlev)
 210       REAL, intent(out) :: eup(nlev), edn(nlev), edr(nlev)
 211
 212 !-----------------------------------------------------------------------------
 213 !     ... local variables
 214 !-----------------------------------------------------------------------------
 215       REAL, PARAMETER    :: eps = 1.E-3
 216 !-----------------------------------------------------------------------------
 217 ! initial conditions:  pi*solar flux = 1;  diffuse incidence = 0
 218 !-----------------------------------------------------------------------------
 219       REAL, PARAMETER    :: pifs = 1.
 220       REAL, PARAMETER    :: fdn0 = 0.
 221
 222       REAL :: mu, sum
 223       REAL :: tausla(0:nlyr)
 224       REAL :: tauc(0:nlyr)
 225       REAL :: mu2(0:nlyr)
 226
 227 !-----------------------------------------------------------------------------
 228 ! internal coefficients and matrix
 229 !-----------------------------------------------------------------------------
 230       INTEGER :: row, nlyrm1
 231       REAL :: lam(nlyr), taun(nlyr), bgam(nlyr)
 232       REAL :: e1(nlyr), e2(nlyr), e3(nlyr), e4(nlyr)
 233       REAL :: cup(nlyr), cdn(nlyr), cuptn(nlyr), cdntn(nlyr)
 234       REAL :: mu1(nlyr)
 235       REAL :: a(2*nlyr), b(2*nlyr), d(2*nlyr), e(2*nlyr), y(2*nlyr)
 236
 237       REAL :: f, g, om, tmpg
 238       REAL :: gam1, gam2, gam3, gam4
 239       REAL :: gi(nlyr), omi(nlyr)
 240
 241 !-----------------------------------------------------------------------------
 242 ! For calculations of Associated Legendre Polynomials for GAMA1,2,3,4
 243 ! in delta-function, modified quadrature, hemispheric constant,
 244 ! Hybrid modified Eddington-delta function metods, p633,Table1.
 245 ! W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
 246 ! W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440,
 247 ! uncomment the following two lines and the appropriate statements further
 248 ! down.
 249 !-----------------------------------------------------------------------------
 250       INTEGER :: mrows, mrowsm1, mrowsm2
 251       REAL    :: expon, expon0, expon1, divisr, tmp, up, dn
 252       REAL    :: ssfc
 253       REAL    :: wrk, wrk1
 254
 255       INTEGER :: i, im1, j, k
 256
 257 !-----------------------------------------------------------------------------
 258 ! MU = cosine of solar zenith angle
 259 ! RSFC = surface albedo
 260 ! TAUU =  unscaled optical depth of each layer
 261 ! OMU  =  unscaled single scattering albedo
 262 ! GU   =  unscaled asymmetry factor
 263 ! KLEV = max dimension of number of layers in atmosphere
 264 ! NLYR = number of layers in the atmosphere
 265 ! NLEVEL = nlyr + 1 = number of levels
 266 !-----------------------------------------------------------------------------
 267
 268       mu = COS( zen*pi/180. )
 269 !-----------------------------------------------------------------------------
 270 !************* compute coefficients for each layer:
 271 ! GAM1 - GAM4 = 2-stream coefficients, different for different approximations
 272 ! EXPON0 = calculation of e when TAU is zero
 273 ! EXPON1 = calculation of e when TAU is TAUN
 274 ! CUP and CDN = calculation when TAU is zero
 275 ! CUPTN and CDNTN = calc. when TAU is TAUN
 276 ! DIVISR = prevents division by zero
 277 !-----------------------------------------------------------------------------
 278       tauc(0:nlyr)   = 0.
 279       tausla(0:nlyr) = 0.
 280       mu2(0:nlyr)    = 1./SQRT(largest)
 281
 282       IF( .NOT. delta ) THEN
 283         gi(1:nlyr)   = gu(1:nlyr)
 284         omi(1:nlyr)  = omu(1:nlyr)
 285         taun(1:nlyr) = tauu(1:nlyr)
 286       ELSE
 287 !-----------------------------------------------------------------------------
 288 ! delta-scaling. Have to be done for delta-Eddington approximation,
 289 ! delta discrete ordinate, Practical Improved Flux Method, delta function,
 290 ! and Hybrid modified Eddington-delta function methods approximations
 291 !-----------------------------------------------------------------------------
 292         DO k = 1, nlyr
 293           f       = gu(k)*gu(k)
 294           wrk     = 1. - f
 295           wrk1    = 1. - omu(k)*f
 296           gi(k)   = (gu(k) - f)/wrk
 297           omi(k)  = wrk*omu(k)/wrk1
 298           taun(k) = wrk1*tauu(k)
 299         ENDDO
 300       END IF
 301
 302 !-----------------------------------------------------------------------------
 303 ! calculate slant optical depth at the top of the atmosphere when zen>90.
 304 ! in this case, higher altitude of the top layer is recommended which can
 305 ! be easily changed in gridz.f.
 306 !-----------------------------------------------------------------------------
 307       IF( zen > 90.0 ) THEN
 308         IF(nid(0) < 0) THEN
 309           tausla(0) = largest
 310         ELSE
 311           sum = 0.0
 312           DO j = 1, nid(0)
 313             sum = sum + 2.*taun(j)*dsdh(0,j)
 314           END DO
 315           tausla(0) = sum
 316         END IF
 317       END IF
 318
 319 layer_loop : &
 320       DO i = 1, nlyr
 321         im1 = i - 1
 322         g  = gi(i)
 323         om = omi(i)
 324         tauc(i) = tauc(im1) + taun(i)
 325
 326 !-----------------------------------------------------------------------------
 327 ! stay away from 1 by precision.  For g, also stay away from -1
 328 !-----------------------------------------------------------------------------
 329         tmpg = MIN( abs(g),1. - precis )
 330         g    = SIGN( tmpg,g )
 331         om   = MIN( om,1. - precis )
 332
 333 !-----------------------------------------------------------------------------
 334 ! calculate slant optical depth
 335 !-----------------------------------------------------------------------------
 336         IF(nid(i) < 0) THEN
 337           tausla(i) = largest
 338         ELSE
 339           sum = 0.0
 340           DO j = 1, MIN(nid(i),i)
 341              sum = sum + taun(j)*dsdh(i,j)
 342           ENDDO
 343           DO j = MIN(nid(i),i)+1,nid(i)
 344              sum = sum + 2.*taun(j)*dsdh(i,j)
 345           ENDDO
 346           tausla(i) = sum
 347           IF(tausla(i) == tausla(im1)) THEN
 348             mu2(i) = SQRT(largest)
 349           ELSE
 350             mu2(i) = (tauc(i) - tauc(im1))/(tausla(i) - tausla(im1))
 351             mu2(i) = SIGN( MAX( ABS(mu2(i)),1./SQRT(largest) ), mu2(i) )
 352           END IF
 353         END IF
 354
 355 !-----------------------------------------------------------------------------
 356 !** the following gamma equations are from pg 16,289, Table 1
 357 !** save mu1 for each approx. for use in converting irradiance to actinic flux
 358 ! Eddington approximation(Joseph et al., 1976, JAS, 33, 2452):
 359 !-----------------------------------------------------------------------------
 360         gam1 =  .25*(7. - om*(4. + 3.*g))
 361         gam2 = -.25*(1. - om*(4. - 3.*g))
 362         gam3 = .25*(2. - 3.*g*mu)
 363         gam4 = 1. - gam3
 364         mu1(i) = 0.5
 365
 366 !-----------------------------------------------------------------------------
 367 ! lambda = pg 16,290 equation 21
 368 ! big gamma = pg 16,290 equation 22
 369 ! checked limit for gam2/gam1 <<1:  bgam -> (1/2)*gma2/gam1
 370 ! so if if gam2 = 0., then bgam = 0.
 371 !-----------------------------------------------------------------------------
 372         lam(i) = sqrt(gam1*gam1 - gam2*gam2)
 373
 374         IF( gam2 /= 0.) THEN
 375           bgam(i) = (gam1 - lam(i))/gam2
 376         ELSE
 377           bgam(i) = 0.
 378         ENDIF
 379
 380         expon = EXP( -lam(i)*taun(i) )
 381
 382 !-----------------------------------------------------------------------------
 383 ! e1 - e4 = pg 16,292 equation 44
 384 !-----------------------------------------------------------------------------
 385         e1(i) = 1. + bgam(i)*expon
 386         e2(i) = 1. - bgam(i)*expon
 387         e3(i) = bgam(i) + expon
 388         e4(i) = bgam(i) - expon
 389
 390 !-----------------------------------------------------------------------------
 391 ! the following sets up for the C equations 23, and 24
 392 ! found on page 16,290
 393 ! prevent division by zero (if LAMBDA=1/MU, shift 1/MU^2 by EPS = 1.E-3
 394 ! which is approx equiv to shifting MU by 0.5*EPS* (MU)**3
 395 !-----------------------------------------------------------------------------
 396         expon0 = EXP( -tausla(im1) )
 397         expon1 = EXP( -tausla(i) )
 398
 399         divisr = lam(i)*lam(i) - 1./(mu2(i)*mu2(i))
 400         tmp    = MAX( eps,abs(divisr) )
 401         divisr = SIGN( tmp,divisr )
 402
 403         up = om*pifs*((gam1 - 1./mu2(i))*gam3 + gam4*gam2)/divisr
 404         dn = om*pifs*((gam1 + 1./mu2(i))*gam4 + gam2*gam3)/divisr
 405
 406 !-----------------------------------------------------------------------------
 407 ! cup and cdn are when tau is equal to zero
 408 ! cuptn and cdntn are when tau is equal to taun
 409 !-----------------------------------------------------------------------------
 410         cup(i) = up*expon0
 411         cdn(i) = dn*expon0
 412         cuptn(i) = up*expon1
 413         cdntn(i) = dn*expon1
 414       end do layer_loop
 415
 416 !-----------------------------------------------------------------------------
 417 !**************** set up matrix ******
 418 ! ssfc = pg 16,292 equation 37  where pi Fs is one (unity).
 419 !-----------------------------------------------------------------------------
 420       ssfc = rsfc*mu*EXP( -tausla(nlyr) )*pifs
 421
 422 !-----------------------------------------------------------------------------
 423 ! MROWS = the number of rows in the matrix
 424 !-----------------------------------------------------------------------------
 425       mrows   = 2*nlyr
 426       mrowsm1 = mrows - 1
 427       mrowsm2 = mrows - 2
 428       nlyrm1  = nlyr - 1
 429
 430 !-----------------------------------------------------------------------------
 431 ! the following are from pg 16,292  equations 39 - 43.
 432 ! set up first row of matrix:
 433 !-----------------------------------------------------------------------------
 434       a(1) = 0.
 435       b(1) = e1(1)
 436       d(1) = -e2(1)
 437       e(1) = fdn0 - cdn(1)
 438
 439 !-----------------------------------------------------------------------------
 440 ! set up odd rows 3 thru (MROWS - 1):
 441 !-----------------------------------------------------------------------------
 442       a(3:mrowsm1:2) = e2(1:nlyrm1)*e3(1:nlyrm1) - e4(1:nlyrm1)*e1(1:nlyrm1)
 443       b(3:mrowsm1:2) = e1(1:nlyrm1)*e1(2:nlyr) - e3(1:nlyrm1)*e3(2:nlyr)
 444       d(3:mrowsm1:2) = e3(1:nlyrm1)*e4(2:nlyr) - e1(1:nlyrm1)*e2(2:nlyr)
 445       e(3:mrowsm1:2) = e3(1:nlyrm1)*(cup(2:nlyr) - cuptn(1:nlyrm1))  &
 446                      + e1(1:nlyrm1)*(cdntn(1:nlyrm1) - cdn(2:nlyr))
 447
 448 !-----------------------------------------------------------------------------
 449 ! set up even rows 2 thru (MROWS - 2):
 450 !-----------------------------------------------------------------------------
 451       a(2:mrowsm2:2) = e2(2:nlyr)*e1(1:nlyrm1) - e3(1:nlyrm1)*e4(2:nlyr)
 452       b(2:mrowsm2:2) = e2(1:nlyrm1)*e2(2:nlyr) - e4(1:nlyrm1)*e4(2:nlyr)
 453       d(2:mrowsm2:2) = e1(2:nlyr)*e4(2:nlyr) - e2(2:nlyr)*e3(2:nlyr)
 454       e(2:mrowsm2:2) = (cup(2:nlyr) - cuptn(1:nlyrm1))*e2(2:nlyr) &
 455                      - (cdn(2:nlyr) - cdntn(1:nlyrm1))*e4(2:nlyr)
 456
 457 !-----------------------------------------------------------------------------
 458 ! set up last row of matrix at MROWS:
 459 !-----------------------------------------------------------------------------
 460       a(mrows) = e1(nlyr) - rsfc*e3(nlyr)
 461       b(mrows) = e2(nlyr) - rsfc*e4(nlyr)
 462       d(mrows) = 0.
 463       e(mrows) = ssfc - cuptn(nlyr) + rsfc*cdntn(nlyr)
 464
 465 !-----------------------------------------------------------------------------
 466 ! solve tri-diagonal matrix:
 467 !-----------------------------------------------------------------------------
 468       CALL tridag( a, b, d, e, y, mrows )
 469
 470 !-----------------------------------------------------------------------------
 471 !*** unfold solution of matrix, compute output fluxes:
 472 !-----------------------------------------------------------------------------
 473 ! the following equations are from pg 16,291  equations 31 & 32
 474 !-----------------------------------------------------------------------------
 475       fdr(1) = EXP( -tausla(0) )
 476       edr(1) = mu * fdr(1)
 477       edn(1) = fdn0
 478       eup(1) = y(1)*e3(1) - y(2)*e4(1) + cup(1)
 479       fdn(1) = edn(1)/mu1(1)
 480       fup(1) = eup(1)/mu1(1)
 481
 482       fdr(2:nlev) = EXP( -tausla(1:nlyr) )
 483       edr(2:nlev) = mu *fdr(2:nlev)
 484       edn(2:nlev) = y(1:mrowsm1:2)*e3(1:nlyr) + y(2:mrows:2)*e4(1:nlyr) + cdntn(1:nlyr)
 485       eup(2:nlev) = y(1:mrowsm1:2)*e1(1:nlyr) + y(2:mrows:2)*e2(1:nlyr) + cuptn(1:nlyr)
 486       fdn(2:nlev) = edn(2:nlev)/mu1(1:nlyr)
 487       fup(2:nlev) = eup(2:nlev)/mu1(1:nlyr)
 488
 489       END SUBROUTINE ps2str
 490
 491       SUBROUTINE tridag( a, b, c, r, u, n)
 492 !-----------------------------------------------------------------------------
 493 ! solve tridiagonal system.  From Numerical Recipies, p. 40
 494 !-----------------------------------------------------------------------------
 495
 496 !-----------------------------------------------------------------------------
 497 !     ... dummy arguments
 498 !-----------------------------------------------------------------------------
 499       INTEGER, intent(in) :: n
 500       REAL,    intent(in) :: a(n), b(n), c(n), r(n)
 501       REAL, intent(out)   :: u(n)
 502
 503 !-----------------------------------------------------------------------------
 504 !     ... local variables
 505 !-----------------------------------------------------------------------------
 506       INTEGER :: j, jm1, jp1
 507       REAL    :: bet
 508       REAL    :: gam(n)
 509       CHARACTER(len=64) :: err_msg
 510
 511       IF (b(1) == 0.) then
 512         call wrf_error_fatal( 'tridag: zero pivot @ n == 1' )
 513       ENDIF
 514       bet  = b(1)
 515       u(1) = r(1)/bet
 516       DO j = 2, n
 517          jm1 = j - 1
 518          gam(j) = c(jm1)/bet
 519          bet    = b(j) - a(j)*gam(j)
 520          IF (bet == 0.) then
 521            write(err_msg,'(''tridag: zero pivot @ n = '',i4)') j
 522            call wrf_error_fatal( trim(err_msg) )
 523          ENDIF
 524          u(j) = (r(j) - a(j)*u(jm1))/bet
 525       END DO
 526
 527       DO j = n - 1, 1, -1
 528          jp1 = j + 1
 529          u(j) = u(j) - gam(jp1)*u(jp1)
 530       END DO
 531
 532       END SUBROUTINE tridag
 533
 534       END MODULE rad_trans