! $Id$ module tr_bpsd use bpsd type(bpsd_device_type), private,save :: device type(bpsd_species_type), private,save :: species type(bpsd_equ1D_type), private,save :: equ1D type(bpsd_metric1D_type),private,save :: metric1D type(bpsd_plasmaf_type), private,save :: plasmaf logical, private, save :: tr_bpsd_init_flag = .TRUE. public contains !======================================================================= subroutine tr_bpsd_init !======================================================================= use trcomm ! local variables integer(4) :: ns,nr,ierr real(8) :: temp(nrmp,nsm,3) !======================================================================= if(tr_bpsd_init_flag) then species%nsmax=0 plasmaf%nsmax=0 plasmaf%nrmax=0 tr_bpsd_init_flag=.FALSE. endif ! write(6,*) 'top of tr_bpsd_init:qp' ! write(6,'(1P5E12.4)') (qp(nr),nr=1,nrmax) ! write(6,*) 'top of tr_bpsd_init:rt' ! write(6,'(1P5E12.4)') (rt(nr,1),nr=1,nrmax) ! pause device%rr=RR device%zz=0.d0 device%ra=RA device%rb=RA+0.1D0 device%bb=BB device%ip=RIP device%elip=RKAP device%trig=RDLT call bpsd_set_data(device,ierr) if(species%nsmax.ne.nsmax) then if(associated(species%data)) then deallocate(species%data) endif species%nsmax=nsmax allocate(species%data(species%nsmax)) endif do ns=1,species%nsmax species%data(ns)%pa=pa(ns) species%data(ns)%pz=pz(ns) species%data(ns)%pz0=pz(ns) enddo call bpsd_set_data(species,ierr) if((equ1D%nrmax.ne.nrmax+1)) then if(associated(equ1D%rho)) then deallocate(equ1D%rho) endif if(associated(equ1D%data)) then deallocate(equ1D%data) endif equ1D%nrmax=nrmax+1 allocate(equ1D%rho(equ1D%nrmax)) allocate(equ1D%data(equ1D%nrmax)) endif equ1D%time=0.d0 equ1D%rho(1)=0.d0 do nr=1,nrmax equ1D%rho(nr+1)=rg(nr) enddo if((metric1D%nrmax.ne.nrmax+1)) then if(associated(metric1D%rho)) then deallocate(metric1D%rho) endif if(associated(metric1D%data)) then deallocate(metric1D%data) endif metric1D%nrmax=nrmax+1 allocate(metric1D%rho(metric1D%nrmax)) allocate(metric1D%data(metric1D%nrmax)) endif metric1D%time=0.d0 metric1D%rho(1)=0.d0 do nr=1,nrmax metric1D%rho(nr+1)=rg(nr) enddo if((plasmaf%nsmax.ne.nsmax).or. & & (plasmaf%nrmax.ne.nrmax+1)) then if(associated(plasmaf%rho)) then deallocate(plasmaf%rho) endif if(associated(plasmaf%data)) then deallocate(plasmaf%data) endif if(associated(plasmaf%qinv)) then deallocate(plasmaf%qinv) endif plasmaf%nsmax=nsmax plasmaf%nrmax=nrmax+1 allocate(plasmaf%rho(plasmaf%nrmax)) allocate(plasmaf%data(plasmaf%nrmax,plasmaf%nsmax)) allocate(plasmaf%qinv(plasmaf%nrmax)) endif plasmaf%time=0.d0 plasmaf%rho(1)=0.d0 do nr=1,nrmax plasmaf%rho(nr+1)=rg(nr) enddo do ns=1,nsmax call mesh_convert_mtog(rn(1:nrmax,ns),temp(1:nrmax,ns,1),nrmax) call mesh_convert_mtog(rt(1:nrmax,ns),temp(1:nrmax,ns,2),nrmax) call mesh_convert_mtog(ru(1:nrmax,ns),temp(1:nrmax,ns,3),nrmax) enddo do nr=1,plasmaf%nrmax do ns=1,plasmaf%nsmax plasmaf%data(nr,ns)%pn=temp(nr,ns,1)*1.d20 plasmaf%data(nr,ns)%pt=temp(nr,ns,2)*1.D3 plasmaf%data(nr,ns)%ptpr=temp(nr,ns,2)*1.D3 plasmaf%data(nr,ns)%ptpp=temp(nr,ns,2)*1.D3 plasmaf%data(nr,ns)%pu=temp(nr,ns,3) enddo enddo !! QP(1:NRMAX)=TTRHOG(1:NRMAX)*ARRHOG(1:NRMAX)*DVRHOG(1:NRMAX) & !! & /(4.D0*PI**2*RDP(1:NRMAX)) QP(1:NRMAX)=TTRHOG(1:NRMAX)*ARRHOG(1:NRMAX)/(4.D0*PI**2*RDPVRHOG(1:NRMAX)) do nr=2,plasmaf%nrmax plasmaf%qinv(nr)=1.d0/qp(nr-1) enddo plasmaf%qinv(1)=((plasmaf%rho(3))**2*plasmaf%qinv(2) & & -(plasmaf%rho(2))**2*plasmaf%qinv(3)) & & /((plasmaf%rho(3))**2-(plasmaf%rho(2))**2) call bpsd_set_data(plasmaf,ierr) return end subroutine tr_bpsd_init !======================================================================= subroutine tr_bpsd_set(ierr) !======================================================================= use trcomm integer(4) :: ierr ! local variables integer(4) :: ns,nr real(8) :: temp(nrmp,nsm,3) !======================================================================= ! write(6,*) 'top of tr_bpsd_set:qp' ! write(6,'(1P5E12.4)') (qp(nr),nr=1,nrmax) ! write(6,*) 'top of tr_bpsd_set: rt' ! write(6,'(1P5E12.4)') (rt(nr,1),nr=1,nrmax) ! pause device%rr=RR device%zz=0.d0 device%ra=RA device%rb=RA+0.1D0 device%bb=BB device%ip=RIP device%elip=RKAP device%trig=RDLT call bpsd_set_data(device,ierr) plasmaf%time=t do ns=1,nsmax call mesh_convert_mtog(rn(1:nrmax,ns),temp(1:nrmax,ns,1),nrmax) call mesh_convert_mtog(rt(1:nrmax,ns),temp(1:nrmax,ns,2),nrmax) call mesh_convert_mtog(ru(1:nrmax,ns),temp(1:nrmax,ns,3),nrmax) enddo do nr=1,plasmaf%nrmax if(nr.eq.1) then plasmaf%rho(nr)=0.d0 else plasmaf%rho(nr)=rg(nr-1) endif do ns=1,plasmaf%nsmax plasmaf%data(nr,ns)%pn=temp(nr,ns,1)*1.d20 plasmaf%data(nr,ns)%pt=temp(nr,ns,2)*1.D3 plasmaf%data(nr,ns)%ptpr=temp(nr,ns,2)*1.D3 plasmaf%data(nr,ns)%ptpp=temp(nr,ns,2)*1.D3 plasmaf%data(nr,ns)%pu=temp(nr,ns,3) enddo enddo ! write(6,*) 'RDP_set:',(rdp(nr),nr=1,nrmax) do nr=1,nrmax !! plasmaf%qinv(nr+1) & !! & =(4.D0*PI**2*RDP(nr))/(TTRHOG(nr)*ARRHOG(nr)*DVRHOG(nr)) plasmaf%qinv(nr+1)=(4.D0*PI**2*RDPVRHOG(nr))/(TTRHOG(nr)*ARRHOG(nr)) enddo plasmaf%qinv(1)=((plasmaf%rho(3))**2*plasmaf%qinv(2) & & -(plasmaf%rho(2))**2*plasmaf%qinv(3)) & & /((plasmaf%rho(3))**2-(plasmaf%rho(2))**2) ! write(6,*) 'end of tr_bpsd_set:qp' ! write(6,'(1P5E12.4)') (qp(nr),nr=1,nrmax) ! write(6,*) 'end of tr_bpsd_set:rt' ! write(6,'(1P5E12.4)') (rt(nr,1),nr=1,nrmax) ! pause call bpsd_set_data(plasmaf,ierr) return end subroutine tr_bpsd_set !======================================================================= subroutine tr_bpsd_get(ierr) !======================================================================= use trcomm integer(4),intent(out) :: ierr ! local variables integer(4) :: ns,nr real(8) :: temp(nrmp,nsm,3) real(8) :: tempx(nrmp,17),psita,dpsitdrho,dvdrho,rgl REAL(8) :: FACTOR0, FACTORM, FACTORP !======================================================================= ! write(6,*) 'top of tr_bpsd_get: qp' ! write(6,'(1P5E12.4)') (qp(nr),nr=1,nrmax) ! write(6,*) 'top of tr_bpsd_get: rt' ! write(6,'(1P5E12.4)') (rt(nr,1),nr=1,nrmax) ! pause call bpsd_get_data(device,ierr) RR=device%rr RA=device%ra BB=device%bb RIP=device%ip RKAP=device%elip RDLT=device%trig if(modelg.eq.3.or.modelg.eq.5) then RIPS=RIP RIPE=RIP endif call bpsd_get_data(plasmaf,ierr) do ns=1,plasmaf%nsmax do nr=1,plasmaf%nrmax temp(nr,ns,1)=plasmaf%data(nr,ns)%pn*1.d-20 temp(nr,ns,2)=plasmaf%data(nr,ns)%pt*1.D-3 temp(nr,ns,3)=plasmaf%data(nr,ns)%pu enddo enddo do nr=2,plasmaf%nrmax qp(nr-1)=1.d0/plasmaf%qinv(nr) enddo Q0=2.d0*qp(1)-qp(2) do ns=1,nsmax call mesh_convert_gtom(temp(1:nrmax,ns,1),rn(1:nrmax,ns),nrmax) call mesh_convert_gtom(temp(1:nrmax,ns,2),rt(1:nrmax,ns),nrmax) call mesh_convert_gtom(temp(1:nrmax,ns,3),ru(1:nrmax,ns),nrmax) enddo if(modelg.eq.3.or.modelg.eq.5.or.modelg.eq.8.or.modelg.eq.9) then equ1D%nrmax=nrmax+1 equ1D%rho(1)=0.d0 do nr=2,nrmax+1 equ1D%rho(nr)=rg(nr-1) enddo call bpsd_get_data(equ1D,ierr) do nr=1,equ1D%nrmax tempx(nr,1)=equ1D%data(nr)%psit tempx(nr,2)=equ1D%data(nr)%psip tempx(nr,3)=equ1D%data(nr)%ppp tempx(nr,4)=equ1D%data(nr)%piq tempx(nr,5)=equ1D%data(nr)%pip*rmu0/(2.d0*pi) tempx(nr,6)=equ1D%data(nr)%pit enddo call data_interpolate_gtom_full(tempx(1,1),PSITRHO,nrmax) call data_interpolate_gtom_full(tempx(1,2),PSIPRHO,nrmax) call data_interpolate_gtom_full(tempx(1,3),PPPRHO,nrmax) call data_interpolate_gtom_full(tempx(1,4),PIQRHO,nrmax) call data_interpolate_gtom_full(tempx(1,5),TTRHO,nrmax) call data_interpolate_gtom_full(tempx(1,6),PIRHO,nrmax) TTRHOG(1:nrmax)=tempx(2:nrmax+1,5) psita=equ1D%data(equ1D%nrmax)%psit metric1D%nrmax=nrmax+1 metric1D%rho(1)=0.d0 do nr=2,nrmax+1 metric1D%rho(nr)=rg(nr-1) enddo call bpsd_get_data(metric1D,ierr) do nr=2,metric1D%nrmax ! metric1D%nrmax = nrmax + 1 rgl=rg(nr-1) ! rgl is equivalent to metric1D%rho(nr) dpsitdrho=2.D0*psita*rgl dvdrho=metric1D%data(nr)%dvpsit*dpsitdrho tempx(nr,1)=dvdrho tempx(nr,2)=metric1D%data(nr)%avegrr2 tempx(nr,3)=metric1D%data(nr)%aver2i tempx(nr,4)=metric1D%data(nr)%avegr tempx(nr,5)=metric1D%data(nr)%avegr2 tempx(nr,6)=metric1D%data(nr)%aveb2 tempx(nr,7)=metric1D%data(nr)%aveb2i tempx(nr,8)=metric1D%data(nr)%avegvr2/dvdrho**2 tempx(nr,9)=metric1D%data(nr)%rs/metric1D%data(nr)%rr tempx(nr,10)=metric1D%data(nr)%rr tempx(nr,11)=metric1D%data(nr)%rs tempx(nr,12)=metric1D%data(nr)%elip tempx(nr,13)=1.d0/metric1D%data(nr)%dvpsip/(2.d0*pi) tempx(nr,14)=metric1D%data(nr)%avegvr2 tempx(nr,15)=metric1D%data(nr)%pvol tempx(nr,16)=metric1D%data(nr)%psur tempx(nr,17)=metric1D%data(nr)%aveb enddo nr=1 ! equivalent to "rho = 0" tempx(nr,1)=0.d0 ! definition tempx(nr,2)=metric1D%data(nr)%avegrr2 tempx(nr,3)=metric1D%data(nr)%aver2i tempx(nr,4)=metric1D%data(nr)%avegr tempx(nr,5)=metric1D%data(nr)%avegr2 tempx(nr,6)=metric1D%data(nr)%aveb2 tempx(nr,7)=metric1D%data(nr)%aveb2i tempx(nr,8)=tempx(2,8) ! dummy because dvdrho=0. tempx(nr,9)=metric1D%data(nr)%rs/metric1D%data(nr)%rr tempx(nr,10)=metric1D%data(nr)%rr tempx(nr,11)=metric1D%data(nr)%rs tempx(nr,12)=metric1D%data(nr)%elip tempx(nr,13)=1.d0/metric1D%data(nr)%dvpsip/(2.d0*pi) tempx(nr,14)=metric1D%data(nr)%avegvr2 tempx(nr,15)=metric1D%data(nr)%pvol tempx(nr,16)=metric1D%data(nr)%psur tempx(nr,17)=metric1D%data(nr)%aveb call data_interpolate_gtom_full(tempx(1,1), DVRHO, nrmax) call data_interpolate_gtom_full(tempx(1,2), ABRHO, nrmax) call data_interpolate_gtom_full(tempx(1,3), ARRHO, nrmax) call data_interpolate_gtom_full(tempx(1,4), AR1RHO,nrmax) call data_interpolate_gtom_full(tempx(1,5), AR2RHO,nrmax) call data_interpolate_gtom_full(tempx(1,12),RKPRHO,nrmax) call data_interpolate_gtom_full(tempx(1,14),ABVRHO,nrmax) call data_interpolate_gtom_full(tempx(1,17),ABB1RHO,nrmax) DVRHOG (1:nrmax)=tempx(2:nrmax+1,1) ! dV/drho ABRHOG (1:nrmax)=tempx(2:nrmax+1,2) ! <|grad rho|^2/R^2>; avegrr2 ARRHOG (1:nrmax)=tempx(2:nrmax+1,3) ! <1/R^2>; aver2i AR1RHOG (1:nrmax)=tempx(2:nrmax+1,4) ! <|grad rho|> AR2RHOG (1:nrmax)=tempx(2:nrmax+1,5) ! <|grad rho|^2> ABB2RHOG(1:nrmax)=tempx(2:nrmax+1,6) ! AIB2RHOG(1:nrmax)=tempx(2:nrmax+1,7) ! <1/B^2> ARHBRHOG(1:nrmax)=tempx(2:nrmax+1,8) ! not used EPSRHO (1:nrmax)=tempx(2:nrmax+1,9) RMJRHO (1:nrmax)=tempx(2:nrmax+1,10) ! local R RMNRHO (1:nrmax)=tempx(2:nrmax+1,11) ! local r RKPRHOG (1:nrmax)=tempx(2:nrmax+1,12) ! local kappa RDPVRHOG(1:nrmax)=tempx(2:nrmax+1,13) ! dpsi/dV ABVRHOG (1:nrmax)=tempx(2:nrmax+1,14) ! <|grad V|^2/R^2> PVOLRHOG(1:nrmax)=tempx(2:nrmax+1,15) ! Plasma volume PSURRHOG(1:nrmax)=tempx(2:nrmax+1,16) ! Plasma surface do nr=1,nrmax ! RDP(nr)=TTRHOG(nr)*ARRHOG(nr)*DVRHOG(nr)/(4.D0*PI**2*QP(nr)) RDP(nr)=DVRHOG(nr)*RDPVRHOG(nr) enddo ! BP(1:NRMAX) =AR1RHOG(1:NRMAX)*RDP(1:NRMAX)/RR ! write(6,*) 'RDP_get:',(rdp(nr),nr=1,nrmax) NR=1 FACTOR0=TTRHO(NR)**2/(RMU0*BB*DVRHO(NR)) FACTORP=ABVRHOG(NR )/TTRHOG(NR ) AJ(NR) =FACTOR0*FACTORP*RDPVRHOG(NR)/DR AJOH(NR)=AJ(NR) DO NR=2,NRMAX FACTOR0=TTRHO(NR)**2/(RMU0*BB*DVRHO(NR)) FACTORM=ABVRHOG(NR-1)/TTRHOG(NR-1) FACTORP=ABVRHOG(NR )/TTRHOG(NR ) AJ(NR) =FACTOR0*(FACTORP*RDPVRHOG(NR)-FACTORM*RDPVRHOG(NR-1))/DR AJOH(NR)=AJ(NR) ENDDO ! write(6,*) 'end of tr_bpsd_get: aj' ! write(6,'(1P5E12.4)') (aj(nr),nr=1,nrmax) ! write(6,*) 'end of tr_bpsd_get: qp' ! write(6,'(1P5E12.4)') (qp(nr),nr=1,nrmax) ! write(6,*) 'end of tr_bpsd_get: rt' ! write(6,'(1P5E12.4)') (rt(nr,1),nr=1,nrmax) ! pause ! Calibration in order to keep consistency between metrics and plasma current if(modelg.eq.3.or.modelg.eq.5) then RIP = ABVRHOG(NRMAX)*RDPVRHOG(NRMAX)/(2.D0*PI*RMU0)*1.D-6 RIPS = RIP RIPE = RIP endif endif return end subroutine tr_bpsd_get ! ----- convert half mesh to origin + grid mesh ----- ! ! Suppose that rho derivative of data be zero ! at the axis and the boundary ! datag(1) at rho = 0 and datag(nrmax+1) at rho = 1 subroutine mesh_convert_mtog(datam,datag,nrmax) implicit none integer(4), intent(in) :: nrmax real(8), intent(in) :: datam(nrmax) real(8), intent(out) :: datag(nrmax+1) datag(1) = (9.d0*datam(1)-datam(2))/8.d0 datag(2:nrmax) = 0.5d0 * (datam(1:nrmax-1) + datam(2:nrmax)) datag(nrmax+1) = (4.d0*datam(nrmax)-datam(nrmax-1))/3.d0 return end subroutine mesh_convert_mtog ! ----- convert origin + grid mesh to half mesh ----- ! ! just invert mesh_convert_gtom !!! CAUTION !!!=======================================================! ! This routine should be used only in case that one reconverts ! ! data converted by "mesh_convert_mtog" routine. ! ! In any other case, one should use "data_interpolate_gtom" routine. ! !=====================================================================! subroutine mesh_convert_gtom(datag,datam,nrmax) implicit none integer(4), intent(in) :: nrmax real(8), intent(in) :: datag(nrmax+1) real(8), intent(out) :: datam(nrmax) real(8) :: c11=9.d0/8.d0,c12=-1.d0/8.d0,c21=0.5d0,c22=0.5d0 real(8) :: det,a11,a12,a21,a22 det=c11*c22-c12*c21 a11= c22/det a12=-c12/det a21=-c21/det a22= c11/det datam(1)=a11*datag(1)+a12*datag(2) datam(2)=a21*datag(1)+a22*datag(2) datam(3:nrmax) = 2.d0 * datag(3:nrmax) - datam(2:nrmax-1) return end subroutine mesh_convert_gtom ! ----- convert half mesh to origin + grid mesh ----- ! ! Suppose that data be zero at the axis ! and rho derivative of data be zero at the boundary subroutine mesh_convert_mtog0(datam,datag,nrmax) implicit none integer(4), intent(in) :: nrmax real(8), intent(in) :: datam(nrmax) real(8), intent(out) :: datag(nrmax+1) datag(1) = 0.d0 datag(2:nrmax) = 0.5d0 * (datam(1:nrmax-1) + datam(2:nrmax)) datag(nrmax+1) = (4.d0*datam(nrmax)-datam(nrmax-1))/3.d0 return end subroutine mesh_convert_mtog0 ! ----- convert origin + grid mesh to half mesh ----- ! ! just invert mesh_convert_mtog0 !!! CAUTION !!!========================================================! ! This routine should be used only in case that one reconverts ! ! data converted by "mesh_convert_mtog" routine. ! ! In any other case, one should use "data_interpolate_gtom0" routine. ! !======================================================================! subroutine mesh_convert_gtom0(datag,datam,nrmax) implicit none integer(4), intent(in) :: nrmax real(8), intent(in) :: datag(nrmax+1) real(8), intent(out) :: datam(nrmax) real(8) :: c11=9.d0/8.d0,c12=-1.d0/8.d0,c21=0.5d0,c22=0.5d0 real(8) :: det,a11,a12,a21,a22 det=c11*c22-c12*c21 a11= c22/det a12=-c12/det a21=-c21/det a22= c11/det datam(1)=0.5d0*datag(2) datam(2:nrmax) = 2.d0 * datag(2:nrmax) - datam(1:nrmax-1) return end subroutine mesh_convert_gtom0 ! ----- interpolate data on half mesh from full data on grid ----- subroutine data_interpolate_gtom_full(datag,datam,nrmax) implicit none integer(4), intent(in) :: nrmax real(8), intent(in) :: datag(nrmax+1) real(8), intent(out) :: datam(nrmax) datam(1:nrmax) = 0.5d0 * (datag(1:nrmax) + datag(2:nrmax+1)) return end subroutine data_interpolate_gtom_full ! ----- interpolate data on half mesh ! from data on grid except the axis ----- subroutine data_interpolate_gtom(datag,datam,nrmax) implicit none integer(4), intent(in) :: nrmax real(8), intent(in) :: datag(nrmax+1) real(8), intent(out) :: datam(nrmax) ! linear extrapolation datam(1) = 1.5d0 * datag(2) - 0.5d0 * datag(3) datam(2:nrmax) = 0.5d0 * (datag(2:nrmax) + datag(3:nrmax+1)) return end subroutine data_interpolate_gtom end module tr_bpsd