pkg/radtrans/radtrans_radmod.F

C $Header$
C $Name$

#include "RADTRANS_OPTIONS.h"

CBOP
C     !ROUTINE: RADTRANS_RADMOD

C     !INTERFACE: ======================================================
      subroutine radtrans_radmod(zd,Edtop,Estop,
     I                           rmud,rmus,rmuu,a,bt,bb,Dmax,
     O                           Edz,Esz,Euz,Eutop)

C     !DESCRIPTION:

C     !USES: ===========================================================
      IMPLICIT NONE

c SLIGHTLY MODIFIED WG CODE PREVIOUSLY CALLED aasack.f
c
c  Model of irradiance in the water column.  Accounts for three 
c  irradiance streams:
c
c  Edz = direct downwelling irradiance
c  Esz = diffuse downwelling irradiance
c  Euz = diffuse upwelling irradiance
c
c  Uses Ackelson's (1994, JGR) mod's to the Aas (1987, AO) 
c  two-stream model.
c
c  Propagation is done in energy units, tests are done in quanta,
c  final is quanta for phytoplankton growth.
c
c  Commented out terms produce a max error of 
c  0.8% in Esz for a > 0.004 and bb > 0.0001 and
c  3.9% in Euz for a > 0.004 and bb > 0.00063
c
C     !INPUT PARAMETERS: ===============================================
C     Edtop        :: direct downward irradiance at top (incl.cos)
C     Estop        :: diffuse downward irradiance at top (incl.cos)
C     a            :: attenuation coefficient
C     bt           :: total scattering coefficient
C     bb           :: backward scattering coefficient (bt = bf + bb)
C     Dmax         :: depth at which Eu = 0 (not used!)
      _RL zd
      _RL Edtop,Estop
      _RL rmud
      _RL rmus,rmuu
      _RL a,bt,bb
      _RL Dmax
c     INTEGER myThid

C     !OUTPUT PARAMETERS: ==============================================
c     Edz          :: direct downwelling irradiance (incl.cos)
c     Esz          :: diffuse downwelling irradiance (incl.cos)
c     Euz          :: diffuse upwelling irradiance (incl.cos)
      _RL Esz,Euz,Edz,Eutop

C     !LOCAL VARIABLES: ================================================
      _RL cd,au,Bu,Cu
      _RL as,Bs,Cs,Bd,Fd
      _RL bquad,cquad,sqarg
      _RL a1,a2,S,SEdz,a2ma1
      _RL rM,rN
      _RL c2,Ta2z,Eutmp
c note - have left WG's assignment of these params so as to keep 
c his code as similar as possible to what he provided, 
c but could assign elsewhere...

      _RL rbot, rd, ru
      data rbot /0.0/ !bottom reflectance
      data rd /1.5/   !these are taken from Ackleson, et al. 1994 (JGR)
      data ru /3.0/
c
c  Downwelling irradiance: Edz, Esz
c  Compute irradiance components at depth
c     direct part
      cd = (a+bt)*rmud
      Edz = Edtop*exp(-cd*zd)
c
      au = a*rmuu
      Bu = ru*bb*rmuu
      Cu = au+Bu
      as = a*rmus
      Bs = rd*bb*rmus
      Cs = as+Bs
      Bd = bb*rmud
      Fd = (bt-bb)*rmud
      bquad = Cs - Cu
      cquad = Bs*Bu - Cs*Cu
      sqarg = bquad*bquad - 4.0*cquad
      a1 = 0.5*(-bquad + sqrt(sqarg))  ! K of Aas
      a2 = 0.5*(-bquad - sqrt(sqarg))
      S = -(Bu*Bd + Cu*Fd)
      SEdz = S*Edz
      a2ma1 = a2 - a1
      rM = SEdz/(a1*a2ma1)
      rN = SEdz/(a2*a2ma1)
c      ea2Dmax = exp(a2ma1*Dmax)
c      c1 = (rN-rM)*exp(-a1*Dmax) - (Estop-rM+rN)*ea2Dmax
c     *                             /(1.0-ea2Dmax)
c      c2 = Estop - rM + rN - c1
      c2 = Estop - rM + rN 
c      a1arg = a1*zd
c      a1arg = min(a1arg,82.8)
c      Ta1z = exp(a1arg)
      Ta2z = exp(a2*zd)
c      Esz = c1*Ta1z + c2*Ta2z + rM - rN
      Esz = c2*Ta2z + rM - rN
      Esz = max(Esz,0.0)
c      Eutmp = ((a1+Cs)*c1)*Ta1z + ((a2+Cs)*c2)*Ta2z + Cs*rM
c     *             - Cs*rN - Fd*Edz
      Eutmp = ((a2+Cs)*c2)*Ta2z + Cs*rM
     *             - Cs*rN - Fd*Edz
      Euz = Eutmp/Bu
      Euz = max(Euz,0.0)
c
c Eu at top of layer
      rM = S*Edtop/(a1*a2ma1)
      rN = S*Edtop/(a2*a2ma1)
      Eutmp = (a2+Cs)*c2 + Cs*rM - Cs*rN - Fd*Edtop
      Eutop = Eutmp/Bu
      Eutop = max(Eutop,0.0)
c
      return
      end
1	jahn	1.1	C $Header$
2			C $Name$
3
4			#include "RADTRANS_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: RADTRANS_RADMOD
8
9			C !INTERFACE: ======================================================
10			subroutine radtrans_radmod(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 Commented out terms produce a max error of
35			c 0.8% in Esz for a > 0.004 and bb > 0.0001 and
36			c 3.9% in Euz for a > 0.004 and bb > 0.00063
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			_RL Esz,Euz,Edz,Eutop
58
59			C !LOCAL VARIABLES: ================================================
60			_RL cd,au,Bu,Cu
61			_RL as,Bs,Cs,Bd,Fd
62			_RL bquad,cquad,sqarg
63			_RL a1,a2,S,SEdz,a2ma1
64			_RL rM,rN
65			_RL c2,Ta2z,Eutmp
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 = Edtopexp(-cdzd)
80			c
81			au = a*rmuu
82			Bu = rubbrmuu
83			Cu = au+Bu
84			as = a*rmus
85			Bs = rdbbrmus
86			Cs = as+Bs
87			Bd = bb*rmud
88			Fd = (bt-bb)*rmud
89			bquad = Cs - Cu
90			cquad = BsBu - CsCu
91			sqarg = bquadbquad - 4.0cquad
92			a1 = 0.5*(-bquad + sqrt(sqarg)) ! K of Aas
93			a2 = 0.5*(-bquad - sqrt(sqarg))
94			S = -(BuBd + CuFd)
95			SEdz = S*Edz
96			a2ma1 = a2 - a1
97			rM = SEdz/(a1*a2ma1)
98			rN = SEdz/(a2*a2ma1)
99			c ea2Dmax = exp(a2ma1*Dmax)
100			c c1 = (rN-rM)exp(-a1Dmax) - (Estop-rM+rN)*ea2Dmax
101			c * /(1.0-ea2Dmax)
102			c c2 = Estop - rM + rN - c1
103			c2 = Estop - rM + rN
104			c a1arg = a1*zd
105			c a1arg = min(a1arg,82.8)
106			c Ta1z = exp(a1arg)
107			Ta2z = exp(a2*zd)
108			c Esz = c1Ta1z + c2Ta2z + rM - rN
109			Esz = c2*Ta2z + rM - rN
110			Esz = max(Esz,0.0)
111			c Eutmp = ((a1+Cs)c1)Ta1z + ((a2+Cs)c2)Ta2z + Cs*rM
112			c * - CsrN - FdEdz
113			Eutmp = ((a2+Cs)c2)Ta2z + Cs*rM
114			* - CsrN - FdEdz
115			Euz = Eutmp/Bu
116			Euz = max(Euz,0.0)
117			c
118			c Eu at top of layer
119			rM = SEdtop/(a1a2ma1)
120			rN = SEdtop/(a2a2ma1)
121			Eutmp = (a2+Cs)c2 + CsrM - CsrN - FdEdtop
122			Eutop = Eutmp/Bu
123			Eutop = max(Eutop,0.0)
124			c
125			return
126			end