program main ! ! Astrid Maute am 6/04 ! for SUN,LINUX & IBM maschines ! calculate the magnetic perturbation at the ground & equivalent current ! process different station which are given in file ./stat.inp ! to change time point which is read in: change start_time ! specify start_time,end_time in iface.F , ! fileapex, *.nc input-file and path in param.F ! reading data from TIEGCM secondary hist output (path + flnm) ! calculate the magnetic perturbation at the ground ! please read the Readme for more information ! use param_module,only: start_time,end_time,abr integer :: ns,num_stat,ios ! start_time = 2 ! 27 ! (user given) index of start time end_time = 9 ! 27 ! (user given) index of end time ! ! get station abreviation open(30,file='./stat.inp',status='OLD',err=99) read(30,*) num_stat ! read number of stations ! do ns = 1,num_stat read(30,'(a3)' ,iostat=ios) abr call iface enddo close(30) ! stop 99 write(6,*) 'Error opening a file' ! end !--------------------------------------------------------------------------- subroutine iface use param_module,only: ndims, | mlon,mlat,mlev, | mslon,mslat,mslev, | mclon,mclat,mclev,mctime, | start_time,end_time, | calparam,r,rad,abr,sglon,sglat, | work,flnm,mtime use eqvcur_module,only: vallat,vallon, | kqlam_org,kqphi_org, | kqlam,kqphi use nc_module,only: noid,start3,rm_start3,ncdefpdim,dim3,count3 use coef_module,only: icm3e ! implicit none ! #ifdef SUN #include "/opt/local/include/netcdf.inc" #else #include "netcdf.inc" #endif ! integer :: idmtime,dmtime,imtime,idnmtime,itime integer :: dtime,dlon,dlat integer :: idtime,idday,idlon,idlat,idec1,ios,iday ! integer :: ncid,istat,id,tid,lonid,latid integer :: i,it,j,k,id1,id2,ilon,ilat,iyr,len_fln integer :: start4(ndims),count4(ndims) integer :: start2(2),count2(2),imagpot,imagprt,iplctr integer :: start1,count1,stride(ndims-1) character(len=120) :: path character(len=120) :: flnmout character :: cdum*3,output*11,file*7 real :: vals(mclon,mclat,mctime) real :: time ! data stride/1,1,1/ ! ! get station abreviation open(20,file='./statloc',status='OLD',err=99) do i = 1,3 read(20,*) enddo 88 read(20,'(a3,29x,f6.2,f7.2)' ,iostat=ios) cdum,sglat,sglon if(cdum.eq.abr) goto 77 if(ios.ge.0) goto 88 77 close(20) if(sglon > 180.) sglon = sglon-360. write(6,*) 'stat ',cdum,sglat,sglon ! if(cdum.ne.abr) then ! abreviation not found write(6,*) 'station ',abr,' not found' stop endif ! ! open output file for magnetic perturbation at a station output = './output.d/' file(1:3) = abr file(4:7) = '.out' path = output//flnm(1:7)//'_'//file open(10,file=path,status='old',iostat=ios,position='append') if(ios.eq.0) write(6,*) 'output in old file ',path if(ios.gt.0) then open(10,file=path,status='new',iostat=ios,ERR=99) write(6,*) 'output in new file ',path endif ! ! define start and count arrays start4(1) = mslon start4(2) = mslat start4(3) = mslev start4(4) = start_time count4(1) = mclon count4(2) = mclat count4(3) = mclev count4(4) = mctime start3(1) = mslon start3(2) = mslat start3(3) = start_time rm_start3(1) = 1 rm_start3(2) = 1 rm_start3(3) = start_time start2(1) = 1 start2(2) = start_time count2(1) = 3 count2(2) = 1 ! name of datafile path = work//flnm path = trim(path) write(6,*) 'opening nc-file:',path ! specify year of run is used for the apex-files (has to be same year) ! year is read from the netcdf file ! open datafile istat = nf_open(path,NF_NOWRITE,ncid) if (istat /= NF_NOERR) then call check_err(istat,'error w/ open') stop endif write(6,*) 'data read from file: ',flnm ! get dimensions id & length istat = nf_inq_dimid(ncid,'time',idtime) if (istat /= NF_NOERR) call check_err(istat,'id time') istat = nf_inq_dimid(ncid,'mlon',idlon) if (istat /= NF_NOERR) call check_err(istat,'id mlon') istat = nf_inq_dimid(ncid,'mlat',idlat) if (istat /= NF_NOERR) call check_err(istat,'id mlat') istat = nf_inq_dimlen(ncid,idtime,dtime) if (istat /= NF_NOERR) call check_err(istat,'dim time') istat = nf_inq_dimlen(ncid,idlon,dlon) if (istat /= NF_NOERR) call check_err(istat,'dim mlon') istat = nf_inq_dimlen(ncid,idlat,dlat) if (istat /= NF_NOERR) call check_err(istat,'dim mlat') ! check global dimensions dmtime = end_time-start_time + 1 if(dmtime > dtime) write(6,*) * 'error in dim of time: data_t ', * dtime,' time =',dmtime if(mlon /= dlon) write(6,*)'error in dim of lon: data_lon ', * dlon, ' mlon =',mlon if(mlat /= dlat) write(6,*)'error in dim of lat: data_lat ', * dlat, ' mlat =',mlat ! check local dimensions with start value if(mctime+dmtime-1 > dtime) write(6,*) * 'error in dim of time: ', * ' chosen time =',mctime+dmtime-1 if(mclon+mslon-1 > dlon) write(6,*)'error in dim of lon ', * 'chosen lon =',mclon+mslon if(mclat+mslat-1 > dlat) write(6,*)'error in dim of lat ', * 'chosen lat =',mclat+mslat ! create output-file len_fln = len(flnm)-3 flnmout(5:4+len_fln) = flnm(1:len_fln) flnmout(1:4) = 'mag_' flnmout(5+len_fln:12+len_fln)= '_qnm.nc' flnmout(13+len_fln:) = ' ' path = work//"maggrd.d/"//flnmout write(6,*) 'output nc-file created:',path istat = nf_create(path,NF_CLOBBER,noid) if (istat /= NF_NOERR) call check_err(istat,'create nc-file ') istat = nf_def_dim(noid,'time',NF_UNLIMITED,idtime) if (istat /= NF_NOERR) call check_err(istat,'defdim time ') istat = nf_def_dim(noid,'mlon',dlon,idlon) if (istat /= NF_NOERR) call check_err(istat,'defdim mlon ') istat = nf_def_dim(noid,'mlat',dlat,idlat) if (istat /= NF_NOERR) call check_err(istat,'defdim mlat ') istat = nf_def_dim(noid,'mtime',3,idnmtime) if (istat /= NF_NOERR) call check_err(istat,'defdim mtime ') ! ids of 3dim array on irregular grid dim3(1) = idlon dim3(2) = idlat dim3(3) = idtime ! istat = nf_def_var(noid,'time',NF_REAL,1,idtime,itime) if (istat /= NF_NOERR) call check_err(istat,'defvar mlon ') ! istat = nf_def_var(noid,'mlon',NF_REAL,1,idlon,ilon) if (istat /= NF_NOERR) call check_err(istat,'defvar mlon ') ! istat = nf_def_var(noid,'mlat',NF_REAL,1,idlat,ilat) if (istat /= NF_NOERR) call check_err(istat,'defvar mlat ') ! istat = nf_def_var(noid,'mtime',NF_REAL,2,(/idnmtime,idtime/), | imtime) if (istat /= NF_NOERR) call check_err(istat,'defvar mtime ') ! istat = nf_enddef(noid) if (istat /= NF_NOERR) call check_err(istat,'enddef ') ! ! get latitudes istat = nf_inq_varid(ncid,'mlat',id1) if (istat /= NF_NOERR) call check_err(istat,'id mlat ') start1 = 1 count1 = mlat istat = nf_get_vara_double(ncid,id1,start1,count1,vallat) if (istat /= NF_NOERR) call check_err(istat,'vara mlat ') ! put latitude in new file istat = nf_put_var_double(noid,ilat,vallat) if (istat /= NF_NOERR) call check_err(istat,'putvar mlat ') ! get longitudes istat = nf_inq_varid(ncid,'mlon',id2) if (istat /= NF_NOERR) call check_err(istat,'id mlon ') start1 = 1 count1 = mlon istat = nf_get_vara_double(ncid,id2,start1,count1,vallon) if (istat /= NF_NOERR) call check_err(istat,'vara mlon ') ! put longitude in new file istat = nf_put_var_double(noid,ilon,vallon) if (istat /= NF_NOERR) call check_err(istat,'putvar mlon ') ! get mtime id istat = nf_inq_varid(ncid,'mtime',idmtime) if (istat /= NF_NOERR) call check_err(istat,'id mtime ') ! get day id istat = nf_inq_varid(ncid,'day',idday) if (istat /= NF_NOERR) call check_err(istat,'id day ') ! get time id istat = nf_inq_varid(ncid,'time',idtime) if (istat /= NF_NOERR) call check_err(istat,'id time ') ! copy attributes to new file for istat = nf_redef(noid) if (istat /= NF_NOERR) call check_err(istat,'redef ') istat = nf_copy_att(ncid,idtime,'units',noid,itime) if (istat /= NF_NOERR) call check_err(istat,'copy att_time unit') istat = nf_copy_att(ncid,idtime,'long_name',noid,itime) if (istat /= NF_NOERR) call check_err(istat,'copy att_time name') istat = nf_copy_att(ncid,idtime,'initial_mtime',noid,itime) if (istat /= NF_NOERR) call check_err(istat,'copy att_mtime name') istat = nf_copy_att(ncid,idtime,'initial_day',noid,itime) if (istat /= NF_NOERR) call check_err(istat,'copy att_day name') istat = nf_copy_att(ncid,idtime,'initial_year',noid,itime) if (istat /= NF_NOERR) call check_err(istat,'copy att_year name') istat = nf_get_att_int(ncid,idtime,'initial_year',iyr) ! put year for apex if (istat /= NF_NOERR) call check_err(istat,'get att_year name') istat = nf_put_att_text(noid,ilon,'units',21, | 'degrees (south north)') if (istat /= NF_NOERR) call check_err(istat,'put att_lon unit') istat = nf_put_att_text(noid,ilat,'units',19, | 'degrees (west east)') if (istat /= NF_NOERR) call check_err(istat,'put att_lon unit') istat = nf_copy_att(ncid,id1,'long_name',noid,ilat) if (istat /= NF_NOERR) call check_err(istat,'copy att_lat name') istat = nf_copy_att(ncid,id2,'long_name',noid,ilon) if (istat /= NF_NOERR) call check_err(istat,'copy att_lon name') istat = nf_copy_att(ncid,idmtime,'units',noid,imtime) if (istat /= NF_NOERR) call check_err(istat, | 'copy att_mtime units') istat = nf_copy_att(ncid,idmtime,'long_name',noid,imtime) if (istat /= NF_NOERR) call check_err(istat, | 'copy att_mtime name') istat = nf_enddef(noid) if (istat /= NF_NOERR) call check_err(istat,'enddef ') ! prepare parameters and output file call calparam ! calculate parameters call ncdefpdim ! nc-output: define plot dimensions & attributes if(icm3e == 1) call sphere ! cm3e relates internal to external component ! reads data for one timestep (time=1) and level(mslev=1) ! if more timesteps a loop is needed count 1 = 1 do it = start_time,end_time ! update start arrays start2(2) = it start3(3) = it rm_start3(3) = it start4(4) = it ! get model time from data file istat = nf_get_vara_int(ncid,idmtime,start2,count2,mtime) if (istat /= NF_NOERR) call check_err(istat,'vara mtime ') ! model day is not the calender day (e.g. run the same day- mtime changing) start1 = it istat = nf_get_vara_int(ncid,idday,start1,count1,iday) if (istat /= NF_NOERR) call check_err(istat,'vara day ') mtime(1) = iday ! put mtime in new file istat = nf_put_vara_int(noid,idnmtime,(/1,it-start_time+1/), | count2,mtime) if (istat /= NF_NOERR) call check_err(istat,'putvar mtime ') ! get time from data file istat = nf_get_vara_double(ncid,idtime,it,1,time) if (istat /= NF_NOERR) call check_err(istat,'vara time ') ! put time in new file istat = nf_put_vara_double(noid,itime,it-start_time+1, | 1,time) if (istat /= NF_NOERR) call check_err(istat,'putvar time ') ! ! K_{q,lam}: horizontal height integrated current density ! stored in secondary history field KQLAM ! get id & datavalues ! northward component positive (Richmond 1994) istat = nf_inq_varid(ncid,'KQLAM',id) if (istat /= NF_NOERR) call check_err(istat,'id KQLAM ') istat = nf_get_vara_double(ncid,id,start3,count3,vals) if (istat /= NF_NOERR) call check_err(istat,'vara KQLAM ') do i = 1,mclon do j = 1,mclat kqlam_org(i,j) = vals(i,j,1) enddo enddo ! K_{q,phi}: horizontal height integrated current density ! stored in secondary history field KQPHI ! get id & datavalues istat = nf_inq_varid(ncid,'KQPHI',id) if (istat /= NF_NOERR) call check_err(istat,'id KQPHI ') istat = nf_get_vara_double(ncid,id,start3,count3,vals) if (istat /= NF_NOERR) call check_err(istat,'vara KQPHI ') do i = 1,mclon do j = 1,mclat kqphi_org(i,j) = vals(i,j,1) enddo enddo if (istat /= NF_NOERR) write(6,*) 'error with closing' ! interpolate kqlam_org, kqphi_org,jmr_org from the irregular magnetic grid to ! the regular magnetic grid ! output: kqlam, kqphi, jmr call mapirr2reg ! calculate mag. longitude of subsolar point for model time ! year must be specified explicitly ! input : mtime,iyr ! calculate : xmlon,date iyr = 2000 call subsollon(mtime,iyr) start3(3) = it-start_time+1 rm_start3(3) = it-start_time+1 ! calculate the magnetic perturbation at the ground call maggrd enddo ! for time loop ! close datafile istat = nf_close(ncid) if (istat /= NF_NOERR) call check_err(istat,'closing data file ') ! close the new created file istat = nf_close(noid) if (istat /= NF_NOERR) call check_err(istat,'closing new file ') ! close(10) ! for mag.perturbation at station ! return 99 write(6,*) 'Error opening a file' end subroutine iface !-------------------------------------------------------------------------- subroutine check_err(status,name) ! implicit none ! #ifdef SUN #include "/opt/local/include/netcdf.inc" #else #include "netcdf.inc" #endif ! integer,intent(in) :: status character(len=*),intent(in) :: name write(6,*) 'error in ',name write(6,*) NF_STRERROR(status) write(6,*) ' ' stop end subroutine check_err !--------------------------------------------------------------------------