#include "dims.h" ! am_02/02: option to change electron-neutral collision frequency SUBROUTINE LAMDAS use cons_module,only: len1,len2,kmax,kmaxm1,kmaxp1,imax,imaxp4, | rmassinv,expz,gask,grav,avo,p0,boltz,dipmin implicit none C **** C **** COMPUTE ION DRAG COEFFICIENTS C **** C **** MAJOR SPECIES AND IONS ARE NUMBERED THUS C **** O2 - 1 O2+ - 1 C **** O - 2 O+ - 2 C **** N2 - 3 N2+ - 3 C **** #include "params.h" #include "fcom.h" #include "index.h" #include "buff.h" #include "phys.h" #include "dynamo.h" #include "trgm.h" #include "extrafields.h" real :: bxm,bx,bxp,by,byp,bz,bzp,bmod,bmodp COMMON/MAGFLD/BXM(ZIMXP,2),BX(ZIMXP,ZJMX),BXP(ZIMXP,4), 1 BY(ZIMXP,ZJMX),BYP(ZIMXP,4),BZ(ZIMXP,ZJMX),BZP(ZIMXP,4), 2 BMOD(ZIMXP,ZJMX),BMODP(ZIMXP,2) real :: rm1,rm2,rm3,v1,v2,v3,v4,v5,v6,t1,t2,t3,t4,t5,t6,t7 ! COMMON/VSCR/RM1(ZIMXP,ZKMXP,3,1),RM2(ZIMXP,ZKMXP,3,1), 1RM3(ZIMXP,ZKMXP,3,1),V1(ZIMXP,ZKMXP,1),V2(ZIMXP,ZKMXP,1), 2V3(ZIMXP,ZKMXP,1),V4(ZIMXP,ZKMXP,1),V5(ZIMXP,ZKMXP,1), 3V6(ZIMXP,ZKMXP,1),T1(ZIMXP),T2(ZIMXP),T3(ZIMXP),T4(ZIMXP), 4T5(ZIMXP),T6(ZIMXP),T7(ZIMXP) !$OMP THREADPRIVATE (/vscr/) !DIR$ TASKCOMMON vscr ! ! Local: real :: RMS(3),qe,v4tmp integer :: i,nps1k,nps2k,ntik,ntk,l,m,no2pk,nopk,nnopk,nmsk,k, | ntek,nlxxk,nlyyk,nlxyk,nzk,nuk,nvk,nwk,nuik,nvik real :: tm1(zimxp),tm2(zimxp),tm3(zimxp),sndip(zimxp) real :: sumfaus, sumfavs, sumfwus, sumfwvs, sumfaun, sumfavn, | sumfwun, sumfwvn, qwind0, qamie0, wtot0, fwindu0, fwindv0, | famieu0, famiev0, work0 ! DATA RMS/32.,16.,30./,QE/1.602E-19/ C **** C **** T1 = BGAUSS C **** T2 = 9.6489E3 * BGAUS ( FOR OMEGA(I) ) C **** T3 = 1.7588028E7 * BGAUSS ( FOR OMEGA(E) ) C **** T4 = SIN(DIP) C **** DO 1 I=1,LEN1 T1(I)=BMOD(I,J) T2(I)=9.6489E3*T1(I) T3(I)=1.7588028E7*T1(I) ! DIPMIN is in trgm.h, initialized in types.f: if (abs(DIPMAG(I,J)) >= dipmin) then t4(i) = sin(dipmag(i,j)) else t4(i) = sin(dipmin) endif T5(I) = T4(I)**2 1 CONTINUE C **** C **** V1 = PSI(O2)/M(O2), V2 = PSI(O)/M(O), V1 = PSI2)/M(N2) C **** AT (K+1/2) C **** NPS1K=NJ+NPS NPS2K=NJ+NPS2 NTIK=NJ+NTI NTK=NJ+NT DO 2 I=1,LEN2 V1(I,1,1)=F(I,NPS1K)*rmassinv(1) V2(I,1,1)=F(I,NPS2K)*rmassinv(2) V3(I,1,1)=(1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv(3) if (v3(i,1,1) < 1.e-20) v3(i,1,1) = 1.e-20 C **** C **** V6 = TR = 0.5*(TI+TN) ( K+1/2 ) C **** V6(I,1,1)=0.5*(F(I,NTIK)+F(I,NTK)) 2 CONTINUE C **** C **** RM1 = NEW(L,M) C **** C **** GENERATE NUMERICAL FACTORS CALL NEW(V6,RM1) C **** MULTIPLY BY MAJOR SPECIES NUMBER DENSITIES DO 3 L=1,3 DO 3 M=1,3 DO 3 I=1,LEN2 RM1(I,1,L,M)=RM1(I,1,L,M)*V1(I,1,M) if (RM1(I,1,L,M) .lt. 0.) | print *,'LAMDAS: negative RM1 for J,I,M,L = ', | J,I,M,L,RM1(I,1,L,M),V6(I,1,1),F(I,NTIK),F(I,NTK) 3 CONTINUE C **** C **** V1(L) = SUM(M=1,3) RM1(L,M) ( K+1/2 ) C **** DO 4 L=1,3 DO 5 I=1,LEN2 V1(I,1,L)=0. 5 CONTINUE DO 4 M=1,3 DO 4 I=1,LEN2 V1(I,1,L)=V1(I,1,L)+RM1(I,1,L,M) 4 CONTINUE C **** C **** RM1(1) = N(O2+), RM1(2) = N(O+), RM1(1) = N(NO+) C **** ( AT (K+1/2) ) C **** NO2PK=NJ+NO2P NOPK=NJ+NOP NNOPK=NNOP DO 6 I=1,LEN2 RM1(I,1,1,1)=F(I,NO2PK) if (RM1(I,1,1,1) .lt. 0.) | print *,'LAMDAS: negative O2+',I,J,F(I,NO2PK) RM1(I,1,2,1)=F(I,NOPK) if (RM1(I,1,2,1) .lt. 0.) | print *,'LAMDAS: negative O+',I,J,F(I,NOPK) RM1(I,1,3,1)=F(I,NNOPK) if (RM1(I,1,3,1) .lt. 0.) | print *,'LAMDAS: negative NO+',I,J,F(I,NNOPK) 6 CONTINUE C **** C **** V6 = N * MBAR ( K+1/2 ) C **** NTK=NJ+NT-1 NMSK=NJ+NMS-1 DO 7 K=1,KMAX NTK=NTK+1 NMSK=NMSK+1 DO 7 I=1,LEN1 V6(I,K,1)=p0*expz(K)*.5*(F(I,NMSK)+F(I,NMSK+1))/(boltz*F(I,NTK)) 7 CONTINUE C **** C **** V1(L) = V1(L)*(N*MBAR)/OMEGAI(L) ( L=1,3 ) ( K=1/2 ) C **** DO 8 L=1,3 DO 8 K=1,KMAX DO 8 I=1,LEN1 V1(I,K,L)=V1(I,K,L)*V6(I,K,1)*RMS(L)/T2(I) 8 CONTINUE C **** C **** V4 = NU/OMEGA (E) C **** NPS1K=NJ+NPS NPS2K=NJ+NPS2 NTEK=NJ+NTE DO 9 I=1,LEN2 V5(I,1,1)=SQRT(F(I,NTEK)) if ((1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv(3) >= 1.E-20) then v4tmp =(1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv(3) else v4tmp = 1.e-20 endif ! am 2001-10-2 modifying the electron-neutral collision frequency ! paper: Atmos. Electrodynamics ! V4(I,1,1)=V6(I,1,1)*(7.2E-15*v4tmp*(F(I,NTEK)/300.)**0.95 ! | +5.2E-15*F(I,NPS1K)*(F(I,NTEK)/300.)**0.79*rmassinv(1) ! | +1.9E-15*F(I,NPS2K)*rmassinv(2)*(F(I,NTEK)/300.)**0.85) ! | *1E6 ! convertion from m^3 to cm^3 ! ! original Banks & Kockards V4(I,1,1)=V6(I,1,1)*(2.33E-11*v4tmp* | F(I,NTEK)*(1.-1.21E-4*F(I,NTEK))+1.82E-10*F(I,NPS1K)* | rmassinv(1)*V5(I,1,1)*(1.+3.6E-2*V5(I,1,1))+8.9E-11* | F(I,NPS2K)*rmassinv(2)*V5(I,1,1)*(1.+5.7E-4*F(I,NTEK))) ! end am 2001-10-2 9 CONTINUE C **** C **** V4 = NU(EN)/OMEGAE C **** DO 10 K=1,KMAX DO 10 I=1,LEN1 ! am 2001-10-17 multiply collision frequency by factor of 4. ! V4(I,K,1)=4.*V4(I,K,1)/T3(I) V4(I,K,1)=V4(I,K,1)/T3(I) 10 CONTINUE C **** C **** RM2 = PEDERSEN FACTORS C **** RM3 = HALL FACTORS C **** DO 11 L=1,3 DO 11 I=1,LEN2 RM2(I,1,L,1)=V4(I,1,1)/(1.+V4(I,1,1)**2)+V1(I,1,L)/(1.+V1(I,1,L) 1**2) RM3(I,1,L,1)=1./(1.+V4(I,1,1)**2)-1./(1.+V1(I,1,L)**2) 11 CONTINUE C **** C **** V1 = SUM(L=1,3) PEDERSEN FACTORS * ION DENSITIES C **** V2 = SUM(L=1,3) HALL FACTORS * ION DENSITIES C **** DO 12 I=1,LEN2 V1(I,1,1)=0. V2(I,1,1)=0. 12 CONTINUE DO 13 L=1,3 DO 13 I=1,LEN2 V1(I,1,1)=V1(I,1,1)+RM2(I,1,L,1)*RM1(I,1,L,1) V2(I,1,1)=V2(I,1,1)+RM3(I,1,L,1)*RM1(I,1,L,1) 13 CONTINUE C **** C **** SET LXX, LYY, LXY C **** DO 14 K=1,KMAX DO 14 I=1,LEN1 V3(I,K,1)=V1(I,K,1)*T1(I)*1.E-1*QE*avo/V6(I,K,1) V4(I,K,1)=V3(I,K,1)*T5(I) V5(I,K,1)=V2(I,K,1)*T1(I)*1.E-1*QE*avo/V6(I,K,1) 14 CONTINUE NLXXK=NLXX+1 NLYYK=NLYY+1 NLXYK=NLXY+1 DO 15 I=1,LEN2-LEN1 if (V3(I,1,1)*V3(I,2,1) .lt. 0.) print *, | 'lamdas: I,J,V3(I,1,1),V3(I,2,1) = ',I,J,V3(I,1,1),V3(I,2,1) if (V3(I,1,1) .lt. 0.) V3(I,1,1) = 0. if (V3(I,2,1) .lt. 0.) V3(I,2,1) = 0. F(I,NLXXK)=SQRT(V3(I,1,1)*V3(I,2,1)) if (V4(I,1,1)*V4(I,2,1) .lt. 0.) print *, | 'lamdas: I,J,V4(I,1,1),V4(I,2,1) = ',I,J,V4(I,1,1),V4(I,2,1) if (V4(I,1,1) .lt. 0.) V4(I,1,1) = 0. if (V4(I,2,1) .lt. 0.) V4(I,2,1) = 0. F(I,NLYYK)=SQRT(V4(I,1,1)*V4(I,2,1)) if (V5(I,1,1)*V5(I,2,1) .lt. 0.) print *, | 'lamdas: I,J,V5(I,1,1),V5(I,2,1) = ',I,J,V5(I,1,1),V5(I,2,1) if (V5(I,1,1) .lt. 0.) V5(I,1,1) = 0. if (V5(I,2,1) .lt. 0.) V5(I,2,1) = 0. F(I,NLXYK)=SQRT(V5(I,1,1)*V5(I,2,1)) 15 CONTINUE C **** C **** BOUNDARIES C **** NLXXK=NLXX NLYYK=NLYY NLXYK=NLXY DO 16 I=1,LEN1 F(I,NLXXK)=SQRT(V3(I,1,1)**3/V3(I,2,1)) F(I,NLYYK)=SQRT(V4(I,1,1)**3/V4(I,2,1)) F(I,NLXYK)=SQRT(V5(I,1,1)**3/V5(I,2,1)) F(I,NLXXK+KMAX)=SQRT(V3(I,KMAX,1)**3/V3(I,KMAXM1,1)) F(I,NLYYK+KMAX)=SQRT(V4(I,KMAX,1)**3/V4(I,KMAXM1,1)) F(I,NLXYK+KMAX)=SQRT(V5(I,KMAX,1)**3/V5(I,KMAXM1,1)) C **** C **** C **** LEVELS 1 AND KMAXP1 C F(I,NLXXK)= 1.5*V3(I,1,1)-0.5*V3(I,2,1) C F(I,NLYYK)= 1.5*V4(I,1,1)-0.5*V4(I,2,1) C F(I,NLXYK)= 1.5*V5(I,1,1)-0.5*V5(I,2,1) C F(I,NLXXK+KMAX)=1.5*V3(I,KMAX,1)-0.5*V3(I,KMAXM1,1) C F(I,NLYYK+KMAX)=1.5*V4(I,KMAX,1)-0.5*V4(I,KMAXM1,1) C F(I,NLXYK+KMAX)=1.5*V5(I,KMAX,1)-0.5*V5(I,KMAXM1,1) C **** C **** 16 CONTINUE NLXYK = NLXY-1 DO 19 K =1,KMAXP1 NLXYK = NLXYK+1 DO 19 I = 1,IMAXP4 F(I,NLXYK) = F(I,NLXYK)*T4(I) 19 CONTINUE C **** C **** COPY SIGMA1, SIGMA2, Z, H, U, V AND W TO COMMON/DYNAMO/ C **** NZK = NJ+NZ-1 NTK = NJ+NT-1 NMSK = NJ+NMS-1 NUK = NJ+NU-1 NVK = NJ+NV-1 NWK = NJ+NW-1 DO 17 K = 1,KMAX NZK = NZK+1 NTK = NTK+1 NMSK = NMSK+1 NUK = NUK+1 NVK = NVK+1 NWK = NWK+1 DO 17 I = 1,IMAX+1 SIGMA1(I,J,K) = V1(I+2,K,1)*1.E10*QE/T1(I+2) SIGMA2(I,J,K) = V2(I+2,K,1)*1.E10*QE/T1(I+2) Z(I,J,K) = F(I+2,NZK) H(I,J,K) = gask*F(I+2,NTK)/(.5*(F(I+2,NMSK)+F(I+2,NMSK+1)) 1 *grav) U(I,J,K) = F(I+2,NUK) V(I,J,K) = F(I+2,NVK) W(I,J,K) = .5*(F(I+2,NWK)+F(I+2,NWK+1))*H(I,J,K) 17 CONTINUE NZK = NZK+1 DO 18 I = 1,IMAX+1 Z(I,J,KMAXP1) = F(I+2,NZK) 18 CONTINUE ! ! add additional fields calculated here ! gl - 20/07/2002 C Get ht integral of Ped, Hall conductance C write (6,"('LAMDAS: ht integration of sigp and sigh for C | j = ',i3)") j do i=1,imax+1 sigp(i,j) = 0. sigh(i,j) = 0. enddo nzk = nj+nz-1 do k=1,kmax nzk = nzk+1 do i=1,imax+1 C Times ht integ by 1.e-2 because appears off for siemens (1.e-2cm/m?) sigp(i,j) = sigp(i,j) + + (f(i,nzk+1)-f(i,nzk)) * sigma1(i,j,k) * 1.e-2 sigh(i,j) = sigh(i,j) + + (f(i,nzk+1)-f(i,nzk)) * sigma2(i,j,k) * 1.e-2 enddo enddo C calculate the Joule heating due to wind and ion drift C qwind = Joule Heting due to wind C qamie = Joule Heating due to ion drift C work = the mechanical work C wtot = total electromagnetic power do i=1,len1 C sndip = sin(dipmag) sndip(i) = sin(dipmag(i,j)) ! if (abs(sndip(i)) .lt. 1.e-3) ! | print *, 'sndip is zero at i= ',i,sndip(i) C term1=(bx**2 + bz**2) tm1(i) = sn2dec(i,j)+(1.-sn2dec(i,j))*sndip(i)**2 C term2=(by**2 + bz**2) tm2(i) = (1.-sn2dec(i,j))+sn2dec(i,j)*sndip(i)**2 C term3=bxby tm3(i) = sndec(i,j)*csdec(i,j) *(1.-sndip(i)**2) enddo C write (6,"('LAMDAS: ht integration of qwind and qamie')") qwind(:,j) = 0. qamie(:,j) = 0. wtot(:,j) = 0. work(:,j) = 0. fwindu(:,j) = 0. fwindv(:,j) = 0. famieu(:,j) = 0. famiev(:,j) = 0. nuik=ndj+nui-1 nvik=ndj+nvi-1 nuk=nj+nu-1 nvk=nj+nv-1 nzk = nj+nz-1 do k=1,kmax nuik = nuik + 1 nvik = nvik + 1 nuk = nuk + 1 nvk = nvk + 1 nzk = nzk+1 do i=1,imax+1 if (abs(sndip(i)) .gt. 0.2) then C write (6,"('LAMDAS: integration of qwind and qamie at C | i= ',i3)") i C qwind is the Joule heating associated with nuetral wind terms C qamie=sigp{(Bx^2+Bz^2)*Wx^2 + (By^2+Bz^2)*Wy^2-2.*Bx*By*Wx*Wy C -2.*B^2*Vx*Wx - 2.*B^2*Vy*Wy} qwind0 = sigma1(i,j,k)*bmod(i,j)**2* + (tm2(i)*f(i,nuk)**2+tm1(i)* + f(i,nvk)**2-2.*tm3(i)*f(i,nuk)*f(i,nvk) + -2.*(f(i,nuk)*f(i,nuik)+f(i,nvk)*f(i,nvik))) C qamie is the Joule heating withou nuetral wind C qamie=sigp{(Bx^2+Bz^2)*Vx^2 + (By^2+Bz^2)*Vy^2+2.*Bx*By*Vx*Vy}/sin(dip)^2 qamie0 = sigma1(i,j,k)*(tm1(i)*f(i,nuik)**2+tm2(i)* + f(i,nvik)**2+2.*tm3(i)*f(i,nuik)*f(i,nvik)) qamie0 = qamie0*bmod(i,j)**2 / sndip(i)**2 C wtot = total electric power C wtot = sigp (E + UxB).E = Qamie + sigp*{(UxB).E} wtot0 = - sigma1(i,j,k)*bmod(i,j)**2*(f(i,nuk)*f(i,nuik) + +f(i,nvk)*f(i,nvik)) + sigma2(i,j,k)*bmod(i,j)**2* + (tm3(i)*(f(i,nuk)* + f(i,nuik)-f(i,nvk)*f(i,nvik))-tm1(i)*f(i,nvk)*f(i,nuik) + +tm2(i)*f(i,nuk)*f(i,nvik))/sndip(i) wtot0 = qamie0 + wtot0 C Calculate the mechanical work done by wind work0 = wtot0 - qwind0 - qamie0 qwind_sec(i,j,k) = qwind0 * 1.e-9 qamie_sec(i,j,k) = qamie0 * 1.e-9 work_sec(i,j,k) = work0 * 1.e-9 wtot_sec(i,j,k) = wtot0 * 1.e-9 C To convert to mW/m^2 by mutipling by e-11 (since sigp (mho/m),height(cm) C B(in Gauss), and velocity (cm/s) qwind(i,j) = qwind(i,j) + + (f(i,nzk+1)-f(i,nzk)) * qwind0 * 1.e-11 qamie(i,j) = qamie(i,j) + + (f(i,nzk+1)-f(i,nzk)) * qamie0 * 1.e-11 wtot(i,j) = wtot(i,j) + + (f(i,nzk+1)-f(i,nzk)) * wtot0 * 1.e-11 work(i,j) = work(i,j) + + (f(i,nzk+1)-f(i,nzk)) * work0 * 1.e-11 C Calculate the horizontal currents fwindu0 = -tm1(i)*f(i,nvk)*sigma1(i,j,k) + + tm3(i)*f(i,nuk)*sigma1(i,j,k) + + sndip(i)*f(i,nuk)*sigma2(i,j,k) fwindv0 = tm2(i)*f(i,nuk)*sigma1(i,j,k) + + sndip(i)*f(i,nvk)*sigma2(i,j,k) + - tm3(i)*f(i,nvk)*sigma1(i,j,k) C Convert wind velocity from cm/s to m/s by multipling 1.e-2 and C Convert B from Gauss to T multipling 1.e-4 and 1.e-2 from cm to m fwindu0 = fwindu0*bmod(i,j)*1.e-8/sndip(i) fwindv0 = fwindv0*bmod(i,j)*1.e-8/sndip(i) famieu0 = -tm1(i)*sigma2(i,j,k)* + f(i,nuik) + sndip(i)*sigma1(i,j,k)*f(i,nvik) + - tm3(i)*sigma2(i,j,k)*f(i,nvik) famiev0 = -tm2(i)*sigma2(i,j,k)* + f(i,nvik) - sndip(i)*sigma1(i,j,k)*f(i,nuik) + - tm3(i)*sigma2(i,j,k)*f(i,nuik) C Convert wind velocity from cm/s to m/s by multipling 1.e-2 and C Convert B from Gauss to T multipling 1.e-4 and 1.e-2 from cm to m famieu0 = famieu0 * bmod(i,j)*1.e-8 / sndip(i)**2 famiev0 = famiev0 * bmod(i,j)*1.e-8 / sndip(i)**2 fwindu(i,j) = fwindu(i,j) + + (f(i,nzk+1)-f(i,nzk)) * fwindu0 fwindv(i,j) = fwindv(i,j) + + (f(i,nzk+1)-f(i,nzk)) * fwindv0 famieu(i,j) = famieu(i,j) + + (f(i,nzk+1)-f(i,nzk)) * famieu0 famiev(i,j) = famiev(i,j) + + (f(i,nzk+1)-f(i,nzk)) * famiev0 fwindu_sec(i,j,k) = fwindu0 fwindv_sec(i,j,k) = fwindv0 famieu_sec(i,j,k) = famieu0 famiev_sec(i,j,k) = famiev0 endif enddo enddo qwind_sec(:,j,zkmxp) = qwind(:,j) qamie_sec(:,j,zkmxp) = qamie(:,j) work_sec(:,j,zkmxp) = work(:,j) wtot_sec(:,j,zkmxp) = wtot(:,j) fwindu_sec(:,j,zkmxp) = fwindu(:,j) fwindv_sec(:,j,zkmxp) = fwindv(:,j) famieu_sec(:,j,zkmxp) = famieu(:,j) famiev_sec(:,j,zkmxp) = famiev(:,j) ! print *,'LAMDAS: j,famiev(:,j) = ',j, famiev(:,j) C Set the values at the poles sumfaus = 0. sumfavs = 0. sumfwus = 0. sumfwvs = 0. sumfaun= 0. sumfavn= 0. sumfwun= 0. sumfwvn= 0. do i=1,imax sumfwus = sumfwus + fwindu(i,1) sumfwun = sumfwun + fwindu(i,36) sumfwvs = sumfwvs + fwindv(i,1) sumfwvn = sumfwvn + fwindv(i,36) sumfaus = sumfaus + famieu(i,1) sumfaun = sumfaun + famieu(i,36) sumfavs = sumfavs + famiev(i,1) sumfavn = sumfavn + famiev(i,36) enddo do i=1,imax+1 fwindu(i,0) = sumfwus/72. fwindu(i,37) = sumfwun/72. fwindv(i,0) = sumfwvs/72. fwindv(i,37) = sumfwvn/72. famieu(i,0) = sumfaus/72. famieu(i,37) = sumfaun/72. famiev(i,0) = sumfavs/72. famiev(i,37) = sumfavn/72. enddo ! end of calculation of additional fields ! RETURN END C