#include "dims.h" SUBROUTINE MINOR2(NX,NXNM,RMX,PHIX,ALFAX,IBND,IBNDB,SMALL,XY,NPDH) use cons_module,only: len1,len2,len3,kmax,kmaxp1,imax,imaxp2,avo, | imaxp4,rmass,expz,boltz,p0,expzmid,expzmid_inv,dtx2inv,grav, | dtsmooth,dtsmooth_div2,rmassinv_o2,rmassinv_o,rmassinv_n2 use input_module,only: difhor 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" #include "diffk.h" ! ! Args: integer,intent(in) :: nx,nxnm,ibnd,ibndb,npdh real,intent(in) :: SMALL(ZKMXP),PHIX(3),rmx,alfax,xy ! ! Local: integer :: i,k,k1,k2,n integer :: nps1k,nps2k,ntk,nxjp2k,nxjp1k,nxjnmk,nxjm1k,nxjm2k,nwk, | npdhk,nxnpk,nxjk,nxnmk,nxk,nxmnk,nmxjm1k,nmxjk,nmxjp1k real :: PHI(2,3),tau,t00,alfa12,alfa21,alfax1,alfax2 ! PHI(:,1)=(/0.,0.673/) PHI(:,2)=(/1.35,0./) ; PHI(:,3)=(/1.11,0.769/) TAU=1.86E+3 T00=273. ! SMAL=1.e-20 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_o2+T5(I)*rmassinv_o+(1.-T4(I)-T5(I))* | rmassinv_n2) 2 CONTINUE C **** S12 = MBAR(K+1/2) DO 3 I=1,LEN3 S12(I,1)=1./(F(I,NPS1K)*rmassinv_o2+F(I,NPS2K)*rmassinv_o+ + (1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv_n2) 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*expz(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 NXJNMK = NJ+NXNM DO 11 I=1,LEN3 S9(I,1)=S10(I,1) S10(I,1) = F(I,NXJNMK) 11 CONTINUE C **** S1 = P, S2 = Q, S3 = R, S4= F NWK=NJNP+NW-1 NPDHK = NPDH-1 DO 13 K=1,KMAX NWK=NWK+1 NPDHK = NPDHK+1 DO 13 I=1,LEN1 S1(I,K)=S11(I,K)/dz*(1./dz+.5*S12(I,K))-expz(K)*(expzmid_inv* 1DIFKK(I,K) 2*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* 1DIFKK(I,K+1) 2*T7(I)*(1./dz+.5*S9(I,K+1))-0.25*(F(I,NWK)+F(I,NWK+1))) 2/dz S2(I,K)=-(S11(I,K)/dz*(1./dz-.5*S12(I,K))+S11(I,K+1)/dz*(1./ 1dz+.5*S12(I,K+1)))+expz(K)*((expzmid_inv*DIFKK(I,K)*(1./dz+.5* 2S9(I,K))+expzmid* 3DIFKK(I,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)+F(I,NPDHK)+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)/DIFKK(I,1) T3(I)=T3(I)*grav*RMX/(DIFKK(I,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 **** 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 C **** INSURE NON-NEGATIVE NX NXK = NJNP+NX-1 NXNMK = NJNP+NXNM-1 DO 20 K=1,KMAXP1 NXK = NXK+1 NXNMK = NXNMK+1 DO 20 I=1,LEN1 ! F(I,NXK) = merge(F(I,NXK),SMALL(K)*XY,F(I,NXK)- ! 1 SMALL(K)*XY>=0.) ! F(I,NXNMK) = merge(F(I,NXNMK),SMALL(K)*XY,F(I,NXNMK)- ! 1 SMALL(K)*XY>=0.) if (f(i,nxk) < small(k)*xy) f(i,nxk) = small(k)*xy if (f(i,nxnmk) < small(k)*xy) f(i,nxnmk) = small(k)*xy 20 CONTINUE C **** FILTER NX NXNPK=NJNP+NX CALL FILTER2(NXNPK,KMAXP1) C **** TIME SMOOTHING NXNMK=NJ+NXNM NXK=NJ+NX NXNPK=NJNP+NX NXMNK=NJNP+NXNM DO 17 I=1,LEN3 F(I,NXNMK) = F(I,NXMNK) F(I,NXMNK)=dtsmooth*F(I,NXK)+dtsmooth_div2* | (F(I,NXNMK)+F(I,NXNPK)) 17 CONTINUE RETURN END C