Context Navigation

readaerosol_interp.F90 @ 1279

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

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

Download in other formats:

Original Format