Context Navigation

readaerosol_interp.F90 @ 3639

Last change on this file since 3639 was 2840, checked in by oboucher, 8 years ago
Correcting the case when aerosol is zero (eg NO3 when missing from file)
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 23.3 KB

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

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

Download in other formats:

Original Format