Context Navigation

mo_simple_plumes.F90 @ 5077

Last change on this file since 5077 was 5075, checked in by abarral, 5 months ago

[continued & end] replace netcdf by lmdz_netcdf.F90 wrapper
"use netcdf" is now only used in lmdz_netcdf.F90 (except ecrad and obsolete/)
<include "netcdf.inc"> is now likewise only used in lmdz_netcdf.F90.

systematically specify explicitely <USE lmdz_netcdf, ONLY:> (probably left some missing, to correct later on)

Further replacement of nf_put_* by nf90_put_* (same for _get_)

[minor] replace deprecated boolean operators along the way

File size: 22.6 KB

Line
1	!>
2	!!
3	!! @brief Module MO_SIMPLE_PLUMES: provides anthropogenic aerosol optical properties as a function of lat, lon
4	!! height, time, and wavelength
5	!!
6	!! @remarks
7	!!
8	!! @author Bjorn Stevens, Stephanie Fiedler and Karsten Peters MPI-Met, Hamburg (v1 release 2016-11-10)
9	!!
10	!! @change-log:
11	!! - 2016-12-05: beta release (BS, SF and KP, MPI-Met)
12	!! - 2016-09-28: revised representation of Twomey effect (SF, MPI-Met)
13	!! - 2015-09-28: bug fixes (SF, MPI-Met)
14	!! - 2016-10-12: revised maximum longitudinal extent of European plume (KP, SF, MPI-Met)
15	!! $ID: n/a$
16	!!
17	!! @par Origin
18	!! Based on code originally developed at the MPI-Met by Karsten Peters, Bjorn Stevens, Stephanie Fiedler
19	!! and Stefan Kinne with input from Thorsten Mauritsen and Robert Pincus
20	!!
21	!! @par Copyright
22	!!
23	!
24	MODULE MO_SIMPLE_PLUMES
25
26	USE lmdz_netcdf, ONLY:nf90_get_var,nf90_close,nf90_inq_varid,nf90_inq_dimid,&
27	nf90_inquire_dimension,nf90_noerr,nf90_nowrite,nf90_open
28
29	IMPLICIT NONE
30
31	INTEGER, PARAMETER :: &
32	nplumes = 9 ,& !< Number of plumes
33	nfeatures = 2 ,& !< Number of features per plume
34	ntimes = 52 ,& !< Number of times resolved per year (52 => weekly resolution)
35	nyears = 251 !< Number of years of available forcing
36
37	LOGICAL, SAVE :: &
38	sp_initialized = .FALSE. !< parameter determining whether input needs to be read
39	!$OMP THREADPRIVATE(sp_initialized)
40
41	REAL, SAVE :: &
42	plume_lat (nplumes) ,& !< latitude of plume center (AOD maximum)
43	plume_lon (nplumes) ,& !< longitude of plume center (AOD maximum)
44	beta_a (nplumes) ,& !< parameter a for beta function vertical profile
45	beta_b (nplumes) ,& !< parameter b for beta function vertical profile
46	aod_spmx (nplumes) ,& !< anthropogenic AOD maximum at 550 for plumes
47	aod_fmbg (nplumes) ,& !< anthropogenic AOD at 550 for fine-mode natural background (idealized to mimic Twomey effect)
48	asy550 (nplumes) ,& !< asymmetry parameter at 550nm for plume
49	ssa550 (nplumes) ,& !< single scattering albedo at 550nm for plume
50	angstrom (nplumes) ,& !< Angstrom parameter for plume
51	sig_lon_E (nfeatures,nplumes) ,& !< Eastward extent of plume feature
52	sig_lon_W (nfeatures,nplumes) ,& !< Westward extent of plume feature
53	sig_lat_E (nfeatures,nplumes) ,& !< Southward extent of plume feature
54	sig_lat_W (nfeatures,nplumes) ,& !< Northward extent of plume feature
55	theta (nfeatures,nplumes) ,& !< Rotation angle of plume feature
56	ftr_weight (nfeatures,nplumes) ,& !< Feature weights
57	year_weight (nyears,nplumes) ,& !< Yearly weight for plume
58	ann_cycle (nfeatures,ntimes,nplumes) !< annual cycle for plume feature
59	!$OMP THREADPRIVATE(plume_lat,plume_lon,beta_a,beta_b,aod_spmx,aod_fmbg,asy550,ssa550,angstrom)
60	!$OMP THREADPRIVATE(sig_lon_E,sig_lon_W,sig_lat_E,sig_lat_W,theta,ftr_weight,year_weight,ann_cycle)
61
62	REAL :: &
63	time_weight (nfeatures,nplumes) ,& !< Time weights
64	time_weight_bg (nfeatures,nplumes) !< as time_weight but for natural background in Twomey effect
65
66	PUBLIC sp_aop_profile
67
68	CONTAINS
69	!
70	! ------------------------------------------------------------------------------------------------------------------------
71	! SP_SETUP: This subroutine should be called at initialization to read the netcdf data that describes the simple plume
72	! climatology. The information needs to be either read by each processor or distributed to processors.
73	!
74	SUBROUTINE sp_setup
75	!
76	USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
77	USE mod_phys_lmdz_omp_data, ONLY: is_omp_root
78	USE mod_phys_lmdz_transfert_para, ONLY: bcast
79	!
80	! ----------
81	!
82	INTEGER :: iret, ncid, DimID, VarID, xdmy
83	CHARACTER (len = 50) :: modname = 'mo_simple_plumes.sp_setup'
84	CHARACTER (len = 80) :: abort_message
85	!
86	! ----------
87	!--only one processor reads the input data
88	IF (is_mpi_root.AND.is_omp_root) THEN
89	!
90	iret = nf90_open("MACv2.0-SP_v1.nc", NF90_NOWRITE, ncid)
91	IF (iret /= NF90_NOERR) THEN
92	abort_message='NetCDF File not opened'
93	CALL abort_physic(modname,abort_message,1)
94	ENDIF
95	!
96	! read dimensions and make sure file conforms to expected size
97	!
98	iret = nf90_inq_dimid(ncid, "plume_number" , DimId)
99	iret = nf90_inquire_dimension(ncid, DimId, len = xdmy)
100	IF (xdmy /= nplumes) THEN
101	abort_message='NetCDF improperly dimensioned -- plume_number'
102	CALL abort_physic(modname,abort_message,1)
103	ENDIF
104	!
105	iret = nf90_inq_dimid(ncid, "plume_feature", DimId)
106	iret = nf90_inquire_dimension(ncid, DimId, len = xdmy)
107	IF (xdmy /= nfeatures) THEN
108	abort_message='NetCDF improperly dimensioned -- plume_feature'
109	CALL abort_physic(modname,abort_message,1)
110	ENDIF
111	!
112	iret = nf90_inq_dimid(ncid, "year_fr" , DimId)
113	iret = nf90_inquire_dimension(ncid, DimID, len = xdmy)
114	IF (xdmy /= ntimes) THEN
115	abort_message='NetCDF improperly dimensioned -- year_fr'
116	CALL abort_physic(modname,abort_message,1)
117	ENDIF
118	!
119	iret = nf90_inq_dimid(ncid, "years" , DimId)
120	iret = nf90_inquire_dimension(ncid, DimID, len = xdmy)
121	IF (xdmy /= nyears) THEN
122	abort_message='NetCDF improperly dimensioned -- years'
123	CALL abort_physic(modname,abort_message,1)
124	ENDIF
125	!
126	! read variables that define the simple plume climatology
127	!
128	iret = nf90_inq_varid(ncid, "plume_lat", VarId)
129	iret = nf90_get_var(ncid, VarID, plume_lat(:), start=(/1/),count=(/nplumes/))
130	IF (iret /= NF90_NOERR) THEN
131	abort_message='NetCDF Error reading plume_lat'
132	CALL abort_physic(modname,abort_message,1)
133	ENDIF
134	!
135	iret = nf90_inq_varid(ncid, "plume_lon", VarId)
136	iret = nf90_get_var(ncid, VarID, plume_lon(:), start=(/1/),count=(/nplumes/))
137	IF (iret /= NF90_NOERR) THEN
138	abort_message='NetCDF Error reading plume_lon'
139	CALL abort_physic(modname,abort_message,1)
140	ENDIF
141	!
142	iret = nf90_inq_varid(ncid, "beta_a" , VarId)
143	iret = nf90_get_var(ncid, VarID, beta_a(:) , start=(/1/),count=(/nplumes/))
144	IF (iret /= NF90_NOERR) THEN
145	abort_message='NetCDF Error reading beta_a'
146	CALL abort_physic(modname,abort_message,1)
147	ENDIF
148	!
149	iret = nf90_inq_varid(ncid, "beta_b" , VarId)
150	iret = nf90_get_var(ncid, VarID, beta_b(:) , start=(/1/),count=(/nplumes/))
151	IF (iret /= NF90_NOERR) THEN
152	abort_message='NetCDF Error reading beta_b'
153	CALL abort_physic(modname,abort_message,1)
154	ENDIF
155	!
156	iret = nf90_inq_varid(ncid, "aod_spmx" , VarId)
157	iret = nf90_get_var(ncid, VarID, aod_spmx(:) , start=(/1/),count=(/nplumes/))
158	IF (iret /= NF90_NOERR) THEN
159	abort_message='NetCDF Error reading aod_spmx'
160	CALL abort_physic(modname,abort_message,1)
161	ENDIF
162	!
163	iret = nf90_inq_varid(ncid, "aod_fmbg" , VarId)
164	iret = nf90_get_var(ncid, VarID, aod_fmbg(:) , start=(/1/),count=(/nplumes/))
165	IF (iret /= NF90_NOERR) THEN
166	abort_message='NetCDF Error reading aod_fmbg'
167	CALL abort_physic(modname,abort_message,1)
168	ENDIF
169	!
170	iret = nf90_inq_varid(ncid, "ssa550" , VarId)
171	iret = nf90_get_var(ncid, VarID, ssa550(:) , start=(/1/),count=(/nplumes/))
172	IF (iret /= NF90_NOERR) THEN
173	abort_message='NetCDF Error reading ssa550'
174	CALL abort_physic(modname,abort_message,1)
175	ENDIF
176	!
177	iret = nf90_inq_varid(ncid, "asy550" , VarId)
178	iret = nf90_get_var(ncid, VarID, asy550(:) , start=(/1/),count=(/nplumes/))
179	IF (iret /= NF90_NOERR) THEN
180	abort_message='NetCDF Error reading asy550'
181	CALL abort_physic(modname,abort_message,1)
182	ENDIF
183	!
184	iret = nf90_inq_varid(ncid, "angstrom" , VarId)
185	iret = nf90_get_var(ncid, VarID, angstrom(:), start=(/1/),count=(/nplumes/))
186	IF (iret /= NF90_NOERR) THEN
187	abort_message='NetCDF Error reading angstrom'
188	CALL abort_physic(modname,abort_message,1)
189	ENDIF
190	!
191	iret = nf90_inq_varid(ncid, "sig_lat_W" , VarId)
192	iret = nf90_get_var(ncid, VarID, sig_lat_W(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
193	IF (iret /= NF90_NOERR) THEN
194	abort_message='NetCDF Error reading sig_lat_W'
195	CALL abort_physic(modname,abort_message,1)
196	ENDIF
197	!
198	iret = nf90_inq_varid(ncid, "sig_lat_E" , VarId)
199	iret = nf90_get_var(ncid, VarID, sig_lat_E(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
200	IF (iret /= NF90_NOERR) THEN
201	abort_message='NetCDF Error reading sig_lat_E'
202	CALL abort_physic(modname,abort_message,1)
203	ENDIF
204	!
205	iret = nf90_inq_varid(ncid, "sig_lon_E" , VarId)
206	iret = nf90_get_var(ncid, VarID, sig_lon_E(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
207	IF (iret /= NF90_NOERR) THEN
208	abort_message='NetCDF Error reading sig_lon_E'
209	CALL abort_physic(modname,abort_message,1)
210	ENDIF
211	!
212	iret = nf90_inq_varid(ncid, "sig_lon_W" , VarId)
213	iret = nf90_get_var(ncid, VarID, sig_lon_W(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
214	IF (iret /= NF90_NOERR) THEN
215	abort_message='NetCDF Error reading sig_lon_W'
216	CALL abort_physic(modname,abort_message,1)
217	ENDIF
218	!
219	iret = nf90_inq_varid(ncid, "theta" , VarId)
220	iret = nf90_get_var(ncid, VarID, theta(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
221	IF (iret /= NF90_NOERR) THEN
222	abort_message='NetCDF Error reading theta'
223	CALL abort_physic(modname,abort_message,1)
224	ENDIF
225	!
226	iret = nf90_inq_varid(ncid, "ftr_weight" , VarId)
227	iret = nf90_get_var(ncid, VarID, ftr_weight(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
228	IF (iret /= NF90_NOERR) THEN
229	abort_message='NetCDF Error reading plume_lat'
230	CALL abort_physic(modname,abort_message,1)
231	ENDIF
232	!
233	iret = nf90_inq_varid(ncid, "year_weight" , VarId)
234	iret = nf90_get_var(ncid, VarID, year_weight(:,:) , start=(/1,1/),count=(/nyears,nplumes /))
235	IF (iret /= NF90_NOERR) THEN
236	abort_message='NetCDF Error reading year_weight'
237	CALL abort_physic(modname,abort_message,1)
238	ENDIF
239	!
240	iret = nf90_inq_varid(ncid, "ann_cycle" , VarId)
241	iret = nf90_get_var(ncid, VarID, ann_cycle(:,:,:) , start=(/1,1,1/),count=(/nfeatures,ntimes,nplumes/))
242	IF (iret /= NF90_NOERR) THEN
243	abort_message='NetCDF Error reading ann_cycle'
244	CALL abort_physic(modname,abort_message,1)
245	ENDIF
246	!
247	iret = nf90_close(ncid)
248	!
249	ENDIF !--root processor
250	!
251	CALL bcast(plume_lat)
252	CALL bcast(plume_lon)
253	CALL bcast(beta_a)
254	CALL bcast(beta_b)
255	CALL bcast(aod_spmx)
256	CALL bcast(aod_fmbg)
257	CALL bcast(asy550)
258	CALL bcast(ssa550)
259	CALL bcast(angstrom)
260	CALL bcast(sig_lon_E)
261	CALL bcast(sig_lon_W)
262	CALL bcast(sig_lat_E)
263	CALL bcast(sig_lat_W)
264	CALL bcast(theta)
265	CALL bcast(ftr_weight)
266	CALL bcast(year_weight)
267	CALL bcast(ann_cycle)
268	!
269	sp_initialized = .TRUE.
270	!
271	RETURN
272	!
273	END SUBROUTINE sp_setup
274	!
275	! ------------------------------------------------------------------------------------------------------------------------
276	! SET_TIME_WEIGHT: The simple plume model assumes that meteorology constrains plume shape and that only source strength
277	! influences the amplitude of a plume associated with a given source region. This routine retrieves the temporal weights
278	! for the plumes. Each plume feature has its own temporal weights which varies yearly. The annual cycle is indexed by
279	! week in the year and superimposed on the yearly mean value of the weight.
280	!
281	SUBROUTINE set_time_weight(year_fr)
282	!
283	! ----------
284	!
285	REAL, INTENT(IN) :: &
286	year_fr !< Fractional Year (1850.0 - 2100.99)
287
288	INTEGER :: &
289	iyear ,& !< Integer year values between 1 and 156 (1850-2100)
290	iweek ,& !< Integer index (between 1 and ntimes); for ntimes=52 this corresponds to weeks (roughly)
291	iplume ! plume number
292	!
293	! ----------
294	!
295	iyear = FLOOR(year_fr) - 1849
296	iweek = FLOOR((year_fr - FLOOR(year_fr)) * ntimes) + 1
297
298	IF ((iweek > ntimes) .OR. (iweek < 1) .OR. (iyear > nyears) .OR. (iyear < 1)) THEN
299	CALL abort_physic('set_time_weight','Time out of bounds',1)
300	ENDIF
301
302	DO iplume=1,nplumes
303	time_weight(1,iplume) = year_weight(iyear,iplume) * ann_cycle(1,iweek,iplume)
304	time_weight(2,iplume) = year_weight(iyear,iplume) * ann_cycle(2,iweek,iplume)
305	time_weight_bg(1,iplume) = ann_cycle(1,iweek,iplume)
306	time_weight_bg(2,iplume) = ann_cycle(2,iweek,iplume)
307	ENDDO
308
309	RETURN
310	END SUBROUTINE set_time_weight
311	!
312	! ------------------------------------------------------------------------------------------------------------------------
313	! SP_AOP_PROFILE: This subroutine calculates the simple plume aerosol and cloud active optical properties based on the
314	! the simple plume fit to the MPI Aerosol Climatology (Version 2). It sums over nplumes to provide a profile of aerosol
315	! optical properties on a host models vertical grid.
316	!
317	SUBROUTINE sp_aop_profile ( &
318	nlevels ,ncol ,lambda ,oro ,lon ,lat , &
319	year_fr ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , &
320	asy_prof )
321	!
322	! ----------
323	!
324	INTEGER, INTENT(IN) :: &
325	nlevels, & !< number of levels
326	ncol !< number of columns
327
328	REAL, INTENT(IN) :: &
329	lambda, & !< wavelength
330	year_fr, & !< Fractional Year (1903.0 is the 0Z on the first of January 1903, Gregorian)
331	oro(ncol), & !< orographic height (m)
332	lon(ncol), & !< longitude
333	lat(ncol), & !< latitude
334	z (ncol,nlevels), & !< height above sea-level (m)
335	dz(ncol,nlevels) !< level thickness (difference between half levels) (m)
336
337	REAL, INTENT(OUT) :: &
338	dNovrN(ncol) , & !< anthropogenic increase in cloud drop number concentration (factor)
339	aod_prof(ncol,nlevels) , & !< profile of aerosol optical depth
340	ssa_prof(ncol,nlevels) , & !< profile of single scattering albedo
341	asy_prof(ncol,nlevels) !< profile of asymmetry parameter
342
343	INTEGER :: iplume, icol, k
344
345	REAL :: &
346	eta(ncol,nlevels), & !< normalized height (by 15 km)
347	z_beta(ncol,nlevels), & !< profile for scaling column optical depth
348	prof(ncol,nlevels), & !< scaled profile (by beta function)
349	beta_sum(ncol), & !< vertical sum of beta function
350	ssa(ncol), & !< single scattering albedo
351	asy(ncol), & !< asymmetry parameter
352	cw_an(ncol), & !< column weight for simple plume (anthropogenic) AOD at 550 nm
353	cw_bg(ncol), & !< column weight for fine-mode natural background AOD at 550 nm
354	caod_sp(ncol), & !< column simple plume anthropogenic AOD at 550 nm
355	caod_bg(ncol), & !< column fine-mode natural background AOD at 550 nm
356	a_plume1, & !< gaussian longitude factor for feature 1
357	a_plume2, & !< gaussian longitude factor for feature 2
358	b_plume1, & !< gaussian latitude factor for feature 1
359	b_plume2, & !< gaussian latitude factor for feature 2
360	delta_lat, & !< latitude offset
361	delta_lon, & !< longitude offset
362	delta_lon_t, & !< threshold for maximum longitudinal plume extent used in transition from 360 to 0 degrees
363	lon1, & !< rotated longitude for feature 1
364	lat1, & !< rotated latitude for feature 2
365	lon2, & !< rotated longitude for feature 1
366	lat2, & !< rotated latitude for feature 2
367	f1, & !< contribution from feature 1
368	f2, & !< contribution from feature 2
369	f3, & !< contribution from feature 1 in natural background of Twomey effect
370	f4, & !< contribution from feature 2 in natural background of Twomey effect
371	aod_550, & !< aerosol optical depth at 550nm
372	aod_lmd, & !< aerosol optical depth at input wavelength
373	lfactor !< factor to compute wavelength dependence of optical properties
374	!
375	! ----------
376	!
377	! initialize input data (by calling setup at first instance)
378	!
379	IF (.NOT.sp_initialized) CALL sp_setup
380	!
381	! get time weights
382	!
383	CALL set_time_weight(year_fr)
384	!
385	! initialize variables, including output
386	!
387	DO k=1,nlevels
388	DO icol=1,ncol
389	aod_prof(icol,k) = 0.0
390	ssa_prof(icol,k) = 0.0
391	asy_prof(icol,k) = 0.0
392	z_beta(icol,k) = MERGE(1.0, 0.0, z(icol,k) >= oro(icol))
393	eta(icol,k) = MAX(0.0,MIN(1.0,z(icol,k)/15000.))
394	ENDDO
395	ENDDO
396	DO icol=1,ncol
397	dNovrN(icol) = 1.0
398	caod_sp(icol) = 0.0
399	caod_bg(icol) = 0.02
400	ENDDO
401	!
402	! sum contribution from plumes to construct composite profiles of aerosol optical properties
403	!
404	DO iplume=1,nplumes
405	!
406	! calculate vertical distribution function from parameters of beta distribution
407	!
408	DO icol=1,ncol
409	beta_sum(icol) = 0.
410	ENDDO
411	DO k=1,nlevels
412	DO icol=1,ncol
413	prof(icol,k) = (eta(icol,k)*(beta_a(iplume)-1.) (1.-eta(icol,k))*(beta_b(iplume)-1.)) dz(icol,k)
414	beta_sum(icol) = beta_sum(icol) + prof(icol,k)
415	ENDDO
416	ENDDO
417	DO k=1,nlevels
418	DO icol=1,ncol
419	prof(icol,k) = ( prof(icol,k) / beta_sum(icol) ) * z_beta(icol,k)
420	ENDDO
421	ENDDO
422	!
423	! calculate plume weights
424	!
425	DO icol=1,ncol
426	!
427	! get plume-center relative spatial parameters for specifying amplitude of plume at given lat and lon
428	!
429	delta_lat = lat(icol) - plume_lat(iplume)
430	delta_lon = lon(icol) - plume_lon(iplume)
431	delta_lon_t = MERGE (260., 180., iplume == 1)
432	delta_lon = MERGE ( delta_lon-SIGN(360.,delta_lon) , delta_lon , ABS(delta_lon) > delta_lon_t)
433
434	a_plume1 = 0.5 / (MERGE(sig_lon_E(1,iplume), sig_lon_W(1,iplume), delta_lon > 0)**2)
435	b_plume1 = 0.5 / (MERGE(sig_lat_E(1,iplume), sig_lat_W(1,iplume), delta_lon > 0)**2)
436	a_plume2 = 0.5 / (MERGE(sig_lon_E(2,iplume), sig_lon_W(2,iplume), delta_lon > 0)**2)
437	b_plume2 = 0.5 / (MERGE(sig_lat_E(2,iplume), sig_lat_W(2,iplume), delta_lon > 0)**2)
438	!
439	! adjust for a plume specific rotation which helps match plume state to climatology.
440	!
441	lon1 = COS(theta(1,iplume))(delta_lon) + SIN(theta(1,iplume))(delta_lat)
442	lat1 = - SIN(theta(1,iplume))(delta_lon) + COS(theta(1,iplume))(delta_lat)
443	lon2 = COS(theta(2,iplume))(delta_lon) + SIN(theta(2,iplume))(delta_lat)
444	lat2 = - SIN(theta(2,iplume))(delta_lon) + COS(theta(2,iplume))(delta_lat)
445	!
446	! calculate contribution to plume from its different features, to get a column weight for the anthropogenic
447	! (cw_an) and the fine-mode natural background aerosol (cw_bg)
448	!
449	f1 = time_weight(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)*2) + (b_plume1 ((lat1)**2))))
450	f2 = time_weight(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)*2) + (b_plume2 ((lat2)**2))))
451	f3 = time_weight_bg(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)*2) + (b_plume1 ((lat1)**2))))
452	f4 = time_weight_bg(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)*2) + (b_plume2 ((lat2)**2))))
453
454	cw_an(icol) = f1 * aod_spmx(iplume) + f2 * aod_spmx(iplume)
455	cw_bg(icol) = f3 * aod_fmbg(iplume) + f4 * aod_fmbg(iplume)
456	!
457	! calculate wavelength-dependent scattering properties
458	!
459	lfactor = MIN(1.0,700.0/lambda)
460	ssa(icol) = (ssa550(iplume) * lfactor*4) / ((ssa550(iplume) lfactor*4) + ((1-ssa550(iplume)) lfactor))
461	asy(icol) = asy550(iplume) * SQRT(lfactor)
462	ENDDO
463	!
464	! distribute plume optical properties across its vertical profile weighting by optical depth and scaling for
465	! wavelength using the angstrom parameter.
466	!
467	lfactor = EXP(-angstrom(iplume) * LOG(lambda/550.0))
468	DO k=1,nlevels
469	DO icol = 1,ncol
470	aod_550 = prof(icol,k) * cw_an(icol)
471	aod_lmd = aod_550 * lfactor
472	caod_sp(icol) = caod_sp(icol) + aod_550
473	caod_bg(icol) = caod_bg(icol) + prof(icol,k) * cw_bg(icol)
474	asy_prof(icol,k) = asy_prof(icol,k) + aod_lmd * ssa(icol) * asy(icol)
475	ssa_prof(icol,k) = ssa_prof(icol,k) + aod_lmd * ssa(icol)
476	aod_prof(icol,k) = aod_prof(icol,k) + aod_lmd
477	ENDDO
478	ENDDO
479	ENDDO
480	!
481	! complete optical depth weighting
482	!
483	DO k=1,nlevels
484	DO icol = 1,ncol
485	asy_prof(icol,k) = MERGE(asy_prof(icol,k)/ssa_prof(icol,k), 0.0, ssa_prof(icol,k) > TINY(1.))
486	ssa_prof(icol,k) = MERGE(ssa_prof(icol,k)/aod_prof(icol,k), 1.0, aod_prof(icol,k) > TINY(1.))
487	ENDDO
488	ENDDO
489	!
490	! calculate effective radius normalization (divisor) factor
491	!
492	DO icol=1,ncol
493	dNovrN(icol) = LOG((1000.0 * (caod_sp(icol) + caod_bg(icol))) + 1.0)/LOG((1000.0 * caod_bg(icol)) + 1.0)
494	ENDDO
495
496	RETURN
497	END SUBROUTINE sp_aop_profile
498
499	END MODULE MO_SIMPLE_PLUMES

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format