      SUBROUTINE CON
C     ****                   DEFINES THE MODEL CONSTANTS
      include "params.h"
      include "blnk.h"
      include "vscr.h"
      include "buff.h"
      include "rfft.h"
      include "index.h"
      include "cons.h"
      include "strt.h"
      include "unit.h"
      include "types.h"
      include "kut.h"
      include "sechis.h"
      DATA ISAVE/3/
C     ****
C     ****     SET PRANDTL NUMBER
C     ****
      DATA PRNDTL/1./
      character*8 nflds_lab(nzflds), nphys_lab(nzphys)
      data nflds_lab
     + /'NT      ','NU      ','NV      ','NPS     ','NPS2    ',
     +  'NPN4S   ','NPNO    ','NOP     ','NPN2D   ','NTI     ',
     +  'NTE     ','NE      ','NO2P    ','NW      ','NZ      ',
     +  'NPHI    ','NPHIH   ','NTNM    ','NUNM    ','NVNM    ',
     +  'NPSNM   ','NPS2NM  ','NN4SNM  ','NPNONM  ','NOPNM   ',
     +  'NMS     ','NVC     ','NKLDU   ','NKLDV   ','NKLDT   ',
     +  'NKLDPS  ','NKLDP2  '/
c
      data nphys_lab
     + /'NRH     ','NFLH    ','NFPH    ','NQDH    ','NKMH    ',
     +  'NPSDH   ','NPSDH2  ','NLXX    ','NLYY    ','NLXY    ',
     +  'NLYX    ','NKT     ','NKM     ','NCP     ','NWTI    ',
     +  'NWTE    ','NNO2    ','NNO     ','NNN2    ','NPDUM   ',
     +  'NNVO2   ','NNVO    ','NNVN2   ','NQO2P   ','NQOP    ',
     +  'NQN2P   ','NQNOP   ','NQNP    ','NQOP2D  ','NQOP2P  ',
     +  'NQTEF   ','NN2P    ','NNP     ','NNOP    ','NQ      ',
     +  'NRJ     ','NUI     ','NVI     ','NWI     '/
c
C     ****                   UNIT NUMBERS
      IHIST = 12
      ISORC = 11
      IUNIN=8
      IUNOT=9
      IMAG=17
!
! bf 5/7/97: changed values of IUSAMI and IUNAMI to avoid conflict
!  with secondary history direct access unit.
!
      IUSAMI = 18
      IUNAMI = 19
C     ****                   IF FORECAST START THEN SET FLAG TO 0
      MODEFC=-1
      IF(MODEST.EQ.3) MODEFC=0
C     ****     INITIALIZE CONTROL PARAMETERS
      IFRST=0
      NERR=0
C     ****                   RESOLUTION PARAMETERS
      IMAX=ZIMX
      JMAX=ZJMX
      KMAX=ZKMX
C     ****                   DERIVED CONSTANTS
      PI=3.14159265358979
      DLAMDA=2.*PI/IMAX
      DPHI=PI/JMAX
      IF(NHEMI.NE.0) DPHI=PI/(2*JMAX)
      DS=(ZST-ZSB)/KMAX
      IMAXH=IMAX/2
      IMAXP2=IMAX+2
      IMAXP4=IMAX+4
      JMAXM1=JMAX-1
      JMAXM2=JMAX-2
      KMAXM1=KMAX-1
      KMAXP1=KMAX+1
      LAST=IMAXP4*KMAX-2
      NPROG=3*KMAXP1
      NPROGP=NPROG+1
C     ****                   BUFFER INFORMATION
C     ****                   COMPUTE NUMBER OF FIELDS ON TAPE
C     ****                   AND INITIALIZE POINTERS
      LEN1=IMAXP4
      LEN2=IMAXP4*KMAX
      LEN3=IMAXP4*KMAXP1
      LONG=LEN2-4
      LBDSK=0
      LBSCM=0
      LBTAP=1
      NCOLS=0
      NCBUF = 0
      DO 15 N=1,NFLDS
      IF(KFLDS(N)-2) 11,12,13
   11 K=1
      LEN=LEN1
      GO TO 14
   12 K=KMAX
      LEN=LEN2
      GO TO 14
   13 K=KMAXP1
      LEN=LEN3
   14 NDEXA(N+1) = NCBUF
      LBSCM=LBSCM+LEN
      NCBUF = NCBUF+K
      IF(N.LE.NDISK) LBDSK=LBDSK+LEN
      IF(N.LE.NTAPE) NCOLS=NCOLS+K
      IF(N.LE.NTAPE) LBTAP=LBTAP+IMAXP2*K
   15 CONTINUE
c
c lbsech = j-slice buffer length for secondary histories
c (dependent on fields selected in secflds input)
c
      lbsech = 1        ! secondary hist buffer length per lat
      if (lusech.gt.0) then
        do n=1,mxfsech
          if (len_trim(secflds(n)).gt.0) lbsech = lbsech + imaxp2*kmaxp1
        enddo
        write(6,"('Buffer length for secondary ',
     +    'histories (lbsech) = ',i6,' (per latitude)')") lbsech
      endif
      LBDSKI = ((LBDSK+511)/512)*512
      LBPHY=0
      NCPHYS=0
      DO 20 N=1,NPHYS
      IF(KPHYS(N)-2) 16,17,18
   16 LEN=LEN1
      K=1
      GO TO 19
   17 LEN=LEN2
      K=KMAX
      GO TO 19
   18 LEN=LEN3
      K=KMAXP1
   19 NDEXB(N+1) = NCPHYS
      NCPHYS = NCPHYS+K
      LBPHY=LBPHY+LEN
   20 CONTINUE
C     ****                   COMPUTE COLUMN POINTERS FOR SCM BUFFERS
      NBADDR(1)=1
      DO 23 I=2,8
   23 NBADDR(I)=NBADDR(I-1)+NCBUF
      NBAPHY=NBADDR(8)+NCBUF
      DO 24 N=1,NPHYS
   24 NDEXB(N+1) = NDEXB(N+1)+NBAPHY
C     ****                   OUTPUT POINTER INFORMATION
      LIM=NFLDS+NPHYS+2
      write(6,"(/72('-'))")
      write(6,"('POINTER VALUES FOR EACH FIELD IN COMMON BLOCK ',
     +    'NDEX:')")
      write(6,"('KMAX=',i3,' KMAXP1=',i3,' ZSB=',f8.3)")
     +  kmax,kmaxp1,zsb
      write(6,"('NDEXA: nflds=',i2)") nflds
      write(6,"('IFLD  NAME      VALUE  KMAX(P1) KPOLE KEQTR')")
      write(6,"('-------------------------------------------')")
      do ifld=1,nflds
        if (kflds(ifld).eq.3) then
          write(6,"(i4,2x,a8,i7,2x,'KMAXP1',2x,2i6)") ifld,
     +      nflds_lab(ifld),ndexa(ifld+1),kpole(ifld),keqtr(ifld)
        elseif (kflds(ifld).eq.2) then
          write(6,"(i4,2x,a8,i7,2x,'KMAX  ',2x,2i6)") ifld,
     +      nflds_lab(ifld),ndexa(ifld+1),kpole(ifld),keqtr(ifld)
        else
          write(6,"(i4,2x,a8,i7,2x,'1     ',2x,2i6)") ifld,
     +      nflds_lab(ifld),ndexa(ifld+1),kpole(ifld),keqtr(ifld)
        endif
      enddo
      write(6,"(/'NDEXB: nphys=',i2)") nphys
      write(6,"('IFLD  NAME      VALUE  KMAX(P1)')")
      write(6,"('-------------------------------')")
      do ifld=1,nphys
        if (kphys(ifld).eq.3) then
          write(6,"(i4,2x,a8,i7,2x,'KMAXP1')") ifld,
     +      nphys_lab(ifld),ndexb(ifld+1)
        elseif (kphys(ifld).eq.2) then
          write(6,"(i4,2x,a8,i7,2x,'KMAX  ')") ifld,
     +      nphys_lab(ifld),ndexb(ifld+1)
        else
          write(6,"(i4,2x,a8,i7,2x,'1     ')") ifld,
     +      nphys_lab(ifld),ndexb(ifld+1)
        endif
      enddo
      write(6,"(72('-')/)")
C     ****                   CHECK IF F DIMENSION IS APPROPRIATE
      MAXFLD=ZFLDX
      MAXPHY=ZPHYX
      WRITE(6,1010)MAXFLD,MAXPHY,NCBUF,NCPHYS
 1010 FORMAT(//* F ARRAY DIMENSIONED FOR FLDX=*I5,2X,*PHYX=*I5,
     A*.*/* MODEL NEEDS NCBUF=*I5,2X,*NCPHYS=*I5,
     B*.*/)
      IF(NCBUF.LE.MAXFLD.AND.NCPHYS.LE.MAXPHY)
     AGO TO 35
      WRITE(6,1020)NCBUF,NCPHYS
 1020 FORMAT(//* EXIT CALLED BY CON*/* NEED FLDX.GE.*I5,2X,*PHYX.GE.
     A*I5)
      STOP
   35 CONTINUE
C     ****                   BUFFER INFORMATION
C     ****                   COMPUTE NUMBER OF FIELDS ON TAPE
C     ****                   AND INITIALIZE POINTERS
      LLEN1=IMAXP4
      LLEN2=IMAXP4*(KMAX-4)
      LLEN3=IMAXP4*(KMAXP1-4)
      LLONG=LLEN2-4
      LLBDSK=0
      LLBSCM=0
      LLBTAP=1
      NNCOLS=0
      NNCBUF=0
      DO N=1,NFLDS
        IF (KFLDS(N).EQ.1) THEN
          K=1
          LEN=LLEN1
        ELSE IF (KFLDS(N).EQ.2) THEN
          K=KMAX-4
          LEN=LLEN2
        ELSE IF (KFLDS(N).EQ.3) THEN
          K=KMAXP1-4
          LEN=LLEN3
        ENDIF
        NNDEXA(N+1) = NNCBUF
        LLBSCM=LLBSCM+LEN
        NNCBUF = NNCBUF+K
        IF(N.LE.NDISK) LLBDSK=LLBDSK+LEN
        IF(N.LE.NTAPE) NNCOLS=NNCOLS+K
        IF(N.LE.NTAPE) LLBTAP=LLBTAP+IMAXP2*K
      ENDDO
      LLBDSKI = ((LLBDSK+511)/512)*512
      LLBPHY=0
      NNCPHYS=0
      DO N=1,NPHYS
        IF(KPHYS(N).EQ.1) THEN
          LEN=LLEN1
          K=1
        ELSE IF (KPHYS(N).EQ.2) THEN
          LEN=LLEN2
          K=KMAX-4
        ELSE IF (KPHYS(N).EQ.3) THEN
          LEN=LLEN3
          K=KMAXP1-4
        ENDIF
        NNDEXB(N+1) = NNCPHYS
        NNCPHYS = NNCPHYS+K
        LLBPHY=LLBPHY+LEN
      ENDDO
C     ****                   INITIALIZE FFT ROUTINE
      CALL SET99(TRIGS,IFAX,IMAX)
C     ****                   COMPUTE PHI VARIATION ARRAYS
      DELT=86400./FLOAT(IMAX)
      JMAXH=JMAX/2
      OMEGA=7.292E-5
      JS=-JMAXH
      IF(NHEMI.NE.0) JS=0
      DO 40 J=1, JMAX
      PHI=(J+JS-.5)*DPHI
      CS(J)=COS(PHI)
      SN(J)=SIN(PHI)
      TN(J)=TAN(PHI)
      COR(J)=2.*OMEGA*SN(J)
   40 CONTINUE
      CSSP(1)=-CS(2)
      CSSP(2)=-CS(1)
C     ****                   IF HEMI CASE, COS HAS SAME SIGN ACROSS EQUA
      IF(NHEMI.NE.0) CSSP(1)=CS(2)
      IF(NHEMI.NE.0) CSSP(2)=CS(1)
      CSNP(1)=-CS(JMAX)
      CSNP(2)=-CS(JMAXM1)
      KUTSP(1) = KUT(2)
      KUTSP(2) = KUT(1)
      KUTNP(1) = KUT(JMAX)
      KUTNP(2) = KUT(JMAXM1)
C     ****     CALCULATE SUN'S DECLINATION
C     ****     (WILL HAVE TO RECALCULATE THIS IF DAY IS ADVANCED)
      DELTA=ATAN(TAN(23.5*PI/180.)*SIN(2.*PI*FLOAT(IIDAY-80)/365.))
C     ****                   GRID CONSTANTS
      C(1)=DLAMDA
      C(2)=DPHI
      C(3)=DS
      C(6)=2.*C(4)
      C(7)=1./C(6)
C     ****                   FINITE DIFFERENCE OPERATOR CONSTANTS
      C(10)=2./(3.*DLAMDA)
      C(11)=1./(12.*DLAMDA)
      C(12)=2./(3.*DPHI)
      C(13)=1./(12.*DPHI)
      C(14)=1./(2.*DLAMDA)
      C(15)=1./(2.*DPHI)
C     ****                   DIFFUSION COEFFICIENT
      C(16)=0.
      C(16)=C(16)/6.
      C(17)=0.
C     ****     INTEGRATION CONSTANTS
      C(18)=C(3)/(1.+.5*C(3))
      C(19)=(1.-.5*C(3))/(1.+.5*C(3))
      C(20)=1./DS
      C(22)=COS(PI/3.)
      C(23)=4.*PI/(24.*60.*60.)
      C(24)=28.9
C     ****     ANNUAL FREQUENCY
      C(25) = C(23)/(2.*365.25)
C     ****     SHAPIRO SMOOTHER CONSTANT
      C(26)=3.0E-2
C     ****     EDDY DIFFUSION
C     C(27)=2.5E-6
      C(27)=5.0E-6
C     ****     EDDY THERMAL CONDUCTION
      C(28)=C(27)/PRNDTL
C     ****                   TIME SMOOTHING CONSTANT
      ALPHA=.95
      C(30)=ALPHA
      C(31)=.5*(1.-ALPHA)
C     ****     RATE CONSTANT FOR ATOMIC OXYGEN RECOMBINATION
      C(44)=0.
C     *****    SURPLUS HEAT PER EVENT CONVERTED FROM E.V. TO ERGS
      C(45)=5.11*1.602E-12
C     ****                   PHYSICAL CONSTANTS
C     ****
C                                 EARTHS RADIUS
      C(51)=6.37122E8
      C(52)=1./C(51)
      C(53)=C(51)*C(51)
C                                 GRAVITATIONAL CONSTANT
      C(54)=870.
      C(55)=1./C(54)
C                                 EARTHS ANGULAR VELOCITY
      C(56)=OMEGA
C                                 GAS CONSTANT
      C(57)=8.314E7
C     ****     H*C
      C(60)=1.9845E-16
C     ****     F107, F107A
C     C(61) = F107
C     C(62) = F107A
C     ****                        G/R
      C(65)=C(54)/C(57)
C     ****     STANDARD PRESSURE
      C(81)=5.0E-4
C     ****     EDDY VISCOSITY
      C(82)=0.
      C(83)=C(81)*C(82)/C(57)*EXP(C(3)/2.)
C     ****     BOLTZMANN CONST
      C(84)=1.38E-16
C     ****     AVOGRADO NUMBER
      C(85)=6.023E23
C     ****     EXP(-.5*DS)
      C(86)=EXP(-.5*DS)
C     ****     EXP(.5*DS)
      C(87)=1./C(86)
C     ****     SUN'S DECLINATION:
C     ****     (WILL HAVE TO RECALCULATE IF DAY IS ADVANCED)
      C(95)=SIN(DELTA)
      C(96)=COS(DELTA)
C     ****     HEIGHT OF BOTTOM OF MODEL
      C(97)=97.E5
C     ****    O+ EDDY DIFFUSION
      C(100)=0.
C     ****                   KO  COEFFICIENT OF EDDY VISCOSITY
      C(109)=0.2
      C(110) = PI
C     ****         SAVE ORIGINAL DT VALUE
      C(111) = C(4)
      C(112)=SQRT(PI)
      JTOTS=0.
      JTOTN=0.
C     ****                   DUMP CONSTANTS
      WRITE(6,70) IMAX,JMAX,KMAX,NFLDS,NPHYS,C(4),DS
   70 FORMAT(*1IMAX=*I3,2X,*JMAX=*I2,2X,*KMAX=*I2,2X,*NFLDS=*I3,2X,
     A*NPHYS=*I3,2X,*DT=*E10.3,2X,*DS=*E10.3)
      WRITE(6,80) LBDSK,LBSCM,LBPHY,LBTAP
   80 FORMAT(* LBDSK=*I6,2X,*LBSCM=*I6,2X,*LBPHY=*I6,2X,*LBTAP=*I6)
      WRITE(6,95) KUT
   95 FORMAT(1X,*KUT = *20I5/(7X,20I5))
      WRITE(6,90) CSSP
   90 FORMAT(*0 J*,4X,*CS(J)*,7X,*F (J)*7X,*SN(J)*,7X,*TN(J)*//
     A(3X,E12.3))
      WRITE(6,100) (J,CS(J),COR(J),SN(J),TN(J),J=1,JMAX)
  100 FORMAT(1H ,I2,4E12.3)
      WRITE(6,105) CSNP
  105 FORMAT(3X,E12.3)
      WRITE(6,130) C
  130 FORMAT(*1C=*10E12.3/(3X,10E12.3))
C     ****     CALCULATE EXPS AND A ARRAYS
      S=ZSB+.5*DS
      EXPS(1)=EXP(-S)
      EXPDS=EXP(-DS)
      A(1)=0.
      DIFK(1)=C(27)
      DIFT(1)=C(28)
      XMUE(1) = C(27)
      DO 1 K=2,KMAX
      EXPS(K)=EXPS(K-1)*EXPDS
      A(K)=A(K-1)*EXPDS
      DIFK(K)=DIFK(K-1)*EXPDS
      DIFT(K) = DIFT(K-1)*EXPDS
      XMUE(K) = DIFK(K)
    1 CONTINUE
      DIFK(KMAXP1)=DIFK(KMAX)*EXPDS
      DIFT(KMAXP1) = DIFT(KMAX)*EXPDS
      XMUE(KMAXP1) = DIFK(KMAXP1)
C     ****
      WRITE(6,135)(K,EXPS(K),A(K),DIFK(K),DIFT(K),K=1,KMAX)
  135 FORMAT(*0 K*,5X,*EXPS(K)*,8X,*A(K)*,5X,*DIFK(K)*,5X,*DIFT(K)*//
     A(1H ,I2,4E12.3))
      RETURN
      END
C
