Context Navigation

readaerosol_interp.F90 @ 1200

Last change on this file since 1200 was 1183, checked in by jghattas, 15 years ago

Ajoute test de coherence entre les longitudes et latitudes des fichiers de forcage des aerosols et le model.

Corection pour la variable mass_solu_aero : changement du nom, commentaire et contenu(il y avait l'aerosol dust en trop)

Changement du nom de module car dans INCA le meme nom existe deja : aerosol_mod => aero_mod

File size: 18.3 KB

Line
1	! $Id$
2	!
3	SUBROUTINE readaerosol_interp(id_aero, r_day, first, pplay, paprs, mass_out, pi_mass_out)
4	!
5	! This routine will return the mass concentration at actual day(mass_out) and
6	! the pre-industrial values(pi_mass_out) for aerosol corresponding to "id_aero".
7	! The mass concentrations for all aerosols are saved in this routine but each
8	! call to this routine only treats the aerosol "id_aero".
9	!
10	! 1) Read in data for the whole year, only at first time step
11	! 2) Interpolate to the actual day, only at new day
12	! 3) Interpolate to the model vertical grid (target grid), only at new day
13	! 4) Test for negative mass values
14
15	USE dimphy, ONLY : klev,klon
16	USE mod_phys_lmdz_para, ONLY : mpi_rank
17	USE readaerosol_mod
18	USE aero_mod, ONLY : naero_spc, name_aero
19	USE write_field_phy
20
21	IMPLICIT NONE
22
23	INCLUDE "YOMCST.h"
24	INCLUDE "chem.h"
25	INCLUDE "temps.h"
26	INCLUDE "clesphys.h"
27	INCLUDE "iniprint.h"
28	INCLUDE "dimensions.h"
29	INCLUDE "comvert.h"
30	!
31	! Input:
32	!****************************************************************************************
33	INTEGER, INTENT(IN) :: id_aero! Identity number for the aerosol to treat
34	REAL, INTENT(IN) :: r_day ! Day of integration
35	LOGICAL, INTENT(IN) :: first ! First model timestep
36	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! pression at model mid-layers
37	REAL, DIMENSION(klon,klev+1),INTENT(IN):: paprs ! pression between model layers
38	!
39	! Output:
40	!****************************************************************************************
41	REAL, INTENT(OUT) :: mass_out(klon,klev) ! Mass of aerosol (monthly mean data,from file) [ug AIBCM/m3]
42	REAL, INTENT(OUT) :: pi_mass_out(klon,klev) ! Mass of preindustrial aerosol (monthly mean data,from file) [ug AIBCM/m3]
43	!
44	! Local Variables:
45	!****************************************************************************************
46	INTEGER :: i, k, ierr
47	INTEGER :: iday, iyr
48	INTEGER :: im, day1, day2, im2
49	INTEGER :: pi_klev_src ! Only for testing purpose
50	INTEGER, SAVE :: klev_src ! Number of vertical levles in source field
51	!$OMP THREADPRIVATE(klev_src)
52
53	REAL, DIMENSION(klon) :: psurf_day, pi_psurf_day
54	REAL, DIMENSION(klon) :: load_src, pi_load_src ! Mass load at source grid
55	REAL, DIMENSION(klon) :: load_tgt, load_tgt_test
56	REAL, DIMENSION(klon,klev) :: delp ! pressure difference in each model layer
57
58	REAL, ALLOCATABLE, DIMENSION(:,:) :: pplay_src ! pression mid-layer at source levels
59	REAL, ALLOCATABLE, DIMENSION(:,:) :: tmp1, tmp2 ! Temporary variables
60	REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: var_year ! VAR in right dimension for the total year
61	REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: pi_var_year ! pre-industrial VAR, -"-
62	!$OMP THREADPRIVATE(var_year,pi_var_year)
63	REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: var_day ! VAR interpolated to the actual day and model grid
64	REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: pi_var_day ! pre-industrial VAR, -"-
65	!$OMP THREADPRIVATE(var_day,pi_var_day)
66	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: psurf_year, pi_psurf_year ! surface pressure for the total year
67	!$OMP THREADPRIVATE(psurf_year, pi_psurf_year)
68	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: load_year, pi_load_year ! load in the column for the total year
69	!$OMP THREADPRIVATE(load_year, pi_load_year)
70
71	REAL, DIMENSION(:,:,:), POINTER :: pt_tmp ! Pointer allocated in readaerosol
72	REAL, POINTER, DIMENSION(:), SAVE :: pt_ap, pt_b ! Pointer for describing the vertical levels
73	!$OMP THREADPRIVATE(pt_ap, pt_b)
74
75	LOGICAL :: lnewday ! Indicates if first time step at a new day
76	LOGICAL,SAVE :: vert_interp ! Indicates if vertical interpolation will be done
77	LOGICAL,SAVE :: debug=.FALSE.! Debugging in this subroutine
78	!$OMP THREADPRIVATE(vert_interp, debug)
79
80
81	!****************************************************************************************
82	! Initialization
83	!
84	!****************************************************************************************
85
86	! Calculation to find if it is a new day
87	iday = INT(r_day)
88	iyr = iday/360
89	iday = iday-iyr*360 ! day of the actual year
90	iyr = iyr + annee_ref ! year of the run
91	im = iday/30 +1 ! the actual month
92
93	! 0.02 is about 0.5/24, namly less than half an hour
94	lnewday = (r_day-FLOAT(iday) < 0.02)
95
96	IF (.NOT. ALLOCATED(var_day)) THEN
97	ALLOCATE( var_day(klon, klev, naero_spc), stat=ierr)
98	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 1',1)
99	ALLOCATE( pi_var_day(klon, klev, naero_spc), stat=ierr)
100	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 2',1)
101
102	ALLOCATE( psurf_year(klon, 12, naero_spc), pi_psurf_year(klon, 12, naero_spc), stat=ierr)
103	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 3',1)
104
105	ALLOCATE( load_year(klon, 12, naero_spc), pi_load_year(klon, 12, naero_spc), stat=ierr)
106	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 4',1)
107
108	lnewday=.TRUE.
109
110	NULLIFY(pt_ap)
111	NULLIFY(pt_b)
112	END IF
113
114	!****************************************************************************************
115	! 1) Read in data : corresponding to the actual year and preindustrial data.
116	! Only for the first day of the year.
117	!
118	!****************************************************************************************
119	IF ( (first .OR. iday==1.) .AND. lnewday ) THEN
120	NULLIFY(pt_tmp)
121
122	! Reading values corresponding to the closest year taking into count the choice of aer_type.
123	! For aer_type=scenario interpolation between 2 data sets is done in readaerosol.
124	CALL readaerosol(name_aero(id_aero), aer_type, iyr, klev_src, pt_ap, pt_b, pt_tmp, &
125	psurf_year(:,:,id_aero), load_year(:,:,id_aero))
126	IF (.NOT. ALLOCATED(var_year)) THEN
127	ALLOCATE(var_year(klon, klev_src, 12, naero_spc), stat=ierr)
128	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 5',1)
129	END IF
130	var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
131
132	! Reading values corresponding to the preindustrial concentrations.
133	CALL readaerosol(name_aero(id_aero), 'preind', iyr, pi_klev_src, pt_ap, pt_b, pt_tmp, &
134	pi_psurf_year(:,:,id_aero), pi_load_year(:,:,id_aero))
135
136	! klev_src must be the same in both files.
137	! Also supposing pt_ap and pt_b to be the same in the 2 files without testing.
138	IF (pi_klev_src /= klev_src) THEN
139	WRITE(lunout,*) 'Error! All forcing files for the same aerosol must have the same vertical dimension'
140	WRITE(lunout,*) 'Aerosol : ', name_aero(id_aero)
141	CALL abort_gcm('readaerosol_interp','Differnt vertical axes in aerosol forcing files',1)
142	END IF
143
144	IF (.NOT. ALLOCATED(pi_var_year)) THEN
145	ALLOCATE(pi_var_year(klon, klev_src, 12, naero_spc), stat=ierr)
146	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 6',1)
147	END IF
148	pi_var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
149
150	IF (debug) THEN
151	CALL writefield_phy('psurf_src',psurf_year(:,:,id_aero),1)
152	CALL writefield_phy('pi_psurf_src',pi_psurf_year(:,:,id_aero),1)
153	CALL writefield_phy('load_year_src',load_year(:,:,id_aero),1)
154	CALL writefield_phy('pi_load_year_src',pi_load_year(:,:,id_aero),1)
155	END IF
156
157	! Pointer no more useful, deallocate.
158	DEALLOCATE(pt_tmp)
159
160	! Test if vertical interpolation will be needed.
161	IF (psurf_year(1,1,id_aero)==not_valid .OR. pi_psurf_year(1,1,id_aero)==not_valid ) THEN
162	! Pressure=not_valid indicates old file format, see module readaerosol
163	vert_interp = .FALSE.
164
165	! If old file format, both psurf_year and pi_psurf_year must be not_valid
166	IF ( psurf_year(1,1,id_aero) /= pi_psurf_year(1,1,id_aero) ) THEN
167	WRITE(lunout,*) 'Warning! All forcing files for the same aerosol must have the same structure'
168	CALL abort_gcm('readaerosol_interp', 'The aerosol files have not the same format',1)
169	END IF
170
171	IF (klev /= klev_src) THEN
172	WRITE(lunout,*) 'Old format of aerosol file do not allowed vertical interpolation'
173	CALL abort_gcm('readaerosol_interp', 'Old aerosol file not possible',1)
174	END IF
175
176	ELSE
177	vert_interp = .TRUE.
178	END IF
179
180	END IF ! IF ( (first .OR. iday==1.) .AND. lnewday ) THEN
181
182	!****************************************************************************************
183	! - 2) Interpolate to the actual day.
184	! - 3) Interpolate to the model vertical grid.
185	!
186	!****************************************************************************************
187
188	IF (lnewday) THEN ! only if new day
189	!****************************************************************************************
190	! 2) Interpolate to the actual day
191	!
192	!****************************************************************************************
193	! Find which months and days to use for time interpolation
194	IF (iday < im*30-15) THEN
195	! in the first half of the month use month before and actual month
196	im2=im-1
197	day2 = im2*30-15
198	day1 = im2*30+15
199	IF (im2 <= 0) THEN
200	! the month is january, thus the month before december
201	im2=12
202	END IF
203	ELSE
204	! the second half of the month
205	im2=im+1
206	IF (im2 > 12) THEN
207	! the month is december, the following thus january
208	im2=1
209	ENDIF
210	day2 = im*30-15
211	day1 = im*30+15
212	END IF
213
214	! Time interpolation, still on vertical source grid
215	ALLOCATE(tmp1(klon,klev_src), tmp2(klon,klev_src),stat=ierr)
216	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 7',1)
217
218	ALLOCATE(pplay_src(klon,klev_src), stat=ierr)
219	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 8',1)
220
221
222	DO k=1,klev_src
223	DO i=1,klon
224	tmp1(i,k) = &
225	var_year(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
226	(var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero))
227
228	tmp2(i,k) = &
229	pi_var_year(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
230	(pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero))
231	END DO
232	END DO
233
234	! Time interpolation for pressure at surface, still on vertical source grid
235	DO i=1,klon
236	psurf_day(i) = &
237	psurf_year(i,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
238	(psurf_year(i,im2,id_aero) - psurf_year(i,im,id_aero))
239
240	pi_psurf_day(i) = &
241	pi_psurf_year(i,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
242	(pi_psurf_year(i,im2,id_aero) - pi_psurf_year(i,im,id_aero))
243	END DO
244
245	! Time interpolation for the load, still on vertical source grid
246	DO i=1,klon
247	load_src(i) = &
248	load_year(i,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
249	(load_year(i,im2,id_aero) - load_year(i,im,id_aero))
250
251	pi_load_src(i) = &
252	pi_load_year(i,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
253	(pi_load_year(i,im2,id_aero) - pi_load_year(i,im,id_aero))
254	END DO
255
256	!****************************************************************************************
257	! 3) Interpolate to the model vertical grid (target grid)
258	!
259	!****************************************************************************************
260
261	IF (vert_interp) THEN
262
263	! - Interpolate variable tmp1 (on source grid) to var_day (on target grid)
264	!********************************************************************************
265	! a) calculate pression at vertical levels for the source grid using the
266	! hybrid-sigma coordinates ap and b and the surface pressure, variables from file.
267	DO k = 1, klev_src
268	DO i = 1, klon
269	pplay_src(i,k)= pt_ap(k) + pt_b(k)*psurf_day(i)
270	END DO
271	END DO
272
273	IF (debug) THEN
274	CALL writefield_phy('psurf_day_src',psurf_day(:),1)
275	CALL writefield_phy('pplay_src',pplay_src(:,:),klev_src)
276	CALL writefield_phy('pplay',pplay(:,:),klev)
277	CALL writefield_phy('day_src',tmp1,klev_src)
278	CALL writefield_phy('pi_day_src',tmp2,klev_src)
279	END IF
280
281	! b) vertical interpolation on pressure leveles
282	CALL pres2lev(tmp1(:,:), var_day(:,:,id_aero), klev_src, klev, pplay_src, pplay, &
283	1, klon, .FALSE.)
284
285	IF (debug) CALL writefield_phy('day_tgt',var_day(:,:,id_aero),klev)
286
287	! c) adjust to conserve total aerosol mass load in the vertical pillar
288	! Calculate the load in the actual pillar and compare with the load
289	! read from aerosol file.
290
291	! Find the pressure difference in each model layer
292	DO k = 1, klev
293	DO i = 1, klon
294	delp(i,k) = paprs(i,k) - paprs (i,k+1)
295	END DO
296	END DO
297
298	! Find the mass load in the actual pillar, on target grid
299	load_tgt(:) = 0.
300	DO k= 1, klev
301	DO i = 1, klon
302	load_tgt(i) = load_tgt(i) + 1/RG * var_day(i,k,id_aero)*delp(i,k)
303	END DO
304	END DO
305
306	! Adjust, uniform
307	DO k = 1, klev
308	DO i = 1, klon
309	var_day(i,k,id_aero) = var_day(i,k,id_aero)*load_src(i)/load_tgt(i)
310	END DO
311	END DO
312
313	IF (debug) THEN
314	load_tgt_test(:) = 0.
315	DO k= 1, klev
316	DO i = 1, klon
317	load_tgt_test(i) = load_tgt_test(i) + 1/RG * var_day(i,k,id_aero)*delp(i,k)
318	END DO
319	END DO
320
321	CALL writefield_phy('day_tgt2',var_day(:,:,id_aero),klev)
322	CALL writefield_phy('load_tgt',load_tgt(:),1)
323	CALL writefield_phy('load_tgt_test',load_tgt_test(:),1)
324	CALL writefield_phy('load_src',load_src(:),1)
325	END IF
326
327	! - Interpolate variable tmp2 (source grid) to pi_var_day (target grid)
328	!********************************************************************************
329	! a) calculate pression at vertical levels at source grid
330	DO k = 1, klev_src
331	DO i = 1, klon
332	pplay_src(i,k)= pt_ap(k) + pt_b(k)*pi_psurf_day(i)
333	END DO
334	END DO
335
336	IF (debug) THEN
337	CALL writefield_phy('pi_psurf_day_src',pi_psurf_day(:),1)
338	CALL writefield_phy('pi_pplay_src',pplay_src(:,:),klev_src)
339	END IF
340
341	! b) vertical interpolation on pressure leveles
342	CALL pres2lev(tmp2(:,:), pi_var_day(:,:,id_aero), klev_src, klev, pplay_src, pplay, &
343	1, klon, .FALSE.)
344
345	IF (debug) CALL writefield_phy('pi_day_tgt',pi_var_day(:,:,id_aero),klev)
346
347	! c) adjust to conserve total aerosol mass load in the vertical pillar
348	! Calculate the load in the actual pillar and compare with the load
349	! read from aerosol file.
350
351	! Find the load in the actual pillar, on target grid
352	load_tgt(:) = 0.
353	DO k = 1, klev
354	DO i = 1, klon
355	load_tgt(i) = load_tgt(i) + 1/RG * pi_var_day(i,k,id_aero)*delp(i,k)
356	END DO
357	END DO
358
359	DO k = 1, klev
360	DO i = 1, klon
361	pi_var_day(i,k,id_aero) = pi_var_day(i,k,id_aero)*pi_load_src(i)/load_tgt(i)
362	END DO
363	END DO
364
365	IF (debug) THEN
366	load_tgt_test(:) = 0.
367	DO k = 1, klev
368	DO i = 1, klon
369	load_tgt_test(i) = load_tgt_test(i) + 1/RG * pi_var_day(i,k,id_aero)*delp(i,k)
370	END DO
371	END DO
372	CALL writefield_phy('pi_day_tgt2',pi_var_day(:,:,id_aero),klev)
373	CALL writefield_phy('pi_load_tgt',load_tgt(:),1)
374	CALL writefield_phy('pi_load_tgt_test',load_tgt_test(:),1)
375	CALL writefield_phy('pi_load_src',pi_load_src(:),1)
376	END IF
377
378
379	ELSE ! No vertical interpolation done
380
381	var_day(:,:,id_aero) = tmp1(:,:)
382	pi_var_day(:,:,id_aero) = tmp2(:,:)
383
384	END IF ! vert_interp
385
386
387	! Deallocation
388	DEALLOCATE(tmp1, tmp2, pplay_src, stat=ierr)
389
390	!****************************************************************************************
391	! 4) Test for negative mass values
392	!
393	!****************************************************************************************
394	IF (MINVAL(var_day(:,:,id_aero)) < 0.) THEN
395	DO k=1,klev
396	DO i=1,klon
397	! Test for var_day
398	IF (var_day(i,k,id_aero) < 0.) THEN
399	IF (iday-day2 < 0.) WRITE(lunout,*) 'iday-day2=',iday-day2
400	IF (var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero) < 0.) THEN
401	WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', &
402	trim(name_aero(id_aero)),'(i,k,im)=', &
403	var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero)
404	END IF
405
406	WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
407	CALL abort_gcm('readaerosol_interp','Error in interpolation 1',1)
408	END IF
409	END DO
410	END DO
411	END IF
412
413	IF (MINVAL(pi_var_day(:,:,id_aero)) < 0. ) THEN
414	DO k=1, klev
415	DO i=1,klon
416	! Test for pi_var_day
417	IF (pi_var_day(i,k,id_aero) < 0.) THEN
418	IF (iday-day2 < 0.) WRITE(lunout,*) 'iday-day2=',iday-day2
419	IF (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero) < 0.) THEN
420	WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', &
421	trim(name_aero(id_aero)),'(i,k,im)=', &
422	pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero)
423	END IF
424
425	WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
426	CALL abort_gcm('readaerosol_interp','Error in interpolation 2',1)
427	END IF
428	END DO
429	END DO
430	END IF
431
432	END IF ! lnewday
433
434	!****************************************************************************************
435	! Copy output from saved variables
436	!
437	!****************************************************************************************
438
439	mass_out(:,:) = var_day(:,:,id_aero)
440	pi_mass_out(:,:) = pi_var_day(:,:,id_aero)
441
442	END SUBROUTINE readaerosol_interp

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

Download in other formats:

Original Format