1 |
|
|
! $Id$ |
2 |
|
|
! |
3 |
|
|
SUBROUTINE 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 |
|
|
ENDIF |
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 |
|
|
ENDIF |
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 |
|
|
ENDIF |
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 |
|
|
ENDIF |
208 |
|
|
ELSE |
209 |
|
|
CALL abort_physic('readaerosol_interp', 'this aer_type not supported',1) |
210 |
|
|
ENDIF |
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 |
|
|
ENDIF |
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 |
|
|
ENDIF |
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 |
|
|
ENDIF |
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 |
|
|
ENDIF |
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 |
|
|
ENDIF |
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 |
|
|
ENDIF |
266 |
|
|
|
267 |
|
|
ELSE |
268 |
|
|
vert_interp = .TRUE. |
269 |
|
|
ENDIF |
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 |
|
|
ENDIF ! 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 |
|
|
ENDIF |
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 |
|
|
ENDIF |
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 |
|
|
ENDIF |
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 |
|
|
ENDDO |
361 |
|
|
ENDDO |
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 |
|
|
ENDDO |
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 |
|
|
ENDDO |
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 |
|
|
ENDDO |
400 |
|
|
ENDDO |
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 |
|
|
ENDIF |
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 |
|
|
ENDDO |
425 |
|
|
ENDDO |
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) + volm *delp(i,k)/RG |
434 |
|
|
ENDDO |
435 |
|
|
ENDDO |
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)/max(1.e-30,load_tgt(i)) |
441 |
|
|
ENDDO |
442 |
|
|
ENDDO |
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) + volm*delp(i,k)/RG |
451 |
|
|
ENDDO |
452 |
|
|
ENDDO |
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 |
|
|
ENDIF |
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 |
|
|
ENDDO |
467 |
|
|
ENDDO |
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 |
|
|
ENDIF |
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) + volm*delp(i,k)/RG |
491 |
|
|
ENDDO |
492 |
|
|
ENDDO |
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)/max(1.e-30,load_tgt(i)) |
497 |
|
|
ENDDO |
498 |
|
|
ENDDO |
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) + volm*delp(i,k)/RG |
507 |
|
|
ENDDO |
508 |
|
|
ENDDO |
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 |
|
|
ENDIF |
514 |
|
|
|
515 |
|
|
|
516 |
|
|
ELSE ! No vertical interpolation done |
517 |
|
|
|
518 |
|
|
var_day(:,:,id_aero) = tmp1(:,:) |
519 |
|
|
pi_var_day(:,:,id_aero) = tmp2(:,:) |
520 |
|
|
|
521 |
|
|
ENDIF ! 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 |
|
|
ENDIF |
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 |
|
|
ENDIF |
546 |
|
|
ENDDO |
547 |
|
|
ENDDO |
548 |
|
|
ENDIF |
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 |
|
|
ENDIF |
561 |
|
|
|
562 |
|
|
WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero) |
563 |
|
|
CALL abort_physic('readaerosol_interp','Error in interpolation 2',1) |
564 |
|
|
ENDIF |
565 |
|
|
ENDDO |
566 |
|
|
ENDDO |
567 |
|
|
ENDIF |
568 |
|
|
|
569 |
|
|
ENDIF ! 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 |
|
|
|
579 |
|
|
END SUBROUTINE readaerosol_interp |