1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
|
*DECK PKIRHO
SUBROUTINE PKIRHO(IPMAP,NALPHA,T,H,PARAMI,PARAMB,RHO)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the reactivity during a Runge-Kutta time step taking into
* account feedback effects.
*
*Copyright:
* Copyright (C) 2017 Ecole Polytechnique de Montreal
* This library is free software; you can redIribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version
*
*Author(s):
* A. Hebert
*
*Parameters: input
* IPMAP pointer to the point kinetic directory
* NALPHA number of feedback parameters
* T time at beggining of step
* H time-step duration
* PARAMI initial values of the global parameters corresponding to
* RHO=0
* PARAMB values of global parameters at beginning of stage
*
*Parameters: ouput
* PARAMB values of global parameters at end of Runge-Kutta time step
* RHO reactivity during Runge-Kutta time step
*
*-----------------------------------------------------------------------
*
USE GANLIB
*----
* Subroutine arguments
*----
TYPE(C_PTR) IPMAP
INTEGER NALPHA
REAL PARAMI(NALPHA),PARAMB(NALPHA)
DOUBLE PRECISION T,H,RHO(3)
*----
* Local variables
*----
TYPE(C_PTR) JPPAR,KPPAR
TYPE(C_PTR) X_PTR,Y_PTR
DOUBLE PRECISION TS(3),DSUM
LOGICAL LCUBIC
*----
* Allocatable arrays
*----
REAL, POINTER, DIMENSION(:) :: X,Y
REAL, ALLOCATABLE, DIMENSION(:) :: TERP,GAR
REAL, ALLOCATABLE, DIMENSION(:,:) :: PARAM
*----
* Compute the values of the global parameters during the time step
*----
TS(1)=T
TS(2)=T+H/2.0D0
TS(3)=T+H
ALLOCATE(PARAM(3,NALPHA))
DO IAL=1,NALPHA
PARAM(:3,IAL)=PARAMB(IAL)
ENDDO
JPPAR=LCMGID(IPMAP,'ALPHA')
DO IAL=1,NALPHA
KPPAR=LCMGIL(JPPAR,IAL)
CALL LCMLEN(KPPAR,'TIME-LAW-T',NXY,ITYLCM)
IF(NXY.NE.0) THEN
ALLOCATE(TERP(NXY))
CALL LCMGPD(KPPAR,'TIME-LAW-T',X_PTR)
CALL C_F_POINTER(X_PTR,X,(/ NXY /))
CALL LCMGPD(KPPAR,'TIME-LAW-P',Y_PTR)
CALL C_F_POINTER(Y_PTR,Y,(/ NXY /))
CALL LCMGET(KPPAR,'TIME-LAW-I',LCUBIC)
DO I=1,3
CALL ALTERP(LCUBIC,NXY,X,REAL(TS(I)),.FALSE.,TERP)
DSUM=0.0D0
DO J=1,NXY
DSUM=DSUM+TERP(J)*Y(J)
ENDDO
PARAM(I,IAL)=REAL(DSUM)
ENDDO
PARAMB(IAL)=PARAM(3,IAL)
DEALLOCATE(TERP)
ENDIF
ENDDO
*----
* Compute the reactivity
*----
RHO(:3)=0.0D0
JPPAR=LCMGID(IPMAP,'ALPHA')
DO IAL=1,NALPHA
KPPAR=LCMGIL(JPPAR,IAL)
CALL LCMLEN(KPPAR,'ALPHA-LAW-P',NXY,ITYLCM)
IF(NXY.NE.0) THEN
ALLOCATE(TERP(NXY),GAR(NXY))
CALL LCMGPD(KPPAR,'ALPHA-LAW-P',X_PTR)
CALL C_F_POINTER(X_PTR,X,(/ NXY /))
CALL LCMGPD(KPPAR,'ALPHA-LAW-R',Y_PTR)
CALL C_F_POINTER(Y_PTR,Y,(/ NXY /))
CALL LCMGET(KPPAR,'ALPHA-LAW-T',ITYPE)
CALL LCMGET(KPPAR,'ALPHA-LAW-I',LCUBIC)
DO I=1,3
IF(ITYPE.EQ.1) THEN
CALL ALTERP(LCUBIC,NXY,X,PARAM(I,IAL),.FALSE.,TERP)
CALL ALTERP(LCUBIC,NXY,X,PARAMI(IAL),.FALSE.,GAR)
DSUM=0.0D0
DO J=1,NXY
DSUM=DSUM+(TERP(J)-GAR(J))*Y(J)
ENDDO
ELSE IF((ITYPE.EQ.2).AND.(PARAMI(IAL).LT.PARAM(I,IAL))) THEN
CALL ALTERI(LCUBIC,NXY,X,PARAMI(IAL),PARAM(I,IAL),TERP)
DSUM=0.0D0
DO J=1,NXY
DSUM=DSUM+TERP(J)*Y(J)
ENDDO
ELSE IF((ITYPE.EQ.2).AND.(PARAMI(IAL).GT.PARAM(I,IAL))) THEN
CALL ALTERI(LCUBIC,NXY,X,PARAM(I,IAL),PARAMI(IAL),TERP)
DSUM=0.0D0
DO J=1,NXY
DSUM=DSUM-TERP(J)*Y(J)
ENDDO
ELSE IF(ITYPE.EQ.3) THEN
DO J=1,NXY
GAR(J)=SQRT(X(J))
ENDDO
GAR1=SQRT(PARAMI(IAL))
GAR2=SQRT(PARAM(I,IAL))
IF(GAR1.LT.GAR2) THEN
CALL ALTERI(LCUBIC,NXY,GAR,GAR1,GAR2,TERP)
DSUM=0.0D0
DO J=1,NXY
DSUM=DSUM+TERP(J)*Y(J)
ENDDO
ELSE IF(GAR2.LT.GAR1) THEN
CALL ALTERI(LCUBIC,NXY,GAR,GAR2,GAR1,TERP)
DSUM=0.0D0
DO J=1,NXY
DSUM=DSUM-TERP(J)*Y(J)
ENDDO
ELSE
CYCLE
ENDIF
ELSE
CYCLE
ENDIF
RHO(I)=RHO(I)+DSUM
ENDDO
DEALLOCATE(GAR,TERP)
ENDIF
ENDDO
DEALLOCATE(PARAM)
RETURN
END
|
