source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/readaerosol_interp.F90

Last change on this file was 3990, checked in by millour, 8 years ago

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