#include "dims.h" SUBROUTINE CMPHOX use cons_module,only: len1,len2,kmax,expz,p0,expzmid,boltz,re implicit none C **** C **** ADVANCES PSI(HOX) BY ONE TIME STEP C **** #include "params.h" #include "fcom.h" #include "vscr.h" #include "index.h" #include "buff.h" #include "crates_const.h" #include "crates_tdep.h" #include "lowbnd.h" #include "phys.h" #include "compcom.h" #include "mwt.h" #include "cmpdat.h" real :: fmin,fmax ! ! Local: real :: PHIHOX(3)=(/0.146, 0.243, 0.162/) real :: small,alfa,xyhox integer :: i,k,NTK,NMSK,NPS1K,NPS2K,NOPK,NPH2OK,NPOHK,NPHO2K, | NPH2K,NPHK,NPHOXK,NPO1K,NPO3K,NPNOK,NPCOK,NPN4SK,NPCH4K, | ibnd,ibndb,njh2otk,njh2olk,njch4tk,nzk C **** C **** NUMBER DENSITY MIXING RATIO OF HOX AT LOWER BOUNDARY C **** C SMALL=1.E-8 SMALL=1.E-12 C **** C **** BOUNDARIES C **** C **** UPPER BOUNDARY CONDITION IS AN UPWARD FLUX OF HOX C **** NTK = NJ+NT+KMAX-1 NZK = NJ+NZ+KMAX NPHOXK = NJ+NPHOX+KMAX-1 NMSK = NJ+NMS+KMAX DO 11 I=1,LEN1 C **** C **** T8 = N*MBAR AT THE UPPER BOUNDARY C **** T8(I) = p0*expz(KMAX)*expzmid*F(I,NMSK)/(boltz*F(I,NTK)) T9(I) = 0.5*(F(I,NPHOXK)+F(I,NPHOXK+1)) T5(I) = (re+F(I,NZK))/(9.555E+4*F(I,NTK)) T6(I) = 2.6183E+3*T9(I)/RMHOX*T8(I)*SQRT(F(I,NTK)) 1 *(1.+T5(I))*EXP(-T5(I)) T7(I) = 2.03E-7*F(I,NTK)*T9(I)/RMHOX*T8(I)*T6(I) C **** REDUCE CHARGE EXCHANGE FLUX TO PREVENT H HOLES C T7(I) = T7(I)*0.5 C T7(I) = T7(I)*0.05 T7(I) = T7(I)*0.01 T4(I) = T6(I)+T7(I) 11 CONTINUE NJH2OTK = NJH2OT-1 NJH2OLK = NJH2OL-1 NJCH4TK = NJCH4T-1 DO 2 K=1,KMAX NJH2OTK = NJH2OTK+1 NJH2OLK = NJH2OLK+1 NJCH4TK = NJCH4TK+1 DO 2 I=1,LEN1 S4(I,K) = 0.5*(F(I,NJH2OLK)+F(I,NJH2OLK+1)) S5(I,K) = 0.5*(F(I,NJH2OTK)+F(I,NJH2OTK+1)) S6(I,K) = S5(I,K)-0.12*S4(I,K) S14(I,K) = 0.5*(XJH2O2(I,K)+XJH2O2(I,K+1)) S15(I,K) = 0.5*(F(I,NJCH4TK)+F(I,NJCH4TK+1)) CC S15(I,K) = 0.5*(XJCH4A(I,K)+XJCH4A(I,K)) 2 CONTINUE C **** C **** SOURCES C **** NTK = NJ+NT-1 NMSK = NJ+NMS-1 NPS1K = NJ+NPS-1 NPS2K = NJ+NPS2-1 NOPK = NJ+NOP-1 NPH2OK = NJ+NPH2O-1 NPOHK = NJ+NPOH-1 NPHO2K = NJ+NPHO2-1 NPH2K = NJ+NPH2-1 NPHK = NJ+NPH-1 NPHOXK = NJ+NPHOX-1 NPO1K = NJ+NPO1-1 NPO3K = NJ+NPO3-1 NPNOK = NJ+NPNO-1 NPCOK = NJ+NPCO-1 NPN4SK = NJ+NPN4S-1 NPCH4K = NJ+NPCH4-1 DO 5 K=1,KMAX NTK = NTK+1 NMSK = NMSK+1 NPS1K = NPS1K+1 NPS2K = NPS2K+1 NOPK = NOPK+1 NPH2OK = NPH2OK+1 NPOHK = NPOHK+1 NPHO2K = NPHO2K+1 NPH2K = NPH2K+1 NPHK = NPHK+1 NPHOXK = NPHOXK+1 NPO1K = NPO1K+1 NPO3K = NPO3K+1 NPNOK = NPNOK+1 NPCOK = NPCOK+1 NPN4SK = NPN4SK+1 NPCH4K = NPCH4K+1 DO 5 I=1,LEN1 C **** C **** S3 = N*MBAR (K+1/2) C **** S3(I,K) = XNMBAR(I,K) S7(I,K) = S3(I,K) S11(I,K) = .5*(F(I,NMSK)+F(I,NMSK+1)) C **** C **** C **** S2 = NUMBER DENSITY PRODUCTION OF HOX*RMHOX(TRUE)/ C **** RMHOX C **** CC S1(I,K) = S3(I,K)**2*((2.*RKM3*F(I,NPH2OK)/RMH2O+ CC 1 2.*RKM4*F(I,NPH2K)/RMH2)*XNO1D(I,K)/RMO1 CC 2 +2.*RKM28(I,K)*F(I,NPO1K)/RMO1*F(I,NPH2K)/RMH2 CC 5 +2.*RKM27(I,K)*F(I,NPO1K)/RMO1*XNH2O2(I,K)/RMH2O2) CC 6 +2.*S6(I,K)*F(I,NPH2OK)/RMH2O*S3(I,K) CC 7 +2.*S14(I,K)*XNH2O2(I,K)/RMH2O2*S3(I,K) CC 8 +PHOXIC(I,K)*F(I,NPH2OK)/RMH2O*S3(I,K) CC 9 +RK11(I,K)*XNHP(I,K)*F(I,NPO1K)/RMO1*S3(I,K) CC 1 +S15(I,K)*F(I,NPCH4K)/RMCH4*S3(I,K)+GAM4 CC 2 *CH3(I,K)*F(I,NPO1K)/RMO1*S3(I,K)+XJCH2OA(I,K) CC 3 *CH2O(I,K)+GAM2A*F(I,NPCH4K)/RMCH4*XNO1D(I,K) CC 4 /RMO1*S3(I,K)*S3(I,K) CC 5 +GAM11(I,K)*CH2O(I,K)*F(I,NPO1K)/RMO1*S3(I,K) CC 6 +GAM9(I,K)*CH3O(I,K)*F(I,NPS1K)*rmassinv_o2*S3(I,K) CC 7 +GAM12(I,K)*CHO(I,K)*F(I,NPS1K)*rmassinv_o2*S3(I,K) CC 8 +GAM13(I,K)*F(I,NPO1K)/RMO1*F(I,NPCH4K)/RMCH4* CC 9 S3(I,K)*S3(I,K)+RK14*F(I,NOPK)*F(I,NPH2K) CC 1 /RMH2*S3(I,K) S1(I,K) = S3(I,K)**2*((2.*RKM3*F(I,NPH2OK)/RMH2O+ 1 2.*RKM4*F(I,NPH2K)/RMH2)*XNO1D(I,K)/RMO1 2 +2.*RKM28(I,K)*F(I,NPO1K)/RMO1*F(I,NPH2K)/RMH2 3 +2.*RKM27(I,K)*F(I,NPO1K)/RMO1*XNH2O2(I,K)/RMH2O2) 4 +2.*S6(I,K)*F(I,NPH2OK)/RMH2O*S3(I,K) 5 +2.*S14(I,K)*XNH2O2(I,K)/RMH2O2*S3(I,K) 6 +PHOXIC(I,K)*F(I,NPH2OK)/RMH2O*S3(I,K) 7 +RK11(I,K)*XNHP(I,K)*F(I,NPO1K)/RMO1*S3(I,K) 8 +S15(I,K)*F(I,NPCH4K)/RMCH4*S3(I,K) 9 +RK14*F(I,NOPK)*F(I,NPH2K)/RMH2*S3(I,K) 1 +(2.*RKM44(I,K)*F(I,NPO1K)/RMO1+2.*RKM45* 2 XNO1D(I,K)/RMO1)*S3(I,K)*F(I,NPCH4K)/RMCH4*S3(I,K) S2(I,K) = (2.*RKM38*RATIO2(I,K)+2.*RKM40(I,K)* 1 RATIO2(I,K)+2.*RKM30(I,K)*RATIO1(I,K)**2+2.*RKM31(I,K) 2 *RATIO1(I,K)*RATIO2(I,K)+2.*RKM41(I,K)*S3(I,K)/S11(I,K) 3 +2.*RKM35(I,K)*RATIO2(I,K)**2)*RATIO3(I,K)**2 CC S3(I,K) = RK12*F(I,NOPK)*RATIO3(I,K) CC 1 +(GAM1(I,K)*F(I,NPCH4K)/RMCH4*S3(I,K)+GAM8(I,K)* CC 2 CH3OOH(I,K)+GAM10*CH2O(I,K))*RATIO1(I,K)* CC 3 RATIO3(I,K)+GAM6(I,K)*CH3O2(I,K)*RATIO2(I,K)* CC 4 RATIO3(I,K) S3(I,K) = RK12*F(I,NOPK)*RATIO3(I,K) ! S1(I,K) = merge(S1(I,K),1.E-20,S1(I,K)-1.E-20>=0.) ! S2(I,K) = merge(S2(I,K),1.E-20,S2(I,K)-1.E-20>=0.) ! S3(I,K) = merge(S3(I,K),1.E-20,S3(I,K)-1.E-20>=0.) if (s1(i,k) < 1.e-20) s1(i,k) = 1.e-20 if (s2(i,k) < 1.e-20) s2(i,k) = 1.e-20 if (s3(i,k) < 1.e-20) s3(i,k) = 1.e-20 5 CONTINUE NMSK = NJ+NMS DO 1 I=1,LEN1 C **** C **** VALUE AT BOTTOM GIVEN BY SPECIFIED NUMBER DENSITY C **** MIXING RATIO XHOXLB C **** T1(I) = 0. T2(I) = 1. C **** T6(I)=(-S3(I,1)+SQRT(S3(I,1)**2+4.*S2(I,1)*S1(I,1))) 1 /(2.*S2(I,1)) T3(I)=-T6(I)*RMTRU(I,1)/S7(I,1) C T3(I) = -SQRT(S1(I,1)/S2(I,1))*RMTRU(I,1)/S7(I,1) 1 CONTINUE ! write(6,"('cmphox: j=',i2,' t3=',/,(6e12.4))") j,t3 ! write(6,"('cmphox: j=',i2,' t4=',/,(6e12.4))") j,t4 IBND = 1 IBNDB = 0 ALFA = -0.38 XYHOX = 1.E-8 CALL HOXCHEM DO 100 I=1,LEN2 S1(I,1)=0. S2(I,1)=0. 100 CONTINUE CALL MINOR2(NPHOX,NPHOXNM,RMHOX,PHIHOX,ALFA,IBND,IBNDB,WHOX, 1XYHOX,NPDHHOX) C **** C **** PARTITION HOX INTO HO2, OH AND H C **** NPHO2K = NJNP+NPHO2-1 NPHOXK = NJNP+NPHOX-1 NPOHK = NJNP+NPOH-1 NPHK = NJNP+NPH-1 DO 7 K=1,KMAX NPHO2K = NPHO2K+1 NPOHK = NPOHK+1 NPHK = NPHK+1 NPHOXK = NPHOXK+1 DO 7 I=1,LEN1 F(I,NPHK) = RATIO3(I,K)*F(I,NPHOXK)*RMH/RMTRU(I,K) F(I,NPOHK) = RATIO1(I,K)*RATIO3(I,K)*F(I,NPHOXK)*RMOH/ 1 RMTRU(I,K) F(I,NPHO2K) = RATIO2(I,K)*RATIO3(I,K)*F(I,NPHOXK)*RMHO2/ 1 RMTRU(I,K) S11(I,K) = SMALL*WHOX(K) ! F(I,NPHK) = merge(F(I,NPHK),S11(I,K),F(I,NPHK)-S11(I,K)>=0.) ! F(I,NPOHK) = merge(F(I,NPOHK),S11(I,K),F(I,NPOHK)-S11(I,K)>=0.) ! F(I,NPHO2K) = merge(F(I,NPHO2K),S11(I,K),F(I,NPHO2K)-S11(I,K)>=0.) if (f(i,nphk) < s11(i,k)) f(i,nphk) = s11(i,k) if (f(i,npohk) < s11(i,k)) f(i,npohk) = s11(i,k) if (f(i,npho2k) < s11(i,k)) f(i,npho2k) = s11(i,k) 7 CONTINUE RETURN END C