Context Navigation

readaerosol_interp.F90 @ 5185

Last change on this file since 5185 was 2840, checked in by oboucher, 8 years ago
Correcting the case when aerosol is zero (eg NO3 when missing from file)
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 23.3 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	USE print_control_mod, ONLY: lunout
24
25	IMPLICIT NONE
26
27	INCLUDE "YOMCST.h"
28	INCLUDE "chem.h"
29	INCLUDE "clesphys.h"
30
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	ENDIF
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	ENDIF
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	ENDIF
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	ENDIF
208	ELSE
209	CALL abort_physic('readaerosol_interp', 'this aer_type not supported',1)
210	ENDIF
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	ENDIF
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	ENDIF
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	ENDIF
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	ENDIF
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	ENDIF
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	ENDIF
266
267	ELSE
268	vert_interp = .TRUE.
269	ENDIF
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	ENDIF ! 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	ENDIF
312	ELSE
313	! the second half of the month
314	im2=im+1
315	day1 = month_mid(im+1)
316	day2 = month_mid(im2+1)
317	IF (im2 > 12) THEN
318	! the month is december, the following thus january
319	im2=1
320	ENDIF
321	ENDIF
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	day1 = month_mid(im)
333	day2 = month_mid(im2)
334	ENDIF
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	ENDDO
361	ENDDO
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	ENDDO
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	ENDDO
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	ENDDO
400	ENDDO
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	ENDIF
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	ENDDO
425	ENDDO
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) + volm *delp(i,k)/RG
434	ENDDO
435	ENDDO
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)/max(1.e-30,load_tgt(i))
441	ENDDO
442	ENDDO
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) + volm*delp(i,k)/RG
451	ENDDO
452	ENDDO
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	ENDIF
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	ENDDO
467	ENDDO
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	ENDIF
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) + volm*delp(i,k)/RG
491	ENDDO
492	ENDDO
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)/max(1.e-30,load_tgt(i))
497	ENDDO
498	ENDDO
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) + volm*delp(i,k)/RG
507	ENDDO
508	ENDDO
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	ENDIF
514
515
516	ELSE ! No vertical interpolation done
517
518	var_day(:,:,id_aero) = tmp1(:,:)
519	pi_var_day(:,:,id_aero) = tmp2(:,:)
520
521	ENDIF ! 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	ENDIF
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	ENDIF
546	ENDDO
547	ENDDO
548	ENDIF
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	ENDIF
561
562	WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
563	CALL abort_physic('readaerosol_interp','Error in interpolation 2',1)
564	ENDIF
565	ENDDO
566	ENDDO
567	ENDIF
568
569	ENDIF ! 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