! module opdif_solve use params,only: nlev,mmx implicit none real,dimension(mmx) :: tp,da,shi,atp,daz,shiz,atpz,atpzz,bntr, | cntr,qqop,xlop,sha,shaz real :: brtr contains !----------------------------------------------------------------------- subroutine opdifsv(ut,dt,colfac) use cons,only: boltz use flds_ionatm,only: te,ti,xiop2p,xiop2d,xne,xinp,xihp,xiop4s use flds_atmos,only: tn,xno,xnn2,xno2,sht,xn2d,xnh,xnco2,xnh2, | xnh2o use flds_ionzrt,only: qop4s ! diag only use flds_modelz,only: gz use flds_ratecoef,only: rk1,rk2,rk8,rk10,rk11,rk12,rk13,rk14,rk15, | rk18,rk19,rk21,rk25,rk26,rk28 use flds_modelz,only: zp implicit none ! ! Args: real,intent(in) :: ut,dt,colfac ! ! Common /syspar/ is retained from old code for pdedif solver: real :: uround integer :: min,mout common/syspar/ uround,min,mout ! ! Local: integer,parameter :: npde=1 integer,parameter :: nwkcn=30000, niwcn=400 real :: usoln(npde,mmx),zpts(mmx),wkcn(nwkcn),tm,ddt integer :: mblft(npde),mbrht(npde),iwcn(niwcn) real,parameter :: eps=1.e-2 integer,parameter :: itype=2, nbug=1 real,parameter :: dzit=0.5, dziu=0.25, zmas=16.*1.66e-24 integer :: k,m,m1,mmx1,index,nnans real :: tr,a,tend real,parameter :: spval=1.e36 ! mblft(1) = 1 mbrht(1) = 2 tm = 0. uround = 7.e-15 ! /syspar/ for pdedif min = 5 ! /syspar/ for pdedif mout = 6 ! /syspar/ for pdedif ddt = 1.e-2 index = 1 m1 = nlev-mmx do k=1,mmx ! was do 1 m = m1+k tp(k) = 0.5*(te(m)+ti(m)) tr=0.5*(ti(m)+tn(m)) a=alog10(tr) da(k)=3.02e+17*tp(k)/(xno(m)*sqrt(tp(k))* | (1.08-0.139*a+4.5e-3*a*a)*colfac+19.9*xnn2(m)+19.5*xno2(m)) shi(k)=zmas*gz(m)/(2.*boltz*tp(k)) sha(k)=sht(m) atp(k)=alog(tp(k)) enddo ! k=1,mmx mmx1 = mmx-1 do k=2,mmx1 ! was do 2 daz(k)=(da(k+1)-da(k-1))/dzit shiz(k)=(shi(k+1)-shi(k-1))/dzit atpz(k)=(atp(k+1)-atp(k-1))/dzit atpzz(k)=(atp(k+1)+atp(k-1)-2.*atp(k))/(dzit*dzit) shaz(k)=(sha(k+1)-sha(k-1))/dzit enddo ! k=2,mmx1 daz(1)=(da(2)-da(1))/dziu atpz(1)=(atp(2)-atp(1))/dziu shiz(1)=(shi(2)-shi(1))/dziu shaz(1)=(sha(2)-sha(1))/dziu atpzz(1)=atpzz(2) daz(mmx)=(da(mmx)-da(mmx-1))/dziu atpz(mmx)=(atp(mmx)-atp(mmx-1))/dziu shiz(mmx)=(shi(mmx)-shi(mmx-1))/dziu shaz(mmx)=(sha(mmx)-sha(mmx-1))/dziu atpzz(mmx)=atpzz(mmx-1) do k=1,mmx ! was do 4 m=m1+k bntr(k)=da(k)*(atpz(k)/sha(k)+shi(k)) | +daz(k)/sha(k) cntr(k)=daz(k) | /sha(k)*(atpz(k)/sha(k)+shi(k))+da(k)*(atpzz(k) | /(sha(k)*sha(k))+shiz(k)/sha(k)-shaz(k)*atpz(k)/(sha(k)**3)) qqop(k)=(rk18(m)*xiop2p(m)*xno(m)+rk19(m)*xiop2p(m) | *xne(m)+rk21(m)*xiop2p(m)+rk25(m)*xiop2d(m)*xno(m)+rk26(m)* | xiop2d(m)*xne(m)+rk28(m)*xiop2d(m) | +rk8(m)*xinp(m)*xno(m)+rk11(m)*xihp(m)*xno(m)) xlop(k)=(rk1(m)*xno2(m)+rk2(m)*xnn2(m)+rk10(m)*xn2d(m)+rk12(m) | *xnh(m)+rk13(m)*xnco2(m)+rk14(m)*xnh2(m)+rk15(m)*xnh2o(m)) enddo ! k=1,mmx brtr=-(atpz(mmx)/sha(mmx)+shi(mmx)) tend=dt do k=1,mmx ! was do 3 m=m1+k zpts(k)=zp(m) usoln(1,k)=xiop4s(m) enddo ! k=1,mmx ! call fprint8('opdifsv for pdedif:',mmx, ! | da , sha , bntr , shaz , qqop , cntr, xlop , qop4s, ! | 'da','sha','bntr','shaz','qqop','cntr','xlop','qop4s') call pdedif(func,bnbdy,npde,mmx,tend,eps,zpts,itype,mblft, | mbrht,nbug,tm,ddt,index,usoln,wkcn,iwcn,nwkcn,niwcn) ! write(6,"('opdifsv after pdedif: tm=',e12.4,' ddt=',e12.4, ! | ' usoln=',/,(6e12.4))") tm,ddt,usoln ! call check_nans(usoln,npde,mmx,1,'usoln',nnans,0,spval,1,1) do k=1,mmx ! was do 24 m=m1+k xiop4s(m)=max(usoln(1,k),1.e-20) enddo ! write(6,"('opdifsv: xiop4s=',/,(6e12.4))") xiop4s end subroutine opdifsv !----------------------------------------------------------------------- subroutine func(npde,nzp2,zn,tm,um,uzn,an,bn,cn) ! ! User-provided subroutine for pdedif solver (see pdedif.F). ! use flds_ionzrt,only: qop4s implicit none ! ! Args: integer,intent(in) :: npde,nzp2 real,intent(in) :: zn(nzp2),um(npde,nzp2),uzn(npde,nzp2),tm real,intent(out) :: an(npde,nzp2),bn(npde,nzp2),cn(npde,nzp2) ! ! Local: integer :: m,mmx1,mmx2,m1,kk,k ! mmx1=mmx+1 mmx2=mmx+2 m1=nlev-mmx do kk=2,mmx1 ! was do 1 k=kk-1 m=m1+k an(1,kk)=da(k)/(sha(k)*sha(k)) bn(1,kk)=bntr(k)/sha(k)-da(k)*shaz(k)/(sha(k)**3) cn(1,kk)=qop4s(m)+qqop(k)+(cntr(k)-xlop(k))*um(1,kk) enddo an(1,1)=an(1,2) bn(1,1)=bn(1,2) cn(1,1)=cn(1,2) an(1,mmx2)=an(1,mmx1) bn(1,mmx2)=bn(1,mmx1) cn(1,mmx2)=cn(1,mmx1) ! write(6,"('func: an=',/,(6e12.4))") an ! write(6,"('func: bn=',/,(6e12.4))") bn ! write(6,"('func: cn=',/,(6e12.4))") cn end subroutine func !----------------------------------------------------------------------- subroutine bnbdy(npde,tm,zl,zr,al,ar,bl,br) ! ! User-provided subroutine for pdedif solver (see pdedif.F). ! use flds_ionatm,only: xiop4s ! ! Args: integer,intent(in) :: npde real,intent(in) :: tm,zl,zr real,intent(in) :: al(npde),ar(npde) real,intent(out) :: bl(npde),br(npde) ! ! Local: integer :: m real,parameter :: fluxe=0. ! old model ! real,parameter :: fluxe=-1.e+8 ! real,parameter :: fluxe=+1.e+8 ! m=nlev-mmx+1 bl(1)=xiop4s(m) br(1)=(brtr*ar(1)-fluxe/da(mmx))*sha(mmx) ! ar(1) goes to NaNQ on 6th call: ! write(6,"('bnbdy: brtr=',e12.4,' ar=',e12.4,' da=',e12.4, ! | ' sha=',e12.4)") brtr,ar(1),da(mmx),sha(mmx) end subroutine bnbdy !----------------------------------------------------------------------- end module opdif_solve