Context Navigation

readaerosol_interp.F90 @ 357

Last change on this file since 357 was 66, checked in by emillour, 14 years ago

EM: Mise a niveau par rapport a la version terrestre (LMDZ5V2.0-dev, rev 1487)

Mise a jour des scripts (terrestres) 'makegcm' et 'create_make_gcm'
Ajout du script 'makelmdz' (version "amelioree", en Bash, de makegcm)
Mise a jour des routines dans phylmd (sauf regr_lat_time_climoz_m.F)

disvert (dans dyn3d et dyn3dpar): passage au Fortran 90
parallel.F90 (dyn3dpar): correction bug
etat0_netcdf.F90 (dyn3d et dyn3dpar) : mise a jour mineure
ce0l.F90 (dyn3dpar) : correction bug
abort_gcm.F (dyn3dpar) : correction bug
ugeostr.F90 (dyn3d et dyn3dpar) : passage au Fortran 90
fluxstokenc_p.F (dyn3dpar) : correction bug
iniacademic.F90 (dyn3d et dyn3dpar) : passage au Fortran 90
friction_p.F (dyn3dpar) : correction bug
infotrac.F90 (dyn3d et dyn3dpar) : correction bug mineur sur lecture traceurs
caladvtrac.F (dyn3d) : modifications cosmetiques

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

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

Download in other formats:

Original Format