pzo_mod.f90

   1 MODULE PZO_CALC
   2
   3   IMPLICIT NONE
   4
   5 CONTAINS
   6
   7    SUBROUTINE PZOISO(POTSP,POTPZ)
   8         Use Input_Data
   9      IMPLICIT NONE
  10      REAL, INTENT(  OUT) :: POTSP,POTPZ
  11
  12      REAL :: DPZ_1,DPZ_S
  13
  14      !!! COMMON
  15      REAL :: RHO,ZETA,ETA
  16      COMMON /QAGON/RHO,ZETA,ETA
  17
  18      ETA=1.E0
  19      CALL PREST(PZ_1=DPZ_1,PZ_S=DPZ_S)
  20
  21      POTSP=CPZ*CSP*RC*DPZ_1
  22      POTPZ=CPZ*RC*((CPZ_Q1+CPZ_Q2/2.E0)*DPZ_1+CPZ_Q2*DPZ_S)
  23
  24      RETURN
  25
  26    END SUBROUTINE PZOISO
  27
  28    SUBROUTINE PZOANISO(POTSP,POTPZ)
  29         Use Input_Data
  30      IMPLICIT NONE
  31      REAL, INTENT(  OUT) :: POTSP,POTPZ
  32
  33      REAL :: DPZ_1,DPZ_ETA1,DPZ_ETA2
  34
  35      !!! COMMON
  36      REAL :: RHO,ZETA,ETA
  37      COMMON /QAGON/RHO,ZETA,ETA
  38
  39      ETA=1.E0
  40      CALL PREST(PZ_1=DPZ_1,PZ_ETA1=DPZ_ETA1,PZ_ETA2=DPZ_ETA2)
  41
  42      POTSP=CPZ*CSP*RC*DPZ_1
  43      POTPZ=CPZ*RC*( CI1*DPZ_1 +&
  44                         CI2*( CNP31_1*DPZ_ETA2 - CNP31_2*DPZ_ETA1 )/ETA_DIF +&
  45                         CI3/((1.E0-ETA1)*(1.E0-ETA2)) *&
  46                             ( (CNE_2-CNE_1)*DPZ_1 +&
  47                               (CNE2_1*DPZ_ETA2 - CNE2_2*DPZ_ETA1)/ETA_DIF ) )
  48      RETURN
  49
  50    END SUBROUTINE PZOANISO
  51
  52    SUBROUTINE PREST(PZ_1,PZ_S,PZ_ETA1,PZ_ETA2)
  53
  54       Use Input_Data, ONLY: RD,RC,ZC,ETA1,ETA2
  55
  56    IMPLICIT NONE
  57
  58 !! 'dummy' and local variables !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  59
  60    REAL PZ_1
  61    REAL,OPTIONAL :: PZ_S,PZ_ETA1,PZ_ETA2
  62
  63    DOUBLE PRECISION DPZ_1,DPZ_S,DPZ_ETA1,DPZ_ETA2
  64 !! GAG Variables  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  65
  66    INTEGER, PARAMETER :: LIMITE=100
  67    DOUBLE PRECISION AL,AU,ABSERR,EPSABS,EPSREL,WORK(4*LIMITE),RO,ZE
  68    INTEGER IER,IWORK(LIMITE),KEY,LAST,LENW,LIMIT,NEVAL
  69
  70    DOUBLE PRECISION, EXTERNAL :: I_FI00,I_FDA00
  71
  72    !!! COMMON
  73    REAL :: RHO,ZETA,ETA
  74    COMMON /QAGON/RHO,ZETA,ETA
  75
  76 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  77    DOUBLE PRECISION RHOP
  78
  79     RO=RHO
  80     ZE=ZETA
  81     AL=0.D0
  82     AU=DBLE(RD)
  83
  84     EPSABS=0.D0
  85     EPSREL=1.D-4
  86     KEY=6
  87     LIMIT=LIMITE
  88     LENW=LIMIT*4
  89
  90 !   WRITE(6,*)RHO*RC,ZETA*ZC
  91 !   ETA=1.E0
  92 !   DO IER=1,1000
  93 !    RHOP=DBLE(IER-1)*RD/DBLE(999)
  94 !    WRITE(55,'(10(E15.8,1X))')RHOP*RC,I_FDA00(RHOP)
  95 !   END DO
  96 !   STOP
  97
  98     ETA=1.E0
  99     CALL DQAG(I_FI00,AL,AU,EPSABS,EPSREL,KEY,DPZ_1,ABSERR,NEVAL, &
 100              IER,LIMIT,LENW,LAST,IWORK,WORK)
 101     IF (IER.NE.0) &
 102        CALL DQAG_CHECK("FI00","PZ_1",DPZ_1,ABSERR,RO*RC,ZE*ZC,NEVAL,IER,LAST)
 103
 104     PZ_1=REAL(DPZ_1)
 105     IF(PRESENT(PZ_S)) THEN
 106         ETA=1.E0
 107         CALL DQAG(I_FDA00,AL,AU,EPSABS,EPSREL,KEY,DPZ_S,ABSERR,NEVAL, &
 108                  IER,LIMIT,LENW,LAST,IWORK,WORK)
 109         IF (IER.NE.0) &
 110            CALL DQAG_CHECK("FDA00","PZ_S",DPZ_S,ABSERR,RO*RC,ZE*ZC,NEVAL,IER,LAST)
 111
 112         PZ_S=REAL(DPZ_S)
 113     END IF
 114     IF(PRESENT(PZ_ETA1).AND.PRESENT(PZ_ETA2)) THEN
 115         ETA=ETA1
 116         CALL DQAG(I_FI00,AL,AU,EPSABS,EPSREL,KEY,DPZ_ETA1,ABSERR,NEVAL, &
 117                  IER,LIMIT,LENW,LAST,IWORK,WORK)
 118         IF (IER.NE.0) &
 119            CALL DQAG_CHECK("FI00","PZ_ETA1",DPZ_ETA1,ABSERR,RO*RC,ZE*ZC,NEVAL,IER,LAST)
 120
 121         ETA=ETA2
 122         CALL DQAG(I_FI00,AL,AU,EPSABS,EPSREL,KEY,DPZ_ETA2,ABSERR,NEVAL, &
 123                  IER,LIMIT,LENW,LAST,IWORK,WORK)
 124         IF (IER.NE.0) &
 125            CALL DQAG_CHECK("FI00","PZ_ETA2",DPZ_ETA2,ABSERR,RO*RC,ZE*ZC,NEVAL,IER,LAST)
 126         PZ_ETA1=REAL(DPZ_ETA1); PZ_ETA2=REAL(DPZ_ETA2)
 127     END IF
 128
 129    RETURN
 130    END SUBROUTINE PREST
 131
 132    SUBROUTINE DQAG_CHECK(FUNC,RESNAME,RES,ABSERR,RO,ZE,NEVAL,IER,LAST)
 133    IMPLICIT NONE
 134    CHARACTER(LEN=*)    :: FUNC,RESNAME
 135    DOUBLE PRECISION    :: RES,ABSERR,RO,ZE
 136    INTEGER             :: NEVAL,IER,LAST
 137      WRITE(6,*)"ERROR AT: RHO: ",RO," ZETA: ",ZE
 138      WRITE(6,*)"QAG ",TRIM(FUNC)," ",TRIM(RESNAME)," RESULT=",RES
 139      WRITE(6,*)"QAG ABSERR=",ABSERR
 140      WRITE(6,*)"QAG NEVAL=",NEVAL
 141      WRITE(6,*)"QAG IER=",IER
 142      WRITE(6,*)"QAG LAST=",LAST
 143
 144      RETURN
 145    END SUBROUTINE DQAG_CHECK
 146
 147   FUNCTION FI00(RHO,ZETA,RHOP,ZM,ETA) RESULT(I00)
 148         Use Input_Data, ONLY: RC,ZC
 149         Use ELIPT_MOD
 150
 151     IMPLICIT NONE
 152     DOUBLE PRECISION, INTENT(IN) :: RHO,ZETA,RHOP,ZM,ETA
 153     DOUBLE PRECISION :: I00
 154
 155     DOUBLE PRECISION ::  PI,DPI
 156     DOUBLE PRECISION :: DELTAM,DELTA0,&
 157                         KM,K0,CKM,CK0,KKM,KK0
 158
 159
 160       PI  = DATAN(1.D0)*4.D0
 161       DPI = 2.D0*PI
 162
 163       DELTAM = SQRT(ETA)*DABS(ZETA-ZM)*ZC/RC
 164       DELTA0 = SQRT(ETA)*DABS(ZETA)*ZC/RC
 165
 166       IF(RHO.EQ.0.D0) THEN
 167        I00=(1.D0/DSQRT(DELTAM**2+RHOP**2)) - &
 168            (1.D0/DSQRT(DELTA0**2+RHOP**2))
 169        I00 = I00/2.D0
 170        RETURN
 171       END IF
 172
 173       KM = 2.D0*DSQRT( RHO*RHOP/((RHO+RHOP)**2+DELTAM**2) )
 174       K0 = 2.D0*DSQRT( RHO*RHOP/((RHO+RHOP)**2+DELTA0**2) )
 175
 176       CKM = (1.D0-KM**2)
 177       CK0 = (1.D0-K0**2)
 178
 179       CALL ELIPTK0(KM,KKM)
 180       CALL ELIPTK0(K0,KK0)
 181
 182 !!!!! I00 expressed in terms of RHO,RHOP,DELTA !!!!!!!!!!!!!!!!!!!
 183
 184       I00 = 2.D0/PI*(1.D0/DSQRT(DELTAM**2+(RHO+RHOP)**2)*KKM - &
 185            1.D0/DSQRT(DELTA0**2+(RHO+RHOP)**2)*KK0)
 186
 187       I00=I00/2.D0
 188
 189
 190       RETURN
 191
 192
 193   END FUNCTION FI00
 194
 195   FUNCTION FDA00(RHO,ZETA,RHOP,ZM,ETA) RESULT(DA00)
 196         Use Input_Data, ONLY: RC,ZC
 197         Use ELIPT_MOD
 198
 199     IMPLICIT NONE
 200     DOUBLE PRECISION, INTENT(IN) :: RHO,ZETA,RHOP,ZM,ETA
 201     DOUBLE PRECISION :: DA00
 202
 203     DOUBLE PRECISION ::  PI,DPI
 204     DOUBLE PRECISION :: SRRHOP,RRHOP15,DELTAM,DELTA0,&
 205                         KM,K0,CKM,CK0,KKM,KK0,EKM,EK0
 206
 207       PI  = DATAN(1.D0)*4.D0
 208       DPI = 2.D0*PI
 209
 210       SRRHOP=SQRT(RHO*RHOP)
 211       RRHOP15=SQRT(RHO*RHOP)**3
 212
 213       DELTAM = DSQRT(ETA)*DABS(ZETA-ZM)*ZC/RC
 214       DELTA0 = DSQRT(ETA)*DABS(ZETA)*ZC/RC
 215
 216       KM = 2.D0*DSQRT( RHO*RHOP/((RHO+RHOP)**2+DELTAM**2) )
 217       K0 = 2.D0*DSQRT( RHO*RHOP/((RHO+RHOP)**2+DELTA0**2) )
 218
 219       CKM = (1.D0-KM**2)
 220       CK0 = (1.D0-K0**2)
 221
 222       CALL ELIPT0(KM,KKM,EKM)
 223       CALL ELIPT0(K0,KK0,EK0)
 224
 225        IF (RHO .EQ. 0.D0) THEN
 226           DA00 = ( DELTAM**2/SQRT(RHOP**2+DELTAM**2)**3) &
 227                -(DELTA0**2/SQRT(RHOP**2+DELTA0**2)**3 )
 228           DA00 = DA00/4.D0
 229           RETURN
 230        END IF
 231
 232       DA00 = DELTAM**2*KM**3/(2.D0*DPI*CKM*RRHOP15)*EKM &
 233             -DELTA0**2*K0**3/(2.D0*DPI*CK0*RRHOP15)*EK0
 234
 235       DA00 = DA00/4.D0
 236
 237       RETURN
 238
 239   END FUNCTION FDA00
 240
 241
 242 END MODULE PZO_CALC