| 1 |
jahn |
1.3 |
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 |
jahn |
1.2 |
C $Name: $ |
| 3 |
jahn |
1.1 |
|
| 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 |
jahn |
1.2 |
O c1out, c2out, |
| 15 |
jahn |
1.1 |
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 |
jahn |
1.2 |
_RL c1out(tlam,Nr), c2out(tlam,Nr) |
| 79 |
jahn |
1.1 |
|
| 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 |
jahn |
1.3 |
kbot = MIN(kmax, Nr) |
| 108 |
jahn |
1.1 |
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 |
jahn |
1.2 |
c1out(nl,k) = 0.0 |
| 121 |
|
|
c2out(nl,k) = 0.0 |
| 122 |
jahn |
1.1 |
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 |
jahn |
1.2 |
c1out(nl,k) = 0. |
| 168 |
|
|
c2out(nl,k) = c2 |
| 169 |
jahn |
1.1 |
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 |
jahn |
1.2 |
c1out(nl,kbot) = 0. |
| 174 |
|
|
c2out(nl,kbot) = c2 |
| 175 |
jahn |
1.1 |
|
| 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 |
jahn |
1.2 |
c1out(nl,k) = c1 |
| 202 |
|
|
c2out(nl,k) = c2 |
| 203 |
jahn |
1.1 |
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 |
jahn |
1.2 |
c1out(nl,kbot) = 0. |
| 208 |
|
|
c2out(nl,kbot) = c2 |
| 209 |
jahn |
1.1 |
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