Context Navigation

readaerosol_interp.F90 @ 3241

Last change on this file since 3241 was 1273, checked in by Laurent Fairhead, 15 years ago
Correction d'un bug sur le calcul de la date courante pour la lecture des aerosols Attention, les resultats sont modifies suite a un decalage d'une journee dans l'ancien calcul de la date LF
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