#include "defs.h" ! module amie_module ! ! Module used to read data from the AMIE outputs (POT,mean energy, ! and energy flux). ! use params_module,only: nlon,nlat,nlev,nlonp1,nlonp4,nlatp1, | nlevp1,nmlon,nmlonp1,nmlat use addfld_module,only: addfld implicit none ! ! Define parameters for AMIE input data file: integer,parameter :: | jamietyp = 2, ! jamietyp = 1 for Gang Lu and = 2 for Aaron Ridley | mxtimes = 5881, ! maximum number of times of AMIE data per file (Dec 13-16 1min) ! Aaron Ridley's amie files: | mxgdays = 1, ! maximum number of days of AMIE data (10 for Gang, 1 for Aaron) | lonmx = 24, ! maximum number of longitudes of AMIE data (1 less than grid) ! colath.F restricts crit(1) between 15-30 colat so crit(2) is between 30 and 45) ! Will set defaults at crad=15, crit(1)=20, crit(2)=35 so AMIE irrel equatorwards of 55mlat ! Ridley AMIE is 2deg 90-44 or 24 lats + 15 lats extended to 14 mlat where potential is zero | ithtrns = 24, ! Aaron: corresponding to trans lat 42 mlat (ie, 90 to 44 dmlat=2) | ithmx = 46, ! max number of AMIE mlats (90/2=45+1) to equator ! Gang Lu's AMIE files: ! | mxgdays = 10, ! maximum number of days of AMIE data ! | lonmx = 36, ! maximum number of longitudes of AMIE data ! | ithtrns = 30, ! corresponding to trans lat 40-deg ! | ithmx = 55, ! maximum number of latitudes of AMIE data ! Both set jmxm the same | jmxm = 2*ithmx-1 ! maximum number of global latitudes (only one 0mlat) ! lonp1 and latp1 are AMIE dim of lon (37 or lonmx+1) and lat (31 or ithmx+1 where 30 is 90 to 45.67 in dmlat=1.67) integer :: lonp1,latp1, | year(mxtimes),month(mxtimes),day(mxtimes),jday(mxtimes) ! Define AMIE output fields real :: | tiepot(nmlonp1,nmlat),tieekv(nmlonp1,nmlat), | tieefx(nmlonp1,nmlat),potg(nlonp1,0:nlatp1), | efxg(nlonp1,0:nlatp1),ekvg(nlonp1,0:nlatp1) ! defined output AMIE fields in TGCM geographic grid real,dimension(nlonp4,nlat) :: | potg_sech, ekvg_sech, efxg_sech real,dimension(nmlonp1,-2:nlevp1) :: tiepot_sech ! ! Define fields for AMIE input data file: ! electric potential in Volt ! mean energy in KeV ! energy flux in W/m^2 (wrong for U MI AMIE - is mW/m2 or ergs/cm2-s) ! amie_cusplat_nh(sh) and amie_cuspmlt_nh(sh) are ! AMIE cusp latitude and MLT in NH and SH ! amie_hpi_nh(sh) are AMIE hemi-integrated power ! amie_pcp_nh(sh) are AMIE polar-cap potential drop ! Saved AMIE outputs with suffix _amie ! ! real,allocatable,dimension(:,:,:) :: ! (lonp1,latp1,ntimes) ! | amie_pot_nh, amie_pot_sh, amie_ekv_nh, amie_ekv_sh, ! | amie_efx_nh, amie_efx_sh ! real,allocatable,dimension(:,:) :: ! (lonp1,latp1) real, dimension(ithtrns+1,lonmx+1) :: | pot_nh_amie,pot_sh_amie, ekv_nh_amie,ekv_sh_amie, | efx_nh_amie,efx_sh_amie real :: | cusplat_nh_amie, cuspmlt_nh_amie, cusplat_sh_amie, | cuspmlt_sh_amie, hpi_sh_amie, hpi_nh_amie, pcp_sh_amie, | pcp_nh_amie, crad(2), phida(2) ! ! 6/12: added by Emery for binary Ridley AMIE real :: Re,areatot real :: pii,dtr,htr real :: del,xmlt,dmlat,dlatm,dlonm,dmltm,rot real, dimension(ithtrns+1) :: mlatdeg,da real, dimension(lonmx+1) :: mltnh real, dimension(lonmx+1,ithtrns*2+2) :: mlts, lats, | TempPotential, AveEOut, EfluxOut contains !----------------------------------------------------------------------- subroutine init_amie ! ! Called from tgcm.F ! (this is not in init.F to avoid circular dependencies) ! use input_module,only: amienh,amiesh ! ! 6/12 added by Emery for Ridley binary AMIE files character (len=100), dimension(100) :: Lines integer :: iError, i,j ! 6/12: Set AMIE NH or SH grid dimensions lonp1 and latp1 (to go 90 to 40 mlat in dmlat=2) latp1 = ithtrns + 1 lonp1 = lonmx + 1 pii = 4.*atan(1.) dtr = pii/180. htr = pii/12. ! Re=radius of the earth in m Re = 6371.e+3 C **** SET GRID SPACING DLATM, DLONG, DLONM C DMLAT=lat spacing in degrees of AMIE apex grid C 6/12: from amie.nc.out in ~ganglu/amie, have 31 lats w dmlat=1.67 from 90 to 40 mlat C and 37 mlts w dmltm=0.67h from 0 to 24MLT, where lonmx=36 (NOT 37!) C 18*3=54 so 90/54=1.67=dmlat (ithtrns=30 to 40 mlat, ithmx=55 to equator, jmxm=2*ithmx-1) ! dtr = pi/180. DMLAT = 180. / FLOAT(jmxm-1) ! AMIE grid spacing in deg (should be 2 deg) DLATM = DMLAT*dtr DLONM = 2.*pii/FLOAT(lonmx) DMLTM = 24./FLOAT(lonmx) ! AMIE grid spacing in MLT (should be 1h) if (len_trim(amienh) > 0) then write(6,"('Reading AMIENH file ',a)") trim(amienh) Lines(2) = trim(amienh) endif if (len_trim(amiesh) > 0) then write(6,"('Reading AMIESH file ',a)") trim(amiesh) Lines(3) = trim(amiesh) else Lines(3) = "mirror" endif Lines(1) = "#AMIEFILES" ! Lines(2) = "/home/emery/archive/amie_umich/b20061213-16.wmds" ! Lines(3) = "mirror" Lines(4) = "" Lines(5) = "" Lines(6) = "" Lines(7) = "#DEBUG" Lines(8) = "2" Lines(9) = "0" Lines(10) = "" Lines(11) = "#END" ! amienh = trim(Lines(2)) ! amiesh = '' write(*,*) "Calling IE_set_inputs" call IE_set_inputs(Lines) write(*,*) "=> Setting up Ionospheric Electrodynamics" call IE_Initialize(iError) write(*,*) "Setting nmlts and nlats" ! Aaron Ridley says the original AMIE comes as 1min NH arrays of 25 MLTs (0-24 - NOT 1-25!!!) ! every 2 deg from 90N to 44N (24!)- However, the code interpolates to whatever grid is ! set below as lats and mlts for input with dimensions set in IO_SetnMLT(LAT)s ! Can find the dimensions of AMIE in readAMIEoutput.f90 call IO_SetnMLTs(lonmx+1) call IO_SetnLats(ithtrns*2+2) areatot = 0. write(*,*) "dmlat dmlon=",dmlat,dlonm/dtr do j=1,latp1 ! In 'mirror', the min/max potentials are flipped from NH to SH ! Want lats from 90 to 40 mlat in both hem mlatdeg(j) = 90.0 - (j-1)*dmlat ! cos(90deg)=0 so pole is ignored with da(1)=0 da(j) = Re*cos(mlatdeg(j)*dtr)*dlonm*(Re*dmlat*dtr) areatot = areatot + da(j)*lonmx do i=1,lonp1 lats(i,j+latp1) = 90.0 - (j-1)*dmlat lats(i,j) = -90.0 + (j-1)*dmlat mltnh(i) = i-1 mlts(i,j) = i-1 mlts(i,j+latp1) = i-1 enddo enddo write (*,*) "areatot da(j)=",areatot,da write(*,*) "Setting grid" call IO_SetGrid(mlts,lats, iError) end subroutine init_amie !----------------------------------------------------------------------- !----------------------------------------------------------------------- subroutine getamie(iyear,iday,iutsec,amie_ibkg,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 ModKind ! 6/12: for AMIE U MI use hist_module,only: modeltime use input_module,only: amiesh,amienh use cons_module,only: pi,ylonm,ylatm use magfield_module,only: im,jm,dim,djm ! use dynamo_module,only: mag2geo ! 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 input_module,only: ! from user input | byimf ! By component of IMF (nT) (e.g., 0. - used in default for phid,phin) ! ! Args: integer,intent(in) :: iyear,iday,iutsec,amie_ibkg,iprint ! ! Local: real :: | potm(lonp1,jmxm),efxm(lonp1,jmxm),ekvm(lonp1,jmxm), | alat(jmxm),alon(lonp1),alatm(jmxm),alonm(lonp1) 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, amie_debug,amieday integer :: iboxcar,inpt real :: model_ut, denoma, f1, f2, fmin, fmax, maxefx ! real :: del,xmlt,dmlat,dlatm,dlonm,dmltm,rot,dtr ! 6/12: added by Emery for Ridley binary AMIE files integer, dimension(7) :: itime real (Real8_) :: rtime ! real :: rtime integer :: iError, imlt,ilat,imltnsv,imltxsv,ilatnsv,ilatxsv,jl real :: value, maxkV, minkV, hpi real :: cps,cpn,latnx,latnn,latsx,latsn,potsx,potsn,potnx,potnn real :: mltnx,mltnn,mltsx,mltsn ! 08/12 bae: added for calccloc real :: model_ctpoten(2),mltd,mltn integer :: ifbad(2),byloc real, dimension(lonmx+1,ithtrns+1) :: potS,potN ! ! 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/ ! ! 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. ! if (iprint > 0) then write(6,"(/,72('-'))") write(6,"('GETAMIE:')") write(6,"('Modeltime(1:4)= ',4i5)") modeltime(1:4) write(6,"('Initial requested iyear=',i4,' iday=',i3,' iutsec=', | i10)") iyear,iday,iutsec endif maxefx=50. ! ! Check for leap year nonleap = int(mod(real(iyear),4.)) ! nonleap year /= 0 ! if (iyear.eq.2000) nonleap = 1 ! ! Check times if (modeltime(3).lt.0.or.modeltime(3).gt.59) then write(6,"('>>> getamie: bad imin=',3i4)") modeltime(3) stop 'getamie' endif if (modeltime(4).lt.0.or.modeltime(4).gt.59) then write(6,"('>>> getamie: bad isec=',3i4)") modeltime(4) stop 'getamie' endif ! model_ut = (modeltime(1)-amieday)*24. + float(iutsec)/3600. model_ut = float(iutsec)/3600. ! 6/12: added by Emery for Ridley binary AMIE ! Assume AMIE is same as model time where modeltime(1) = daynumber of year! itime(1) = iyear do i=1,12 if (nonleap == 0) then if (ndmon_lp(i) <= modeltime(1)) then itime(2) = i+1 endif else if (ndmon_nl(i) <= modeltime(1)) then itime(2) = i+1 endif endif enddo if (nonleap == 0) then itime(3) = modeltime(1)-ndmon_lp(itime(2)-1) else itime(3) = modeltime(1)-ndmon_nl(itime(2)-1) endif itime(4) = modeltime(2) itime(5) = modeltime(3) itime(6) = modeltime(4) itime(7) = 0 ! Calling time conversion for Ridley AMIE ! write(*,*) nonleap,itime call time_int_to_real(itime, rtime) ! write(*,*) "rtime =",rtime call IO_SetTime(rtime) ! write(*,*) "after SetTime rtime =",rtime call get_AMIE_values(rtime) ! write(*,*) "after get_AMIE rtime =",rtime ! ! AveEOut in keV (from ~1-9 keV) call IO_GetAveE(AveEOut, iError) if (iError /= 0) then write(*,*) "Error : ",iError else ! write(*,*) "AveEOut itime=",itime ! write(*,*) AveEOut endif ! EfluxOut in mW/m2? (~24 mW/m2 max 06349 0030UT - has to be mW/m2!) ! 6/12: Emery - U MI AMIE EFluxOut is in erg/cm2-s or mW/m2 (NOT in W/m2!) call IO_GetEFlux(EfluxOut, iError) if (iError /= 0) then write(*,*) "Error : ",iError else ! write(*,*) "EfluxOut itime=",itime ! write(*,*) EfluxOut endif call IO_GetPotential(TempPotential, iError) if (iError /= 0) then write(*,*) "Error : ",iError else if (itime(4).eq.0 .and. itime(5).eq.2 .and. itime(3).eq.13) then write(*,*) "i,rtime=",itime,rtime ! mlts = 10x 1-25 (same as 25 IEi_HavenMLTS) ! write(*,*) "mlts=",mlts ! lats = 25x 86,82,78,74,70,66,62,58,54,50 (10 mlats - diff from 39 IEi_HavenLats - expected 88,87-50) ! write(*,*) "lats=",lats ! TempPotential in Volts write(*,*) "SH pot in V, -90 to -42 mlat" write(*,*) ((TempPotential(i,j),i=1,lonmx),j=1,latp1) write(*,*) "NH pot in V, 90 to 42 mlat" write(*,*) ((TempPotential(i,j),i=1,lonmx),j=latp1+1,latp1*2) endif ! get rid of outragous flux values > maxefx ! Calculate min/max CP and hpi_s,nh_amie in GW do j=1,2 jl = (j-1)*latp1 hpi = 0. maxkV = 0. minkV = 0. do imlt=1,lonp1 do ilat=1,latp1 ! write(*,*) jl,imlt,ilat if (jl < 1) then potS(imlt,ilat) = TempPotential(imlt,ilat+jl) pot_sh_amie(imlt,ilat) = TempPotential(imlt,ilat+jl) ekv_sh_amie(imlt,ilat) = AveEOut(imlt,ilat+jl) efx_sh_amie(imlt,ilat) = EfluxOut(imlt,ilat+jl) if (imlt maxkV) then maxkV = TempPotential(imlt,ilat+jl) ilatxsv = ilat+jl imltxsv = imlt endif if (TempPotential(imlt,ilat+jl) < minkV) then minkV = TempPotential(imlt,ilat+jl) ilatnsv = ilat+jl imltnsv = imlt endif enddo enddo if (jl < 1) then cps = (maxkV-minkV)*1.e-3 potsx = maxkV*1.e-3 potsn = minkV*1.e-3 latsx = mlatdeg(ilatxsv) latsn = mlatdeg(ilatnsv) mltsx = mlts(imltxsv,1) mltsn = mlts(imltnsv,1) hpi_sh_amie = hpi*1.e-12 else cpn = (maxkV-minkV)*1.e-3 potnx = maxkV*1.e-3 potnn = minkV*1.e-3 latnx = lats(1,ilatxsv) latnn = lats(1,ilatnsv) mltnx = mlts(imltxsv,1) mltnn = mlts(imltnsv,1) ! Convert from mW/m2 * m2(da) * 1.e-12GW/mW hpi_nh_amie = hpi*1.e-12 endif ! write(*,*) "min,max tot CP = ",minkV*1.e-3,mlts(imltnsv,ilatnsv),lats(imltnsv,ilatnsv), & ! maxkV*1.e-3,mlts(imltxsv,ilatxsv),lats(imltnsv,ilatnsv), (maxkV-minkV)*1.e-3 write(*,*) "min,max tot CP HP= ",minkV*1.e-3,maxkV*1.e-3, | (maxkV-minkV)*1.e-3,hpi*1.e-12 enddo ! for j=1,2 and jl endif ! for if can read potential w/o error ! Assign pcp_s,nh_amie as cps,n in kV of Polar Cap Potential drop pcp_sh_amie = cps pcp_nh_amie = cpn ! Calculate theta0 (cusplat) and set MLT for cusp also from calccloc? ! Set defaults at crad=15deg, crit1=crad+5 and crit2=crad+20 in colath.F cusplat_nh_amie = 75. cusplat_sh_amie = 75. cuspmlt_nh_amie = 11.75 cuspmlt_sh_amie = 11.75 ! ! The OLTMAX latitude also defines the co-latitude theta0, which in ! turn determines crit1(+2.5deg) and crit2(-12.5deg) which are used ! in TIE-GCM as the boundaries of the polar cap and the region of ! influence of the high-lat potential versus the low-lat dynamo potential ! Define this latitude to be between 70 and 77.5 degrees (then 60-75 mlat!) crad(1) = (90.-cusplat_sh_amie)*pi/180. crad(2) = (90.-cusplat_nh_amie)*pi/180. phida(1) = (cuspmlt_sh_amie - 12.) * pi / 12. phida(2) = (cuspmlt_nh_amie - 12.) * pi / 12. ! 08/12 bae: Recalculate convection (and auroral) parameters from calccloc call calccloc (mlatdeg,mlatdeg,ithtrns+1,mltnh,mltnh, | lonmx+1,potS,potN,model_ctpoten,ifbad) ! 08/12 bae: Skip resetting crad,phida since not used in aurora.F anymore pcp_sh_amie = model_ctpoten(1) pcp_nh_amie = model_ctpoten(2) ! 08/12 bae: If both hemispheres cannot be calculated (ifbad>0 in both), set defaults if (ifbad(1)>0 .and. ifbad(2)>0) then offc(1) = 0. offc(2) = 0. offa(1) = 0. offa(2) = 0. dskofc(1) = 0. dskofc(2) = 0. dskofa(1) = -2.5*pi/180. dskofa(2) = -2.5*pi/180. theta0(1) = 12.*pi/180. theta0(2) = 12.*pi/180. rrad(1) = theta0(1) + 2.*pi/180. rrad(2) = theta0(2) + 2.*pi/180. ! Use parameterization defaults for phid (phid(MLT)=9.39 +/- 0.21By - 12) ! and phin (phin(MLT)=23.50 +/- 0.15By - 12) mltd = 9.39 + 0.21*byloc phid(1) = (mltd-12.) * pi/12. mltd = 9.39 - 0.21*byloc phid(2) = (mltd-12.) * pi/12. mltn = 23.50 + 0.15*byloc phin(1) = (mltn-12.) * pi/12. mltn = 23.50 - 0.15*byloc phin(2) = (mltn-12.) * pi/12. endif ! ! gl - 14/07/2002 changed from Emery's original amieterp.f ! ! Change mag coords so latitude extends to the equator in into the opposite ! pole, and so MLT is converted to apex longitues between -180 and +180. ! Put EPOTDUP, EMKEV and WPM2 into S. Hem and N. Hem alike in EPOTM, EKEVM, ! and EFLXM, dimensioned (IMXMP,JMNH) where MLAT goes from -90 to ~-2 ! for ihsoln=1, and from +90 to ~+2 for ihsoln=2 since GRDSET is for ! -90 to the equator. (Will work for +90 to the equator also.) ! AMIE grid goes from -90 to ~-40 for SH and +90 to ~+40 for NH (6/12: w dmlat=1.67deg)) ! TGCM grid goes from -87.5 to -2.5 for SH and +2.5 to +87.5 for NH (6/12: w dglat=5deg) ! Hence, at very end, need to reverse order of NH in whole globe arrays. ! ! do i=1,lonp1 ! Initialize arrays around equator ! do j=latp1+1,ithmx ! potm(i,j) = 0. ! potm(i,jmxm+1-j) = 0. ! ekvm(i,j) = ekv_sh_amie(i,latp1) ! ekvm(i,jmxm+1-j) = ekv_nh_amie(i,latp1) ! efxm(i,j) = 0. ! efxm(i,jmxm+1-j) = 0. ! enddo ! Put in AMIE arrays from pole to JMFT ! do j=1,latp1 ! potm(i,j) = pot_sh_amie(i,j) ! potm(i,latp1+1-j) = pot_nh_amie(i,j) ! ekvm(i,j) = ekv_sh_amie(i,j) ! ekvm(i,latp1+1-j) = ekv_nh_amie(i,j) ! efxm(i,j) = efx_sh_amie(i,j) ! efxm(i,latp1+1-j) = efx_nh_amie(i,j) ! enddo ! enddo ! mlongitude starts from 180 degree rot = 250./15. - iutsec/3600. dmltm = 24./float(lonmx) do i=1,lonp1 xmlt = (i-1)*dmltm - rot + 24. xmlt = AMOD(xmlt,24.) m = IFIX(xmlt/dmltm + 1.01) mp1 = m + 1 IF (mp1 > lonp1) mp1 = 2 del = xmlt - (m-1)*dmltm C Initialize arrays around equator do j=latp1+1,ithmx potm(i,j) = 0. potm(i,jmxm+1-j) = 0. ekvm(i,j) = (1.-del)*ekv_sh_amie(m,latp1) + | del*ekv_sh_amie(mp1,latp1) ekvm(i,jmxm+1-j) = (1.-del)*ekv_nh_amie(m,latp1) + | del*ekv_nh_amie(mp1,latp1) efxm(i,j) = 0. efxm(i,jmxm+1-j) = 0. enddo C Put in AMIE arrays from pole to latp1 do j=1,latp1 potm(i,j) = (1.-del)*pot_sh_amie(m,j) + del*pot_sh_amie(mp1,j) potm(i,jmxm+1-j) = (1.-del)*pot_nh_amie(m,j) + | del*pot_nh_amie(mp1,j) ekvm(i,j) = (1.-del)*ekv_sh_amie(m,j) + del*ekv_sh_amie(mp1,j) ekvm(i,jmxm+1-j) = (1.-del)*ekv_nh_amie(m,j) + | del*ekv_nh_amie(mp1,j) efxm(i,j) = (1.-del)*efx_sh_amie(m,j) + del*efx_sh_amie(mp1,j) efxm(i,jmxm+1-j) = (1.-del)*efx_nh_amie(m,j) + | del*efx_nh_amie(mp1,j) enddo enddo C Set up coeffs to go between EPOTM(IMXMP,JMNH) and TIEPOT(IMAXM,JMAXMH) C C **** SET GRID SPACING DLATM, DLONG, DLONM C DMLAT=lat spacing in degrees of AMIE apex grid C 6/12: from amie.nc.out in ~ganglu/amie, have 31 lats w dmlat=1.67 from 90 to 40 mlat C and 37 mlts w dmltm=0.67h from 0 to 24MLT, where lonmx=36 (NOT 37!) C 18*3=54 so 90/54=1.67=dmlat (ithtrns=30 to 40 mlat, ithmx=55 to equator, jmxm=2*ithmx-1) ! dtr = pi/180. ! DMLAT = 180. / FLOAT(jmxm-1) ! DLATM = DMLAT*dtr ! DLONM = 2.*PI/FLOAT(lonmx) ! DMLTM = 24./FLOAT(lonmx) C **** C **** SET ARRAY YLATM (LATITUDE VALUES FOR GEOMAGNETIC GRID C **** alatm(1) = -PI/2. alat(1) = -90. alatm(jmxm) = PI/2. alat(jmxm) = 90. do i = 2,ithmx alat(i) = alat(i-1)+dlatm/dtr alat(jmxm+1-i) = alat(jmxm+2-i)-dlatm/dtr alatm(i) = alatm(i-1)+dlatm alatm(jmxm+1-i) = alatm(jmxm+2-i)-dlatm enddo alon(1) = -pi/dtr alonm(1) = -pi do i=2,lonp1 alon(i) = alon(i-1) + dlonm/dtr alonm(i) = alonm(i-1) + dlonm enddo ! ylatm and ylonm are arrays of latitudes and longitudes of the ! distored magnetic grids in radian - from consdyn.h ! Convert from apex magnetic grid to distorted magnetic grid ! ! Allocate workspace for regrid routine rgrd2.f: lw = nmlonp1+nmlat+2*nmlonp1 if (.not. allocated(w)) allocate(w(lw),stat=ier) if (ier /= 0) write(6,"('>>> horizontal_interp: error allocating', | ' w(lw): lw=',i6,' ier=',i4)") lw,ier liw = nmlonp1 + nmlat if (.not. allocated(iw)) allocate(iw(liw),stat=ier) if (ier /= 0) write(6,"('>>> horzontal_interp: error allocating', | ' iw(liw): liw=',i6,' ier=',i4)") liw,ier intpol(:) = 1 ! linear (not cubic) interp in both dimensions if (alatm(1) > ylatm(1)) alatm(1) = ylatm(1) if (alatm(jmxm) < ylatm(nmlat)) alatm(jmxm) = ylatm(nmlat) if (alonm(1) > ylonm(1)) alonm(1) = ylonm(1) if (alonm(lonp1) < ylonm(nmlonp1)) alonm(lonp1) = ylonm(nmlonp1) call rgrd2(lonp1,jmxm,alonm,alatm,potm,nmlonp1,nmlat, | ylonm,ylatm,tiepot,intpol,w,lw,iw,liw,ier) call rgrd2(lonp1,jmxm,alonm,alatm,ekvm,nmlonp1,nmlat, | ylonm,ylatm,tieekv,intpol,w,lw,iw,liw,ier) call rgrd2(lonp1,jmxm,alonm,alatm,efxm,nmlonp1,nmlat, | ylonm,ylatm,tieefx,intpol,w,lw,iw,liw,ier) C Convert from TGCM distorted magnetic grid to geographic one ! CALL GRDTRP(tiepot(1,1),potg(1,0),im(1,0),jm(1,0), ! | dim(1,0),djm(1,0),IMAXGP,IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM) ! CALL GRDTRP(tieekv(1,1),ekvg(1,0),im(1,0),jm(1,0), ! | DIM(1,0),DJM(1,0),IMAXGP,IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM) ! CALL GRDTRP(tieefx(1,1),efxg(1,0),im(1,0),jm(1,0), ! | DIM(1,0),DJM(1,0),IMAXGP,IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM) call mag2geo(tiepot(1,1),potg(1,0),im(1,0),jm(1,0), | dim(1,0),djm(1,0),nlonp1,nmlonp1,nlon,nlat+2,nmlon,nmlat) call mag2geo(tieekv(1,1),ekvg(1,0),im(1,0),jm(1,0), | dim(1,0),djm(1,0),nlonp1,nmlonp1,nlon,nlat+2,nmlon,nmlat) call mag2geo(tieefx(1,1),efxg(1,0),im(1,0),jm(1,0), | dim(1,0),djm(1,0),nlonp1,nmlonp1,nlon,nlat+2,nmlon,nmlat) C **** C **** INSERT PERIODIC POINTS C **** DO j = 1,nlat ekvg(nlonp1,j) = ekvg(1,j) efxg(nlonp1,j) = efxg(1,j) potg(nlonp1,j) = potg(1,j) ENDDO ! temp - save the AMIE output in geographic grid amie_debug = 1 if (amie_debug > 0) then ! do j=1,nmlat ! do i=1,nmlonp1 ! tiepot_sech(i,:) = tiepot(i,j) ! enddo ! call addfld('TIEPOT',' ',' ',tiepot_sech(:,1:nlevp1), ! | 'mlon',1,nmlonp1,'mlev',1,nlevp1,j) ! enddo ! tiepot(nmlonp1,nmlat) ! call addfsech_ij('TIEPOT',' ',' ',tiepot,1,nmlonp1,1,nmlat) ! call addfsech_ij('TIEEKV',' ',' ',tieekv,1,nmlonp1,1,nmlat) ! call addfsech_ij('TIEEFX',' ',' ',tieefx,1,nmlonp1,1,nmlat) call addfld('TIEPOT',' ',' ',tiepot,'mlon',1,nmlonp1, | 'mlat',1,nmlat,0) call addfld('TIEEKV',' ',' ',tieekv,'mlon',1,nmlonp1, | 'mlat',1,nmlat,0) call addfld('TIEEFX',' ',' ',tieefx,'mlon',1,nmlonp1, | 'mlat',1,nmlat,0) do i=1,nlon potg_sech(i+2,:) = potg(i,1:nlat) ekvg_sech(i+2,:) = ekvg(i,1:nlat) efxg_sech(i+2,:) = efxg(i,1:nlat) enddo do i=1,2 potg_sech(i,:) = potg_sech(nlon+i,:) potg_sech(nlon+i+2,:) = potg_sech(i,:) ekvg_sech(i,:) = ekvg_sech(nlon+i,:) ekvg_sech(nlon+i+2,:) = ekvg_sech(i,:) efxg_sech(i,:) = efxg_sech(nlon+i,:) efxg_sech(nlon+i+2,:) = efxg_sech(i,:) enddo ! call addfsech_ij('POTG',' ',' ',potg_sech,1,nlonp4,1,nlat) ! call addfsech_ij('EKVG',' ',' ',ekvg_sech,1,nlonp4,1,nlat) ! call addfsech_ij('EFXG',' ',' ',efxg_sech,1,nlonp4,1,nlat) call addfld('POTG',' ',' ',potg_sech,'lon',1,nlonp4, | 'lat',1,nlat,0) call addfld('EKVG',' ',' ',ekvg_sech,'lon',1,nlonp4, | 'lat',1,nlat,0) call addfld('EFXG',' ',' ',efxg_sech,'lon',1,nlonp4, | 'lat',1,nlat,0) endif ! if (iprint > 0) then write(6,"('getamie: AMIE data interpolated to date and time')") write(6,"(72('-'),/)") endif end subroutine getamie !------------------------------------------------------------------- subroutine boxcar_ave(x,y,lon,lat,mtime,itime,ibox) ! ! perform boxcar average ! ! Args: integer,intent(in) :: lon,lat,mtime,itime,ibox real,dimension(lon,lat,mtime), intent(in) :: x real,dimension(lon,lat), intent(out) :: y ! Local: integer :: i, iset, iset1 ! iset = itime - ibox/2 if (iset < 1) iset = 1 iset1 = iset + ibox if (iset1 > mtime) then iset1 = mtime iset = iset1 - ibox endif ! write(6,"('boxcar_ave: mtime,itime,ibox',3i5)") ! | mtime,itime,ibox ! y(:,:) = 0. do i=iset,iset1 y(:,:) = y(:,:) + x(:,:,i) enddo if (ibox > 0) y(:,:) = y(:,:)/ibox ! end subroutine boxcar_ave !----------------------------------------------------------------------- ! BOP ! !IROUTINE: mag2geo ! !INTERFACE: subroutine mag2geo(am,ag,im,jm,dim,djm,lg,lm,nlong,nlatg,nlonm, | nlatm) ! !USES: ! !DESCRIPTION: ! Transform field am on geomagnetic grid to ag on the geographic grid using ! indices im,jm and weights am. Return field ag on geographic grid. ! ! !PARAMETERS: integer,intent(in) :: lg,lm,nlong,nlatg,nlonm,nlatm integer,intent(in) :: im(lg,*),jm(lg,*) real,intent(in) :: am(lm,*),dim(lg,*),djm(lg,*) ! !RETURN VALUE: real,intent(out) :: ag(lg,*) ! ! !REVISION HISTORY: ! 05.03.11 ! ! EOP ! ! Local: integer :: ig,jg ! do jg=1,nlatg do ig=1,nlong ag(ig,jg) = | am(im(ig,jg) ,jm(ig,jg)) *(1.-dim(ig,jg))*(1.-djm(ig,jg))+ | am(im(ig,jg)+1,jm(ig,jg)) * dim(ig,jg) *(1.-djm(ig,jg))+ | am(im(ig,jg) ,jm(ig,jg)+1)*(1.-dim(ig,jg))*djm(ig,jg)+ | am(im(ig,jg)+1,jm(ig,jg)+1)* dim(ig,jg) *djm(ig,jg) enddo ! ig=1,nlong enddo ! jg=1,nlatg end subroutine mag2geo !----------------------------------------------------------------------- end module amie_module