Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

readaerosol_interp.F90

Last change on this file was 1347, checked in by Laurent Fairhead, 15 years ago

Additions to aerosol outputs for CMIP5 exercise
(Needed because of chageset r1346 LF)

Additions aux sorties aérosols pour l'exercice CMIP5
(Nécessaires suite au changeset r1346 LF)

Michael, Anne

File size: 22.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format