13,18c13 < #ifdef SUN < #include "netcdf3.inc" < #else < #include "netcdf.inc" < #endif < ! --- > 21,26c16,34 < | mxgdays = 10, ! maximum number of days of AMIE data < | mxtimes = 5881, ! maximum number of times of AMIE data per day < | ithtrns = 30, ! corresponding to trans lat 40-deg < | ithmx = 55, ! maximum number of latitudes of AMIE data < | jmxm = 2*ithmx-1, ! maximum number of global latitudes < | lonmx = 36 ! maximum number of longitudes of AMIE data --- > | 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) 42c50 < ! energy flux in W/m^2 --- > ! energy flux in W/m^2 (wrong for U MI AMIE - is mW/m2 or ergs/cm2-s) 49,52c57,61 < 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,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) :: 55,59d63 < real, allocatable,dimension(:) :: ! (ntimes) < | amie_cusplat_nh, amie_cuspmlt_nh, amie_hpi_nh, < | amie_pcp_nh, amie_nh_ut, < | amie_cusplat_sh, amie_cuspmlt_sh, amie_hpi_sh, < | amie_pcp_sh, amie_sh_ut 64a69,78 > ! 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 > > 73a88,111 > ! 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) > > 76c114 < call rdamie_nh --- > Lines(2) = trim(amienh) 80,111c118 < call rdamie_sh < endif < end subroutine init_amie < !----------------------------------------------------------------------- < subroutine rdamie_nh < ! < ! Read AMIE data for the northern hemisphere from amienh < ! < use input_module,only: amienh < ! < ! 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_pot,idv_ekv, < | idv_efx,idv_cusplat,idv_cuspmlt,idv_hpi,idv_pcp < integer :: i ! TEMP for write-out of UT,HP,CP < ! < ! Acquire mss file: < dskfile = ' ' < call getfile(amienh,dskfile) < ! < write(6,"(/,72('-'))") < write(6,"('RDAMIE_NH: read AMIE data for northern hemisphere:')") < ! < ! Open netcdf file: < istat = nf_open(dskfile,NF_NOWRITE,ncid) < if (istat /= NF_NOERR) then < write(6,"(/,'>>> rdamie_nh: error opening netcdf NH AMIE data ', < | 'diurnal file ',a)") trim(dskfile) < stop 'rdamie_nh' --- > Lines(3) = trim(amiesh) 113c120 < write(6,"('rdamie_nh: opened file ',a)") trim(dskfile) --- > Lines(3) = "mirror" 115,287c122,178 < ! < ! Get AMIE grid dimension: < istat = nf_inq_dimid(ncid,'lon',id_lon) < istat = nf_inq_dimlen(ncid,id_lon,lonp1) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting AMIE longitude dimension') < < istat = nf_inq_dimid(ncid,'lat',id_lat) < istat = nf_inq_dimlen(ncid,id_lat,latp1) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting AMIE latitude dimension') < 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, < | 'rdamie_nh: Error getting time dimension') < write(6,"('rdamie_nh: ntimes=',i3)") ntimes < ! < ! Search for requested AMIE output fields < istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim) < ! < ! Get 1-D AMIE 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: < ! real,allocatable,dimension(:,:,:) :: ! (ntimes) < ! | amie_cusplat_nh, amie_cuspmlt_nh, amie_hpi_nh, < ! | amie_pcp_nh, amie_nh_ut, < ! < if (.not. allocated(amie_nh_ut)) < | allocate(amie_nh_ut(ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_nh: error allocating', < | ' amie_nh_ut: ntimes=',i3)")ntimes < if (.not. allocated(amie_cusplat_nh)) < | allocate(amie_cusplat_nh(ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_nh: error allocating', < | ' amie_cusplat_nh: ntimes=',i3)")ntimes < if (.not. allocated(amie_cuspmlt_nh)) < | allocate(amie_cuspmlt_nh(ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_nh: error allocating', < | ' amie_cuspmlt_nh: ntimes=',i3)")ntimes < if (.not. allocated(amie_hpi_nh)) < | allocate(amie_hpi_nh(ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_nh: error allocating', < | ' amie_hpi_nh: ntimes=',i3)")ntimes < if (.not. allocated(amie_pcp_nh)) < | allocate(amie_pcp_nh(ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_nh: error allocating', < | ' amie_pcp_nh: 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 NH AMIE UT id') < istat = nf_get_var_double(ncid,idv_ut,amie_nh_ut) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE variable ut') < ! < ! Get HPI < istat = nf_inq_varid(ncid,'hpi',idv_hpi) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE hpi id') < istat = nf_get_var_double(ncid,idv_hpi,amie_hpi_nh) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE variable hpi') < ! < ! Get PCP < istat = nf_inq_varid(ncid,'pcp',idv_pcp) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE pcp id') < istat = nf_get_var_double(ncid,idv_pcp,amie_pcp_nh) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE variable pcp') < ! 8/12 TEMP bae: < write (6,"(1x,'i,ut(hr,min),hpi,pcp=',i3,4e12.3)") < | (i,amie_nh_ut(i),amie_nh_ut(i)*60.,amie_hpi_nh(i), < | amie_pcp_nh(i),i=1,30) < ! < ! Get cusplat < istat = nf_inq_varid(ncid,'cusplat',idv_cusplat) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE cusplat id') < istat = nf_get_var_double(ncid,idv_cusplat,amie_cusplat_nh) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE variable cusplat') < ! < ! Get cuspmlt < istat = nf_inq_varid(ncid,'cuspmlt',idv_cuspmlt) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE cusplat id') < istat = nf_get_var_double(ncid,idv_cuspmlt,amie_cuspmlt_nh) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE variable cuspmlt') < ! < ! Allocate 2-d fields: < ! real,allocatable,dimension(:,:,:) :: ! (lonp1,latp1) < ! | pot_nh_amie, ekv_nh_amie, efx_nh_amie < ! < if (.not. allocated(pot_nh_amie)) < | allocate(pot_nh_amie(lonp1,latp1),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_sh: error allocating', < | ' pot_nh_amie: lonp1=',i3,' latp1=',i3)")lonp1,latp1 < if (.not. allocated(ekv_nh_amie)) < | allocate(ekv_nh_amie(lonp1,latp1),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_sh: error allocating', < | ' ekv_nh_amie: lonp1=',i3,' latp1=',i3)")lonp1,latp1 < if (.not. allocated(efx_nh_amie)) < | allocate(efx_nh_amie(lonp1,latp1),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_sh: error allocating', < | ' efx_nh_amie: lonp1=',i3,' latp1=',i3)")lonp1,latp1 < ! < ! Allocate 3-d fields: < ! real,allocatable,dimension(:,:,:) :: ! (lonp1,latp1,ntimes) < ! | amie_pot_nh, amie_pot_sh, amie_ekv_nh, amie_ekv_sh, < ! | amie_efx_nh, amie_efx_sh < ! < if (.not. allocated(amie_pot_nh)) < | allocate(amie_pot_nh(lonp1,latp1,ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_nh: error allocating', < | ' amie_pot_nh: lonp1=',i3,' latp1=',i3,' ntimes=',i3)") < | lonp1,latp1,ntimes < if (.not. allocated(amie_ekv_nh)) < | allocate(amie_ekv_nh(lonp1,latp1,ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_nh: error allocating', < | ' amie_ekv_nh: lonp1=',i3,' latp1=',i3,' ntimes=',i3)") < | lonp1,latp1,ntimes < if (.not. allocated(amie_efx_nh)) < | allocate(amie_efx_nh(lonp1,latp1,ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_nh: error allocating', < | ' amie_efx_nh: lonp1=',i3,' latp1=',i3,' ntimes=',i3)") < | lonp1,latp1,ntimes < ! < ! Get 3-D AMIE fields (lon,lat,ntimes) < ! < ! AMIE electric potential < istat = nf_inq_varid(ncid,'pot',idv_pot) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE electric potential id') < istat = nf_get_var_double(ncid,idv_pot,amie_pot_nh) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE variable pot') < ! < ! AMIE mean energy < istat = nf_inq_varid(ncid,'ekv',idv_ekv) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE mean energy id') < istat = nf_get_var_double(ncid,idv_ekv,amie_ekv_nh) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE variable ekv') < ! < ! AMIE energy flux < istat = nf_inq_varid(ncid,'efx',idv_efx) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE energy flux id') < istat = nf_get_var_double(ncid,idv_efx,amie_efx_nh) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_nh: Error getting NH AMIE variable efx') < ! < ! Close the file: < istat = nf_close(ncid) < write(6,"('Completed read from NH AMIE data file ',a)") < | trim(dskfile) < write(6,"(72('-'),/)") < end subroutine rdamie_nh --- > > 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 289,483c180 < subroutine rdamie_sh < ! < ! Read AMIE data for the southern hemisphere from amiesh < ! < use input_module,only: amiesh < ! < ! 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_pot,idv_ekv, < | idv_efx,idv_cusplat,idv_cuspmlt,idv_hpi,idv_pcp < ! < ! Acquire mss file: < dskfile = ' ' < call getfile(amiesh,dskfile) < ! < write(6,"(/,72('-'))") < write(6,"('RDAMIE_SH: read AMIE data for southern hemisphere:')") < ! < ! Open netcdf file: < istat = nf_open(dskfile,NF_NOWRITE,ncid) < if (istat /= NF_NOERR) then < write(6,"(/,'>>> rdamie_nh: error opening netcdf SH AMIE data ', < | 'diurnal file ',a)") trim(dskfile) < stop 'rdamie_sh' < endif < ! < ! Get AMIE grid dimension: < istat = nf_inq_dimid(ncid,'lon',id_lon) < istat = nf_inq_dimlen(ncid,id_lon,lonp1) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting AMIE longitude dimension') < istat = nf_inq_dimid(ncid,'lat',id_lat) < istat = nf_inq_dimlen(ncid,id_lat,latp1) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting AMIE latitude dimension') < ! < ! 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, < | 'rdamie_sh: Error getting time dimension') < ! < ! Search for requested AMIE output fields < istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim) < ! < ! Get 1-D AMIE 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: < ! real,allocatable,dimension(:,:,:) :: ! (ntimes) < ! | amie_cusplat_sh, amie_cuspmlt_sh, amie_hpi_sh, < ! | amie_pcp_sh, amie_sh_ut, < ! < if (.not. allocated(amie_sh_ut)) < | allocate(amie_sh_ut(ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_sh: error allocating', < | ' amie_sh_ut: ntimes=',i3)")ntimes < if (.not. allocated(amie_cusplat_sh)) < | allocate(amie_cusplat_sh(ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_sh: error allocating', < | ' amie_cusplat_sh: ntimes=',i3)")ntimes < if (.not. allocated(amie_cuspmlt_sh)) < | allocate(amie_cuspmlt_sh(ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_sh: error allocating', < | ' amie_cuspmlt_sh: ntimes=',i3)")ntimes < if (.not. allocated(amie_hpi_sh)) < | allocate(amie_hpi_sh(ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_sh: error allocating', < | ' amie_hpi_sh: ntimes=',i3)")ntimes < if (.not. allocated(amie_pcp_sh)) < | allocate(amie_pcp_sh(ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_sh: error allocating', < | ' amie_pcp_sh: ntimes=',i3)")ntimes < ! < ! Get ut < istat = nf_inq_varid(ncid,'ut',idv_ut) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE UT id') < istat = nf_get_var_double(ncid,idv_ut,amie_sh_ut) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE variable ut') < ! < ! Get HPI < istat = nf_inq_varid(ncid,'hpi',idv_hpi) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE hpi id') < istat = nf_get_var_double(ncid,idv_hpi,amie_hpi_sh) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE variable hpi') < ! < ! Get PCP < istat = nf_inq_varid(ncid,'pcp',idv_pcp) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE pcp id') < istat = nf_get_var_double(ncid,idv_pcp,amie_pcp_sh) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE variable hpi') < ! < ! Get cusplat < istat = nf_inq_varid(ncid,'cusplat',idv_cusplat) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE cusplat id') < istat = nf_get_var_double(ncid,idv_cusplat,amie_cusplat_sh) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE variable cusplat') < ! < ! Get cuspmlt < istat = nf_inq_varid(ncid,'cuspmlt',idv_cuspmlt) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE cusplat id') < istat = nf_get_var_double(ncid,idv_cuspmlt,amie_cuspmlt_sh) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE variable cuspmlt') < ! < ! Allocate 2-d fields: < ! real,allocatable,dimension(:,:,:) :: ! (lonp1,latp1) < ! | pot_sh_amie, ekv_sh_amie, efx_sh_amie < ! < if (.not. allocated(pot_sh_amie)) < | allocate(pot_sh_amie(lonp1,latp1),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_sh: error allocating', < | ' pot_sh_amie: lonp1=',i3,' latp1=',i3)")lonp1,latp1 < if (.not. allocated(ekv_sh_amie)) < | allocate(ekv_sh_amie(lonp1,latp1),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_sh: error allocating', < | ' ekv_sh_amie: lonp1=',i3,' latp1=',i3)")lonp1,latp1 < if (.not. allocated(efx_sh_amie)) < | allocate(efx_sh_amie(lonp1,latp1),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_sh: error allocating', < | ' efx_sh_amie: lonp1=',i3,' latp1=',i3)")lonp1,latp1 < < ! < ! Allocate 3-d fields: < ! real,allocatable,dimension(:,:,:) :: ! (lonp1,latp1,ntimes) < ! | amie_pot_nh, amie_pot_sh, amie_ekv_nh, amie_ekv_sh, < ! | amie_efx_nh, amie_efx_sh < ! < if (.not. allocated(amie_pot_sh)) < | allocate(amie_pot_sh(lonp1,latp1,ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_sh: error allocating', < | ' amie_pot_nh: lonp1=',i3,' latp1=',i3,' ntimes=',i3)") < | lonp1,latp1,ntimes < if (.not. allocated(amie_ekv_sh)) < | allocate(amie_ekv_sh(lonp1,latp1,ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_sh: error allocating', < | ' amie_ekv_nh: lonp1=',i3,' latp1=',i3,' ntimes=',i3)") < | lonp1,latp1,ntimes < if (.not. allocated(amie_efx_sh)) < | allocate(amie_efx_sh(lonp1,latp1,ntimes),stat=ier) < if (ier /= 0) write(6,"('>>> rdamie_sh: error allocating', < | ' amie_efx_nh: lonp1=',i3,' latp1=',i3,' ntimes=',i3)") < | lonp1,latp1,ntimes < ! < ! Get 3-D AMIE fields (lon,lat,ntimes) < ! < ! AMIE electric potential < istat = nf_inq_varid(ncid,'pot',idv_pot) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE electric potential id') < istat = nf_get_var_double(ncid,idv_pot,amie_pot_sh) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE variable pot') < ! < ! AMIE mean energy < istat = nf_inq_varid(ncid,'ekv',idv_ekv) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE mean energy id') < istat = nf_get_var_double(ncid,idv_ekv,amie_ekv_sh) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE variable ekv') < ! < ! AMIE energy flux < istat = nf_inq_varid(ncid,'efx',idv_efx) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE energy flux id') < istat = nf_get_var_double(ncid,idv_efx,amie_efx_sh) < if (istat /= NF_NOERR) call rpt_ncerr(istat, < | 'rdamie_sh: Error getting SH AMIE variable efx') < ! < ! Close the file: < istat = nf_close(ncid) < write(6,"('Completed read from SH AMIE data file ',a)") < | trim(dskfile) < write(6,"(72('-'),/)") < end subroutine rdamie_sh --- > 491a189 > use ModKind ! 6/12: for AMIE U MI 496a195,199 > ! 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) 512c215,228 < real :: del,xmlt,dmlat,dlatm,dlonm,dmltm,rot,dtr --- > ! 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 > 522,577c238,245 < ! Reading Charley Barth's auroral modification factor < integer, parameter:: nxtimes=47 < integer :: nxyear,nxday(nxtimes) < real :: fxfact_nh, fxfact_sh,fxtime,fxtime1,fxtime2 < real, dimension(nxtimes):: ut_nh,ut_sh,fx_nh,fx_sh < ! < ! External: < real,external :: finterp < ! < nxyear = 1998 < nxday(1:nxtimes)=(/267, < | 267, 267, 267, 267, 267, 267, 267, < | 267, 267, 267, 267, 267, 267, 267, < | 267, 268, 268, 268, 268, 268, 268, < | 268, 268, 268, 268, 268, 268, 268, < | 268, 268, 269, 269, 269, 269, 269, < | 269, 269, 269, 269, 269, 269, 269, < | 269, 269, 269, 270/) < ut_nh(1:nxtimes)=(/0.0, < | 0.68872, 2.28569, 3.88347, 5.47933, 7.07654, 8.67133, < | 10.2672, 11.8640, 13.4619, 15.0575, 16.6542, 18.2520, < | 19.8473, 21.4437, 23.0409, 0.63913, 2.23508, 3.83221, < | 5.42711, 7.02325, 8.62048, 10.2171, 11.8147, 13.4118, < | 15.0064, 16.6022, 18.1989, 19.7967, 21.3921, 22.9884, < | 0.58565, 2.18393, 3.78000, 5.37563, 6.97244, 8.56865, < | 10.1659, 11.7654, 13.3620, 14.9582, 16.5557, 18.1509, < | 19.7470, 21.3442, 22.9422, 0./) < fx_nh(1:nxtimes)=(/1.0, < | 0.954679, 1.01677, 0.805276, 0.662938, 0.680216, 1.09095, < | 0.883288,0.972496, 0.720709, 0.579112, 0.818075, 1.28912, < | 1.35952, 1.36740, 1.03170, 1.34414, 1.53157, 0.90830, < | 0.824193, 1.10184, 1.35307, 0.999406, 0.788716, 0.76273, < | 0.647832, 0.57241, 0.52249, 1.48004, 2.36089, 1.86270, < | 1.87322, 1.80098, 1.57992, 1.33901, 1.21373, 1.42878, < | 1.28782, 1.2731, 1.31435, 1.13873, 0.984014, 1.04712, < | 1.48459, 1.75276, 1.94827, 1./) < < ut_sh(1:nxtimes)=(/0.0, < | 0.104409, 1.70090, 3.29847, 4.89378, 6.49047, 8.08820, < | 9.68356, 11.2801, 12.8777, 14.4745, 16.0726, 17.6666, < | 19.2649, 20.8609, 22.4578, 0.05568, 1.65123, 3.24788, < | 4.84236, 6.44144, 8.03825, 9.63280, 11.2301, 12.8267, < | 14.4227, 16.0197, 17.6161, 19.2119, 20.8103, 22.4062, < | 0.02522, 1.59759, 3.19667, 4.79183, 6.38823, 7.98401, < | 9.58092, 11.1789, 12.7793, 14.3752, 15.9723, 17.5671, < | 19.1629, 20.7597, 22.3572, 0./) < < fx_sh(1:nxtimes)=(/1.0, < | 0.544252, 0.487370, 0.572375, 0.605849, 0.605344, 0.776822, < | 0.895618, 0.641604, 0.590224, 0.638560, 0.739044, 0.863948, < | 1.22977, 1.16256, 0.881014, 0.669191, 0.617876, 0.498846, < | 0.998803, 0.826644, 0.899562, 1.41258, 1.05342, 0.972168, < | 1.00838, 1.50718, 2.00560, 2.01770, 1.85613, 1.31682, < | 0.920968, 1.64630, 2.33071, 1.96759, 1.83485, 1.79565, < | 1.74479, 1.52049, 1.76034, 1.85011, 1.32154, 1.19636, < | 1.28152, 1.42212, 1.49638, 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. 585a254,255 > > maxefx=50. 591,604c261,262 < ! Check times: < ! < if (year(1) /= iyear) then < write(6,"('>>> getamie: wrong year,iyear=',2i4)") < | year(1),iyear < stop 'getamie' < endif < amieday = jday(1) < if (amie_ibkg > 0) amieday = jday(1) - 1 < ! if (amieday /= iday .and. amie_ibkg == 0) then < ! write(6,"('>>> getamie: wrong amieday,modelday=',2i4)") < ! | amieday,iday < ! stop 'getamie' < ! endif --- > > ! Check times 616,727c274,286 < maxefx=50. < if (amieday==314 .and. abs(model_ut-7.333) < 0.1) < | maxefx=21. < if (amieday==314 .and. abs(model_ut-8.833) < 0.1) < | maxefx=14. < if (amieday==318) maxefx=10. < if (amieday==318 .and. model_ut < 5.75) maxefx=3. < if (amieday==318 .and. abs(model_ut-4.75) < 0.1) maxefx=1. < if (amieday==318 .and. abs(model_ut-21.0) < 0.1) maxefx=2. < < ! Find out the multiplication factor from Charley Barth < fxfact_nh = 1. < fxfact_sh = 1. < if (nxyear == iyear) then < ! For the southern hemisphere < iset = nxtimes < iset1 = nxtimes < do i=1,nxtimes < fxtime=(nxday(i)-iday)*24.+ut_nh(i) < if (fxtime < model_ut) iset=i < enddo < iset1 = iset + 1 < if (iset==nxtimes) iset1=iset < fxtime1=(iday-nxday(iset))*24.+ut_nh(iset) < fxtime2=(iday-nxday(iset1))*24.+ut_nh(iset1) < denoma = fxtime2 - fxtime1 < if (denoma == 0.) then < f1 = 1. < f2 = 0. < else < ! f1 = (fxtime2-model_ut)/denoma < ! f2 = (model_ut-fxtime1)/denoma < f1 = (fxtime2-(model_ut+(modeltime(1)-jday(i))*24.))/denoma < f2 = (model_ut+(modeltime(1)-jday(i))*24.-fxtime1)/denoma < endif < fxfact_nh = f1*fx_nh(iset)+f2*fx_nh(iset1) < ! For the southern hemisphere < iset = nxtimes < iset1 = nxtimes < do i=1,nxtimes < fxtime=(nxday(i)-iday)*24.+ut_sh(i) < ! if (fxtime < model_ut) iset=i < if (fxtime < model_ut+(modeltime(1)-jday(i))*24.) iset=i < enddo < iset1 = iset + 1 < if (iset==nxtimes) iset1=iset < fxtime1=(iday-nxday(iset))*24.+ut_sh(iset) < fxtime2=(iday-nxday(iset1))*24.+ut_sh(iset1) < denoma = fxtime2 - fxtime1 < if (denoma == 0.) then < f1 = 1. < f2 = 0. < else < ! f1 = (fxtime2-model_ut)/denoma < ! f2 = (model_ut-fxtime1)/denoma < f1 = (fxtime2-(model_ut+(modeltime(1)-jday(i))*24.))/denoma < f2 = (model_ut+(modeltime(1)-jday(i))*24.-fxtime1)/denoma < endif < fxfact_sh = f1*fx_sh(iset)+f2*fx_sh(iset1) < endif < ! < ! Interpolate AMIE 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 < pot_sh_amie(:,:) = 0. < ekv_sh_amie(:,:) = 0. < efx_sh_amie(:,:) = 0. < cusplat_sh_amie = 0. < cuspmlt_sh_amie = 0. < hpi_sh_amie = 0. < pcp_sh_amie = 0. < ! < n = len_trim(amiesh) < iboxcar = 0 < if (n > 0) then < nn = size(amie_sh_ut) < if (amie_ibkg == 0) then < iset = nn < iset1 = nn < do i=1,nn < ! if (amie_sh_ut(i) < model_ut) iset = i < if (amie_sh_ut(i) < model_ut+(modeltime(1)-jday(i))*24.) < | iset = i < enddo < write(6,"('getamie: AMIE SH Data nn,iset,day1,day2=',4i5)") < | nn,iset,jday(1),jday(nn) < iset1 = iset + 1 < if (iset == nn) iset1 = iset < ! get rid of outragous flux values < inpt = 0 < do j=1,latp1 < do i=1,lonp1 < if (amie_efx_sh(i,j,iset) > maxefx) then < if (i > 1) then < amie_efx_sh(i,j,iset)=amie_efx_sh(i-1,j,iset) < else < amie_efx_sh(i,j,iset)=amie_efx_sh(i+1,j,iset) < endif < if (j < 12) amie_efx_sh(i,j,iset) = 1. < inpt = inpt + 1 < endif < if (iset1 /= iset .and. amie_efx_sh(i,j,iset1) > maxefx) then < if (i > 1) then < amie_efx_sh(i,j,iset1)=amie_efx_sh(i-1,j,iset1) < else < amie_efx_sh(i,j,iset1)=amie_efx_sh(i+1,j,iset1) < endif < if (j < 12) amie_efx_sh(i,j,iset1) = 1. < inpt = inpt + 1 < endif < enddo --- > ! 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 729,730c288,313 < if (inpt > 0) write(6,"('maxefx = ',f4.1, < | ' Bad efx points SH= ',i3)")maxefx,inpt --- > 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 732,805c315,322 < denoma = amie_sh_ut(iset1) - amie_sh_ut(iset) < if (denoma == 0.) then < f1 = 1. < f2 = 0. < else < ! f1 = (amie_sh_ut(iset1) - model_ut)/denoma < ! f2 = (model_ut - amie_sh_ut(iset))/denoma < f1 = (amie_sh_ut(iset1) - (model_ut+(modeltime(1)- < | jday(iset1))*24.))/denoma < f2 = (model_ut+(modeltime(1)-jday(iset1))*24. - < | amie_sh_ut(iset))/denoma < endif < ! write(6,"('getamie: AMIE SH Data n,iset,modeltime,f1,f2 =', < ! | 4i5,2f5.2)")n,iset,modeltime(1),jday(iset1),f1,f2 < write(6,"('getamie: AMIE SH Data model_day,model_ut,amie_day,', < | 'amie_ut,f1,f2 =',i4,f7.1,i4,f7.1,2f5.2)")modeltime(1), < | model_ut,jday(iset),amie_sh_ut(iset1),f1,f2 < cusplat_sh_amie = (f1*amie_cusplat_sh(iset1) + < | f2*amie_cusplat_sh(iset)) < cuspmlt_sh_amie = (f1*amie_cuspmlt_sh(iset1) + < | f2*amie_cuspmlt_sh(iset)) < hpi_sh_amie = (f1*amie_hpi_sh(iset1) + f2*amie_hpi_sh(iset)) < pcp_sh_amie = (f1*amie_pcp_sh(iset1) + f2*amie_pcp_sh(iset)) < if (iboxcar == 0) then < pot_sh_amie(:,:) = (f1*amie_pot_sh(:,:,iset1) + < | f2*amie_pot_sh(:,:,iset)) < ekv_sh_amie(:,:) = (f1*amie_ekv_sh(:,:,iset1) + < | f2*amie_ekv_sh(:,:,iset)) < efx_sh_amie(:,:) = (f1*amie_efx_sh(:,:,iset1) + < | f2*amie_efx_sh(:,:,iset)) < else < call boxcar_ave(amie_pot_sh,pot_sh_amie,lonp1,latp1, < | nn,iset,iboxcar) < call boxcar_ave(amie_efx_sh,efx_sh_amie,lonp1,latp1, < | nn,iset,iboxcar) < call boxcar_ave(amie_ekv_sh,ekv_sh_amie,lonp1,latp1, < | nn,iset,iboxcar) < endif < else < if (amie_ibkg == 1) then < pot_sh_amie(:,:) = amie_pot_sh(:,:,1) < ekv_sh_amie(:,:) = amie_ekv_sh(:,:,1) < efx_sh_amie(:,:) = amie_efx_sh(:,:,1) < cusplat_sh_amie = amie_cusplat_sh(1) < cuspmlt_sh_amie = amie_cuspmlt_sh(1) < hpi_sh_amie = amie_hpi_sh(1) < pcp_sh_amie = amie_pcp_sh(1) < else if (amie_ibkg == 3) then < pot_sh_amie(:,:) = amie_pot_sh(:,:,241) < ekv_sh_amie(:,:) = amie_ekv_sh(:,:,241) < efx_sh_amie(:,:) = amie_efx_sh(:,:,241) < cusplat_sh_amie = amie_cusplat_sh(241) < cuspmlt_sh_amie = amie_cuspmlt_sh(241) < hpi_sh_amie = amie_hpi_sh(241) < pcp_sh_amie = amie_pcp_sh(241) < else < do i=1,nn < pot_sh_amie(:,:) = pot_sh_amie(:,:) + amie_pot_sh(:,:,1) < ekv_sh_amie(:,:) = ekv_sh_amie(:,:) + amie_ekv_sh(:,:,1) < efx_sh_amie(:,:) = efx_sh_amie(:,:) + amie_efx_sh(:,:,1) < cusplat_sh_amie = cusplat_sh_amie + amie_cusplat_sh(1) < cuspmlt_sh_amie = cuspmlt_sh_amie + amie_cuspmlt_sh(1) < hpi_sh_amie = hpi_sh_amie + amie_hpi_sh(1) < pcp_sh_amie = pcp_sh_amie + amie_pcp_sh(1) < enddo < pot_sh_amie(:,:) = pot_sh_amie(:,:)/nn < ekv_sh_amie(:,:) = ekv_sh_amie(:,:)/nn < efx_sh_amie(:,:) = efx_sh_amie(:,:)/nn < cusplat_sh_amie = cusplat_sh_amie/nn < cuspmlt_sh_amie = cuspmlt_sh_amie/nn < hpi_sh_amie = hpi_sh_amie/nn < pcp_sh_amie = pcp_sh_amie/nn < endif < 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 807,820c324,380 < ! get rid of outragous flux values < inpt = 0 < do j=1,latp1 < do i=1,lonp1 < if (efx_sh_amie(i,j) > maxefx) then < ! efx_sh_amie(i,j) = 1. < if (i > 1) then < efx_sh_amie(i,j)=efx_sh_amie(i-1,j) < else < efx_sh_amie(i,j)=efx_sh_amie(i+1,j) < endif < inpt = inpt + 1 < endif < enddo --- > > 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 hpi = hpi + da(ilat)*efx_sh_amie(imlt,ilat) > endif > else > potN(imlt,ilat) = TempPotential(imlt,ilat+jl) > pot_nh_amie(imlt,ilat) = TempPotential(imlt,ilat+jl) > ekv_nh_amie(imlt,ilat) = AveEOut(imlt,ilat+jl) > efx_nh_amie(imlt,ilat) = EfluxOut(imlt,ilat+jl) > if (imlt hpi = hpi + da(ilat)*efx_nh_amie(imlt,ilat) > endif > endif > if (TempPotential(imlt,ilat+jl) > 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 822,869d381 < if (inpt > 0) write(6,"('Total bad efx points SH= ',i3)")inpt < ! < ! get NH AMIE data < pot_nh_amie(:,:) = 0. < ekv_nh_amie(:,:) = 0. < efx_nh_amie(:,:) = 0. < cusplat_nh_amie = 0. < cuspmlt_nh_amie = 0. < hpi_nh_amie = 0. < pcp_nh_amie = 0. < < n = len_trim(amienh) < iboxcar = 0 < if (n > 0) then < nn = size(amie_nh_ut) < ! write(6,"('getamie: Interpolate AMIE NH Data nn=',i3)")nn < if (amie_ibkg == 0) then < iset = nn < iset1 = nn < do i=1,nn < if (amie_nh_ut(i) < model_ut+(modeltime(1)-jday(i))*24.) < | iset = i < enddo < iset1 = iset + 1 < if (iset == nn) iset1 = iset < ! get rid of outragous flux values < inpt = 0 < do j=1,latp1 < do i=1,lonp1 < if (amie_efx_nh(i,j,iset) > maxefx) then < if (i > 1) then < amie_efx_nh(i,j,iset)=amie_efx_nh(i-1,j,iset) < else < amie_efx_nh(i,j,iset)=amie_efx_nh(i+1,j,iset) < endif < if (j < 12) amie_efx_nh(i,j,iset) = 1. < inpt = inpt + 1 < endif < if (iset1 /= iset .and. amie_efx_nh(i,j,iset1) > maxefx) then < if (i > 1) then < amie_efx_nh(i,j,iset1)=amie_efx_nh(i-1,j,iset1) < else < amie_efx_nh(i,j,iset1)=amie_efx_nh(i+1,j,iset1) < endif < if (j < 12) amie_efx_nh(i,j,iset1) = 1. < inpt = inpt + 1 < endif < enddo 871,872c383,419 < if (inpt > 0) write(6,"('maxefx = ',f4.1, < | ' Bad efx points NH= ',i3)")maxefx,inpt --- > 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 874,994d420 < denoma = amie_nh_ut(iset1) - amie_nh_ut(iset) < if (denoma == 0.) then < f1 = 1. < f2 = 0. < else < ! f1 = (amie_nh_ut(iset1) - model_ut)/denoma < ! f2 = (model_ut - amie_nh_ut(iset))/denoma < f1 = (amie_nh_ut(iset1) - (model_ut+(modeltime(1)- < | jday(iset1))*24.))/denoma < f2 = (model_ut+(modeltime(1)-jday(iset1))*24. - < | amie_nh_ut(iset))/denoma < endif < cusplat_nh_amie = (f1*amie_cusplat_nh(iset1) + < | f2*amie_cusplat_nh(iset)) < cuspmlt_nh_amie = (f1*amie_cuspmlt_nh(iset1) + < | f2*amie_cuspmlt_nh(iset)) < hpi_nh_amie = (f1*amie_hpi_nh(iset1) + f2*amie_hpi_nh(iset)) < pcp_nh_amie = (f1*amie_pcp_nh(iset1) + f2*amie_pcp_nh(iset)) < if (iboxcar == 0) then < pot_nh_amie(:,:) = (f1*amie_pot_nh(:,:,iset1) + < | f2*amie_pot_nh(:,:,iset)) < ekv_nh_amie(:,:) = (f1*amie_ekv_nh(:,:,iset1) + < | f2*amie_ekv_nh(:,:,iset)) < efx_nh_amie(:,:) = (f1*amie_efx_nh(:,:,iset1) + < | f2*amie_efx_nh(:,:,iset)) < ! write(6,"('ekv_nh_amie min, max = ',2e12.4)") < ! | minval(ekv_nh_amie),maxval(ekv_nh_amie) < else < call boxcar_ave(amie_pot_nh,pot_nh_amie,lonp1,latp1, < | nn,iset,iboxcar) < ! call fminmax(amie_pot_nh(:,:,iset),lonp1*latp1,fmin,fmax) < ! write(6,"('AMIE pot max,min = ',2f8.0)")fmax,fmin < ! call fminmax(pot_nh_amie(:,:),lonp1*latp1,fmin,fmax) < ! write(6,"('boxcar_ave AMIE pot max,min= ',2f8.0)")fmax,fmin < call boxcar_ave(amie_efx_nh,efx_nh_amie,lonp1,latp1, < | nn,iset,iboxcar) < ! call fminmax(amie_efx_nh(:,:,iset),lonp1*latp1,fmin,fmax) < ! write(6,"('AMIE efx max,min = ',2f8.0)")fmax,fmin < ! call fminmax(efx_nh_amie(:,:),lonp1*latp1,fmin,fmax) < ! write(6,"('boxcar_ave AMIE efx max,min= ',2f8.0)")fmax,fmin < call boxcar_ave(amie_ekv_nh,ekv_nh_amie,lonp1,latp1, < | nn,iset,iboxcar) < ! call fminmax(amie_ekv_nh(:,:,iset),lonp1*latp1,fmin,fmax) < ! write(6,"('AMIE ekv max,min = ',2f8.0)")fmax,fmin < ! call fminmax(ekv_nh_amie(:,:),lonp1*latp1,fmin,fmax) < ! write(6,"('boxcar_ave AMIE ekv max,min= ',2f8.0)")fmax,fmin < endif < else < if (amie_ibkg == 1) then < pot_nh_amie(:,:) = amie_pot_nh(:,:,1) < ekv_nh_amie(:,:) = amie_ekv_nh(:,:,1) < efx_nh_amie(:,:) = amie_efx_nh(:,:,1) < cusplat_nh_amie = amie_cusplat_nh(1) < cuspmlt_nh_amie = amie_cuspmlt_nh(1) < hpi_nh_amie = amie_hpi_nh(1) < pcp_nh_amie = amie_pcp_nh(1) < else if (amie_ibkg == 3) then < pot_nh_amie(:,:) = amie_pot_nh(:,:,241) < ekv_nh_amie(:,:) = amie_ekv_nh(:,:,241) < efx_nh_amie(:,:) = amie_efx_nh(:,:,241) < cusplat_nh_amie = amie_cusplat_nh(241) < cuspmlt_nh_amie = amie_cuspmlt_nh(241) < hpi_nh_amie = amie_hpi_nh(241) < pcp_nh_amie = amie_pcp_nh(241) < else < do i=1,nn < pot_nh_amie(:,:) = pot_nh_amie(:,:) + amie_pot_nh(:,:,1) < ekv_nh_amie(:,:) = ekv_nh_amie(:,:) + amie_ekv_nh(:,:,1) < efx_nh_amie(:,:) = efx_nh_amie(:,:) + amie_efx_nh(:,:,1) < cusplat_nh_amie = cusplat_nh_amie + amie_cusplat_nh(1) < cuspmlt_nh_amie = cuspmlt_nh_amie + amie_cuspmlt_nh(1) < hpi_nh_amie = hpi_nh_amie + amie_hpi_nh(1) < pcp_nh_amie = pcp_nh_amie + amie_pcp_nh(1) < enddo < pot_nh_amie(:,:) = pot_nh_amie(:,:)/nn < ekv_nh_amie(:,:) = ekv_nh_amie(:,:)/nn < efx_nh_amie(:,:) = efx_nh_amie(:,:)/nn < cusplat_nh_amie = cusplat_nh_amie/nn < cuspmlt_nh_amie = cuspmlt_nh_amie/nn < hpi_nh_amie = hpi_nh_amie/nn < pcp_nh_amie = pcp_nh_amie/nn < endif < endif < endif < ! get rid of outragous flux values < inpt = 0 < do j=1,latp1 < do i=1,lonp1 < if (efx_nh_amie(i,j) > maxefx) then < ! efx_nh_amie(i,j) = 1. < if (i > 1) then < efx_nh_amie(i,j)=efx_nh_amie(i-1,j) < else < efx_nh_amie(i,j)=efx_nh_amie(i+1,j) < endif < inpt = inpt + 1 < endif < ! Reduce auroral flux intensity for HPI < 80 GW for Sep. 2005 event < ! if (hpi_nh_amie < 70.) efx_nh_amie(i,j) = efx_nh_amie(i,j)/2. < enddo < enddo < if (inpt > 0) write(6,"('Total bad efx points NH = ',i3)")inpt < ! Multiply Barth's factor for auroral fluxes < if (nxyear==iyear) then < do j=1,19 < do i=1,lonp1 < if (ekv_nh_amie(i,j) > 2.5) < | ekv_nh_amie(i,j) = 1.5*ekv_nh_amie(i,j) < if (ekv_sh_amie(i,j) > 2.5) < | ekv_sh_amie(i,j) = 1.5*ekv_sh_amie(i,j) < enddo < enddo < endif < fxfact_nh = 1. < fxfact_sh = 1. < if (fxfact_nh > 1. .or. fxfact_sh > 1.) then < write(6,"('Barth Multiplication factors = ',2i5,2f7.3)") < | iyear, nxyear, fxfact_nh, fxfact_sh < efx_nh_amie(:,:) = efx_nh_amie(:,:)*fxfact_nh < efx_sh_amie(:,:) = efx_sh_amie(:,:)*fxfact_sh < endif 1000,1014c426 < ! Define this latitude to be between 70 and 77.5 degrees < ! < ! if (cusplat_sh_amie > 65.0) then < ! cusplat_sh_amie = 65.0 < ! cuspmlt_sh_amie = 11. < ! endif < if (cusplat_sh_amie > 75.0) then < cusplat_sh_amie = 75.0 < cuspmlt_sh_amie = 11. < endif < if (cusplat_sh_amie < 60.0) then < cusplat_sh_amie = 60.0 < cuspmlt_sh_amie = 11. < endif < ! cusplat_sh_amie = amin1(65.0,cusplat_sh_amie) --- > ! Define this latitude to be between 70 and 77.5 degrees (then 60-75 mlat!) 1016,1032d427 < ! if (cusplat_nh_amie > 65.0) then < ! cusplat_nh_amie = 65.0 < ! cuspmlt_nh_amie = 11. < ! endif < if (cusplat_nh_amie > 75.0) then < cusplat_nh_amie = 75.0 < cuspmlt_nh_amie = 11. < endif < if (cusplat_nh_amie < 60.0) then < cusplat_nh_amie = 60.0 < cuspmlt_nh_amie = 11. < endif < ! cusplat_nh_amie = amin1(65.0,cusplat_nh_amie) < if (cuspmlt_sh_amie > 12.5) cuspmlt_sh_amie = 12.5 < if (cuspmlt_sh_amie < 11.0) cuspmlt_sh_amie = 11.0 < if (cuspmlt_nh_amie > 12.5) cuspmlt_nh_amie = 12.5 < if (cuspmlt_nh_amie < 11.0) cuspmlt_nh_amie = 11.0 1036a432,464 > > ! 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 > 1046,1047c474,475 < ! AMIE grid goes from -90 to ~-40 for SH and +90 to ~+40 for NH. < ! TGCM grid goes from -87.5 to -2.5 for SH and +2.5 to +87.5 for NH. --- > ! 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) 1110,1115c538,546 < C DMLAT=lat spacing in degrees of AMIE apex grid < dtr = pi/180. < DMLAT = 180. / FLOAT(jmxm-1) < DLATM = DMLAT*dtr < DLONM = 2.*PI/FLOAT(lonmx) < DMLTM = 24./FLOAT(lonmx) --- > 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) 1269,1270c700,701 < subroutine mag2geo (am,ag,im,jm,dim,djm,lg,lm,nlong,nlatg, < | nlonm,nlatm) --- > subroutine mag2geo(am,ag,im,jm,dim,djm,lg,lm,nlong,nlatg,nlonm, > | nlatm) 1302,1317c733 < !------------------------------------------------------------------- < 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 --- >