source: LMDZ6/trunk/libf/phylmd/readaerosol_interp.f90 @ 5363

Last change on this file since 5363 was 5292, checked in by abarral, 7 weeks ago

Move academic.h chem.h chem_spla.h to module

  • 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: 23.3 KB
RevLine 
[1179]1! $Id$
2!
[1337]3SUBROUTINE readaerosol_interp(id_aero, itap, pdtphys, r_day, first, pplay, paprs, t_seri, mass_out, pi_mass_out, load_src)
[1179]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
[5292]15USE chem_mod_h
16    USE clesphys_mod_h
[1265]17  USE ioipsl
[1179]18  USE dimphy, ONLY : klev,klon
[5282]19  USE mod_phys_lmdz_para, ONLY : mpi_rank
[1179]20  USE readaerosol_mod
[1183]21  USE aero_mod, ONLY : naero_spc, name_aero
[1179]22  USE write_field_phy
[1237]23  USE phys_cal_mod
[1716]24  USE pres2lev_mod
[2311]25  USE print_control_mod, ONLY: lunout
[1179]26
[5285]27  USE yomcst_mod_h
[5274]28IMPLICIT NONE
[1179]29
[5274]30
[2344]31
[1179]32!
33! Input:
34!****************************************************************************************
35  INTEGER, INTENT(IN)                    :: id_aero! Identity number for the aerosol to treat
[1237]36  INTEGER, INTENT(IN)                    :: itap   ! Physic step count
37  REAL, INTENT(IN)                       :: pdtphys! Physic day step
[1179]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
[1216]42  REAL, DIMENSION(klon,klev), INTENT(IN) :: t_seri ! air temperature
[1179]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]
[1337]48  REAL, INTENT(OUT) :: load_src(klon) ! Load of aerosol (monthly mean data,from file) [kg/m3]
[1179]49!     
50! Local Variables:
51!****************************************************************************************
52  INTEGER                         :: i, k, ierr
[1237]53  INTEGER                         :: iday, iyr, lmt_pas
[1265]54!  INTEGER                         :: im, day1, day2, im2
55  INTEGER                         :: im, im2
56  REAL                            :: day1, day2
[1179]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
[1216]61  REAL                              :: zrho      ! Air density [kg/m3]
62  REAL                              :: volm      ! Volyme de melange [kg/kg]
[1179]63  REAL, DIMENSION(klon)             :: psurf_day, pi_psurf_day
[1337]64  REAL, DIMENSION(klon)             :: pi_load_src  ! Mass load at source grid
[1179]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)
[1265]84  INTEGER, SAVE                     :: nbr_tsteps ! number of time steps in file read
85  REAL, DIMENSION(14), SAVE         :: month_len, month_start, month_mid
[1270]86!$OMP THREADPRIVATE(nbr_tsteps, month_len, month_start, month_mid)
[1265]87  REAL                              :: jDay
[1179]88
89  LOGICAL            :: lnewday      ! Indicates if first time step at a new day
[1237]90  LOGICAL            :: OLDNEWDAY
[1179]91  LOGICAL,SAVE       :: vert_interp  ! Indicates if vertical interpolation will be done
[1181]92  LOGICAL,SAVE       :: debug=.FALSE.! Debugging in this subroutine
[1179]93!$OMP THREADPRIVATE(vert_interp, debug)
[1492]94  CHARACTER(len=8)      :: type
95  CHARACTER(len=8)      :: filename
[1179]96
97
98!****************************************************************************************
99! Initialization
100!
101!****************************************************************************************
102
103! Calculation to find if it is a new day
[1237]104
[1249]105  IF(mpi_rank == 0 .AND. debug )then
[1237]106     PRINT*,'CONTROL PANEL REGARDING TIME STEPING'
107  ENDIF
108
109  ! Use phys_cal_mod
[1273]110  iday= day_cur
111  iyr = year_cur
112  im  = mth_cur
[1237]113
[1273]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
[1265]119  CALL ymds2ju(iyr, im, iday, 0., jDay)
[1273]120!   CALL ymds2ju(iyr, im, iday-(im-1)*30, 0., jDay)
[1179]121
[1273]122
[1237]123  IF(MOD(itap-1,NINT(86400./pdtphys)) == 0)THEN
124     lnewday=.TRUE.
[1249]125  ELSE
126     lnewday=.FALSE.
[1237]127  ENDIF
[1179]128
[1249]129  IF(mpi_rank == 0 .AND. debug)then
[1237]130     ! 0.02 is about 0.5/24, namly less than half an hour
[1403]131     OLDNEWDAY = (r_day-REAL(iday) < 0.02)
[1237]132     ! Once per day, update aerosol fields
133     lmt_pas = NINT(86400./pdtphys)
[1403]134     PRINT*,'r_day-REAL(iday) =',r_day-REAL(iday)
[1237]135     PRINT*,'itap =',itap
136     PRINT*,'pdtphys =',pdtphys
137     PRINT*,'lmt_pas =',lmt_pas
138     PRINT*,'iday =',iday
139     PRINT*,'r_day =',r_day
[1246]140     PRINT*,'day_cur =',day_cur
141     PRINT*,'mth_cur =',mth_cur
142     PRINT*,'year_cur =',year_cur
[1237]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
[1179]149  IF (.NOT. ALLOCATED(var_day)) THEN
[1181]150     ALLOCATE( var_day(klon, klev, naero_spc), stat=ierr)
[2311]151     IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 1',1)
[1181]152     ALLOCATE( pi_var_day(klon, klev, naero_spc), stat=ierr)
[2311]153     IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 2',1)
[1179]154
[1270]155     ALLOCATE( psurf_year(klon, 12, naero_spc), pi_psurf_year(klon, 12, naero_spc), stat=ierr)
[2311]156     IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 3',1)
[1179]157
[1270]158     ALLOCATE( load_year(klon, 12, naero_spc), pi_load_year(klon, 12, naero_spc), stat=ierr)
[2311]159     IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 4',1)
[1270]160
[1179]161     lnewday=.TRUE.
162
163     NULLIFY(pt_ap)
164     NULLIFY(pt_b)
[2840]165  ENDIF
[1179]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!****************************************************************************************
[1246]172  IF ( (first .OR. iday==0) .AND. lnewday ) THEN
[1179]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.
[1582]177     ! If aer_type=mix1, mix2 or mix3, the run type and file name depends on the aerosol.
[1492]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'
[2840]190        ENDIF
[1492]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'
[2840]199        ENDIF
[1582]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
[1584]203           filename='aerosols'
[1582]204           type='annuel'
205        ELSE
206           filename='aerosols'
207           type='preind'
[2840]208        ENDIF
[1492]209     ELSE
[2311]210        CALL abort_physic('readaerosol_interp', 'this aer_type not supported',1)
[2840]211     ENDIF
[1492]212
213     CALL readaerosol(name_aero(id_aero), type, filename, iyr, klev_src, pt_ap, pt_b, pt_tmp, &
[1270]214          psurf_year(:,:,id_aero), load_year(:,:,id_aero))
[1179]215     IF (.NOT. ALLOCATED(var_year)) THEN
[1270]216        ALLOCATE(var_year(klon, klev_src, 12, naero_spc), stat=ierr)
[2311]217        IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 5',1)
[2840]218     ENDIF
[1179]219     var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
220
221     ! Reading values corresponding to the preindustrial concentrations.
[1492]222     type='preind'
223     CALL readaerosol(name_aero(id_aero), type, filename, iyr, pi_klev_src, pt_ap, pt_b, pt_tmp, &
[1270]224          pi_psurf_year(:,:,id_aero), pi_load_year(:,:,id_aero))
[1179]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)
[2311]231        CALL abort_physic('readaerosol_interp','Differnt vertical axes in aerosol forcing files',1)
[2840]232     ENDIF
[1179]233
234     IF (.NOT. ALLOCATED(pi_var_year)) THEN
[1270]235        ALLOCATE(pi_var_year(klon, klev_src, 12, naero_spc), stat=ierr)
[2311]236        IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 6',1)
[2840]237     ENDIF
[1179]238     pi_var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
239   
240     IF (debug) THEN
[1270]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)
[1179]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)
[2840]247     ENDIF
[1179]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'
[2311]260           CALL abort_physic('readaerosol_interp', 'The aerosol files have not the same format',1)
[2840]261        ENDIF
[1179]262       
263        IF (klev /= klev_src) THEN
264           WRITE(lunout,*) 'Old format of aerosol file do not allowed vertical interpolation'
[2311]265           CALL abort_physic('readaerosol_interp', 'Old aerosol file not possible',1)
[2840]266        ENDIF
[1179]267
268     ELSE
269        vert_interp = .TRUE.
[2840]270     ENDIF
[1179]271
[1265]272!    Calendar initialisation
273!
274     DO i = 2, 13
[1403]275       month_len(i) = REAL(ioget_mon_len(year_cur, i-1))
[1265]276       CALL ymds2ju(year_cur, i-1, 1, 0.0, month_start(i))
277     ENDDO
[1403]278     month_len(1) = REAL(ioget_mon_len(year_cur-1, 12))
[1265]279     CALL ymds2ju(year_cur-1, 12, 1, 0.0, month_start(1))
[1403]280     month_len(14) = REAL(ioget_mon_len(year_cur+1, 1))
[1265]281     CALL ymds2ju(year_cur+1, 1, 1, 0.0, month_start(14))
282     month_mid(:) = month_start (:) + month_len(:)/2.
[1273]283     
284     if (debug) then
285       write(lunout,*)' month_len = ',month_len
286       write(lunout,*)' month_mid = ',month_mid
287     endif
[1265]288
[2840]289  ENDIF  ! IF ( (first .OR. iday==0) .AND. lnewday ) THEN
[1179]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!****************************************************************************************
[1265]302    ! Find which months and days to use for time interpolation
[1270]303     nbr_tsteps = 12
[1265]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
[2840]312          ENDIF
[1265]313       ELSE
314          ! the second half of the month
315          im2=im+1
[2523]316          day1 = month_mid(im+1)
317          day2 = month_mid(im2+1)
[1265]318          IF (im2 > 12) THEN
319             ! the month is december, the following thus january
320             im2=1
321          ENDIF
[2840]322       ENDIF
[1265]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
[2523]333          day1 = month_mid(im)
334          day2 = month_mid(im2)
[2840]335       ENDIF
[1179]336     ELSE
[2311]337       CALL abort_physic('readaerosol_interp', 'number of months undefined',1)
[1265]338     ENDIF
[1273]339     if (debug) then
340       write(lunout,*)' jDay, day1, day2, im, im2 = ', jDay, day1, day2, im, im2
341     endif
[1265]342
343 
[1179]344     ! Time interpolation, still on vertical source grid
345     ALLOCATE(tmp1(klon,klev_src), tmp2(klon,klev_src),stat=ierr)
[2311]346     IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 7',1)
[1179]347
348     ALLOCATE(pplay_src(klon,klev_src), stat=ierr)
[2311]349     IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 8',1)
[1179]350     
351
352     DO k=1,klev_src
353        DO i=1,klon
354           tmp1(i,k) = &
[1269]355                var_year(i,k,im2,id_aero) - (jDay-day2)/(day1-day2) * &
[1179]356                (var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero))
357           
358           tmp2(i,k) = &
[1269]359                pi_var_year(i,k,im2,id_aero) - (jDay-day2)/(day1-day2) * &
[1179]360                (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero))
[2840]361        ENDDO
362     ENDDO
[1179]363
364     ! Time interpolation for pressure at surface, still on vertical source grid
365     DO i=1,klon
366        psurf_day(i) = &
[1269]367             psurf_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
[1179]368             (psurf_year(i,im2,id_aero) - psurf_year(i,im,id_aero))
369       
370        pi_psurf_day(i) = &
[1269]371             pi_psurf_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
[1179]372             (pi_psurf_year(i,im2,id_aero) - pi_psurf_year(i,im,id_aero))
[2840]373     ENDDO
[1179]374
375     ! Time interpolation for the load, still on vertical source grid
376     DO i=1,klon
377        load_src(i) = &
[1269]378             load_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
[1179]379             (load_year(i,im2,id_aero) - load_year(i,im,id_aero))
380       
381        pi_load_src(i) = &
[1269]382             pi_load_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
[1179]383             (pi_load_year(i,im2,id_aero) - pi_load_year(i,im,id_aero))
[2840]384     ENDDO
[1179]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)
[2840]400           ENDDO
401        ENDDO
[1179]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)
[2840]409        ENDIF
[1179]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)
[2840]425           ENDDO
426        ENDDO
[1179]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
[1216]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]
[2840]434              load_tgt(i) = load_tgt(i) + volm *delp(i,k)/RG
435           ENDDO
436        ENDDO
[1179]437       
438        ! Adjust, uniform
439        DO k = 1, klev
440           DO i = 1, klon
[2840]441              var_day(i,k,id_aero) = var_day(i,k,id_aero)*load_src(i)/max(1.e-30,load_tgt(i))
442           ENDDO
443        ENDDO
[1179]444       
445        IF (debug) THEN
446           load_tgt_test(:) = 0.
447           DO k= 1, klev
448              DO i = 1, klon
[1216]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]
[2840]451                 load_tgt_test(i) = load_tgt_test(i) + volm*delp(i,k)/RG
452              ENDDO
453           ENDDO
[1179]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)
[2840]459        ENDIF
[1179]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)
[2840]467           ENDDO
468        ENDDO
[1179]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)
[2840]473        ENDIF
[1179]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
[1216]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]
[2840]491              load_tgt(i) = load_tgt(i) + volm*delp(i,k)/RG
492           ENDDO
493        ENDDO
[1179]494
495        DO k = 1, klev
496           DO i = 1, klon
[2840]497              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))
498           ENDDO
499        ENDDO
[1179]500
501        IF (debug) THEN
502           load_tgt_test(:) = 0.
503           DO k = 1, klev
504              DO i = 1, klon
[1216]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]
[2840]507                 load_tgt_test(i) = load_tgt_test(i) + volm*delp(i,k)/RG
508              ENDDO
509           ENDDO
[1179]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)
[2840]514        ENDIF
[1179]515
516
517     ELSE   ! No vertical interpolation done
518
519        var_day(:,:,id_aero)    = tmp1(:,:)
520        pi_var_day(:,:,id_aero) = tmp2(:,:)
521
[2840]522     ENDIF ! vert_interp
[1179]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
[1265]537                 IF (jDay-day2 < 0.) WRITE(lunout,*) 'jDay-day2=',jDay-day2
[1179]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)
[2840]542                 ENDIF
[1179]543                 WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
[1273]544                 WRITE(lunout,*) 'day1, day2, jDay = ', day1, day2, jDay
[2311]545                 CALL abort_physic('readaerosol_interp','Error in interpolation 1',1)
[2840]546              ENDIF
547           ENDDO
548        ENDDO
549     ENDIF
[1179]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
[1265]556                 IF (jDay-day2 < 0.) WRITE(lunout,*) 'jDay-day2=',jDay-day2
[1179]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)
[2840]561                 ENDIF
[1179]562                 
563                 WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
[2311]564                 CALL abort_physic('readaerosol_interp','Error in interpolation 2',1)
[2840]565              ENDIF
566           ENDDO
567        ENDDO
568     ENDIF
[1179]569
[2840]570  ENDIF ! lnewday
[1179]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.