/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/calc_surface.f90
Go to the documentation of this file.00001 subroutine calc_surface(rsnew,emd) 00002 use phys_constant, only : long 00003 use grid_parameter, only : nrf, ntf, npf 00004 use coordinate_grav_r, only : drg, rg 00005 use def_matter, only : rs 00006 ! use def_matter_parameter, only : emdc 00007 implicit none 00008 real(long), pointer :: emd(:,:,:) 00009 real(long), pointer :: rsnew(:,:) 00010 real(long) :: delta 00011 integer :: itf, ipf 00012 !quad real(long) :: sgg(1:4), gg(1:4,1:4) 00013 !quad real(long) :: r1,r2,r3,q1,q2,q3, aa,bb,cc,det,rs1,rs2,rs2a,rs2b 00014 ! 00015 do ipf = 0, npf 00016 do itf = 0, ntf 00017 delta = drg(nrf)*(0.0d0 - emd(nrf ,itf,ipf)) & 00018 & /(emd(nrf,itf,ipf) - emd(nrf-1,itf,ipf)) 00019 !robust delta = (drg(nrf)+drg(nrf-1)) & 00020 !robust & *(0.0d0 - emd(nrf ,itf,ipf)) & 00021 !robust & /(emd(nrf,itf,ipf) - emd(nrf-2,itf,ipf)) 00022 rsnew(itf,ipf) = rs(itf,ipf)*(1.0d0 + delta) 00023 ! 00024 ! =============== 3 point interpolation ============================ 00025 !quad q1 = emd(nrf ,itf,ipf); r1 = rs(itf,ipf)*rg(nrf ) 00026 !quad q2 = emd(nrf-1,itf,ipf); r2 = rs(itf,ipf)*rg(nrf-1) 00027 !quad q3 = emd(nrf-2,itf,ipf); r3 = rs(itf,ipf)*rg(nrf-2) 00028 !quad 00029 !quad sgg(1) = q1 00030 !quad sgg(2) = q2 00031 !quad sgg(3) = q3 00032 !quad 00033 !quad gg(1,1) = r1*r1 00034 !quad gg(1,2) = r1 00035 !quad gg(1,3) = 1.0d0 00036 !quad 00037 !quad gg(2,1) = r2*r2 00038 !quad gg(2,2) = r2 00039 !quad gg(2,3) = 1.0d0 00040 !quad 00041 !quad gg(3,1) = r3*r3 00042 !quad gg(3,2) = r3 00043 !quad gg(3,3) = 1.0d0 00044 !quad 00045 !quad call minv(gg,sgg,3,4) 00046 !quad aa = sgg(1) 00047 !quad bb = sgg(2) 00048 !quad cc = sgg(3) 00049 !quad det = bb*bb-4.0d0*aa*cc 00050 !quad 00051 !quad rs1 = rs(itf,ipf)*(1.0d0 + delta) 00052 !quad 00053 !quad rs2 = 0.0d0 00054 !quad !rs2a = (-bb+sqrt(det))/2.0d0/aa 00055 !quad!uryu if( det >= 0 ) then 00056 !quad if( det >= 0 .and. aa > 0 ) then 00057 !quad rs2 = (-bb-sqrt(det))/2.0d0/aa ! Correct solution 00058 !quad!uryu if ((itf==(ntf/2).or.itf==(ntf/4)).and.ipf<=2) & 00059 !quad!uryu write(6,'(2i4,1p,10e15.5)') & 00060 !quad!uryu itf,ipf,rs1,rs2,r1,r2,r3,aa,bb,cc,det 00061 !quad else 00062 !quad!uryu rs2 = -bb/2.0d0/aa ! Approximate solution. Connect to minimum. 00063 !quad rs2 = rsnew(itf,ipf) 00064 !quad!uryu if ((itf==(ntf/2).or.itf==(ntf/4)).and.ipf<=2) & 00065 !quad!uryu write(6,'(2i4,1p,e15.5,1p,10e15.5)') & 00066 !quad!uryu itf,ipf, rs1, rs2, r1,r2,r3, aa,bb,cc,det 00067 !quad end if 00068 !quad rsnew(itf,ipf) = rs2 00069 ! ================================================================== 00070 ! 00071 end do 00072 end do 00073 rsnew(0,1:npf) = rsnew(0,0) 00074 end subroutine calc_surface
