source: LMDZ6/branches/Amaury_dev/libf/phylmd/readaerosol_interp.F90 @ 5157

Last change on this file since 5157 was 5144, checked in by abarral, 8 weeks ago

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