*DK ALFPAC                                                                ALFPAC.1
      SUBROUTINE LFK (N,M,CP,W)                                           ALFPAC.2
      INCLUDE 'param.h'                                                   ALFPAC.3
      SAVE                                                                ALFPAC.3
C     ****                                                                ALFPAC.4
C     ****     FILE ALFPAC DOCUMENTATION ,DEC 1976 , VERSION 0            ALFPAC.5
C     ****     PAUL N SWARZTRAUBER , NATIONAL CENTER FOR ATMOSPHERIC      ALFPAC.6
C     ****     RESEARCH BOULDER COLORADO 80303 , WHICH IS SPONSORED BY    ALFPAC.7
C     ****     THE NATIONAL SCIENCE FOUNDATION                            ALFPAC.8
C     ****                                                                ALFPAC.9
C     ****     THIS FILE CONTAINS THE DOCUMENTATION FOR FILE ALFPAC       ALFPAC.10
C     ****     WHICH CONTAINS SEVERAL PROGRAMS FOR COMPUTING THE          ALFPAC.11
C     ****     NORMALIZED ASSOCIATED LEGENDRE FUNCTIONS PBAR(N,M,THETA)   ALFPAC.12
C     ****     WHERE THETA IS COLATITUDE.  THEY ARE NORMALIZED SO THAT    ALFPAC.13
C     ****     THE INTEGRAL OF PBAR(N,M,THETA)**2*SIN(THETA) ON THE       ALFPAC.14
C     ****     INTERVAL THETA=0 TO THETA=PI IS EQUAL TO 1.  THE RELATION  ALFPAC.15
C     ****     BETWEEN PBAR(N,M,THETA) AND THE ASSOCIATED LEGENDRE        ALFPAC.16
C     ****     FUNCTION P(N,M,THETA) IS GIVEN BY -                        ALFPAC.17
C     ****                                                                ALFPAC.18
C     ****     PBAR(N,M,THETA)=SQRT((N+.5)*FACTORIAL(N-M)/FACTORIAL(N+M)) ALFPAC.19
C     ****                             *P(N,M,THETA)                      ALFPAC.20
C     ****                                                                ALFPAC.21
C     ****     FOR MORE DETAIL SEE BELOUSOV - TABLES OF NORMALIZED        ALFPAC.22
C     ****     ASSOCIATED LEGENDRE POLYNOMIALS, MACMILLAN,1962            ALFPAC.23
C     ****                                                                ALFPAC.24
C     ****     A BRIEF DISCRIPTION OF EACH SUBROUTINE FOLLOWS WITH THE    ALFPAC.25
C     ****     DETAILS OF THE ARGUMENTS GIVEN FURTHER BELOW               ALFPAC.26
C     ****                                                                ALFPAC.27
C     ****     SUBROUTINE LFK                                             ALFPAC.28
C     ****                                                                ALFPAC.29
C     ****     SUBROUTINE LFK COMPUTES THE FOURIER COEFFICIENTS IN THE    ALFPAC.30
C     ****     TRIGNOMETRIC SERIES REPRESENTATION OF PBAR(N,M,THETA)      ALFPAC.31
C     ****                                                                ALFPAC.32
C     ****     SUBROUTINE LFP                                             ALFPAC.33
C     ****                                                                ALFPAC.34
C     ****     SUBROUTINE LFP USES THE COEFFICIENTS COMPUTED BY           ALFPAC.35
C     ****     SUBROUTINE LFK AND TABULATES PBAR(N,M,THETA) AT THE        ALFPAC.36
C     ****     COLATITUDES THETA(I)=(I-1)*PI/(L-1) I=1,...L.              ALFPAC.37
C     ****                                                                ALFPAC.38
C     ****     SUBROUTINE LFPT                                            ALFPAC.39
C     ****                                                                ALFPAC.40
C     ****     SUBROUTINE LFPT USES THE COEFFICIENTS COMPUTED BY          ALFPAC.41
C     ****     SUBROUTINE LFK TO COMPUTE PBAR(N,M,THETA) AT A FIXED       ALFPAC.42
C     ****     COLATITUDE THETA.                                          ALFPAC.43
C     ****                                                                ALFPAC.44
C     ****     SUBROUTINE LFMA                                            ALFPAC.45
C     ****                                                                ALFPAC.46
C     ****     GIVEN N,I AND L THEN SUBROUTINE LFM TABULATES              ALFPAC.47
C     ****     PBAR(N,M,THETA) FOR M=0,...,N AT COLATITUDE                ALFPAC.48
C     ****     THETA=(I-1)*PI/(L-1)                                       ALFPAC.49
C     ****                                                                ALFPAC.50
C     ****     SUBROUTINE LFMB                                            ALFPAC.51
C     ****                                                                ALFPAC.52
C     ****     LFMB IS THE SAME AS LFMA EXCEPT THAT IT WILL BE MORE       ALFPAC.53
C     ****     EFFICIENT IN CERTIAN APPLICATIONS. SEE THE DISCRIPTIONS    ALFPAC.54
C     ****     LFMA AND LFMB GIVEN BELOW                                  ALFPAC.55
C     ****                                                                ALFPAC.56
C     ****     SUBROUTINE LFNA                                            ALFPAC.57
C     ****                                                                ALFPAC.58
C     ****     GIVEN M, I AND L THEN SUBROUTINE LFNA TABULATES            ALFPAC.59
C     ****     PBAR(N,M,THETA) FOR N=M,...,L-1 AT THETA=(I-1)*PI/(L-1)    ALFPAC.60
C     ****                                                                ALFPAC.61
C     ****     SUBROUTINE LFNB                                            ALFPAC.62
C     ****                                                                ALFPAC.63
C     ****     LFNB IS THE SAME AS LFNA EXCEPT THAT IT WILL BE MORE       ALFPAC.64
C     ****     EFFICIENT IN CERTIAN APPLICATIONS. SEE THE DISCRIPTIONS    ALFPAC.65
C     ****     OF LFNA AND LFNB GIVEN BELOW                               ALFPAC.66
C     ****                                                                ALFPAC.67
C     ****     SUBROUTINE BELSOV                                          ALFPAC.68
C     ****                                                                ALFPAC.69
C     ****     GIVEN N,L AND I, SUBROUTINE BELSOV USES THE METHOD OF      ALFPAC.70
C     ****     BELOUSOV TO TABULATE PBAR(N,M,THETA) FOR M=0,...,N AT      ALFPAC.71
C     ****     THETA=(I-1)*PI/(L-1).                                      ALFPAC.72
C     ****                                                                ALFPAC.73
C     ****     THE ACCURACY OF THE SUBROUTINES WAS CHECKED USING DOUBLE   ALFPAC.74
C     ****     PRECISION VERSIONS OF EACH PROGRAM WHICH WITH THE          ALFPAC.75
C     ****     EXCEPTION OF BELSOV ARE ALSO ON THIS FILE. THESE PROGRAMS  ALFPAC.76
C     ****     CAN BE ACCESSED BY PRECEEDING THE NAME WITH D.  I.E. DLFK  ALFPAC.77
C     ****     CALLS THE DOUBLE VERSION OF LFK.                           ALFPAC.78
C     ****                                                                ALFPAC.79
C     ****                 TABLE OF MAXIMUM ABSOLUTE ERRORS               ALFPAC.80
C     ****                 ON THE CONTROL DATA 7600 FOR                   ALFPAC.81
C     ****                 N=0,...,L  M=0,...,N AND I=1,...,L             ALFPAC.82
C     ****                                                                ALFPAC.83
C     ****                 THE 7600 HAS A MACHINE EPSILON OF              ALFPAC.84
C     ****                 2.**(-48)=3.6D-15                              ALFPAC.85
C     ****                                                                ALFPAC.86
C     ****             LFK     LFMA    LFMB    LFNA   LFNB   BELSOV       ALFPAC.87
C     ****                                                                ALFPAC.88
C     ****     L= 51 7.5D-12 7.5D-12 2.8D-12 3.6D-11 2.2D-11 9.6D-13      ALFPAC.89
C     ****                                                                ALFPAC.90
C     ****     L=101 8.0D-11 2.2D-11 7.8D-12 2.9D-10 1.2D-10 3.4D-12      ALFPAC.91
C     ****                                                                ALFPAC.92
C     ****     L=201 6.3D-10 4.8D-11 1.6D-11 3.3D-09 1.2D-09 7.8D-12      ALFPAC.93
C     ****                                                                ALFPAC.94
C     ****                 TIMES IN MILLISECONDS FOR N=0,...,L            ALFPAC.95
C     ****                 M=0,...,N AND I=1,...,L ON THE                 ALFPAC.96
C     ****                 CONTROL DATA 7600.                             ALFPAC.97
C     ****                                                                ALFPAC.98
C     ****                                                                ALFPAC.99
C     ****             LFK     LFMA    LFMB    LFNA    LFNB   BELSOV      ALFPAC.100
C     ****                                                                ALFPAC.101
C     ****     L= 51   264      722    1279     870    1246    1752       ALFPAC.102
C     ****                                                                ALFPAC.103
C     ****     L=101  1787     4789    8118    5761    8837   12391       ALFPAC.104
C     ****                                                                ALFPAC.105
C     ****     L=201 12773    33965   56123   40993   65499   90998       ALFPAC.106
C     ****                                                                ALFPAC.107
C     ****     SUBROUTINE LFP TAKES APPROXIMATELY L TIMES SUBROUTINE LFK. ALFPAC.108
C     ****     WHEN THE FAST FOURIER TRANSFORM IS INCORPORATED INTO LFP   ALFPAC.109
C     ****     IT WILL TAKE ABOUT LOG2(L) TIMES SUBROUTINE LFK.           ALFPAC.110
C     ****                                                                ALFPAC.111
C     ****     SUBROUTINE LFK (N,M,CP,W)                                  ALFPAC.112
C     ****                                                                ALFPAC.113
C     ****     SUBROUTINE LFK COMPUTES THE FOURIER COEFFICIENTS IN THE    ALFPAC.114
C     ****     TRIGNOMETRIC SERIES REPRESENTATION OF PBAR(N,M,THETA).     ALFPAC.115
C     ****     LFK IS INTENDED TO BE FOLLOWED BY A CALL OF SUBROUTINE     ALFPAC.116
C     ****     LFP WHICH TABULATES THE FUNCTION. THE FORM OF THE          ALFPAC.117
C     ****     TRIGNOMETRIC SERIES IS GIVEN BELOW IN THE DOCUMENTATION    ALFPAC.118
C     ****     FOR SUBROUTINE LFP                                         ALFPAC.119
C     ****                                                                ALFPAC.120
C     ****     IN ADDITION TO COMPUTING THE ASSOCIATED LEGENDRE FUNCTIONS ALFPAC.121
C     ****     THESE COEFFICIENTS CAN BE USED TO OBTAIN ACCURATE VALUES   ALFPAC.122
C     ****     FOR EXPRESSIONS INVOLVING DERIVATIVES OR INTEGRALS OF THE  ALFPAC.123
C     ****     FUNCTIONS                                                  ALFPAC.124
C     ****                                                                ALFPAC.125
C     ****     INPUT PARAMETERS                                           ALFPAC.126
C     ****                                                                ALFPAC.127
C     ****     N      THE DEGREE (SUBSCRIPT) OF PBAR(N,M,THETA)           ALFPAC.128
C     ****                                                                ALFPAC.129
C     ****     M      THE ORDER (SUPERSCRIPT) OF PBAR(N,M,THETA)          ALFPAC.130
C     ****                                                                ALFPAC.131
C     ****     NOTE-  LFK CAN BE CALLED IN ANY ORDER OF N AND M, HOWEVER  ALFPAC.132
C     ****            IT IS SOMEWHAT MORE EFFICIENT IF M CHANGES MORE     ALFPAC.133
C     ****            RAPIDLY THAN N                                      ALFPAC.134
C     ****                                                                ALFPAC.135
C     ****     W      A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST     ALFPAC.136
C     ****            INT(2.5*N)+4. IF M=0 OR M=1 THIS ARRAY IS NOT USED  ALFPAC.137
C     ****            AND W MAY BE A DUMMY VARIABLE.                      ALFPAC.138
C     ****                                                                ALFPAC.139
C     ****     OUTPUT PARAMETERS                                          ALFPAC.140
C     ****                                                                ALFPAC.141
C     ****     CP     AN ARRAY OF LENGTH AT LEAST INT(N/2)+1 LOCATIONS    ALFPAC.142
C     ****            WHICH CONTAINS THE COEFFICIENTS                     ALFPAC.143
C     ****                                                                ALFPAC.144
C     ****     W      CONTAINS RESULTS WHICH MUST NOT BE DESTROYED IF LFK ALFPAC.145
C     ****            WILL BE CALLED AGAIN WITH N UNCHANGED. THIS DOES    ALFPAC.146
C     ****            NOT APPLY IF M=0 OR M=1 IN WHICH CASE W MAY BE A    ALFPAC.147
C     ****            DUMMY VARIABLE.                                     ALFPAC.148
C     ****                                                                ALFPAC.149
      DIMENSION CP(1),W(1)                                                ALFPAC.150
      IDW=N/2+1                                                           ALFPAC.151
      IW1=N+2                                                             ALFPAC.152
      IW2=IW1+IDW                                                         ALFPAC.153
      IW3=IW2+IDW                                                         ALFPAC.154
      CALL LFK1 (N,M,CP,W,W(IW1),W(IW2),W(IW3))                           ALFPAC.155
      RETURN                                                              ALFPAC.156
      END                                                                 ALFPAC.157
C                                                                         ALFPAC.158
*DK LFP                                                                   ALFPAC.159
      SUBROUTINE LFP (N,M,L,CP,P)                                         ALFPAC.160
      INCLUDE 'param.h'                                                   ALFPAC.161
      SAVE                ! This save is needed for computation           ALFPAC.162
C     ****                                                                ALFPAC.163
C     ****     SUBROUTINE LFP (N,M,L,CP,P)                                ALFPAC.164
C     ****                                                                ALFPAC.165
C     ****     SUBROUTINE LFP USES THE COEFFICIENTS CP GENERATED BY       ALFPAC.166
C     ****     SUBROUTINE LFK AND TABULATES PBAR(N,M,THETA) AT            ALFPAC.167
C     ****     THETA=(I-1)*PI/(L-1) FOR I=1,...,L.                        ALFPAC.168
C     ****                                                                ALFPAC.169
C     ****     INPUT PARAMETERS                                           ALFPAC.170
C     ****                                                                ALFPAC.171
C     ****     N      THE DEGREE (SUBSCRIPT) OF PBAR(N,M,THETA)           ALFPAC.172
C     ****                                                                ALFPAC.173
C     ****     M      THE ORDER (SUPERSCRIPT) OF PBAR(N,M,THETA)          ALFPAC.174
C     ****                                                                ALFPAC.175
C     ****     CP     AN ARRAY OF LENGTH AT LEAST INT(N/2)+1 LOCATIONS    ALFPAC.176
C     ****            WHICH CONTAINS THE COEFFICIENTS COMPUTED BY         ALFPAC.177
C     ****            SUBROUTINE LFK.                                     ALFPAC.178
C     ****                                                                ALFPAC.179
C     ****     L      THE NUMBER OF COLATITUDES THETA(I)=(I-1)*PI/(L-1)   ALFPAC.180
C     ****            FOR I=1,...,L.                                      ALFPAC.181
C     ****                                                                ALFPAC.182
C     ****     OUTPUT PARAMETER                                           ALFPAC.183
C     ****                                                                ALFPAC.184
C     ****     P      AN ARRAY OF LENGTH L WHICH CONTAINS A TABULATION OF ALFPAC.185
C     ****            PBAR(N,M,THETA).   THE TABULATION IS COMPUTED FROM  ALFPAC.186
C     ****            ONE OF THE FOLLOWING SERIES DEPENDING ON THE PARITY ALFPAC.187
C     ****            OF N AND M                                          ALFPAC.188
C     ****                                                                ALFPAC.189
C     ****        FOR N EVEN,M EVEN AND FOR I=1,...L                      ALFPAC.190
C     ****                                                                ALFPAC.191
C     ****            P(I)=.5*CP(1) PLUS THE SUM FROM K=1 TO K=N/2 OF     ALFPAC.192
C     ****                 CP(K)*COS(2*K*THETA(I))                        ALFPAC.193
C     ****                                                                ALFPAC.194
C     ****        FOR N EVEN,M ODD AND FOR I=1,...L                       ALFPAC.195
C     ****                                                                ALFPAC.196
C     ****            P(I)=THE SUM FROM K=1 TO K=N/2 OF                   ALFPAC.197
C     ****                                                                ALFPAC.198
C     ****                 CP(K)*SIN(2*K*THETA(I))                        ALFPAC.199
C     ****                                                                ALFPAC.200
C     ****        FOR N ODD,M EVEN AND FOR I=1,...L                       ALFPAC.201
C     ****                                                                ALFPAC.202
C     ****            P(I)=THE SUM FROM K=1 TO K=(N+1)/2 OF               ALFPAC.203
C     ****                                                                ALFPAC.204
C     ****                 CP(K)*COS((2*K-1)*THETA(I))                    ALFPAC.205
C     ****                                                                ALFPAC.206
C     ****        FOR N ODD,M ODD AND FOR I=1,...L                        ALFPAC.207
C     ****                                                                ALFPAC.208
C     ****            P(I)=THE SUM FROM K=1 TO K=(N+1)/2 OF               ALFPAC.209
C     ****                                                                ALFPAC.210
C     ****                 CP(K)*SIN((2*K-1)*THETA(I))                    ALFPAC.211
C     ****                                                                ALFPAC.212
      DIMENSION CP(1),P(1)                                                ALFPAC.213
      IF(N)10,10,40                                                       ALFPAC.214
   10 IF(M)20,20,40                                                       ALFPAC.215
   20 DO 30 I=1,L                                                         ALFPAC.216
       P(I)=.707106781186548                                              ALFPAC.217
   30 CONTINUE                                                            ALFPAC.218
      GO TO 240                                                           ALFPAC.219
   40 LS2=(L+1)/2                                                         ALFPAC.220
      LM1=L-1                                                             ALFPAC.221
      NP1=N+1                                                             ALFPAC.222
      PI=3.14159265358979                                                 ALFPAC.223
      DT=PI/LM1                                                           ALFPAC.224
      NMOD=MOD(N,2)                                                       ALFPAC.225
      MMOD=MOD(M,2)                                                       ALFPAC.226
      IF(NMOD)50,50,120                                                   ALFPAC.227
   50 IF(MMOD)60,60,90                                                    ALFPAC.228
   60 KDO=N/2+1                                                           ALFPAC.229
      DO 80 I=1,LS2                                                       ALFPAC.230
       TH=DBLE((I-1))*DT                                                  ALFPAC.231
       CDT=DCOS(TH+TH)                                                    ALFPAC.232
       SDT=DSIN(TH+TH)                                                    ALFPAC.233
       CT=1.                                                              ALFPAC.234
       ST=0.                                                              ALFPAC.235
       SUM=.5*CP(1)                                                       ALFPAC.236
       DO 70 KP1=2,KDO                                                    ALFPAC.237
C       ****                                                              ALFPAC.238
C       ****     K=KP1-1                                                  ALFPAC.239
C       ****                                                              ALFPAC.240
        CTH=CDT*CT-SDT*ST                                                 ALFPAC.241
        ST=SDT*CT+CDT*ST                                                  ALFPAC.242
        CT=CTH                                                            ALFPAC.243
        SUM=SUM+CP(KP1)*CT                                                ALFPAC.244
C       ****                                                              ALFPAC.245
C       ****     SUM=SUM+CP(KP1)*COS((K+K)*TH)                            ALFPAC.246
C       ****                                                              ALFPAC.247
   70  CONTINUE                                                           ALFPAC.248
       P(I)=SUM                                                           ALFPAC.249
   80 CONTINUE                                                            ALFPAC.250
      GO TO 190                                                           ALFPAC.251
   90 KDO=N/2                                                             ALFPAC.252
      DO 110 I=1,LS2                                                      ALFPAC.253
       TH=DBLE((I-1))*DT                                                  ALFPAC.254
       CDT=DCOS(TH+TH)                                                    ALFPAC.255
       SDT=DSIN(TH+TH)                                                    ALFPAC.256
       CT=1.                                                              ALFPAC.257
       ST=0.                                                              ALFPAC.258
       SUM=0.                                                             ALFPAC.259
       DO 100 K=1,KDO                                                     ALFPAC.260
       CTH=CDT*CT-SDT*ST                                                  ALFPAC.261
       ST=SDT*CT+CDT*ST                                                   ALFPAC.262
       CT=CTH                                                             ALFPAC.263
       SUM=SUM+CP(K)*ST                                                   ALFPAC.264
C      ****                                                               ALFPAC.265
C      ****     SUM=SUM+CP(K)*SIN((K+K)*TH)                               ALFPAC.266
C      ****                                                               ALFPAC.267
  100  CONTINUE                                                           ALFPAC.268
       P(I)=SUM                                                           ALFPAC.269
  110 CONTINUE                                                            ALFPAC.270
      GO TO 190                                                           ALFPAC.271
  120 KDO=(N+1)/2                                                         ALFPAC.272
      IF(MMOD)130,130,160                                                 ALFPAC.273
  130 DO 150 I=1,LS2                                                      ALFPAC.274
       TH=DBLE((I-1))*DT                                                  ALFPAC.275
       CDT=DCOS(TH+TH)                                                    ALFPAC.276
       SDT=DSIN(TH+TH)                                                    ALFPAC.277
       CT=DCOS(TH)                                                        ALFPAC.278
       ST=-DSIN(TH)                                                       ALFPAC.279
       SUM=0.                                                             ALFPAC.280
       DO 140 K=1,KDO                                                     ALFPAC.281
        CTH=CDT*CT-SDT*ST                                                 ALFPAC.282
        ST=SDT*CT+CDT*ST                                                  ALFPAC.283
        CT=CTH                                                            ALFPAC.284
        SUM=SUM+CP(K)*CT                                                  ALFPAC.285
C       ****                                                              ALFPAC.286
C       ****     SUM=SUM+CP(K)*COS((K+K-1)*TH)                            ALFPAC.287
C       ****                                                              ALFPAC.288
  140  CONTINUE                                                           ALFPAC.289
       P(I)=SUM                                                           ALFPAC.290
  150 CONTINUE                                                            ALFPAC.291
      GO TO 190                                                           ALFPAC.292
  160 DO 180 I=1,LS2                                                      ALFPAC.293
      TH=DBLE((I-1))*DT                                                   ALFPAC.294
       CDT=DCOS(TH+TH)                                                    ALFPAC.295
       SDT=DSIN(TH+TH)                                                    ALFPAC.296
       CT=DCOS(TH)                                                        ALFPAC.297
       ST=-DSIN(TH)                                                       ALFPAC.298
       SUM=0.                                                             ALFPAC.299
       DO 170 K=1,KDO                                                     ALFPAC.300
        CTH=CDT*CT-SDT*ST                                                 ALFPAC.301
        ST=SDT*CT+CDT*ST                                                  ALFPAC.302
        CT=CTH                                                            ALFPAC.303
        SUM=SUM+CP(K)*ST                                                  ALFPAC.304
C       ****                                                              ALFPAC.305
C       ****     SUM=SUM+CP(K)*SIN((K+K-1)*TH)                            ALFPAC.306
C       ****                                                              ALFPAC.307
  170  CONTINUE                                                           ALFPAC.308
       P(I)=SUM                                                           ALFPAC.309
  180 CONTINUE                                                            ALFPAC.310
  190 IF(MOD(N-M,2))220,200,220                                           ALFPAC.311
  200 DO 210 I=1,LS2                                                      ALFPAC.312
       LMI=L-I                                                            ALFPAC.313
       P(LMI+1)=P(I)                                                      ALFPAC.314
  210 CONTINUE                                                            ALFPAC.315
      GO TO 240                                                           ALFPAC.316
  220 DO 230 I=1,LS2                                                      ALFPAC.317
       LMI=L-I                                                            ALFPAC.318
       P(LMI+1)=-P(I)                                                     ALFPAC.319
  230 CONTINUE                                                            ALFPAC.320
  240 RETURN                                                              ALFPAC.321
      END                                                                 ALFPAC.322
C                                                                         ALFPAC.323
*DK LFK1                                                                  ALFPAC.324
      SUBROUTINE LFK1 (N,M,CP,D,A,B,C)                                    ALFPAC.325
      INCLUDE 'param.h'                                                   ALFPAC.326
      SAVE                                                                ALFPAC.327
      DIMENSION CP(1),D(1),A(1),B(1),C(1)                                 ALFPAC.328
      DATA NLAST /-1/                                                     ALFPAC.329
      DATA IINN/2/                                                        ALFPAC.330
      DATA IINM/3/                                                        ALFPAC.331
      IF(M-1)10,10,130                                                    ALFPAC.332
C     ****                                                                ALFPAC.333
C     ****     COEFFICIENTS FOR M=0                                       ALFPAC.334
C     ****                                                                ALFPAC.335
   10 NMOD=MOD(N,2)                                                       ALFPAC.336
      FN=DBLE(N)                                                          ALFPAC.337
      IF(NMOD)20,20,30                                                    ALFPAC.338
   20 K=-1                                                                ALFPAC.339
      CP(1)=2.*DSQRT(FN+.5)*QP(N-1)**2                                    ALFPAC.340
      GO TO  40                                                           ALFPAC.341
   30 K=0                                                                 ALFPAC.342
      CP(1)=2*DSQRT(FN+.5d0)*(FN+1)/FN*QP(N)**2                           ALFPAC.343
   40 KS=1                                                                ALFPAC.344
      FK=K                                                                ALFPAC.345
      FNPK=FN+FK                                                          ALFPAC.346
      FNMK=FN-FK                                                          ALFPAC.347
   50 K=K+2                                                               ALFPAC.348
      IF(N-K)70,70,60                                                     ALFPAC.349
   60 KS=KS+1                                                             ALFPAC.350
      FNPK=FNPK+2.                                                        ALFPAC.351
      FNMK=FNMK-2.                                                        ALFPAC.352
      CP(KS)=FNPK*(FNMK+1.)/FNMK*CP(KS-1)/(FNPK+1.)                       ALFPAC.353
      GO TO 50                                                            ALFPAC.354
   70 IF(M)340,340, 80                                                    ALFPAC.355
C     ****                                                                ALFPAC.356
C     ****     COEFFICIENTS FOR M=1                                       ALFPAC.357
C     ****                                                                ALFPAC.358
   80 ENM=1./DSQRT(FN*(FN+1.))                                            ALFPAC.359
      IF(NMOD)90,90,110                                                   ALFPAC.360
   90 KSM=KS-1                                                            ALFPAC.361
      TK=0.                                                               ALFPAC.362
      DO 100 K=1,KSM                                                      ALFPAC.363
       TK=TK+2.                                                           ALFPAC.364
       CP(K)=ENM*TK*CP(K+1)                                               ALFPAC.365
  100 CONTINUE                                                            ALFPAC.366
      GO TO 340                                                           ALFPAC.367
  110 TK=-1.                                                              ALFPAC.368
      DO 120 K=1,KS                                                       ALFPAC.369
       TK=TK+2.                                                           ALFPAC.370
       CP(K)=ENM*TK*CP(K)                                                 ALFPAC.371
  120 CONTINUE                                                            ALFPAC.372
      GO TO 340                                                           ALFPAC.373
  130 NP1=N+1                                                             ALFPAC.374
      MSQ=M*M                                                             ALFPAC.375
      TMSQ=MSQ+MSQ                                                        ALFPAC.376
      MMOD=MOD(M,2)                                                       ALFPAC.377
      IF(N-NLAST)140,250,140                                              ALFPAC.378
  140 NLAST=N                                                             ALFPAC.379
      NNP1=N*(N+1)                                                        ALFPAC.380
      NMOD=MOD(N,2)                                                       ALFPAC.381
      FN=DBLE(N)                                                          ALFPAC.382
      FNNP1=FN*(FN+1.)                                                    ALFPAC.383
      FNNP12=FNNP1-2.                                                     ALFPAC.384
      FNNP16=FNNP1-6.                                                     ALFPAC.385
      HFNM2=FNNP12/2.                                                     ALFPAC.386
      TN=DBLE(N+N)                                                        ALFPAC.387
      IF(NMOD)200,150,200                                                 ALFPAC.388
C     ****                                                                ALFPAC.389
C     ****     STARTING COEFFICIENTS FOR N EVEN                           ALFPAC.390
C     ****                                                                ALFPAC.391
150   D(1)=2.*DSQRT(FN+.5)*QP(N-1)**2                                     ALFPAC.392
      IP=0                                                                ALFPAC.393
      MH=0                                                                ALFPAC.394
      FM=0.                                                               ALFPAC.395
      FNMM=FN-FM                                                          ALFPAC.396
      FNPM=FN+FM                                                          ALFPAC.397
  160 MH=MH+1                                                             ALFPAC.398
      IP=1-IP                                                             ALFPAC.399
      IF(N-MH)250,170,170                                                 ALFPAC.400
  170 FM=FM+1.                                                            ALFPAC.401
      FNMM=FNMM-1.                                                        ALFPAC.402
      FNPM=FNPM+1.                                                        ALFPAC.403
      ALPHA=DSQRT(FNPM*(FNMM+1.))                                         ALFPAC.404
      IF(IP)180,180,190                                                   ALFPAC.405
  180 D(MH+1)=HFNM2/(FM-1.)*D(MH)/ALPHA                                   ALFPAC.406
      GO TO 160                                                           ALFPAC.407
  190 D(MH+1)=ALPHA/HFNM2*FM*D(MH)                                        ALFPAC.408
      GO TO 160                                                           ALFPAC.409
C     ****                                                                ALFPAC.410
C     ****     STARTING COEFFICIENTS FOR N ODD                            ALFPAC.411
C     ****                                                                ALFPAC.412
  200 D(1)=2.*DSQRT(FN+.5)*(FN+1.)/FN*QP(N)**2                            ALFPAC.413
      IP=0                                                                ALFPAC.414
      MH=0                                                                ALFPAC.415
      FM=0.                                                               ALFPAC.416
      FNMM=FN-FM                                                          ALFPAC.417
      FNPM=FN+FM                                                          ALFPAC.418
  210 MH=MH+1                                                             ALFPAC.419
      IP=1-IP                                                             ALFPAC.420
      IF(N-MH)250,220,220                                                 ALFPAC.421
  220 FM=FM+1.                                                            ALFPAC.422
      FNMM=FNMM-1.                                                        ALFPAC.423
      FNPM=FNPM+1.                                                        ALFPAC.424
      ALPHA=DSQRT(FNPM*(FNMM+1.))                                         ALFPAC.425
      IF(IP)230,230,240                                                   ALFPAC.426
  230 D(MH+1)=ALPHA/(FM-1.)*D(MH)                                         ALFPAC.427
      GO TO 210                                                           ALFPAC.428
  240 D(MH+1)=FM/ALPHA*D(MH)                                              ALFPAC.429
      GO TO 210                                                           ALFPAC.430
  250 IF(NMOD)260,260,300                                                 ALFPAC.431
C     ****                                                                ALFPAC.432
C     ****     N IS EVEN                                                  ALFPAC.433
C     ****                                                                ALFPAC.434
  260 NH=N/2                                                              ALFPAC.435
      NHDO=NH+1                                                           ALFPAC.436
      TI=-2.                                                              ALFPAC.437
      DO 270 IP1=1,NHDO                                                   ALFPAC.438
       TI=TI+2.                                                           ALFPAC.439
       A(IP1)=(TI-1.)*(TI-2.)-FNNP1                                       ALFPAC.440
       B(IP1)=2.*(FNNP1-TI*TI-TMSQ)                                       ALFPAC.441
       C(IP1)=(TI+1.)*(TI+2.)-FNNP1                                       ALFPAC.442
  270 CONTINUE                                                            ALFPAC.443
      IF(MMOD)280,280,290                                                 ALFPAC.444
C     ****                                                                ALFPAC.445
C     ****     N EVEN, M EVEN                                             ALFPAC.446
C     ****                                                                ALFPAC.447
  280 B(1)=FNNP1-TMSQ                                                     ALFPAC.448
      C(1)=2.-FNNP1                                                       ALFPAC.449
      CP(1)=D(M+1)                                                        ALFPAC.450
      CP(IINN)=-A(IINN)*D(M+1)                                            ALFPAC.451
      CALL TRIVEC (NH,A(IINN),B(IINN),C(IINN),CP(IINN))                   ALFPAC.452
      RETURN                                                              ALFPAC.453
C     ****                                                                ALFPAC.454
C     *****     N EVEN, M ODD                                             ALFPAC.455
C     ****                                                                ALFPAC.456
  290 B(IINN)=2.*FNNP1-TMSQ-4.                                            ALFPAC.457
      C(IINN)=12.-FNNP1                                                   ALFPAC.458
      CP(1)=D(M+1)                                                        ALFPAC.459
      CP(IINN)=-A(IINM)*D(M+1)                                            ALFPAC.460
      CALL TRIVEC(NH-1,A(IINM),B(IINM),C(IINM),CP(IINN))                  ALFPAC.461
      RETURN                                                              ALFPAC.462
C     ****                                                                ALFPAC.463
C     ****     N ODD                                                      ALFPAC.464
C     ****                                                                ALFPAC.465
  300 NH=(N+1)/2                                                          ALFPAC.466
      TI=0.                                                               ALFPAC.467
      DO 310 I=1,NH                                                       ALFPAC.468
       TI=TI+2.                                                           ALFPAC.469
       A(I)=(TI-2.)*(TI-3.)-FNNP1                                         ALFPAC.470
       B(I)=2.*(FNNP1-(TI-1.)**2-TMSQ)                                    ALFPAC.471
       C(I)=TI*(TI+1.)-FNNP1                                              ALFPAC.472
  310 CONTINUE                                                            ALFPAC.473
      IF(MMOD)320,320,330                                                 ALFPAC.474
C     ****                                                                ALFPAC.475
C     ****     N ODD, M EVEN                                              ALFPAC.476
C     ****                                                                ALFPAC.477
  320 B(1)=FNNP1-TMSQ-TMSQ-2.                                             ALFPAC.478
      C(1)=6.-FNNP1                                                       ALFPAC.479
      CP(1)=D(M+1)                                                        ALFPAC.480
      CP(IINN)=-A(IINN)*D(M+1)                                            ALFPAC.481
      CALL TRIVEC(NH-1,A(IINN),B(IINN),C(IINN),CP(IINN))                  ALFPAC.482
      RETURN                                                              ALFPAC.483
C     ****                                                                ALFPAC.484
C     ****     N ODD, M ODD                                               ALFPAC.485
C     ****                                                                ALFPAC.486
  330 B(1)=3.*FNNP1-TMSQ-TMSQ-2.                                          ALFPAC.487
      C(1)=6.-FNNP1                                                       ALFPAC.488
      CP(1)=D(M+1)                                                        ALFPAC.489
      CP(IINN)=-A(IINN)*D(M+1)                                            ALFPAC.490
      CALL TRIVEC(NH-1,A(IINN),B(IINN),C(IINN),CP(IINN))                  ALFPAC.491
  340 RETURN                                                              ALFPAC.492
      END                                                                 ALFPAC.493
C                                                                         ALFPAC.494
*DK QP                                                                    ALFPAC.495
      FUNCTION QP(L)                                                      ALFPAC.496
      INCLUDE 'param.h'                                                   ALFPAC.497
      SAVE                                                                ALFPAC.498
      QP=1.                                                               ALFPAC.499
      FL=DBLE(L)                                                          ALFPAC.500
   10 IF(FL)30,30,20                                                      ALFPAC.501
   20 QP=FL*QP/(FL+1)                                                     ALFPAC.502
      FL=FL-2                                                             ALFPAC.503
      GO TO 10                                                            ALFPAC.504
   30 RETURN                                                              ALFPAC.505
      END                                                                 ALFPAC.506
C                                                                         ALFPAC.507
*DK                                                                       BLOCKD.1
      BLOCK DATA CONSTS                                                   BLOCKD.2
      INCLUDE 'param.h'                                                   BLOCKD.3
      INCLUDE 'const.h'                                                   BLOCKD.4
      INCLUDE 'start.h'                                                   BLOCKD.5
      DATA NFLDS/32/,NHIST/24/,ITAP/1/                                    BLOCKD.6
      DATA KPOLE/1,-1,-1,1,1,1,1,1,0,1,1,0,0,1,1,1,1,1,-1,-1,1,1,1,1,1,   BLOCKD.7
     11,1,0,0,0,0,0/                                                      BLOCKD.8
      DATA KFLDS/16*3,1,10*3,5*2/                                         BLOCKD.9
      DATA RMASS/32.,16.,28./                                             BLOCKD.10
      DATA KUT/1,2,3,5,6,7,9,10,11,13,14,15,17,10*17,17,15,14,13,11,      BLOCKD.11
     110,9,7,6,5,3,2,1/                                                   BLOCKD.12
      DATA NESTIN/11,12,13,14/,NESTOT/21,22,23,24/
      DATA NESTRD/'n1in','n2in','n3in','n4in'/
      DATA NESTWT/'n1ot','n2ot','n3ot','n4ot'/
      END                                                                 BLOCKD.13
C                                                                         BLOCKD.14
*DK BNDRY                                                                  BNDRY.1
      SUBROUTINE BNDRY                                                     BNDRY.2
C     ****                                                                 BNDRY.3
C     ****     SPECIFY THE LOWER BOUNDARY CONDITIONS FOR                   BNDRY.4
C     ****     NEUTRAL VELOCITY: U AND V, AND NEURAL TEMPERATURE           BNDRY.5
C     ****                                                                 BNDRY.6
      INCLUDE 'param.h'                                                    BNDRY.7
      INCLUDE 'const.h'                                                    BNDRY.8
      INCLUDE 'start.h'                                                    BNDRY.9
      INCLUDE 'bndry.h'                                                    BNDRY.10
      COMMON/TIDES/TIDE(10),TIDE2(2)                                       BNDRY.11
      DIMENSION B(5,19),RL(5),BHOUR(5)                                     BNDRY.12
      COMPLEX*16 ZEE(5),CL(5),EXPDLM                                       BNDRY.13
      DATA B/                                                              BNDRY.14
     2 0.969152, 0.0     , 0.216046, 0.0     , 0.093838,                   BNDRY.15
     3 0.0     , 0.909763, 0.0     , 0.342113, 0.0     ,                   BNDRY.16
     4-0.245226, 0.0     , 0.798445, 0.0     , 0.421218,                   BNDRY.17
     5 0.0     ,-0.408934, 0.0     , 0.645517, 0.0     ,                   BNDRY.18
     6 0.024633, 0.0     ,-0.543993, 0.0     , 0.464159,                   BNDRY.19
     7 0.0     , 0.071127, 0.0     ,-0.643189, 0.0     ,                   BNDRY.20
     8-0.001292, 0.0     , 0.139613, 0.0     ,-0.699495,                   BNDRY.21
     9 0.0     ,-0.006673, 0.0     , 0.225090, 0.0     ,                   BNDRY.22
     1 0.000042, 0.0     ,-0.019654, 0.0     , 0.320141,                   BNDRY.23
     1 0.0     , 0.000394, 0.0     ,-0.043345, 0.0     ,                   BNDRY.24
     2-0.000001, 0.0     , 0.001772, 0.0     ,-0.079831,                   BNDRY.25
     3 0.0     ,-0.000016, 0.0     , 0.005401, 0.0     ,                   BNDRY.26
     4 0.0     , 0.0     ,-0.000112, 0.0     , 0.012932,                   BNDRY.27
     5 0.0     , 0.0     , 0.0     ,-0.000476, 0.0     ,                   BNDRY.28
     6 0.0     , 0.0     , 0.000005, 0.0     ,-0.001490,                   BNDRY.29
     7 0.0     , 0.0     , 0.0     , 0.000031, 0.0     ,                   BNDRY.30
     8 0.0     , 0.0     , 0.0     , 0.0     , 0.000129,                   BNDRY.31
     9 0.0     , 0.0     , 0.0     ,-0.000002, 0.0     ,                   BNDRY.32
     1 0.0     , 0.0     , 0.0     , 0.0     ,-0.000009/                   BNDRY.33
      DATA RL/                                                             BNDRY.34
     1 7.8519D5, 3.6665D5, 2.1098D5, 1.3671D5, 0.9565D5/                   BNDRY.35
      DATA T2/184./,T1/181./                                               BNDRY.36
      DO I=1,5                                                             BNDRY.37
       BHOUR(I)=TIDE(I+5)                                                  BNDRY.38
      ENDDO                                                                BNDRY.39
      CI=(0.,1.)                                                           BNDRY.40
      IF(ITIDE.EQ.0)GO TO 13                                               BNDRY.41
      DO N=1,5                                                             BNDRY.42
       ZEE(N)=TIDE(N)*CDEXP(CI*C(110)*BHOUR(N)/6.)                         BNDRY.43
       CL(N)=CDSQRT(DCMPLX(C(57)/(C(24)*C(54)*RL(N))*                      BNDRY.44
     1       (T1*2./7.+(T2-T1)/C(3))-.25))-.5*CI                           BNDRY.45
      ENDDO                                                                BNDRY.46
      JM=2*JMAX+1                                                          BNDRY.47
C     ****                                                                 BNDRY.48
C     ****     SET UP HOUGH FUNCTIONS                                      BNDRY.49
C     ****                                                                 BNDRY.50
      DO L=2,20                                                            BNDRY.51
       LM1=L-1                                                             BNDRY.52
       DO M=2,3                                                            BNDRY.53
        MM1=M-1                                                            BNDRY.54
        CALL LFK(L,M,CP,W)                                                 BNDRY.55
        CALL LFP(L,M,JM,CP,PT)                                             BNDRY.56
        DO J=1,JMAX                                                        BNDRY.57
         P(J,LM1,MM1)=PT(2*(JMAX+1-J))                                     BNDRY.58
        ENDDO                                                              BNDRY.59
       ENDDO                                                               BNDRY.60
       DO J=1,JMAX                                                         BNDRY.61
        P(J,LM1,2)=DSQRT(DBLE(L*(L+1)-6))*P(J,LM1,2)-2*TN(J)*P(J,LM1,1)    BNDRY.62
       ENDDO                                                               BNDRY.63
      ENDDO                                                                BNDRY.64
      DO L=1,5                                                             BNDRY.65
       DO LD=1,2                                                           BNDRY.66
        DO J=1,JMAX                                                        BNDRY.67
         CALL SDOT(HOUGH(J,L,LD),19,P(J,1,LD),JMAX,B(L,1),5)               BNDRY.68
        ENDDO                                                              BNDRY.69
       ENDDO                                                               BNDRY.70
      ENDDO                                                                BNDRY.71
C     ****                                                                 BNDRY.72
C     ****      GENERATE ZB, UB, VB AND TB                                 BNDRY.73
C     ****                                                                 BNDRY.74
      DO J=1,JMAX                                                          BNDRY.75
       TB(J)=0.                                                            BNDRY.76
       ZB(J)=0.                                                            BNDRY.77
       DZB(J)=0.                                                           BNDRY.78
      ENDDO                                                                BNDRY.79
      DO L=1,5                                                             BNDRY.80
       DO J=1,JMAX                                                         BNDRY.81
        ZB(J)=ZB(J)+ZEE(L)*HOUGH(J,L,1)                                    BNDRY.82
        DZB(J)=DZB(J)+ZEE(L)*HOUGH(J,L,2)                                  BNDRY.83
        TB(J)=TB(J)+CI*C(24)*C(54)/C(57)*ZEE(L)*CL(L)*HOUGH(J,L,1)         BNDRY.84	
       ENDDO                                                               BNDRY.85
      ENDDO                                                                BNDRY.86
      DO J=1,JMAX                                                          BNDRY.87
       UB(J)=C(23)*C(51)*(1.-(COR(J)/C(23))**2)                            BNDRY.88
       VB(J)=CI*C(54)*(DZB(J)-2.*COR(J)/(C(23)*CS(J))*ZB(J))/UB(J)         BNDRY.89
       UB(J)=C(54)*(COR(J)/C(23)*DZB(J)-2./CS(J)*ZB(J))/UB(J)              BNDRY.90
      ENDDO                                                                BNDRY.91
      GO TO 11                                                             BNDRY.92
   13 CONTINUE                                                             BNDRY.93
C     ****                                                                 BNDRY.94
C     ****      ZERO BOUNDARY CONDITION                                    BNDRY.95
C     ****                                                                 BNDRY.96
      DO J=1,JMAX                                                          BNDRY.97
       ZB(J)=0.                                                            BNDRY.98
       TB(J)=0.                                                            BNDRY.99
       UB(J)=0.                                                            BNDRY.100
       VB(J)=0.                                                            BNDRY.101
      ENDDO                                                                BNDRY.102
   11 CONTINUE                                                             BNDRY.103
C     ****                                                                 BNDRY.104
C     ****      CALCULATE LONGITUDINAL STRUCTURE                           BNDRY.105
C     ****                                                                 BNDRY.106
      RLAMDA=-2.*C(1)                                                      BNDRY.107
      BND(1)=CDEXP(CI*2.*RLAMDA)                                           BNDRY.108
      EXPDLM=CDEXP(CI*2.*C(1))                                             BNDRY.109
      DO I=2,IMAXP4                                                        BNDRY.110
       BND(I)=BND(I-1)*EXPDLM                                              BNDRY.112
      ENDDO                                                                BNDRY.113
      RETURN                                                               BNDRY.114
      END                                                                  BNDRY.115
C                                                                          BNDRY.116
*DK BNDRY2                                                                BNDRY2.1
      SUBROUTINE BNDRY2                                                   BNDRY2.2
C     ****                                                                BNDRY2.3
C     ****     TIDAL BOUNDARY CONDITION FOR DIURNAL GRAVITATIONAL MODE    BNDRY2.4
C     ****     (1,1)                                                      BNDRY2.5
C     ****                                                                BNDRY2.6
      INCLUDE 'param.h'                                                   BNDRY2.7
      INCLUDE 'const.h'                                                   BNDRY2.8
      INCLUDE 'start.h'                                                   BNDRY2.9
      INCLUDE 'bndry.h'                                                   BNDRY2.10
      COMMON/TIDES/TIDE(10),TIDE2(2)                                      BNDRY2.11
      DIMENSION B(1,19),RL (1),BHOUR2(1)                                  BNDRY2.12
      COMPLEX*16 ZEE(1),CL(1),EXPDLM                                      BNDRY2.13
      DATA B/                                                             BNDRY2.14
     2 0.282710,                                                          BNDRY2.15
     3 0.0     ,                                                          BNDRY2.16
     4-0.638229,                                                          BNDRY2.17
     5 0.0     ,                                                          BNDRY2.18
     6 0.620521,                                                          BNDRY2.19
     7 0.0     ,                                                          BNDRY2.20
     8-0.336408,                                                          BNDRY2.21
     9 0.0     ,                                                          BNDRY2.22
     1 0.117021,                                                          BNDRY2.23
     1 0.0     ,                                                          BNDRY2.24
     2-0.028332,                                                          BNDRY2.25
     3 0.0     ,                                                          BNDRY2.26
     4 0.005042,                                                          BNDRY2.27
     5 0.0     ,                                                          BNDRY2.28
     6-0.000686,                                                          BNDRY2.29
     7 0.0     ,                                                          BNDRY2.30
     8 0.000074,                                                          BNDRY2.31
     9 0.0     ,                                                          BNDRY2.32
     1-0.000006/                                                          BNDRY2.33
      DATA RL/                                                            BNDRY2.34
     1 0.6909D5/                                                          BNDRY2.35
      DATA T2/184./,T1/181./                                              BNDRY2.36
      DO I=1,1                                                            BNDRY2.37
       BHOUR2(I)=TIDE2(I+1)                                               BNDRY2.38
      ENDDO                                                               BNDRY2.39
      CI=(0.,1.)                                                          BNDRY2.40
      IF(ITIDE2.EQ.0)GO TO 13                                             BNDRY2.41
      DO N=1,1                                                            BNDRY2.42
       ZEE(N)=TIDE2(N)*CDEXP(CI*C(110)*BHOUR2(N)/12.)                     BNDRY2.43
       CL(N)=CDSQRT(DCMPLX(C(57)/(C(24)*C(54)*RL(N))*                     BNDRY2.44
     1       (T1*2./7.+(T2-T1)/C(3))-.25))-.5*CI                          BNDRY2.45
      ENDDO                                                               BNDRY2.46
      JM=2*JMAX+1                                                         BNDRY2.47
C     ****                                                                BNDRY2.48
C     ****     SET UP HOUGH FUNCTIONS                                     BNDRY2.49
C     ****                                                                BNDRY2.50
      DO L=1,19                                                           BNDRY2.51
       DO M=1,2                                                           BNDRY2.52
        CALL LFK(L,M,CP,W)                                                BNDRY2.53
        CALL LFP(L,M,JM,CP,PT)                                            BNDRY2.54
        DO J=1,JMAX                                                       BNDRY2.55
         P(J,L,M)=PT(2*(JMAX+1-J))                                        BNDRY2.56
        ENDDO                                                             BNDRY2.57
       ENDDO                                                              BNDRY2.58
       DO J=1,JMAX                                                        BNDRY2.59
        P(J,L,2)=DSQRT(DBLE(L*(L+1)-2))*P(J,L,2)-TN(J)*P(J,L,1)           BNDRY2.60
       ENDDO                                                              BNDRY2.61
      ENDDO                                                               BNDRY2.62
      DO L=1,1                                                            BNDRY2.63
       DO LD=1,2                                                          BNDRY2.64
        DO J=1,JMAX                                                       BNDRY2.65
         CALL SDOT(HOUGH(J,L,LD),19,P(J,1,LD),JMAX,B(L,1),1)              BNDRY2.66
        ENDDO                                                             BNDRY2.67
       ENDDO                                                              BNDRY2.68
      ENDDO                                                               BNDRY2.69
C     ****                                                                BNDRY2.70
C     ****      GENERATE ZB2, UB2, VB2 AND TB2                            BNDRY2.71
C     ****                                                                BNDRY2.72
      DO J=1,JMAX                                                         BNDRY2.73
       TB2(J)=0.                                                          BNDRY2.74
       ZB2(J)=0.                                                          BNDRY2.75
       DZB(J)=0.                                                          BNDRY2.76
      ENDDO                                                               BNDRY2.77
      DO L=1,1                                                            BNDRY2.78
       DO J=1,JMAX                                                        BNDRY2.79
        ZB2(J)=ZB2(J)+ZEE(L)*HOUGH(J,L,1)                                 BNDRY2.80
        DZB(J)=DZB(J)+ZEE(L)*HOUGH(J,L,2)                                 BNDRY2.81
        TB2(J)=TB2(J)+CI*C(24)*C(54)/C(57)*ZEE(L)*CL(L)*HOUGH(J,L,1)      BNDRY2.82
       ENDDO                                                              BNDRY2.83
      ENDDO                                                               BNDRY2.84
      DO J=1,JMAX                                                         BNDRY2.85
       UB2(J)=.5*C(23)*C(51)*(1.-(COR(J)/(.5*C(23)))**2)                  BNDRY2.86
       VB2(J)=CI*C(54)*(DZB(J)-COR(J)/(.5*C(23)*CS(J))*ZB2(J))/UB2(J)     BNDRY2.87
       UB2(J)=C(54)*(COR(J)/(.5*C(23))*DZB(J)-1./CS(J)*ZB2(J))/UB2(J)     BNDRY2.88
      ENDDO                                                               BNDRY2.89
      GO TO 11                                                            BNDRY2.90
   13 CONTINUE                                                            BNDRY2.91
C     ****                                                                BNDRY2.92
C     ****      ZERO BOUNDARY CONDITION                                   BNDRY2.93
C     ****                                                                BNDRY2.94
      DO J=1,JMAX                                                         BNDRY2.95
       ZB2(J)=0.                                                          BNDRY2.96
       TB2(J)=0.                                                          BNDRY2.97
       UB2(J)=0.                                                          BNDRY2.98
       VB2(J)=0.                                                          BNDRY2.99
      ENDDO                                                               BNDRY2.100
   11 CONTINUE                                                            BNDRY2.101
C     ****                                                                BNDRY2.102
C     ****      CALCULATE LONGITUDINAL STRUCTURE                          BNDRY2.103
C     ****                                                                BNDRY2.104
      RLAMDA=-2.*C(1)                                                     BNDRY2.105
      BND2(1)=CDEXP(CI*RLAMDA)                                            BNDRY2.106
      EXPDLM=CDEXP(CI*C(1))                                               BNDRY2.107
      DO I=2,IMAXP4                                                       BNDRY2.108
       BND2(I)=BND2(I-1)*EXPDLM                                           BNDRY2.109
      ENDDO                                                               BNDRY2.110
      RETURN                                                              BNDRY2.111
      END                                                                 BNDRY2.112
C                                                                         BNDRY2.113
*DK BNDRYA                                                                BNDRYA.1
      SUBROUTINE BNDRYA                                                   BNDRYA.2
C     ****                                                                BNDRYA.3
C     ****     TIDAL BOUNDARY CONDITION FOR ANNUAL MODE                   BNDRYA.4
C     ****                                                                BNDRYA.5
      INCLUDE "param.h"                                                   BNDRYA.6
      INCLUDE "const.h"                                                   BNDRYA.7
      INCLUDE "start.h"                                                   BNDRYA.8
      DIMENSION B(6,24),RL(0:6)                                           BNDRYA.9
      COMPLEX*16 CC(0:6,0:6),CL(0:6),EXPT,DZB,ZZB                         BNDRYA.10
      COMMON/CBNDRYA/P(ZJMX,24,2),HOUGH(ZJMX,0:6,2),PT(4*ZJMX+1),CP(20),  BNDRYA.11
     1W(70)                                                               BNDRYA.12
      COMMON/CBNDRYA1/DZB(ZJMX),ZZB(ZJMX)                                 BNDRYA.13
      DATA ((B(I,J),I = 1,6),J = 1,12)/                                   BNDRYA.14
     1-0.882922, 0.000000,-0.345087, 0.000000,-0.202228, 0.000000,        BNDRYA.15
     2 0.000000,-0.930826, 0.000000,-0.301357, 0.000000,-0.152720,        BNDRYA.16
     3-0.466226, 0.000000, 0.567457, 0.000000, 0.407114, 0.000000,        BNDRYA.17
     4 0.000000,-0.362673, 0.000000, 0.694431, 0.000000, 0.438014,        BNDRYA.18
     5-0.055436, 0.000000, 0.711847, 0.000000,-0.163050, 0.000000,        BNDRYA.19
     6 0.000000,-0.044983, 0.000000, 0.625545, 0.000000,-0.325772,        BNDRYA.20
     7-0.002909, 0.000000, 0.225723, 0.000000,-0.749160, 0.000000,        BNDRYA.21
     8 0.000000,-0.002773, 0.000000, 0.186467, 0.000000,-0.723674,        BNDRYA.22
     9-0.000086, 0.000000, 0.034940, 0.000000,-0.435919, 0.000000,        BNDRYA.23
     * 0.000000,-0.000103, 0.000000, 0.029425, 0.000000,-0.379254,        BNDRYA.24
     1-0.000002, 0.000000, 0.003267, 0.000000,-0.122687, 0.000000,        BNDRYA.25
     2 0.000000,-0.000003, 0.000000, 0.002928, 0.000000,-0.104008/        BNDRYA.26
      DATA ((B(I,J),I = 1,6),J = 13,24)/                                  BNDRYA.27
     3 0.0     , 0.000000, 0.000206, 0.000000,-0.021267, 0.000000,        BNDRYA.28
     4 0.0     , 0.0     , 0.000000, 0.000202, 0.000000,-0.018228,        BNDRYA.29
     5 0.0     , 0.0     , 0.000009, 0.000000,-0.002540, 0.000000,        BNDRYA.30
     6 0.0     , 0.0     , 0.000000, 0.000010, 0.000000,-0.002252,        BNDRYA.31
     7 0.0     , 0.0     , 0.0     , 0.000000,-0.000223, 0.000000,        BNDRYA.32
     8 0.0     , 0.0     , 0.0     , 0.0     , 0.000000,-0.000208,        BNDRYA.33
     9 0.0     , 0.0     , 0.0     , 0.0     ,-0.000015, 0.000000,        BNDRYA.34
     * 0.0     , 0.0     , 0.0     , 0.0     , 0.000000,-0.000015,        BNDRYA.35
     1 0.0     , 0.0     , 0.0     , 0.0     ,-0.000001, 0.000000,        BNDRYA.36
     2 0.0     , 0.0     , 0.0     , 0.0     , 0.000000,-0.000001,        BNDRYA.37
     3 0.0     , 0.0     , 0.0     , 0.0     , 0.0     , 0.000000,        BNDRYA.38
     4 0.0     , 0.0     , 0.0     , 0.0     , 0.0     , 0.0     /        BNDRYA.39
      DATA (RL(N),N=1,6)/                                                 BNDRYA.40
     1  -10.8409D5,-7.0243D5,-2.4874D5,-1.9696D5,-1.0694D5,               BNDRYA.41
     2  -0.9119D5/                                                        BNDRYA.42
      DATA T2/184./,T1/181./                                              BNDRYA.43
C     ****                                                                BNDRYA.44
C     ****     ANNUAL BOUNDARY COEFFICIENTS FOR LOWER BOUNDARY AT         BNDRYA.45
C     ****       Z = -7.                                                  BNDRYA.46
C     ****                                                                BNDRYA.47
      DATA((CC(K,N),N=0,6),K=0,3)/                                        BNDRYA.48
     1  ( 0.136291D+03, 0.000000D+00),(-0.200536D-01, 0.000000D+00),      BNDRYA.49
     2  ( 0.423456D+00, 0.000000D+00),(-0.143623D-02, 0.000000D+00),      BNDRYA.50
     3  ( 0.262889D+00, 0.000000D+00),( 0.365122D-02, 0.000000D+00),      BNDRYA.51
     4  ( 0.102716D+00, 0.000000D+00),( 0.826209D+00,-0.116536D-01),      BNDRYA.52
     5  (-0.455993D+00,-0.311380D-02),( 0.271258D-01,-0.373420D-01),      BNDRYA.53
     6  (-0.433095D-01, 0.165546D-01),( 0.713386D-02,-0.112538D-01),      BNDRYA.54
     7  (-0.153849D-01, 0.103063D-01),( 0.214366D-02,-0.570878D-02),      BNDRYA.55
     8  (-0.316298D+00, 0.115053D+00),(-0.159072D-01,-0.245495D-02),      BNDRYA.56
     9  ( 0.302211D+00,-0.132446D-01),( 0.230750D-02, 0.170566D-03),      BNDRYA.57
     *  ( 0.100434D+00,-0.299227D-02),( 0.264555D-02,-0.137723D-03),      BNDRYA.58
     1  ( 0.499098D-01,-0.110255D-02),(-0.520584D-03,-0.114124D-02),      BNDRYA.59
     2  ( 0.178599D-01, 0.561092D-02),( 0.557591D-03, 0.176165D-02),      BNDRYA.60
     3  (-0.148151D-02, 0.749397D-03),(-0.617325D-03, 0.778743D-03),      BNDRYA.61
     4  (-0.530835D-03, 0.641767D-03),(-0.964206D-03, 0.551394D-03)/      BNDRYA.62
      DATA((CC(K,N),N=0,6),K=4,6)/                                        BNDRYA.63
     1  (-0.137927D-02, 0.866386D-03),(-0.242825D-02, 0.441184D-03),      BNDRYA.64
     2  ( 0.120715D-02,-0.136729D-02),( 0.122657D-03, 0.316213D-04),      BNDRYA.65
     3  ( 0.390769D-03,-0.162978D-03),( 0.378377D-03,-0.195668D-04),      BNDRYA.66
     4  ( 0.366912D-03,-0.681579D-04),(-0.470068D-03,-0.118650D-05),      BNDRYA.67
     5  ( 0.120025D-02,-0.797459D-03),( 0.622700D-03,-0.424648D-04),      BNDRYA.68
     6  (-0.537275D-03, 0.101658D-03),( 0.222407D-03,-0.828812D-05),      BNDRYA.69
     7  (-0.209097D-03, 0.828365D-04),( 0.945940D-04, 0.317248D-04),      BNDRYA.70
     8  ( 0.341903D-03, 0.192246D-04),( 0.129833D-03, 0.247156D-04),      BNDRYA.71
     9  (-0.610206D-03, 0.591081D-06),(-0.102160D-03,-0.434110D-04),      BNDRYA.72
     *  ( 0.196672D-04,-0.305687D-04),(-0.905354D-04,-0.813929D-04),      BNDRYA.73
     1  ( 0.569460D-05,-0.116661D-03)/                                    BNDRYA.74
      DATA SCALE/1./                                                      BNDRYA.75
      RT2=DSQRT(2.D0)                                                     BNDRYA.76
      ITIDEA=0                                                            BNDRYA.77
      CI=(0.,1.)                                                          BNDRYA.78
      IF(ITIDEA.EQ.0)GO TO 13                                             BNDRYA.79
      JM=2*JMAX+1                                                         BNDRYA.80
C     ****                                                                BNDRYA.81
C     ****     HEIGHT VARIATION                                           BNDRYA.82
C     ****                                                                BNDRYA.83
      DO N=1,6                                                            BNDRYA.84
       CL(N)=-CDSQRT(DCMPLX(C(57)/(C(24)*C(54)*RL(N))*                    BNDRYA.85
     1       (T1*2./7.+(T2-T1)/C(3))-.25))-.5*CI                          BNDRYA.86
      ENDDO                                                               BNDRYA.87
      CL(0)=0.                                                            BNDRYA.88
C     ****                                                                BNDRYA.89
C     ****     SET UP LEGENDRE POLYNOMIALS                                BNDRYA.90
C     ****                                                                BNDRYA.91
      DO L=1,24                                                           BNDRYA.92
       DO M=0,1                                                           BNDRYA.93
        MP1=M+1                                                           BNDRYA.94
        CALL LFK(L,M,CP,W)                                                BNDRYA.95
        CALL LFP(L,M,JM,CP,PT)                                            BNDRYA.96
        DO J=1,JMAX                                                       BNDRYA.97
         P(J,L,MP1)=PT(2*(JMAX+1-J))                                      BNDRYA.98
        ENDDO                                                             BNDRYA.99
       ENDDO                                                              BNDRYA.100
       M=0                                                                BNDRYA.101
       RM=DBLE(M)                                                         BNDRYA.102
       DO J=1,JMAX                                                        BNDRYA.103
        P(J,L,2)=DSQRT(DBLE(L*(L+1)-M*(M+1)))*P(J,L,2)-RM*TN(J)*          BNDRYA.104
     1           P(J,L,1)                                                 BNDRYA.105
       ENDDO                                                              BNDRYA.106
      ENDDO                                                               BNDRYA.107
C     ****                                                                BNDRYA.108
C     ****     NOW EVALUATE HOUGH FUNCTIONS                               BNDRYA.109
C     ****                                                                BNDRYA.110
      DO L=1,6                                                            BNDRYA.111
       DO LD=1,2                                                          BNDRYA.112
        DO J=1,JMAX                                                       BNDRYA.113
         CALL SDOT(HOUGH(J,L,LD),24,P(J,1,LD),JMAX,B(L,1),6)              BNDRYA.114
        ENDDO                                                             BNDRYA.115
       ENDDO                                                              BNDRYA.116
      ENDDO                                                               BNDRYA.117
C     ****                                                                BNDRYA.118
C     ****     HOUGH FUNCTION OF ORDE ZERO`                               BNDRYA.119
C     ****                                                                BNDRYA.120
      DO J=1,JMAX                                                         BNDRYA.121
       HOUGH(J,0,1)=1./RT2                                                BNDRYA.122
       HOUGH(J,0,2)=0.                                                    BNDRYA.123
      ENDDO                                                               BNDRYA.124
C     ****                                                                BNDRYA.125
C     ****     GENERATE ZBA, TBA, UBA, VBA                                BNDRYA.126
C     ****                                                                BNDRYA.127
      DO J=1,JMAX                                                         BNDRYA.128
       ZBA(J)=0.                                                          BNDRYA.129
       TBA(J)=0.                                                          BNDRYA.130
       UBA(J)=0.                                                          BNDRYA.131
       VBA(J)=0.                                                          BNDRYA.132
      ENDDO                                                               BNDRYA.133
C     ****                                                                BNDRYA.134
C     ****     SUMMATION OVER FREQUENCY, K, CALCULATION OF PHASE FACTOR   BNDRYA.135
C     ****                                                                BNDRYA.136
      DO K=0,6                                                            BNDRYA.137
       EXPT=CDEXP(CI*DBLE(K*(IIDAY-1))*86400.*C(25))                      BNDRYA.138
C      ****                                                               BNDRYA.139
C      ****     SUMMATION OVER ORDER, N                                   BNDRYA.140
C      ****                                                               BNDRYA.141
       DO N=0,6                                                           BNDRYA.142
        FACTOR=SCALE                                                      BNDRYA.143
        IF(K.EQ.0.AND.N.EQ.0)FACTOR=1.                                    BNDRYA.144
        DO J=1,JMAX                                                       BNDRYA.145
         ZZB(J)=CC(K,N)*HOUGH(J,N,1)*EXPT*1.D5*FACTOR                     BNDRYA.146
         DZB(J)=CC(K,N)*HOUGH(J,N,2)*EXPT*1.D5*FACTOR                     BNDRYA.147
         W(J)=C(54)/(C(51)*((DBLE(K)*C(25))**2-COR(J)**2))                BNDRYA.148
         ZBA(J)=ZBA(J)+ZZB(J)                                             BNDRYA.149
         TBA(J)=TBA(J)+CI*C(24)*C(54)/C(57)*CL(N)*ZZB(J)                  BNDRYA.150
         UBA(J)=UBA(J)+W(J)*(COR(J)*DZB(J)-RM*DBLE(K)*C(25)/              BNDRYA.151
     1          CS(J)*ZZB(J))                                             BNDRYA.152
         VBA(J)=VBA(J)+CI*W(J)*(DBLE(K)*C(25)*DZB(J)-RM*COR(J)/           BNDRYA.153
     1          CS(J)*ZZB(J))                                             BNDRYA.154
         ENDDO                                                            BNDRYA.155
        ENDDO                                                             BNDRYA.156
       ENDDO                                                              BNDRYA.157
      GO TO 11                                                            BNDRYA.158
   13 CONTINUE                                                            BNDRYA.159
C     ****                                                                BNDRYA.160
C     ****     ZERO BOUNDARY CONDITION                                    BNDRYA.161
C     ****                                                                BNDRYA.162
      DO J=1,JMAX                                                         BNDRYA.163
       ZBA(J)=CC(0,0)/RT2*1.D5                                            BNDRYA.164
       TBA(J)=0.                                                          BNDRYA.165
       UBA(J)=0.                                                          BNDRYA.166
       VBA(J)=0.                                                          BNDRYA.167
      ENDDO                                                               BNDRYA.168
   11 CONTINUE                                                            BNDRYA.169
C     ****                                                                BNDRYA.170
C     ****     LONGITUDINAL STRUCTURE                                     BNDRYA.171
C     ****                                                                BNDRYA.172
      DO I=1,IMAXP4                                                       BNDRYA.173
       BNDA(I)=1.                                                         BNDRYA.174
      ENDDO                                                               BNDRYA.175
      RETURN                                                              BNDRYA.176
      END                                                                 BNDRYA.177
C                                                                         BNDRYA.178
*DK CNSTNT                                                                CNSTNT.1
      SUBROUTINE CNSTNT                                                   CNSTNT.2
C     ****                                                                CNSTNT.3
C     ****     DEFINES AND CALCULATES CONSTANTS NEEDED BY PROGRAM         CNSTNT.4
C     ****                                                                CNSTNT.5
      INCLUDE 'param.h'                                                   CNSTNT.6
      INCLUDE 'epara.h'                                                   CNSTNT.7
      INCLUDE 'econs.h'                                                   CNSTNT.8
      DATA ALFA/1.668/                                                    CNSTNT.9
      DIMENSION TABL2(181,3:5),THETA0(JMAXM)                               CNSTNT.11
C     ****                                                                CNSTNT.13
C     ****    NOW SET CONSTANTS                                           CNSTNT.14
C     ****                                                                CNSTNT.15
      RE=6.378165D8                                                       CNSTNT.16
      H0=9.0D6                                                            CNSTNT.17
      H00=9.7D6                                                           CNSTNT.18
      R1=1.06D7                                                           CNSTNT.24
      R0=RE+H0                                                            CNSTNT.25
      R00=RE+H00                                                          CNSTNT.26
      PI=4*DATAN(1.D0)                                                    CNSTNT.28
      DTR=PI/180.                                                         CNSTNT.29
      RTD=180./PI                                                         CNSTNT.30
      DLATM=PI/DBLE(JMAXM-1)                                              CNSTNT.31
      DLONM=2.*PI/DBLE(IMAXM)                                             CNSTNT.32
C     ****                                                                CNSTNT.33
C     ****     FILL ARRAY YLATM (EQUALLY SPACED IN THETAO BUT HOLDS       CNSTNT.34
C     ****       CORRESPONDING VALUE OF THETAS)                           CNSTNT.35
C     ****                                                                CNSTNT.36
      DO J=1,JMAXM                                                        CNSTNT.37
C      ****                                                               CNSTNT.38
C      ****    THETA0 = EQUALLY SPACED GRID VALUES                        CNSTNT.39
C      ****                                                               CNSTNT.40
       THETA0(J)=-PI/2.+DBLE(J-1)*DLATM                                   CNSTNT.41
      ENDDO                                                               CNSTNT.42
      DO J=2,JMAXM-1                                                      CNSTNT.43
C      ****                                                               CNSTNT.44
C      ****     TANTH0 = ABS(TAN(THETA0))                                 CNSTNT.46
C      ****                                                               CNSTNT.47
       TANTH0=DABS(DTAN(THETA0(J)))                                       CNSTNT.48
C      ****                                                               CNSTNT.49
C      ****     HAMH0 = HA-H0                                             CNSTNT.50
C      ****                                                               CNSTNT.51
       HAMH0=R1*TANTH0+R0*TANTH0**(2+2*ALFA)/                             CNSTNT.52
     1          (1+TANTH0**2)**ALFA                                       CNSTNT.53
C      ****                                                               CNSTNT.54
C      ****     TANTHS = ABS(TAN(THETAS))                                 CNSTNT.55
C      ****                                                               CNSTNT.56
       TANTHS=DSQRT(HAMH0/R0)                                             CNSTNT.57
C      ****                                                               CNSTNT.58
C      ****     YLATM = TANTHS                                            CNSTNT.59
C      ****                                                               CNSTNT.60
       YLATM(J)=DSIGN(DATAN(TANTHS),THETA0(J))                            CNSTNT.61
C      ****                                                               CNSTNT.62
C      ****     RCOS0S = COS(THETA0)/COS(THETAS)                          CNSTNT.63
C      ****                                                               CNSTNT.64
       RCOS0S(J)=DSQRT((1+TANTHS**2)/(1+TANTH0**2))                       CNSTNT.65
C      ****                                                               CNSTNT.66
C      ****     DT0DTS = D(THETA0)/D(THETAS)                              CNSTNT.67
C      ****                                                               CNSTNT.68
C      ****     DT1DTS = DT0DTS DIVIDED BY ABS[SIN(Im)]. REMIANS FINITE   CNSTNT.69
C      ****              AND NONZERO AY MAGNETIC EQUATOR.                 CNSTNT.70
C      ****                                                               CNSTNT.71
       TANTHS2=TANTHS**2                                                  CNSTNT.72
       DT1DTS(J)=(R0*DSQRT(1+4*TANTHS2)*(1+TANTHS2))/                     CNSTNT.73
     1           (R1*(1+TANTH0**2)+2*R0*TANTH0**(2*ALFA+1)*               CNSTNT.74
     2           (1+ALFA+TANTH0**2)/(1+TANTH0**2)**ALFA)                  CNSTNT.75
       DT0DTS(J)=DT1DTS(J)*2*TANTHS/DSQRT(1+4*TANTHS2)                    CNSTNT.76
      ENDDO                                                               CNSTNT.77
C     ****                                                                CNSTNT.78
C     ****     NOW DO POLES                                               CNSTNT.79
C     ****                                                                CNSTNT.80
      YLATM(1)=THETA0(1)                                                  CNSTNT.81
      YLATM(JMAXM)=THETA0(JMAXM)                                          CNSTNT.82
      RCOS0S(1)=1.                                                        CNSTNT.83
      RCOS0S(JMAXM)=1.                                                    CNSTNT.84
      DT0DTS(1)=1.                                                        CNSTNT.85
      DT0DTS(JMAXM)=1.                                                    CNSTNT.86
      DT1DTS(1)=1.                                                        CNSTNT.87
      DT1DTS(JMAXM)=1.                                                    CNSTNT.88
C     ****                                                                CNSTNT.89
C     ****     FILL YLONM                                                 CNSTNT.90
C     ****                                                                CNSTNT.91
      DO I=1,IMAXMP                                                       CNSTNT.92
       YLONM(I)=-PI+DBLE(I-1)*DLONM                                       CNSTNT.93
      ENDDO                                                               CNSTNT.94
      DTHETA=PI/(2.*180.)                                                 CNSTNT.95
      TABLE(1,1)=0.                                                       CNSTNT.96
      TABLE(1,2)=0.                                                       CNSTNT.97
      DO I=2,181                                                          CNSTNT.98
       TABLE(I,1)=TABLE(I-1,1)+DTHETA                                     CNSTNT.99
      ENDDO                                                               CNSTNT.100
      DO I=2,180                                                          CNSTNT.101
       TABL2(I,4)=DTAN(TABLE(I,1))                                        CNSTNT.102
       TABLE(I,2)=TABLE(I,1)                                              CNSTNT.103
      ENDDO                                                               CNSTNT.104
      DO N=1,7                                                            CNSTNT.105
       DO I=2,180                                                         CNSTNT.106
        TABL2(I,3)=TABLE(I,2)                                             CNSTNT.107
        TABLE(I,2)=DTAN(TABL2(I,3))                                       CNSTNT.108
        TABL2(I,5)=DSQRT(R1/R0*TABLE(I,2)+TABLE(I,2)**                    CNSTNT.109
     1             (2*(1+ALFA))/(1+TABLE(I,2)**2)**ALFA)                  CNSTNT.110
        TABLE(I,2)=TABL2(I,3)-(TABL2(I,5)-TABL2(I,4))*2*TABL2(I,5)/       CNSTNT.111
     2             (R1/R0*(1+TABLE(I,2)**2)+2*TABLE(I,2)**(2*ALFA+1)*     CNSTNT.112
     3             (1+ALFA+TABLE(I,2)**2)/(1+TABLE(I,2)**2)**ALFA)        CNSTNT.113
       ENDDO                                                              CNSTNT.114
      ENDDO                                                               CNSTNT.115
      RETURN                                                              CNSTNT.116
      END                                                                 CNSTNT.117
C                                                                         CNSTNT.118
*DK CON                                                                      CON.1
      SUBROUTINE CON                                                         CON.2
C     ****                                                                   CON.3
C     ****          DEFINES THE MODEL CONSTANTS                              CON.4
C     ****                                                                   CON.5
      INCLUDE 'param.h'                                                      CON.6
      INCLUDE 'buffe.h'                                                      CON.7
      INCLUDE 'crfft.h'                                                      CON.8
      INCLUDE 'const.h'                                                      CON.9
      INCLUDE 'index.h'                                                      CON.10
      INCLUDE 'start.h'                                                      CON.11
      INCLUDE 'field.h'                                                      CON.12
      DATA NPHYS/39/                                                         CON.13
C     ****                                                                   CON.14
C     ****          RESOLUTION PARAMETERS                                    CON.15
C     ****                                                                   CON.16
      IMAX=ZIMX                                                              CON.17
      JMAX=ZJMX                                                              CON.18
      KMAX=ZKMX                                                              CON.19
C     ****                                                                   CON.20
C     ****          DERIVED CONSTANTS                                        CON.21
C     ****                                                                   CON.22
      PI=3.14159265358979                                                    CON.23
      DLAMDA=2*PI/IMAX                                                       CON.24
      DPHI=PI/JMAX                                                           CON.25
      DS=(ZST-ZSB)/KMAX                                                      CON.26
      IMAXH=IMAX/2                                                           CON.27
      IMAXP2=IMAX+2                                                          CON.28
      IMAXP4=IMAX+4                                                          CON.29
      JMAXM1=JMAX-1                                                          CON.31
      JMAXM2=JMAX-2                                                          CON.32
      KMAXM1=KMAX-1                                                          CON.33
      KMAXP1=KMAX+1                                                          CON.34
C     ****                                                                   CON.35
C     ****          BUFFER INFORMATION AND INITIALIZE POINTERS               CON.36
C     ****                                                                   CON.37
      LEN1=IMAXP4                                                            CON.38
      LEN2=IMAXP4*KMAX                                                       CON.39
      LEN3=IMAXP4*KMAXP1                                                     CON.40
      NCBUF=0                                                                CON.41
      DO N=1,NFLDS                                                           CON.42
       IF(N.EQ.17)THEN                                                       CON.43
        K=1                                                                  CON.44
       ELSE                                                                  CON.45
        K=KMAXP1                                                             CON.46
       ENDIF                                                                 CON.47
       NDEXA(N+1)=NCBUF                                                      CON.48
       NCBUF=NCBUF+K                                                         CON.49
      ENDDO                                                                  CON.50
      NCPHYS=0                                                               CON.51
      DO N=1,NPHYS                                                           CON.52
       K=KMAXP1                                                              CON.53
       NDEXB(N+1)=NCPHYS                                                     CON.54
       NCPHYS=NCPHYS+K                                                       CON.55
      ENDDO                                                                  CON.56
C     ****                                                                   CON.57
C     ****          COMPUTE COLUMN POINTERS FOR SCM BUFFERS                  CON.58
      NBADDR(1)=1                                                            CON.59
      DO I=2,8                                                               CON.60
       NBADDR(I)=NBADDR(I-1)+NCBUF                                           CON.61
      ENDDO                                                                  CON.62
      NBAPHY=NBADDR(8)+NCBUF                                                 CON.63
      DO N=1,NPHYS                                                           CON.64
       NDEXB(N+1)=NDEXB(N+1)+NBAPHY                                          CON.65
      ENDDO                                                                  CON.66
C     ****                                                                   CON.67
C     ****          OUTPUT POINTER INFORMATION                               CON.68
C     ****                                                                   CON.69
      LIM=NFLDS+NPHYS+2                                                      CON.70
      WRITE(6,1005) (NDEXA(N),N=1,LIM)                                       CON.71
1005  FORMAT(//'  POINTER VALUES FOR EACH FIELD IN COMMON BLOCK NDEX',       CON.72
     A       //(10X,5I5))                                                    CON.73
C     ****                                                                   CON.74
C     ****          INITIALIZE FFT ROUTINE                                   CON.75
C     ****                                                                   CON.76
      CALL SET99(TRIGS,IFAX,IMAX)                                            CON.77
C     ****                                                                   CON.78
C     ****          COMPUTE PHI VARIATION ARRAYS                             CON.79
C     ****                                                                   CON.80
      DELT=86400/DBLE(IMAX)                                                  CON.81
      JMAXH=JMAX/2                                                           CON.82
      OMEGA=7.292D-5                                                         CON.83
      DO J=1,JMAX                                                            CON.84
       PHI=(J-JMAXH-0.5)*DPHI                                                CON.85
       CS(J)=DCOS(PHI)                                                       CON.86
       SN(J)=DSIN(PHI)                                                       CON.87
       TN(J)=DTAN(PHI)                                                       CON.88
       COR(J)=2.*OMEGA*SN(J)                                                 CON.89
      ENDDO                                                                  CON.90
      CSSP(1)=-CS(2)                                                         CON.91
      CSSP(2)=-CS(1)                                                         CON.92
      CSNP(1)=-CS(JMAX)                                                      CON.93
      CSNP(2)=-CS(JMAXM1)                                                    CON.94
      KUTSP(1)=KUT(2)                                                        CON.95
      KUTSP(2)=KUT(1)                                                        CON.96
      KUTNP(1)=KUT(JMAX)                                                     CON.97
      KUTNP(2)=KUT(JMAXM1)                                                   CON.98
C     ****                                                                   CON.99
C     ****          CALCULATE SUN'S DECLINATION                              CON.100
C     ****          (WILL HAVE TO RECALCULATE THIS IF DAY IS ADVANCED)       CON.101
C     ****                                                                   CON.102
      DELTA=DATAN(DTAN(23.5*PI/180.)*DSIN(2.*PI*DBLE(IIDAY-80)/365.))        CON.103
C     ****                                                                   CON.104
C     ****          GRID CONSTANTS                                           CON.105
C     ****                                                                   CON.106
      C(1)=DLAMDA                                                            CON.107
      C(2)=DPHI                                                              CON.108
      C(3)=DS                                                                CON.109
      C(6)=2.*C(4)                                                           CON.110
      C(7)=1./C(6)                                                           CON.111
C     ****                                                                   CON.112
C     ****          FINITE DIFFERENCE OPERATOR CONSTANTS                     CON.113
C     ****                                                                   CON.114
      C(10)=2./(3.*DLAMDA)                                                   CON.115
      C(11)=1./(12.*DLAMDA)                                                  CON.116
      C(12)=2./(3.*DPHI)                                                     CON.117
      C(13)=1./(12.*DPHI)                                                    CON.118
      C(14)=1./(2.*DLAMDA)                                                   CON.119
      C(15)=1./(2.*DPHI)                                                     CON.120
C     ****                                                                   CON.121
C     ****          INTEGRATION CONSTANTS                                    CON.122
C     ****                                                                   CON.123
      C(18)=C(3)/(1.+.5*C(3))                                                CON.124
      C(19)=(1.-.5*C(3))/(.1+.5*C(3))                                        CON.125
      C(20)=1./DS                                                            CON.126
      C(22)=DCOS(PI/3.)                                                      CON.127
      C(23)=4.*PI/(24.*60.*60.)                                              CON.128
      C(24)=28.9                                                             CON.129
C     ****                                                                   CON.130
C     ****     ANNUAL FREQUENCY                                              CON.131
C     ****                                                                   CON.132
      C(25)=C(23)/(2.*365.25)                                                CON.133
C     ****                                                                   CON.134
C     ****          SHAPIRO SMOOTHER CONSTANT                                CON.135
C     ****                                                                   CON.136
      C(26)=3.0D-2                                                           CON.137
C     ****                                                                   CON.138
C     ****          EDDY DIFFUSION                                           CON.139
C     ****                                                                   CON.140
      C(27)=5.D-6                                                            CON.141
C     ****                                                                   CON.142
C     ****          EDDY THERMAL CONDUCTION                                  CON.143
C     ****                                                                   CON.144
      C(28)=C(27)                                                            CON.145
C     ****                                                                   CON.146
C     ****          TIME SMOOTHING CONSTANT                                  CON.147
C     ****                                                                   CON.148
      ALPHA=.95                                                              CON.149
      C(30)=ALPHA                                                            CON.150
      C(31)=.5*(1.-ALPHA)                                                    CON.151
C     ****                                                                   CON.152
C     *****         SURPLUS HEAT PER EVENT CONVERTED FROM E.V. TO ERGS       CON.153
C     ****                                                                   CON.154
      C(45)=5.11*1.602D-12                                                   CON.155
C     ****                                                                   CON.156
C     ****          PHYSICAL CONSTANTS                                       CON.157
C     ****                                                                   CON.158
C     ****          EARTHS RADIUS                                            CON.159
C     ****                                                                   CON.160
      C(51)=6.37122D8                                                        CON.161
      C(52)=1./C(51)                                                         CON.162
      C(53)=C(51)*C(51)                                                      CON.163
C     ****                                                                   CON.164
C     ****          GRAVITATIONAL CONSTANT                                   CON.165
C     ****                                                                   CON.166
      C(54)=870.                                                             CON.167
      C(55)=1./C(54)                                                         CON.168
C     ****                                                                   CON.169
C     ****          EARTHS ANGULAR VELOCITY                                  CON.170
C     ****                                                                   CON.171
      C(56)=OMEGA                                                            CON.172
C     ****                                                                   CON.173
C     ****          GAS CONSTANT                                             CON.174
C     ****                                                                   CON.175
      C(57)=8.314D7                                                          CON.176
C     ****                                                                   CON.177
C     ****          H*C                                                      CON.178
C     ****                                                                   CON.179
      C(60)=1.9845D-16                                                       CON.180
C     ****                                                                   CON.181
C     ****          G/R                                                      CON.182
C     ****                                                                   CON.183
      C(65)=C(54)/C(57)                                                      CON.184
C     ****                                                                   CON.185
C     ****          STANDARD PRESSURE                                        CON.186
C     ****                                                                   CON.187
      C(81)=5.0D-4                                                           CON.188
C     ****                                                                   CON.189
C     ****          BOLTZMANN CONST                                          CON.190
C     ****                                                                   CON.191
      C(84)=1.38D-16                                                         CON.192
C     ****                                                                   CON.193
C     ****          AVOGRADO NUMBER                                          CON.194
C     ****                                                                   CON.195
      C(85)=6.023D23                                                         CON.196
C     ****                                                                   CON.197
C     ****          EXP(-.5*DS)                                              CON.198
C     ****                                                                   CON.199
      C(86)=DEXP(-.5*DS)                                                     CON.200
C     ****                                                                   CON.201
C     ****          EXP(.5*DS)                                               CON.202
C     ****                                                                   CON.203
      C(87)=1./C(86)                                                         CON.204
C     ****                                                                   CON.205
C     ****          SUN'S DECLINATION:                                       CON.206
C     ****          (WILL HAVE TO RECALCULATE IF DAY IS ADVANCED)            CON.207
C     ****                                                                   CON.208
      C(95)=DSIN(DELTA)                                                      CON.209
      C(96)=DCOS(DELTA)                                                      CON.210
C     ****                                                                   CON.211
C     ****          KO  COEFFICIENT OF EDDY VISCOSITY                        CON.212
C     ****                                                                   CON.213
      C(109)=0.2                                                             CON.214
      C(110)=PI                                                              CON.215
C     ****                                                                   CON.216
C     ****        SAVE ORIGINAL DT VALUE                                     CON.217
C     ****                                                                   CON.218
      C(112)=DSQRT(PI)                                                       CON.219
C     ****                                                                   CON.220
C     ****          CALCULATE EXPS AND A ARRAYS                              CON.221
C     ****                                                                   CON.222
      S=ZSB+.5*DS                                                            CON.223
      EXPS(1)=DEXP(-S)                                                       CON.224
      EXPDS=DEXP(-DS)                                                        CON.225
      A(1)=1.D-6*DEXP(-7.-S)                                                 CON.226
      DIFK(1)=C(27)                                                          CON.227
      DIFT(1)=C(28)                                                          CON.228
      XMUE(1)=C(27)                                                          CON.229
      DO 1 K=2,KMAX                                                          CON.230
      EXPS(K)=EXPS(K-1)*EXPDS                                                CON.231
      DIFK(K)=DIFK(K-1)*EXPDS                                                CON.232
      DIFT(K)=DIFT(K-1)*EXPDS                                                CON.233
      XMUE(K)=DIFK(K)                                                        CON.234
    1 CONTINUE                                                               CON.235
      DIFK(KMAXP1)=DIFK(KMAX)*EXPDS                                          CON.236
      DIFT(KMAXP1)=DIFT(KMAX)*EXPDS                                          CON.237
      XMUE(KMAXP1)=DIFK(KMAXP1)                                              CON.238
      RETURN                                                                 CON.239
      END                                                                    CON.240
C                                                                            CON.241
*DK EINPUT                                                                EINPUT.1
      SUBROUTINE EINPUT                                                   EINPUT.2
      INCLUDE 'param.h'                                                   EINPUT.3
      INCLUDE 'const.h'                                                   EINPUT.4
      INCLUDE 'epara.h'                                                   EINPUT.5
      INCLUDE 'econs.h'                                                   EINPUT.6
      INCLUDE 'efild.h'                                                   EINPUT.8
      INCLUDE 'eterp.h'                                                   EINPUT.9
      INCLUDE 'edphi.h'                                                   EINPUT.10
C     ****                                                                EINPUT.11
C     ****     READ IN GEOMAGNETIC COORDINATE TRANSFORM METRIX            EINPUT.12
C     ****                                                                EINPUT.13
      IMAG=23                                                             EINPUT.14
      OPEN(IMAG,FILE='magnet.dat',FORM='UNFORMATTED')                     EINPUT.15
      DO I=1,IMAXGP*(JMAXGP+1)                                            EINPUT.16
       READ(IMAG)ALATM(I,0),ALONM(I,0),XB(I,0),YB(I,0),ZZB(I,0),          EINPUT.17
     1           BMOD(I,0)                                                EINPUT.18
      ENDDO                                                               EINPUT.19
      DO I=1,IMAXGP*(JMAXGP+1)*4                                          EINPUT.20
       READ(IMAG)RJAC(I,0,1,1)                                            EINPUT.21
      ENDDO                                                               EINPUT.22
      DO I=1,IMAXGP*(JMAXGP+1)*6                                          EINPUT.23
       READ(IMAG)AV(I,0,1,1)                                              EINPUT.24
      ENDDO                                                               EINPUT.25
      DO J=1,JMAXM                                                        EINPUT.29
       DO I=1,IMAXMP                                                      EINPUT.30
        READ(IMAG)IG(I,J),JG(I,J)                                         EINPUT.31
       ENDDO                                                              EINPUT.32
      ENDDO                                                               EINPUT.33
      DO J=1,JMAXM                                                        EINPUT.34
       DO I=1,IMAXMP                                                      EINPUT.35
        READ(IMAG)WT(1,I,J),WT(2,I,J),WT(3,I,J),WT(4,I,J)                 EINPUT.36
       ENDDO                                                              EINPUT.37
      ENDDO                                                               EINPUT.38
      DO I=1,IMAXGP*JMAXG*6                                               EINPUT.51
       READ(IMAG)DVEC(I,1,1,1)                                            EINPUT.52
      ENDDO                                                               EINPUT.53
      DO I=1,IMAXGP*JMAXG                                                 EINPUT.54
       READ(IMAG)DDDARR(I,1)                                              EINPUT.55
      ENDDO                                                               EINPUT.56
      DO I=1,IMAXGP*JMAXG                                                 EINPUT.57
       READ(IMAG)BE3ARR(I,1)                                              EINPUT.58
      ENDDO                                                               EINPUT.59
      CLOSE(IMAG)                                                         EINPUT.60
C     ****                                                                EINPUT.84
C     ****     CALCULATE MAGNETIC FIELDS                                  EINPUT.85
C     ****                                                                EINPUT.86
      CALL MAGDYN                                                         EINPUT.87
      RETURN                                                              EINPUT.106
      END                                                                 EINPUT.107
C                                                                         EINPUT.108
*DK FGET                                                                    FGET.1
      FUNCTION FGET(STR,IER)                                                FGET.2
      IMPLICIT REAL*8 (A-H,O-Z)                                             FGET.3
      PARAMETER(TOOBIG=.3402499962D+39,TOOSMALL=-.3402499962D+39)           FGET.4
      CHARACTER STR*(*)                                                     FGET.5
C     ****                                                                  FGET.6
C     ****     GET A REAL FROM STDIN:                                       FGET.7
C     ****     IF IER=0 ON RETURN, ALL'S WELL                               FGET.8
C     ****     IF IER=1 ON RETURN, FLOAT NUMBER WAS TOO BIG OR TOO SMALL    FGET.9
C     ****     IF IER=2 ON RETURN, A DAD STRING WAS READ (E.G. ALPHAS       FGET.10
C     ****                         INSTEAD OF DIGITS)                       FGET.11
C     ****                                                                  FGET.12
      IER=0                                                                 FGET.13
      READ(STR,*,ERR=900)R                                                  FGET.14
      IF(R.GT.TOOBIG.OR.R.LT.TOOSMALL)THEN                                  FGET.15
       IER=1                                                                FGET.16
       RETURN                                                               FGET.17
      ENDIF                                                                 FGET.18
      FGET=R                                                                FGET.19
      RETURN                                                                FGET.20
C     ****                                                                  FGET.21
C     ****     BAD STRING                                                   FGET.22
C     ****                                                                  FGET.23
900   IRE=2                                                                 FGET.24
      RETURN                                                                FGET.25
      END                                                                   FGET.26
C                                                                           FGET.27
*DK FOGCM                                                                  FOGCM.1
      PROGRAM FOGCM                                                        FOGCM.2
C     ****                                                                 FOGCM.3
C     ****          WORKSTATION GENERATION OF THE TIGCM (VERSION 13)       FOGCM.4
C     ****          JULY, 1999          WENBIN WANG                        FOGCM.5
C     ****                                                                 FOGCM.6
      INCLUDE 'param.h'                                                    FOGCM.7
      INCLUDE 'start.h'                                                    FOGCM.8
      INCLUDE 'nestcnt.h'                                                  FOGCM.9
C     ****                                                                 FOGCM.10
C     ****          INITIALIZE THE MODEL, ASSIGN FORTRAN UNITS             FOGCM.11
C     ****                                                                 FOGCM.12
      IHIST=3                                                              FOGCM.13
      IOTPT=7                                                              FOGCM.14
      IUNIN=8                                                              FOGCM.15
      IUNOT=9                                                              FOGCM.16
      INEST(1)=96                                                          FOGCM.17
      INEST(2)=97                                                          FOGCM.18
      INEST(3)=98                                                          FOGCM.19
      INEST(4)=99                                                          FOGCM.20
      CALL START                                                           FOGCM.21
C     ****                                                                 FOGCM.22
C     ****          OPEN OUTPUT FILE OF THE MODEL. INEST(I), I=1,2,3,4     FOGCM.23
C     ****          ARE UNITS FOR FIRST, SECOND, THIRD AND FOURTH LEVEL    FOGCM.24
C     ****          LEVEL NESTED GRID OUTPUTS, RESPECTIVELY.               FOGCM.25
C     ****                                                                 FOGCM.26
      OPEN(UNIT=IOTPT,FILE=LABHIS(1)(1:LNBLNK(LABHIS(1))),                 FOGCM.27
     +     FORM='UNFORMATTED')                                             FOGCM.28
      DO I=1,NESTL                                                         FOGCM.29
       OPEN(UNIT=INEST(I),FILE=LABHISN(1,I)(1:LNBLNK(LABHISN(1,I))),       FOGCM.30
     +     FORM='UNFORMATTED')                                             FOGCM.31
c       OPEN(UNIT=NESTIN(I),FILE=NESTRD(I)(1:LNBLNK(NESTRD(I))),
c     +     FORM='UNFORMATTED')
c       OPEN(UNIT=NESTOT(I),FILE=NESTWT(I)(1:LNBLNK(NESTWT(I))),
c     +     FORM='UNFORMATTED')
      ENDDO                                                                FOGCM.32
      OPEN(11,FILE='Pot-Height.dat',FORM='UNFORMATTED')
      OPEN(12,FILE='Pre-Height.dat',FORM='UNFORMATTED')
      OPEN(13,FILE='Sig-Height.dat',FORM='UNFORMATTED')
      OPEN(14,FILE='Hgt-Height.dat',FORM='UNFORMATTED')   
C     ****                                                                 FOGCM.34
C     ****          ADVANCE THE MODEL IN TIME                              FOGCM.35
C     ****                                                                 FOGCM.36
      DO WHILE (ITER.LT.NSTP)                                              FOGCM.37
       CALL ADVNCE                                                         FOGCM.38
C      ****                                                                FOGCM.39
C      ****          CHECK FOR HISTORY FILE WRITE                          FOGCM.40
C      ****                                                                FOGCM.41
       IF(IHIS.EQ.0)CALL OUTPUT                                            FOGCM.42
      ENDDO                                                                FOGCM.43
      WRITE(6,"('NORMAL EXIT AT ITER = ',I6)")ITER                         FOGCM.44
      STOP                                                                 FOGCM.46
      END                                                                  FOGCM.47
C                                                                          FOGCM.48
*DK GETANS                                                                GETANS.1
      SUBROUTINE GETANS(ANS,NU)                                           GETANS.2
      CHARACTER*160 ANS                                                   GETANS.3
C     ****                                                                GETANS.4
C     ****     READ CHARACTER ANS FROM LOGICAL UNIT NU                    GETANS.5
C     ****                                                                GETANS.6
      READ(NU,"(A)")ANS                                                   GETANS.7
      RETURN                                                              GETANS.8
      END                                                                 GETANS.9
*DK GRIDINT                                                               GRDINT.1
      SUBROUTINE GRDINT(AM,AG,IG,JG,WT,LG,LM,IMAXG,JMAXG,IMAXM,JMAXM,     GRDINT.2
     1JM)                                                                 GRDINT.3
      IMPLICIT REAL*8 (A-H,O-Z)                                           GRDINT.4
C     ****                                                                GRDINT.5
C     ****     TRANSFORM SCALAR FIELD GIVEN ON GEOGRAPHIC GRID TO         GRDINT.6
C     ****       GEOMAGNETIC GRID USING INDICES AND WEIGHTS GENERATED     GRDINT.7
C     ****       BY GRDSTM.                                               GRDINT.8
C     ****                                                                GRDINT.9
C     ****     INPUT:                                                     GRDINT.10
C     ****      AG(LG,1) IS 2-DIM SCALAR FIELD TO BE TRANSFORMED TO       GRDINT.11
C     ****        GEOMAGNETIC GRID.                                       GRDINT.12
C     ****        NB: PERIODIC POINT IS REPEATED SO THAT                  GRDINT.13
C     ****        AG(IMAXG+1,JG) = AG(1,JG)                               GRDINT.14
C     ****        (LM.GE.IMAXG+1)                                         GRDINT.15
C     ****      IMAXG, JMAXG  DIMENSIONS OF GEOGRAPHIC GRID               GRDINT.16
C     ****      IMAXM, JMAXM  DIMENSIONS OF GEOMAGNETIC GRID              GRDINT.17
C     ****      LG IS FIRST DIMENSION OF AG IN CALLING PROGRAM            GRDINT.18
C     ****        (LG.GE.IMAXG+1)                                         GRDINT.19
C     ****      IG(LM,1) GIVES ROUNDED DOWN GEOGRAPHIC LONGITUDE INDEX    GRDINT.20
C     ****        FOR EACH GEOMAGNETIC GRID POINT                         GRDINT.21
C     ****      JG(LM,1) GIVES ROUNDED DOWN GEOGRAPHIC LATITUDE INDEX     GRDINT.22
C     ****        FOR EACH GEOMAGNETIC GRID POINT                         GRDINT.23
C     ****      WT(4,LM,1) WEIGHTS FOR 4 CORNERS OF EACH GEOMAGNETIC      GRDINT.24
C     ****        GRID ELEMENT                                            GRDINT.25
C     ****      LM IS FIRST DIMENSION OF ARRAYS AM, IG, JG, WT IN         GRDINT.26
C     ****        CALLING PROGRAM                                         GRDINT.27
C     ****        (LM.GE.IMAXM)                                           GRDINT.28
C     ****      IEM IS FLAG INDICATING TYPE OF TRANSFORMATION             GRDINT.29
C     ****        IEM = 1 FOR OUPUT FIELDS AT DIP EQUATOR                 GRDINT.30
C     ****        IEM = 0 FOR GLOBAL OUTPUT FIELDS                        GRDINT.31
C     ****      NOTE:  IG, JG, WT ARE PRODUCED BY PREVIOUS CALL TO        GRDINT.32
C     ****        GRDSTM.                                                 GRDINT.33
C     ****                                                                GRDINT.34
C     ****     OUTPUT:                                                    GRDINT.35
C     ****      AM(LM,1) 2-DIM SCALAR FIELD TRANSFORMED TO GEOMAGNETIC    GRDINT.36
C     ****        LATITUDE/LONGITUDE GRID                                 GRDINT.37
C     ****                                                                GRDINT.38
      DIMENSION AM(LM,1),AG(LG,1),IG(LM,1),JG(LM,1),WT(4,LM,1)            GRDINT.39
C     ****                                                                GRDINT.40
C     ****     CARRY OUT INTERPOLATION                                    GRDINT.41
C     ****                                                                GRDINT.42
      DO IM=1,IMAXM                                                       GRDINT.43
        AM(IM,1)=                                                         GRDINT.44
     1    AG(IG(IM,JM),JG(IM,JM))*WT(1,IM,JM)+                            GRDINT.45
     2    AG(IG(IM,JM)+1,JG(IM,JM))*WT(2,IM,JM)+                          GRDINT.46
     3    AG(IG(IM,JM)+1,JG(IM,JM)+1)*WT(3,IM,JM)+                        GRDINT.47
     4    AG(IG(IM,JM),JG(IM,JM)+1)*WT(4,IM,JM)                           GRDINT.48
      ENDDO                                                               GRDINT.49
      RETURN                                                              GRDINT.50
      END                                                                 GRDINT.51
C                                                                         GRDINT.52
*DK ICHK                                                                    ICHK.1
      LOGICAL FUNCTION ICHK(ITIME)                                          ICHK.2
C     ****                                                                  ICHK.3
C     ****     CHECK FOR TIME OUT OF RANGE                                  ICHK.4
C     ****                                                                  ICHK.5
      DIMENSION ITIME(3)                                                    ICHK.6
      ICHK=ITIME(1).LT.0.OR.ITIME(2).LT.0.OR.ITIME(2).GE.24.OR.             ICHK.7
     1     ITIME(3).LT.0.OR.ITIME(3).GE.60                                  ICHK.8
      RETURN                                                                ICHK.9
      END                                                                   ICHK.10
C                                                                           ICHK.11
*DK IGET                                                                    IGET.1
      INTEGER FUNCTION IGET(STR,IER)                                        IGET.2
      PARAMETER(TOOBIG=2147483647.,TOOSMALL=-2147483648.)                   IGET.3
      CHARACTER STR*(*)                                                     IGET.4
C     ****                                                                  IGET.5
C     ****     GET A INTEGER FROM STDIN:                                    IGET.6
C     ****     IF IER=0 ON RETURN, ALL'S WELL                               IGET.7
C     ****     IF IER=1 ON RETURN, INTEGER NUMBER WAS TOO BIG OR TOO SMALL  IGET.8
C     ****     IF IER=2 ON RETURN, A DAD STRING WAS READ (E.G. ALPHAS       IGET.9
C     ****                         INSTEAD OF DIGITS)                       IGET.10
C     ****                                                                  IGET.11
      IER=0                                                                 IGET.12
      READ(STR,*,ERR=900)R                                                  IGET.13
      IF(R.GT.TOOBIG.OR.R.LT.TOOSMALL)THEN                                  IGET.14
       IER=1                                                                IGET.15
       RETURN                                                               IGET.16
      ENDIF                                                                 IGET.17
      IGET=IFIX(R)                                                          IGET.18
      RETURN                                                                IGET.19
C     ****                                                                  IGET.20
C     ****     BAD STRING                                                   IGET.21
C     ****                                                                  IGET.22
900   IRE=2                                                                 IGET.23
      RETURN                                                                IGET.24
      END                                                                   IGET.25
C                                                                           IGET.26
*DK INIT                                                                    INIT.1
      SUBROUTINE INIT                                                       INIT.2
      INCLUDE 'param.h'                                                     INIT.3
      INCLUDE 'buffe.h'                                                     INIT.4
C     ****                                                                  INIT.5
C     ****                   INITIALIZE BUFFER INDICIES                     INIT.6
C     ****                                                                  INIT.7
      DO M=1,MODULO                                                         INIT.8
       NJTMP(M+1)=NBADDR(M)                                                 INIT.9
      ENDDO                                                                 INIT.10
      RETURN                                                                INIT.11
      END                                                                   INIT.12
C                                                                           INIT.13
*DK INPUT                                                                  INPUT.1
      SUBROUTINE INPUT                                                     INPUT.2
C     ****                                                                 INPUT.3
C     ****    THIS IS THE MAIN SUBROUTINE FOR INTERACTIVE INPUT OF CONTROL INPUT.4
C     ****    PARAMETERS FOR TIENG RUN ON WORKSTATIONS                     INPUT.5
C     ****    NOV., 1999       WENBIN WANG                                 INPUT.6
C     ****                                                                 INPUT.7
C     ****     LABEL='LABEL' (UP TO 160 CHARACTERS)                        INPUT.8
C     ****           LABEL FOR THIS RUN                                    INPUT.9
C     ****     DATE=IYEAR,IDAY                           (INTEGERS)        INPUT.10
C     ****          IYEAR AND IDAY ARE THE YEAR (1979, SAY) AND DAYNUMBER  INPUT.11
C     ****           OF THE YEAR FOR WHICH THE MODEL IS TO BE APPLIED.     INPUT.12
C     ****     SOURCE='PATH/FILENAME'                                      INPUT.13
C     ****          IDENTIFIES SOURCE VOLUME AND ITS PATH                  INPUT.14
C     ****     OUTPUT='PATH/FILE1',                                        INPUT.15
C     ****            'PATH/FILE2', - - - -                                INPUT.16
C     ****          IDENTIFICATION OF UP TO 10 OUTPUT VOLUMES TO BE USED   INPUT.17
C     ****          SEQUENTIALLY.  LIMIT OF 24 HISTORIES PER VOLUME.       INPUT.18
C     ****          EACH VOLUME IDENTIFIED AS DESCRIBED UNDER 'DATA'.      INPUT.19
C     ****     START=IDAYST,IHRST,MINST   (INTEGERS)                       INPUT.20
C     ****           STARTING TIME: DAY,HOUR,MINUTE                        INPUT.21
C     ****     STOP=IDAYND1,IHRND1,MINND1  (INTEGERS)                      INPUT.22
C     ****           STOPPING TIME                                         INPUT.23
C     ****     HIST=IDAYH1,IHRH1,MINH1   (INTEGERS)                        INPUT.24
C     ****           INTERVAL AT WHICH HISTORIES TO BE WRITTEN             INPUT.25
C     ****     STEP=ISTEP       (INTEGER)                                  INPUT.26
C     ****           STEP IN SECONDS                                       INPUT.27
C     ****     DIFHOR=DIFHOR           (INTEGER)                           INPUT.28
C     ****           DIFHOR = 1 IF HORIZONTAL VARIATION OF EDDY DIFFUSIO   INPUT.29
C     ****           IS TO BE CALCULATED BY CALLING USER PROVIDED          INPUT.30
C     ****           SUBROUTINE DFACT(T,J,SECS,IMX)                        INPUT.31
C     ****           DIFHOR = 0 FOR NO HORIZONTAL VARIATION OF EDDY DIFF   INPUT.32
C     ****     TIDE = Z22,Z23,Z24,Z25,Z26,PH22,PH23,PH24,PH25,PH26         INPUT.33
C     ****           (FLOATING)                                            INPUT.34
C     ****           Z22,Z23,Z24,Z25,Z26 ARE AMPLITUDES IN CM.,            INPUT.35
C     ****           PH22,PH23,PH24,PH25,PH26 ARE PHASES IN HRS.           INPUT.36
C     ****           OF 22,23,24,25,26 COMPONENTS.                         INPUT.37
C     ****           NOTE:  WHEN TIDAL BOUNDARY CONTITION IS USED, BOUND   INPUT.38
C     ****           FIELDS FOR T, U, V AND Z ARE PLOTTED.  TO LOOK AT T   INPUT.39
C     ****           PLOTS, A SUITABLE DISPOSE OF $PLT SHOULD BE INCLUDE   INPUT.40
C     ****           YOUR JCL.                                             INPUT.41
C     ****     TIDE2 = Z11, PH11     (FLOATING)                            INPUT.42
C     ****          Z11 = AMPLITUDE OF (1,1) TIDAL COMPONENT  (CM)         INPUT.43
C     ****          PH11 = PHASE OF (1,1) TIDAL COMPONENT  (HOURS)         INPUT.44
C     ****               -12..LE.PH11.LT.12.                               INPUT.45
C     ****                                                                 INPUT.46
C     ****     THE FOLLOWING PARAMETER READS ARE IN A SEPARATE READ FILE,  INPUT.47
C     ****     AND CAN BE TIME VARYING.                                    INPUT.48
C     ****                                                                 INPUT.49
C     ****          EACH INPUT IS DEFINED BY A TABLE OF VALUES FOR         INPUT.50
C     ****          MONOTONICALLY INCREASING TIMES.  VALUE FOR GIVEN       INPUT.51
C     ****          TIME IS FOUND BY LINEAR INTERPOLATION.  IF TIME IS     INPUT.52
C     ****          OUTSIDE TABLE, LINEAR EXTRAPOLATION OCCURS.  READFILE  INPUT.53
C     ****          ENTRY IS OF FOLLOWING FORM:                            INPUT.54
C     ****          PARAM=ID1,IH1,IM1,VAL1,ID2,IH2,IM2,VAL2, - - - - -     INPUT.55
C     ****          ,IDN,IHN,IMN,VALN                                      INPUT.56
C     ****          WHERE ID1=DAY1, IH1=HOUR1, IM1=MINUTE1                 INPUT.57
C     ****          ARE INTEGERS GIVING UNIVERSAL TIME, VAL1 IS            INPUT.58
C     ****          CORRESPONDING FLOATING POINT INPUT VALUE.  FOR         INPUT.59
C     ****          THIS FORM OF INPUT, N.GE.2, N.LE.50.                   INPUT.60
C     ****          PARAM=VAL DEFINES A CONSTANT VALUE.                    INPUT.61
C     ****                                                                 INPUT.62
C     ****     F107 = F107        (FLOATING)                               INPUT.63
C     ****          F107 IS THE DAILY 10.7 CM SOLAR FLUX VALUE, WHICH      INPUT.64
C     ****           IF VARYING, SHOULD BE CENTERED AT 12 UT.              INPUT.65
C     ****     F107A = F107A      (FLOATING)                               INPUT.66
C     ****          F107A IS THE 3 SOLAR ROTATION (81 DAY) AVERAGE 10.7    INPUT.67
C     ****           CM SOLAR FLUX.  IF VARYING, SHOULD BE CENTERED AT     INPUT.68
C     ****           12 UT.                                                INPUT.69
C     ****     AURORAL INPUTS                                              INPUT.70
C     ****                                                                 INPUT.71
C     ****      GENERAL INDICES OF AURORAL POWER, CONVECTION STRENGTH AND  INPUT.72
C     ****      CONVECTION ASSYMETRY, UPON WHICH ALL THE OTHER CONVECTION  INPUT.73
C     ****      AND OVAL CHARACTERISTICS CAN BE PARAMETERIZED IN TAIL.     INPUT.74
C     ****      1)  POWER      (FLOATING)                                  INPUT.75
C     ****           POWER IS THE HEMISPHERIC POWER IN GW, USUALLY TAKEN   INPUT.76
C     ****           FROM DAVE EVANS.                                      INPUT.77
C     ****      2)  CTPOTEN    (FLOATING)                                  INPUT.78
C     ****           CTPOTEN IS THE CROSS-TAIL POTENTIAL IN KV, USUALLY    INPUT.79
C     ****           INFERRED FROM THE REIFF AND LUHMANN (1986) FORMULA    INPUT.80
C     ****           OF VBSIN3(THETA/2).                                   INPUT.81
C     ****      3)  BYIMF      (FLOATING)                                  INPUT.82
C     ****           BYIMF IS THE GSM BY COMPONENT IN NT OF THE SOLAR IMF  INPUT.83
C     ****           MAGNETIC FIELD, USUALLY TAKEN FROM THE IMF TAPE.      INPUT.84
C     ****                                                                 INPUT.85
C     ****     START INPUT                                                 INPUT.86
C     ****                                                                 INPUT.87
      INCLUDE 'param.h'                                                    INPUT.88
      INCLUDE 'const.h'                                                    INPUT.89
      INCLUDE 'start.h'                                                    INPUT.90
      INCLUDE 'nestcnt.h'                                                  INPUT.91
C     ****                                                                 INPUT.92
C     ****     READ PARAMETERS FOR THE AURORAL INPUT                       INPUT.93
C     ****                                                                 INPUT.94
      PARAMETER (NAURPX=5,NLEX=600)                                        INPUT.95
      DIMENSION DUM(4,NLEX,NAURPX),ICNT(5)                                 INPUT.96
      CHARACTER*8 NAMEHA(NAURPX),X                                         INPUT.97
      REAL*8 MAG(2,2)                                                      INPUT.98
      COMMON/AURDAT0/AURPS(NLEX,2,NAURPX),PARAMV(NAURPX)                   INPUT.99
      COMMON/AURDAT1/NAURP,NPTS(NAURPX)                                    INPUT.100
C     ****                                                                 INPUT.101
C     ****     HEELIS CASE PLUS AURORA (ASYMMETRICAL)                      INPUT.102
C     ****                                                                 INPUT.103
      DATA NAMEHA/'    F107','   F107A',' H POWER',' C POTEN',             INPUT.104
     1            '  IMF BY'/                                              INPUT.105
C     ****                                                                 INPUT.106
C     ****   SET UP TIDES                                                  INPUT.107
C     ****                                                                 INPUT.108
      DIMENSION TIDEA(5,12),TIDEP(5,12)                                    INPUT.109
      CHARACTER*3 TEMPMONTH                                                INPUT.110
      DATA TIDEA/2.8E+4,0.7E+4,1.0E+4,1.5E+4,0.6E+4,                       INPUT.111
     1            2.2E+4,0.7E+4,1.1E+4,1.3E+4,0.3E+4,                      INPUT.112
     2            1.5E+4,0.1E+4,2.8E+4,0.1E+4,0.2E+4,                      INPUT.113
     3            1.8E+4,0.3E+4,2.2E+4,0.6E+4,0.3E+4,                      INPUT.114
     4            1.9E+4,0.3E+4,0.5E+4,0.9E+4,0.5E+4,                      INPUT.115
     5            2.3E+4,0.5E+4,0.4E+4,0.8E+4,0.4E+4,                      INPUT.116
     6            2.3E+4,0.8E+4,0.7E+4,1.1E+4,0.6E+4,                      INPUT.117
     7            2.2E+4,1.1E+4,0.5E+4,1.3E+4,0.6E+4,                      INPUT.118
     8            1.5E+4,0.4E+4,1.4E+4,0.9E+4,0.6E+4,                      INPUT.119
     9            1.7E+4,0.4E+4,1.8E+4,0.2E+4,0.4E+4,                      INPUT.120
     1            2.4E+4,0.4E+4,1.2E+4,1.0E+4,0.7E+4,                      INPUT.121
     1            2.7E+4,0.7E+4,0.9E+4,1.2E+4,0.5E+4/                      INPUT.122
      DATA TIDEP/10.4,1.2,7.8,9.0,8.3,                                     INPUT.123
     1            10.2,0.8,8.1,8.8,7.8,                                    INPUT.124
     2            10.3,5.8,7.6,8.3,8.4,                                    INPUT.125
     3            10.2,5.6,7.7,1.4,6.2,                                    INPUT.126
     4            10.3,8.5,7.7,2.2,7.4,                                    INPUT.127
     5            10.5,8.6,6.6,3.0,8.3,                                    INPUT.128
     6            10.2,10.5,3.8,3.2,8.1,                                   INPUT.129
     7            10.2,10.7,2.1,2.7,7.4,                                   INPUT.130
     8            10.4,1.6,8.2,1.3,0.8,                                    INPUT.131
     9            10.3,3.4,7.7,4.0,2.9,                                    INPUT.132
     1            10.5,1.8,8.0,8.5,6.8,                                    INPUT.133
     1            10.5,0.9,7.4,8.8,8.1/                                    INPUT.134
      COMMON/TIDES/TIDE(10),TIDE2(2)                                       INPUT.135
      CHARACTER*160 LABEL,ANS,PNAME                                        INPUT.136
      LOGICAL ICHK,ISTPCK                                                  INPUT.137
      CHARACTER*60 SOURCE,OUTPUT(20)                                       INPUT.138
      INTEGER START(3),STOP(3),HIST(3),STEP,DIFHOR                         INPUT.139
C     ****                                                                 INPUT.140
C     ****    SHOW THE VERSION OF TIGCM                                    INPUT.141
C     ****                                                                 INPUT.142
C      CALL VERSION                                                        INPUT.143
      NU=89                                                                INPUT.144
      OPEN(NU,FILE='inputdata',STATUS='OLD')                               INPUT.145
C     ****                                                                 INPUT.146
C     ****      INITIALIZE                                                 INPUT.147
C     ****                                                                 INPUT.148
      DO I=1,160                                                           INPUT.149
       LABEL(I:I)=' '                                                      INPUT.150
      ENDDO                                                                INPUT.151
      DO I=1,5                                                             INPUT.152
       ICNT(I)=0                                                           INPUT.153
      ENDDO                                                                INPUT.154
      DO I=1,3                                                             INPUT.155
       ISOURCE(I)=0                                                        INPUT.156
      ENDDO                                                                INPUT.157
      DO I=1,NAURPX*NLEX*2                                                 INPUT.159
       AURPS(I,1,1)=0.                                                     INPUT.160
      ENDDO                                                                INPUT.161
      ICHAGN=0                                                             INPUT.162
      ICHVOL=0                                                             INPUT.163
      MFIL=24                                                              INPUT.164
C     ****                                                                 INPUT.165
C     ****    OPEN FILE TO STORE INPUT PARAMETERS                          INPUT.166
C     ****                                                                 INPUT.167
7001  WRITE(6,"(///,'****   File Name Input   ****')")                     INPUT.168
      WRITE(6,"('Please input a file name to hold input parameters',       INPUT.169
     +            ' like ddmmyyuser01 ')")                                 INPUT.170
      WRITE(6,"(' Enter your file name: ',$)")                             INPUT.171
      CALL GETANS(ANS,NU)                                                  INPUT.172
      NPNAME=LNBLNK(ANS)                                                   INPUT.173
      DO I=1,160                                                           INPUT.174
       PNAME(I:I)=' '                                                      INPUT.175
      ENDDO                                                                INPUT.176
      DO I=1,NPNAME                                                        INPUT.177
       PNAME(I:I)=ANS(I:I)                                                 INPUT.178
      ENDDO                                                                INPUT.179
      WRITE(6,"(/,'Your input parameters wrill be stored in file: ',       INPUT.180
     +           A)")PNAME(1:NPNAME)                                       INPUT.181
      WRITE(6,"(' Type C or c to continue, Type R or r to rename the ',    INPUT.182
     +            'file: ',$)")                                            INPUT.183
      CALL GETANS(ANS,NU)                                                  INPUT.184
      IF(ANS(1:1).EQ.'R'.OR.ANS(1:1).EQ.'r')GOTO 7001                      INPUT.185
      IF(ANS(1:1).EQ.' '.OR.ANS(1:1).EQ.'C'.OR.ANS(1:1).EQ.'c')GOTO 7002   INPUT.186
C     ****                                                                 INPUT.187
C     ****     LABEL INPUT                                                 INPUT.188
C     ****                                                                 INPUT.189
7002  WRITE(6,"(///,'****    Label Input   ****')")                        INPUT 190
      WRITE(6,"(' Please input the LABEL of the run, up to 160',           INPUT.191
     +         ' characters')")                                            INPUT.192
      CALL GETANS(ANS,NU)                                                  INPUT.193
      NLABEL=LNBLNK(ANS)                                                   INPUT.194
      DO I=1,160                                                           INPUT.195
       LABEL(I:I)=' '                                                      INPUT.196
      ENDDO                                                                INPUT.197
      DO I=1,NLABEL                                                        INPUT.198
       LABEL(I:I)=ANS(I:I)                                                 INPUT.199
      ENDDO                                                                INPUT.200
      IF(NLABEL.GT.160)THEN                                                INUPT.201
       WRITE(6,"(' INPUT -- LABEL: MORE THAN 160 CHARACTERS')")            INPUT.202
       GOTO 7002                                                           INPUT.203
      ENDIF                                                                INPUT.204
      IF(ICHAGN.EQ.1)THEN                                                  INPUT.205
       GOTO 7030                                                           INPUT.206
      ENDIF                                                                INPUT.207
C     ****                                                                 INPUT.208
C     ****     DATE                                                        INPUT.209
C     ****                                                                 INPUT.210
7003  WRITE(6,"(///,'****   Date Input   ****')")                          INPUT.211
      WRITE(6,"(' Please enter the YEAR for which the model is to be',     INPUT.212
     +' applied: ',$)")                                                    INPUT.213
      CALL GETANS(ANS,NU)                                                  INPUT.214
      IYEAR=IGET(ANS,IER)                                                  INPUT.215
      WRITE(6,"(' Please enter the DAY for which the model is to be',      INPUT.216
     +' applied: ',$)")                                                    INPUT.217
      CALL GETANS(ANS,NU)                                                  INPUT.218
      IIDAY=IGET(ANS,IER)                                                  INPUT.219
      IF(IYEAR.LT.1900.OR.IYEAR.GT.3000.OR.IIDAY.LT.1.OR.IIDAY.GT.366)     INPUT.220
     1  THEN                                                               INPUT.221
       WRITE(6,"(' INPUT -- DATE: PARAMETER OUT OF RANGE')")               INPUT.222
       GOTO 7003                                                           INPUT.223
      ENDIF                                                                INPUT.224
      IF(ICHAGN.EQ.1)THEN                                                  INPUT.225
       GOTO 7030                                                           INPUT.226
      ENDIF                                                                INPUT.227
C     ****                                                                 INPUT.228
C     ****     STEP                                                        INPUT.229
C     ****                                                                 INPUT.230
7006  WRITE(6,"(///,'****   Time Step Input   ****')")                     INPUT.231
      WRITE(6,"(' Please Enter STEP of this run in seconds :',$)")         INPUT.232
      CALL GETANS(ANS,NU)                                                  INPUT.233
      STEP=IGET(ANS,IER)                                                   INPUT.234
      IF(STEP.LE.0)THEN                                                    INPUT.235
       WRITE(6,"(' INPUT -- STEP: STEP .LE. 0')")                          INPUT.236
       GOTO 7006                                                           INPUT.237
      ENDIF                                                                INPUT.248
      C(4)=STEP                                                            INPUT.249
      IF(ICHAGN.EQ.1.OR.ICHVOL.EQ.1)THEN                                   INPUT.250
       GOTO 7008                                                           INPUT.251
      ENDIF                                                                INPUT.252
C     ****                                                                 INPUT.253
C     ****     START                                                       INPUT.254
C     ****                                                                 INPUT.255
7007  WRITE(6,"(///,'****    Start Time Input   ****')")                   INPUT.256
      WRITE(6,"(' Please enter START time of the run: ',                   INPUT.257
     + /'-- (seperate day,hour and minute with comma): ',$)")              INPUT.258
      CALL GETANS(ANS,NU)                                                  INPUT.259
      IF(ANS(1:1).EQ.' ')THEN                                              INPUT.260
       WRITE(6,"('There is no default -- Please try again')")              INPUT.261
       GOTO 7007                                                           INPUT.262
      ENDIF                                                                INPUT.263
C     ****                                                                 INPUT.264
C     ****     Find Commas                                                 INPUT.265
C     ****                                                                 INPUT.266
      ICOMMA1=0                                                            INPUT.267
      ICOMMA2=0                                                            INPUT.268
      DO I=1,160                                                           INPUT.269
       IF(ANS(I:I).EQ.',')THEN                                             INPUT.270
        IF(ICOMMA1.EQ.0) ICOMMA1=I                                         INPUT.271
        IF(ICOMMA1.GT.0) ICOMMA2=I                                         INPUT.272
       ENDIF                                                               INPUT.273
      ENDDO                                                                INPUT.274
      IF(ICOMMA1.EQ.O.OR.ICOMMA2.EQ.0)THEN                                 INPUT.275
       WRITE(6,"('>>>Please seperate day, hour and minute by ',            INPUT.276
     +       'commas <<<')")                                               INPUT.277
       GOTO 7007                                                           INPUT.278
      ENDIF                                                                INPUT.279
C     ****                                                                 INPUT.280
C     ****     Get Day, Hour, Minute                                       INPUT.281
C     ****                                                                 INPUT.282
      NSTART=LNBLNK(ANS)                                                   INPUT.283
      START(1)=IGET(ANS(1:ICOMMA1-1),IER)                                  INPUT.284
      START(2)=IGET(ANS(ICOMMA1+1:ICOMMA2-1),IER)                          INPUT.285
      START(3)=IGET(ANS(ICOMMA2+1:NSTART),IER)                             INPUT.286
      IF(START(1).GT.365.OR.START(1).LT.0.OR.START(2).GT.24.               INPUT.287
     +   OR.START(2).LT.0.OR.START(3).GT.60.OR.START(3).LT.0)THEN          INPUT.288
       WRITE(6,"('>>> Invalid input of day or hour or minute <<<')")       INPUT.289
       WRITE(6,"('Please try again')")                                     INPUT.290
       GOTO 7007                                                           INPUT.291
      ENDIF                                                                INPUT.292
      IF(ICHK(START))THEN                                                  INPUT.293
       WRITE(6,"(' INPUT -- START: TIME OUT OF RANGE')")                   INPUT.294
       GOTO 7007                                                           INPUT.295
      ENDIF                                                                INPUT.296
      IF(ISTPCK(START,STEP))THEN                                           INPUT.297
       WRITE(6,"(' INPUT -- START: TIME NOT A MULTIPLE OF STEP')")         INPUT.298
       GOTO 7007                                                           INPUT.299
      ENDIF                                                                INPUT.300
7008  ITER=ITERA(START,STEP)                                               INPUT.301
      IF(ICHAGN.EQ.1.OR.ICHVOL.EQ.1)THEN                                   INPUT.302
       GOTO 7009                                                           INPUT.303
      ENDIF                                                                INPUT.304
C     ****                                                                 INPUT.305
C     ****     SOURCE                                                      INPUT.306
C     ****                                                                 INPUT.307
7009  WRITE(6,"(///,'****   Source Input   ****')")                        INPUT.308
7010  CONTINUE                                                             INPUT.309
      DO I=1,60                                                            INPUT.310
       SOURCE(I:I)=' '                                                     INPUT.311
      ENDDO                                                                INPUT.312
      WRITE(6,"(' Do you want to start with your source file(y/n)? :'      INPUT.313
     +,$)")                                                                INPUT.314
      CALL GETANS(ANS,NU)                                                  INPUT.315
      IF(ANS(1:1).EQ.'Y'.OR.ANS(1:1).EQ.'y')THEN                           INPUT.316
       ISTRT=1                                                             INPUT.317
       WRITE(6,"(' Please enter your SOURCE filename: ',$)")               INPUT.318
       CALL GETANS(ANS,NU)                                                 INPUT.319
       NSOURCE=LNBLNK(ANS)                                                 INPUT.320
       SOURCE(1:NSOURCE)=ANS(1:NSOURCE)                                    INPUT.321
      ELSE                                                                 INPUT.322
       ISTRT=0                                                             INPUT.323
       WRITE(6,"(' Please enter HISTORY filename: ',$)")                   INPUT.324
       CALL GETANS(ANS,NU)                                                 INPUT.325
       NSOURCE=LNBLNK(ANS)                                                 INPUT.326
       SOURCE(1:NSOURCE)=ANS(1:NSOURCE)                                    INPUT.327
      ENDIF                                                                INPUT.328
      NSOURCE=LNBLNK(SOURCE)                                               INPUT.329
      DO I=1,NSOURCE                                                       INPUT.330
       LABSRC(I:I)=SOURCE(I:I)                                             INPUT.331
      ENDDO                                                                INPUT.332
      DO I=1,3                                                             INPUT.333
       ISOURCE(I)=START(I)                                                 INPUT.334
      ENDDO                                                                INPUT.335
      IF(ICHAGN.EQ.1.OR.ICHVOL.EQ.1)THEN                                   INPUT.336
       GOTO 7013                                                           INPUT.337
      ENDIF                                                                INPUT.338
C     ****                                                                 INPUT.339
C     ****     OUTPUT                                                      INPUT.340
C     ****                                                                 INPUT.341
7011  WRITE(6,"(///,'****    Output File Input   ****')")                  INPUT.342
      WRITE(6,"(' Please enter your OUTPUT filename: ',$)")                INPUT.343
      CALL GETANS(ANS,NU)                                                  INPUT.344
      NOUTPUT=LNBLNK(ANS)                                                  INPUT.345
      DO J=1,20                                                            INPUT.346
       DO I=1,60                                                           INPUT.347
        OUTPUT(J)(I:I)=' '                                                 INPUT.348
       ENDDO                                                               INPUT.349
      ENDDO                                                                INPUT.350
      OUTPUT(1)(1:NOUTPUT)=ANS(1:NOUTPUT)                                  INPUT.351
      OUTPUT(1)(NOUTPUT+1:NOUTPUT+6)='XX.dat'                              INPUT.352
      WRITE(6,"('Your OUTPUT file will be : ',a)")OUTPUT(1)                INPUT.353
      WRITE(6,"('XX is the volume number with maximum of 20')")            INPUT.354
      NOUTPUT=LNBLNK(OUTPUT(1))                                            INPUT.355
      DO I=1,NOUTPUT                                                       INPUT.356
       LABHIS(1)(I:I)=OUTPUT(1)(I:I)                                       INPUT.357
      ENDDO                                                                INPUT.358
      NOUTPUTN=LNBLNK(ANS)                                                 INPUT.359
      LABHISN(1,1)(1:2)='n1'                                               INPUT.360
      LABHISN(1,1)(3:NOUTPUTN+2)=ANS(1:NOUTPUT)                            INPUT.361
      LABHISN(1,1)(3+NOUTPUTN:NOUTPUTN+8)='XX.dat'                         INPUT.362
      NOUTPUTN=NOUTPUTN+8                                                  INPUT.363
      IF(ICHAGN.EQ.1)THEN                                                  INPUT.364
       GOTO 7029                                                           INPUT.365
      ENDIF                                                                INPUT.366
C     ****                                                                 INPUT.367
C     ****     STOP                                                        INPUT.368
C     ****                                                                 INPUT.369
7012  WRITE(6,"(///,'****    Stop Time Input   ****')")                    INPUT.370
      WRITE(6,"(' Please enter STOP time of the run: ',                    INPUT.371
     + /'-- (seperate day,hour and minute with comma): ',$)")              INPUT.372
      CALL GETANS(ANS,NU)                                                  INPUT.373
      IF(ANS(1:1).EQ.' ')THEN                                              INPUT.374
       WRITE(6,"('There is no default -- Please try again')")              INPUT.375
       GOTO 7012                                                           INPUT.376
      ENDIF                                                                INPUT.377
C     ****                                                                 INPUT.378
C     ****     Find Commas                                                 INPUT.379
C     ****                                                                 INPUT.380
      ICOMMA1=0                                                            INPUT.381
      ICOMMA2=0                                                            INPUT.382
      DO I=1,160                                                           INPUT.383
       IF(ANS(I:I).EQ.',')THEN                                             INPUT.384
        IF(ICOMMA1.EQ.0) ICOMMA1=I                                         INPUT.385
        IF(ICOMMA1.GT.0) ICOMMA2=I                                         INPUT.386
       ENDIF                                                               INPUT.387
      ENDDO                                                                INPUT.388
      IF(ICOMMA1.EQ.O.OR.ICOMMA2.EQ.0)THEN                                 INPUT.389
       WRITE(6,"('>>>Please seperate day, hour and minute by ',            INPUT.390
     +       'commas <<<')")                                               INPUT.391
       GOTO 7012                                                           INPUT.392
      ENDIF                                                                INPUT.393
C     ****                                                                 INPUT.394
C     ****     Get Day, Hour, Minute                                       INPUT.395
C     ****                                                                 INPUT.396
      NSTOP=LNBLNK(ANS)                                                    INPUT.397
      STOP(1)=IGET(ANS(1:ICOMMA1-1),IER)                                   INPUT.398
      STOP(2)=IGET(ANS(ICOMMA1+1:ICOMMA2-1),IER)                           INPUT.399
      STOP(3)=IGET(ANS(ICOMMA2+1:NSTOP),IER)                               INPUT.400
      IF(STOP(1).GT.365.OR.STOP(1).Lt.0.OR.STOP(2).GT.24.                  INPUT.401
     +   OR.STOP(2).LT.0.OR.STOP(3).GT.60.OR.STOP(3).LT.0)THEN             INPUT.402
       WRITE(6,"('>>> Invalid input of day or hour or minute <<<')")       INPUT.403
       WRITE(6,"('Please try again')")                                     INPUT.404
       GOTO 7012                                                           INPUT.405
      ENDIF                                                                INPUT.406
      IF(ICHK(STOP))THEN                                                   INPUT.407
       WRITE(6,"(' INPUT -- STOP: TIME OUT OF RANGE')")                    INPUT.408
       GOTO 7012                                                           INPUT.409
      ENDIF                                                                INPUT.410
7013  IF(ISTPCK(STOP,STEP))THEN                                            INPUT.421
       WRITE(6,"(' INPUT -- STOP: TIME',I1,' NOT A MULTIPLE OF STEP')")    INPUT.412
       GOTO 7012                                                           INPUT.413
      ENDIF                                                                INPUT.414
      NNSTP=ITERA(STOP,STEP)                                               INPUT.415
      NSTP=NNSTP                                                           INPUT.416
C     ****                                                                 INPUT.417
C     ****     ITER.LT.NNSTP                                               INPUT.418
C     ****                                                                 INPUT.419
      IF(ITER.GE.NNSTP)THEN                                                INPUT.420
       WRITE(6,"(' INPUT -- MUST HAVE START TIME .LT. STOP TIME')")        INPUT.421
      GOTO 7012                                                            INPUT.422
      ENDIF                                                                INPUT.423
      IF(I2CHAGN.EQ.1.OR.ICHVOL.EQ.1)THEN                                  INPUT.424
       GOTO 7015                                                           INPUT.425
      ENDIF                                                                INPUT.426
C     ****                                                                 INPUT.427
C     ****     HIST                                                        INPUT.428
C     ****                                                                 INPUT.429
7014  WRITE(6,"(///,'****    History Time Input   ****')")                 INPUT.430
      WRITE(6,"(' Please enter HISTORY time of the run: ',                 INPUT.431
     + /'-- (seperate day,hour and minute with comma): ',$)")              INPUT.432
      CALL GETANS(ANS,NU)                                                  INPUT.433
      IF(ANS(1:1).EQ.' ')THEN                                              INPUT.434
       WRITE(6,"('There is no default -- Please try again')")              INPUT.435
       GOTO 7014                                                           INPUT.436
      ENDIF                                                                INPUT.437
C     ****                                                                 INPUT.438
C     ****     Find Commas                                                 INPUT.439
C     ****                                                                 INPUT.440
      ICOMMA1=0                                                            INPUT.441
      ICOMMA2=0                                                            INPUT.442
      DO I=1,160                                                           INPUT.443
       IF(ANS(I:I).EQ.',')THEN                                             INPUT.444
        IF(ICOMMA1.EQ.0) ICOMMA1=I                                         INPUT.445
        IF(ICOMMA1.GT.0) ICOMMA2=I                                         INPUT.446
       ENDIF                                                               INPUT.447
      ENDDO                                                                INPUT.448
      IF(ICOMMA1.EQ.O.OR.ICOMMA2.EQ.0)THEN                                 INPUT.449
       WRITE(6,"('>>>Please seperate day, hour and minute by ',            INPUT.450
     +       'commas <<<')")                                               INPUT.451
       GOTO 7014                                                           INPUT.452
      ENDIF                                                                INPUT.453
C     ****                                                                 INPUT.454
C     ****     Get Day, Hour, Minute                                       INPUT.455
C     ****                                                                 INPUT.456
      NHIST1=LNBLNK(ANS)                                                   INPUT.457
      HIST(1)=IGET(ANS(1:ICOMMA1-1),IER)                                   INPUT.458
      HIST(2)=IGET(ANS(ICOMMA1+1:ICOMMA2-1),IER)                           INPUT.459
      HIST(3)=IGET(ANS(ICOMMA2+1:NHIST1),IER)                              INPUT.460
      IF(HIST(1).GT.365.OR.HIST(1).LT.0.OR.HIST(2).GT.24.                  INPUT.461
     +   OR.HIST(2).LT.0.OR.HIST(3).GT.60.OR.HIST(3).LT.0)THEN             INPUT.462
       WRITE(6,"('>>> Invalid input of day or hour or minute <<<')")       INPUT.463
       WRITE(6,"('Please try again')")                                     INPUT.464
       GOTO 7014                                                           INPUT.465
      ENDIF                                                                INPUT.466
      IF(ICHK(HIST))THEN                                                   INPUT.467
       WRITE(6,"(' INPUT -- HIST: TIME OUT OF RANGE')")                    INPUT.468
       GOTO 7014                                                           INPUT.469
      ENDIF                                                                INPUT.470
7015  IF(ISTPCK(HIST,STEP))THEN                                            INPUT.471
       WRITE(6,"(' INPUT -- HIST: TIME',I1,' NOT A MULTIPLE OF STEP')")    INPUT.472
       GOTO 7014                                                           INPUT.473
      ENDIF                                                                INPUT.474
      MHIS=ITERA(HIST,STEP)                                                INPUT.475
      IF(ICHAGN.EQ.1.OR.ICHVOL.EQ.1)THEN                                   INPUT.476
       GOTO 7028                                                           INPUT.477
      ENDIF                                                                INPUT.478
C     ****                                                                 INPUT.479
C     ****     DIFHOR                                                      INPUT.480
C     ****                                                                 INPUT.481
7018  WRITE(6,"(///,'****   Horizontal Diffusion Flag   ****')")           INPUT.482
      WRITE(6,"(' Please Enter Horizontal Diffusion Flag  ',               INPUT.483
     +/'-- 0--> no varaiation',/'-- 1--> variation',                       INPUT.484
     +/'-- (default is 1): ',$)")                                          INPUT.485
      CALL GETANS(ANS,NU)                                                  INPUT.486
      IF(ANS(1:1).EQ.' ')THEN                                              INPUT.487
       DIFHOR=1                                                            INPUT.488
      ELSE                                                                 INPUT.489
       DIFHOR=IGET(ANS,IER)                                                INPUT.490
      ENDIF                                                                INPUT.491
      IF(DIFHOR.NE.0.AND.DIFHOR.NE.1)THEN                                  INPUT.492
       WRITE(6,"(' INPUT -- DIFHOR: MUST BE 0 OR 1')")                     INPUT.493
       GOTO 7018                                                           INPUT.494
      ENDIF                                                                INPUT.495
      ICOMP=DIFHOR                                                         INPUT.496
      IF(ICHAGN.EQ.1)THEN                                                  INPUT.497
       GOTO 7030                                                           INPUT.498
      ENDIF                                                                INPUT.499
C     ****                                                                 INPUT.500
C     ****     TIDE                                                        INPUT.501
C     ****                                                                 INPUT.502
      OPEN (7,FILE='dateUT.dat')                                           INPUT.503
      READ(7,7700)TEMPMONTH                                                INPUT.504
      CLOSE(7)                                                             INPUT.505
7700  FORMAT(4X,A3,21X)                                                    INPUT.506
      IF(TEMPMONTH.EQ.'Jan')THEN                                           INPUT.507
       DO N=1,5                                                            INPUT.508
        TIDE(N)=TIDEA(N,1)                                                 INPUT.509
        TIDE(N+5)=TIDEP(N,1)                                               INPUT.510
       ENDDO                                                               INPUT.511
      ELSEIF(TEMPMONTH.EQ.'Feb')THEN                                       INPUT.512
       DO N=1,5                                                            INPUT.513
        TIDE(N)=TIDEA(N,2)                                                 INPUT.514
        TIDE(N+5)=TIDEP(N,2)                                               INPUT.515
       ENDDO                                                               INPUT.516
      ELSEIF(TEMPMONTH.EQ.'Mar')THEN                                       INPUT.517
       DO N=1,5                                                            INPUT.518
        TIDE(N)=TIDEA(N,3)                                                 INPUT.519
        TIDE(N+5)=TIDEP(N,3)                                               INPUT.520
       ENDDO                                                               INPUT.521
      ELSEIF(TEMPMONTH.EQ.'Apr')THEN                                       INPUT.522
       DO N=1,5                                                            INPUT.523
        TIDE(N)=TIDEA(N,4)                                                 INPUT.524
        TIDE(N+5)=TIDEP(N,4)                                               INPUT.525
       ENDDO                                                               INPUT.526
      ELSEIF(TEMPMONTH.EQ.'May')THEN                                       INPUT.527
       DO N=1,5                                                            INPUT.528
        TIDE(N)=TIDEA(N,5)                                                 INPUT.529
        TIDE(N+5)=TIDEP(N,5)                                               INPUT.530
       ENDDO                                                               INPUT.531
      ELSEIF(TEMPMONTH.EQ.'Jun')THEN                                       INPUT.532
       DO N=1,5                                                            INPUT.533
        TIDE(N)=TIDEA(N,6)                                                 INPUT.534
        TIDE(N+5)=TIDEP(N,6)                                               INPUT.535
       ENDDO                                                               INPUT.536
      ELSEIF(TEMPMONTH.EQ.'Jul')THEN                                       INPUT.537
       DO N=1,5                                                            INPUT.538
        TIDE(N)=TIDEA(N,7)                                                 INPUT.539
        TIDE(N+5)=TIDEP(N,7)                                               INPUT.540
       ENDDO                                                               INPUT.541
      ELSEIF(TEMPMONTH.EQ.'Aug')THEN                                       INPUT.542
       DO N=1,5                                                            INPUT.543
        TIDE(N)=TIDEA(N,8)                                                 INPUT.544
        TIDE(N+5)=TIDEP(N,8)                                               INPUT.545
       ENDDO                                                               INPUT.546
      ELSEIF(TEMPMONTH.EQ.'Sep')THEN                                       INPUT.547
       DO N=1,5                                                            INPUT.548
        TIDE(N)=TIDEA(N,9)                                                 INPUT.549
        TIDE(N+5)=TIDEP(N,9)                                               INPUT.550
       ENDDO                                                               INPUT.551
      ELSEIF(TEMPMONTH.EQ.'Oct')THEN                                       INPUT.552
       DO N=1,5                                                            INPUT.553
        TIDE(N)=TIDEA(N,10)                                                INPUT.554
        TIDE(N+5)=TIDEP(N,10)                                              INPUT.555
       ENDDO                                                               INPUT.556
      ELSEIF(TEMPMONTH.EQ.'Nov')THEN                                       INPUT.557
       DO N=1,5                                                            INPUT.558
        TIDE(N)=TIDEA(N,11)                                                INPUT.559
        TIDE(N+5)=TIDEP(N,11)                                              INPUT.560
       ENDDO                                                               INPUT.561
      ELSEIF(TEMPMONTH.EQ.'Dec')THEN                                       INPUT.562
       DO N=1,5                                                            INPUT.563
        TIDE(N)=TIDEA(N,12)                                                INPUT.564
        TIDE(N+5)=TIDEP(N,12)                                              INPUT.565
       ENDDO                                                               INPUT.566
      ENDIF                                                                INPUT.567
      ITIDE=0                                                              INPUT.568
      DO I=1,5                                                             INPUT.569
       IF(TIDE(I).GT.0)ITIDE=1                                             INPUT.570
      ENDDO                                                                INPUT.571
C     ****                                                                 INPUT.572
C     ****     TIDE2                                                       INPUT.573
C     ****                                                                 INPUT.574
      TIDE2(1)=20000.                                                      INPUT.575
      TIDE2(1)=18.                                                         INPUT.576
      ITIDE2=0                                                             INPUT.577
      DO N=1,1                                                             INPUT.578
       IF(TIDE2(N).GT.0.)ITIDE2=1                                          INPUT.579
      ENDDO                                                                INPUT.580
C     ****                                                                 INPUT.581
C     ****     COLFAC                                                      INPUT.582
C     ****                                                                 INPUT.583
      COLFAC=1.5                                                           INPUT.584
C     ****                                                                 INPUT.576
C     ****     NEST GRID                                                   INPUT.577
C     ****                                                                 INPUT.578
7024  WRITE(6,"(///,'****   Nested Grid   ****')")                         INPUT.579
      NEST=0                                                               INPUT.580
      WRITE(6,"(' NESTED GRID Function Enabled (default is 0',
     1            ' not enabled, 1 -->enabled <RETURN>) :',$)")            INPUT.582
      CALL GETANS(ANS,NU)                                                  INPUT.583
      NEST=IGET(ANS,IER)
      IF(NEST.GT.1.OR.NEST.LT.0)THEN                                       INPUT.584
       WRITE(6,*)'INPUT SHOULD BE EITHER 1 OR 0, SELECT AGAIN'             INPUT.585
       GOTO 7024
      ENDIF                                                                INPUT.586
      IF(ICHAGN.EQ.1)THEN                                                  INPUT.587
       GOTO 7030                                                           INPUT.588
      ENDIF                                                                INPUT.589
C     ****                                                                 INPUT.585
C     ****     MAG                                                         INPUT.586
C     ****                                                                 INPUT.587
      MAG(1,1)=-74.5                                                       INPUT.588
      MAG(2,1)=127.                                                        INPUT.589
      MAG(1,2)=79.                                                         INPUT.590
      MAG(2,2)=-70.                                                        INPUT.591
      DO I=1,2                                                             INPUT.592
       DIPOLE(I)=90.0+(3-2*I)*MAG(1,I)                                     INPUT.593
       PMLONG(I)=MAG(2,I)+(I-2)*180.0                                      INPUT.594
      ENDDO                                                                INPUT.595
C     ****                                                                 START.120
C     ****     SPECIFY NESTED GRID DOMAIN PARAMETERS                       START.121
C     ****                                                                 START.122
      IF(NEST.EQ.1)THEN
       CALL NESTREAD
      ENDIF
C     ****                                                                 INPUT.596
C     ****     CHECK THAT SUFFICIENT VOLUMES FOR OUTPUT                    INPUT.597
C     ****                                                                 INPUT.598
7028  NVOL=(NNSTP-ITER)/MHIS                                               INPUT.599
      NVOL=(NVOL+MFIL)/MFIL                                                INPUT.600
      ICHVOL=0                                                             INPUT.601
      IF(NVOL.GT.10)THEN                                                   INPUT.602
       WRITE(6,"(/////,' Output volume exceeds 10')")                      INPUT.603
       WRITE(6,"('--- Please select the parameter you want to change :'    INPUT.604
     +      ,$)")                                                          INPUT.605
       WRITE(6,*)                                                          INPUT.606
       WRITE(6,"('1). STOP TIME')")                                        INPUT.607
       WRITE(6,"('2). START TIME')")                                       INPUT.608
       WRITE(6,"('3). STEP')")                                             INPUT.609
       WRITE(6,"('4). HISTORY TIME')")                                     INPUT.610
       WRITE(6,"('Please enter your choice :',$)")                         INPUT.611
       CALL GETANS(ANS,NU)                                                 INPUT.612
       IN=IGET(ANS,IER)                                                    INPUT.613
       IF(IN.LT.1.OR.IN.GT.4)THEN                                          INPUT.614
        WRITE(6,"(' Not a valid selection, try again')")                   INPUT.615
        GOTO 7028                                                          INPUT.616
       ELSE IF(IN.EQ.1)THEN                                                INPUT.617
        ICHVOL=1                                                           INPUT.618
        GOTO 7012                                                          INPUT.619
       ELSE IF(IN.EQ.2)THEN                                                INPUT.620
        ICHVOL=1                                                           INPUT.621
        GOTO 7007                                                          INPUT.622
       ELSE IF(IN.EQ.3)THEN                                                INPUT.623
        ICHVOL=1                                                           INPUT.624
        GOTO 7006                                                          INPUT.625
       ELSE IF(IN.EQ.4)THEN                                                INPUT.626
        ICHVOL=1                                                           INPUT.627
        GOTO 7014                                                          INPUT.628
       ENDIF                                                               INPUT.629
      ENDIF                                                                INPUT.630
C     ****                                                                 INPUT.631
C     ****     DEFINE OUTPUT NAME                                          INPUT.632
C     ****                                                                 INPUT.633
7029  DO II=2,NVOL                                                         INPUT.634
       DO I=1,NOUTPUT                                                      INPUT.635
        LABHIS(II)(I:I)=LABHIS(1)(I:I)                                     INPUT.636
       ENDDO                                                               INPUT.637
       DO I=1,NOUTPUTN                                                     INPUT.638
        DO L=1,NESTL                                                       INPUT.639
         LABHISN(II,L)(I:I)=LABHISN(1,1)(I:I)                              INPUT.640
         WRITE(X,8001)L                                                    INPUT.641
         LABHISN(II,L)(2:2)=X                                              INPUT.642
        ENDDO                                                              INPUT.643
       ENDDO                                                               INPUT.644
       IF(II.LE.9)THEN                                                     INPUT.645
        WRITE(X,8001)II                                                    INPUT.646
        LABHIS(II)(NOUTPUT-5:NOUTPUT-5)='0'                                INPUT.647
        LABHIS(II)(NOUTPUT-4:NOUTPUT-4)=X                                  INPUT.648
        DO L=1,NESTL                                                       INPUT.649
         LABHISN(II,L)(NOUTPUTN-5:NOUTPUTN-5)='0'                          INPUT.650
         LABHISN(II,L)(NOUTPUTN-4:NOUTPUTN-4)=X                            INPUT.651
        ENDDO                                                              INPUT.652
       ELSE                                                                INPUT.653
        WRITE(X,8002)II                                                    INPUT.654
        LABHIS(II)(NOUTPUT-5:NOUTPUT-4)=X                                  INPUT.655
        DO L=1,NESTL                                                       INPUT.656
         LABHISN(II,L)(NOUTPUTN-5:NOUTPUTN-4)=X                            INPUT.657
        ENDDO                                                              INPUT.658
       ENDIF                                                               INPUT.659
      ENDDO                                                                INPUT.661
      LABHIS(1)(NOUTPUT-5:NOUTPUT-5)='0'                                   INPUT.662
      LABHIS(1)(NOUTPUT-4:NOUTPUT-4)='1'                                   INPUT.663
      LABHISN(1,1)(NOUTPUTN-5:NOUTPUTN-5)='0'                              INPUT.664
      LABHISN(1,1)(NOUTPUTN-4:NOUTPUTN-4)='1'                              INPUT.665
      DO L=2,NESTL
       LABHISN(1,L)=LABHISN(1,1)
      ENDDO
      LABHISN(1,2)(2:2)='2'
      LABHISN(1,3)(2:2)='3'
      LABHISN(1,4)(2:2)='4'
8001  FORMAT(I1)                                                           INPUT.666
8002  FORMAT(I2)                                                           INPUT.667
      OPEN(7,FILE=PNAME,STATUS='UNKNOWN')                                  INPUT.668
7030  CONTINUE                                                             INPUT.669
      WRITE(6,"(///////,'----- Those Are Your Input Parameters For ',      INPUT.670
     +          'This Run -----')")                                        INPUT.671
      WRITE(6,*)                                                           INPUT.672
      WRITE(6,"('1). LABEL   = ',A)")LABEL(1:LNBLNK(LABEL))                INPUT.673
      WRITE(6,"('2). DATE    = ',I4,'(YEAR), ',I3,'(DAT)')")IYEAR,IIDAY    INPUT.674
      WRITE(6,"('4). STEP    = ',I3)")STEP                                 INPUT.675
      DO I=1,NVOL                                                          INPUT.676
       WRITE(6,"('5). OUTPUT  = ',A)")LABHIS(I)(1:LNBLNK(OUTPUT))          INPUT.677
      ENDDO                                                                INPUT.678
      DO I=1,NVOL                                                          INPUT.679
       DO L=1,NESTL                                                        INPUT.680
        WRITE(6,"('5). OUTPUT  = ',A)")                                    INPUT.681
     1            LABHISN(I,L)(1:LNBLNK(LABHISN(I,L)))                     INPUT.682
       ENDDO                                                               INPUT.683
      ENDDO                                                                INPUT.684
      WRITE(6,"('6). SOURCE  = ',A)")SOURCE(1:LNBLNK(SOURCE))              INPUT.685
      WRITE(6,"('8). START   = ',I3,'(DAYS)',I2,'(HOURS)',I2,'(MINUTES)'   INPUT.686
     +     )")START(1),START(2),START(3)                                   INPUT.687
      WRITE(6,"('9). STOP    = ',I3,'(DAYS)',I2,'(HOURS)',I2,'(MINUTES)'   INPUT.688
     +     )")STOP(1),STOP(2),STOP(3)                                      INPUT.689
      WRITE(6,"('10). HIST   = ',I3,'(DAYS)',I2,'(HOURS)',I2,'(MINUTES)'   INPUT.690
     +     )")HIST(1),HIST(2),HIST(3)                                      INPUT.691
      WRITE(6,"('12). DIFHOR = ',I1)")DIFHOR                               INPUT.692
      WRITE(6,"('14). TIDES  = ',5F10.2,'(AMPLITUDES)')")TIDE(1),          INPUT.693
     +     TIDE(2),TIDE(3),TIDE(4),TIDE(5)                                 INPUT.694
      WRITE(6,"('              ',5F10.2,'(PHASES)')")TIDE(6),TIDE(7),      INPUT.695
     +     TIDE(8),TIDE(9),TIDE(10)                                        INPUT.696
      WRITE(6,"('15). TIDE2  = ', 2F10.2)")TIDE2(1),TIDE2(2)               INPUT.697
      WRITE(6,"('16). COLFAC = ',F10.2)")COLFAC                            INPUT.698
      WRITE(6,*)                                                           INPUT.699
      WRITE(6,"(' Do You Want To Change Any of Those Parameters? :',$)")   INPUT.700
      CALL GETANS(ANS,NU)                                                  INPUT.701
      IF(ANS(1:1).EQ.'Y'.OR.ANS(1:1).EQ.'y')THEN                           INPUT.702
       WRITE(6,"(' Enter the No. of the Parameter You Want to Change :',   INPUT.703
     +      $)")                                                           INPUT.704
       CALL GETANS(ANS,NU)                                                 INPUT.705
       ICHECK=IGET(ANS,IER)                                                INPUT.706
       ICHAGN=1                                                            INPUT.707
       IF(ICHECK.EQ.1)THEN                                                 INPUT.708
        GOTO 7002                                                          INPUT.709
       ELSE IF(ICHECK.EQ.2)THEN                                            INPUT.710
        GOTO 7003                                                          INPUT.711
       ELSE IF(ICHECK.EQ.4)THEN                                            INPUT.712
        GOTO 7006                                                          INPUT.713
       ELSE IF(ICHECK.EQ.5)THEN                                            INPUT.714
        GOTO 7011                                                          INPUT.715
       ELSE IF(ICHECK.EQ.6)THEN                                            INPUT.716
        GOTO 7009                                                          INPUT.717
       ELSE IF(ICHECK.EQ.8)THEN                                            INPUT.718
        GOTO 7007                                                          INPUT.719
       ELSE IF(ICHECK.EQ.9)THEN                                            INPUT.720
        GOTO 7012                                                          INPUT.721
       ELSE IF(ICHECK.EQ.10)THEN                                           INPUT.722
        GOTO 7014                                                          INPUT.723
       ELSE IF(ICHECK.EQ.12)THEN                                           INPUT.724
        GOTO 7018                                                          INPUT.725
       ENDIF                                                               INPUT.726
      ENDIF                                                                INPUT.727
      WRITE(6,"('Saving the parameters to ',A)")PNAME(1:LNBLNK(PNAME))     INPUT.728
      WRITE(6,"('****************************************')")              INPUT.729
      WRITE(7,"('1). LABEL   = ',A)")LABEL(1:LNBLNK(LABEL))                INPUT.730
      WRITE(7,"('2). DATE    = ',I4,'(YEAR), ',I3,'(DAT)')")IYEAR,IIDAY    INPUT.731
      WRITE(7,"('4). STEP    = ',I4)")STEP                                 INPUT.732
      DO I=1,NVOL                                                          INPUT.733
       WRITE(7,"('5). OUTPUT  = ',A)")LABHIS(I)(1:LNBLNK(OUTPUT))          INPUT.734
      ENDDO                                                                INPUT.735
      DO I=1,NVOL                                                          INPUT.736
       DO L=1,NESTL                                                        INPUT.737
        WRITE(7,"('5). OUTPUT  = ',A)")                                    INPUT.738
     1             LABHISN(I,L)(1:LNBLNK(LABHISN(I,L)))                    INPUT.739
       ENDDO                                                               INPUT.740
      ENDDO                                                                INPUT.741
      WRITE(7,"('6). SOURCE  = ',A)")SOURCE(1:LNBLNK(SOURCE))              INPUT.742
      WRITE(7,"('8). START   = ',I3,'(DAYS)',I2,'(HOURS)',I2,'(MINUTES)'   INPUT.743
     +     )")START(1),START(2),START(3)                                   INPUT.744
      WRITE(7,"('9). STOP    = ',I3,'(DAYS)',I2,'(HOURS)',I2,'(MINUTES)'   INPUT.745
     +     )")STOP(1),STOP(2),STOP(3)                                      INPUT.746
      WRITE(7,"('10). HIST   = ',I3,'(DAYS)',I2,'(HOURS)',I2,'(MINUTES)'   INPUT.757
     +     )")HIST(1),HIST(2),HIST(3)                                      INPUT.748
      WRITE(7,"('12). DIFHOR = ',I1)")DIFHOR                               INPUT.749
      WRITE(7,"('14). TIDES  = ',5F10.2,'(AMPLITUDES)')")TIDE(1),          INPUT.750
     +     TIDE(2),TIDE(3),TIDE(4),TIDE(5)                                 INPUT.751
      WRITE(7,"('             ',5F10.2,'(PHASES)')")TIDE(6),TIDE(7),       INPUT.752
     +     TIDE(8),TIDE(9),TIDE(10)                                        INPUT.753
      WRITE(7,"('15). TIDE2  = ', 2F10.2)")TIDE2(1),TIDE2(2)               INPUT.754
      WRITE(7,"('16). COLFAC = ',F7.2)")COLFAC                             INPUT.755
      CLOSE(7)                                                             INPUT.756
      CLOSE(NU)                                                            INPUT.757
C     ****                                                                 INPUT.758
C     ****     RESET READ DICTIONARY                                       INPUT.759
C     ****                                                                 INPUT.760
      NAURP=NAURPX                                                         INPUT.761
      DO I=1,NAURP                                                         INPUT.762
       ICNT(I)=0                                                           INPUT.763
      ENDDO                                                                INPUT.764
      LUNIT=7                                                              INPUT.775
      OPEN(LUNIT,FILE='dinput.dat',STATUS='old')                           INPUT.766
C     ****                                                                 INPUT.767
C     ****     READ IN SOLAR EMISSION AND AURORAL PARAMETERS               INPUT.768
C     ****                                                                 INPUT.769
      DO J=1,NAURP                                                         INPUT.770
       READ(LUNIT,"(I7)")ICNT(J)                                           INPUT.771
       IF(ICNT(J).EQ.1)THEN                                                INPUT.772
        READ(LUNIT,"(I7,F6.1,59X,A18)")JJJ,DUM(1,1,J),COMMENT              INPUT.773
       ELSE                                                                INPUT.774
        DO K=1,ICNT(J)/4                                                   INPUT.775
         READ(LUNIT,"(I7,4F6.1,31X,A18)")JJJ,DUM(1,K,J),DUM(2,K,J)         INPUT.776
     +        ,DUM(3,K,J),DUM(4,K,J),COMMENT                               INPUT.777
        ENDDO                                                              INPUT.778
       ENDIF                                                               INPUT.779
      ENDDO                                                                INPUT.780
      CLOSE(LUNIT)                                                         INPUT.781
C     ****                                                                 INPUT.782
C     ****     INTERPOLATION OF AURORAL PARAMETERS IN TIME                 INPUT.783
C     ****                                                                 INPUT.784
100   DO N=1,NAURP                                                         INPUT.785
       NPTS(N)=(ICNT(N)+3)/4                                               INPUT.786
       IF (NPTS(N) .GT. NLEX) THEN                                         INPUT.787
        WRITE (6,"(1X,'INPUT -- ',A8,' NPTS =',I3,                         INPUT.788
     1      ' WHICH IS .GT. NLEX')")NAMEHA(N),NPTS(N)                      INPUT.789
        STOP                                                               INPUT.790
       ENDIF                                                               INPUT.791
      ENDDO                                                                INPUT.792
C     ****                                                                 INPUT.793
C     ****     CHECK NUMBER OF DATA POINTS                                 INPUT.794
C     ****                                                                 INPUT.795
      DO N=1,NAURP                                                         INPUT.796
       II=ICNT(N)                                                          INPUT.797
       IF (II.GT.1.AND.II.LT.8.AND.MOD(II,4).NE.0) THEN                    INPUT.798
        WRITE(6,"(1X,'INPUT -- ',A8,' NO. OF ENTRIES MUST BE 1 OR ',       INPUT.799
     1   'MULTIPLE OF 4 AND GE 8.  ICNT NPTS =', 2I5)")NAMEHA(N),          INPUT.800
     2   ICNT(N),NPTS(N)                                                   INPUT.801
        IERROR=1                                                           INPUT.802
       ENDIF                                                               INPUT.803
      ENDDO                                                                INPUT.804
C     ****                                                                 INPUT.805
C     ****     ASSIGN DATA POINTS INTO COMMON BLLOCKS                      INPUT.806
C     ****                                                                 INPUT.807
      DO N=1,NAURP                                                         INPUT.808
       IF (NPTS(N).EQ.1)THEN                                               INPUT.809
        AURPS(1,2,N)=DUM(1,1,N)                                            INPUT.800
       ELSE                                                                INPUT.811
        DO M=1,NPTS(N)                                                     INPUT.812
         AURPS(M,1,N)=DUM(1,M,N)*24.+DUM(2,M,N)+DUM(3,M,N)/60.             INPUT.813
         AURPS(M,2,N)=DUM(4,M,N)                                           INPUT.814
        ENDDO                                                              INPUT.815
       ENDIF                                                               INPUT.816
      ENDDO                                                                INPUT.817
C     ****                                                                 INPUT.818
C     ****     TIMES FOR AURORAL PARAMETERS MUST BE MONOTONIC              INPUT.819
C     ****                                                                 INPUT.820
      DO 240 N=1,NAURP                                                     INPUT.821
       IF (NPTS(N).EQ.1)GO TO 240                                          INPUT.822
       DO M=2,NPTS(N)                                                      INPUT.823
        IF (AURPS(M,1,N).LE.AURPS(M-1,1,N))THEN                            INPUT.824
         WRITE (6,"(1X,'INPUT --',A8,':  TABLE OF TIMES NOT ',             INPUT.825
     1    'MONOTONIC')")NAMEHA(N)                                          INPUT.826
         WRITE (6,"(1X,'NPTS N,N+1 =', 3I4, '  D H M (N,N+1)', 6I3)")      INPUT.827
     1    NPTS(N),M-1,M,(DUM(K1,M-1,N),K1=1,3),(DUM(K2,M,N),K2=1,3)        INPUT.828
         IERROR = 1                                                        INPUT.829
        ENDIF                                                              INPUT.830
       ENDDO                                                               INPUT.831
  240 CONTINUE                                                             INPUT.832
C     ****                                                                 INPUT.833
C     ****   TERMINATE IF ERROR FLAG SET                                   INPUT.834
C     ****                                                                 INPUT.835
      IF (IERROR.EQ.1) STOP 'INPUT'                                        INPUT.836
      RETURN                                                               INPUT.837
      END                                                                  INPUT.838
C                                                                          INPUT.839
*DK ISTPCK                                                                ISTPCK.1
      LOGICAL FUNCTION ISTPCK(ITIME,ISTEP)                                ISTPCK.2
C     ****                                                                ISTPCK.3
C     ****     CHECK FOR TIME NOT A MULTIPLE OF ISTEP                     ISTPCK.4
C     ****                                                                ISTPCK.5
      DIMENSION ITIME(3)                                                  ISTPCK.6
      ISTPCK=MOD(((ITIME(1)*24+ITIME(2))*60+ITIME(3))*60,ISTEP).          ISTPCK.7
     1       NE.0                                                         ISTPCK.8
      RETURN                                                              ISTPCK.9
      END                                                                 ISTPCK.10
C                                                                         ISTPCK.11
*DK ITERA                                                                  ITERA.1
      FUNCTION ITERA(ITIME,ISTEP)                                          ITERA.2
C     ****                                                                 ITERA.3
C     ****     CALCULATES TIME STEP NUMBER (ITERATION)                     ITERA.4
C     ****                                                                 ITERA.5
      DIMENSION ITIME(3)                                                   ITERA.6
      ITERA=((ITIME(1)*24+ITIME(2))*60+ITIME(3))*60/ISTEP                  ITERA.7
      RETURN                                                               ITERA.8
      END                                                                  ITERA.9
C                                                                          ITERA.10
*DK LNBLNK                                                                LNBLNK.1
      FUNCTION LNBLNK(STR)                                                LNBLNK.2
      CHARACTER*(*) STR                                                   LNBLNK.3
C     ****                                                                LNBLNK.4
C     ****     RETURN NUMBER OF CHARACTERS UP TO FIRST OCCURRENCE OF A    LNBLNK.5
C     ****     BLANK                                                      LNBLNK.6
C     ****                                                                LNBLNK.7
      LNBLNK=0                                                            LNBLNK.8
      LENSTR=LEN(STR)                                                     LNBLNK.9
      DO I=1,LENSTR                                                       LNBLNK.10
       IF(STR(I:I+1).EQ.' ')THEN                                          LNBLNK.11
        LNBLNK=I-1                                                        LNBLNK.12
        RETURN                                                            LNBLNK.13
       ENDIF                                                              LNBLNK.14
      ENDDO                                                               LNBLNK.15
      LNBLNK=LENSTR                                                       LNBLNK.16
      RETURN                                                              LNBLNK.17
      END                                                                 LNBLNK.18
*DK MAGDYN                                                                MAGDYN.1
      SUBROUTINE MAGDYN                                                   MAGDYN.2
C     ****                                                                MAGDYN.3
C     ****     SET UP VERSIONS OF /TRGM/ AND /MAGFLD/ CALCULATED          MAGDYN.4
C     ****       FROM TRUE FIELD AVAILABLE IN /FIELD/                     MAGDYN.5
C     ****                                                                MAGDYN.6
      INCLUDE 'param.h'                                                   MAGDYN.7
      INCLUDE 'const.h'                                                   MAGDYN.8
      INCLUDE 'econs.h'                                                   MAGDYN.9
      INCLUDE 'epara.h'
      include 'efild.h'                                                   MAGDYN.10
      INTEGER ZIMXP2, ZJMXP2                                              MAGDYN.11
      PARAMETER(ZIMXP2=ZIMX+2,ZJMXP2=ZJMX+2)                              MAGDYN.12
      COMMON/MAGFLD/BX(-1:ZIMXP2,-1:ZJMXP2),BY(-1:ZIMXP2,-1:ZJMXP2),      MAGDYN.13
     1  BZ(-1:ZIMXP2,-1:ZJMXP2),BMOD2(-1:ZIMXP2,-1:ZJMXP2)                MAGDYN.14
      COMMON/TRGM/RLATM(ZIMXP,ZJMX),RLONM(ZIMXP,ZJMX),                    MAGDYN.15
     1  DIPMAG(ZIMXP,ZJMX),DECMAG(ZIMXP,ZJMX),SNDEC(ZIMXP,ZJMX),          MAGDYN.16
     2  CSDEC(ZIMXP,ZJMX),SN2DEC(ZIMXP,ZJMX),SNCSDC(ZIMXP,ZJMX),          MAGDYN.17
     3  DUMDUM(ZJMX,3)                                                    MAGDYN.18
C     ****                                                                MAGDYN.19
C     ****     FILL /TGCM/ AND /MAGFLD/                                   MAGDYN.20
C     ****                                                                MAGDYN.21
      SIN10=1.D-6                                                         MAGDYN.22
      COS10=DSQRT(1.-SIN10**2)                                            MAGDYN.23
      DIPMIN=3.0                                                          MAGDYN.24
      SINDIP=DSIN(DIPMIN*DTR)                                             MAGDYN.25
      DO J=1,JMAX                                                         MAGDYN.26
       JJ=5+9*(J-1)
       DO I=1,IMAX                                                        MAGDYN.27
        II=(I-1)*5+1
        RLATM(I+2,J)=ALATM(II,JJ)                                           MAGDYN.28
        RLONM(I+2,J)=ALONM(II,JJ)                                           MAGDYN.29
        DIPMAG(I+2,J)=DATAN(ZZB(II,JJ)/DSQRT(XB(II,JJ)**2+YB(II,JJ)**2))        MAGDYN.30
C       ****                                                              MAGDYN.31
C       ****     AlLTER DIPMAG SO THAT THE SIN(DIPMAG) NEVER GOES BELOW   MAGDYN.32
C       ****     0.01 (I.E. MINIMUM DIPMAG IS 0.573 DEGREES -- AFFECTS    MAGDYN.33
C       ****     6 LOCATIONS --12.5 TO 12.5 DEGREES). SINDIP IS USED IN   MAGDYN.34
C       ****     LAMDAS, AND SETTEI.                                      MAGDYN.35
C       ****                                                              MAGDYN.36
        SINDIP=DSIN(DIPMAG(I+2,J))                                        MAGDYN.37
        IF(ABS(DIPMAG(I+2,J))/DTR.LT.DIPMIN)THEN                          MAGDYN.38
         IF(DIPMAG(I+2,J).GE.0.)DIPMAG(I+2,J)=DIPMIN*DTR                  MAGDYN.39
         IF(DIPMAG(I+2,J).LT.0.)DIPMAG(I+2,J)=-DIPMIN*DTR                 MAGDYN.40
        ENDIF                                                             MAGDYN.41
        DECMAG(I+2,J)=-DATAN2(YB(II,JJ),XB(II,JJ))                            MAGDYN.42
        SNDEC(I+2,J)=DSIN(DECMAG(I+2,J))                                  MAGDYN.43
        CSDEC(I+2,J)=DCOS(DECMAG(I+2,J))                                  MAGDYN.44
        SN2DEC(I+2,J)=SNDEC(I+2,J)**2                                     MAGDYN.45
        SNCSDC(I+2,J)=SNDEC(I+2,J)*CSDEC(I+2,J)                           MAGDYN.46
        BX(I,J)=YB(II,JJ)/BMOD(II,JJ)                                         MAGDYN.47
        BY(I,J)=XB(II,JJ)/BMOD(II,JJ)                                         MAGDYN.48
        BZ(I,J)=-ZZB(II,JJ)/BMOD(II,JJ)                                       MAGDYN.49
        BMOD2(I,J)=BMOD(II,JJ)                                              MAGDYN.50
C       ****                                                              MAGDYN.51
C       ****          SET MINIMUM DIP TO 10 DEGREES                       MAGDYN.52
C       ****                                                              MAGDYN.53
        DIF=DABS(BZ(I,J))-SIN10                                           MAGDYN.54
        IF(DIF.LT.0.)THEN                                                 MAGDYN.55
         BX(I,J)=BX(I,J)*COS10/DSQRT(1.-BZ(I,J)**2)                       MAGDYN.56
         BY(I,J)=BY(I,J)*COS10/DSQRT(1.-BZ(I,J)**2)                       MAGDYN.57
         BZ(I,J)=DSIGN(SIN10,BZ(I,J))                                     MAGDYN.58
        ENDIF                                                             MAGDYN.59
       ENDDO                                                              MAGDYN.60
      ENDDO                                                               MAGDYN.61
C     ****                                                                MAGDYN.62
C     ****     VALUES AT J = -1, 0, ZJMAXP1, ZJNMP2                       MAGDYN.63
C     ****                                                                MAGDYN.64
      DO J=1,2                                                            MAGDYN.65
       DO I=1,IMAX                                                        MAGDYN.66
        BX(I,J-2)=-BX(1+MOD(I-1+IMAX/2,IMAX),3-J)                         MAGDYN.67
        BY(I,J-2)=-BY(1+MOD(I-1+IMAX/2,IMAX),3-J)                         MAGDYN.68
        BZ(I,J-2)=BZ(1+MOD(I-1+IMAX/2,IMAX),3-J)                          MAGDYN.69
        BMOD2(I,J-2)=BMOD2(1+MOD(I-1+ZIMX/2,ZIMX),3-J)                    MAGDYN.70
        BX(I,JMAX+J)=-BX(1+MOD(I-1+IMAX/2,IMAX),JMAX+1-J)                 MAGDYN.71
        BY(I,JMAX+J)=-BY(1+MOD(I-1+IMAX/2,IMAX),JMAX+1-J)                 MAGDYN.72
        BZ(I,JMAX+J)=BZ(1+MOD(I-1+IMAX/2,IMAX),JMAX+1-J)                  MAGDYN.73
        BMOD2(I,ZJMX+J)=BMOD2(1+MOD(I-1+ZIMX/2,ZIMX),ZJMX+1-J)            MAGDYN.74
       ENDDO                                                              MAGDYN.75
      ENDDO                                                               MAGDYN.76
C     ****                                                                MAGDYN.77
C     ****     PERIODIC POINTS                                            MAGDYN.78
C     ****                                                                MAGDYN.79
      DO I=1,2                                                            MAGDYN.80
       DO J=1,8*JMAX                                                      MAGDYN.81
        RLATM(I,J)=RLATM(I+IMAX,J)                                        MAGDYN.82
        RLATM(I+IMAXP2,J)=RLATM(I+2,J)                                    MAGDYN.83
       ENDDO                                                              MAGDYN.84
      ENDDO                                                               MAGDYN.85
      DO I=1,2                                                            MAGDYN.86
       DO J=-1,4*(JMAX+4)-2                                               MAGDYN.87
        BX(I-2,J)=BX(I-2+IMAX,J)                                          MAGDYN.88
        BX(I+IMAX,J)=BX(I,J)                                              MAGDYN.89
       ENDDO                                                              MAGDYN.90
      ENDDO                                                               MAGDYN.91
      RETURN                                                              MAGDYN.92
      END                                                                 MAGDYN.93
C                                                                         MAGDYN.94
*DK NESTREAD                                                            NESTREAD.1
      SUBROUTINE NESTREAD                                               NESTREAD.2
      IMPLICIT REAL*8 (A-H,O-Z)                                         NESTREAD.3
      INCLUDE 'nestcnt.h'                                               NESTREAD.4
      INCLUDE 'start.h'                                                 NESTREAD.5
      DIMENSION NTIME(4)                                                NESTREAD.6
      OPEN(12,FILE='nestsrt.dat')                                       NESTREAD.7
      READ(12,1001)NLEVEL                                               NESTREAD.8
      NESTL=NLEVEL                                                      NESTREAD.9
      READ(12,1001)ISTEP                                                NESTREAD.10
1001  FORMAT(I6)                                                        NESTREAD.11
1003  FORMAT(4I6)                                                       NESTREAD.12
1002  FORMAT(5I6)                                                       NESTREAD.13
      DO N=1,NLEVEL                                                     NESTREAD.14
       READ(12,1002)NESTRI(N),NESTJ1(N),NESTJM(N),NESTI1(N),NESTIM(N)   NESTREAD.15
       NTIME(1)=1                                                       NESTREAD.16
       IF(N.GT.1)THEN                                                   NESTREAD.17
        DO I=2,N                                                        NESTREAD.18
         NTIME(N)=NTIME(N-1)*NESTRI(N-1)                                NESTREAD.19
        ENDDO                                                           NESTREAD.20
       ENDIF                                                            NESTREAD.21
       READ(12,1003)JJ,MSDAY(N),MSHOR(N),MSMIN(N)                       NESTREAD.22
       READ(12,1003)JJ,MEDAY(N),MEHOR(N),MEMIN(N)                       NESTREAD.23
       NES(N)=((MSDAY(N)*24+MSHOR(N))*60+MSMIN(N))*60*NTIME(N)/ISTEP    NESTREAD.24
       NEE(N)=((MEDAY(N)*24+MEHOR(N))*60+MEMIN(N))*60*NTIME(N)/ISTEP    NESTREAD.25
      ENDDO                                                             NESTREAD.26
      READ(12,1003)MHIST(1),MHIST(2),MHIST(3),MHIST(4)                  NESTREAD.27
      PHIS(1)=-87.5+DBLE(NESTJ1(1)-1)*5.                                NESTREAD.28
      SLAM(1)=-180.+DBLE(NESTI1(1)-3)*5.                                NESTREAD.29
      IF(NLEVEL.GE.2)THEN                                               NESTREAD.30
       DO N=2,NLEVEL                                                    NESTREAD.31
        RATIO=1.                                                        NESTREAD.32
        DO I=1,N-1                                                      NESTREAD.33
         RATIO=RATIO*NESTRI(I)                                          NESTREAD.34
        ENDDO                                                           NESTREAD.35
        PHIS(N)=PHIS(N-1)+DBLE(NESTJ1(N)-1)*5./RATIO                    NESTREAD.36
        SLAM(N)=SLAM(N-1)+DBLE(NESTI1(N)-1)*5./RATIO                    NESTREAD.37
       ENDDO                                                            NESTREAD.38
      ENDIF                                                             NESTREAD.39
      READ(12,1003)NFDBK(1),NFDBK(2),NFDBK(3),NFDBK(4)
      CLOSE(12)                                                         NESTREAD.40
      RETURN                                                            NESTREAD.41
      END                                                               NESTREAD.42
*DK OUTPUT                                                                OUTPUT.1
      SUBROUTINE OUTPUT                                                   OUTPUT.2
C     ****                                                                OUTPUT.3
C     ****     WRITES A HISTORY FILE                                      OUTPUT.4
C     ****                                                                OUTPUT.5
      INCLUDE 'param.h'                                                   OUTPUT.6
      INCLUDE 'index.h'                                                   OUTPUT.7
      INCLUDE 'buffe.h'                                                   OUTPUT.8
      INCLUDE 'const.h'                                                   OUTPUT.9
      INCLUDE 'start.h'                                                   OUTPUT.10
      INCLUDE 'epara.h'                                                   OUTPUT.11
      INCLUDE 'edphi.h'                                                   OUTPUT.12
      PARAMETER (NAURPX=5,NLEX=600)                                       OUTPUT.13
      COMMON/AURDAT0/AURPS(NLEX,2,NAURPX),PARAMV(NAURPX)                  OUTPUT.14
      COMMON/AURDAT1/NAURP,NPTS(NAURPX)                                   OUTPUT.15
      COMMON /IONCR/ ISTAR,IHEM,THETA0(2),OFFA(2),OFFC(2),DSKOFA(2),      OUTPUT.16
     1DSKOFC(2),PHID(2),PHIN(2),PHIDP0(2),PHIDM0(2),PHINP0(2),PHINM0(2),  OUTPUT.17
     2 PSIM(2),PSIE(2),PCEN(2),R1(2),DDUUMM(IMAXMP*JMX0+IMX0*JMX0+2)      OUTPUT.18
      COMMON /OVALR/ RRAD(2),H0,RH,RROTH,E0,REE,RROTE,FC,ALFAC,FD,ALFAD   OUTPUT.19
     1, ALFK,ALF6,ALF21,RROT6,RROT21,RD6,RD6V,RD21,RH6,RH21,RT6,RT21      OUTPUT.20
     2, ALFA0,RALFA,ALFA20,RALFA2,E20,RE2                                 OUTPUT.21
      COMMON/MASS/RMN4S,RMN2D,RMNO,BRN2D,COLFAC                           OUTPUT.22
      COMMON/INSOL/ F107,F107A                                            OUTPUT.23
      REAL*4 IHEAD(108)                                                   OUTPUT.24
      DATA IHEAD/108*0./                                                  OUTPUT.25
      DATA NH/1/                                                          OUTPUT.26
C     ****                                                                OUTPUT.27
C     ****     LOAD HEADER                                                OUTPUT.28
C     ****                                                                OUTPUT.29
      IHEAD(1)=DBLE(ITER)                                                 OUTPUT.30
      X=ITER*C(4)+1.D-4                                                   OUTPUT.31
      IHEAD(2)=IDINT(X/86400.)                                            OUTPUT.32
      X=X-IHEAD(2)*86400.                                                 OUTPUT.33
      IHEAD(3)=IDINT(X/3600.)                                             OUTPUT.34
      X=X-IHEAD(3)*3600.                                                  OUTPUT.35
      IHEAD(4)=IDINT(X/60.)                                               OUTPUT.36
      NPAR1=4+1                                                           OUTPUT.37
      IHEAD(NPAR1)=NAURP                                                  OUTPUT.38
      DO N=1,NAURP                                                        OUTPUT.39
        IHEAD(NPAR1+N)=PARAMV(N)                                          OUTPUT.40
      ENDDO                                                               OUTPUT.41
      NTH=30                                                              OUTPUT.42
      DO N=1,NTH                                                          OUTPUT.43
       IHEAD(NPAR1+NAURP+N)=THETA0(N)                                     OUTPUT.44
      ENDDO                                                               OUTPUT.45
      NROT0=NPAR1+NAURP+NTH                                               OUTPUT.46
      NROT=30                                                             OUTPUT.47
      DO N=1,NROT                                                         OUTPUT.48
       IHEAD(NROT0+N)=RRAD(N)                                             OUTPUT.49
      ENDDO                                                               OUTPUT.50
C     ****                                                                OUTPUT.51
C     ****     SAVE MODEL NUMBER, O+O COLLISION FACTOR AND F107 FLUX IN   OUTPUT.52
C     ****     THE HEADER.                                                OUTPUT.53
C     ****                                                                OUTPUT.54
      NEX1=NROT0+NROT+1                                                   OUTPUT.55
      IHEAD(NEX1)=4.40                                                    OUTPUT.56
      IHEAD(NEX1+1)=COLFAC                                                OUTPUT.57
      IHEAD(NEX1+2)=F107                                                  OUTPUT.58
      IHEAD(NEX1+3)=F107A                                                 OUTPUT.59
      IHEAD(108)=IFIL                                                     OUTPUT.60
C     ****                                                                OUTPUT.61
C     ****     WRITE HISTORY FILE                                         OUTPUT.62
C     *****                                                               OUTPUT.63
      IFIL=IFIL+1                                                         OUTPUT.64
      WRITE(IOTPT)(IHEAD(I),I=1,108)                                      OUTPUT.65
C     ****                                                                OUTPUT.66
C     ****     GENERAL LINE LOOP                                          OUTPUT.67
C     ****                                                                OUTPUT.68
      DO J=1,JMAX                                                         OUTPUT.69
       WRITE(IOTPT)((G(I,K,J),I=1,IMAXP4),K=1,KMAXP1*NHIST+1)             OUTPUT.89
      ENDDO                                                               OUTPUT.90
      IHIS=MHIS                                                           OUTPUT.91
C     ****                                                                OUTPUT.92
C     ****     CHECK FOR FULL HISTORY VOLUME                              OUTPUT.93
C     ****                                                                OUTPUT.94
      IF(IFIL.LT.MFIL)THEN                                                OUTPUT.95
       RETURN                                                             OUTPUT.96
C      ****                                                               OUTPUT.97
C      ****     CHECK FOR END OF RUN                                      OUTPUT.98
C      ****                                                               OUTPUT.99
      ELSE IF(ITER.GE.NSTP)THEN                                           OUTPUT.100
       RETURN                                                             OUTPUT.101
      ELSE                                                                OUTPUT.102
       CLOSE(IOTPT)                                                       OUTPUT.103
       ITAP=ITAP+1                                                        OUTPUT.104
       OPEN(UNIT=IOTPT,FILE=LABHIS(ITAP)(1:LNBLNK(LABHIS(ITAP))),         OUTPUT.105
     +      STATUS='NEW',FORM='UNFORMATTED')                              OUTPUT.106
       IFIL=0                                                             OUTPUT.107
       RETURN                                                             OUTPUT.108
      ENDIF                                                               OUTPUT.109
      END                                                                 OUTPUT.110
C                                                                         OUTPUT.111
*DK RUNMDS                                                                RUNMDS.1
      BLOCK DATA RUNTIME                                                  RUNMDS.2
C     ****                                                                RUNMDS.3
C     ****     THIS BLOCK DATA CARRIES RAY'S RUN-TIME MODS                RUNMDS.4
C     ****                                                                RUNMDS.5
      INCLUDE 'param.h'                                                   RUNMDS.6
      REAL*8 NOB                                                          RUNMDS.7
      COMMON/RUNT/TTBOUND,NOB(ZJMX),AVTO,HOR(ZJMX)                        RUNMDS.8
      COMMON/MASS/RMN4S,RMN2D,RMNO,BRN2D,CEE                              RUNMDS.9
      DATA RMN4S,RMN2D,RMNO/14.,14.,30./                                  RUNMDS.10
C     ****                                                                RUNMDS.11
C     ****     VALUE OF T AT LOWER BOUNDARY                               RUNMDS.12
C     ****                                                                RUNMDS.13
      DATA TTBOUND/181./                                                  RUNMDS.14
C     ****                                                                RUNMDS.15
C     ****     SET BRN2D FOR CMPN2D                                       RUNMDS.16
C     ****                                                                RUNMDS.17
      DATA BRN2D/0.6/                                                     RUNMDS.18
C     ****                                                                RUNMDS.19
C     ****     N(NO) AT LOWER BOUNDARY, A FUNCTION OF LATITUDE            RUNMDS.20
C     ****       -- 36 VALUES                                             RUNMDS.21
C     ****                                                                RUNMDS.22
      DATA NOB/ZJMX*4.D6/                                                 RUNMDS.23
C     ****                                                                RUNMDS.24
C     ****     SET AVTO                                                   RUNMDS.25
C     ****                                                                RUNMDS.26
      DATA AVTO/1.D-12/                                                   RUNMDS.27
C     ****                                                                RUNMDS.28
C     ****     SET ARRAY GIVING HORIZONTAL VARIATION OF EDDY DIFFUSION    RUNMDS.29
C     ****     AND EDDY THERMAL CONDUCTIVITY.  IF VALUES ARE UNITY,       RUNMDS.30
C     ****     VALUE OF KE AT BOTTOM OF MODEL IS 5.D-6                    RUNMDS.31
C     ****                                                                RUNMDS.32
      DATA HOR/36*.5/                                                     RUNMDS.33
      END                                                                 RUNMDS.34
C                                                                         RUNMDS.35
*DK SDOT                                                                    SDOT.1
      SUBROUTINE SDOT(DOT,N,SX,INCX,SY,INCY)                                SDOT.2
      IMPLICIT REAL*8 (A-H,O-Z)                                             SDOT.3
C     ****     BEGIN PROLOGUE  SDOT                                         SDOT.4
C     ****     DATE WRITTEN   791001   (YYMMDD)                             SDOT.5
C     ****     REVISION DATE  820801   (YYMMDD)                             SDOT.6
C     ****     CATEGORY NO.  D1A4                                           SDOT.7
C     ****     KEYWORDS  BLAS,INNER PRODUCT,LINEAR ALGEBRA,VECTOR           SDOT.8
C     ****     AUTHOR  LAWSON, C. L., (JPL)                                 SDOT.9
C     ****     HANSON, R. J., (SNLA)                                        SDOT.10
C     ****     KINCAID, D. R., (U. OF TEXAS)                                SDOT.11
C     ****     KROGH, F. T., (JPL)                                          SDOT.12
C     ****     PURPOSE  S.P. inner product of s.p. vectors                  SDOT.13
C     ****     DESCRIPTION                                                  SDOT.14
C     ****                                                                  SDOT.15
C     ****           B L A S  Subprogram                                    SDOT.16
C     ****     Description of Parameters                                    SDOT.17
C     ****                                                                  SDOT.18
C     ****     --Input--                                                    SDOT.19
C     ****        N  number of elements in input vector(s)                  SDOT.20
C     ****       SX  single precision vector with N elements                SDOT.21
C     ****     INCX  storage spacing between elements of SX                 SDOT.22
C     ****       SY  single precision vector with N elements                SDOT.23
C     ****     INCY  storage spacing between elements of SY                 SDOT.24
C     ****                                                                  SDOT.25
C     ****     --Output--                                                   SDOT.26
C     ****     SDOT  single precision dot product (zero if N .LE. 0)        SDOT.27
C     ****                                                                  SDOT.28
C     ****     Returns the dot product of single precision SX and SY.       SDOT.29
C     ****     SDOT = sum for I = 0 to N-1 of                               SDOT.30
C     ****     SX(LX+I*INCX) * SY(LY+I*INCY),                               SDOT.31
C     ****     where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N, and LY is  SDOT.32
C     ****     defined in a similar way using INCY.                         SDOT.33
C     ****     ***REFERENCES                                                SDOT.34
C     ****        LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T.,       SDOT.35
C     ****        *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*,     SDOT.36
C     ****         ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL          SDOT.37
C     ****         SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323    SDOT.38
C     ****     ***ROUTINES CALLED  (NONE)                                   SDOT.39
C     ****     ***END PROLOGUE  SDOT                                        SDOT.40
C     ****                                                                  SDOT.41
C     ****      EXTERNAL SX,SY                                              SDOT.42
C     ****     ***FIRST EXECUTABLE STATEMENT  SDOT                          SDOT.43
      DIMENSION SX(N),SY(N)                                                 SDOT.44
      DOT=0.0                                                               SDOT.45
      IF(N.LE.0)RETURN                                                      SDOT.46
      IF(INCX.EQ.INCY) IF(INCX-1)5,20,60                                    SDOT.47
    5 CONTINUE                                                              SDOT.48
C     ****                                                                  SDOT.49
C     ****        CODE FOR UNEQUAL INCREMENTS OR NONPOSITIVE INCREMENTS.    SDOT.50
C     ****                                                                  SDOT.51
      IX=1                                                                  SDOT.52
      IY=1                                                                  SDOT.53
      IF(INCX.LT.0)IX=(-N+1)*INCX+1                                         SDOT.54
      IF(INCY.LT.0)IY=(-N+1)*INCY+1                                         SDOT.55
      DO I=1,N                                                              SDOT.56
       DOT=DOT+SX(IX)*SY(IY)                                                SDOT.57
       IX=IX+INCX                                                           SDOT.58
       IY=IY+INCY                                                           SDOT.59
      ENDDO                                                                 SDOT.60
      RETURN                                                                SDOT.61
C     ****                                                                  SDOT.62
C     ****     CODE FOR BOTH INCREMENTS EQUAL TO 1                          SDOT.63
C     ****                                                                  SDOT.64
C     ****     CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5. SDOT.65
C     ****                                                                  SDOT.66
   20 M=MOD(N,5)                                                            SDOT.67
      IF(M.EQ.0)GO TO 40                                                    SDOT.68
      DO I=1,M                                                              SDOT.69
       DOT=DOT+SX(I)*SY(I)                                                  SDOT.70
      ENDDO                                                                 SDOT.71
      IF(N.LT.5)RETURN                                                      SDOT.72
   40 MP1=M+1                                                               SDOT.73
      DO I=MP1,N,5                                                          SDOT.74
       DOT=DOT+SX(I)*SY(I)+SX(I+1)*SY(I+1)+                                 SDOT.75
     1     SX(I+2)*SY(I+2)+SX(I+3)*SY(I+3)+SX(I+4)*SY(I+4)                  SDOT.76
      ENDDO                                                                 SDOT.77
      RETURN                                                                SDOT.78
C     ****                                                                  SDOT.79
C     ****        CODE FOR POSITIVE EQUAL INCREMENTS .NE.1.                 SDOT.80
C     ****                                                                  SDOT.81
   60 CONTINUE                                                              SDOT.82
      NS=N*INCX                                                             SDOT.83
      DO I=1,NS,INCX                                                        SDOT.84
       DOT=DOT+SX(I)*SY(I)                                                  SDOT.85
      ENDDO                                                                 SDOT.86
      RETURN                                                                SDOT.87
      END                                                                   SDOT.88
C                                                                           SDOT.89
*DK START                                                                  START.1
      SUBROUTINE START                                                     START.2
C     ****                                                                 START.3
C     ****          STARTS THE MODEL FROM HISTORY FILE                     START.4
C     ****                                                                 START.5
      INCLUDE 'param.h'                                                    START.6
      INCLUDE 'buffe.h'                                                    START.7
      INCLUDE 'const.h'                                                    START.8
      INCLUDE 'index.h'                                                    START.9
      INCLUDE 'start.h'                                                    START.10
      LOGICAL FOUND                                                        START.13
      REAL*4 IHEAD(108)                                                    START.14
C     ****                                                                 START.15
C     ****          READ CONTROL PARAMETERS                                START.16
C     ****                                                                 START.17
      CALL INPUT                                                           START.18
C     ****                                                                 START.19
C     ****          DEFINE MODEL CONSTANTS                                 START.20
C     ****                                                                 START.21
      CALL CON                                                             START.22
C     ****                                                                 START.23
C     ****          CALCULATE DYNAMO PARAMETERS                            START.24
C     ****                                                                 START.25
      CALL CNSTNT                                                          START.26
C     ****                                                                 START.27
C     ****          ZERO SCM BUFFER AREA                                   START.28
C     ****                                                                 START.29
      NCBUFK=8*NCBUF+NCPHYS                                                START.30
      DO K=1,NCBUFK                                                        START.31
       DO I=1,IMAXP4                                                       START.32
        F(I,K)=0.0                                                         START.33
       ENDDO                                                               START.34
      ENDDO                                                                START.35
      IF(ISTRT.EQ.1)THEN                                                   START.36
       OPEN(UNIT=IHIST,FILE='ehist',FORM='UNFORMATTED')                    START.37
      ELSE                                                                 START.38
       OPEN(UNIT=IHIST,FILE=LABSRC(1:LNBLNK(LABSRC)),FORM='UNFORMATTED')   START.39
      ENDIF                                                                START.40
      IHIS=MHIS                                                            START.41
      IF(ISTRT.EQ.0)THEN                                                   START.42
       FOUND=.FALSE.                                                       START.43
C      ****                                                                START.44
C      ****    UNBLOCK THE HEADER                                          START.45
C      ****                                                                START.46
1000   CONTINUE                                                            START.47
       READ(IHIST,END=900)(IHEAD(I),I=1,108)                               START.48
C      ****                                                                START.49
C      ****    CONVERT HISTORY FILE HEADER TO COMPUTATION TIME             START.50
C      ****                                                                START.51
       NDAY=IFIX(IHEAD(2))                                                 START.52
       NHOR=IFIX(IHEAD(3))                                                 START.53
       NMIN=IFIX(IHEAD(4))                                                 START.54
C      ****                                                                START.55
C      ****    CHECK FOR DESIRED TIME                                      START.56
C      ****                                                                START.57
       IF(NDAY.EQ.ISOURCE(1).AND.NHOR.EQ.ISOURCE(2).AND.NMIN.              START.58
     1    EQ.ISOURCE(3))FOUND=.TRUE.                                       START.59
       IF(.NOT.FOUND)THEN                                                  START.60
        WRITE(6,"(' Searching for history ',i2,':',i2,':',i2,              START.61
     1    ' Found history ',i2,':',i2,':',i2,' -- still searching')")      START.62
     2    ISOURCE(1),ISOURCE(2),ISOURCE(3),NDAY,NHOR,NMIN                  START.63
c       *****                                                              START.64
c       *****     UNBLOCK THE DATA                                         START.65
c       *****                                                              START.66
        DO J=1,JMAX                                                        START.67
         READ(IHIST,END=900)((F(I,K),I=1,IMAXP4),K=1,KMAXP1*NHIST+1)       START.68
        ENDDO                                                              START.69
        GOTO 1000                                                          START.70
       ELSE                                                                START.71
        WRITE(6,"(' Searching for ut ',i2,':',i2,':',i2,                   START.72
     1    ' Found history ',i2,':',i2,':',i2,' -- Success! ')")            START.73
     2    ISOURCE(1),ISOURCE(2),ISOURCE(3),NDAY,NHOR,NMIN                  START.74
c       *****                                                              START.75
c       *****     READ IN THE DATA                                         START.76
c       *****                                                              START.77
        DO J=1,JMAX                                                        START.78
         READ(IHIST)((F(I,K),I=1,IMAXP4),K=1,KMAXP1*NHIST+1)               START.79
C        ****                                                              START.80
C        ****    STORE F ARRAY INTO G ARRRAY                               START.81
C        ****                                                              START.82
         DO K=1,KMAXP1*NHIST+1                                             START.83
          DO I=1,IMAXP4                                                    START.84
           G(I,K,J)=F(I,K)                                                 START.85
          ENDDO                                                            START.86
         ENDDO                                                             START.87
        ENDDO                                                              START.95
       ENDIF                                                               START.96
      ELSE                                                                 START.97
c      *****                                                               START.98
c      *****     READ IN THE DATA                                          START.99
c      *****                                                               START.100
       DO J=1,JMAX                                                         START.101
        READ(IHIST)((F(I,K),I=1,IMAXP4),K=1,KMAXP1*NHIST+1)                START.102
C       ****                                                               START.103
C       ****    STORE F ARRAY INTO G ARRRAY                                START.104
C       ****                                                               START.105
        DO K=1,KMAXP1*NHIST+1                                              START.106
         DO I=1,IMAXP4                                                     START.107
          G(I,K,J)=F(I,K)                                                  START.108
         ENDDO                                                             START.109
        ENDDO                                                              START.110
       ENDDO                                                               START.118
      ENDIF                                                                START.119
      CLOSE(IHIST)                                                         START.120
C     ****                                                                 START.121
C     ****     CALCULATE TN UN PURTURBATIONS BY TIDES AT LOWER BOUNDARY    START.122
C     ****                                                                 START.123
      CALL BNDRY                                                           START.124
      CALL BNDRY2                                                          START.125
      CALL BNDRYA                                                          START.126
C     ****                                                                 START.127
C     ****     READ IN MAGNETIC COORDINATE TRANSFORMATION METRIX           START.128
C     ****                                                                 START.129
      CALL EINPUT                                                          START.130
      RETURN                                                               START.131
900   WRITE(6,*)'DESIRED TIME OUT OF THE RANGE OF HISTORY FILE, PLEASE     START.132
     1TRY ANOTHER HISTORY FILE'                                            START.133
      STOP                                                                 START.134
      END                                                                  START.135
C                                                                          START.136
*DK TRIVEC                                                                TRIVEC.1
      SUBROUTINE TRIVEC (N,A,B,C,E)                                       TRIVEC.2
      INCLUDE 'param.h'                                                   TRIVEC.3
      DIMENSION A(1),B(1),C(1),E(1)                                       TRIVEC.4
      DATA IINN/2/                                                        TRIVEC.5
      IF(N)120,120,10                                                     TRIVEC.6
   10 R=E(1)                                                              TRIVEC.7
      IF(N-2)20,30,40                                                     TRIVEC.8
   20 E(1)=R/B(1)                                                         TRIVEC.9
      RETURN                                                              TRIVEC.10
   30 QIH=A(IINN)                                                         TRIVEC.11
      BIH=B(IINN)                                                         TRIVEC.12
      GO TO 70                                                            TRIVEC.13
   40 QIH=A(N)                                                            TRIVEC.14
      BIH=B(N)                                                            TRIVEC.15
      DO 60 IDO=3,N                                                       TRIVEC.16
       I=N-IDO+2                                                          TRIVEC.17
       IF(DABS(BIH).LT.DABS(C(I)))GO TO 50                                TRIVEC.18
       RATIO=C(I)/BIH                                                     TRIVEC.19
       C(I)=0.                                                            TRIVEC.20
       B(I+1)=QIH/BIH                                                     TRIVEC.21
C      ****                                                               TRIVEC.22
C      ****     PB(I+1)=0.                                                TRIVEC.23
C      ****     Q(I+1)=QIH/BIH                                            TRIVEC.24
C      ****                                                               TRIVEC.25
       BIH=B(I)-RATIO*QIH                                                 TRIVEC.26
       QIH=A(I)                                                           TRIVEC.27
       GO TO 60                                                           TRIVEC.28
   50  B(I+1)=B(I)/C(I)                                                   TRIVEC.29
       C(I)=A(I)/C(I)                                                     TRIVEC.30
       BIH1=QIH-BIH*B(I+1)                                                TRIVEC.31
       QIH=-BIH*C(I)                                                      TRIVEC.32
       BIH=BIH1                                                           TRIVEC.33
C      ****                                                               TRIVEC.34
C      ****     PB(I+1)=A(I)/C(I)                                         TRIVEC.35
C      ****     Q(I+1)=B(I)/C(I)                                          TRIVEC.36
C      ****                                                               TRIVEC.37
   60 CONTINUE                                                            TRIVEC.38
   70 IF(DABS(BIH).LT.DABS(C(1)))GO TO 80                                 TRIVEC.39
      Q2=QIH/BIH                                                          TRIVEC.40
C     ****                                                                TRIVEC.41
C     ****     Q(2)=QIH/BIH                                               TRIVEC.42
C     ****                                                                TRIVEC.43
      BIH = B(1)-C(1)/BIH*QIH                                             TRIVEC.44
      E(1) = R/BIH                                                        TRIVEC.45
      E(IINN) = -Q2*E(1)                                                  TRIVEC.46
C     ****                                                                TRIVEC.47
C     ****     E(2)=-Q(2)*E(1)                                            TRIVEC.48
C     ****                                                                TRIVEC.49
      GO TO 90                                                            TRIVEC.50
   80 RATIO=BIH/C(1)                                                      TRIVEC.51
      BIH=QIH-RATIO*B(1)                                                  TRIVEC.52
      RIH=-RATIO*R                                                        TRIVEC.53
      E(1)=RIH/BIH                                                        TRIVEC.54
      E(IINN)=(R-B(1)*E(1))/C(1)                                          TRIVEC.55
   90 IF (N-3) 120,100,100                                                TRIVEC.56
  100 DO 110 I=3,N                                                        TRIVEC.57
       E(I)=-B(I)*E(I-1)-C(I-1)*E(I-2)                                    TRIVEC.58
C      ****                                                               TRIVEC.59
C      ****    E(I)=-Q(I)*E(I-1)-PB(I)*E(I-2)                             TRIVEC.60
C      ****                                                               TRIVEC.61
  110 CONTINUE                                                            TRIVEC.62
  120 RETURN                                                              TRIVEC.63
      END                                                                 TRIVEC.64
C                                                                         TRIVEC.65
