c
c          Program ----- veljp3d.f -----
c
c    This program is used to calculate 3-D P-wave velocity
c    distribution beneath the Japan Islands which is obtained 
c    by a simultaneous inversion of arrival time data from local,
c    regional and teleseismic events.  For details, see "Deep
c    structure of the Japan subduction zone" by Abdelwahed & Zhao,
c    Phys. Earth Planet. Inter. 162, 32-52, 2007.
c
c       The meaningful range of this model is as follows:
c        latitude : 28 - 45.6 N
c        longitude: 128-146 E
c        depth    : 0  - 700 km
c
c                            Prof. Dapeng Zhao 
c                            Department of Geophysics
c                            Tohoku University
c                            Sendai 980-8578
c                            Japan
c                            zhao@aob.geophys.tohoku.ac.jp
c                                            
      OPEN(2,FILE="JP3D2007")
      OPEN(3,FILE="datadis")
      PID  = 0.017453
      PRINT*," Now reading the 3-D velocity data"
      PRINT*," Please WAIT a moment ......"
      CALL INPUTJP
      CALL INPUT1
      CALL INPUT2
100   PRINT*,"Input Latitude, Longitude & Depth:"
      PRINT*,"  e.g., 36.9 141.3 100.4 ....."
      READ(*,*) PHI,RAM,DEP
      PE     = (90.0-PHI)*PID
      RE     = RAM*PID
      HE     = DEP
c   calculate depths of the Conrad, the Moho and
c   the plate boundary beneath the location (PHI,RAM)
      CALL HLAY(PE,RE,H1,1)
      CALL HLAY(PE,RE,H2,2)
      CALL HLAY(PE,RE,H3,3)
c   when LAY = 1, the focus is in the upper crust;
c   when LAY = 2, the focus is in the lower crust;
c   when LAY = 3, the focus is in the mantle wedge;
c   when LAY = 4, the focus is beneath the plate boundary.
      IF(HE.LE.H1)                   THEN
         LAY = 1
      ELSE IF(HE.GT.H1.AND.HE.LE.H2) THEN
         LAY = 2
      ELSE IF(HE.GT.H2.AND.HE.LE.H3) THEN
         LAY = 3
      ELSE
         LAY = 4
      END IF
      CALL VEL1D(HE,V1,LAY,1)
      CALL VEL3(PE,RE,HE,DV,LAY,1)
      VP    = V1*(1.0+DV*0.01)
      PRINT*," V1d = ",V1,"km/s"
      PRINT*," V3d = ",VP,"km/s"
      PRINT*," Vp-perturbation = ",DV,"%"
      IF(PHI.NE.0.0) GO TO 100
      STOP
      END

      SUBROUTINE INPUT1

	 PARAMETER (MPA=47,MRA=48,MHA=20)
      COMMON/VMOD3D/NPA,NRA,NHA,PNA(MPA),RNA(MRA),HNA(MHA),
     &              VELAP(MPA,MRA,MHA) 
	rewind(2)
100   FORMAT(3I3)
      READ(2,100)  NPA,NRA,NHA
      CALL PUT1(NPA,NRA,NHA,PNA,RNA,HNA,VELAP)
      CALL BLDMAP
      RETURN
      END

      SUBROUTINE PUT1(NPX,NRX,NHX,PNX,RNX,HNX,VELXP)
      DIMENSION  VELXP(NPX,NRX,NHX),
     &           PNX(NPX),RNX(NRX),HNX(NHX)
      READ(2,110) (PNX(I),I=1,NPX)
      READ(2,110) (RNX(I),I=1,NRX)
      READ(2,120) (HNX(I),I=1,NHX)
      DO 1 K = 1,NHX !-1
      DO 1 I = 1,NPX
1     READ(2,140) (VELXP(I,J,K),J=1,NRX)
110   FORMAT(6(9F7.2/))
120   FORMAT(3(8F7.2/))
140   FORMAT(4(14F5.2/))
      RETURN
      END

      SUBROUTINE INPUT2
      COMMON/DISCON/PN(51),RN(63),DEPA(51,63),DEPB(51,63),DEPC(51,63)
      READ(3,100)  NP,NR
      READ(3,110) (PN(I),I=1,NP)
      READ(3,120) (RN(I),I=1,NR)
      DO 1  I = NP,1,-1
      READ(3,130) (DEPA(I,J),J=1,NR)
1     CONTINUE
      DO 2  I = NP,1,-1
      READ(3,130) (DEPB(I,J),J=1,NR)
2     CONTINUE
      DO 3  I = NP,1,-1
      READ(3,130) (DEPC(I,J),J=1,NR)
3     CONTINUE
100   FORMAT(2I6)
110   FORMAT(5(10F7.2/),F7.2)
120   FORMAT(6(10F7.2/),3F7.2)
130   FORMAT(6(10F7.1/),3F7.1)
      RETURN
      END

      SUBROUTINE BLDMAP
      PARAMETER (MKA=3301,MKB=3301)

	PARAMETER (MPA=47,MRA=48,MHA=20) 
      COMMON/VMOD3D/NPA,NRA,NHA,PNA(MPA),RNA(MRA),HNA(MHA),
     &              VELAP(MPA,MRA,MHA)
      COMMON/LOCATE/PLA,RLA,HLA,IPLOCA(MKA),IRLOCA(MKA),IHLOCA(MKA)

      CALL LOCX(PNA,RNA,HNA,NPA,NRA,NHA,MKA,
     &     PLA,RLA,HLA,IPLOCA,IRLOCA,IHLOCA)

      RETURN
      END

      SUBROUTINE LOCX(PNX,RNX,HNX,NPX,NRX,NHX,MKX,
     &           PLX,RLX,HLX,IPLOCX,IRLOCX,IHLOCX)
      DIMENSION  IPLOCX(MKX),IRLOCX(MKX),IHLOCX(MKX),
     &           PNX(NPX),RNX(NRX),HNX(NHX)
      PLX      = 1.0-PNX(1)*100.0
      IPMAX    = IFIX(PNX(NPX)*100.0+PLX)
      IP       = 1
      DO 10 I  = 1,IPMAX
      IP1      = IP+1
      PNOW     = (FLOAT(I)-PLX)/100.0
      IF(PNOW.GE.PNX(IP1))   IP = IP1
      IPLOCX(I)= IP
10    CONTINUE
      RLX      = 1.0-RNX(1)*100.0
      IRMAX    = IFIX(RNX(NRX)*100.0+RLX)
      IR       = 1
      DO 20 I  = 1,IRMAX
      IR1      = IR+1
      RNOW     = (FLOAT(I)-RLX)/100.0
      IF(RNOW.GE.RNX(IR1))   IR = IR1
      IRLOCX(I)= IR
20    CONTINUE
      HLX      = 1.0-HNX(1)
      IHMAX    = IFIX(HNX(NHX)+HLX)
      IH       = 1
      DO 30 I  = 1,IHMAX
      IH1      = IH+1
      HNOW     = FLOAT(I)-HLX
      IF(HNOW.GE.HNX(IH1))   IH = IH1
      IHLOCX(I)= IH
30    CONTINUE
      RETURN
      END

      SUBROUTINE VEL3(PE,RE,HE,V,LAY,IPS)
      PARAMETER (MKA=3301,MKB=3301)

      PARAMETER (MPA=47,MRA=48,MHA=20)
      COMMON/VMOD3D/NPA,NRA,NHA,PNA(MPA),RNA(MRA),HNA(MHA),
     &              VELAP(MPA,MRA,MHA)
      COMMON/LOCATE/PLA,RLA,HLA,IPLOCA(MKA),IRLOCA(MKA),IHLOCA(MKA)

      COMMON/WEIGHT/WV(8),IP,JP,KP,IP1,JP1,KP1
      COMMON/PRHFD/P,R,H,PF,RF,HF,PF1,RF1,HF1,PD,RD,HD
      PIDEG = 0.017453
      P     = 90.0-PE/PIDEG
      R     = RE/PIDEG
      H     = HE
c      IF(LAY.LE.3)       THEN
         CALL PRHF(IPLOCA,IRLOCA,IHLOCA,PLA,RLA,HLA,
     &             PNA,RNA,HNA,MPA,MRA,MHA,MKA)
c      ELSE IF(LAY.EQ.4)  THEN
c         CALL PRHF(IPLOCB,IRLOCB,IHLOCB,PLB,RLB,HLB,
c     &             PNB,RNB,HNB,MPB,MRB,MHB,MKB)
c      ELSE
c      END IF
      WV(1) = PF1*RF1*HF1
      WV(2) = PF*RF1*HF1
      WV(3) = PF1*RF*HF1
      WV(4) = PF*RF*HF1
      WV(5) = PF1*RF1*HF
      WV(6) = PF*RF1*HF
      WV(7) = PF1*RF*HF
      WV(8) = PF*RF*HF
c   calculate velocity
c      IF(LAY.LE.3)      THEN
         CALL VABPS(MPA,MRA,MHA,VELAP,V)
c      ELSE IF(LAY.EQ.4) THEN
c         CALL VABPS(MPB,MRB,MHB,VELBP,V)
c      ELSE
c      END IF
      RETURN
      END

      SUBROUTINE VABPS(MP,MR,MH,V,VEL)
      COMMON/WEIGHT/WV(8),IP,JP,KP,IP1,JP1,KP1
      DIMENSION  V(MP,MR,MH)
      VEL = WV(1)*V(IP,JP,KP)  + WV(2)*V(IP1,JP,KP)
     &    + WV(3)*V(IP,JP1,KP) + WV(4)*V(IP1,JP1,KP)
     &    + WV(5)*V(IP,JP,KP1) + WV(6)*V(IP1,JP,KP1)
     &    + WV(7)*V(IP,JP1,KP1)+ WV(8)*V(IP1,JP1,KP1)
      RETURN
      END

      SUBROUTINE INTMAP(R,IRLOC,NR,RL,IR)
      DIMENSION IRLOC(NR)
      IS      = IFIX(R+RL)
      IR      = IRLOC(IS)
      RETURN
      END

      SUBROUTINE PRHF(IPLOCX,IRLOCX,IHLOCX,PLX,RLX,HLX,
     &                PNX,RNX,HNX,MPX,MRX,MHX,MKX)
      COMMON/WEIGHT/WV(8),IP,JP,KP,IP1,JP1,KP1
      COMMON/PRHFD/P,R,H,PF,RF,HF,PF1,RF1,HF1,PD,RD,HD
      DIMENSION  IPLOCX(MKX),IRLOCX(MKX),IHLOCX(MKX),
     &           PNX(MPX),RNX(MRX),HNX(MHX)
      CALL LIMIT(PNX(1),PNX(MPX),P)
      CALL LIMIT(RNX(1),RNX(MRX),R)
      CALL LIMIT(HNX(1),HNX(MHX),H)

      CALL INTMAP(P*100.0,IPLOCX,MKX,PLX,IP)
      CALL INTMAP(R*100.0,IRLOCX,MKX,RLX,JP)
      CALL INTMAP(H,IHLOCX,MKX,HLX,KP)
      IP1   = IP+1
      JP1   = JP+1
      KP1   = KP+1
      PD    = PNX(IP1)-PNX(IP)
      RD    = RNX(JP1)-RNX(JP)
      HD    = HNX(KP1)-HNX(KP)
      PF    = (P-PNX(IP))/PD
      RF    = (R-RNX(JP))/RD
      HF    = (H-HNX(KP))/HD
      PF1   = 1.0-PF
      RF1   = 1.0-RF
      HF1   = 1.0-HF
      RETURN
      END

      SUBROUTINE HLAY(PE,RE,HE,IJK)
      COMMON/DISCON/PN(51),RN(63),DEPA(51,63),DEPB(51,63),DEPC(51,63)
      PID    = 0.017453
      P      = 90.0-PE/PID
      R      = RE/PID
      CALL LIMIT(PN(1),PN(51),P)
      CALL LIMIT(RN(1),RN(63),R)
      DO 1 I = 1,50
      I1     = I+1
      IF(P.GE.PN(I).AND.P.LT.PN(I1)) GO TO 11
1     CONTINUE
11    CONTINUE
      DO 2 J = 1,62
      J1     = J+1
      IF(R.GE.RN(J).AND.R.LT.RN(J1)) GO TO 22
2     CONTINUE
22    CONTINUE
      PF    = (P-PN(I))/(PN(I1)-PN(I))
      RF    = (R-RN(J))/(RN(J1)-RN(J))
      PF1   = 1.0-PF
      RF1   = 1.0-RF
      WV1   = PF1*RF1
      WV2   = PF*RF1
      WV3   = PF1*RF
      WV4   = PF*RF
      IF(IJK.EQ.1)       THEN
        HE  = WV1*DEPA(I,J)  + WV2*DEPA(I1,J)
     &      + WV3*DEPA(I,J1) + WV4*DEPA(I1,J1)
      ELSE IF(IJK.EQ.2)  THEN
        HE  = WV1*DEPB(I,J)  + WV2*DEPB(I1,J)
     &      + WV3*DEPB(I,J1) + WV4*DEPB(I1,J1)
      ELSE IF(IJK.EQ.3)  THEN
        HE  = WV1*DEPC(I,J)  + WV2*DEPC(I1,J)
     &      + WV3*DEPC(I,J1) + WV4*DEPC(I1,J1)
      ELSE
      END IF
      RETURN
      END

      SUBROUTINE LIMIT(C1,C2,C)
      A1    = AMIN1(C1,C2)
      A2    = AMAX1(C1,C2)
      IF(C.LT.A1)   C = A1
      IF(C.GT.A2)   C = A2
      RETURN
      END

      SUBROUTINE VEL1D(HE,V,LAY,IPS)
      IF(LAY.EQ.1)      THEN
        V    = 6.0
        IF(IPS.EQ.2)    V = 3.5
      ELSE IF(LAY.EQ.2) THEN
        V    = 6.7
        IF(IPS.EQ.2)    V = 3.8
      ELSE IF(LAY.GE.3) THEN
        HM   = 40.0
        IF(HE.LT.HM)    THEN
          CALL JPMODEL(IPS,HM,VM)
          V  = VM-(HM-HE)*0.003
        ELSE 
          CALL JPMODEL(IPS,HE,V)
        END IF
      ELSE
      END IF
      RETURN
      END

      SUBROUTINE INPUTJP
      COMMON/JPMODV/VP(29),VS(29),RA(29),DEPJ(29)
      DIMENSION VP1(29),VS1(29),RA1(29)
      DATA VP1/7.75, 7.94, 8.13, 8.33, 8.54, 8.75, 8.97,
     &         9.50, 9.91,10.26,10.55,10.99,11.29,11.50, 
     &        11.67,11.85,12.03,12.20,12.37,12.54,12.71,  
     &        12.87,13.02,13.16,13.32,13.46,13.60,13.64,13.64/
      DATA VS1/4.353,4.444,4.539,4.638,4.741,4.850,4.962, 
     &         5.227,5.463,5.670,5.850,6.125,6.295,6.395,  
     &         6.483,6.564,6.637,6.706,6.770,6.833,6.893,  
     &         6.953,7.012,7.074,7.137,7.199,7.258,7.314,7.304/
      DATA RA1/1.00,0.99,0.98,0.97,0.96,0.95,0.94,0.93, 
     &         0.92,0.91,0.90,0.88,0.86,0.84,0.82,0.80,  
     &         0.78,0.76,0.74,0.72,0.70,0.68,0.66,0.64,  
     &         0.62,0.60,0.58,0.56,0.55/
      DO 1 L  = 1,29
      VP(L)   = VP1(L)
      VS(L)   = VS1(L)
      RA(L)   = RA1(L)
      DEPJ(L) = 40.0+6325.59*(1.0-RA1(L))
1     CONTINUE
      RETURN
      END

      SUBROUTINE JPMODEL(IPS,H,V)
      COMMON/JPMODV/VP(29),VS(29),RA(29),DEPJ(29)
      DO 2 K = 1,28
      K1     = K+1
      H1     = DEPJ(K)
      H2     = DEPJ(K1)
      IF(H.GE.H1.AND.H.LT.H2) GO TO 3
2     CONTINUE
3     CONTINUE
      H12    = (H-H1)/(H2-H1)
      IF(IPS.EQ.1)  THEN
         V   = (VP(K1)-VP(K))*H12+VP(K)
      ELSE
         V   = (VS(K1)-VS(K))*H12+VS(K)
      END IF
      RETURN
      END