| 1 |
C $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/monod/monod_radtrans_direct.F,v 1.2 2012/08/24 19:45:36 jahn Exp $ |
| 2 |
C $Name: $ |
| 3 |
|
| 4 |
#include "DARWIN_OPTIONS.h" |
| 5 |
|
| 6 |
CBOP |
| 7 |
C !ROUTINE: MONOD_RADTRANS_DIRECT |
| 8 |
|
| 9 |
C !INTERFACE: ========================================================== |
| 10 |
subroutine MONOD_RADTRANS_DIRECT( |
| 11 |
I H,rmud,Edsf,Essf,a_k,bt_k,bb_k,kmax, |
| 12 |
O Edbot,Esbot,Eubot,Estop,Eutop, |
| 13 |
O tirrq,tirrwq, |
| 14 |
O amp1, amp2, |
| 15 |
I myThid) |
| 16 |
|
| 17 |
C !DESCRIPTION: |
| 18 |
c |
| 19 |
c Model of irradiance in the water column. Accounts for three |
| 20 |
c irradiance streams [Ackleson, Balch, Holligan, JGR, 1994], |
| 21 |
c |
| 22 |
c Edbot = direct downwelling irradiance in W/m2 per waveband |
| 23 |
c Esbot = diffuse downwelling irradiance in W/m2 per waveband |
| 24 |
c Eubot = diffuse upwelling irradiance in W/m2 per waveband |
| 25 |
c |
| 26 |
c all defined at the bottom of each layer. Also computed are Estop, |
| 27 |
c Eutop at the top of each layer which should be very close to Esbot, |
| 28 |
c Eubot of the layer above. |
| 29 |
c |
| 30 |
c The Ed equation is integrated exactly, Es and Eu are computed by |
| 31 |
c solving a set of linear equation for the amplitudes in the exact |
| 32 |
c solution [see, e.g., Kylling, Stamnes, Tsay, JAC, 1995]. |
| 33 |
c The boundary condition in the deepest wet layer is |
| 34 |
c downward-decreasing modes only (i.e., zero irradiance at infinite |
| 35 |
c depth, assuming the optical properties of the last layer). |
| 36 |
c |
| 37 |
c Also computed are scalar radiance and PAR at the grid cell center |
| 38 |
c (both in uEin/m2/s). |
| 39 |
c |
| 40 |
C !USES: =============================================================== |
| 41 |
IMPLICIT NONE |
| 42 |
#include "SIZE.h" /* Nr */ |
| 43 |
#include "EEPARAMS.h" |
| 44 |
#include "MONOD_SIZE.h" |
| 45 |
#include "SPECTRAL_SIZE.h" /* tlam */ |
| 46 |
#include "SPECTRAL.h" /* WtouEin */ |
| 47 |
#include "WAVEBANDS_PARAMS.h" /* darwin_PAR_ilamLo/Hi |
| 48 |
darwin_radmodThresh |
| 49 |
darwin_rmus darwin_rmuu */ |
| 50 |
|
| 51 |
C !INPUT PARAMETERS: =================================================== |
| 52 |
C H :: layer thickness (including hFacC!) |
| 53 |
C rmud :: inv.cosine of direct (underwater solar) zenith angle |
| 54 |
C Edsf :: direct downwelling irradiance below surface per waveband |
| 55 |
C Essf :: diffuse downwelling irradiance below surface per waveband |
| 56 |
C a_k :: absorption coefficient per level and waveband (1/m) |
| 57 |
C bt_k :: total scattering coefficient per level and waveband (1/m) |
| 58 |
C = forward + back scattering coefficient |
| 59 |
C bb_k :: backscattering coefficient per level and waveband (1/m) |
| 60 |
C kmax :: maximum number of layers to compute |
| 61 |
_RL H(Nr) |
| 62 |
_RL rmud |
| 63 |
_RL Edsf(tlam), Essf(tlam) |
| 64 |
_RL a_k(Nr,tlam), bt_k(Nr,tlam), bb_k(Nr,tlam) |
| 65 |
INTEGER kmax |
| 66 |
INTEGER myThid |
| 67 |
|
| 68 |
C !OUTPUT PARAMETERS: ================================================== |
| 69 |
C Edbot :: direct downwelling irradiance at bottom of layer |
| 70 |
C Esbot :: diffuse downwelling irradiance at bottom of layer |
| 71 |
C Eubot :: diffuse upwelling irradiance at bottom of layer |
| 72 |
C Estop :: diffuse downwelling irradiance at top of layer |
| 73 |
C Eutop :: diffuse upwelling irradiance at top of layer |
| 74 |
C tirrq :: total scalar irradiance at cell center (uEin/m2/s) |
| 75 |
C tirrwq :: total scalar irradiance at cell center per waveband |
| 76 |
C amp1 :: amplitude of downward increasing mode |
| 77 |
C amp2 :: amplitude of downward decreasing mode |
| 78 |
_RL Edbot(tlam,Nr),Esbot(tlam,Nr),Eubot(tlam,Nr) |
| 79 |
_RL Estop(tlam,Nr),Eutop(tlam,Nr) |
| 80 |
_RL tirrq(Nr) |
| 81 |
_RL tirrwq(tlam,Nr) |
| 82 |
_RL amp1(tlam,Nr), amp2(tlam,Nr) |
| 83 |
CEOP |
| 84 |
|
| 85 |
#ifdef DAR_RADTRANS |
| 86 |
|
| 87 |
C !LOCAL VARIABLES: ==================================================== |
| 88 |
INTEGER k, nl, kbot |
| 89 |
_RL Edtop(tlam,Nr) |
| 90 |
_RL Etopwq, Ebotwq |
| 91 |
_RL zd |
| 92 |
_RL rmus,rmuu |
| 93 |
_RL cd,au,Bu,Cu |
| 94 |
_RL as,Bs,Cs,Bd,Fd |
| 95 |
_RL bquad,D |
| 96 |
_RL kappa1,kappa2,denom |
| 97 |
_RL c1,c2 |
| 98 |
_RL r2(Nr),r1(Nr),x(Nr),y(Nr) |
| 99 |
_RL ed(Nr),e2(Nr),e1(Nr) |
| 100 |
_RL a3d(2*Nr), b3d(2*Nr), c3d(2*Nr), y3d(2*Nr) |
| 101 |
_RL rd, ru |
| 102 |
data rd /1.5 _d 0/ !these are taken from Ackleson, et al. 1994 (JGR) |
| 103 |
data ru /3.0 _d 0/ |
| 104 |
|
| 105 |
rmus = darwin_rmus |
| 106 |
rmuu = darwin_rmuu |
| 107 |
|
| 108 |
c find deepest wet layer |
| 109 |
kbot = MIN(kmax, Nr) |
| 110 |
DO WHILE (H(kbot).EQ.0 .AND. kbot.GT.1) |
| 111 |
kbot = kbot - 1 |
| 112 |
ENDDO |
| 113 |
IF (H(kbot).EQ.0) kbot = kbot - 1 |
| 114 |
|
| 115 |
DO nl = 1,tlam |
| 116 |
DO k=1,Nr |
| 117 |
Edtop(nl,k) = 0.0 |
| 118 |
Estop(nl,k) = 0.0 |
| 119 |
Eutop(nl,k) = 0.0 |
| 120 |
Edbot(nl,k) = 0.0 |
| 121 |
Esbot(nl,k) = 0.0 |
| 122 |
Eubot(nl,k) = 0.0 |
| 123 |
amp1(nl,k) = 0.0 |
| 124 |
amp2(nl,k) = 0.0 |
| 125 |
ENDDO |
| 126 |
ENDDO |
| 127 |
IF (kbot.GT.0) THEN |
| 128 |
DO nl=1,tlam |
| 129 |
IF (Edsf(nl) .GE. darwin_radmodThresh .OR. |
| 130 |
& Essf(nl) .GE. darwin_radmodThresh) THEN |
| 131 |
DO k=1,kbot |
| 132 |
zd = H(k) |
| 133 |
cd = (a_k(k,nl)+bt_k(k,nl))*rmud |
| 134 |
au = a_k(k,nl)*rmuu |
| 135 |
Bu = ru*bb_k(k,nl)*rmuu |
| 136 |
Cu = au+Bu |
| 137 |
as = a_k(k,nl)*rmus |
| 138 |
Bs = rd*bb_k(k,nl)*rmus |
| 139 |
Cs = as+Bs |
| 140 |
Bd = bb_k(k,nl)*rmud |
| 141 |
Fd = (bt_k(k,nl)-bb_k(k,nl))*rmud |
| 142 |
bquad = Cs + Cu |
| 143 |
D = 0.5*(bquad + SQRT(bquad*bquad - 4.0*Bs*Bu)) |
| 144 |
kappa1 = D - Cs |
| 145 |
kappa2 = Cs - Bs*Bu/D ! == D - Cu |
| 146 |
r1(k) = Bu/D |
| 147 |
r2(k) = Bs/D |
| 148 |
denom = (cd-Cs)*(cd+Cu) + Bs*Bu |
| 149 |
x(k) = -((cd+Cu)*Fd+Bu*Bd)/denom |
| 150 |
y(k) = (-Bs*Fd+(cd-Cs)*Bd)/denom |
| 151 |
ed(k) = EXP(-cd*zd) |
| 152 |
e1(k) = EXP(-kappa1*zd) |
| 153 |
e2(k) = EXP(-kappa2*zd) |
| 154 |
ENDDO |
| 155 |
|
| 156 |
C integrate Ed equation first |
| 157 |
Edtop(nl,1) = Edsf(nl) |
| 158 |
DO k=1,kbot-1 |
| 159 |
Edbot(nl,k) = Edtop(nl,k)*ed(k) |
| 160 |
Edtop(nl,k+1) = Edbot(nl,k) |
| 161 |
ENDDO |
| 162 |
Edbot(nl,kbot) = Edtop(nl,kbot)*ed(kbot) |
| 163 |
|
| 164 |
C setup tridiagonal matrix of continuity/boundary conditions |
| 165 |
C variables: c2(1), c1(1), c2(2), ..., c1(kbot) |
| 166 |
C a3d,b3d,c3d: lower, main and upper diagonal |
| 167 |
C y3d: right-hand side |
| 168 |
C |
| 169 |
C top b.c.: c2(1) + e1(1)*r1(1)*c1(1) = Essf - x(1)*Edsf |
| 170 |
a3d(1) = 0. _d 0 ! not used |
| 171 |
b3d(1) = 1. ! A(1,1)*c2(1) |
| 172 |
c3d(1) = e1(1)*r1(1) ! A(1,2)*c1(1) |
| 173 |
y3d(1) = Essf(nl) - x(1)*Edsf(nl) |
| 174 |
C continuity at layer boundaries |
| 175 |
DO k=1, kbot-1 |
| 176 |
a3d(2*k) = (1. - r2(k)*r1(k+1))*e2(k) ! A(2k,2k-1)*c2(k) |
| 177 |
b3d(2*k) = r1(k) - r1(k+1) ! + A(2k,2k )*c1(k) |
| 178 |
c3d(2*k) = -1. + r2(k+1)*r1(k+1) ! + A(2k,2k+1)*c2(k+1) |
| 179 |
y3d(2*k)= (x(k+1) - x(k) - r1(k+1)*(y(k+1)-y(k)))*Edbot(nl,k) |
| 180 |
a3d(2*k+1) = 1 - r1(k)*r2(k) ! A(2k+1,2k )*c1(k) |
| 181 |
b3d(2*k+1) = r2(k) - r2(k+1) ! + A(2k+1,2k+1)*c2(k+1) |
| 182 |
c3d(2*k+1) = (-1. + r1(k+1)*r2(k))*e1(k+1) ! + A(2k+1,2k+2)*c1(k+1) |
| 183 |
y3d(2*k+1)= (y(k+1) - y(k) - r2(k)*(x(k+1)-x(k)))*Edbot(nl,k) |
| 184 |
ENDDO |
| 185 |
c bottom boundary condition: c1 = 0 |
| 186 |
a3d(2*kbot) = 0. _d 0 ! A(2*kbot,2*kbot-1)*c2(kbot) |
| 187 |
b3d(2*kbot) = 1. _d 0 ! + A(2*kbot,2*kbot )*c1(kbot) |
| 188 |
c3d(2*kbot) = 0. _d 0 ! not used |
| 189 |
y3d(2*kbot) = 0. _d 0 ! = 0 |
| 190 |
|
| 191 |
CALL SOLVE_TRIDIAGONAL_PIVOT(a3d,b3d,c3d,y3d,2*kbot,myThid) |
| 192 |
|
| 193 |
C compute irradiances |
| 194 |
DO k=1,kbot |
| 195 |
c2 = y3d(2*k-1) |
| 196 |
c1 = y3d(2*k) |
| 197 |
Estop(nl,k) = c2 + r1(k)*e1(k)*c1 + x(k)*Edtop(nl,k) |
| 198 |
Esbot(nl,k) = e2(k)*c2 + r1(k)*c1 + x(k)*Edbot(nl,k) |
| 199 |
Eutop(nl,k) = r2(k)*c2 + e1(k)*c1 + y(k)*Edtop(nl,k) |
| 200 |
Eubot(nl,k) = r2(k)*e2(k)*c2 + c1 + y(k)*Edbot(nl,k) |
| 201 |
amp1(nl,k) = c1 |
| 202 |
amp2(nl,k) = c2 |
| 203 |
ENDDO |
| 204 |
IF (kbot .LT. Nr) THEN |
| 205 |
Estop(nl,kbot+1) = Esbot(nl,kbot) |
| 206 |
Eutop(nl,kbot+1) = Eubot(nl,kbot) |
| 207 |
ENDIF |
| 208 |
|
| 209 |
C endif thresh |
| 210 |
ENDIF |
| 211 |
|
| 212 |
DO k = 1,Nr |
| 213 |
C convert to scalar irradiance in quanta |
| 214 |
#ifdef DAR_RADTRANS_RMUS_PAR |
| 215 |
C use rmus for all 3 components (?) |
| 216 |
Etopwq = (Edtop(nl,k)+Estop(nl,k)+Eutop(nl,k))*WtouEins(nl) |
| 217 |
Ebotwq = (Edbot(nl,k)+Esbot(nl,k)+Eubot(nl,k))*WtouEins(nl) |
| 218 |
tirrwq(nl,k) = SQRT(Etopwq*Ebotwq)*rmus |
| 219 |
#else |
| 220 |
C use appropriate average cosine for each component |
| 221 |
Etopwq = (rmud*Edtop(nl,k)+rmus*Estop(nl,k)+rmuu*Eutop(nl,k)) |
| 222 |
& *WtouEins(nl) |
| 223 |
Ebotwq = (rmud*Edbot(nl,k)+rmus*Esbot(nl,k)+rmuu*Eubot(nl,k)) |
| 224 |
& *WtouEins(nl) |
| 225 |
C and interpolate |
| 226 |
tirrwq(nl,k) = SQRT(Etopwq*Ebotwq) |
| 227 |
#endif |
| 228 |
ENDDO |
| 229 |
|
| 230 |
C enddo nl |
| 231 |
ENDDO |
| 232 |
C endif kbot.gt.0 |
| 233 |
ENDIF |
| 234 |
|
| 235 |
DO k = 1,Nr |
| 236 |
C sum PAR range |
| 237 |
tirrq(k) = 0.0 |
| 238 |
DO nl = darwin_PAR_ilamLo,darwin_PAR_ilamHi |
| 239 |
tirrq(k) = tirrq(k) + tirrwq(nl,k) |
| 240 |
ENDDO |
| 241 |
ENDDO |
| 242 |
c |
| 243 |
#endif /* DAR_RADTRANS */ |
| 244 |
|
| 245 |
return |
| 246 |
end |
| 247 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph