Context Navigation

readreleases.f90 @ 3481cc1

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug

Last change on this file since 3481cc1 was 3481cc1, checked in by Ignacio Pisso <ip@…>, 5 years ago
move license from headers to a different file
Property mode set to `100644`
File size: 24.6 KB

Rev	Line
[e200b7a]	1	subroutine readreleases
	2
	3	!*****************************************************************************
	4	! *
	5	! This routine reads the release point specifications for the current *
	6	! model run. Several release points can be used at the same time. *
	7	! *
	8	! Author: A. Stohl *
	9	! *
	10	! 18 May 1996 *
	11	! *
	12	! Update: 29 January 2001 *
	13	! Release altitude can be either in magl or masl *
[b4d29ce]	14	! HSO, 12 August 2013
	15	! Added optional namelist input
[e200b7a]	16	! *
	17	!*****************************************************************************
	18	! *
	19	! Variables: *
	20	! decay decay constant of species *
	21	! dquer [um] mean particle diameters *
	22	! dsigma e.g. dsigma=10 or dsigma=0.1 means that 68% of the mass*
	23	! are between 0.1dquer and 10dquer *
	24	! ireleasestart, ireleaseend [s] starting time and ending time of each *
	25	! release *
	26	! kindz 1: zpoint is in m agl, 2: zpoint is in m asl, 3: zpoint*
	27	! is in hPa *
	28	! npart number of particles to be released *
	29	! nspec number of species to be released *
	30	! density [kg/m3] density of the particles *
	31	! rm [s/m] Mesophyll resistance *
	32	! species name of species *
	33	! xmass total mass of each species *
	34	! xpoint1,ypoint1 geograf. coordinates of lower left corner of release *
	35	! area *
	36	! xpoint2,ypoint2 geograf. coordinates of upper right corner of release *
	37	! area *
[341f4b7]	38	! weta_gas, wetb_gas parameters for below-cloud scavenging (gas) *
	39	! crain_aero, csnow_aero parameters for below-cloud scavenging (aerosol) *
	40	! ccn_aero, in_aero parameters for in-cloud scavenging (aerosol) *
[e200b7a]	41	! zpoint1,zpoint2 height range, over which release takes place *
[f13406c]	42	! num_min_discrete if less, release cannot be randomized and happens at *
	43	! time mid-point of release interval *
[8a65cb0]	44	! lroot true if serial version, or if MPI and root process *
[e200b7a]	45	! *
	46	!*****************************************************************************
	47
	48	use point_mod
	49	use xmass_mod
	50	use par_mod
	51	use com_mod
	52
	53	implicit none
	54
[5184a7c]	55	integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat,irel,ispc,nsettle
[f13406c]	56	integer,parameter :: num_min_discrete=100
[e200b7a]	57	real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
[f13406c]	58	real(kind=dp) :: jul1,jul2,julm,juldate
[5184a7c]	59	real,parameter :: eps2=1.e-9
[e200b7a]	60	character(len=50) :: line
	61	logical :: old
	62
[b4d29ce]	63	! help variables for namelist reading
	64	integer :: numpoints, parts, readerror
	65	integer*2 :: zkind
	66	integer :: idate1, itime1, idate2, itime2
	67	real :: lon1,lon2,lat1,lat2,z1,z2
	68	character*40 :: comment
	69	integer,parameter :: unitreleasesout=2
	70	real,allocatable, dimension (:) :: mass
	71	integer,allocatable, dimension (:) :: specnum_rel,specnum_rel2
	72
	73	! declare namelists
	74	namelist /releases_ctrl/ &
[5184a7c]	75	nspec, &
	76	specnum_rel
[b4d29ce]	77
	78	namelist /release/ &
[5184a7c]	79	idate1, itime1, &
	80	idate2, itime2, &
	81	lon1, lon2, &
	82	lat1, lat2, &
	83	z1, z2, &
	84	zkind, &
	85	mass, &
	86	parts, &
	87	comment
[e200b7a]	88
[b4d29ce]	89	numpoint=0
[e200b7a]	90
[b4d29ce]	91	! allocate with maxspec for first input loop
	92	allocate(mass(maxspec),stat=stat)
	93	if (stat.ne.0) write(,)'ERROR: could not allocate mass'
	94	allocate(specnum_rel(maxspec),stat=stat)
	95	if (stat.ne.0) write(,)'ERROR: could not allocate specnum_rel'
[e200b7a]	96
[b4d29ce]	97	! presetting namelist releases_ctrl
	98	nspec = -1 ! use negative value to determine failed namelist input
	99	specnum_rel = 0
	100
	101	!sec, read release to find how many releasepoints should be allocated
	102	open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',form='formatted',err=999)
	103
	104	! check if namelist input provided
	105	read(unitreleases,releases_ctrl,iostat=readerror)
	106
	107	! prepare namelist output if requested
[8a65cb0]	108	if (nmlout.and.lroot) then
	109	open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',access='append',status='replace',err=1000)
[e200b7a]	110	endif
	111
[b4d29ce]	112	if ((readerror.ne.0).or.(nspec.lt.0)) then
[e200b7a]	113
[5184a7c]	114	! no namelist format, close file and allow reopening in old format
[b4d29ce]	115	close(unitreleases)
	116	open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',err=999)
[e200b7a]	117
[b4d29ce]	118	readerror=1 ! indicates old format
[e200b7a]	119
[5184a7c]	120	! Check the format of the RELEASES file (either in free format,
	121	! or using a formatted mask)
	122	! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
	123	!**************************************************************************
[e200b7a]	124
[b4d29ce]	125	call skplin(12,unitreleases)
	126	read (unitreleases,901) line
	127	901 format (a)
	128	if (index(line,'Total') .eq. 0) then
	129	old = .false.
	130	else
	131	old = .true.
	132	endif
	133	rewind(unitreleases)
[e200b7a]	134
[b4d29ce]	135
[5184a7c]	136	! Skip first 11 lines (file header)
	137	!**********************************
[b4d29ce]	138
	139	call skplin(11,unitreleases)
	140
	141	read(unitreleases,*,err=998) nspec
	142	if (old) call skplin(2,unitreleases)
	143	do i=1,nspec
	144	read(unitreleases,*,err=998) specnum_rel(i)
	145	if (old) call skplin(2,unitreleases)
	146	end do
	147
	148	numpoint=0
	149	100 numpoint=numpoint+1
	150	read(unitreleases,*,end=25)
	151	read(unitreleases,*,err=998,end=25) idum,idum
	152	if (old) call skplin(2,unitreleases)
	153	read(unitreleases,*,err=998) idum,idum
	154	if (old) call skplin(2,unitreleases)
	155	read(unitreleases,*,err=998) xdum
	156	if (old) call skplin(2,unitreleases)
	157	read(unitreleases,*,err=998) xdum
	158	if (old) call skplin(2,unitreleases)
[e200b7a]	159	read(unitreleases,*,err=998) xdum
	160	if (old) call skplin(2,unitreleases)
[b4d29ce]	161	read(unitreleases,*,err=998) xdum
	162	if (old) call skplin(2,unitreleases)
	163	read(unitreleases,*,err=998) idum
	164	if (old) call skplin(2,unitreleases)
	165	read(unitreleases,*,err=998) xdum
	166	if (old) call skplin(2,unitreleases)
	167	read(unitreleases,*,err=998) xdum
	168	if (old) call skplin(2,unitreleases)
	169	read(unitreleases,*,err=998) idum
	170	if (old) call skplin(2,unitreleases)
	171	do i=1,nspec
	172	read(unitreleases,*,err=998) xdum
	173	if (old) call skplin(2,unitreleases)
	174	end do
[5184a7c]	175	!save compoint only for the first 1000 release points
[b4d29ce]	176	read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
	177	if (old) call skplin(1,unitreleases)
	178
	179	goto 100
	180
	181	25 numpoint=numpoint-1
	182
	183	else
	184
	185	readerror=0
	186	do while (readerror.eq.0)
	187	idate1=-1
	188	read(unitreleases,release,iostat=readerror)
	189	if ((idate1.lt.0).or.(readerror.ne.0)) then
	190	readerror=1
	191	else
	192	numpoint=numpoint+1
	193	endif
	194	end do
	195	readerror=0
	196	endif ! if namelist input
	197
	198	rewind(unitreleases)
	199
[4c64400]	200	if (nspec.gt.maxspec) goto 994
	201
[b4d29ce]	202	! allocate arrays of matching size for number of species (namelist output)
	203	deallocate(mass)
	204	allocate(mass(nspec),stat=stat)
	205	if (stat.ne.0) write(,)'ERROR: could not allocate mass'
	206	allocate(specnum_rel2(nspec),stat=stat)
	207	if (stat.ne.0) write(,)'ERROR: could not allocate specnum_rel2'
	208	specnum_rel2=specnum_rel(1:nspec)
[b0434e1]	209	deallocate(specnum_rel)
[5184a7c]	210	! eso: BUG, crashes here for nspec=12 and maxspec=6,
	211	! TODO: catch error and exit
[b4d29ce]	212	allocate(specnum_rel(nspec),stat=stat)
	213	if (stat.ne.0) write(,)'ERROR: could not allocate specnum_rel'
	214	specnum_rel=specnum_rel2
	215	deallocate(specnum_rel2)
	216
	217	!allocate memory for numpoint releaspoints
	218	allocate(ireleasestart(numpoint),stat=stat)
	219	if (stat.ne.0) write(,)'ERROR: could not allocate ireleasestart'
	220	allocate(ireleaseend(numpoint),stat=stat)
	221	if (stat.ne.0) write(,)'ERROR: could not allocate ireleaseend'
	222	allocate(xpoint1(numpoint),stat=stat)
	223	if (stat.ne.0) write(,)'ERROR: could not allocate xpoint1'
	224	allocate(xpoint2(numpoint),stat=stat)
	225	if (stat.ne.0) write(,)'ERROR: could not allocate xpoint2'
	226	allocate(ypoint1(numpoint),stat=stat)
	227	if (stat.ne.0) write(,)'ERROR: could not allocate ypoint1'
	228	allocate(ypoint2(numpoint),stat=stat)
	229	if (stat.ne.0) write(,)'ERROR: could not allocate ypoint2'
	230	allocate(zpoint1(numpoint),stat=stat)
	231	if (stat.ne.0) write(,)'ERROR: could not allocate zpoint1'
	232	allocate(zpoint2(numpoint),stat=stat)
	233	if (stat.ne.0) write(,)'ERROR: could not allocate zpoint2'
	234	allocate(kindz(numpoint),stat=stat)
	235	if (stat.ne.0) write(,)'ERROR: could not allocate kindz'
	236	allocate(xmass(numpoint,maxspec),stat=stat)
	237	if (stat.ne.0) write(,)'ERROR: could not allocate xmass'
	238	allocate(rho_rel(numpoint),stat=stat)
	239	if (stat.ne.0) write(,)'ERROR: could not allocate rho_rel'
	240	allocate(npart(numpoint),stat=stat)
	241	if (stat.ne.0) write(,)'ERROR: could not allocate npart'
	242	allocate(xmasssave(numpoint),stat=stat)
	243	if (stat.ne.0) write(,)'ERROR: could not allocate xmasssave'
	244
[8a65cb0]	245	if (lroot) write (,) 'Releasepoints : ', numpoint
[b4d29ce]	246
	247	do i=1,numpoint
	248	xmasssave(i)=0.
[e200b7a]	249	end do
	250
	251	!now save the information
	252	DEP=.false.
	253	DRYDEP=.false.
	254	WETDEP=.false.
	255	OHREA=.false.
	256	do i=1,maxspec
	257	DRYDEPSPEC(i)=.false.
[c8fc724]	258	WETDEPSPEC(i)=.false.
[e200b7a]	259	end do
	260
[b4d29ce]	261	if (readerror.ne.0) then
[5184a7c]	262	! Skip first 11 lines (file header)
	263	!**********************************
[e200b7a]	264
[b4d29ce]	265	call skplin(11,unitreleases)
[e200b7a]	266
[5184a7c]	267	! Assign species-specific parameters needed for physical processes
	268	!*************************************************************************
[e200b7a]	269
[b4d29ce]	270	read(unitreleases,*,err=998) nspec
	271	if (nspec.gt.maxspec) goto 994
	272	if (old) call skplin(2,unitreleases)
	273	endif
[e200b7a]	274
[b4d29ce]	275	! namelist output
[8a65cb0]	276	if (nmlout.and.lroot) then
[b4d29ce]	277	write(unitreleasesout,nml=releases_ctrl)
	278	endif
[e200b7a]	279
	280	do i=1,nspec
[b4d29ce]	281	if (readerror.ne.0) then
	282	read(unitreleases,*,err=998) specnum_rel(i)
	283	if (old) call skplin(2,unitreleases)
	284	call readspecies(specnum_rel(i),i)
	285	else
	286	call readspecies(specnum_rel(i),i)
	287	endif
[e200b7a]	288
[5184a7c]	289	! For backward runs, only 1 species is allowed
	290	!*********************************************
[e200b7a]	291
[5184a7c]	292	!if ((ldirect.lt.0).and.(nspec.gt.1)) then
	293	!write(,) '#####################################################'
	294	!write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
	295	!write(,) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
	296	!write(,) '#####################################################'
	297	! stop
	298	!endif
[e200b7a]	299
[5184a7c]	300	! Molecular weight
	301	!*****************
[e200b7a]	302
[b4d29ce]	303	if (((iout.eq.2).or.(iout.eq.3)).and.(weightmolar(i).lt.0.)) then
[e200b7a]	304	write(,) 'For mixing ratio output, valid molar weight'
	305	write(,) 'must be specified for all simulated species.'
	306	write(,) 'Check table SPECIES or choose concentration'
	307	write(,) 'output instead if molar weight is not known.'
	308	stop
	309	endif
	310
[5184a7c]	311	! Radioactive decay
	312	!******************
[e200b7a]	313
	314	decay(i)=0.693147/decay(i) !conversion half life to decay constant
	315
	316
[5184a7c]	317	! Dry deposition of gases
	318	!************************
[e200b7a]	319
[b4d29ce]	320	if (reldiff(i).gt.0.) rm(i)=1./(henry(i)/3000.+100.*f0(i)) ! mesophyll resistance
[e200b7a]	321
[5184a7c]	322	! Dry deposition of particles
	323	!****************************
[e200b7a]	324
	325	vsetaver(i)=0.
	326	cunningham(i)=0.
	327	dquer(i)=dquer(i)*1000000. ! Conversion m to um
[b4d29ce]	328	if (density(i).gt.0.) then ! Additional parameters
	329	call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
[e200b7a]	330	do j=1,ni
	331	fract(i,j)=fracth(j)
	332	schmi(i,j)=schmih(j)
	333	vset(i,j)=vsh(j)
	334	cunningham(i)=cunningham(i)+cun*fract(i,j)
	335	vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
	336	end do
[8a65cb0]	337	if (lroot) write(,) 'Average settling velocity: ',i,vsetaver(i)
[e200b7a]	338	endif
	339
[5184a7c]	340	! Dry deposition for constant deposition velocity
	341	!************************************************
[e200b7a]	342
	343	dryvel(i)=dryvel(i)*0.01 ! conversion to m/s
	344
[8a65cb0]	345	! Check if wet deposition or OH reaction shall be calculated
	346	!***********************************************************
[05cf28d]	347
	348	! ESO 04.2016 check for below-cloud scavenging (gas or aerosol)
	349	if ((dquer(i).le.0..and.(weta_gas(i).gt.0. .or. wetb_gas(i).gt.0.)) .or. &
	350	&(dquer(i).gt.0. .and. (crain_aero(i) .gt. 0. .or. csnow_aero(i).gt.0.))) then
[e200b7a]	351	WETDEP=.true.
[c8fc724]	352	WETDEPSPEC(i)=.true.
[8a65cb0]	353	if (lroot) then
[05cf28d]	354	write (,) ' Below-cloud scavenging: ON'
[5184a7c]	355	! write (,) 'Below-cloud scavenging coefficients: ',weta(i),i
[8a65cb0]	356	end if
[4fbe7a5]	357	else
[05cf28d]	358	if (lroot) write (,) ' Below-cloud scavenging: OFF'
[e200b7a]	359	endif
[5184a7c]	360
	361	! NIK 31.01.2013 + 10.12.2013 + 15.02.2015
[05cf28d]	362	if (dquer(i).gt.0..and.(ccn_aero(i).gt.0. .or. in_aero(i).gt.0.)) then
[4fbe7a5]	363	WETDEP=.true.
[c8fc724]	364	WETDEPSPEC(i)=.true.
[8a65cb0]	365	if (lroot) then
[05cf28d]	366	write (,) ' In-cloud scavenging: ON'
[5184a7c]	367	! write (,) 'In-cloud scavenging coefficients: ',&
	368	! &ccn_aero(i),in_aero(i),i !,wetc_in(i), wetd_in(i),i
[8a65cb0]	369	end if
	370	else
[05cf28d]	371	if (lroot) write (,) ' In-cloud scavenging: OFF'
[4fbe7a5]	372	endif
	373
[8a65cb0]	374	if (ohcconst(i).gt.0.) then
[e200b7a]	375	OHREA=.true.
[05cf28d]	376	if (lroot) write (,) ' OHreaction switched on: ',ohcconst(i),i
[e200b7a]	377	endif
	378
[b4d29ce]	379	if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then
[e200b7a]	380	DRYDEP=.true.
	381	DRYDEPSPEC(i)=.true.
	382	endif
	383
[c8fc724]	384	end do ! end loop over species
[e200b7a]	385
[b4d29ce]	386	if (WETDEP.or.DRYDEP) DEP=.true.
[e200b7a]	387
	388	! Read specifications for each release point
	389	!*******************************************
[b4d29ce]	390	numpoints=numpoint
[e200b7a]	391	numpoint=0
	392	numpartmax=0
	393	releaserate=0.
[5184a7c]	394	101 numpoint=numpoint+1
[b4d29ce]	395
	396	if (readerror.lt.1) then ! reading namelist format
	397
	398	if (numpoint.gt.numpoints) goto 250
	399	zkind = 1
	400	mass = 0
	401	parts = 0
	402	comment = ' '
	403	read(unitreleases,release,iostat=readerror)
	404	id1=idate1
	405	it1=itime1
	406	id2=idate2
	407	it2=itime2
	408	xpoint1(numpoint)=lon1
	409	xpoint2(numpoint)=lon2
	410	ypoint1(numpoint)=lat1
	411	ypoint2(numpoint)=lat2
	412	zpoint1(numpoint)=z1
	413	zpoint2(numpoint)=z2
	414	kindz(numpoint)=zkind
	415	do i=1,nspec
	416	xmass(numpoint,i)=mass(i)
	417	end do
	418	npart(numpoint)=parts
	419	compoint(min(1001,numpoint))=comment
	420
[5184a7c]	421	! namelist output
[8a65cb0]	422	if (nmlout.and.lroot) then
[b4d29ce]	423	write(unitreleasesout,nml=release)
	424	endif
	425
[e200b7a]	426	else
	427
[b4d29ce]	428	read(unitreleases,*,end=250)
	429	read(unitreleases,*,err=998,end=250) id1,it1
	430	if (old) call skplin(2,unitreleases)
	431	read(unitreleases,*,err=998) id2,it2
	432	if (old) call skplin(2,unitreleases)
	433	read(unitreleases,*,err=998) xpoint1(numpoint)
	434	if (old) call skplin(2,unitreleases)
	435	read(unitreleases,*,err=998) ypoint1(numpoint)
	436	if (old) call skplin(2,unitreleases)
	437	read(unitreleases,*,err=998) xpoint2(numpoint)
	438	if (old) call skplin(2,unitreleases)
	439	read(unitreleases,*,err=998) ypoint2(numpoint)
	440	if (old) call skplin(2,unitreleases)
	441	read(unitreleases,*,err=998) kindz(numpoint)
	442	if (old) call skplin(2,unitreleases)
	443	read(unitreleases,*,err=998) zpoint1(numpoint)
	444	if (old) call skplin(2,unitreleases)
	445	read(unitreleases,*,err=998) zpoint2(numpoint)
	446	if (old) call skplin(2,unitreleases)
	447	read(unitreleases,*,err=998) npart(numpoint)
	448	if (old) call skplin(2,unitreleases)
	449	do i=1,nspec
	450	read(unitreleases,*,err=998) xmass(numpoint,i)
	451	if (old) call skplin(2,unitreleases)
	452	mass(i)=xmass(numpoint,i)
	453	end do
[5184a7c]	454	!save compoint only for the first 1000 release points
[b4d29ce]	455	if (numpoint.le.1000) then
	456	read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
	457	comment=compoint(numpoint)(1:40)
	458	else
	459	read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
	460	comment=compoint(1001)(1:40)
	461	endif
	462	if (old) call skplin(1,unitreleases)
	463
[5184a7c]	464	! namelist output
[8a65cb0]	465	if (nmlout.and.lroot) then
[b4d29ce]	466	idate1=id1
	467	itime1=it1
	468	idate2=id2
	469	itime2=it2
	470	lon1=xpoint1(numpoint)
	471	lon2=xpoint2(numpoint)
	472	lat1=ypoint1(numpoint)
	473	lat2=ypoint2(numpoint)
	474	z1=zpoint1(numpoint)
	475	z2=zpoint2(numpoint)
	476	zkind=kindz(numpoint)
	477	parts=npart(numpoint)
	478	write(unitreleasesout,nml=release)
	479	endif
	480
	481	if (numpoint.le.1000) then
	482	if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
	483	(xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
	484	(compoint(numpoint)(1:8).eq.' ')) goto 250
	485	else
	486	if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
	487	(xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
	488	endif
	489
	490	endif ! if namelist format
[e200b7a]	491
	492	! If a release point contains no particles, stop and issue error message
	493	!***********************************************************************
	494
	495	if (npart(numpoint).eq.0) then
	496	write(,) 'FLEXPART MODEL ERROR'
	497	write(,) 'RELEASES file is corrupt.'
	498	write(,) 'At least for one release point, there are zero'
	499	write(,) 'particles released. Make changes to RELEASES.'
	500	stop
	501	endif
	502
[8ee24a5]	503	! If FLEXPART is run for backward deposition force zpoint
	504	!*********************************************************************
	505	if (WETBKDEP) then
[5184a7c]	506	zpoint1(numpoint)=0.
	507	zpoint2(numpoint)=20000.
	508	kindz(numpoint)=1
[8ee24a5]	509	endif
	510	if (DRYBKDEP) then
[5184a7c]	511	zpoint1(numpoint)=0.
	512	zpoint2(numpoint)=2.*href
	513	kindz(numpoint)=1
[8ee24a5]	514	endif
	515
	516
[e200b7a]	517	! Check whether x coordinates of release point are within model domain
	518	!*********************************************************************
	519
[5184a7c]	520	if (xpoint1(numpoint).lt.xlon0) &
	521	xpoint1(numpoint)=xpoint1(numpoint)+360.
	522	if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) &
	523	xpoint1(numpoint)=xpoint1(numpoint)-360.
	524	if (xpoint2(numpoint).lt.xlon0) &
	525	xpoint2(numpoint)=xpoint2(numpoint)+360.
	526	if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) &
	527	xpoint2(numpoint)=xpoint2(numpoint)-360.
[e200b7a]	528
	529	! Determine relative beginning and ending times of particle release
	530	!******************************************************************
	531
	532	jul1=juldate(id1,it1)
	533	jul2=juldate(id2,it2)
[f13406c]	534	julm=(jul1+jul2)/2.
[e200b7a]	535	if (jul1.gt.jul2) then
	536	write(,) 'FLEXPART MODEL ERROR'
	537	write(,) 'Release stops before it begins.'
	538	write(,) 'Make changes to file RELEASES.'
	539	stop
	540	endif
	541	if (mdomainfill.eq.0) then ! no domain filling
	542	if (ldirect.eq.1) then
	543	if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
	544	write(,) 'FLEXPART MODEL ERROR'
	545	write(,) 'Release starts before simulation begins or ends'
	546	write(,) 'after simulation stops.'
	547	write(,) 'Make files COMMAND and RELEASES consistent.'
	548	stop
	549	endif
[f13406c]	550	if (npart(numpoint).gt.num_min_discrete) then
	551	ireleasestart(numpoint)=int((jul1-bdate)*86400.)
	552	ireleaseend(numpoint)=int((jul2-bdate)*86400.)
	553	else
	554	ireleasestart(numpoint)=int((julm-bdate)*86400.)
	555	ireleaseend(numpoint)=int((julm-bdate)*86400.)
	556	endif
[e200b7a]	557	else if (ldirect.eq.-1) then
	558	if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
	559	write(,) 'FLEXPART MODEL ERROR'
	560	write(,) 'Release starts before simulation begins or ends'
	561	write(,) 'after simulation stops.'
	562	write(,) 'Make files COMMAND and RELEASES consistent.'
	563	stop
	564	endif
[f13406c]	565	if (npart(numpoint).gt.num_min_discrete) then
	566	ireleasestart(numpoint)=int((jul1-bdate)*86400.)
	567	ireleaseend(numpoint)=int((jul2-bdate)*86400.)
	568	else
	569	ireleasestart(numpoint)=int((julm-bdate)*86400.)
	570	ireleaseend(numpoint)=int((julm-bdate)*86400.)
	571	endif
[e200b7a]	572	endif
	573	endif
	574
[4c64400]	575
[e200b7a]	576	! Determine the release rate (particles per second) and total number
	577	! of particles released during the simulation
	578	!*******************************************************************
	579
	580	if (ireleasestart(numpoint).ne.ireleaseend(numpoint)) then
	581	releaserate=releaserate+real(npart(numpoint))/ &
	582	real(ireleaseend(numpoint)-ireleasestart(numpoint))
	583	else
	584	releaserate=99999999
	585	endif
	586	numpartmax=numpartmax+npart(numpoint)
[b4d29ce]	587	goto 101
[e200b7a]	588
[5184a7c]	589	250 close(unitreleases)
[e200b7a]	590
[8a65cb0]	591	if (nmlout.and.lroot) then
[b4d29ce]	592	close(unitreleasesout)
	593	endif
	594
[8a65cb0]	595	if (lroot) write (,) 'Particles allocated (maxpart) : ',maxpart
	596	if (lroot) write (,) 'Particles released (numpartmax): ',numpartmax
[e200b7a]	597	numpoint=numpoint-1
	598
	599	if (ioutputforeachrelease.eq.1) then
	600	maxpointspec_act=numpoint
	601	else
	602	maxpointspec_act=1
	603	endif
	604
[5184a7c]	605	! Disable settling if more than 1 species at any release point
	606	! or if MQUASILAG and more than one species
	607	!*************************************************************
	608
	609	if (mquasilag.ne.0) then
	610	if (nspec.gt.1) lsettling=.false.
	611	else
	612	do irel=1,numpoint
	613	nsettle=0
	614	do ispc=1,nspec
	615	if (xmass(irel,ispc).gt.eps2) nsettle=nsettle+1
	616	end do
	617	if (nsettle.gt.1) lsettling=.false.
	618	end do
	619	end if
	620
	621	if (lroot) then
	622	if (.not.lsettling) then
	623	write(,) 'WARNING: more than 1 species per release point, settling &
	624	&disabled'
	625	end if
	626	end if
	627
[b4d29ce]	628	! Check, whether the total number of particles may exceed totally allowed
	629	! number of particles at some time during the simulation
	630	!************************************************************************
	631
[e200b7a]	632	if (releaserate.gt. &
	633	0.99*real(maxpart)/real(lage(nageclass))) then
[8a65cb0]	634	if (numpartmax.gt.maxpart.and.lroot) then
[5184a7c]	635	write(,) '#####################################################'
	636	write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
	637	write(,) '#### ####'
	638	write(,) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####'
	639	write(,) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####'
	640	write(,) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED ####'
	641	write(,) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####'
	642	write(,) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####'
	643	write(,) '#### HOWEVER, FLEXPART MAY HAVE TO STOP ####'
	644	write(,) '#### AT SOME TIME DURING THE SIMULATION. PLEASE ####'
	645	write(,) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT. ####'
	646	write(,) '#####################################################'
[e200b7a]	647	write(,) 'Maximum release rate may be: ',releaserate, &
	648	' particles per second'
	649	write(,) 'Maximum allowed release rate is: ', &
	650	real(maxpart)/real(lage(nageclass)),' particles per second'
	651	write(,) &
	652	'Total number of particles released during the simulation is: ', &
	653	numpartmax
	654	write(,) 'Maximum allowed number of particles is: ',maxpart
	655	endif
	656	endif
	657
[5184a7c]	658
[4c64400]	659	if (lroot) then
	660	write(*,FMT='(A,ES14.7)') ' Total mass released:', sum(xmass(1:numpoint,1:nspec))
	661	end if
	662
[e200b7a]	663	return
	664
[5184a7c]	665	994 write(,) '#####################################################'
[e200b7a]	666	write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
	667	write(,) '#### ####'
	668	write(,) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS####'
	669	write(,) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####'
	670	write(,) '#####################################################'
	671	stop
	672
[5184a7c]	673	998 write(,) '#####################################################'
[e200b7a]	674	write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
	675	write(,) '#### ####'
	676	write(,) '#### FATAL ERROR - FILE "RELEASES" IS ####'
	677	write(,) '#### CORRUPT. PLEASE CHECK YOUR INPUTS FOR ####'
	678	write(,) '#### MISTAKES OR GET A NEW "RELEASES"- ####'
	679	write(,) '#### FILE ... ####'
	680	write(,) '#####################################################'
	681	stop
	682
	683
[5184a7c]	684	999 write(,) '#####################################################'
[e200b7a]	685	write(,) ' FLEXPART MODEL SUBROUTINE READRELEASES: '
	686	write(,)
	687	write(,) 'FATAL ERROR - FILE CONTAINING PARTICLE RELEASE POINTS'
	688	write(,) 'POINTS IS NOT AVAILABLE OR YOU ARE NOT'
	689	write(,) 'PERMITTED FOR ANY ACCESS'
	690	write(,) '#####################################################'
	691	stop
	692
[b4d29ce]	693	1000 write(,) ' #### FLEXPART MODEL ERROR! FILE "RELEASES" #### '
	694	write(,) ' #### CANNOT BE OPENED IN THE DIRECTORY #### '
	695	write(*,'(a)') path(2)(1:length(2))
	696	stop
	697
[e200b7a]	698	end subroutine readreleases

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

Download in other formats:

Original Format