#include "dims.h" SUBROUTINE ELDEN use cons_module,only: len1,len2,len3,kmax,kmaxm1,imaxp4,rmass, | expz,p0,boltz,rmassinv_o2,rmassinv_o,rmassinv_n2 implicit none C **** C **** CALCULATES ELECTRON DENSITY AND ION NUMBER DENSITIES C **** #include "params.h" #include "fcom.h" #include "vscr.h" #include "index.h" #include "buff.h" #include "phys.h" #include "crates_const.h" #include "crates_tdep.h" #include "mwt.h" ! ! Local: integer :: i,k,ntk,nmsk,nqnopk,nqo2pk,nqn2pk,nopk,nnpk,nps1k, | nps2k,npo1k,npnok,npn4sk,npco2k,nn2pk,no2pk,nnopk,nek C **** C **** S12 = N*MBAR =PO*EXP(-Z)*MBAR/(K*T) C **** NTK=NJ+NT-1 NMSK=NJ+NMS-1 DO 2 K=1,KMAX NTK=NTK+1 NMSK=NMSK+1 DO 2 I=1,LEN1 S12(I,K)=p0*expz(K)*.5*(F(I,NMSK)+F(I,NMSK+1))/(boltz*F(I,NTK)) 2 CONTINUE C **** C **** CALCULATE A, B, C, D, E, (F+G) AND H FOR QUARTIC C **** NQNOPK=NQNOP NQO2PK=NQO2P NQN2PK=NQN2P NOPK=NJ+NOP NNPK=NNP NPS1K=NJ+NPS NPS2K=NJ+NPS2 NPO1K=NJ+NPO1 NPNOK=NJ+NPNO NPN4SK=NJ+NPN4S NPCO2K=NJ+NPCO2 DO 5 I=1,LEN2 C **** C **** S10 = A = QI(NO+)+K2*N(O+)*N(N2)+K7*N(N+)*N(02) C **** S10(I,1)=.5*(F(I,NQNOPK)+F(I,NQNOPK+1))+ 1S12(I,1)*(RK2(I,1)*F(I,NOPK)*(1.-F(I,NPS1K)- 1F(I,NPS2K))*rmassinv_n2+RK7*F(I,NNPK)*F(I,NPS1K)*rmassinv_o2) C **** C **** S9 = B = QI(O2+)+K1*N(O+)*N(O2)+K6*N(N+)*N(02) C +K13*N(CO2)*N(O+)+K23*N(O2)*N(O+(2P)) C **** +K27*N(O+(2D))*N(O2) C **** S9(I,1)=.5*(F(I,NQO2PK)+F(I,NQO2PK+1))+ 1S12(I,1)*((RK1(I,1)*F(I,NOPK)+RK6* 2F(I,NNPK))*F(I,NPS1K)*rmassinv_o2+RK27* 3XIOP2D(I,1)*F(I,NPS1K)*rmassinv_o2+RK13*F(I,NOPK)* 4F(I,NPCO2K)/RMCO2+RK23*XIOP2P(I,1)*F(I,NPS1K)*rmassinv_o2) C **** C **** S8 = C = K4*N(N4S)+K5*N(NO) C **** S8(I,1)=S12(I,1)*(RK4*F(I,NPN4SK)/RMN4S+RK5*F(I,NPNOK)/ 1RMNO) C **** C **** S7 = D = QI(N2+)+K16*N(N2)*N(O+(2P))+K24*N(N2)*N(O+(2D)) C **** S7(I,1)=.5*(F(I,NQN2PK)+F(I,NQN2PK+1))+(RK16*XIOP2P(I,1) 1 +RK24*XIOP2D(I,1))*(1.-F(I,NPS1K)-F(I,NPS2K)) 2 *rmassinv_n2*S12(I,1) C **** C **** S6 = E =K3*N(O)+K9*N(O2) C **** S6(I,1)=S12(I,1)*(RK3(I,1)*F(I,NPO1K)*rmassinv_o+RK9*F(I,NPS1K) 1*rmassinv_o2) C **** C **** S5 = F+G = N(O+)+N(N+) C **** S5(I,1)=F(I,NOPK)+F(I,NNPK) C **** C **** S4 = H = K9*N(02) C **** S4(I,1)=S12(I,1)*RK9*F(I,NPS1K)*rmassinv_o2 C **** C **** ! S10(I,1) = merge(S10(I,1),1.E-20,S10(I,1)-1.E-20>=0.) ! S9(I,1) = merge(S9(I,1),1.E-20,S9(I,1)-1.E-20>=0.) ! S8(I,1) = merge(S8(I,1),1.E-20,S8(I,1)-1.E-20>=0.) ! S7(I,1) = merge(S7(I,1),1.E-20,S7(I,1)-1.E-20>=0.) ! S6(I,1) = merge(S6(I,1),1.E-20,S6(I,1)-1.E-20>=0.) ! S5(I,1) = merge(S5(I,1),1.E-20,S5(I,1)-1.E-20>=0.) ! S4(I,1) = merge(S4(I,1),1.E-20,S4(I,1)-1.E-20>=0.) ! if (s10(i,1) < 1.e-20) s10(i,1) = 1.e-20 if (s9(i,1) < 1.e-20) s9(i,1) = 1.e-20 if (s8(i,1) < 1.e-20) s8(i,1) = 1.e-20 if (s7(i,1) < 1.e-20) s7(i,1) = 1.e-20 if (s6(i,1) < 1.e-20) s6(i,1) = 1.e-20 if (s5(i,1) < 1.e-20) s5(i,1) = 1.e-20 if (s4(i,1) < 1.e-20) s4(i,1) = 1.e-20 C **** C **** CALCULATE QUARTIC COEFFICIENTS C **** C **** S15 = A0 C **** S15(I,1)=-S6(I,1)*S8(I,1)*(S10(I,1)+S9(I,1)+S7(I,1)) C **** C **** S14 = A1/4 C **** S14(I,1)=-(RA1(I,1)*(S6(I,1)*(S8(I,1)*S5(I,1)+S9(I,1))+S7(I,1)* 1(S8(I,1)+S4(I,1)))+RA2(I,1)*(S6(I,1)*(S10(I,1)+S7(I,1))-S4(I,1)* 2S7(I,1))+RA3(I,1)*S8(I,1)*(S10(I,1)+S9(I,1)))/4. C **** C **** S13 = A2/6 C **** S13(I,1)=(RA1(I,1)*(S6(I,1)*S8(I,1)-(RA2(I,1)*S6(I,1)+RA3(I,1)* 1S8(I,1))*S5(I,1)-RA2(I,1)*S7(I,1)-RA3(I,1)*S9(I,1))-RA2(I,1)* 2RA3(I,1)*S10(I,1))/6. C **** C **** S12 = A3/4 C **** S12(I,1)=(RA1(I,1)*(RA2(I,1)*S6(I,1)+RA3(I,1)*S8(I,1)-RA2(I,1)* 1RA3(I,1)*S5(I,1)))/4. C **** C **** S11 = A4 C **** S11(I,1)=RA1(I,1)*RA2(I,1)*RA3(I,1) 5 CONTINUE C **** C **** SOLVE QUARTIC ( S5 = N(E) ) C **** ! Vquart uses 1st arg to reference s15,s14,s13,s12,s11 in that order. ! (s1,s2,s3 are workspace) ! CALL VQUART(S15,S5,S3,S2,S1,LEN3,1,LEN2) C **** C **** ENSURE NON NEGATIVE N(E) C **** C **** C **** CALCULATE N(N2+), N(02+) AND N(NO+) C **** NN2PK=NN2P NO2PK=NJNP+NO2P NNOPK=NNOP DO 7 I=1,LEN2 c 1/21/97: ensure non-negative Ne: ! 11/30/98 change minimum S5 according to modsrc.kibo: ! S5(I,1) = merge(S5(I,1),1.E+2,S5(I,1)-1.E+2>=0.) ! S5(I,1) = merge(S5(I,1),3.1E+3,S5(I,1)-3.1E+3>=0.) if (s5(i,1) < 3.1e+3) s5(i,1) = 3.1e+3 F(I,NN2PK)=S7(I,1)/(S6(I,1)+RA3(I,1)*S5(I,1)) F(I,NO2PK)=(S9(I,1)+S4(I,1)*S7(I,1)/(S6(I,1)+RA3(I,1)*S5(I,1)))/ | (S8(I,1)+RA2(I,1)*S5(I,1)) F(I,NNOPK)=(S10(I,1)+S7(I,1)*(S6(I,1)-S4(I,1))/(S6(I,1)+RA3(I,1)* | S5(I,1))+S8(I,1)*(S9(I,1)+S4(I,1)*S7(I,1)/(S6(I,1)+RA3(I,1)* | S5(I,1)))/(S8(I,1)+RA2(I,1)*S5(I,1)))/(RA1(I,1)*S5(I,1)) ! F(I,NO2PK) = merge(F(I,NO2PK),1.E-20,F(I,NO2PK)-1.E-20>=0.) ! F(I,NN2PK) = merge(F(I,NN2PK),1.E-20,F(I,NN2PK)-1.E-20>=0.) ! F(I,NNOPK) = merge(F(I,NNOPK),1.E-20,F(I,NNOPK)-1.E-20>=0.) if (f(i,no2pk) < 1.e-20) f(i,no2pk) = 1.e-20 if (f(i,nn2pk) < 1.e-20) f(i,nn2pk) = 1.e-20 if (f(i,nnopk) < 1.e-20) f(i,nnopk) = 1.e-20 7 CONTINUE ! ! 6/26/98 btf: Define top level O2+ for utask/mtask: ! no2pk = njnp+no2p+kmax do i=1,imaxp4 f(i,no2pk) = f(i,no2pk-1) enddo C **** C **** C **** PLACE N(E) IN SLOT IN F-ARRAY C **** NEK=NJNP+NE DO 8 I=1,LEN2-LEN1 F(I,NEK+1)=SQRT(S5(I,1)*S5(I,2)) 8 CONTINUE DO 9 I=1,LEN1 F(I,NEK)=SQRT(S5(I,1)**3/S5(I,2)) F(I,NEK+KMAX)=SQRT(S5(I,KMAX)**3/S5(I,KMAXM1)) 9 CONTINUE C call addfsech('NOP' ,' ',' ',F(1,NNOP) ,zimxp,zkmxp,zkmx,j) C RETURN END C