| 1 |
C $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/monod/monod_radtrans_iter.F,v 1.2 2012/07/30 15:21:51 jahn Exp $ |
| 2 |
C $Name: $ |
| 3 |
|
| 4 |
#include "DARWIN_OPTIONS.h" |
| 5 |
|
| 6 |
CBOP |
| 7 |
C !ROUTINE: MONOD_RADTRANS_ITER |
| 8 |
|
| 9 |
C !INTERFACE: ========================================================== |
| 10 |
subroutine MONOD_RADTRANS_ITER( |
| 11 |
I H,rmud,Edsf,Essf,a_k,bt_k,bb_k,kmax,niter, |
| 12 |
O Edbot,Esbot,Eubot,Eutop, |
| 13 |
O tirrq,tirrwq, |
| 14 |
O c1out, c2out, |
| 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: |
| 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 Propagation is done in energy units, tests are done in quanta, |
| 27 |
c final is quanta for phytoplankton growth. |
| 28 |
c |
| 29 |
c The Ed equation is integrated exactly. |
| 30 |
c Es and Eu are first computed using a truncation to downward- |
| 31 |
c decreasing modes a la Aas that makes Es continuous. |
| 32 |
c Then niter alternating upward and downward integrations are performed, |
| 33 |
c each time using Es at the top and Eu at the bottom of each layer as a |
| 34 |
c boundary condition. The boundary condition in the deepest wet layer |
| 35 |
c is always downward-decreasing modes only. |
| 36 |
c During upward integrations, Eu is made continuous, during downward |
| 37 |
c integrations, Es. |
| 38 |
c At the end, Ed and Es are continuous, but Eu is so only approximately. |
| 39 |
c |
| 40 |
C !USES: =============================================================== |
| 41 |
IMPLICIT NONE |
| 42 |
#include "SIZE.h" /* Nr */ |
| 43 |
C#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 |
C niter :: number of up-down iterations after initial Aas integration |
| 62 |
_RL H(Nr) |
| 63 |
_RL rmud |
| 64 |
_RL Edsf(tlam), Essf(tlam) |
| 65 |
_RL a_k(Nr,tlam), bt_k(Nr,tlam), bb_k(Nr,tlam) |
| 66 |
INTEGER kmax,niter |
| 67 |
INTEGER myThid |
| 68 |
|
| 69 |
C !OUTPUT PARAMETERS: ================================================== |
| 70 |
C Edbot :: direct downwelling irradiance at bottom of layer |
| 71 |
C Esbot :: diffuse downwelling irradiance at bottom of layer |
| 72 |
C Eubot :: diffuse upwelling irradiance at bottom of layer |
| 73 |
C tirrq :: total scalar irradiance at cell center (uEin/m2/s) |
| 74 |
C tirrwq :: total scalar irradiance at cell center per waveband |
| 75 |
_RL Edbot(tlam,Nr),Esbot(tlam,Nr),Eubot(tlam,Nr),Eutop(tlam,Nr) |
| 76 |
_RL tirrq(Nr) |
| 77 |
_RL tirrwq(tlam,Nr) |
| 78 |
_RL c1out(tlam,Nr), c2out(tlam,Nr) |
| 79 |
|
| 80 |
#ifdef DAR_RADTRANS |
| 81 |
|
| 82 |
C !LOCAL VARIABLES: ==================================================== |
| 83 |
INTEGER k, nl, iter, kbot |
| 84 |
_RL Edtop(tlam,Nr),Estop(tlam,Nr) |
| 85 |
_RL Etopwq, Ebotwq |
| 86 |
_RL zd |
| 87 |
_RL rmus,rmuu |
| 88 |
|
| 89 |
C !LOCAL VARIABLES: ================================================ |
| 90 |
_RL cd,au,Bu,Cu |
| 91 |
_RL as,Bs,Cs,Bd,Fd |
| 92 |
_RL bquad,cquad,sqarg |
| 93 |
_RL a1,a2,denom |
| 94 |
_RL c1,c2,tmp,Esnew,Eunew |
| 95 |
_RL R2(Nr),R1(Nr),x(Nr),y(Nr) |
| 96 |
_RL expAddr(Nr),expAsdr(Nr),expmAudr(Nr),idenom(Nr) |
| 97 |
c |
| 98 |
_RL rbot, rd, ru |
| 99 |
data rbot /0.0/ !bottom reflectance (not used) |
| 100 |
data rd /1.5/ !these are taken from Ackleson, et al. 1994 (JGR) |
| 101 |
data ru /3.0/ |
| 102 |
CEOP |
| 103 |
rmus = darwin_rmus |
| 104 |
rmuu = darwin_rmuu |
| 105 |
|
| 106 |
c find deepest wet layer |
| 107 |
kbot = MIN(kmax, Nr) |
| 108 |
DO WHILE (H(kbot).EQ.0 .AND. kbot.GT.1) |
| 109 |
kbot = kbot - 1 |
| 110 |
ENDDO |
| 111 |
|
| 112 |
DO nl = 1,tlam |
| 113 |
DO k=1,Nr |
| 114 |
Edtop(nl,k) = 0.0 |
| 115 |
Estop(nl,k) = 0.0 |
| 116 |
Eutop(nl,k) = 0.0 |
| 117 |
Edbot(nl,k) = 0.0 |
| 118 |
Esbot(nl,k) = 0.0 |
| 119 |
Eubot(nl,k) = 0.0 |
| 120 |
c1out(nl,k) = 0.0 |
| 121 |
c2out(nl,k) = 0.0 |
| 122 |
ENDDO |
| 123 |
IF (Edsf(nl) .GE. darwin_radmodThresh .OR. |
| 124 |
& Essf(nl) .GE. darwin_radmodThresh) THEN |
| 125 |
DO k=1,kbot |
| 126 |
zd = H(k) |
| 127 |
cd = (a_k(k,nl)+bt_k(k,nl))*rmud |
| 128 |
au = a_k(k,nl)*rmuu |
| 129 |
Bu = ru*bb_k(k,nl)*rmuu |
| 130 |
Cu = au+Bu |
| 131 |
as = a_k(k,nl)*rmus |
| 132 |
Bs = rd*bb_k(k,nl)*rmus |
| 133 |
Cs = as+Bs |
| 134 |
Bd = bb_k(k,nl)*rmud |
| 135 |
Fd = (bt_k(k,nl)-bb_k(k,nl))*rmud |
| 136 |
bquad = Cs - Cu |
| 137 |
cquad = Bs*Bu - Cs*Cu |
| 138 |
sqarg = bquad*bquad - 4.0*cquad |
| 139 |
a1 = 0.5*(-bquad + sqrt(sqarg)) |
| 140 |
a2 = 0.5*(-bquad - sqrt(sqarg)) ! K of Aas |
| 141 |
R1(k) = (a1+Cs)/Bu |
| 142 |
R2(k) = (a2+Cs)/Bu |
| 143 |
denom = (cd-Cs)*(cd+Cu) + Bs*Bu |
| 144 |
x(k) = -((cd+Cu)*Fd+Bu*Bd)/denom |
| 145 |
y(k) = (-Bs*Fd+(cd-Cs)*Bd)/denom |
| 146 |
expAddr(k) = exp(-cd*zd) |
| 147 |
expmAudr(k) = exp(-a1*zd) |
| 148 |
expAsdr(k) = exp(a2*zd) |
| 149 |
idenom(k) = 1./(R1(k)-R2(k)*expAsdr(k)*expmAudr(k)) |
| 150 |
ENDDO |
| 151 |
|
| 152 |
C integrate Ed equation first |
| 153 |
Edtop(nl,1) = Edsf(nl) |
| 154 |
DO k=1,kbot-1 |
| 155 |
Edbot(nl,k) = Edtop(nl,k)*expAddr(k) |
| 156 |
Edtop(nl,k+1) = Edbot(nl,k) |
| 157 |
ENDDO |
| 158 |
Edbot(nl,kbot) = Edtop(nl,kbot)*expAddr(kbot) |
| 159 |
|
| 160 |
C start with Aas solution (no increasing mode) |
| 161 |
Estop(nl,1) = Essf(nl) |
| 162 |
DO k=1,kbot-1 |
| 163 |
c2 = Estop(nl,k) - x(k)*Edtop(nl,k) |
| 164 |
Estop(nl,k+1) = MAX(0., c2*expAsdr(k) + x(k)*Edbot(nl,k)) |
| 165 |
Eubot(nl,k) = MAX(0., R2(k)*c2*expAsdr(k) + y(k)*Edbot(nl,k)) |
| 166 |
Eutop(nl,k) = R2(k)*c2 + y(k)*Edtop(nl,k) |
| 167 |
c1out(nl,k) = 0. |
| 168 |
c2out(nl,k) = c2 |
| 169 |
ENDDO |
| 170 |
C Aas b.c. in bottom layer |
| 171 |
c2 = Estop(nl,kbot) - x(kbot)*Edtop(nl,kbot) |
| 172 |
Eutop(nl,kbot) = R2(kbot)*c2 + y(kbot)*Edtop(nl,kbot) |
| 173 |
c1out(nl,kbot) = 0. |
| 174 |
c2out(nl,kbot) = c2 |
| 175 |
|
| 176 |
c improve solution iteratively |
| 177 |
DO iter=1,niter |
| 178 |
c bottom boundary condition |
| 179 |
Eubot(nl,kbot-1) = Eutop(nl,kbot) |
| 180 |
|
| 181 |
DO k=kbot-1,2,-1 |
| 182 |
c compute Eubot(k-1) from Estop(k) and Eubot(k) |
| 183 |
tmp = Estop(nl,k)-x(k)*Edtop(nl,k) |
| 184 |
c1 = (Eubot(nl,k)-R2(k)*expAsdr(k)*tmp-y(k)*Edbot(nl,k)) |
| 185 |
& *idenom(k) |
| 186 |
c2 = (R1(k)*tmp + y(k)*expmAudr(k)*Edbot(nl,k) |
| 187 |
& - expmAudr(k)*Eubot(nl,k))*idenom(k) |
| 188 |
Eunew = R2(k)*c2 + R1(k)*expmAudr(k)*c1 + y(k)*Edtop(nl,k) |
| 189 |
Eubot(nl,k-1) = MAX(0., Eunew) |
| 190 |
ENDDO |
| 191 |
DO k=1,kbot-1 |
| 192 |
c compute Estop(k+1) from Estop(k) and Eubot(k) |
| 193 |
tmp = Estop(nl,k) - x(k)*Edtop(nl,k) |
| 194 |
c1 = (Eubot(nl,k)-R2(k)*expAsdr(k)*tmp-y(k)*Edbot(nl,k)) |
| 195 |
& *idenom(k) |
| 196 |
c2 = (R1(k)*tmp + y(k)*expmAudr(k)*Edbot(nl,k) |
| 197 |
& - expmAudr(k)*Eubot(nl,k))*idenom(k) |
| 198 |
Esnew = expAsdr(k)*c2 + c1 + x(k)*Edbot(nl,k) |
| 199 |
Estop(nl,k+1) = MAX(0., Esnew) |
| 200 |
Eutop(nl,k) = R2(k)*c2+R1(k)*expmAudr(k)*c1+y(k)*Edtop(nl,k) |
| 201 |
c1out(nl,k) = c1 |
| 202 |
c2out(nl,k) = c2 |
| 203 |
ENDDO |
| 204 |
C Aas b.c. in bottom layer |
| 205 |
c2 = Estop(nl,kbot) - x(kbot)*Edtop(nl,kbot) |
| 206 |
Eutop(nl,kbot) = R2(kbot)*c2 + y(kbot)*Edtop(nl,kbot) |
| 207 |
c1out(nl,kbot) = 0. |
| 208 |
c2out(nl,kbot) = c2 |
| 209 |
C enddo iter |
| 210 |
ENDDO |
| 211 |
|
| 212 |
c compute missing fields |
| 213 |
C uses c2 from previous iteration! |
| 214 |
Esbot(nl,kbot) = c2*expAsdr(kbot) + x(kbot)*Edbot(nl,kbot) |
| 215 |
Eubot(nl,kbot) = R2(kbot)*c2*expAsdr(kbot) |
| 216 |
& + y(kbot)*Edbot(nl,kbot) |
| 217 |
|
| 218 |
C Es is continuous now (unless negative...) |
| 219 |
DO k=1,kbot-1 |
| 220 |
Esbot(nl,k) = Estop(nl,k+1) |
| 221 |
ENDDO |
| 222 |
C endif thresh |
| 223 |
ENDIF |
| 224 |
|
| 225 |
DO k = 1,Nr |
| 226 |
#ifdef DAR_RADTRANS_RMUS_PAR |
| 227 |
Etopwq = (Edtop(nl,k)+Estop(nl,k)+Eutop(nl,k))*WtouEins(nl) |
| 228 |
Ebotwq = (Edbot(nl,k)+Esbot(nl,k)+Eubot(nl,k))*WtouEins(nl) |
| 229 |
C interpolate and convert to scalar using rmus only!? |
| 230 |
tirrwq(nl,k) = sqrt(Etopwq*Ebotwq)*rmus |
| 231 |
#else |
| 232 |
C convert to scalar irradiance in quanta |
| 233 |
Etopwq = (rmud*Edtop(nl,k)+rmus*Estop(nl,k)+rmuu*Eutop(nl,k)) |
| 234 |
& *WtouEins(nl) |
| 235 |
Ebotwq = (rmud*Edbot(nl,k)+rmus*Esbot(nl,k)+rmuu*Eubot(nl,k)) |
| 236 |
& *WtouEins(nl) |
| 237 |
C and interpolate |
| 238 |
tirrwq(nl,k) = sqrt(Etopwq*Ebotwq) |
| 239 |
#endif |
| 240 |
ENDDO |
| 241 |
|
| 242 |
C enddo nl |
| 243 |
ENDDO |
| 244 |
|
| 245 |
DO k = 1,Nr |
| 246 |
C sum PAR range |
| 247 |
tirrq(k) = 0.0 |
| 248 |
DO nl = darwin_PAR_ilamLo,darwin_PAR_ilamHi |
| 249 |
tirrq(k) = tirrq(k) + tirrwq(nl,k) |
| 250 |
ENDDO |
| 251 |
ENDDO |
| 252 |
c |
| 253 |
#endif /* DAR_RADTRANS */ |
| 254 |
|
| 255 |
return |
| 256 |
end |
| 257 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph