Context Navigation

readaerosolstrato2_rrtm.F90 @ 3863

Last change on this file since 3863 was 3413, checked in by Laurent Fairhead, 6 years ago
Inclusion of Yann's latest (summer/fall 2018) modifications for convergence of DYNAMICO/LMDZ physics YM/LF
File size: 15.7 KB

Line
1	!
2	! $Id: readaerosolstrato2_rrtm.F90 2526 2016-05-26 22:13:40Z oboucher $
3	!
4
5	SUBROUTINE readaerosolstrato2_rrtm(debut)
6
7	USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_dimid, &
8	nf95_inq_varid, nf95_open
9	USE netcdf, ONLY: nf90_get_var, nf90_noerr, nf90_nowrite
10
11	USE phys_cal_mod, ONLY : mth_cur
12	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, grid2dTo1d_glo, grid_type, unstructured
13	USE mod_phys_lmdz_mpi_data
14	USE mod_phys_lmdz_omp_data
15	USE mod_phys_lmdz_para
16	USE phys_state_var_mod
17	USE phys_local_var_mod
18	USE aero_mod
19	USE dimphy
20	USE YOERAD, ONLY : NLW
21	USE YOMCST
22	USE xios
23	IMPLICIT NONE
24
25	INCLUDE "clesphys.h"
26
27	CHARACTER (len = 80) :: abort_message
28	CHARACTER (LEN=20) :: modname = 'readaerosolstrato2'
29
30	! Variable input
31	LOGICAL, INTENT(IN) :: debut
32
33	! Variables locales
34	INTEGER n_lat ! number of latitudes in the input data
35	INTEGER n_lon ! number of longitudes
36	INTEGER n_lev ! number of levels in the input data
37	INTEGER n_month ! number of months in the input data
38	INTEGER n_wav ! number of wavelengths in the input data
39	REAL, POINTER:: latitude(:)
40	REAL, POINTER:: time(:)
41	REAL, POINTER:: lev(:)
42	REAL, POINTER:: wav(:)
43	INTEGER i,k,wave,band
44	INTEGER, SAVE :: mth_pre=1
45	!$OMP THREADPRIVATE(mth_pre)
46
47	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: tau_aer_strat
48	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: piz_aer_strat
49	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cg_aer_strat
50	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: taulw_aer_strat
51	!$OMP THREADPRIVATE(tau_aer_strat,piz_aer_strat,cg_aer_strat,taulw_aer_strat)
52
53	! Champs reconstitues
54	REAL, ALLOCATABLE:: tauaerstrat(:, :, :, :)
55	REAL, ALLOCATABLE:: pizaerstrat(:, :, :, :)
56	REAL, ALLOCATABLE:: cgaerstrat(:, :, :, :)
57	REAL, ALLOCATABLE:: taulwaerstrat(:, :, :, :)
58
59	REAL, ALLOCATABLE:: tauaerstrat_mois(:, :, :, :)
60	REAL, ALLOCATABLE:: pizaerstrat_mois(:, :, :, :)
61	REAL, ALLOCATABLE:: cgaerstrat_mois(:, :, :, :)
62	REAL, ALLOCATABLE:: taulwaerstrat_mois(:, :, :, :)
63
64	REAL, ALLOCATABLE:: tauaerstrat_mois_glo(:, :, :)
65	REAL, ALLOCATABLE:: pizaerstrat_mois_glo(:, :, :)
66	REAL, ALLOCATABLE:: cgaerstrat_mois_glo(:, :, :)
67	REAL, ALLOCATABLE:: taulwaerstrat_mois_glo(:, :, :)
68	REAL, ALLOCATABLE:: tauaerstrat_mpi(:, :, :)
69	REAL, ALLOCATABLE:: pizaerstrat_mpi(:, :, :)
70	REAL, ALLOCATABLE:: cgaerstrat_mpi(:, :, :)
71	REAL, ALLOCATABLE:: taulwaerstrat_mpi(:, :, :)
72
73	! For NetCDF:
74	INTEGER ncid_in ! IDs for input files
75	INTEGER varid, ncerr
76
77	!--------------------------------------------------------
78
79	IF (.not.ALLOCATED(tau_aer_strat)) ALLOCATE(tau_aer_strat(klon,klev,NSW))
80	IF (.not.ALLOCATED(piz_aer_strat)) ALLOCATE(piz_aer_strat(klon,klev,NSW))
81	IF (.not.ALLOCATED(cg_aer_strat)) ALLOCATE(cg_aer_strat(klon,klev,NSW))
82
83	IF (.not.ALLOCATED(taulw_aer_strat)) ALLOCATE(taulw_aer_strat(klon,klev,NLW))
84
85	!--we only read monthly strat aerosol data
86	IF (debut.OR.mth_cur.NE.mth_pre) THEN
87
88	!--only root reads the data
89	IF (is_mpi_root.AND.is_omp_root) THEN
90
91	!--check mth_cur
92	IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
93	print *,'probleme avec le mois dans readaerosolstrat =', mth_cur
94	ENDIF
95
96	!--initialize n_lon as input data is 2D (lat-alt) only
97	n_lon = nbp_lon
98
99	!--Starts with SW optical properties
100
101	CALL nf95_open("tauswstrat.2D.nc", nf90_nowrite, ncid_in)
102
103	CALL nf95_inq_varid(ncid_in, "LEV", varid)
104	CALL nf95_gw_var(ncid_in, varid, lev)
105	n_lev = size(lev)
106	IF (n_lev.NE.klev) THEN
107	abort_message='Le nombre de niveaux n est pas egal a klev'
108	CALL abort_physic(modname,abort_message,1)
109	ENDIF
110
111	CALL nf95_inq_varid(ncid_in, "LAT", varid)
112	CALL nf95_gw_var(ncid_in, varid, latitude)
113	n_lat = size(latitude)
114
115	IF (grid_type/=unstructured) THEN
116	IF (n_lat.NE.nbp_lat) THEN
117	print *, 'latitude=', n_lat, nbp_lat
118	abort_message='Le nombre de lat n est pas egal a nbp_lat'
119	CALL abort_physic(modname,abort_message,1)
120	ENDIF
121	ENDIF
122
123	CALL nf95_inq_varid(ncid_in, "TIME", varid)
124	CALL nf95_gw_var(ncid_in, varid, time)
125	n_month = size(time)
126	IF (n_month.NE.12) THEN
127	abort_message='Le nombre de month n est pas egal a 12'
128	CALL abort_physic(modname,abort_message,1)
129	ENDIF
130
131	CALL nf95_inq_varid(ncid_in, "WAV", varid)
132	CALL nf95_gw_var(ncid_in, varid, wav)
133	n_wav = size(wav)
134	print *, 'WAV aerosol strato=', n_wav, wav
135	IF (n_wav.NE.NSW) THEN
136	abort_message='Le nombre de wav n est pas egal a NSW'
137	CALL abort_physic(modname,abort_message,1)
138	ENDIF
139
140	ALLOCATE(tauaerstrat(n_lat, n_lev, n_wav, n_month))
141	ALLOCATE(pizaerstrat(n_lat, n_lev, n_wav, n_month))
142	ALLOCATE(cgaerstrat(n_lat, n_lev, n_wav, n_month))
143
144	!--reading stratospheric aerosol tau per layer
145	CALL nf95_inq_varid(ncid_in, "TAU_SUN", varid)
146	ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)
147	print *,'code erreur readaerosolstrato=', ncerr, varid
148
149	!--reading stratospheric aerosol omega per layer
150	CALL nf95_inq_varid(ncid_in, "OME_SUN", varid)
151	ncerr = nf90_get_var(ncid_in, varid, pizaerstrat)
152	print *,'code erreur readaerosolstrato=', ncerr, varid
153
154	!--reading stratospheric aerosol g per layer
155	CALL nf95_inq_varid(ncid_in, "GGG_SUN", varid)
156	ncerr = nf90_get_var(ncid_in, varid, cgaerstrat)
157	print *,'code erreur readaerosolstrato sw=', ncerr, varid
158
159	CALL nf95_close(ncid_in)
160
161
162	IF (grid_type/=unstructured) THEN
163	ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
164	ALLOCATE(pizaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
165	ALLOCATE(cgaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
166
167	ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev, n_wav))
168	ALLOCATE(pizaerstrat_mois_glo(klon_glo, n_lev, n_wav))
169	ALLOCATE(cgaerstrat_mois_glo(klon_glo, n_lev, n_wav))
170	!--select the correct month
171	!--and copy into 1st longitude
172	tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur)
173	pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur)
174	cgaerstrat_mois(1,:,:,:) = cgaerstrat(:,:,:,mth_cur)
175
176	!--copy longitudes
177	DO i=2, n_lon
178	tauaerstrat_mois(i,:,:,:) = tauaerstrat_mois(1,:,:,:)
179	pizaerstrat_mois(i,:,:,:) = pizaerstrat_mois(1,:,:,:)
180	cgaerstrat_mois(i,:,:,:) = cgaerstrat_mois(1,:,:,:)
181	ENDDO
182
183	!---reduce to a klon_glo grid
184	DO band=1, NSW
185	CALL grid2dTo1d_glo(tauaerstrat_mois(:,:,:,band),tauaerstrat_mois_glo(:,:,band))
186	CALL grid2dTo1d_glo(pizaerstrat_mois(:,:,:,band),pizaerstrat_mois_glo(:,:,band))
187	CALL grid2dTo1d_glo(cgaerstrat_mois(:,:,:,band),cgaerstrat_mois_glo(:,:,band))
188	ENDDO
189	ENDIF
190	!--Now LW optical properties
191	!
192
193	CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in)
194
195	CALL nf95_inq_varid(ncid_in, "LEV", varid)
196	CALL nf95_gw_var(ncid_in, varid, lev)
197	n_lev = size(lev)
198	IF (n_lev.NE.klev) THEN
199	abort_message='Le nombre de niveaux n est pas egal a klev'
200	CALL abort_physic(modname,abort_message,1)
201	ENDIF
202
203	CALL nf95_inq_varid(ncid_in, "LAT", varid)
204	CALL nf95_gw_var(ncid_in, varid, latitude)
205	n_lat = size(latitude)
206
207	IF (grid_type/=unstructured) THEN
208	IF (n_lat.NE.nbp_lat) THEN
209	abort_message='Le nombre de lat n est pas egal a nbp_lat'
210	CALL abort_physic(modname,abort_message,1)
211	ENDIF
212	ENDIF
213
214	CALL nf95_inq_varid(ncid_in, "TIME", varid)
215	CALL nf95_gw_var(ncid_in, varid, time)
216	n_month = size(time)
217	IF (n_month.NE.12) THEN
218	abort_message='Le nombre de month n est pas egal a 12'
219	CALL abort_physic(modname,abort_message,1)
220	ENDIF
221
222	CALL nf95_inq_varid(ncid_in, "WAV", varid)
223	CALL nf95_gw_var(ncid_in, varid, wav)
224	n_wav = size(wav)
225	print *, 'WAV aerosol strato=', n_wav, wav
226	IF (n_wav.NE.NLW) THEN
227	abort_message='Le nombre de wav n est pas egal a NLW'
228	CALL abort_physic(modname,abort_message,1)
229	ENDIF
230
231	ALLOCATE(taulwaerstrat(n_lat, n_lev, n_wav, n_month))
232
233	!--reading stratospheric aerosol lw tau per layer
234	CALL nf95_inq_varid(ncid_in, "TAU_EAR", varid)
235	ncerr = nf90_get_var(ncid_in, varid, taulwaerstrat)
236	print *,'code erreur readaerosolstrato lw=', ncerr, varid
237
238	CALL nf95_close(ncid_in)
239
240	IF (grid_type/=unstructured) THEN
241
242	ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
243	ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
244
245	!--select the correct month
246	!--and copy into 1st longitude
247	taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur)
248	!--copy longitudes
249	DO i=2, n_lon
250	taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:)
251	ENDDO
252
253	!---reduce to a klon_glo grid
254	DO band=1, NLW
255	CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band))
256	ENDDO
257	ENDIF
258
259	ELSE !--proc other than mpi_root and omp_root
260	!--dummy allocation needed for debug mode
261
262	ALLOCATE(tauaerstrat_mois_glo(1,1,1))
263	ALLOCATE(pizaerstrat_mois_glo(1,1,1))
264	ALLOCATE(cgaerstrat_mois_glo(1,1,1))
265	ALLOCATE(taulwaerstrat_mois_glo(1,1,1))
266
267	ALLOCATE(tauaerstrat(0,0,0,12))
268	ALLOCATE(pizaerstrat(0,0,0,12))
269	ALLOCATE(cgaerstrat(0,0,0,12))
270	ALLOCATE(taulwaerstrat(0,0,0,12))
271
272
273	ENDIF !--is_mpi_root and is_omp_root
274
275	!$OMP BARRIER
276
277	!--keep memory of previous month
278	mth_pre=mth_cur
279
280	IF (grid_type==unstructured) THEN
281
282	IF (is_omp_master) THEN
283	ALLOCATE(tauaerstrat_mpi(klon_mpi, klev, NSW))
284	ALLOCATE(pizaerstrat_mpi(klon_mpi, klev, NSW))
285	ALLOCATE(cgaerstrat_mpi(klon_mpi, klev, NSW))
286	ALLOCATE(taulwaerstrat_mpi(klon_mpi, klev, NLW))
287
288	CALL xios_send_field("tauaerstrat_in",SPREAD(tauaerstrat(:,:,:,mth_cur),1,8))
289	CALL xios_recv_field("tauaerstrat_out",tauaerstrat_mpi)
290	CALL xios_send_field("pizaerstrat_in",SPREAD(pizaerstrat(:,:,:,mth_cur),1,8))
291	CALL xios_recv_field("pizaerstrat_out",pizaerstrat_mpi)
292	CALL xios_send_field("cgaerstrat_in",SPREAD(cgaerstrat(:,:,:,mth_cur),1,8))
293	CALL xios_recv_field("cgaerstrat_out",cgaerstrat_mpi)
294	CALL xios_send_field("taulwaerstrat_in",SPREAD(taulwaerstrat(:,:,:,mth_cur),1,8))
295	CALL xios_recv_field("taulwaerstrat_out",taulwaerstrat_mpi)
296	ELSE
297	ALLOCATE(tauaerstrat_mpi(0, 0, 0))
298	ALLOCATE(pizaerstrat_mpi(0, 0, 0))
299	ALLOCATE(cgaerstrat_mpi(0, 0, 0))
300	ALLOCATE(taulwaerstrat_mpi(0, 0, 0))
301	ENDIF
302
303	CALL scatter_omp(tauaerstrat_mpi,tau_aer_strat)
304	CALL scatter_omp(pizaerstrat_mpi,piz_aer_strat)
305	CALL scatter_omp(cgaerstrat_mpi,cg_aer_strat)
306	CALL scatter_omp(taulwaerstrat_mpi,taulw_aer_strat)
307	ELSE
308
309	!--scatter on all proc
310	CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
311	CALL scatter(pizaerstrat_mois_glo,piz_aer_strat)
312	CALL scatter(cgaerstrat_mois_glo,cg_aer_strat)
313	CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat)
314
315	DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois)
316	DEALLOCATE(taulwaerstrat_mois)
317	ENDIF
318
319	IF (is_mpi_root.AND.is_omp_root) THEN
320	DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat,taulwaerstrat)
321	ENDIF
322
323
324	!$OMP BARRIER
325
326	ENDIF !--debut ou nouveau mois
327
328	!--total vertical aod at the 5 SW wavelengths
329	!--for now use band 3 AOD into all 5 wavelengths
330	!--it is only a reasonable approximation for 550 nm (wave=2)
331	band=3
332	DO i=1, klon
333	DO k=1, klev
334	IF (stratomask(i,k).GT.0.999999) THEN
335	DO wave=1, nwave_sw
336	tausum_aero(i,wave,id_STRAT_phy)=tausum_aero(i,wave,id_STRAT_phy)+tau_aer_strat(i,k,band)
337	ENDDO
338	ENDIF
339	ENDDO
340	ENDDO
341
342	!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
343	DO band=1, NSW
344	WHERE (stratomask.GT.0.999999)
345	!--anthropogenic aerosols bands 1 to NSW
346	cg_aero_sw_rrtm(:,:,2,band) = ( cg_aero_sw_rrtm(:,:,2,band)piz_aero_sw_rrtm(:,:,2,band)tau_aero_sw_rrtm(:,:,2,band) + &
347	cg_aer_strat(:,:,band)piz_aer_strat(:,:,band)tau_aer_strat(:,:,band) ) / &
348	MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
349	piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
350	piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
351	piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / &
352	MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 )
353	tau_aero_sw_rrtm(:,:,2,band) = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band)
354	!--natural aerosols bands 1 to NSW
355	cg_aero_sw_rrtm(:,:,1,band) = ( cg_aero_sw_rrtm(:,:,1,band)piz_aero_sw_rrtm(:,:,1,band)tau_aero_sw_rrtm(:,:,1,band) + &
356	cg_aer_strat(:,:,band)piz_aer_strat(:,:,band)tau_aer_strat(:,:,band) ) / &
357	MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
358	piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
359	piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
360	piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / &
361	MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 )
362	tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band)
363	!--no stratospheric aerosol in index 1 for these tests
364	! cg_aero_sw_rrtm(:,:,1,band) = cg_aero_sw_rrtm(:,:,1,band)
365	! piz_aero_sw_rrtm(:,:,1,band) = piz_aero_sw_rrtm(:,:,1,band)
366	! tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band)
367	ENDWHERE
368	ENDDO
369
370	!--total vertical aod at 10 um
371	!--this is approximated from band 7 of RRTM
372	band=7
373	DO i=1, klon
374	DO k=1, klev
375	IF (stratomask(i,k).GT.0.999999) THEN
376	DO wave=1, nwave_lw
377	tausum_aero(i,nwave_sw+wave,id_STRAT_phy)=tausum_aero(i,nwave_sw+wave,id_STRAT_phy)+taulw_aer_strat(i,k,band)
378	ENDDO
379	ENDIF
380	ENDDO
381	ENDDO
382
383	DO band=1, NLW
384	WHERE (stratomask.GT.0.999999)
385	tau_aero_lw_rrtm(:,:,2,band) = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)
386	tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band)
387	!--no stratospheric aerosols in index 1 for these tests
388	! tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,1,band)
389	ENDWHERE
390	ENDDO
391
392	!--default SSA value if there is no aerosol
393	!--to avoid 0 values that seems to cause some problem to RRTM
394	WHERE (tau_aero_sw_rrtm.LT.1.e-14)
395	piz_aero_sw_rrtm = 1.0
396	ENDWHERE
397
398	!--in principle this should not be necessary
399	!--as these variables have min values already but just in case
400	!--put 1e-15 min value to both SW and LW AOD
401	tau_aero_sw_rrtm = MAX(tau_aero_sw_rrtm,1.e-15)
402	tau_aero_lw_rrtm = MAX(tau_aero_lw_rrtm,1.e-15)
403
404	END SUBROUTINE readaerosolstrato2_rrtm

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

Download in other formats:

Original Format