Context Navigation

readaerosol_interp.F90 @ 3872

Last change on this file since 3872 was 1583, checked in by jghattas, 13 years ago
Modified filename for specific case aer_type=mix3. Now same filename used for SO4 and the rest of aerosols.
File size: 23.3 KB

Rev	Line
[1179]	1	! $Id$
	2	!
[1337]	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]
[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)
[1483]	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
	131	OLDNEWDAY = (r_day-FLOAT(iday) < 0.02)
	132	! Once per day, update aerosol fields
	133	lmt_pas = NINT(86400./pdtphys)
	134	PRINT*,'r_day-FLOAT(iday) =',r_day-FLOAT(iday)
	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)
[1179]	151	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 1',1)
[1181]	152	ALLOCATE( pi_var_day(klon, klev, naero_spc), stat=ierr)
[1179]	153	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 2',1)
	154
[1270]	155	ALLOCATE( psurf_year(klon, 12, naero_spc), pi_psurf_year(klon, 12, naero_spc), stat=ierr)
	156	IF (ierr /= 0) CALL abort_gcm('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)
	159	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 4',1)
	160
[1179]	161	lnewday=.TRUE.
	162
	163	NULLIFY(pt_ap)
	164	NULLIFY(pt_b)
	165	END IF
	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.
[1483]	177	! If aer_type=mix1, the run type and file name depends on the aerosol.
	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'
	190	END IF
	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'
	199	END IF
[1581]	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
[1583]	203	filename='aerosols'
[1581]	204	type='annuel'
	205	ELSE
	206	filename='aerosols'
	207	type='preind'
	208	END IF
[1483]	209	ELSE
	210	CALL abort_gcm('readaerosol_interp', 'this aer_type not supported',1)
	211	END IF
	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)
	217	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 5',1)
[1179]	218	END IF
	219	var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
	220
	221	! Reading values corresponding to the preindustrial concentrations.
[1483]	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)
	231	CALL abort_gcm('readaerosol_interp','Differnt vertical axes in aerosol forcing files',1)
	232	END IF
	233
	234	IF (.NOT. ALLOCATED(pi_var_year)) THEN
[1270]	235	ALLOCATE(pi_var_year(klon, klev_src, 12, naero_spc), stat=ierr)
	236	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 6',1)
[1179]	237	END IF
	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)
	247	END IF
	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'
	260	CALL abort_gcm('readaerosol_interp', 'The aerosol files have not the same format',1)
	261	END IF
	262
	263	IF (klev /= klev_src) THEN
	264	WRITE(lunout,*) 'Old format of aerosol file do not allowed vertical interpolation'
	265	CALL abort_gcm('readaerosol_interp', 'Old aerosol file not possible',1)
	266	END IF
	267
	268	ELSE
	269	vert_interp = .TRUE.
	270	END IF
	271
[1265]	272	! Calendar initialisation
	273	!
	274	DO i = 2, 13
	275	month_len(i) = float(ioget_mon_len(year_cur, i-1))
	276	CALL ymds2ju(year_cur, i-1, 1, 0.0, month_start(i))
	277	ENDDO
	278	month_len(1) = float(ioget_mon_len(year_cur-1, 12))
	279	CALL ymds2ju(year_cur-1, 12, 1, 0.0, month_start(1))
	280	month_len(14) = float(ioget_mon_len(year_cur+1, 1))
	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
[1246]	289	END IF ! 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
	312	END IF
	313	ELSE
	314	! the second half of the month
	315	im2=im+1
	316	day2 = month_mid(im+1)
	317	day1 = month_mid(im2+1)
	318	IF (im2 > 12) THEN
	319	! the month is december, the following thus january
	320	im2=1
	321	ENDIF
	322	END IF
	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
	333	day2 = month_mid(im)
	334	day1 = month_mid(im2)
	335	END IF
[1179]	336	ELSE
[1265]	337	CALL abort_gcm('readaerosol_interp', 'number of months undefined',1)
	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)
	346	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 7',1)
	347
	348	ALLOCATE(pplay_src(klon,klev_src), stat=ierr)
	349	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 8',1)
	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))
	361	END DO
	362	END DO
	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))
	373	END DO
	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))
	384	END DO
	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)
	400	END DO
	401	END DO
	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)
	409	END IF
	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)
	425	END DO
	426	END DO
	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]
	434	load_tgt(i) = load_tgt(i) + 1/RG * volm *delp(i,k)
[1179]	435	END DO
	436	END DO
	437
	438	! Adjust, uniform
	439	DO k = 1, klev
	440	DO i = 1, klon
	441	var_day(i,k,id_aero) = var_day(i,k,id_aero)*load_src(i)/load_tgt(i)
	442	END DO
	443	END DO
	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]
	451	load_tgt_test(i) = load_tgt_test(i) + 1/RG * volm*delp(i,k)
[1179]	452	END DO
	453	END DO
	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)
	459	END IF
	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)
	467	END DO
	468	END DO
	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)
	473	END IF
	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]
	491	load_tgt(i) = load_tgt(i) + 1/RG * volm * delp(i,k)
[1179]	492	END DO
	493	END DO
	494
	495	DO k = 1, klev
	496	DO i = 1, klon
	497	pi_var_day(i,k,id_aero) = pi_var_day(i,k,id_aero)*pi_load_src(i)/load_tgt(i)
	498	END DO
	499	END DO
	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]
	507	load_tgt_test(i) = load_tgt_test(i) + 1/RG * volm * delp(i,k)
[1179]	508	END DO
	509	END DO
	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)
	514	END IF
	515
	516
	517	ELSE ! No vertical interpolation done
	518
	519	var_day(:,:,id_aero) = tmp1(:,:)
	520	pi_var_day(:,:,id_aero) = tmp2(:,:)
	521
	522	END IF ! vert_interp
	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)
	542	END IF
	543	WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
[1273]	544	WRITE(lunout,*) 'day1, day2, jDay = ', day1, day2, jDay
[1179]	545	CALL abort_gcm('readaerosol_interp','Error in interpolation 1',1)
	546	END IF
	547	END DO
	548	END DO
	549	END IF
	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)
	561	END IF
	562
	563	WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
	564	CALL abort_gcm('readaerosol_interp','Error in interpolation 2',1)
	565	END IF
	566	END DO
	567	END DO
	568	END IF
	569
	570	END IF ! lnewday
	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
	580	END SUBROUTINE readaerosol_interp

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

Download in other formats:

Original Format