Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

mo_simple_plumes.F90 @ 4688

Last change on this file since 4688 was 3318, checked in by jghattas, 7 years ago
Integration des developpements fait sur la trunk par O. Boucher concernant le MACspV2 aerosol plume climatology(nouveau option flag_aersol=7). Les commits suivants fait sur le trunk sont ici merge : [3274], [3279], [3287], [3288], [3290], [3295], [3296], [3297]. Tout les modifications dans newmicro.f90 ne sont pas retenu mais les changemnets lie au flag_aerosol=7 sont prise.
File size: 20.2 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 netcdf
27
28	IMPLICIT NONE
29
30	INTEGER, PARAMETER :: &
31	nplumes = 9 ,& !< Number of plumes
32	nfeatures = 2 ,& !< Number of features per plume
33	ntimes = 52 ,& !< Number of times resolved per year (52 => weekly resolution)
34	nyears = 251 !< Number of years of available forcing
35
36	LOGICAL, SAVE :: &
37	sp_initialized = .FALSE. !< parameter determining whether input needs to be read
38	!$OMP THREADPRIVATE(sp_initialized)
39
40	REAL :: &
41	plume_lat (nplumes) ,& !< latitude of plume center (AOD maximum)
42	plume_lon (nplumes) ,& !< longitude of plume center (AOD maximum)
43	beta_a (nplumes) ,& !< parameter a for beta function vertical profile
44	beta_b (nplumes) ,& !< parameter b for beta function vertical profile
45	aod_spmx (nplumes) ,& !< anthropogenic AOD maximum at 550 for plumes
46	aod_fmbg (nplumes) ,& !< anthropogenic AOD at 550 for fine-mode natural background (idealized to mimic Twomey effect)
47	asy550 (nplumes) ,& !< asymmetry parameter at 550nm for plume
48	ssa550 (nplumes) ,& !< single scattering albedo at 550nm for plume
49	angstrom (nplumes) ,& !< Angstrom parameter for plume
50	sig_lon_E (nfeatures,nplumes) ,& !< Eastward extent of plume feature
51	sig_lon_W (nfeatures,nplumes) ,& !< Westward extent of plume feature
52	sig_lat_E (nfeatures,nplumes) ,& !< Southward extent of plume feature
53	sig_lat_W (nfeatures,nplumes) ,& !< Northward extent of plume feature
54	theta (nfeatures,nplumes) ,& !< Rotation angle of plume feature
55	ftr_weight (nfeatures,nplumes) ,& !< Feature weights
56	time_weight (nfeatures,nplumes) ,& !< Time weights
57	time_weight_bg (nfeatures,nplumes) ,& !< as time_weight but for natural background in Twomey effect
58	year_weight (nyears,nplumes) ,& !< Yearly weight for plume
59	ann_cycle (nfeatures,ntimes,nplumes) !< annual cycle for plume feature
60	!$OMP THREADPRIVATE(plume_lat,plume_lon,beta_a,beta_b,aod_spmx,aod_fmbg,asy550,ssa550,angstrom)
61	!$OMP THREADPRIVATE(sig_lon_E,sig_lon_W,sig_lat_E,sig_lat_W,theta,ftr_weight,year_weight,ann_cycle)
62
63	PUBLIC sp_aop_profile
64
65	CONTAINS
66	!
67	! ------------------------------------------------------------------------------------------------------------------------
68	! SP_SETUP: This subroutine should be called at initialization to read the netcdf data that describes the simple plume
69	! climatology. The information needs to be either read by each processor or distributed to processors.
70	!
71	SUBROUTINE sp_setup
72	!
73	USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
74	USE mod_phys_lmdz_omp_data, ONLY: is_omp_root
75	USE mod_phys_lmdz_transfert_para, ONLY: bcast
76	!
77	! ----------
78	!
79	INTEGER :: iret, ncid, DimID, VarID, xdmy
80	!
81	! ----------
82	!--only one processor reads the input data
83	IF (is_mpi_root.AND.is_omp_root) THEN
84	!
85	iret = nf90_open("MACv2.0-SP_v1.nc", NF90_NOWRITE, ncid)
86	IF (iret /= NF90_NOERR) STOP 'NetCDF File not opened'
87	!
88	! read dimensions and make sure file conforms to expected size
89	!
90	iret = nf90_inq_dimid(ncid, "plume_number" , DimId)
91	iret = nf90_inquire_dimension(ncid, DimId, len = xdmy)
92	IF (xdmy /= nplumes) STOP 'NetCDF improperly dimensioned -- plume_number'
93	!
94	iret = nf90_inq_dimid(ncid, "plume_feature", DimId)
95	iret = nf90_inquire_dimension(ncid, DimId, len = xdmy)
96	IF (xdmy /= nfeatures) STOP 'NetCDF improperly dimensioned -- plume_feature'
97	!
98	iret = nf90_inq_dimid(ncid, "year_fr" , DimId)
99	iret = nf90_inquire_dimension(ncid, DimID, len = xdmy)
100	IF (xdmy /= ntimes) STOP 'NetCDF improperly dimensioned -- year_fr'
101	!
102	iret = nf90_inq_dimid(ncid, "years" , DimId)
103	iret = nf90_inquire_dimension(ncid, DimID, len = xdmy)
104	IF (xdmy /= nyears) STOP 'NetCDF improperly dimensioned -- years'
105	!
106	! read variables that define the simple plume climatology
107	!
108	iret = nf90_inq_varid(ncid, "plume_lat", VarId)
109	iret = nf90_get_var(ncid, VarID, plume_lat(:), start=(/1/),count=(/nplumes/))
110	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lat'
111	iret = nf90_inq_varid(ncid, "plume_lon", VarId)
112	iret = nf90_get_var(ncid, VarID, plume_lon(:), start=(/1/),count=(/nplumes/))
113	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lon'
114	iret = nf90_inq_varid(ncid, "beta_a" , VarId)
115	iret = nf90_get_var(ncid, VarID, beta_a(:) , start=(/1/),count=(/nplumes/))
116	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading beta_a'
117	iret = nf90_inq_varid(ncid, "beta_b" , VarId)
118	iret = nf90_get_var(ncid, VarID, beta_b(:) , start=(/1/),count=(/nplumes/))
119	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading beta_b'
120	iret = nf90_inq_varid(ncid, "aod_spmx" , VarId)
121	iret = nf90_get_var(ncid, VarID, aod_spmx(:) , start=(/1/),count=(/nplumes/))
122	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading aod_spmx'
123	iret = nf90_inq_varid(ncid, "aod_fmbg" , VarId)
124	iret = nf90_get_var(ncid, VarID, aod_fmbg(:) , start=(/1/),count=(/nplumes/))
125	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading aod_fmbg'
126	iret = nf90_inq_varid(ncid, "ssa550" , VarId)
127	iret = nf90_get_var(ncid, VarID, ssa550(:) , start=(/1/),count=(/nplumes/))
128	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading ssa550'
129	iret = nf90_inq_varid(ncid, "asy550" , VarId)
130	iret = nf90_get_var(ncid, VarID, asy550(:) , start=(/1/),count=(/nplumes/))
131	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading asy550'
132	iret = nf90_inq_varid(ncid, "angstrom" , VarId)
133	iret = nf90_get_var(ncid, VarID, angstrom(:), start=(/1/),count=(/nplumes/))
134	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading angstrom'
135	!
136	iret = nf90_inq_varid(ncid, "sig_lat_W" , VarId)
137	iret = nf90_get_var(ncid, VarID, sig_lat_W(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
138	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lat_W'
139	iret = nf90_inq_varid(ncid, "sig_lat_E" , VarId)
140	iret = nf90_get_var(ncid, VarID, sig_lat_E(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
141	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lat_E'
142	iret = nf90_inq_varid(ncid, "sig_lon_E" , VarId)
143	iret = nf90_get_var(ncid, VarID, sig_lon_E(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
144	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lon_E'
145	iret = nf90_inq_varid(ncid, "sig_lon_W" , VarId)
146	iret = nf90_get_var(ncid, VarID, sig_lon_W(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
147	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lon_W'
148	iret = nf90_inq_varid(ncid, "theta" , VarId)
149	iret = nf90_get_var(ncid, VarID, theta(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
150	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading theta'
151	iret = nf90_inq_varid(ncid, "ftr_weight" , VarId)
152	iret = nf90_get_var(ncid, VarID, ftr_weight(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
153	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lat'
154	iret = nf90_inq_varid(ncid, "year_weight" , VarId)
155	iret = nf90_get_var(ncid, VarID, year_weight(:,:) , start=(/1,1/),count=(/nyears,nplumes /))
156	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading year_weight'
157	iret = nf90_inq_varid(ncid, "ann_cycle" , VarId)
158	iret = nf90_get_var(ncid, VarID, ann_cycle(:,:,:) , start=(/1,1,1/),count=(/nfeatures,ntimes,nplumes/))
159	IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading ann_cycle'
160	!
161	iret = nf90_close(ncid)
162	!
163	ENDIF !--root processor
164	!
165	CALL bcast(plume_lat)
166	CALL bcast(plume_lon)
167	CALL bcast(beta_a)
168	CALL bcast(beta_b)
169	CALL bcast(aod_spmx)
170	CALL bcast(aod_fmbg)
171	CALL bcast(asy550)
172	CALL bcast(ssa550)
173	CALL bcast(angstrom)
174	CALL bcast(sig_lon_E)
175	CALL bcast(sig_lon_W)
176	CALL bcast(sig_lat_E)
177	CALL bcast(sig_lat_W)
178	CALL bcast(theta)
179	CALL bcast(ftr_weight)
180	CALL bcast(year_weight)
181	CALL bcast(ann_cycle)
182	!
183	sp_initialized = .TRUE.
184	!
185	RETURN
186	!
187	END SUBROUTINE sp_setup
188	!
189	! ------------------------------------------------------------------------------------------------------------------------
190	! SET_TIME_WEIGHT: The simple plume model assumes that meteorology constrains plume shape and that only source strength
191	! influences the amplitude of a plume associated with a given source region. This routine retrieves the temporal weights
192	! for the plumes. Each plume feature has its own temporal weights which varies yearly. The annual cycle is indexed by
193	! week in the year and superimposed on the yearly mean value of the weight.
194	!
195	SUBROUTINE set_time_weight(year_fr)
196	!
197	! ----------
198	!
199	REAL, INTENT(IN) :: &
200	year_fr !< Fractional Year (1850.0 - 2100.99)
201
202	INTEGER :: &
203	iyear ,& !< Integer year values between 1 and 156 (1850-2100)
204	iweek ,& !< Integer index (between 1 and ntimes); for ntimes=52 this corresponds to weeks (roughly)
205	iplume ! plume number
206	!
207	! ----------
208	!
209	iyear = FLOOR(year_fr) - 1849
210	iweek = FLOOR((year_fr - FLOOR(year_fr)) * ntimes) + 1
211
212	IF ((iweek > ntimes) .OR. (iweek < 1) .OR. (iyear > nyears) .OR. (iyear < 1)) THEN
213	CALL abort_physic('set_time_weight','Time out of bounds')
214	ENDIF
215
216	DO iplume=1,nplumes
217	time_weight(1,iplume) = year_weight(iyear,iplume) * ann_cycle(1,iweek,iplume)
218	time_weight(2,iplume) = year_weight(iyear,iplume) * ann_cycle(2,iweek,iplume)
219	time_weight_bg(1,iplume) = ann_cycle(1,iweek,iplume)
220	time_weight_bg(2,iplume) = ann_cycle(2,iweek,iplume)
221	END DO
222
223	RETURN
224	END SUBROUTINE set_time_weight
225	!
226	! ------------------------------------------------------------------------------------------------------------------------
227	! SP_AOP_PROFILE: This subroutine calculates the simple plume aerosol and cloud active optical properties based on the
228	! the simple plume fit to the MPI Aerosol Climatology (Version 2). It sums over nplumes to provide a profile of aerosol
229	! optical properties on a host models vertical grid.
230	!
231	SUBROUTINE sp_aop_profile ( &
232	nlevels ,ncol ,lambda ,oro ,lon ,lat , &
233	year_fr ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , &
234	asy_prof )
235	!
236	! ----------
237	!
238	INTEGER, INTENT(IN) :: &
239	nlevels, & !< number of levels
240	ncol !< number of columns
241
242	REAL, INTENT(IN) :: &
243	lambda, & !< wavelength
244	year_fr, & !< Fractional Year (1903.0 is the 0Z on the first of January 1903, Gregorian)
245	oro(ncol), & !< orographic height (m)
246	lon(ncol), & !< longitude
247	lat(ncol), & !< latitude
248	z (ncol,nlevels), & !< height above sea-level (m)
249	dz(ncol,nlevels) !< level thickness (difference between half levels) (m)
250
251	REAL, INTENT(OUT) :: &
252	dNovrN(ncol) , & !< anthropogenic increase in cloud drop number concentration (factor)
253	aod_prof(ncol,nlevels) , & !< profile of aerosol optical depth
254	ssa_prof(ncol,nlevels) , & !< profile of single scattering albedo
255	asy_prof(ncol,nlevels) !< profile of asymmetry parameter
256
257	INTEGER :: iplume, icol, k
258
259	REAL :: &
260	eta(ncol,nlevels), & !< normalized height (by 15 km)
261	z_beta(ncol,nlevels), & !< profile for scaling column optical depth
262	prof(ncol,nlevels), & !< scaled profile (by beta function)
263	beta_sum(ncol), & !< vertical sum of beta function
264	ssa(ncol), & !< single scattering albedo
265	asy(ncol), & !< asymmetry parameter
266	cw_an(ncol), & !< column weight for simple plume (anthropogenic) AOD at 550 nm
267	cw_bg(ncol), & !< column weight for fine-mode natural background AOD at 550 nm
268	caod_sp(ncol), & !< column simple plume anthropogenic AOD at 550 nm
269	caod_bg(ncol), & !< column fine-mode natural background AOD at 550 nm
270	a_plume1, & !< gaussian longitude factor for feature 1
271	a_plume2, & !< gaussian longitude factor for feature 2
272	b_plume1, & !< gaussian latitude factor for feature 1
273	b_plume2, & !< gaussian latitude factor for feature 2
274	delta_lat, & !< latitude offset
275	delta_lon, & !< longitude offset
276	delta_lon_t, & !< threshold for maximum longitudinal plume extent used in transition from 360 to 0 degrees
277	lon1, & !< rotated longitude for feature 1
278	lat1, & !< rotated latitude for feature 2
279	lon2, & !< rotated longitude for feature 1
280	lat2, & !< rotated latitude for feature 2
281	f1, & !< contribution from feature 1
282	f2, & !< contribution from feature 2
283	f3, & !< contribution from feature 1 in natural background of Twomey effect
284	f4, & !< contribution from feature 2 in natural background of Twomey effect
285	aod_550, & !< aerosol optical depth at 550nm
286	aod_lmd, & !< aerosol optical depth at input wavelength
287	lfactor !< factor to compute wavelength dependence of optical properties
288	!
289	! ----------
290	!
291	! initialize input data (by calling setup at first instance)
292	!
293	IF (.NOT.sp_initialized) CALL sp_setup
294	!
295	! get time weights
296	!
297	CALL set_time_weight(year_fr)
298	!
299	! initialize variables, including output
300	!
301	DO k=1,nlevels
302	DO icol=1,ncol
303	aod_prof(icol,k) = 0.0
304	ssa_prof(icol,k) = 0.0
305	asy_prof(icol,k) = 0.0
306	z_beta(icol,k) = MERGE(1.0, 0.0, z(icol,k) >= oro(icol))
307	eta(icol,k) = MAX(0.0,MIN(1.0,z(icol,k)/15000.))
308	END DO
309	END DO
310	DO icol=1,ncol
311	dNovrN(icol) = 1.0
312	caod_sp(icol) = 0.0
313	caod_bg(icol) = 0.02
314	END DO
315	!
316	! sum contribution from plumes to construct composite profiles of aerosol optical properties
317	!
318	DO iplume=1,nplumes
319	!
320	! calculate vertical distribution function from parameters of beta distribution
321	!
322	DO icol=1,ncol
323	beta_sum(icol) = 0.
324	END DO
325	DO k=1,nlevels
326	DO icol=1,ncol
327	prof(icol,k) = (eta(icol,k)*(beta_a(iplume)-1.) (1.-eta(icol,k))*(beta_b(iplume)-1.)) dz(icol,k)
328	beta_sum(icol) = beta_sum(icol) + prof(icol,k)
329	END DO
330	END DO
331	DO k=1,nlevels
332	DO icol=1,ncol
333	prof(icol,k) = ( prof(icol,k) / beta_sum(icol) ) * z_beta(icol,k)
334	END DO
335	END DO
336	!
337	! calculate plume weights
338	!
339	DO icol=1,ncol
340	!
341	! get plume-center relative spatial parameters for specifying amplitude of plume at given lat and lon
342	!
343	delta_lat = lat(icol) - plume_lat(iplume)
344	delta_lon = lon(icol) - plume_lon(iplume)
345	delta_lon_t = MERGE (260., 180., iplume == 1)
346	delta_lon = MERGE ( delta_lon-SIGN(360.,delta_lon) , delta_lon , ABS(delta_lon) > delta_lon_t)
347
348	a_plume1 = 0.5 / (MERGE(sig_lon_E(1,iplume), sig_lon_W(1,iplume), delta_lon > 0)**2)
349	b_plume1 = 0.5 / (MERGE(sig_lat_E(1,iplume), sig_lat_W(1,iplume), delta_lon > 0)**2)
350	a_plume2 = 0.5 / (MERGE(sig_lon_E(2,iplume), sig_lon_W(2,iplume), delta_lon > 0)**2)
351	b_plume2 = 0.5 / (MERGE(sig_lat_E(2,iplume), sig_lat_W(2,iplume), delta_lon > 0)**2)
352	!
353	! adjust for a plume specific rotation which helps match plume state to climatology.
354	!
355	lon1 = COS(theta(1,iplume))(delta_lon) + SIN(theta(1,iplume))(delta_lat)
356	lat1 = - SIN(theta(1,iplume))(delta_lon) + COS(theta(1,iplume))(delta_lat)
357	lon2 = COS(theta(2,iplume))(delta_lon) + SIN(theta(2,iplume))(delta_lat)
358	lat2 = - SIN(theta(2,iplume))(delta_lon) + COS(theta(2,iplume))(delta_lat)
359	!
360	! calculate contribution to plume from its different features, to get a column weight for the anthropogenic
361	! (cw_an) and the fine-mode natural background aerosol (cw_bg)
362	!
363	f1 = time_weight(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)*2) + (b_plume1 ((lat1)**2))))
364	f2 = time_weight(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)*2) + (b_plume2 ((lat2)**2))))
365	f3 = time_weight_bg(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)*2) + (b_plume1 ((lat1)**2))))
366	f4 = time_weight_bg(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)*2) + (b_plume2 ((lat2)**2))))
367
368	cw_an(icol) = f1 * aod_spmx(iplume) + f2 * aod_spmx(iplume)
369	cw_bg(icol) = f3 * aod_fmbg(iplume) + f4 * aod_fmbg(iplume)
370	!
371	! calculate wavelength-dependent scattering properties
372	!
373	lfactor = MIN(1.0,700.0/lambda)
374	ssa(icol) = (ssa550(iplume) * lfactor*4) / ((ssa550(iplume) lfactor*4) + ((1-ssa550(iplume)) lfactor))
375	asy(icol) = asy550(iplume) * SQRT(lfactor)
376	END DO
377	!
378	! distribute plume optical properties across its vertical profile weighting by optical depth and scaling for
379	! wavelength using the angstrom parameter.
380	!
381	lfactor = EXP(-angstrom(iplume) * LOG(lambda/550.0))
382	DO k=1,nlevels
383	DO icol = 1,ncol
384	aod_550 = prof(icol,k) * cw_an(icol)
385	aod_lmd = aod_550 * lfactor
386	caod_sp(icol) = caod_sp(icol) + aod_550
387	caod_bg(icol) = caod_bg(icol) + prof(icol,k) * cw_bg(icol)
388	asy_prof(icol,k) = asy_prof(icol,k) + aod_lmd * ssa(icol) * asy(icol)
389	ssa_prof(icol,k) = ssa_prof(icol,k) + aod_lmd * ssa(icol)
390	aod_prof(icol,k) = aod_prof(icol,k) + aod_lmd
391	END DO
392	END DO
393	END DO
394	!
395	! complete optical depth weighting
396	!
397	DO k=1,nlevels
398	DO icol = 1,ncol
399	asy_prof(icol,k) = MERGE(asy_prof(icol,k)/ssa_prof(icol,k), 0.0, ssa_prof(icol,k) > TINY(1.))
400	ssa_prof(icol,k) = MERGE(ssa_prof(icol,k)/aod_prof(icol,k), 1.0, aod_prof(icol,k) > TINY(1.))
401	END DO
402	END DO
403	!
404	! calculate effective radius normalization (divisor) factor
405	!
406	DO icol=1,ncol
407	dNovrN(icol) = LOG((1000.0 * (caod_sp(icol) + caod_bg(icol))) + 1.0)/LOG((1000.0 * caod_bg(icol)) + 1.0)
408	END DO
409
410	RETURN
411	END SUBROUTINE sp_aop_profile
412
413	END MODULE MO_SIMPLE_PLUMES

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

Download in other formats:

Original Format