| blob | 0a3c1f6982bd3c7286f63b3bc376e235493d875d |
1 module module_ftuv_driver
3 use module_wave_data, only : nw
5 implicit none
7 save
9 integer, parameter :: dp = selected_real_kind(14,300)
10 !----------------------------------------------------
11 ! private variables
12 !----------------------------------------------------
13 integer, private, parameter :: nref = 51
14 real(dp), private, parameter :: m2s = 60._dp
15 real(dp), private, parameter :: Pa2hPa = .01_dp
17 integer, private :: nz, kcon
18 integer, private :: next, last
19 integer, private :: curjulday = 0
21 !++alma 2012-12-01 modis landuse from sw
22 integer, private, allocatable :: luse2usgs(:)
23 !--
25 real(dp), private :: esfact = 1.0_dp
26 real(dp), private :: dobson
28 real(dp), private :: dels
30 real(dp), private :: zref(nref), tref(nref), airref(nref), o3ref(nref)
31 real(dp), private :: albedo(nw-1,2)
33 type column_density
34 integer :: ncoldens_levs
35 integer :: ndays_of_year
36 real(dp), pointer :: col_levs(:)
37 real(dp), pointer :: day_of_year(:)
38 real(dp), pointer :: o3_col_dens(:,:,:,:)
39 real(dp), pointer :: o2_col_dens(:,:,:,:)
40 logical :: is_allocated
41 end type column_density
43 type(column_density), private, allocatable :: col_dens(:)
44 logical, private :: photo_inti_initialized = .false.
47 data zref/ 00.0, 01.0, 02.0, 03.0, 04.0, 05.0, 06.0, 07.0, 08.0, 09.0, &
48 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, &
49 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0, &
50 30.0, 31.0, 32.0, 33.0, 34.0, 35.0, 36.0, 37.0, 38.0, 39.0, &
51 40.0, 41.0, 42.0, 43.0, 44.0, 45.0, 46.0, 47.0, 48.0, 49.0, &
52 50.0 /
54 data tref/ 288.150, 281.651, 275.154, 268.659, 262.166, 255.676, 249.187, &
55 242.700, 236.215, 229.733, 223.252, 216.774, 216.650, 216.650, &
56 216.650, 216.650, 216.650, 216.650, 216.650, 216.650, 216.650, &
57 217.581, 218.574, 219.567, 220.560, 221.552, 222.544, 223.536, &
58 224.527, 225.518, 226.509, 227.500, 228.490, 230.973, 233.743, &
59 236.513, 239.282, 242.050, 244.818, 247.584, 250.350, 253.114, &
60 255.878, 258.641, 261.403, 264.164, 266.925, 269.684, 270.650, &
61 270.650, 270.650 /
63 data airref/ 2.55E+19, 2.31E+19, 2.09E+19, 1.89E+19, 1.70E+19, 1.53E+19, &
64 1.37E+19, 1.23E+19, 1.09E+19, 9.71E+18, 8.60E+18, 7.59E+18, &
65 6.49E+18, 5.54E+18, 4.74E+18, 4.05E+18, 3.46E+18, 2.96E+18, &
66 2.53E+18, 2.16E+18, 1.85E+18, 1.57E+18, 1.34E+18, 1.14E+18, &
67 9.76E+17, 8.33E+17, 7.12E+17, 6.09E+17, 5.21E+17, 4.47E+17, &
68 3.83E+17, 3.28E+17, 2.82E+17, 2.41E+17, 2.06E+17, 1.76E+17, &
69 1.51E+17, 1.30E+17, 1.12E+17, 9.62E+16, 8.31E+16, 7.19E+16, &
70 6.23E+16, 5.40E+16, 4.70E+16, 4.09E+16, 3.56E+16, 3.11E+16, &
71 2.74E+16, 2.42E+16, 2.14E+16 /
73 data o3ref/ 8.00E+11, 7.39E+11, 6.90E+11, 6.22E+11, 5.80E+11, 5.67E+11, &
74 5.70E+11, 5.86E+11, 6.50E+11, 8.23E+11, 1.13E+12, 1.57E+12, &
75 2.02E+12, 2.24E+12, 2.35E+12, 2.57E+12, 2.95E+12, 3.47E+12, &
76 4.04E+12, 4.49E+12, 4.77E+12, 4.88E+12, 4.86E+12, 4.73E+12, &
77 4.54E+12, 4.32E+12, 4.03E+12, 3.65E+12, 3.24E+12, 2.85E+12, &
78 2.52E+12, 2.26E+12, 2.03E+12, 1.80E+12, 1.58E+12, 1.40E+12, &
79 1.22E+12, 1.04E+12, 8.73E+11, 7.31E+11, 6.07E+11, 4.95E+11, &
80 3.98E+11, 3.18E+11, 2.54E+11, 2.03E+11, 1.62E+11, 1.30E+11, &
81 1.04E+11, 8.27E+10, 6.61E+10/
84 !----------------------------------------------------
85 ! private parameters
86 !----------------------------------------------------
87 integer, private, parameter :: pid_no2 = 4
88 integer, private, parameter :: pid_n2o = 9
89 integer, private, parameter :: pid_n2o5 = 8
90 integer, private, parameter :: pid_o31d = 2
91 integer, private, parameter :: pid_o33p = 3
92 integer, private, parameter :: pid_hono = 12
93 integer, private, parameter :: pid_hno3 = 13
94 integer, private, parameter :: pid_hno4 = 14
95 integer, private, parameter :: pid_no3o2 = 5
96 integer, private, parameter :: pid_no3o = 6
97 integer, private, parameter :: pid_h2o2 = 11
98 integer, private, parameter :: pid_ch2om = 16
99 integer, private, parameter :: pid_ch2or = 15
100 integer, private, parameter :: pid_ald = 17
101 integer, private, parameter :: pid_ch3choa = 17
102 integer, private, parameter :: pid_ch3chob = 18
103 integer, private, parameter :: pid_ch3choc = 19
104 integer, private, parameter :: pid_op1 = 24
105 integer, private, parameter :: pid_op2 = 1 ! unknown
106 integer, private, parameter :: pid_paa = 16
107 integer, private, parameter :: pid_ket = 23
108 integer, private, parameter :: pid_glya = 21
109 integer, private, parameter :: pid_glyb = 21
110 integer, private, parameter :: pid_mgly = 22
111 integer, private, parameter :: pid_dcb = 1 ! unknown
112 integer, private, parameter :: pid_onit = 25
113 integer, private, parameter :: pid_pan = 26
114 integer, private, parameter :: pid_mvk = 27
115 integer, private, parameter :: pid_macr = 28
116 integer, private, parameter :: pid_hyac = 29
117 integer, private, parameter :: pid_glyald= 30
119 ! added for CB05CL
120 ! These numbers are referred to from module_wave_data.F (sjref(ii,XX))
121 integer, private, parameter :: pid_cl2 = 56
122 integer, private, parameter :: pid_hocl = 57
123 integer, private, parameter :: pid_fmcl = 58
125 contains
127 subroutine ftuv_driver( id, curr_secs, dtstep, config_flags, &
128 gmt, julday, &
129 p_phy, t_phy, rho_phy, p8w, t8w, &
130 xlat, xlong, z_at_w, &
131 moist, chem,gd_cloud,gd_cloud2, &
132 ph_no2,ph_o31d,ph_o33p,ph_hono, &
133 ph_hno3,ph_hno4,ph_no3o2,ph_no3o, &
134 ph_h2o2,ph_ch2om,ph_ch2or,ph_ald, &
135 ph_op1,ph_op2,ph_paa,ph_ket,ph_glya, &
136 ph_glyb,ph_mgly,ph_dcb,ph_onit, &
137 ph_macr,ph_ch3coc2h5,ph_n2o, &
138 ph_pan,ph_mpan,ph_acetol,ph_gly, &
139 ph_bigald,ph_mek,ph_c2h5ooh,ph_c3h7ooh, &
140 ph_pooh,ph_rooh,ph_xooh,ph_isopooh, &
141 ph_alkooh,ph_mekooh,ph_tolooh, &
142 ph_terpooh,ph_n2o5,ph_mvk,ph_glyald, &
143 ph_hyac, &
144 ph_cl2,ph_hocl,ph_fmcl, &
145 ivgtyp, &
146 ntuv_lev, ntuv_bin, ntuv_pht, &
147 ph_radfld, ph_adjcoe, ph_prate, &
148 wc_x, zref_x, &
149 tauaer1, tauaer2, tauaer3, tauaer4, & !rajesh
150 waer1, waer2, waer3, waer4, & !rajesh
151 gaer1, gaer2, gaer3, gaer4, & !rajesh
152 ids,ide, jds,jde, kds,kde, &
153 ims,ime, jms,jme, kms,kme, &
154 its,ite, jts,jte, kts,kte )
156 use module_configure
157 use module_state_description
158 use module_model_constants
159 use module_data_sorgam, only : mw_so4_aer
161 use module_wave_data, only : tuv_jmax, nw, wc
162 use module_ftuv_subs, only : calc_zenith
164 implicit none
166 !----------------------------------------------------
167 ! input arguments
168 !----------------------------------------------------
169 integer, intent(in ) :: id, julday, &
170 ids,ide, jds,jde, kds,kde, &
171 ims,ime, jms,jme, kms,kme, &
172 its,ite, jts,jte, kts,kte
173 integer, intent(in ) :: ntuv_lev, ntuv_bin, ntuv_pht
174 integer, intent(in) :: ivgtyp(ims:ime,jms:jme)
176 real, intent(in ) :: dtstep, gmt
177 real(dp), intent(in ) :: curr_secs
179 type(grid_config_rec_type), intent(in ) :: config_flags
181 real, intent(in ) :: p_phy(ims:ime,kms:kme,jms:jme), &
182 t_phy(ims:ime,kms:kme,jms:jme), &
183 rho_phy(ims:ime,kms:kme,jms:jme), &
184 p8w(ims:ime,kms:kme,jms:jme), &
185 t8w(ims:ime,kms:kme,jms:jme), &
186 z_at_w(ims:ime,kms:kme,jms:jme)
188 !----------------------------------------------------
189 ! these arrays are for cloudwater/ice coming from
190 ! grell convection, optional
191 !----------------------------------------------------
192 real, optional, &
193 intent(in ) :: gd_cloud(ims:ime,kms:kme,jms:jme), &
194 gd_cloud2(ims:ime,kms:kme,jms:jme)
195 real, intent(in ) :: xlat(ims:ime,jms:jme), &
196 xlong(ims:ime,jms:jme)
197 real, intent(in ) :: moist(ims:ime,kms:kme,jms:jme,num_moist)
198 real, intent(in ) :: chem(ims:ime,kms:kme,jms:jme,num_chem)
199 !rajesh: add arrays for aerosol optical properties
200 real, dimension( ims:ime, kms:kme, jms:jme ), &
201 intent(in ) :: tauaer1, tauaer2, tauaer3, tauaer4, &
202 waer1, waer2, waer3, waer4, &
203 gaer1, gaer2, gaer3, gaer4
205 !----------------------------------------------------
206 ! output arguments
207 !----------------------------------------------------
208 real, dimension( ims:ime, kms:kme, jms:jme ), &
209 intent(out ) :: ph_no2, ph_o31d, ph_o33p, ph_hono, &
210 ph_hno3, ph_hno4, ph_no3o2, ph_no3o, &
211 ph_h2o2, ph_ch2om, ph_ch2or, ph_ald, &
212 ph_op1, ph_op2, ph_paa, ph_ket, ph_glya, &
213 ph_glyb, ph_mgly, ph_dcb, &
214 ph_onit, ph_macr, ph_ch3coc2h5, ph_n2o, &
215 ph_pan, ph_mpan, ph_acetol, ph_gly, &
216 ph_bigald, ph_mek, ph_c2h5ooh, ph_c3h7ooh,&
217 ph_pooh,ph_rooh,ph_xooh,ph_isopooh, &
218 ph_alkooh,ph_mekooh,ph_tolooh,ph_terpooh, &
219 ph_n2o5,ph_mvk,ph_glyald,ph_hyac, &
220 ph_cl2,ph_hocl,ph_fmcl
222 real, dimension( ims:ime, ntuv_lev, jms:jme, ntuv_bin ), &
223 intent(out ) :: ph_radfld
224 real, dimension( ims:ime, ntuv_lev, jms:jme, ntuv_pht ), &
225 intent(out ) :: ph_adjcoe
226 real, dimension( ims:ime, ntuv_lev, jms:jme, ntuv_pht ), &
227 intent(out ) :: ph_prate
228 real, dimension(ntuv_bin), &
229 intent(out ) :: wc_x
230 real, dimension(ntuv_lev), &
231 intent(out ) :: zref_x
233 !----------------------------------------------------
234 ! local parameters
235 !----------------------------------------------------
236 real(dp), parameter :: xair_kg = 1.0e3_dp/28.97_dp*6.023e23_dp*1.0e-6_dp
237 real(dp), parameter :: ph_no3o_factor = 1.1236_dp
238 real(dp), parameter :: kg_per_amu = 1.65979e-27_dp
240 !----------------------------------------------------
241 ! local scalars
242 !----------------------------------------------------
243 integer :: i, j, k, kp1, m
244 integer :: astat, istat
245 integer :: isorg
246 integer :: lu
247 real(dp) :: xtime, xhour, xmin, gmtp
248 real(dp) :: wrk, qs, es
249 real(dp) :: alat, along, azim, zen
250 real(dp) :: o3top
251 real(dp) :: o2top
252 real(dp) :: atm_mass_den
253 logical :: chm_is_mozart
255 !----------------------------------------------------
256 ! local arrays
257 !----------------------------------------------------
258 real(dp) :: zen_angle(its:ite,jts:jte)
259 real(dp) :: temp(kts:kte+1), o3(kts:kte+1), air(kts:kte+1)
260 real(dp) :: zh(kts:kte+1), rh(kts:kte+1), xlwc(kts:kte+1)
262 real(dp) :: acb1(kts:kte+1), acb2(kts:kte+1), aoc1(kts:kte+1), aoc2(kts:kte+1)
263 real(dp) :: aant(kts:kte+1), aso4(kts:kte+1), asal(kts:kte+1)
264 ! rajesh: add arrays
265 real(dp) :: tauaer300(kts:kte), tauaer400(kts:kte+1), &
266 tauaer600(kts:kte), tauaer999(kts:kte+1)
267 real(dp) :: waer300(kts:kte), waer400(kts:kte), &
268 waer600(kts:kte), waer999(kts:kte)
269 real(dp) :: gaer300(kts:kte), gaer400(kts:kte), &
270 gaer600(kts:kte), gaer999(kts:kte)
272 real(dp) :: p_jtop(its:ite,jts:jte)
273 real(dp) :: o2_exo_col(its:ite,jts:jte)
274 real(dp) :: o3_exo_col(its:ite,jts:jte)
275 real(dp) :: prate(kts:kte+1,tuv_jmax)
276 real(dp) :: radfld(nz,nw-1)
277 real(dp) :: adjcoe(nz,tuv_jmax)
278 real(dp) :: prate0(nz,tuv_jmax)
280 chm_is_mozart = config_flags%chem_opt == MOZART_KPP .or. &
281 config_flags%chem_opt == MOZCART_KPP .or. &
282 config_flags%chem_opt == T1_MOZCART_KPP .or. &
283 config_flags%chem_opt == MOZART_MOSAIC_4BIN_KPP .or. &
284 config_flags%chem_opt == MOZART_MOSAIC_4BIN_AQ_KPP
285 !----------------------------------------------------
286 ! set o3,o2 column densities at model top
287 !----------------------------------------------------
288 if( chm_is_mozart ) then
289 p_jtop(its:ite,jts:jte) = Pa2hPa * p_phy(its:ite,kte,jts:jte)
290 call p_interp( o2_exo_col, o3_exo_col, p_jtop, &
291 id, its, ite, jts, jte )
292 else
293 o2top = 3.60906e+21_dp
294 o3top = 2.97450e+16_dp
295 endif
297 isorg=0
298 aer_select: SELECT CASE(config_flags%chem_opt)
299 CASE (RADM2SORG,RADM2SORG_KPP,RACMSORG_KPP,RADM2SORG_AQCHEM,RACM_ESRLSORG_KPP,RACMSORG_AQ,RACMSORG_AQCHEM_KPP, &
300 RACM_ESRLSORG_AQCHEM_KPP,CBMZSORG,CBMZSORG_AQ, &
301 CB05_SORG_AQ_KPP,CB05_SORG_VBS_AQ_KPP)
302 isorg=1
303 CALL wrf_debug(15,'SORGAM aerosols initialization ')
304 CASE (MOZART_MOSAIC_4BIN_KPP, MOZART_MOSAIC_4BIN_AQ_KPP)
305 CALL wrf_debug(15,'MOZART_MOSAIC_4BIN_KPP, MOZART_MOSAIC_4BIN_AQ_KPP aerosols initialization ')
306 CASE DEFAULT
307 CALL wrf_debug(15,'no aerosols initialization yet')
308 CALL wrf_message('no aerosols initialization yet')
309 END SELECT aer_select
311 !----------------------------------------------------
312 ! calculate the GMT
313 !----------------------------------------------------
314 xtime = curr_secs/60._dp
315 xhour = real( int(gmt + .01_dp) + int(xtime/m2s),kind=dp )
316 xmin = m2s*gmt + ( xtime - xhour*m2s )
317 gmtp = mod( xhour, 24.0_dp )
318 gmtp = gmtp + xmin/m2s
320 xlwc(:) = 0._dp !initialize in case no cloud water is present
322 !----------------------------------------------------
323 ! set photolysis rates
324 !----------------------------------------------------
325 J_TILE_LOOP : do j = jts, jte
326 I_TILE_LOOP : do i = its, ite
328 !----------------------------------------------------
329 ! set land use type
330 !----------------------------------------------------
331 if( chm_is_mozart ) then
332 !++alma 2012-12-01 modis landuse from sw
333 ! if( ivgtyp(i,j) /= 16 ) then
334 if( luse2usgs( ivgtyp(i,j) ) /= 16 ) then
335 !--
337 lu = 1
338 else
339 lu = 2
340 endif
341 else
342 lu = 1
343 endif
345 if( chm_is_mozart ) then
346 o3top = o3_exo_col(i,j)
347 o2top = o2_exo_col(i,j)
348 endif
350 level_loop : &
351 do k = kts, kte
352 kp1 = k + 1
353 temp(kp1) = t_phy(i,k,j) ! temperature(K)
354 air(kp1) = xair_kg*rho_phy(i,k,j) ! air num.(molecules/cm^3)
355 o3(kp1) = chem(i,k,j,p_o3)*1.0e-6_dp ! ozone(VMR)
357 wrk = t_phy(i,k,j) - 273._dp
358 es = 6.11_dp*10.0_dp**(7.63_dp*wrk/(241.9_dp + wrk)) ! Magnus EQ
359 qs = 0.622_dp*es/(p_phy(i,k,j)*0.01_dp) ! sat mass mix (H2O)
360 wrk = moist(i,k,j,p_qv)/qs
361 rh(kp1) = max( 0.0_dp, min( 1.0_dp, wrk) ) ! relative huminity
363 if( p_qc > 1 ) then
364 xlwc(kp1) = 1.0e3_dp*moist(i,k,j,p_qc)*rho_phy(i,k,j) ! cloud water(g/m^3)
365 end if
366 if( xlwc(kp1) < 1.0e-6_dp ) then
367 xlwc(kp1) = 0.0_dp
368 end if
370 zh(kp1) = .5_dp*(z_at_w(i,kp1,j)+z_at_w(i,k,j))*0.001_dp - z_at_w(i,kts,j)*0.001_dp ! height(km)
371 ! zh(k+1) = z(i,k,j)*0.001 - z_at_w(i,1,j)*0.001 ! height(km)
373 !===============================================================================
374 ! XUEXI
375 ! ADD the aerosol effect on J
376 ! chem is in the unit of ug/m3
377 ! cb1 = p_ecj * 1.0
378 ! cb2 = p_ecj * 0.0
379 ! oc1 = p_orgpaj ( Primary organic conc. Aitken mode)
380 ! + p_p25aj ( Unknown PM25 conc. Acc. mode )
381 ! ant = p_no3aj ( Nitrate conc. Acc. mode )
382 ! + p_nh4aj ( Ammonium conc. Acc. mode )
383 ! so4 = p_so4aj ( Sulfate conc. Acc. mode )
384 ! sa1 = p_seas ( Coarse soil-derived aero sols )
385 !
386 !=======================================================
388 acb1(kp1) = 0._dp
389 acb2(kp1) = 0._dp
390 aoc1(kp1) = 0._dp
391 aoc2(kp1) = 0._dp
392 aant(kp1) = 0._dp
393 aso4(kp1) = 0._dp
394 asal(kp1) = 0._dp
395 ! rajesh: initialize aerosol optical properties to zero
396 tauaer300(k) = 0._dp
397 tauaer400(k) = 0._dp
398 tauaer600(k) = 0._dp
399 tauaer999(k) = 0._dp
400 waer300(k) = 0._dp
401 waer400(k) = 0._dp
402 waer600(k) = 0._dp
403 waer999(k) = 0._dp
404 gaer300(k) = 0._dp
405 gaer400(k) = 0._dp
406 gaer600(k) = 0._dp
407 gaer999(k) = 0._dp
409 if( isorg == 1 ) then
410 acb1(kp1) = chem(i,k,j,p_ecj)
411 ! acb2(kp1) = 0.0_dp
412 ! aoc1(kp1) = chem(i,k,j,p_orgpaj) + chem(i,k,j,p_orgbaj) + chem(i,k,j,p_orgaj)
414 if(config_flags%chem_opt == CB05_SORG_VBS_AQ_KPP) then
415 aoc1(kp1) = chem(i,k,j,p_orgpaj)+chem(i,k,j,p_asoa1j)+chem(i,k,j,p_asoa2j) &
416 + chem(i,k,j,p_asoa3j)+chem(i,k,j,p_asoa4j) &
417 + chem(i,k,j,p_bsoa1j)+chem(i,k,j,p_bsoa2j) &
418 + chem(i,k,j,p_bsoa3j)+chem(i,k,j,p_bsoa4j)
419 else
420 aoc1(kp1) = chem(i,k,j,p_orgpaj)
421 endif
422 ! aoc1(kp1) = chem(i,k,j,p_orgpaj)
423 ! aoc2(kp1) = 0.0_dp
424 aant(kp1) = chem(i,k,j,p_no3aj) + chem(i,k,j,p_nh4aj)
425 aso4(kp1) = chem(i,k,j,p_so4aj)
426 ! asal(kp1) = chem(i,k,j,p_seas)
427 asal(kp1) = chem(i,k,j,p_seas) + chem(i,k,j,p_naaj) + chem(i,k,j,p_claj)
428 elseif( config_flags%chem_opt == MOZART_KPP .or. &
429 config_flags%chem_opt == MOZCART_KPP .or. &
430 config_flags%chem_opt == T1_MOZCART_KPP ) then
431 aso4(kp1) = chem(i,k,j,p_sulf)*air(kp1)*real(mw_so4_aer,kind=dp)*kg_per_amu*1.e-9_dp
432 if( config_flags%chem_opt == MOZCART_KPP .or. &
433 config_flags%chem_opt == T1_MOZCART_KPP ) then
434 atm_mass_den = rho_phy(i,k,j)
435 acb1(kp1) = chem(i,k,j,p_bc1)*atm_mass_den
436 acb2(kp1) = chem(i,k,j,p_bc2)*atm_mass_den
437 aoc1(kp1) = chem(i,k,j,p_oc1)*atm_mass_den
438 aoc2(kp1) = chem(i,k,j,p_oc2)*atm_mass_den
439 ! aant(kp1) = 0._dp
440 asal(kp1) = (chem(i,k,j,p_seas_1) + chem(i,k,j,p_seas_2) &
441 + chem(i,k,j,p_seas_3) + chem(i,k,j,p_seas_4))*atm_mass_den
442 endif
443 endif
444 ! rajesh: Extract column of aerosol optical properties
445 if(config_flags%aer_ra_feedback == 1) then
446 tauaer300(k) = tauaer1(i,k,j)
447 tauaer400(k) = tauaer2(i,k,j)
448 tauaer600(k) = tauaer3(i,k,j)
449 tauaer999(k) = tauaer4(i,k,j)
450 waer300(k) = waer1(i,k,j)
451 waer400(k) = waer2(i,k,j)
452 waer600(k) = waer3(i,k,j)
453 waer999(k) = waer4(i,k,j)
454 gaer300(k) = gaer1(i,k,j)
455 gaer400(k) = gaer2(i,k,j)
456 gaer600(k) = gaer3(i,k,j)
457 gaer999(k) = gaer4(i,k,j)
458 endif
459 ! else
460 ! acb1(kp1) = 0._dp
461 ! acb2(kp1) = 0._dp
462 ! aoc1(kp1) = 0._dp
463 ! aoc2(kp1) = 0._dp
464 ! aant(kp1) = 0._dp
465 ! aso4(kp1) = 0._dp
466 ! asal(kp1) = 0._dp
467 enddo level_loop
469 temp(1) = t8w(i,kts,j)
470 air(1) = xair_kg*( p8w(i,kts,j)/t8w(i,kts,j)/r_d )
471 o3(1) = o3(2)
472 rh(1) = rh(2)
473 xlwc(1) = 0._dp
474 zh(1) = 0._dp
475 acb1(1) = acb1(2)
476 acb2(1) = acb2(2)
477 aoc1(1) = aoc1(2)
478 aoc2(1) = aoc2(2)
479 aant(1) = aant(2)
480 aso4(1) = aso4(2)
481 asal(1) = asal(2)
482 ! rajesh
483 ! tauaer300(1) = tauaer300(2)
484 ! tauaer400(1) = tauaer400(2)
485 ! tauaer600(1) = tauaer600(2)
486 ! tauaer999(1) = tauaer999(2)
487 ! waer300(1) = waer300(2)
488 ! waer400(1) = waer400(2)
489 ! waer600(1) = waer600(2)
490 ! waer999(1) = waer999(2)
491 ! gaer300(1) = gaer300(2)
492 ! gaer400(1) = gaer400(2)
493 ! gaer600(1) = gaer600(2)
494 ! gaer999(1) = gaer999(2)
496 !----------------------------------------------------
497 ! smooth air density ...
498 !----------------------------------------------------
499 do k = kts+1, kte-1
500 air(k) = .5_dp*air(k) + .25_dp*( air(k-1) + air(k+1) )
501 enddo
503 alat = real(xlat(i,j),kind=dp)
504 along = -real(xlong(i,j),kind=dp)
505 call calc_zenith( alat, along, julday, gmtp, azim, zen )
506 zen_angle(i,j) = zen
508 call photo( config_flags%chem_opt,&
509 kte+1, &
510 tuv_jmax, &
511 julday, &
512 gmtp, &
513 alat, along, &
514 o3top, o2top, dobson, &
515 lu, &
516 zh, &
517 temp, &
518 air, &
519 rh, &
520 xlwc, &
521 o3, &
522 acb1, &
523 acb2, &
524 aoc1, &
525 aoc2, &
526 aant, &
527 aso4, &
528 asal, &
529 ! rajesh: Add aerosol optical property columns as arguments
530 tauaer300, tauaer400, &
531 tauaer600, tauaer999, &
532 waer300, waer400, &
533 waer600, waer999, &
534 gaer300, gaer400, &
535 gaer600, gaer999, &
536 config_flags%aer_ra_feedback, &
537 prate, &
538 radfld, &
539 adjcoe, &
540 prate0 )
542 do k = kts, kte
543 ph_no2(i,k,j) = max( 0._dp,prate(k,pid_no2)*m2s )
544 ph_o31d(i,k,j) = max( 0._dp,prate(k,pid_o31d)*m2s )
545 ph_o33p(i,k,j) = max( 0._dp,prate(k,pid_o33p)*m2s )
546 ph_hono(i,k,j) = max( 0._dp,prate(k,pid_hono)*m2s )
547 ph_hno3(i,k,j) = max( 0._dp,prate(k,pid_hno3)*m2s )
548 ph_hno4(i,k,j) = max( 0._dp,prate(k,pid_hno4)*m2s )
549 ph_no3o2(i,k,j) = max( 0._dp,prate(k,pid_no3o2)*m2s )
550 ph_no3o(i,k,j) = max( 0._dp,prate(k,pid_no3o)*m2s*ph_no3o_factor )
551 ph_h2o2(i,k,j) = max( 0._dp,prate(k,pid_h2o2)*m2s )
552 ph_ch2om(i,k,j) = max( 0._dp,prate(k,pid_ch2om)*m2s )
553 ph_ch2or(i,k,j) = max( 0._dp,prate(k,pid_ch2or)*m2s )
554 ph_ald(i,k,j) = max( 0._dp,prate(k,pid_ald)*m2s )
555 ph_op1(i,k,j) = max( 0._dp,prate(k,pid_op1)*m2s )
556 ph_op2(i,k,j) = max( 0._dp,prate(k,pid_op2)*m2s )
557 ph_paa(i,k,j) = max( 0._dp,prate(k,pid_paa)*m2s )
558 ph_ket(i,k,j) = max( 0._dp,prate(k,pid_ket)*m2s )
559 ph_glya(i,k,j) = max( 0._dp,prate(k,pid_glya)*m2s )
560 ph_glyb(i,k,j) = max( 0._dp,prate(k,pid_glyb)*m2s )
561 ph_mgly(i,k,j) = max( 0._dp,prate(k,pid_mgly)*m2s )
562 ph_dcb(i,k,j) = max( 0._dp,prate(k,pid_dcb)*m2s )
563 ph_onit(i,k,j) = max( 0._dp,prate(k,pid_onit)*m2s )
564 ph_macr(i,k,j) = max( 0._dp,prate(k,pid_macr) *m2s )
565 ph_ch3coc2h5(i,k,j) = max( 0._dp,prate(k,23)*m2s )
566 ph_cl2(i,k,j) = max( 0._dp,prate(k,pid_cl2)*m2s )
567 ph_hocl(i,k,j) = max( 0._dp,prate(k,pid_hocl)*m2s )
568 ph_fmcl(i,k,j) = max( 0._dp,prate(k,pid_fmcl)*m2s )
569 ph_pan(i,k,j) = max( 0._dp,prate(k,pid_pan)*m2s )
570 ph_n2o5(i,k,j) = max( 0._dp,prate(k,pid_n2o5)*m2s )
571 enddo
573 wc_x(1:min(ntuv_bin,nw-1)) = wc(1:min(ntuv_bin,nw-1))
574 zref_x(1:min(nz,ntuv_lev)) = zref(1:min(nz,ntuv_lev))
576 if( chm_is_mozart ) then
577 do k = kts, kte
578 ph_n2o(i,k,j) = max( 0._dp,prate(k,pid_n2o)*m2s )
579 ph_n2o5(i,k,j) = max( 0._dp,prate(k,pid_n2o5)*m2s )
580 ph_mvk(i,k,j) = max( 0._dp,prate(k,pid_mvk)*m2s )
581 ph_pan(i,k,j) = max( 0._dp,prate(k,pid_pan)*m2s )
582 ph_ald(i,k,j) = max( 0._dp,(prate(k,pid_ch3choa)+prate(k,pid_ch3chob)+prate(k,pid_ch3choc))*m2s )
583 ph_gly(i,k,j) = ph_mgly(i,k,j)
584 ph_mek(i,k,j) = ph_ket(i,k,j)
585 ph_acetol(i,k,j) = ph_op1(i,k,j)
586 ph_pooh(i,k,j) = ph_op1(i,k,j)
587 ph_glyald(i,k,j) = max( 0._dp,prate(k,pid_glyald) *m2s )
588 ph_hyac(i,k,j) = max( 0._dp,2._dp*prate(k,pid_hyac) *m2s )
589 ph_hno4(i,k,j) = ph_hno4(i,k,j) + 1.e-5_dp*m2s
590 enddo
591 end if
593 enddo I_TILE_LOOP
594 enddo J_TILE_LOOP
596 end subroutine ftuv_driver
598 !++alma 2012-12-01 modis landse from
599 ! subroutine ftuv_init( id,ips, ipe, jps, jpe, kte, &
600 ! ide, jde, config_flags )
601 subroutine ftuv_init( id,ips, ipe, jps, jpe, kte, &
602 ide, jde, config_flags, num_land_cat, mminlu_loc)
603 !--
605 !---------------------------------------------------------------------
606 ! ... new initialization routine for ftuv
607 !---------------------------------------------------------------------
609 use module_state_description, only : mozart_kpp, mozcart_kpp, &
610 t1_mozcart_kpp,&
611 mozart_mosaic_4bin_kpp,&
612 mozart_mosaic_4bin_aq_kpp
613 use module_ftuv_subs, only : aer_init
614 use module_configure, only : grid_config_rec_type
616 implicit none
618 !---------------------------------------------------------------------
619 ! ... dummy arguments
620 !---------------------------------------------------------------------
621 integer, intent(in) :: id
622 integer, intent(in) :: ips, ipe, jps, jpe, kte
623 integer, intent(in) :: ide, jde
625 !++alma 2012-12-01 modis landuse from sw
626 integer, intent(in) :: num_land_cat
627 character(len=*), intent(in) :: mminlu_loc
628 type(grid_config_rec_type), intent(in ) :: config_flags
629 !--
632 !---------------------------------------------------------------------
633 ! ... local variables
634 !---------------------------------------------------------------------
635 integer :: astat
636 integer :: ncid
637 integer :: dimid
638 integer :: varid
639 integer :: max_dom
640 integer :: cpos
641 integer :: iend, jend
642 integer :: lon_e, lat_e
643 integer :: ncoldens_levs
644 integer :: ndays_of_year
645 integer :: i
646 real, allocatable :: coldens(:,:,:,:)
647 character(len=128) :: err_msg
648 character(len=64) :: filename
649 character(len=3) :: id_num
651 LOGICAL , EXTERNAL :: wrf_dm_on_monitor
653 #ifdef NETCDF
654 include 'netcdf.inc'
655 #define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes ( A , size ( A ) * DWORDSIZE )
656 #else
657 if( config_flags%chem_opt == MOZART_KPP .or. &
658 config_flags%chem_opt == MOZCART_KPP .or. &
659 config_flags%chem_opt == T1_MOZCART_KPP .or. &
660 config_flags%chem_opt == MOZART_MOSAIC_4BIN_KPP .or. &
661 config_flags%chem_opt == MOZART_MOSAIC_4BIN_AQ_KPP ) then
662 call wrf_message( 'ftuv_init: mozart,mozcart,mozart_mosaic_* chem option requires netcdf' )
663 call wrf_abort
664 end if
665 #endif
667 !---------------------------------------------------------------------
668 ! ... general initializer
669 !---------------------------------------------------------------------
670 if( id == 1 .and. .not. photo_inti_initialized ) then
671 write(err_msg,*) 'ftuv_init: calling photo_inti for id = ',id
672 call wrf_message( trim(err_msg) )
673 call photo_inti( config_flags%chem_opt, kte, 20._dp )
674 endif
675 #ifdef NETCDF
676 is_mozart : &
677 if( config_flags%chem_opt == MOZART_KPP .or. &
678 config_flags%chem_opt == MOZCART_KPP .or. &
679 config_flags%chem_opt == T1_MOZCART_KPP .or. &
680 config_flags%chem_opt == MOZART_MOSAIC_4BIN_KPP .or. &
681 config_flags%chem_opt == MOZART_MOSAIC_4BIN_AQ_KPP ) then
682 !---------------------------------------------------------------------
683 ! ... allocate column_density type
684 !---------------------------------------------------------------------
685 if( id == 1 .and. .not. allocated(col_dens) ) then
686 CALL nl_get_max_dom( 1,max_dom )
687 allocate( col_dens(max_dom),stat=astat )
688 if( astat /= 0 ) then
689 CALL wrf_message( 'ftuv_init: failed to allocate column_density type col_dens' )
690 CALL wrf_abort
691 end if
692 write(err_msg,*) 'ftuv_init: intializing ',max_dom,' domains'
693 call wrf_message( trim(err_msg) )
694 col_dens(:)%is_allocated = .false.
695 endif
696 !---------------------------------------------------------------------
697 ! ... open column density netcdf file
698 !---------------------------------------------------------------------
699 col_dens_allocated : &
700 if( .not. col_dens(id)%is_allocated ) then
701 if( wrf_dm_on_monitor() ) then
702 cpos = index( config_flags%exo_coldens_inname, '<domain>' )
703 if( cpos > 0 ) then
704 write(id_num,'(i3)') 100+id
705 filename = config_flags%exo_coldens_inname(:cpos-1) // 'd' // id_num(2:3)
706 else
707 filename = trim( config_flags%exo_coldens_inname )
708 endif
709 err_msg = 'ftuv_init: intializing domain ' // id_num(2:3)
710 call wrf_message( trim(err_msg) )
711 err_msg = 'ftuv_init: failed to open file ' // trim(filename)
712 call handle_ncerr( nf_open( trim(filename), nf_noclobber, ncid ), trim(err_msg) )
713 !---------------------------------------------------------------------
714 ! ... get dimensions
715 !---------------------------------------------------------------------
716 err_msg = 'ftuv_init: failed to get col_dens levels id'
717 call handle_ncerr( nf_inq_dimid( ncid, 'coldens_levs', dimid ), trim(err_msg) )
718 err_msg = 'ftuv_init: failed to get col_dens levels'
719 call handle_ncerr( nf_inq_dimlen( ncid, dimid, ncoldens_levs ), trim(err_msg) )
720 err_msg = 'ftuv_init: failed to get number of days in year id'
721 call handle_ncerr( nf_inq_dimid( ncid, 'ndays_of_year', dimid ), trim(err_msg) )
722 err_msg = 'ftuv_init: failed to get number of days in year'
723 call handle_ncerr( nf_inq_dimlen( ncid, dimid, ndays_of_year ), trim(err_msg) )
724 err_msg = 'ftuv_init: failed to get west_east id'
725 call handle_ncerr( nf_inq_dimid( ncid, 'west_east', dimid ), trim(err_msg) )
726 err_msg = 'ftuv_init: failed to get west_east'
727 call handle_ncerr( nf_inq_dimlen( ncid, dimid, lon_e ), trim(err_msg) )
728 err_msg = 'ftuv_init: failed to get south_north id'
729 call handle_ncerr( nf_inq_dimid( ncid, 'south_north', dimid ), trim(err_msg) )
730 err_msg = 'ftuv_init: failed to get south_north'
731 call handle_ncerr( nf_inq_dimlen( ncid, dimid, lat_e ), trim(err_msg) )
732 end IF
733 !---------------------------------------------------------------------
734 ! ... bcast the dimensions
735 !---------------------------------------------------------------------
736 CALL wrf_dm_bcast_bytes ( ncoldens_levs , IWORDSIZE )
737 CALL wrf_dm_bcast_bytes ( ndays_of_year , IWORDSIZE )
738 CALL wrf_dm_bcast_bytes ( lon_e , IWORDSIZE )
739 CALL wrf_dm_bcast_bytes ( lat_e , IWORDSIZE )
740 !---------------------------------------------------------------------
741 ! ... allocate local arrays
742 !---------------------------------------------------------------------
743 iend = min( ipe,ide-1 )
744 jend = min( jpe,jde-1 )
745 allocate( coldens(lon_e,lat_e,ncoldens_levs,ndays_of_year), stat=astat )
746 if( astat /= 0 ) then
747 call wrf_message( 'ftuv_init: failed to allocate coldens' )
748 call wrf_abort
749 end if
750 !---------------------------------------------------------------------
751 ! ... allocate column_density type component arrays
752 !---------------------------------------------------------------------
753 col_dens(id)%ncoldens_levs = ncoldens_levs
754 col_dens(id)%ndays_of_year = ndays_of_year
755 allocate( col_dens(id)%col_levs(ncoldens_levs), &
756 col_dens(id)%day_of_year(ndays_of_year), stat=astat )
757 if( astat /= 0 ) then
758 call wrf_message( 'ftuv_init: failed to allocate col_levs,day_of_year' )
759 call wrf_abort
760 end if
761 allocate( col_dens(id)%o3_col_dens(ips:iend,jps:jend,ncoldens_levs,ndays_of_year), &
762 col_dens(id)%o2_col_dens(ips:iend,jps:jend,ncoldens_levs,ndays_of_year), stat=astat )
763 if( astat /= 0 ) then
764 call wrf_message( 'ftuv_init: failed to allocate o3_col_dens,o2_col_dens' )
765 call wrf_abort
766 end if
767 col_dens(id)%is_allocated = .true.
768 !---------------------------------------------------------------------
769 ! ... read arrays
770 !---------------------------------------------------------------------
771 IF ( wrf_dm_on_monitor() ) THEN
772 err_msg = 'ftuv_init: failed to get col_levs variable id'
773 call handle_ncerr( nf_inq_varid( ncid, 'coldens_levs', varid ), trim(err_msg) )
774 err_msg = 'ftuv_init: failed to read col_levs variable'
775 call handle_ncerr( nf_get_var_double( ncid, varid, col_dens(id)%col_levs ), trim(err_msg) )
776 err_msg = 'ftuv_init: failed to get days_of_year variable id'
777 call handle_ncerr( nf_inq_varid( ncid, 'days_of_year', varid ), trim(err_msg) )
778 err_msg = 'ftuv_init: failed to read days_of_year variable'
779 call handle_ncerr( nf_get_var_double( ncid, varid, col_dens(id)%day_of_year ), trim(err_msg) )
780 err_msg = 'ftuv_init: failed to col_dens variable id'
781 call handle_ncerr( nf_inq_varid( ncid, 'o3_column_density', varid ), trim(err_msg) )
782 err_msg = 'ftuv_init: failed to read col_dens variable'
783 call handle_ncerr( nf_get_var_real( ncid, varid, coldens ), trim(err_msg) )
784 end if
786 DM_BCAST_MACRO(col_dens(id)%col_levs)
787 write(*,*) 'ftuv_init: bcast col_levs'
788 DM_BCAST_MACRO(col_dens(id)%day_of_year)
789 write(*,*) 'ftuv_init: bcast day_of_year'
790 CALL wrf_dm_bcast_bytes ( coldens, size(coldens)*RWORDSIZE )
791 write(*,*) 'ftuv_init: bcast o3_col_dens'
793 col_dens(id)%o3_col_dens(ips:iend,jps:jend,:ncoldens_levs,:ndays_of_year) = &
794 coldens(ips:iend,jps:jend,:ncoldens_levs,:ndays_of_year)
796 IF ( wrf_dm_on_monitor() ) THEN
797 err_msg = 'ftuv_init: failed to col_dens variable id'
798 call handle_ncerr( nf_inq_varid( ncid, 'o2_column_density', varid ), trim(err_msg) )
799 err_msg = 'ftuv_init: failed to read col_dens variable'
800 call handle_ncerr( nf_get_var_real( ncid, varid, coldens ), trim(err_msg) )
801 !---------------------------------------------------------------------
802 ! ... close column density netcdf file
803 !---------------------------------------------------------------------
804 err_msg = 'ftuv_init: failed to close file ' // trim(filename)
805 call handle_ncerr( nf_close( ncid ), trim(err_msg) )
806 end if
808 CALL wrf_dm_bcast_bytes ( coldens, size(coldens)*RWORDSIZE )
809 write(*,*) 'ftuv_init: bcast o2_col_dens'
811 col_dens(id)%o2_col_dens(ips:iend,jps:jend,:ncoldens_levs,:ndays_of_year) = &
812 coldens(ips:iend,jps:jend,:ncoldens_levs,:ndays_of_year)
814 deallocate( coldens )
816 write(*,*) ' '
817 write(*,*) 'ftuv_init: coldens variables for domain ',id
818 write(*,*) 'ftuv_init: days_of_year'
819 write(*,'(1p,5g15.7)') col_dens(id)%day_of_year(:)
820 write(*,*) 'ftuv_init: coldens levels'
821 write(*,'(1p,5g15.7)') col_dens(id)%col_levs(:)
822 write(*,*) ' '
823 endif col_dens_allocated
825 !++alma 2012-12-01 modis landuse from sw
826 ! if( id == 1 ) then
827 if( id == 1 .and. .not. allocated(luse2usgs) ) then
828 !allocate( luse2usgs(config_flags%num_land_cat),stat=astat )
829 print*,"num_land_cat: ", num_land_cat
830 allocate( luse2usgs(num_land_cat),stat=astat )
831 if( astat /= 0 ) then
832 CALL wrf_message( 'ftuv_init: failed to allocate luse2usgs')
833 CALL wrf_abort
834 end if
835 if( trim(mminlu_loc) == 'USGS' ) then
836 !luse2usgs(:) = (/ (i,i=1,config_flags%num_land_cat) /)
837 luse2usgs(:) = (/ (i,i=1,num_land_cat) /)
838 elseif( trim(mminlu_loc) == 'MODIFIED_IGBP_MODIS_NOAH' ) then
839 luse2usgs(:) = (/ 14,13,12,11,15,8,9,10,10,7, &
840 17,4,1,5,24,19,16,21,22,23 /)
841 endif
842 endif
843 !--
846 endif is_mozart
847 #endif
848 !---------------------------------------------------------------------
849 ! ... initialize aerosol routine
850 !---------------------------------------------------------------------
851 if( id == 1 .and. .not. photo_inti_initialized ) then
852 write(*,*) 'ftuv_init: initializing aerosol routine'
853 call aer_init
854 write(*,*) 'ftuv_init: finished'
855 photo_inti_initialized = .true.
856 endif
858 end subroutine ftuv_init
860 #ifdef NETCDF
861 subroutine handle_ncerr( ret, mes )
862 !---------------------------------------------------------------------
863 ! ... netcdf error handling routine
864 !---------------------------------------------------------------------
866 implicit none
868 !---------------------------------------------------------------------
869 ! ... dummy arguments
870 !---------------------------------------------------------------------
871 integer, intent(in) :: ret
872 character(len=*), intent(in) :: mes
874 include 'netcdf.inc'
876 if( ret /= nf_noerr ) then
877 call wrf_message( trim(mes) )
878 call wrf_message( trim(nf_strerror(ret)) )
879 call wrf_abort
880 end if
882 end subroutine handle_ncerr
883 #endif
885 subroutine photo_inti( chem_opt, nlev, ztopin )
887 use module_wave_data, only : nw, wl, wave_data_inti
888 use module_ftuv_subs, only : nwint, wlint, schu_inti, inter_inti
889 use module_state_description, only : mozart_kpp, mozcart_kpp, &
890 t1_mozcart_kpp,&
891 mozart_mosaic_4bin_kpp,&
892 mozart_mosaic_4bin_aq_kpp
894 implicit none
896 !---------------------------------------------------------------------
897 ! input arguments
898 !---------------------------------------------------------------------
899 integer, intent(in) :: chem_opt
900 integer, intent(in) :: nlev
901 real(dp), intent(in) :: ztopin
903 !---------------------------------------------------------------------
904 ! local arguments
905 !---------------------------------------------------------------------
906 integer :: k, kdis, nabv, iw
907 real(dp) :: ztop
909 if( chem_opt /= MOZART_KPP .and. &
910 chem_opt /= MOZCART_KPP .and. &
911 chem_opt /= T1_MOZCART_KPP .and. &
912 chem_opt /= MOZART_MOSAIC_4BIN_KPP .and. &
913 chem_opt /= MOZART_MOSAIC_4BIN_AQ_KPP ) then
914 !---------------------------------------------------------------------
915 ! change unit from m to km
916 !---------------------------------------------------------------------
917 ztop = ztopin
918 !---------------------------------------------------------------------
919 ! set the top exo model level
920 !---------------------------------------------------------------------
921 do k = 1,nref
922 if( ztop < zref(k) ) then
923 exit
924 end if
925 enddo
926 if( k == 1 .or. k > nref ) then
927 call wrf_message( 'photo_inti: ztop is not in zref range' )
928 call wrf_abort
929 endif
931 kcon = k + 1
932 kdis = nref - kcon
933 nabv = kdis/2 + 1
934 if( mod(kdis,2) /= 0 ) then
935 nabv = nabv + 1
936 endif
937 nz = nlev + nabv
938 write(*,*) '******************************************'
939 write(*,*) 'photo_inti: kcon,kdis,nabv,nlev,nz = ',kcon,kdis,nabv,nlev,nz
940 write(*,*) '******************************************'
941 dobson = 265.0_dp
942 else
943 nz = nlev
944 dobson = 0._dp
945 endif
947 !---------------------------------------------------------------------
948 ! set albedo for land
949 !---------------------------------------------------------------------
950 do iw = 1, nw-1
951 if( wl(iw)<400.0_dp ) then
952 albedo(iw,1) = 0.05_dp
953 else if( (wl(iw)>=400.0_dp ) .and. (wl(iw)<450.0_dp) ) then
954 albedo(iw,1) = 0.06_dp
955 else if( (wl(iw)>=450.0_dp ) .and. (wl(iw)<500.0_dp) ) then
956 albedo(iw,1) = 0.08_dp
957 else if( (wl(iw)>=500.0_dp ) .and. (wl(iw)<550.0_dp) ) then
958 albedo(iw,1) = 0.10_dp
959 else if( (wl(iw)>=550.0_dp ) .and. (wl(iw)<600.0_dp) ) then
960 albedo(iw,1) = 0.11_dp
961 else if( (wl(iw)>=600.0_dp ) .and. (wl(iw)<640.0_dp) ) then
962 albedo(iw,1) = 0.12_dp
963 else if( (wl(iw)>=640.0_dp ) .and. (wl(iw)<660.0_dp) ) then
964 albedo(iw,1) = 0.135_dp
965 else if( wl(iw)>=660.0_dp ) then
966 albedo(iw,1) = 0.15_dp
967 end if
968 end do
970 !---------------------------------------------------------------------
971 ! set albedo for sea water
972 !
973 ! For mozart_kpp these are interpolated from reflectance values
974 ! for IGBP type 17 (ocean water). Source: NASA LARC
975 ! (http://www-surf.larc.nasa.gov/surf/pages/explan.html)
976 !---------------------------------------------------------------------
977 if( chem_opt == MOZART_KPP .or. &
978 chem_opt == MOZCART_KPP .or. &
979 chem_opt == T1_MOZCART_KPP .or. &
980 chem_opt == MOZART_MOSAIC_4BIN_KPP .or. &
981 chem_opt == MOZART_MOSAIC_4BIN_AQ_KPP) then
982 albedo(1:17,2) = (/ 0.0747_dp, 0.0755_dp, 0.0767_dp, 0.0783_dp, 0.0802_dp, &
983 0.0825_dp, 0.0852_dp, 0.0882_dp, 0.0914_dp, 0.0908_dp, &
984 0.0763_dp, 0.0725_dp, 0.0689_dp, 0.0632_dp, 0.0570_dp, &
985 0.0585_dp, 0.0535_dp /)
986 else
987 albedo(1:17,2) = (/ 0.8228_dp, 0.8140_dp, 0.8000_dp, 0.7820_dp, 0.7600_dp, &
988 0.7340_dp, 0.7040_dp, 0.6700_dp, 0.6340_dp, 0.4782_dp, &
989 0.3492_dp, 0.3000_dp, 0.2700_dp, 0.2400_dp, 0.1700_dp, &
990 0.0800_dp, 0.0600_dp /)
991 endif
993 !---------------------------------------------------------------------
994 ! intialize other modules
995 !---------------------------------------------------------------------
996 call schu_inti
997 call wave_data_inti( nz )
998 call inter_inti( nw, wl, nwint, wlint )
1000 end subroutine photo_inti
1002 subroutine photo( chem_opt, nlev, njout, julday, gmtp, alat, &
1003 along, o3top, o2top, o3toms, lu, zin, tlevin, &
1004 airlevin, rhin, xlwcin, o3in, acb1in, &
1005 acb2in, aoc1in, aoc2in, aantin, aso4in, &
1006 asalin, tauaer300in,tauaer400in,tauaer600in,tauaer999in, &
1007 waer300in, waer400in, waer600in, waer999in, &
1008 gaer300in, gaer400in, gaer600in, gaer999in, &
1009 aer_ra_feedback, prate, radfld, adjcoe, prate0 )
1011 use module_ftuv_subs, only : calc_zenith, photoin
1013 implicit none
1015 !----------------------------------------------------
1016 ! input arguments
1017 !----------------------------------------------------
1018 integer, intent(in) :: chem_opt
1019 integer, intent(in) :: lu
1020 integer, intent(in) :: nlev, njout
1021 integer, intent(in) :: julday
1022 real(dp), intent(in) :: gmtp
1023 real(dp), intent(in) :: alat, along
1024 real(dp), intent(in) :: o3toms
1025 real(dp), intent(in) :: o3top
1026 real(dp), intent(in) :: o2top
1028 real(dp), intent(in) :: zin(nlev)
1029 real(dp), intent(in) :: tlevin(nlev)
1030 real(dp), intent(in) :: airlevin(nlev)
1031 real(dp), intent(in) :: rhin(nlev)
1032 real(dp), intent(in) :: xlwcin(nlev)
1033 real(dp), intent(in) :: o3in(nlev)
1034 real(dp), intent(in) :: acb1in(nlev)
1035 real(dp), intent(in) :: acb2in(nlev)
1036 real(dp), intent(in) :: aoc1in(nlev)
1037 real(dp), intent(in) :: aoc2in(nlev)
1038 real(dp), intent(in) :: aantin(nlev)
1039 real(dp), intent(in) :: aso4in(nlev)
1040 real(dp), intent(in) :: asalin(nlev)
1042 ! rajesh: declare dust variables
1043 real(dp), intent(in) :: tauaer300in(nlev-1)
1044 real(dp), intent(in) :: tauaer400in(nlev-1)
1045 real(dp), intent(in) :: tauaer600in(nlev-1)
1046 real(dp), intent(in) :: tauaer999in(nlev-1)
1047 real(dp), intent(in) :: waer300in(nlev-1)
1048 real(dp), intent(in) :: waer400in(nlev-1)
1049 real(dp), intent(in) :: waer600in(nlev-1)
1050 real(dp), intent(in) :: waer999in(nlev-1)
1051 real(dp), intent(in) :: gaer300in(nlev-1)
1052 real(dp), intent(in) :: gaer400in(nlev-1)
1053 real(dp), intent(in) :: gaer600in(nlev-1)
1054 real(dp), intent(in) :: gaer999in(nlev-1)
1055 INTEGER, INTENT(IN) :: aer_ra_feedback
1057 !----------------------------------------------------
1058 ! output arguments
1059 !----------------------------------------------------
1060 real(dp), intent(out) :: prate(nlev,njout)
1061 real(dp), intent(out) :: radfld(nz,nw-1)
1062 real(dp), intent(out) :: adjcoe(nz,njout)
1063 real(dp), intent(out) :: prate0(nz,njout)
1065 !----------------------------------------------------
1066 ! local arugments
1067 !----------------------------------------------------
1068 integer :: astat, istat
1069 integer :: k, n
1070 real(dp) :: azim, zen
1071 real(dp), allocatable :: z_ph(:), tlev_ph(:), tlay_ph(:), airlev_ph(:)
1072 real(dp), allocatable :: rh_ph(:), xlwc_ph(:)
1073 real(dp), allocatable :: o3_ph(:)
1074 real(dp), allocatable :: acb1_ph(:), acb2_ph(:)
1075 real(dp), allocatable :: aoc1_ph(:), aoc2_ph(:)
1076 real(dp), allocatable :: aant_ph(:), aso4_ph(:), asal_ph(:)
1078 ! rajesh: declare allocatable variables for aerosol optical properties
1079 real(dp), allocatable :: tauaer300_ph(:), tauaer400_ph(:), &
1080 tauaer600_ph(:), tauaer999_ph(:), &
1081 waer300_ph(:), waer400_ph(:), &
1082 waer600_ph(:), waer999_ph(:), &
1083 gaer300_ph(:), gaer400_ph(:), &
1084 gaer600_ph(:), gaer999_ph(:)
1085 real(dp), allocatable :: ftuv(:,:)
1087 !-----------------------------------------------------------------------------
1088 ! Now xx-ftuv only considers the following 28 photolysis reactions ...
1089 !-----------------------------------------------------------------------------
1090 ! 1 o2 + hv -> o + o
1091 ! 2 o3 -> o2 + o(1d)
1092 ! 3 o3 -> o2 + o(3p)
1093 ! 4 no2 -> no + o(3p)
1094 ! 5 no3 -> no + o2
1095 ! 6 no3 -> no2 + o(3p)
1096 ! 7 n2o5 -> no3 + no + o(3p)
1097 ! 8 n2o5 -> no3 + no2
1098 ! 9 n2o + hv -> n2 + o(1d)
1099 ! 10 ho2 + hv -> oh + o
1100 ! 11 h2o2 -> 2 oh
1101 ! 12 hno2 -> oh + no
1102 ! 13 hno3 -> oh + no2
1103 ! 14 hno4 -> ho2 + no2
1104 ! 15 ch2o -> h + hco
1105 ! 16 ch2o -> h2 + co
1106 ! 17 ch3cho -> ch3 + hco
1107 ! 18 ch3cho -> ch4 + co
1108 ! 19 ch3cho -> ch3co + h
1109 ! 20 c2h5cho -> c2h5 + hco
1110 ! 21 chocho -> products
1111 ! 22 ch3cocho -> products
1112 ! 23 ch3coch3
1113 ! 24 ch3ooh -> ch3o + oh
1114 ! 25 ch3ono2 -> ch3o+no2
1115 ! 26 pan + hv -> products
1116 ! 27 mvk + hv -> products
1117 ! 28 macr + hv -> products
1118 ! 29 hyac + hv -> products
1119 ! 30 glyald + hv -> products
1120 !-----------------------------------------------------------------------------
1122 !-----------------------------------------------------------------------------
1123 ! solar zenith angle calculation
1124 !-----------------------------------------------------------------------------
1125 call calc_zenith( alat, along, julday, gmtp, azim, zen )
1126 if( zen == 90.0_dp ) then
1127 zen = 89.0_dp
1128 elseif( zen >= 95.0_dp ) then
1129 prate(:,:) = 0.0_dp
1130 return
1131 endif
1133 !---------------------------------------------------------------------
1134 ! allocate memory space
1135 !---------------------------------------------------------------------
1136 astat = 0
1137 allocate( ftuv(nz,njout), stat=istat )
1138 astat = astat + istat
1139 allocate( z_ph(nz), tlev_ph(nz), tlay_ph(nz-1), airlev_ph(nz), stat=istat )
1140 astat = astat + istat
1141 allocate( rh_ph(nz), xlwc_ph(nz), stat=istat )
1142 astat = astat + istat
1143 allocate( o3_ph(nz), stat=istat )
1144 astat = astat + istat
1145 allocate( acb1_ph(nz), acb2_ph(nz), aoc1_ph(nz), aoc2_ph(nz), stat=istat )
1146 astat = astat + istat
1147 allocate( aant_ph(nz), aso4_ph(nz), asal_ph(nz), stat=istat )
1148 astat = astat + istat
1150 ! rajesh: allocate memory space to aerosol optical property variables
1151 allocate( tauaer300_ph(nz-1), tauaer400_ph(nz-1), tauaer600_ph(nz-1), &
1152 tauaer999_ph(nz-1), stat=istat)
1153 astat = astat + istat
1154 allocate( waer300_ph(nz-1), waer400_ph(nz-1), waer600_ph(nz-1), &
1155 waer999_ph(nz-1), stat=istat)
1156 astat = astat + istat
1157 allocate( gaer300_ph(nz-1), gaer400_ph(nz-1), gaer600_ph(nz-1), &
1158 gaer999_ph(nz-1), stat=istat)
1159 astat = astat + istat
1161 if( astat /= 0 ) then
1162 call wrf_message( 'ftuv_driver: failed to allocate _ph arrays' )
1163 call wrf_abort
1164 endif
1166 rh_ph(:) = 0.0_dp
1167 xlwc_ph(:) = 0.0_dp
1168 acb1_ph(:) = 0.0_dp
1169 acb2_ph(:) = 0.0_dp
1170 aoc1_ph(:) = 0.0_dp
1171 aoc2_ph(:) = 0.0_dp
1172 aant_ph(:) = 0.0_dp
1173 aso4_ph(:) = 0.0_dp
1174 asal_ph(:) = 0.0_dp
1176 ! rajesh: Initialize variables to zero
1177 tauaer300_ph(:) = 0.0_dp
1178 tauaer400_ph(:) = 0.0_dp
1179 tauaer600_ph(:) = 0.0_dp
1180 tauaer999_ph(:) = 0.0_dp
1181 waer300_ph(:) = 0.0_dp
1182 waer400_ph(:) = 0.0_dp
1183 waer600_ph(:) = 0.0_dp
1184 waer999_ph(:) = 0.0_dp
1185 gaer300_ph(:) = 0.0_dp
1186 gaer400_ph(:) = 0.0_dp
1187 gaer600_ph(:) = 0.0_dp
1188 gaer999_ph(:) = 0.0_dp
1190 !----------------------------------------------------
1191 ! assgin vertical profiles
1192 ! model levels first
1193 !----------------------------------------------------
1194 z_ph(1:nlev) = zin(1:nlev)
1195 tlev_ph(1:nlev) = tlevin(1:nlev)
1196 airlev_ph(1:nlev) = airlevin(1:nlev)
1197 rh_ph(1:nlev) = rhin(1:nlev)
1198 xlwc_ph(1:nlev) = xlwcin(1:nlev)
1199 o3_ph(1:nlev) = o3in(1:nlev)
1200 acb1_ph(1:nlev) = acb1in(1:nlev)
1201 acb2_ph(1:nlev) = acb2in(1:nlev)
1202 aoc1_ph(1:nlev) = aoc1in(1:nlev)
1203 aoc2_ph(1:nlev) = aoc2in(1:nlev)
1204 aant_ph(1:nlev) = aantin(1:nlev)
1205 aso4_ph(1:nlev) = aso4in(1:nlev)
1206 asal_ph(1:nlev) = asalin(1:nlev)
1208 ! rajesh:
1209 tauaer300_ph(1:nlev-1) = tauaer300in(1:nlev-1)
1210 tauaer400_ph(1:nlev-1) = tauaer400in(1:nlev-1)
1211 tauaer600_ph(1:nlev-1) = tauaer600in(1:nlev-1)
1212 tauaer999_ph(1:nlev-1) = tauaer999in(1:nlev-1)
1213 waer300_ph(1:nlev-1) = waer300in(1:nlev-1)
1214 waer400_ph(1:nlev-1) = waer400in(1:nlev-1)
1215 waer600_ph(1:nlev-1) = waer600in(1:nlev-1)
1216 waer999_ph(1:nlev-1) = waer999in(1:nlev-1)
1217 gaer300_ph(1:nlev-1) = gaer300in(1:nlev-1)
1218 gaer400_ph(1:nlev-1) = gaer400in(1:nlev-1)
1219 gaer600_ph(1:nlev-1) = gaer600in(1:nlev-1)
1220 gaer999_ph(1:nlev-1) = gaer999in(1:nlev-1)
1222 tlay_ph(1:nlev-1) = 0.5_dp*(tlev_ph(1:nlev-1) + tlev_ph(2:nlev))
1224 !----------------------------------------------------
1225 ! exo model levels
1226 !----------------------------------------------------
1227 if( nz > nlev ) then
1228 z_ph(nlev+1:nz-1) = zref(kcon:nref:2)
1229 tlev_ph(nlev+1:nz-1) = tref(kcon:nref:2)
1230 airlev_ph(nlev+1:nz-1) = airref(kcon:nref:2)
1231 o3_ph(nlev+1:nz-1) = o3ref(kcon:nref:2)/airref(kcon:nref:2) ! Changed to vmr
1233 z_ph(nz) = zref(nref)
1234 tlev_ph(nz) = tref(nref)
1235 airlev_ph(nz) = airref(nref)
1236 o3_ph(nz) = o3ref(nref)/airref(nref)
1237 tlay_ph(nlev) = 0.5_dp*(tlev_ph(nlev) + tref(kcon))
1238 tlay_ph(nlev+1:) = 0.5_dp*( tlev_ph(nlev+1:nz-1) + tlev_ph(nlev+2:) )
1239 end if
1242 call photoin( chem_opt, nz, zen, o3toms, esfact, &
1243 o3top, o2top, albedo(:,lu), z_ph, tlev_ph, &
1244 tlay_ph, airlev_ph, rh_ph, xlwc_ph, o3_ph, &
1245 acb1_ph, acb2_ph, aoc1_ph, aoc2_ph, aant_ph, &
1246 aso4_ph, asal_ph, &
1247 tauaer300_ph, tauaer400_ph, tauaer600_ph, &
1248 tauaer999_ph, waer300_ph, waer400_ph, &
1249 waer600_ph, waer999_ph, gaer300_ph, &
1250 gaer400_ph, gaer600_ph, gaer999_ph, &
1251 aer_ra_feedback, ftuv, adjcoe, radfld )
1253 do n = 1,njout
1254 prate(1:nlev-1,n) = ftuv(2:nlev,n)
1255 prate0(1:nz,n) = ftuv(1:nz,n)
1256 enddo
1258 deallocate( ftuv, z_ph, tlev_ph, tlay_ph, airlev_ph, rh_ph, xlwc_ph, &
1259 o3_ph, acb1_ph, acb2_ph, aoc1_ph, aoc2_ph, aant_ph, aso4_ph, &
1260 asal_ph, tauaer300_ph, tauaer400_ph,tauaer600_ph, tauaer999_ph, &
1261 waer300_ph, waer400_ph, waer600_ph, waer999_ph, &
1262 gaer300_ph, gaer400_ph, gaer600_ph, gaer999_ph )
1264 end subroutine photo
1266 subroutine ftuv_timestep_init( id, julday )
1267 !-----------------------------------------------------------------------------
1268 ! ... setup the time interpolation
1269 !-----------------------------------------------------------------------------
1271 use module_ftuv_subs, only : sundis
1273 implicit none
1275 !-----------------------------------------------------------------------------
1276 ! ... dummy arguments
1277 !-----------------------------------------------------------------------------
1278 integer, intent(in) :: id ! domain id
1279 integer, intent(in) :: julday ! day of year at present time step
1281 !-----------------------------------------------------------------------------
1282 ! ... local variables
1283 !-----------------------------------------------------------------------------
1284 integer :: m
1285 real(dp) :: calday
1287 calday = real( julday,kind=dp)
1288 if( calday < col_dens(id)%day_of_year(1) ) then
1289 next = 1
1290 last = 12
1291 dels = (365._dp + calday - col_dens(id)%day_of_year(12)) &
1292 / (365._dp + col_dens(id)%day_of_year(1) - col_dens(id)%day_of_year(12))
1293 else if( calday >= col_dens(id)%day_of_year(12) ) then
1294 next = 1
1295 last = 12
1296 dels = (calday - col_dens(id)%day_of_year(12)) &
1297 / (365. + col_dens(id)%day_of_year(1) - col_dens(id)%day_of_year(12))
1298 else
1299 do m = 11,1,-1
1300 if( calday >= col_dens(id)%day_of_year(m) ) then
1301 exit
1302 end if
1303 end do
1304 last = m
1305 next = m + 1
1306 dels = (calday - col_dens(id)%day_of_year(m)) / (col_dens(id)%day_of_year(m+1) - col_dens(id)%day_of_year(m))
1307 end if
1308 !-----------------------------------------------------------------------------
1309 ! set solar distance factor
1310 !-----------------------------------------------------------------------------
1311 if( curjulday /= julday ) then
1312 curjulday = julday
1313 call sundis( curjulday, esfact )
1314 endif
1316 end subroutine ftuv_timestep_init
1318 subroutine p_interp( o2_exo_col, o3_exo_col, ptop, &
1319 id, its, ite, jts, jte )
1320 !---------------------------------------------------------------
1321 ! ... pressure interpolation for exo col density
1322 !---------------------------------------------------------------
1324 implicit none
1326 !---------------------------------------------------------------
1327 ! ... dummy arguments
1328 !---------------------------------------------------------------
1329 integer, intent(in) :: id
1330 integer, intent(in) :: its, ite
1331 integer, intent(in) :: jts, jte
1332 real(dp), intent(in) :: ptop(its:ite,jts:jte) ! pressure at photolysis top (hPa)
1333 real(dp), intent(out) :: o2_exo_col(its:ite,jts:jte) ! exo model o2 column density (molecules/cm^2)
1334 real(dp), intent(out) :: o3_exo_col(its:ite,jts:jte) ! exo model o3 column density (molecules/cm^2)
1336 !---------------------------------------------------------------
1337 ! ... local variables
1338 !---------------------------------------------------------------
1339 integer :: i, j, k, ku, kl
1340 real(dp) :: pinterp
1341 real(dp) :: delp
1342 real(dp) :: tint_vals(2)
1344 lat_loop : &
1345 do j = jts,jte
1346 long_loop : &
1347 do i = its,ite
1348 pinterp = ptop(i,j)
1349 if( pinterp < col_dens(id)%col_levs(1) ) then
1350 ku = 1
1351 kl = 1
1352 delp = 0._dp
1353 else
1354 do ku = 2,col_dens(id)%ncoldens_levs
1355 if( pinterp <= col_dens(id)%col_levs(ku) ) then
1356 kl = ku - 1
1357 delp = log( pinterp/col_dens(id)%col_levs(kl) )/log( col_dens(id)%col_levs(ku)/col_dens(id)%col_levs(kl) )
1358 exit
1359 end if
1360 end do
1361 end if
1362 tint_vals(1) = col_dens(id)%o2_col_dens(i,j,kl,last) &
1363 + delp * (col_dens(id)%o2_col_dens(i,j,ku,last) &
1364 - col_dens(id)%o2_col_dens(i,j,kl,last))
1365 tint_vals(2) = col_dens(id)%o2_col_dens(i,j,kl,next) &
1366 + delp * (col_dens(id)%o2_col_dens(i,j,ku,next) &
1367 - col_dens(id)%o2_col_dens(i,j,kl,next))
1368 o2_exo_col(i,j) = tint_vals(1) + dels * (tint_vals(2) - tint_vals(1))
1369 tint_vals(1) = col_dens(id)%o3_col_dens(i,j,kl,last) &
1370 + delp * (col_dens(id)%o3_col_dens(i,j,ku,last) &
1371 - col_dens(id)%o3_col_dens(i,j,kl,last))
1372 tint_vals(2) = col_dens(id)%o3_col_dens(i,j,kl,next) &
1373 + delp * (col_dens(id)%o3_col_dens(i,j,ku,next) &
1374 - col_dens(id)%o3_col_dens(i,j,kl,next))
1375 o3_exo_col(i,j) = tint_vals(1) + dels * (tint_vals(2) - tint_vals(1))
1376 end do long_loop
1377 end do lat_loop
1379 end subroutine p_interp
1381 end module module_ftuv_driver