Context Navigation

readaerosol_interp.F90 @ 5442

Last change on this file since 5442 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

IPSLCM5A configuration (comparison between trunk and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

IPSLCM5A configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

MH07 configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

pour une config IPSLCM5A (comparaison tronc/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux (à part sur RN et Pb)
- fichiers histoire égaux

pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux
- fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

pour une config IPSLCM5A:
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

pour une config MH07
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

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

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

Download in other formats:

Original Format