module meped ! ! Import electron data from MEPED (Medium Energy Proton and Electron ! Detector) instrument on NOAA POES (TIROS) satellite. ! implicit none real,allocatable :: elec_pw(:), elec_char(:) integer,allocatable :: meped_year(:),meped_day(:),meped_sec(:) real,allocatable :: time_meped(:) real :: missing_data integer :: ntime contains !----------------------------------------------------------------------- subroutine rdmeped(meped_file) use nchist_module,only: handle_ncerr ! ! Read MEPED data file (this is called once per run from init, if ! namelist read meped_file is non-blank). ! implicit none #include ! ! Args: character(len=*),intent(in) :: meped_file ! ! Local: integer :: istat,ncid,iatt character(len=120) :: char120,name integer :: id_time,idv_time,id_year,id_day,id_sec integer :: id_elecpw,id_elecchar,id_missing ! ! Open file: write(6,"(/,72('-'))") istat = nf_open(meped_file,NF_NOWRITE,ncid) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_open for netcdf ', | 'file ',a,' status=',a)") trim(meped_file),istat call handle_ncerr(istat,char120) endif ! ! From ncdump (only electron vars): ! ! netcdf meped { ! dimensions: ! time = 6992 ; ! variables: ! double time(time) ; ! int year(time) ; ! int day(time) ; ! int sec(time) ; ! double elec_pw(time) ; ! elec_pw:long_name = "Hemispheric Power for electrons" ; ! elec_pw:units = "GW" ; ! double el_char(time) ; ! el_char:long_name = "Characteristic energy for electons" ; ! el_char:units = "KeV" ; ! ! Get time dimension: istat = nf_inq_dimid(ncid,'time',id_time) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting time dimension id') istat = nf_inq_dim(ncid,id_time,name,ntime) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting time dimension length') write(6,"('rdmeped: id_time=',i3,' ntime=',i5)") id_time,ntime ! ! Allocate data arrays: allocate(elec_pw(ntime),stat=istat) if (istat /= 0) | write(6,"('>>> rdmeped: error allocating elec_pw: ntime=',i5)") | ntime allocate(elec_char(ntime),stat=istat) if (istat /= 0) | write(6,"('>>> rdmeped: error allocating elec_char: ntime=', | i5)") ntime allocate(time_meped(ntime),stat=istat) if (istat /= 0) | write(6,"('>>> rdmeped: error allocating time_meped: ntime=', | i5)") ntime allocate(meped_year(ntime),stat=istat) if (istat /= 0) | write(6,"('>>> rdmeped: error allocating meped_year: ntime=', | i5)") ntime allocate(meped_day(ntime),stat=istat) if (istat /= 0) | write(6,"('>>> rdmeped: error allocating meped_day: ntime=', | i5)") ntime allocate(meped_sec(ntime),stat=istat) if (istat /= 0) | write(6,"('>>> rdmeped: error allocating meped_sec: ntime=', | i5)") ntime ! ! Get variable id's and read data (electron and time data only): ! ! Electron power (GW): istat = nf_inq_varid(ncid,'elec_pw',id_elecpw) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting elec_pw var id') istat = nf_get_var_double(ncid,id_elecpw,elec_pw) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting elec_pw var') write(6,"('rdmeped: elec_pw min,max=',2e12.4)") | minval(elec_pw),maxval(elec_pw) ! ! Electron characteristic energy (KeV): istat = nf_inq_varid(ncid,'el_char',id_elecchar) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting el_char var id') istat = nf_get_var_double(ncid,id_elecchar,elec_char) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting el_char var') write(6,"('rdmeped: elec_char min,max=',2e12.4)") | minval(elec_char),maxval(elec_char) ! ! MEPED time (real total secs): istat = nf_inq_varid(ncid,'time',idv_time) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting time var id') istat = nf_get_var_double(ncid,idv_time,time_meped) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting time var') ! ! MEPED year: istat = nf_inq_varid(ncid,'year',id_year) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting year var id') istat = nf_get_var_int(ncid,id_year,meped_year) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting year var') ! ! MEPED day: istat = nf_inq_varid(ncid,'day',id_day) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting day var id') istat = nf_get_var_int(ncid,id_day,meped_day) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting day var') ! ! MEPED seconds: istat = nf_inq_varid(ncid,'sec',id_sec) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting sec var id') istat = nf_get_var_int(ncid,id_sec,meped_sec) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting sec var') ! ! Missing data value: istat = nf_inq_attid(ncid,id_missing,'missing_value',iatt) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdmeped: Error getting miising_data att id') istat = nf_get_att_double(ncid,NF_GLOBAL,id_missing, | name,missing_data) write(6,"('rdmeped: missing_data=',e12.4)") missing_data ! ! write(6,"('rdmeped: elec_char=',/,(8e10.2))") elec_char ! write(6,"('rdmeped: elec_pw =',/,(8e10.2))") elec_pw ! write(6,"('rdmeped: meped_year=',/,(12i5))") meped_year ! write(6,"('rdmeped: meped_day =',/,(12i5))") meped_day ! write(6,"('rdmeped: meped_sec =',/,(10i7))") meped_sec ! write(6,"(72('-'),/)") end subroutine rdmeped !----------------------------------------------------------------------- subroutine get_meped(model_year,model_day,rmodel_secs,epower, | echar) ! ! Interpolate MEPED data to current model time. ! This is called at every timestep from aurora.F, if namelist read ! meped_file is non-blank (data was read at run initialization ! by sub rdmeped). ! implicit none ! ! Args: integer,intent(in) :: model_year,model_day real,intent(in) :: rmodel_secs real,intent(out) :: epower,echar ! ! Local: integer :: i,i0,i1,ndpy,nspy real :: total_model_secs integer(kind=8) :: meped_total_sec0,meped_total_sec1, | model_total_secs ! ! External: real,external :: finterp_bigints ! ! Total model seconds (since Jan 1, year 0000) ndpy = 365 if (mod(model_year,4)==0) ndpy = 366 ! leap year nspy = 3600*24*ndpy ! number of seconds per year total_model_secs = float(nspy)*float(model_year)+ | 3600.*24.*float(model_day)+rmodel_secs ! ! Bracket model time: i0 = 0 do i=2,ntime if (total_model_secs <= time_meped(i) .and. | total_model_secs >= time_meped(i-1)) then i0 = i-1 i1 = i endif enddo ! i=1,ntime if (i0==0) then write(6,"(/,'>>> get_meped: could not bracket model ', | 'time: year=',i4,' day=',i3,' secs=',f10.2)") | model_year,model_day,rmodel_secs write(6,"('MEPED data file contains data for ', | 'year,day ',i4,', ',i3,' to ',i4,', ',i3)") | meped_year(1),meped_day(1),meped_year(ntime), | meped_day(ntime) call shutdown('get_meped') else ! write(6,"('get_meped: bracket model time: i0=',i5,' i1=',i5, ! | ' time_meped(i0)=',f15.2,' time_meped(i1)=',f15.2)") ! | i0,i1,time_meped(i0),time_meped(i1) endif ! ! Linear interpolation (8-byte ints for total seconds): ! meped_total_sec0 = int(time_meped(i0),8) meped_total_sec1 = int(time_meped(i1),8) model_total_secs = int(total_model_secs,8) epower = finterp_bigints(elec_pw(i0),elec_pw(i1),meped_total_sec0, | meped_total_sec1,model_total_secs) echar = finterp_bigints(elec_char(i0),elec_char(i1), | meped_total_sec0,meped_total_sec1,model_total_secs) ! write(6,"('get_meped: model year,day,secs=',i4,', ',i3,', ',i6, ! | ' model total secs=',i12)") model_year,model_day, ! | int(rmodel_secs),model_total_secs ! write(6,"('get_meped: meped year0,day0,sec0=',i4,', ',i3,', ', ! | i6,' meped total secs0=',i15)") meped_year(i0),meped_day(i0), ! | meped_sec(i0),meped_total_sec0 ! write(6,"('get_meped: meped year1,day1,sec1=',i4,', ',i3,', ', ! | i6,' meped total secs1=',i15)") meped_year(i1),meped_day(i1), ! | meped_sec(i1),meped_total_sec1 ! write(6,"('elec_pw0,1 =',2e12.4,' epower=',e12.4)") elec_pw(i0), ! | elec_pw(i1),epower ! write(6,"('elec_char0,1=',2e12.4,' echar =',e12.4)") ! | elec_char(i0),elec_char(i1),echar end subroutine get_meped !----------------------------------------------------------------------- end module meped