Context Navigation

readaerosol_interp.F90 @ 1265

Last change on this file since 1265 was 1265, checked in by Laurent Fairhead, 15 years ago
Modifications pour utiliser le calendrier réaliste et les nouveaux fichiers a 14 valeurs mensuelles
File size: 23.5 KB

Line
1	! $Id$
2	!
3	SUBROUTINE readaerosol_interp(id_aero, itap, pdtphys, r_day, first, pplay, paprs, t_seri, 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 ioipsl
16	USE dimphy, ONLY : klev,klon
17	USE mod_phys_lmdz_para, ONLY : mpi_rank
18	USE readaerosol_mod
19	USE aero_mod, ONLY : naero_spc, name_aero
20	USE write_field_phy
21	USE phys_cal_mod
22
23	IMPLICIT NONE
24
25	INCLUDE "YOMCST.h"
26	INCLUDE "chem.h"
27	INCLUDE "temps.h"
28	INCLUDE "clesphys.h"
29	INCLUDE "iniprint.h"
30	INCLUDE "dimensions.h"
31	INCLUDE "comvert.h"
32	!
33	! Input:
34	!****************************************************************************************
35	INTEGER, INTENT(IN) :: id_aero! Identity number for the aerosol to treat
36	INTEGER, INTENT(IN) :: itap ! Physic step count
37	REAL, INTENT(IN) :: pdtphys! Physic day step
38	REAL, INTENT(IN) :: r_day ! Day of integration
39	LOGICAL, INTENT(IN) :: first ! First model timestep
40	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! pression at model mid-layers
41	REAL, DIMENSION(klon,klev+1),INTENT(IN):: paprs ! pression between model layers
42	REAL, DIMENSION(klon,klev), INTENT(IN) :: t_seri ! air temperature
43	!
44	! Output:
45	!****************************************************************************************
46	REAL, INTENT(OUT) :: mass_out(klon,klev) ! Mass of aerosol (monthly mean data,from file) [ug AIBCM/m3]
47	REAL, INTENT(OUT) :: pi_mass_out(klon,klev) ! Mass of preindustrial aerosol (monthly mean data,from file) [ug AIBCM/m3]
48	!
49	! Local Variables:
50	!****************************************************************************************
51	INTEGER :: i, k, ierr
52	INTEGER :: iday, iyr, lmt_pas
53	! INTEGER :: im, day1, day2, im2
54	INTEGER :: im, im2
55	REAL :: day1, day2
56	INTEGER :: pi_klev_src ! Only for testing purpose
57	INTEGER, SAVE :: klev_src ! Number of vertical levles in source field
58	!$OMP THREADPRIVATE(klev_src)
59
60	REAL :: zrho ! Air density [kg/m3]
61	REAL :: volm ! Volyme de melange [kg/kg]
62	REAL, DIMENSION(klon) :: psurf_day, pi_psurf_day
63	REAL, DIMENSION(klon) :: load_src, pi_load_src ! Mass load at source grid
64	REAL, DIMENSION(klon) :: load_tgt, load_tgt_test
65	REAL, DIMENSION(klon,klev) :: delp ! pressure difference in each model layer
66
67	REAL, ALLOCATABLE, DIMENSION(:,:) :: pplay_src ! pression mid-layer at source levels
68	REAL, ALLOCATABLE, DIMENSION(:,:) :: tmp1, tmp2 ! Temporary variables
69	REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: var_year ! VAR in right dimension for the total year
70	REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: pi_var_year ! pre-industrial VAR, -"-
71	!$OMP THREADPRIVATE(var_year,pi_var_year)
72	REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: var_day ! VAR interpolated to the actual day and model grid
73	REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: pi_var_day ! pre-industrial VAR, -"-
74	!$OMP THREADPRIVATE(var_day,pi_var_day)
75	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: psurf_year, pi_psurf_year ! surface pressure for the total year
76	!$OMP THREADPRIVATE(psurf_year, pi_psurf_year)
77	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: load_year, pi_load_year ! load in the column for the total year
78	!$OMP THREADPRIVATE(load_year, pi_load_year)
79
80	REAL, DIMENSION(:,:,:), POINTER :: pt_tmp ! Pointer allocated in readaerosol
81	REAL, DIMENSION(:,:), POINTER :: pt_tmp_surf ! Pointer allocated in readaerosol
82	REAL, DIMENSION(:,:), POINTER :: pt_tmp_load ! Pointer allocated in readaerosol
83	REAL, DIMENSION(:,:), POINTER :: pt_tmp_pi_surf ! Pointer allocated in readaerosol
84	REAL, DIMENSION(:,:), POINTER :: pt_tmp_pi_load ! Pointer allocated in readaerosol
85	REAL, POINTER, DIMENSION(:), SAVE :: pt_ap, pt_b ! Pointer for describing the vertical levels
86	!$OMP THREADPRIVATE(pt_ap, pt_b)
87	INTEGER, SAVE :: nbr_tsteps ! number of time steps in file read
88	REAL, DIMENSION(14), SAVE :: month_len, month_start, month_mid
89	!$OMP THREADPRIVATE(month_len, month_start, month_mid)
90	REAL :: jDay
91
92	!$OMP THREADPRIVATE(nbr_tsteps)
93	LOGICAL :: lnewday ! Indicates if first time step at a new day
94	LOGICAL :: OLDNEWDAY
95	LOGICAL,SAVE :: vert_interp ! Indicates if vertical interpolation will be done
96	LOGICAL,SAVE :: debug=.FALSE.! Debugging in this subroutine
97	!$OMP THREADPRIVATE(vert_interp, debug)
98
99
100	!****************************************************************************************
101	! Initialization
102	!
103	!****************************************************************************************
104
105	! Calculation to find if it is a new day
106
107	IF(mpi_rank == 0 .AND. debug )then
108	PRINT*,'CONTROL PANEL REGARDING TIME STEPING'
109	ENDIF
110
111	! Use phys_cal_mod
112	! iday= day_cur
113	! iyr = year_cur
114	! im = mth_cur
115
116	iday = INT(r_day)
117	iyr = iday/360
118	iday = iday-iyr*360 ! day of the actual year
119	iyr = iyr + annee_ref ! year of the run
120	im = iday/30 +1 ! the actual month
121	CALL ymds2ju(iyr, im, iday, 0., jDay)
122
123	IF(MOD(itap-1,NINT(86400./pdtphys)) == 0)THEN
124	lnewday=.TRUE.
125	ELSE
126	lnewday=.FALSE.
127	ENDIF
128
129	IF(mpi_rank == 0 .AND. debug)then
130	! 0.02 is about 0.5/24, namly less than half an hour
131	OLDNEWDAY = (r_day-FLOAT(iday) < 0.02)
132	! Once per day, update aerosol fields
133	lmt_pas = NINT(86400./pdtphys)
134	PRINT*,'r_day-FLOAT(iday) =',r_day-FLOAT(iday)
135	PRINT*,'itap =',itap
136	PRINT*,'pdtphys =',pdtphys
137	PRINT*,'lmt_pas =',lmt_pas
138	PRINT*,'iday =',iday
139	PRINT*,'r_day =',r_day
140	PRINT*,'day_cur =',day_cur
141	PRINT*,'mth_cur =',mth_cur
142	PRINT*,'year_cur =',year_cur
143	PRINT*,'NINT(86400./pdtphys) =',NINT(86400./pdtphys)
144	PRINT*,'MOD(0,1) =',MOD(0,1)
145	PRINT*,'lnewday =',lnewday
146	PRINT*,'OLDNEWDAY =',OLDNEWDAY
147	ENDIF
148
149	IF (.NOT. ALLOCATED(var_day)) THEN
150	ALLOCATE( var_day(klon, klev, naero_spc), stat=ierr)
151	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 1',1)
152	ALLOCATE( pi_var_day(klon, klev, naero_spc), stat=ierr)
153	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 2',1)
154
155	! ALLOCATE( psurf_year(klon, 12, naero_spc), pi_psurf_year(klon, 12, naero_spc), stat=ierr)
156	! IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 3',1)
157	!
158	! ALLOCATE( load_year(klon, 12, naero_spc), pi_load_year(klon, 12, naero_spc), stat=ierr)
159	! IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 4',1)
160
161	lnewday=.TRUE.
162
163	NULLIFY(pt_ap)
164	NULLIFY(pt_b)
165	END IF
166
167	!****************************************************************************************
168	! 1) Read in data : corresponding to the actual year and preindustrial data.
169	! Only for the first day of the year.
170	!
171	!****************************************************************************************
172	IF ( (first .OR. iday==0) .AND. lnewday ) THEN
173	NULLIFY(pt_tmp)
174	NULLIFY(pt_tmp_surf)
175	NULLIFY(pt_tmp_pi_surf)
176	NULLIFY(pt_tmp_load)
177	NULLIFY(pt_tmp_pi_load)
178
179	! Reading values corresponding to the closest year taking into count the choice of aer_type.
180	! For aer_type=scenario interpolation between 2 data sets is done in readaerosol.
181	CALL readaerosol(name_aero(id_aero), aer_type, iyr, klev_src, pt_ap, pt_b, pt_tmp, &
182	pt_tmp_surf,pt_tmp_load,nbr_tsteps)
183	IF (.NOT. ALLOCATED(var_year)) THEN
184	ALLOCATE(var_year(klon, klev_src, nbr_tsteps, naero_spc), stat=ierr)
185	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 5.1',1)
186	ALLOCATE( psurf_year(klon, nbr_tsteps, naero_spc), load_year(klon, nbr_tsteps, naero_spc), stat=ierr)
187	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 5.2',1)
188	END IF
189	var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
190	psurf_year(:,:,id_aero) = pt_tmp_surf(:,:)
191	load_year(:,:,id_aero) = pt_tmp_load(:,:)
192
193	! Reading values corresponding to the preindustrial concentrations.
194	CALL readaerosol(name_aero(id_aero), 'preind', iyr, pi_klev_src, pt_ap, pt_b, pt_tmp, &
195	pt_tmp_surf,pt_tmp_load,nbr_tsteps)
196
197	! klev_src must be the same in both files.
198	! Also supposing pt_ap and pt_b to be the same in the 2 files without testing.
199	IF (pi_klev_src /= klev_src) THEN
200	WRITE(lunout,*) 'Error! All forcing files for the same aerosol must have the same vertical dimension'
201	WRITE(lunout,*) 'Aerosol : ', name_aero(id_aero)
202	CALL abort_gcm('readaerosol_interp','Differnt vertical axes in aerosol forcing files',1)
203	END IF
204
205	IF (.NOT. ALLOCATED(pi_var_year)) THEN
206	ALLOCATE(pi_var_year(klon, klev_src, nbr_tsteps, naero_spc), stat=ierr)
207	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 6.1',1)
208	ALLOCATE( pi_psurf_year(klon, nbr_tsteps, naero_spc), pi_load_year(klon, nbr_tsteps, naero_spc), stat=ierr)
209	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 6.2',1)
210	END IF
211	pi_var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
212	pi_psurf_year(:,:,id_aero) = pt_tmp_surf(:,:)
213	pi_load_year(:,:,id_aero) = pt_tmp_load(:,:)
214
215	IF (debug) THEN
216	CALL writefield_phy('var_year_first',var_year(:,:,1,id_aero),klev_src)
217	CALL writefield_phy('var_year_last',var_year(:,:,nbr_tsteps,id_aero),klev_src)
218	CALL writefield_phy('psurf_src',psurf_year(:,:,id_aero),1)
219	CALL writefield_phy('pi_psurf_src',pi_psurf_year(:,:,id_aero),1)
220	CALL writefield_phy('load_year_src',load_year(:,:,id_aero),1)
221	CALL writefield_phy('pi_load_year_src',pi_load_year(:,:,id_aero),1)
222	END IF
223
224	! Pointer no more useful, deallocate.
225	DEALLOCATE(pt_tmp)
226	DEALLOCATE(pt_tmp_surf)
227	DEALLOCATE(pt_tmp_load)
228
229	! Test if vertical interpolation will be needed.
230	IF (psurf_year(1,1,id_aero)==not_valid .OR. pi_psurf_year(1,1,id_aero)==not_valid ) THEN
231	! Pressure=not_valid indicates old file format, see module readaerosol
232	vert_interp = .FALSE.
233
234	! If old file format, both psurf_year and pi_psurf_year must be not_valid
235	IF ( psurf_year(1,1,id_aero) /= pi_psurf_year(1,1,id_aero) ) THEN
236	WRITE(lunout,*) 'Warning! All forcing files for the same aerosol must have the same structure'
237	CALL abort_gcm('readaerosol_interp', 'The aerosol files have not the same format',1)
238	END IF
239
240	IF (klev /= klev_src) THEN
241	WRITE(lunout,*) 'Old format of aerosol file do not allowed vertical interpolation'
242	CALL abort_gcm('readaerosol_interp', 'Old aerosol file not possible',1)
243	END IF
244
245	ELSE
246	vert_interp = .TRUE.
247	END IF
248
249	! Calendar initialisation
250	!
251	DO i = 2, 13
252	month_len(i) = float(ioget_mon_len(year_cur, i-1))
253	CALL ymds2ju(year_cur, i-1, 1, 0.0, month_start(i))
254	ENDDO
255	month_len(1) = float(ioget_mon_len(year_cur-1, 12))
256	CALL ymds2ju(year_cur-1, 12, 1, 0.0, month_start(1))
257	month_len(14) = float(ioget_mon_len(year_cur+1, 1))
258	CALL ymds2ju(year_cur+1, 1, 1, 0.0, month_start(14))
259	month_mid(:) = month_start (:) + month_len(:)/2.
260
261	write(55,*)'month_len = ',month_len
262	write(55,*)'month_start = ',month_start
263	write(55,*)'month_mid = ',month_mid
264
265	END IF ! IF ( (first .OR. iday==0) .AND. lnewday ) THEN
266
267	!****************************************************************************************
268	! - 2) Interpolate to the actual day.
269	! - 3) Interpolate to the model vertical grid.
270	!
271	!****************************************************************************************
272
273	IF (lnewday) THEN ! only if new day
274	!****************************************************************************************
275	! 2) Interpolate to the actual day
276	!
277	!****************************************************************************************
278	! Find which months and days to use for time interpolation
279	IF (nbr_tsteps == 12) then
280	IF (jDay < month_mid(im+1)) THEN
281	im2=im-1
282	day2 = month_mid(im2+1)
283	day1 = month_mid(im+1)
284	IF (im2 <= 0) THEN
285	! the month is january, thus the month before december
286	im2=12
287	END IF
288	ELSE
289	! the second half of the month
290	im2=im+1
291	day2 = month_mid(im+1)
292	day1 = month_mid(im2+1)
293	IF (im2 > 12) THEN
294	! the month is december, the following thus january
295	im2=1
296	ENDIF
297	END IF
298	ELSE IF (nbr_tsteps == 14) then
299	im = im + 1
300	IF (jDay < month_mid(im)) THEN
301	! in the first half of the month use month before and actual month
302	im2=im-1
303	day2 = month_mid(im2)
304	day1 = month_mid(im)
305	ELSE
306	! the second half of the month
307	im2=im+1
308	day2 = month_mid(im)
309	day1 = month_mid(im2)
310	END IF
311	ELSE
312	CALL abort_gcm('readaerosol_interp', 'number of months undefined',1)
313	ENDIF
314
315	! ! Find which months and days to use for time interpolation
316	! IF (iday < im*30-15) THEN
317	! ! in the first half of the month use month before and actual month
318	! im2=im-1
319	! day2 = im2*30-15
320	! day1 = im2*30+15
321	! IF (im2 <= 0) THEN
322	! ! the month is january, thus the month before december
323	! im2=12
324	! END IF
325	! ELSE
326	! ! the second half of the month
327	! im2=im+1
328	! IF (im2 > 12) THEN
329	! ! the month is december, the following thus january
330	! im2=1
331	! ENDIF
332	! day2 = im*30-15
333	! day1 = im*30+15
334	! END IF
335	! jDay = jDay+1
336	! write(55,*)'iday, jDay, im, im2, day1, day2 = ',iday, jDay, im, im2, day1, day2, nbr_tsteps
337
338	! Time interpolation, still on vertical source grid
339	ALLOCATE(tmp1(klon,klev_src), tmp2(klon,klev_src),stat=ierr)
340	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 7',1)
341
342	ALLOCATE(pplay_src(klon,klev_src), stat=ierr)
343	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 8',1)
344
345
346	DO k=1,klev_src
347	DO i=1,klon
348	tmp1(i,k) = &
349	var_year(i,k,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * &
350	(var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero))
351
352	tmp2(i,k) = &
353	pi_var_year(i,k,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * &
354	(pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero))
355	END DO
356	END DO
357
358	! Time interpolation for pressure at surface, still on vertical source grid
359	DO i=1,klon
360	psurf_day(i) = &
361	psurf_year(i,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * &
362	(psurf_year(i,im2,id_aero) - psurf_year(i,im,id_aero))
363
364	pi_psurf_day(i) = &
365	pi_psurf_year(i,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * &
366	(pi_psurf_year(i,im2,id_aero) - pi_psurf_year(i,im,id_aero))
367	END DO
368
369	! Time interpolation for the load, still on vertical source grid
370	DO i=1,klon
371	load_src(i) = &
372	load_year(i,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * &
373	(load_year(i,im2,id_aero) - load_year(i,im,id_aero))
374
375	pi_load_src(i) = &
376	pi_load_year(i,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * &
377	(pi_load_year(i,im2,id_aero) - pi_load_year(i,im,id_aero))
378	END DO
379
380	!****************************************************************************************
381	! 3) Interpolate to the model vertical grid (target grid)
382	!
383	!****************************************************************************************
384
385	IF (vert_interp) THEN
386
387	! - Interpolate variable tmp1 (on source grid) to var_day (on target grid)
388	!********************************************************************************
389	! a) calculate pression at vertical levels for the source grid using the
390	! hybrid-sigma coordinates ap and b and the surface pressure, variables from file.
391	DO k = 1, klev_src
392	DO i = 1, klon
393	pplay_src(i,k)= pt_ap(k) + pt_b(k)*psurf_day(i)
394	END DO
395	END DO
396
397	IF (debug) THEN
398	CALL writefield_phy('psurf_day_src',psurf_day(:),1)
399	CALL writefield_phy('pplay_src',pplay_src(:,:),klev_src)
400	CALL writefield_phy('pplay',pplay(:,:),klev)
401	CALL writefield_phy('day_src',tmp1,klev_src)
402	CALL writefield_phy('pi_day_src',tmp2,klev_src)
403	END IF
404
405	! b) vertical interpolation on pressure leveles
406	CALL pres2lev(tmp1(:,:), var_day(:,:,id_aero), klev_src, klev, pplay_src, pplay, &
407	1, klon, .FALSE.)
408
409	IF (debug) CALL writefield_phy('day_tgt',var_day(:,:,id_aero),klev)
410
411	! c) adjust to conserve total aerosol mass load in the vertical pillar
412	! Calculate the load in the actual pillar and compare with the load
413	! read from aerosol file.
414
415	! Find the pressure difference in each model layer
416	DO k = 1, klev
417	DO i = 1, klon
418	delp(i,k) = paprs(i,k) - paprs (i,k+1)
419	END DO
420	END DO
421
422	! Find the mass load in the actual pillar, on target grid
423	load_tgt(:) = 0.
424	DO k= 1, klev
425	DO i = 1, klon
426	zrho = pplay(i,k)/t_seri(i,k)/RD ! [kg/m3]
427	volm = var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
428	load_tgt(i) = load_tgt(i) + 1/RG * volm *delp(i,k)
429	END DO
430	END DO
431
432	! Adjust, uniform
433	DO k = 1, klev
434	DO i = 1, klon
435	var_day(i,k,id_aero) = var_day(i,k,id_aero)*load_src(i)/load_tgt(i)
436	END DO
437	END DO
438
439	IF (debug) THEN
440	load_tgt_test(:) = 0.
441	DO k= 1, klev
442	DO i = 1, klon
443	zrho = pplay(i,k)/t_seri(i,k)/RD ! [kg/m3]
444	volm = var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
445	load_tgt_test(i) = load_tgt_test(i) + 1/RG * volm*delp(i,k)
446	END DO
447	END DO
448
449	CALL writefield_phy('day_tgt2',var_day(:,:,id_aero),klev)
450	CALL writefield_phy('load_tgt',load_tgt(:),1)
451	CALL writefield_phy('load_tgt_test',load_tgt_test(:),1)
452	CALL writefield_phy('load_src',load_src(:),1)
453	END IF
454
455	! - Interpolate variable tmp2 (source grid) to pi_var_day (target grid)
456	!********************************************************************************
457	! a) calculate pression at vertical levels at source grid
458	DO k = 1, klev_src
459	DO i = 1, klon
460	pplay_src(i,k)= pt_ap(k) + pt_b(k)*pi_psurf_day(i)
461	END DO
462	END DO
463
464	IF (debug) THEN
465	CALL writefield_phy('pi_psurf_day_src',pi_psurf_day(:),1)
466	CALL writefield_phy('pi_pplay_src',pplay_src(:,:),klev_src)
467	END IF
468
469	! b) vertical interpolation on pressure leveles
470	CALL pres2lev(tmp2(:,:), pi_var_day(:,:,id_aero), klev_src, klev, pplay_src, pplay, &
471	1, klon, .FALSE.)
472
473	IF (debug) CALL writefield_phy('pi_day_tgt',pi_var_day(:,:,id_aero),klev)
474
475	! c) adjust to conserve total aerosol mass load in the vertical pillar
476	! Calculate the load in the actual pillar and compare with the load
477	! read from aerosol file.
478
479	! Find the load in the actual pillar, on target grid
480	load_tgt(:) = 0.
481	DO k = 1, klev
482	DO i = 1, klon
483	zrho = pplay(i,k)/t_seri(i,k)/RD ! [kg/m3]
484	volm = pi_var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
485	load_tgt(i) = load_tgt(i) + 1/RG * volm * delp(i,k)
486	END DO
487	END DO
488
489	DO k = 1, klev
490	DO i = 1, klon
491	pi_var_day(i,k,id_aero) = pi_var_day(i,k,id_aero)*pi_load_src(i)/load_tgt(i)
492	END DO
493	END DO
494
495	IF (debug) THEN
496	load_tgt_test(:) = 0.
497	DO k = 1, klev
498	DO i = 1, klon
499	zrho = pplay(i,k)/t_seri(i,k)/RD ! [kg/m3]
500	volm = pi_var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
501	load_tgt_test(i) = load_tgt_test(i) + 1/RG * volm * delp(i,k)
502	END DO
503	END DO
504	CALL writefield_phy('pi_day_tgt2',pi_var_day(:,:,id_aero),klev)
505	CALL writefield_phy('pi_load_tgt',load_tgt(:),1)
506	CALL writefield_phy('pi_load_tgt_test',load_tgt_test(:),1)
507	CALL writefield_phy('pi_load_src',pi_load_src(:),1)
508	END IF
509
510
511	ELSE ! No vertical interpolation done
512
513	var_day(:,:,id_aero) = tmp1(:,:)
514	pi_var_day(:,:,id_aero) = tmp2(:,:)
515
516	END IF ! vert_interp
517
518
519	! Deallocation
520	DEALLOCATE(tmp1, tmp2, pplay_src, stat=ierr)
521
522	!****************************************************************************************
523	! 4) Test for negative mass values
524	!
525	!****************************************************************************************
526	IF (MINVAL(var_day(:,:,id_aero)) < 0.) THEN
527	DO k=1,klev
528	DO i=1,klon
529	! Test for var_day
530	IF (var_day(i,k,id_aero) < 0.) THEN
531	IF (jDay-day2 < 0.) WRITE(lunout,*) 'jDay-day2=',jDay-day2
532	IF (var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero) < 0.) THEN
533	WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', &
534	trim(name_aero(id_aero)),'(i,k,im)=', &
535	var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero)
536	END IF
537
538	WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
539	CALL abort_gcm('readaerosol_interp','Error in interpolation 1',1)
540	END IF
541	END DO
542	END DO
543	END IF
544
545	IF (MINVAL(pi_var_day(:,:,id_aero)) < 0. ) THEN
546	DO k=1, klev
547	DO i=1,klon
548	! Test for pi_var_day
549	IF (pi_var_day(i,k,id_aero) < 0.) THEN
550	IF (jDay-day2 < 0.) WRITE(lunout,*) 'jDay-day2=',jDay-day2
551	IF (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero) < 0.) THEN
552	WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', &
553	trim(name_aero(id_aero)),'(i,k,im)=', &
554	pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero)
555	END IF
556
557	WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
558	CALL abort_gcm('readaerosol_interp','Error in interpolation 2',1)
559	END IF
560	END DO
561	END DO
562	END IF
563
564	END IF ! lnewday
565
566	!****************************************************************************************
567	! Copy output from saved variables
568	!
569	!****************************************************************************************
570
571	mass_out(:,:) = var_day(:,:,id_aero)
572	pi_mass_out(:,:) = pi_var_day(:,:,id_aero)
573
574	END SUBROUTINE readaerosol_interp

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

Download in other formats:

Original Format