Line data Source code
1 : module ebert_curry_ice_optics
2 :
3 :
4 : use shr_kind_mod, only: r8 => shr_kind_r8
5 : use physconst, only: gravit
6 : use ppgrid, only: pcols, pver
7 : use physics_types, only: physics_state
8 : use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_get_field, pbuf_old_tim_idx
9 : use constituents, only: cnst_get_ind
10 : use radconstants, only: nswbands, nlwbands, get_sw_spectral_boundaries
11 : use cam_abortutils, only: endrun
12 :
13 : implicit none
14 : private
15 : save
16 :
17 : public :: &
18 : ec_rad_props_init, &
19 : ec_ice_optics_sw, &
20 : ec_ice_get_rad_props_lw
21 :
22 :
23 : real(r8), public, parameter:: scalefactor = 1._r8 !500._r8/917._r8
24 :
25 : ! indices into pbuf
26 : integer :: iciwp_idx = 0
27 : integer :: iclwp_idx = 0
28 : integer :: cld_idx = 0
29 : integer :: rei_idx = 0
30 :
31 : ! indices into constituents for old optics
32 : integer :: ixcldice ! cloud ice water index
33 : integer :: ixcldliq ! cloud liquid water index
34 :
35 :
36 : !==============================================================================
37 : contains
38 : !==============================================================================
39 :
40 1536 : subroutine ec_rad_props_init()
41 :
42 : integer :: err
43 :
44 1536 : iciwp_idx = pbuf_get_index('ICIWP',errcode=err)
45 1536 : iclwp_idx = pbuf_get_index('ICLWP',errcode=err)
46 1536 : cld_idx = pbuf_get_index('CLD')
47 1536 : rei_idx = pbuf_get_index('REI')
48 :
49 : ! old optics
50 1536 : call cnst_get_ind('CLDICE', ixcldice)
51 1536 : call cnst_get_ind('CLDLIQ', ixcldliq)
52 :
53 1536 : end subroutine ec_rad_props_init
54 :
55 : !==============================================================================
56 :
57 0 : subroutine ec_ice_optics_sw (state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f, oldicewp)
58 :
59 : type(physics_state), intent(in) :: state
60 : type(physics_buffer_desc), pointer :: pbuf(:)
61 :
62 : real(r8),intent(out) :: ice_tau (nswbands,pcols,pver) ! extinction optical depth
63 : real(r8),intent(out) :: ice_tau_w (nswbands,pcols,pver) ! single scattering albedo * tau
64 : real(r8),intent(out) :: ice_tau_w_g(nswbands,pcols,pver) ! asymmetry parameter * tau * w
65 : real(r8),intent(out) :: ice_tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w
66 : logical, intent(in) :: oldicewp
67 :
68 0 : real(r8), pointer, dimension(:,:) :: rei
69 0 : real(r8), pointer, dimension(:,:) :: cldn
70 0 : real(r8), pointer, dimension(:,:) :: tmpptr
71 : real(r8), dimension(pcols,pver) :: cicewp
72 : real(r8), dimension(nswbands) :: wavmin
73 : real(r8), dimension(nswbands) :: wavmax
74 : !
75 : ! ice water coefficients (Ebert and Curry,1992, JGR, 97, 3831-3836)
76 : real(r8) :: abari(4) = & ! a coefficient for extinction optical depth
77 : (/ 3.448e-03_r8, 3.448e-03_r8,3.448e-03_r8,3.448e-03_r8/)
78 : real(r8) :: bbari(4) = & ! b coefficient for extinction optical depth
79 : (/ 2.431_r8 , 2.431_r8 ,2.431_r8 ,2.431_r8 /)
80 : real(r8) :: cbari(4) = & ! c coefficient for single scat albedo
81 : (/ 1.00e-05_r8 , 1.10e-04_r8 ,1.861e-02_r8,.46658_r8 /)
82 : real(r8) :: dbari(4) = & ! d coefficient for single scat albedo
83 : (/ 0.0_r8 , 1.405e-05_r8,8.328e-04_r8,2.05e-05_r8 /)
84 : real(r8) :: ebari(4) = & ! e coefficient for asymmetry parameter
85 : (/ 0.7661_r8 , 0.7730_r8 ,0.794_r8 ,0.9595_r8 /)
86 : real(r8) :: fbari(4) = & ! f coefficient for asymmetry parameter
87 : (/ 5.851e-04_r8, 5.665e-04_r8,7.267e-04_r8,1.076e-04_r8/)
88 :
89 : real(r8) :: abarii ! A coefficient for current spectral band
90 : real(r8) :: bbarii ! B coefficient for current spectral band
91 : real(r8) :: cbarii ! C coefficient for current spectral band
92 : real(r8) :: dbarii ! D coefficient for current spectral band
93 : real(r8) :: ebarii ! E coefficient for current spectral band
94 : real(r8) :: fbarii ! F coefficient for current spectral band
95 :
96 : ! Minimum cloud amount (as a fraction of the grid-box area) to
97 : ! distinguish from clear sky
98 : real(r8), parameter :: cldmin = 1.0e-80_r8
99 :
100 : ! Decimal precision of cloud amount (0 -> preserve full resolution;
101 : ! 10^-n -> preserve n digits of cloud amount)
102 : real(r8), parameter :: cldeps = 0.0_r8
103 :
104 : integer :: ns, i, k, indxsl, lchnk, Nday
105 : integer :: itim_old
106 : real(r8) :: tmp1i, tmp2i, tmp3i, g
107 :
108 0 : Nday = state%ncol
109 0 : lchnk = state%lchnk
110 :
111 0 : itim_old = pbuf_old_tim_idx()
112 0 : call pbuf_get_field(pbuf, cld_idx,cldn, start=(/1,1,itim_old/), kount=(/pcols,pver,1/))
113 0 : call pbuf_get_field(pbuf, rei_idx,rei)
114 :
115 0 : if(oldicewp) then
116 0 : do k=1,pver
117 0 : do i = 1,Nday
118 0 : cicewp(i,k) = 1000.0_r8*state%q(i,k,ixcldice)*state%pdel(i,k) /(gravit* max(0.01_r8,cldn(i,k)))
119 : end do
120 : end do
121 : else
122 0 : if (iciwp_idx<=0) then
123 0 : call endrun('ec_ice_optics_sw: oldicewp must be set to true since ICIWP was not found in pbuf')
124 : endif
125 0 : call pbuf_get_field(pbuf, iciwp_idx, tmpptr)
126 0 : cicewp(1:pcols,1:pver) = 1000.0_r8*tmpptr(1:pcols,1:pver)
127 : endif
128 :
129 0 : call get_sw_spectral_boundaries(wavmin,wavmax,'microns')
130 :
131 0 : do ns = 1, nswbands
132 :
133 0 : if(wavmax(ns) <= 0.7_r8) then
134 : indxsl = 1
135 0 : else if(wavmax(ns) <= 1.25_r8) then
136 : indxsl = 2
137 0 : else if(wavmax(ns) <= 2.38_r8) then
138 : indxsl = 3
139 0 : else if(wavmax(ns) > 2.38_r8) then
140 : indxsl = 4
141 : end if
142 :
143 0 : abarii = abari(indxsl)
144 0 : bbarii = bbari(indxsl)
145 0 : cbarii = cbari(indxsl)
146 0 : dbarii = dbari(indxsl)
147 0 : ebarii = ebari(indxsl)
148 0 : fbarii = fbari(indxsl)
149 :
150 0 : do k=1,pver
151 0 : do i=1,Nday
152 :
153 : ! note that optical properties for ice valid only
154 : ! in range of 13 > rei > 130 micron (Ebert and Curry 92)
155 0 : if (cldn(i,k) >= cldmin .and. cldn(i,k) >= cldeps) then
156 0 : tmp1i = abarii + bbarii/max(13._r8,min(scalefactor*rei(i,k),130._r8))
157 0 : ice_tau(ns,i,k) = cicewp(i,k)*tmp1i
158 : else
159 0 : ice_tau(ns,i,k) = 0.0_r8
160 : endif
161 :
162 0 : tmp2i = 1._r8 - cbarii - dbarii*min(max(13._r8,scalefactor*rei(i,k)),130._r8)
163 0 : tmp3i = fbarii*min(max(13._r8,scalefactor*rei(i,k)),130._r8)
164 : ! Do not let single scatter albedo be 1. Delta-eddington solution
165 : ! for non-conservative case has different analytic form from solution
166 : ! for conservative case, and raddedmx is written for non-conservative case.
167 0 : ice_tau_w(ns,i,k) = ice_tau(ns,i,k) * min(tmp2i,.999999_r8)
168 0 : g = ebarii + tmp3i
169 0 : ice_tau_w_g(ns,i,k) = ice_tau_w(ns,i,k) * g
170 0 : ice_tau_w_f(ns,i,k) = ice_tau_w(ns,i,k) * g * g
171 :
172 : end do ! End do i=1,Nday
173 : end do ! End do k=1,pver
174 : end do ! nswbands
175 :
176 0 : end subroutine ec_ice_optics_sw
177 :
178 : !==============================================================================
179 :
180 0 : subroutine ec_ice_get_rad_props_lw(state, pbuf, abs_od, oldicewp)
181 :
182 : type(physics_state), intent(in) :: state
183 : type(physics_buffer_desc),pointer :: pbuf(:)
184 : real(r8), intent(out) :: abs_od(nlwbands,pcols,pver)
185 : logical, intent(in) :: oldicewp
186 :
187 : real(r8) :: gicewp(pcols,pver)
188 : real(r8) :: gliqwp(pcols,pver)
189 : real(r8) :: cicewp(pcols,pver)
190 : real(r8) :: cliqwp(pcols,pver)
191 : real(r8) :: ficemr(pcols,pver)
192 : real(r8) :: cwp(pcols,pver)
193 : real(r8) :: cldtau(pcols,pver)
194 :
195 0 : real(r8), pointer, dimension(:,:) :: cldn
196 0 : real(r8), pointer, dimension(:,:) :: rei
197 : integer :: ncol, itim_old, lwband, i, k, lchnk
198 :
199 : real(r8) :: kabs, kabsi
200 :
201 : real(r8) kabsl ! longwave liquid absorption coeff (m**2/g)
202 : parameter (kabsl = 0.090361_r8)
203 :
204 0 : real(r8), pointer, dimension(:,:) :: iclwpth, iciwpth
205 :
206 :
207 0 : ncol = state%ncol
208 0 : lchnk = state%lchnk
209 :
210 0 : itim_old = pbuf_old_tim_idx()
211 0 : call pbuf_get_field(pbuf, rei_idx, rei)
212 0 : call pbuf_get_field(pbuf, cld_idx, cldn, start=(/1,1,itim_old/), kount=(/pcols,pver,1/))
213 :
214 :
215 0 : if(oldicewp) then
216 0 : do k=1,pver
217 0 : do i = 1,ncol
218 0 : gicewp(i,k) = state%q(i,k,ixcldice)*state%pdel(i,k)/gravit*1000.0_r8 ! Grid box ice water path.
219 0 : gliqwp(i,k) = state%q(i,k,ixcldliq)*state%pdel(i,k)/gravit*1000.0_r8 ! Grid box liquid water path.
220 0 : cicewp(i,k) = gicewp(i,k) / max(0.01_r8,cldn(i,k)) ! In-cloud ice water path.
221 0 : cliqwp(i,k) = gliqwp(i,k) / max(0.01_r8,cldn(i,k)) ! In-cloud liquid water path.
222 : ficemr(i,k) = state%q(i,k,ixcldice) / &
223 0 : max(1.e-10_r8,(state%q(i,k,ixcldice)+state%q(i,k,ixcldliq)))
224 : end do
225 : end do
226 0 : cwp(:ncol,:pver) = cicewp(:ncol,:pver) + cliqwp(:ncol,:pver)
227 : else
228 0 : if (iclwp_idx<=0 .or. iciwp_idx<=0) then
229 0 : call endrun('ec_ice_get_rad_props_lw: oldicewp must be set to true since ICIWP and/or ICLWP were not found in pbuf')
230 : endif
231 0 : call pbuf_get_field(pbuf, iclwp_idx, iclwpth)
232 0 : call pbuf_get_field(pbuf, iciwp_idx, iciwpth)
233 0 : do k=1,pver
234 0 : do i = 1,ncol
235 0 : cwp(i,k) = 1000.0_r8 *iciwpth(i,k) + 1000.0_r8 *iclwpth(i,k)
236 0 : ficemr(i,k) = 1000.0_r8*iciwpth(i,k)/(max(1.e-18_r8,cwp(i,k)))
237 : end do
238 : end do
239 : endif
240 :
241 0 : do k=1,pver
242 0 : do i=1,ncol
243 :
244 : ! Note from Andrew Conley:
245 : ! Optics for RK no longer supported, This is constructed to get
246 : ! close to bit for bit. Otherwise we could simply use ice water path
247 : !note that optical properties for ice valid only
248 : !in range of 13 > rei > 130 micron (Ebert and Curry 92)
249 0 : kabsi = 0.005_r8 + 1._r8/min(max(13._r8,scalefactor*rei(i,k)),130._r8)
250 0 : kabs = kabsi*ficemr(i,k) ! kabsl*(1._r8-ficemr(i,k)) + kabsi*ficemr(i,k)
251 : !emis(i,k) = 1._r8 - exp(-1.66_r8*kabs*clwp(i,k))
252 0 : cldtau(i,k) = kabs*cwp(i,k)
253 : end do
254 : end do
255 :
256 0 : do lwband = 1,nlwbands
257 0 : abs_od(lwband,1:ncol,1:pver)=cldtau(1:ncol,1:pver)
258 : enddo
259 :
260 0 : end subroutine ec_ice_get_rad_props_lw
261 :
262 : !==============================================================================
263 :
264 : end module ebert_curry_ice_optics
|
