C	this is a "standard shell" only,   fsubmodel & fsubaniso

C	Sven Ivansson 1990




	SUBROUTINE FSUBMODEL ( IFSMTYPLOK,ZLOK, UA2LOK,UB2LOK,RHOLOK )
C	input:    ifsmtyplok  model type, parameter continuity within each 
C                             model type is assumed
C	          zlok        depth in km
C	output:   ua2         "squared P-slowness" in (s/km)**2
C	          ub2         "squared S-slowness" in (s/km)**2
C	          rho         density in "10**12 kg"/(km**3)

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C	nol <= noldim

	INCLUDE 'INMODEL.INC'
C	possible checks could be (using the usual layer index id) : 
C	  ifsmtyplok = ifsmtyp(id) > 0
C	  zd(id-1) <= zlok <= zd(id)
C	  nss < id <= nsw  --->  ub2lok := (0.0d0,0.0d0)

	INTEGER IFSMTYPLOK
	REAL*8 ZLOK

	REAL*8 RHOLOK
	COMPLEX*16 UA2LOK,UB2LOK

C	real*8 qfakt, aa,bb,dbla,dblb, slask
C	complex*16 cslask
C	qfakt = 20 * 3.141592653589793 * dlog10(dexp(1.0d0))
C
C	slask = (zlok-z0) / (z1-z0)
C	aa = (1-slask)*aa0 + slask*aa1
C	dbla = (1-slask)*dbla0 + slask*dbla1
C	cslask = dcmplx ( 1.0d0 , 0.5d0*dbla/qfakt ) / aa
C	ua2lok = cslask * cslask
C	bb = (1-slask)*bb0 + slask*bb1
C	dblb = (1-slask)*dblb0 + slask*dblb1
C	cslask = dcmplx ( 1.0d0 , 0.5d0*dblb/qfakt ) / bb
C	ub2lok = cslask * cslask
C	rholok = (1-slask)*rholok0 + slask*rholok1
C

	IF (IFSMTYPLOK.EQ.1) THEN
	  UA2LOK = DCMPLX ( 1/2.25D0 , 0.0D0 )    ! 1/1.5**2
	  UB2LOK = (0.0D0,0.0D0)
	  RHOLOK = 1.0D0
	ELSE IF (IFSMTYPLOK.EQ.2) THEN
	  UA2LOK = DCMPLX ( 1/9.0D0 , 0.0D0 )     ! 1/3.0**2
	  UB2LOK = DCMPLX ( 1/3.0D0 , 0.0D0 )     ! 1/(3.0/SQRT(3.0))**2
	  RHOLOK = 2.0D0
	ELSE
	  STOP 760
	END IF

	CALL U2AZIM ( UA2LOK, UA2LOK )
	CALL U2AZIM ( UB2LOK, UB2LOK )

ccc	%%% NOTE !!! NO TEST WITH U2TEST IS DONE !!! BEWARE OF THIS !!!
ccc	             TO BE "SAFE": RUN "RMPLOT" WHERE THIS TEST IS DONE

	RETURN
	END



	SUBROUTINE FSUBANISO ( IFSMTYPLOK,ZLOK, UA2LOK,UB2LOK,RHOLOK, 
     $	UC2LOK,UF2LOK )
C	input:    ifsmtyplok  model type, parameter continuity within each 
C                             model type is assumed
C	          zlok        depth in km
C	output:   ua2         "squared P-slowness" in (s/km)**2
C	          ub2         "squared S-slowness" in (s/km)**2
C	          rho         density in "10**12 kg"/(km**3)
C	          uc2         "squared C-slowness" in (s/km)**2
C	          uf2         "squared F-slowness" in (s/km)**2
C	( see White "Underground Sound" )

	IMPLICIT NONE

	INCLUDE 'RPPARAMETER.INC'
C	nol <= noldim

	INCLUDE 'INMODEL.INC'
C	possible checks could be (using the usual layer index id) : 
C	  ifsmtyplok = ifsmtyp(id) > 0
C	  zd(id-1) <= zlok <= zd(id)
C	  nss < id <= nsw  --->  ub2lok := (0.0d0,0.0d0)

	INTEGER IFSMTYPLOK
	REAL*8 ZLOK

	REAL*8 RHOLOK
	COMPLEX*16 UA2LOK,UB2LOK, UC2LOK,UF2LOK

	IF (IFSMTYPLOK.EQ.1) THEN
	  UA2LOK = DCMPLX ( 1/2.25D0 , 0.0D0 )    ! 1/1.5**2
	  UB2LOK = (0.0D0,0.0D0)
	  RHOLOK = 1.0D0
	  UC2LOK = UA2LOK
	  UF2LOK = UA2LOK
	ELSE IF (IFSMTYPLOK.EQ.2) THEN
	  UA2LOK = DCMPLX ( 1/9.0D0 , 0.0D0 )     ! 1/3.0**2
	  UB2LOK = DCMPLX ( 1/3.0D0 , 0.0D0 )     ! 1/(3.0/SQRT(3.0))**2
	  RHOLOK = 2.0D0
	  UC2LOK = UA2LOK
	  UF2LOK = 1/UA2LOK - 2/UB2LOK
	ELSE
	  STOP 760
	END IF

	CALL U2AZIM ( UA2LOK, UA2LOK )
	CALL U2AZIM ( UB2LOK, UB2LOK )
	CALL U2AZIM ( UC2LOK, UC2LOK )
	CALL U2AZIM ( UF2LOK, UF2LOK )

ccc	%%% NOTE !!! NO TEST WITH U2TEST IS DONE !!! BEWARE OF THIS !!!
ccc	             TO BE "SAFE": RUN "RMPLOT" WHERE THIS TEST IS DONE

	RETURN
	END