! copyright (c) 2008 Swedish Defence Research Agency (FOI) 
! may not be copied or distributed without written permission from FOI

C***********************************************************************
C  SHf: "SH with equivalent fluid density" (se achenbach baksid 58)
C    * only for fluid (no solid) homogeneous (no airy) layers
C    * not for "local reaction"
C    * DBLA as "-Im(c**2) / [2*Re(c**2)]" (se zhang et al GJI22,19)
C        hence UA2 = 1 / [AA*AA*(1-2*i*DBLA)]
C                           < NOT  UA = [1+i*DBLA/(2*QFAKT)] / AA >
C    * RHO as "densities for equivalent fluid case with DBLA=0"
C        hence in computations RHO -> RHO * "c**2 for DBLA=0"/c**2
C                                   = RHO * UA2*Re[1/UA2]
C                                       < NOT  RHO * UA2/[Re(UA)]**2 >
C***********************************************************************

C	subroutines rmodel,u2azim,paramz  for rpress.f etc

C	Sven Ivansson  1990
C	option "nsourcfw<-1" for "beam source array" added may 2020 (for rmodpfw)



	SUBROUTINE RMODEL ( LUNIN,LUNOUT, MHQEFF )
C	lunin,lunout, mhqeff  are input
C	  ( if  lunin = 0  the "lunin data" are taken from a previous call )
C	  ( the "nsh" are multiplied with mhqeff whenever  "iimhqani" = 1,2 )
C	  ( special "iimhqani" = -3,-4 option for sponge-rfl minimization )
C	common /inmodel/ is updated
C	  ( if  nsourcfw > 0  variables in "sourcfw.inc" are updated and 
C	                   hydrophones in the solid are allowed )
C	input-file records (apart from initial for "zsourc","..fw","top") : 
C	     " zz  iifsmtyp iimhqani "
C	         for  iifsmtyp = 0 
C	             " aa bb rrho  dbla dblb   nsh nh "
C	             " cc ff       dblc  dblf "        for "aniso" (extra line)
C	         for  iifsmtyp nonzero 
C	             "    bb                   nsh nh "
C	             except for the very first case when i=0
C	     note:  "anisotropic layer" characterized by  iimhqani = -2,2
C	                   iisolid  determined from  bb & iimhqani
C	            "fsubmodel layer"   characterized by  iifsmtyp nonzero 
C	                   then  iabs(iifsmtyp)  will be the "fsubmodel 
C	                      (fsubaniso) type"
C	                   iifsmtyp < 0  will imply approximating discrete 
C	                      layers
C	                   iifsmtyp > 0  will imply  "fsubmodel (fsubaniso) 
C	                      continuation"
C	            "fluid layer"       characterized by  nonzero bb
C	            "rigid bottom,top"  characterized by  rrho < 0.0q0
C	            "free  bottom,top"  characterized by  rrho = 0.0q0
C	             nsh = "parts (with a - for gradient)"
C	                   ( 0 for "new start" (only for iifsmtyp=0) )
C	             nh < 0  will imply "1-nh" hydrophones with one at the 
C	                   upper boundary (which must then be the upper bounda-
C	                   ry for the fluid and for the first hydrophone)
C	zsourc := matching depth (also source depth if nsourcfw=0)
C	ua2sourc  := "ua2 at the matching depth"
C	rhosourc  := "rho at the matching depth"
C	nss := "nr of non-halfspace solid layers above the fluid"
C	nso := "nr of non-halfspace layers above the matching depth (source)"
C	nsw := "nr of upper solid layers and fluid layers disregarding 
C	       the half-spaces"
C	nol := "number of layers including the lower half-space but not the 
C	       upper half-space"
C	  % one must always have  " 0 <= nss <= nso <= nsw < nol " 
C	    and the fluid layers must not be interlaced with solid layers
C	  % when "nss=nso" or "nso=nsw", such a source is nevertheless "on the 
C	    fluid side" when "nsourcfw=0"
C	  % it is possible to have fluid throughout ("nss=0" and "nol=nsw+1", 
C	    use boundaries that are not homogeneous solid half-spaces)
C	  % it is possible to have solid throughout ("nss=nso=nsw=0", take 
C	    "ns=nshc=0" to indicate this, and use boundaries that are not 
C	    homogeneous fluid half-spaces)
C	zd(0:nol)  := lower layer depths (km),  zd(nol):=-1.0q0
C	dz(0:nol)  := layer thicknesses (km),  dz(0):=-1.0q0 dz(nol):=-1.0q0
C	isolid(0:nol) 
C	           :=  0 for fluid layer (which is always isotropic) and 
C	                     'homogeneous fluid half-space or free or rigid 
C	                      boundary (isolid(0) isolid(nol))', 
C	               1 for isotropic solid layer or 'homogeneous solid 
C	                      half-space boundary (isolid(0) isolid(nol))', 
C	              -1 for transversely isotropic solid layer with vertical 
C	                     symmetry 
C	                     ( see  White (1983) "Underground Sound", 
C	                            Takeuchi & Saito in Bolt (ed.) "Methods in 
C	                               Computational Physics" vol. 11, 
C	                            Mandal & Mitchell in J Geophys 59 (1986), 
C	                            Levin in Geophysics 44 (1979) )
C	ifsmtyp(0:nol) 
C	           := -7,-6,-5,-4 for a "usual" layer with constant parameters 
C	                   below "like" & "different" refer to "dz,ua2,ub2,rho"
C	                   then  -7  for  "like the layer below" 
C	                         -6  for  "like both adjacent layers" 
C	                         -5  for  "like the layer above" 
C	                         -4  for  "different from adjacent layers" 
C	              -3,-2,-1,0  for a "usual" layer with "ua2 gradient"
C	                   only allowed in the fluid between the half-spaces 
C	                   ( i.e. for  nss < id <= nsw  )
C	                   below "like" & "different" refer to "dz,ub2,rho"
C	                   then  -3  for  "like the layer below" 
C	                         -2  for  "like both adjacent layers" 
C	                         -1  for  "like the layer above" 
C	                          0  for  "different from adjacent layers" 
C	              > 0         for a "fsubmodel/cont" layer
C	                   then  ifsmtyp  gives the model-type identification 
C	                             for calls to fsubmodel (or fsubaniso)
C	ua2(0:nol) := squared complex P-slownesses ( (s/km)**2 ) such that 
C	            the "ua" have positive real parts and nonnegative 
C	            imaginary parts  ( "ua" = cdsqrt("ua2") )  
C	            ( except when  "ifsmtyp" > 0 , then use "fsubmodel/cont" )
C	            ( set to zero when "irrelevant" )
C	            ( "exp(-i*omeg*t)" with Re(omeg)>0.0 is behind )
C	            note: when  " ifsmtyp = -3,-2,-1,0 "  "ua2" is squared 
C	            P-slowness at the middle of the layer
C	ub2(0:nol) := squared complex S-slownesses ( (s/km)**2 ) such that 
C	            the "ub" have positive real parts and nonnegative 
C	            imaginary parts  ( "ub" = cdsqrt("ub2") )
C	            ( except when  "ifsmtyp" > 0 , then use "fsubmodel/cont" )
C	            ( set to zero when "irrelevant" )
C	            ( "exp(-i*omeg*t)" with Re(omeg)>0.0 is behind )
C	            however, when  " ifsmtyp = -3,-2,-1,0 "  "ub2" is "squared 
C	            P-slowness" difference ( " bottom - top " ) for the layer 
C	rho(0:nol) := densities ( "10**12 kg"/(km**3), or kg/(dm**3) )  
C	            ( except when  "ifsmtyp" > 0 , then use "fsubmodel/cont" )
C	            ( set to zero when "irrelevant" )
C	            note: a rigid half-space is characterized by  
C	                  " rho(nol) < 0.0q0 "   " rho(0) < 0.0q0 "
C	                  EXCEPT " rho(0) = -99.0q0 " for "multi adaptation"
C	                         " rho(nol) = -99.0q0 " for "local reaction"
C	            note: a vacuum half-space is characterized by  
C	                  " rho(nol) = 0.0q0 "   " rho(0) = 0.0q0 "
C	twovb2(0:nol) := 2/ub2(0:nol)      ( but zero when "irrelevant" )
C	uc2(0:nol) := (for anisotropic layer) squared complex "C"-slownesses 
C	            ( (s/km)**2 ) such that the "uc" have positive real parts 
C	            and nonnegative imaginary parts  ( "uc" = cdsqrt("uc2") )  
C	            ( except when  "ifsmtyp" > 0 , then use "fsubaniso/cont" )
C	            ( set to zero when "irrelevant" )
C	            ( "exp(-i*omeg*t)" with Re(omeg)>0.0 is behind )
C	uf2(0:nol) := (for anisotropic layer) squared complex "F"-slownesses 
C	            ( (s/km)**2 ) such that the "uf" have positive real parts 
C	            and nonnegative imaginary parts  ( "uf" = cdsqrt("uf2") )  
C	            ( except when  "ifsmtyp" > 0 , then use "fsubaniso/cont" )
C	            ( set to zero when "irrelevant" )
C	            ( "exp(-i*omeg*t)" with Re(omeg)>0.0 is behind )
C	nhyd := "number of hydrophones"
C	idhyd(1:nhyd) := "interface numbers for the hydrophones" (increasing)
C	  % for rpress,  " nss <= idhyd(1) < ... < idhyd(nhyd) <= nsw "  needed
C	ihdd(0:nol) := "hydrophone nr:s for the interfaces" ( ihdd(nol)=0 )
C	            ( 0 for an interface without a hydrophone )

	USE INMODELa

	IMPLICIT NONE

	INCLUDE '../RPPARAMETER.INC'
C	nol <= noldim , nhyd <= nhyddim
	INTEGER NPKDIM
	REAL*16 TOLHOM
	PARAMETER ( NPKDIM = 2200, TOLHOM = 1.0Q-14 )

	INCLUDE '../SOURCFW.INC'
	INCLUDE 'BEAMARR.INC'

	INTEGER LUNIN,LUNOUT, MHQEFF

	INTEGER NPKT, ISOURC, NS,NSHC, ID, ID0,IDX, IH, I,K,M, NN, 
     $  NSHLOK, II0, NSOURCFW_TMP, 
     $	IISOLID(0:npkdim), IIFSMTYP(-1:npkdim), IIMHQANI(0:npkdim), 
     $	NSH(-1:npkdim), NH(0:npkdim), ISOLID0,IFSMTYP0
	REAL*16 ZSOURC0,ZSOURC1, AA,BB,DBLA,DBLB, CC,FF,DBLC,DBLF, 
     $	RHOOLD,RHONEW, QFAKT, SLASK, SLASK1, DD, 
     $	ZZ(-1:npkdim),RRHO(0:npkdim), DZ0,RHO0
	COMPLEX*32 UA2OLD,UB2OLD,UC2OLD,UF2OLD, UA2NEW,UB2NEW,UC2NEW,
     $	UF2NEW, UUA2SOURC, CSLASK, QCMPLX, 
     $	UUA2(-1:npkdim),UUB2(-1:npkdim), UUC2(0:npkdim),UUF2(0:npkdim)
	LOGICAL DIMDONE

	SAVE NPKT, ISOURC, NS, NSOURCFW_TMP, ZZ,IISOLID,IIFSMTYP,IIMHQANI, 
     &	 UUA2SOURC, UUA2,UUB2,RRHO,NSH,NH, UUC2,UUF2, 
     & 	 ISOLID0,IFSMTYP0, DZ0,RHO0, DIMDONE
	DATA DIMDONE/.FALSE./

	QFAKT =20 * 3.14159265358979323846264338327950280Q0 * LOG10(EXP(1.0Q0))

	IF (MHQEFF.LE.0) STOP 570

C	------------------------------------------------

	IF (LUNIN.NE.0) THEN

	  READ(LUNIN,*) ZSOURC0,ZSOURC1, NS,NSHC
	  IF ( ZSOURC0.GT.ZSOURC1 .OR. NSHC.LT.0 ) STOP 640
	  IF ( NS.LT.0 .OR. NSHC.LT.NS ) STOP 640
	  READ(LUNIN,*) NSOURCFW
	  IF (NSOURCFW.EQ.0) THEN
C	    source at the matching depth
	    NSOURCFW_TMP = NSOURCFW
	    NSARR = 0  ! nsarr now usable as indicator
	  ELSE IF (NSOURCFW.GT.0) THEN
C	    source array
	    IF (NSOURCFW.GE.nsourcfwdim) STOP 3617  ! ".ge.": paramz special
	    DO I = 1, NSOURCFW
	      READ(LUNIN,*) NSOFW0(I), AMSOURCFW(I)
	    END DO
	    IF (NSOFW0(1).LE.0) STOP 4162
	    DO I = 2, NSOURCFW
	      IF ( NSOFW0(I) .LE. NSOFW0(I-1) ) STOP 4163
	    END DO
	    NSOURCFW_TMP = NSOURCFW
	    NSARR = 0  ! nsarr now usable as indicator
	  ELSE IF (NSOURCFW.EQ.-1) THEN
	    NSOURCFW_TMP = NSOURCFW
C	    beam array (see "24/8 07")
C	      nsarr:  nr of beam sources, pos/neg for below/above "nso"
C	      angarr: in degrees, pos/neg for downwards/upwards beam
C	      cutarr: weight at outmost beam sources rel center
	    READ(LUNIN,*) NSARR, ANGARR,CUTARR
	    IF (NSARR.GT.0) THEN
	      IDOWNARR = +1  ! beam array just below "nso"
	    ELSE
	      IDOWNARR = -1  ! beam array just above "nso"
	      NSARR = -NSARR
	    END IF	    
	    IF (NSARR.LT.2) STOP '%%% rmodel, nsarr not needed'
	    IF (NSARR.GT.nsourcfwdim) STOP '%%% rmodel, nsarr dim'
	    IF ( ABS(ANGARR) .GT. 90.0Q0 ) STOP '%%% rmodel, angarr'
	    IF (CUTARR.LE.0.0Q0) STOP '%%% rmodel, cutarr'
	  ELSE
C	    beam source array
	    IF (-NSOURCFW.GE.nsourcfwdim) STOP 3617  ! ".ge.": paramz special
	    READ(LUNIN,*) NSOFW0(1)  ! amsourcfw(1) to be set later
	    IF (NSOFW0(1).LE.0) STOP 4162
	    DO I = 2, -NSOURCFW
	      NSOFW0(I) = NSOFW0(I-1) + 1  ! amsourcfw(i) to be set later
	    END DO
	    NSOURCFW_TMP = -NSOURCFW
	    NSARR = 0  ! nsarr now usable as indicator
C           nr of beam sources = "-nsourcfw" >=2
C	      angarr: in degrees, pos/neg for downwards/upwards beam
C	      cutarr: weight at outmost beam sources rel center
	    READ(LUNIN,*) ANGARR,CUTARR
	    IF ( ABS(ANGARR) .GT. 90.0Q0 ) STOP '%%% rmodel, angarr'
	    IF (CUTARR.LE.0.0Q0) STOP '%%% rmodel, cutarr'
	  END IF

	  READ(LUNIN,*) AA,BB,RHO0,DBLA,DBLB
	  IF (RHO0.GT.0.0Q0) THEN
	    UUA2(-1) = 1 / ( AA * AA * QCMPLX (1.0Q0,-2*DBLA) )
	    IF (BB.LT.0.0Q0) BB = AA/SQRT(3.0Q0)
	    IF (BB.NE.0.0Q0) THEN
	      ISOLID0 = 1
	      IF (DBLB.LT.0.0Q0) DBLB = DBLA*(9.0Q0/4.0Q0)
	      CSLASK = QCMPLX ( 1.0Q0 , 0.5Q0*DBLB/QFAKT ) / BB
	      UUB2(-1) = CSLASK * CSLASK
	    ELSE
	      ISOLID0 = 0
	      IF (DBLB.NE.0.0Q0) STOP 5440
	      UUB2(-1) = (0.0Q0,0.0Q0)
	    END IF
	  ELSE
	    ISOLID0 = 0
	    UUA2(-1) = (0.0Q0,0.0Q0)
	    UUB2(-1) = (0.0Q0,0.0Q0)
	    IF (RHO0.EQ.-99.0Q0) THEN
C	      "multi" adaptation above
	       READ(LUNIN,*) ICONTCUT0
	    ELSE IF (RHO0.LT.0.0Q0) THEN
	      RHO0 = -1.0Q0
	    END IF
	  END IF
	  DZ0 = -1.0Q0
	  IFSMTYP0 = -4

	  ISOURC = 0
	  IIFSMSOME = 0
	  FSMDZMIN = -1.0Q0
	  IIMHQSOME = 0
	  I = -1

	  DO WHILE (.TRUE.)

	    I = I + 1
	    IF (I.GT.npkdim) STOP 4501
	    READ(LUNIN,*,END=21) ZZ(I),IIFSMTYP(I),IIMHQANI(I)
	    IF (I.EQ.0) THEN
	      IF (IIFSMTYP(I).EQ.0) THEN
		ZZ(-1) = ZZ(0)
	        IIFSMTYP(-1) = 0
	        II0 = 0
	      ELSE
	        II0 = 1
	        GOTO 11
	      END IF
	    ELSE IF (I.EQ.1) THEN
	      IF (II0.EQ.1) THEN
	        IF (IIFSMTYP(I).NE.IIFSMTYP(I-1)) STOP 3761
	        IF (IIMHQANI(I).NE.IIMHQANI(I-1)) STOP 3761
	      END IF
	    END IF
	    IF (IIFSMTYP(I).EQ.0) THEN
	      READ(LUNIN,*) AA,BB,RRHO(I),DBLA,DBLB,
     $	      NSH(I),NH(I)   !  nsh = "pre-parts (with a - for gradient)"
	      IF (RRHO(I).EQ.-99.0Q0) THEN
C               local reaction below
	        IF ( I.LE.1 .OR. IISOLID(I-1).NE.0 )
     &	          STOP '%%% rmodel, local-reaction problem'
	        IF ( AA.NE.0.0Q0 .OR. BB.NE.0.0Q0 .OR. DBLB.NE.0.0Q0 ) 
     &	         STOP '%%% rmodel, local-reaction problem'
	        READ(LUNIN,*) ZLOCREACTNOL  ! se pierce ovan s 107
	        IF (DBLA.EQ.0.0Q0) STOP '%%% rmodel, local-reaction problem'
	        IF (DBLA.LT.0.0Q0) THEN
	          CSLASK = SQRT(UUA2(I-1))
		  IF ( REAL(CSLASK).LE.0.0Q0 .OR. AIMAG(CSLASK).LT.0.0Q0 )
     &	           STOP '%%% rmodel, uua2 problem'
	          ZLOCREACTNOL = ZLOCREACTNOL * RRHO(I-1)/CSLASK !given normalized
	        END IF
	        DBLA = 0.0Q0  ! indicator no longer needed
	      ELSE IF (RRHO(I).LE.0.0Q0) THEN
	        AA = 0.0Q0
	        BB = 0.0Q0
	        DBLA = 0.0Q0
	        DBLB = 0.0Q0
	        IF (RRHO(I).LT.0.0Q0) RRHO(I) = -1.0Q0
	      END IF
	      IF (I.EQ.0) THEN
	        IF (NSH(0).NE.0) STOP '%%% rmodel, nsh(0).ne.0'
	        IF (NH(0).NE.0) STOP '%%% rmodel, nh(0).ne.0'
	      END IF
	    ELSE
	      READ(LUNIN,*) BB,NSH(I),NH(I)   !  nsh as above
	      IF (IIFSMTYP(I).GT.0) THEN
	        IF (IIFSMSOME.EQ.0) THEN
	          IIFSMSOME = 1
	          FSMDZMIN = ( ZZ(I) - ZZ(I-1) ) / NSH(I)
	        ELSE
	          FSMDZMIN = MIN ( (ZZ(I)-ZZ(I-1))/NSH(I), FSMDZMIN )
	        END IF
	      END IF
	    END IF
	    IF (IIMHQANI(I).GT.0) THEN
	      IF (IIMHQANI(I).GT.2) STOP 2403
	      IF ( IIFSMTYP(I).GT.0 .OR. NSH(I).EQ.0 ) STOP 1405
	      IIMHQSOME = 1
	    END IF
	    IF (IABS(IIMHQANI(I)).NE.2) THEN
	      IF ( IIMHQANI(I).NE.0 .AND. IIMHQANI(I).NE.1 .AND. 
     &	       IIMHQANI(I).NE.-3 .AND. IIMHQANI(I).NE.-4 ) 
     &	       STOP '%%% rmodel, iimhqani'
	      IF (BB.EQ.0.0Q0) THEN
	        IISOLID(I) = 0
	      ELSE
	        IISOLID(I) = 1
	      END IF
	    ELSE
	      IF (IIFSMTYP(I).EQ.0) THEN
	        READ(LUNIN,*) CC,FF,DBLC,DBLF
	        IF ( DBLC.LT.0.0Q0 .OR. DBLF.LT.0.0Q0 ) STOP 3400
	      END IF
	      IISOLID(I) = -1
	      IF (BB.EQ.0.0Q0) STOP 2300 ! no fluid anisotropic layer allowed
	    END IF
	    IF (NSH(I).EQ.0) THEN
	      IF (ZZ(I).NE.ZZ(I-1)) STOP 2810
	    ELSE
	      IF (ZZ(I).LT.ZZ(I-1)) STOP 2811
	    END IF
	    IF (IIFSMTYP(I).EQ.0) THEN
	      IF (NSH(I).NE.0) THEN
		IF (I.EQ.0) STOP 6570
	        IF (IISOLID(I).NE.IISOLID(I-1)) STOP 6571
	      END IF
	      IF (IIFSMTYP(I-1).NE.0) THEN
	        IF (NSH(I).NE.0) STOP 3001
	      END IF
	    ELSE
	      IF (NSH(I).EQ.0) STOP 2821
	      IF (BB.GT.0.0Q0) STOP 3012
	      IF ( IIFSMTYP(I).GT.0 .AND. NSH(I).LE.0 ) STOP 3023
	      IF (IIFSMTYP(I-1).EQ.0) THEN
	        IF (NSH(I-1).EQ.0) STOP 3034
	      END IF
	    END IF

	    IF (IISOLID(I).EQ.0) THEN
C	      fluid now
	      IF ( IIFSMTYP(I).EQ.0 .AND. DBLB.NE.0.0Q0 ) STOP 5440
	    ELSE
C	      solid now
	      IF (NSH(I).LT.0) STOP 5060   ! gradient allowed only in fluid
	      IF (IIFSMTYP(I).EQ.0) THEN
	        IF (BB.LT.0.0Q0) BB = AA/SQRT(3.0Q0)
	        IF (DBLB.LT.0.0Q0) DBLB = DBLA*(9.0Q0/4.0Q0)
	      END IF
	    END IF

	    IF (IIFSMTYP(I).EQ.0) THEN
	      IF (AA.NE.0.0Q0) THEN
	        UUA2(I) = 1 / ( AA * AA * QCMPLX (1.0Q0,-2*DBLA) )
	      ELSE
	        UUA2(I) = (0.0Q0,0.0Q0)
	      END IF
	      IF (IISOLID(I).NE.0) THEN
	        CSLASK = QCMPLX ( 1.0Q0 , 0.5Q0*DBLB/QFAKT ) / BB
	        UUB2(I) = CSLASK * CSLASK
	      ELSE
	        UUB2(I) = (0.0Q0,0.0Q0)
	      END IF
	      IF (IISOLID(I).EQ.-1) THEN
	        CSLASK = QCMPLX ( 1.0Q0 , 0.5Q0*DBLC/QFAKT ) / CC
	        UUC2(I) = CSLASK * CSLASK
	        CSLASK = QCMPLX ( 1.0Q0 , 0.5Q0*DBLF/QFAKT ) / FF
	        UUF2(I) = CSLASK * CSLASK
	      END IF
C	      note that the complex slownesses are the primary entities
	    END IF

	    IF (I.GT.0) THEN
	      IF (NSH(I).NE.0) THEN
	        IF (ZZ(I).EQ.ZSOURC1) THEN
	          IF (ZZ(I-1).EQ.ZSOURC0) THEN
	            IF ( IISOLID(I).EQ.0 .AND. IABS(NSH(I)).EQ.NSHC ) 
     $	              THEN
	              IF ( 0.LT.NS .AND. NS.LT.NSHC ) THEN
	                SLASK1 = DBLE(NS) / NSHC
	                ZSOURC = (1-SLASK1)*ZSOURC0 + SLASK1*ZSOURC1
	              ELSE IF (NS.EQ.0) THEN
	                SLASK1 = 0.0Q0
	                ZSOURC = ZSOURC0
	              ELSE
	                SLASK1 = 1.0Q0
	                ZSOURC = ZSOURC1
	              END IF
	              ISOURC = I
	              IF ( IIMHQANI(I).EQ.-3 .OR. IIMHQANI(I).EQ.-4 ) 
     &	               STOP '%%% rmodel, iimhqani'
	              IF (IIFSMTYP(I).EQ.0) THEN
	                UUA2SOURC = (1-SLASK1)*UUA2(I-1) + SLASK1*UUA2(I)
	                RHOSOURC = (1-SLASK1)*RRHO(I-1) + SLASK1*RRHO(I)
	              ELSE
	                STOP '%%% rmodel, SHf'  ! avoid fsubmodel
	              END IF
	            END IF
	          END IF
	        END IF
	      END IF
	    END IF

 11	    CONTINUE
	  END DO

 21	  IF (ISOURC.EQ.0) THEN
	    IF ( ZSOURC0.NE.ZZ(0) .OR. ZSOURC1.NE.ZZ(0) ) STOP 4504
	    IF ( NS.NE.0 .OR. NSHC.NE.0 ) STOP 4507 !"solid below" tested later
	    ZSOURC = ZZ(0)
	    UUA2SOURC = (1/1.5Q0)**2   ! dummy
	    RHOSOURC = 1.0Q0           ! dummy
	  END IF
	  NPKT = I - 1
	  IF (NPKT.LT.0) STOP 1010
	  IF ( IISOLID(NPKT).LT.0 .OR. IIFSMTYP(NPKT).NE.0 .OR. 
     $	  IIMHQANI(NPKT).NE.0 .OR. NSH(NPKT).NE.0 .OR. NH(NPKT).NE.0 ) 
     $	  STOP 5555

	  IF ( (NSOURCFW.GT.0.OR.NSOURCFW.LT.-1) .AND. IIMHQSOME.NE.0 ) THEN
	    DO I = 1, NPKT
	      IF (NSH(I).GT.0) THEN
	        IF (IIMHQANI(I).LE.0) STOP 3516  !"nsofw() mhqeff" modif
	      END IF
	    END DO
	  END IF

	END IF

C	------------------------------------------------

	IF (.NOT.DIMDONE) THEN
	  NOLDIM = 0
	  NHYDDIM = 0
	  DO I = 1, NPKT
	    IF (IIMHQANI(I).LE.0) THEN
	      NSHLOK = NSH(I)
	    ELSE
	      NSHLOK = MHQEFF * NSH(I)
	    END IF
	    IF (IISOLID(I).EQ.0) THEN
	      IF (NH(I).LT.0) NHYDDIM = 1  ! fluid now
	    END IF
	    NHYDDIM = NHYDDIM + IABS(NH(I))
	    NOLDIM = NOLDIM + IABS(NSHLOK)
	  END DO
	  NOLDIM = NOLDIM + 1
	  IF (NHYDDIM.EQ.0) NHYDDIM = 1
	  ALLOCATE ( ZD(0:noldim),DZ(0:noldim), 
     $	   ISOLID(0:noldim),IFSMTYP(0:noldim), 
     $     UA2(0:noldim),UB2(0:noldim),RHO(0:noldim), TWOVB2(0:noldim), 
     $	   UC2(0:noldim),UF2(0:noldim), 
     $	   IDHYD(nhyddim),IHDD(0:noldim) )  ! INMODELa
C	  index 0 for upper half-space, index nol for lower half-space
	  DIMDONE = .TRUE.
	END IF

	ISOLID(0) = ISOLID0
	IFSMTYP(0) = IFSMTYP0
	DZ(0) = DZ0
	RHO(0) = RHO0
	CALL U2AZIM ( UUA2(-1), UA2(0) )
	CALL U2AZIM ( UUB2(-1), UB2(0) )
	CALL U2TEST (QFAKT,UA2(0),UB2(0))
	IF (UB2(0).NE.(0.0Q0,0.0Q0)) THEN
	  TWOVB2(0) = 2 / UB2(0)
	ELSE
	  TWOVB2(0) = (0.0Q0,0.0Q0)
	END IF
	UC2(0) = (0.0Q0,0.0Q0)
	UF2(0) = (0.0Q0,0.0Q0)

	CALL U2AZIM ( UUA2SOURC, UA2SOURC )

	NSS = 0
	NSW = 0
	NOL = 0
	NOLHOML = 0
	NSO = 0     ! not set below if isourc.eq.0
	NHYD = 0

	IF (RRHO(NPKT).GT.0.0Q0) THEN
	  DO II0 = NPKT-1, 1, -1
	    IF (IIFSMTYP(II0).NE.0) GOTO 140
	    IF (IABS(IIMHQANI(II0)).EQ.2) GOTO 140
	    IF ( ABS(RRHO(II0)-RRHO(NPKT))/RRHO(NPKT) .GT. 
     $	    TOLHOM ) GOTO 140
	    IF (IISOLID(II0).NE.IISOLID(NPKT)) GOTO 140
	    IF ( ABS(UUA2(II0)-UUA2(NPKT))/ABS(UUA2(NPKT)) .GT. 
     $	    TOLHOM ) GOTO 140
	    IF (IISOLID(II0).EQ.1) THEN
	      IF ( ABS(UUB2(II0)-UUB2(NPKT))/ABS(UUB2(NPKT)) .GT. 
     $	      TOLHOM ) GOTO 140
	    END IF
	  END DO
	  II0 = 0
 140	  CONTINUE
	ELSE
	  II0 = NPKT - 1
	END IF

	DO I = 1, NPKT

	  IF (IIMHQANI(I).LE.0) THEN
	    NSHLOK = NSH(I)
	  ELSE
	    NSHLOK = MHQEFF * NSH(I)
	  END IF

	  IF (IISOLID(I).EQ.0) THEN
C	    fluid now
	    IF ( NSW.NE.NOL .AND. I.LT.NPKT ) STOP 5433
	    IF (I.EQ.ISOURC) THEN
	      IF (IIMHQANI(I).LE.0) THEN
	        NSO = NOL + NS
	      ELSE
	        NSO = NOL + MHQEFF*NS
	      END IF
	    END IF
	    IF (NH(I).LT.0) THEN
	      IF ( NSS.NE.NSW .OR. NHYD.NE.0 ) STOP 6060
	      NHYD = 1
	      IDHYD(NHYD) = NSS
	    END IF
	    NSW = NSW + IABS(NSHLOK)
	  ELSE
C	    solid now
	    IF (ISOURC.GT.0) THEN
	      IF (NSS.EQ.NSW) THEN
	        NSS = NSS + NSHLOK
	        NSW = NSS
	      END IF
	    END IF
	    IF (NH(I).LT.0) STOP 43167
	  END IF
	  IF (NH(I).NE.0) THEN
	    IF (NSHLOK.EQ.0) STOP 1332
	    NN = IABS(NSHLOK)/IABS(NH(I))
	    IF ( NN.EQ.0 .OR. (NN*IABS(NH(I))).NE.IABS(NSHLOK) ) STOP 1321
	    DO K = 1, IABS(NH(I))
	      NHYD = NHYD + 1
	      IF (NHYD.GT.nhyddim) STOP 7013
	      IDHYD(NHYD) = NOL + K*NN
	    END DO
	  END IF
	  NOL = NOL + IABS(NSHLOK)
	  IF (I.EQ.II0) NOLHOML = NOL

	END DO

	NOL = NOL + 1
	IF (NOL.GT.noldim) STOP 5556
	NOLHOML = NOLHOML + 1
	IF ( NSS.LT.0 .OR. NSO.LT.NSS .OR. NSW.LT.NSO .OR. NOL.LE.NSW ) THEN
	  WRITE(*,*) '%%% nss nso nsw nol: ',NSS,NSO,NSW,NOL
	  STOP 4000
	END IF
	IF ( IIMHQSOME.EQ.0 .OR. MHQEFF.EQ.1 ) THEN
	  DO I = 1, NSOURCFW_TMP
	    NSOFW(I) = NSOFW0(I)
	    IF (NSOFW(I).GE.NOL) STOP 4162
	  END DO
	ELSE
C	  "partial implementation only"
	  DO I = 1, NSOURCFW_TMP
	    NSOFW(I) = MHQEFF * NSOFW0(I)
	    IF (NSOFW(I).GE.NOL) STOP 4163
	  END DO
	  IF (NSOURCFW.EQ.-1) STOP '%%% rmodel, beam array & mhqeff'
	END IF
	IF (NSOURCFW.EQ.-1) THEN
	  IF (IDOWNARR.GT.0) THEN
	    IF (NSO.GE.NSW) STOP '%%% rmodel, beam array & idownarr'
	  ELSE
	    IF (NSO.LE.NSS) STOP '%%% rmodel, beam array & idownarr'
	  END IF
        END IF

	ID = 0
	ZD(0) = ZZ(0)

	DO I = 1, NPKT

	  IF (IIMHQANI(I).LE.0) THEN
	    NSHLOK = NSH(I)
	  ELSE
	    NSHLOK = MHQEFF * NSH(I)
	  END IF
	  ID0 = ID
	  IDX = IABS(NSHLOK)
	  ID = ID0 + IDX
	  IF (IDX.GT.0) THEN
	    DD = (ZZ(I)-ZZ(I-1)) / IDX
	    ZD(ID) = ZZ(I)
	    DZ(ID) = DD
	    ISOLID(ID) = IISOLID(I)
	    DO K = IDX-1, 1, -1
	      SLASK1 = DBLE(K) / IDX
	      SLASK = 1 - SLASK1
	      M = ID0 + K
	      ZD(M) = SLASK*ZZ(I-1) + SLASK1*ZZ(I)
	      DZ(M) = DD
	      ISOLID(M) = IISOLID(I)
	    END DO
	  END IF

	  IF (IIFSMTYP(I).EQ.0) THEN

	    IF (NSHLOK.GT.0) THEN

	      IF (IDX.EQ.1) THEN
	        IFSMTYP(ID) = -4
	      ELSE IF ( UUA2(I-1).NE.UUA2(I) .OR. UUB2(I-1).NE.UUB2(I) 
     $	      .OR. RRHO(I-1).NE.RRHO(I) ) THEN
	        DO M = ID0+1, ID
	          IFSMTYP(M) = -4
	        END DO
	      ELSE IF (IISOLID(I).GE.0) THEN
	        IFSMTYP(ID0+1) = -7
	        IFSMTYP(ID) = -5
	        DO M = ID0+2, ID-1
	          IFSMTYP(M) = -6
	        END DO
	      ELSE
	        IF ( UUC2(I-1).NE.UUC2(I) .OR. UUF2(I-1).NE.UUF2(I) ) 
     $	        THEN
	          DO M = ID0+1, ID
	            IFSMTYP(M) = -4
	          END DO
	        ELSE
	          IFSMTYP(ID0+1) = -7
	          IFSMTYP(ID) = -5
	          DO M = ID0+2, ID-1
	            IFSMTYP(M) = -6
	          END DO
	        END IF
	      END IF
	      CALL U2AZIM ( UUA2(I-1), UA2OLD )
	      CALL U2AZIM ( UUA2(I), UA2NEW )
	      CALL U2AZIM ( UUB2(I-1), UB2OLD )
	      CALL U2AZIM ( UUB2(I), UB2NEW )
	      DO K = 1, IDX
	        SLASK1 = (K-0.5Q0) / IDX
	        IF ( IIMHQANI(I) .EQ. -3 ) THEN
	          SLASK1 = SLASK1*SLASK1*SLASK1  ! for special sponge layer
	        ELSE IF ( IIMHQANI(I) .EQ. -4 ) THEN
	          SLASK1 = SLASK1*SLASK1         ! for special sponge layer
	        END IF
	        SLASK = 1 - SLASK1
	        M = ID0 + K
	        UA2(M) = SLASK*UA2OLD + SLASK1*UA2NEW
	        UB2(M) = SLASK*UB2OLD + SLASK1*UB2NEW
	        CALL U2TEST (QFAKT,UA2(M),UB2(M))
	        RHO(M) = SLASK*RRHO(I-1) + SLASK1*RRHO(I)
	        IF (UB2(M).NE.(0.0Q0,0.0Q0)) THEN
	          TWOVB2(M) = 2 / UB2(M)
	        ELSE
	          TWOVB2(M) = (0.0Q0,0.0Q0)
	        END IF
	      END DO
	      IF (IISOLID(I).EQ.-1) THEN
	        CALL U2AZIM ( UUC2(I-1), UC2OLD )
	        CALL U2AZIM ( UUC2(I), UC2NEW )
	        CALL U2AZIM ( UUF2(I-1), UF2OLD )
	        CALL U2AZIM ( UUF2(I), UF2NEW )
	        DO K = 1, IDX
	          SLASK1 = (K-0.5Q0) / IDX
	          SLASK = 1 - SLASK1
	          M = ID0 + K
	          UC2(M) = SLASK*UC2OLD + SLASK1*UC2NEW
	          UF2(M) = SLASK*UF2OLD + SLASK1*UF2NEW
	        END DO
	      ELSE
	        DO K = 1, IDX
	          M = ID0 + K
	          UC2(M) = (0.0Q0,0.0Q0)
	          UF2(M) = (0.0Q0,0.0Q0)
		END DO
	      END IF

	    ELSE IF (NSHLOK.LT.0) THEN

	      IF (IDX.EQ.1) THEN
	        IFSMTYP(ID) = 0
	      ELSE IF (RRHO(I-1).NE.RRHO(I)) THEN
	        DO M = ID0+1, ID
	          IFSMTYP(M) = 0
	        END DO
	      ELSE
	        IFSMTYP(ID0+1) = -3
	        IFSMTYP(ID) = -1
	        DO M = ID0+2, ID-1
	          IFSMTYP(M) = -2
	        END DO
	      END IF
	      CALL U2AZIM ( UUA2(I-1), UA2OLD )
	      CALL U2TEST (QFAKT,UA2OLD,(0.0Q0,0.0Q0))
	      CALL U2AZIM ( UUA2(I), UA2NEW )
	      CALL U2TEST (QFAKT,UA2NEW,(0.0Q0,0.0Q0))
	      CSLASK = (UA2NEW-UA2OLD) / IDX
	      IF (CSLASK.EQ.(0.0Q0,0.0Q0)) STOP 4507  ! gradient unnecessary
	      DO K = 1, IDX
	        SLASK1 = (K-0.5Q0) / IDX
	        SLASK = 1 - SLASK1
	        M = ID0 + K
	        UA2(M) = SLASK*UA2OLD + SLASK1*UA2NEW
	        UB2(M) = CSLASK
	        RHO(M) = SLASK*RRHO(I-1) + SLASK1*RRHO(I)
	        TWOVB2(M) = (0.0Q0,0.0Q0)
	        UC2(M) = (0.0Q0,0.0Q0)
	        UF2(M) = (0.0Q0,0.0Q0)
	      END DO

	    ELSE

	      IF (ZZ(I).NE.ZZ(I-1)) STOP 3174

	    END IF

	  ELSE IF (IIFSMTYP(I).LT.0) THEN

	    IF (NSHLOK.GT.0) THEN

	      IF (IISOLID(I).GE.0) THEN
	        STOP '%%% rmodel, SHf'  ! avoid fsubmodel
	      ELSE
	        STOP '%%% rmodel, SHf'  ! avoid fsubaniso
	      END IF
	      IF (IISOLID(I).EQ.0) THEN
	        IF (UB2OLD.NE.(0.0Q0,0.0Q0)) STOP 5701
	      ELSE
	        IF (UB2OLD.EQ.(0.0Q0,0.0Q0)) STOP 5702
	      END IF
	      DO M = ID0+1, ID
	        IFSMTYP(M) = -4
	        IF (IISOLID(I).GE.0) THEN
	          STOP '%%% rmodel, SHf'  ! avoid fsubmodel
	        ELSE
	          STOP '%%% rmodel, SHf'  ! avoid fsubaniso
	        END IF
	        IF (IISOLID(I).EQ.0) THEN
	          IF (UB2NEW.NE.(0.0Q0,0.0Q0)) STOP 5703
	        ELSE
	          IF (UB2NEW.EQ.(0.0Q0,0.0Q0)) STOP 5704
	        END IF
	        UA2(M) = 0.5Q0 * ( UA2OLD + UA2NEW )
	        UB2(M) = 0.5Q0 * ( UB2OLD + UB2NEW )
	        CALL U2TEST (QFAKT,UA2(M),UB2(M))
	        RHO(M) = 0.5Q0 * ( RHOOLD + RHONEW )
	        IF (UB2(M).NE.(0.0Q0,0.0Q0)) THEN
	          TWOVB2(M) = 2 / UB2(M)
	        ELSE
	          TWOVB2(M) = (0.0Q0,0.0Q0)
	        END IF
	        UA2OLD = UA2NEW
	        UB2OLD = UB2NEW
	        RHOOLD = RHONEW
	        IF (IISOLID(I).EQ.-1) THEN
	          UC2(M) = 0.5Q0 * ( UC2OLD + UC2NEW )
	          UF2(M) = 0.5Q0 * ( UF2OLD + UF2NEW )
	          UC2OLD = UC2NEW
	          UF2OLD = UF2NEW
		ELSE
	          UC2(M) = (0.0Q0,0.0Q0)
	          UF2(M) = (0.0Q0,0.0Q0)
	        END IF
	      END DO

	    ELSE

	      IF (IISOLID(I).NE.0) STOP 4507

	      STOP '%%% rmodel, SHf'  ! avoid fsubmodel

	    END IF

	  ELSE

	    DO M = ID0+1, ID
	      IFSMTYP(M) = IIFSMTYP(I)
	      UA2(M) = (0.0Q0,0.0Q0)
	      UB2(M) = (0.0Q0,0.0Q0)
	      RHO(M) = 0.0Q0
	      TWOVB2(M) = (0.0Q0,0.0Q0)
	      UC2(M) = (0.0Q0,0.0Q0)
	      UF2(M) = (0.0Q0,0.0Q0)
	    END DO

	  END IF

	END DO

	ZD(NSO) = ZSOURC    ! be careful and avoid a numerical difference
	DO I = 1, NSOURCFW_TMP
	  IF (ISOLID(NSOFW(I)).EQ.-1) STOP 3416   ! "aniso" not "fw" impl
	  CALL PARAMZ ( ZD(NSOFW(I)), NN, UA2SOURCFW(I),
     $	  UB2SOURCFW(I),RHOSOURCFW(I) )  ! note: special usage, see "paramz"
	  IF (NN.NE.2) STOP 3416      ! only "both and equal" allowed
	END DO
	IF (NSOURCFW.EQ.-1) THEN
C	  beam array
	  IF ( UA2(NSO+IDOWNARR) .NE. UA2SOURC ) STOP '%% beam arr, ua2'
	  IF ( RHO(NSO+IDOWNARR) .NE. RHOSOURC ) STOP '%% beam arr, rho'
	END IF

	IF (RHO(0).EQ.-99.0Q0) THEN
	   IF ( NSS.NE.0 .OR. NSW.LT.1 ) STOP '%rmodel, "multi" probl'
	END IF

	IF ((ID+1).NE.NOL) STOP 5103
	ZD(NOL) = ZD(NOL-1)      ! artificial
      	DZ(NOL) = -1.0Q0         ! artificial
	ISOLID(NOL) = IISOLID(NPKT)
	IFSMTYP(NOL) = -4
	CALL U2AZIM ( UUA2(NPKT), UA2(NOL) )
	CALL U2AZIM ( UUB2(NPKT), UB2(NOL) )
	CALL U2TEST (QFAKT,UA2(NOL),UB2(NOL))
	RHO(NOL) = RRHO(NPKT)
	IF (RHO(NOL).EQ.-99.0Q0) THEN
	  IF ( NSS.EQ.NSW .OR. NOL.NE.NSW+1 ) STOP '%rmodel, loc-react probl'
	END IF
	IF (UB2(NOL).NE.(0.0Q0,0.0Q0)) THEN
	  TWOVB2(NOL) = 2 / UB2(NOL)
	ELSE
	  TWOVB2(NOL) = (0.0Q0,0.0Q0)
	END IF
	UC2(NOL) = (0.0Q0,0.0Q0)
	UF2(NOL) = (0.0Q0,0.0Q0)

	DO ID = 0, NOL
	  IHDD(ID) = 0
	END DO
	DO IH = 1, NHYD
	  IHDD(IDHYD(IH)) = IH
	END DO

	IF (NSOURCFW.EQ.-1) THEN
	  IF (IDOWNARR.GT.0) THEN
C	    receiver depths only allowed below beam source depths
	    DO IH = 1, NHYD
	      IF ( IDHYD(IH) .LE. NSO ) STOP '%%% beam array, idhyd'
	    END DO
	  ELSE
C	    receiver depths only allowed above beam source depths
	    DO IH = 1, NHYD
	      IF ( IDHYD(IH) .GE. NSO ) STOP '%%% beam array, idhyd'
	    END DO
	  END IF
	END IF

	DO ID = 1, NSS
	  IF (ISOLID(ID).EQ.0) STOP 5617
	END DO
	DO ID = NSS+1, NSW
	  IF (ISOLID(ID).NE.0) STOP 5617
	END DO
	DO ID = NSW+1, NOL-1
	  IF (ISOLID(ID).EQ.0) STOP 5617
	END DO
	IF ( RHO(0).LE.0.0Q0 .AND. ISOLID(0).NE.0 ) STOP 5616
	IF ( RHO(NOL).LE.0.0Q0 .AND. ISOLID(NOL).NE.0 ) STOP 5616
	IF (RHO(NOL).LE.0.0Q0) THEN
	  IF (NOLHOML.NE.NOL) STOP 5615
	ELSE
	  IF ( ISOLID(NOLHOML) .NE. ISOLID(NOL) ) STOP 5615
	  IF (ISOLID(NOL).NE.0) THEN
	    IF (NOLHOML.LE.NSW) STOP 5615
	  END IF
	END IF
	IF (NSS.EQ.NSW) THEN
	  IF (NSW.NE.0) STOP 5614
	  IF (ISOLID(1).NE.1) STOP 5613  ! "solid below" required when 'nsw=0'
	END IF

C	------------------------------------------------

	IF (LUNOUT.GT.0) THEN

	  WRITE(LUNOUT,*)
	  WRITE(LUNOUT,*)
	  WRITE(LUNOUT,*) ' IIFSMSOME,IIMHQSOME are ',IIFSMSOME,
     $	  IIMHQSOME
	  WRITE(LUNOUT,*) ' IAZIMI(temp),MHQEFF now: ',IAZIMI,MHQEFF
	  WRITE(LUNOUT,*)
	  WRITE(LUNOUT,*) ' ZSOURC   is ',ZSOURC, '  ! matching depth'
	  WRITE(LUNOUT,*) ' UASOURC  is ',SQRT(UA2SOURC)
	  WRITE(LUNOUT,*) ' RHOSOURC is ',RHOSOURC
	  CALL U2TOREAL ( QFAKT,UA2SOURC, AA,DBLA )
	  WRITE(LUNOUT,*) ' "1/Re(uasourc) dblasourc" are ',SNGL(AA),
     $	  SNGL(DBLA)
	  WRITE(LUNOUT,*)
	  IF (NSOURCFW.EQ.0) THEN
	    WRITE(LUNOUT,*) ' SOURCE AT THE MATCHING DEPTH'
	  ELSE IF (NSOURCFW.GT.0) THEN 
	    WRITE(LUNOUT,*) ' SOURCE ARRAY (NSOURCFW=',NSOURCFW,' )  ',
     $      ' nr  depth  Re(strength) Im(strength) :'
	    DO I = 1, NSOURCFW
	      WRITE(LUNOUT,13) I, SNGL(ZD(NSOFW(I))), 
     $        SNGL(REAL(AMSOURCFW(I))),SNGL(AIMAG(AMSOURCFW(I)))
 13	      FORMAT(1X,1X,I4,2X,F12.6,3X,F12.6,1X,F12.6)
	    END DO
	  ELSE IF (NSOURCFW.EQ.-1) THEN 
	    WRITE(LUNOUT,*) ' BEAM ARRAY AT THE MATCHING DEPTH'
	    WRITE(LUNOUT,*) ' nsarr,idownarr=',NSARR,IDOWNARR
	    WRITE(LUNOUT,*) ' angarr=',ANGARR
	    WRITE(LUNOUT,*) ' cutarr=',CUTARR
          ELSE
	    WRITE(LUNOUT,*) ' BEAM SOURCE ARRAY (|NSOURCFW|=',-NSOURCFW,' )  ',
     $      ' nr  depth :'
	    DO I = 1, -NSOURCFW
	      WRITE(LUNOUT,'(1X,1X,I4,2X,F12.6)') I, SNGL(ZD(NSOFW(I)))
	    END DO
	    WRITE(LUNOUT,*) ' angarr=',ANGARR
	    WRITE(LUNOUT,*) ' cutarr=',CUTARR
	  END IF
	  WRITE(LUNOUT,*)
	  WRITE(LUNOUT,*) ' GIVEN DEPTH MODEL :'
	  WRITE(LUNOUT,*) ' IIMHQANI ZZ   1/RE(UA) 1/RE(UB)   ',
     $	  'RHO   DBLA DBLB   NSH IH'
	  CALL U2TOREAL ( QFAKT,UA2(0), AA,DBLA )
	  CALL U2TOREAL ( QFAKT,UB2(0), BB,DBLB )
	  WRITE(LUNOUT,14) SNGL(AA),SNGL(BB),SNGL(RHO(0)),
     $	  SNGL(DBLA),SNGL(DBLB)
 14	  FORMAT(1X,1X,2X,1X,12X,3F10.6,1X,2F10.6)
	  DO I = 0, NPKT
	    IF (IIMHQANI(I).LE.0) THEN
	      NSHLOK = NSH(I)
	    ELSE
	      NSHLOK = MHQEFF * NSH(I)
	    END IF
	    IF (IIFSMTYP(I).EQ.0) THEN
	      CALL U2AZIM ( UUA2(I), UA2NEW )
	      CALL U2AZIM ( UUB2(I), UB2NEW )
	      CALL U2TOREAL ( QFAKT,UA2NEW, AA,DBLA )
	      CALL U2TOREAL ( QFAKT,UB2NEW, BB,DBLB )
	      WRITE(LUNOUT,23) IIMHQANI(I),SNGL(ZZ(I)),SNGL(AA),
     $	      SNGL(BB),SNGL(RRHO(I)),SNGL(DBLA),SNGL(DBLB),NSHLOK,NH(I)
 23	      FORMAT(1X,1X,I2,1X,F12.6,3F10.6,1X,2F10.6,1X,I5,1X,I3)
	      IF (IISOLID(I).EQ.-1) THEN
	        CALL U2AZIM ( UUC2(I), UC2NEW )
	        CALL U2AZIM ( UUF2(I), UF2NEW )
	        CALL U2TOREAL ( QFAKT,UC2NEW, CC,DBLC )
	        CALL U2TOREAL ( QFAKT,UF2NEW, FF,DBLF )
	        WRITE(LUNOUT,27) '      aniso ', 
     $	        SNGL(CC),SNGL(FF), SNGL(DBLC),SNGL(DBLF)
 27	        FORMAT(1X,2X,1X,1X,A12,2F10.6,10X,1X,2F10.6)
	      END IF
	    ELSE IF (I.GT.0) THEN
	      WRITE(LUNOUT,33) IIMHQANI(I),SNGL(ZZ(I)),
     $	      '     "fsubmodel"',' (lower limit) ( iifsmtyp =',
     $	      IIFSMTYP(I),' )    ',NH(I)
 33	      FORMAT(1X,1X,I2,1X,F12.6,A16,A27,I5,A6,3X,1X,I3)
	      IF (IISOLID(I).EQ.-1) WRITE(LUNOUT,27) '      aniso '
	    ELSE
	      WRITE(LUNOUT,33) IIMHQANI(I),SNGL(ZZ(I)),
     $	      '     "fsubmodel"',' top                       '
	    END IF
	  END DO

	  WRITE(LUNOUT,*)
	  WRITE(LUNOUT,*) ' NHYD: ',NHYD
	  WRITE(LUNOUT,*) ' NSS,NSO,NSW,NOL: ',NSS,NSO,NSW,NOL
	  WRITE(LUNOUT,*) ' NOLHOML: ',NOLHOML

	  WRITE(LUNOUT,*)
	  WRITE(LUNOUT,*) ' LAYER MODEL (with current mhqeff ',
     $	  'applied for iimhqani=1,2) :'
	  IF (NSS.EQ.0) WRITE(LUNOUT,*)
	  DO ID = 0, NOL
	    IF (IFSMTYP(ID).LE.-4) THEN
	      CALL U2TOREAL ( QFAKT,UA2(ID), AA,DBLA )
	      CALL U2TOREAL ( QFAKT,UB2(ID), BB,DBLB )
	      IF (IHDD(ID).NE.0) THEN
	        WRITE(LUNOUT,43) SNGL(DZ(ID)),SNGL(AA),SNGL(BB),
     $	        SNGL(RHO(ID)),SNGL(DBLA),SNGL(DBLB),IHDD(ID)
	      ELSE
	        WRITE(LUNOUT,43) SNGL(DZ(ID)),SNGL(AA),SNGL(BB),
     $	        SNGL(RHO(ID)),SNGL(DBLA),SNGL(DBLB)
	      END IF
 43	      FORMAT(1X,2X,1X,1X,F12.6,3F10.6,1X,2F10.6,1X,5X,1X,I3)
	      IF (ISOLID(ID).EQ.-1) THEN
	        CALL U2TOREAL ( QFAKT,UC2(ID), CC,DBLC )
	        CALL U2TOREAL ( QFAKT,UF2(ID), FF,DBLF )
	        WRITE(LUNOUT,47) '      aniso ',SNGL(CC),SNGL(FF), 
     $	        SNGL(DBLC),SNGL(DBLF)
 47	        FORMAT(1X,2X,1X,1X,A12,2F10.6,10X,1X,2F10.6)
	      END IF
	    ELSE IF (IFSMTYP(ID).LE.0) THEN
	      CALL U2TOREAL ( QFAKT,UA2(ID), AA,DBLA )
	      IF (IHDD(ID).NE.0) THEN
	        WRITE(LUNOUT,53) SNGL(DZ(ID)),SNGL(AA),' gradient ',
     $	        SNGL(RHO(ID)),SNGL(DBLA),' (middle) ',IHDD(ID)
	      ELSE
	        WRITE(LUNOUT,53) SNGL(DZ(ID)),SNGL(AA),' gradient ',
     $	        SNGL(RHO(ID)),SNGL(DBLA),' (middle) ',IHDD(ID)
	      END IF
 53	      FORMAT
     &	      (1X,2X,1X,1X,F12.6,F10.6,A10,F10.6,1X,F10.6,A10,1X,5X,1X,I3)
	    ELSE
	      IF (IHDD(ID).NE.0) THEN
	        WRITE(LUNOUT,63) SNGL(DZ(ID)),'     "fsubmodel/cont" ',
     $	        'layer ( ifsmtyp = ',IFSMTYP(ID),' ) ',IHDD(ID)
	      ELSE
	        WRITE(LUNOUT,63) SNGL(DZ(ID)),'     "fsubmodel/cont" ',
     $	        'layer ( ifsmtyp = ',IFSMTYP(ID),' ) '
	      END IF
 63	      FORMAT(1X,2X,1X,1X,F12.6,A22,A18,I5,A3,3X,1X,5X,1X,I3)
	      IF (ISOLID(ID).EQ.-1) WRITE(LUNOUT,47) '      aniso '
	    END IF
	    IF (ID.EQ.NSO) WRITE(LUNOUT,*) '           *     zsourc=',
     $	    ZSOURC, '  ! matching depth'
	    IF ( ID.EQ.NSS .OR. ID.EQ.NSW .OR. ID.EQ.NOL ) 
     $	    WRITE(LUNOUT,*)
	  END DO
	  IF (RHO(NOL).GT.0.0Q0) THEN
	    WRITE(LUNOUT,*)
	    WRITE(LUNOUT,*) ' BRANCH POINT(S) IN "0.0<=arg(u)<pi":   ',
     $	    '( u = "up rotated by -arg(omeg)" )'
	    CSLASK = SQRT(UA2(NOL))
	    WRITE(LUNOUT,*) '  for P   ',CSLASK
	    WRITE(LUNOUT,*) '    ',REAL(CSLASK)*AIMAG(CSLASK),
     $	    AIMAG(CSLASK)-REAL(CSLASK)
	    IF (ISOLID(NOL).NE.0) THEN
	      CSLASK = SQRT(UB2(NOL))
	      WRITE(LUNOUT,*) '  for S   ',CSLASK
	      WRITE(LUNOUT,*) '  ',REAL(CSLASK)*AIMAG(CSLASK),
     $	      AIMAG(CSLASK)-REAL(CSLASK)
	    END IF
	  ELSE IF (RHO(NOL).EQ.-99.0Q0) THEN
	    WRITE(LUNOUT,*)
	    WRITE(LUNOUT,*) ' ZLOCREACTNOL: ', ZLOCREACTNOL
	  END IF
	  IF (RHO(0).GT.0.0Q0) THEN
	    WRITE(LUNOUT,*)
	    WRITE(LUNOUT,*) ' BRANCH POINT(S) IN "0.0<=arg(u)<pi":   ',
     $	    '( u = "up rotated by -arg(omeg)" )'
	    CSLASK = SQRT(UA2(0))
	    WRITE(LUNOUT,*) '  for surf-P   ',CSLASK
	    WRITE(LUNOUT,*) '    ',REAL(CSLASK)*AIMAG(CSLASK),
     $	    AIMAG(CSLASK)-REAL(CSLASK)
	    IF (ISOLID(0).NE.0) THEN
	      CSLASK = SQRT(UB2(0))
	      WRITE(LUNOUT,*) '  for surf-S   ',CSLASK
	      WRITE(LUNOUT,*) '  ',REAL(CSLASK)*AIMAG(CSLASK),
     $	      AIMAG(CSLASK)-REAL(CSLASK)
	    END IF
	  ELSE IF (RHO(0).EQ.-99.0Q0) THEN
	    WRITE(LUNOUT,*)
	    WRITE(LUNOUT,*) ' ICONTCUT0: ', ICONTCUT0
	  END IF

	END IF

	RETURN
	END







	SUBROUTINE U2TOREAL ( QFAKT,U2, VELPHASE,DBL )

	IMPLICIT NONE

	REAL*16 QFAKT
	COMPLEX*32 U2

	REAL*16 VELPHASE,DBL

	IF (U2.NE.(0.0Q0,0.0Q0)) THEN
	  VELPHASE = SQRT ( REAL(1/U2) )
	  DBL = -0.5Q0 * AIMAG(1/U2) / REAL(1/U2)
	ELSE
	  VELPHASE = 0.0Q0
	  DBL = 0.0Q0
	END IF

	RETURN
	END






	SUBROUTINE U2AZIM ( U2IN, U2UT )

	USE INMODELa

	IMPLICIT NONE

	COMPLEX*32 U2IN, U2UT

	COMPLEX*32 ULOK, EIMULT

	IF (IAZIMI.EQ.0) THEN
	  U2UT = U2IN
	ELSE IF (U2IN.EQ.(0.0Q0,0.0Q0)) THEN
	  U2UT = (0.0Q0,0.0Q0)
	ELSE
	  ULOK = SQRT(U2IN)
	  ULOK = ULOK + FREAZIM*AIMAG(ULOK)/REAL(ULOK) + 
     $	  EIMULT(-FIMAZIM*AIMAG(ULOK))
	  U2UT = ULOK*ULOK
	END IF

	RETURN
	END








	SUBROUTINE PARAMZ ( Z, INDAB, UA2LOK,UB2LOK,RHOLOK )
C	to get parameters at depth z
C	output:   indab = 0 for "none"
C	                 -1 for "above" only
C	                  1 for "below" only
C	                 -2 for "both" and "unequal"
C	                  2 for "both" and "equal"
C	          ua2lok(1),ub2lok(1),rholok(1)  for "above" if present
C	          ua2lok(2),ub2lok(2),rholok(2)  for "below" if present
ccc NOT YET IMPLEMENTED FOR ANISOTROPIC LAYER !!!

	USE INMODELa

	IMPLICIT NONE

	INCLUDE '../RPPARAMETER.INC'
C	nol <= noldim , nhyd <= nhyddim

	REAL*16 Z

	INTEGER INDAB
	REAL*16 RHOLOK(1:2)
	COMPLEX*32 UA2LOK(1:2),UB2LOK(1:2)

	INTEGER ID, I
	REAL*16 SLASK
	COMPLEX*32 CSLASK

	IF (Z.LT.ZD(0)) THEN

	  IF (RHO(0).LE.0.0Q0) THEN
	    INDAB = 0
	    RETURN
	  ELSE
	    INDAB = 2
	    UA2LOK(1) = UA2(0)
	    UB2LOK(1) = UB2(0)
	    RHOLOK(1) = RHO(0)
	    UA2LOK(2) = UA2LOK(1)
	    UB2LOK(2) = UB2LOK(1)
	    RHOLOK(2) = RHOLOK(1)
	  END IF

	ELSE IF (Z.GT.ZD(NOL-1)) THEN

	  IF (RHO(NOL).LE.0.0Q0) THEN
	    INDAB = 0
	    RETURN
	  ELSE
	    INDAB = 2
	    UA2LOK(1) = UA2(NOL)
	    UB2LOK(1) = UB2(NOL)
	    RHOLOK(1) = RHO(NOL)
	    UA2LOK(2) = UA2LOK(1)
	    UB2LOK(2) = UB2LOK(1)
	    RHOLOK(2) = RHOLOK(1)
	  END IF

	ELSE

	  ID = 1
	  DO WHILE (Z.GT.ZD(ID))
	    ID = ID + 1
	  END DO
C	  now  1 <= id <= nol-1  ,  zd(id-1) < z <= zd(id)  except when z=zd(0)

	  DO I = 1, 2

	    IF (I.EQ.1) THEN
	      IF (Z.EQ.ZD(0)) THEN
	        IF (RHO(0).LE.0.0Q0) THEN
	          INDAB = 0
	        ELSE
	          INDAB = -1   ! not necessarily the final value
	          UA2LOK(I) = UA2(0)
	          UB2LOK(I) = UB2(0)
	          RHOLOK(I) = RHO(0)
	        END IF
	        GOTO 100
	      ELSE
	        INDAB = -1   ! not necessarily the final value
	      END IF
	    ELSE
	      IF (Z.EQ.ZD(0)) THEN
C	        now  id = 1
		IF (INDAB.EQ.0) THEN
	          INDAB = 1
	        ELSE
	          INDAB = -2
	        END IF
	      ELSE IF (Z.EQ.ZD(ID)) THEN
	        IF (ID.EQ.(NOL-1)) THEN
	          IF (RHO(NOL).LE.0.0Q0) GOTO 100
	        END IF
	        ID = ID + 1
C	        now  2 <= id <= nol  ,  zd(id-1) = z < zd(id)
	        IF (INDAB.EQ.0) THEN
	          INDAB = 1
	        ELSE
	          INDAB = -2
	        END IF
	      ELSE
	        IF (INDAB.NE.0) THEN
	          INDAB = 2
	          UA2LOK(2) = UA2LOK(1)
	          UB2LOK(2) = UB2LOK(1)
	          RHOLOK(2) = RHOLOK(1)
	          GOTO 100
	        END IF
	      END IF
	    END IF

	    IF (IFSMTYP(ID).LE.-4) THEN
	      UA2LOK(I) = UA2(ID)
	      UB2LOK(I) = UB2(ID)
	      RHOLOK(I) = RHO(ID)
	    ELSE IF (IFSMTYP(ID).LE.0) THEN
	      SLASK =  ( Z - 0.5Q0*(ZD(ID-1)+ZD(ID)) ) / DZ(ID)
	      UA2LOK(I) = UA2(ID) + SLASK*UB2(ID)
	      UB2LOK(I) = (0.0Q0,0.0Q0)
	      RHOLOK(I) = RHO(ID)
	    ELSE
	      STOP '%%% rmodel, SHf'  ! avoid fsubmodel
	    END IF

 100	    CONTINUE

	  END DO

	  IF (INDAB.EQ.-2) THEN
ccc	    if (ua2lok(2).ne.ua2lok(1)) return
ccc	    if (ub2lok(2).ne.ub2lok(1)) return
ccc	    if (rholok(2).ne.rholok(1)) return
	    CSLASK = UA2LOK(2) - UA2LOK(1)
	    IF (ABS(REAL(CSLASK)).GT.1.0Q-14) RETURN
	    IF (ABS(AIMAG(CSLASK)).GT.1.0Q-14) RETURN
	    CSLASK = UB2LOK(2) - UB2LOK(1)
	    IF (ABS(REAL(CSLASK)).GT.1.0Q-14) RETURN
	    IF (ABS(AIMAG(CSLASK)).GT.1.0Q-14) RETURN
	    IF (ABS(RHOLOK(2)-RHOLOK(1)).GT.1.0Q-14) RETURN
	    UA2LOK(2) = UA2LOK(1)
	    UB2LOK(2) = UB2LOK(1)
	    RHOLOK(2) = RHOLOK(1)
	    INDAB = 2
	  END IF

	END IF

	RETURN
	END
