#include "dims.h" SUBROUTINE ELDEN use cons_module,only: len1,len2,len3,kmax,kmaxm1,rmassinv,expz, | boltz,p0 use crates_module,only: rk4,rk5,rk6,rk7,rk8,rk9,rk10,rk16,rk23, | rk26 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_tdep.h" #include "mwt.h" ! ! Local: integer :: i,k integer :: ntk,nmsk,nnpk,nqnpk,nopk,npn2dk,nps1k,nps2k,nqnopk, | nqo2pk,nqn2pk,npnok,npn4sk,nno2k,nn2pk,no2pk,nnopk,nek C **** C **** S12 = N*MBAR =PO*EXP(-Z)*MBAR/(K*T) C **** ! call addfsech('OP_ITP',' ',' ',f(1,nj+nop) ,zimxp,zkmxp,zkmx,j) ! call addfsech('OP_ITC',' ',' ',f(1,njnp+nop ),zimxp,zkmxp,zkmx,j) 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 ! call addfsech('PKT_BARM',' ',' ',s12,zimxp,zkmxp,zkmx,j) C **** C **** CALCULATE N(N+) C **** NNPK=NNP ! NQNPK=NQNP NOPK=NJ+NOP NPN2DK=NJ+NPN2D NPS1K=NJ+NPS NPS2K=NJ+NPS2 DO 3 I=1,LEN2 S11(I,1)=(.5*(F(I,NQNPK)+F(I,NQNPK+1)) | +RK10*F(I,NOPK)*F(I,NPN2DK)*S12(I,1)/ | RMN2D)/(S12(I,1)*((RK6+RK7)*F(I,NPS1K)*rmassinv(1)+ | RK8*F(I,NPS2K)*rmassinv(2))) 3 CONTINUE DO 4 I=1,LEN2 F(I,NNPK)=S11(I,1) 4 CONTINUE ! call addfsech('NPLUS',' ',' ',f(1,nnp),zimxp,zkmxp,zkmx,j) C **** C **** CALCULATE A, B, C, D, E, (F+G) AND H FOR QUARTIC C **** NQNOPK=NQNOP NQO2PK=NQO2P NQN2PK=NQN2P NOPK=NJNP+NOP NNPK=NNP NPS1K=NJ+NPS NPS2K=NJ+NPS2 NPNOK=NJ+NPNO NPN4SK=NJ+NPN4S NNO2K=NNO2 DO 5 I=1,LEN2 C **** C **** S10 = A = QI(NO+)+K2*N(O+)*N(N2)+K7*N(N+)*N(02)+B9*N(NO) C **** S10(I,1)=.5*(F(I,NQNOPK)+F(I,NQNOPK+1))+ | S12(I,1)*(RK2(I,1)*F(I,NOPK)*(1.-F(I,NPS1K)- | F(I,NPS2K))*rmassinv(3)+RK7*F(I,NNPK)*F(I,NPS1K)* | rmassinv(1)+.5*(BETA9(I,1)+BETA9(I+LEN1,1))*F(I,NPNOK)/RMNO) C **** C **** S9 = B = QI(O2+)+K1*N(O+)*N(O2)+K6*N(N+)*N(02) C **** S9(I,1)=.5*(F(I,NQO2PK)+F(I,NQO2PK+1))+ | S12(I,1)*(RK1(I,1)*F(I,NOPK)+RK6* | F(I,NNPK))*F(I,NPS1K)*rmassinv(1)+RK26* | XIOP2D(I,1)*F(I,NPS1K)*rmassinv(1) 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)/RMNO) C **** C **** S7 = D = QI(N2+) C **** S7(I,1)=.5*(F(I,NQN2PK)+F(I,NQN2PK+1))+(RK16*XIOP2P(I,1) 1 +RK23*XIOP2D(I,1))*(1.-F(I,NPS1K)-F(I,NPS2K)) 2 *rmassinv(3) C **** C **** S6 = E =K3*N(O)+K9*N(O2) C **** S6(I,1)=S12(I,1)*(RK3(I,1)*F(I,NPS2K)*rmassinv(2)+RK9* | F(I,NPS1K)*rmassinv(1)) 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(1) 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)* | (S8(I,1)+S4(I,1)))+RA2(I,1)*(S6(I,1)*(S10(I,1)+S7(I,1))-S4(I,1)* | S7(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)* | S8(I,1))*S5(I,1)-RA2(I,1)*S7(I,1)-RA3(I,1)*S9(I,1))-RA2(I,1)* | RA3(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)* | RA3(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 ! call addfsech('A_COEF',' ',' ',s10,zimxp,zkmxp,zkmx,j) ! call addfsech('B_COEF',' ',' ',s9 ,zimxp,zkmxp,zkmx,j) ! call addfsech('C_COEF',' ',' ',s8 ,zimxp,zkmxp,zkmx,j) ! call addfsech('D_COEF',' ',' ',s7 ,zimxp,zkmxp,zkmx,j) ! call addfsech('E_COEF',' ',' ',s6 ,zimxp,zkmxp,zkmx,j) ! call addfsech('FG_COEF',' ',' ',s5 ,zimxp,zkmxp,zkmx,j) ! call addfsech('H_COEF',' ',' ',s4 ,zimxp,zkmxp,zkmx,j) ! call addfsech('A0' ,' ',' ',s15,zimxp,zkmxp,zkmx,j) ! call addfsech('A1' ,' ',' ',s14,zimxp,zkmxp,zkmx,j) ! call addfsech('A2' ,' ',' ',s13,zimxp,zkmxp,zkmx,j) ! call addfsech('A3' ,' ',' ',s12,zimxp,zkmxp,zkmx,j) ! call addfsech('A4' ,' ',' ',s11,zimxp,zkmxp,zkmx,j) C **** C **** SOLVE QUARTIC ( S5 = N(E) ) C **** ! ! 2/00: changed args to vquart (s.a. vquart.f): ! old call: CALL VQUART(S15,S5,S3,S2,S1,LEN3,1,LEN2,j) CALL VQUART(S15,S5,S3,S2,S1,LEN3,1,LEN2,j) ! ! s15,s14,s13,s12,s11 are coefficients for the quartic solver. ! Positive roots of quartic are returned in s5: ! ! call vquart(s15,s14,s13,s12,s11, s5,imaxp4,kmax,kmaxp1,j) C **** C **** ENSURE NON NEGATIVE N(E) C **** C **** C **** CALCULATE N(N2+), N(02+) AND N(NO+) C **** ! call addfsech('ROOT' ,' ',' ',s5,zimxp,zkmxp,zkmx,j) NN2PK=NN2P NO2PK=NJNP+NO2P NNOPK=NNOP DO 7 I=1,LEN2 c 5/7/98: ensure non-negative Ne: C S5(I,1) = merge(S5(I,1),3.3E+3,S5(I,1)-3.3E+3>=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 C S5(I,1) = merge(S5(I,1),2.8E+3,S5(I,1)-2.8E+3>=0.) C S5(I,1) = merge(S5(I,1),2.5E+3,S5(I,1)-2.5E+3>=0.) 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)) 7 CONTINUE call addfsech('NPLUSb',' ',' ' ,f(1,nnp),zimxp,zkmxp,zkmx,j) call addfsech('N2P_ELD',' ',' ',f(1,nn2p),zimxp,zkmxp,zkmx,j) call addfsech('O2P_ELD',' ',' ',f(1,njnp+no2p),zimxp,zkmxp, | zkmx,j) call addfsech('NOP_ELD',' ',' ',f(1,nnop),zimxp,zkmxp,zkmx,j) 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)) C F(I,NEK)= 1.5*S5(I,1)-0.5*S5(I,2) C F(I,NEK+KMAX)=1.5*S5(I,KMAX)-0.5*S5(I,KMAXM1) 9 CONTINUE ! call addfsech('NE_ELDEN' ,' ',' ',f(1,njnp+ne), ! | zimxp,zkmxp,zkmx,j) RETURN END C