C Calculate conductance from EUV from solcond (PEDSolar), x1fitall and QJ from fini in AMIE C plus Robinson formula using mean electron energy and energy flux for aurora (PEDAurora) C Calculate the estimated Joule heating (QJ) from DMSP dawn-dusk passes using IDM Vy and Ped(A+S) C assuming that RPA Vx ram direction is not as important, and ignoring the neutral wind Vn C 6/11 Barbara Emery: This code has not been tested and won't work. It needs C to be revised to work using DMSP satellite data. common/solap/sdate,clatp,polon,sbsllat,sbsllon,wkap(lwk) COMMON/FXFY/TWOPI,ROT C LONMX = number of longitudes (+1?) in your mlon array (not same for PII=4.*ATAN(1.) RAD = 180./PII TWOPI=2.*PII IYD = (IYEAR-1900)*1000 + IDAYN UT=FLOAT(IHR)+(FLOAT(IMIN)/60.) TYME(NTIME) = FLOAT(IDA-IDA1)*24. + UT C ROT IS MAGNETIC LONGITUDE OF LOCAL MIDNIGHT IN RADIANS. ROT = 15.*IHR + IMIN/4.E0 - 70. write (6,"(1x,'rot =', f7.2)") -rot ROT = -ROT/RAD C Find subsolar point from subsol call subsol (iyear,idayn,ihr,imin,0.,sbsllat,sbsllon) write (6,"(1x,'subsolar point from subsol(lat,lon) =',2f8.2)") | sbsllat,sbsllon C Find alon at 0 MLT from mlt2alon in magloctm.f and revise ROT call mlt2alon (0.,sbsllat,sbsllon,clatp,polon,alonx) write (6,"(1x,'mag lon of 0 MLT from mlt2alon =',f8.2)") alonx write (6,"(1x,'rot =', f7.2)") alonx rot = alonx/rad C c-------------------------------------------------------------------- C CONDUCT (and AURCOND) assumes a colatitude between 0 and 90 degrees C (N.H.), so IDAYMN and OFFMAG are different for different IHSOLN. C Save both N.H. and S.H. CCHI and SCHI for use with total conductance C (which is radar and magnetometer conductance). C IDAYMN = day number of minimum solar radiance IDAYMNN = -11 SDN = -SIN(23.5/RAD)*COS((IDAYN-IDAYMNN)*TWOPI/365.24) CDN = SQRT(1. - SDN*SDN) IDAYMNS = 172 SDS = -SIN(23.5/RAD)*COS((IDAYN-IDAYMNS)*TWOPI/365.24) CDS = SQRT(1. - SDS*SDS) C OFFMAG = offset of the magnetic pole from the geographic pole OFFMAGN = 8.7 SEN = SIN(OFFMAGN/RAD) CEN = SQRT(1. - SEN*SEN) OFFMAGS = 16.0 SES = SIN(OFFMAGS/RAD) CES = SQRT(1. - SES*SES) C CCHI,SCHI = cos and sin of the mag colat of the subsolar point CCHIN = CEN*SDN + SEN*CDN*COS((-83. + 15.*(UT-12.))/RAD) SCHIN = SQRT(1. - CCHIN*CCHIN) CCHIS = CES*SDS + SES*CDS*COS((127. + 15.*(UT-12.))/RAD) SCHIS = SQRT(1. - CCHIS*CCHIS) IF (IHSOLN .EQ. 1) THEN CCHI = CCHIN SCHI = SCHIN ELSE CCHI = CCHIS SCHI = SCHIS ENDIF WRITE (6,"(1X/1X,I3,1X,I5,4I3,' IDAYN=',I5/ | 1X,'MAG LON OF 24 MLT, MAG COLAT OF SUBSOLAR PT =', 2F8.2)") | NTIME,IYEAR,IMON,IDA,IHR,IMIN,IDAYN, | ROT*RAD,ATAN2(SCHI,CCHI)*RAD c-------------------------------------------------------------------- C Calculate the Pedersen and Hall conductance using Robinson's formulae C Robinson et al, JGR, 2565-2569, 1987. C See also Rich et al, Ann. Geophys., 527-534, 1987. Please note that C there is an error in eq 14. Do not divide by 2! -- This reference C is how the DMSP people do it. They only look at energies above 450 eV C usually. DO 1540 ITH=0,ITHMX DO 1540 LON=0,LONMX C Z= COLATITUDE in radians C Y=LONGITUDE FROM MIDNIGHT in radians Z = ... Y = ... C EFLX2D is the electron energy flux in mW/m2 (ergs/cm2) C EKEV2D is the mean electron energy in keV C SIGP(edersen)A(aurora) and SIGPS(olar) are in ??, as is the HALL conductance (SIGH) CALL SOLCOND(Z,Y,CCHI,SCHI,SPED,SHAL) SIGPA(LON,ITH) = SQRT(EFLX2D(LON,ITH)) * 40. * EKEV2D(LON,ITH) | / (16. + EKEV2D(LON,ITH)**2) SIGHA(LON,ITH) = 0.45 * (EKEV2D(LON,ITH)**0.85) * SIGPA(LON,ITH) SIGP(LON,ITH) = SQRT(SIGPA(LON,ITH)**2 + SIGPS(LON,ITH)**2) SIGH(LON,ITH) = SQRT(SIGHA(LON,ITH)**2 + SIGHS(LON,ITH)**2) 1540 CONTINUE C Calculate the estimated Joule heating in mW/m2(?) C QJ(?) ~ E(mV/m?)*Ped(A+S)(mho?) (need to work out units etc) C Vi = ExB/B^2 so estimate E=Vy(m/s?)*B(nT?) C CALCULATE JOULE HEATING (STORED IN SMQJ in units??) DO 5130 LON=1,LONMX 5130 CURT(LON) = 0. DO 5137 ITH=0,ITHTRNS DO 5135 LON=1,LONMX C EE and EN are the electric field in the east and north directions in units?? SMQJ(LON,ITH) = SIGP(LON,ITH)*(EE(LON,ITH)*EE(LON,ITH) + 1 EN(LON,ITH)*EN(LON,ITH)) IF (ITH .GT. ITHTRNS) GO TO 5135 C Only calculate hemispheric average to ITHTRNS CURT(LON) = CURT(LON) + SMQJ(LON,ITH)*ST(ITH) 5135 CONTINUE SMQJ(0,ITH) = SMQJ(LONMX,ITH) 5137 CONTINUE TOTJUL = 0. DO 5138 LON=1,LONMX 5138 TOTJUL = TOTJUL + CURT(LON) SMPLJL(NTIME) = TOTJUL*RI*RI*DLON*DTH TJH(NTIME) = 0. stop end SUBROUTINE SOLCOND(Z,Y,CCHI,SCHI,SGP,SGH) SAVE C include "Headers/fxfy.h" COMMON/FXFY/TWOPI,ROT C Z= COLATITUDE C Y=LONGITUDE FROM MIDNIGHT PII = TWOPI/2.E0 SZ=SIN(Z) CZ=COS(Z) C COSK=CZ*CCHI + SZ*COS(Y-PII-ROT)*SCHI COSK=CZ*CCHI + SZ*COS(Y-PII )*SCHI FCOSK=0.03+EXP(1.803*TANH(3.833*COSK)+0.5*COSK-2.332) FCOSK=FCOSK+0.03 BBP=SQRT(1.0-0.99524*SZ*SZ)*(1.0+0.3*CZ*CZ) BBH=SQRT(1.0-0.01504*(1.0-CZ)-0.97986*SZ*SZ)*(1.0+0.5*CZ*CZ) SGP=12.5*FCOSK/BBP SGH=17.0*FCOSK/BBH RETURN END