Context Navigation

readaerosolstrato_ecrad.F90 @ 4790

Last change on this file since 4790 was 4790, checked in by idelkadi, 4 months ago
Continued work on implementing the Ecrad code in LMDZ (Integration of aerosols):
File size: 24.6 KB

Line
1	!
2	! $Id: readaerosolstrato_ecrad.F90 tlurton $
3	!
4	SUBROUTINE readaerosolstrato_ecrad(config, debut, ok_volcan)
5
6	USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_dimid, &
7	nf95_inq_varid, nf95_open
8	USE netcdf, ONLY: nf90_get_var, nf90_noerr, nf90_nowrite
9
10	USE phys_cal_mod, ONLY : mth_cur
11	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, grid2dTo1d_glo, grid_type, unstructured
12	USE mod_phys_lmdz_mpi_data
13	USE mod_phys_lmdz_omp_data
14	USE mod_phys_lmdz_para
15	USE phys_state_var_mod
16	USE phys_local_var_mod
17	! USE physiq_mod, ONLY : ok_volcan
18	USE aero_mod
19	USE dimphy
20	! USE YOERAD
21	USE radiation_config, ONLY : config_type
22	USE YOMCST
23	#ifdef CPP_XIOS
24	USE xios
25	#endif
26
27	IMPLICIT NONE
28
29	INCLUDE "clesphys.h"
30
31	CHARACTER (len = 80) :: abort_message
32	CHARACTER (LEN=20) :: modname = 'readaerosolstrato_ecrad'
33
34	! Variable input
35	LOGICAL, INTENT(IN) :: debut
36	LOGICAL, INTENT(IN) :: ok_volcan !activate volcanic diags
37
38	! Variables locales
39	INTEGER n_lat ! number of latitudes in the input data
40	INTEGER n_lon ! number of longitudes
41	INTEGER n_lev ! number of levels in the input data
42	INTEGER n_month ! number of months in the input data
43	INTEGER n_wav ! number of wavelengths in the input data
44	REAL, ALLOCATABLE:: latitude(:)
45	REAL, ALLOCATABLE:: time(:)
46	REAL, ALLOCATABLE:: lev(:)
47	REAL, ALLOCATABLE:: wav(:)
48	INTEGER i,k,wave,band
49	INTEGER, SAVE :: mth_pre=1
50	!$OMP THREADPRIVATE(mth_pre)
51
52	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: tau_aer_strat
53	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: piz_aer_strat
54	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cg_aer_strat
55	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: taulw_aer_strat
56	! add ThL
57	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: pizlw_aer_strat
58	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cglw_aer_strat
59	! end add ThL
60	! mod ThL
61	!$OMP THREADPRIVATE(tau_aer_strat,piz_aer_strat,cg_aer_strat,taulw_aer_strat,pizlw_aer_strat,cglw_aer_strat)
62	! end mod ThL
63
64	! Champs reconstitues
65	REAL, ALLOCATABLE:: tauaerstrat(:, :, :, :)
66	REAL, ALLOCATABLE:: pizaerstrat(:, :, :, :)
67	REAL, ALLOCATABLE:: cgaerstrat(:, :, :, :)
68	REAL, ALLOCATABLE:: taulwaerstrat(:, :, :, :)
69	! add ThL
70	REAL, ALLOCATABLE:: pizlwaerstrat(:, :, :, :)
71	REAL, ALLOCATABLE:: cglwaerstrat(:, :, :, :)
72	! end add ThL
73
74	REAL, ALLOCATABLE:: tauaerstrat_mois(:, :, :, :)
75	REAL, ALLOCATABLE:: pizaerstrat_mois(:, :, :, :)
76	REAL, ALLOCATABLE:: cgaerstrat_mois(:, :, :, :)
77	REAL, ALLOCATABLE:: taulwaerstrat_mois(:, :, :, :)
78	! add ThL
79	REAL, ALLOCATABLE:: pizlwaerstrat_mois(:, :, :, :)
80	REAL, ALLOCATABLE:: cglwaerstrat_mois(:, :, :, :)
81	! end add ThL
82
83	REAL, ALLOCATABLE:: tauaerstrat_mois_glo(:, :, :)
84	REAL, ALLOCATABLE:: pizaerstrat_mois_glo(:, :, :)
85	REAL, ALLOCATABLE:: cgaerstrat_mois_glo(:, :, :)
86	REAL, ALLOCATABLE:: taulwaerstrat_mois_glo(:, :, :)
87	! add ThL
88	REAL, ALLOCATABLE:: pizlwaerstrat_mois_glo(:, :, :)
89	REAL, ALLOCATABLE:: cglwaerstrat_mois_glo(:, :, :)
90	! end add ThL
91	REAL, ALLOCATABLE:: tauaerstrat_mpi(:, :, :)
92	REAL, ALLOCATABLE:: pizaerstrat_mpi(:, :, :)
93	REAL, ALLOCATABLE:: cgaerstrat_mpi(:, :, :)
94	REAL, ALLOCATABLE:: taulwaerstrat_mpi(:, :, :)
95	! add ThL
96	REAL, ALLOCATABLE:: pizlwaerstrat_mpi(:, :, :)
97	REAL, ALLOCATABLE:: cglwaerstrat_mpi(:, :, :)
98	! end add ThL
99
100	! For NetCDF:
101	INTEGER ncid_in ! IDs for input files
102	INTEGER varid, ncerr
103
104	TYPE(config_type), INTENT(IN) :: config
105
106	!--------------------------------------------------------
107
108	IF (.not.ALLOCATED(tau_aer_strat)) ALLOCATE(tau_aer_strat(klon,klev,config%n_bands_sw))
109	IF (.not.ALLOCATED(piz_aer_strat)) ALLOCATE(piz_aer_strat(klon,klev,config%n_bands_sw))
110	IF (.not.ALLOCATED(cg_aer_strat)) ALLOCATE(cg_aer_strat(klon,klev,config%n_bands_sw))
111
112	IF (.not.ALLOCATED(taulw_aer_strat)) ALLOCATE(taulw_aer_strat(klon,klev,config%n_bands_lw))
113	! add ThL
114	IF (.not.ALLOCATED(pizlw_aer_strat)) ALLOCATE(pizlw_aer_strat(klon,klev,config%n_bands_lw))
115	IF (.not.ALLOCATED(cglw_aer_strat)) ALLOCATE(cglw_aer_strat(klon,klev,config%n_bands_lw))
116	! end add ThL
117
118	!--we only read monthly strat aerosol data
119	IF (debut.OR.mth_cur.NE.mth_pre) THEN
120
121	!--only root reads the data
122	IF (is_mpi_root.AND.is_omp_root) THEN
123
124	!--check mth_cur
125	IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
126	print *,'probleme avec le mois dans readaerosolstrat =', mth_cur
127	ENDIF
128
129	!--initialize n_lon as input data is 2D (lat-alt) only
130	n_lon = nbp_lon
131
132	!--Starts with SW optical properties
133
134	CALL nf95_open("tauswstrat.2D.nc", nf90_nowrite, ncid_in)
135
136	CALL nf95_inq_varid(ncid_in, "LEV", varid)
137	CALL nf95_gw_var(ncid_in, varid, lev)
138	n_lev = size(lev)
139	IF (n_lev.NE.klev) THEN
140	abort_message='Le nombre de niveaux n est pas egal a klev'
141	CALL abort_physic(modname,abort_message,1)
142	ENDIF
143
144	CALL nf95_inq_varid(ncid_in, "LAT", varid)
145	CALL nf95_gw_var(ncid_in, varid, latitude)
146	n_lat = size(latitude)
147
148	IF (grid_type/=unstructured) THEN
149	IF (n_lat.NE.nbp_lat) THEN
150	print *, 'latitude=', n_lat, nbp_lat
151	abort_message='Le nombre de lat n est pas egal a nbp_lat'
152	CALL abort_physic(modname,abort_message,1)
153	ENDIF
154	ENDIF
155
156	CALL nf95_inq_varid(ncid_in, "TIME", varid)
157	CALL nf95_gw_var(ncid_in, varid, time)
158	n_month = size(time)
159	IF (n_month.NE.12) THEN
160	abort_message='Le nombre de month n est pas egal a 12'
161	CALL abort_physic(modname,abort_message,1)
162	ENDIF
163
164	CALL nf95_inq_varid(ncid_in, "WAV", varid)
165	CALL nf95_gw_var(ncid_in, varid, wav)
166	n_wav = size(wav)
167	print *, 'WAV aerosol strato=', n_wav, wav
168	! IF (n_wav.NE.config%n_bands_sw) THEN
169	! abort_message='Le nombre de wav n est pas egal a config%n_bands_sw'
170	! CALL abort_physic(modname,abort_message,1)
171	! ENDIF
172
173	ALLOCATE(tauaerstrat(n_lat, n_lev, n_wav, n_month))
174	ALLOCATE(pizaerstrat(n_lat, n_lev, n_wav, n_month))
175	ALLOCATE(cgaerstrat(n_lat, n_lev, n_wav, n_month))
176
177	!--reading stratospheric aerosol tau per layer
178	CALL nf95_inq_varid(ncid_in, "TAU_SUN", varid)
179	ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)
180	print *,'code erreur readaerosolstrato=', ncerr, varid
181
182	!--reading stratospheric aerosol omega per layer
183	CALL nf95_inq_varid(ncid_in, "OME_SUN", varid)
184	ncerr = nf90_get_var(ncid_in, varid, pizaerstrat)
185	print *,'code erreur readaerosolstrato=', ncerr, varid
186
187	!--reading stratospheric aerosol g per layer
188	CALL nf95_inq_varid(ncid_in, "GGG_SUN", varid)
189	ncerr = nf90_get_var(ncid_in, varid, cgaerstrat)
190	print *,'code erreur readaerosolstrato sw=', ncerr, varid
191
192	CALL nf95_close(ncid_in)
193
194
195	IF (grid_type/=unstructured) THEN
196	ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
197	ALLOCATE(pizaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
198	ALLOCATE(cgaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
199
200	ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev, n_wav))
201	ALLOCATE(pizaerstrat_mois_glo(klon_glo, n_lev, n_wav))
202	ALLOCATE(cgaerstrat_mois_glo(klon_glo, n_lev, n_wav))
203	!--select the correct month
204	!--and copy into 1st longitude
205	tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur)
206	pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur)
207	cgaerstrat_mois(1,:,:,:) = cgaerstrat(:,:,:,mth_cur)
208
209	!--copy longitudes
210	DO i=2, n_lon
211	tauaerstrat_mois(i,:,:,:) = tauaerstrat_mois(1,:,:,:)
212	pizaerstrat_mois(i,:,:,:) = pizaerstrat_mois(1,:,:,:)
213	cgaerstrat_mois(i,:,:,:) = cgaerstrat_mois(1,:,:,:)
214	ENDDO
215
216	!---reduce to a klon_glo grid
217	DO band=1, config%n_bands_sw
218	CALL grid2dTo1d_glo(tauaerstrat_mois(:,:,:,band),tauaerstrat_mois_glo(:,:,band))
219	CALL grid2dTo1d_glo(pizaerstrat_mois(:,:,:,band),pizaerstrat_mois_glo(:,:,band))
220	CALL grid2dTo1d_glo(cgaerstrat_mois(:,:,:,band),cgaerstrat_mois_glo(:,:,band))
221	ENDDO
222	ENDIF
223	!--Now LW optical properties
224	!
225
226	CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in)
227
228	CALL nf95_inq_varid(ncid_in, "LEV", varid)
229	CALL nf95_gw_var(ncid_in, varid, lev)
230	n_lev = size(lev)
231	IF (n_lev.NE.klev) THEN
232	abort_message='Le nombre de niveaux n est pas egal a klev'
233	CALL abort_physic(modname,abort_message,1)
234	ENDIF
235
236	CALL nf95_inq_varid(ncid_in, "LAT", varid)
237	CALL nf95_gw_var(ncid_in, varid, latitude)
238	n_lat = size(latitude)
239
240	IF (grid_type/=unstructured) THEN
241	IF (n_lat.NE.nbp_lat) THEN
242	abort_message='Le nombre de lat n est pas egal a nbp_lat'
243	CALL abort_physic(modname,abort_message,1)
244	ENDIF
245	ENDIF
246
247	CALL nf95_inq_varid(ncid_in, "TIME", varid)
248	CALL nf95_gw_var(ncid_in, varid, time)
249	n_month = size(time)
250	IF (n_month.NE.12) THEN
251	abort_message='Le nombre de month n est pas egal a 12'
252	CALL abort_physic(modname,abort_message,1)
253	ENDIF
254
255	CALL nf95_inq_varid(ncid_in, "WAV", varid)
256	CALL nf95_gw_var(ncid_in, varid, wav)
257	n_wav = size(wav)
258	print *, 'WAV aerosol strato=', n_wav, wav
259	! IF (n_wav.NE.config%n_bands_lw) THEN
260	! abort_message='Le nombre de wav n est pas egal a config%n_bands_lw'
261	! CALL abort_physic(modname,abort_message,1)
262	! ENDIF
263
264	ALLOCATE(taulwaerstrat(n_lat, n_lev, n_wav, n_month))
265	! add ThL
266	ALLOCATE(pizlwaerstrat(n_lat, n_lev, n_wav, n_month))
267	ALLOCATE(cglwaerstrat(n_lat, n_lev, n_wav, n_month))
268	! end add ThL
269
270	!--reading stratospheric aerosol lw tau per layer
271	CALL nf95_inq_varid(ncid_in, "TAU_EAR", varid)
272	ncerr = nf90_get_var(ncid_in, varid, taulwaerstrat)
273	print *,'code erreur readaerosolstrato lw=', ncerr, varid
274
275	! add/mod ThL
276	!--with ECRAd, we also read omega and g:
277	!--reading stratospheric aerosol omega per layer
278	CALL nf95_inq_varid(ncid_in, "OME_EAR", varid)
279	ncerr = nf90_get_var(ncid_in, varid, pizlwaerstrat)
280	print *,'code erreur readaerosolstrato lw=', ncerr, varid
281
282	!--reading stratospheric aerosol g per layer
283	CALL nf95_inq_varid(ncid_in, "GGG_EAR", varid)
284	ncerr = nf90_get_var(ncid_in, varid, cglwaerstrat)
285	print *,'code erreur readaerosolstrato lw=', ncerr, varid
286
287	CALL nf95_close(ncid_in)
288
289	IF (grid_type/=unstructured) THEN
290
291	ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
292	ALLOCATE(pizlwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
293	ALLOCATE(cglwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
294
295	ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
296	ALLOCATE(pizlwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
297	ALLOCATE(cglwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
298	!--select the correct month
299	!--and copy into 1st longitude
300	taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur)
301	pizlwaerstrat_mois(1,:,:,:) = pizlwaerstrat(:,:,:,mth_cur)
302	cglwaerstrat_mois(1,:,:,:) = cglwaerstrat(:,:,:,mth_cur)
303	!--copy longitudes
304	DO i=2, n_lon
305	taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:)
306	pizlwaerstrat_mois(i,:,:,:) = pizlwaerstrat_mois(1,:,:,:)
307	cglwaerstrat_mois(i,:,:,:) = cglwaerstrat_mois(1,:,:,:)
308	ENDDO
309
310	!---reduce to a klon_glo grid
311	DO band=1, config%n_bands_lw
312	CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band))
313	CALL grid2dTo1d_glo(pizlwaerstrat_mois(:,:,:,band),pizlwaerstrat_mois_glo(:,:,band))
314	CALL grid2dTo1d_glo(cglwaerstrat_mois(:,:,:,band),cglwaerstrat_mois_glo(:,:,band))
315	ENDDO
316	ENDIF
317
318	ELSE !--proc other than mpi_root and omp_root
319	!--dummy allocation needed for debug mode
320
321	ALLOCATE(tauaerstrat_mois_glo(1,1,1))
322	ALLOCATE(pizaerstrat_mois_glo(1,1,1))
323	ALLOCATE(cgaerstrat_mois_glo(1,1,1))
324	ALLOCATE(taulwaerstrat_mois_glo(1,1,1))
325	ALLOCATE(pizlwaerstrat_mois_glo(1,1,1))
326	ALLOCATE(cglwaerstrat_mois_glo(1,1,1))
327
328	ALLOCATE(tauaerstrat(0,0,0,12))
329	ALLOCATE(pizaerstrat(0,0,0,12))
330	ALLOCATE(cgaerstrat(0,0,0,12))
331	ALLOCATE(taulwaerstrat(0,0,0,12))
332	ALLOCATE(pizlwaerstrat(0,0,0,12))
333	ALLOCATE(cglwaerstrat(0,0,0,12))
334	! end add/mod ThL
335
336	ENDIF !--is_mpi_root and is_omp_root
337
338	!$OMP BARRIER
339
340	!--keep memory of previous month
341	mth_pre=mth_cur
342
343	IF (grid_type==unstructured) THEN
344
345	#ifdef CPP_XIOS
346
347	IF (is_omp_master) THEN
348	ALLOCATE(tauaerstrat_mpi(klon_mpi, klev, config%n_bands_sw))
349	ALLOCATE(pizaerstrat_mpi(klon_mpi, klev, config%n_bands_sw))
350	ALLOCATE(cgaerstrat_mpi(klon_mpi, klev, config%n_bands_sw))
351	ALLOCATE(taulwaerstrat_mpi(klon_mpi, klev, config%n_bands_lw))
352	! add ThL
353	ALLOCATE(pizlwaerstrat_mpi(klon_mpi, klev, config%n_bands_lw))
354	ALLOCATE(cglwaerstrat_mpi(klon_mpi, klev, config%n_bands_lw))
355	! end add ThL
356
357	CALL xios_send_field("tauaerstrat_in",SPREAD(tauaerstrat(:,:,:,mth_cur),1,8))
358	CALL xios_recv_field("tauaerstrat_out",tauaerstrat_mpi)
359	CALL xios_send_field("pizaerstrat_in",SPREAD(pizaerstrat(:,:,:,mth_cur),1,8))
360	CALL xios_recv_field("pizaerstrat_out",pizaerstrat_mpi)
361	CALL xios_send_field("cgaerstrat_in",SPREAD(cgaerstrat(:,:,:,mth_cur),1,8))
362	CALL xios_recv_field("cgaerstrat_out",cgaerstrat_mpi)
363	CALL xios_send_field("taulwaerstrat_in",SPREAD(taulwaerstrat(:,:,:,mth_cur),1,8))
364	CALL xios_recv_field("taulwaerstrat_out",taulwaerstrat_mpi)
365	! add ThL
366	CALL xios_send_field("pizlwaerstrat_in",SPREAD(pizlwaerstrat(:,:,:,mth_cur),1,8))
367	CALL xios_recv_field("pizlwaerstrat_out",pizlwaerstrat_mpi)
368	CALL xios_send_field("cglwaerstrat_in",SPREAD(cglwaerstrat(:,:,:,mth_cur),1,8))
369	CALL xios_recv_field("cglwaerstrat_out",cglwaerstrat_mpi)
370	! end add ThL
371	ELSE
372	ALLOCATE(tauaerstrat_mpi(0, 0, 0))
373	ALLOCATE(pizaerstrat_mpi(0, 0, 0))
374	ALLOCATE(cgaerstrat_mpi(0, 0, 0))
375	ALLOCATE(taulwaerstrat_mpi(0, 0, 0))
376	! add ThL
377	ALLOCATE(pizlwaerstrat_mpi(0, 0, 0))
378	ALLOCATE(cglwaerstrat_mpi(0, 0, 0))
379	! end add ThL
380	ENDIF
381
382	CALL scatter_omp(tauaerstrat_mpi,tau_aer_strat)
383	CALL scatter_omp(pizaerstrat_mpi,piz_aer_strat)
384	CALL scatter_omp(cgaerstrat_mpi,cg_aer_strat)
385	CALL scatter_omp(taulwaerstrat_mpi,taulw_aer_strat)
386	! add ThL
387	CALL scatter_omp(pizlwaerstrat_mpi,pizlw_aer_strat)
388	CALL scatter_omp(cglwaerstrat_mpi,cglw_aer_strat)
389	! end add ThL
390	#endif
391	ELSE
392
393	!--scatter on all proc
394	CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
395	CALL scatter(pizaerstrat_mois_glo,piz_aer_strat)
396	CALL scatter(cgaerstrat_mois_glo,cg_aer_strat)
397	CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat)
398	! add ThL
399	CALL scatter(pizlwaerstrat_mois_glo,pizlw_aer_strat)
400	CALL scatter(cglwaerstrat_mois_glo,cglw_aer_strat)
401	! end add ThL
402	IF (is_mpi_root.AND.is_omp_root) DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois, taulwaerstrat_mois, pizlwaerstrat_mois, cglwaerstrat_mois) ! mod ThL
403
404	ENDIF
405
406	IF (is_mpi_root.AND.is_omp_root) THEN
407	DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat, taulwaerstrat, pizlwaerstrat, cglwaerstrat) ! mod ThL
408	ENDIF
409
410
411	!$OMP BARRIER
412
413	ENDIF !--debut ou nouveau mois
414
415	!--total vertical aod at the SW wavelengths
416	!--for now use ECRAD band #10 AOD into all wavelengths
417	!--it is only a reasonable approximation for 550 nm
418	!--WAV(10)=0.53 µm
419	!--<!!!> This is considering that wavebands are in their original order,
420	!--and have NOT been re-ordered monotonously, i.e. with
421	!--WAV(1) = 3.46 µm
422	!--WAV(2) = 2.79 µm
423	!--(...)
424	!--WAV(last-1) = 0.23 µm
425	!--WAV(last)= 8.02 µm
426	!--If bands were to be re-ordered, i.e. last band be put in 1st position,
427	!--then we should consider band=11 to fetch WAV=0.53 µm.
428	band=10
429	DO i=1, klon
430	DO k=1, klev
431	IF (stratomask(i,k).GT.0.999999) THEN
432	DO wave=1, config%n_bands_sw
433	tausum_aero(i,wave,id_STRAT_phy)=tausum_aero(i,wave,id_STRAT_phy)+tau_aer_strat(i,k,band)
434	ENDDO
435	ENDIF
436	ENDDO
437	ENDDO
438
439	IF (.NOT. ok_volcan) THEN
440	!
441	!--this is the default case
442	!--stratospheric aerosols are added to both index 2 and 1 for double radiation calls
443	!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
444	DO band=1, config%n_bands_sw
445	WHERE (stratomask.GT.0.999999)
446	!--strat aerosols are added to index 2: natural and anthropogenic aerosols for bands 1 to config%n_bands_sw
447	cg_aero_sw_rrtm(:,:,2,band) = ( cg_aero_sw_rrtm(:,:,2,band)piz_aero_sw_rrtm(:,:,2,band)tau_aero_sw_rrtm(:,:,2,band) + &
448	cg_aer_strat(:,:,band)piz_aer_strat(:,:,band)tau_aer_strat(:,:,band) ) / &
449	MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
450	piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
451	piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
452	piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / &
453	MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 )
454	tau_aero_sw_rrtm(:,:,2,band) = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band)
455	!--strat aerosols are added to index 1: natural aerosols only for bands 1 to config%n_bands_sw
456	cg_aero_sw_rrtm(:,:,1,band) = ( cg_aero_sw_rrtm(:,:,1,band)piz_aero_sw_rrtm(:,:,1,band)tau_aero_sw_rrtm(:,:,1,band) + &
457	cg_aer_strat(:,:,band)piz_aer_strat(:,:,band)tau_aer_strat(:,:,band) ) / &
458	MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
459	piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
460	piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
461	piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / &
462	MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 )
463	tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band)
464	ENDWHERE
465	ENDDO
466	!
467	ELSE
468	!
469	!--this is the VOLMIP case
470	!--stratospheric aerosols are only added to index 2 in this case
471	!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
472	DO band=1, config%n_bands_sw
473	WHERE (stratomask.GT.0.999999)
474	!--strat aerosols are added to index 2 : natural and anthropogenic aerosols for bands 1 to config%n_bands_sw
475	cg_aero_sw_rrtm(:,:,2,band) = ( cg_aero_sw_rrtm(:,:,2,band)piz_aero_sw_rrtm(:,:,2,band)tau_aero_sw_rrtm(:,:,2,band) + &
476	cg_aer_strat(:,:,band)piz_aer_strat(:,:,band)tau_aer_strat(:,:,band) ) / &
477	MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
478	piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
479	piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
480	piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / &
481	MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 )
482	tau_aero_sw_rrtm(:,:,2,band) = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band)
483	ENDWHERE
484	ENDDO
485	ENDIF
486
487	!--total vertical aod at 10 um
488	!--this is approximated from band 7 of ECRAD (9.73 µm)
489	band=7
490	DO i=1, klon
491	DO k=1, klev
492	IF (stratomask(i,k).GT.0.999999) THEN
493	DO wave=1, config%n_bands_lw
494	tausum_aero(i,config%n_bands_sw+wave,id_STRAT_phy)=tausum_aero(i,config%n_bands_sw+wave,id_STRAT_phy)+taulw_aer_strat(i,k,band)
495	ENDDO
496	ENDIF
497	ENDDO
498	ENDDO
499
500	IF (.NOT. ok_volcan) THEN
501	!--this is the default case
502	!--stratospheric aerosols are added to both index 2 and 1
503	DO band=1, config%n_bands_lw
504	WHERE (stratomask.GT.0.999999)
505	! add ThL
506	cg_aero_lw_rrtm(:,:,2,band) = ( cg_aero_lw_rrtm(:,:,2,band)piz_aero_lw_rrtm(:,:,2,band)tau_aero_lw_rrtm(:,:,2,band) + &
507	cglw_aer_strat(:,:,band)pizlw_aer_strat(:,:,band)taulw_aer_strat(:,:,band) ) / &
508	MAX( piz_aero_lw_rrtm(:,:,2,band)*tau_aero_lw_rrtm(:,:,2,band) + &
509	pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band), 1.e-15 )
510	piz_aero_lw_rrtm(:,:,2,band) = ( piz_aero_lw_rrtm(:,:,2,band)*tau_aero_lw_rrtm(:,:,2,band) + &
511	pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band) ) / &
512	MAX( tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band), 1.e-15 )
513	cg_aero_lw_rrtm(:,:,1,band) = ( cg_aero_lw_rrtm(:,:,1,band)piz_aero_lw_rrtm(:,:,1,band)tau_aero_lw_rrtm(:,:,1,band) + &
514	cglw_aer_strat(:,:,band)pizlw_aer_strat(:,:,band)taulw_aer_strat(:,:,band) ) / &
515	MAX( piz_aero_lw_rrtm(:,:,1,band)*tau_aero_lw_rrtm(:,:,1,band) + &
516	pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band), 1.e-15 )
517	piz_aero_lw_rrtm(:,:,1,band) = ( piz_aero_lw_rrtm(:,:,1,band)*tau_aero_lw_rrtm(:,:,1,band) + &
518	pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band) ) / &
519	MAX( tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band), 1.e-15 )
520	! end add ThL
521	tau_aero_lw_rrtm(:,:,2,band) = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)
522	tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band)
523	ENDWHERE
524	ENDDO
525	!
526	ELSE
527	!
528	! mod ThL
529	!--this is the VOLMIP case
530	DO band=1, config%n_bands_lw
531	WHERE (stratomask.GT.0.999999)
532	!--stratospheric aerosols are not added to index 1
533	!--and we copy index 2 in index 1 because we want the same dust aerosol LW properties as above
534	cg_aero_lw_rrtm(:,:,1,band) = cg_aero_lw_rrtm(:,:,2,band)
535	piz_aero_lw_rrtm(:,:,1,band) = piz_aero_lw_rrtm(:,:,2,band)
536	tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,2,band)
537
538	!--stratospheric aerosols are only added to index 2
539	cg_aero_lw_rrtm(:,:,2,band) = (cg_aero_lw_rrtm(:,:,2,band)piz_aero_lw_rrtm(:,:,2,band)tau_aero_lw_rrtm(:,:,2,band) + &
540	cglw_aer_strat(:,:,band)pizlw_aer_strat(:,:,band)taulw_aer_strat(:,:,band)) / &
541	MAX( piz_aero_lw_rrtm(:,:,2,band)*tau_aero_lw_rrtm(:,:,2,band) + &
542	pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band), 1.e-15 )
543	piz_aero_lw_rrtm(:,:,2,band) = (piz_aero_lw_rrtm(:,:,2,band)*tau_aero_lw_rrtm(:,:,2,band) + &
544	pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band) ) / &
545	MAX( tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band), 1.e-15 )
546	tau_aero_lw_rrtm(:,:,2,band) = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)
547	ENDWHERE
548	ENDDO
549	ENDIF
550	! end mod ThL
551
552	!--default SSA value if there is no aerosol
553	!--to avoid 0 values that seems to cause some problem to RRTM
554	WHERE (tau_aero_sw_rrtm.LT.1.e-14)
555	piz_aero_sw_rrtm = 1.0
556	ENDWHERE
557
558	!--in principle this should not be necessary
559	!--as these variables have min values already but just in case
560	!--put 1e-15 min value to both SW and LW AOD
561	tau_aero_sw_rrtm = MAX(tau_aero_sw_rrtm,1.e-15)
562	tau_aero_lw_rrtm = MAX(tau_aero_lw_rrtm,1.e-15)
563
564	END SUBROUTINE readaerosolstrato_ecrad

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

Download in other formats:

Original Format