Cocal: /home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/printout_emfield

00001 subroutine printout_emfield_strength(iseq)
00002   use phys_constant,  only : long, Gauss, statV
00003   use grid_parameter, only : nrg, ntg, npg, ntgeq, nrf
00004   use def_matter_parameter, only : radi
00005   use def_faraday_tensor, only : fxd_grid, fyd_grid, fzd_grid, fijd_grid
00006   implicit none
00007   integer :: irg, itg, ipg, iseq, ii, jj
00008   integer :: E_max_gp(3,3), B_max_gp(3,3)
00009   integer :: Ep_max_gp(3), Bp_max_gp(3)
00010   integer :: Et_max_gp(3), Bt_max_gp(3)
00011   real(long) :: E_max(3),B_max(3)
00012   real(long) :: Ep_max, Et_max, Bp_max, Bt_max  ! Poloidal and Toroidal
00013   real(long) :: E_cens(4,3), B_cens(4,3)  ! 1 = Center,     2 = Equator
00014   real(long) :: Ep_cens(4), Et_cens(4)    ! 3 = North Pole, 4 = South Pole
00015   real(long) :: Bp_cens(4), Bt_cens(4)
00016   real(long) :: tmpE(3), tmpB(3), tmpEt, tmpEp, tmpBt, tmpBp
00017   character(len=31) char_31
00018 !
00019 ! --- Find maximum value
00020 !
00021   E_max(1:3) = 0.0d0 ; Ep_max = 0.0d0 ; Et_max = 0.0d0
00022   B_max(1:3) = 0.0d0 ; Bp_max = 0.0d0 ; Bt_max = 0.0d0
00023   ipg = 0
00024   do itg = 0, ntg
00025     do irg = 0, nrg
00026 !
00027       tmpE(1) = fxd_grid(irg,itg,ipg)/radi
00028       tmpE(2) = fyd_grid(irg,itg,ipg)/radi
00029       tmpE(3) = fzd_grid(irg,itg,ipg)/radi
00030       tmpB(1) =   fijd_grid(irg,itg,ipg,3)/radi
00031       tmpB(2) = - fijd_grid(irg,itg,ipg,2)/radi
00032       tmpB(3) =   fijd_grid(irg,itg,ipg,1)/radi
00033 !
00034       tmpEp = sqrt(tmpE(1)**2+tmpE(3)**2)
00035       tmpBp = sqrt(tmpB(1)**2+tmpB(3)**2)
00036       tmpEt = dabs(tmpE(2))
00037       tmpBt = dabs(tmpB(2))
00038 !
00039       do ii = 1, 3
00040         if (dabs(tmpE(ii)).ge.dabs(E_max(ii))) then 
00041           E_max(ii) = tmpE(ii)
00042           E_max_gp(ii,1) = irg ; E_max_gp(ii,2) = itg ; E_max_gp(ii,3) = ipg
00043         end if
00044         if (dabs(tmpB(ii)).ge.dabs(B_max(ii))) then
00045           B_max(ii) = tmpB(ii)
00046           B_max_gp(ii,1) = irg ; B_max_gp(ii,2) = itg ; B_max_gp(ii,3) = ipg
00047         end if
00048       end do
00049       if (tmpEp.ge.Ep_max) then
00050         Ep_max = tmpEp
00051         Ep_max_gp(1) = irg ; Ep_max_gp(2) = itg ; Ep_max_gp(3) = ipg
00052       end if
00053       if (tmpEt.ge.Et_max) then
00054         Et_max = tmpEt
00055         Et_max_gp(1) = irg ; Et_max_gp(2) = itg ; Et_max_gp(3) = ipg
00056       end if
00057       if (tmpBp.ge.Bp_max) then
00058         Bp_max = tmpBp
00059         Bp_max_gp(1) = irg ; Bp_max_gp(2) = itg ; Bp_max_gp(3) = ipg
00060       end if
00061       if (tmpBt.ge.Bt_max) then
00062         Bt_max = tmpBt
00063         Bt_max_gp(1) = irg ; Bt_max_gp(2) = itg ; Bt_max_gp(3) = ipg
00064       end if
00065     end do
00066   end do
00067 !
00068   do ii = 1, 4
00069     if (ii.eq.1) then ; irg = 0   ; itg = 0     ; ipg = 0 ; end if
00070     if (ii.eq.2) then ; irg = nrf ; itg = ntgeq ; ipg = 0 ; end if
00071     if (ii.eq.3) then ; irg = nrf ; itg = 0     ; ipg = 0 ; end if
00072     if (ii.eq.4) then ; irg = nrf ; itg = ntg   ; ipg = 0 ; end if
00073     E_cens(ii,1) = fxd_grid(irg,itg,ipg)/radi
00074     E_cens(ii,2) = fyd_grid(irg,itg,ipg)/radi
00075     E_cens(ii,3) = fzd_grid(irg,itg,ipg)/radi
00076     B_cens(ii,1) =   fijd_grid(irg,itg,ipg,3)/radi
00077     B_cens(ii,2) = - fijd_grid(irg,itg,ipg,2)/radi
00078     B_cens(ii,3) =   fijd_grid(irg,itg,ipg,1)/radi
00079     Ep_cens(ii) = sqrt(E_cens(ii,1)**2 + E_cens(ii,3)**2)
00080     Bp_cens(ii) = sqrt(B_cens(ii,1)**2 + B_cens(ii,3)**2)
00081     Et_cens(ii) = dabs(E_cens(ii,2))
00082     Bt_cens(ii) = dabs(B_cens(ii,2))
00083   end do
00084 !
00085   if (iseq.eq.1) then 
00086     open(190,file='rnsEMFseq.dat',status='unknown')
00087   end if
00088   write(190,*) '== Sequence number == ', iseq
00089   write(190,*) '## EM fields in G = c = M = 1/4pi epsi = 1 unit ##'
00090 !
00091   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Ex   max  =', E_max(1), &
00092   &             '(',E_max_gp(1,1),',',E_max_gp(1,2),',',E_max_gp(1,3),')'
00093   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Ey   max  =', E_max(2), &
00094   &             '(',E_max_gp(2,1),',',E_max_gp(2,2),',',E_max_gp(2,3),')'
00095   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Ez   max  =', E_max(3), &
00096   &             '(',E_max_gp(3,1),',',E_max_gp(3,2),',',E_max_gp(3,3),')'
00097   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Epol max  =', Ep_max, &
00098   &                   '(',Ep_max_gp(1),',',Ep_max_gp(2),',',Ep_max_gp(3),')'
00099   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Etor max  =', Et_max, &
00100   &                   '(',Et_max_gp(1),',',Et_max_gp(2),',',Et_max_gp(3),')'
00101   do ii = 1, 4
00102     if (ii.eq.1) then ; irg = 0   ; itg = 0     ; ipg = 0 ; end if
00103     if (ii.eq.2) then ; irg = nrf ; itg = ntgeq ; ipg = 0 ; end if
00104     if (ii.eq.3) then ; irg = nrf ; itg = 0     ; ipg = 0 ; end if
00105     if (ii.eq.4) then ; irg = nrf ; itg = ntg   ; ipg = 0 ; end if
00106     if (ii.eq.1) char_31 = ' -- Value at the center     -- '
00107     if (ii.eq.2) char_31 = ' -- Value at the equator    -- '
00108     if (ii.eq.3) char_31 = ' -- Value at the north pole -- '
00109     if (ii.eq.4) char_31 = ' -- Value at the south pole -- '
00110     write(190,'((a31,3x,3(a1,i3),a1))') char_31,'(',irg,',',itg,',',ipg,')'
00111     write(190,'(a12,1p,3e22.14)') '(Ex,Ey,Ez) =', (E_cens(ii,jj),jj=1,3)
00112     write(190,'(a12,1p,2e22.14)') '(Epol,Etor)=', Ep_cens(ii), Et_cens(ii)
00113   end do
00114 !
00115   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Bx   max  =', B_max(1), &
00116   &             '(',B_max_gp(1,1),',',B_max_gp(1,2),',',B_max_gp(1,3),')'
00117   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' By   max  =', B_max(2), &
00118   &             '(',B_max_gp(2,1),',',B_max_gp(2,2),',',B_max_gp(2,3),')'
00119   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Bz   max  =', B_max(3), &
00120   &             '(',B_max_gp(3,1),',',B_max_gp(3,2),',',B_max_gp(3,3),')'
00121   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Bpol max  =', Bp_max, &
00122   &                   '(',Bp_max_gp(1),',',Bp_max_gp(2),',',Bp_max_gp(3),')'
00123   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Btor max  =', Bt_max, &
00124   &                   '(',Bt_max_gp(1),',',Bt_max_gp(2),',',Bt_max_gp(3),')'
00125   do ii = 1, 4
00126     if (ii.eq.1) then ; irg = 0   ; itg = 0     ; ipg = 0 ; end if
00127     if (ii.eq.2) then ; irg = nrf ; itg = ntgeq ; ipg = 0 ; end if
00128     if (ii.eq.3) then ; irg = nrf ; itg = 0     ; ipg = 0 ; end if
00129     if (ii.eq.4) then ; irg = nrf ; itg = ntg   ; ipg = 0 ; end if
00130     if (ii.eq.1) char_31 = ' -- Value at the center     -- '
00131     if (ii.eq.2) char_31 = ' -- Value at the equator    -- '
00132     if (ii.eq.3) char_31 = ' -- Value at the north pole -- '
00133     if (ii.eq.4) char_31 = ' -- Value at the south pole -- '
00134     write(190,'((a31,3x,3(a1,i3),a1))') char_31,'(',irg,',',itg,',',ipg,')'
00135     write(190,'(a12,1p,3e22.14)') '(Bx,By,Bz) =', (B_cens(ii,jj),jj=1,3)
00136     write(190,'(a12,1p,2e22.14)') '(Bpol,Btor)=', Bp_cens(ii), Bt_cens(ii)
00137   end do
00138 !
00139   write(190,*) '## EM fields in Gauss ##' 
00140 !
00141   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Ex   max  =',E_max(1)*statV,&
00142   &             '(',E_max_gp(1,1),',',E_max_gp(1,2),',',E_max_gp(1,3),')'
00143   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Ey   max  =',E_max(2)*statV,&
00144   &             '(',E_max_gp(2,1),',',E_max_gp(2,2),',',E_max_gp(2,3),')'
00145   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Ez   max  =',E_max(3)*statV,&
00146   &             '(',E_max_gp(3,1),',',E_max_gp(3,2),',',E_max_gp(3,3),')'
00147   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Epol max  =',Ep_max*statV,&
00148   &                   '(',Ep_max_gp(1),',',Ep_max_gp(2),',',Ep_max_gp(3),')'
00149   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Etor max  =',Et_max*statV,&
00150   &                   '(',Et_max_gp(1),',',Et_max_gp(2),',',Et_max_gp(3),')'
00151   do ii = 1, 4
00152     if (ii.eq.1) then ; irg = 0   ; itg = 0     ; ipg = 0 ; end if
00153     if (ii.eq.2) then ; irg = nrf ; itg = ntgeq ; ipg = 0 ; end if
00154     if (ii.eq.3) then ; irg = nrf ; itg = 0     ; ipg = 0 ; end if
00155     if (ii.eq.4) then ; irg = nrf ; itg = ntg   ; ipg = 0 ; end if
00156     if (ii.eq.1) char_31 = ' -- Value at the center     -- '
00157     if (ii.eq.2) char_31 = ' -- Value at the equator    -- '
00158     if (ii.eq.3) char_31 = ' -- Value at the north pole -- '
00159     if (ii.eq.4) char_31 = ' -- Value at the south pole -- '
00160     write(190,'((a31,3x,3(a1,i3),a1))') char_31,'(',irg,',',itg,',',ipg,')'
00161     write(190,'(a12,1p,3e22.14)') '(Ex,Ey,Ez) =', (E_cens(ii,jj)*statV,jj=1,3)
00162     write(190,'(a12,1p,2e22.14)') '(Epol,Etor)=', Ep_cens(ii)*statV, &
00163     &                                             Et_cens(ii)*statV
00164   end do
00165 !
00166   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Bx   max  =',B_max(1)*Gauss,&
00167   &             '(',B_max_gp(1,1),',',B_max_gp(1,2),',',B_max_gp(1,3),')'
00168   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' By   max  =',B_max(2)*Gauss,&
00169   &             '(',B_max_gp(2,1),',',B_max_gp(2,2),',',B_max_gp(2,3),')'
00170   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Bz   max  =',B_max(3)*Gauss,&
00171   &             '(',B_max_gp(3,1),',',B_max_gp(3,2),',',B_max_gp(3,3),')'
00172   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Bpol max  =',Bp_max*Gauss,&
00173   &                   '(',Bp_max_gp(1),',',Bp_max_gp(2),',',Bp_max_gp(3),')'
00174   write(190,'(a12,1p,1e22.14,3x,3(a1,i3),a1)') ' Btor max  =',Bt_max*Gauss,&
00175   &                   '(',Bt_max_gp(1),',',Bt_max_gp(2),',',Bt_max_gp(3),')'
00176   do ii = 1, 4
00177     if (ii.eq.1) then ; irg = 0   ; itg = 0     ; ipg = 0 ; end if
00178     if (ii.eq.2) then ; irg = nrf ; itg = ntgeq ; ipg = 0 ; end if
00179     if (ii.eq.3) then ; irg = nrf ; itg = 0     ; ipg = 0 ; end if
00180     if (ii.eq.4) then ; irg = nrf ; itg = ntg   ; ipg = 0 ; end if
00181     if (ii.eq.1) char_31 = ' -- Value at the center     -- '
00182     if (ii.eq.2) char_31 = ' -- Value at the equator    -- '
00183     if (ii.eq.3) char_31 = ' -- Value at the north pole -- '
00184     if (ii.eq.4) char_31 = ' -- Value at the south pole -- '
00185     write(190,'((a31,3x,3(a1,i3),a1))') char_31,'(',irg,',',itg,',',ipg,')'
00186     write(190,'(a12,1p,3e22.14)') '(Bx,By,Bz) =', (B_cens(ii,jj)*Gauss,jj=1,3)
00187     write(190,'(a12,1p,2e22.14)') '(Bpol,Btor)=', Bp_cens(ii)*Gauss, &
00188     &                                             Bt_cens(ii)*Gauss
00189   end do
00190 !
00191 end subroutine printout_emfield_strength
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/printout_emfield_strength.f90
