cat > param.grid << '\eof' c +++param.grid+++ cc major dimensions parameter(nx=51,ny=51,nz=51) parameter(nxyz=nx*ny*nz) '\eof' cat > param.sds << '\eof' c super droplets (SDs): only in cloudy volumes... c nppg - typically, twice the number of SDs per gridbox parameter(nppg=30) parameter(np=(nx-1)*(ny-1)*nz*nppg) '\eof' cat > aerosol.inc << '\eof' PARAMETER(nca=30,ncap=nca+1) common /drops/ conc(nca),qcc(nca),radc(nca) REAL rad0,mas0,mas0s,rho_a COMMON /aerosol1/ rad0(nca),mas0(nca),mas0s(ncap),rho_a COMMON /aerosol2/ iter(nca,nca),frac(nca,nca) COMMON /aerosol3/ s_activ(nca),r_activ(nca) common /chemistry/ rho_w,amol_w,rho_s,amol_s,vanhof common /kinetic/ alpha_c,alpha_t,delta_v,delta_t '\eof' cat > src.for << '\eof' program moist_convec ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Simple 3d periodic domain flow model with stretched vertical grid c c Author: W. Grabowski (grabow@ucar.edu) with routines taken c from early 1-PE version of the EULAG model from Piotr... c c with super-droplets and detailed droplet growth including activation c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cc include 'param.grid' include 'param.sds' include 'aerosol.inc' cc number of time steps: parameter(ntime0=1800*40) ! 30 min cc grid (for regular dzeta isosurfaces) data dx,dy,dz,dt /0.04,0.04,0.02,.025/ common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi dimension xx(nx),yy(ny),zz(nz) cc random noise variables (in pbl) dimension thn(nx,ny,nz),qvn(nx,ny,nz) common /rnd_forc/ thn_a,qvn_a,nz_noise CC MODEL VARIABLES cc thermodynamics dimension theta(nx,ny,nz),qv(nx,ny,nz),qc(nx,ny,nz),qr(nx,ny,nz) cc dynamics dimension ux(nx,ny,nz),uy(nx,ny,nz),uz(nx,ny,nz) ! current time level dimension uxp(nx,ny,nz),uyp(nx,ny,nz),uzp(nx,ny,nz) ! previous time level cc forces for model variables dimension ft(nx,ny,nz),fx(nx,ny,nz),fy(nx,ny,nz),fz(nx,ny,nz), * fqv(nx,ny,nz),fqc(nx,ny,nz),fqr(nx,ny,nz) cc advective velocities: dimension uxa(nx+1,ny,nz),uya(nx,ny+1,nz),uza(nx,ny,nz+1) dimension sup(nx,ny,nz) cc profiles common /strtch/ height(nz),gac(nz) common /prof_d/ rho0(nz),th0(nz),th_e(nz),ux_e(nz),uy_e(nz) common /prof_m/ qv_e(nz),tm_e(nz) common /prof_a/ tau(nz) cc required variables: dimension den(nx,ny,nz),p(nx,ny,nz) dimension pfx(nx,ny,nz),pfy(nx,ny,nz),pfz(nx,ny,nz) dimension scr1(nx,ny,nz),scr2(nx,ny,nz),scr3(nx,ny,nz) common /boundary/ thsrf,qvsrf,thtop,qvtop,thwal,qvwal CCC SD arrays c SD atributes: two-time level c SD atributes: position dimension xp(np),yp(np),zp(np) dimension xpn(np),ypn(np),zpn(np) c number of SDs, cloud water, multiplicity, radius dimension qcp(np),plic(np),radp(np) real radad(np),masad(np) c common /activ/ rad_sdb,dcon cc constants common /const/ gg,cp,rg,rv common/reference/ tt0,ee0 common/latent/hlatv,hlats ccc SDs: function used in critical radius calculation: fun(rr,rrae)=a_over_3b*(rr**3-rrae**3) - rr**2 funp(rr)=3.*a_over_3b*rr**2 - 2.*rr cc grid: time=0. dxi=1./dx dyi=1./dy dzi=1./dz dti=1./dt do i=1,nx xx(i)=float(i-1)*dx enddo do j=1,ny yy(j)=float(j-1)*dy enddo cc zz is regular (ustretched) dzeta grid do k=1,nz zz(k)=float(k-1)*dz enddo call zstrtch(zz,nz,dz) cc initialize moisture parameters: call moist_init cc initialize model profiles: call prof_init do i=1,nx do j=1,ny do k=1,nz den(i,j,k)=rho0(k)*gac(k) enddo enddo enddo cc top and bottom and wall conditions: call top_bot_wall cc initial fields call noise(thn,nx,ny,nz,nz-1) a=rg/rv c=hlatv/cp b=hlatv/(rv*tt0) d=hlatv/rv e=-cp/rg del=0.1 do k=1,nz thetme=th_e(k)/tm_e(k) prek=1.e5*thetme**e do i=1,nx do j=1,ny theta(i,j,k)=th_e(k) + del*thn(i,j,k) qv(i,j,k)=qv_e(k) qc(i,j,k)=0. qr(i,j,k)=0. ux(i,j,k)=ux_e(k) uy(i,j,k)=uy_e(k) uz(i,j,k)=0. uxp(i,j,k)=ux_e(k) uyp(i,j,k)=uy_e(k) uzp(i,j,k)=0. ft(i,j,k)=0. fx(i,j,k)=0. fy(i,j,k)=0. fz(i,j,k)=0. fqv(i,j,k)=0. fqc(i,j,k)=0. fqr(i,j,k)=0. p(i,j,k)=0. thi=1./theta(i,j,k) y=b*thetme*tt0*thi ees=ee0*exp(b-y) qvs=a*ees/(prek-ees) sup(i,j,k)=qv(i,j,k)/qvs - 1. enddo enddo enddo ccc initial CCN calls: call micro_set cc surface tension, vapor diffusivity, thermal conductivity: cc HERE we assume that temp and pre variability is not important temp=tm_e(1) sigma=76.1-0.155*(temp-273.16) sigma=sigma*1.e-3 ! in N/m difvp=2.11 * (temp/273.16)**1.94 * 101325./prek difvp=difvp*1.e-5 ! in m**2/s tcond=1.5e-11*temp**3 - 4.8e-8*temp**2 + 1.0e-4*temp - 3.9e-4 tcond=tcond ! in W/mK print*,'-- sigma,difvp,tcond: ',sigma,difvp,tcond ccc set initial condition: print*,'--- nca= ',nca rho_w=1.e3 pi=4.*atan(1.) do ic=1,nca cc s_active and r_activ: rs0=rad0(ic) sm0=mas0(ic) aaa=2*sigma/(rv*rho_w*temp) bbb=vanhof*sm0*amol_w/amol_s/(4./3.*pi*rho_w) a_over_3b = sqrt(aaa/3./bbb) print*,aaa,bbb cc drop: r_dr= 1./a_over_3b ! initial drop size... term1=aaa/r_dr r3=max(0.001*rs0**3, r_dr**3 - rs0**3) term2=bbb/r3 s_eq=exp(term1-term2) - 1. do it=1,5 r_dr=r_dr - fun(r_dr,rs0)/funp(r_dr) term1=aaa/r_dr r3=max(0.001*rs0**3, r_dr**3 - rs0**3) term2=bbb/r3 s_eq=exp(term1-term2) - 1. enddo c print*,' ' print 103,ic,rs0*1.e6,r_dr*1.e6,s_eq*1.e2 103 format(1x,i4,3e14.4,' ic,r0,ract,sact ') r_activ(ic)=r_dr s_activ(ic)=s_eq ccc eqlibrium with assumed initial RH: r_dr=r_activ(ic) supers=sup(1,1,1) ! only for S=const do it=1,15 ccc=alog(1.+supers) fun11=r_dr*(r_dr**3-rs0**3)*ccc - aaa*(r_dr**3-rs0**3) 1 + bbb*r_dr fun11p=(4.*r_dr**3-rs0**3)*ccc - 3.*aaa*r_dr**2 + bbb r_dr=r_dr - fun11/fun11p term1=aaa/r_dr r3=r_dr**3 - rs0**3 term2=bbb/r3 s_eq=exp(term1-term2) - 1. enddo print 104,supers*1.e2,rs0*1.e6,r_activ(ic)*1.e6, 1 s_activ(ic)*1.e2,r_dr*1.e6 104 format(1x,5e14.4,' sup,r0,r_act,s_act,r_qe ') radc(ic)=r_dr qcc(ic)=4./3.*pi*rho_w*(radc(ic)**3 - rad0(ic)**3) * conc(ic) print*,radc(ic),conc(ic),qcc(ic),' init rad,icon,qc' enddo ! loop over bins pi=4.*atan(1.) rho_w=1.e3 fact=4./3.*pi*rho_w npx=0 dxh=.5*dx dyh=.5*dy dzh=.5*dz do i=1,nx-1 do j=1,ny-1 do k=1,nz c SDS initialization: do ipp=npx+1,npx+nppg ic=ipp-npx cccc random position rand1=min(.999,max(.001,rand()))-.5 rand2=min(.999,max(.001,rand()))-.5 rand3=min(.999,max(.001,rand()))-.5 if(k.eq.1 .or. k.eq.nz) rand3=min(.999,max(.001,rand())) cc xp(ipp)=(i-1)*dx + rand1*dx yp(ipp)=(j-1)*dy + rand2*dy zp(ipp)=(k-1)*dz + rand3*dz if(k.eq.1) zp(ipp)= rand3*dz/2. if(k.eq.nz) zp(ipp)=(k-1)*dz - rand3*dz/2. plic(ipp)=conc(ic)*rho0(k) radp(ipp)=radc(ic) qcp(ipp)=qcc(ic) radad(ipp)=rad0(ic) masad(ipp)=mas0(ic) cccccccccc c if(ipp.le.nppg) then c if(ipp.le.10*nppg) then c print*,' i,j,k: ',i,j,k c print(284),ic,radp(ipp),plic(ipp),qcp(ipp) c print(285),radad(ipp),masad(ipp) c endif 284 format(2x,i5,3d16.5,' ic,radp,plic,qcp') 285 format(2x,2d16.5,' radad,masad') cccccccccc enddo npx=npx+nppg cc end do end do end do n=nx m=ny l=nz print*,'check npx: ',npx,(n-1)*(m-1)*l*nppg ccc check: c do ip=1,npx c i=(xp(ip)+dxh)/dx + 1 c j=(yp(ip)+dyh)/dy + 1 c k=(zp(ip)+dzh)/dz + 1 c c i=max(1,min(nx,i)) c j=max(1,min(ny,j)) c k=max(1,min(nz,k)) c cc if(i.eq.5 .and. j.eq.5 .and. k.eq.5) then c if(i.eq.15 .and. j.eq.15 .and. k.eq.15) then c print*,ip,radp(ip) c endif c enddo cc test cloud water: dxh=.5*dx dyh=.5*dy dzh=.5*dz do ip=1,npx i=(xp(ip)+dxh)/dx + 1 j=(yp(ip)+dyh)/dy + 1 k=(zp(ip)+dzh)/dz + 1 i=max(1,min(n,i)) j=max(1,min(m,j)) k=max(1,min(l,k)) qc(i,j,k)=qc(i,j,k)+qcp(ip) enddo cc periodicity do j=1,m-1 do k=1,l qc(n,j,k)=qc(1,j,k) enddo enddo do i=1,n do k=1,l qc(i,m,k)=qc(i,1,k) enddo enddo an=1.e20 ax=-1.e20 do i=1,n do j=1,m do k=1,l an=min(an,qc(i,j,k)) ax=max(ax,qc(i,j,k)) enddo enddo enddo print*,'--- initial min,max qc: ',an,ax Ccc plot initial fields: call diagno_1(ux,uy,uz,theta,nx,ny,nz,scr1,scr2,den) call diagno_2(qv,qc,qr,nx,ny,nz,scr1,scr2,den) call diagno_4(npx,xp,yp,zp,qcp,plic,radp) ccc save initial data: c write(17) time,nx,nz c write(17) ux,uy,uz,theta,qv,qc c write(18) time,npx,xp,yp,zp,qcp,plic,radp CCCC MARCH FORWARD IN TIME: ntime=ntime0 do itime=1,ntime ! TIME LOOP print*,'*** itime, time: ',itime,time cc extrapolate in time to get advective momentums: call velprd_1(ux,uxp,uxa,uy,uyp,uya,uz,uzp,uza,nx,ny,nz,den) cc transport SDs: call sd_adv_sed(ux,uxp,uy,uyp,uz,uzp,xp,yp,zp,radp,npx) c call sd_adv(ux,uxp,uy,uyp,uz,uzp,xp,yp,zp,npx) ccc eliminate SDs that fell out... c call remove(xp,yp,zp,qcp,plic,radp,npx) cc save previous velocities: do i=1,nxyz uxp(i,1,1)=ux(i,1,1) uyp(i,1,1)=uy(i,1,1) uzp(i,1,1)=uz(i,1,1) enddo cc add half of the force: do i=1,nxyz ux(i,1,1) = ux(i,1,1)+.5*dt*fx(i,1,1) uy(i,1,1) = uy(i,1,1)+.5*dt*fy(i,1,1) uz(i,1,1) = uz(i,1,1)+.5*dt*fz(i,1,1) cc theta and qv: 1st order theta(i,1,1)=theta(i,1,1)+dt*ft(i,1,1) qv(i,1,1) = qv(i,1,1)+dt*fqv(i,1,1) enddo CC ADVECTION: c liner: 1-iord=1, 0-iord prescribed inside mpdata liner=0 if(itime/6*6.eq.itime) liner=1 c if(itime/3*3.eq.itime) liner=1 c if(itime/1*1.eq.itime) liner=1 call mpdata_3d(uxa,uya,uza,ux,den,1,liner) call mpdata_3d(uxa,uya,uza,uy,den,2,liner) call mpdata_3d(uxa,uya,uza,uz,den,3,liner) cc get new thermodynamic variables: call mpdata_3d(uxa,uya,uza,theta,den,4,liner) call mpdata_3d(uxa,uya,uza,qv,den,5,liner) cc save velocities after advection into advective velocities: cc (used as scratch) do k=1,nz do j=1,ny do i=1,nx uxa(i,j,k)=ux(i,j,k) uya(i,j,k)=uy(i,j,k) uza(i,j,k)=uz(i,j,k) enddo enddo enddo ccc 1st, finish dynamics... cc derive qc for buoyancy: do i=1,nxyz qc(i,1,1)=0. enddo dxh=.5*dx dyh=.5*dy dzh=.5*dz do ip=1,npx i=(xp(ip)+dxh)/dx + 1 j=(yp(ip)+dyh)/dy + 1 k=(zp(ip)+dzh)/dz + 1 i=max(1,min(nx,i)) j=max(1,min(ny,j)) k=max(1,min(nz,k)) qc(i,j,k)=qc(i,j,k)+qcp(ip) enddo cc periodicity do j=1,ny-1 do k=1,nz qc(nx,j,k)=qc(1,j,k) enddo enddo do i=1,nx do k=1,nz qc(i,ny,k)=qc(i,1,k) enddo enddo cc add buoyancy epsb=rv/rg-1. do k=1,nz do j=1,ny do i=1,nx fx(i,j,k)=gg*( (theta(i,j,k)-th_e(k))/th0(k) * + epsb*(qv(i,j,k)-qv_e(k))-qc(i,j,k) ) enddo enddo enddo call integz(fx,fy,nx,ny,nz) cc apply do i=1,nxyz uz(i,1,1) = uz(i,1,1) +.5*dt*fx(i,1,1) enddo cc calculate pressure gradient force: c epp=1.e-6 epp=1.e-7 itp=100 call gcrk_2(p,fx,fy,fz,ux,uy,uz,nx,ny,nz,itp,epp) call prforc_2(p,fx,fy,fz,ux,uy,uz,nx,ny,nz) do i=1,nxyz ux(i,1,1)=fx(i,1,1) uy(i,1,1)=fy(i,1,1) uz(i,1,1)=fz(i,1,1) enddo cc calculate velocity forces (using saved velocities after advection): do k=1,nz do j=1,ny do i=1,nx fx(i,j,k)=(ux(i,j,k)-uxa(i,j,k)) *2./dt fy(i,j,k)=(uy(i,j,k)-uya(i,j,k)) *2./dt fz(i,j,k)=(uz(i,j,k)-uza(i,j,k)) *2./dt enddo enddo enddo ccc 2nd, derive new forces for thermodynamics cc derive new S: do k=1,nz thetme=th_e(k)/tm_e(k) prek=1.e5*thetme**e do i=1,nx do j=1,ny thi=1./theta(i,j,k) y=b*thetme*tt0*thi ees=ee0*exp(b-y) qvs=a*ees/(prek-ees) sup(i,j,k)=qv(i,j,k)/qvs-1. enddo enddo enddo call extr11(sup,nxyz,an,ax) print*,' sup min max: ',an,ax cc set forces to zero: do i=1,nxyz ft(i,1,1)=0. fqv(i,1,1)=0. enddo cc SDs thermodynamics; NOTE: theta and qv tendencies are brought back, cc SD properties are imediately updated... c call sd_phy_old(theta,qv,sup,ft,fqv, c 1 npx,xp,yp,zp,qcp,plic,radp) call sd_phy(theta,sup,ft,fqv,np,xp,yp,zp,qcp, 1 plic,radp,radad,masad) cc chamber boundaries: cc bottom and top: call boundry(theta,qv) cc walls call walls(theta,qv,ft,fqv,ux,uy,uz,fx,fy,fz) cc update clock (in minutes...) time = float(itime)*dt/60. CCcc output and plot: dtout=.2 ! in min dtape=1. ! in min if(amod(time+.1*dt/60.,dtout).lt.0.5*dt/60.) then c if(itime/20*20.eq.itime) then call diagno_1(ux,uy,uz,theta,nx,ny,nz,scr1,scr2,den) call diagno_2(qv,qc,qr,nx,ny,nz,scr1,scr2,den) call diagno_4(npx,xp,yp,zp,qcp,plic,radp) write(17) time,nx,nz write(17) ux,uy,uz,theta,qv,qc print*,'wrote bulk tape for t = ',time endif if(amod(time+.1*dt/60.,dtape).lt.0.5*dt/60.) then write(18) time,npx,xp,yp,zp,qcp,plic,radp print*,'wrote SD tape. npx,t = ',npx,time endif enddo ! TIME LOOP cc finished... stop end c spectrum: subroutine micro_set include 'aerosol.inc' c c MASS GRID (x) AND INTERFACE (o) STRUCTURE: c c 1 2 3 4 5 c o o o o o class boundaries c |--x--|--x--|--x--|--x--| classes c 1 2 3 4 c call chem_kin pi=4.*atan(1.) rho_a=rho_s open (10,file="new.CCN.data",status='old') read(10,*) do ic=1,nca read(10,*) ic1,rad0(ic),a1,a2,conc(ic) rad0(ic)=rad0(ic)*1.e-6 conc(ic)=conc(ic)*1.e6 mas0(ic)=4./3.*pi*rho_a*rad0(ic)**3 print 104,ic,rad0(ic)*1.e6,conc(ic)*1.e-6 104 format(2x,'++ bin,rada,conc: ',i4,2e16.6) enddo sum=0. do ic=1,nca sum=sum+conc(ic)*1.e-6 enddo print*,' ---- total conc (cm**-3): ',sum return end subroutine micro_set_old include 'aerosol.inc' c c MASS GRID (x) AND INTERFACE (o) STRUCTURE: c c 1 2 3 4 5 c o o o o o class boundaries c |--x--|--x--|--x--|--x--| classes c 1 2 3 4 c ccc lognormal aerosol distribution: flogn(r,rg,sigma,totc)=totc/(sqrt(2.*pi)*r*log(sigma)) 1 * exp(-log(r/rg)**2/(2.*log(sigma)**2) ) cc cc parameters for aerosol distribution as in MG2007: ccc one mode, as MG07: c totc1=100. *1.e6 ! total conc in m**-3 PRISTINE cC totc1=1000.*1.e6 ! total conc in m**-3 POLLUTED c rg1 =0.05 *1.e-6 ! mean radius in m c sigma1=2. ! std dev, no units ccc one mode - larger mode taken from VOCALS cc note totc2 replaced with totc1 below... c totc1=100. *1.e6 ! total conc in m**-3 c rg1 =0.075 *1.e-6 ! mean radius in m c sigma1=1.6 ! std dev, no units ccc two modes, as AR paper, table 5: ccc PRISTINE: c totc1=125. *1.e6 ! total conc in m**-3 c rg1 =0.011 *1.e-6 ! mean radius in m c sigma1=1.2 ! std dev, no units c totc2=65. *1.e6 ! total conc in m**-3 c rg2 =0.06 *1.e-6 ! mean radius in m c sigma2=1.7 ! std dev, no units ccc POLLUTED: c totc1=160. *1.e6 ! total conc in m**-3 c rg1 =0.029 *1.e-6 ! mean radius in m c sigma1=1.36 ! std dev, no units c totc2=380. *1.e6 ! total conc in m**-3 c rg2 =0.071 *1.e-6 ! mean radius in m c sigma2=1.57 ! std dev, no units cc corrected VOCALS: c totc1=60. *1.e6 ! total conc in m**-3 c rg1 =0.020 *1.e-6 ! mean radius in m c sigma1=1.4 ! std dev, no units c totc2=40. *1.e6 ! total conc in m**-3 c rg2 =0.075 *1.e-6 ! mean radius in m c sigma2=1.6 ! std dev, no units c modified VOCALS to have 200: totc1=120. *1.e6 ! total conc in m**-3 rg1 =0.020 *1.e-6 ! mean radius in m sigma1=1.4 ! std dev, no units totc2=80. *1.e6 ! total conc in m**-3 rg2 =0.075 *1.e-6 ! mean radius in m sigma2=1.6 ! std dev, no units call chem_kin pi=4.*atan(1.) rho_a=rho_s amn=2.e-9 amx=200.e-9 cc log: np=nca del=(alog10(amx)-alog10(amn))/(np-1) print*,del print*,'-- log:' do ic=1,np rad0(ic)=10.**(alog10(amn) + (ic-1)*del) print*,rad0(ic) mas0(ic)=4./3.*pi*rho_a*rad0(ic)**3 enddo do ic=2,np rads=(rad0(ic)+rad0(ic-1))/2. mas0s(ic)=4./3.*pi*rho_a*rads**3 enddo rads=rad0(1)-.5*(rad0(2)-rad0(1)) mas0s(1)=4./3.*pi*rho_a*rads**3 rads=rad0(nca)+.5*(rad0(nca)-rad0(nca-1)) mas0s(ncap)=4./3.*pi*rho_a*rads**3 print*,'mas0s: ',(mas0s(ic),ic=1,ncap) cc cccc AEROSOL CONCENTRATIONS IN BINS: mint=1000 ! subdivisions of each bin for integration of distribution do ic=1,nca rleft=(3.*mas0s(ic)/(4.*pi*rho_a))**(1./3.) rrigh=(3.*mas0s(ic+1)/(4.*pi*rho_a))**(1./3.) dra=(rrigh-rleft)/float(mint) sum=0. do iint=1,mint rrr=rleft + .5*dra + (iint-1)*dra cc 1 mode: c sum=sum + dra*flogn(rrr,rg1,sigma1,totc1) cc 2 modes: sum=sum + dra*flogn(rrr,rg1,sigma1,totc1) 1 + dra*flogn(rrr,rg2,sigma2,totc2) enddo conc(ic)=sum ccccccccccccc limit for class above..... c if(ic.ge.63) conc(ic)=1.e-18 ! 130 cccccccccccccccccccccccccccccccccccccccccccccccc print 104,ic,rad0(ic)*1.e6,conc(ic)*1.e-6 104 format(2x,'++ bin,rada,conc: ',i4,2e16.6) cc cloud water and drop radius within this class enddo sum=0. do ic=1,nca sum=sum+conc(ic)*1.e-6 enddo print*,' ---- total conc (cm**-3): ',sum cc adjustment: coe=200./sum sum=0. do ic=1,nca conc(ic)=conc(ic)*coe sum=sum+conc(ic)*1.e-6 enddo print*,' ---- adjusted total conc (cm**-3): ',sum cc end of adjustment do ic=1,nca print 105,ic,rad0(ic)*1.e6,conc(ic)*1.e-6 enddo 105 format(2x,i4,2e16.6,' ++ bin,rada,conc: ') return end subroutine chem_kin common /chemistry/ rho_w,amol_w,rho_s,amol_s,vanhof common /kinetic/ alpha_c,alpha_t,delta_v,delta_t cc water: amol_w=18.02 rho_w=1.e3 cc salt: c amonium sulfate: amol_s=132.14 ! amonium sulfate rho_s=1.77e3 vanhof=3. c sodium chloride c amol_s=56.44 ! sodium chloride c rho_s=2.17e3 c vanhof=2. cc kinetic coefficients: alpha_c=0.036 alpha_t=0.7 delta_v=2.16e-7 delta_t=1.04e-7 return end subroutine vel_check(ux,uy,uz,den,scr1,nx,ny,nz) dimension ux(nx,ny,nz),uy(nx,ny,nz),uz(nx,ny,nz), 1 den(nx,ny,nz),scr1(nx,ny,nz) call rhsdiv_2(ux,uy,uz,den,scr1,nx,ny,nz,1) divmx=-1.e15 divmn= 1.e15 do k=1,nz do j=1,ny do i=1,nx divmx=amax1(divmx,scr1(i,j,k)) divmn=amin1(divmn,scr1(i,j,k)) enddo enddo enddo print 251,divmn,divmx 251 format(1x,' --> min, max div from rhsdiv: ',2e12.4) return end subroutine noise(ff,nx,ny,nz,nzn) dimension ff(nx,ny,nz) double precision rand do i=1,nx do j=1,ny do k=1,nz ff(i,j,k)=0. enddo enddo enddo ccc different sequence... c do i=1,10000 do i=1,1990 aaa=2.*(rand()-.5) enddo do j=1,ny-1 do i=1,nx-1 do k=2,nzn ff(i,j,k)=2.*(rand()-.5) enddo enddo enddo do j=1,ny-1 do k=1,nzn ff(nx,j,k)=ff(1,j,k) enddo enddo do i=1,nx do k=1,nzn ff(i,ny,k)=ff(i,1,k) enddo enddo return end subroutine div_adv(ux,uy,uz,n1,n2,n3) dimension ux(n1+1,n2,n3),uy(n1,n2+1,n3),uz(n1,n2,n3+1) include 'param.grid' dimension div(nx,ny,nz) do i=1,nx do j=1,ny do k=1,nz div(i,j,k)= 1 ux(i+1,j,k)-ux(i,j,k) 2 +uy(i,j+1,k)-uy(i,j,k) 3 +uz(i,j,k+1)-uz(i,j,k) enddo enddo enddo call minmax(div,nxyz,amn,amx) print 201,amn,amx 201 format(1x,' --> min, max div: ',2e12.4) return end subroutine velprd_1(ux,uxp,uxa,uy,uyp,uya,uz,uzp,uza,nx,ny,nz,rho) dimension ux(nx,ny,nz),uy(nx,ny,nz),uz(nx,ny,nz) dimension uxp(nx,ny,nz),uyp(nx,ny,nz),uzp(nx,ny,nz) dimension uxa(nx+1,ny,nz),uya(nx,ny+1,nz),uza(nx,ny,nz+1), 1 rho(nx,ny,nz) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi do k=1,nz do j=1,ny do i=2,nx uxa(i,j,k) = . (0.75*(ux(i-1,j,k)*rho(i-1,j,k)+ux(i,j,k)*rho(i,j,k)) - . 0.25*(uxp(i-1,j,k)*rho(i-1,j,k)+uxp(i,j,k)*rho(i,j,k)))*dt/dx enddo cc cyclic in horizontal uxa(1,j,k) = uxa(nx,j,k) uxa(nx+1,j,k) = uxa(2,j,k) enddo enddo do k=1,nz do i=1,nx do j=2,ny uya(i,j,k) = . (0.75*(uy(i,j-1,k)*rho(i,j-1,k)+uy(i,j,k)*rho(i,j,k)) - . 0.25*(uyp(i,j-1,k)*rho(i,j-1,k)+uyp(i,j,k)*rho(i,j,k)))*dt/dy enddo cc cyclic in horizontal uya(i,1,k) = uya(i,ny,k) uya(i,ny+1,k) = uya(i,2,k) enddo enddo do i=1,nx do j=1,ny do k=2,nz uza(i,j,k) = . (0.75*(uz(i,j,k-1)*rho(i,j,k-1)+uz(i,j,k)*rho(i,j,k)) - . 0.25*(uzp(i,j,k-1)*rho(i,j,k-1)+uzp(i,j,k)*rho(i,j,k)) )*dt/dz enddo cc zero flux in vertical uza(i,j,1) = - uza(i,j,2) uza(i,j,nz+1) = - uza(i,j,nz) enddo enddo return end subroutine mpdata_3d(u1,u2,u3,x,h,iflg,liner) include 'param.grid' parameter(n=nx,m=ny,l=nz) parameter(n1=n+1,n2=m+1,n3=l+1) parameter(n1m=n1-1,n2m=n2-1,n3m=n3-1) dimension u1(n1,n2m,n3m),u2(n1m,n2,n3m),u3(n1m,n2m,n3), * x(n1m,n2m,n3m),h(n1m,n2m,n3m) dimension v1(n1,n2m,n3m),v2(n1m,n2,n3m),v3(n1m,n2m,n3), * f1(n1,n2m,n3m),f2(n1m,n2,n3m),f3(n1m,n2m,n3), * cp(n1m,n2m,n3m),cn(n1m,n2m,n3m),mx(n1m,n2m,n3m), * mn(n1m,n2m,n3m),bcx(n2m,n3m,2),bcy(n1m,n3m,2) real mx,mn parameter(iord0=2,isor=1,nonos0=1,idiv=0) data ep/1.e-10/ c pp(y)= amax1(0.,y) pn(y)=-amin1(0.,y) donor(y1,y2,a)=pp(a)*y1-pn(a)*y2 vdyf(x1,x2,a,r)=(abs(a)-a**2/r)*(abs(x2)-abs(x1)) 1 /(abs(x2)+abs(x1)+ep) vcorr(a,b,y1,y2,r)=-0.125*a*b*y1/(y2*r) vcor31(a,x0,x1,x2,x3,r)= -(a -3.*abs(a)*a/r+2.*a**3/r**2)/3. 1 *(abs(x0)+abs(x3)-abs(x1)-abs(x2)) 2 /(abs(x0)+abs(x3)+abs(x1)+abs(x2)+ep) vcor32(a,b,y1,y2,r)=0.25*b/r*(abs(a)-2.*a**2/r)*y1/y2 vcor33(a,b,c,y1,y2,r)=-a*b*c/(24.**r**2)*y1/y2 vdiv1(a1,a2,a3,r)=0.25*a2*(a3-a1)/r vdiv2(a,b1,b2,b3,b4,r)=0.25*a*(b1+b2-b3-b4)/r c ibcx=1 ibcy=1 iord=iord0 if(liner.eq.1) iord=1 nonos=nonos0 if(liner.eq.1) nonos=1 if(isor.eq.3) iord=max0(iord,3) do j=1,n2m do i=1,n1m do k=1,n3 v3(i,j,k) = u3(i,j,k) enddo end do end do do k=1,n3m do j=1,n2m do i=1,n1 v1(i,j,k) = u1(i,j,k) end do end do do i=1,n1m do j=1,n2 v2(i,j,k) = u2(i,j,k) end do end do end do if (nonos.eq.1) then do k=2,n3m-1 do j=2,n2m-1 do i=2,n1m-1 mx(i,j,k)=amax1(x(i-1,j,k),x(i,j,k),x(i+1,j,k), . x(i,j-1,k),x(i,j+1,k),x(i,j,k+1),x(i,j,k-1)) mn(i,j,k)=amin1(x(i-1,j,k),x(i,j,k),x(i+1,j,k), . x(i,j-1,k),x(i,j+1,k),x(i,j,k+1),x(i,j,k-1)) end do end do end do endif do 30 itr=1,iord do 331 k=1,n3m do 331 j=1,n2m do 331 i=2,n1m 331 f1(i,j,k)=donor(x(i-1,j,k),x(i,j,k),v1(i,j,k)) do k=1,n3m do j=1,n2m f1(1 ,j,k)= f1(n1m,j,k) f1(n1,j,k)= f1(2,j,k) enddo enddo do 332 k=1,n3m do 332 j=2,n2m do 332 i=1,n1m 332 f2(i,j,k)=donor(x(i,j-1,k),x(i,j,k),v2(i,j,k)) do k=1,n3m do i=1,n1m f2(i,1 ,k)= f2(i,n2m,k) f2(i,n2,k)= f2(i,2,k) enddo enddo do 333 k=2,n3m do 333 j=1,n2m do 333 i=1,n1m 333 f3(i,j,k)=donor(x(i,j,k-1),x(i,j,k),v3(i,j,k)) if (iflg.eq.7) then do j=1,n2m do i=1,n1m f3(i,j, 1)=donor(x(i,j,1),x(i,j,1),v3(i,j,1)) f3(i,j,n3)=donor(x(i,j,n3m),x(i,j,n3m),v3(i,j,n3)) end do end do else do j=1,n2m do i=1,n1m f3(i,j, 1)=-f3(i,j, 2) f3(i,j,n3)=-f3(i,j,n3m) enddo enddo end if do 334 k=1,n3m do 334 j=1,n2m do 334 i=1,n1m 334 x(i,j,k)=x(i,j,k)-( f1(i+1,j,k)-f1(i,j,k) . +f2(i,j+1,k)-f2(i,j,k) . +f3(i,j,k+1)-f3(i,j,k) )/h(i,j,k) if(itr.eq.iord) go to 6 do 48 k=1,n3m do 48 j=1,n2m do 48 i=1,n1 f1(i,j,k)=v1(i,j,k) 48 v1(i,j,k)=0. do 49 k=1,n3m do 49 j=1,n2 do 49 i=1,n1m f2(i,j,k)=v2(i,j,k) 49 v2(i,j,k)=0. do 50 k=1,n3 do 50 j=1,n2m do 50 i=1,n1m f3(i,j,k)=v3(i,j,k) 50 v3(i,j,k)=0. compute antidiffusive velocities in x direction do 51 k=2,n3-2 do 51 j=2,n2-2 do 51 i=2,n1m 51 v1(i,j,k)=vdyf(x(i-1,j,k),x(i,j,k),f1(i,j,k), * .5*(h(i-1,j,k)+h(i,j,k))) * +vcorr(f1(i,j,k), * f2(i-1,j,k)+f2(i-1,j+1,k)+f2(i,j+1,k)+f2(i,j,k), * abs(x(i-1,j+1,k))+abs(x(i,j+1,k)) * -abs(x(i-1,j-1,k))-abs(x(i,j-1,k)), * abs(x(i-1,j+1,k))+abs(x(i,j+1,k)) * +abs(x(i-1,j-1,k))+abs(x(i,j-1,k))+ep, * .5*(h(i-1,j,k)+h(i,j,k))) * +vcorr(f1(i,j,k), * f3(i-1,j,k)+f3(i-1,j,k+1)+f3(i,j,k+1)+f3(i,j,k), * abs(x(i-1,j,k+1))+abs(x(i,j,k+1)) * -abs(x(i-1,j,k-1))-abs(x(i,j,k-1)), * abs(x(i-1,j,k+1))+abs(x(i,j,k+1)) * +abs(x(i-1,j,k-1))+abs(x(i,j,k-1))+ep, * .5*(h(i-1,j,k)+h(i,j,k))) if(ibcy.eq.1) then do k=2,n3-2 do i=2,n1m v1(i,1,k)=vdyf(x(i-1,1,k),x(i,1,k),f1(i,1,k), * .5*(h(i-1,1,k)+h(i,1,k))) * +vcorr(f1(i,1,k), * f2(i-1,1,k)+f2(i-1,2,k)+f2(i,2,k)+f2(i,1,k), * abs(x(i-1,2,k))+abs(x(i,2,k)) * -abs(x(i-1,n2-2,k))-abs(x(i,n2-2,k)), * abs(x(i-1,2,k))+abs(x(i,2,k)) * +abs(x(i-1,n2-2,k))+abs(x(i,n2-2,k))+ep, * .5*(h(i-1,1,k)+h(i,1,k))) * +vcorr(f1(i,1,k), * f3(i-1,1,k)+f3(i-1,1,k+1)+f3(i,1,k+1)+f3(i,1,k), * abs(x(i-1,1,k+1))+abs(x(i,1,k+1)) * -abs(x(i-1,1,k-1))-abs(x(i,1,k-1)), * abs(x(i-1,1,k+1))+abs(x(i,1,k+1)) * +abs(x(i-1,1,k-1))+abs(x(i,1,k-1))+ep, * .5*(h(i-1,1,k)+h(i,1,k))) v1(i,n2-1,k)=v1(i,1,k) enddo enddo endif if(idiv.eq.1) then do 511 k=2,n3-2 do 511 j=2-ibcy,n2-2+ibcy do 511 i=2,n1m 511 v1(i,j,k)=v1(i,j,k) * -vdiv1(f1(i-1,j,k),f1(i,j,k),f1(i+1,j,k), * .5*(h(i-1,j,k)+h(i,j,k))) * -vdiv2(f1(i,j,k),f2(i-1,j+1,k),f2(i,j+1,k),f2(i-1,j,k), * f2(i,j,k), .5*(h(i-1,j,k)+h(i,j,k))) * -vdiv2(f1(i,j,k),f3(i-1,j,k+1),f3(i,j,k+1),f3(i-1,j,k), * f3(i,j,k), .5*(h(i-1,j,k)+h(i,j,k))) endif if(isor.eq.3) then do 61 k=2,n3-2 do 61 j=2-ibcy,n2-2+ibcy do 61 i=3,n1-2 61 v1(i,j,k)=v1(i,j,k) +vcor31(f1(i,j,k), 1 x(i-2,j,k),x(i-1,j,k),x(i,j,k),x(i+1,j,k), 1 .5*(h(i-1,j,k)+h(i,j,k))) if(ibcx.eq.1) then do k=2,n3-2 do j=2,n2-2 v1(2,j,k)=v1(2,j,k) +vcor31(f1(2,j,k), 1 x(n1-2,j,k),x(1,j,k),x(2,j,k),x(3,j,k), 1 .5*(h(1,j,k)+h(2,j,k))) v1(n1-1,j,k)=v1(n1-1,j,k) +vcor31(f1(n1-1,j,k),x(n1-3,j,k), 1 x(n1-2,j,k),x(n1-1,j,k),x(2,j,k), 1 .5*(h(n1-2,j,k)+h(n1-1,j,k))) enddo enddo endif do 62 k=2,n3-2 do 62 j=2,n2-2 do 62 i=3-ibcx,n1-2+ibcx 62 v1(i,j,k)=v1(i,j,k) 1 +vcor32(f1(i,j,k), * f2(i-1,j,k)+f2(i-1,j+1,k)+f2(i,j+1,k)+f2(i,j,k), * abs(x(i,j+1,k))-abs(x(i,j-1,k))-abs(x(i-1,j+1,k)) * +abs(x(i-1,j-1,k)), abs(x(i,j+1,k))+abs(x(i,j-1,k)) * +abs(x(i-1,j+1,k))+abs(x(i-1,j-1,k))+ep, * .5*(h(i-1,j,k)+h(i,j,k))) 1 +vcor32(f1(i,j,k), * f3(i-1,j,k)+f3(i-1,j,k+1)+f3(i,j,k+1)+f3(i,j,k), * abs(x(i,j,k+1))-abs(x(i,j,k-1))-abs(x(i-1,j,k+1)) * +abs(x(i-1,j,k-1)), abs(x(i,j,k+1))+abs(x(i,j,k-1)) * +abs(x(i-1,j,k+1))+abs(x(i-1,j,k-1))+ep, * .5*(h(i-1,j,k)+h(i,j,k))) do 63 k=2,n3-2 do 63 j=2,n2-2 do 63 i=3-ibcx,n1-2+ibcx 63 v1(i,j,k)=v1(i,j,k) + vcor33(f1(i,j,k), * f2(i-1,j,k)+f2(i-1,j+1,k)+f2(i,j+1,k)+f2(i,j,k), * f3(i-1,j,k)+f3(i-1,j,k+1)+f3(i,j,k+1)+f3(i,j,k), * abs(x(i-1,j+1,k+1))-abs(x(i-1,j-1,k+1)) * -abs(x(i-1,j+1,k-1))+abs(x(i-1,j-1,k-1)) * +abs(x(i,j+1,k+1))-abs(x(i,j-1,k+1)) * -abs(x(i,j+1,k-1))+abs(x(i,j-1,k-1)), * abs(x(i-1,j+1,k+1))+abs(x(i-1,j-1,k+1)) * +abs(x(i-1,j+1,k-1))+abs(x(i-1,j-1,k-1)) * +abs(x(i,j+1,k+1))+abs(x(i,j-1,k+1)) * +abs(x(i,j+1,k-1))+abs(x(i,j-1,k-1))+ep, * .5*(h(i-1,j,k)+h(i,j,k))) endif compute antidiffusive velocities in y direction do 52 k=2,n3-2 do 52 j=2,n2m do 52 i=2,n1-2 52 v2(i,j,k)=vdyf(x(i,j-1,k),x(i,j,k),f2(i,j,k), * .5*(h(i,j-1,k)+h(i,j,k))) * +vcorr(f2(i,j,k), * f1(i,j-1,k)+f1(i,j,k)+f1(i+1,j,k)+f1(i+1,j-1,k), * abs(x(i+1,j-1,k))+abs(x(i+1,j,k)) * -abs(x(i-1,j-1,k))-abs(x(i-1,j,k)), * abs(x(i+1,j-1,k))+abs(x(i+1,j,k)) * +abs(x(i-1,j-1,k))+abs(x(i-1,j,k))+ep, * .5*(h(i,j-1,k)+h(i,j,k))) * +vcorr(f2(i,j,k), * f3(i,j-1,k)+f3(i,j,k)+f3(i,j,k+1)+f3(i,j-1,k+1), * abs(x(i,j-1,k+1))+abs(x(i,j,k+1)) * -abs(x(i,j-1,k-1))-abs(x(i,j,k-1)), * abs(x(i,j-1,k+1))+abs(x(i,j,k+1)) * +abs(x(i,j-1,k-1))+abs(x(i,j,k-1))+ep, * .5*(h(i,j-1,k)+h(i,j,k))) if(ibcx.eq.1) then do k=2,n3-2 do j=2,n2m v2(1,j,k)=vdyf(x(1,j-1,k),x(1,j,k),f2(1,j,k), * .5*(h(1,j-1,k)+h(1,j,k))) * +vcorr(f2(1,j,k), * f1(1,j-1,k)+f1(1,j,k)+f1(2,j,k)+f1(2,j-1,k), * abs(x(2,j-1,k))+abs(x(2,j,k)) * -abs(x(n1-2,j-1,k))-abs(x(n1-2,j,k)), * abs(x(2,j-1,k))+abs(x(2,j,k)) * +abs(x(n1-2,j-1,k))+abs(x(n1-2,j,k))+ep, * .5*(h(1,j-1,k)+h(1,j,k))) * +vcorr(f2(1,j,k), * f3(1,j-1,k)+f3(1,j,k)+f3(1,j,k+1)+f3(1,j-1,k+1), * abs(x(1,j-1,k+1))+abs(x(1,j,k+1)) * -abs(x(1,j-1,k-1))-abs(x(1,j,k-1)), * abs(x(1,j-1,k+1))+abs(x(1,j,k+1)) * +abs(x(1,j-1,k-1))+abs(x(1,j,k-1))+ep, * .5*(h(1,j-1,k)+h(1,j,k))) v2(n1-1,j,k)=v2(1,j,k) enddo enddo endif if(idiv.eq.1) then do 521 k=2,n3-2 do 521 j=2,n2m do 521 i=2-ibcx,n1-2+ibcx 521 v2(i,j,k)=v2(i,j,k) * -vdiv1(f2(i,j-1,k),f2(i,j,k),f2(i,j+1,k), * .5*(h(i,j-1,k)+h(i,j,k))) * -vdiv2(f2(i,j,k),f1(i+1,j-1,k),f1(i+1,j,k),f1(i,j-1,k), * f1(i,j,k), .5*(h(i,j-1,k)+h(i,j,k))) * -vdiv2(f2(i,j,k),f3(i,j-1,k+1),f3(i,j,k+1),f3(i,j-1,k), * f3(i,j,k), .5*(h(i,j-1,k)+h(i,j,k))) endif if(isor.eq.3) then do 71 k=2,n3-2 do 71 j=3,n2-2 do 71 i=2-ibcx,n1-2+ibcx 71 v2(i,j,k)=v2(i,j,k) +vcor31(f2(i,j,k), 1 x(i,j-2,k),x(i,j-1,k),x(i,j,k),x(i,j+1,k), 1 .5*(h(i,j-1,k)+h(i,j,k))) if(ibcy.eq.1) then do k=2,n3-2 do i=2,n1-2 v2(i,2,k)=v2(i,2,k) +vcor31(f2(i,2,k), 1 x(i,n2-2,k),x(i,1,k),x(i,2,k),x(i,3,k), 1 .5*(h(i,1,k)+h(i,2,k))) v2(i,n2-1,k)=v2(i,n2-1,k) +vcor31(f2(i,n2-1,k),x(i,n2-3,k), 1 x(i,n2-2,k),x(i,n2-1,k),x(i,2,k), 1 .5*(h(i,n2-2,k)+h(i,n2-1,k))) enddo enddo endif do 72 k=2,n3-2 do 72 j=3-ibcy,n2-2+ibcy do 72 i=2,n1-2 72 v2(i,j,k)=v2(i,j,k) 1 +vcor32(f2(i,j,k), * f1(i,j-1,k)+f1(i+1,j-1,k)+f1(i+1,j,k)+f1(i,j,k), * abs(x(i+1,j,k))-abs(x(i-1,j,k))-abs(x(i+1,j-1,k)) * +abs(x(i-1,j-1,k)), abs(x(i+1,j,k))+abs(x(i-1,j,k)) * +abs(x(i+1,j-1,k))+abs(x(i-1,j-1,k))+ep, * .5*(h(i,j-1,k)+h(i,j,k))) 1 +vcor32(f2(i,j,k), * f3(i,j-1,k)+f3(i,j-1,k+1)+f3(i,j,k+1)+f3(i,j,k), * abs(x(i,j,k+1))-abs(x(i,j,k-1))-abs(x(i,j-1,k+1)) * +abs(x(i,j-1,k-1)), abs(x(i,j,k+1))+abs(x(i,j,k-1)) * +abs(x(i,j-1,k+1))+abs(x(i,j-1,k-1))+ep, * .5*(h(i,j-1,k)+h(i,j,k))) do 73 k=2,n3-2 do 73 j=3-ibcy,n2-2+ibcy do 73 i=2,n1-2 73 v2(i,j,k)=v2(i,j,k) + vcor33(f2(i,j,k), * f1(i,j-1,k)+f1(i+1,j-1,k)+f1(i+1,j,k)+f1(i,j,k), * f3(i,j-1,k)+f3(i,j-1,k+1)+f3(i,j,k+1)+f3(i,j,k), * abs(x(i+1,j-1,k+1))-abs(x(i-1,j-1,k+1)) * -abs(x(i+1,j-1,k-1))+abs(x(i-1,j-1,k-1)) * +abs(x(i+1,j,k+1))-abs(x(i-1,j,k+1)) * -abs(x(i+1,j,k-1))+abs(x(i-1,j,k-1)), * abs(x(i+1,j-1,k+1))+abs(x(i-1,j-1,k+1)) * +abs(x(i+1,j-1,k-1))+abs(x(i-1,j-1,k-1)) * +abs(x(i+1,j,k+1))+abs(x(i-1,j,k+1)) * +abs(x(i+1,j,k-1))+abs(x(i-1,j,k-1))+ep, * .5*(h(i,j-1,k)+h(i,j,k))) endif compute antidiffusive velocities in z direction do 53 k=2,n3m do 53 j=2,n2-2 do 53 i=2,n1-2 53 v3(i,j,k)=vdyf(x(i,j,k-1),x(i,j,k),f3(i,j,k), * .5*(h(i,j,k-1)+h(i,j,k))) * +vcorr(f3(i,j,k), * f1(i,j,k-1)+f1(i,j,k)+f1(i+1,j,k)+f1(i+1,j,k-1), * abs(x(i+1,j,k-1))+abs(x(i+1,j,k)) * -abs(x(i-1,j,k-1))-abs(x(i-1,j,k)), * abs(x(i+1,j,k-1))+abs(x(i+1,j,k)) * +abs(x(i-1,j,k-1))+abs(x(i-1,j,k))+ep, * .5*(h(i,j,k-1)+h(i,j,k))) * +vcorr(f3(i,j,k), * f2(i,j,k-1)+f2(i,j+1,k-1)+f2(i,j+1,k)+f2(i,j,k), * abs(x(i,j+1,k-1))+abs(x(i,j+1,k)) * -abs(x(i,j-1,k-1))-abs(x(i,j-1,k)), * abs(x(i,j+1,k-1))+abs(x(i,j+1,k)) * +abs(x(i,j-1,k-1))+abs(x(i,j-1,k))+ep, * .5*(h(i,j,k-1)+h(i,j,k))) if(ibcx.eq.1) then do k=2,n3m do j=2,n2-2 v3(1,j,k)=vdyf(x(1,j,k-1),x(1,j,k),f3(1,j,k), * .5*(h(1,j,k-1)+h(1,j,k))) * +vcorr(f3(1,j,k), * f1(1,j,k-1)+f1(1,j,k)+f1(2,j,k)+f1(2,j,k-1), * abs(x(2,j,k-1))+abs(x(2,j,k)) * -abs(x(n1-2,j,k-1))-abs(x(n1-2,j,k)), * abs(x(2,j,k-1))+abs(x(2,j,k)) * +abs(x(n1-2,j,k-1))+abs(x(n1-2,j,k))+ep, * .5*(h(1,j,k-1)+h(1,j,k))) * +vcorr(f3(1,j,k), * f2(1,j,k-1)+f2(1,j+1,k-1)+f2(1,j+1,k)+f2(1,j,k), * abs(x(1,j+1,k-1))+abs(x(1,j+1,k)) * -abs(x(1,j-1,k-1))-abs(x(1,j-1,k)), * abs(x(1,j+1,k-1))+abs(x(1,j+1,k)) * +abs(x(1,j-1,k-1))+abs(x(1,j-1,k))+ep, * .5*(h(1,j,k-1)+h(1,j,k))) v3(n1-1,j,k)=v3(1,j,k) enddo enddo endif if(ibcy.eq.1) then do k=2,n3m do i=2,n1-2 v3(i,1,k)=vdyf(x(i,1,k-1),x(i,1,k),f3(i,1,k), * .5*(h(i,1,k-1)+h(i,1,k))) * +vcorr(f3(i,1,k), * f1(i,1,k-1)+f1(i,1,k)+f1(i+1,1,k)+f1(i+1,1,k-1), * abs(x(i+1,1,k-1))+abs(x(i+1,1,k)) * -abs(x(i-1,1,k-1))-abs(x(i-1,1,k)), * abs(x(i+1,1,k-1))+abs(x(i+1,1,k)) * +abs(x(i-1,1,k-1))+abs(x(i-1,1,k))+ep, * .5*(h(i,1,k-1)+h(i,1,k))) * +vcorr(f3(i,1,k), * f2(i,1,k-1)+f2(i,2,k-1)+f2(i,2,k)+f2(i,1,k), * abs(x(i,2,k-1))+abs(x(i,2,k)) * -abs(x(i,n2-2,k-1))-abs(x(i,n2-2,k)), * abs(x(i,2,k-1))+abs(x(i,2,k)) * +abs(x(i,n2-2,k-1))+abs(x(i,n2-2,k))+ep, * .5*(h(i,1,k-1)+h(i,1,k))) v3(i,n2-1,k)=v3(i,1,k) enddo enddo if(ibcx.eq.1) then do k=2,n3m v3(1,1,k)=vdyf(x(1,1,k-1),x(1,1,k),f3(1,1,k), * .5*(h(1,1,k-1)+h(1,1,k))) * +vcorr(f3(1,1,k), * f1(1,1,k-1)+f1(1,1,k)+f1(2,1,k)+f1(2,1,k-1), * abs(x(2,1,k-1))+abs(x(2,1,k)) * -abs(x(n1-2,1,k-1))-abs(x(n1-2,1,k)), * abs(x(2,1,k-1))+abs(x(2,1,k)) * +abs(x(n1-2,1,k-1))+abs(x(n1-2,1,k))+ep, * .5*(h(1,1,k-1)+h(1,1,k))) * +vcorr(f3(1,1,k), * f2(1,1,k-1)+f2(1,2,k-1)+f2(1,2,k)+f2(1,1,k), * abs(x(1,2,k-1))+abs(x(1,2,k)) * -abs(x(1,n2-2,k-1))-abs(x(1,n2-2,k)), * abs(x(1,2,k-1))+abs(x(1,2,k)) * +abs(x(1,n2-2,k-1))+abs(x(1,n2-2,k))+ep, * .5*(h(1,1,k-1)+h(1,1,k))) v3(n1-1,n2-1,k)=v3(1,1,k) v3(n1-1,1,k)=v3(1,1,k) v3(1,n2-1,k)=v3(1,1,k) enddo endif endif if(idiv.eq.1) then do 531 k=2,n3m do 531 j=2-ibcy,n2-2+ibcy do 531 i=2-ibcx,n1-2+ibcx 531 v3(i,j,k)=v3(i,j,k) * -vdiv1(f3(i,j,k-1),f3(i,j,k),f3(i,j,k+1), * .5*(h(i,j,k-1)+h(i,j,k))) * -vdiv2(f3(i,j,k),f1(i+1,j,k-1),f1(i+1,j,k),f1(i,j,k-1), * f1(i,j,k), .5*(h(i,j,k-1)+h(i,j,k))) * -vdiv2(f3(i,j,k),f2(i,j+1,k-1),f2(i,j+1,k),f2(i,j,k-1), * f2(i,j,k), .5*(h(i,j,k-1)+h(i,j,k))) endif if(isor.eq.3) then do 81 k=3,n3-2 do 81 j=2-ibcy,n2-2+ibcy do 81 i=2-ibcx,n1-2+ibcx 81 v3(i,j,k)=v3(i,j,k) +vcor31(f3(i,j,k), 1 x(i,j,k-2),x(i,j,k-1),x(i,j,k),x(i,j,k+1), 1 .5*(h(i,j,k-1)+h(i,j,k))) do 82 k=3,n3-2 do 82 j=2,n2-2 do 82 i=2,n1-2 82 v3(i,j,k)=v3(i,j,k) 1 +vcor32(f3(i,j,k), * f2(i,j,k-1)+f2(i,j+1,k-1)+f2(i,j+1,k)+f2(i,j,k), * abs(x(i,j+1,k))-abs(x(i,j-1,k))-abs(x(i,j+1,k-1)) * +abs(x(i,j-1,k-1)), abs(x(i,j+1,k))+abs(x(i,j-1,k)) * +abs(x(i,j+1,k-1))+abs(x(i,j-1,k-1))+ep, * .5*(h(i,j,k-1)+h(i,j,k))) 1 +vcor32(f3(i,j,k), * f1(i,j,k-1)+f1(i+1,j,k-1)+f1(i+1,j,k)+f1(i,j,k), * abs(x(i+1,j,k))-abs(x(i-1,j,k))-abs(x(i+1,j,k-1)) * +abs(x(i-1,j,k-1)), abs(x(i+1,j,k))+abs(x(i-1,j,k)) * +abs(x(i+1,j,k-1))+abs(x(i-1,j,k-1))+ep, * .5*(h(i,j,k-1)+h(i,j,k))) do 83 k=2,n3-2 do 83 j=2,n2-2 do 83 i=2,n1-2 83 v3(i,j,k)=v3(i,j,k) + vcor33(f3(i,j,k), * f2(i,j,k-1)+f2(i,j+1,k-1)+f2(i,j+1,k)+f2(i,j,k), * f1(i,j,k-1)+f1(i+1,j,k-1)+f1(i+1,j,k)+f1(i,j,k), * abs(x(i+1,j+1,k-1))-abs(x(i+1,j-1,k-1)) * -abs(x(i-1,j+1,k-1))+abs(x(i-1,j-1,k-1)) * +abs(x(i+1,j+1,k))-abs(x(i+1,j-1,k)) * -abs(x(i-1,j+1,k))+abs(x(i-1,j-1,k)), * abs(x(i+1,j+1,k-1))+abs(x(i+1,j-1,k-1)) * +abs(x(i-1,j+1,k-1))+abs(x(i-1,j-1,k-1)) * +abs(x(i+1,j+1,k))+abs(x(i+1,j-1,k)) * +abs(x(i-1,j+1,k))+abs(x(i-1,j-1,k))+ep, * .5*(h(i,j,k-1)+h(i,j,k))) endif if(nonos.eq.1) then c non-osscilatory option do 401 k=2,n3m-1 do 401 j=2,n2m-1 do 401 i=2,n1m-1 mx(i,j,k)=amax1(x(i-1,j,k),x(i,j,k),x(i+1,j,k),mx(i,j,k), . x(i,j-1,k),x(i,j+1,k),x(i,j,k+1),x(i,j,k-1)) 401 mn(i,j,k)=amin1(x(i-1,j,k),x(i,j,k),x(i+1,j,k),mn(i,j,k), . x(i,j-1,k),x(i,j+1,k),x(i,j,k+1),x(i,j,k-1)) do 402 k=2,n3m-1 do 402 j=2,n2m-1 do 402 i=2,n1m 402 f1(i,j,k)=donor(x(i-1,j,k),x(i,j,k),v1(i,j,k)) do 403 k=2,n3m-1 do 403 j=2,n2m do 403 i=2,n1m-1 403 f2(i,j,k)=donor(x(i,j-1,k),x(i,j,k),v2(i,j,k)) do 4033 k=2,n3m do 4033 j=2,n2m-1 do 4033 i=2,n1m-1 4033 f3(i,j,k)=donor(x(i,j,k-1),x(i,j,k),v3(i,j,k)) do 404 k=2,n3m-1 do 404 j=2,n2m-1 do 404 i=2,n1m-1 cp(i,j,k)=(mx(i,j,k)-x(i,j,k))*h(i,j,k)/ 1( pn(f1(i+1,j,k))+pp(f1(i,j,k)) 2 +pn(f2(i,j+1,k))+pp(f2(i,j,k)) 3 +pn(f3(i,j,k+1))+pp(f3(i,j,k))+ep) cn(i,j,k)=(x(i,j,k)-mn(i,j,k))*h(i,j,k)/ 1( pp(f1(i+1,j,k))+pn(f1(i,j,k)) 2 +pp(f2(i,j+1,k))+pn(f2(i,j,k)) 3 +pp(f3(i,j,k+1))+pn(f3(i,j,k))+ep) 404 continue do 405 k=3,n3m-1 do 405 j=3,n2m-1 do 405 i=3,n1m-1 v1(i,j,k)=pp(v1(i,j,k))* 1 ( amin1(1.,cp(i,j,k),cn(i-1,j,k))*pp(sign(1., x(i-1,j,k))) 1 +amin1(1.,cp(i-1,j,k),cn(i,j,k))*pp(sign(1.,-x(i-1,j,k))) ) 2 -pn(v1(i,j,k))* 2 ( amin1(1.,cp(i-1,j,k),cn(i,j,k))*pp(sign(1., x(i ,j,k ))) 2 +amin1(1.,cp(i,j,k),cn(i-1,j,k))*pp(sign(1.,-x(i ,j,k ))) ) v2(i,j,k)=pp(v2(i,j,k))* 1 ( amin1(1.,cp(i,j,k),cn(i,j-1,k))*pp(sign(1., x(i,j-1,k))) 1 +amin1(1.,cp(i,j-1,k),cn(i,j,k))*pp(sign(1.,-x(i,j-1,k))) ) 1 -pn(v2(i,j,k))* 2 ( amin1(1.,cp(i,j-1,k),cn(i,j,k))*pp(sign(1., x(i ,j,k ))) 2 +amin1(1.,cp(i,j,k),cn(i,j-1,k))*pp(sign(1.,-x(i ,j,k ))) ) 405 v3(i,j,k)=pp(v3(i,j,k))* 1 ( amin1(1.,cp(i,j,k),cn(i,j,k-1))*pp(sign(1., x(i,j,k-1))) 1 +amin1(1.,cp(i,j,k-1),cn(i,j,k))*pp(sign(1.,-x(i,j,k-1))) ) 1 -pn(v3(i,j,k))* 2 ( amin1(1.,cp(i,j,k-1),cn(i,j,k))*pp(sign(1., x(i ,j,k ))) 2 +amin1(1.,cp(i,j,k),cn(i,j,k-1))*pp(sign(1.,-x(i ,j,k ))) ) endif if (ibcx.eq.1) then do k=1,n3m do j=1,n2m v1(1 ,j,k)=v1(n1m,j,k) v1(n1,j,k)=v1(2 ,j,k) end do end do end if if (ibcy.eq.1) then do k=1,n3m do i=1,n1m v2(i, 1,k)=v2(i,n2m,k) v2(i,n2,k)=v2(i,2 ,k) end do end do end if 30 continue 6 continue return end c I tried to convert back pressure solver from 2d babyE into 3D c problem in precon - different codes in babyE and 3D from Piotr... c TESTS WITHOUT PRECOND IN BABY EULAG: OK without it.... c subroutine gcrk_2(p,pfx,pfy,pfz,u,v,w,n1,n2,n3,itr,eps0) real p(*),pfx(*),pfy(*),pfz(*),u(*),v(*),w(*) include 'param.grid' parameter(n=nx,m=ny,l=nz) parameter(nn=nxyz,nl=nxyz) common// r(nn),qr(nn),ar(nn) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi common /prof_d/ rho0(nz),th0(nz),th_e(nz),ux_e(nz),uy_e(nz) common /strtch/ height(nz),gac(nz) common/itero/ niter,nitsm,icount,eer,eem parameter (lord=3) dimension x(nn,lord),ax(nn,lord),ax2(lord),axar(lord),del(lord) dimension rho2d(nx,ny,nz) cconvergence test modes ************************************************** logical ctest * data ctest/.false./ * cc data ctest/.true./ * c parameter (nplt=100) * c dimension err(0:nplt),xitr(0:nplt) * c if(ctest) then * c itr=6000/lord * c ner=60 * c snorm=1./float(n*l) * c eps0=1.e-15 * c endif * cconvergence test modes ************************************************** do k=1,l do j=1,m do i=1,n rho2d(i,j,k)=rho0(k)*gac(k) enddo enddo enddo eps=eps0*dti epa=1.e-30 nlc=0 itmn=5 cc iprc 0-no preconditioning, 1-with precon iprc=0 if(iprc.eq.1) STOP 'preconditioner not reday...' call precon_1(r,qr,ar,pfx,pfy,pfz,rho2d,iprc,0) do k=1,nl r(k)=0. ar(k)=0. qr(k)=0. enddo do i=1,lord do k=1,nl x(k,i)=0. ax(k,i)=0. enddo enddo call prforc_2(p,pfx,pfy,pfz,u,v,w,n1,n2,n3) call rhsdiv_2(pfx,pfy,pfz,rho2d,r,n1,n2,n3,-1) call precon_1(r,qr,ar,pfx,pfy,pfz,rho2d,iprc,1) eer0=0. eem0=-1.e15 rl20=0. do k=1,nl eer0=eer0+qr(k)**2 eem0=amax1(eem0,abs(qr(k))) rl20=rl20+r(k)**2 enddo eer0=amax1(eer0,epa) eem0=amax1(eem0,epa) rl20=amax1(rl20,epa) cconvergence test modes ************************************************** c if(ctest) then * c do ier=0,nplt * c err(ier)=eps * c enddo * c eer=-1.e15 * c do 3 k=1,nl * c 3 eer=amax1(eer,abs(r(k))) * c err(0)=eer * c print 300, err(0) * c 300 format(4x,e11.4,' residual error at it=1') * c endif * cconvergence test modes ************************************************** do k=1,nl x(k,1)=qr(k) enddo call laplc_2(x(1,1),ax(1,1),pfx,pfy,pfz,n1,n2,n3) do 100 it=1,itr do i=1,lord ax2(i)=0. rax=0. do k=1,nl rax=rax+r(k)*ax(k,i) ax2(i)=ax2(i)+ax(k,i)*ax(k,i) enddo ax2(i)=amax1(epa,ax2(i)) beta=-rax/ax2(i) dvmx=-1.e15 rl2=0. do k=1,nl p(k)=p(k)+beta* x(k,i) r(k)=r(k)+beta*ax(k,i) dvmx=amax1(dvmx,abs(r(k))) rl2=rl2+r(k)*r(k) enddo if(dvmx.le.eps.and.it.ge.itmn) go to 200 if(rl2.ge.rl20.and..not.ctest) go to 200 rl20=amax1(rl2,epa) call precon_1(r,qr,ar,pfx,pfy,pfz,rho2d,iprc,1) call laplc_2(qr,ar,pfx,pfy,pfz,n1,n2,n3) nlc=nlc+1 do ii=1,i axar(ii)=0. do k=1,nl axar(ii)=axar(ii)+ax(k,ii)*ar(k) enddo del(ii)=-axar(ii)/ax2(ii) enddo if(i.lt.lord) then do k=1,nl x(k,i+1)=qr(k) ax(k,i+1)=ar(k) enddo do ii=1,i do k=1,nl x(k,i+1)= x(k,i+1)+del(ii)* x(k,ii) ax(k,i+1)=ax(k,i+1)+del(ii)*ax(k,ii) enddo enddo else do k=1,nl x(k,1)=qr(k)+del(1)* x(k,1) ax(k,1)=ar(k)+del(1)*ax(k,1) enddo do ii=2,i do k=1,nl x(k,1 )= x(k,1)+del(ii)* x(k,ii) x(k,ii)=0. ax(k,1 )=ax(k,1)+del(ii)*ax(k,ii) ax(k,ii)=0. enddo enddo endif cconvergence test modes ************************************************** c if(ctest) then * c if(nlc/ner*ner.eq.nlc) then * c ier=nlc/ner * c eer=-1.e15 * c do 50 k=1,nl * c 50 eer=amax1(eer,abs(r(k))) * c err(ier)=eer * c endif * c endif * cconvergence test modes ************************************************** enddo 100 continue 200 continue eer=0. eem=-1.e15 do k=1,nl eer=eer+qr(k)**2 eem=amax1(eem,abs(qr(k))) enddo eer=eer/eer0 eem=eem/eem0 niter=nlc nitsm=nitsm+niter icount=icount+1 cconvergence test modes ************************************************** c if(ctest) then * c print 301, (err(ier),ier=1,nplt,1) * c 301 format(4x,5e11.4) * c do 400 ier=0,nplt * c xitr(ier)=ier*ner * c 400 err(ier)=alog10(err(ier)*dt ) * c plmx=float(itr*lord) * c call set(.1,.9,.1,.9,0.,plmx,-10.,0.,1) * c call labmod('(i4)','(f5.0)',4,4,2,2,20,20,0) * c call periml(4,10,5,2) * c call dashdc('$$$$$$$$$$$$$$$$$$$$',10,12) * c call curved(xitr,err,nplt+1) * c i1=int(102.4+409.6) * c call wtstr(cpux(i1),cpuy(50),'niter',3,0,0) * c call wtstr(cpux(17),cpuy(i1),'log e',3,90,0) * c call frame * c endif * cconvergence test modes ************************************************** return end subroutine precon_1(rhs,p,r,c11,c22,c33,d,iflg,jfl) include 'param.grid' parameter(n=nx,m=ny,l=nz,nl=nxyz) dimension rhs(n,m,l),p(n,m,l),r(n,m,l), . c11(n,m,l),c22(n,m,l),c33(n,m,l),d(n,m,l) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi compute available storage dimension e(n,m,0:l-1),f(n,m,0:l-1) . ,px(n,m,l),dgh(n,m,l),po(n,m,l) data beta/-1.e15/ data itr,line/2,1/ if(iflg.eq.0) then do i=1,nl p(i,1,1)=rhs(i,1,1) enddo return endif c omg=.7 c oms=1.-omg c dxi2=0.25*dxi*dxi c dzi2=0.25*dzi*dzi c do i=1,nl c c33(i,1)=d(i,1)*dzi2 c c11(i,1)=d(i,1)*dxi2 c dgh(i,1)=0. c po(i,1)=0. c p(i,1)=0. c r(i,1)=0. c enddo c if(line.eq.1) then c do k=1,l c do i=2,n-1 c dgh(i,k)=c11(i+1,k)+c11(i-1,k) c enddo c dgh(1,k)=c11(2,k)+c11(n-1,k) c dgh(n,k)=c11(2,k)+c11(n-1,k) c enddo c endif c c if(jfl.eq.0) then c if(line.eq.0) then c beta=-1.e15 c do i=1,nl c beta=amax1(beta,abs(c11(i,1))/d(i,1)) c enddo c beta=0.5/beta c else c beta=1. c endif c return c endif c beti=1./beta*(1-line) c c do 100 it=1,itr c do i=1,nl c r(i,1)=r(i,1)+d(i,1)*(beti*p(i,1)-rhs(i,1)) c . +dgh(i,1)*p(i,1) c enddo c do i=1,n c e(i,0)=1. c f(i,0)=0. c dn=d(i,1)*beti+2.*c33(i,2)+dgh(i,1) c e(i,1)=2.*c33(i,2)/dn c f(i,1)= r(i,1)/dn c enddo c do k=2,l-1 c do i=1,n c dn=c33(i,k+1)+c33(i,k-1)*(1.-e(i,k-2))+d(i,k)*beti c . + dgh(i,k) c e(i,k)= c33(i,k+1)/dn c f(i,k)=(c33(i,k-1)*f(i,k-2)+r(i,k))/dn c enddo c enddo c do i=1,n c dn=d(i,l)*beti+2.*(1.-e(i,l-2))*c33(i,l-1) c . + dgh(i,l) c p(i,l)=(r(i,l)+2.*f(i,l-2)*c33(i,l-1))/dn c p(i,l-1)=f(i,l-1)/(1.-e(i,l-1)) c enddo c do k=l-2,1,-1 c do i=1,n c p(i,k)=e(i,k)*p(i,k+2)+f(i,k) c enddo c enddo c c c if(line.eq.1) then c do i=1,nl c p(i,1)=oms*po(i,1)+omg*p(i,1) c po(i,1)= p(i,1) c enddo c endif c c if(it.eq.itr) go to 101 c do k=1,l c do i=2,n-1 c px(i,k)=c11(i,k)*(p(i+1,k)-p(i-1,k)) c enddo c px(1,k)=c11(1,k)*(p(2,k)-p(n-1,k)) c px(n,k)=c11(n,k)*(p(2,k)-p(n-1,k)) c enddo c c do k=1,l c do 91 i=2,n-1 c 91 r(i,k)=px(i+1,k)-px(i-1,k) c r(1,k)=(px(2,k)-px(n-1,k)) c r(n,k)=(px(2,k)-px(n-1,k)) c enddo c c 100 continue c 101 continue c return end subroutine laplc_2(p,r,u,v,w,n1,m1,l1) dimension p(n1,m1,l1),r(n1,m1,l1),u(n1,m1,l1),v(n1,m1,l1), 1 w(n1,m1,l1) include 'param.grid' parameter(n=nx,m=ny,l=nz) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi common /strtch/ height(nz),gac(nz) common /prof_d/ rho0(nz),th0(nz),th_e(nz),ux_e(nz),uy_e(nz) dimension px(n,m,l),py(n,m,l),pz(n,m,l),pe(n+1,m+1,l+1) nm=n*m ml=m*l nml=n*m*l dxil=.5*dxi dyil=.5*dyi dzil=.5*dzi compute pressure derivatives everywhere do 1 i=2,n-1 do 1 jk=1,ml 1 px(i,jk,1)= dxil*(p(i+1,jk,1)-p(i-1,jk,1)) do 18 jk=1,ml px(1,jk,1)=dxil*(p(2,jk,1)-p(n-1,jk,1)) px(n,jk,1)=dxil*(p(2,jk,1)-p(n-1,jk,1)) 18 continue j3=1 do 2 k=1,l do 2 i=1,n do 28 j=1+j3,m-j3 28 py(i,j,k)= dyil*(p(i,j+j3,k)-p(i,j-j3,k)) py(i,1,k)=dyil*(p(i,1+j3,k)-p(i,m-j3,k)) py(i,m,k)=dyil*(p(i,1+j3,k)-p(i,m-j3,k)) 2 continue do 3 k=2,l-1 do 3 ij=1,nm 3 pz(ij,1,k)=dzil*(p(ij,1,k+1)-p(ij,1,k-1)) do 38 ij=1,nm pz(ij,1,1)= dzi*(p(ij,1,2)-p(ij,1,1)) 38 pz(ij,1,l)= dzi*(p(ij,1,l)-p(ij,1,l-1)) compute interior pressure forces do 10 k=2,l-1 do 10 ij=1,nm u(ij,1,k)=-px(ij,1,k) v(ij,1,k)=-py(ij,1,k) 10 w(ij,1,k)=-pz(ij,1,k) compute pressure forces at the boundaries do 21 k=1,l,l-1 do 21 j=1,m do 21 i=1,n w(i,j,k)=0. u(i,j,k)=-px(i,j,k) 21 v(i,j,k)=-py(i,j,k) do 99 k=1,l coef=rho0(k)*gac(k) do 99 ij=1,nm u(ij,1,k)=coef*u(ij,1,k) v(ij,1,k)=coef*v(ij,1,k) 99 w(ij,1,k)=coef*w(ij,1,k) compute laplacian do 91 i=2,n-1 do 91 jk=1,ml 91 r(i,jk,1)=dxil*(u(i+1,jk,1)-u(i-1,jk,1)) do 911 jk=1,ml r(1,jk,1)=dxil*(u(2,jk,1)-u(n-1,jk,1)) r(n,jk,1)=dxil*(u(2,jk,1)-u(n-1,jk,1)) 911 continue j3=1 do 92 k=1,l do 92 i=1,n do 921 j=1+j3,m-j3 921 r(i,j,k)=r(i,j,k)+dyil*(v(i,j+j3,k)-v(i,j-j3,k)) r(i,1,k)=r(i,1,k)+dyil*(v(i,1+j3,k)-v(i,m-j3,k)) r(i,m,k)=r(i,m,k)+dyil*(v(i,1+j3,k)-v(i,m-j3,k)) 92 continue do 93 k=2,l-1 do 93 ij=1,nm 93 r(ij,1,k)=r(ij,1,k)+dzil*(w(ij,1,k+1)-w(ij,1,k-1)) do 931 ij=1,nm r(ij,1,1)=r(ij,1,1)+dzi*(w(ij,1,2)-w(ij,1,1)) r(ij,1,l)=r(ij,1,l)+dzi*(w(ij,1,l)-w(ij,1,l-1)) 931 continue do 94 k=1,l do 94 ij=1,nm 94 r(ij,1,k)=-r(ij,1,k)/(rho0(k)*gac(k)) return end subroutine prforc_2(p,pfx,pfy,pfz,u,v,w,n1,n2,n3) dimension p(n1,n2,n3) dimension u(n1,n2,n3),v(n1,n2,n3),w(n1,n2,n3) dimension pfx(n1,n2,n3),pfy(n1,n2,n3),pfz(n1,n2,n3) include 'param.grid' parameter(n=nx,m=ny,l=nz) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi dimension px(n,m,l),py(n,m,l),pz(n,m,l) dxil=.5*dxi dyil=.5*dyi dzil=.5*dzi compute pressure derivatives everywhere ml=m*l do 1 i=2,n-1 do 1 jk=1,ml 1 px(i,jk,1)=dxil*(p(i+1,jk,1)-p(i-1,jk,1)) do 18 jk=1,ml px(1,jk,1)=dxil*(p(2,jk,1)-p(n-1,jk,1)) px(n,jk,1)=dxil*(p(2,jk,1)-p(n-1,jk,1)) 18 continue j3=1 do 2 k=1,l do 2 i=1,n do 28 j=1+j3,m-j3 28 py(i,j,k)= dyil*(p(i,j+j3,k)-p(i,j-j3,k)) py(i,1,k)=dyil*(p(i,1+j3,k)-p(i,m-j3,k)) py(i,m,k)=dyil*(p(i,1+j3,k)-p(i,m-j3,k)) 2 continue nm=n*m do 3 k=2,l-1 do 3 ij=1,nm 3 pz(ij,1,k)=dzil*(p(ij,1,k+1)-p(ij,1,k-1)) do 38 ij=1,nm pz(ij,1,1)= dzi*(p(ij,1,2)-p(ij,1,1)) 38 pz(ij,1,l)= dzi*(p(ij,1,l)-p(ij,1,l-1)) compute interior pressure forces do 21 i=1,n do 21 j=1,m do 10 k=2,l-1 pfx(i,j,k)=u(i,j,k)-px(i,j,k) pfy(i,j,k)=v(i,j,k)-py(i,j,k) 10 pfz(i,j,k)=w(i,j,k)-pz(i,j,k) pfz(i,j,1)=0. pfz(i,j,l)=0. pfx(i,j,1)=u(i,j,1)-px(i,j,1) pfx(i,j,l)=u(i,j,l)-px(i,j,l) pfy(i,j,1)=v(i,j,1)-py(i,j,1) 21 pfy(i,j,l)=v(i,j,l)-py(i,j,l) return end subroutine rhsdiv_2(u,v,w,d,r,n1,m1,l1,iflg) dimension u(n1,m1,l1),v(n1,m1,l1),w(n1,m1,l1), 1 d(n1,m1,l1),r(n1,m1,l1) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi include 'param.grid' parameter(n=nx,m=ny,l=nz) nm=n*m ml=m*l nml=n*m*l do 200 i=1,nml 200 r(i,1,1)=0. dxil=.5*dxi dyil=.5*dyi dzil=.5*dzi do 1 i=2,n-1 do 1 j=1,ml 1 r(i,j,1)=dxil*(u(i+1,j,1)*d(i+1,j,1)-u(i-1,j,1)*d(i-1,j,1)) do 11 j=1,ml r(1,j,1)=dxil*(u(2,j,1)*d(2,j,1)-u(n-1,j,1)*d(n-1,j,1)) r(n,j,1)=dxil*(u(2,j,1)*d(2,j,1)-u(n-1,j,1)*d(n-1,j,1)) 11 continue j3=1 do 2 k=1,l do 2 i=1,n do 12 j=1+j3,m-j3 12 r(i,j,k)=r(i,j,k) 2 +dyil*(v(i,j+j3,k)*d(i,j+j3,k)-v(i,j-j3,k)*d(i,j-j3,k)) r(i,1,k)=r(i,1,k) 2 +dyil*(v(i,1+j3,k)*d(i,1+j3,k)-v(i,m-j3,k)*d(i,m-j3,k)) r(i,m,k)=r(i,m,k) 2 +dyil*(v(i,1+j3,k)*d(i,1+j3,k)-v(i,m-j3,k)*d(i,m-j3,k)) 2 continue do 3 k=2,l-1 do 3 i=1,nm 3 r(i,1,k)=r(i,1,k) 3 +dzil*(w(i,1,k+1)*d(i,1,k+1)-w(i,1,k-1)*d(i,1,k-1)) do 13 i=1,nm r(i,1,1)=r(i,1,1)+dzi*(w(i,1,2)*d(i,1,2)-w(i,1,1)*d(i,1,1)) 13 r(i,1,l)=r(i,1,l)+dzi*(w(i,1,l)*d(i,1,l)-w(i,1,l-1)*d(i,1,l-1)) if(iflg.ne.0) then do 4 i=1,nml 4 r(i,1,1)=iflg*r(i,1,1)/d(i,1,1) endif return end subroutine minmax(a,n,an,ax) dimension a(n) an= 1.e15 ax=-1.e15 do i=1,n an=amin1(a(i),an) ax=amax1(a(i),ax) enddo return end subroutine diagno_1(ux,uy,uz,th,nx,ny,nz,scr1,scr2,rho) dimension ux(nx,ny,nz),uy(nx,ny,nz),uz(nx,ny,nz), . th(nx,ny,nz),scr1(nx,ny,nz),scr2(nx,ny,nz) dimension rho(nx,ny,nz) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi nxyz=nx*ny*nz do i=1,nx do j=1,ny do k=1,nz scr2(i,j,k)=rho(i,j,k) enddo enddo enddo print 200, time 200 format(1x,' ****** analysis for time (min): ',f8.2) call minmax(ux,nxyz,amn,amx) print 201,amn,amx 201 format(1x,' --> min, max ux: ',2e12.4) call minmax(uy,nxyz,amn,amx) print 301,amn,amx 301 format(1x,' --> min, max uy: ',2e12.4) call minmax(uz,nxyz,amn,amx) print 202,amn,amx 202 format(1x,' --> min, max uz: ',2e12.4) cour=0. do i=1,nxyz cour=amax1(cour,abs(ux(i,1,1))*dt/dx+abs(uy(i,1,1))*dt/dy . +abs(uz(i,1,1))*dt/dz) enddo print 302,cour 302 format(1x,' --> max courant: ',e12.4) call minmax(th,nxyz,amn,amx) print 203,amn,amx 203 format(1x,' --> min, max th: ',2e12.4) if(cour.gt..9) then stop 'courant' endif return end subroutine diagno_2(th,ux,uz,nx,ny,nz,scr1,scr2,rho) dimension ux(nx,ny,nz),uz(nx,ny,nz), . th(nx,ny,nz),scr1(nx,ny,nz),scr2(nx,ny,nz) dimension rho(nx,ny,nz) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi nxyz=nx*ny*nz do i=1,nx do j=1,ny do k=1,nz scr2(i,j,k)=rho(i,j,k) enddo enddo enddo call minmax(th,nxyz,amn,amx) print 201,amn,amx 201 format(1x,' --> min, max qv: ',2e12.4) call minmax(ux,nxyz,amn,amx) print 301,amn,amx 301 format(1x,' --> min, max qc: ',2e12.4) call minmax(uz,nxyz,amn,amx) print 202,amn,amx 202 format(1x,' --> min, max qr: ',2e12.4) return end subroutine diagno_4(npx,xp,yp,zp,qcp,plic,radp) include 'param.grid' include 'param.sds' dimension qc(nx,ny,nz),con(nx,ny,nz),rad(nx,ny,nz) common /prof_d/ rho0(nz),th0(nz),th_e(nz),ux_e(nz),uy_e(nz) common /prof_m/ qv_e(nz),tm_e(nz) dimension xp(np),yp(np),zp(np) dimension qcp(np),plic(np),radp(np) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi dimension ply(nx,nz),plx(ny,nz),zz(nz) character*50 lhead do k=1,nz zz(k)=float(k-1)*dz enddo do i=1,nxyz rad(i,1,1)=0. con(i,1,1)=0. enddo dxh=.5*dx dyh=.5*dy dzh=.5*dz do ip=1,npx i=(xp(ip)+dxh)/dx + 1 j=(yp(ip)+dyh)/dy + 1 k=(zp(ip)+dzh)/dz + 1 i=max(1,min(nx,i)) j=max(1,min(ny,j)) k=max(1,min(nz,k)) con1=plic(ip) con(i,j,k)=con(i,j,k) + con1 rad(i,j,k)=rad(i,j,k) + radp(ip)*con1 enddo do i=1,nx-1 do j=1,ny-1 do k=1,nz if(con(i,j,k).gt.1.) then rad(i,j,k)=rad(i,j,k)/con(i,j,k) else rad(i,j,k)=0. endif con(i,j,k)=con(i,j,k)*1.e-6 ! to cm**-3 rad(i,j,k)=rad(i,j,k)*1.e6 ! to micron enddo enddo enddo cc cyclicity do j=1,ny-1 do k=1,nz con(nx,j,k)=con(1,j,k) rad(nx,j,k)=rad(1,j,k) enddo enddo do i=1,nx do k=1,nz con(i,ny,k)=con(i,1,k) rad(i,ny,k)=rad(i,1,k) enddo enddo call minmax(con,nxyz,amn,amx) print 202,amn,amx 202 format(1x,' --> min, max conc (cm**3): ',2e12.4) call minmax(rad,nxyz,amn,amx) print 203,amn,amx 203 format(1x,' --> min, max rad (micron): ',2e12.4) return end subroutine moist_init create parametrs for the model: cc mass, terminal velocity, diameter data ar,br,cr,dr /5.2e2,3.,130.,0.5/ data as,bs,cs,ds /2.5e-2,2.,4.,0.25/ cc collection ef., alpha, beta data er,alphr,betr /0.8, 1., 2./ data es,alphs,bets /0.2, .3, 3./ cc No data anor,anos /2*1.e7/ cc latent heats: data hlatv,hlats /2.53e6,2.84e6/ cc cloud droplet concentration (per cc) data dconc /70./ ! must be between 50 and 2000 cc limiting temperatures data tup,tdn /168.,153./ cc gammas: data gamb1r,gambd1r /6.0,11.7/ data gamb1s,gambd1s /2.0,2.56/ cc reference temperature and saturated vapor pressure: c data tt0,ee0 /273.16,611./ c data tt0,ee0 /283.16,1228./ data tt0,ee0 /293.16,2337.3/ common/rain_p0/ ar,br,cr,dr,er,alphr,betr,gamb1r,gambd1r,anor common/rain_p1/ dconc,ddisp common/snow_p0/ as,bs,cs,ds,es,alphs,bets,gamb1s,gambd1s,anos common/temp_p/ tup,tdn common/latent/hlatv,hlats common/reference/ tt0,ee0 common /const/ gg,cp,rg,rv data gg,cp,rg,rv /9.72,1005.,287.,461./ cc check consistency if(dconc.lt.50..or.dconc.gt.2000.) then print*,' *** inconsistent droplet concentration. stop.' stop 'dconc' endif cc calculate relative dispersion for Berry's autoconversion: ddisp=0.146-5.964e-2*alog(dconc/2000.) c print 2075,anor,anos,dconc,ddisp c2075 format(1x,' N0 in raindrop distr.: ',e15.3/ c 1 1x,' N0 in snowflake distr.: ',e15.3/ c 1 1x,' Berry parameters of cloud droplet spectrum:'/ c 1 1x,' droplet conc., relative disp.: ',2f12.4) return end subroutine integz(a,b,n1,n2,n3) dimension a(n1,n2,n3),b(n1,n2,n3) nm=n1*n2 nml=n1*n2*n3 do k=2,n3-1 do i=1,nm b(i,1,k)=0.25*(a(i,1,k+1)+2.*a(i,1,k)+a(i,1,k-1)) enddo enddo do i=1,nm b(i,1,1 )=0.5*(a(i,1, 2)+a(i,1, 1)) b(i,1,n3)=0.5*(a(i,1,n3)+a(i,1,n3-1)) enddo do i=1,nml a(i,1,1)=b(i,1,1) enddo return end subroutine prof_init include 'param.grid' common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi common /strtch/ height(nz),gac(nz) dimension xx(nx),zz(nz) common /prof_d/ rho0(nz),th0(nz),th_e(nz),ux_e(nz),uy_e(nz) common /prof_m/ qv_e(nz),tm_e(nz) common/temp_p/ tup,tdn common /const/ gg,cp,rg,rv common/reference/ tt0,ee0 common/latent/hlatv,hlats a=rg/rv c=hlatv/cp b=hlats/rv d=hlatv/rv e=-cp/rg tu=280. tl=299. do k=1,nz zz(k)=height(k) c th0(k) =tl + (tu-tl)*height(k)/height(nz) th0(k) =(tu+tl)/2. rho0(k)=1. th_e(k)=th0(k) tm_e(k)=th_e(k) ux_e(k)=0. uy_e(k)=0. pre=1.e5 tt=th0(k) delt=(tt-tt0)/(tt*tt0) esw=ee0*exp(d * delt) qvsw=a * esw /(pre-esw) qvs=qvsw qv_e(k)=qvs print*,'k,t,qv: ',k,th_e(k),qv_e(k) enddo return end subroutine zstrtch(zz,nz1,dz) cc define vertical grid for stretched coordinate, derive Jacobian dimension zz(nz1) include 'param.grid' common /strtch/ height(nz),gac(nz) dimension zb(500) if(nz.gt.500) stop 'grid formulation' top=float(nz-1)*dz cc exponent of stretching function: c ex1=5./4. c ex1=4.1/2.9 ex1=1. aa=top / top**ex1 cc vertical coordinate: pi=4.*atan(1.) do k=1,nz cc no stretched zb(k)=aa*zz(k)**ex1 cc stretched c zb(k)=top * (sin(0.5*pi*zz(k)/top))**2 height(k)=zb(k) enddo cc jacobian: do k=1,nz if(k.eq.1 ) gac(k)=(zb(2)-zb(1))/dz if(k.eq.nz) gac(k)=(zb(nz)-zb(nz-1))/dz if(k.ne.1 .and. k.ne.nz) gac(k)=(zb(k+1)-zb(k-1))/(2.*dz) print*,k,zb(k),gac(k) enddo cc check consistency (int gac ddzeta = H) sum1=.5*(gac(1)*dz + gac(nz)*dz) do i=2,nz-1 sum1=sum1+gac(i)*dz enddo print*,'int Jacobian before adj: ',sum1 cc adjust numerical jacobian: coe=float(nz-1)*dz/sum1 do i=1,nz gac(i)=gac(i)*coe enddo cc check: sum1=.5*(gac(1)*dz + gac(nz)*dz) do i=2,nz-1 sum1=sum1+gac(i)*dz enddo print*,'int Jacobian after adj: ',sum1 return end function cvmgm(a,ab,ac) if(ac.lt.0) then cvmgm=a else cvmgm=ab endif return end subroutine interp(ar,pl,nx,ny,nz,zz,hh) dimension ar(nx,ny,nz),pl(nx,ny,nz),zz(nz),hh(nz) do i=1,nx do j=1,ny pl(i,j,1)=ar(i,j,1) pl(i,j,nz)=ar(i,j,nz) do k=2,nz-1 do kk=2,nz iisn=kk-1 if(hh(kk).ge.zz(k)) go to 665 enddo print*,'should not be here, stop' stop 'loop 665' 665 continue coe2=(zz(k)-hh(iisn))/(hh(iisn+1)-hh(iisn)) pl(i,j,k)=coe2*ar(i,j,iisn+1) + (1.-coe2)*ar(i,j,iisn) enddo enddo enddo return end subroutine extr11(qc,n,an,ax) dimension qc(n) an=1.e20 ax=-1.e20 do i=1,n an=min(an,qc(i)) ax=max(ax,qc(i)) enddo return end subroutine extr_loc(ara,nx,ny,nz, 1 amn,imn,jmn,kmn,amx,imx,jmx,kmx) dimension ara(nx,ny,nz) amn=1.e15 amx=-1.e15 do i=1,nx do j=1,ny do k=1,nz amn=min(amn,ara(i,j,k) ) if(ara(i,j,k).eq.amn) then imn=i jmn=j kmn=k endif amx=max(amx,ara(i,j,k) ) if(ara(i,j,k).eq.amx) then imx=i jmx=j kmx=k endif enddo enddo enddo return end subroutine sd_adv_sed(ux,uxp,uy,uyp,uz,uzp,xp,yp,zp,radp,npx) include 'param.grid' include 'param.sds' COMMON /aerosol3/ s_activ(30),r_activ(30) dimension ux(nx,ny,nz),uxp(nx,ny,nz), 1 uy(nx,ny,nz),uyp(nx,ny,nz), 1 uz(nx,ny,nz),uzp(nx,ny,nz) dimension xp(np),yp(np),zp(np),radp(np) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi cc centered-in-time trajectory xxl=(nx-1)*dx yyl=(ny-1)*dy dxh=.5*dx dyh=.5*dy dzh=.5*dz do ip=1,npx cc predictor: i=(xp(ip)+dxh)/dx + 1 i=max(1,min(i,nx)) ip1=i+1 im1=i-1 if(i.eq.nx) ip1=2 if(i.eq. 1) im1=nx-1 j=(yp(ip)+dyh)/dy + 1 j=max(1,min(j,ny)) jp1=j+1 jm1=j-1 if(j.eq.ny) jp1=2 if(j.eq. 1) jm1=ny-1 k=(zp(ip)+dzh)/dz + 1 uup=.5*(ux(ip1,j,k)+ux(i,j,k)) uum=.5*(ux(im1,j,k)+ux(i,j,k)) delx=(xp(ip) + dxh - (i-1)*dx )/dx delx=max(0.,min(1.,delx)) uxc=(1.-delx)*uum + delx*uup xxf=xp(ip)+uxc*dt if(xxf.ge.xxl-dxh) xxf=xxf-xxl if(xxf.le. 0.-dxh) xxf=xxf+xxl vvp=.5*(uy(i,jp1,k)+uy(i,j,k)) vvm=.5*(uy(i,jm1,k)+uy(i,j,k)) dely=(yp(ip) + dyh - (j-1)*dy )/dy dely=max(0.,min(1.,dely)) uyc=(1.-dely)*vvm + dely*vvp yyf=yp(ip)+uyc*dt if(yyf.ge.yyl-dyh) yyf=yyf-yyl if(yyf.le. 0.-dyh) yyf=yyf+yyl if(k.gt.1 .and. k.lt.nz) then wwp=.5*(uz(i,j,k)+uz(i,j,k+1)) wwm=.5*(uz(i,j,k)+uz(i,j,k-1)) endif if(k.eq.1) then wwp=.5*(uz(i,j,k)+uz(i,j,k+1)) wwm=-wwp endif if(k.eq.nz) then wwm=.5*(uz(i,j,k)+uz(i,j,k-1)) wwp=-wwm endif delz=(zp(ip) + dzh - (k-1)*dz )/dz delz=max(0.,min(1.,delz)) uzc=(1.-delz)*wwm + delz*wwp zzf=zp(ip)+uzc*dt cc corrector: c do iter=1,0 ! no iteration do iter=1,1 ! one iteration c do iter=1,2 ! two iterations i=(xxf+dxh)/dx + 1 i=max(1,min(nx,i)) ip1=i+1 im1=i-1 if(i.eq.nx) ip1=2 if(i.eq. 1) im1=nx-1 j=(yyf+dyh)/dy + 1 j=max(1,min(ny,j)) jp1=j+1 jm1=j-1 if(j.eq.ny) jp1=2 if(j.eq. 1) jm1=ny-1 k=(zzf+dzh)/dz + 1 k=max(1,min(nz,k)) uup=.5*(ux(ip1,j,k)+ux(i,j,k)) uum=.5*(ux(im1,j,k)+ux(i,j,k)) uup1=.5*(uxp(ip1,j,k)+uxp(i,j,k)) uum1=.5*(uxp(im1,j,k)+uxp(i,j,k)) uup2=2.*uup-uup1 uum2=2.*uum-uum1 vvp=.5*(uy(i,jp1,k)+uy(i,j,k)) vvm=.5*(uy(i,jm1,k)+uy(i,j,k)) vvp1=.5*(uyp(i,jp1,k)+uyp(i,j,k)) vvm1=.5*(uyp(i,jm1,k)+uyp(i,j,k)) vvp2=2.*vvp-vvp1 vvm2=2.*vvm-vvm1 if(k.gt.1 .and. k.lt.nz) then wwp=.5*(uz(i,j,k)+uz(i,j,k+1)) wwm=.5*(uz(i,j,k)+uz(i,j,k-1)) wwp1=.5*(uzp(i,j,k)+uzp(i,j,k+1)) wwm1=.5*(uzp(i,j,k)+uzp(i,j,k-1)) endif if(k.eq.1) then wwp=.5*(uz(i,j,k)+uz(i,j,k+1)) wwm=-wwp wwp1=.5*(uzp(i,j,k)+uzp(i,j,k+1)) wwm1=-wwp1 endif if(k.eq.nz) then wwm=.5*(uz(i,j,k)+uz(i,j,k-1)) wwp=-wwm wwm1=.5*(uzp(i,j,k)+uzp(i,j,k-1)) wwp1=-wwm1 endif wwp2=2.*wwp-wwp1 wwm2=2.*wwm-wwm1 delx=(xxf + dxh - (i-1)*dx )/dx delx=max(0.,min(1.,delx)) uxc2=(1.-delx)*uum2 + delx*uup2 xxf=xp(ip)+(uxc+uxc2)*dt/2. c if(xxf.ge.xxl-dxh) xxf=xxf-xxl if(xxf.le. 0.-dxh) xxf=xxf+xxl dely=(yyf + dyh - (j-1)*dy )/dy dely=max(0.,min(1.,dely)) uyc2=(1.-dely)*vvm2 + dely*vvp2 yyf=yp(ip)+(uyc+uyc2)*dt/2. c if(yyf.ge.yyl-dyh) yyf=yyf-yyl if(yyf.le. 0.-dyh) yyf=yyf+yyl delz=(zzf + dzh - (k-1)*dz )/dz delz=max(0.,min(1.,delz)) uzc2=(1.-delz)*wwm2 + delz*wwp2 zzf=zp(ip)+(uzc+uzc2)*dt/2. enddo ! iterative corrector xp(ip)=xxf yp(ip)=yyf zp(ip)=zzf enddo ! do ip=1,npx cc sedimentation: pi=4.*atan(1.) rho_w=1.e3 do ip=1,npx cc setup fall velocity: qc0=4./3.*radp(ip)**3. * pi * rho_w dia=2.e6*radp(ip) amg=qc0 * 1.e3 ! in g if(dia.le.134.43) then vt=4.5795e5*amg**(2./3.) go to 101 endif if(dia.lt.1511.64) then vt=4.962e3*amg**(1./3.) go to 101 endif if(dia.lt.3477.84) then vt=1.732e3*amg**(1./6.) go to 101 endif vt=917. 101 continue vt=vt*1.e-2 zp(ip)=zp(ip)-vt*dt cc move SD to a random location with size of 1/2 activation: if(zp(ip).le.0.) then iclas=mod(ip,nppg) if(iclas.eq.0) iclas=nppg radp(ip)=0.5*r_activ(iclas) rand1=min(.999,max(.001,rand())) rand2=min(.999,max(.001,rand())) rand3=min(.999,max(.001,rand())) xp(ip)=rand1*xxl yp(ip)=rand2*yyl zp(ip)=rand3*zzl c print*,'moved:',radp(ip),xp(ip),yp(ip),zp(ip) endif enddo ! do ip=1,npx return end subroutine sd_adv(ux,uxp,uy,uyp,uz,uzp,xp,yp,zp,npx) include 'param.grid' include 'param.sds' dimension ux(nx,ny,nz),uxp(nx,ny,nz), 1 uy(nx,ny,nz),uyp(nx,ny,nz), 1 uz(nx,ny,nz),uzp(nx,ny,nz) dimension xp(np),yp(np),zp(np) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi cc centered-in-time trajectory xxl=(nx-1)*dx yyl=(ny-1)*dy dxh=.5*dx dyh=.5*dy dzh=.5*dz do ip=1,npx cc predictor: i=(xp(ip)+dxh)/dx + 1 i=max(1,min(i,nx)) ip1=i+1 im1=i-1 if(i.eq.nx) ip1=2 if(i.eq. 1) im1=nx-1 j=(yp(ip)+dyh)/dy + 1 j=max(1,min(j,ny)) jp1=j+1 jm1=j-1 if(j.eq.ny) jp1=2 if(j.eq. 1) jm1=ny-1 k=(zp(ip)+dzh)/dz + 1 uup=.5*(ux(ip1,j,k)+ux(i,j,k)) uum=.5*(ux(im1,j,k)+ux(i,j,k)) delx=(xp(ip) + dxh - (i-1)*dx )/dx delx=max(0.,min(1.,delx)) uxc=(1.-delx)*uum + delx*uup xxf=xp(ip)+uxc*dt if(xxf.ge.xxl-dxh) xxf=xxf-xxl if(xxf.le. 0.-dxh) xxf=xxf+xxl vvp=.5*(uy(i,jp1,k)+uy(i,j,k)) vvm=.5*(uy(i,jm1,k)+uy(i,j,k)) dely=(yp(ip) + dyh - (j-1)*dy )/dy dely=max(0.,min(1.,dely)) uyc=(1.-dely)*vvm + dely*vvp yyf=yp(ip)+uyc*dt if(yyf.ge.yyl-dyh) yyf=yyf-yyl if(yyf.le. 0.-dyh) yyf=yyf+yyl if(k.gt.1 .and. k.lt.nz) then wwp=.5*(uz(i,j,k)+uz(i,j,k+1)) wwm=.5*(uz(i,j,k)+uz(i,j,k-1)) endif if(k.eq.1) then wwp=.5*(uz(i,j,k)+uz(i,j,k+1)) wwm=-wwp endif if(k.eq.nz) then wwm=.5*(uz(i,j,k)+uz(i,j,k-1)) wwp=-wwm endif delz=(zp(ip) + dzh - (k-1)*dz )/dz delz=max(0.,min(1.,delz)) uzc=(1.-delz)*wwm + delz*wwp zzf=zp(ip)+uzc*dt cc corrector: c do iter=1,0 ! no iteration do iter=1,1 ! one iteration c do iter=1,2 ! two iterations i=(xxf+dxh)/dx + 1 i=max(1,min(nx,i)) ip1=i+1 im1=i-1 if(i.eq.nx) ip1=2 if(i.eq. 1) im1=nx-1 j=(yyf+dyh)/dy + 1 j=max(1,min(ny,j)) jp1=j+1 jm1=j-1 if(j.eq.ny) jp1=2 if(j.eq. 1) jm1=ny-1 k=(zzf+dzh)/dz + 1 k=max(1,min(nz,k)) uup=.5*(ux(ip1,j,k)+ux(i,j,k)) uum=.5*(ux(im1,j,k)+ux(i,j,k)) uup1=.5*(uxp(ip1,j,k)+uxp(i,j,k)) uum1=.5*(uxp(im1,j,k)+uxp(i,j,k)) uup2=2.*uup-uup1 uum2=2.*uum-uum1 vvp=.5*(uy(i,jp1,k)+uy(i,j,k)) vvm=.5*(uy(i,jm1,k)+uy(i,j,k)) vvp1=.5*(uyp(i,jp1,k)+uyp(i,j,k)) vvm1=.5*(uyp(i,jm1,k)+uyp(i,j,k)) vvp2=2.*vvp-vvp1 vvm2=2.*vvm-vvm1 if(k.gt.1 .and. k.lt.nz) then wwp=.5*(uz(i,j,k)+uz(i,j,k+1)) wwm=.5*(uz(i,j,k)+uz(i,j,k-1)) wwp1=.5*(uzp(i,j,k)+uzp(i,j,k+1)) wwm1=.5*(uzp(i,j,k)+uzp(i,j,k-1)) endif if(k.eq.1) then wwp=.5*(uz(i,j,k)+uz(i,j,k+1)) wwm=-wwp wwp1=.5*(uzp(i,j,k)+uzp(i,j,k+1)) wwm1=-wwp1 endif if(k.eq.nz) then wwm=.5*(uz(i,j,k)+uz(i,j,k-1)) wwp=-wwm wwm1=.5*(uzp(i,j,k)+uzp(i,j,k-1)) wwp1=-wwm1 endif wwp2=2.*wwp-wwp1 wwm2=2.*wwm-wwm1 delx=(xxf + dxh - (i-1)*dx )/dx delx=max(0.,min(1.,delx)) uxc2=(1.-delx)*uum2 + delx*uup2 xxf=xp(ip)+(uxc+uxc2)*dt/2. c if(xxf.ge.xxl-dxh) xxf=xxf-xxl if(xxf.le. 0.-dxh) xxf=xxf+xxl dely=(yyf + dyh - (j-1)*dy )/dy dely=max(0.,min(1.,dely)) uyc2=(1.-dely)*vvm2 + dely*vvp2 yyf=yp(ip)+(uyc+uyc2)*dt/2. c if(yyf.ge.yyl-dyh) yyf=yyf-yyl if(yyf.le. 0.-dyh) yyf=yyf+yyl delz=(zzf + dzh - (k-1)*dz )/dz delz=max(0.,min(1.,delz)) uzc2=(1.-delz)*wwm2 + delz*wwp2 zzf=zp(ip)+(uzc+uzc2)*dt/2. enddo ! iterative corrector xp(ip)=xxf yp(ip)=yyf zp(ip)=zzf enddo ! do ip=1,npx return end subroutine sd_phy(th,sup,ftt,fqv,nsd,xp,yp,zp,qcp, 1 plic,radp,radad,masad) include 'param.grid' include 'param.sds' include 'aerosol.inc' parameter(n=nx,m=ny,l=nz) dimension th(n,m,l),sup(n,m,l),ftt(n,m,l),fqv(n,m,l) c SDs: dimension xp(np),yp(np),zp(np) dimension qcp(np),plic(np),radp(np) real radad(np),masad(np) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi common /const/ gg,cp,rg,rv common/reference/ tt0,ee0 common/latent/hlat,hlats common /prof_m/ qv_e(nz),tm_e(nz) ccc function describing droplet growth rate drdt(r17)=gfact/r17 1 * (supers - exp(aaa/r17-bbb/(r17**3 -rs0**3)) + 1.) t00=tt0 t00i=1./t00 pi=4.*atan(1.) rho_w=1.e3 fact=4./3.*pi*rho_w a=rg/rv b=hlat/(rv*tt0) c=hlat/cp d=hlat/rv e=-cp/rg xxlv=hlat rho0=1. dxh=.5*dx dyh=.5*dy dzh=.5*dz do ip=1,nsd i=(xp(ip)+dxh)/dx + 1 j=(yp(ip)+dyh)/dy + 1 k=(zp(ip)+dzh)/dz + 1 i=max(1,min(n,i)) j=max(1,min(m,j)) k=max(1,min(l,k)) cc no interpolation supers=sup(i,j,k) rad_o=radp(ip) temp=th(i,j,k) prek=1.e5 temp=tm_e(1) sigma=76.1-0.155*(temp-273.16) sigma=sigma*1.e-3 ! in N/m difvp=2.11 * (temp/273.16)**1.94 * 101325./prek difvp=difvp*1.e-5 ! in m**2/s tcond=1.5e-11*temp**3 - 4.8e-8*temp**2 + 1.0e-4*temp - 3.9e-4 tcond=tcond ! in W/mK c print*,'-- in sd_phys: sigma,difvp,tcond: ',sigma,difvp,tcond rs0=radad(ip) sm0=masad(ip) r_dr=radp(ip) aaa=2*sigma/(rv*rho_w*temp) bbb=vanhof*sm0*amol_w/amol_s/(4./3.*pi*rho_w) c print*,aaa,bbb,' ---aaa,bbb' cc kinetic: coev=sqrt(2.*pi/rv/temp) coet=sqrt(2.*pi/rg/temp) coe1=r_dr/(r_dr+delta_v) + difvp/(r_dr*alpha_c)*coev difvpp=difvp / coe1 coe2=r_dr/(r_dr+delta_t) + tcond/(r_dr*alpha_t*rho0*cp)*coet tcondp=tcond / coe2 c print*,tcond,coe2,tcondp,' ---tcond,coe2,tcondp' thi=1./temp ess=ee0*exp( hlat/rv * (t00i - thi)) cc termodynamic factor with kinetic effects: gfact1=rho_w*rv*temp/ess/difvpp gfact2=rho_w*hlat/tcondp/temp * (hlat/rv/temp - 1.) gfact = 1./(gfact1+gfact2) c print*,ess,' ---ess' c print*,gfact1,' ---gfact1' c print*,gfact2,' ---gfact2' c print*,gfact,' ---gfact' r_old=r_dr cc do explicit, but substep if needed:: eps=0.05 rr=r_dr tt=0. it=0 100 continue it=it+1 tau=rr/abs(drdt(rr)) c print*,drdt(rr) dts=eps*tau if(tt+dts.ge.dt) dts=dt-tt ak1=dts*drdt(rr) rr = rr + ak1 tt = tt + dts if(dt-tt.lt.1.e-5) go to 701 go to 100 701 continue c print*,'== i,j,k,it,rr: ',i,j,k,it,rr r_dr=rr radp(ip)=max(0.,r_dr) qcp(ip)=fact * radp(ip)**3. *plic(ip)/rho0 dqc=fact * (radp(ip)**3.-rad_o**3.) 1 * plic(ip)/rho0 / dt fqv(i,j,k)=fqv(i,j,k)-dqc ftt(i,j,k)=ftt(i,j,k)+dqc*hlat/cp enddo ! ip loop cc apply cyclicity for forces: do j=1,n-1 do k=1,l ftt(n,j,k)=ftt(1,j,k) fqv(n,j,k)=fqv(1,j,k) enddo enddo do i=1,n do k=1,l ftt(i,m,k)=ftt(i,1,k) fqv(i,m,k)=fqv(i,1,k) enddo enddo return end subroutine sd_phy_old(theta,qv,sup,ftt,fqv, 1 nsd,xp,yp,zp,qcp,plic,radp) include 'param.grid' include 'param.sds' common /prof_d/ rho0(nz),th0(nz),th_e(nz),ux_e(nz),uy_e(nz) common /prof_m/ qv_e(nz),tm_e(nz) dimension sup(nx,ny,nz),ftt(nx,ny,nz),fqv(nx,ny,nz) dimension theta(nx,ny,nz),qv(nx,ny,nz) c SDs: dimension xp(np),yp(np),zp(np) dimension qcp(np),plic(np),radp(np) c SDs parameters assigned to gridboxes: real npxz(nx,ny,nz),coxz(nx,ny,nz),tausqe(nx,ny,nz) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi common /const/ gg,cp,rg,rv common/reference/ tt0,ee0 common/latent/hlat,hlats c parameters for activation: common /activ/ rad_sdb,dcon c flag for activation: dimension iflg(nx,ny,nz) bl(f00,f10,f01,f11,x,y) = 1 f00*(1.-x)*(1.-y) + f10*x*(1.-y) + f01*(1.-x)*y + f11*x*y cc activation and condensation growth: aa1=0.9152e-10 ! constant in droplet growth equation r_sh=1.86e-6 ! including kinetic effects c r_sh=0. ! excluding kinetic effects pi=4.*atan(1.) rho_w=1.e3 fact=4./3.*pi*rho_w a=rg/rv b=hlat/(rv*tt0) c=hlat/cp d=hlat/rv e=-cp/rg xxlv=hlat cc assign number of SDs and local concentartion to gridboxes do i=1,nxyz npxz(i,1,1)=0. coxz(i,1,1)=0. c tausqe(i,1,1)=0. enddo dxh=.5*dx dyh=.5*dy dzh=.5*dz do ip=1,nsd i=(xp(ip)+dxh)/dx + 1 j=(yp(ip)+dyh)/dy + 1 k=(zp(ip)+dzh)/dz + 1 i=max(1,min(nx,i)) j=max(1,min(ny,j)) k=max(1,min(nz,k)) npxz(i,j,k)=npxz(i,j,k)+1. coe=1. if(k.eq.1 .or. k.eq.nz) coe=2. con=plic(ip) coxz(i,j,k)=coxz(i,j,k)+con*coe c tausqe(i,j,k)=tausqe(i,j,k) + radp(ip)*con enddo cccc constant coefficient: c do i=1,nxyz c tausqe(i,1,1)=max(tausqe(i,1,1),.1) c tausqe(i,1,1)=1./(2.8e-4*tausqe(i,1,1)) c enddo ccc variable coefficient: c do k=1,nz c thetme=th_e(k)/tm_e(k) c prek=1.e5*thetme**e c do j=1,ny c do i=1,nx c tt=theta(i,j,k)/thetme c dv = 8.794e-5*tt**1.81/prek c coef=4.*pi*dv c tausqe(i,j,k)=max(tausqe(i,j,k),.1) c tausqe(i,j,k)=1./(coef*tausqe(i,j,k)) c enddo c enddo c enddo cc activation: do i=1,nx-1 do j=1,ny-1 do k=1,nz cc flag for activation: iflg(i,j,k) = 0 if(sup(i,j,k).lt.sact) go to 388 ! no activation needed cc local concentration conc=coxz(i,j,k) 222 continue if(conc .ge. ccnmx) go to 388 ! no need to activate... conc=conc+dcon cc new or additional activation: nsd=nsd+1 ipp=nsd if(nsd.gt.np) then print*,'after activation',i,j,k,np,nsd STOP 'stop nsd.gt.np' endif cccc random position at activation rand1=rand()-.5 rand2=rand()-.5 rand3=rand()-.5 if(k.eq.1 .or. k.eq.nz) rand3=rand() cc xp(ipp)=(i-1)*dx + rand1*dx yp(ipp)=(j-1)*dy + rand2*dy zp(ipp)=(k-1)*dz + rand3*dz if(k.eq.1) zp(ipp)= rand3*dz/2. if(k.eq.nz) zp(ipp)=(k-1)*dz - rand3*dz/2. plic(ipp)=dcon*rho0(k) radp(ipp)=rad_sdb + 0.2e-6*2.*(rand()-.5) qcp(ipp)=fact * radp(ipp)**3.*dcon dqc=qcp(ipp)/dt fqv(i,j,k)=fqv(i,j,k)-dqc ftt(i,j,k)=ftt(i,j,k)+dqc*hlat/cp*th_e(k)/tm_e(k) cc flag for activation: iflg(i,j,k) = 1 go to 222 388 continue ! no need to activate... enddo ! k loop enddo ! j loop enddo ! i loop cc reposition droplets after activation: call sd_rep_act(xp,yp,zp,nsd,iflg) cc condensation/evaporation: cc interp=0/1 - no/yes interpolation to SD position interp=0 c interp=1 dxh=.5*dx dyh=.5*dy dzh=.5*dz do ip=1,nsd i=(xp(ip)+dxh)/dx + 1 j=(yp(ip)+dyh)/dy + 1 k=(zp(ip)+dzh)/dz + 1 i=max(1,min(nx,i)) j=max(1,min(ny,j)) k=max(1,min(nz,k)) cc no interpolation supers=sup(i,j,k) rad_o=radp(ip) cc constant: aaa=aa1 cc variable c thetme=th_e(k)/tm_e(k) c prek=1.e5*thetme**e c tt=theta(i,j,k)/thetme c dv = 8.794e-5*tt**1.81/prek c thi=1./theta(i,j,k) c y=b*thetme*tt0*thi c ees=ee0*exp(b-y) c qvs=a*ees/(prek-ees) c dqsdt=xxlv*qvs/(rv*tt**2) * prek/(prek-ees) ! dqs/dT cc dqsdt = xxlv*qvs/(rv*tt**2) ! dqs/dT c ab = 1.+dqsdt*xxlv/cp c aaa=qvs*dv/rho_w/ab radn2=amax1(0., (rad_o+r_sh)**2 + 2.*aaa*supers*dt) rad_n=sqrt(radn2) - r_sh radp(ip)=max(0.,rad_n) qcp(ip)=fact * radp(ip)**3. *plic(ip)/rho0(k) dqc=fact * (radp(ip)**3.-rad_o**3.) 1 * plic(ip)/rho0(k) / dt fqv(i,j,k)=fqv(i,j,k)-dqc ftt(i,j,k)=ftt(i,j,k)+dqc*hlat/cp*th_e(k)/tm_e(k) cc logic for total evaporation: if(rad_n.lt.1.e-8) then radp(ip)=0. endif enddo ! ip loop cc eliminate empty SDs call remove(xp,yp,zp,qcp,plic,radp,nsd) cc apply cyclicity for forces: do j=1,ny-1 do k=1,nz ftt(nx,j,k)=ftt(1,j,k) fqv(nx,j,k)=fqv(1,j,k) enddo enddo do i=1,nx do k=1,nz ftt(i,ny,k)=ftt(i,1,k) fqv(i,ny,k)=fqv(i,1,k) enddo enddo return end subroutine sd_periodic(xp,yp,zp,qcp,plic,radp,nsd) include 'param.grid' include 'param.sds' dimension xp(np),yp(np),zp(np) dimension qcp(np),plic(np),radp(np) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi dxh=.5*dx dyh=.5*dy cc eliminate SDs from i=nx and j=ny planes: do ip=1,nsd i=(xp(ip)+dxh)/dx + 1 j=(yp(ip)+dyh)/dy + 1 i=max(1,min(nx,i)) j=max(1,min(ny,j)) if(i.eq.nx .or. j.eq.ny) then radp(ip)=0. endif enddo call remove(xp,yp,zp,qcp,plic,radp,nsd) cc add SDs from i=1 to i=nx nsd0=nsd xxl=(nx-1)*dx do ip=1,nsd0 i=(xp(ip)+dxh)/dx + 1 if(i.eq.1) then cc add SDs to the nx column: nsd=nsd+1 xp(nsd)=xp(ip) + xxl yp(nsd)=yp(ip) zp(nsd)=zp(ip) qcp(nsd)=qcp(ip) plic(nsd)=plic(ip) radp(nsd)=radp(ip) endif enddo cc add SDs from j=1 to j=ny nsd0=nsd yyl=(ny-1)*dy do ip=1,nsd0 j=(yp(ip)+dyh)/dy + 1 if(j.eq.1) then cc add SDs to the nx column: nsd=nsd+1 xp(nsd)=xp(ip) yp(nsd)=yp(ip) + yyl zp(nsd)=zp(ip) qcp(nsd)=qcp(ip) plic(nsd)=plic(ip) radp(nsd)=radp(ip) endif enddo return end subroutine remove(xp,yp,zp,qcp,plic,radp,nsd) include 'param.grid' include 'param.sds' cc this routine removes evaporated super-droplets dimension xp(np),yp(np),zp(np) dimension qcp(np),plic(np),radp(np) dimension xpn(np),ypn(np),zpn(np) dimension qcpn(np),plicn(np),radpn(np) cc zero new locations: do ip=1,np xpn(ip)=0. ypn(ip)=0. zpn(ip)=0. qcpn(ip)=0. plicn(ip)=0. radpn(ip)=0. enddo nsdn=0 do ip=1,nsd if(radp(ip).gt.1.e-7) then nsdn=nsdn+1 xpn(nsdn)=xp(ip) ypn(nsdn)=yp(ip) zpn(nsdn)=zp(ip) qcpn(nsdn)=qcp(ip) plicn(nsdn)=plic(ip) radpn(nsdn)=radp(ip) endif enddo nsd=nsdn do ip=1,nsd xp(ip)=xpn(ip) yp(ip)=ypn(ip) zp(ip)=zpn(ip) qcp(ip)=qcpn(ip) plic(ip)=plicn(ip) radp(ip)=radpn(ip) enddo do ip=nsd+1,np xp(ip)=0. yp(ip)=0. zp(ip)=0. qcp(ip)=0. plic(ip)=0. radp(ip)=0. enddo return end subroutine sd_rep_act(xp,yp,zp,npx,iflg) include 'param.grid' include 'param.sds' dimension iflg(nx,ny,nz) dimension xp(np),yp(np),zp(np) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi dxh=.5*dx dyh=.5*dy dzh=.5*dz do ip=1,npx i=(xp(ip)+dxh)/dx + 1 j=(yp(ip)+dyh)/dy + 1 k=(zp(ip)+dzh)/dz + 1 i=max(1,min(nx,i)) j=max(1,min(ny,j)) k=max(1,min(nz,k)) if(iflg(i,j,k).eq.1) then rand1=rand()-.5 rand2=rand()-.5 rand3=rand()-.5 if(k.eq.1 .or. k.eq.nz) rand3=rand() cc xp(ip)=(i-1)*dx + rand1*dx yp(ip)=(j-1)*dy + rand2*dy zp(ip)=(k-1)*dz + rand3*dz if(k.eq.1) zp(ip)= rand3*dz/2. if(k.eq.nz) zp(ip)=(k-1)*dz - rand3*dz/2. endif enddo return end subroutine walls(th,qv,fth,fqv,ux,uy,uz,fx,fy,fz) cc dumping at lateral walls include 'param.grid' dimension th(nx,ny,nz),qv(nx,ny,nz),ux(nx,ny,nz),uy(nx,ny,nz) dimension fth(nx,ny,nz),fqv(nx,ny,nz),fx(nx,ny,nz),fy(nx,ny,nz) dimension uz(nx,ny,nz),fz(nx,ny,nz) common /prof_d/ rho0(nz),th0(nz),th_e(nz),ux_e(nz),uy_e(nz) common /prof_m/ qv_e(nz),tm_e(nz) common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi common /strtch/ height(nz),gac(nz) common /boundary/ thsrf,qvsrf,thtop,qvtop,thwal,qvwal dimension ss1(nx,ny,nz),ss2(nx,ny,nz),ss3(nx,ny,nz),ss4(nx,ny,nz) dimension thsfl(nx,ny),qvsfl(nx,ny),usfl(nx,ny),vsfl(nx,ny) common /const/ gg,cp,rg,rv common/reference/ tt0,ee0 common/latent/hlatv,hlats cc cc relaxation time scale: tau=0.1 cc wall location: ii=10 jj=10 ccc NOTE: no 2 in qv and th forces: uncentered for SDs... do j=1,ny do k=1,nz do i=-1,1 fth(ii-i,j,k)=fth(ii-i,j,k)-(th(ii-i,j,k)-thwal)/tau fqv(ii-i,j,k)=fqv(ii-i,j,k)-(qv(ii-i,j,k)-qvwal)/tau fx(ii-i,j,k)= fx(ii-i,j,k)-(ux(ii-i,j,k))/tau * 2. fy(ii-i,j,k)= fy(ii-i,j,k)-(uy(ii-i,j,k))/tau * 2. fz(ii-i,j,k)= fz(ii-i,j,k)-(uz(ii-i,j,k))/tau * 2. enddo enddo enddo do i=1,nx do k=1,nz do j=-1,1 fth(i,jj-j,k)=fth(i,jj-j,k)-(th(i,jj-j,k)-thwal)/tau fqv(i,jj-j,k)=fqv(i,jj-j,k)-(qv(i,jj-j,k)-qvwal)/tau fx(i,jj-j,k)= fx(i,jj-j,k)-(ux(i,jj-j,k))/tau * 2. fy(i,jj-j,k)= fy(i,jj-j,k)-(uy(i,jj-j,k))/tau * 2. fz(i,jj-j,k)= fz(i,jj-j,k)-(uz(i,jj-j,k))/tau * 2. enddo enddo enddo return end subroutine top_bot_wall include 'param.grid' common/grid/ time,dt,dx,dy,dz,dti,dxi,dyi,dzi common/temp_p/ tup,tdn common /const/ gg,cp,rg,rv common/reference/ tt0,ee0 common/latent/hlatv,hlats common /boundary/ thsrf,qvsrf,thtop,qvtop,thwal,qvwal a=rg/rv c=hlatv/cp b=hlats/rv d=hlatv/rv e=-cp/rg temp_l=299. temp_u=280. pre=1.e5 thsrf=temp_l thtop=temp_u thwal=(temp_u+temp_l)/2. tt=temp_l delt=(tt-tt0)/(tt*tt0) esw=ee0*exp(d * delt) qvsrf=a * esw /(pre-esw) *0.95 c tt=temp_u delt=(tt-tt0)/(tt*tt0) esw=ee0*exp(d * delt) qvtop=a * esw /(pre-esw) *0.95 tt=thwal delt=(tt-tt0)/(tt*tt0) esw=ee0*exp(d * delt) qvwal=a * esw /(pre-esw) *0.95 print*,'--thsrf,qvsrf,thtop,qvtop: ',thsrf,qvsrf,thtop,qvtop print*,'--thwal,qvwal: ',thwal,qvwal return end subroutine boundry(theta,qv) include 'param.grid' dimension theta(nx,ny,nz),qv(nx,ny,nz) common /boundary/ thsrf,qvsrf,thtop,qvtop,thwal,qvwal do i=1,nx do j=1,ny theta(i,j, 1)=thsrf theta(i,j,nz)=thtop qv(i,j,1)=qvsrf qv(i,j,nz)=qvtop enddo enddo return end '\eof' pgf90 -r8 src.for