Context Navigation

readaerosol_interp.F90 @ 2176

Last change on this file since 2176 was 1910, checked in by Laurent Fairhead, 11 years ago
Merged trunk changes r1860:1909 into testing branch
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.4 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
[1750]	22	USE pres2lev_mod
[1179]	23
	24	IMPLICIT NONE
	25
	26	INCLUDE "YOMCST.h"
	27	INCLUDE "chem.h"
	28	INCLUDE "temps.h"
	29	INCLUDE "clesphys.h"
	30	INCLUDE "iniprint.h"
	31	INCLUDE "dimensions.h"
	32	INCLUDE "comvert.h"
	33	!
	34	! Input:
	35	!****************************************************************************************
	36	INTEGER, INTENT(IN) :: id_aero! Identity number for the aerosol to treat
[1237]	37	INTEGER, INTENT(IN) :: itap ! Physic step count
	38	REAL, INTENT(IN) :: pdtphys! Physic day step
[1179]	39	REAL, INTENT(IN) :: r_day ! Day of integration
	40	LOGICAL, INTENT(IN) :: first ! First model timestep
	41	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! pression at model mid-layers
	42	REAL, DIMENSION(klon,klev+1),INTENT(IN):: paprs ! pression between model layers
[1216]	43	REAL, DIMENSION(klon,klev), INTENT(IN) :: t_seri ! air temperature
[1179]	44	!
	45	! Output:
	46	!****************************************************************************************
	47	REAL, INTENT(OUT) :: mass_out(klon,klev) ! Mass of aerosol (monthly mean data,from file) [ug AIBCM/m3]
	48	REAL, INTENT(OUT) :: pi_mass_out(klon,klev) ! Mass of preindustrial aerosol (monthly mean data,from file) [ug AIBCM/m3]
[1337]	49	REAL, INTENT(OUT) :: load_src(klon) ! Load of aerosol (monthly mean data,from file) [kg/m3]
[1179]	50	!
	51	! Local Variables:
	52	!****************************************************************************************
	53	INTEGER :: i, k, ierr
[1237]	54	INTEGER :: iday, iyr, lmt_pas
[1265]	55	! INTEGER :: im, day1, day2, im2
	56	INTEGER :: im, im2
	57	REAL :: day1, day2
[1179]	58	INTEGER :: pi_klev_src ! Only for testing purpose
	59	INTEGER, SAVE :: klev_src ! Number of vertical levles in source field
	60	!$OMP THREADPRIVATE(klev_src)
	61
[1216]	62	REAL :: zrho ! Air density [kg/m3]
	63	REAL :: volm ! Volyme de melange [kg/kg]
[1179]	64	REAL, DIMENSION(klon) :: psurf_day, pi_psurf_day
[1337]	65	REAL, DIMENSION(klon) :: pi_load_src ! Mass load at source grid
[1179]	66	REAL, DIMENSION(klon) :: load_tgt, load_tgt_test
	67	REAL, DIMENSION(klon,klev) :: delp ! pressure difference in each model layer
	68
	69	REAL, ALLOCATABLE, DIMENSION(:,:) :: pplay_src ! pression mid-layer at source levels
	70	REAL, ALLOCATABLE, DIMENSION(:,:) :: tmp1, tmp2 ! Temporary variables
	71	REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: var_year ! VAR in right dimension for the total year
	72	REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: pi_var_year ! pre-industrial VAR, -"-
	73	!$OMP THREADPRIVATE(var_year,pi_var_year)
	74	REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: var_day ! VAR interpolated to the actual day and model grid
	75	REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: pi_var_day ! pre-industrial VAR, -"-
	76	!$OMP THREADPRIVATE(var_day,pi_var_day)
	77	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: psurf_year, pi_psurf_year ! surface pressure for the total year
	78	!$OMP THREADPRIVATE(psurf_year, pi_psurf_year)
	79	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: load_year, pi_load_year ! load in the column for the total year
	80	!$OMP THREADPRIVATE(load_year, pi_load_year)
	81
	82	REAL, DIMENSION(:,:,:), POINTER :: pt_tmp ! Pointer allocated in readaerosol
	83	REAL, POINTER, DIMENSION(:), SAVE :: pt_ap, pt_b ! Pointer for describing the vertical levels
	84	!$OMP THREADPRIVATE(pt_ap, pt_b)
[1265]	85	INTEGER, SAVE :: nbr_tsteps ! number of time steps in file read
	86	REAL, DIMENSION(14), SAVE :: month_len, month_start, month_mid
[1270]	87	!$OMP THREADPRIVATE(nbr_tsteps, month_len, month_start, month_mid)
[1265]	88	REAL :: jDay
[1179]	89
	90	LOGICAL :: lnewday ! Indicates if first time step at a new day
[1237]	91	LOGICAL :: OLDNEWDAY
[1179]	92	LOGICAL,SAVE :: vert_interp ! Indicates if vertical interpolation will be done
[1181]	93	LOGICAL,SAVE :: debug=.FALSE.! Debugging in this subroutine
[1179]	94	!$OMP THREADPRIVATE(vert_interp, debug)
[1492]	95	CHARACTER(len=8) :: type
	96	CHARACTER(len=8) :: filename
[1179]	97
	98
	99	!****************************************************************************************
	100	! Initialization
	101	!
	102	!****************************************************************************************
	103
	104	! Calculation to find if it is a new day
[1237]	105
[1249]	106	IF(mpi_rank == 0 .AND. debug )then
[1237]	107	PRINT*,'CONTROL PANEL REGARDING TIME STEPING'
	108	ENDIF
	109
	110	! Use phys_cal_mod
[1273]	111	iday= day_cur
	112	iyr = year_cur
	113	im = mth_cur
[1237]	114
[1273]	115	! iday = INT(r_day)
	116	! iyr = iday/360
	117	! iday = iday-iyr*360 ! day of the actual year
	118	! iyr = iyr + annee_ref ! year of the run
	119	! im = iday/30 +1 ! the actual month
[1265]	120	CALL ymds2ju(iyr, im, iday, 0., jDay)
[1273]	121	! CALL ymds2ju(iyr, im, iday-(im-1)*30, 0., jDay)
[1179]	122
[1273]	123
[1237]	124	IF(MOD(itap-1,NINT(86400./pdtphys)) == 0)THEN
	125	lnewday=.TRUE.
[1249]	126	ELSE
	127	lnewday=.FALSE.
[1237]	128	ENDIF
[1179]	129
[1249]	130	IF(mpi_rank == 0 .AND. debug)then
[1237]	131	! 0.02 is about 0.5/24, namly less than half an hour
[1403]	132	OLDNEWDAY = (r_day-REAL(iday) < 0.02)
[1237]	133	! Once per day, update aerosol fields
	134	lmt_pas = NINT(86400./pdtphys)
[1403]	135	PRINT*,'r_day-REAL(iday) =',r_day-REAL(iday)
[1237]	136	PRINT*,'itap =',itap
	137	PRINT*,'pdtphys =',pdtphys
	138	PRINT*,'lmt_pas =',lmt_pas
	139	PRINT*,'iday =',iday
	140	PRINT*,'r_day =',r_day
[1246]	141	PRINT*,'day_cur =',day_cur
	142	PRINT*,'mth_cur =',mth_cur
	143	PRINT*,'year_cur =',year_cur
[1237]	144	PRINT*,'NINT(86400./pdtphys) =',NINT(86400./pdtphys)
	145	PRINT*,'MOD(0,1) =',MOD(0,1)
	146	PRINT*,'lnewday =',lnewday
	147	PRINT*,'OLDNEWDAY =',OLDNEWDAY
	148	ENDIF
	149
[1179]	150	IF (.NOT. ALLOCATED(var_day)) THEN
[1181]	151	ALLOCATE( var_day(klon, klev, naero_spc), stat=ierr)
[1179]	152	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 1',1)
[1181]	153	ALLOCATE( pi_var_day(klon, klev, naero_spc), stat=ierr)
[1179]	154	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 2',1)
	155
[1270]	156	ALLOCATE( psurf_year(klon, 12, naero_spc), pi_psurf_year(klon, 12, naero_spc), stat=ierr)
	157	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 3',1)
[1179]	158
[1270]	159	ALLOCATE( load_year(klon, 12, naero_spc), pi_load_year(klon, 12, naero_spc), stat=ierr)
	160	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 4',1)
	161
[1179]	162	lnewday=.TRUE.
	163
	164	NULLIFY(pt_ap)
	165	NULLIFY(pt_b)
	166	END IF
	167
	168	!****************************************************************************************
	169	! 1) Read in data : corresponding to the actual year and preindustrial data.
	170	! Only for the first day of the year.
	171	!
	172	!****************************************************************************************
[1246]	173	IF ( (first .OR. iday==0) .AND. lnewday ) THEN
[1179]	174	NULLIFY(pt_tmp)
	175
	176	! Reading values corresponding to the closest year taking into count the choice of aer_type.
	177	! For aer_type=scenario interpolation between 2 data sets is done in readaerosol.
[1665]	178	! If aer_type=mix1, mix2 or mix3, the run type and file name depends on the aerosol.
[1492]	179	IF (aer_type=='preind' .OR. aer_type=='actuel' .OR. aer_type=='annuel' .OR. aer_type=='scenario') THEN
	180	! Standard case
	181	filename='aerosols'
	182	type=aer_type
	183	ELSE IF (aer_type == 'mix1') THEN
	184	! Special case using a mix of decenal sulfate file and annual aerosols(all aerosols except sulfate)
	185	IF (name_aero(id_aero) == 'SO4') THEN
	186	filename='so4.run '
	187	type='scenario'
	188	ELSE
	189	filename='aerosols'
	190	type='annuel'
	191	END IF
	192	ELSE IF (aer_type == 'mix2') THEN
	193	! Special case using a mix of decenal sulfate file and natrual aerosols
	194	IF (name_aero(id_aero) == 'SO4') THEN
	195	filename='so4.run '
	196	type='scenario'
	197	ELSE
	198	filename='aerosols'
	199	type='preind'
	200	END IF
[1665]	201	ELSE IF (aer_type == 'mix3') THEN
	202	! Special case using a mix of annual sulfate file and natrual aerosols
	203	IF (name_aero(id_aero) == 'SO4') THEN
	204	filename='aerosols'
	205	type='annuel'
	206	ELSE
	207	filename='aerosols'
	208	type='preind'
	209	END IF
[1492]	210	ELSE
	211	CALL abort_gcm('readaerosol_interp', 'this aer_type not supported',1)
	212	END IF
	213
	214	CALL readaerosol(name_aero(id_aero), type, filename, iyr, klev_src, pt_ap, pt_b, pt_tmp, &
[1270]	215	psurf_year(:,:,id_aero), load_year(:,:,id_aero))
[1179]	216	IF (.NOT. ALLOCATED(var_year)) THEN
[1270]	217	ALLOCATE(var_year(klon, klev_src, 12, naero_spc), stat=ierr)
	218	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 5',1)
[1179]	219	END IF
	220	var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
	221
	222	! Reading values corresponding to the preindustrial concentrations.
[1492]	223	type='preind'
	224	CALL readaerosol(name_aero(id_aero), type, filename, iyr, pi_klev_src, pt_ap, pt_b, pt_tmp, &
[1270]	225	pi_psurf_year(:,:,id_aero), pi_load_year(:,:,id_aero))
[1179]	226
	227	! klev_src must be the same in both files.
	228	! Also supposing pt_ap and pt_b to be the same in the 2 files without testing.
	229	IF (pi_klev_src /= klev_src) THEN
	230	WRITE(lunout,*) 'Error! All forcing files for the same aerosol must have the same vertical dimension'
	231	WRITE(lunout,*) 'Aerosol : ', name_aero(id_aero)
	232	CALL abort_gcm('readaerosol_interp','Differnt vertical axes in aerosol forcing files',1)
	233	END IF
	234
	235	IF (.NOT. ALLOCATED(pi_var_year)) THEN
[1270]	236	ALLOCATE(pi_var_year(klon, klev_src, 12, naero_spc), stat=ierr)
	237	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 6',1)
[1179]	238	END IF
	239	pi_var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
	240
	241	IF (debug) THEN
[1270]	242	CALL writefield_phy('var_year_jan',var_year(:,:,1,id_aero),klev_src)
	243	CALL writefield_phy('var_year_dec',var_year(:,:,12,id_aero),klev_src)
[1179]	244	CALL writefield_phy('psurf_src',psurf_year(:,:,id_aero),1)
	245	CALL writefield_phy('pi_psurf_src',pi_psurf_year(:,:,id_aero),1)
	246	CALL writefield_phy('load_year_src',load_year(:,:,id_aero),1)
	247	CALL writefield_phy('pi_load_year_src',pi_load_year(:,:,id_aero),1)
	248	END IF
	249
	250	! Pointer no more useful, deallocate.
	251	DEALLOCATE(pt_tmp)
	252
	253	! Test if vertical interpolation will be needed.
	254	IF (psurf_year(1,1,id_aero)==not_valid .OR. pi_psurf_year(1,1,id_aero)==not_valid ) THEN
	255	! Pressure=not_valid indicates old file format, see module readaerosol
	256	vert_interp = .FALSE.
	257
	258	! If old file format, both psurf_year and pi_psurf_year must be not_valid
	259	IF ( psurf_year(1,1,id_aero) /= pi_psurf_year(1,1,id_aero) ) THEN
	260	WRITE(lunout,*) 'Warning! All forcing files for the same aerosol must have the same structure'
	261	CALL abort_gcm('readaerosol_interp', 'The aerosol files have not the same format',1)
	262	END IF
	263
	264	IF (klev /= klev_src) THEN
	265	WRITE(lunout,*) 'Old format of aerosol file do not allowed vertical interpolation'
	266	CALL abort_gcm('readaerosol_interp', 'Old aerosol file not possible',1)
	267	END IF
	268
	269	ELSE
	270	vert_interp = .TRUE.
	271	END IF
	272
[1265]	273	! Calendar initialisation
	274	!
	275	DO i = 2, 13
[1403]	276	month_len(i) = REAL(ioget_mon_len(year_cur, i-1))
[1265]	277	CALL ymds2ju(year_cur, i-1, 1, 0.0, month_start(i))
	278	ENDDO
[1403]	279	month_len(1) = REAL(ioget_mon_len(year_cur-1, 12))
[1265]	280	CALL ymds2ju(year_cur-1, 12, 1, 0.0, month_start(1))
[1403]	281	month_len(14) = REAL(ioget_mon_len(year_cur+1, 1))
[1265]	282	CALL ymds2ju(year_cur+1, 1, 1, 0.0, month_start(14))
	283	month_mid(:) = month_start (:) + month_len(:)/2.
[1273]	284
	285	if (debug) then
	286	write(lunout,*)' month_len = ',month_len
	287	write(lunout,*)' month_mid = ',month_mid
	288	endif
[1265]	289
[1246]	290	END IF ! IF ( (first .OR. iday==0) .AND. lnewday ) THEN
[1179]	291
	292	!****************************************************************************************
	293	! - 2) Interpolate to the actual day.
	294	! - 3) Interpolate to the model vertical grid.
	295	!
	296	!****************************************************************************************
	297
	298	IF (lnewday) THEN ! only if new day
	299	!****************************************************************************************
	300	! 2) Interpolate to the actual day
	301	!
	302	!****************************************************************************************
[1265]	303	! Find which months and days to use for time interpolation
[1270]	304	nbr_tsteps = 12
[1265]	305	IF (nbr_tsteps == 12) then
	306	IF (jDay < month_mid(im+1)) THEN
	307	im2=im-1
	308	day2 = month_mid(im2+1)
	309	day1 = month_mid(im+1)
	310	IF (im2 <= 0) THEN
	311	! the month is january, thus the month before december
	312	im2=12
	313	END IF
	314	ELSE
	315	! the second half of the month
	316	im2=im+1
	317	day2 = month_mid(im+1)
	318	day1 = month_mid(im2+1)
	319	IF (im2 > 12) THEN
	320	! the month is december, the following thus january
	321	im2=1
	322	ENDIF
	323	END IF
	324	ELSE IF (nbr_tsteps == 14) then
	325	im = im + 1
	326	IF (jDay < month_mid(im)) THEN
	327	! in the first half of the month use month before and actual month
	328	im2=im-1
	329	day2 = month_mid(im2)
	330	day1 = month_mid(im)
	331	ELSE
	332	! the second half of the month
	333	im2=im+1
	334	day2 = month_mid(im)
	335	day1 = month_mid(im2)
	336	END IF
[1179]	337	ELSE
[1265]	338	CALL abort_gcm('readaerosol_interp', 'number of months undefined',1)
	339	ENDIF
[1273]	340	if (debug) then
	341	write(lunout,*)' jDay, day1, day2, im, im2 = ', jDay, day1, day2, im, im2
	342	endif
[1265]	343
	344
[1179]	345	! Time interpolation, still on vertical source grid
	346	ALLOCATE(tmp1(klon,klev_src), tmp2(klon,klev_src),stat=ierr)
	347	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 7',1)
	348
	349	ALLOCATE(pplay_src(klon,klev_src), stat=ierr)
	350	IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 8',1)
	351
	352
	353	DO k=1,klev_src
	354	DO i=1,klon
	355	tmp1(i,k) = &
[1269]	356	var_year(i,k,im2,id_aero) - (jDay-day2)/(day1-day2) * &
[1179]	357	(var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero))
	358
	359	tmp2(i,k) = &
[1269]	360	pi_var_year(i,k,im2,id_aero) - (jDay-day2)/(day1-day2) * &
[1179]	361	(pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero))
	362	END DO
	363	END DO
	364
	365	! Time interpolation for pressure at surface, still on vertical source grid
	366	DO i=1,klon
	367	psurf_day(i) = &
[1269]	368	psurf_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
[1179]	369	(psurf_year(i,im2,id_aero) - psurf_year(i,im,id_aero))
	370
	371	pi_psurf_day(i) = &
[1269]	372	pi_psurf_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
[1179]	373	(pi_psurf_year(i,im2,id_aero) - pi_psurf_year(i,im,id_aero))
	374	END DO
	375
	376	! Time interpolation for the load, still on vertical source grid
	377	DO i=1,klon
	378	load_src(i) = &
[1269]	379	load_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
[1179]	380	(load_year(i,im2,id_aero) - load_year(i,im,id_aero))
	381
	382	pi_load_src(i) = &
[1269]	383	pi_load_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
[1179]	384	(pi_load_year(i,im2,id_aero) - pi_load_year(i,im,id_aero))
	385	END DO
	386
	387	!****************************************************************************************
	388	! 3) Interpolate to the model vertical grid (target grid)
	389	!
	390	!****************************************************************************************
	391
	392	IF (vert_interp) THEN
	393
	394	! - Interpolate variable tmp1 (on source grid) to var_day (on target grid)
	395	!********************************************************************************
	396	! a) calculate pression at vertical levels for the source grid using the
	397	! hybrid-sigma coordinates ap and b and the surface pressure, variables from file.
	398	DO k = 1, klev_src
	399	DO i = 1, klon
	400	pplay_src(i,k)= pt_ap(k) + pt_b(k)*psurf_day(i)
	401	END DO
	402	END DO
	403
	404	IF (debug) THEN
	405	CALL writefield_phy('psurf_day_src',psurf_day(:),1)
	406	CALL writefield_phy('pplay_src',pplay_src(:,:),klev_src)
	407	CALL writefield_phy('pplay',pplay(:,:),klev)
	408	CALL writefield_phy('day_src',tmp1,klev_src)
	409	CALL writefield_phy('pi_day_src',tmp2,klev_src)
	410	END IF
	411
	412	! b) vertical interpolation on pressure leveles
	413	CALL pres2lev(tmp1(:,:), var_day(:,:,id_aero), klev_src, klev, pplay_src, pplay, &
	414	1, klon, .FALSE.)
	415
	416	IF (debug) CALL writefield_phy('day_tgt',var_day(:,:,id_aero),klev)
	417
	418	! c) adjust to conserve total aerosol mass load in the vertical pillar
	419	! Calculate the load in the actual pillar and compare with the load
	420	! read from aerosol file.
	421
	422	! Find the pressure difference in each model layer
	423	DO k = 1, klev
	424	DO i = 1, klon
	425	delp(i,k) = paprs(i,k) - paprs (i,k+1)
	426	END DO
	427	END DO
	428
	429	! Find the mass load in the actual pillar, on target grid
	430	load_tgt(:) = 0.
	431	DO k= 1, klev
	432	DO i = 1, klon
[1216]	433	zrho = pplay(i,k)/t_seri(i,k)/RD ! [kg/m3]
	434	volm = var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
	435	load_tgt(i) = load_tgt(i) + 1/RG * volm *delp(i,k)
[1179]	436	END DO
	437	END DO
	438
	439	! Adjust, uniform
	440	DO k = 1, klev
	441	DO i = 1, klon
	442	var_day(i,k,id_aero) = var_day(i,k,id_aero)*load_src(i)/load_tgt(i)
	443	END DO
	444	END DO
	445
	446	IF (debug) THEN
	447	load_tgt_test(:) = 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 = var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
	452	load_tgt_test(i) = load_tgt_test(i) + 1/RG * volm*delp(i,k)
[1179]	453	END DO
	454	END DO
	455
	456	CALL writefield_phy('day_tgt2',var_day(:,:,id_aero),klev)
	457	CALL writefield_phy('load_tgt',load_tgt(:),1)
	458	CALL writefield_phy('load_tgt_test',load_tgt_test(:),1)
	459	CALL writefield_phy('load_src',load_src(:),1)
	460	END IF
	461
	462	! - Interpolate variable tmp2 (source grid) to pi_var_day (target grid)
	463	!********************************************************************************
	464	! a) calculate pression at vertical levels at source grid
	465	DO k = 1, klev_src
	466	DO i = 1, klon
	467	pplay_src(i,k)= pt_ap(k) + pt_b(k)*pi_psurf_day(i)
	468	END DO
	469	END DO
	470
	471	IF (debug) THEN
	472	CALL writefield_phy('pi_psurf_day_src',pi_psurf_day(:),1)
	473	CALL writefield_phy('pi_pplay_src',pplay_src(:,:),klev_src)
	474	END IF
	475
	476	! b) vertical interpolation on pressure leveles
	477	CALL pres2lev(tmp2(:,:), pi_var_day(:,:,id_aero), klev_src, klev, pplay_src, pplay, &
	478	1, klon, .FALSE.)
	479
	480	IF (debug) CALL writefield_phy('pi_day_tgt',pi_var_day(:,:,id_aero),klev)
	481
	482	! c) adjust to conserve total aerosol mass load in the vertical pillar
	483	! Calculate the load in the actual pillar and compare with the load
	484	! read from aerosol file.
	485
	486	! Find the load in the actual pillar, on target grid
	487	load_tgt(:) = 0.
	488	DO k = 1, klev
	489	DO i = 1, klon
[1216]	490	zrho = pplay(i,k)/t_seri(i,k)/RD ! [kg/m3]
	491	volm = pi_var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
	492	load_tgt(i) = load_tgt(i) + 1/RG * volm * delp(i,k)
[1179]	493	END DO
	494	END DO
	495
	496	DO k = 1, klev
	497	DO i = 1, klon
	498	pi_var_day(i,k,id_aero) = pi_var_day(i,k,id_aero)*pi_load_src(i)/load_tgt(i)
	499	END DO
	500	END DO
	501
	502	IF (debug) THEN
	503	load_tgt_test(:) = 0.
	504	DO k = 1, klev
	505	DO i = 1, klon
[1216]	506	zrho = pplay(i,k)/t_seri(i,k)/RD ! [kg/m3]
	507	volm = pi_var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
	508	load_tgt_test(i) = load_tgt_test(i) + 1/RG * volm * delp(i,k)
[1179]	509	END DO
	510	END DO
	511	CALL writefield_phy('pi_day_tgt2',pi_var_day(:,:,id_aero),klev)
	512	CALL writefield_phy('pi_load_tgt',load_tgt(:),1)
	513	CALL writefield_phy('pi_load_tgt_test',load_tgt_test(:),1)
	514	CALL writefield_phy('pi_load_src',pi_load_src(:),1)
	515	END IF
	516
	517
	518	ELSE ! No vertical interpolation done
	519
	520	var_day(:,:,id_aero) = tmp1(:,:)
	521	pi_var_day(:,:,id_aero) = tmp2(:,:)
	522
	523	END IF ! vert_interp
	524
	525
	526	! Deallocation
	527	DEALLOCATE(tmp1, tmp2, pplay_src, stat=ierr)
	528
	529	!****************************************************************************************
	530	! 4) Test for negative mass values
	531	!
	532	!****************************************************************************************
	533	IF (MINVAL(var_day(:,:,id_aero)) < 0.) THEN
	534	DO k=1,klev
	535	DO i=1,klon
	536	! Test for var_day
	537	IF (var_day(i,k,id_aero) < 0.) THEN
[1265]	538	IF (jDay-day2 < 0.) WRITE(lunout,*) 'jDay-day2=',jDay-day2
[1179]	539	IF (var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero) < 0.) THEN
	540	WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', &
	541	trim(name_aero(id_aero)),'(i,k,im)=', &
	542	var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero)
	543	END IF
	544	WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
[1273]	545	WRITE(lunout,*) 'day1, day2, jDay = ', day1, day2, jDay
[1179]	546	CALL abort_gcm('readaerosol_interp','Error in interpolation 1',1)
	547	END IF
	548	END DO
	549	END DO
	550	END IF
	551
	552	IF (MINVAL(pi_var_day(:,:,id_aero)) < 0. ) THEN
	553	DO k=1, klev
	554	DO i=1,klon
	555	! Test for pi_var_day
	556	IF (pi_var_day(i,k,id_aero) < 0.) THEN
[1265]	557	IF (jDay-day2 < 0.) WRITE(lunout,*) 'jDay-day2=',jDay-day2
[1179]	558	IF (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero) < 0.) THEN
	559	WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', &
	560	trim(name_aero(id_aero)),'(i,k,im)=', &
	561	pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero)
	562	END IF
	563
	564	WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
	565	CALL abort_gcm('readaerosol_interp','Error in interpolation 2',1)
	566	END IF
	567	END DO
	568	END DO
	569	END IF
	570
	571	END IF ! lnewday
	572
	573	!****************************************************************************************
	574	! Copy output from saved variables
	575	!
	576	!****************************************************************************************
	577
	578	mass_out(:,:) = var_day(:,:,id_aero)
	579	pi_mass_out(:,:) = pi_var_day(:,:,id_aero)
	580
	581	END SUBROUTINE readaerosol_interp

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

Download in other formats:

Original Format