Context Navigation

readaerosol_interp.F90 @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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

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

Download in other formats:

Original Format