CFPRIME DSTAB1                                                          DSTAB1      1
C***********************************************************************RIGHTE      1
C*                                                                     *RIGHTE      2
C*        Extended Transient-Midterm Stability Program (ETMSP)         *RIGHTE      3
C*        ====================================================         *RIGHTE      4
C*                                                                     *RIGHTE      5
C*       Copyright (C) (1993) Electric Power Research Institute        *RIGHTE      6
C*                                                                     *RIGHTE      7
C* EPRI reserves all rights in the software and associated             *RIGHTE      8
C* documentation as delivered. The software and its associated         *RIGHTE      9
C* documention or any portion thereof may not be reproduced in any form*RIGHTE     10
C* whatsoever without the consent of EPRI, except as explicity provided*RIGHTE     11
C* for by license. Licensee shall submit a detailed description to EPRI*RIGHTE     12
C* of any changes or modifications to the software within two (2)      *RIGHTE     13
C* months after such changes or modifications are made. At EPRI's sole *RIGHTE     14
C* discretion, the software, as changed or modified, shall be given a  *RIGHTE     15
C* new designation sufficiently different from its then current        *RIGHTE     16
C* designation as to prevent mistake,confusion, or deception as between*RIGHTE     17
C* the current software as changed or modified.                        *RIGHTE     18
C*                                                                     *RIGHTE     19
C* A license under EPRI's rights in the software can be obtained       *RIGHTE     20
C* directly from EPRI.                                                 *RIGHTE     21
C*                                                                     *RIGHTE     22
C* This program and supporting materials may be independently available*RIGHTE     23
C* from Electric Power Software centre (EPSC) for an appropriate       *RIGHTE     24
C* distribution fee. The EPSC is a service of                          *RIGHTE     25
C*                                                                     *RIGHTE     26
C*                    Power Computing Company                          *RIGHTE     27
C*                    1930 Hi Line                                     *RIGHTE     28
C*                    Dallas, Tx 75207                                 *RIGHTE     29
C*                                                                     *RIGHTE     30
C* NEITHER EPRI, ANY MEMBER OF EPRI, NOR ANY PERSON OR  ORGANIZATION   *RIGHTE     31
C* ACTING ON BEHAVE OF THEM:                                           *RIGHTE     32
C*                                                                     *RIGHTE     33
C* 1. MAKE ANY WARRANTY OR REPRESENTATION WHATSOEVER, EXPRESS OR       *RIGHTE     34
C*    IMPLIED, INCLUDING ANY WARRANTY OF MERCHANTABILITY OR FITNESS OR *RIGHTE     35
C*    ANY PURPOSE WITH RESPECT TO THE SOFTWARE OR ITS ASSOCIATED       *RIGHTE     36
C*    DOCUMENTATION; OR                                                *RIGHTE     37
C*                                                                     *RIGHTE     38
C* 2. ASSUMES ANY LIABILITY WHATSOEVER WITH RESPECT TO ANY USE OF THE  *RIGHTE     39
C*    SOFTWARE OR ITS ASSOCIATED DOCUMENTATION OR ANY PORTION THEREOF  *RIGHTE     40
C*    OR WITH RESPECT TO ANY DAMAGES WHICH MAY RESULT FROM SUCH USE.   *RIGHTE     41
C*                                                                     *RIGHTE     42
C* THIS NOTICE MAY NOT BE REMOVED FROM THE SOFTWARE BY ANY USER        *RIGHTE     43
C* THEREOF.                                                            *RIGHTE     44
C*                                                                     *RIGHTE     45
C***********************************************************************RIGHTE     46
      SUBROUTINE DSTAB1(NXS,NWSS)                                       DSTAB1      3
CD                                                                      DSTAB1      4
CD    TO CALCULATE THE DERIVATIVES OF THE STATES OF TYPES 1, 2 AND 3    DSTAB1      5
CD    STABILIZER, USED BY BOTH RUNGE-KUTTA AND TRAPEZOIDAL METHOD.      DSTAB1      6
CD                                                                      DSTAB1      7
CD    THE DATA STRUCTURES OF THESE TYPES OF STABILIZERS ARE GIVEN BELOW.DSTAB1      8
CD                                                                      DSTAB1      9
CD   4.1 Stabilizer type 1                                              DSTAB1     10
CD                                                                      DSTAB1     11
CD     Stabilizer data                                                  DSTAB1     12
CD                                                                      DSTAB1     13
CD     IDDAT(JS+1)    = BUSR         : sensing bus number               DSTAB1     14
CD     IDDAT(JS+2)    = TYPE         : type of input signal:            DSTAB1     15
CD                                     = 1 rotor speed                  DSTAB1     16
CD                                     = 2 acceleration power           DSTAB1     17
CD                                     = 3 bus frequency                DSTAB1     18
CD                                     = 4 electrical power             DSTAB1     19
CD                                     = 5 bus voltage magnitude        DSTAB1     20
CD                                     = 6 voltage deviation from steadyDSTAB1     21
CD                                        state value                   DSTAB1     22
CD     DDAT(JS+ 3)    = KS                                              DSTAB1     23
CD     DDAT(JS+ 4)    = T5                                              DSTAB1     24
CD     DDAT(JS+ 5)    = T1                                              DSTAB1     25
CD     DDAT(JS+ 6)    = T2                                              DSTAB1     26
CD     DDAT(JS+ 7)    = T3                                              DSTAB1     27
CD     DDAT(JS+ 8)    = T4                                              DSTAB1     28
CD     DDAT(JS+ 9)    = VSMAX                                           DSTAB1     29
CD     DDAT(JS+10)    = VSMIN                                           DSTAB1     30
CD     DDAT(JS+11)    = TF                                              DSTAB1     31
CD     DDAT(JS+12)    = RC                                              DSTAB1     32
CD     DDAT(JS+13)    = XC                                              DSTAB1     33
CD     DDAT(JS+14)    = PTHR                                            DSTAB1     34
CD                                                                      DSTAB1     35
CD     Stabilizer states                                                DSTAB1     36
CD                                                                      DSTAB1     37
CD     XSTATE(JSS+1)  = TF block if input is bus frequency              DSTAB1     38
CD                                   : optional                         DSTAB1     39
CD     XSTATE(JSS+2)  = T5 block     : compulsory                       DSTAB1     40
CD     XSTATE(JSS+3)  = T2 block     : optional                         DSTAB1     41
CD     XSTATE(JSS+4)  = T4 block     : optional                         DSTAB1     42
CD                                                                      DSTAB1     43
CD     Stabilizer Y-variables                                           DSTAB1     44
CD                                                                      DSTAB1     45
CD     VSS(JYS+1)     = stabilizer output                               DSTAB1     46
CD                                                                      DSTAB1     47
CD   4.2 Stabilizer type 2                                              DSTAB1     48
CD                                                                      DSTAB1     49
CD     Stabilizer data                                                  DSTAB1     50
CD                                                                      DSTAB1     51
CD     IDDAT(JS+1)    = BUSR         : sensing bus number               DSTAB1     52
CD     IDDAT(JS+2)    = TYPE         : type of input signal:            DSTAB1     53
CD                                     = 1 rotor speed                  DSTAB1     54
CD                                     = 2 acceleration power           DSTAB1     55
CD                                     = 3 bus frequency                DSTAB1     56
CD                                     = 4 electrical power             DSTAB1     57
CD                                     = 5 bus voltage magnitude        DSTAB1     58
CD                                     = 6 voltage deviation from steadyDSTAB1     59
CD                                        state value                   DSTAB1     60
CD     DDAT(JS+ 3)    = KS                                              DSTAB1     61
CD     DDAT(JS+ 4)    = T5                                              DSTAB1     62
CD     DDAT(JS+ 5)    = T1                                              DSTAB1     63
CD     DDAT(JS+ 6)    = T2                                              DSTAB1     64
CD     DDAT(JS+ 7)    = T3                                              DSTAB1     65
CD     DDAT(JS+ 8)    = T4                                              DSTAB1     66
CD     DDAT(JS+ 9)    = VSMAX                                           DSTAB1     67
CD     DDAT(JS+10)    = VSMIN                                           DSTAB1     68
CD     DDAT(JS+11)    = TF                                              DSTAB1     69
CD     DDAT(JS+12)    = RC                                              DSTAB1     70
CD     DDAT(JS+13)    = XC                                              DSTAB1     71
CD     DDAT(JS+14)    = PTHR                                            DSTAB1     72
CD     DDAT(JS+15)    = KANG                                            DSTAB1     73
CD     DDAT(JS+16)    = TANG                                            DSTAB1     74
CD     DDAT(JS+17)    = DWTHR                                           DSTAB1     75
CD     DDAT(JS+18)    = ETHR                                            DSTAB1     76
CD     DDAT(JS+19)    = VS3MX                                           DSTAB1     77
CD     DDAT(JS+20)    = VS3MN                                           DSTAB1     78
CD     DDAT(JS+21)    = VANGM                                           DSTAB1     79
CD     DDAT(JS+22)    = VANGMX                                          DSTAB1     80
CD                                                                      DSTAB1     81
CD     Stabilizer states                                                DSTAB1     82
CD                                                                      DSTAB1     83
CD     XSTATE(JSS+1)  = TF block if input is bus frequency              DSTAB1     84
CD                                   : optional                         DSTAB1     85
CD     XSTATE(JSS+2)  = T5 block     : compulsory                       DSTAB1     86
CD     XSTATE(JSS+3)  = T2 block     : compulsory                       DSTAB1     87
CD     XSTATE(JSS+4)  = T4 block     : optional                         DSTAB1     88
CD     XSTATE(JSS+5)  = TANG block        : compulsory                  DSTAB1     89
CD                                                                      DSTAB1     90
CD     Stabilizer Y-variables                                           DSTAB1     91
CD                                                                      DSTAB1     92
CD     VSS(JYS+1)     = stabilizer output                               DSTAB1     93
CD                                                                      DSTAB1     94
CD   4.3 Stabilizer type 3                                              DSTAB1     95
CD                                                                      DSTAB1     96
CD     Stabilizer data                                                  DSTAB1     97
CD                                                                      DSTAB1     98
CD     IDDAT(JS+1)    = BUSR         : sensing bus number               DSTAB1     99
CD     IDDAT(JS+2)    = TYPE         : type of input signal:            DSTAB1    100
CD                                     = 1 rotor speed                  DSTAB1    101
CD                                     = 2 acceleration power           DSTAB1    102
CD                                     = 3 bus frequency                DSTAB1    103
CD                                     = 4 electrical power             DSTAB1    104
CD                                     = 5 bus voltage magnitude        DSTAB1    105
CD                                     = 6 voltage deviation from steadyDSTAB1    106
CD                                        state value                   DSTAB1    107
CD     DDAT(JS+ 3)    = KS                                              DSTAB1    108
CD     DDAT(JS+ 4)    = T5                                              DSTAB1    109
CD     DDAT(JS+ 5)    = T1                                              DSTAB1    110
CD     DDAT(JS+ 6)    = T2                                              DSTAB1    111
CD     DDAT(JS+ 7)    = T3                                              DSTAB1    112
CD     DDAT(JS+ 8)    = T4                                              DSTAB1    113
CD     DDAT(JS+ 9)    = VSMAX                                           DSTAB1    114
CD     DDAT(JS+10)    = VSMIN                                           DSTAB1    115
CD     DDAT(JS+11)    = TF                                              DSTAB1    116
CD     DDAT(JS+12)    = RC                                              DSTAB1    117
CD     DDAT(JS+13)    = XC                                              DSTAB1    118
CD     DDAT(JS+14)    = PTHR                                            DSTAB1    119
CD     DDAT(JS+15)    = KANG                                            DSTAB1    120
CD     DDAT(JS+16)    = TANG                                            DSTAB1    121
CD     DDAT(JS+17)    = DWTHR                                           DSTAB1    122
CD     DDAT(JS+18)    = EIHR                                            DSTAB1    123
CD     DDAT(JS+19)    = VS3MX                                           DSTAB1    124
CD     DDAT(JS+20)    = VS3MN                                           DSTAB1    125
CD     DDAT(JS+21)    = VANGM                                           DSTAB1    126
CD     DDAT(JS+22)    = VANGMX                                          DSTAB1    127
CD                                                                      DSTAB1    128
CD     Stabilizer states                                                DSTAB1    129
CD                                                                      DSTAB1    130
CD     XSTATE(JSS+1)  = TF block if input is bus frequency              DSTAB1    131
CD                                   : optional                         DSTAB1    132
CD     XSTATE(JSS+2)  = T5 block     : compulsory                       DSTAB1    133
CD     XSTATE(JSS+3)  = T2 block     : compulsory                       DSTAB1    134
CD     XSTATE(JSS+4)  = T4 block     : optional                         DSTAB1    135
CD     XSTATE(JSS+5)  = TANG block        : compulsory                  DSTAB1    136
CD                                                                      DSTAB1    137
CD     Stabilizer Y-variables                                           DSTAB1    138
CD                                                                      DSTAB1    139
CD     VSS(JYS+1)     = stabilizer output                               DSTAB1    140
CD                                                                      DSTAB1    141
      INCLUDE 'PARX.INS'                                                DSTAB1    142
      INCLUDE 'CNMDAT.INS'                                              DSTAB1    143
      INCLUDE 'CNEXC.INS'                                               DSTAB1    144
      INCLUDE 'CNXSTA.INS'                                              DSTAB1    145
      INCLUDE 'CNQSTA.INS'                                              DSTAB1    146
      INCLUDE 'CNYVAR.INS'                                              DSTAB1    147
      INCLUDE 'CNSGAL.INS'                                              DSTAB1    148
      INCLUDE 'CNFLAG.INS'                                              DSTAB1    149
      INCLUDE 'CNUNIT.INS'                                              DSTAB1    150
      INCLUDE 'CNB9.INS'                                                DSTAB1    151
C                                                                       DSTAB1    152
*IF DEBUG                                                               DSTAB1    153
      LDEB=(JDEB(84).NE.0)                                              DSTAB1    154
*ENDIF                                                                  DSTAB1    155
      ICOUNT=ICOUNT+1                                                   DSTAB1    156
C                                                                       DSTAB1    157
C                                                                       DSTAB1    158
*IF DEBUG                                                               DSTAB1    159
      IF((LDEB).AND.ICOUNT.EQ.4)WRITE(IUD6,1100) JS,LSTATE,             DSTAB1    160
     &  (IDDAT(JS+LK),LK=1,2),(DDAT(JS+LK),LK=3,20)                     DSTAB1    161
*ENDIF                                                                  DSTAB1    162
C                                                                       DSTAB1    163
C   ISTABN = 1                                                          DSTAB1    164
C                                                                       DSTAB1    165
      IF (ISTABN.EQ.1) THEN                                             DSTAB1    166
C                                                                       DSTAB1    167
         KSTATE=KSTATE+1                                                DSTAB1    168
         DSTATE(LSTATE+KSTATE)=(DDAT(JS+3)*VSI1-XSTATE(JSTATE+KSTATE))  DSTAB1    169
     *                         /DDAT(JS+4)                              DSTAB1    170
         VSX=DDAT(JS+4)*DSTATE(LSTATE+KSTATE)                           DSTAB1    171
C                                                                       DSTAB1    172
         IF (IDDAT(JS+6).NE.0) THEN                                     DSTAB1    173
           KSTATE=KSTATE+1                                              DSTAB1    174
           DSTATE(LSTATE+KSTATE)=(VSX-XSTATE(JSTATE+KSTATE))/DDAT(JS+6) DSTAB1    175
           VSX=XSTATE(JSTATE+KSTATE)+DDAT(JS+5)*DSTATE(LSTATE+KSTATE)   DSTAB1    176
*IF DEBUG                                                               DSTAB1    177
          IF((LDEB).AND.ICOUNT.EQ.4)WRITE(IUD6,1250)JSTATE,KSTATE,VSI1, DSTAB1    178
     &     DSTATE(LSTATE+KSTATE),XSTATE(JSTATE+KSTATE),VSX              DSTAB1    179
 1250     FORMAT('TYPE 1 BLOCK #2 : JSTATE,KSTATE,VSI1,XDOT,X,VSX ',    DSTAB1    180
     &    2I5,4F12.5)                                                   DSTAB1    181
*ENDIF                                                                  DSTAB1    182
         ENDIF                                                          DSTAB1    183
         IF (IDDAT(JS+8).NE.0) THEN                                     DSTAB1    184
           KSTATE=KSTATE+1                                              DSTAB1    185
           DSTATE(LSTATE+KSTATE)=(VSX-XSTATE(JSTATE+KSTATE))/DDAT(JS+8) DSTAB1    186
           VSX=XSTATE(JSTATE+KSTATE)+DDAT(JS+7)*DSTATE(LSTATE+KSTATE)   DSTAB1    187
*IF DEBUG                                                               DSTAB1    188
          IF((LDEB).AND.ICOUNT.EQ.4)WRITE(IUD6,1300)JSTATE,KSTATE,VSI1, DSTAB1    189
     &     DSTATE(LSTATE+KSTATE),XSTATE(JSTATE+KSTATE),VSX              DSTAB1    190
 1300     FORMAT('TYPE 1 BLOCK #3 : JSTATE,KSTATE,VSI1,XDOT,X,VSX ',    DSTAB1    191
     &    2I5,4F12.5)                                                   DSTAB1    192
*ENDIF                                                                  DSTAB1    193
         ENDIF                                                          DSTAB1    194
C                                                                       DSTAB1    195
          IF(VSX.GT.DDAT(JS+9)) VSX=DDAT(JS+9)                          DSTAB1    196
          IF(VSX.LT.DDAT(JS+10)) VSX=DDAT(JS+10)                        DSTAB1    197
          VSS(JYS+1)=VSX                                                DSTAB1    198
C                                                                       DSTAB1    199
        ELSE                                                            DSTAB1    200
C                                                                       DSTAB1    201
C   ISTAB=2,3                                                           DSTAB1    202
C                                                                       DSTAB1    203
          KSTATE=KSTATE+1                                               DSTAB1    204
          DSTATE(LSTATE+KSTATE)=(DDAT(JS+3)*VSI1-XSTATE(JSTATE+KSTATE)) DSTAB1    205
     *                          /DDAT(JS+4)                             DSTAB1    206
          VS1=DDAT(JS+4)*DSTATE(LSTATE+KSTATE)                          DSTAB1    207
C                                                                       DSTAB1    208
          KSTATE=KSTATE+1                                               DSTAB1    209
          DSTATE(LSTATE+KSTATE)=(VS1-XSTATE(JSTATE+KSTATE))/DDAT(JS+6)  DSTAB1    210
          VSX=DDAT(JS+5)*DSTATE(LSTATE+KSTATE)+XSTATE(JSTATE+KSTATE)    DSTAB1    211
C                                                                       DSTAB1    212
C        SECOND PHASE LEAD - OPTIONAL                                   DSTAB1    213
C                                                                       DSTAB1    214
          IF (IDDAT(JS+8).NE.0)THEN                                     DSTAB1    215
          KSTATE=KSTATE+1                                               DSTAB1    216
          DSTATE(LSTATE+KSTATE)=(VSX-XSTATE(JSTATE+KSTATE))/DDAT(JS+8)  DSTAB1    217
          VSX=DDAT(JS+7)*DSTATE(LSTATE+KSTATE)+XSTATE(JSTATE+KSTATE)    DSTAB1    218
          ENDIF                                                         DSTAB1    219
C                                                                       DSTAB1    220
C       S1-ANG BLOCK-S2 PATH                                            DSTAB1    221
C       IDDAT(JS+23) IS THE STATUS OF S1 : = 0  ORIGINALLY OPEN         DSTAB1    222
C                                          = 1  CLOSE                   DSTAB1    223
C                                          = 2  REOPEN AND STAY OPENED  DSTAB1    224
C                                                                       DSTAB1    225
          XANG=0.0                                                      DSTAB1    226
          IF (IDDAT(JS+23).NE.2)THEN                                    DSTAB1    227
             IF (ET.LT.DDAT(JS+18).AND.XVARB.GE.REFMX.AND.VSI1.GE.      DSTAB1    228
     +           DDAT(JS+17)) IDDAT(JS+23)=1                            DSTAB1    229
             IF (IDDAT(JS+23).EQ.1) XANG=VSI1*DDAT(JS+15)               DSTAB1    230
          ENDIF                                                         DSTAB1    231
          KSTATE=KSTATE+1                                               DSTAB1    232
          DSTATE(LSTATE+KSTATE)=(XANG-XSTATE(JSTATE+KSTATE))/DDAT(JS+16)DSTAB1    233
          VANG=XSTATE(JSTATE+KSTATE)                                    DSTAB1    234
          IF (VANG.GE.DDAT(JS+21)) THEN                                 DSTAB1    235
             VANG=DDAT(JS+21)                                           DSTAB1    236
             IF (DSTATE(LSTATE+KSTATE).GT.0.0) THEN                     DSTAB1    237
                DSTATE(LSTATE+KSTATE)=0.0                               DSTAB1    238
                QSTATE(JSTATE+KSTATE)=0.0                               DSTAB1    239
             END IF                                                     DSTAB1    240
          END IF                                                        DSTAB1    241
          IF (ISTABN.EQ.3) THEN                                         DSTAB1    242
             IF (VSI1.LT.0.0) THEN                                      DSTAB1    243
                VANG=0.0                                                DSTAB1    244
                IDDAT(JS+24)=0                                          DSTAB1    245
             ELSE                                                       DSTAB1    246
                IDDAT(JS+24)=1                                          DSTAB1    247
             END IF                                                     DSTAB1    248
          END IF                                                        DSTAB1    249
C                                                                       DSTAB1    250
C      LIMITERS ON VANG,VSX; SUMIMG JUNCTION; LIMITER ON VSS            DSTAB1    251
C                                                                       DSTAB1    252
          IF (VANG.GT.DDAT(JS+22)) VANG=DDAT(JS+22)                     DSTAB1    253
          IF (VANG.LT.0.0) VANG=0.0                                     DSTAB1    254
          IF (VSX.GT.DDAT(JS+19)) VSX=DDAT(JS+19)                       DSTAB1    255
          IF (VSX.LT.DDAT(JS+20)) VSX=DDAT(JS+20)                       DSTAB1    256
          VSSUM=VSX+VANG                                                DSTAB1    257
          IF (VSSUM.GT.DDAT(JS+9)) VSSUM=DDAT(JS+9)                     DSTAB1    258
          IF (VSSUM.LT.DDAT(JS+10)) VSSUM=DDAT(JS+10)                   DSTAB1    259
C                                                                       DSTAB1    260
           VSS(JYS+1)=VSSUM                                             DSTAB1    261
C                                                                       DSTAB1    262
C     CHECK IF SW1 SHOULD REOPEN                                        DSTAB1    263
C                                                                       DSTAB1    264
          IF (LIR4.AND.IDDAT(JS+23).EQ.1) THEN                          DSTAB1    265
             IF(XVARB.LT.REFMX.OR.VSI1.LT.DDAT(JS+17)) IDDAT(JS+23)=2   DSTAB1    266
          ENDIF                                                         DSTAB1    267
        ENDIF                                                           DSTAB1    268
C                                                                       DSTAB1    269
      LSTATE=LSTATE+KSTATE                                              DSTAB1    270
      JSTATE=JSTATE+KSTATE                                              DSTAB1    271
      NWSS=JSTATE                                                       DSTAB1    272
C                                                                       DSTAB1    273
*IF DEBUG                                                               DSTAB1    274
      IF((LDEB).AND.ICOUNT.EQ.4)WRITE(IUD6,1100) JS,LSTATE,             DSTAB1    275
     &  (IDDAT(JS+LK),LK=1,2),(DDAT(JS+LK),LK=3,17)                     DSTAB1    276
*ENDIF                                                                  DSTAB1    277
       IF(ICOUNT.EQ.4)ICOUNT=0                                          DSTAB1    278
C                                                                       DSTAB1    279
      RETURN                                                            DSTAB1    280
C                                                                       DSTAB1    281
 1100 FORMAT(1X,'DB-DSTAB:JS,LSTATE,DDAT(1:17) ',2I5,2X,2I5,/,          DSTAB1    282
     & (8X,10G12.5))                                                    DSTAB1    283
C                                                                       DSTAB1    284
      END                                                               DSTAB1    285
