00001 subroutine interpo_patch_to_active_patch_mpt(fnc,cfn,rc,thc,phic)
00002 use phys_constant, only : long, pi
00003 use grid_parameter_interpo, only : nrg_itp, ntg_itp, npg_itp
00004 use coordinate_grav_extended_interpo
00005 implicit none
00006 real(long), pointer :: fnc(:,:,:)
00007 real(long), intent(out) :: cfn
00008 real(long), intent(in) :: rc, thc, phic
00009 real(long) :: r4(4), th4(4), phi4(4), fr4(4), ft4(4), fp4(4)
00010 integer :: irg, itg, ipg, irgex, itgex, ipgex
00011 integer :: ir0, it0, ip0, irg0 , itg0 , ipg0, ii, jj, kk
00012 real(long), external :: lagint_4th
00013
00014
00015
00016
00017 cfn = 0.0d0
00018
00019 do irg = 0, nrg_itp+1
00020 if (rc.lt.rgex_itp(irg).and.rc.ge.rgex_itp(irg-1)) &
00021 & ir0 = min0(irg-2,nrg_itp-3)
00022 end do
00023 do itg = 0, ntg_itp+1
00024 if (thc.lt.thgex_itp(itg).and.thc.ge.thgex_itp(itg-1)) it0 = itg-2
00025 end do
00026 do ipg = 0, npg_itp+1
00027 if (phic.lt.phigex_itp(ipg).and.phic.ge.phigex_itp(ipg-1)) ip0 = ipg-2
00028 end do
00029
00030 do ii = 1, 4
00031 irg0 = ir0 + ii - 1
00032 itg0 = it0 + ii - 1
00033 ipg0 = ip0 + ii - 1
00034 r4(ii) = rgex_itp(irg0)
00035 th4(ii) = thgex_itp(itg0)
00036 phi4(ii) = phigex_itp(ipg0)
00037 end do
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047 do kk = 1, 4
00048 ipg0 = ip0 + kk - 1
00049 do jj = 1, 4
00050 itg0 = it0 + jj - 1
00051 do ii = 1, 4
00052 irg0 = ir0 + ii - 1
00053 irgex = irgex_r_itp(irg0)
00054 itgex = itgex_r_itp(itgex_th_itp(itg0),irg0)
00055 ipgex = ipgex_r_itp(ipgex_th_itp(ipgex_phi_itp(ipg0),itg0),irg0)
00056 fr4(ii) = fnc(irgex,itgex,ipgex)
00057 end do
00058 ft4(jj) = lagint_4th(r4,fr4,rc)
00059 end do
00060 fp4(kk) = lagint_4th(th4,ft4,thc)
00061 end do
00062 cfn = lagint_4th(phi4,fp4,phic)
00063
00064 end subroutine interpo_patch_to_active_patch_mpt