Context Navigation

readaerosol_interp.F90 @ 1373

Last change on this file since 1373 was 1347, checked in by Laurent Fairhead, 15 years ago

Additions to aerosol outputs for CMIP5 exercise
(Needed because of chageset r1346 LF)

Additions aux sorties aérosols pour l'exercice CMIP5
(Nécessaires suite au changeset r1346 LF)

Michael, Anne

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

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

Download in other formats:

Original Format