pkg/radtrans/radtrans_solve.F

b55e95f1ff Oliv*0001 #include "RADTRANS_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: RADTRANS_DIRECT
                0005
                0006 C !INTERFACE: ==========================================================
                0007       subroutine RADTRANS_SOLVE(
                0008      I                   H,rmud,Edsf,Essf,a_k,bt_k,bb_k,kbot,
                0009      O                   Edbot,Esbot,Eubot,Estop,Eutop,
                0010      O                   amp1, amp2, x, y,
                0011      O                   r1, r2, kappa1, kappa2,
                0012      I                   myThid)
                0013
                0014 C !DESCRIPTION:
                0015 c
                0016 c  Model of irradiance in the water column.  Accounts for three
                0017 c  irradiance streams [Ackleson, Balch, Holligan, JGR, 1994],
                0018 c
                0019 c  Edbot = direct downwelling irradiance in W/m2 per waveband
                0020 c  Esbot = diffuse downwelling irradiance in W/m2 per waveband
                0021 c  Eubot = diffuse upwelling irradiance in W/m2 per waveband
                0022 c
                0023 c  all defined at the bottom of each layer.  Also computed are Estop,
                0024 c  Eutop at the top of each layer which should be very close to Esbot,
                0025 c  Eubot of the layer above.
                0026 c
                0027 c  The Ed equation is integrated exactly, Es and Eu are computed by
                0028 c  solving a set of linear equation for the amplitudes in the exact
                0029 c  solution [see, e.g., Kylling, Stamnes, Tsay, JAC, 1995].
                0030 c  The boundary condition in the deepest wet layer is
                0031 c  downward-decreasing modes only (i.e., zero irradiance at infinite
                0032 c  depth, assuming the optical properties of the last layer).
                0033 c
                0034 c  Also computed are scalar radiance and PAR at the grid cell center
                0035 c  (both in uEin/m2/s).
                0036 c
                0037 C !USES: ===============================================================
                0038       IMPLICIT NONE
                0039 #include "SIZE.h"
                0040 #include "EEPARAMS.h"
                0041 #include "RADTRANS_SIZE.h"
                0042 #include "RADTRANS_PARAMS.h"
                0043
                0044 C !INPUT PARAMETERS: ===================================================
                0045 C     H     :: layer thickness (including hFacC!)
                0046 C     rmud  :: inv.cosine of direct (underwater solar) zenith angle
                0047 C     Edsf  :: direct downwelling irradiance below surface per waveband
                0048 C     Essf  :: diffuse downwelling irradiance below surface per waveband
                0049 C     a_k   :: absorption coefficient per level and waveband (1/m)
                0050 C     bt_k  :: total scattering coefficient per level and waveband (1/m)
                0051 C              = forward + back scattering coefficient
                0052 C     bb_k  :: backscattering coefficient per level and waveband (1/m)
                0053 C     kbot  :: maximum number of layers to compute
                0054       _RL H(Nr)
                0055       _RL rmud
                0056       _RL Edsf, Essf
                0057       _RL a_k(Nr), bt_k(Nr), bb_k(Nr)
                0058       INTEGER kbot
                0059       INTEGER myThid
                0060
                0061 C !OUTPUT PARAMETERS: ==================================================
                0062 C     Edbot  :: direct downwelling irradiance at bottom of layer
                0063 C     Esbot  :: diffuse downwelling irradiance at bottom of layer
                0064 C     Eubot  :: diffuse upwelling irradiance at bottom of layer
                0065 C     Estop  :: diffuse downwelling irradiance at top of layer
                0066 C     Eutop  :: diffuse upwelling irradiance at top of layer
                0067 C     amp1   :: amplitude of downward increasing mode
                0068 C     amp2   :: amplitude of downward decreasing mode
                0069       _RL Edbot(Nr),Esbot(Nr),Eubot(Nr)
                0070       _RL Estop(Nr),Eutop(Nr)
                0071       _RL amp1(Nr), amp2(Nr)
                0072       _RL x(Nr), y(Nr)
                0073       _RL kappa1(Nr),kappa2(Nr)
                0074       _RL r2(Nr),r1(Nr)
                0075 CEOP
                0076
                0077 #ifdef ALLOW_RADTRANS
                0078
                0079 C !LOCAL VARIABLES: ====================================================
                0080       INTEGER k
                0081       _RL Edtop(Nr)
                0082       _RL Etopwq, Ebotwq
                0083       _RL zd
                0084       _RL rmus,rmuu
                0085       _RL cd,au,Bu,Cu
                0086       _RL as,Bs,Cs,Bd,Fd
                0087       _RL bquad,D
                0088       _RL denom
                0089       _RL c1,c2
                0090       _RL ed(Nr),e2(Nr),e1(Nr)
                0091       _RL a3d(2*Nr), b3d(2*Nr), c3d(2*Nr), y3d(2*Nr)
                0092
                0093       rmus = RT_rmus
                0094       rmuu = RT_rmuu
                0095
                0096       DO k=1,Nr
                0097         Edtop(k) = 0.0
                0098         Estop(k) = 0.0
                0099         Eutop(k) = 0.0
                0100         Edbot(k) = 0.0
                0101         Esbot(k) = 0.0
                0102         Eubot(k) = 0.0
                0103         amp1(k) = 0.0
                0104         amp2(k) = 0.0
                0105         kappa1(k) = 0.0
                0106         kappa2(k) = 0.0
                0107         r1(k) = 0.0
                0108         r2(k) = 0.0
                0109         x(k) = 0.0
                0110         y(k) = 0.0
                0111       ENDDO
                0112       IF (kbot.GT.0 .AND.
                0113      &    (Edsf.GE.RT_sfcIrrThresh .OR. Essf.GE.RT_sfcIrrThresh)) THEN
                0114         DO k=1,kbot
                0115           zd = H(k)
                0116           cd = (a_k(k)+bt_k(k))*rmud
                0117           au = a_k(k)*rmuu
                0118           Bu = RT_ru*bb_k(k)*rmuu
                0119           Cu = au+Bu
                0120           as = a_k(k)*rmus
                0121           Bs = RT_rd*bb_k(k)*rmus
                0122           Cs = as+Bs
                0123           Bd = bb_k(k)*rmud
                0124           Fd = (bt_k(k)-bb_k(k))*rmud
                0125           bquad = Cs + Cu
                0126           D = 0.5*(bquad + SQRT(bquad*bquad - 4.0*Bs*Bu))
                0127           kappa1(k) = D - Cs
                0128           kappa2(k) = Cs - Bs*Bu/D  ! == D - Cu
                0129           r1(k) = Bu/D
                0130           r2(k) = Bs/D
                0131           denom = (cd-Cs)*(cd+Cu) + Bs*Bu
                0132           x(k) = -((cd+Cu)*Fd+Bu*Bd)/denom
                0133           y(k) = (-Bs*Fd+(cd-Cs)*Bd)/denom
                0134           ed(k) = EXP(-cd*zd)
                0135           e1(k) = EXP(-kappa1(k)*zd)
                0136           e2(k) = EXP(-kappa2(k)*zd)
                0137         ENDDO
                0138
                0139 C integrate Ed equation first
                0140         Edtop(1) = Edsf
                0141         DO k=1,kbot-1
                0142           Edbot(k) = Edtop(k)*ed(k)
                0143           Edtop(k+1) = Edbot(k)
                0144         ENDDO
                0145         Edbot(kbot) = Edtop(kbot)*ed(kbot)
                0146
                0147 C setup tridiagonal matrix of continuity/boundary conditions
                0148 C variables: c2(1), c1(1), c2(2), ..., c1(kbot)
                0149 C a3d,b3d,c3d: lower, main and upper diagonal
                0150 C y3d: right-hand side
                0151 C
                0152 C top b.c.: c2(1) + e1(1)*r1(1)*c1(1) = Essf - x(1)*Edsf
                0153         a3d(1) = 0. _d 0  ! not used
                0154         b3d(1) = 1.           ! A(1,1)*c2(1)
                0155         c3d(1) = e1(1)*r1(1)  ! A(1,2)*c1(1)
                0156         y3d(1) = Essf - x(1)*Edsf
                0157 C continuity at layer boundaries
                0158         DO k=1, kbot-1
                0159           a3d(2*k) = (1. - r2(k)*r1(k+1))*e2(k)  !   A(2k,2k-1)*c2(k)
                0160           b3d(2*k) = r1(k) - r1(k+1)             ! + A(2k,2k  )*c1(k)
                0161           c3d(2*k) = -1. + r2(k+1)*r1(k+1)       ! + A(2k,2k+1)*c2(k+1)
                0162           y3d(2*k)=(x(k+1) - x(k) - r1(k+1)*(y(k+1)-y(k)))
                0163      &             *Edbot(k)
                0164           a3d(2*k+1) = 1 - r1(k)*r2(k)                !   A(2k+1,2k  )*c1(k)
                0165           b3d(2*k+1) = r2(k) - r2(k+1)                ! + A(2k+1,2k+1)*c2(k+1)
                0166           c3d(2*k+1) = (-1. + r1(k+1)*r2(k))*e1(k+1)  ! + A(2k+1,2k+2)*c1(k+1)
                0167           y3d(2*k+1)=(y(k+1) - y(k) - r2(k)*(x(k+1)-x(k)))
                0168      &               *Edbot(k)
                0169         ENDDO
                0170 c bottom boundary condition: c1 = 0
                0171         a3d(2*kbot) = 0. _d 0  !   A(2*kbot,2*kbot-1)*c2(kbot)
                0172         b3d(2*kbot) = 1. _d 0  ! + A(2*kbot,2*kbot  )*c1(kbot)
                0173         c3d(2*kbot) = 0. _d 0  ! not used
                0174         y3d(2*kbot) = 0. _d 0  ! = 0
                0175
                0176         CALL RADTRANS_SOLVE_TRIDIAG(a3d,b3d,c3d,y3d,2*kbot,myThid)
                0177
                0178 C compute irradiances
                0179         DO k=1,kbot
                0180           c2 = y3d(2*k-1)
                0181           c1 = y3d(2*k)
                0182           Estop(k) = c2 + r1(k)*e1(k)*c1 + x(k)*Edtop(k)
                0183           Esbot(k) = e2(k)*c2 + r1(k)*c1 + x(k)*Edbot(k)
                0184           Eutop(k) = r2(k)*c2 + e1(k)*c1 + y(k)*Edtop(k)
                0185           Eubot(k) = r2(k)*e2(k)*c2 + c1 + y(k)*Edbot(k)
                0186           amp1(k) = c1
                0187           amp2(k) = c2
                0188         ENDDO
                0189         IF (kbot .LT. Nr) THEN
                0190           Estop(kbot+1) = Esbot(kbot)
                0191           Eutop(kbot+1) = Eubot(kbot)
                0192         ENDIF
                0193
                0194 C     endif kbot.gt.0
                0195       ENDIF
                0196
                0197 #endif /* ALLOW_RADTRANS */
                0198
                0199       RETURN
                0200       END
                0201