#include "dims.h" SUBROUTINE OPLUS use cons_module,only: kmax,kmaxm1,kmaxp1,imax,imaxh,imaxp2, | len1,len2,rmassinv,expz,cs,p0,dtx2inv,boltz,dtsmooth, | dtsmooth_div2,shapiro,grav,gask,re,dphi,dlamda use crates_module,only: rk10,rk16,rk17,rk18,rk21,rk22,rk23,rk24, | rk26,rk27 implicit none C **** C **** ADVANCES N(O+) BY ONE TIME STEP C **** C **** ON ENTRY, T7 = UPWARD NUMBER FLUX OF O+ AT TOP C **** #include "params.h" #include "fcom.h" #include "vscr.h" #include "index.h" #include "buff.h" #include "phys.h" #include "crates_tdep.h" #include "mwt.h" C **** COMMON BLOCK DESCRIBING MAGNETIC FIELD real :: bxm,bx,bxp,by,byp,bz,bzp,bmod,bmodp COMMON/MAGFLD/BXM(ZIMXP,2),BX(ZIMXP,ZJMX),BXP(ZIMXP,4), 1 BY(ZIMXP,ZJMX),BYP(ZIMXP,4),BZ(ZIMXP,ZJMX),BZP(ZIMXP,4), 2 BMOD(ZIMXP,ZJMX),BMODP(ZIMXP,2) C **** MOLECULAR WEIGHTS OF O+ AND N2D ! ! Local: integer :: KUTT(ZJMX),ntjm1k,ntk,ntjp1k,ntem1d,ntek,ntep1k,ntim1k, | ntik,ntip1k,i,nmsm1k,nmsk,ntem1k,nmsp1k,nps1mk,nps1k,nps1pk, | nps2mk,nps2k,nps2pk,nopm1k,nopk,nopp1k,k,ntem2k,ntep2k, | nopm2k,nopp2k,ntim2k,ntip2k,nuk,nvk,nwk,nopmk,noppk,nuik,nvik, | nojp2k,nojp1k,nojnmk,nojm1k,nojm2k,nwik,npn2dk,nqop2pk,nqop2dk, | nek,nqopk,nopnmk,nopnpk,nopmnk,k1,k2,n real :: rmop,explic,gmr,fmin,fmax character(len=80) :: title ! DATA KUTT/1,2,3,5,6,7,9,10,11,13,14,15,17,10*17,17,15,14,13,11,10, 19,7,6,5,3,2,1/ DATA RMOP/16./ DATA EXPLIC/1./ C **** C **** T1 = DIV(B) C **** CALL DIVB(T1,J) C **** C **** S4 = TP(J-1), S5 = TP(J), S6 = TP(J+1) C **** NTJM1K = NJM1+NT NTK = NJ+NT NTJP1K = NJP1+NT NTEM1K = NJM1+NTE NTEK = NJ+NTE NTEP1K = NJP1+NTE NTIM1K = NJM1+NTI NTIK = NJ+NTI NTIP1K = NJP1+NTI DO 40 I = 1,LEN2 S4(I,1) = .5*(F(I,NTEM1K)+F(I,NTIM1K)) S5(I,1) = .5*(F(I,NTEK)+F(I,NTIK)) S6(I,1) = .5*(F(I,NTEP1K)+F(I,NTIP1K)) 40 CONTINUE C **** C **** S13 = D(J-1), S14 = D(J), S15 = D(J+1) C **** NMSM1K = NJM1+NMS NMSK = NJ+NMS NMSP1K = NJP1+NMS NPS1MK = NJM1+NPS NPS1K = NJ+NPS NPS1PK = NJP1+NPS NPS2MK = NJM1+NPS2 NPS2K = NJ+NPS2 NPS2PK = NJP1+NPS2 CALL RRK(F(1,NTJM1K),F(1,NMSM1K),F(1,NPS1MK),F(1,NPS2MK),S4,S13,j) CALL RRK(F(1,NTK),F(1,NMSK),F(1,NPS1K),F(1,NPS2K),S5,S14,j) CALL RRK(F(1,NTJP1K),F(1,NMSP1K),F(1,NPS1PK),F(1,NPS2PK),S6,S15,j) C **** C **** S4 = 2.*TP(J-1), S5 = 2.*TP(J), S6 = 2.*TP(J+1) C **** DO 41 I = 1,LEN2 S4(I,1) = 2.*S4(I,1) S5(I,1) = 2.*S5(I,1) S6(I,1) = 2.*S6(I,1) 41 CONTINUE C **** C **** S10 = H(J-1), S11 = H(J), S12 = H(J+1) C **** DO 42 I = 1,LEN2 S10(I,1) = gask*F(I,NTJM1K)/(.5*(F(I,NMSM1K)+F(I,NMSM1K+1))* 1 grav) S11(I,1) = gask*F(I,NTK)/(.5*(F(I,NMSK)+F(I,NMSK+1))* 1 grav) S12(I,1) = gask*F(I,NTJP1K)/(.5*(F(I,NMSP1K)+F(I,NMSP1K+1))* 1 grav) 42 CONTINUE C **** C **** S7 = (D/(H*DZ)*2.*TP+M*G/R)*N(O+) (J-1) C **** S8 = (D/(H*DZ)*2.*TP+M*G/R)*N(O+) (J) C **** S9 = (D/(H*DZ)*2.*TP+M*G/R)*N(O+) (J+1) C **** NOPM1K = NJM1+NOP NOPK = NJ+NOP NOPP1K = NJP1+NOP CALL DIFFUS(S4,F(1,NOPM1K),S10,S7) CALL DIFFUS(S5,F(1,NOPK),S11,S8) CALL DIFFUS(S6,F(1,NOPP1K),S12,S9) C **** C **** S9 = (B(H).DEL(H))*(D/(H*DZ)*TP+M*G/R)*N(O+) (J) C **** CALL BDOTDH(S7,S8,S9,J,S9) C **** C **** S9 = D*BZ*S9 (J) C **** DO 43 K = 1,KMAX DO 44 I = 1,LEN1 S9(I,K) = S14(I,K)*BZ(I,J)*S9(I,K) 44 CONTINUE 43 CONTINUE C **** C **** S3 = 2.*TP*N(O+) (J-2) C **** S4 = 2.*TP*N(O+) (J-1) C **** S5 = 2.*TP*N(O+) (J) C **** S6 = 2.*TP*N(O+) (J+1) C **** S7 = 2.*TP*N(O+) (J+2) C **** NTEM2K = NJM2+NTE NTEP2K = NJP2+NTE NOPM2K = NJM2+NOP NOPP2K = NJP2+NOP NTIM2K = NJM2+NTI NTIP2K = NJP2+NTI DO 45 I = 1,LEN2 S3(I,1) = F(I,NOPM2K)*(F(I,NTEM2K)+F(I,NTIM2K)) S4(I,1) = S4(I,1)*F(I,NOPM1K) S5(I,1) = S5(I,1)*F(I,NOPK) S6(I,1) = S6(I,1)*F(I,NOPP1K) S7(I,1) = F(I,NOPP2K)*(F(I,NTEP2K)+F(I,NTIP2K)) 45 CONTINUE C **** C **** S1 = (B(H).DEL(H))*2.*TP*N(O+) (J-1) C **** S2 = (B(H).DEL(H))*2.*TP*N(O+) (J) C **** S3 = (B(H).DEL(H))*2.*TP*N(O+) (J+1) C **** CALL BDOTDH(S3,S4,S5,J-1,S1) CALL BDOTDH(S4,S5,S6,J,S2) CALL BDOTDH(S5,S6,S7,J+1,S3) C **** C **** S1 = S1*D(J-1), S2 = S2*D(J), S3 = S3*D(J+1) C **** DO 46 I = 1,LEN2 S1(I,1) = S1(I,1)*S13(I,1) S2(I,1) = S2(I,1)*S14(I,1) S3(I,1) = S3(I,1)*S15(I,1) 46 CONTINUE C **** C **** S8 = (B(H).DEL(H))*D*(B(H).DEL(H))*2.*TP*N(O+) (J) C **** CALL BDOTDH(S1,S2,S3,J,S8) C **** C **** S7 = (BZ*D/(H*DZ)+DIV(*B))*S2 C **** CALL BDZDVB(S2,T1,S11,J,S7) C **** C **** COLLECT NEW EXPLICIT TERMS IN S4 C **** ( RHS OF TRIDIAGONALSYSTEM ) C **** DO 47 I = 1,LEN2 S4(I,1) = -EXPLIC*(S7(I,1)+S8(I,1)+S9(I,1)) 47 CONTINUE C **** C **** S9 = H(J-1), S10 = H(J), S11 = H(J+1) C **** DO 48 I=1,LEN2 S9(I,1) = S10(I,1) S10(I,1) = S11(I,1) S11(I,1) = S12(I,1) 48 CONTINUE C **** C **** S6 = (B.U)*N(O+) (J-1) C **** S7 = (B.U)*N(O+) (J) C **** S8 = (B.U)*N(O+) (J+1) C **** NUK = NJM1+NU-1 NVK = NJM1+NV-1 NWK = NJM1+NW-1 NOPK = NJM1+NOP-1 DO 49 K = 1,KMAX NUK = NUK+1 NVK = NVK+1 NWK = NWK+1 NOPK = NOPK+1 DO 50 I = 1,LEN1 S6(I,K)=(BX(I,J-1)*F(I,NUK)+BY(I,J-1)*F(I,NVK)+S9(I,K)* 1 BZ(I,J-1)*.5*(F(I,NWK)+F(I,NWK+1)))*F(I,NOPK) 50 CONTINUE 49 CONTINUE NUK = NJ+NU-1 NVK = NJ+NV-1 NWK = NJ+NW-1 NOPK = NJ+NOP-1 DO 51 K = 1,KMAX NUK = NUK+1 NVK = NVK+1 NWK = NWK+1 NOPK = NOPK+1 DO 52 I = 1,LEN1 S7(I,K)=(BX(I,J)*F(I,NUK)+BY(I,J)*F(I,NVK)+S10(I,K)* 1 BZ(I,J)*.5*(F(I,NWK)+F(I,NWK+1)))*F(I,NOPK) 52 CONTINUE 51 CONTINUE NUK = NJP1+NU-1 NVK = NJP1+NV-1 NWK = NJP1+NW-1 NOPK = NJP1+NOP-1 DO 53 K = 1,KMAX NUK = NUK+1 NVK = NVK+1 NWK = NWK+1 NOPK = NOPK+1 DO 54 I = 1,LEN1 S8(I,K)=(BX(I,J+1)*F(I,NUK)+BY(I,J+1)*F(I,NVK)+S11(I,K)* 1 BZ(I,J+1)*.5*(F(I,NWK)+F(I,NWK+1)))*F(I,NOPK) 54 CONTINUE 53 CONTINUE C **** C **** NOW ADD TERMS TO F IN S4 C **** NOPMK=NJM1+NOP-1 NOPK=NJ+NOP-1 NOPPK=NJP1+NOP-1 NUIK=NUI-1 NVIK=NVI-1 DO 24 K=1,KMAX NOPMK=NOPMK+1 NOPK=NOPK+1 NOPPK=NOPPK+1 NUIK=NUIK+1 NVIK=NVIK+1 DO 24 I = 3,IMAXP2 S4(I,K)=S4(I,K)+1./(2.*re)*(1./(CS(J)*dlamda)* 1(BX(I,J)*(S7(I+1,K)-S7(I-1,K))+.5*(F(I,NUIK)+F(I,NUIK+1))* 2BMOD(I,J)**2*(F(I+1,NOPK)/BMOD(I+1,J)**2-F(I-1,NOPK)/BMOD(I-1,J) 3**2))+1./dphi*(BY(I,J)*(S8(I,K)-S6(I,K))+.5*(F(I,NVIK)+F(I,NVIK+1) 4)*BMOD(I,J)**2*(F(I,NOPPK)/BMOD(I,J+1)**2-F(I,NOPMK)/BMOD(I,J-1)** 52))) 24 CONTINUE C **** PERIODIC POINTS DO 33 I=1,2 DO 33 K=1,KMAX S4(I,K)=S4(I+IMAX,K) S4(I+IMAXP2,K)=S4(I+2,K) 33 CONTINUE C **** C **** T1 = DIV(B)/BZ C **** DO 1 I=1,LEN1 T1(I)=T1(I)/BZ(I,J) 1 CONTINUE C **** C **** S7=D(J-1), S8=D(J), S9=D(J+1) C **** DO 2 I=1,LEN2 S7(I,1)=S13(I,1) S8(I,1)=S14(I,1) S9(I,1)=S15(I,1) 2 CONTINUE C **** C **** S15 = 1./(H*DZ), S14 = TP = .5*(TI+TE) C **** NTIK=NJ+NTI NTEK=NJ+NTE DO 3 I=1,LEN2 S15(I,1)=1./(S10(I,1)*dz) S14(I,1)=.5*(F(I,NTIK)+F(I,NTEK)) 3 CONTINUE C **** C **** S13 = 2.*TP/(H(-)*DZ)+GMR (K = 1,KMAXP1) C **** S12 = 2.*TP/(H(+)*DZ)-GMR (K = 0,KMAX) C **** C **** GMR = G*M(0+)/(2.*R) GMR=grav*RMOP/(2.*gask) DO 4 I=1,LEN2-LEN1 S13(I,2)=2.*S14(I,2)*(.5*(S15(I,1)+S15(I,2)))+GMR S12(I,2)=2.*S14(I,1)*(.5*(S15(I,1)+S15(I,2)))-GMR 4 CONTINUE C **** UPPER AND LOWER BOUNDARIES DO 5 I=1,LEN1 S13(I,1)=2.*S14(I,1)*(1.5*S15(I,1)-0.5*S15(I,2))+GMR S13(I,KMAXP1)=2.*(2.*S14(I,KMAX)-S14(I,KMAXM1))*(1.5*S15(I,KMAX) 1-0.5*S15(I,KMAXM1))+GMR S12(I,1)=2.*(2.*S14(I,1)-S14(I,2))*(1.5*S15(I,1)-0.5*S15(I,2))-GMR S12(I,KMAXP1)=2.*S14(I,KMAX)*(1.5*S15(I,KMAX)-0.5*S15(I,KMAXM1)) 1-GMR 5 CONTINUE C **** C **** S11 = D(K+1/2) C **** DO 8 I=1,LEN2-LEN1 S11 (I,2)=.5*(S8(I,1)+S8(I,2)) 8 CONTINUE C **** UPPER AND LOWER BOUNDARIES DO 9 I=1,LEN1 S11(I,1)=(1.5*S8(I,1)-0.5*S8(I,2)) S11(I,KMAXP1)=(1.5*S8(I,KMAX)-.5*S8(I,KMAXM1)) 9 CONTINUE C **** C **** S7 = (DIV(B)+(DH*D*BZ)/(D*BZ) C **** C **** PERIODIC POINTS C **** DO 91 I = 1,2 DO 92 K = 1,3*KMAXP1 S9(I,K) = S9(I+IMAX,K) S9(I+IMAXP2,K) = S9(I+2,K) 92 CONTINUE 91 CONTINUE DO 93 K = 1,KMAX DO 94 I = 3,IMAXP2 S7(I,K) = T1(I)+1./(re*S8(I,K)*BZ(I,J)**2)*(BX(I,J)/CS(J)* | (S8(I+1,K)*BZ(I+1,J)-S8(I-1,K)*BZ(I-1,J))/(2.*dlamda)+ | BY(I,J)*(S9(I,K)*BZ(I,J+1)-S7(I,K)*BZ(I,J-1))/(2.*dphi)) 94 CONTINUE 93 CONTINUE C **** C **** PERIODIC POINTS C **** DO 96 I = 1,2 DO 97 K = 1,KMAXP1 S7(I,K) = S7(I+IMAX,K) S7(I+IMAXP2,K) = S7(I+2,K) 97 CONTINUE 96 CONTINUE C **** C **** S10 = (1./(H*DZ)-(DIV(B)+DH*D*BZ/(D*BZ))/(2*BZ))*BZ**2 C **** S9 = (1./(H*DZ)+(DIV(B)+DH*D*BZ/(D*BZ))/(2*BZ))*BZ**2 C **** DO 10 K=1,KMAX DO 10 I=1,LEN1 S10(I,K)=(S15(I,K)-.5*S7(I,K))*BZ(I,J)**2 S9(I,K)=(S15(I,K)+.5*S7(I,K))*BZ(I,J)**2 10 CONTINUE C **** C **** SHAPIRO SMOOTHER C **** SET S4 TO N(O+)/(2*DT) (N-1) C **** NOJP2K=NJP2+NOPNM NOJP1K=NJP1+NOPNM NOJNMK=NJ+NOPNM NOJM1K=NJM1+NOPNM NOJM2K=NJM2+NOPNM DO 36 I=1,LEN2 S1(I,1)=F(I,NOJNMK)-shapiro*(F(I,NOJP2K)+F(I,NOJM2K)-4.* | (F(I,NOJP1K)+F(I,NOJM1K))+6.*F(I,NOJNMK)) 36 CONTINUE DO 37 I=3,LEN2-2 S4(I,1)=S4(I,1)-(S1(I,1)-shapiro*(S1(I+2,1)+S1(I-2,1)-4.* | (S1(I+1,1)+S1(I-1,1))+6.*S1(I,1)))*dtx2inv 37 CONTINUE C **** C **** BEGIN COEFFS P,Q,R IN S1,S2,S3 C **** DO 11 I=1,LEN2 S1(I,1)=S10(I,1)*S11(I,1)*S12(I,1) S2(I,1)=-(S9(I,1)*S11(I,2)*S12(I,2)+S10(I,1)*S11(I,1)*S13(I,1)) S3(I,1)=S9(I,1)*S11(I,2)*S13(I,2) 11 CONTINUE C **** C **** S7 = (B.U) C **** DO 12 I=1,LEN2 C **** S6 = H S6(I,1)=1./(S15(I,1)*dz) 12 CONTINUE NUK=NJ+NU-1 NVK=NJ+NV-1 NWK=NJ+NW-1 DO 13 K=1,KMAX NUK=NUK+1 NVK=NVK+1 NWK=NWK+1 DO 13 I=1,LEN1 C **** S7 = B.U S7(I,K)=BX(I,J)*F(I,NUK)+BY(I,J)*F(I,NVK)+S6(I,K)*BZ(I,J)*.5* 1(F(I,NWK)+F(I,NWK+1)) 13 CONTINUE C **** C **** FINISH P AND R IN S1 AND S3 C **** NWIK=NWI-1 DO 14 K=1,KMAXM1 NWIK=NWIK+1 DO 14 I=1,LEN1 S1(I,K+1)=S1(I,K+1)+(BZ(I,J)*S7(I,K)+.5*(F(I,NWIK+1)+F(I,NWIK+2))) 1*.5*S15(I,K+1) S2(I,K)=S2(I,K)-.5*(F(I,NWIK)+F(I,NWIK+1))*6./re S3(I,K)=S3(I,K)-(BZ(I,J)*S7(I,K+1)+.5*(F(I,NWIK)+F(I,NWIK+1)))*.5* 1S15(I,K) 14 CONTINUE C **** UPPER AND LOWER BOUNDARIES NWIK=NWI DO 15 I=1,LEN1 S1(I,1)=S1(I,1)+(BZ(I,J)*(2.*S7(I,1)-S7(I,2))+.5*(F(I,NWIK)+ 1F(I,NWIK+1)))*.5*S15(I,1) S2(I,KMAX)=S2(I,KMAX)-.5*(F(I,NWIK+KMAX)+F(I,NWIK+KMAXM1))*6./re S3(I,KMAX)=S3(I,KMAX)-(BZ(I,J)*(2.*S7(I,KMAX)-S7(I,KMAXM1))+.5* 1(F(I,NWIK+KMAX)+F(I,NWIK+KMAXM1)))*.5*S15(I,KMAX) 15 CONTINUE C **** C **** ADDITIONS TO Q IN S2 C **** DO 16 K=1,KMAX DO 16 I=1,LEN1 S2(I,K)=S2(I,K)-S7(I,K)*T1(I)*BZ(I,J)-dtx2inv 16 CONTINUE C **** C **** UPPER BOUNDARY CONDITION C **** C **** T2 = A, T3= B NWIK=NWI+KMAX DO 17 I=1,LEN1 C **** C **** MODS FOR OPLUS -- UPPER BOUNDARY CONDITION C **** C T2(I) = .5*(2.*S7(I,KMAX)-S7(I,KMAXM1))*BZ(I,J)+.5*F(I,NWIK) C T2(I)=0. IS DIFFUSIVE EQUILIBRIUM T2(I) = 0. C T2(I) = 0.5*F(I,NWIK) IS OPFLUX AND WIND COMPONENT C T2(I) = 0.5*F(I,NWIK) T3(I) = -BZ(I,J)**2*S11(I,KMAXP1)*S12(I,KMAXP1)-T2(I) T2(I) = -BZ(I,J)**2*S11(I,KMAXP1)*S13(I,KMAXP1)+T2(I) C **** S2 = Q = Q+B/A*R S2(I,KMAX)=S2(I,KMAX)+T3(I)/T2(I)*S3(I,KMAX) C **** S4 = F = F -R/A*PHI S4(I,KMAX)=S4(I,KMAX)-T7(I)*S3(I,KMAX)/T2(I) S3(I,KMAX)=0. 17 CONTINUE C **** C **** C **** SOURCES AND SINKS C **** NTK=NJ+NT-1 NMSK=NJ+NMS-1 NPS1K=NJ+NPS-1 NPS2K=NJ+NPS2-1 NPN2DK=NJ+NPN2D-1 NQOP2PK=NQOP2P-1 NQOP2DK=NQOP2D-1 NEK=NJ+NE-1 NOPK=NJ+NOP-1 DO 26 K=1,KMAX NTK=NTK+1 NMSK=NMSK+1 NPS1K=NPS1K+1 NPS2K=NPS2K+1 NPN2DK=NPN2DK+1 NQOP2PK=NQOP2PK+1 NQOP2DK=NQOP2DK+1 NEK=NEK+1 NOPK=NOPK+1 DO 26 I=1,LEN1 C **** CALCULATE LOSS IN S13, SUBSTRACT FROM S2 XIOP2P(I,K) = .5*(F(I,NQOP2PK)+F(I,NQOP2PK+1))/((p0*expz(K)* | .5*(F(I,NMSK)+F(I,NMSK+1))/(boltz*F(I,NTK))*((RK16+ | RK17)*(1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv(3)+RK18* | F(I,NPS2K)*rmassinv(2)))+(RK19(I,K)+RK20(I,K))*F(I,NEK)+ | RK21+RK22) XIOP2D(I,K) = (.5*(F(I,NQOP2DK)+F(I,NQOP2DK+1))+(RK20(I,K)* | F(I,NEK)+RK22)*XIOP2P(I,K))/((p0*expz(K)*.5*(F(I,NMSK)+ | F(I,NMSK+1))/(boltz*F(I,NTK))*(RK23*(1.-F(I,NPS1K)- | F(I,NPS2K))*rmassinv(3)+RK24*F(I,NPS2K)*rmassinv(2)+ | RK26*F(I,NPS1K)*rmassinv(1)))+RK25(I,K)*F(I,NEK)+RK27) S13(I,K)=p0*expz(K)*.5*(F(I,NMSK)+F(I,NMSK+1))/(boltz* | F(I,NTK))*(RK1(I,K)*F(I,NPS1K)*rmassinv(1)+RK2(I,K)* | (1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv(3)+RK10*F(I,NPN2DK)/ | RMN2D) S2(I,K)=S2(I,K)-S13(I,K) 26 CONTINUE C **** ADD SOURCE TERM TO RHS (F) IN S4 NQOPK=NQOP-1 NEK=NJ+NE-1 NPS2K=NJ+NPS2-1 NMSK=NJ+NMS-1 NTK=NJ+NT-1 DO 27 K=1,KMAX NQOPK=NQOPK+1 NEK=NEK+1 NPS2K=NPS2K+1 NMSK=NMSK+1 NTK=NTK+1 ! ! 9/99: Bug fix from tgcm13: ! Remove * from *F(I,NEK) in following line: ! 1 *F(I,NEK)+RK21)*XIOP2P(I,K)-(RK25(I,K)*F(I,NEK)+ ! so RK19(I,K) is not raised to power F(I,NEK) ! DO 27 I=1,LEN1 S4(I,K)=S4(I,K)-.5*(F(I,NQOPK)+F(I,NQOPK+1))-(RK19(I,K)* 1 F(I,NEK)+RK21)*XIOP2P(I,K)-(RK25(I,K)*F(I,NEK)+ 2 RK27)*XIOP2D(I,K)-(RK18*XIOP2P(I,K)+RK24* 3 XIOP2D(I,K))* 4 F(I,NPS2K)*rmassinv(2)* 5 p0*expz(K)*.5*(F(I,NMSK)+F(I,NMSK+1))/(boltz* 6 F(I,NTK)) 27 CONTINUE C **** LOWER BOUNDARY CONDITION C **** C **** N(O+) = Q/L NQOPK=NQOP DO 28 I=1,LEN1 S2(I,1)=S2(I,1)-S1(I,1) S4(I,1)=S4(I,1)-2.*S1(I,1)*F(I,NQOPK)/(1.5*S13(I,1)-.5*S13(I,2)) S1(I,1)=0. 28 CONTINUE C **** C **** SOLVE TRIDIAGONAL SYSTEM C **** NOPK=NJNP+NOP CALL TRSOLV(S1,S2,S3,S4,S7,F(1,NOPK),LEN1,3,LEN1-2,KMAXP1,1,KMAX, 11) C **** C **** FILTER N(O+) C **** IF(KUTT(J).LT.IMAXH)THEN NOPK=NJNP+NOP CALL FILTER2(NOPK,KMAX,KUTT(J)) ENDIF C **** TIME SMOOTHING NOPNMK=NJ+NOPNM NOPK=NJ+NOP NOPNPK=NJNP+NOP NOPMNK=NJNP+NOPNM DO 29 I=1,LEN2 F(I,NOPMNK)=dtsmooth*F(I,NOPK)+dtsmooth_div2* | (F(I,NOPNMK)+F(I,NOPNPK)) 29 CONTINUE C **** PERIODIC POINTS K1=NJNP+NOP K2=K1+KMAXM1 DO 30 N=1,2 DO 31 I=1,2 DO 31 K=K1,K2 F(I,K)=F(I+IMAX,K) F(I+IMAXP2,K)=F(I+2,K) 31 CONTINUE K1=NJNP+NOPNM K2=K1+KMAXM1 30 CONTINUE ! ! KMAXP1: ! nopk = njnp+nop+kmax nopnmk = njnp+nopnm+kmax do i=1,len1 f(i,nopk) = 0. f(i,nopnmk) = 0. enddo C **** INSURE NON-NEGATIVE N(O+) NOPK=NJNP+NOP NOPNMK=NJNP+NOPNM DO 32 I=1,LEN2 if (F(I,NOPK) < 1.E-5) f(i,nopk) = 1.e-5 if (F(I,NOPNMK) < 1.E-5) f(i,nopnmk) = 1.e-5 32 CONTINUE RETURN END C