source: LMDZ4/branches/LMDZ4-dev-20091210/libf/phylmd/readaerosol_interp.F90 @ 3241

Last change on this file since 3241 was 1273, checked in by Laurent Fairhead, 15 years ago

Correction d'un bug sur le calcul de la date courante pour la lecture
des aerosols
Attention, les resultats sont modifies suite a un decalage d'une journee
dans l'ancien calcul de la date
LF

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