!----------------------------------------------------------------------- subroutine cism_input_analysis ! ! 08/12 bae: Call cism_input_analysis here to: ! 1) calculate MIX energy flux geflux from gflx and geng ! 2) calculate the Hemispheric power for each hemisphere and define power as the average ! 3) convert geflux from geographic to magnetic coordinates if ifarad>0 ! 4) convert gpot from geographic to magnetic coordinates ! 5) call calccloc to find the Convection Location (theta0=convection radius) ! and arad (ifarad>0) ! 6) calculate ctpoten from the average of the NH and SH MIX potential drops ! 7) make any default corrections (offc=-0.4deg towards noon always; crad=Rc=12 deg ! when both hemispheres have ifbad>0) ! 8) colath called in heelis or in advance for NONE since .inp = 'HEELIS' or 'NONE' use cism_coupling_module,only: gpot,geng,gflx !CISM use input_module,only: byimf,ctpoten,power,kp ! byimf ! BY component of IMF (nT) (e.g., 0.) ! ctpoten ! cross-cap potential (kV) (e.g., 45.) ! power ! hemispheric power (GW) (e.g., 16.) ! kp ! interpolated kp from 3-h real Kp (e.g. 0 to 9) use magfield_module,only: sunlons ! sun's longitude in dipole coordinates in radians ! sunlons(geog nlat all same) from sub sunloc use cons_module,only: ! see cons.F | ylonm,ylatm, ! magnetic grid longitudes (nmlonp1) and latitudes (nmlat) in radians | dphi, ! delta latitude geographic (radians) | dtr, ! degrees to radians (pi/180) | rtd, ! radians to degrees (180/pi) | pi, ! 4.*atan(1.) | crit ! critical co-latitudes (dimensioned 2) in radians ! 08/12 bae: added for calccloc when ifbad>0 for both hemispheres (reset defaults) use aurora_module,only: theta0,offa,dskofa,phid,phin,rrad,offc, | dskofc use dynamo_module,only: phihm use params_module,only: | nlon, ! number of geographic longitudes (at 5 deg, nlon==72) | nlonp1, ! nlon+1 | nlonp4, ! nlon+4 | nlonp2, ! nlon+2 | nlev, ! number of pressure levels (e.g., 28) | nlevp1, ! nlev+1 | nlat, ! number of geographic latitudes (at 5 deg, nlat==36) | nlatp1, ! nlat+1 | nlatp2, ! nlat+2 | nmlon, ! number of geomagnetic grid longitudes | nmlonp1,! nmlon+1 | nmlat, ! number of geomagnetic grid latitudes | nmlatp1,! nmlat+1 | nmlath, ! (nmlat+1)/2 (index to magnetic equator) | nmlev, ! number of geomagnetic pressure levels (nmlev==nlevp1+3) | nmlevp1,! number of geomagnetic midpoint levels | nimlevp1,! number of geomagnetic interface levels | spval use magfield_module,only: ig,jg,wt ! 09/12 bae: Need geo2mag for arad calc from eflux use dynamo_module,only: geo2mag ! 11/08 EMERY added eflux=flux*fac_p2e*alfa (fac_p2e=2*1.602e-9, so 2*alfa in this factor) real, parameter:: fac_p2e = 3.204e-9 ! convert from particle to energy flux real, dimension(nlonp1,0:nlatp1) :: geflux ! energy flux from CMIT in geog coords, but with poles at 0 and nlatp1 (REQUIRED for geo2mag!!) real, dimension(nmlonp1,nmlat) :: meflux ! energy flux from CMIT in mag coords real :: Re,areatot,dglat,dglon,ghpi_sh_mix,ghpi_nh_mix real, dimension(nlat) :: da,glatdeg integer :: nl4,nml4,i,jj,j,ifcalccloc,ifarad real :: pcp_sh_mix,pcp_nh_mix ! cross-tail potential drop in kV for SH and NH from MIX real :: pcp_sh_min,pcp_nh_min,pcp_sh_max,pcp_nh_max ! min/max kV potentials real :: ghpi_sh_max,ghpi_nh_max ! 08/12 bae: added for calccloc real :: model_ctpoten(2),mltd,mltn,byloc integer :: ifbad(2) real, dimension(nmlonp1) :: mltSh,mltNh real, dimension(nmlat/4-1) :: magSdeg,magNdeg ! 09/12 changed dimension to nmlonp1 to correspond to dimension in calccloc! real, dimension(nmlonp1,nmlat/4-1) :: potS,potN,efxS,efxN ! 09/12 bae: real :: maxefx(2) integer :: i180 ! TEST turning off all calccloc: ifcalccloc=1,0 if do,not call calccloc ifcalccloc = 1 ! ifarad = 1,0 if do,not calculate the auroral radius (arad) from efxS,N ! ifarad = 2 (do 2 things:) calc Ra and use Rc=Ra-4 deg instead of ! the default 12 deg for Rc when both ifbad>0 ifarad = 1 ! Limit size of byimf in phin and phid calculations (as in aurora.F) ! NOTE: This byloc is assymetric in hemisphere, which WAS correct ! according to Emery et al, 2012 (NCAR TN#491, ! http://nldr.library.ucar.edu/repository/collections/TECH-NOTE-000-000-000-856.) byloc = byimf if (byloc .gt. 7.) byloc = 7. if (byloc .lt. -11.) byloc = -11. ! Print-out testing: ! write (6,"(1x,'cism_input_analysis byimf,loc,power,ctpoten,kp=', ! | 2f7.2,3e12.4)") byimf,byloc,power,ctpoten,kp ! Should calculate hemispheric power from MIX energy flux in GW for both hem and store in power ! Calculate min/max CP and hpi_s,nh_amie in GW ! Can calculate HP in geog coord instead of in magnetic coords ! Re = radius of the Earth in m at the level of the ionosphere (use 6370e+3 in MIX) Re = 6371.e+3 areatot = 0. dglat = 180./nlat dglon = 360./nlon do j=1,nlat glatdeg(j) = -90.0 + dglat/2. + (j-1)*dglat ! NOTE: The HP from this calc at the cos(lat) of the data is ~12% less than the ! MIX calcFaceArea python calculations. Can get good agreement in the MIX grid ! if change da(j) to be for the j-1 latitude equatorwards of the j data. (Won't do) da(j) = Re*cos(glatdeg(j)*dtr)*dglon*dtr*(Re*dglat*dtr) areatot = areatot + da(j)*nlon enddo ! geng from MIX has been divided by 2 to go from mean (av) energy to Characteristic (or Maxwellian) or alpha energy in keV ! 11/08 EMERY added eflux=flux*fac_p2e*alfa (fac_p2e=2*1.602e-9, so 2*alfa in this factor) ! write (6,"(1x,'areatot,4pi fac_p2e gflx(1,1),geng(1,1)=',5e12.4)") ! | areatot,areatot/(Re*Re),fac_p2e,gflx(1,1),geng(1,1) maxefx(1) = 0. maxefx(2) = 0. do j=1,nlat ! Found geflux(nlonp1,j) not same as geflux(1,j), so go only to nlon and set for nlonp1! do i=1,nlon ! NOTE: nlonp4 is 2 points before i=1 and 2 points after i=nlon geflux(i,j) = gflx(i+2,j)*fac_p2e*geng(i+2,j) ! Find peak fluxes in SH (1) and NH (2) if (j .lt. nlat/2) maxefx(1) = max(geflux(i,j),maxefx(1)) if (j .ge. nlat/2) maxefx(2) = max(geflux(i,j),maxefx(2)) enddo geflux(nlonp1,j) = geflux(1,j) ! periodic point enddo ! Now calculate geflux(i,0) and geflux(i,nlatp1) at across the SH and NH poles for geo2mag ! (NOTE: 0 never used since ig goes from 1-36 and geo2mag uses 1-37) do i=1,nlon i180 = i + nlon/2 if (i180 .gt. nlonp1) i180 = i180 - nlon geflux(i,0) = geflux(i180,1) geflux(i,nlatp1) = geflux(i180,nlat) enddo ! write (6,"(1x,'maxefx gS,N mW/m2 =',2e12.4)") maxefx ! Calculate Hemispheric power (HP) for each hemisphere from the pole to ~45 deg (nlat/4) ! SH nl4 = nlat/4+1 ghpi_sh_mix = 0. ghpi_sh_max = 0. do i=1,nlon do j=1,nl4 ghpi_sh_mix = ghpi_sh_mix + da(j)*geflux(i,j) ghpi_sh_max = max(geflux(i,j),ghpi_sh_max) enddo enddo ghpi_nh_mix = 0. ghpi_nh_max = 0. do i=1,nlon do jj=1,nl4 j = nlat - jj + 1 ghpi_nh_mix = ghpi_nh_mix + da(j)*geflux(i,j) ghpi_nh_max = max(geflux(i,j),ghpi_nh_max) enddo enddo ! Convert from mW/m2 * m2(da) * 1.e-12GW/mW to GW ghpi_sh_mix = ghpi_sh_mix*1.e-12 ghpi_nh_mix = ghpi_nh_mix*1.e-12 power = 0.5*(ghpi_sh_mix + ghpi_nh_mix) ! Print-out testing: write (6,"(1x,'MIX HP power (S,N,b), maxS,N mW/m2=',5e12.4)") | ghpi_sh_mix,ghpi_nh_mix,power,ghpi_sh_max,ghpi_nh_max ! This section on finding min/max from gpot can be deleted unless ! it is desired to use ctpoten from gpot instead of from phihm ! Find gpot min/max potentials pcp_sh_mix = 0. pcp_nh_mix = 0. pcp_sh_min = 0. pcp_nh_min = 0. pcp_sh_max = 0. pcp_nh_max = 0. do i=1,nlon ! NOTE: nlonp4 is 2 points before i=1 and 2 points after i=nlon do j=1,nl4 pcp_sh_min = min(gpot(i+2,j),pcp_sh_min) pcp_sh_max = max(gpot(i+2,j),pcp_sh_max) enddo do jj=1,nl4 j = nlat - jj + 1 pcp_nh_min = min(gpot(i+2,j),pcp_nh_min) pcp_nh_max = max(gpot(i+2,j),pcp_nh_max) enddo enddo ! Convert from V to kV pcp_sh_min = pcp_sh_min*1.e-3 pcp_nh_min = pcp_nh_min*1.e-3 pcp_sh_max = pcp_sh_max*1.e-3 pcp_nh_max = pcp_nh_max*1.e-3 pcp_sh_mix = pcp_sh_max - pcp_sh_min pcp_nh_mix = pcp_nh_max - pcp_nh_min ! Print-out testing: ! write (6,"(1x,'MIX gpot model_ctpoten (S,N,av) =',3e12.4)") ! | pcp_sh_mix,pcp_nh_mix,0.5*(pcp_sh_mix+pcp_nh_mix) ! Found 5 deg grid gives ~+1kV CP at 06348 2201 UT: 163S+199N~181kV vs 162S+198N~180kV ! vs MIX values of 165.0S+200.4N~182.7kV or ~3/180 or 2% less ! For lower kV, 5-10kV less ~50-75kV, or ~10-15% less CP ! Set MLT (mlat=ylatm*rtd) ! Find MLT from sunlons (sunlons(1-nlat) are all the same value so use the first value) ! 09/12 changed nmlon to nmlonp1 do i=1,nmlonp1 mltSh(i) = (ylonm(i)-sunlons(1)) * rtd / 15. + 12. if (mltSh(i) .gt. 24.) mltSh(i) = mltSh(i) - 24. if (mltSh(i) .lt. 0.) mltSh(i) = mltSh(i) + 24. mltNh(i) = mltSh(i) enddo ! write (6,"(1x,'nmlonp1 sunlons MltS,N,ylonm(1,p1) =',i3, ! | 7e12.4)") nmlonp1,sunlons(1)*rtd,mltSh(1),mltNh(1), ! | ylonm(1)*rtd,mltSh(nmlonp1),mltNh(nmlonp1),ylonm(nmlonp1)*rtd if (ifarad .gt. 0) then ! Convert MIX geflux from geographic to geomag coords and set efxS and N do j=1,nmlat call geo2mag(meflux(1,j),geflux(1:nlonp1,:),ig,jg,wt,nlonp1, | nmlonp1,nmlon,nmlat,j) meflux(nmlonp1,j) = meflux(1,j) ! periodic point enddo endif ! ifarad>0 ! write (6,"(1x,'ifcalccloc =',i1)") ifcalccloc if (ifcalccloc .eq. 1) then ! ! Convert LFM geographic potential gpotm(nlonp4,0:nlatp1) to geomagnetic coordinates phihm(nmlonp1,nmlat) call cism_pot2mag ! write (6,"(1x,'cism_pot2mag phihm(1,1)=',e12.4)") phihm(1,1) ! get ~1/4 of globe (more points in equatorial region than in polar regions so subtract 1) ! Get points from poles to about +/-35 mlat ! nml4=23 from -90 about every 2 deg to -70, every 3 to -30: ! 23 magS = -0.9000E+02 -0.8812E+02 -0.8624E+02 -0.8433E+02 -0.8240E+02 ! -0.8043E+02 -0.7841E+02 -0.7633E+02 -0.7419E+02 -0.7197E+02 ! -0.6967E+02 -0.6728E+02 -0.6480E+02 -0.6222E+02 -0.5955E+02 ! -0.5678E+02 -0.5393E+02 -0.5100E+02 -0.4801E+02 -0.4498E+02 ! -0.4192E+02 -0.3885E+02 -0.3582E+02 nml4 = nmlat/4-1 ! Save about half nmlat per hemisphere (nml4) do jj=1,nml4 magSdeg(jj) = ylatm(jj)*rtd ! 09/12 changed nmlon here and below to nmlonp1 (nmlonp1=nmlon+1) do i=1,nmlonp1 if (ifarad>0) efxS(i,jj) = meflux(i,jj) potS(i,jj) = phihm(i,jj) enddo enddo do j=1,nml4 jj = nmlat - j + 1 magNdeg(j) = ylatm(jj)*rtd do i=1,nmlonp1 if (ifarad>0) efxN(i,j) = meflux(i,jj) potN(i,j) = phihm(i,jj) enddo enddo ! write (6,"(1x,'magS =',10e12.4)") magSdeg ! write (6,"(1x,'magN =',10e12.4)") magNdeg ! write (6,"(1x,'nml4 magS,N potS,N =',i3,4e12.4)") ! | nml4,magSdeg(1),magNdeg(1),potS(1,1),potN(1,1) ! 08/12 bae: Call calccloc here to get revised variables for aurora_cons ! - theta0 (and phid for cusp location) are used later in crit, drizzle, and cusp ! drizzle and cusp are found in aurora.F which is called in dynamics ! calccloc sets phid and phin to be no more than +/2h from defaults with byloc (usually 0) ! 09/12 bae: efxS and efxN are used to find the auroral radius instead of arad=crad+2 call calccloc (magSdeg,magNdeg,nml4,mltSh,mltNh, | nmlonp1,potS,potN,efxS,efxN,ifarad,model_ctpoten,ifbad) ! 09/12 changed nmlon to nmlonp1 in call to calccloc ctpoten = 0.5 * (model_ctpoten(1)+model_ctpoten(2)) ! Print-out testing: write (6,"(1x,'MIX ifbad model_ctpoten (S,N,av) =',2i2,3e12.4)") | ifbad,model_ctpoten,ctpoten ! Choose if pcp_s,nh_mix from gpot or phihm - but is not used anyway, just calculated! ! pcp_sh_mix = model_ctpoten(1) ! pcp_nh_mix = model_ctpoten(2) ! Print-out testing: write (6,"(1x,'CalcCLoc S,N: offc,dskofc,theta0,phid,n=', | 10e12.4)") offc*rtd,dskofc*rtd,theta0*rtd, | phid*rtd/15.-12.,phin*rtd/15.-12. ! Set center of convection to always be in the center (and use subsequent defaults for the auroral center) ! Now set rrad and theta0 and others to defaults ! 9/12 bae: cannot have offc=dskofc=0; offc=-0.4deg towards noon from MIX average offc(1) = -0.4/rtd offc(2) = -0.4/rtd ! TEMP TEST center at pole - no dynamo (pfrac=0 everywhere) problem fixed in colath.F ! offc(1) = -0.0/rtd ! offc(2) = -0.0/rtd offa(1) = offc(1) offa(2) = offc(2) dskofc(1) = 0./rtd dskofc(2) = 0./rtd dskofa(1) = dskofc(1)-2.5/rtd dskofa(2) = dskofc(2)-2.5/rtd ! Set phid (and phin) to the defaults ! phid is used in the cusp location which will be removed or replaced by the Univ of NH aurora including cusp ! Use parameterization defaults for phid (phid(MLT)=9.39 +/- 0.21By - 12) ! and phin (phin(MLT)=23.50 +/- 0.15By - 12) ! NOTE: byloc and byimf are usually zero for CMIT runs so phid=9.39 and phin=23.50 MLT mltd = 9.39 + 0.21*byloc phid(1) = (mltd-12.) * 15.*dtr mltd = 9.39 - 0.21*byloc phid(2) = (mltd-12.) * 15.*dtr mltn = 23.50 + 0.15*byloc phin(1) = (mltn-12.) * 15.*dtr mltn = 23.50 - 0.15*byloc phin(2) = (mltn-12.) * 15.*dtr ! 08/12 bae: If both hems cannot be calculated (ifbad>0 in both), set theta0 defaults ! 09/12 bae: The auroral radius can be used to calculate theta0 for ifbad>0 if (ifbad(1)>0 .and. ifbad(2)>0) then theta0(1) = 12./rtd theta0(2) = 12./rtd if (ifarad .eq. 0) then rrad(1) = theta0(1) + 2./rtd rrad(2) = theta0(2) + 2./rtd endif endif ! Print-out testing: ! write (6,"(1x,'After Rev S,N: offc,dskofc,theta0,rrad,phid,n=', ! | 12e12.4)") offc*rtd,dskofc*rtd,theta0*rtd,rrad*rtd, ! | phid*rtd/15.-12.,phin*rtd/15.-12. endif ! for ifcalccloc>0 (test turning it on or off here and in aurora.F) end subroutine cism_input_analysis !-----------------------------------------------------------------------