C------------------------------------------------------------------------------
C     nprho.for
C------------------------------------------------------------------------------
! 30 Apr 2010 - Slight change in setting model labeling
! 31 Mar 2010 - Alternate way to define a structure for a value plus gradients
C 24 Mar 2010 - Changed subroutine name to NPRHO
C 24 Mar 2010 - Added ENTRY PDENSITYa
C 24 Mar 2010 - Added ENTRY IPDENSITYa
C     19 Jul 2006 - ipdensity (no argument) and ipdensityg (w. argument)
C     26 Jan 2006 - 
C     25 Jan 2006 - version for gravity waves
C     GENERIC DENSITY PERTURBATION MODEL
C------------------------------------------------------------------------------
      SUBROUTINE nprho
C------------------------------------------------------------------------------

      implicit double precision (a-h,o-z)

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

      TYPE :: values_and_derivs_1_2
         REAL(KIND=8) :: value
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: ptmtm,prtm,pthtm,pphtm
         REAL(KIND=8) :: ptmr,prr,pthr,pphr
         REAL(KIND=8) :: ptmth,prth,pthth,pphth
         REAL(KIND=8) :: ptmph,prph,pthph,pphph
      END TYPE 
      TYPE (values_and_derivs_1_2) :: rho, potentialrho

C------------------------------------------------------------------------------
! An alternate way to define a structure for a value plus first and second 
!     derivatives:                                                ! 31 Mar 2010

      TYPE :: gradient
         DOUBLE PRECISION :: tm,r,th,ph
      END TYPE gradient

      TYPE :: gradient2
         TYPE (gradient) :: tm,r,th,ph
      END TYPE gradient2

      TYPE :: value_and_gradients_1_2
         DOUBLE PRECISION :: value
         TYPE (gradient) :: p
         TYPE (gradient2) :: pp
      END TYPE value_and_gradients_1_2

      TYPE (value_and_gradients_1_2) rho2, potentialrho2

!      rho2%value = 1
!      rho2%p%r = 2
!      rho2%pp%r%th = 23
!      rho2%pp%r%ph = 24
!      rho2%pp%ph%th = 43
!      rho2%pp%ph%ph = 44

!      PRINT*, 'rho2 = ', rho2                                 ! testing types

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

      ENTRY PDENSITY(rho,potentialrho)

      ENTRY PDENSITYg(rho,potentialrho)

      RETURN

      ENTRY PDENSITYa                                             ! 24 Mar 2010

      RETURN

C------------------------------------------------------------------------------
      END SUBROUTINE nprho
C------------------------------------------------------------------------------


C------------------------------------------------------------------------------
      SUBROUTINE SETPDENSITY
C------------------------------------------------------------------------------
C
C     This subroutine is called once.
C
      implicit double precision (a-h,o-z)

C     COMMON DECK "B16" INSERTED HERE
      INTEGER DTMX,DTNTBL,DTITBL,DTFRMTB,IDSDT(10)
      COMMON/B16/DTMX,DTNTBL(10),DTITBL(10),DTFRMTB(10),DTGP(10)
      EQUIVALENCE (DTGP,IDSDT)
C
      INTEGER DWMX ,DWNTBL(10),DWITBL(10),DWFRMTB(10)
C
      DATA DWMX/1/
      DATA DWNTBL/1,11,8*0/
      DATA DWITBL/1,9*0/
      DATA DWFRMTB/1,9*0/
C
      DTMX=DWMX
      CALL IMOVE(DTNTBL,DWNTBL,10)
      CALL IMOVE(DTITBL,DWITBL,10)
      CALL IMOVE(DTFRMTB,DWFRMTB,10)
C
      RETURN
C
C------------------------------------------------------------------------------
      END SUBROUTINE SETPDENSITY
C------------------------------------------------------------------------------


C------------------------------------------------------------------------------
      SUBROUTINE IPDENSITY  ! 19 Jul 2006
C------------------------------------------------------------------------------
C     This subroutine is called once for each new W array.

      implicit double precision (a-h,o-z)

      TYPE :: model
         CHARACTER :: model_name
         REAL(KIND=8) :: data_id
         CHARACTER :: pert_name
         REAL(KIND=8) :: pert_id
      END TYPE model
      TYPE (model) :: moddg

C     COMMON DECK "DS" INSERTED HERE
      DOUBLE PRECISION MODDS
      COMMON/DS/MODDS(4), ds,Pdst,PdsR,PdsTH,PdsPH
C     COMMON DECK "WW" INSERTED HERE
      PARAMETER (NWARSZ=1000)
      COMMON/WW/dum(10), MAXW,W(NWARSZ)
C
C     DELTA DENSITY    675-799
      EQUIVALENCE (W(676),DDFORM)
C
C------------------------------------------------------------------------------
      ENTRY IPDENSITYa                                            ! 24 Mar 2010
C------------------------------------------------------------------------------
C     This subroutine is called once for each new W array.

      IF(W(675) .NE. 0.0) THEN
           PRINT*,' Model check number mismatch: '
           STOP ' W(675) should be 0 if you are using NPRHO'
      ENDIF
C
      MODDS(3) = 'NPRHO'
      MODDS(4) = W(677)
      RETURN

C------------------------------------------------------------------------------
      ENTRY IPDENSITYg(moddg)
C------------------------------------------------------------------------------
C     This subroutine is called once for each new W array.

      IF(W(675) .NE. 0.0) THEN
           PRINT*,' Model check number mismatch: '
           STOP ' W(675) should be 0 if you are using NPRHO'
      ENDIF
C
      moddg%pert_name = 'NPRHO'
      moddg%pert_id = W(677)
!      MODDS(3) = moddg%pert_name                                 ! 30 Apr 2010
!      MODDS(4) = moddg%pert_id                                   ! 30 Apr 2010
      MODDS(3) = 'NPRHO'                                          ! 30 Apr 2010
      MODDS(4) = W(677)                                           ! 30 Apr 2010

      RETURN

C------------------------------------------------------------------------------
      END SUBROUTINE IPDENSITY
C------------------------------------------------------------------------------