/[MITgcm]/MITgcm_contrib/darwin2/pkg/radtrans/radtrans_radmod.F
ViewVC logotype

Contents of /MITgcm_contrib/darwin2/pkg/radtrans/radtrans_radmod.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.1 - (show annotations) (download)
Wed Apr 13 18:56:26 2011 UTC (14 years, 3 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt64k_20130723, ctrb_darwin2_ckpt65w_20160512, ctrb_darwin2_ckpt65j_20150225, ctrb_darwin2_ckpt63l_20120405, ctrb_darwin2_ckpt66g_20170424, ctrb_darwin2_ckpt64h_20130528, ctrb_darwin2_ckpt66k_20171025, ctrb_darwin2_ckpt66n_20180118, ctrb_darwin2_ckpt62v_20110413, ctrb_darwin2_ckpt65v_20160409, ctrb_darwin2_ckpt65s_20160114, ctrb_darwin2_ckpt65_20140718, ctrb_darwin2_ckpt64m_20130820, ctrb_darwin2_ckpt66d_20170214, ctrb_darwin2_ckpt64r_20131210, ctrb_darwin2_ckpt65m_20150615, ctrb_darwin2_ckpt65q_20151118, ctrb_darwin2_ckpt65o_20150914, ctrb_darwin2_ckpt64f_20130405, ctrb_darwin2_ckpt63f_20111201, ctrb_darwin2_ckpt64a_20121116, ctrb_darwin2_ckpt65p_20151023, ctrb_darwin2_ckpt64n_20130826, ctrb_darwin2_ckpt65e_20140929, ctrb_darwin2_ckpt64o_20131024, ctrb_darwin2_ckpt64v_20140411, ctrb_darwin2_ckpt64z_20140711, ctrb_darwin2_ckpt65l_20150504, ctrb_darwin2_ckpt65z_20160929, ctrb_darwin2_ckpt65n_20150729, ctrb_darwin2_ckpt62y_20110526, ctrb_darwin2_ckpt64y_20140622, ctrb_darwin2_ckpt65d_20140915, ctrb_darwin2_ckpt64t_20140202, ctrb_darwin2_ckpt66h_20170602, ctrb_darwin2_ckpt64i_20130622, ctrb_darwin2_ckpt62x_20110513, ctrb_darwin2_ckpt64s_20140105, ctrb_darwin2_ckpt62w_20110426, ctrb_darwin2_ckpt64x_20140524, ctrb_darwin2_ckpt63o_20120629, ctrb_darwin2_ckpt64e_20130305, ctrb_darwin2_ckpt65x_20160612, ctrb_darwin2_ckpt66f_20170407, ctrb_darwin2_ckpt63c_20111011, ctrb_darwin2_ckpt63i_20120124, ctrb_darwin2_ckpt65g_20141120, ctrb_darwin2_ckpt63m_20120506, ctrb_darwin2_ckpt63s_20120908, ctrb_darwin2_ckpt65k_20150402, ctrb_darwin2_ckpt63e_20111107, ctrb_darwin2_ckpt64w_20140502, ctrb_darwin2_ckpt63b_20110830, ctrb_darwin2_ckpt63j_20120217, ctrb_darwin2_ckpt66a_20161020, ctrb_darwin2_ckpt63r_20120817, ctrb_darwin2_ckpt64g_20130503, ctrb_darwin2_ckpt64l_20130806, ctrb_darwin2_ckpt63g_20111220, ctrb_darwin2_ckpt65f_20141014, ctrb_darwin2_ckpt64c_20130120, ctrb_darwin2_ckpt63a_20110804, ctrb_darwin2_ckpt66b_20161219, ctrb_darwin2_ckpt64u_20140308, ctrb_darwin2_ckpt64j_20130704, ctrb_darwin2_ckpt65i_20150123, ctrb_darwin2_ckpt66j_20170815, ctrb_darwin2_ckpt65y_20160801, ctrb_darwin2_ckpt63h_20111230, ctrb_darwin2_ckpt63p_20120707, ctrb_darwin2_ckpt66c_20170121, ctrb_darwin2_ckpt65a_20140728, ctrb_darwin2_ckpt65b_20140812, ctrb_darwin2_ckpt65t_20160221, ctrb_darwin2_ckpt64p_20131118, ctrb_darwin2_ckpt63d_20111107, ctrb_darwin2_ckpt63q_20120731, ctrb_darwin2_ckpt63_20110728, ctrb_darwin2_ckpt64b_20121224, ctrb_darwin2_ckpt64d_20130219, ctrb_darwin2_ckpt66o_20180209, ctrb_darwin2_ckpt66e_20170314, ctrb_darwin2_ckpt64_20121012, ctrb_darwin2_ckpt64q_20131118, ctrb_darwin2_baseline, ctrb_darwin2_ckpt64p_20131024, ctrb_darwin2_ckpt65u_20160315, ctrb_darwin2_ckpt65r_20151221, ctrb_darwin2_ckpt66i_20170718, ctrb_darwin2_ckpt63n_20120604, ctrb_darwin2_ckpt63k_20120317, ctrb_darwin2_ckpt65c_20140830, ctrb_darwin2_ckpt62z_20110622, ctrb_darwin2_ckpt66l_20171025, ctrb_darwin2_ckpt65h_20141217, ctrb_darwin2_ckpt66m_20171213, HEAD 
darwin2 initial checkin
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 = 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 a1 = 0.5*(-bquad + sqrt(sqarg)) ! K of Aas
93 a2 = 0.5*(-bquad - sqrt(sqarg))
94 S = -(Bu*Bd + Cu*Fd)
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(-a1*Dmax) - (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 = c1*Ta1z + c2*Ta2z + 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 * - Cs*rN - Fd*Edz
113 Eutmp = ((a2+Cs)*c2)*Ta2z + Cs*rM
114 * - Cs*rN - Fd*Edz
115 Euz = Eutmp/Bu
116 Euz = max(Euz,0.0)
117 c
118 c Eu at top of layer
119 rM = S*Edtop/(a1*a2ma1)
120 rN = S*Edtop/(a2*a2ma1)
121 Eutmp = (a2+Cs)*c2 + Cs*rM - Cs*rN - Fd*Edtop
122 Eutop = Eutmp/Bu
123 Eutop = max(Eutop,0.0)
124 c
125 return
126 end
  ViewVC Help
Powered by ViewVC 1.1.22