00001 subroutine sourceterm_MWtemp_CF(sou)
00002   use phys_constant, only : long
00003   use grid_parameter, only : nrg, ntg, npg
00004   use def_metric, only : alph, psi
00005   use def_faraday_tensor, only : fxd, fyd, fzd
00006   use def_emfield_derivatives, only : Lie_bAxd_grid, Lie_bAyd_grid, &
00007   &                                   Lie_bAzd_grid
00008   use interface_interpo_linear_type0
00009   use interface_grgrad_midpoint
00010   use make_array_3d
00011   implicit none
00012   real(long), pointer :: sou(:,:,:)
00013   real(long), pointer :: fnc2(:,:,:)
00014   real(long), pointer :: grad2x(:,:,:), grad2y(:,:,:), grad2z(:,:,:)
00015   real(long), pointer :: dxLbAx(:,:,:), dyLbAx(:,:,:), dzLbAx(:,:,:)
00016   real(long), pointer :: dxLbAy(:,:,:), dyLbAy(:,:,:), dzLbAy(:,:,:)
00017   real(long), pointer :: dxLbAz(:,:,:), dyLbAz(:,:,:), dzLbAz(:,:,:)
00018   real(long) :: psigc, alphgc, a2divp2, dxLbAxgc, dyLbAygc, dzLbAzgc
00019   real(long) :: fxdgc, fydgc, fzdgc
00020   real(long) :: dxps2al, dyps2al, dzps2al, dps2alfa
00021   integer    :: irg, itg, ipg
00022 
00023 
00024 
00025 
00026   call alloc_array3d(fnc2,0,nrg,0,ntg,0,npg)
00027   call alloc_array3d(grad2x,1,nrg,1,ntg,1,npg)
00028   call alloc_array3d(grad2y,1,nrg,1,ntg,1,npg)
00029   call alloc_array3d(grad2z,1,nrg,1,ntg,1,npg)
00030   call alloc_array3d(dxLbAx,1,nrg,1,ntg,1,npg)
00031   call alloc_array3d(dyLbAx,1,nrg,1,ntg,1,npg)
00032   call alloc_array3d(dzLbAx,1,nrg,1,ntg,1,npg)
00033   call alloc_array3d(dxLbAy,1,nrg,1,ntg,1,npg)
00034   call alloc_array3d(dyLbAy,1,nrg,1,ntg,1,npg)
00035   call alloc_array3d(dzLbAy,1,nrg,1,ntg,1,npg)
00036   call alloc_array3d(dxLbAz,1,nrg,1,ntg,1,npg)
00037   call alloc_array3d(dyLbAz,1,nrg,1,ntg,1,npg)
00038   call alloc_array3d(dzLbAz,1,nrg,1,ntg,1,npg)
00039   fnc2(0:nrg,0:ntg,0:npg) = psi(0:nrg,0:ntg,0:npg)**2/alph(0:nrg,0:ntg,0:npg)
00040   call grgrad_midpoint(fnc2,grad2x,grad2y,grad2z)
00041   call grgrad_midpoint(Lie_bAxd_grid,dxLbAx,dyLbAx,dzLbAx)
00042   call grgrad_midpoint(Lie_bAyd_grid,dxLbAy,dyLbAy,dzLbAy)
00043   call grgrad_midpoint(Lie_bAzd_grid,dxLbAz,dyLbAz,dzLbAz)
00044 
00045   do ipg = 1, npg
00046     do itg = 1, ntg
00047       do irg = 1, nrg
00048         call interpo_linear_type0(psigc,psi,irg,itg,ipg)
00049         call interpo_linear_type0(alphgc,alph,irg,itg,ipg)
00050         a2divp2 = alphgc**2/psigc**2
00051         dxps2al = grad2x(irg,itg,ipg)
00052         dyps2al = grad2y(irg,itg,ipg)
00053         dzps2al = grad2z(irg,itg,ipg)
00054         fxdgc = fxd(irg,itg,ipg)
00055         fydgc = fyd(irg,itg,ipg)
00056         fzdgc = fzd(irg,itg,ipg)
00057         dxLbAxgc = dxLbAx(irg,itg,ipg)
00058         dyLbAygc = dyLbAy(irg,itg,ipg)
00059         dzLbAzgc = dzLbAz(irg,itg,ipg)
00060 
00061         dps2alfa = a2divp2*(dxps2al*fxdgc+dyps2al*fydgc+dzps2al*fzdgc)
00062 
00063         sou(irg,itg,ipg) = dps2alfa + dxLbAxgc + dyLbAygc + dzLbAzgc
00064 
00065       end do
00066     end do
00067   end do
00068 
00069   deallocate(fnc2)
00070   deallocate(grad2x)
00071   deallocate(grad2y)
00072   deallocate(grad2z)
00073   deallocate(dxLbAx)
00074   deallocate(dyLbAx)
00075   deallocate(dzLbAx)
00076   deallocate(dxLbAy)
00077   deallocate(dyLbAy)
00078   deallocate(dzLbAy)
00079   deallocate(dxLbAz)
00080   deallocate(dyLbAz)
00081   deallocate(dzLbAz)
00082 
00083 end subroutine sourceterm_MWtemp_CF
00084