| 1 |
jahn |
1.1 |
C $Header$ |
| 2 |
|
|
C $Name$ |
| 3 |
|
|
|
| 4 |
|
|
#include "RADTRANS_OPTIONS.h" |
| 5 |
|
|
|
| 6 |
|
|
CBOP |
| 7 |
|
|
C !ROUTINE: RADTRANS_RADMOD_DECR |
| 8 |
|
|
|
| 9 |
|
|
C !INTERFACE: ====================================================== |
| 10 |
|
|
subroutine radtrans_radmod_decr(zd,Edtop,Estop, |
| 11 |
|
|
I rmud,rmus,rmuu,a,bt,bb,Dmax, |
| 12 |
|
|
O Edz,Esz,Euz,Eutop) |
| 13 |
|
|
|
| 14 |
|
|
C !DESCRIPTION: |
| 15 |
|
|
|
| 16 |
|
|
C !USES: =========================================================== |
| 17 |
|
|
IMPLICIT NONE |
| 18 |
|
|
|
| 19 |
|
|
c SLIGHTLY MODIFIED WG CODE PREVIOUSLY CALLED aasack.f |
| 20 |
|
|
c |
| 21 |
|
|
c Model of irradiance in the water column. Accounts for three |
| 22 |
|
|
c irradiance streams: |
| 23 |
|
|
c |
| 24 |
|
|
c Edz = direct downwelling irradiance |
| 25 |
|
|
c Esz = diffuse downwelling irradiance |
| 26 |
|
|
c Euz = diffuse upwelling irradiance |
| 27 |
|
|
c |
| 28 |
|
|
c Uses Ackelson's (1994, JGR) mod's to the Aas (1987, AO) |
| 29 |
|
|
c two-stream model. |
| 30 |
|
|
c |
| 31 |
|
|
c Propagation is done in energy units, tests are done in quanta, |
| 32 |
|
|
c final is quanta for phytoplankton growth. |
| 33 |
|
|
c |
| 34 |
|
|
c In the spirit of Aas, the solution in each layer is truncated to the |
| 35 |
|
|
c 2 downward decreasing modes. Ed and Es are matched between layers, |
| 36 |
|
|
c Eu is discontinuous. |
| 37 |
|
|
c |
| 38 |
|
|
C !INPUT PARAMETERS: =============================================== |
| 39 |
|
|
C Edtop :: direct downward irradiance at top (incl.cos) |
| 40 |
|
|
C Estop :: diffuse downward irradiance at top (incl.cos) |
| 41 |
|
|
C a :: attenuation coefficient |
| 42 |
|
|
C bt :: total scattering coefficient |
| 43 |
|
|
C bb :: backward scattering coefficient (bt = bf + bb) |
| 44 |
|
|
C Dmax :: depth at which Eu = 0 (not used!) |
| 45 |
|
|
_RL zd |
| 46 |
|
|
_RL Edtop,Estop |
| 47 |
|
|
_RL rmud |
| 48 |
|
|
_RL rmus,rmuu |
| 49 |
|
|
_RL a,bt,bb |
| 50 |
|
|
_RL Dmax |
| 51 |
|
|
c INTEGER myThid |
| 52 |
|
|
|
| 53 |
|
|
C !OUTPUT PARAMETERS: ============================================== |
| 54 |
|
|
c Edz :: direct downwelling irradiance (incl.cos) |
| 55 |
|
|
c Esz :: diffuse downwelling irradiance (incl.cos) |
| 56 |
|
|
c Euz :: diffuse upwelling irradiance (incl.cos) |
| 57 |
|
|
c Eutop :: diffuse upwelling irradiance at top of layer |
| 58 |
|
|
_RL Esz,Euz,Edz,Eutop |
| 59 |
|
|
|
| 60 |
|
|
C !LOCAL VARIABLES: ================================================ |
| 61 |
|
|
_RL cd,au,Bu,Cu |
| 62 |
|
|
_RL as,Bs,Cs,Bd,Fd |
| 63 |
|
|
_RL bquad,cquad,sqarg |
| 64 |
|
|
_RL a2,denom,x,y |
| 65 |
|
|
_RL c2,Ta2z |
| 66 |
|
|
c note - have left WG's assignment of these params so as to keep |
| 67 |
|
|
c his code as similar as possible to what he provided, |
| 68 |
|
|
c but could assign elsewhere... |
| 69 |
|
|
|
| 70 |
|
|
_RL rbot, rd, ru |
| 71 |
|
|
data rbot /0.0/ !bottom reflectance |
| 72 |
|
|
data rd /1.5/ !these are taken from Ackleson, et al. 1994 (JGR) |
| 73 |
|
|
data ru /3.0/ |
| 74 |
|
|
c |
| 75 |
|
|
c Downwelling irradiance: Edz, Esz |
| 76 |
|
|
c Compute irradiance components at depth |
| 77 |
|
|
c direct part |
| 78 |
|
|
cd = (a+bt)*rmud |
| 79 |
|
|
Edz = Edtop*exp(-cd*zd) |
| 80 |
|
|
c |
| 81 |
|
|
au = a*rmuu |
| 82 |
|
|
Bu = ru*bb*rmuu |
| 83 |
|
|
Cu = au+Bu |
| 84 |
|
|
as = a*rmus |
| 85 |
|
|
Bs = rd*bb*rmus |
| 86 |
|
|
Cs = as+Bs |
| 87 |
|
|
Bd = bb*rmud |
| 88 |
|
|
Fd = (bt-bb)*rmud |
| 89 |
|
|
bquad = Cs - Cu |
| 90 |
|
|
cquad = Bs*Bu - Cs*Cu |
| 91 |
|
|
sqarg = bquad*bquad - 4.0*cquad |
| 92 |
|
|
c a1 = 0.5*(-bquad + sqrt(sqarg)) |
| 93 |
|
|
a2 = 0.5*(-bquad - sqrt(sqarg)) ! K of Aas |
| 94 |
|
|
denom = (cd-Cs)*(cd+Cu) + Bs*Bu |
| 95 |
|
|
x = -((cd+Cu)*Fd+Bu*Bd)/denom |
| 96 |
|
|
y = (-Bs*Fd+(cd-Cs)*Bd)/denom |
| 97 |
|
|
c2 = Estop - x*Edtop |
| 98 |
|
|
Ta2z = exp(a2*zd) |
| 99 |
|
|
Esz = c2*Ta2z + x*Edz |
| 100 |
|
|
Esz = max(Esz,0.0) |
| 101 |
|
|
Euz = ((a2+Cs)*c2)*Ta2z/Bu + y*Edz |
| 102 |
|
|
Euz = max(Euz,0.0) |
| 103 |
|
|
c |
| 104 |
|
|
c Eu at top of layer |
| 105 |
|
|
Eutop = (a2+Cs)*c2/Bu + y*Edtop |
| 106 |
|
|
Eutop = max(Eutop,0.0) |
| 107 |
|
|
c |
| 108 |
|
|
return |
| 109 |
|
|
end |
Parent Directory
| 
Revision Log
| 
Revision Graph