#include "dims.h" SUBROUTINE MINOR(NX,NXNM,RMX,PHIX,ALFAX,IBND,IBNDB,idebug) 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,idebug 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,idebug) ! if (idebug > 0) ! | call addfsech('HADVEC',' ',' ',s3,zimxp,zkmxp,zkmx,j) DO 24 I=1,LEN3 S13(I,1)=S3(I,1) ! hadvec 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 if (idebug > 0) | call addfsech('HADVEC',' ',' ',s13,zimxp,zkmxp,zkmx,j) 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 if (idebug > 0) then call addfsech('SLOSS',' ',' ',s14,zimxp,zkmxp,zkmx,j) call addfsech('SPROD',' ',' ',s15,zimxp,zkmxp,zkmx,j) endif 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)) ! if (idebug > 0) ! | write(6,"('minor: lat=',i2,' i=',i2,' bo2=',e12.4,' bo1=', ! | e12.4,' xmbari=',e12.4)") j,i,t4(i),t5(i),t6(i) 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)) ! xmbar_k S10(I,1)=(S12(I,1)-T6(I))/dz ! dmdz S9(I,1)=.5*(T5(I)+F(I,NPS2K)) ! pso1 S8(I,1)=.5*(T4(I)+F(I,NPS1K)) ! pso2 S7(I,1)=(F(I,NPS2K)-T5(I))/dz ! do1dz S6(I,1)=(F(I,NPS1K)-T4(I))/dz ! do2dz 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 if (idebug > 0) then call addfsech('XMBAR_KH',' ',' ',s12,zimxp,zkmxp,zkmx,j) call addfsech('XMBAR_K' ,' ',' ',s11,zimxp,zkmxp,zkmx,j) call addfsech('DMDZ0' ,' ',' ',s10,zimxp,zkmxp,zkmx,j) call addfsech('PSO1' ,' ',' ',s9 ,zimxp,zkmxp,zkmx,j) call addfsech('PSO2' ,' ',' ',s8 ,zimxp,zkmxp,zkmx,j) call addfsech('DO1DZ' ,' ',' ',s7 ,zimxp,zkmxp,zkmx,j) call addfsech('DO2DZ' ,' ',' ',s6 ,zimxp,zkmxp,zkmx,j) call addfsech('TNI' ,' ',' ',s5 ,zimxp,zkmxp,zkmx,j) endif 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* ! s0prod 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)) ! alfa11 S2(I,1)=ALFA12*S8(I,1) ! alfa12 S3(I,1)=ALFA21*S9(I,1) ! alfa21 S4(I,1)=-(PHI(2,3)+ALFA21*S8(I,1)) ! alfa22 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 if (idebug > 0) then call addfsech('S0PROD',' ',' ',s15,zimxp,zkmxp,zkmx ,j) call addfsech('ALFA11',' ',' ',s1 ,zimxp,zkmxp,zkmxp,j) call addfsech('ALFA12',' ',' ',s2 ,zimxp,zkmxp,zkmxp,j) call addfsech('ALFA21',' ',' ',s3 ,zimxp,zkmxp,zkmxp,j) call addfsech('ALFA22',' ',' ',s4 ,zimxp,zkmxp,zkmxp,j) call addfsech('EX' ,' ',' ',s12,zimxp,zkmxp,zkmxp,j) call addfsech('DMDZ1' ,' ',' ',s10,zimxp,zkmxp,zkmxp,j) call addfsech('AX0' ,' ',' ',s11,zimxp,zkmxp,zkmxp,j) endif 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 if (idebug > 0) then call addfsech('THDIFF0',' ',' ',s12,zimxp,zkmxp,zkmxp,j) endif 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 ! if (idebug > 0) then ! do i=1,len1 ! write(6,"('minor: lat=',i2,' i=',i2,' difk=',/,(6e12.4))") ! | j,i,difk(:) ! enddo ! write(6,"('minor: lat=',i2,' dfactor=',/,(6e12.4))") ! | j,t7(:) ! endif 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 if (idebug > 0) then call addfsech('AX1' ,' ',' ',s11,zimxp,zkmxp,zkmxp,j) call addfsech('DMDZ2' ,' ',' ',s9 ,zimxp,zkmxp,zkmxp,j) call addfsech('THDIFF1',' ',' ',s12,zimxp,zkmxp,zkmxp,j) call addfsech('FSMOOTH',' ',' ',s10,zimxp,zkmxp,zkmx,j) call addfsech('W_OMEGA',' ',' ',f(1,njnp+nw),zimxp,zkmxp,zkmx,j) ! ! call addfsech('P_COEF' ,' ',' ',s1 ,zimxp,zkmxp,zkmx,j) ! call addfsech('R_COEF' ,' ',' ',s3 ,zimxp,zkmxp,zkmx,j) ! call addfsech('Q_COEF' ,' ',' ',s2 ,zimxp,zkmxp,zkmx,j) ! call addfsech('F_RHS' ,' ',' ',s4 ,zimxp,zkmxp,zkmx,j) endif 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 if (idebug > 0) then call addfsech('P_COEF' ,' ',' ',s1 ,zimxp,zkmxp,zkmxp,j) call addfsech('R_COEF' ,' ',' ',s3 ,zimxp,zkmxp,zkmxp,j) call addfsech('Q_COEF' ,' ',' ',s2 ,zimxp,zkmxp,zkmxp,j) call addfsech('F_RHS' ,' ',' ',s4 ,zimxp,zkmxp,zkmxp,j) endif 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,j,idebug) if (idebug > 0) then call addfsech('MNR_SOLV' ,' ',' ',f(1,njnp+nx),zimxp,zkmxp, | zkmxp,j) endif C **** FILTER NX IF(KUTT(J).LT.IMAXH)THEN NXNPK=NJNP+NX CALL FILTER2(NXNPK,KMAXP1,KUTT(J),0) ENDIF if (idebug > 0) then call addfsech('MNR_FILT' ,' ',' ',f(1,njnp+nx),zimxp,zkmxp, | zkmxp,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 if (idebug > 0) then call addfsech('MNR_SMOO' ,' ',' ',f(1,njnp+nxnm),zimxp,zkmxp, | zkmxp,j) endif 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 if (idebug > 0) then call addfsech('MNR_OUT' ,' ',' ',f(1,njnp+nx),zimxp,zkmxp, | zkmxp,j) call addfsech('MNR_TM1' ,' ',' ',f(1,njnp+nxnm),zimxp,zkmxp, | zkmxp,j) endif RETURN END C