54 integer :: j,ix,jy,kz,ncolumn,numparttot
55 real :: gridarea(0:nymax-1),pp(nzmax),ylat,ylatp,ylatm,hzone,
ran1
56 real :: cosfactm,cosfactp,deltacol,dz1,dz2,dz,pnew,fractus
57 real,
parameter :: pih=pi/180.
58 real :: colmass(0:nxmax-1,0:nymax-1),colmasstotal,zposition
60 integer :: ixm,ixp,jym,jyp,indzm,indzp,in,indzh,i,jj
61 real :: pvpart,ddx,ddy,rddx,rddy,p1,p2,p3,p4,y1(2)
63 integer :: idummy = -11
71 nx_we(1)=max(int(xpoint1(1)),0)
72 nx_we(2)=min((int(xpoint2(1))+1),nxmin1)
73 ny_sn(1)=max(int(ypoint1(1)),0)
74 ny_sn(2)=min((int(ypoint2(1))+1),nymin1)
79 if (xglobal.and.sglobal.and.nglobal)
then
80 if ((nx_we(1).eq.0).and.(nx_we(2).eq.nxmin1).and. &
81 (ny_sn(1).eq.0).and.(ny_sn(2).eq.nymin1))
then
91 if (xglobal) nx_we(2)=min(nx_we(2),nx-2)
98 do jy=ny_sn(1),ny_sn(2)
99 ylat=ylat0+
real(jy)*dy
102 if ((ylatm.lt.0).and.(ylatp.gt.0.))
then
105 cosfactp=cos(ylatp*pih)*r_earth
106 cosfactm=cos(ylatm*pih)*r_earth
107 if (cosfactp.lt.cosfactm)
then
108 hzone=sqrt(r_earth**2-cosfactp**2)- &
109 sqrt(r_earth**2-cosfactm**2)
111 hzone=sqrt(r_earth**2-cosfactm**2)- &
112 sqrt(r_earth**2-cosfactp**2)
115 gridarea(jy)=2.*pi*r_earth*hzone*dx/360.
125 cosfactp=cos(ylatp*pih)*r_earth
126 hzone=sqrt(r_earth**2-cosfactm**2)- &
127 sqrt(r_earth**2-cosfactp**2)
128 gridarea(0)=2.*pi*r_earth*hzone*dx/360.
134 ylat=ylat0+
real(nymin1)*dy
138 cosfactm=cos(ylatm*pih)*r_earth
139 hzone=sqrt(r_earth**2-cosfactp**2)- &
140 sqrt(r_earth**2-cosfactm**2)
141 gridarea(nymin1)=2.*pi*r_earth*hzone*dx/360.
149 do jy=ny_sn(1),ny_sn(2)
150 do ix=nx_we(1),nx_we(2)
151 pp(1)=rho(ix,jy,1,1)*r_air*tt(ix,jy,1,1)
152 pp(nz)=rho(ix,jy,nz,1)*r_air*tt(ix,jy,nz,1)
153 colmass(ix,jy)=(pp(1)-pp(nz))/ga*gridarea(jy)
154 colmasstotal=colmasstotal+colmass(ix,jy)
158 write(*,*)
'Atm. mass: ',colmasstotal
161 if (ipin.eq.0) numpart=0
168 do jy=ny_sn(1),ny_sn(2)
169 ylat=ylat0+
real(jy)*dy
170 do ix=nx_we(1),nx_we(2)
171 ncolumn=nint(0.999*
real(npart(1))*colmass(ix,jy)/ &
173 if (ncolumn.eq.0) goto 30
174 if (ncolumn.gt.numcolumn) numcolumn=ncolumn
181 pp(kz)=rho(ix,jy,kz,1)*r_air*tt(ix,jy,kz,1)
185 deltacol=(pp(1)-pp(nz))/
real(ncolumn)
186 pnew=pp(1)+deltacol/2.
198 if (ncolumn.gt.20)
then
201 pnew=pp(1)-
ran1(idummy)*(pp(1)-pp(nz))
205 if ((pp(kz).ge.pnew).and.(pp(kz+1).lt.pnew))
then
215 xtra1(numpart+jj)=
real(ix)-0.5+
ran1(idummy)
216 if (ix.eq.0) xtra1(numpart+jj)=
ran1(idummy)
217 if (ix.eq.nxmin1) xtra1(numpart+jj)= &
218 real(nxmin1)-
ran1(idummy)
219 ytra1(numpart+jj)=
real(jy)-0.5+
ran1(idummy)
220 ztra1(numpart+jj)=(height(kz)*dz2+height(kz+1)*dz1)*dz
221 if (ztra1(numpart+jj).gt.height(nz)-0.5) &
222 ztra1(numpart+jj)=height(nz)-0.5
227 ixm=int(xtra1(numpart+jj))
228 jym=int(ytra1(numpart+jj))
231 ddx=xtra1(numpart+jj)-
real(ixm)
232 ddy=ytra1(numpart+jj)-
real(jym)
240 if (height(i).gt.ztra1(numpart+jj))
then
247 dz1=ztra1(numpart+jj)-height(indzm)
248 dz2=height(indzp)-ztra1(numpart+jj)
252 y1(in)=p1*pv(ixm,jym,indzh,1) &
253 +p2*pv(ixp,jym,indzh,1) &
254 +p3*pv(ixm,jyp,indzh,1) &
255 +p4*pv(ixp,jyp,indzh,1)
257 pvpart=(dz2*y1(1)+dz1*y1(2))*dz
258 if (ylat.lt.0.) pvpart=-1.*pvpart
264 if (((ztra1(numpart+jj).gt.3000.).and. &
265 (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1))
then
269 nclass(numpart+jj)=min(int(
ran1(idummy)* &
270 real(nclassunc))+1,nclassunc)
271 numparticlecount=numparticlecount+1
272 npoint(numpart+jj)=numparticlecount
273 idt(numpart+jj)=mintime
275 itramem(numpart+jj)=0
276 itrasplit(numpart+jj)=itra1(numpart+jj)+ldirect* &
278 xmass1(numpart+jj,1)=colmass(ix,jy)/
real(ncolumn)
279 if (mdomainfill.eq.2) xmass1(numpart+jj,1)= &
280 xmass1(numpart+jj,1)*pvpart*48./29.*ozonescale/10.**9
288 numparttot=numparttot+ncolumn
289 if (ipin.eq.0) numpart=numpart+jj
298 if (numpart.gt.maxpart)
then
299 write(*,*)
'numpart too large: change source in init_atm_mass.f'
300 write(*,*)
'numpart: ',numpart,
' maxpart: ',maxpart
304 xmassperparticle=colmasstotal/
real(numparttot)
311 if ((xtra1(j).lt.0.).or.(xtra1(j).ge.
real(nxmin1)).or. &
312 (ytra1(j).lt.0.).or.(ytra1(j).ge.
real(nymin1))) then
330 fractus=
real(numcolumn)/
real(nz)
331 write(*,*)
'Total number of particles at model start: ',numpart
332 write(*,*)
'Maximum number of particles per column: ',numcolumn
333 write(*,*)
'If ',fractus,
' <1, better use more particles'
334 fractus=sqrt(max(fractus,1.))/2.
336 do jy=ny_sn(1),ny_sn(2)
337 do ix=nx_we(1),nx_we(2)
338 ncolumn=nint(0.999/fractus*
real(npart(1))*colmass(ix,jy) &
340 if (ncolumn.gt.maxcolumn) stop
'maxcolumn too small'
341 if (ncolumn.eq.0) goto 80
349 if (ix.eq.nx_we(1)) numcolumn_we(1,jy)=ncolumn
350 if (ix.eq.nx_we(2)) numcolumn_we(2,jy)=ncolumn
351 if (jy.eq.ny_sn(1)) numcolumn_sn(1,ix)=ncolumn
352 if (jy.eq.ny_sn(2)) numcolumn_sn(2,ix)=ncolumn
359 pp(kz)=rho(ix,jy,kz,1)*r_air*tt(ix,jy,kz,1)
365 deltacol=(pp(1)-pp(nz))/
real(ncolumn)
366 pnew=pp(1)+deltacol/2.
370 if ((pp(kz).ge.pnew).and.(pp(kz+1).lt.pnew))
then
374 zposition=(height(kz)*dz2+height(kz+1)*dz1)*dz
375 if (zposition.gt.height(nz)-0.5) zposition=height(nz)-0.5
381 if (ix.eq.nx_we(1)) zcolumn_we(1,jy,j)=zposition
382 if (ix.eq.nx_we(2)) zcolumn_we(2,jy,j)=zposition
383 if (jy.eq.ny_sn(1)) zcolumn_sn(1,ix,j)=zposition
384 if (jy.eq.ny_sn(2)) zcolumn_sn(2,ix,j)=zposition
389 acc_mass_we(1,jy,j)=0.
390 acc_mass_we(2,jy,j)=0.
391 acc_mass_sn(1,jy,j)=0.
392 acc_mass_sn(2,jy,j)=0.
406 open(unitboundcond,file=path(2)(1:length(2))//
'boundcond.bin', &
408 read(unitboundcond) numcolumn_we,numcolumn_sn, &
409 zcolumn_we,zcolumn_sn,acc_mass_we,acc_mass_sn
subroutine init_domainfill