#include "dims.h" SUBROUTINE MINOR(NX,NXNM,RMX,PHIX,ALFAX,IBND,IBNDB) use input_module,only: difhor use cons_module,only: expz,rmassinv,rmass,imaxh,difk,imaxp4, | imax,kmax,len1,len2,len3,kmaxp1,imaxp2,grav,dtx2inv,avo, | dtsmooth,dtsmooth_div2,p0,boltz,shapiro,expzmid,expzmid_inv implicit none C **** C **** ADVANCES MINOR SPECIES BY ONE TIME STEP C **** NX = INDEX OF MINOR SPECIES FOR TIME T(N) C **** NXNM = INDEX OF MINOR SPECIES FOR TIME T(N-1) C **** RMX = MOLECULAR WEIGHT OF MINOR SPECIES C **** PHIX = DIFFUSION VECTOR, PHI(X,N),N=1,3 C **** ALFAX = THERMAL DIFFUSION COEFFICIENT C **** C **** LOWER BOUNDARY C **** IF IBNDB = 0, C **** T1 = A, T2 = B, T3 = C DEFINE LOWER BOUNDARY CONDITION C **** WHERE A*DPSX/DZ + B*PSX + C = 0. C **** IF IBNDB = 1, C **** T3 = TOTAL UPWARD NUMBER FLUX OF X AT LOWER BOUNDARY C **** UPPER BOUNDARY C **** IF IBND =0, C **** T4 = DIFFUSIVE UPWARD NUMBER FLUX OF X AT UPPER C **** BOUNDARY C **** IF IBND =1, C **** T4 = TOTAL UPWARD NUMBER FLUX OF X AT UPPER C **** BOUNDARY C **** IF IBND =2, C **** T4 = VALUE OF PSX AT UPPER BOUNDARY C **** C **** SOURCES C **** S1 = SX/N(X), WHERE SX IS PORTION OF NUMBER DENSITY C **** SOURCE PROPORTIONAL TO N(X), THE MINOR SPECIES NUMBER C **** DENSITY (LEVEL K+1/2) C **** S2 = S0, PORTION OF NUMBER DENSITY SOURCE INDEPENDENT C **** OF N(X). (LEVEL K+1/2) C **** #include "params.h" #include "fcom.h" #include "vscr.h" #include "index.h" #include "buff.h" #include "phys.h" #include "cmpbnd.h" ! ! Args: integer,intent(in) :: nx,nxnm,ibnd,ibndb real,intent(in) :: rmx,alfax ! ! Local: integer :: KUTT(ZJMX),i,k,nps1k,nps2k,ntk,nxjp2k,nxjp1k,nxjnmk, | nxjm1k,nxjm2k,nwk,nxnpk,nxnmk,nxk,nxmnk,nxjk,nmxjm1k,nmxjk, | nmxjp1k,k1,k2,n real :: PHIX(3),PHI(2,3),small,tau,t00,alfa12,alfa21,alfax1,alfax2 ! 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 PHI/0.,0.673,1.35,0.,1.11,0.769/,TAU/1.86E+3/,T00/273./, C **** C **** VALUE OF 'SMALL' FOR MINOR C **** 1 SMALL/1.E-12/ C **** S13 = HORIZONAL ADVECTION CALL ADVEC(NX,S3) DO 24 I=1,LEN3 S13(I,1)=S3(I,1) 24 CONTINUE C **** PERIODIC POINTS DO 42 I=1,2 DO 42 K=1,KMAX S13(I,K)=S13(I+IMAX,K) S13(I+IMAXP2,K)=S13(I+2,K) 42 CONTINUE C **** S14=SX/N(X)(K+1/2), S15=S0(K+1/2) DO 1 I=1,LEN2 S15(I,1)=S2(I,1) S14(I,1)=S1(I,1) 1 CONTINUE C **** T4 = PS1(-1/2), T5 = PS2(-1/2), T6 = MBAR(-1/2) NPS1K=NJ+NPS NPS2K=NPS1K+KMAXP1 DO 2 I=1,LEN1 T7(I)=T4(I) T4(I) = B(I,1,1)*F(I,NPS1K)+B(I,1,2)*F(I,NPS2K)+FB(I,1) T5(I) = B(I,2,1)*F(I,NPS1K)+B(I,2,2)*F(I,NPS2K)+FB(I,2) T6(I)=1./(T4(I)*rmassinv(1)+T5(I)*rmassinv(2)+(1.-T4(I)-T5(I))* | rmassinv(3)) 2 CONTINUE C **** S12 = MBAR(K+1/2) DO 3 I=1,LEN3 S12(I,1)=1./(F(I,NPS1K)*rmassinv(1)+F(I,NPS2K)*rmassinv(2)+ | (1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv(3)) 3 CONTINUE C **** S11 = MBAR(K), S10 = DM/DZ(K), S9 = PS2(K), C **** S8 = PS1(K), S7 = DPS2/DZ(K), S6 = DPS1/DZ(K) C **** LOWER BOUNDARY DO 4 I=1,LEN1 S11(I,1)=.5*(T6(I)+S12(I,1)) S10(I,1)=(S12(I,1)-T6(I))/dz S9(I,1)=.5*(T5(I)+F(I,NPS2K)) S8(I,1)=.5*(T4(I)+F(I,NPS1K)) S7(I,1)=(F(I,NPS2K)-T5(I))/dz S6(I,1)=(F(I,NPS1K)-T4(I))/dz T4(I)=T7(I) 4 CONTINUE C **** LEVELS K=2,KMAXP1 DO 5 I=LEN1+1,LEN3 S11(I,1)=.5*(S12(I,1)+S12(I-LEN1,1)) S10(I,1)=(S12(I,1)-S12(I-LEN1,1))/dz S9(I,1)=.5*(F(I,NPS2K)+F(I,NPS2K-1)) S8(I,1)=.5*(F(I,NPS1K)+F(I,NPS1K-1)) S7(I,1)=(F(I,NPS2K)-F(I,NPS2K-1))/dz S6(I,1)=(F(I,NPS1K)-F(I,NPS1K-1))/dz 5 CONTINUE C **** S5=T(TOT,K) C **** LOWER AND UPPER BOUNDARIES NTK=NJ+NT+KMAX DO 6 I=1,LEN1 S5(I,1)=F(I,NTK) S5(I,KMAXP1)=F(I,NTK-1) 6 CONTINUE C **** LEVELS K=2,KMAX NTK=NJ+NT DO 7 K=2,KMAX NTK=NTK+1 DO 7 I=1,LEN1 S5(I,K)=.5*(F(I,NTK)+F(I,NTK-1)) 7 CONTINUE C **** EVALUATE ALFA12, ALFA21, ALFAM1, ALFAM2 ALFA12=PHI(1,2)-PHI(1,3) ALFA21=PHI(2,1)-PHI(2,3) ALFAX1=PHIX(1)-PHIX(3) ALFAX2=PHIX(2)-PHIX(3) C **** S15=(S0*MX/NMBAR)(K+1/2) NTK=NJ+NT-1 DO 8 K=1,KMAX NTK=NTK+1 DO 8 I=1,LEN1 S15(I,K)=S15(I,K)*RMX*boltz*F(I,NTK)/(p0* 1expz(K)*S12(I,K)) 8 CONTINUE C **** S1 = ALFA(1,1,K), S2 = ALFA(1,2,K), C **** S2 = ALFA(2,1,K), S4 = ALFA(2,2,K) DO 9 I=1,LEN3 S1(I,1)=-(PHI(1,3)+ALFA12*S9(I,1)) S2(I,1)=ALFA12*S8(I,1) S3(I,1)=ALFA21*S9(I,1) S4(I,1)=-(PHI(2,3)+ALFA21*S8(I,1)) C **** S12=EX(K) S12(I,1)=((ALFAX1*S4(I,1)-ALFAX2*S3(I,1))*(S6(I,1)-(1.-(RMASS(1)+ | S10(I,1))/S11(I,1))*S8(I,1))+(ALFAX2*S1(I,1)-ALFAX1*S2(I,1))* | (S7(I,1)-(1.-(RMASS(2)+S10(I,1))/S11(I,1))*S9(I,1)))/(S1(I,1)* | S4(I,1)-S2(I,1)*S3(I,1))+1.-(RMX+S10(I,1))/S11(I,1) C **** S10 = (DM/DZ)/MBAR S10(I,1)=S10(I,1)/S11(I,1) C **** S11 = AX(K) S11(I,1)=-S11(I,1)/(TAU*RMASS(3))*(T00/S5(I,1))**0.25/(PHIX(3)+ | ALFAX1*S8(I,1)+ALFAX2*S9(I,1)) 9 CONTINUE C **** S12=EX-ALFAX*D/DS(LN(T(TOT)) (THERMAL DIFFUSION TERM) DO 21 I=LEN1+1,LEN2 S12(I,1)=S12(I,1)-ALFAX*(S5(I+LEN1,1)-S5(I-LEN1,1))/(2.*dz* | S5(I,1)) 21 CONTINUE C **** LOWER BOUNDARY AND UPPER BOUNDARY DO 22 I=1,LEN1 S12(I,1)=S12(I,1)-ALFAX*(S5(I+LEN1,1)-S5(I,1))/(dz*S5(I,1)) 22 CONTINUE DO 23 I=1,LEN1 S12(I+LEN2,1)=S12(I+LEN2,1)-ALFAX*(S5(I+LEN2,1)- | S5(I+LEN2-LEN1,1))/(dz*S5(I+LEN2,1)) 23 CONTINUE C **** T7=DFACT DO 10 I=1,LEN1 T7(I)=1. 10 CONTINUE IF(difhor.EQ.1)CALL DFACT(T7,J,IMAXP4) C **** SHAPIRO SMOOTHER, S10=NXNM NXJP2K=NJP2+NXNM NXJP1K=NJP1+NXNM NXJNMK=NJ+NXNM NXJM1K=NJM1+NXNM NXJM2K=NJM2+NXNM DO 11 I=1,LEN3 S9(I,1)=S10(I,1) S8(I,1)=F(I,NXJNMK)-shapiro*(F(I,NXJP2K)+F(I,NXJM2K)-4.* | (F(I,NXJP1K)+F(I,NXJM1K))+6.*F(I,NXJNMK)) 11 CONTINUE DO 12 I=3,LEN3-2 S10(I,1)=S8(I,1)-shapiro*(S8(I+2,1)+S8(I-2,1)-4.*(S8(I+1,1)+ | S8(I-1,1))+6.*S8(I,1)) 12 CONTINUE C **** S1 = P, S2 = Q, S3 = R, S4= F NWK=NJNP+NW-1 DO 13 K=1,KMAX NWK=NWK+1 DO 13 I=1,LEN1 S1(I,K)=S11(I,K)/dz*(1./dz+.5*S12(I,K))-expz(K)*(expzmid_inv* | DIFK(K)*T7(I)*(1./dz-.5*S9(I,K))+0.25*(F(I,NWK)+F(I,NWK+1)))/dz S3(I,K)=S11(I,K+1)/dz*(1./dz-.5*S12(I,K+1))-expz(K)*(expzmid* | DIFK(K+1)*T7(I)*(1./dz+.5*S9(I,K+1))-0.25*(F(I,NWK)+F(I,NWK+1))) | /dz S2(I,K)=-(S11(I,K)/dz*(1./dz-.5*S12(I,K))+S11(I,K+1)/dz*(1./ | dz+.5*S12(I,K+1)))+expz(K)*((expzmid_inv*DIFK(K)*(1./dz+.5* | S9(I,K))+expzmid*DIFK(K+1)*(1./dz-.5*S9(I,K+1)))*T7(I)/dz- | S14(I,K)+dtx2inv) S4(I,K)=expz(K)*(S10(I,K)*dtx2inv-S13(I,K)+S15(I,K)) 13 CONTINUE C **** IF IBNDB = 1, SET FLUX BOUNDARY CONDITION AT BOTTOM IF(IBNDB.EQ.1)THEN NWK=NJNP+NW DO 26 I=1,LEN1 T1(I)=1. T2(I)=S9(I,1)-F(I,NWK)/DIFK(1) T3(I)=T3(I)*grav*RMX/(DIFK(1)*p0*expz(1)*expzmid_inv*avo) 26 CONTINUE ENDIF C **** BOUNDARIES C **** MODIFY EX(KMAXP1) IF IBND=1 IF(IBND.EQ.1)THEN NWK=NJNP+NW+KMAX DO 25 I=1,LEN1 T7(I)=expz(KMAX)*expzmid*F(I,NWK)/S11(I,KMAXP1) S12(I,KMAXP1)=S12(I,KMAXP1)-T7(I) 25 CONTINUE ENDIF DO 14 I=1,LEN1 C **** LOWER BOUNDARY S2(I,1)=S2(I,1)+S1(I,1)*(T1(I)+.5*T2(I)*dz)/(T1(I)-.5*T2(I)*dz 1) S4(I,1)=S4(I,1)-S1(I,1)*T3(I)*dz/(T1(I)-.5*T2(I)*dz) S1(I,1)=0. C **** UPPER BOUNDARY IF(IBND.EQ.0.OR.IBND.EQ.1)THEN C **** UPWARD NUMBER FLUX IS IN T4 S1(I,KMAXP1)=1.+.5*dz*S12(I,KMAXP1) S2(I,KMAXP1)=S1(I,KMAXP1)-2. S3(I,KMAXP1)=0. S4(I,KMAXP1)=-grav*RMX*T4(I)*dz/(p0*S11(I,KMAXP1)*avo) ELSE C **** VALUE AT UPPER BNDRY IN T4 S1(I,KMAXP1)=.5 S2(I,KMAXP1)=.5 S3(I,KMAXP1)=0. S4(I,KMAXP1)=T4(I) ENDIF 14 CONTINUE C **** SOLVE TRIDIAGONAL SYSTEM NXNPK=NJNP+NX CALL TRSOLV(S1,S2,S3,S4,S7,F(1,NXNPK),LEN1,3,LEN1-2,KMAXP1,1, 1KMAXP1,1) C **** FILTER NX IF(KUTT(J).LT.IMAXH)THEN NXNPK=NJNP+NX CALL FILTER2(NXNPK,KMAXP1,KUTT(J)) ENDIF C **** TIME SMOOTHING NXNMK=NJ+NXNM-1 NXK=NJ+NX-1 NXNPK=NJNP+NX-1 NXMNK=NJNP+NXNM-1 DO 17 K = 1,KMAXP1 NXNMK=NXNMK+1 NXK=NXK+1 NXNPK=NXNPK+1 NXMNK=NXMNK+1 DO 17 I=3,IMAX+2 F(I,NXMNK)=dtsmooth*F(I,NXK)+dtsmooth_div2* | (F(I,NXNMK)+F(I,NXNPK)) 17 CONTINUE C **** PERIODIC POINTS NXJM1K = NJM1 + NX - 1 NXJK = NJ + NX - 1 NXJP1K = NJP1 + NX - 1 NMXJM1K = NJM1 + NXNM - 1 NMXJK = NJ + NXNM - 1 NMXJP1K = NJP1 + NXNM - 1 DO 203 K = 1,KMAX + 1 NXJM1K = NXJM1K + 1 NXJK= NXJK + 1 NXJP1K = NXJP1K + 1 NMXJM1K = NMXJM1K + 1 NMXJK = NMXJK + 1 NMXJP1K = NMXJP1K + 1 DO 203 I = 1,2 F(I,NXJM1K) = F(I+IMAX,NXJM1K) F(I+IMAX+2,NXJM1K) = F(I+2,NXJM1K) F(I,NXJK) = F(I+IMAX,NXJK) F(I+IMAX+2,NXJK) = F(I+2,NXJK) F(I,NXJP1K) = F(I+IMAX,NXJP1K) F(I+IMAX+2,NXJP1K) = F(I+2,NXJP1K) F(I,NMXJM1K) = F(I+IMAX,NMXJM1K) F(I+IMAX+2,NMXJM1K) = F(I+2,NMXJM1K) F(I,NMXJK) = F(I+IMAX,NMXJK) F(I+IMAX+2,NMXJK) = F(I+2,NMXJK) F(I,NMXJP1K) = F(I+IMAX,NMXJP1K) F(I+IMAX+2,NMXJP1K) = F(I+2,NMXJP1K) 203 CONTINUE C **** INSURE NON-NEGATIVE NX NXJM1K = NJM1 + NX NXJK = NJ + NX NXJP1K = NJP1 + NX NMXJM1K = NJM1 + NXNM NMXJK = NJ + NXNM NMXJP1K = NJP1 + NXNM DO 201 K = 2,KMAX NXJM1K = NXJM1K + 1 NXJK = NXJK + 1 NXJP1K = NXJP1K + 1 NMXJM1K = NMXJM1K + 1 NMXJK = NMXJK + 1 NMXJP1K = NMXJP1K + 1 DO 201 I = 3,IMAX+2 S11(I,K) = 0.1*(F(I+1,NXJK)+F(I-1,NXJK)+F(I,NXJM1K)+ 1 F(I,NXJP1K)+F(I,NXJK-1)+F(I,NXJK+1)+F(I-1,NXJM1K)+ 2 F(I-1,NXJP1K)+F(I+1,NXJM1K)+F(I+1,NXJP1K)) S12(I,K) = 0.1*(F(I+1,NMXJK)+F(I-1,NMXJK)+F(I,NMXJM1K)+ 1 F(I,NMXJP1K)+F(I,NMXJK-1)+F(I,NMXJK+1)+F(I-1,NMXJM1K)+ 2 F(I-1,NMXJP1K)+F(I+1,NMXJM1K)+F(I+1,NMXJP1K)) 201 CONTINUE DO 202 I=3,IMAX+2 S11(I,1) = SMALL S11(I,KMAXP1) = SMALL S12(I,1) = SMALL S12(I,KMAXP1) = SMALL 202 CONTINUE NXK = NJNP+NX-1 NXNMK = NJNP+NXNM-1 DO 20 K=1,KMAXP1 NXK = NXK+1 NXNMK = NXNMK+1 DO 20 I=3,IMAX+2 if (f(i,nxk) < small) f(i,nxk) = s11(i,k) if (f(i,nxnmk) < small) f(i,nxnmk) = s12(i,k) 20 CONTINUE C **** INSERT PERIODIC POINTS K1=NJNP+NX K2=K1+KMAX DO 18 N=1,2 DO 19 I=1,2 DO 19 K=K1,K2 F(I,K)=F(I+IMAX,K) F(I+IMAXP2,K)=F(I+2,K) 19 CONTINUE K1=NJNP+NXNM K2=K1+KMAX 18 CONTINUE RETURN END C