! module mk_drifts use input,only: tmpdir,ncfile_dipdec,ncfile_magfield use hist,only: history use den_convert,only: denconv,barm,pkt implicit none include "netcdf.inc" ! ! Module to calculate ion velocities. ! ! (ionvel==5 added, and others corrected by Ganglu, 9/07) ! ionvel = 0 -> read ion velocities from the (secondary) history ! ionvel = 1 -> ExB velocity ! ionvel = 2 -> Ion velocity (ExB+unvn) ! ionvel = 3 -> Neutral wind component of ion velocity (efield=0) ! ionvel = 4 -> Electric field component of ion velocity (unvn=0) ! ionvel = 5 -> Vertical ion drift due to ion diffusion (UI,VI as in ionvel=2) ! ! Dipdec and neutral atmos are used only when ionvel > 1 ! real,allocatable,save :: + bx(:,:),by(:,:),bz(:,:), ! mag field + epot(:,:,:), ! epot + hts(:,:,:), ! heights + phi(:,:,:), ! epot w/ poles and wrap pts + dipdec(:,:,:), ! dip and declination (radians) + un(:,:,:),vn(:,:,:),wn(:,:,:), ! neutral velocities + xo2(:,:,:),xo1(:,:,:),xn2(:,:,:), ! neutral densities + te(:,:,:),ti(:,:,:),tn(:,:,:), ! temperatures + electron(:,:,:) ! electron density integer,external :: ixfindc ! real,parameter :: xmaso=32., wto2=32.,wto1=16.,wtn2=28. ! btf 6/23/04, as per Art: real,parameter :: xmaso=31.5, wto2=32.,wto1=16.,wtn2=28. character(len=NF_MAX_NAME) :: name contains !------------------------------------------------------------------- subroutine savephi(flat,h,imx,kmx,jmx,lat,ionvel) ! ! Save phi(imxp3,kmx,jmxp2) for current lat. ! Also save neutrals u,v,w,o2,o,n2 if ionvel > 1. ! (flat must be intent(inout) only for dummy arg in denconv -- ! it is not actually changed) ! ! Args: integer,intent(in) :: imx,kmx,jmx,lat,ionvel type(history),intent(in) :: h real,intent(inout) :: flat(imx,kmx,h%nflds) ! ! Locals: integer :: ier,ixpoten,ixz,i,ncalls=0,j integer :: ixun,ixvn,ixwn,ixo2,ixo1,ixte,ixti,ixtn,ixne integer :: jmxm1,jmxp2,imxm1,imxp1,imxp3 character(len=8) :: dunits real :: fmin,fmax ! ncalls = ncalls+1 jmxm1 = jmx-1 jmxp2 = jmx+2 imxm1 = imx-1 imxp1 = imx+1 imxp3 = imx+3 if (.not.allocated(phi)) then allocate(phi(imxp3,kmx,jmxp2),stat=ier) if (ier /= 0) call allocerr(ier, + "allocate 3d phi from savephi.") endif if (.not.allocated(epot)) then allocate(epot(imx,kmx,jmx),stat=ier) if (ier /= 0) call allocerr(ier, + "allocate 3d epot from savephi.") endif if (.not.allocated(hts)) then allocate(hts(imx,kmx,jmx),stat=ier) if (ier /= 0) call allocerr(ier, + "allocate 3d hts from savephi.") endif ! ! Save 3d epot and hts: ixpoten = ixfindc(h%fnames,h%nflds,'POTEN ') if (ixpoten <= 0) then write(6,"('>>> savephi: need POTEN for ion drifts.')") stop 'savephi' endif ixz = ixfindc(h%fnames,h%nflds,'Z ') if (ixz <= 0) then write(6,"('>>> savephi: need Z for ion drifts.')") stop 'savephi' endif epot(:,:,lat) = flat(:,:,ixpoten) hts(:,:,lat) = flat(:,:,ixz) ! if (ionvel /= 3) then phi(3:imxp1,:,lat+1) = flat(1:imxm1,:,ixpoten) if (lat == 1) then do i=1,imxm1 phi(i+2,:,lat) = flat(1+mod(i+jmxm1,imxm1),:,ixpoten) enddo elseif (lat == jmx) then do i=1,imxm1 phi(i+2,:,jmxp2) = flat(1+mod(i+jmxm1,imxm1),:,ixpoten) enddo endif ! ! Periodic points: 1,2<-73,74 and 75,76<-3,4 if (ncalls == jmx) then do j=1,jmxp2 phi(1:2,:,j) = phi(imxm1+1:imxm1+2,:,j) phi(imxp1+1:imxp1+2,:,j) = phi(3:4,:,j) enddo endif else ! ionvel==3 -> set phi==0 phi = 0. ! whole array expression endif if (ncalls == jmx) ncalls = 0 ! reset for next history ! ! Save neutrals if necessary: if (ionvel > 1) then if (.not.allocated(un)) then allocate(un(imx,kmx,jmx),stat=ier) if (ier /= 0) call allocerr(ier, + "allocate 3d un from savephi.") endif if (.not.allocated(vn)) then allocate(vn(imx,kmx,jmx),stat=ier) if (ier /= 0) call allocerr(ier, + "allocate 3d vn from savephi.") endif if (.not.allocated(wn)) then allocate(wn(imx,kmx,jmx),stat=ier) if (ier /= 0) call allocerr(ier, + "allocate 3d wn from savephi.") endif if (.not.allocated(te)) then allocate(te(imx,kmx,jmx),stat=ier) if (ier /= 0) call allocerr(ier, + "allocate 3d te from savephi.") endif if (.not.allocated(ti)) then allocate(ti(imx,kmx,jmx),stat=ier) if (ier /= 0) call allocerr(ier, + "allocate 3d ti from savephi.") endif if (.not.allocated(tn)) then allocate(tn(imx,kmx,jmx),stat=ier) if (ier /= 0) call allocerr(ier, + "allocate 3d tn from savephi.") endif if (.not.allocated(electron)) then allocate(electron(imx,kmx,jmx),stat=ier) if (ier /= 0) call allocerr(ier, + "allocate 3d NE from savephi.") endif if (.not.allocated(xo2)) then allocate(xo2(imx,kmx,jmx),stat=ier) if (ier /= 0) call allocerr(ier, + "allocate 3d xo2 from savephi.") endif if (.not.allocated(xo1)) then allocate(xo1(imx,kmx,jmx),stat=ier) if (ier /= 0) call allocerr(ier, + "allocate 3d xo from savephi.") endif if (.not.allocated(xn2)) then allocate(xn2(imx,kmx,jmx),stat=ier) if (ier /= 0) call allocerr(ier, + "allocate 3d xn2 from savephi.") endif ixun = ixfindc(h%fnames,h%nflds,'UN ') ixvn = ixfindc(h%fnames,h%nflds,'VN ') ixwn = ixfindc(h%fnames,h%nflds,'W ') if (ixwn==0) ixwn = ixfindc(h%fnames,h%nflds,'OMEGA ') ixo2 = ixfindc(h%fnames,h%nflds,'O2 ') ixo1 = ixfindc(h%fnames,h%nflds,'O1 ') ixte = ixfindc(h%fnames,h%nflds,'TE ') ixti = ixfindc(h%fnames,h%nflds,'TI ') ixtn = ixfindc(h%fnames,h%nflds,'TN ') ixne = ixfindc(h%fnames,h%nflds,'NE ') if (ixun==0.or.ixvn==0.or.ixwn==0.or.ixo2==0.or.ixo1==0 .or. + ixte==0.or.ixti==0.or.ixtn==0 .or. ixne==0) then write(6,"('>>> isavephi: need neutrals for ionvel=',i2, + ': ixun,vn,wn,o2,o=',9i5)") ionvel,ixun,ixvn,ixwn,ixo2,ixo1, + ixte,ixti,ixtn,ixne stop 'ionvel' endif ! ! Get neutral winds. If ionvel==4, set neutral winds = 0. ! if (ionvel /= 4) then un(:,:,lat) = flat(:,:,ixun) vn(:,:,lat) = flat(:,:,ixvn) wn(:,:,lat) = flat(:,:,ixwn) else un = 0. vn = 0. wn = 0. endif te(:,:,lat) = flat(:,:,ixte) ti(:,:,lat) = flat(:,:,ixti) tn(:,:,lat) = flat(:,:,ixtn) electron(:,:,lat) = flat(:,:,ixne) ! call fminmax(un(:,:,lat),imx*kmx,fmin,fmax) ! write(6,"('UN: min,max=',2e12.4)") fmin,fmax ! call fminmax(wn(:,:,lat),imx*kmx,fmin,fmax) ! write(6,"('WN: min,max=',2e12.4)") fmin,fmax ! call fminmax(te(:,:,lat),imx*kmx,fmin,fmax) ! write(6,"('TE: min,max=',2e12.4)") fmin,fmax ! call fminmax(tn(:,:,lat),imx*kmx,fmin,fmax) ! write(6,"('TN: min,max=',2e12.4)") fmin,fmax ! call fminmax(electron(:,:,lat),imx*kmx,fmin,fmax) ! write(6,"('NE: min,max=',2e12.4)") fmin,fmax ! ! Get neutral densities in cm3 ! (note barm and pkt have already been defined by mkdenparms) ! (note flat not changed by denconv) ! call denconv(flat(:,:,ixo2),xo2(:,:,lat), + 1,barm,pkt,wto2,imx,kmx,-1,'MMR','CM3',dunits) call denconv(flat(:,:,ixo1),xo1(:,:,lat), + 1,barm,pkt,wto1,imx,kmx,-1,'MMR','CM3',dunits) xn2(:,:,lat) = max(.00001,(1.-flat(:,:,ixo2)-flat(:,:,ixo1))) call denconv(xn2(:,:,lat),xn2(:,:,lat), + 0,barm,pkt,wtn2,imx,kmx,-1,'MMR','CM3',dunits) endif ! ionvel > 1 return end subroutine savephi !------------------------------------------------------------------- subroutine mkdrifts(ui,mkui,vi,mkvi,wi,mkwi, | h,imx,kmx,jmx,lat,ionvel) use proc,only: dlat,dlon,gcmlat,spval ! ! Calculate ion drift velocities at current latitude: ! (3d arrays phi, epot, and hts in this module were saved by savephi) ! ! Args: real,intent(out) :: ui(imx,kmx),vi(imx,kmx),wi(imx,kmx) ! real,intent(in) :: epot(imx,kmx) type(history),intent(in) :: h integer,intent(in) :: imx,kmx,jmx,lat,ionvel logical,intent(in) :: mkui,mkvi,mkwi ! ! Locals: real,parameter :: pi=3.14159, dtr=pi/180., rs=6.475165e+8, + scale=1.e6, bgaus=.4, colfac=1.5, grav=870. integer :: ncalls=0,ixpoten,i,k real :: fmin0,fmax0, fmin1,fmax1, fmin2,fmax2 real :: ex(kmx),ey(kmx),ez(kmx),b(imx) integer,external :: ixfindc real :: hh,cosd,sind,si,ci,ve,vu,uik,vik,wik, + tun,tvn,twn,to1,to2,tn2,vinp,rhoi,coi,unmag,vnmag, + rnu_o2p_o2,rnu_op_o2,rnu_nop_o2,rnu_nop_o,rnu_op_o, + rnu_o2p_o,rnu_nop_n2,rnu_o2p_n2,rnu_op_n2,nop_n2, + rnu_o2p,rnu_op,rnu_nop,tnti,tite,H_p,D_a, + dndz,dtdz real :: wi_ion(imx,kmx) ! ncalls = ncalls+1 ! ! Get mag field. Also dip and declination if necessary: ! 5/20/02: read netcdf file for mag field data. ! if (ncalls == 1) then ! ! Sub get magnc reads netcdf files for either 2.5 or 5.0 res. call getmagnc(imx,jmx) if (ionvel > 1) call getdipdec(imx,jmx) endif ! ! Calculate electric field: ! (epot is electric potential at current latitude, ! phi is 3d epot defined earlier w/ poles and periodic points) ! do i=1,imx if (ionvel /= 3) then ex(:) = -(phi(i+3,:,lat+1) - phi(i+1,:,lat+1)) / + (2.*dlon*dtr*(rs+hts(i,:,lat)*1.e5)*cos(gcmlat(lat)*dtr)) ey(:) = -(phi(i+2,:,lat+2) - phi(i+2,:,lat)) / + (2.*dlat*dtr*(rs+hts(i,:,lat)*1.e5)) ez(2:kmx-1) = -((epot(i,3:kmx,lat)-epot(i,2:kmx-1,lat)) / + ((hts(i,3:kmx,lat)-hts(i,2:kmx-1,lat))*1.e5) + + (epot(i,2:kmx-1,lat)-epot(i,1:kmx-2,lat)) / + ((hts(i,2:kmx-1,lat)-hts(i,1:kmx-2,lat))*1.e5)) / 2. ez(1) = ez(2)*2.-ez(3) ez(kmx) = ez(kmx-1)*2.-ez(kmx-2) ! ! ionvel == 3 -> set electric field 0 (would be calculated to 0 ! anyway because phi==0, but whole array ops are faster): else ex = 0. ey = 0. ez = 0. endif ! ! Define all 3 components of ion velocities: ! b(i) = (bx(i,lat)**2+by(i,lat)**2+bz(i,lat)**2)**0.5 ! ! ExB only: if (ionvel == 1) then if (mkui) | ui(i,:) = (ey(:)*bz(i,lat)-ez(:)*by(i,lat))*scale/b(i)**2 if (mkvi) | vi(i,:) = (ez(:)*bx(i,lat)-ex(:)*bz(i,lat))*scale/b(i)**2 if (mkwi) | wi(i,:) = (ex(:)*by(i,lat)-ey(:)*bx(i,lat))*scale/b(i)**2 ! ! Ion velocity w/ neutrals (3d neutrals were saved by savephi): ! (ionvel may be 2, 3, or 4) else hh = sqrt(bx(i,lat)**2+by(i,lat)**2) cosd = cos(dipdec(i,lat,2)) sind = sin(dipdec(i,lat,2)) si = sin(dipdec(i,lat,1)) ci = cos(dipdec(i,lat,1)) do k=1,kmx uik = (ey(k)*bz(i,lat)-ez(k)*by(i,lat))*scale/b(i)**2 vik = (ez(k)*bx(i,lat)-ex(k)*bz(i,lat))*scale/b(i)**2 wik = (ex(k)*by(i,lat)-ey(k)*bx(i,lat))*scale/b(i)**2 ! ve is the esatward component of perpendicular ExB drift ve = (uik*by(i,lat)-vik*bx(i,lat)) / hh ! vn is the vertical component of perpendicular ExB drift ! vu = -((uik*bx(i,lat)+vik*by(i,lat))/(hh*b(i))*bz(i,lat))+ ! + hh/b(i)*wik vu = b(i)*wik/hh tun = un(i,k,lat) tvn = vn(i,k,lat) twn = wn(i,k,lat) to1= xo1(i,k,lat) to2 = xo2(i,k,lat) tn2 = xn2(i,k,lat) if (to1.ne.spval.and.tn2.ne.spval.and.to2.ne.spval.and. + tun.ne.spval.and.tvn.ne.spval.and.twn.ne.spval) then ! vinp is the ion-neutral collision frequency ! vinp = 3.e-10 * (to1+tn2+to2) ! using the same calculation at in lamdas.F ! O2 collision frequencies: tnti = 0.5*(ti(i,k,lat)+tn(i,k,lat)) rnu_o2p_o2 = 2.59E-11*sqrt(tnti)* | (1.-0.073*alog10(tnti))**2 rnu_op_o2 = 6.64E-10 rnu_nop_o2 = 4.27E-10 ! O collision frequencies: rnu_o2p_o = 2.31E-10 rnu_op_o = 3.67e-11*sqrt(tnti)* ! O+ ~ O (resonant) | (1.-0.064*alog10(tnti))**2*colfac rnu_nop_o = 2.44E-10 ! N2 collision frequencies: rnu_o2p_n2 = 4.13E-10 ! O2+ ~ N2 rnu_op_n2 = 6.82E-10 ! O+ ~ N2 rnu_nop_n2 = 4.34E-10 rnu_o2p = rnu_o2p_o2*xo2(i,k,lat)+rnu_o2p_o*xo1(i,k,lat)+ + rnu_o2p_n2*xn2(i,k,lat) rnu_op = rnu_op_o2*xo2(i,k,lat)+rnu_op_o*xo1(i,k,lat)+ + rnu_op_n2*xn2(i,k,lat) rnu_nop = rnu_nop_o2*xo2(i,k,lat)+rnu_nop_o*xo1(i,k,lat)+ + rnu_nop_n2*xn2(i,k,lat) vinp = rnu_o2p+rnu_op+rnu_nop tite = 0.5*(ti(i,k,lat)+te(i,k,lat)) ! Hp is the plasma scale height and using O+ mass H_p = 2.*1.38e-16*tite/(16.*1.673e-24*grav) ! D_a is the ion diffusion coefficient D_a = 2.*1.38e-16*tite/(16.*1.763e-24*vinp) ! wi_ion is vertical ion drift due to ion diffusion if (k > 1 .and. k < kmx) then dndz = (electron(i,k+1,lat)-electron(i,k-1,lat))/ + (hts(i,k+1,lat)-hts(i,k-1,lat))/1.e5 dtdz = 0.5*(ti(i,k+1,lat)+te(i,k+1,lat)- + ti(i,k-1,lat)-te(i,k-1,lat))/ + (hts(i,k+1,lat)-hts(i,k-1,lat))/1.e5 endif ! Z (hts) is in km - convert to cm by a factor of 1.e5 if (k == 1) then dndz = (electron(i,2,lat)-electron(i,1,lat))/ + (hts(i,2,lat)-hts(i,1,lat))/1.e5 dtdz = 0.5*(ti(i,2,lat)+te(i,2,lat)- + ti(i,1,lat)-te(i,1,lat))/ + (hts(i,2,lat)-hts(i,1,lat))/1.e5 endif if (k == kmx) then dndz = (electron(i,kmx,lat)-electron(i,kmx-1,lat))/ + (hts(i,kmx,lat)-hts(i,kmx-1,lat))/1.e5 dtdz = 0.5*(ti(i,kmx,lat)+te(i,kmx,lat)- + ti(i,kmx-1,lat)-te(i,kmx-1,lat))/ + (hts(i,kmx,lat)-hts(i,kmx-1,lat))/1.e5 endif if (k < kmx) then ! tite==0 at kmx wi_ion(i,k) = -si*si*D_a*(dndz/electron(i,k,lat) + + dtdz/tite + 1./H_p) ! convert wi_ion from cm/s to m/s wi_ion(i,k) = wi_ion(i,k)/100. endif ! if ((i==10.and.k==kmx).or.(i==10.and.k==kmx-1)) ! | write(6,"('i,k,lat=',3i4,' Z,tite,barm,vinp,H_p,D_a,dndz,dtdz', ! | ' = ',8e12.3,' wi_ion=',d12.4)") i,k,lat,hts(i,k,lat),tite, ! | barm(i,k),vinp,H_p,D_a,dndz,dtdz,wi_ion(i,k) ! and rhoi is the ratio of ion-neutral collision freq and ion gyro freq ! rhoi = vinp / (9.57946e+3 * b(i) / xmaso) rhoi = vinp / (9.57946e+3 * b(i) / barm(i,k)) coi = 1. / (1. + rhoi*rhoi) unmag = tun*cosd - tvn*sind vnmag = tun*sind + tvn*cosd if (mkui) | ui(i,k) =coi*(rhoi*rhoi*unmag+rhoi*(-vnmag*si-twn*ci))+ | coi*(ve+rhoi*vu) if (mkvi) | vi(i,k) = coi*(rhoi*rhoi*vnmag+rhoi*unmag*si+ | vnmag*ci*ci-twn*si*ci)+coi*(vu*si-rhoi*ve*si) if (mkwi) then wi(i,k) = coi*(rhoi*rhoi*twn+rhoi*unmag*ci-vnmag*si*ci+ | twn*si*si)+coi*(vu*ci-rhoi*ve*ci) if (k < kmx) then if (ionvel==2) wi(i,k)=wi(i,k)+wi_ion(i,k) if (ionvel==5) wi(i,k)=wi_ion(i,k) else ! wi(i,k) = spval wi(i,k) = wi(i,k-1) endif endif else ! a term is spval if (mkui) ui(i,k) = spval if (mkvi) vi(i,k) = spval if (mkwi) wi(i,k) = spval endif enddo ! k=1,kmx ! ! 11/24/11 btf: double-res runs have problems w/ top 2 k values of wi: ! (minimum is often -inf at mid to high latitudes) ! do k=kmx-1,kmx wi(i,k) = wi(i,kmx-2) enddo ! write(6,"('lat=',i4,' i=',i4,' kmx=',i4,' wi(i,:)=',/, ! | (6e12.4))") lat,i,kmx,wi(i,:) endif ! ionvel enddo ! i=1,imx ! ! Deallocate per-history arrays used in ion velocities: ! (reset ncalls for next history) ! if (ncalls == jmx) then deallocate(phi) deallocate(epot) deallocate(hts) if (ionvel > 1) then deallocate(un) deallocate(vn) deallocate(wn) deallocate(xo1) deallocate(xo2) deallocate(xn2) deallocate(te) deallocate(ti) deallocate(tn) deallocate(electron) endif ncalls = 0 endif return end subroutine mkdrifts !------------------------------------------------------------------- subroutine getdipdec(imx,jmx) include "netcdf.inc" ! ! Read dipdec data from mss netcdf file: ! ! Args: integer,intent(in) :: imx,jmx ! ! Local: character(len=1024) :: dskfile character(len=120) :: char120 integer :: ier,ncid,istat,id_lon,id_lat,nlon,nlat,idv_dip,idv_dec, | icount2(2),istart2(2),ids2(2),i,j real :: fmin,fmax ! if (allocated(dipdec)) return allocate(dipdec(imx,jmx,2),stat=ier) if (ier /= 0) call allocerr(ier, | "allocate dipdec from getdipdec.") ! ! Get dipdec field from mss file: if (len_trim(ncfile_dipdec)==0) then if (jmx==36) then ! 5-deg res ncfile_dipdec = '$TGCMDATA/dipdec_5.0h.nc' elseif (jmx==72) then ! 2.5-deg res ncfile_dipdec = '$TGCMDATA/dipdec_2.5h.nc' else write(6,"('>>> getmagnc: unknown jmx=',i4,' -- need', | ' ncfile_dipdec')") jmx stop 'ncfile_dipdec' endif endif dskfile = ' ' call expand_path(ncfile_dipdec) call getfile(ncfile_dipdec,dskfile) ! ! Open dipdec file for reading: istat = nf_open(dskfile,NF_NOWRITE,ncid) if (istat /= NF_NOERR) then write(char120,"('>>> getdipdec: Error return from nf_open for ', | 'netcdf file ',a,' status=',a)") trim(dskfile),istat call handle_ncerr(istat,char120) stop 'dipdec.nc open' else write(6,"('getdipdec: Opened dipdec data file ',a, | ' ncid=',i3)") trim(ncfile_dipdec),ncid endif ! ! Read and check dimensions: istat = nf_inq_dimid(ncid,'lon',id_lon) if (istat /= NF_NOERR) then write(char120,"('>>> getdipdec: error obtaining lon ', | 'dim id from netcdf file ',a,' status=',a)") trim(dskfile), | istat call handle_ncerr(istat,char120) stop 'magfield file' endif istat = nf_inq_dim(ncid,id_lon,name,nlon) if (nlon /= imx) then write(6,"(/,'>>> getdipdec: length of lon read from', | ' history does not match parameter nlon.')") write(6,"(' nlon (read)=',i4,' imx (param)=',i4)") nlon,imx stop 'getdipdec' endif istat = nf_inq_dimid(ncid,'lat',id_lat) if (istat /= NF_NOERR) then write(char120,"('>>> getdipdec: error obtaining lat ', | 'dim id from netcdf file ',a,' status=',a)") trim(dskfile), | istat call handle_ncerr(istat,char120) stop 'magfield file' endif istat = nf_inq_dim(ncid,id_lat,name,nlat) if (nlat /= jmx) then write(6,"(/,'>>> getdipdec: length of lat read from', | ' history does not match parameter nlat.')") write(6,"(' nlat (read)=',i4,' jmx (param)=',i4)") nlat,jmx stop 'getdipdec' endif ! ! Read dip: istat = nf_inq_varid(ncid,'DIP',idv_dip) istart2(:) = 1 icount2(1) = imx icount2(2) = jmx istat = nf_get_vara_real(ncid,idv_dip,istart2,icount2, | dipdec(1,1,1)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error return from nf_get_vara_real for var dip') ! do j=1,jmx ! write(6,"('getdipdec: j=',i2,' dipdec(:,j,1)=',/,(6e12.4))") ! | j,dipdec(:,j,1) ! enddo call fminmax(dipdec,imx*jmx,fmin,fmax) write(6,"('getdipdec: dip min,max=',2e12.4)") fmin,fmax ! ! Read dec: istat = nf_inq_varid(ncid,'DEC',idv_dec) istart2(:) = 1 icount2(1) = imx icount2(2) = jmx istat = nf_get_vara_real(ncid,idv_dec,istart2,icount2, | dipdec(1,1,2)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error return from nf_get_vara_real for var dec') ! do j=1,jmx ! write(6,"('getdipdec: j=',i2,' dipdec(:,j,2)=',/,(6e12.4))") ! | j,dipdec(:,j,2) ! enddo call fminmax(dipdec(1,1,2),imx*jmx,fmin,fmax) write(6,"('getdipdec: dec min,max=',2e12.4,/)") fmin,fmax ! end subroutine getdipdec !------------------------------------------------------------------- subroutine getmagnc(imx,jmx) include "netcdf.inc" ! ! Read magnetic field data from external file: ! (this replaces sub getmag, which reads old cray fortran binary file) ! ! Args: integer,intent(in) :: imx,jmx ! ! Local: character(len=1024) :: dskfile character(len=120) :: char120 integer :: ier,istat,count_2d(2),start_2d(2),idv_bx,idv_by,idv_bz integer :: id_lon,id_lat,nlon,nlat integer :: i,j,ncid real :: fmin,fmax real,dimension(imx-1,jmx) :: rdbx,rdby,rdbz ! ! If bx,by,bz are already allocated (previous history), just return ! if (allocated(bx).and.allocated(by).and.allocated(bz)) | return allocate(bx(imx,jmx),stat=ier) allocate(by(imx,jmx),stat=ier) allocate(bz(imx,jmx),stat=ier) if (ier /= 0) call allocerr(ier, | "allocate bx,by,bz from getmagnc.") ! ! Get magnetic field from mss file: ! If it was not provided in namelist input, try looking in $TGCMDATA. if (len_trim(ncfile_magfield)==0) then if (jmx==36) then ! 5-deg res ncfile_magfield = '$TGCMDATA/magfield_5.0h.nc' elseif (jmx==72) then ! 2.5-deg res ncfile_magfield = '$TGCMDATA/magfield_2.5h.nc' else write(6,"('>>> getmagnc: unknown jmx=',i4,' -- need', | ' ncfile_magfield')") jmx stop 'ncfile_magfield' endif endif dskfile = ' ' call expand_path(ncfile_magfield) call getfile(ncfile_magfield,dskfile) ! ! Open magnetic field data file for reading: istat = nf_open(dskfile,NF_NOWRITE,ncid) if (istat /= NF_NOERR) then write(char120,"('>>> getmagnc: Error return from nf_open for ', | 'netcdf file ',a,' status=',a)") trim(dskfile),istat call handle_ncerr(istat,char120) stop 'magfield file' else write(6,"('getmagnc: Opened magfield data file ',a, | ' ncid=',i3)") trim(ncfile_magfield),ncid endif ! ! Read and check dimensions: ! ! netcdf magfield_5.0h { ! dimensions: ! lon = 73 ; ! lat = 36 ; ! istat = nf_inq_dimid(ncid,'lon',id_lon) if (istat /= NF_NOERR) then write(char120,"('>>> getmagnc: error obtaining lon ', | 'dim id from netcdf file ',a,' status=',a)") trim(dskfile), | istat call handle_ncerr(istat,char120) stop 'magfield file' endif istat = nf_inq_dim(ncid,id_lon,name,nlon) ! ! 9/20/11 btf: magfield_2.5h.nc lon dim=144, but imx=145 ! 11/22/11 btf: This is a problem only w/ 2.5-deg file. ! Really should rewrite magfield_2.5h with lon=145, ! i.e., with periodic point. ! if (nlon == 145 .and. nlon /= imx-1) then write(6,"(/,'>>> getmagnc: length of lon read from', | ' history does not match parameter nlon.')") write(6,"(' nlon (read)=',i4,' imx (param)=',i4)") nlon,imx stop 'getmagnc' endif istat = nf_inq_dimid(ncid,'lat',id_lat) if (istat /= NF_NOERR) then write(char120,"('>>> getmagnc: error obtaining lat ', | 'dim id from netcdf file ',a,' status=',a)") trim(dskfile), | istat call handle_ncerr(istat,char120) stop 'magfield file' endif istat = nf_inq_dim(ncid,id_lat,name,nlat) if (nlat /= jmx) then write(6,"(/,'>>> getmagnc: length of lat read from', | ' history does not match parameter nlat.')") write(6,"(' nlat (read)=',i4,' jmx (param)=',i4)") nlat,jmx stop 'getmagnc' endif ! ! Read by,bx,bz: start_2d(:) = 1 count_2d(1) = imx-1 count_2d(2) = jmx ! ! BX: istat = nf_inq_varid(ncid,"BX",idv_bx) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error return from nf_inq_varid to get var id for bx') istat = nf_get_vara_real(ncid,idv_bx,start_2d,count_2d,rdbx) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error return from nf_get_vara_real for var bx') do j=1,jmx do i=1,imx-1 bx(i,j) = rdbx(i,j) enddo ! i=1,imx-1 bx(imx,j) = bx(1,j) ! periodic point enddo ! j=1,jmx call fminmax(bx,imx*jmx,fmin,fmax) write(6,"('getmagnc: bx min,max=',2e12.4)") fmin,fmax ! ! BY: istat = nf_inq_varid(ncid,"BY",idv_by) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error return from nf_inq_varid to get var id for by') istat = nf_get_vara_real(ncid,idv_by,start_2d,count_2d,rdby) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error return from nf_get_vara_real for var by') do j=1,jmx do i=1,imx-1 by(i,j) = rdby(i,j) enddo ! i=1,imx-1 by(imx,j) = by(1,j) enddo ! j=1,jmx call fminmax(by,imx*jmx,fmin,fmax) write(6,"('getmagnc: by min,max=',2e12.4)") fmin,fmax ! ! BZ: istat = nf_inq_varid(ncid,"BZ",idv_bz) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error return from nf_inq_varid to get var id for bz') istat = nf_get_vara_real(ncid,idv_bz,start_2d,count_2d,rdbz) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error return from nf_get_vara_real for var bz') do j=1,jmx do i=1,imx-1 bz(i,j) = rdbz(i,j) enddo ! i=1,imx-1 bz(imx,j) = bz(1,j) enddo ! j=1,jmx call fminmax(bz,imx*jmx,fmin,fmax) write(6,"('getmagnc: bz min,max=',2e12.4)") fmin,fmax ! ! Close netcdf file: istat = nf_close(ncid) end subroutine getmagnc end module mk_drifts