#include "defs.h" ! module low_proton_module ! ! Module used to read data from the auroral proton data (mean energy, ! and energy flux). ! use params_module,only: nlon,nlat,nlonp1,nlonp4 implicit none ! #ifdef SUN #include "netcdf3.inc" #else #include "netcdf.inc" #endif ! ! Define parameters for auroral proton data file: integer,parameter :: | mxgdays = 1, ! maximum number of days of proton data | mxtimes = 721, ! maximum number of times of proton data per day | ilat = 11, ! number of latitudes of proton data (50/5+1) | ilon = 73 ! number of longitudes of proton data(360/5+1) integer :: lonp1,latp1, | year(mxtimes),month(mxtimes),day(mxtimes),jday(mxtimes) ! Define auroral proton output fields real :: kev_lp(nlonp4,nlat),efx_lp(nlonp4,nlat) ! ! Define fields for auroral proton input data file: ! mean energy in KeV ! energy flux in erg.cm-2.s-1 (or mW/m^2) real, allocatable,dimension(:) :: proton_ut real,allocatable,dimension(:,:,:) :: ! (lonp1,latp1,ntimes) | fuv_kev_lp, fuv_efx_lp ! contains !----------------------------------------------------------------------- subroutine rdproton ! ! Read AMIE data for the northern hemisphere from amienh ! use input_module,only: tempdir,aurora_proton ! ! Local: character(len=80) :: dskfile integer :: istat,ntimes,ncid,ndims,nvars,ngatts,idunlim,ier integer :: id_lon,id_lat,id_time,lon,lat,mlt, | idv_year,idv_mon,idv_day,idv_jday, | idv_iter,idv_ut,idv_lon,idv_lat,idv_kev,idv_efx ! ! Acquire mss file: ! call mkdiskflnm(aurora_proton,dskfile) ! call getms(aurora_proton,dskfile,tempdir,' ') dskfile = ' ' call getfile(aurora_proton,dskfile) ! write(6,"(/,72('-'))") write(6,"('RDPROTON: read proton data')") ! ! Open netcdf file: istat = nf_open(dskfile,NF_NOWRITE,ncid) if (istat /= NF_NOERR) then write(6,"(/,'>>> rdproton: error opening netcdf proton data ', | 'diurnal file ',a)") trim(dskfile) stop 'rdproton' else endif ! ! Get proton grid dimension: istat = nf_inq_dimid(ncid,'dim1_LON',id_lon) istat = nf_inq_dimlen(ncid,id_lon,lonp1) if (istat /= NF_NOERR) call rpt_ncerr(istat, | 'rdproton: Error getting proton longitude dimension') ! lonp1 = 72 istat = nf_inq_dimid(ncid,'dim1_LAT',id_lat) istat = nf_inq_dimlen(ncid,id_lat,latp1) if (istat /= NF_NOERR) call rpt_ncerr(istat, | 'rdamie_nh: Error getting proton latitude dimension') ! latp1 = 10 write(6,"('lonp1=',i3,' latp1=',i3)") lonp1,latp1 ! ! Get time dimension: istat = nf_inq_unlimdim(ncid,id_time) istat = nf_inq_dimlen(ncid,id_time,ntimes) if (istat /= NF_NOERR) call rpt_ncerr(istat, | 'rdproton: Error getting time dimension') ! write(6,"('rdproton: ntimes=',i3)") ntimes ! ! Search for requested proton output fields istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim) ! ! Get 1-D proton fields (ntimes) istat = nf_inq_varid(ncid,'YEAR',idv_year) istat = nf_get_var_int(ncid,idv_year,year) istat = nf_inq_varid(ncid,'MONTH',idv_mon) istat = nf_get_var_int(ncid,idv_mon,month) istat = nf_inq_varid(ncid,'DAY',idv_day) istat = nf_get_var_int(ncid,idv_day,day) istat = nf_inq_varid(ncid,'JDAY',idv_jday) istat = nf_get_var_int(ncid,idv_jday,jday) ! ! Allocate 1-d fields: ! if (.not. allocated(proton_ut)) | allocate(proton_ut(ntimes),stat=ier) if (ier /= 0) write(6,"('>>> rdproton: error allocating', | ' proton_ut: ntimes=',i3)")ntimes ! ! Get ut istat = nf_inq_varid(ncid,'UT',idv_ut) if (istat /= NF_NOERR) call rpt_ncerr(istat, | 'rdamie_nh: Error getting auroral proton UT id') istat = nf_get_var_double(ncid,idv_ut,proton_ut) if (istat /= NF_NOERR) call rpt_ncerr(istat, | 'rdamie_nh: Error getting proton variable ut') ! ! Allocate 3-d fields: ! real,allocatable,dimension(:,:,:) :: ! (lonp1,latp1,ntimes) ! | kev_lp, efx_lp ! if (.not. allocated(fuv_kev_lp)) | allocate(fuv_kev_lp(lonp1,latp1,ntimes),stat=ier) if (ier /= 0) write(6,"('>>> rdproton: error allocating', | ' kev_lp: lonp1=',i3,' latp1=',i3,' ntimes=',i3)") | lonp1,latp1,ntimes if (.not. allocated(fuv_efx_lp)) | allocate(fuv_efx_lp(lonp1,latp1,ntimes),stat=ier) if (ier /= 0) write(6,"('>>> rdamie_nh: error allocating', | ' efx_lp: lonp1=',i3,' latp1=',i3,' ntimes=',i3)") | lonp1,latp1,ntimes ! ! Get 3-D AMIE fields (lon,lat,ntimes) ! ! Proton char. energy istat = nf_inq_varid(ncid,'P_CHAR',idv_kev) if (istat /= NF_NOERR) call rpt_ncerr(istat, | 'rdproton: Error getting proton char. energy id') istat = nf_get_var_double(ncid,idv_kev,fuv_kev_lp) if (istat /= NF_NOERR) call rpt_ncerr(istat, | 'rdproton: Error getting proton char. energy') ! ! Proton energy flux istat = nf_inq_varid(ncid,'P_FLUX',idv_efx) if (istat /= NF_NOERR) call rpt_ncerr(istat, | 'rdproton: Error getting proton energy flux id') istat = nf_get_var_double(ncid,idv_efx,fuv_efx_lp) if (istat /= NF_NOERR) call rpt_ncerr(istat, | 'rdproton: Error getting proton energy flux') ! ! Close the file: istat = nf_close(ncid) write(6,"('Completed read from auroral proton data file ',a)") | trim(dskfile) write(6,"(72('-'),/)") end subroutine rdproton !----------------------------------------------------------------------- subroutine getproton(iyear,iday,iutsec,iprint) ! ! Read AMIE outputs from amie_ncfile file, returning electric potential, ! auroral mean energy and energy flux at current date and time, ! and the data is linearly interpolated to the model time ! gl - 12/07/2002 ! use hist_module,only: modeltime use input_module,only: aurora_proton use cons_module,only: pi,ylonm,ylatm use magfield_module,only: im,jm,dim,djm use dynamo_module,only: mag2geo ! ! Args: integer,intent(in) :: iyear,iday,iutsec,iprint ! ! Local: integer :: ier,lw,liw,isign,intpol(2) integer,allocatable :: iw(:) real,allocatable :: w(:) integer :: i,j,ndmon_nl(13),ndmon_lp(13),nonleap integer :: nn, iset, iset1, m, mp1, n, proton_debug,protonday real :: model_ut, denoma, f1, f2 real :: del,xmlt,dmlat,dlatm,dlonm,dmltm,rot,dtr ! ! Data: GSWM data given at the middle of each month ! -> day of year for this data assuming non-leap year ! J F M A M J J A S O N D J data ndmon_nl ! non_leap | /31,59,90,120,151,181,212,243,273,304,334,365,396/ data ndmon_lp ! leap year | /31,60,91,121,152,182,213,244,274,305,335,366,397/ ! ! External: real,external :: finterp ! if (iprint > 0) then write(6,"(/,72('-'))") write(6,"('GETPROTON:')") write(6,"('Initial requested iyear=',i4,' iday=',i3,' iutsec=', | i10)") iyear,iday,iutsec endif ! ! Check for leap year nonleap = int(mod(real(iyear),4.)) ! nonleap year /= 0 ! if (iyear.eq.2000) nonleap = 1 ! ! Check times: ! if (year(1) /= iyear) then write(6,"('>>> getproton: wrong year,iyear=',i4,2x,i4)") | year(1),iyear stop 'getproton' endif protonday = jday(1) if (protonday /= iday) then write(6,"('>>> getproton: wrong protonday,modelday=',2i4)") | protonday,iday stop 'getproton' endif ! ! if (lonp1 /= nlonp1) then ! write(6,"('>>> getproton: wrong longitude dimension =',2i4)") ! | lonp1,nlonp1 ! stop 'getproton' ! endif ! if (modeltime(3).lt.0.or.modeltime(3).gt.59) then write(6,"('>>> getproton: bad imin=',3i4)") modeltime(3) stop 'getproton' endif if (modeltime(4).lt.0.or.modeltime(4).gt.59) then write(6,"('>>> getproton: bad isec=',3i4)") modeltime(4) stop 'getproton' endif model_ut = (modeltime(1)-protonday)*24. + float(iutsec)/3600. ! ! Interpolate proton data to modeltime iutsec ! amie_ibkg = 0 use real UT AMIE data ! = 1 use the first AMIE volumne as the background ! = 2 use the 24-hr average AMIE volumne as the background kev_lp(:,:) = 0. efx_lp(:,:) = 0. ! n = len_trim(aurora_proton) if (n > 0) then ! write(6,"('getproton: Interpolate NH proton Data')") nn = size(proton_ut) iset = nn iset1 = nn do i=1,nn if (proton_ut(i) < model_ut) iset = i enddo iset1 = iset + 1 if (iset == nn) iset1 = iset denoma = proton_ut(iset1) - proton_ut(iset) if (denoma == 0.) then f1 = 1. f2 = 0. else f1 = (proton_ut(iset1) - model_ut)/denoma f2 = (model_ut - proton_ut(iset))/denoma endif ! ! assign the proton data to TGCM grids do i=1,lonp1 do j=1,latp1 kev_lp(i+2,nlat-j+1) = (f1*fuv_kev_lp(i,j,iset1) + | f2*fuv_kev_lp(i,j,iset)) efx_lp(i+2,nlat-j+1) = (f1*fuv_efx_lp(i,j,iset1) + | f2*fuv_efx_lp(i,j,iset)) if (efx_lp(i+2,nlat-j+1) .gt. 0. .and. | kev_lp(i+2,nlat-j+1) .lt. 1.) kev_lp(i+2,nlat-j+1)=1.0 ! kev_lp(i,nlat+j-10) = (f1*fuv_kev_lp(i,j,iset1) + ! | f2*fuv_kev_lp(i,j,iset)) ! efx_lp(i,nlat+j-10) = (f1*fuv_efx_lp(i,j,iset1) + ! | f2*fuv_efx_lp(i,j,iset)) enddo enddo kev_lp(nlonp1,:) = kev_lp(3,:) efx_lp(nlonp1,:) = efx_lp(3,:) ! do j=31,nlat ! write(6,"('j = ',i3,'efx_lp(:,j) = ',/,(6f12.2))") ! | j,efx_lp(:,j) ! enddo ! do i=1,2 kev_lp(i,:) = kev_lp(nlon+i,:) kev_lp(nlonp4+i-2,:) = kev_lp(i+2,:) efx_lp(i,:) = efx_lp(nlon+i,:) efx_lp(nlonp4+i-2,:) = efx_lp(i+2,:) enddo ! call addfsech_ij('PRO_KEV','Proton Mean Energy','KeV ', | kev_lp,1,nlonp4,1,nlat) call addfsech_ij('PRO_EFX','Proton Energy Flux','mW/m^2', | efx_lp,1,nlonp4,1,nlat) endif ! if (iprint > 0) then write(6,"('getproton: proton data interpolated ', | 'to date and time')") write(6,"(72('-'),/)") endif end subroutine getproton !------------------------------------------------------------------- subroutine rpt_ncerr(istat,msg) ! ! Handle a netcdf lib error: ! integer,intent(in) :: istat character(len=*),intent(in) :: msg ! write(6,"(/72('-'))") write(6,"('>>> Error from netcdf library:')") write(6,"(a)") trim(msg) write(6,"('istat=',i5)") istat write(6,"(a)") nf_strerror(istat) write(6,"(72('-')/)") return end subroutine rpt_ncerr end module low_proton_module !----------------------------------------------------------------- ! ! 5/29/98 btf: Original code from Marina Galand, here modified for ! use in tiegcm (tgcm13) ! 9/9/99 btf: added implicit none for tgcm13mt. ! 10/13/03 G. Lu: modified for tiegcm1 ! subroutine low_proton(z,tn,o1,o2,barm,qia,lev0,lev1,lon0,lon1,lat) ! use cons_module,only: | p0, ! standard pressure (5.0e-4) | expz, ! exp(-zp) where zp is log pressure (ln(p0/p)) | grav, ! gravity | gask, ! gas constant (8.314e7) | boltz, ! boltzman's constant (1.38e-16) | rmass_o2, rmass_o1, rmass_n2, ! molelcular weights of o2,o,n2 | rmassinv_o2,rmassinv_o1,rmassinv_n2, ! 1./rmass | avo ! avogadro number (6.023e23) ! use low_proton_module,only: kev_lp,efx_lp ! implicit none ! ! Determination of electron and ion production rates ! - noted QTI and QIA respectively - ! induced by a proton beam of an average energy of ! ALP keV and an incident energy flux of EFLUX erg.cm-2.s-1. ! Restriction : 2 < ALP < 40 keV ! Note : ALP = 2*E0 (for a Maxwellian) ! ! Args: integer :: lev0,lev1,lon0,lon1,lat real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) :: | z,tn,o2,o1,barm ! ! Local: ! qia is output ion production for N2+, O2+, O+, NO+, N+ real,dimension(lev0:lev1,lon0:lon1,5),intent(out) :: qia ! ! Locals: integer :: i,k,n,m,ntk,npo2k,npo1k,nzk,ii real :: ALPV,EL,KL,RANGE,F_E,R_norm(lon0:lon1), + LAMBDA,ETA,W(6),shape(3),Kd(4),E,QIA_TOT real,dimension(lev0:lev1,lon0:lon1) :: + xnn2, ! number density of n2 + xno2, ! number density of o2 + xno, ! number density of o + rho, p0ez_mbar, qti,n2 ! ! Note: qti = output electron production, not currently used in tiegcm: ! ! Parameters real,parameter :: + EminV=100., + Mean_mass=4.23e-23, + tab_shape(9)=(/1.,0.,0., 0.,0.1,1., 0.,0.7,1./), + tab_EL(6)=(/0.87, 0.84, 0.883, 0.9, 0.92, 0.92/), + tab_KL(6)=(/4.5e-18,8.1e-18,7.0e-18,7.0e-18,7.0e-18,8.e-18/), + tab_E0(6)=(/13.5, 6.5, 3.5, 2.5, 1.5, 0./), + P_W1(7)=(/37.90, 2.95e-1, 2.03e-1, -1.57e-2, 3.44e-4, 0., 0./), + P_W2(7)=(/156.68, -63.86, 17.33, -2.32, 1.641e-1, -5.83e-3, + 8.21e-5/), + P_W3(7)=(/8.01e-1, 2.09, -2.91e-2, 1.58e-4, 0., 0., 0./), + P_W4(7)=(/14.26, 2.73, -8.55e-2, 2.93e-3, 0., 0., 0./), + P_W5(7)=(/28.81, 32.69, -1.83, 4.02e-2, 0., 0., 0./), + P_W6(7)=(/18.94, 1.32, -8.25e-2, 3.80e-3, -7.06e-5, + 0., 0./), + P_Kd1(3)=(/3.85e-2, 3.03e-2, -4.74e-4/), + P_Kd2(3)=(/4.56e-2, 1.27e-2, -1.15e-4/), + P_Kd3(3)=(/8.90e-2, 3.95e-2, -4.21e-4/), + P_Kd4(3)=(/3.15e-1, 2.19e-2, -3.92e-4/), + E_Kd(4)=(/30., 50., 40., 30./) ! ! exp(-z)*p0 1. ! p0ez_mbar = ---------- * ------------------ ! 1.38e-16*tn o2/32.+o/16.+n2/28. ! ! Where: o2, o, n2 are mass mixing ratios. ! Used in number densities calculation below. ! QIA(:,:,:) = 0. do i=lon0,lon1 do k=lev0,lev1 p0ez_mbar(k,i) = p0*expz(k)/(boltz*tn(k,i)* | (o2(k,i)*rmassinv_o2+ | o1(k,i)*rmassinv_o1+ | (1.-o2(k,i)-o1(k,i))*rmassinv_n2)) n2(k,i) = max(.00001,(1.-o1(k,i)-o2(k,i))) ! n2 mmr enddo enddo ! ! Calculate number densities of O2, O, and N2 from mass mixing ! ratios in f-array. Store in locals xno2, xno, xnn2. ! Also sum total number density rho = o2+o+n2 (gm/cm3) ! do i=lon0,lon1 do k=lev0,lev1 rho(k,i) = (o2(k,i)+o1(k,i)+n2(k,i))*p0ez_mbar(k,i)* + 1.66e-24 xno2(k,i) = o2(k,i)*p0ez_mbar(k,i)/rmass_o2 xno(k,i) = o1(k,i)*p0ez_mbar(k,i)/rmass_o1 xnn2(k,i) = n2(k,i) *p0ez_mbar(k,i)/rmass_n2 enddo enddo ! ! Loop down from top: do k=lev1-1,lev0,-1 ! ! Longitude loop: do i=lon0,lon1 if (efx_lp(i,lat) > 0. .and. lat > 25) then ! ! RANGE (g.cm-2) ALPV=kev_lp(i,lat)*1.e3 do ii=6,1,-1 if (kev_lp(i,lat).gt.2.*tab_E0(ii)) then EL=tab_EL(ii) KL=tab_KL(ii) endif enddo RANGE = Mean_mass/KL/(1.-EL)*(ALPV**(1.-EL)-EminV**(1.-EL)) F_E = ALPV**(1.-EL)/KL/(1.-EL)/RANGE*Mean_mass ! ! R_norm : Normalized atmospheric scattering depth (g.cm-2) if (k .eq. lev1-1) then R_norm(i)=0. else R_norm(i)=R_norm(i)+(rho(k+1,i)+rho(k,i))/2* + (z(k+1,i)-z(k,i))/RANGE endif ! ! LAMBDA : Normalized energy deposition function if (abs(1.-EL)/(1.-EL)*(F_E-R_norm(i)).gt.0.) then LAMBDA=KL*RANGE/Mean_mass/(ALPV-EminV)* + (((1.-EL)*KL*RANGE/Mean_mass* + (F_E-R_norm(i)))**(EL/(1.-EL))) else LAMBDA=0. endif ! ! ETA : Energy deposition rate (eV.cm-2.s-1) ETA=efx_lp(i,lat)*1.e-7/(1.602e-19)/RANGE* + rho(k,i)*LAMBDA ! W(:)=0. do ii=1,7 W(1)=W(1)+P_W1(ii)*(kev_lp(i,lat)/2)**(ii-1) W(2)=W(2)+P_W2(ii)*(kev_lp(i,lat)/2)**(ii-1) W(3)=W(3)+P_W3(ii)*(kev_lp(i,lat)/2)**(ii-1) W(4)=W(4)+P_W4(ii)*(kev_lp(i,lat)/2)**(ii-1) W(5)=W(5)+P_W5(ii)*(kev_lp(i,lat)/2)**(ii-1) W(6)=W(6)+P_W6(ii)*(kev_lp(i,lat)/2)**(ii-1) enddo ! ! QTI : Electron production rate (cm-3.s-1) QTI(k,i)=ETA/W(6) ! ! QIA(i,K,n) : Ion production rate (cm-3.s-1) ! n=1,2,3,4,5 -> N2+,O2+,O+,NO+,N+ ! shape(1)=tab_shape(1)*xnn2(k,i)/(tab_shape(1)*xnn2(k,i)+ + tab_shape(2)*xno2(k,i)+tab_shape(3)*xno(k,i)) shape(2)=tab_shape(5)*xno2(k,i)/(tab_shape(4)*xnn2(k,i)+ + tab_shape(5)*xno2(k,i)+tab_shape(6)*xno(k,i)) shape(3)=tab_shape(9)*xno(k,i)/(tab_shape(7)*xnn2(k,i)+ + tab_shape(8)*xno2(k,i)+tab_shape(9)*xno(k,i)) ! do n=1,4 Kd(n)=0. if (kev_lp(i,lat).gt.E_Kd(n)) then E=E_Kd(n) else E=kev_lp(i,lat) endif select case (n) case (1) do m=1,3 Kd(n)=Kd(n)+P_Kd1(m)*E**(m-1) enddo case (2) do m=1,3 Kd(n)=Kd(n)+P_Kd2(m)*E**(m-1) enddo case (3) do m=1,3 Kd(n)=Kd(n)+P_Kd3(m)*E**(m-1) enddo case (4) do m=1,3 Kd(n)=Kd(n)+P_Kd4(m)*E**(m-1) enddo end select enddo ! QIA(k,i,1)=shape(1)*(ETA/W(1)/(1+Kd(1))+ETA/W(2)/(1+Kd(2))) QIA(k,i,5)=shape(1)*(ETA/W(1)*Kd(1)/(1+Kd(1))+ + ETA/W(2)*Kd(2)/(1+Kd(2))) QIA(k,i,2)=shape(2)*(ETA/W(3)/(1+Kd(3))+ETA/W(4)/(1+Kd(4))) QIA(k,i,3)=shape(2)*(ETA/W(3)*Kd(3)/(1+Kd(3))+ + ETA/W(4)*Kd(4)/(1+Kd(4)))+shape(3)*ETA/W(5) QIA_TOT=QIA(k,i,1)+QIA(k,i,2)+QIA(k,i,3)+QIA(k,i,5) ! if (i .eq. lon0) then ! write(6,"('z(k:lev1,i) = ',/,(6e12.2))")z(k:lev1,i) ! write(6,"('rho(k:lev1,i) = ',/,(6e12.2))")rho(k:lev1,i) ! write(6,"('k,i = ',2i4,'RANGE,R_norm(i) = ',2e12.2)") ! + k,i,RANGE,R_norm(i) ! endif ! if (lat > 30 .and. k.eq.14) then ! write(6,"('k,i,shape(1:3) = ',2i4,3e12.3)")k,i,shape(:) ! write(6,"('efx_lp(i,lat),kev_lp(i,lat) = ',2e12.2)") ! + efx_lp(i,lat),kev_lp(i,lat) ! write(6,"('ETA,F_E,RANGE,LAMBDA,R_norm,rho = ',6e12.2)") ! + ETA,F_E,RANGE,LAMBDA,R_norm(i),rho(k,i) !! write(6,"('xno,xno2,xnn2,rho= ',4e12.3)") !! + xno(k,i),xno2(k,i),xnn2(k,i),rho(k,i) ! write(6,"('W(1:6) = ',6e12.3)")W(:) ! write(6,"('Kd(k,i) = ',4e12.3)")Kd(:) ! write(6,"('QTI(k,i) = ',e12.3)")QTI(k,i) ! write(6,"('QIA(k,i,:) = ',5e12.3)")QIA(k,i,:) ! endif if (QIA_TOT.gt.1.e-20) then do n=1,3 QIA(k,i,n)=QIA(k,i,n)*QTI(k,i)/QIA_TOT enddo QIA(k,i,5)=QIA(k,i,5)*QTI(k,i)/QIA_TOT endif endif enddo enddo ! end subroutine low_proton