Context Navigation

readaerosol_interp.F90 @ 3814

Last change on this file since 3814 was 3814, checked in by ymipsl, 10 years ago

remove all dynamic dependency in LMDZ physics except for the include "dimensions.h"

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

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

Download in other formats:

Original Format