      SUBROUTINE ADVNCE
C     ****                   ADVANCES THE MODEL IN TIME
      include "params.h"
      PARAMETER (IMAXG=ZIMX, JMAXG=ZJMX)
      PARAMETER (IMAXGP=IMAXG+1, JMAXGP=JMAXG+1)
      include "amie.h"
      include "blnk.h"
      include "buff.h"
      include "cons.h"
      include "consts.h"
      include "cterp.h"
      include "dynphi.h"
      include "index.h"
      include "phys.h"
      include "sechis.h"
      include "strt.h"
      include "unit.h"
      include "vscr.h"
      common /sight/ sigp(73,0:37),sigh(73,0:37),qjn(73,0:37),
     | tec(73,0:37),qwind(73,0:37),qamie(73,0:37),wtot(73,0:37),
     | work(73,0:37),fwindu(73,0:37),fwindv(73,0:37),
     | famieu(73,0:37),famiev(73,0:37),ftot(73,0:37),
     | famie(73,0:37),fwind(73,0:37)
      COMMON/TRGM/RLATM(ZIMXP,ZJMX),RLONM(ZIMXP,ZJMX),
     1  DIPMAG(ZIMXP,ZJMX),DECMAG(ZIMXP,ZJMX),SNDEC(ZIMXP,ZJMX),
     2  CSDEC(ZIMXP,ZJMX),SN2DEC(ZIMXP,ZJMX),SNCSDC(ZIMXP,ZJMX),
     3  DUMDUM(ZJMX,3)
      dimension sigpm(imaxmp,jmaxm),sighm(imaxmp,jmaxm)
      common /facvne/face(imaxmp,jmaxmh),facvn(imaxmp,jmaxmh),
     |  fact(imaxmp,jmaxmh),fac3(imaxmp,jmaxmh,3)
      dimension factm(imaxmp,jmaxm),factg(imaxgp,0:jmaxgp)
      character*80 uthdr
      character*8 blank
      data blank/'        '/
      dimension gcmlon(73),gcmlat(0:37),st(zjmx)
      integer mtime(3)

C  ekv0 = average mean energy from POLAR each MLT=0,28
      dimension ekv0(29)
      data ekv0/1.61,1.59,1.63,1.58,1.68,1.84,2.04,2.07,1.89,1.62,
     |          1.23,1.13,0.97,0.97,0.80,0.64,0.59,0.64,0.74,0.80,
     |          0.88,1.08,1.21,1.41,1.61,1.59,1.63,1.58,1.68/

      DIMENSION DTPOT1(80,49),DTPOT2(80,49),
     | DPKVA1(73,36),DVIA1(73,36,3),DKEVA1(73,36),DFLXA1(73,36),
     | DPKVA2(73,36),DVIA2(73,36,3),DKEVA2(73,36),DFLXA2(73,36)
C     ****                   INITIALIZE BUFFER INDICES
      MODULO=8
      CALL INIT
      jfirst = 0
C  Set up gcmlon and gcmlat to include wrap-around and poles
      do i=1,73
	gcmlon(i) = -185. + i*5.
      enddo
      do i=1,36
	gcmlat(i) = -92.5 + i*5.
      enddo
	gcmlat(0) = -90.
	gcmlat(37) = 90.

C     ****                   TIME LOOP
      TIME = SECOND(X)
      TIMET = TIMEF(X)
  100 CONTINUE
      jfirst = jfirst + 1
      ITER=ITER+1
      write(6,"('ADVANCE: iter = ',i6)") iter
      JTOTS=0.
      JTOTN=0.
C     ****
C     ****     FIND UT TIME BETWEEN 0-24 HR (MODEL DAY = UT DAY, 11/89)
C     ****
      SECS = AMOD(FLOAT(ITER)*C(4),86400.)
      IF (IADVDA .EQ. 1) THEN
       IDAYIT = ITER*IFIX(C(4))/86400
       IF (IDAYIT*86400 .EQ. ITER*IFIX(C(4))) THEN
C       ****
C       ****     ADVANCE DAY
C       ****
        IDAPR = IIDAY
        IYRPR = IYEAR
        IIDAY = IIDAY + 1
C       ****
C       ****     LPYR = 1(0)  IF IS (NOT) A LEAP YEAR
C       ****
        IYR4 = IYEAR/4
        IYR100 = IYEAR/100
        LPYR = 0
        IF (IYR4*4 .EQ. IYEAR .AND. IYR100*100 .NE. IYEAR) LPYR=1
        IENDA = 365 + LPYR
        IF (IIDAY .GT. IENDA) THEN
         IYEAR = IYEAR + 1
         IIDAY = IIDAY - IENDA
        ENDIF                          !   FOR PAST YEAR'S END
C       ****
C       ****     RECALCULATE SUN'S DECLINATION
C       ****
        DELTA = ATAN(TAN(23.5*PI/180.)*SIN(2.*PI*FLOAT(IIDAY-80)/365.))
        C(95) = SIN(DELTA)
        C(96) = COS(DELTA)
        WRITE(6,"(1X,'ADVNCE:  ADVANCING DAY (PREVIOUS,PRESENT)=',4I5)")
     1   IDAPR,IYRPR,IIDAY,IYEAR
       ENDIF                           !   FOR ADVANCING DAY
      ENDIF                           !   FOR IADVDA FLAG ON
C     ****
C     ****     CALCULATE AURORAL PARAMETERS FOR THIS TIME STEP
C     ****
C     ****     FIND UT TIME  (MODEL DAY = UT DAY, 11/89)
      SEC = FLOAT(ITER)*C(4)
      SECTGCM = FLOAT(ITER)*C(4)
      CALL SUN

C  Use AMIE calculation
      IF (IAMIE.EQ.1) THEN
       DO 125 IH=1,2
       IF (NTST(IH).LT.NTIMS(IH) .AND. SEC.GT.UTSP(2,IH)) THEN
	NTST(IH) = NTST(IH) + 1
	UTSP(1,IH) = UTSP(2,IH)
	UTSA(1,IH) = UTSA(2,IH)
	UTSASN(1,IH) = UTSASN(2,IH)
	HPASN(1,IH) = HPASN(2,IH)
	HPA(1,IH) = HPA(2,IH)
	CPA(1,IH) = CPA(2,IH)
	CUSPTA(1,IH) = CUSPTA(2,IH)
	CUSPLA(1,IH) = CUSPLA(2,IH)
	IUAMI = IUSAMI
	J1 = 1
	J2 = 18
	IF (IH .EQ. 2) THEN
	 J1 = 19
	 J2 = 36
	 IUAMI = IUNAMI
	ENDIF
	DO 105 JA=J1,J2
	  DO 105 IA=1,73
	    POTKVA1(IA,JA) = POTKVA2(IA,JA)
	    EKEVA1(IA,JA) = EKEVA2(IA,JA)
	    svEKEVA1(IA,JA) = EKEVA2(IA,JA)
	    EFLXA1(IA,JA) = EFLXA2(IA,JA)
	    DO 105 KA=1,3
	      VIXYZA1(IA,JA,KA) = VIXYZA2(IA,JA,KA)
  105   CONTINUE
	DO 106 JA=1,36
	  DO 106 IA=1,73
	    EKEVASN(IA,JA,1,IH) = EKEVASN(IA,JA,2,IH)
	    EFLXASN(IA,JA,1,IH) = EFLXASN(IA,JA,2,IH)
  106   CONTINUE
	J2 = 47
	IF (IH .EQ. 2) THEN
	 J1 = 48
	 J2 = 97
	ENDIF
	DO 115 JA=J1,J2
	  DO 115 IA=1,81
	    TIEPOT1(IA,JA) = TIEPOT2(IA,JA)
  115   CONTINUE
C       READ (IUAMI) UTSA(2,IH),HPA(2,IH),CPA(2,IH),CUSPTA(2,IH),
	READ (IUAMI) IUTSA,     HPA(2,IH),CPA(2,IH),CUSPTA(2,IH),
     1    CUSPLA(2,IH),      DPKVA2,DTPOT2,DVIA2,DKEVA2,DFLXA2
C    1    CUSPLA(2,IH),ICONJ,DPKVA2,DTPOT2,DVIA2,DKEVA2,DFLXA2
	UTSA(2,IH) = IUTSA
C  Lummerzheim gave me characteristic energy, not mean energy!!
	fackv = 1.
	WRITE (6,"(1x,'ADVNC read AMIE: IH UTSA HP CP CMLT CLAT =',
     |   I2,F9.0,2F6.1,2F6.2)") IH,UTSA(2,IH),HPA(2,IH),CPA(2,IH),
     |   CUSPTA(2,IH),CUSPLA(2,IH)
	UTSP(2,IH) = UTSA(2,IH)
	UTSASN(2,IH) = UTSA(2,IH)
	HPASN(2,IH) = HPA(2,IH)
C       ****      Set dummy1 in array1, dummy2 in array2 for hemisphere
	DO 116 JA=1,36
	  DO 116 IA=1,73
	    EKEVASN(IA,JA,2,IH) = DKEVA2(IA,JA) * fackv
	    EFLXASN(IA,JA,2,IH) = DFLXA2(IA,JA)
  116   CONTINUE
	JA = 0
	IF (IH .EQ. 2) JA = 18
	DO 62 I=1,73
	DO 62 J=1,18
	POTKVA2(I,J+JA) = DPKVA2(I,J+JA)
	EKEVA2(I,J+JA) = DKEVA2(I,J+JA) * fackv
	svEKEVA2(I,J+JA) = DKEVA2(I,J+JA) * fackv
	EFLXA2(I,J+JA) = DFLXA2(I,J+JA)
	DO 62 K=1,3
   62   VIXYZA2(I,J+JA,K) = DVIA2(I,J+JA,K)
	IF (IH .EQ. 1) THEN
	 DO 64 J=1,49
	 DO 63 I=1,80
   63    TIEPOT2(I,J) = DTPOT2(I,J)
   64    TIEPOT2(81,J) = TIEPOT2(1,J)
	ELSE
	 DO 66 J=1,49
	 DO 65 I=1,80
   65    TIEPOT2(I,97+1-J) = DTPOT2(I,J)
   66    TIEPOT2(81,97+1-J) = TIEPOT2(1,97+1-J)
	ENDIF
	j=1
c       write (6,"(1x,'j=1 kVfxeV',6e12.4)") potkva1(1,j),potkva1(73,j),
c    |   eflxa1(1,j),eflxa1(73,1),ekeva1(1,j),ekeva1(73,j)
c       write (6,"(1x,10e12.4)") potkva2(1,j),potkva2(73,j),
c    |   eflxa2(1,j),eflxa2(73,1),ekeva2(1,j),ekeva2(73,j)
	if (ih .eq. 1) then
c       CALL EZCNTR (tiepot2,81,20)
c       CALL EZCNTR (eflxa2,73,10)
c       CALL EZCNTR (ekeva2,73,10)
	endif
       ENDIF                    !  FOR READING AMIE
  125  CONTINUE
C       ****      If ICONJ=1, find times closest to SECTGCM and
C       ****      substitute values in opposite hemisphere for FLX,KEV
	if (iconj .eq. 1) call fndutsa
C TEMP for May 1996 - since Ne neg in elden.f at j=3, k=1, i=64 where there
C   is a blob of eflux and low kev, multiply keV by factor 2 or less
c        uts = sec - 1641600.
c        uth = uts/3600.
c        nuth = nint(uth) + 1
c        fac = (uth-16.)*2.5/1.5
c        write (6,"(1x,'uth facm fac =',3f8.2)") uth,facm,fac
	 DO 78 I=1,73
c        lt = i/3 + nuth + 12
c        if (lt .gt. 24) lt = lt - 24
c        ekv1 = amax1(ekv0(lt),3.-fac)
	 DO 78 ja=1,36
c        ekeva1(I,Ja) = svekeva1(I,Ja)*fac
c        ekeva2(I,Ja) = svekeva2(I,Ja)*fac
	 ekeva1(I,Ja) = 3.
	 ekeva2(I,Ja) = 3.
c        ekeva1(I,Ja) = 4.
c        ekeva2(I,Ja) = 4.
c        ekeva1(I,Ja) = ekv1
c        ekeva2(I,Ja) = ekv1
c        ekeva1(I,Ja) = amax1(3.,svekeva1(I,Ja))
c        ekeva2(I,Ja) = amax1(3.,svekeva2(I,Ja))
   78    continue
c       if (sec .lt. 1642000) CALL EZCNTR (ekeva1,73,10)
c       if (sec .gt. 1702440) CALL EZCNTR (ekeva1,73,10)
c       CALL EZCNTR (ekeva1,73,10)
c       CALL EZCNTR (ekeva2,73,10)
      ENDIF                     !  FOR IAMIE=1
      call getgpi
      CALL TAIL (SEC)
      IHIS=IHIS-1
      ISAV=ISAV-1
      if (iter.le.isecstop) then
        ihissech = ihissech-1
        isavsech = isavsech-1
      endif
C     ****                   FOR THE FIRST FORWARD TIME STEP OF A
C     ****                   FORECAST START, HALVE DT.
      IF(MODEFC)150,130,140
  130 C(4) = C(111)/2.
      C(6)   = 2.*C(4)
      C(7)   = 1./C(6)
      MODEFC=1
      GO TO 150
C     ****                   RESTORE DT TO ORIGINAL VALUE
  140 C(4) = C(111)
      C(6)   = 2.*C(4)
      C(7)   = 1./C(6)
      MODEFC=-1
  150 CONTINUE
C     ****                   BEGIN SWEEP OF MESH FROM SOUTH TO NORTH
C     ****                   BRING IN LINE 1
      CALL VDRIFT
      CALL FLIP
      CALL UNLCRD(IUNIN,NJIN,LBDSKI,NSIN)
      CALL UNLCCK(IUNIN,NSIN,LEN)
      IF(LEN.EQ.LBDSKI.AND.NSIN.EQ.0) GO TO 155
      NUM=155
      GO TO 4000
C     ****                   BRING IN LINES 2 AND 3
C     ****                   AND REARRANGE LINES 1 AND 2
C     ****                   TO FORM LINES 0 AND -1
  155 DO 270 JP2=1,2
      CALL FLIP
      CALL UNLCRD(IUNIN,NJIN,LBDSKI,NSIN)
C     ****                   ADD DIAGNOSTIC QUANTITIES ONTO LINE JP2
      J=JP2-2
      CALL ADDIAG
C     ****                   REARRANGE LINE JP2 TO GET LINE -JP2+1
      NJPUT=NJP1
      IF(JP2.EQ.2) NJPUT=NJM1
      IF(NHEMI.NE.0) GO TO 200
C     ****                   CONSTRUCT LINES -1,0 FOR GLOBAL CASE
      ILIM=IMAXH+2
      DO 195 N=1,NFLDS
C     ****                   CHECK IF FIELD IS NEEDED ACROSS POLE
      IF(KPOLE(N)) 160,195,160
C     ****                   TRANSFER Z TYPES FIRST ITERATION ONLY
  160 IF(IFRST.EQ.1 .AND. KPOLE(N).EQ.2) GO TO 195
C     ****                   DETERMINE LENGTH OF FIELD
      IF(KFLDS(N)-2) 165,167,170
  165 KLIM=1
      LEN=LEN1
      GO TO 175
  167 KLIM=KMAX
      LEN=LEN2
      GO TO 175
  170 KLIM=KMAXP1
      LEN=LEN3
C     ****                   SET POINTERS TO INPUT AND OUTPUT BUFFERS
  175 NXK   = NJP2+NDEXA(N+1)
      NYK   = NJPUT+NDEXA(N+1)
C     ****                   CHECK IF FIELD IS TO HAVE SIGN CHANGED
      IF(KPOLE(N)) 180,195,185
C     ****                   CHANGE SIGN
  180 DO 184 K=1,KLIM
      ID=IMAXH+1
      DO 182 I=1,ILIM
      F(ID+2,NYK)=-F(I+2,NXK)
      F(I   ,NYK)=-F(ID ,NXK)
      ID=ID+1
  182 CONTINUE
      NXK=NXK+1
      NYK=NYK+1
  184 CONTINUE
      GO TO 195
C     ****                   DO NOT CHANGE SIGN
  185 DO 188 K=1,KLIM
      ID=IMAXH+1
      DO 187 I=1,ILIM
      F(ID+2,NYK)=F(I+2,NXK)
      F(I   ,NYK)=F(ID ,NXK)
      ID=ID+1
  187 CONTINUE
      NXK=NXK+1
      NYK=NYK+1
  188 CONTINUE
  195 CONTINUE
      GO TO 267
C     ****                   CONSTRUCT LINES -1,0 FOR HEMISPHERIC CASE
  200 DO 265 N=1,NFLDS
C     ****                   CHECK IF FIELD IS NEEDED ACROSS EQUATOR
      IF(KEQTR(N)) 205,265,205
C     ****                   TRANSFER Z TYPES FIRST ITERATION ONLY
  205 IF(IFRST.EQ.1 .AND. KEQTR(N).EQ.2) GO TO 265
C     ****                   DETERMINE LENGTH OF FIELD
      IF(KFLDS(N)-2) 210,215,220
  210 KLIM=1
      LEN=LEN1
      GO TO 225
  215 KLIM=KMAX
      LEN=LEN2
      GO TO 225
  220 KLIM=KMAXP1
      LEN=LEN3
C     ****                   SET POINTERS TO FIELDS IN SCM BUFFERS
  225 NJP2K = NJP2+NDEXA(N+1)
      NJPUTK= NJPUT+NDEXA(N+1)
C     ****                   CHECK IF FIELD IS TO HAVE SIGN CHANGED
      IF(KEQTR(N)) 230,265,260
C     ****                   CHANGE SIGN AND STORE IN PROPER SCM BUFFER
  230 DO 235 I=1,LEN
  235 F(I,NJPUTK)=-1.*F(I,NJP2K)
      GO TO 265
C     ****                   DO NOT CHANGE SIGN AND STORE IN PROPER
C     ****                   SCM BUFFER.
  260 DO 263 I=1,LEN
  263 F(I,NJPUTK)=F(I,NJP2K)
  265 CONTINUE
  267 CONTINUE
C     ****                   CHECK FOR COMPLETION OF BUFFERING
      CALL UNLCCK(IUNIN,NSIN,LEN)
      IF(LEN.EQ.LBDSKI.AND.NSIN.EQ.0) GO TO 270
      NUM=270
      GO TO 4000
  270 CONTINUE
C     ****                   CHANGE BUFFER POINTERS SO THAT LINE J IS -2
C     ****                   SO THAT INITIALIZATION CALLS TO HDIF
C     ****                   CAN BE MADE.
      DO 280 J=1,5
  280 CALL FLIP
C     ****                   ITITIALIZE HORIZONTAL DIFFUSION QUANTITIES
      DO 290 JJ=1,3
      CALL FLIP
      J=JJ-3
  290 CALL HDIF
C     ****                   GENERAL LINE LOOP (LATITUDINAL VARIATION)
      DO 3000 J=1,JMAX

       if (iter==6340) write(6,"('advnce enter j loop: j=',i2)") j

      RACS=1./(C(51)*CS(J))
      CALL FLIP
C     ****                   START IN NEXT LINE NORTH
      IF(J.GE.JMAXM2) GO TO 310
      CALL UNLCRD(IUNIN,NJIN,LBDSKI,NSIN)
C     ****                   START OUT PREVIOUS LINE SOUTH
  310 IF(J.EQ.1) GO TO 315
      CALL UNLCWT(IUNOT,NJOT,LBDSKI,NSOT)
C     ****                   ADD DIAGNOSTIC QUANTITIES ONTO LINE J+2
  315 CONTINUE
      IF(J.GT.JMAXM2)GO TO 405
      CALL ADDIAG
  405 CONTINUE
      IF(J.EQ.JMAX) IFRST=1
C     ******************************************************************
C     ****
C     ****                   DYNAMICS
C     ****
C     ******************************************************************
C     ****
C     ****    CALL TO BNDCMP SHOULD LEAD CALLS IN DYNAMICS SECTION OF
C     ****      ADVNCE TO ENSURE THAT MATRICES B(ZIMXP,2,2) AND VECTORS
C     ****      FB(ZIMXP,2) ARE AVAILABLE WHEN NEEDED
C     ****
      CALL BNDCMP
C     ****
C     ****    CALCULATE W AT N*DT
C     ****
      CALL SWDOT
C     ****
C     ****     CALCULATE ION DRIFT VELOCITY
C     ****
      CALL VDRIFT2
C     ****     CALCULATE CP, KT, AND KM
      CALL CPKTKM
C     ****      CLEAR NQO2P, NQOP, NQN2P, NQNOP, NQNP
      NQPK=NQO2P
      DO 1 I=1,5*LEN3
      F(I,NQPK)=0.
    1 CONTINUE
C     ****     CALCULATE COLUMN NUMBER DENSITIES
      CALL CHAPMN
C     ****     CALCULATE SOLAR HEATING
C     ****     CALCULATE REACTION RATES
      CALL RATES
      CALL QRJ
C     ****     CALCULATE BACKGROUND (NIGHT-TIME) IONIZATION
      CALL QINITE
C     ****     ADD IONIZATION DUE TO SOLAR X-RAYS
      CALL XRAY
C     ****     GENERATE AURORAL FIELDS, UI, VI, WI.  FIELDS RETURNED
C     ****     IN T1 THRU T7
      CALL HEELIS
C     ****     CALCULATE AURORAL ADDITIONS TO IONIZATION RATES
      CALL ORORA
C     ****     SOLVE FOR N(O+)
C     ****     SET FLUX AT UPPER BOUNDARY
      CALL OPFLUX
      CALL OPLUS
C     ****     COMPUTE ELECTRON DENSITY, N(E), AND ION NUMBER DENSITIES,
C     ****     N(O2+), N(N2+), N(NO+), N(N+)
      CALL ELDEN

      if (iter==6340) write(6,"('advnce before lamdas: j=',i2)") j

C     ****     ION-DRAG PARAMETERS, LXX, LYY, LXY (IN DIPOLE COORDINATE
C     ****     SYSTEM)
      CALL LAMDAS

      if (iter==6340) write(6,"('advnce after lamdas: j=',i2)") j

C     ****     ROTATE LAMDAS TO GEOGRAPHIC SYSTEM
      CALL ROTATE
C     ****     ELECTRON PRODUCTION OF N4S AND N2D
      CALL QTIEFF
C     ****     ODD NITROGEN CHEMISTRY
      CALL CMPN2D
      CALL CMPN4S
      CALL CMPNO
      CALL CMPO2O
C     ****     HEATING AND O2 DISSOCIATION DUE TO ODD NITROGEN AND ION
C     ****     CHEMISTRY
      CALL QJNNO
      CALL QJION
C     ****     CALCULATE ELECTRON AND ION TEMPERATURES
      CALL SETTEI
C     ****     CALCULATE IMPLICIT AND EXPLICIT COOLING TERMS IN NWTI
C     ****     AND NWTE
      CALL NEWTON
C     ****     ADD IN O3P COOLING
      CALL NEWTO3P
C     ****     CALCULATE FLF, FPH, NQDH, NKMH, NPSDH, NPSDH2
C     ****     AT (N-1)*DT
      CALL HDIF
C     ****     ADVANCE T, U, V BY ONE TIME STEP
      CALL DT
      CALL DUV
C     ****     COMPOSITION OF MAJOR SPECIES
      CALL COMP
C     ******************************************************************
C     ****
C     ****                   END OF DYNAMICS
C     ****
C     ******************************************************************
C     ****                   CHECK IF ON LINE JMAXM2 OR JMAXM1
      IF(J.LT.JMAXM2) GO TO 2980
      IF(J.EQ.JMAX) GO TO 2990
C     ****                   REARRANGE LINE NJGET TO GET LINE NJIN
      NJGET=NJP2
      IF(J.EQ.JMAXM1) NJGET=NJ
      ILIM=IMAXH+2
      DO 2975 N=1,NFLDS
C     ****                   CHECK IF FIELD IS NEEDED ACROSS POLE
      IF(KPOLE(N)) 2916,2975,2916
C     ****                   TRANSFER Z TYPES FIRST ITERATION ONLY
 2916 IF(IFRST.EQ.1.AND.KPOLE(N).EQ.2) GO TO 2975
C     ****                   DETERMINE LENGTH OF FIELD
      IF(KFLDS(N)-2) 2920,2925,2930
 2920 KLIM=1
      LEN=LEN1
      GO TO 2935
 2925 KLIM=KMAX
      LEN=LEN2
      GO TO 2935
 2930 KLIM=KMAXP1
      LEN=LEN3
C     ****                   BRING IN FIELD FROM SCM BUFFER
 2935 NXK   = NJGET+NDEXA(N+1)
      NYK   = NJIN+NDEXA(N+1)
C     ****                   CHECK IF FIELD IS TO HAVE SIGN CHANGED
      IF(KPOLE(N)) 2940,2975,2955
C     ****                   CHANGE SIGN
 2940 DO 2950 K=1,KLIM
      ID=IMAXH+1
      DO 2945 I=1,ILIM
      F(ID+2,NYK)=-F(I+2,NXK)
      F(I   ,NYK)=-F(ID ,NXK)
      ID=ID+1
 2945 CONTINUE
      NXK=NXK+1
      NYK=NYK+1
 2950 CONTINUE
      GO TO 2975
C     ****                   DO NOT CHANGE SIGN
 2955 DO 2965 K=1,KLIM
      ID=IMAXH+1
      DO 2960 I=1,ILIM
      F(ID+2,NYK)=F(I+2,NXK)
      F(I   ,NYK)=F(ID ,NXK)
      ID=ID+1
 2960 CONTINUE
      NXK=NXK+1
      NYK=NYK+1
 2965 CONTINUE
C     ****                   TRANSFER REARRANGED FIELD
C     ****                   TO CORRECT LOCATION IN SCM BUFFER
 2975 CONTINUE
      GO TO 2990
C     ****                   CHECK FOR COMPLETION OF BUFFERING
 2980 CALL UNLCCK(IUNIN,NSIN,LEN)
      IF(LEN.EQ.LBDSKI.AND.NSIN.EQ.0) GO TO 2990
      NUM=2990
      GO TO 4000
 2990 CONTINUE
      IF(J.EQ.1)GO TO 3000
      CALL UNLCCK(IUNOT,NSOT,LEN)
      IF(LEN.EQ.LBDSKI.AND.NSOT.EQ.0) GO TO 3000
      NUM=3000
      GO TO 4000
 3000 CONTINUE
cccccccccccccccccc end of J latitude loop ccccccccccccccccccccccccccccccccc

CCCCCCCCCCCCCCCC calculate the field-aligned currents CCCCCCCCCCCCCCCCCC
        pii = 4.*atan(1.)
        re = 6.481e+6
        dlat = 5.
        dlon = 5.
        dthr = dlat*pii/180.
        dlonr = dlon*pii/180.
C  zjmx goes from south to north
        do j=1,zjmx/2
        colat = 92.5*pii/180. - dthr*j
        st(j) = sin(90.*pii/180. - colat)
        st(zjmx+1-j) = st(j)
        enddo

      do j=1,zjmx
        do i=1,imax+1       
          famie(i,j) = 0.
          fwind(i,j) = 0.
          ftot(i,j) = 0.
        enddo
      enddo

        do 5150 j=2, zjmx-1
C     print *,'ADVNCE: calculate field-aligned current for j=',j
           j1 = zjmx+1-j
C          if (abs(st(j)) .lt. 1.e-20) then
C             write(6,"('advnce: glat,st(j)=',f6.2,f7.4)")
C    +                -92.5+dlat*j,st(j)
C             go to 5150
C           endif
           if (abs(sin(dipmag(i,j))) .lt. 1.e-5) then
              write(6,"('advnce: glat,sndip(i,j)=',f6.2,f7.4)")
     +                -92.5+dlat*j,sin(dipmag(i,j))
              go to 5150
            endif
          do 5145 i=1, imax+1
            lonm = i -1
            if (lonm .eq. 0) lonm = imax+1
            lonp = i + 1
            if (lonp .eq. imax+2) lonp = 1
C  Field-aligned current from wind
            fwind(i,j) = ((fwindu(lonp,j)-fwindu(lonm,j))/
     1        dlonr + (st(j+1)*fwindv(i,j+1)-st(j-1)*
     2        fwindv(i,j-1))/dthr) / (2.*re*st(j)) /
     3        (-sin(dipmag(i,j)))
C  Field-aligned current from ion drift
            famie(i,j) = ((famieu(lonp,j)-famieu(lonm,j))/
     1        dlonr + (st(j+1)*famiev(i,j+1)-st(j-1)*
     2        famiev(i,j-1))/dthr) / (2.*re*st(j)) /
     3        (-sin(dipmag(i,j)))
C  Total field-aligned currents
            ftot(i,j) = famie(i,j) + fwind(i,j)
5145       continue
5150     continue
        do 5155 i=1, imax+1
           fwind(i,1) = (fwind(i,2)-fwind(i,1)) / re/dthr
           fwind(i,zjmx) = (fwind(i,zjmx)-fwind(i,zjmx-1))/re/dthr
           famie(i,1) = (famie(i,2)-famie(i,1)) / re/dthr
           famie(i,zjmx) = (famie(i,zjmx)-famie(i,zjmx-1))/re/dthr
           ftot(i,1) = ftot(i,1) + ftot(i,1)
           ftot(i,zjmx) = ftot(i,zjmx) + ftot(i,zjmx)
5155     continue
CCCCCCCCCCCCCCCCCCC end of FAC calculation CCCCCCCCCCCCCCCCCCCCC



C     ****                   TRANSFER LINE JMAX TO DISK
      CALL FLIP
      CALL UNLCWT(IUNOT,NJOT,LBDSKI,NSOT)
      CALL UNLCCK(IUNOT,NSOT,LEN)
      IF(LEN.EQ.LBDSKI.AND.NSOT.EQ.0) GO TO 3100
      NUM=3100
      GO TO 4000
 3100 ITEMP=IUNIN
      IUNIN=IUNOT
      IUNOT=ITEMP
      REWIND IUNIN
      REWIND IUNOT
c     JTOTN=JTOTN*C(81)*C(55)*C(3)*C(53)*C(1)*C(2)
c     IF(NHEMI.EQ.0)JTOTS=JTOTS*C(81)*C(55)*C(3)*C(53)*C(1)*C(2)
C     ****
C     ****     SET UP AND SOLVE DYNAMO EQUATIONS
C     ****
      CALL DYN

C  Copy over FAC from half globe to full globe in magnetic coords
	do i=1,imaxmp
	do j=1,jmaxmh
	factm(i,j) = fact(i,j)
	enddo
	js = jmaxmh
	do j=jmaxmh+1,jmaxm
	js = js - 1
	factm(i,j) = fact(i,js)
	enddo
	enddo
C  Transform FAC from magnetic to geographic grid (including poles)
      CALL GRDTRP(factm(1,1),factg(1,0),IM(1,0),JM(1,0),
     1  DIM(1,0),DJM(1,0),IMAXGP,IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM)

C  Set values at the poles for Ped/Hall/Joule/TEC ht integration
	      sumts = 0.
	      sumtn = 0.
	      sumjs = 0.
	      sumjn = 0.
	      sumps = 0.
	      sumpn = 0.
	      sumhs = 0.
	      sumhn = 0.
	      sumqwn = 0.
	      sumqws = 0.
	      sumqan = 0.
	      sumqas = 0.
	      sumwkn = 0.
	      sumwks = 0.
	      sumwtn = 0.
	      sumwts = 0.
	      sumwfn = 0.
	      sumwfs = 0.
	      sumwfwn = 0.
	      sumwfws = 0.
	      sumwfan = 0.
	      sumwfas = 0.
	      do i=1,imax
	       sumts = sumts + tec(i,1)
	       sumtn = sumtn + tec(i,36)
	       sumjs = sumjs + qjn(i,1)
	       sumjn = sumjn + qjn(i,36)
	       sumps = sumps + sigp(i,1)
	       sumpn = sumpn + sigp(i,36)
	       sumhs = sumhs + sigh(i,1)
	       sumhn = sumhn + sigh(i,36)
	       sumqws = sumqws + qwind(i,1)
	       sumqwn = sumqwn + qwind(i,36)
	       sumqas = sumqas + qamie(i,1)
	       sumqan = sumqan + qamie(i,36)
	       sumwks = sumwks + work(i,1)
	       sumwkn = sumwkn + work(i,36)
	       sumwts = sumwts + wtot(i,1)
	       sumwtn = sumwtn + wtot(i,36)
	       sumwfs = sumwfs + ftot(i,1)
	       sumwfn = sumwfn + ftot(i,36)
	       sumwfas = sumwfas + famie(i,1)
	       sumwfan = sumwfan + famie(i,36)
	       sumwfws = sumwfws + fwind(i,1)
	       sumwfwn = sumwfwn + fwind(i,36)
              enddo
	      do i=1,imax+1
		tec(i,0) = sumts/72.
		tec(i,37) = sumtn/72.
		qjn(i,0) = sumjs/72.
		qjn(i,37) = sumjn/72.
		sigp(i,0) = sumps/72.
		sigp(i,37) = sumpn/72.
		sigh(i,0) = sumhs/72.
		sigh(i,37) = sumhn/72.
		qwind(i,0) = sumqws/72.
		qwind(i,37) = sumqwn/72.
		qamie(i,0) = sumqas/72.
		qamie(i,37) = sumqan/72.
		work(i,0) = sumwks/72.
		work(i,37) = sumwkn/72.
		wtot(i,0) = sumwts/72.
		wtot(i,37) = sumwtn/72.
		ftot(i,0) = sumwfs/72.
		ftot(i,37) = sumwfn/72.
		famie(i,0) = sumwfas/72.
		famie(i,37) = sumwfan/72.
		fwind(i,0) = sumwfws/72.
		fwind(i,37) = sumwfwn/72.
              enddo

C  Calculate global values of Conductances/Joule heat/TEC/FAC on geog grid
C       sud = re*re*dlon*dlat*dtr*dtr
	sud = c(51)*c(51)*c(1)*c(2)
	sums = 0.
	sumn = 0.
	do i=1,imax
	  do j=1,jmax/2
	    jn = j+jmax/2
	    sums = sums + sud*cs(j)
	    sumn = sumn + sud*cs(jn)
          enddo
        enddo
	      sumfs = 0.
	      sumfn = 0.
	      sumts = 0.
	      sumtn = 0.
	      sumjs = 0.
	      sumjn = 0.
	      sumps = 0.
	      sumpn = 0.
	      sumhs = 0.
	      sumhn = 0.
      do j=1,jmax/2
	jn = j+jmax/2
	abbs = cs(j) * sud / sums
	abbn = cs(jn) * sud / sumn
	      do i=1,imax
	       sumfs = sumfs + factg(i,j)*abbs
	       sumfn = sumfn + factg(i,jn)*abbn
	       sumts = sumts + tec(i,j)*abbs
	       sumtn = sumtn + tec(i,jn)*abbn
	       sumjs = sumjs + qjn(i,j)*abbs
	       sumjn = sumjn + qjn(i,jn)*abbn
	       sumps = sumps + sigp(i,j)*abbs
	       sumpn = sumpn + sigp(i,jn)*abbn
	       sumhs = sumhs + sigh(i,j)*abbs
	       sumhn = sumhn + sigh(i,jn)*abbn
              enddo
      enddo
	write (6,"(1x,'SH FAC/TEC/QJ/SIGP,H=',5e12.4)") sumfs,sumts,
     |   sumjs,sumps,sumhs
	write (6,"(1x,'NH FAC/TEC/QJ/SIGP,H=',5e12.4)") sumfn,sumtn,
     |   sumjn,sumpn,sumhn

C  Plot out height integrated Ped/Hall conductance on geographic grid
C       CALL EZCNTR (sigp(1,1),73,18)
C       CALL EZCNTR (sigh(1,1),73,18)
C  Convert height integrated Ped/Hall conductance from geographic to mag
C  ig(imaxmp,jmaxm), jg(imaxmp,jmaxm), wt(4,imaxmp,jmaxm)
c     do 3010 k=1,jmaxm
c       CALL GRDINT (sigpm(1,k),sigp(1,0),ig,jg,wt,imaxgp,imaxmp,
c    |   imaxg,jmaxg+2,imaxm,jmaxm,k)
c       sigpm(imaxmp,k) = sigpm(1,k)
c       CALL GRDINT (sighm(1,k),sigh(1,0),ig,jg,wt,imaxgp,imaxmp,
c    |   imaxg,jmaxg+2,imaxm,jmaxm,k)
c       sighm(imaxmp,k) = sighm(1,k)
c3010 continue
C  Plot out on magnetic grid
C       CALL EZCNTR (sigpm,imaxmp,jmaxmh)
C       CALL EZCNTR (sighm,imaxmp,jmaxmh)

C  Write out extra parameters for later processing
      ihrut = secs/3600.
      imnut = (secs-float(ihrut)*3600.+29.)/60.
      write (uthdr,"(2i4,1x,i2.2,':',i2.2,' height-integrated')")
     | iyear,iiday,ihrut,imnut
      mtime(1) = iiday
      mtime(2) = ihrut
      mtime(3) = imnut
c
c     subroutine wrxdr(flnm,f,nx,ny,xx,yy,xlab,ylab,lab1,
c    +  lab2,lab3,lab4)
c     subroutine wrxdr(flnm,f,nx,ny,xx,yy,xlab,ylab,lab1,
c    +  lab2,lab3,lab4,mtime,iclose)
c
c
C     ****                   CHECK FOR EXIT
C     ****                   CHECK FOR SAVE TAPE WRITE
c     write(6,"('ADVNCE iter',i5,': ihis=',i2,
c    +  ' isav=',i2,' sechis=',i2,' secsav=',i2,' secstrt,stop=',
c    +  2i5)") iter,ihis,isav,ihissech,isavsech,isecstart,isecstop
 3130 if (isav.eq.0.or.isavsech.eq.0) goto 3150
      if (ihis.ne.0.and.ihissech.ne.0) goto 3140
      if (iter.ge.nstp) then
        isav = 0
        isavsech = 0
      endif
      goto 3150
 3140 if (iter.lt.nstp) goto 100
      isav = 0
      isavsech = 0
 3150 TIME = SECOND(X)-TIME
      TIMET = TIMEF(X)-TIMET

!     call wrxdr('tec.xdr',tec,zimxp1,jmax+2,gcmlon,gcmlat,
!    |  'LONGITUDE (DEG)','LATITUDE (DEG)',uthdr,
!    |  'Total Electron Content (cm-2)',blank,blank,mtime,0)
!     call wrxdr('qjn.xdr',qjn,zimxp1,jmax+2,gcmlon,gcmlat,
!    |  'LONGITUDE (DEG)','LATITUDE (DEG)',uthdr,
!    |  'Joule Heat (ergs/cm2)',blank,blank,mtime,0)
!     call wrxdr('sigp.xdr',sigp,zimxp1,jmax+2,gcmlon,gcmlat,
!    |  'LONGITUDE (DEG)','LATITUDE (DEG)',uthdr,
!    |  'Pedersen Conductance (mhos)',blank,blank,mtime,0)
!     call wrxdr('sigh.xdr',sigh,zimxp1,jmax+2,gcmlon,gcmlat,
!    |  'LONGITUDE (DEG)','LATITUDE (DEG)',uthdr,
!    |  'Hall Conductance (mhos)',blank,blank,mtime,0)
!     call wrxdr('qwind.xdr',qwind,zimxp1,jmax+2,gcmlon,gcmlat,
!    |  'LONGITUDE (DEG)','LATITUDE (DEG)',uthdr,
!    |  'WIND Joule Heat (mW/m2)',blank,blank,mtime,0)
!     call wrxdr('qamie.xdr',qamie,zimxp1,jmax+2,gcmlon,gcmlat,
!    |  'LONGITUDE (DEG)','LATITUDE (DEG)',uthdr,
!    |  'AMIE Joule Heat (mW/m2)',blank,blank,mtime,0)
!     call wrxdr('wtot.xdr',wtot,zimxp1,jmax+2,gcmlon,gcmlat,
!    |  'LONGITUDE (DEG)','LATITUDE (DEG)',uthdr,
!    |  'TOTAL E&M POWER (mW/m2)',blank,blank,mtime,0)
!     call wrxdr('work.xdr',work,zimxp1,jmax+2,gcmlon,gcmlat,
!    |  'LONGITUDE (DEG)','LATITUDE (DEG)',uthdr,
!    |  'MECHANICAL WORK (mW/m2)',blank,blank,mtime,0)
!     call wrxdr('ftot.xdr',ftot,zimxp1,jmax+2,gcmlon,gcmlat,
!    |  'LONGITUDE (DEG)','LATITUDE (DEG)',uthdr,
!    |  'Total Field Aligned Current (A/m2)',blank,blank,mtime,0)
!     call wrxdr('famie.xdr',famie,zimxp1,jmax+2,gcmlon,gcmlat,
!    |  'LONGITUDE (DEG)','LATITUDE (DEG)',uthdr,
!    |  'Field Aligned Current due to Ion Drift(A/m2)',blank,blank,
!    |  mtime,0)
!     call wrxdr('fwind.xdr',fwind,zimxp1,jmax+2,gcmlon,gcmlat,
!    |  'LONGITUDE (DEG)','LATITUDE (DEG)',uthdr,
!    |  'Field Aligned Current due to Wind (A/m2)',blank,blank,mtime,0)

      WRITE(6,3170) TIMET
 3170 format('TIME FROM TIMEF CALL = ',e12.5)
      WRITE(6,3160) ITER,TIME
 3160 format('ADVANCE EXITED AT ITER = ',i6,' TIME IN INNER LOOP = ',
     +  e12.5,' SECS')
      RETURN
 4000 WRITE(6,4005) ITER
 4005 FORMAT(*1ERROR EXIT AT ITER =*I6,5X,*RESTART*)
      WRITE(6,4010)NUM,LEN,NSIN,NSOT,NSIN1,NSIN2,J
 4010 FORMAT(*0ADVANCE- I/O PROBLEM AT NUM=*I5,2X,*LEN=*I5,2X,*NSIN=*I5,
     A2X,*NSOT=*I5,2X,*NSIN1=*I5,2X,*NSIN2=*I5,2X,*J=*I2)
      CALL EXIT
      END
C
