Context Navigation

readaerosol_interp.F90 @ 1271

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

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

Download in other formats:

Original Format