Context Navigation

readaerosol_interp.f90 @ 5282

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

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

Download in other formats:

Original Format