source: LMDZ4/branches/LMDZ4_AR5/libf/phylmd/readaerosol_interp.F90 @ 4289

Last change on this file since 4289 was 1583, checked in by jghattas, 13 years ago

Modified filename for specific case aer_type=mix3. Now same filename used
for SO4 and the rest of aerosols.

File size: 23.3 KB
Line 
1! $Id$
2!
3SUBROUTINE 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-FLOAT(iday) < 0.02)
132     ! Once per day, update aerosol fields
133     lmt_pas = NINT(86400./pdtphys)
134     PRINT*,'r_day-FLOAT(iday) =',r_day-FLOAT(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, 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  IF (aer_type == 'mix3') THEN
201        ! Special case using a mix of annual sulfate file and natrual aerosols
202        IF (name_aero(id_aero) == 'SO4') THEN
203           filename='aerosols'
204           type='annuel'
205        ELSE
206           filename='aerosols'
207           type='preind'
208        END IF
209     ELSE
210        CALL abort_gcm('readaerosol_interp', 'this aer_type not supported',1)
211     END IF
212
213     CALL readaerosol(name_aero(id_aero), type, filename, iyr, klev_src, pt_ap, pt_b, pt_tmp, &
214          psurf_year(:,:,id_aero), load_year(:,:,id_aero))
215     IF (.NOT. ALLOCATED(var_year)) THEN
216        ALLOCATE(var_year(klon, klev_src, 12, naero_spc), stat=ierr)
217        IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 5',1)
218     END IF
219     var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
220
221     ! Reading values corresponding to the preindustrial concentrations.
222     type='preind'
223     CALL readaerosol(name_aero(id_aero), type, filename, iyr, pi_klev_src, pt_ap, pt_b, pt_tmp, &
224          pi_psurf_year(:,:,id_aero), pi_load_year(:,:,id_aero))
225
226     ! klev_src must be the same in both files.
227     ! Also supposing pt_ap and pt_b to be the same in the 2 files without testing.
228     IF (pi_klev_src /= klev_src) THEN
229        WRITE(lunout,*) 'Error! All forcing files for the same aerosol must have the same vertical dimension'
230        WRITE(lunout,*) 'Aerosol : ', name_aero(id_aero)
231        CALL abort_gcm('readaerosol_interp','Differnt vertical axes in aerosol forcing files',1)
232     END IF
233
234     IF (.NOT. ALLOCATED(pi_var_year)) THEN
235        ALLOCATE(pi_var_year(klon, klev_src, 12, naero_spc), stat=ierr)
236        IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 6',1)
237     END IF
238     pi_var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
239   
240     IF (debug) THEN
241        CALL writefield_phy('var_year_jan',var_year(:,:,1,id_aero),klev_src)
242        CALL writefield_phy('var_year_dec',var_year(:,:,12,id_aero),klev_src)
243        CALL writefield_phy('psurf_src',psurf_year(:,:,id_aero),1)
244        CALL writefield_phy('pi_psurf_src',pi_psurf_year(:,:,id_aero),1)
245        CALL writefield_phy('load_year_src',load_year(:,:,id_aero),1)
246        CALL writefield_phy('pi_load_year_src',pi_load_year(:,:,id_aero),1)
247     END IF
248
249     ! Pointer no more useful, deallocate.
250     DEALLOCATE(pt_tmp)
251
252     ! Test if vertical interpolation will be needed.
253     IF (psurf_year(1,1,id_aero)==not_valid .OR. pi_psurf_year(1,1,id_aero)==not_valid ) THEN
254        ! Pressure=not_valid indicates old file format, see module readaerosol
255        vert_interp = .FALSE.
256
257        ! If old file format, both psurf_year and pi_psurf_year must be not_valid
258        IF (  psurf_year(1,1,id_aero) /= pi_psurf_year(1,1,id_aero) ) THEN
259           WRITE(lunout,*) 'Warning! All forcing files for the same aerosol must have the same structure'
260           CALL abort_gcm('readaerosol_interp', 'The aerosol files have not the same format',1)
261        END IF
262       
263        IF (klev /= klev_src) THEN
264           WRITE(lunout,*) 'Old format of aerosol file do not allowed vertical interpolation'
265           CALL abort_gcm('readaerosol_interp', 'Old aerosol file not possible',1)
266        END IF
267
268     ELSE
269        vert_interp = .TRUE.
270     END IF
271
272!    Calendar initialisation
273!
274     DO i = 2, 13
275       month_len(i) = float(ioget_mon_len(year_cur, i-1))
276       CALL ymds2ju(year_cur, i-1, 1, 0.0, month_start(i))
277     ENDDO
278     month_len(1) = float(ioget_mon_len(year_cur-1, 12))
279     CALL ymds2ju(year_cur-1, 12, 1, 0.0, month_start(1))
280     month_len(14) = float(ioget_mon_len(year_cur+1, 1))
281     CALL ymds2ju(year_cur+1, 1, 1, 0.0, month_start(14))
282     month_mid(:) = month_start (:) + month_len(:)/2.
283     
284     if (debug) then
285       write(lunout,*)' month_len = ',month_len
286       write(lunout,*)' month_mid = ',month_mid
287     endif
288
289  END IF  ! IF ( (first .OR. iday==0) .AND. lnewday ) THEN
290 
291!****************************************************************************************
292! - 2) Interpolate to the actual day.
293! - 3) Interpolate to the model vertical grid.
294!
295!****************************************************************************************
296
297  IF (lnewday) THEN ! only if new day
298!****************************************************************************************
299! 2) Interpolate to the actual day
300!
301!****************************************************************************************
302    ! Find which months and days to use for time interpolation
303     nbr_tsteps = 12
304     IF (nbr_tsteps == 12) then
305       IF (jDay < month_mid(im+1)) THEN
306          im2=im-1
307          day2 = month_mid(im2+1)
308          day1 = month_mid(im+1)
309          IF (im2 <= 0) THEN
310             ! the month is january, thus the month before december
311             im2=12
312          END IF
313       ELSE
314          ! the second half of the month
315          im2=im+1
316          day2 = month_mid(im+1)
317          day1 = month_mid(im2+1)
318          IF (im2 > 12) THEN
319             ! the month is december, the following thus january
320             im2=1
321          ENDIF
322       END IF
323     ELSE IF (nbr_tsteps == 14) then
324       im = im + 1
325       IF (jDay < month_mid(im)) THEN
326          ! in the first half of the month use month before and actual month
327          im2=im-1
328          day2 = month_mid(im2)
329          day1 = month_mid(im)
330       ELSE
331          ! the second half of the month
332          im2=im+1
333          day2 = month_mid(im)
334          day1 = month_mid(im2)
335       END IF
336     ELSE
337       CALL abort_gcm('readaerosol_interp', 'number of months undefined',1)
338     ENDIF
339     if (debug) then
340       write(lunout,*)' jDay, day1, day2, im, im2 = ', jDay, day1, day2, im, im2
341     endif
342
343 
344     ! Time interpolation, still on vertical source grid
345     ALLOCATE(tmp1(klon,klev_src), tmp2(klon,klev_src),stat=ierr)
346     IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 7',1)
347
348     ALLOCATE(pplay_src(klon,klev_src), stat=ierr)
349     IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 8',1)
350     
351
352     DO k=1,klev_src
353        DO i=1,klon
354           tmp1(i,k) = &
355                var_year(i,k,im2,id_aero) - (jDay-day2)/(day1-day2) * &
356                (var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero))
357           
358           tmp2(i,k) = &
359                pi_var_year(i,k,im2,id_aero) - (jDay-day2)/(day1-day2) * &
360                (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero))
361        END DO
362     END DO
363
364     ! Time interpolation for pressure at surface, still on vertical source grid
365     DO i=1,klon
366        psurf_day(i) = &
367             psurf_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
368             (psurf_year(i,im2,id_aero) - psurf_year(i,im,id_aero))
369       
370        pi_psurf_day(i) = &
371             pi_psurf_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
372             (pi_psurf_year(i,im2,id_aero) - pi_psurf_year(i,im,id_aero))
373     END DO
374
375     ! Time interpolation for the load, still on vertical source grid
376     DO i=1,klon
377        load_src(i) = &
378             load_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
379             (load_year(i,im2,id_aero) - load_year(i,im,id_aero))
380       
381        pi_load_src(i) = &
382             pi_load_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
383             (pi_load_year(i,im2,id_aero) - pi_load_year(i,im,id_aero))
384     END DO
385
386!****************************************************************************************
387! 3) Interpolate to the model vertical grid (target grid)
388!
389!****************************************************************************************
390
391     IF (vert_interp) THEN
392
393        ! - Interpolate variable tmp1 (on source grid) to var_day (on target grid)
394        !********************************************************************************
395        ! a) calculate pression at vertical levels for the source grid using the
396        !    hybrid-sigma coordinates ap and b and the surface pressure, variables from file.
397        DO k = 1, klev_src
398           DO i = 1, klon
399              pplay_src(i,k)= pt_ap(k) + pt_b(k)*psurf_day(i)
400           END DO
401        END DO
402       
403        IF (debug) THEN
404           CALL writefield_phy('psurf_day_src',psurf_day(:),1)
405           CALL writefield_phy('pplay_src',pplay_src(:,:),klev_src)
406           CALL writefield_phy('pplay',pplay(:,:),klev)
407           CALL writefield_phy('day_src',tmp1,klev_src)
408           CALL writefield_phy('pi_day_src',tmp2,klev_src)
409        END IF
410
411        ! b) vertical interpolation on pressure leveles
412        CALL pres2lev(tmp1(:,:), var_day(:,:,id_aero), klev_src, klev, pplay_src, pplay, &
413             1, klon, .FALSE.)
414       
415        IF (debug) CALL writefield_phy('day_tgt',var_day(:,:,id_aero),klev)
416       
417        ! c) adjust to conserve total aerosol mass load in the vertical pillar
418        !    Calculate the load in the actual pillar and compare with the load
419        !    read from aerosol file.
420       
421        ! Find the pressure difference in each model layer
422        DO k = 1, klev
423           DO i = 1, klon
424              delp(i,k) = paprs(i,k) - paprs (i,k+1)
425           END DO
426        END DO
427
428        ! Find the mass load in the actual pillar, on target grid
429        load_tgt(:) = 0.
430        DO k= 1, klev
431           DO i = 1, klon
432              zrho = pplay(i,k)/t_seri(i,k)/RD       ! [kg/m3]
433              volm = var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
434              load_tgt(i) = load_tgt(i) + 1/RG * volm *delp(i,k)
435           END DO
436        END DO
437       
438        ! Adjust, uniform
439        DO k = 1, klev
440           DO i = 1, klon
441              var_day(i,k,id_aero) = var_day(i,k,id_aero)*load_src(i)/load_tgt(i)
442           END DO
443        END DO
444       
445        IF (debug) THEN
446           load_tgt_test(:) = 0.
447           DO k= 1, klev
448              DO i = 1, klon
449                 zrho = pplay(i,k)/t_seri(i,k)/RD       ! [kg/m3]
450                 volm = var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
451                 load_tgt_test(i) = load_tgt_test(i) + 1/RG * volm*delp(i,k)
452              END DO
453           END DO
454           
455           CALL writefield_phy('day_tgt2',var_day(:,:,id_aero),klev)
456           CALL writefield_phy('load_tgt',load_tgt(:),1)
457           CALL writefield_phy('load_tgt_test',load_tgt_test(:),1)
458           CALL writefield_phy('load_src',load_src(:),1)
459        END IF
460
461        ! - Interpolate variable tmp2 (source grid) to pi_var_day (target grid)
462        !********************************************************************************
463        ! a) calculate pression at vertical levels at source grid   
464        DO k = 1, klev_src
465           DO i = 1, klon
466              pplay_src(i,k)= pt_ap(k) + pt_b(k)*pi_psurf_day(i)
467           END DO
468        END DO
469
470        IF (debug) THEN
471           CALL writefield_phy('pi_psurf_day_src',pi_psurf_day(:),1)
472           CALL writefield_phy('pi_pplay_src',pplay_src(:,:),klev_src)
473        END IF
474
475        ! b) vertical interpolation on pressure leveles
476        CALL pres2lev(tmp2(:,:), pi_var_day(:,:,id_aero), klev_src, klev, pplay_src, pplay, &
477             1, klon, .FALSE.)
478
479        IF (debug) CALL writefield_phy('pi_day_tgt',pi_var_day(:,:,id_aero),klev)
480
481        ! c) adjust to conserve total aerosol mass load in the vertical pillar
482        !    Calculate the load in the actual pillar and compare with the load
483        !    read from aerosol file.
484
485        ! Find the load in the actual pillar, on target grid
486        load_tgt(:) = 0.
487        DO k = 1, klev
488           DO i = 1, klon
489              zrho = pplay(i,k)/t_seri(i,k)/RD          ! [kg/m3]
490              volm = pi_var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
491              load_tgt(i) = load_tgt(i) + 1/RG * volm * delp(i,k)
492           END DO
493        END DO
494
495        DO k = 1, klev
496           DO i = 1, klon
497              pi_var_day(i,k,id_aero) = pi_var_day(i,k,id_aero)*pi_load_src(i)/load_tgt(i)
498           END DO
499        END DO
500
501        IF (debug) THEN
502           load_tgt_test(:) = 0.
503           DO k = 1, klev
504              DO i = 1, klon
505                 zrho = pplay(i,k)/t_seri(i,k)/RD          ! [kg/m3]
506                 volm = pi_var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
507                 load_tgt_test(i) = load_tgt_test(i) + 1/RG * volm * delp(i,k)
508              END DO
509           END DO
510           CALL writefield_phy('pi_day_tgt2',pi_var_day(:,:,id_aero),klev)
511           CALL writefield_phy('pi_load_tgt',load_tgt(:),1)
512           CALL writefield_phy('pi_load_tgt_test',load_tgt_test(:),1)
513           CALL writefield_phy('pi_load_src',pi_load_src(:),1)
514        END IF
515
516
517     ELSE   ! No vertical interpolation done
518
519        var_day(:,:,id_aero)    = tmp1(:,:)
520        pi_var_day(:,:,id_aero) = tmp2(:,:)
521
522     END IF ! vert_interp
523
524
525     ! Deallocation
526     DEALLOCATE(tmp1, tmp2, pplay_src, stat=ierr)
527
528!****************************************************************************************
529! 4) Test for negative mass values
530!
531!****************************************************************************************
532     IF (MINVAL(var_day(:,:,id_aero)) < 0.) THEN
533        DO k=1,klev
534           DO i=1,klon
535              ! Test for var_day
536              IF (var_day(i,k,id_aero) < 0.) THEN
537                 IF (jDay-day2 < 0.) WRITE(lunout,*) 'jDay-day2=',jDay-day2
538                 IF (var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero) < 0.) THEN
539                    WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', &
540                         trim(name_aero(id_aero)),'(i,k,im)=',           &
541                         var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero)
542                 END IF
543                 WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
544                 WRITE(lunout,*) 'day1, day2, jDay = ', day1, day2, jDay
545                 CALL abort_gcm('readaerosol_interp','Error in interpolation 1',1)
546              END IF
547           END DO
548        END DO
549     END IF
550
551     IF (MINVAL(pi_var_day(:,:,id_aero)) < 0. ) THEN
552        DO k=1, klev
553           DO i=1,klon
554              ! Test for pi_var_day
555              IF (pi_var_day(i,k,id_aero) < 0.) THEN
556                 IF (jDay-day2 < 0.) WRITE(lunout,*) 'jDay-day2=',jDay-day2
557                 IF (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero) < 0.) THEN
558                    WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', &
559                         trim(name_aero(id_aero)),'(i,k,im)=',           &
560                         pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero)
561                 END IF
562                 
563                 WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
564                 CALL abort_gcm('readaerosol_interp','Error in interpolation 2',1)
565              END IF
566           END DO
567        END DO
568     END IF
569
570  END IF ! lnewday
571
572!****************************************************************************************
573! Copy output from saved variables
574!
575!****************************************************************************************
576
577  mass_out(:,:)    = var_day(:,:,id_aero)
578  pi_mass_out(:,:) = pi_var_day(:,:,id_aero)
579 
580END SUBROUTINE readaerosol_interp
Note: See TracBrowser for help on using the repository browser.