#include "slgcm.main"
Status: R

#include "slgcm.inout"
#include "slgcm.interp"
#include "slgcm.zangs"
#include "slgcm.misc"
#include "slgcm.dirsol"
#include "slgcm.rdn"
#include "slgcm.ras"
#include "slgcm.ut"
#include "slgcm.fill"
#include "slgcm.geopl"
C
      SUBROUTINE SLGCST(DT, RSFLG, INIT, FIRST, TAU, DIURNL, SEASNL)
C
#include "slgcm.slgprm"
#include "slgcm.slgcon"
#include "slgcm.slgprg"
#include "slgcm.slclok"
#include "slgcm.slgclm"
#include "slgcm.slgout"
#include "slgcm.slgwrk"
C
C
      LOGICAL RSFLG, INIT, CUMFRI, DIURNL, SEASNL, FIRST
      double precision TAU
C
      DIMENSION SLUH(IJM0M0), SLVH(IJM0M0), SPUH(IJM0M0), SPVH(IJM0M0)
     *,         SLAM(IJM0M0), SPHI(IJM0M0), SSIG(IJM0M0)
     *,         CBV(IJM0M0),  SBV(IJM0M0)
C
      DIMENSION  BL(IJM0M0,LM), BT(IJP3M0,LM)
     *,         BLS(IJM0M0,LM), BTS(IJM0M0,LM)
     *,         HK1(IJM0M0,LM), GHT(IJM0M0,LM), DHT(IJM0M0,LM)
     *,        UNSS(IJM0M0,LM), GTP(IJM0M0,LM)
C
     *,         UST(IJP3M0,LM), VST(IJP3M0,LM),   SDS(IJP3M0,LM)
     *,         QST(IJP3M0,LM), YUUPT(IJP3M0,LM), YVVPT(IJP3M0,LM)
      DIMENSION TX1(IJP3M0),      TX2(IJP3M0),      TX3(IJP3M0)
C
      DIMENSION UTS(IJM0M0),    VTS(IJM0M0), STS(IJM0M0)
     *,         UMH(IJM0M0),    VMH(IJM0M0), SMH(IJM0M0)
     *,       AUUPT(IJM0M0),  AVVPT(IJM0M0)
     *,      YUUPTS(IJM0M0), YUVPTS(IJM0M0)
     *,      YVUPTS(IJM0M0), YVVPTS(IJM0M0)
C
      DIMENSION  UPN(IJM0M0,LM),  VPN(IJM0M0,LM),  PPN(IJM0M0,LM)
     *,         UPNO(IJM0M0),     VPNO(IJM0M0),    PPNO(IJM0M0)
C
      DIMENSION  WF1(IJM0M0), EA(IJM0M0)
C
      DIMENSION  IWK(IJM0M0), IPR(IJM0M0), JPR(IJM0M0), KPR(IJM0M0)
      DIMENSION  C1A(IJM0M0), C2A(IJM0M0), C3A(IJM0M0), C4A(IJM0M0)
      DIMENSION  C1B(IJM0M0), C2B(IJM0M0), C3B(IJM0M0), C4B(IJM0M0)
C
      DIMENSION  WKU(IJM0M0), WKV(IJM0M0), WKT(IJM0M0)
      DIMENSION  PPNT(IJM0M0), PPNB(IJM0M0), WKQ(IJM0M0)
C
      LOGICAL    STABLE(IJM0M0)
      DIMENSION  ALF(IJM0M0,LM), BET(IJM0M0,LM), GAM(IJM0M0,LM)
     *,          HOL(IJM0M0,LM), ETA(IJM0M0,LM), CLT(IJM0M0,LM)
     *,          HST(IJM0M0,LM), GMH(IJM0M0,LM), QOL(IJM0M0,LM)
     *,          WFN(IJM0M0),    AKM(IJM0M0), EKB(IJM0M0)
     *,          VTC(IJM0M0),    VTD(IJM0M0), VTE(IJM0M0)
     *,          VTF(IJM0M0),    VTG(IJM0M0)
     *,          IAR(IJM0M0),    I1R(IJM0M0), I2R(IJM0M0)
     *,          UOC(IJM0M0,LM), VOC(IJM0M0,LM)
     *,          PKI(IJM0M0)
C
      DIMENSION  TM1SP(LM),  TM1NP(LM), LMIN(LM), LMAX(LM)
      dimension  RASAL(LMM1)
C
      DIMENSION   RHS(IJM0M0,2), PHI(IJM0M0,2), WSV(13*IM), IX(IM)

C
      EQUIVALENCE (C1A(1),WKU(1)),     (C2A(1),WKV(1)), (C3A(1),WKT(1))
     *,           (C4A(1),PPNT(1)),    (C1B(1),PPNB(1)),(C2B(1),WKQ(1))
      EQUIVALENCE
     *            (TN0(1,1),   BTS(1,1))
     *,           (DIV(1,1),   UNSS(1,1))
C
      EQUIVALENCE
     *            (UMH(1),   PHI(1,1))
     *,           (VMH(1),   PHI(1,2))
     *,           (SMH(1),   RHS(1,1))
     *,           (sluh(1),   RHS(1,2))
     *,           (AUUPT(1), YUUPTS(1))
     *,           (WF1(1),   YUUPTS(1))
     *,           (AVVPT(1), YUVPTS(1))
     *,           (EA(1),    YUVPTS(1))
     *,           (CBV(1),   YVUPTS(1))
     *,           (SBV(1),   YVVPTS(1))
     *,           (sluh(1),   slvh(1))
     *,           (spuh(1),   spvh(1))
C
      EQUIVALENCE (WKX(1,13),          UPN(1,1))
     *,           (WKX(1,13+LM),       VPN(1,1))
     *,           (WKX(1,13+LM+LM),    PPN(1,1))
     *,           (WKX(1,13+LM*3),     BLS (1,1))
     *,           (WKX(1,13+LM*4),     GTP (1,1))
      EQUIVALENCE (WKY(1,1),           UST(1,1))
     *,           (WKY(1,1+LM),        VST(1,1))
     *,           (WKY(1,1+LM+LM),     SDS(1,1))
     *,           (WKY(1,1+LM+LM+LM),  YUUPT(1,1))
     *,           (WKY(1,1+LM*4),      YVVPT(1,1))
     *,           (WKY(1,1+LM*5),      QST(1,1))
      
      EQUIVALENCE (GTP(1,1),   BL(1,1))
     *,           (PPN(1,1),   HK1(1,1))
     *,           (GHT(1,1),   YUUPT(1,1))
     *,           (DHT(1,1),   YVVPT(1,1))
C
      EQUIVALENCE (UPN(1,1),   ALF(1,1))
     *,           (VPN(1,1),   BET(1,1))
     *,           (PPN(1,1),   GAM(1,1))
     *,           (SDS(1,1),   HOL(1,1))
     *,           (YUUPT(1,1), ETA(1,1))
     *,           (YVVPT(1,1), CLT(1,1))
     *,           ( BLS(1,1),  GMH(1,1))
     *,           ( UST(1,1),  HST(1,1))
     *,           ( VST(1,1),  QOL(1,1))
     *,           ( QST(1,1),  UOC(1,1))
     *,           ( GTP(1,1),  VOC(1,1))
C
      EQUIVALENCE (UPNO(1), VTE(1))
     *,           (VPNO(1), VTF(1))
     *,           (PPNO(1), VTG(1))
C
C
      LOGICAL INIT1, INIT2
      DATA INIT1/.TRUE./, INIT2/.TRUE./, CUMFRI/.TRUE.  /
      DATA NSTEP/0/
      data ia/158/
      SAVE
C
      NSTEP = NSTEP + 1
C
      CALL WRAP(1,   IJM0M0, LGP)
      CALL WRAP(1,   IJM0M0, PIK)
C
      DO 2002 L=1,LM
C
      CALL WRAP(I1P1,IJM1M0, UN0(1,L))
      CALL WRAP(1,   IJM1M0, VN0(1,L))
      CALL WRAP(1,   IJM0M0, TN0(1,L))
      CALL WRAP(1,   IJM0M0, QN0(1,L))
      CALL WRAP(1,   IJM0M0, SGD(1,L))
      CALL WRAP(1,   IJM0M0, DIV(1,L))
C
      DO 2002 I=1,IJM0M0
      TN0(I,L) = TN0(I,L) * PIK(I)*(PRH(L)/TBR(L)) - 1.0
 2002 CONTINUE
C
C
      IF (FIRST .OR. INIT) THEN
         PRINT *,' DRCONV=',DRCONV,' DIFT=',DIFT,' SHPFLT=',SHPFLT
     *,          ' DL1=',DL1,' DL2=',DL2,' SURFR=',SURFR
     *,          ' FOURPT=',FOURPT,' TRIQAD=',TRIQAD,' TRICUB=',TRICUB
     *,          ' RSFLG=',RSFLG,' VDIFQ=',VDIFQ
     *,          ' ADIBAT=',ADIBAT,' FORWRD=',FORWRD,' DIGFIL='
     *,            DIGFIL, ' CDA=',CDA,' CDB=',CDB,' ADMASS=',ADMASS
     *,          ' TEXF1=',TEXF1,' TEXF2=',TEXF2
     *,          ' TEXFa=',TEXFa,' TEXFb=',TEXFb
C
         RASAL(:) = RASALF
C
         CALL WRAP(1,   IJM0M0, LGPO)
         DO 2003 I=1,IJM0M0
         WK1(I) = (0.01 * EXP(LGPO(I)+PBR)) ** RKAP
 2003    CONTINUE
C
         DO 5 L=1,LM
         CALL WRAP(I1P1,IJM1M0, UN1(1,L))
         CALL WRAP(1,   IJM1M0, VN1(1,L))
         CALL WRAP(1,   IJM0M0, TN1(1,L))
         CALL WRAP(1,   IJM0M0, SGDO(1,L))
         CALL WRAP(1,   IJM0M0, DIVO(1,L))
C
         CALL FILPU(UN1(1,L))
         CALL FILPV(VN1(1,L))
C
         DO 5 I=1,IJM0M0
         TN1(I,L) = TN1(I,L) * WK1(I)*(PRH(L)/TBR(L)) - 1.0
    5    CONTINUE
C
         IF (.NOT. RSFLG) THEN
C
            DO 3015 I=1,IJM0M0
            WKA(I) = PIB(I) * 100.0
            WK2(I) = 1.0 / PIK(I)
            WK3(I) = 1.0 / WK1(I)
 3015       CONTINUE
            DO 3020 L=1,LM
            DO 3020 I=1,IJM0M0
            TN0(I,L) = (1.0 + TN0(I,L)) * WK2(I) * (TBR(L) / PRH(L))
            TN1(I,L) = (1.0 + TN1(I,L)) * WK3(I) * (TBR(L) / PRH(L))
 3020       CONTINUE
C
            CALL DCADJ2(IJM0M0, IJM0M0, LM, CMP, CPM, AMP, APM
     *,                        QMP, QPM, TN0, QN0, WK5, WK6, DT)
C
            DO 3030 L=1,LM
            DO 3030 I=1,IJM0M0
            TN0(I,L) = TN0(I,L) * PIK(I) * (PRH(L) / TBR(L)) - 1.0
            TN1(I,L) = TN1(I,L) * WK1(I) * (PRH(L) / TBR(L)) - 1.0
 3030       CONTINUE
C
C           PRELIMINARY CALCULATION OF DIVERGENCE AND SIGDOT
C
            DO 16 I=1,IJM0M0
            WK6(I) = 0.0
            WK0(I) = 0.0
   16       CONTINUE
            DO 17 I=1,IJM2M1
            WK6(I1+I) = (WKA(I1+I) + WKA(I1+I+1)) * HALF
   17       CONTINUE
            DO 18 I=I1P1,IJM1M0,IMP3
            WK6(IMP2+I) = WK6(2+I)
   18       CONTINUE
            DO 19 I=1,IJM1M0
            WK7(I) = (WKA(I1+I) + WKA(I)) * HALF
   19       CONTINUE
C
            DO 50 L=1,LM
            DO 20 I=1,IJM1M0
            WK8(I) = UN0(I,L) * WK6(I)
            WK9(I) = VN0(I,L) * WK7(I)
   20       CONTINUE
            DO 21 I=2,IJM2M0
            WK1(I1+I) = UN0(I1+I,L) - UN0(I1+I-1,L)
            WK0(I1+I) = WK8(I1+I  ) - WK8(I1+I-1  )
   21       CONTINUE
            DO 25 I=I1P1,IJM1M0,IMP3
            WK1(I) = WK1(IM+I)
            WK0(I) = WK0(IM+I)
   25       CONTINUE
            DO 30 I=1,IJM1M0
            WK5(I) = VN0(I,L) * CPHV(I)
            WK8(I) = WK9(I  ) * CPHV(I)
   30       CONTINUE
            DO 35 I=I1P1,IJM1M0
            DIV(I,L) = WK1(I) * CDU(I) + (WK5(I) - WK5(I-I1)) * CDV(I)
            SGD(I,L) = WK0(I) * CDU(I) + (WK8(I) - WK8(I-I1)) * CDV(I)
   35       CONTINUE
            call avrx3(div(1,l),sgd(1,l))
C
C     CALCULATE THE DIVERGENCE AND SIGDOT AT THE POLES
C
            DNPOLE = 0.0
            DSPOLE = 0.0
            SNPOLE = 0.0
            SSPOLE = 0.0
            DO 36 I=1,IM
            DNPOLE = DNPOLE + VN0(IJM2M0+I,L)
            DSPOLE = DSPOLE + VN0(I,L)
            SNPOLE = SNPOLE + WK9(IJM2M0+I)
            SSPOLE = SSPOLE + WK9(I)
   36       CONTINUE
            TEM    = 4.0 / (DY*IM)
            DNPOLE = DNPOLE * TEM
            DSPOLE = DSPOLE * TEM
            SNPOLE = SNPOLE * TEM
            SSPOLE = SSPOLE * TEM
C
            DO 40 I=1,IMP3
            DIV(I+IJM1M0,L) = -DNPOLE
            DIV(I,L)        =  DSPOLE
            SGD(I+IJM1M0,L) = -SNPOLE
            SGD(I,L)        =  SSPOLE
   40       CONTINUE
C
            CALL WRAP(1,IJM0M0,DIV(1,L))
C
   50       CONTINUE
C
            DO 60 L=2,LM
            DO 55 I=1,IJM0M0
            SGD(I,L) = SGD(I,L-1) + SGD(I,L)
   55       CONTINUE
   60       CONTINUE
            DO 80 L=1,LM
            DO 65 I=1,IJM0M0
            SGD(I,L) = (SGD(I,LM)*SIG(L+1) - SGD(I,L)) / WKA(I)
   65       CONTINUE
            DO 75 I=1,IJM0M0
            sgdo(I,L) = sgd(I,L)
            divo(I,L) = div(I,L)
   75       CONTINUE
   80       CONTINUE
            sgd (:,LM) = 0.0
            sgdo(:,LM) = 0.0
C
            GBR(:) = 0.0
            CALL GLMNP2(IMP3, JM, PHIS, DX(1,2), DY
     *,                                ARE2, DX(1,1), PHISM)
            GBR(:) = PHISM
C
            PRINT *,' GBR=',GBR, 'PBR=',PBR
C
         ENDIF
         IF (.NOT. INIT) FIRST = .FALSE.
      ENDIF
C
      CALL GETGTP(IJM0M0, LM, DSL, TBR, RGAS, GBR, PHIS, TN0
     *,                                       LGP, GTP, WK1)
C
C
c
C     START OF THE MAIN LOOP
C
C
c     IF (DRCONV) THEN
c        itau = tau
c        if (mod(itau,3600) .eq. 0) then
c        DO 3015 I=1,IJM0M0
c        WK2(I) = (0.01 * EXP(LGP(I)+PBR)) ** RKAP
c        WK1(I) = 1.0 / WK2(I)
c3015    CONTINUE
c        DO 3020 L=1,LM
c        DO 3020 I=1,IJM0M0
c        TN0(I,L) = (1.0 + TN0(I,L)) * WK1(I) * (TBR(L) / PRH(L))
c3020    CONTINUE
C
c        CALL DCADJ(IJM0M0, IJM0M0, LM, CMP, CPM, AMP, APM
c    *,                        TN0, WK5, WK6, DT)
C
C
C        CALL BICK (DSG, PRJ, PRH, CPHU, CPHV, TBR
C    *,                   TN0, WK1, WK3, PPNT, PPNB)
C
C        DO 3029 I=1,IJM0M0
C        TN0(I,1) = TN0(I,1)  + PPNT(I) * DT
C        TN0(I,LM) = TN0(I,LM)  + PPNB(I) * DT
C3029    CONTINUE
C
C
c        DO 3030 L=1,LM
c        DO 3030 I=1,IJM0M0
c        TN0(I,L) = TN0(I,L) * WK2(I) * (PRH(L) / TBR(L)) - 1.0
c3030    CONTINUE
c        ENDIF
c     ENDIF
C
C
      DO 205 J=1,JMM1
      II = (J-1)*IMP3
      DO 205 I=1,IMP3
      WK7(II+I) = FCOR(J,2)
      WK8(II+I) = FCOR(J,1)
  205 CONTINUE
C
C    COMPUTE FOR N(V)
C
      DO 207 I=I1P1,IJM1M0
      VTC(I) = (LGPO(I) - LGPO(I+1)) * CDU(I)
      WK9(I) = (LGP (I) - LGP (I+1)) * CDU(I)
  207 CONTINUE
      tem = 0.5 / DY
      DO 208 I=1,IJM1M0
      VTD(I) = (LGPO(I) - LGPO(I1+I)) * tem
      WK0(I) = (LGP (I) - LGP (I1+I)) * tem
  208 CONTINUE
C
C
      DO 500 L=1,LM
      CALL FILPU(UN0(1,L))
      CALL FILPV(VN0(1,L))
C
      IF (DIFT) THEN
         CALL DELSQR (CPHU, CPHV, UN1(1,L), UPNO)
         CALL DELSQV (CPHU, CPHV, VN1(1,L), VPNO)
         CALL DELSQR (CPHU, CPHV, TN1(1,L), PPNO)
C
         CALL DELSQR (CPHU, CPHV, UN0(1,L), UPN(1,L))
         CALL DELSQV (CPHU, CPHV, VN0(1,L), VPN(1,L))
         CALL DELSQR (CPHU, CPHV, TN0(1,L), PPN(1,L))
 5213    CONTINUE
      ELSE
         UPNO(:) = 0.0
         VPNO(:) = 0.0
         PPNO(:) = 0.0
C
         UPN(:,L) = 0.0
         VPN(:,L) = 0.0
         PPN(:,L) = 0.0
      ENDIF
C
C CALCULATE V AT U-PTS  (SAVED IN WK3)
C
      CALL UPTV(IJM0M0, VN0(1,L), WK3, TX1, TX2, TX3, IJP3M0)

C
C**********************
C CALCULATE YU (IN YUUPT) AND NLAMBDA (IN UPN) AT U-PTS
C**********************
C
      DO 1226 I=1,IJM0M0
      WK5(I) = GTP(I,L) - DIV(I,L)*CDB
 1226 CONTINUE
C
      RTBR  = RGAS * TBR(L)
      RTBRH = RTBR * 0.5
C
      DO 217 I=I1P1,IJM1M1
      UPNO(I) = VTC(I) * (TN1(I,L)+TN1(I+1,L)) * RTBRH
     *        + UPNO(I)
C
      UPN(I,L) = WK9(I) * (TN0(I,L) + TN0(I+1,L)) * RTBRH
     *         + UPN(I,L)
C
      YUUPT(I+I2,L) = UN0(I,L) + (WK7(I)*WK3(I)
     *              - (WK5(I+1)-WK5(I))*CDU(I) + UPN(I,L))  * DL2DT
  217 CONTINUE
C
      CALL WRAP(I1P1,IJM1M0,YUUPT(I2+1,L))
      CALL FILPU(YUUPT(I2+1,L))
      CALL FILLU(YUUPT(1,L))

c     CALL WRTTEM(77, IMP3, JM, IM, yuupt(i2+1,l), WK1)
C
      CALL WRAP(I1P1,IJM1M0,UPN(1,L))
      CALL WRAP(I1P1,IJM1M0,UPNO)
      CALL FILPU(UPN(1,L))
      CALL FILPU(UPNO)
C
C CALCULATE U AT V-PTS (SAVED IN WK3)
C
       CALL VPTU(IJM0M0, UN0(1,L), WK3, TX1, TX2, TX3, IJP3M0)

C
C**********************
C CALCULATE YV (IN YVVPT) AND NTHETA (IN VPN)  AT V-PTS
C**********************

         DO 216 I=1,IJM1M0
         VPNO(I) = VTD(I) * (TN1(I1+I,L)+TN1(I,L)) * RTBR
     *           + VPNO(I)
C
         VPN(I,L) = WK0(I) * (TN0(I1+I,L) + TN0(I,L))* RTBR
     *            + VPN(I,L)
C
         YVVPT(I+I2,L) = VN0(I,L) - (WK8(I) * WK3(I)
     *          + (WK5(I1+I) - WK5(I)) * (1.0/DY) - VPN(I,L)) * DL2DT
  216    CONTINUE
C
      CALL WRAP(1,IJM1M0,YVVPT(I2+1,L))
      CALL FILPV(YVVPT(I2+1,L))
      CALL FILLV(YVVPT(1,L))
C
      CALL WRAP(1,IJM1M0,VPN(1,L))
      CALL WRAP(1,IJM1M0,VPNO)
      CALL FILPV(VPN(1,L))
      CALL FILPV(VPNO)

c     CALL WRTTEM(77, IMP3, JM, IM, yvvpt(i2+1,l), WK1)
C
C
C  CALCULATE BL, BT, NT
C
C
      IF (L .EQ. 1) THEN
         DO 301 I=1,IJM0M0
         TEM1    = DIVO(I,L) + SGDO(I,L)*ODS(L)
         TEM2    = SGDO(I,L) * HALF
         PPNO(I) = (TEM2 * GMA(L) - TEM1 * RKAP) * TN1(I,L)
     *               + PPNO(I)
C
C
         WK1(I  ) = DIV(I,L  ) + SGD(I,1) * ODS(L)
         BL(I,L)  = LGP(I)     - WK1(I) * DL2DT
C
         WK2(I  ) = SGD(I,1) * (GMA(L)*0.5)
C
         PPN(I,L) = (WK2(I) - WK1(I) * RKAP) * TN0(I,L)
     *            + PPN(I,L)
C    *            + PPNT(I)
C
         BT(I+I2,L)  = TN0(I,L) - LGP(I) * RKAP
     *                          + (WK2(I) + PPN(I,L)) * DL2DT
  301    CONTINUE
      ELSEIF (L .LT. LM) THEN
         DO 302 I=1,IJM0M0
         TEM1    = DIVO(I,L) + (SGDO(I,L) - SGDO(I,L-1))*ODS(L)
         TEM2    = (SGDO(I,L) + SGDO(I,L-1)) * HALF
         PPNO(I) = (TEM2 * GMA(L) - TEM1 * RKAP) * TN1(I,L)
     *               + PPNO(I)
C
         WK1(I  ) = DIV(I,L  ) + (SGD(I,L) - SGD(I,L-1)) * ODS(L)
         BL(I,L)  = LGP(I)     - WK1(I) * DL2DT
C
         WK2(I  ) = (SGD(I,L) + SGD(I,L-1)) * (GMA(L)*0.5)
C
         PPN(I,L) = (WK2(I) - WK1(I) * RKAP) * TN0(I,L)
     *            + PPN(I,L)
C
         BT(I+I2,L)  = TN0(I,L) - LGP(I) * RKAP
     *                          + (WK2(I) + PPN(I,L)) * DL2DT
  302    CONTINUE
      ELSE
         DO 303 I=1,IJM0M0
         TEM1    = DIVO(I,L) + (SGDO(I,L) - SGDO(I,L-1))*ODS(L)
         TEM2    = (SGDO(I,L) + SGDO(I,L-1)) * HALF
         PPNO(I) = (TEM2 * GMA(L) - TEM1 * RKAP) * TN1(I,L)
     *               + PPNO(I)
C
C
         WK1(I  ) = DIV(I,L  ) + (SGD(I,L) - SGD(I,L-1)) * ODS(L)
         BL(I,L)  = LGP(I)     - WK1(I) * DL2DT
C
         WK2(I  ) = (SGD(I,L) + SGD(I,L-1)) * (GMA(L)*0.5)
C
         PPN(I,L) = (WK2(I) - WK1(I) * RKAP) * TN0(I,L)
     *            + PPN(I,L)
C        PPN(I,L) = PPN(I,L) + PPNB(I)
C
         BT(I+I2,L)  = TN0(I,L) - LGP(I) * RKAP
     *                          + (WK2(I) + PPN(I,L)) * DL2DT
  303    CONTINUE
      ENDIF

      CALL WRAP(1, IJM0M0, BL(1,L))
      CALL WRAP(1, IJM0M0, BT(I2+1,L))
      CALL FILLP(BT(1,L))
C
C        TIME EXTRAPOLATION
C
C  SAVE U(N+1/2) AND V(N+1/2) IN ARRAYS UST AND VST FOR USE IN
C  INTERPOLATION, ETC
C
C
         DO 310 I=1,JSU1M0
         UST(I+I2,L)  = UN0(I,L) + TEXFa * (UN0(I,L) - UN1(I,L))
         VST(I+I2,L)  = VN0(I,L) + TEXFa * (VN0(I,L) - VN1(I,L))
         SDS(I+I2,L)  = SGD(I,L) + TEXFa * (SGD(I,L) - SGDO(I,L))
         QST(I+I2,L)  = QN0(I,L) 
         TEMT  = PPN(I,L) + TEXFb * (PPN(I,L) - PPNO(I))
         TEMU  = UPN(I,L) + TEXFb * (UPN(I,L) - UPNO(I))
         TEMV  = VPN(I,L) + TEXFb * (VPN(I,L) - VPNO(I))
C
         UN1(I,L)  = UN0(I,L)
         VN1(I,L)  = VN0(I,L)
         SGDO(I,L) = SGD(I,L)
C
         PPN(I,L)  = TEMT
         UPN(I,L)  = TEMU
         VPN(I,L)  = TEMV
  310    CONTINUE
C
         DO 311 II=1,JUP0M0
         I = II + JSU1M0
         UST(I+I2,L)  = UN0(I,L) + TEXF1 * (UN0(I,L) - UN1(I,L))
         VST(I+I2,L)  = VN0(I,L) + TEXF1 * (VN0(I,L) - VN1(I,L))
         SDS(I+I2,L)  = SGD(I,L) + TEXF1 * (SGD(I,L) - SGDO(I,L))
         QST(I+I2,L)  = QN0(I,L) 
         TEMT  = PPN(I,L) + TEXF2 * (PPN(I,L) - PPNO(I))
         TEMU  = UPN(I,L) + TEXF2 * (UPN(I,L) - UPNO(I))
         TEMV  = VPN(I,L) + TEXF2 * (VPN(I,L) - VPNO(I))
C
         UN1(I,L)  = UN0(I,L)
         VN1(I,L)  = VN0(I,L)
         SGDO(I,L) = SGD(I,L)
C
         PPN(I,L)  = TEMT
         UPN(I,L)  = TEMU
         VPN(I,L)  = TEMV
  311    CONTINUE
C
         DO 312 II=1,JNU1M0
         I = II + JNU*IMP3
         UST(I+I2,L)  = UN0(I,L) + TEXFa * (UN0(I,L) - UN1(I,L))
         VST(I+I2,L)  = VN0(I,L) + TEXFa * (VN0(I,L) - VN1(I,L))
         SDS(I+I2,L)  = SGD(I,L) + TEXFa * (SGD(I,L) - SGDO(I,L))
         QST(I+I2,L)  = QN0(I,L) 
         TEMT  = PPN(I,L) + TEXFb * (PPN(I,L) - PPNO(I))
         TEMU  = UPN(I,L) + TEXFb * (UPN(I,L) - UPNO(I))
         TEMV  = VPN(I,L) + TEXFb * (VPN(I,L) - VPNO(I))
C
         UN1(I,L)  = UN0(I,L)
         VN1(I,L)  = VN0(I,L)
         SGDO(I,L) = SGD(I,L)
C
         PPN(I,L)  = TEMT
         UPN(I,L)  = TEMU
         VPN(I,L)  = TEMV
  312    CONTINUE
C
      CALL FILPU(UST(I2+1,L))
      CALL FILPV(VST(I2+1,L))
C
      CALL FILLU(UST(1,L))
      CALL FILLV(VST(1,L))
      CALL FILLP(SDS(1,L))
      CALL FILLP(QST(1,L))
C
  500 CONTINUE
C
       DO 510 L=1,LM
       DO 510 I=1,IJM0M0
       DIVO(I,L) = DIV(I,L)
       TN1(I,L)  = TN0(L,L)
  510  CONTINUE
       DO 520 I=1,IJM0M0
       LGPO(I)   = LGP(I)
  520  CONTINUE
C
C
C  NOW 'CENTER IN SPACE'.
C
         do 530 l=1,lm
         lmin(l) = 1
         lmax(l) = lm
  530    continue   
C
C
      DO 600 L=1,LM

      IF (ROTATE) THEN
         DO 320 I=1,IJM0M0
         VTS(I) = VST(I+I2,L)*DTOA
  320    CONTINUE
         DO 321 I=1,JSU1M0
         UTS(I) = UST(I+I2,L)*DTOA
  321    CONTINUE
         DO 322 I=1,JNU1M0
         UTS(I+JNU*IMP3) = UST(I+JNU*IMP3+I2,L)*DTOA
  322    CONTINUE
         DO 323 I=1,JUP0M0
         UTS(JSU1M0+I) = UST(JSU1M0+I+I2,L)*RH1U(JSU1M0+I)*DT
  323    CONTINUE
      ELSE
         DO 324 I=1,IJM0M0
         VTS(I) = VST(I+I2,L)*DTOA
         UTS(I) = UST(I+I2,L)*RH1U(I)*DT
  324    CONTINUE
      ENDIF
C
      IF (L .EQ. 1) THEN
         DO 325 I=1,IJM0M0
         STS(I) = SDS(I+I2,L) * DTHF
  325    CONTINUE
      ELSE
         DO 326 I=1,IJM0M0
         STS(I) = (SDS(I+I2,L) + SDS(I+I2,L-1)) * DTHF
  326    CONTINUE
      ENDIF
C
C
C*******************************************************************
C THE ZANG* ROUTINES ARE FOR 3-D
C*******************************************************************

C
C CALCULATE COORDS FOR U-PT  USING GEOGRAPHICAL OR ROTATED COORD SYSTEM
C

         CALL ZANGUR(UTS, VTS, STS
     *,              CPHU, SPHU, PHU, XCR(1,2), SGB(L), SIG
     *,              SPUH, SLUH, SSIG, CBV, SBV
     *,              WK1, WK2, WK3, TX1, TX2, TX3)
C
C
         call clkon
c        CALL MINMXL(LM, IJM0M0, I1P1, NSTEP, SGB, SSIG, LMIN, LMAX)
C
         do 556 I=1,IJM0M0
         KPR(I) = LMIN(L) - 1
  556    continue
         DO 559 K=LMIN(L),LMAX(L)
         DO 559 I=I1P1,IJM1M0
         IF (SSIG(I) .GT. SGB(K))  KPR(I) = K
  559    CONTINUE
         call clkoff(kprtim)
C
C
         CALL BILL(UST, C1A, SPUH, SLUH, SSIG
     *,            IWK, IPR, JPR, KPR, WK1, WK2, WK3
     *,            SGX, SGB, 0.5, 1.0, 2, JMM1)
         CALL BILL(VST, C2A, SPUH, SLUH, SSIG
     *,            IWK, IPR, JPR, KPR, WK1, WK2, WK3
     *,            SGX, SGB, 1.0, 0.5, 2, JMM1)
         CALL BILL(SDS, SMH, SPUH, SLUH, SSIG
     *,            IWK, IPR, JPR, KPR, WK1, WK2, WK3
     *,            SGX, SGB, 1.0, 1.0, 2, JMM1)
C
         IF (ROTATE) THEN
            DO 610 I=1,JSU2M0
            UMH(I1+I) = C1A(I1+I) * CBV(I1+I) - C2A(I1+I) * SBV(I1+I)
            VMH(I1+I) = C2A(I1+I) * CBV(I1+I) + C1A(I1+I) * SBV(I1+I)
  610       CONTINUE
C
            DO 612 II=1,JUP0M0
            I = II + JSU1M0
            UMH(I) = C1A(I) / COS(SPUH(I))
            VMH(I) = C2A(I)
  612       CONTINUE
C
            DO 1614 I=1,JNU2M0
            UMH(I+JNU*IMP3) = C1A(I+JNU*IMP3) * CBV(I+JNU*IMP3)
     *                      - C2A(I+JNU*IMP3) * SBV(I+JNU*IMP3)
            VMH(I+JNU*IMP3) = C2A(I+JNU*IMP3) * CBV(I+JNU*IMP3)
     *                      + C1A(I+JNU*IMP3) * SBV(I+JNU*IMP3)
 1614       CONTINUE
C
            DO 614 I=1,IJM2M0
            UMH(I1+I) = UMH(I1+I) * DTOA
            VMH(I1+I) = VMH(I1+I) * DTOA
            SMH(I1+I) = SMH(I1+I) * DT
  614       CONTINUE
         ELSE
            DO 617 I=I1P1,IJM1M0
            UMH(I) = C1A(I) / COS(SPUH(I)) * DTOA
            VMH(I) = C2A(I) * DTOA
            SMH(I) = SMH(I) * DT
  617       CONTINUE
         ENDIF
C
         CALL FILPU(UMH)
         CALL FILPU(VMH)
         CALL FILPU(SMH)
C
         CALL WRAP(I1P1,IJM1M0,UMH)
         CALL WRAP(I1P1,IJM1M0,VMH)
         CALL WRAP(I1P1,IJM1M0,SMH)
C
C
      CALL ZANGU(UMH, VMH, SMH, ONE
     *,          CPHU, SPHU, PHU, XCR(1,2), SGB(L), SIG
     *,          SPHI,   SLAM, SSIG, WK1, WK2)
C
C
      call clkon
c     CALL MINMXL(LM, IJM0M0, I1P1, NSTEP, SGB, SSIG, LMIN, LMAX)
C
      do 756 I=1,IJM0M0
      KPR(I) = LMIN(L) - 1
  756 continue
      DO 759 K=LMIN(L),LMAX(L)
      DO 759 I=I1P1,IJM1M0
      IF (SSIG(I) .GT. SGB(K))  KPR(I) = K
  759 CONTINUE
      call clkoff(kprtim)
C
C
       do 761 I=1,IJM0M0
       yuupts(I) = -999.0
       yvupts(I) = -999.0
  761  CONTINUE
C
         CALL BIC1nu(YUUPT, YUUPTS, SPHI, SLAM, SSIG
     *,            IWK, IPR, JPR, KPR
     *,            C1A, C2A, C3A, C4A, C1B, C2B, C3B, C4B
     *,            WK1, WK2, WK3, WK4
     *,            SGX, SGB, S1A, 0.5, 1.0, 2, JMM1)
C
         CALL BIC1nu(YVVPT, YVUPTS, SPHI, SLAM, SSIG
     *,            IWK, IPR, JPR, KPR
     *,            C1A, C2A, C3A, C4A, C1B, C2B, C3B, C4B
     *,            WK1, WK2, WK3, WK4
     *,            SGX, SGB, S1A, 1.0, 0.5, 2, JMM1)
C
C
C CALCULATE AUUPT(I,J) AT U-POINTS USING THE FULLY CENTERED APPROACH
C
      DO 9602 J=2,JMM1
      II = (J-1)*IMP3
      DO 9602 I=1,IMP3
      WK1(II+I) = BLMIK(I,1) - SLUH(II+I)
      WK2(II+I) = BLMIK(I,1) - SLAM(II+I)
 9602 CONTINUE
C
      DO 765 I=I1P1,IJM1M0
      WK4(I) = COS(SPUH(I))
      WK5(I) = COS(WK1(I))
      WK6(I) = SIN(WK1(I))
      WK7(I) = COS(WK2(I))
      WK8(I) = SIN(WK2(I))
      WK2(I) = SIN(SPHI(I))
C
      WKA(I) = SPHU(I) * SIN(SPUH(I))
      WK0(I) = CPHU(I) * WK4(I)
      WK0(I) =  1.0    / (WKA(I) + WK0(I) * WK5(I))
C
      WK9(I) = CPHU(I) * WK8(I) * WK4(I) * WK6(I) * WK0(I) + WK7(I)
C
      WK2(I) =  WK2 (I) * WK8(I)
     *       + (SPHU(I) * COS(SPHI(I)) - CPHU(I) * WK2(I) * WK7(I))
     *                  * WK6(I) * WK4(I) * WK0(I)
C
      AUUPT(I) =  WK9(I) * YUUPTS(I)  +  WK2(I) * YVUPTS(I)
     *         +  DL1DT  * UPN(I,L)
C
 765  CONTINUE
C
      CALL WRAP(I1P1,IJM1M0,AUUPT)
      CALL FILPU(AUUPT)
C
c     CALL WRTTEM(77, IMP3, JM, IM, yuupts, WK1)
c     CALL WRTTEM(77, IMP3, JM, IM, yvupts, WK1)
c     CALL WRTTEM(77, IMP3, JM, IM, upn(1,l), WK1)
c     CALL WRTTEM(77, IMP3, JM, IM, AUUPT, WK1)
C
C
C CALCULATE COORDS FOR V-PT  USING GEOGRAPHICAL OR ROTATED COORD SYSTEM
C
         CALL ZANGVR(UTS, VTS, STS
     *,              CPHV, SPHV, PHV, XCR(1,1), SGB(L), SIG
     *,              SPVH, SLVH, SSIG, CBV, SBV
     *,              WK1, WK2, WK3, TX1, TX2, TX3)
C
C
         call clkon
c        CALL MINMXL(LM, IJM0M0, I1P1, NSTEP, SGB, SSIG, LMIN, LMAX)
         do 1615 I=1,IJM0M0
         KPR(I) = LMIN(L) - 1
 1615    CONTINUE
C
         DO 1616 K=LMIN(L),LMAX(L)
         DO 1616 I=1,IJM1M0
         IF (SSIG(I) .GT. SGB(K))  KPR(I) = K
 1616    CONTINUE
         call clkoff(kprtim)
C
C
         CALL BILL(UST, C1A, SPVH, SLVH, SSIG
     *,            IWK, IPR, JPR, KPR, WK1, WK2, WK3
     *,            SGX, SGB, 0.5, 1.0, 1, JMM1)
         CALL BILL(VST, C2A, SPVH, SLVH, SSIG
     *,            IWK, IPR, JPR, KPR, WK1, WK2, WK3
     *,            SGX, SGB, 1.0, 0.5, 1, JMM1)
         CALL BILL(SDS, SMH, SPVH, SLVH, SSIG
     *,            IWK, IPR, JPR, KPR, WK1, WK2, WK3
     *,            SGX, SGB, 1.0, 1.0, 1, JMM1)
C
         IF (ROTATE) THEN
            DO 616 I=1,JSV1M0
            UMH(I) =  C1A(I) * CBV(I) - C2A(I) * SBV(I)
            VMH(I) =  C2A(I) * CBV(I) + C1A(I) * SBV(I)
  616       CONTINUE
C
            DO 1618 I=1,JVP0M0
            UMH(I+JSV1M0) = C1A(I+JSV1M0) / COS(SPVH(I+JSV1M0))
            VMH(I+JSV1M0) = C2A(I+JSV1M0)
 1618       CONTINUE
C
            DO 620 I=1,JNV2M0
            UMH(I+JNV*IMP3) = C1A(I+JNV*IMP3) * CBV(I+JNV*IMP3)
     *                      - C2A(I+JNV*IMP3) * SBV(I+JNV*IMP3)
            VMH(I+JNV*IMP3) = C2A(I+JNV*IMP3) * CBV(I+JNV*IMP3)
     *                      + C1A(I+JNV*IMP3) * SBV(I+JNV*IMP3)
  620       CONTINUE
C
            DO 622 I=1,IJM1M0
            UMH(I) =  UMH(I) * DTOA
            VMH(I) =  VMH(I) * DTOA
            SMH(I) =  SMH(I) * DT
  622       CONTINUE
         ELSE
            DO 1619 I=1,IJM1M0
            UMH(I) = C1A(I) / COS(SPVH(I)) * DTOA
            VMH(I) = C2A(I) * DTOA
            SMH(I) = SMH(I) * DT
 1619       CONTINUE
         ENDIF
C
         CALL FILPV(UMH)
         CALL FILPV(VMH)
         CALL FILPV(SMH)
         CALL WRAP(1,IJM1M0,UMH)
         CALL WRAP(1,IJM1M0,VMH)
         CALL WRAP(1,IJM1M0,SMH)
C
C
C  INTERPOLATE TO GET YU*, YV*, NLAMBDA*, NTHETA*  (V-PTS)
C
      CALL ZANGV(UMH, VMH, SMH, ONE
     *,          CPHV, SPHV, PHV, XCR(1,1), SGB(L), SIG
     *,          SPHI,   SLAM, SSIG, WK1, WK2)
C
C
         call clkon
c        CALL MINMXL(LM, IJM0M0, I1P1, NSTEP, SGB, SSIG, LMIN, LMAX)
         DO 776 I=1,IJM0M0
         KPR(I) = LMIN(L) - 1
  776    CONTINUE
C
         DO 779 K=LMIN(L),LMAX(L)
         DO 779 I=1,IJM1M0
         IF (SSIG(I) .GT. SGB(K))  KPR(I) = K
  779    CONTINUE
         call clkoff(kprtim)
C
C
         CALL BIC1nu(YVVPT, YVVPTS, SPHI, SLAM, SSIG
     *,            IWK, IPR, JPR, KPR
     *,            C1A, C2A, C3A, C4A, C1B, C2B, C3B, C4B
     *,            WK1, WK2, WK3, WK4
     *,            SGX,  SGB, S1A, 1.0, 0.5, 1, JMM1)
C
         CALL BIC1nu(YUUPT, YUVPTS, SPHI, SLAM, SSIG
     *,            IWK, IPR, JPR, KPR
     *,            C1A, C2A, C3A, C4A, C1B, C2B, C3B, C4B
     *,            WK1, WK2, WK3, WK4
     *,            SGX,  SGB, S1A, 0.5, 1.0, 1, JMM1)
C
C CALCULATE AVVPT(I,J) AT U-POINTS USING THE FULLY CENTERED APPROACH
C
      DO 9601 J=1,JMM1
      II = (J-1)*IMP3
      DO 9601 I=1,IMP3
      WK1(II+I) = BLMIK(I,2) - SLVH(II+I)
      WK2(II+I) = BLMIK(I,2) - SLAM(II+I)
 9601 CONTINUE
C
      DO 785 I=1,IJM1M0
      WK3(I) = SIN(SPVH(I))
      WK4(I) = COS(SPVH(I))
      WK5(I) = COS(WK1(I))
      WK6(I) = COS(WK2(I))
      WK7(I) = SIN(WK2(I))
      WK8(I) = SIN(SPHI(I))
C
      WK2(I) = CPHV(I)*WK4(I)*WK5(I) +   SPHV(I)*WK3(I)

C
      WK8(I) =  (WK8(I)*WK3(I)*WK6(I) + WK4(I)*COS(SPHI(I))*WK5(I))
      WK9(I) = - WK3(I) * WK7(I)
      AVVPT(I) = (WK9(I) * YUVPTS(I) + WK8(I) * YVVPTS(I))/ WK2(I)
     *         +  DL1DT  * VPN(I,L)
C
  785 CONTINUE
C
      CALL WRAP(1,IJM1M0,AVVPT)
      CALL FILPV(AVVPT)

c     CALL WRTTEM(77, IMP3, JM, IM, AvvPT, WK1)
C
C
C  CALCULATE A TERM IN FORCING FOR HELMHOLTZ EQUATION
C
      CALL UPTV(IJM0M0, AVVPT, WK2, TX1, TX2, TX3, IJP3M0)
C
      DO 814 I=1,IMP3
      WK7(I) = FDTHF(1,1)
c*sp10     WK7(I) = FDTHF(1,1) * 0.25
      WK0(I) = ARR1(1,1)
  814 CONTINUE
      DO 815 J=2,JMM1
      II = (J-1)*IMP3
      DO 815 I=1,IMP3
      WK7(II+I) = FDTHF(J,1) 
      WK8(II+I) = FDTHF(J,2)
      WK9(II+I) = ARR1(J,2)
      WK0(II+I) = ARR1(J,1)
  815 CONTINUE
      DO 817 I=I1P1,IJM1M0
      UN0(I,L) = (WK2(I) * WK8(I) + AUUPT(I)) * WK9(I)
  817 CONTINUE
cdec7 DO 826 I=I1P1,IJM1M0,IMP3
c     UN0(IM+I+2,L) = UN0(I+2,L)
c 826 CONTINUE


      DO 820 I=I1P1,IJM1M1
      WK3(I+1) = UN0(I+1,L) - UN0(I,L)
  820 CONTINUE
      DO 821 I=I1P1,IJM1M0,IMP3
      WK3(I)      = WK3(IM+I)
  821 CONTINUE
c
      CALL VPTU(IJM0M0, AUUPT, WK1, TX1, TX2, TX3, IJP3M0)
c
      DO 823 I=1,IJM1M0
      VN0(I,L) = (AVVPT(I) - WK1(I)*WK7(I)) * WK0(I)
  823 CONTINUE


      DO 824 I=1,IJM1M0
      WK2(I) = VN0(I,L) * CPHV(I)
  824 CONTINUE
      DO 825 I=I1P1,IJM1M0
      UNSS(I,L) = WK3(I) * CDU(I) + (WK2(I) - WK2(I-I1)) * CDV(I)
  825 CONTINUE
C
C
C CALCULATE TERM1 AT THE POLES

C
      TEMA = 0.0
      TEMB = 0.0
      DO 830 I=1,IM
      TEMA = TEMA + VN0(I,L)
      TEMB = TEMB - VN0(IJM2M0+I,L)
  830 CONTINUE
C
      TEMN     = 4.0 / (IM*DY)
      TEMS     = TEMN
      TM1NP(L) = TEMN * TEMB
      TM1SP(L) = TEMS * TEMA
C
C
C
         CALL ZANGPR(UTS, VTS, STS, HALF
     *,            CPHU, SPHU, PHU, XCR(1,1), SGB(L), SIG
     *,            SPHI, SLAM, SSIG, CBV, SBV
     *,            WK1, WK2, WK3,    WK4,  WK5, WK6, WK7)
C
C
  675    CONTINUE
C
         call clkon
c        CALL MINMXL(LM, IJM0M0, I1P1, NSTEP, SGB, SSIG, LMIN, LMAX)
         KPR(:) = LMIN(L) - 1
         DO 670 K=LMIN(L),LMAX(L)
         DO 670 I=1,IJM0M0
         IF (SSIG(I) .GT. SGB(K))  KPR(I) = K
  670    CONTINUE
         call clkoff(kprtim)
C
C
         CALL BILL(UST,  C1A, SPHI, SLAM, SSIG
     *,            IWK, IPR, JPR, KPR, WK1, WK2, WK3
     *,            SGX, SGB, 0.5, 1.0, 2, JMM1)
         CALL BILL(VST,  C2A, SPHI, SLAM, SSIG
     *,            IWK, IPR, JPR, KPR, WK1, WK2, WK3
     *,            SGX, SGB, 1.0, 0.5, 2, JMM1)
         CALL BILL(SDS,  SMH, SPHI, SLAM, SSIG
     *,            IWK, IPR, JPR, KPR, WK1, WK2, WK3
     *,            SGX, SGB, 1.0, 1.0, 2, JMM1)
C
         IF (ROTATE) THEN
            DO 2610 I=1,JSU2M0
            UMH(I1+I) = C1A(I1+I) * CBV(I1+I) - C2A(I1+I) * SBV(I1+I)
            VMH(I1+I) = C2A(I1+I) * CBV(I1+I) + C1A(I1+I) * SBV(I1+I)
 2610       CONTINUE
C
            DO 2612 II=1,JUP0M0
            I = II + JSU1M0
            UMH(I) = C1A(I) / COS(SPHI(I))
            VMH(I) = C2A(I)
 2612       CONTINUE
C
            DO 2614 I=1,JNU2M0
            UMH(I+JNU*IMP3) = C1A(I+JNU*IMP3) * CBV(I+JNU*IMP3)
     *                      - C2A(I+JNU*IMP3) * SBV(I+JNU*IMP3)
            VMH(I+JNU*IMP3) = C2A(I+JNU*IMP3) * CBV(I+JNU*IMP3)
     *                      + C1A(I+JNU*IMP3) * SBV(I+JNU*IMP3)
 2614       CONTINUE
C
            DO 2619 I=1,IJM2M0
            UMH(I1+I) = UMH(I1+I) * DTOA
            VMH(I1+I) = VMH(I1+I) * DTOA
 2619       CONTINUE
         ELSE
            DO 2615 I=I1P1,IJM1M0
            UMH(I) = C1A(I) / COS(SPHI(I)) * DTOA
            VMH(I) = C2A(I) * DTOA
 2615       CONTINUE
         ENDIF
C
      CALL WRAP(1,IJM0M0,UMH)
      CALL WRAP(1,IJM0M0,VMH)
      CALL WRAP(1,IJM0M0,SMH)
C
      DO 676 I=1,IMP3
      SMH(I) = SDS(I+I2,L)
      SMH(I+IJM1M0) = SDS(I+IJP1M0,L)
  676 CONTINUE
      DO 695 I=1,IJM0M0
      SMH(I) = SMH(I) * DT
  695 CONTINUE
C
C
C INTERPOLATE TO GET NT*, BL* AND BT* FOR CALCULATION OF H(K)
C
      DO 706 I=1,IJM0M0
      TX1(I+I2) = BL(I,L)
  706 CONTINUE
      CALL FILLP(TX1)
C
      CALL ZANGP(UMH, VMH, SMH, ONE
     *,          CPHU, SPHU, PHU, XCR(1,1), SGB(L), SIG
     *,          SPHI, SLAM, SSIG, WK1, WK2)
C
c
         call clkon
c        CALL MINMXL(LM, IJM0M0, I1P1, NSTEP, SGB, SSIG, LMIN, LMAX)
         KPR(:) = LMIN(L) - 1
         DO 748 K=LMIN(L),LMAX(L)
         DO 748 I=1,IJM0M0
         IF (SSIG(I) .GT. SGB(K))  KPR(I) = K
  748    CONTINUE
         call clkoff(kprtim)
C
         CALL BICPnu(VST(I2+1,L), BT,  TX1, QST, cbv, sbv, QN0(1,L)
     *,                                  SPHI, SLAM, SSIG
     *,            IWK, IPR, JPR, KPR
     *,            C1A, C2A, C3A, C4A, C1B, C2B, C3B, C4B
     *,            WK1, WK2, WK3, WK4, WK5, WK6, WK7, WK8
     *,            SGX, DTOA, SGB, S1A, 1.0, 1.0, 2, JMM1, UST(I2+1,L))
C
      DO 712 I=1,IJM0M0
      BTS(I,L) = CBV(I) + PPN(I,L)*DL1DT
      BLS(I,L) = SBV(I)
  712 CONTINUE
C
  600 CONTINUE
C
      DO 751 I=1,IJM0M0
      WK1(I) = 0.0
      WK2(I) = 0.0
  751 CONTINUE
C
      DO 750 L=1,LM
C
      CALL WRAP(1,IJM0M0,BLS(1,L))
      CALL WRAP(1,IJM0M0,BTS(1,L))
C
C     COMPUTE H(K) TERM (PG 118 OF NOTES)
C
      DO 752 I=1,IJM0M0
      TEM     = BLS(I,L) * DSG(L)
      PPN(I,L) = BTS(I,L) + (WK1(I) + TEM * 0.5) * GMA(L)
      WK1(I)   = WK1(I)   + TEM
  752 CONTINUE
  750 CONTINUE
C
      DO 755 L=LM,1,-1
      DO 753 I=1,IJM0M0
      TEM      =  PPN(I,L) * (TBR(L)*DSL(L))
      HK1(I,L) = PHIS(I) + (WK2(I) + WK1(I)*HFC(L) + TEM*0.5) * RGAS
      WK2(I  ) =  WK2(I) + TEM
  753 CONTINUE
      CALL WRAP(1,IJM0M0,HK1(1,L))
C
      CALL GLMNP2(IMP3, JM, HK1(1,L), DX(1,2), DY
     *,                                ARE2, DX(1,1), GBR(L))
      DO 754 I=1,IJM0M0
      HK1(I,L) = HK1(I,L) - GBR(L)
  754 CONTINUE
C
  755 CONTINUE
C
C
c     do 1498 l=1,lm
c     CALL WRTTEM(77, IMP3, JM, IM, HK1(1,L), WK1)
c1498 continue
c     do 1499 l=1,lm
c     CALL WRTTEM(77, IMP3, JM, IM, UNSS(1,L), WK1)
c1499 continue
c
c
      DO 1500 L=1,LM
C
      DO 1010  I=1,IJM0M0
      EA(I)  = 0.0
      TX1(I) = 0.0
 1010 CONTINUE
C
C COMPUTE H(HAT) AND A(HAT)
C

      DO 1030 K=1,LM
      DO 1020 I=1,IJM0M0
      TX1(I) = TX1(I) +  HK1(I,K) * EI(L,K)
 1020 CONTINUE
      DO 1025 I=I1P1,IJM1M0
      EA(I)     = EA(I)     + UNSS(I,K) * EI(L,K)
 1025 CONTINUE
 1030 CONTINUE
c
      if (l .lt. lmavrx) then
        call avrx3(tx1,ea)
      endif
c
c     CALL FILLP(TX1)
C
C     CALCULATE THE FORCING TERM FOR THE HELMHOLTZ EQUATION.
C
      CDFAC = HFCDDT*SQMU(L)
      TEMTX = (1.0 - CDFAC) * SQMU(L)
      TEMLK = HFDTLK(L)     * SQMU(L)
C
      DO 1100 I=I1P1,IJM1M0
      II = I2 + I
      WF1(I) = EA(I) * TEMLK - TX1(I) * TEMTX
 1100 CONTINUE
C
C CALCULATE FORCING AT POLES
C
      EANP  = 0.0
      EASP  = 0.0
      DO 1105 K=1,LM
      EASP = EASP + TM1SP(K) * EI(L,K)
      EANP = EANP + TM1NP(K) * EI(L,K)
 1105 CONTINUE
      EASP = EASP * TEMLK
      EANP = EANP * TEMLK
C
      TEMNP  = 0.0
      TEMSP  = 0.0
      DO 1110 I=1,IM
      TEMNP  =  TEMNP  + TX1(I+IJM1M0)
      TEMSP  =  TEMSP  + TX1(I)
 1110 CONTINUE
C
      temnp = temnp / im
      temsp = temsp / im
C
      TM2NP =  - TEMNP * TEMTX
      TM2SP =  - TEMSP * TEMTX
      DO 1115 I=1,IMP3
      WF1(I+IJM1M0) = EANP+TM2NP
      WF1(I)        = EASP+TM2SP

      TX1(I+IJM1M0) = TEMNP
      TX1(I)        = TEMSP
 1115 CONTINUE
c
c     CALL WRTTEM(77, IMP3, JM, IM, wf1, WK1)

C
      CALL WRAP(1,IJM0M0,WF1)
C
C
C    solve the elliptic equation
C
C
      DO 1120 I=1,IJM0M0
      RHS(I,1) = WF1(I)
      RHS(I,2) = 0.0
 1120 CONTINUE
      CALL CLKON
C
      IF (MGSOL) THEN
         DO 4115 J=1,JM
         II = (J-1) * IMP3
         DO 4115 I=1,IMP3
         WK1(II+I) = ARR1(J,2) * WF1(II+I)
 4115    CONTINUE
C*
         CALL GLMNP2(IMP3, JM, WK1, DX(1,2), DY, ARE2, DX(1,1), TEM)
C*
         DO 4116 J=1,JM
         II = (J-1) * IMP3
         DO 4116 I=1,IMP3
         RHS(II+I,1) = RHS(II+I,1) - TEM * (1.0 + FUSQ(J))
 4116    CONTINUE
         TEMPHI = - TEM / SQMU(L)
C
C        DEFINITIONS FOR THE MULTI-GRID SOLVER
C
         NOVLI  = 3
         NOVLJ1 = 0
         NOVLJ2 = 0
         IVERS = 2
         DELT   = DT
         REARTH = AE
         IDMPHI = IJM0M0*2
         HELMCO = ZERO
C

         PHIB = HFDTLK(L) / (DL1DT)
C
         CALL MGSOLV(IMP3,NOVLI,JM,NOVLJ1,NOVLJ2,IVERS,0,OMEGA,
     1              DELT,REARTH,PHIB,CDA,HELMCO,PHI,RHS,IDMPHI,MGERRO)
C
         IF (MGERRO .NE. 0) WRITE(6,*) 'MGERRO=',MGERRO

         PHI(:,1) = PHI(:,1) + TEMPHI
      ELSE
         IF (INIT1) THEN
            ITRG = 3*IMP2 + 1
            DO 1112 J=1,JM
            DSC1(J) = ARR1(J,2)
            DSC2(J) = 0.0
            DSC3(J) = ARR1(J,1)
            DSC4(J) = ARR1(J,2)*ARR1(J,2)*(1.0-FUSQ(J))*(DL1DT*TWOM)
     *                                   *CPH(J,2) / AE
            DSC5(J) = 0.0
 1112       CONTINUE
         ENDIF
         DO 1116 J=1,JM
         DSC6(J) = - SQMU(L)
 1116    CONTINUE
C   
         CALL ELLSOL(IMP3, IM, JM, AE, INIT1, INIT2, .TRUE.
     *,    DSC1, DSC2, DSC3, DSC4, DSC5, DSC6, RHS, PHI, WSV, UST, IX)
C
         INIT1 = .FALSE.
         INIT2 = .FALSE.
c
      ENDIF
C
      CALL CLKOFF(MGSLTIME)
C
C     uvm & vum NOW CONTAIN GHAT AND DHAT.
C
      CALL WRAP(1, IJM0M0, PHI(1,1))
      DO 1125 I=1,IJM0M0
      GHT(I,L) = PHI(I,1) + CDFAC*TX1(I)
      DHT(I,L) = (TX1(I) - GHT(I,L)) * (1.0/HFDTLK(L))
 1125 CONTINUE
 1500 CONTINUE
C
      DO 2000 L=1,LM
c
c     CALL WRTTEM(77, IMP3, JM, IM, ght(1,l), WK1)
C
C     HERE GTP  CONTAINS G AND DIV CONTAINS D
C
      DO 1510  I=1,IJM0M0
      GTP(I,L) = 0.0
      DIV(I,L) = 0.0
 1510 CONTINUE
      DO 1530 K=1,LM
      DO 1520 I=1,IJM0M0
      GTP(I,L) = GTP(I,L) + GHT(I,K) * E(L,K)
      DIV(I,L) = DIV(I,L) + DHT(I,K) * E(L,K)
 1520 CONTINUE
 1530 CONTINUE
C
C     OBTAIN U(T+DT)
C
      DO 1535 I=1,IJM0M0
      WK0(I) =   GTP(I,L) -   DIV(I,L)*CDA
 1535 CONTINUE
      DO 1545 I=1,IJM0M1
      WK3(I) =   WK0(I+1) +   WK0(I)
 1545 CONTINUE
      DO 1550 I=I1P1,IJM1M1
      WK4(I) = WK3(I1+I) - WK3(I-I1)
      WK5(I) = WK0(I+1)  - WK0(I)
 1550 CONTINUE
      DO 1565 J=1,JM
      II = (J-1) * IMP3
      DO 1560 I=1,IMP3
      WK6(II+I) = ARR3(J)
      WK7(II+I) = ARR4(J)
 1560 CONTINUE
 1565 CONTINUE
      DO 1570 I=I1P1,IJM1M0
      UN0(I,L) =  UN0(I,L)
     *         -  WK5(I) * WK6(I) * RDLM
     *         -  WK4(I) * WK7(I) * (RDPH*0.25)
 1570 CONTINUE
C
      CALL WRAP(I1P1,IJM1M0,UN0(1,L))
      CALL FILPU(UN0(1,L))
C
C     OBTAIN V(T+DT)
C
      DO 1590 I=1,IJM0M2
      WK3(I+1) =   WK0(I+2) -   WK0(I)
 1590 CONTINUE
      DO 1580 I=1,IJM0M0,IMP3
      WK3(IMP2+I) =   WK3(I+2)
      WK3(I)      =   WK3(IM+I)
 1580 CONTINUE
C
      DO 1598 J=1,JM
      II = (J-1)*IMP3
      DO 1597 I=1,IMP3
      WK6(II+I) = ARR5(J)
      WK7(II+I) = ARR6(J)
 1597 CONTINUE
 1598 CONTINUE
      DO 1600 I=1,IJM1M0
      VN0(I,L) =  VN0(I,L)
     *         + (WK3(I)    + WK3(I1+I)) * WK6(I) * (RDLM*0.25)
     *         - (WK0(I1+I) - WK0(I))    * WK7(I) * RDPH
 1600 CONTINUE
C
      CALL WRAP(1,IJM1M0,VN0(1,L))
C
 2000 CONTINUE
c     do 2001 l=1,lm
c     CALL WRTTEM(77, IMP3, JM, IM, div(1,l), WK1)
c2001 CONTINUE
C
C     OBTAIN SIGDOT(T+DT)
C
      DO 2410 I=1,IJM0M0
      UPN(I,1) = (BLS(I,1) - DIV(I,1)*DL1DT) * DSG(1)
 2410 CONTINUE
      DO 2421 L=2,LM
      DO 2420 I=1,IJM0M0
      UPN(I,L) = UPN(I,L-1) + (BLS(I,L) - DIV(I,L)*DL1DT) * DSG(L)
 2420 CONTINUE
 2421 CONTINUE
C
      DO 2505 I=1,IJM0M0
      LGP(I)    = UPN(I,LM) * (1.0/(1.0-SIG(1)))
      SGD(I,LM) = 0.0
      WK1(I)    = LGP(I) * RKAP
 2505 CONTINUE
C
C     OBTAIN SIGMADOT AND T(T+DT)
C
      TEM = 1.0 / DL1DT
      DO 2522 L=LMM1,1,-1
      LP1 = L + 1
      DO 2520 I=1,IJM0M0
      SGD(I,L) = (UPN(I,L) - LGP(I) * (SIG(LP1)-SIG(1))) * TEM
C
      TN0(I,LP1) = (SGD(I,L) + SGD(I,LP1)) * (DL1DT*GMA(LP1)*0.5)
     *                       + WK1(I) + BTS(I,LP1)
 2520 CONTINUE
 2522 CONTINUE
      DO 2570 I=1,IJM0M0
      TN0(I,1) = SGD(I,1) * (DL1DT*GMA(1)*0.5) + WK1(I) + BTS(I,1)
 2570 CONTINUE
C
c     CALL GETGTP(IJM0M0, LM, DSL, TBR, RGAS, PHIS, TN0
c    *,                                 LGP,  GHT,  WK1)
C
      IF (SHPFLT) THEN
         DO 2650 L=1,LM
         CALL FILTML(UN0(1,L), WK1, .TRUE., WK3, WK4)
         CALL FILTML(TN0(1,L), WK2, .TRUE., WK3, WK4)
         DO 2574 I=I1P1,IJM1M0
         UN0(I,L) = UN0(I,L) + WK1(I) * DT
         TN0(I,L) = TN0(I,L) + WK2(I) * DT
 2574    CONTINUE
         CALL FILTML(VN0(1,L), WK1, .FALSE., WK3, WK4)
         DO 2575 I=1,IJM1M0
         VN0(I,L) = VN0(I,L) + WK1(I) * DT
 2575    CONTINUE
 2650    CONTINUE
      ENDIF
      DO 2655 L=1,LM
      CALL WRAP(1,IJM0M0,UN0(1,L))
      CALL WRAP(1,IJM0M0,VN0(1,L))
      CALL WRAP(1,IJM0M0,TN0(1,L))
      CALL WRAP(1,IJM0M0,QN0(1,L))
 2655 CONTINUE
      CALL WRAP(1,IJM0M0,LGP)
C
C
      DO 4004 I=1,IJM0M0
      PIB(I) = (0.01 * EXP(LGP(I)+PBR))
 4004 CONTINUE
c
      IF (ADMASS) THEN
         CALL ADJMAS(IMP3, JM, PIB, LGP, DX, DY, ARE2, PBR, SRFPRS
     *,                                                     PCHANG)
      ENDIF
c
      DO 4005 I=1,IJM0M0
      PIK(I) = PIB(I) ** RKAP
      PKI(I) = 1.0 / PIK(I)
 4005 CONTINUE
      DO 4007 L=1,LM
      DO 4007 I=1,IJM0M0
      TN0(I,L) = (1.0 + TN0(I,L)) * TBR(L) * PKI(I) * (1.0 / PRH(L))
 4007 CONTINUE
c
      nblst = DT / 450 + 0.001
      bdt   = dt / float(nblst)
c
      IF (ADIBAT) THEN
         IF (DT .GT. 0) THEN
            do 4008 nbc=1,nblst
            CALL DCADJ2(IJM0M0, IJM0M0, LM, CMP, CPM, AMP, APM
     *,                        QMP, QPM, TN0, QN0, WK5, WK6, BDT)
 4008       continue
         endif
         RETURN
      ENDIF
C
C
C     ************  Eulerian Physics **************
C
c
c     apply radiation increments
c
      CALL SOLAR(DIURNL,SEASNL,DAYOFYR(),SECOFDY(),WK1,SLRIJ)
C
      CALL WRAP(1,IJM0M0,SLRIJ)
C
      do 6666 nbc=1,nblst
c
      DO 4010 L=1,LM
      DO 4010 I=1,IJM0M0
      TN0(I,L) = TN0(I,L) + (FLW(I,L) + FSW(I,L) * SLRIJ(I))  
     *                                * PKI(I) * (BDT/PRH(L))
 4010 CONTINUE
c
c    Vertical diffusion of theta and q
c
      DO 4300 L=1,LM
      IF (L .EQ. 1) THEN
         DO 4301 I=1,IJM0M0
         TEM    = PIK(I) * (TN0(I,L) + TN0(I,L+1))
         TEM    = VDQ(L+1) / (TEM*TEM)
         WKT(I) = TEM * (TN0(I,L+1) - TN0(I,L)) * PRJ(L+1)
         WKQ(I) = TEM * (QN0(I,L+1) - QN0(I,L))
C
         TN0(I,L) = TN0(I,L) + WKT(I) * ODP(L) * BDT
         QN0(I,L) = QN0(I,L) + WKQ(I) * ODS(L) * BDT
 4301    CONTINUE
      ELSEIF (L .LT. LM) THEN
         DO 4302 I=1,IJM0M0
         TEM    = PIK(I) * (TN0(I,L) + TN0(I,L+1))
         TEM    = VDQ(L+1) / (TEM*TEM)
         TEMD   = TEM * (TN0(I,L+1) - TN0(I,L)) * PRJ(L+1)
         TEMQ   = TEM * (QN0(I,L+1) - QN0(I,L))
C
         TN0(I,L) = TN0(I,L) + (TEMD - WKT(I)) * ODP(L) * BDT
         QN0(I,L) = QN0(I,L) + (TEMQ - WKQ(I)) * ODS(L) * BDT
         WKT(I) = TEMD
         WKQ(I) = TEMQ
 4302    CONTINUE
      ELSE
         DO 4303 I=1,IJM0M0
         TN0(I,L) = TN0(I,L) - WKT(I) * ODP(LM) * BDT
         QN0(I,L) = QN0(I,L) - WKQ(I) * ODS(LM) * BDT
 4303    CONTINUE
      ENDIF
 4300 CONTINUE 
C
C
C   SURFACE FLUXES
C
      call filpu(un0(1,lm))
      call filpv(vn0(1,lm))
      DO 4011 I=1,IJM0M0
      WK1(I) = UN0(I,LM) * UN0(I,LM)
      WK2(I) = VN0(I,LM) * VN0(I,LM)
 4011 CONTINUE
      DO 4015 i=I1P1,IJM1M0
      EKB(I+1) = (WK1(I)+WK1(I+1) + WK2(I+1) + wk2(I1+I+1)) * 0.5
 4015 CONTINUE
      DO 4020 I=I1P1,IJM1M0,IMP3
      EKB(I) = EKB(I+IM)
 4020 CONTINUE
      TEM1 = 0.0
      TEM2 = 0.0
      DO 4021 I=1,IM
      TEM1 = TEM1 + EKB(I+I1)
      TEM2 = TEM2 + EKB(I+IJM2M0)
 4021 CONTINUE
      TEM1 = TEM1 / IM
      TEM2 = TEM2 / IM
      DO 4025 I=1,IMP3
      EKB(I) = TEM1
      EKB(I+IJM1M0) = TEM2
 4025 CONTINUE
C
      DO 4030 I=1,IJM0M0
      IF (EKB(I) .LT. V2MIN) EKB(I) = V2MIN
      IF (CAPIJ(I) .EQ. 0.0) TSF(I) = SSTIJ(I)
C
      WK4(I) = TN0(I,LM) - TSF(I)*PKI(I)*(ONE/PRJ(LMP1))
C
      WK5(I) = PIB(I)
      WK8(I) = TSF(I)
4030  CONTINUE
C
C
      CALL VQSAT(WK8,WK5,WK1,WK6,IJM0M0,.TRUE. )
C
      DO 2425 I=1,IJM0M0
      WK5(I) = - WK1(I) + QN0(I,LM)
C
      IF (WK5(I) .LT. ZERO) WK5(I) = WK5(I) * WETIJ(I)
C
      WK1(I)    = WK4(I) / EKB(I) * PIK(I) * HPK(LM)
      RIB(I)    = RIB(I) + WK1(I)
      STABLE(I) = WK1(I) .GT. ZERO
2425  CONTINUE

C
      
      DO 2450 I=1,IJM0M0
      IF (STABLE(I)) THEN
         WK2(I) = SQRT(BSF*WK1(I) + ONE)
         WK3(I) = ONE / WK2(I)
      ELSE
         WK1(I) = - WK1(I)
         WK2(I) = CDSIJ(I)*SQRT(WK1(I)) + ONE
         WK3(I) = ONE / WK2(I)
         WK2(I) = WK3(I)
      ENDIF
C
      WK9(I) = (TWO*BSF)   * WK1(I) * WK3(I) + ONE
      WK0(I) = (THREE*BSF) * WK1(I) * WK2(I) + ONE
C
      WK8(I) = ONE / (WK0(I) * WK9(I))
      WK2(I) = WK9(I)
2450  CONTINUE
      DO 2500 I=1,IJM0M0
      IF (STABLE(I)) THEN
         WK9(I) = WK0(I) * WK8(I)
         WK0(I) = WK2(I) * WK8(I)
      ENDIF
C
      WK9(I) = WK9(I)*CDRIJ(I)
      WK0(I) = WK0(I)*CDRIJ(I)
c     WK9(I) = CDRIJ(I)
c     WK0(I) = CDRIJ(I)
C
      WK8(I) = SQRT(EKB(I))
      CUB(I) = CUB(I) + WK9(I)
      CTB(I) = CTB(I) + WK0(I)
      VIB(I) = VIB(I) + WK8(I)
2500  CONTINUE
C
      DO 2510 I=1,IJM0M0
      WK8(I) = WK8(I) * PIB(I) * (CMB2PA/(RGAS*PRJ(LMP1)))
     *                * PKI(I) / TN0(I,LM)
C
      WK0(I) = WK0(I) * WK8(I)
      WK9(I) = WK9(I) * WK8(I)
C
C   SENSIBLE HEAT FLUX
C
      WK1(I)    = WK4(I) * WK0(I)
      WK3(I)    = WK1(I) * PIK(I) * (CP*PRJ(LMP1))
      SHF(I)    = SHF(I) - WK3(I)
      tem       = 1.0 / pib(i)
      TN0(I,LM) = TN0(I,LM) - WK1(I)*(PRJ(LMP1)*BDT/(PRH(LM)*CNV(LM)))
     *                      * tem
C
C   LATENT HEAT FLUX
C
      WK1(I) = WK5(I) * WK0(I)
      EVP(I) = EVP(I) - WK1(I)*ALHL
C
      WK3(I)    = WK3(I)    + WK1(I)*ALHL
      QN0(I,LM) = QN0(I,LM) - WK1(I)*(BDT/CNV(LM)) * tem
C
2510  CONTINUE
C
      DO 2550 I=1,IJM0M0
      WK1(I)   = SLRIJ(I) * DSW(I)
C
      DLR(I) = DLR(I) + DLW(I)
      OLR(I) = OLR(I) + OLW(I)
      DSR(I) = DSR(I) + WK1(I)
      WK4(I) = RSW(I) * SLRIJ(I)
      RSR(I) = RSR(I) + WK4(I)
      SLI(I) = SLI(I) + SLRIJ(I)
C
      WK3(I) = DLW(I) + WK3(I) + WK1(I)
      WK1(I) = TSF(I) * TSF(I)
      WK1(I) = WK1(I) * WK1(I)
      WK3(I) = WK3(I) - WK1(I)*STFBOL
C     ULR(I) = ULR(I) + VT1(I)*STFBOL
C
      WK1(I) = TSF(I) - TICE
      IF (WK1(I) .GT. ZERO) WK1(I) = ZERO
C
      WK1(I) = WK1(I) * THKIJ(I)
      WK3(I) = WK3(I) - WK1(I)
      FIC(I) = FIC(I) - WK1(I)
C
      WK3(I) = WK3(I) * CAPIJ(I)
      TSF(I) = TSF(I) + WK3(I) * BDT
cariesTSF(I) = TSF(I) + WK3(I)*(HALF*BDT)
      WK1(I) = TSF(I) - TICE
      IF (WK1(I) .GT. ZERO .AND. THKIJ(I) .GT. ZERO) THEN
         TSF(I) = TSF(I) - WK1(I)
         FIC(I) = FIC(I) - WK1(I) / (BDT*CAPIJ(I))
caries   FIC(I) = FIC(I) - WK1(I) / (HALF*BDT*CAPIJ(I))
      ENDIF
      GTM(I) = GTM(I) + TSF(I)
2550  CONTINUE
2551  CONTINUE
C
C
C    SURFACE FRICTION
C
      DO 2595 I=I1P1,IJM1M0
      WK1(I) = (WK9(I+1) + WK9(I)) * UN0(I,LM)
      UDR(I) = UDR(I) + WK1(I)*0.5
C
      TEM = 0.5 / (PIB(I+1) + PIB(I))
      UN0(I,LM) = UN0(I,LM) - WK1(I) * TEM * (BDT/CNV(LM))
2595  CONTINUE


      DO 2620 I=1,IJM1M0
      WK2(I) = (WK9(I+I1) + WK9(I)) * VN0(I,LM)
C
      VDR(I) = VDR(I) + WK2(I)*0.5
C
      TEM = 0.5 / (PIB(I+I1) + PIB(I))
      VN0(I,LM) = VN0(I,LM) - WK2(I) * TEM * (BDT/CNV(LM))
2620  CONTINUE
 4575 CONTINUE
c
C     CALL CLKON
C
C     DRY CONVECTIVE ADJUSTMENT
C
      CALL DCADJ2(IJM0M0, IJM0M0, LM, CMP, CPM, AMP, APM
     *,                        QMP, QPM, TN0, QN0, WK5, WK6, BDT)
 
C
C     PREPARATION TO CALL RAS
C
      DO 3005 I=1,IJM0M0
      PIK(I) = PIB(I) ** RKAP
      WK0(I) = ONE / PIK(I)
3005  CONTINUE
C     CALL PTOKAPA(IJM0M0, WK1, PIK)
C
      DO 3600 L=1,LM
      DO 3595 I=1,IJM0M0
      CLT(I,L) = ZERO
3595  CONTINUE
3600  CONTINUE
C
C     SET ALF AND BET
C
 
      CALL SETAB(IJM0M0, IJM0M0, LM, PRH, SGB, PIB, TN0, PIK, ALF, BET
     *,                              GAM, WK1, WK2, WK3, WK4, WK5, WK6)
C
C     LARGE-SCALE PRECIPITATION
C
      CALL LSP(IJM0M0, IJM0M0, LM,  CNV, PRH, DSG
     *,                PIB,    TN0, QN0, ALF, BET, GAM, WK0
c    *,                CLT,    DTL, PRE, WK1, WK2, BDT, SGB, PIK)
     *,                CLT,    DTL, DQL, PRE, WK1, WK2, BDT)
C
C     INCREMENT CLOUD AMOUNTS
C
      ICLDL = ICLDL + 1
      IMCA  = IMCA  + 1
      ISTEP = ISTEP + 1

      DO 3615 L=1,LM
      DO 3610 I=1,IJM0M0
      IF (CLT(I,L) .GT. ONE) CLT(I,L) = ONE
      CLD(I,L) = CLD(I,L) + CLT(I,L)
      COU(I,L) = COU(I,L) + CLT(I,L)
      CLT(I,L) = ZERO
3610  CONTINUE
3615  CONTINUE
C
C
C     RELAXED ARAKAWA-SCHUBERT
C
C
C     CALL REGULAR RAS
C
      CALL CLKON
C
      IF (.NOT. CUMFRI) THEN
            CALL RAS (IJM0M0, IJM0M0, LM,  BDT
     *,            SGO,    ODS, PRH, PRJ, FPK, HPK, SIG(2)
     *,            PIB,    TN0,    QN0, PIK, WK0
     *,            PCU,    CLT,    TCU, QCU
     *,            ALF,    BET,    GAM, HOL, ETA
     *,            WK1,    WK2,    WK3, WK4, WK5, WK6, WK7, WK8, WK9
     *,            WFN,    AKM)
      ELSE
C
C     CALL RAS CUMLUS FRICTION
C
      DO 5710 L=1,LM
      CALL FILPU(UN0(1,L))
      call wrap(1,ijm0m0,un0(1,l))
      DO 5702 I=1,IMP3
      UOC(I,L)        = 0.0
      UOC(I+IJM1M0,L) = 0.0
      VOC(I,L)        = 0.0
      VOC(I+IJM1M0,L) = 0.0
 5702 CONTINUE
      DO 5710 I=1,IJM2M0
      UOC(I1+I+1,L) = UN0(I1+I+1,L) + UN0(I1+I,L)
      VOC(I+I1,L) = VN0(I,L)        + VN0(I+I1,L)
 5710 CONTINUE
      DO 5711 L=1,LM
      DO 5711 I=I1P1,IJM1M0, IMP3
      UOC(I,L) = UOC(I+IM,L)
 5711 CONTINUE
C
c     do 5719 iras=1,6
c     CALL SETAB(IJM0M0, IJM0M0, LM, PRH, SGB, PIB, TN0, PIK, ALF, BET
c    *,                              GAM, WK1, WK2, WK3, WK4, WK5, WK6)
c
C    call to old RAS
C
c     CALL RASV (IJM0M0, IJM0M0, LM,  BDT, IRAS, NCRND
c    *,         SGB,    ODS, PRH, PRJ, FPK, HPK, SIG(2)
c    *,         PIB,    PIK, WK0, ALF, BET, GAM
c    *,         UOC,    VOC, TN0, QN0
c    *,         PCU,    CLT, TCU, QCU
c    *,         HOL,    QOL, ETA, HST, GMH
c    *,         WK1,    WK2, WK3, WK4, WK5, WK6, WK7, WK8, WK9
c    *,         WFN,    AKM)
c
C
      ISEED = CTIME() + (NBC - 1) * BDT + 0.001
C
      CALL RASV (IJM0M0, IJM0M0, LM,  BDT
     *,          ISEED, NCRND, KRMIN, KRMAX, FRAC, RASAL, .TRUE.
     *,          CP,  ALHL, GRAV, CRTMSF
C
     *,          SGB,    ODS, PRH, PRJ, FPK, HPK, SIG(2)
     *,          PIB,    PIK, WK0, ALF, BET, GAM
     *,          UOC,    VOC, TN0, QN0
     *,          PCU,    CLT, TCU, QCU
C
     *,          HOL,    QOL, ETA, HST, GMH
     *,          WK1,    WK2, WK3, WK4, WK5, WK6, WK7, WK8
     *,          WK9,    WFN, AKM, VTC, VTD, VTE, VTF, VTG
     *,          IAR,    I1R,    I2R)
C
c     call to plug-compatible RAS
c
c     CALL RASV (IJM0M0, IJM0M0, LM,  BDT
c    *,          ISEED, NCRND, KRMIN, KRMAX, FRAC, RASAL, BOTOP
c    *,          CP,  ALHL, RKAP, GRAV, P00, CRTMSF
c    *,         SGB,    ODS, PRH, PRJ, FPK, HPK, SIG(2)
c    *,         PIB,    PIK, WK0, ALF, BET, GAM
c    *,         UOC,    VOC, TN0, QN0
c    *,         PCU,    CLT, TCU, QCU
c    *,         HOL,    QOL, ETA, HST, GMH
c    *,         WK1,    WK2, WK3, WK4, WK5, WK6, WK7, WK8, WK9
c    *,         WFN,    AKM)
c
c5719 continue
C
C
      DO 5750 L=1,LM
      call wrap(1,ijm0m0,uoc(1,l))
      DO 5720 I=I1P1,IJM1M0
      WK1(I) = (UOC(I,L) - UN0(I-1,L)  - UN0(I,L)) * 0.25
      WK2(I) = (VOC(I,L) - VN0(I-I1,L) - VN0(I,L)) * 0.25
 5720 CONTINUE
      DO 5721 I=I1P1,IJM1M0, IMP3
      WK1(I) = WK1(IM+I)
 5721 CONTINUE
      DO 5722 I=1,IMP3
      WK2(I) = 0.0
      WK2(I+IJM1M0) = 0.0
 5722 CONTINUE
      
      DO 5725 I=I1P1,IJM1M0
      UN0(I,L) = UN0(I,L)  + WK1(I+1) + WK1(I)
 5725 CONTINUE
      call wrap(i1p1,ijm1m0,un0(1,l))
C
      DO 5735 I=1,IJM1M0
      VN0(I,L) = VN0(I,L) + WK2(I) + WK2(I+I1)
 5735 CONTINUE
 5750 CONTINUE
C
      ENDIF
      ICLDC = ICLDC + 1
      DO 5760 L=1,LM
      DO 5755 I=1,IJM0M0
      IF (CLT(I,L) .GT. ONE) CLT(I,L) = ONE
      CLN(I,L) = CLN(I,L) + CLT(I,L)
      COC(I,L) = COC(I,L) + CLT(I,L)
 5755 CONTINUE
 5760 CONTINUE
c
      CALL CLKOFF(CONVTIME)

c 
      CALL FILL(IJM0M0, IJM0M0, LM, PIB, QN0, CNV, DSG, FIL
     *,                             WK1, WK2, Wk3, BDT)
C
C
C
c     CALL WRTPBP(IJM0M0, LGP, DIV(1,LM), ISEL, NPT, PBR, NSTEP
c    *,                        PSAV(1,NSTEP), DSAV(1,NSTEP))
C
 6666 continue
C
      RETURN
      END
      SUBROUTINE RESTRD (KRS, SHFTR)
C
#include "slgcm.slgprm"
#include "slgcm.slgcon"
#include "slgcm.slgprg"
#include "slgcm.slgwrk"
C
C*********************************************************************
C********************* SUBROUTINE RESTRD *****************************
C********************** 11 AUGUST 1989  ******************************
C*********************************************************************
C
C
      DOUBLE PRECISION TAU
      LOGICAL SHFTR
      REAL   HEADER(50)
c
      print *,' In RESTRD: Reading unit = ',KRS
C
      READ(KRS,END=100,ERR=100) TAU, ICLDC, ICLDL, HEADER
C
      PBR = HEADER(1)
      DO 10 L=1,LM
      GBR(L) = HEADER(L+1)
   10 CONTINUE
      PRINT *,' FOUND TAU=',TAU,' PBR=',PBR
      PRINT *,' GBR=',(GBR(L),L=1,LM)
      PRINT *,' PBR=',PBR
C
      DO 20 L=1,LM
      READ (KRS) (UN1 (I,L),I=1,IJM0M0)
     *,          (VN1 (I,L),I=1,IJM0M0)
     *,          (TN1 (I,L),I=1,IJM0M0)
     *,          (SGDO(I,L),I=1,IJM0M0)
     *,          (DIVO(I,L),I=1,IJM0M0)
c    *,          (GTPO(I,L),I=1,IJM0M0)
   20 CONTINUE
      READ (KRS) (LGPO(I),I=1,IJM0M0)
C
      DO 30 L=1,LM
      READ (KRS) (UN0(I,L),I=1,IJM0M0)
     *,          (VN0(I,L),I=1,IJM0M0)
     *,          (TN0(I,L),I=1,IJM0M0)
     *,          (SGD(I,L),I=1,IJM0M0)
     *,          (DIV(I,L),I=1,IJM0M0)
c    *,          (GTP(I,L),I=1,IJM0M0)
     *,          (QN0(I,L),I=1,IJM0M0)
   30 CONTINUE
C
      READ (KRS) (LGP(I),I=1,IJM0M0)
C
      DO 40 L=1,LM
      READ (KRS) (CLD(I,L),I=1,IJM0M0)
     *,          (CLN(I,L),I=1,IJM0M0)
   40 CONTINUE
C
      READ (KRS) (TSF(I),I=1,IJM0M0)
C
      READ (KRS) (PIB(I),I=1,IJM0M0)
C
C
      IF (SHFTR) THEN
         UN1(:,:) = UN0(:,:)
         VN1(:,:) = VN0(:,:)
         TN1(:,:) = TN0(:,:)
         SGDO(:,:) = SGD(:,:)
         DIVO(:,:) = DIV(:,:)
c        GTPO(:,:) = GTP(:,:)
         LGPO(:  ) = LGP(:  )
      ENDIF
C
      PIK(:) = PIB(:) ** RKAP
C
      WRITE(6,2000) SETTIME(TAU)
      WRITE(6,3000) YEAR(), MONTH() , DAY()
      WRITE(6,4000) HOUR(), MINUTE(), SECND()

      RETURN
100   CONTINUE
      PRINT *, ' NO RESTART AVAILABLE : STOPPING IN RESTRD   '
      print *,' TAU has value ',tau
      STOP
C
2000  FORMAT(//,' ** READ RESTART:    TAU =',E12.5)
3000  FORMAT(//,'     YEAR=',I5,'  MONTH=',I5,'    DAY=',I5)
4000  FORMAT(//,'     HOUR=',I5,' MINUTE=',I5,' SECOND=',I5)

C
      END
      SUBROUTINE RESTWT(KRS)
C
C
#include "slgcm.slgprm"
#include "slgcm.slgcon"
#include "slgcm.slgprg"
#include "slgcm.slgwrk"
C
C
C*********************************************************************
C*********************** SLAG    MODEL *******************************
C********************* SUBROUTINE RESTWT *****************************
C********************** 11 AUGUST 1989  ******************************
C*********************************************************************
C
C
      DOUBLE PRECISION TAU
      REAL   HEADER(50)
C
      TAU = CTIME()
      HEADER(1) = PBR
      DO 10 L=1,LM
      HEADER(L+1) = GBR(L)
   10 CONTINUE
C
      WRITE(KRS) TAU, ICLDC, ICLDL, HEADER
C
      DO 20 L=1,LM
      WRITE (KRS) (UN1(I,L),I=1,IJM0M0)
     *,           (VN1 (I,L),I=1,IJM0M0)
     *,           (TN1 (I,L),I=1,IJM0M0)
     *,           (SGDO(I,L),I=1,IJM0M0)
     *,           (DIVO(I,L),I=1,IJM0M0)
c    *,           (GTPO(I,L),I=1,IJM0M0)
   20 CONTINUE
      WRITE (KRS) (LGPO(I),I=1,IJM0M0)
C
      DO 30 L=1,LM
      WRITE (KRS) (UN0(I,L),I=1,IJM0M0)
     *,           (VN0(I,L),I=1,IJM0M0)
     *,           (TN0(I,L),I=1,IJM0M0)
     *,           (SGD(I,L),I=1,IJM0M0)
     *,           (DIV(I,L),I=1,IJM0M0)
c    *,           (GTP(I,L),I=1,IJM0M0)
     *,           (QN0(I,L),I=1,IJM0M0)
   30 CONTINUE
C
      WRITE (KRS) (LGP(I),I=1,IJM0M0)
C
      DO 40 L=1,LM
      WRITE (KRS) (CLD(I,L),I=1,IJM0M0)
     *,           (CLN(I,L),I=1,IJM0M0)
   40 CONTINUE
C
      WRITE (KRS) (TSF(I),I=1,IJM0M0)
C
      WRITE (KRS) (PIB(I),I=1,IJM0M0)
C
C
      WRITE(6,2000) CTIME()
      WRITE(6,3000) YEAR(), MONTH(), DAY()
      WRITE(6,4000) HOUR(), MINUTE(),SECND()

C
      RETURN
C
2000  FORMAT(//,' ** WROTE RESTART.    TAU =',E12.5)
3000  FORMAT(//,'     YEAR=',I5,'  MONTH=',I5,'    DAY=',I5)
4000  FORMAT(//,'     HOUR=',I5,' MINUTE=',I5,' SECOND=',I5)

C
      END
      SUBROUTINE MINMXL(LM, IJM0M0, I1P1, NSTEP, SGB, SSIG, LMIN, LMAX)
C
      DIMENSION SSIG(IJM0M0,LM), LMIN(LM), LMAX(LM), SGB(LM)
C
      IJM1M0 = IJM0M0 - I1P1 + 1
      IJM2M0 = IJM0M0 - I1P1 + 1 - I1P1 + 1
      DO 10 L=1,LM
      LMIN(L) = 1
      LMAX(L) = LM
  10  CONTINUE
c
      DO 50 L=1,LM
      CALL WRAP(1,IJM0M0,SSIG(1,L))
c     SGMIN = 1.0
c     SGMAX = 0.0
c
c     DO 20 I=I1P1,IJM1M0
c     SGMIN = AMIN1(SSIG(I,L),SGMIN)
c     SGMAX = AMAX1(SSIG(I,L),SGMAX)
c  20 CONTINUE
C     PRINT *,' SGMIN=',SGMIN,' SGMAX=',SGMAX,' L=',L,' IJM2M0=',IJM2M0
      IMIN = ISMIN(IJM2M0, SSIG(I1P1,L), 1)
      IMAX = ISMAX(IJM2M0, SSIG(I1P1,L), 1)
      SGMIN = SSIG(IMIN+i1p1-1,L)
      SGMAX = SSIG(IMAX+i1p1-1,L)
c     PRINT *,' SGMIN=',SGMIN,' SGMAX=',SGMAX,' IMIN=',IMIN
c    *,       ' IMAX=',IMAX
C
      DO 25 K=1,LM
      IF (SGMIN .GT. SGB(K))  LMIN(L) = K
   25 CONTINUE
c
      DO 50 K=LM,1,-1
      IF (SGMAX .LT. SGB(K))  LMAX(L) = K
   50 CONTINUE
c
c     DO 60 L=1,LM
c     LMIN(L) = LMIN(L) - 2
c     LMAX(L) = LMAX(L) + 2
c     LMIN(L) = MAX0(LMIN(L),1)
c     LMAX(L) = MIN0(LMAX(L),LM)
c     LMIN(L) = 1
c     LMAX(L) = LM
c  60 CONTINUE
C
      RETURN
      END
      SUBROUTINE ADJMAS(IMP3, JM, PIB, LGP, DX, DY, ARE2, PBR, SRFPRS
     *,                                                        PCHANG)
c
      DIMENSION PIB(IMP3,JM), LGP(IMP3,JM), DX(JM,2)
      REAL LGP
c
c
      CALL GLMNP2(IMP3, JM, PIB, DX(1,2), DY, ARE2, DX(1,1), PIBM)
C
      TEM = SRFPRS - PIBM
      PCHANG = PCHANG + TEM
      print *,' tem=',tem,' pibm=',pibm,' srfprs=',srfprs
      DO 10 I=1,IMP3*JM
      PIB(I,1) = PIB(I,1) + TEM
      LGP(I,1) = ALOG(PIB(I,1)*100.0) - PBR
   10 CONTINUE
c
      RETURN
      END
      SUBROUTINE WRTTEM(KWT, IMP3, JM, IM, FLD, TX1)
C
      DIMENSION  FLD(IMP3,JM), TX1(IM,JM)
C
      DO 10 J=1,JM
      DO 10 I=1,IM
      TX1(I,J) = FLD(I,J)
   10 CONTINUE
C
      WRITE(KWT) TX1
C
      RETURN
      END
      SUBROUTINE GETGTP(IM, LM, DSL, TBR, RGAS, GBR, PHIS, TN0
     *,                                         LGP, GTP, WK1)
C
      DIMENSION TBR(LM), DSL(LM), GBR(LM)
      DIMENSION PHIS(IM), TN0(IM,LM), LGP(IM), GTP(IM,LM), WK1(IM)
      REAL LGP
C
C     CALCULATE GTP
C
      DO 10 I=1,IM
      WK1(I) = 0.0
  10  CONTINUE
      DO 20 L=LM,1,-1
      DO 20 I=1,IM
      TEM      = TN0(I,L) * (DSL(L)*TBR(L))
      GTP(I,L) = WK1(I  ) + TEM*0.5
      WK1(I  ) = WK1(I  ) + TEM
  20  CONTINUE
      DO 40 L=1,LM
      DO 40 I=1,IM
      GTP(I,L) = PHIS(I) + (GTP(I,L) + LGP(I)*TBR(L)) * RGAS - GBR(L)
   40 CONTINUE
C
      RETURN
      END
      SUBROUTINE WRAP(IB,IE,XX)
#include "slgcm.slgprm"
      DIMENSION XX(IJM0M0)
C
      DO 10 I=IB,IE,IMP3
      XX(IM+I)   = XX(I)
      XX(IMP1+I) = XX(I+1)
      XX(IMP2+I) = XX(I+2)
   10 CONTINUE
C
      RETURN
      END
*********************************************************************
      SUBROUTINE FIPOLE(WIN,FI,VIN,ANGIN,SPOLE)
#include "slgcm.slgprm"
      DIMENSION  WIN(IMP3,JMP3), FI(IMP3,JM)
C
      LOGICAL SPOLE
C
      GL1ST = ANGIN
      IF (SPOLE) THEN
         PH1ST = -PIO2 + VIN
         IF (PH1ST .LE. -PIO2 + DPH*0.05) THEN
             FIP = WIN(1,3)
            GO TO 101
         ENDIF
      ELSE
         PH1ST =  PIO2 - VIN
         IF (PH1ST .GE. PIO2 - DPH*0.05) THEN
            FIP = WIN(1,JMP2)
            GO TO 101
         ENDIF
      ENDIF
C
      XGL1ST = (GL1ST+DLM)*RDLM
      XPH1ST = (PH1ST+PIO2)*RDPH + 3
C
C     IF (IPR .LE. 1 .OR. IPR .GE. IMP2 .OR. JPR .LE. 1 .OR.
C    *    JPR .GE. JMP2) PRINT *,' IN FIPOLE IPR=',IPR,' JPR=', JPR
C     PRINT *,' IN FIPOLE IPR=',IPR,' JPR=', JPR
C
      IF (.NOT. TRICUB) THEN
         IPR    = XGL1ST + 0.5
         JPR    = XPH1ST + 0.5
         BTPR   = XPH1ST-JPR
         ALPR   = XGL1ST-IPR
         IF (IPR .EQ. 1) IPR = IMP1
         IPRP1  = IPR + 1
         JPRP1  = JPR + 1
         IPRM1  = IPR - 1
         JPRM1  = JPR - 1
         TEM1   = 1.0 - ALPR
         TEM2   = 1.0 + ALPR
         C1A =  - ALPR * TEM1 * 0.5
         C2A =    TEM1 * TEM2
         C3A =    ALPR * TEM2 * 0.5
C
         TEM1   = 1.0 - BTPR
         TEM2   = 1.0 + BTPR
         C1B =  - BTPR * TEM1 * 0.5
         C2B =    TEM1 * TEM2
         C3B =    BTPR * TEM2 * 0.5
C
         FIP =  C1A*(C1B*WIN(IPRM1,JPRM1 ) + C2B*WIN(IPRM1,JPR )
     *             + C3B*WIN(IPRM1,JPRP1 ))
     *       +  C2A*(C1B*WIN(IPR  ,JPRM1 ) + C2B*WIN(IPR  ,JPR )
     *             + C3B*WIN(IPR  ,JPRP1 ))
     *       +  C3A*(C1B*WIN(IPRP1,JPRM1 ) + C2B*WIN(IPRP1,JPR )
     *             + C3B*WIN(IPRP1,JPRP1 ))
      ELSE
         IPR    = XGL1ST
         JPR    = XPH1ST
         BTPR   = XPH1ST-JPR
         ALPR   = XGL1ST-IPR
         IF (IPR .EQ. 1) IPR = IMP1
         IPRP1  = IPR + 1
         JPRP1  = JPR + 1
         IPRM1  = IPR - 1
         JPRM1  = JPR - 1
         ALPRSQ = ALPR*ALPR
         BTPRSQ = BTPR*BTPR
         IPRP2  = IPR + 2
         JPRP2  = JPR + 2
         TWOMAL = (2.-ALPR) * 0.5
         TWOMBT = (2.-BTPR) * 0.5
         ONMASQ = 1. - ALPRSQ
         ONMBSQ = 1. - BTPRSQ
C
C
         C1A = ONMASQ*TWOMAL
         C2A = (ALPRSQ+ALPR)*TWOMAL
         C3A = THIRD*(ALPRSQ-ALPR)*TWOMAL
         C4A = SIXTH*ONMASQ*(-ALPR)
C
         C1B = ONMBSQ*TWOMBT
         C2B = (BTPRSQ+BTPR)*TWOMBT
         C3B = THIRD*(BTPRSQ-BTPR)*TWOMBT
         C4B = SIXTH*ONMBSQ*(-BTPR)
C
         FIP  = C1A*(C1B*WIN(IPR  ,JPR   ) + C2B*WIN(IPR  ,JPRP1 )
     *             + C3B*WIN(IPR  ,JPRM1 ) + C4B*WIN(IPR  ,JPRP2 ))
     *        + C2A*(C1B*WIN(IPRP1,JPR   ) + C2B*WIN(IPRP1,JPRP1 )
     *             + C3B*WIN(IPRP1,JPRM1 ) + C4B*WIN(IPRP1,JPRP2 ))
     *        + C3A*(C1B*WIN(IPRM1,JPR   ) + C2B*WIN(IPRM1,JPRP1 )
     *             + C3B*WIN(IPRM1,JPRM1 ) + C4B*WIN(IPRM1,JPRP2 ))
     *        + C4A*(C1B*WIN(IPRP2,JPR   ) + C2B*WIN(IPRP2,JPRP1 )
     *             + C3B*WIN(IPRP2,JPRM1 ) + C4B*WIN(IPRP2,JPRP2 ))
C
C     ELSE
c        FIP = ALPR*(BTPR*WIN(IPRP1,JPRP1) + (1.0-BTPR)*WIN(IPRP1,JPR))
c    *   + (1-ALPR)*(BTPR*WIN(IPR  ,JPRP1) + (1.0-BTPR)*WIN(IPR,  JPR))
      ENDIF
  101 CONTINUE
C
      IF (SPOLE) THEN
         DO 310 I=1,IMP3
         FI(I,1) = FIP
  310    CONTINUE
      ELSE
         DO 320 I=1,IMP3
         FI(I,JM) = FIP
  320    CONTINUE
      ENDIF
C
      RETURN
      END
C
        SUBROUTINE VNPOLE(VIN,CLAM,SLAM,VOUT,ANGOUT)
#include "slgcm.slgprm"
C
        DIMENSION  VIN(IJM0M0), CLAM(IMP3,2), SLAM(IMP3,2)
     *,            WKA(IM),     WKB(IM)
C
        SUMVCS = 0.0
        SUMVSN = 0.0
        DO 10 I=1,IM
        WKA(I) = VIN(IJM2M0+I) * CLAM(I,2)
        WKB(I) = VIN(IJM2M0+I) * SLAM(I,2)
   10   CONTINUE
        DO 20 I=1,IM
        SUMVCS    = SUMVCS + WKA(I)
        SUMVSN    = SUMVSN + WKB(I)
   20   CONTINUE
        SUMVCS    = SUMVCS * TWOOIM
        SUMVSN    = SUMVSN * TWOOIM
C
C
        IF (SUMVCS .NE. 0.0) THEN
           DEL    = ATAN2(SUMVSN,SUMVCS)
           VOUT   = SQRT(SUMVCS*SUMVCS+SUMVSN*SUMVSN)
        ELSE
           DEL    = 0.0
           VOUT   = 0.0
        ENDIF
C
C  REDEFINE LAM* SUCH THAT DLM < LAM* < (TWOPI+DLM).
C
        RATIO  = (DEL-DLM)/TWOPI
        IRATIO = RATIO
        IF (RATIO .LT. 0) IRATIO = IRATIO-1
        RESULT = (RATIO-IRATIO)*TWOPI
        ANGOUT = RESULT + DLM
C
  100   CONTINUE
C
        RETURN
        END
        SUBROUTINE VSPOLE(VIN,CLAM,SLAM,VOUT,ANGOUT)
#include "slgcm.slgprm"
C
        DIMENSION  VIN(IJM0M0), CLAM(IMP3,2), SLAM(IMP3,2)
     *,            WKA(IM),     WKB(IM)
C
        SUMVCS = 0.0
        SUMVSN = 0.0
        DO 10 I=1,IM
        WKA(I) = VIN(I) * CLAM(I,2)
        WKB(I) = VIN(I) * SLAM(I,2)
   10   CONTINUE
        DO 20 I=1,IM
        SUMVCS    = SUMVCS + WKA(I)
        SUMVSN    = SUMVSN + WKB(I)
   20   CONTINUE
        SUMVCS    = SUMVCS * TWOOIM
        SUMVSN    = SUMVSN * TWOOIM
C
        IF (SUMVCS .NE. 0.0) THEN
           DEL   = ATAN2(SUMVSN,SUMVCS) + PI
           VOUT  = SQRT(SUMVCS*SUMVCS+SUMVSN*SUMVSN)
       ELSE
           DEL    = PI
           VOUT   = 0.0
      ENDIF
C
C  REDEFINE LAM* SUCH THAT DLM < LAM* < (TWOPI+DLM).
C
        RATIO  = (DEL-DLM)/TWOPI
        IRATIO = RATIO
        IF (RATIO .LT. 0) IRATIO = IRATIO-1
        RESULT = (RATIO-IRATIO)*TWOPI
        ANGOUT = RESULT + DLM
  100   CONTINUE
C
        RETURN
        END
        SUBROUTINE FSTFOR(IM, TWOOIM, PIFAC, VIN, CLAM, SLAM, VOUT,DEL)
C
        DIMENSION  VIN(IM), CLAM(IM), SLAM(IM), WKA(IM), WKB(IM)
C
        SUMVCS = 0.0
        SUMVSN = 0.0
        DO 10 I=1,IM
        WKA(I) = VIN(I) * CLAM(I)
        WKB(I) = VIN(I) * SLAM(I)
   10   CONTINUE
        DO 20 I=1,IM
        SUMVCS    = SUMVCS + WKA(I)
        SUMVSN    = SUMVSN + WKB(I)
   20   CONTINUE
        SUMVCS    = SUMVCS * TWOOIM
        SUMVSN    = SUMVSN * TWOOIM
C
        IF (SUMVCS .NE. 0.0) THEN
           DEL   = ATAN2(SUMVSN,SUMVCS) + PIFAC
           VOUT  = SQRT(SUMVCS*SUMVCS+SUMVSN*SUMVSN)
       ELSE
           DEL    = PIFAC
           VOUT   = 0.0
      ENDIF
C
        RETURN
        END
        SUBROUTINE VPOLE (IM, DLM, TWOPI, TWOOIM, PIFAC
     *,                  VIN, CLAM, SLAM, VOUT, ANGOUT)
C
        DIMENSION  VIN(IM), CLAM(IM), SLAM(IM), WKA(IM), WKB(IM)
C
        SUMVCS = 0.0
        SUMVSN = 0.0
        DO 10 I=1,IM
        WKA(I) = VIN(I) * CLAM(I)
        WKB(I) = VIN(I) * SLAM(I)
   10   CONTINUE
        DO 20 I=1,IM
        SUMVCS    = SUMVCS + WKA(I)
        SUMVSN    = SUMVSN + WKB(I)
   20   CONTINUE
        SUMVCS    = SUMVCS * TWOOIM
        SUMVSN    = SUMVSN * TWOOIM
C
C
        IF (SUMVCS .NE. 0.0) THEN
           DEL   = ATAN2(SUMVSN,SUMVCS) + PIFAC
       ELSEIF (SUMVSN .GT. 0.0) THEN
           DEL    = TWOPI * 0.25 + PIFAC
       ELSE
           DEL    = - TWOPI * 0.25 + PIFAC
       ENDIF
       VOUT  = SQRT(SUMVCS*SUMVCS+SUMVSN*SUMVSN)
C
C  REDEFINE LAM* SUCH THAT DLM < LAM* < (TWOPI+DLM).
C
        RATIO  = (DEL-DLM)/TWOPI
        IRATIO = RATIO
        IF (RATIO .LT. 0) IRATIO = IRATIO-1
        RESULT = (RATIO-IRATIO)*TWOPI
        ANGOUT = RESULT + DLM
  100   CONTINUE
C
        RETURN
        END
        SUBROUTINE UPOLE (IM, DLM, TWOPI, TWOOIM, PIFAC
     *,                  UIN, CLAM, SLAM, VOUT, ANGOUT)
C
        DIMENSION  UIN(IM), CLAM(IM), SLAM(IM), WKA(IM), WKB(IM)
C
        SUMVCS = 0.0
        SUMVSN = 0.0
        DO 10 I=1,IM
        WKA(I) = UIN(I) * CLAM(I)
        WKB(I) = UIN(I) * SLAM(I)
   10   CONTINUE
        DO 20 I=1,IM
        SUMVCS    = SUMVCS + WKA(I)
        SUMVSN    = SUMVSN + WKB(I)
   20   CONTINUE
        SUMVCS    = SUMVCS * TWOOIM
        SUMVSN    = SUMVSN * TWOOIM
C
        IF (SUMVSN .NE. 0.0) THEN
           DEL   = -ATAN2(SUMVCS,SUMVSN) + PIFAC
       ELSE
           DEL    = - TWOPI * 0.25 + PIFAC
       ENDIF
       VOUT  = SQRT(SUMVCS*SUMVCS+SUMVSN*SUMVSN)
C
C  REDEFINE LAM* SUCH THAT DLM < LAM* < (TWOPI+DLM).
C
        RATIO  = (DEL-DLM)/TWOPI
        IRATIO = RATIO
        IF (RATIO .LT. 0) IRATIO = IRATIO-1
        RESULT = (RATIO-IRATIO)*TWOPI
        ANGOUT = RESULT + DLM
  100   CONTINUE
C
        RETURN
        END
        SUBROUTINE FILPV(VFILL)
C
C*********************************************************************
C  FILL IN THE 'GHOST' VALUES AT V - POINTS J=JM
C*********************************************************************
C
C
#include "slgcm.slgprm"
C
        DIMENSION VFILL(IJM0M0)
C
        DO 10 I=1,IMB2
        VFILL(IJM1M0+I)      = - VFILL(IJM2M0+I+IMB2)
        VFILL(IJM1M0+I+IMB2) = - VFILL(IJM2M0+I)
   10   CONTINUE
C
        VFILL(IJM1M0+IMP1) = VFILL(IJM1M0+1)
        VFILL(IJM1M0+IMP2) = VFILL(IJM1M0+2)
        VFILL(IJM1M0+IMP3) = VFILL(IJM1M0+3)
C
        RETURN
        END
      SUBROUTINE FILPU(UFILL)
C
C*********************************************************************
C  FILL IN THE 'GHOST' VALUES AT U - POINTS J=1 AND J=JM
C*********************************************************************
C
#include "slgcm.slgprm"
#include "slgcm.slgclm"
C
      DIMENSION UFILL(IJM0M0)
C
c     DO 10 I=1,IMB2
c     UFILL(I)      = (UFILL(I1+I)      - UFILL(I1+I+IMB2)) * 0.5
c     UFILL(I+IMB2) = -UFILL(I)
c     UFILL(IJM1M0+I) = (UFILL(IJM2M0+I)-UFILL(IJM2M0+I+IMB2)) * 0.5
c     UFILL(IJM1M0+I+IMB2) = -UFILL(IJM1M0+I)
c  10 CONTINUE
C
C
      annp = 0.0
      bnnp = 0.0
      ansp = 0.0
      bnsp = 0.0
      DO 20 I=1,IM
      annp = annp + ufill(ijm2m0+i) * clm(i,1)
      bnnp = bnnp + ufill(ijm2m0+i) * slm(i,1)
      ansp = ansp + ufill(i+i1)     * clm(i,1)
      bnsp = bnsp + ufill(i+i1)     * slm(i,1)
 20   CONTINUE
      annp = annp * twooim
      bnnp = bnnp * twooim
      ansp = ansp * twooim
      bnsp = bnsp * twooim
      do 30 i=1,im
      ufill(i)        = ansp * clm(i,1) + bnsp * slm(i,1)
      ufill(i+ijm1m0) = annp * clm(i,1) + bnnp * slm(i,1)
  30  continue
C
      UFILL(IMP1) = UFILL(1)
      UFILL(IMP2) = UFILL(2)
      UFILL(IMP3) = UFILL(3)
      UFILL(IMP1+IJM1M0) = UFILL(1+IJM1M0)
      UFILL(IMP2+IJM1M0) = UFILL(2+IJM1M0)
      UFILL(IMP3+IJM1M0) = UFILL(3+IJM1M0)
C

      RETURN
      END
      SUBROUTINE FILPUN(UFILL, VFS, VFN)
C
C  FILL IN THE 'POLE' VALUES AT THE U-POINT USING FIRST FOURIER COMP.
C
#include "slgcm.slgprm"
#include "slgcm.slgclm"
C
      DIMENSION UFILL(IJM0M0), VFS(IM), VFN(IM)
C
      CALL FSTFOR (IM, TWOOIM, ZERO
     *,           VFN, CLM(1,2), SLM(1,2), CVVN, GLM00N)
      CALL FSTFOR (IM, TWOOIM, PI
     *,           VFS, CLM(1,2), SLM(1,2), CVVS, GLM00S)
C
      DO 10 I=1,IM
      UFILL(I)        = CVVS * SIN((I-1)*DLM-GLM00S)
      UFILL(I+IJM1M0) = CVVN * SIN((I-1)*DLM-GLM00N)
   10 CONTINUE
C
      UFILL(IMP1) = UFILL(1)
      UFILL(IMP2) = UFILL(2)
      UFILL(IMP3) = UFILL(3)
      UFILL(IMP1+IJM1M0) = UFILL(1+IJM1M0)
      UFILL(IMP2+IJM1M0) = UFILL(2+IJM1M0)
      UFILL(IMP3+IJM1M0) = UFILL(3+IJM1M0)
C
      RETURN
      END
        SUBROUTINE FILLV(VFILL)
C
C*********************************************************************
C  FILL IN THE 'GHOST' VALUES AT V - POINTS J=1,2 and J=JMP2,JMP3
C*********************************************************************
C
#include "slgcm.slgprm"
C
        DIMENSION VFILL(IMP3,JMP3)
C
C
        DO 10 I=1,IMB2
        VFILL(I,1   ) = - VFILL(I+IMB2,4   )
        VFILL(I,2   ) = - VFILL(I+IMB2,3   )
        VFILL(I,JMP2) = - VFILL(I+IMB2,JMP1)
        VFILL(I,JMP3) = - VFILL(I+IMB2,JM  )
C
        VFILL(I+IMB2,1   ) = - VFILL(I,4   )
        VFILL(I+IMB2,2   ) = - VFILL(I,3   )
        VFILL(I+IMB2,JMP2) = - VFILL(I,JMP1)
        VFILL(I+IMB2,JMP3) = - VFILL(I,JM  )
   10   CONTINUE
C
        VFILL(IMP1,1   ) = VFILL(1,1   )
        VFILL(IMP2,1   ) = VFILL(2,1   )
        VFILL(IMP3,1   ) = VFILL(3,1   )
        VFILL(IMP1,2   ) = VFILL(1,2   )
        VFILL(IMP2,2   ) = VFILL(2,2   )
        VFILL(IMP3,2   ) = VFILL(3,2   )
        VFILL(IMP1,JMP2) = VFILL(1,JMP2)
        VFILL(IMP2,JMP2) = VFILL(2,JMP2)
        VFILL(IMP3,JMP2) = VFILL(3,JMP2)
        VFILL(IMP1,JMP3) = VFILL(1,JMP3)
        VFILL(IMP2,JMP3) = VFILL(2,JMP3)
        VFILL(IMP3,JMP3) = VFILL(3,JMP3)
C
        RETURN
        END
C
      SUBROUTINE FILLU(UFILL)
C
C*********************************************************************
C  FILL IN THE 'GHOST' VALUES AT U - POINTS J=1,2,3 and J=JMP2,JMP3
C*********************************************************************
C
#include "slgcm.slgprm"
C
      DIMENSION UFILL(IMP3,JMP3)
C
      DO 10 I=1,IMB2
      UFILL(I,   1) = -UFILL(I+IMB2,5)
      UFILL(I,   2) = -UFILL(I+IMB2,4)
      UFILL(I,JMP3) = -UFILL(I+IMB2,JMP1)
c9/17/92  UFILL(I,JMP2) = (UFILL(I,JMP1)-UFILL(I+IMB2,JMP1)) * HALF
C
      UFILL(I+IMB2,   1) = -UFILL(I,5)
      UFILL(I+IMB2,   2) = -UFILL(I,4)
      UFILL(I+IMB2,JMP3) = -UFILL(I,JMP1)
c9/17/92 UFILL(I+IMB2,JMP2) = (UFILL(I+IMB2,JMP1)-UFILL(I,JMP1)) * HALF
   10 CONTINUE
c9/17/92
c     DO 20 I=1,IM
c     UFILL(I,3) = (UFILL(I,2)+UFILL(I,4)) * HALF
c  20 CONTINUE
C
      UFILL(IMP1,1   ) = UFILL(1,1   )
      UFILL(IMP2,1   ) = UFILL(2,1   )
      UFILL(IMP3,1   ) = UFILL(3,1   )
      UFILL(IMP1,2   ) = UFILL(1,2   )
      UFILL(IMP2,2   ) = UFILL(2,2   )
      UFILL(IMP3,2   ) = UFILL(3,2   )
c9/17/92
c     UFILL(IMP1,3   ) = UFILL(1,3   )
c     UFILL(IMP2,3   ) = UFILL(2,3   )
c     UFILL(IMP3,3   ) = UFILL(3,3   )
c     UFILL(IMP1,JMP2) = UFILL(1,JMP2)
c     UFILL(IMP2,JMP2) = UFILL(2,JMP2)
c     UFILL(IMP3,JMP2) = UFILL(3,JMP2)
c
      UFILL(IMP1,JMP3) = UFILL(1,JMP3)
      UFILL(IMP2,JMP3) = UFILL(2,JMP3)
      UFILL(IMP3,JMP3) = UFILL(3,JMP3)
C
      RETURN
      END
C
      SUBROUTINE FILLP(FI)
C
C*********************************************************************
C  FILL IN THE 'GHOST' VALUES AT PHI - POINTS J=1,2 and J=JMP3
C*********************************************************************
C
#include "slgcm.slgprm"
C
      DIMENSION FI(IJP3M0)
C
      DO 10 I=1,IMB2
      FI(I)        = FI(I+IMB2+I4)
      FI(I1+I)     = FI(I+IMB2+I3)
      FI(I+IJP2M0) = FI(I+IMB2+IJM0M0)
C
      FI(I+IMB2)        = FI(I+I4)
      FI(I+IMB2+I1)     = FI(I+I3)
      FI(I+IMB2+IJP2M0) = FI(I+IJM0M0)
   10 CONTINUE
C
      FI(IMP1)     = FI(1   )
      FI(IMP2)     = FI(2   )
      FI(IMP3)     = FI(3   )
      FI(IMP1+I1   )  = FI(1+I1  )
      FI(IMP2+I1   )  = FI(2+I1  )
      FI(IMP3+I1   )  = FI(3+I1  )
      FI(IMP1+IJP2M0) = FI(1+IJP2M0)
      FI(IMP2+IJP2M0) = FI(2+IJP2M0)
      FI(IMP3+IJP2M0) = FI(3+IJP2M0)
C
      RETURN
      END

      SUBROUTINE BILL(WIN, WOU, PHII, BLMI, SGI
     *,               IWK, IPR, JPR, KPR, WK1, WKY, WKZ
     *,               SGX, SGB, FACLAM, FACPHI, JA, JB)
C
C
C  USES A TRILINEAR SCHEME TO INTERPOLATE THE FIELD WIN TO THE POINT
C  (PHIIN,BLMIN) AND PUTS THE INTERPOLATED VALUE IN WOUT.
C  ON A C-GRID, FACLAM AND FACPHI MUST BE SET AS FOLLOWS DEPENDING ON
C  WHETHER THE POINTS BEING USED ARE THE U V OR FI POINTS.
C  U POINTS: FACLAM = 0.5; FACPHI = 1.0
C  V POINTS: FACLAM = 1.0; FACPHI = 0.5
C  H POINTS: FACLAM = 1.0; FACPHI = 1.0
C
#include "slgcm.slgprm"
C
      DIMENSION  WIN(IJP3M0,LM),  WOU(IJM0M0)
     *,         PHII(IJM0M0), BLMI(IJM0M0), SGI(IJM0M0)
     *,          IWK(IJM0M0), IPR(IJM0M0),  JPR(IJM0M0),  KPR(IJM0M0)
     *,          WK1(IJM0M0), WKZ(IJM0M0),  WKY(IJM0M0)
     *,          SGB(LM), SGX(LM)
C    *,     lpr(ijm0m0,lm)
C
      JL  = JB - JA + 1
      IJL = IMP3*JL
      IJA = (JA-1)*IMP3
      IJB = IJA + 1
      IJE = IJA + IJL
C
      FACDLM = (2.-FACLAM)*DLM
C
C     DO 9 I=1,IJM0M0*lm
C     lpr(i,1) = kpr(i,1) + 1
C     if (kpr(i,1) .le. 1) then
C        kpr(i,1) = lpr(i,1)
C     elseif (kpr(i,1) .ge. lm) then
C        lpr(i,1) = kpr(i,1)
C     endif
C   9 CONTINUE
C
C
C  FIRST MAKE SURE DLM < LAM* < (TWOPI+DLM) FOR THE V AND FI POINTS
C  AND THAT 1.5*DLM < LAM* < (TWOPI+1.5*DLM) FOR THE U POINTS
C
C
      DO 120 I=ijb,ije
      WK1(I) = (BLMI(I)-FACDLM) * (ONE/TWOPI)
      IWK(I) =  WK1(I)
      IPR(I) =  IWK(I)
  120 CONTINUE
C
      DO 125 I=ijb,ije
      IF (WK1(I) .LT. 0.0) THEN
         IWK(I) = IPR(I) - 1
      ENDIF
  125 CONTINUE
C
      DO 130 I=ijb,ije
      TEM1   = (WK1(I) - IWK(I))*IM + TWO
      IPR(I) = TEM1
      WKZ(I) = TEM1   - IPR(I)
      IWK(I) = IPR(I)
C
      TEM2   = (PHII(I) + PIO2)*RDPH + FACPHI
      JPR(I) = TEM2
      WKY(I) = TEM2   - JPR(I)
C
  130 CONTINUE
      DO 135 I=ijb,ije
      IF (IWK(I) .EQ. 1)  IPR(I) = IMP1
  135 CONTINUE
C
      DO 210 I=ijb,ije
      IP = IPR(I) + (JPR(I)+1)*IMP3
C
      LP = KPR(I)
C     LP1 = lPR(I,L)
      LP1 = LP + 1
C     IF (LP .EQ. 0)  LP  = 1
C     IF (LP .EQ. LM) LP1 = LM
      IF (LP .LE. 0) then
        LP  = 1
        LP1 = 1
      endif
      IF (LP .gE. lm) then
        LP  = lm
        LP1 = lm
      endif
C
      TEMK = (SGI(I) - SGB(lp)) * SGX(lp)
C
      TEM    = 1.0 - WKY(I)
      TEM1   = TEM*WIN(IP,LP)     + WKY(I)*WIN(IP+I1,LP)
      TEM2   = TEM*WIN(IP+1,LP)   + WKY(I)*WIN(IP+1+I1,LP)
      TEM3   = TEM*WIN(IP,LP1)    + WKY(I)*WIN(IP+I1,LP1)
      TEM4   = TEM*WIN(IP+1,LP1)  + WKY(I)*WIN(IP+1+I1,LP1)
C
      TEMA   = (1.0 - WKZ(I))
C
      WOU(I) = (TEMA * TEM3 + WKZ(I) * TEM4) * TEMK
     *       +   (TEMA * TEM1 + WKZ(I) * TEM2) * (1.0 - TEMK)
C
  210 CONTINUE
C
C
      RETURN
      END
      SUBROUTINE FILTML(A,B,ULATS,C,D)
#include "slgcm.slgprm"
C
C*********************************************************************
C*********************** ARIES   MODEL *******************************
C********************* SUBROUTINE FILTML *****************************
C********************** 9 FEBRUARY 1987 ******************************
C*********************************************************************
C
C
      INTEGER     ORDER, PASSES
      PARAMETER   (ORDER=8,        TAU=1.*3600.)
      PARAMETER   (PASSES=ORDER/2, SHPFAC=(-.25)**PASSES)
C
C   ARGUMENTS
C
      DIMENSION    A(IMP3,JM), B(IMP3,JM), C(IMP3,JM), D(IMP3,JM)
      LOGICAL      ULATS
C
C   ******** ZONAL SMOOTHING **********
C
      CALL WRAP(1,IJM0M0,A)
C
      DO 20 I=1,IJM0M0
      C(I,1) = A(I,1)
20    CONTINUE
C
      DO 100 K=1,PASSES
      DO 30 I=1,IJM0M1
      B(I,1) = C(I+1,1) - C(I,1)
30    CONTINUE
      DO 40 I=3,IJM0M0,IMP3
      B(IM+I,1) = B(I,1)
40    CONTINUE
      DO 50 I=1,IJM0M1
      C(I+1,1)  = B(I+1,1) - B(I,1)
50    CONTINUE
      DO 60 I=1,IJM0M0,IMP3
      C(I,1)    = C(IM+I,1)
60    CONTINUE
100   CONTINUE
C
C   SAVE ZONALLY SMOOTH FIELD IN D
C
      DO 110 I=1,IJM0M0
      D(I,1) = A(I,1) - C(I,1)*SHPFAC
110   CONTINUE
C
C    ************** MERIDIONAL SMOOTHING *************
C    FIELDS AT U AND V LATITUDES OF THE C-GRID ARE TREATED DIFFERENTLY.
C
      DO 120 I=1,IJM0M0
      C(I,1) = D(I,1)
120   CONTINUE
C
C
      IF (ULATS) THEN
         DO 150 K=1,PASSES
         DO 140 I=1,IJM1M0
         B(I,1) = C(I,2) - C(I,1)
140      CONTINUE
         DO 145 I=1,IJM2M0
         C(I,2) = B(I,2) - B(I,1)
145      CONTINUE
         C(:,1)  = ZERO
         C(:,JM) = ZERO
150      CONTINUE
      ELSE
         DO 200 K=1,PASSES
         DO 160 I=1,IMB2
         B(I     , 1) = C(I     , 1) + C(I+IMB2, 1)
         B(I+IMB2, 1) = C(I+IMB2, 1) + C(I     , 1)
         B(I     ,JM) = C(I     ,JM) + C(I+IMB2,JM)
         B(I+IMB2,JM) = C(I+IMB2,JM) + C(I     ,JM)
  160    CONTINUE
         DO 170 I=1,IJM1M0
         B(I,2) = C(I,2) - C(I,1)
170      CONTINUE
         DO 180 I=1,IJM1M0
         C(I,1) = B(I,2) - B(I,1)
180      CONTINUE
200      CONTINUE
      ENDIF
C
C    TAU IS THE DAMPING TIME IN SECONDS.
C
      DO 210 I=1,IJM0M0
      B(I,1) = (D(I,1) - C(I,1)*SHPFAC - A(I,1)) * (ONE/TAU)
210   CONTINUE
C
      CALL WRAP(1,IJM0M0,B)
      RETURN
      END
      SUBROUTINE BICPnu(VST, W1I, W2I, W3I, W1O, W2O, W3O
     *,                 PHII, BLMI,SGI
     *,                 IWK, IPR, JPR, KPR
     *,                 C1A, C2A, C3A, C4A, C1B, C2B, C3B, C4B
     *,                 WK1, WK2, WK3, WK4, WK5, WK6, WK7, WK8
     *,                 SGX, DTOA, SGB, S1A, FACLAM,FACPHI,JA,JB,UST)
C
C
C  USES A BICUBIC SCHEME TO INTERPOLATE THE FIELD W1I TO THE POINT
C  (PHII,BLMI) AND PUTS THE INTERPOLATED VALUE IN W1O.
C  ON A C-GRID, FACLAM AND FACPHI MUST BE SET AS FOLLOWS DEPENDING ON
C  WHETHER THE POINTS BEING USED ARE THE U V OR FI POINTS.
C  U POINTS: FACLAM = 0.5; FACPHI = 1.0
C  V POINTS: FACLAM = 1.0; FACPHI = 0.5
C  H POINTS: FACLAM = 1.0; FACPHI = 1.0
C
#include "slgcm.slgprm"
#include "slgcm.slgclm"
C
      DIMENSION  W1I(IJP3M0,LM),  W2I(IJP3M0), W3I(IJP3M0,LM)
     *,          W1O(IJM0M0),     W2O(IJM0M0), W3O(IJM0M0)
C
      DIMENSION  PHII(IJM0M0),    BLMI(IJM0M0), SGI(IJM0M0)
     *,          VST(IJM0M0),     SGB(LM),      UST(IJM0M0)
C
      DIMENSION  C1A(IJM0M0), C2A(IJM0M0), C3A(IJM0M0), C4A(IJM0M0)
      DIMENSION  C1B(IJM0M0), C2B(IJM0M0), C3B(IJM0M0), C4B(IJM0M0)
     *,          c5a(ijm0m0), c5b(ijm0m0)
      DIMENSION  WK1(IJM0M0), WK2(IJM0M0), WK3(IJM0M0), WK4(IJM0M0)
     *,          WK5(IJM0M0), WK6(IJM0M0), WK7(IJM0M0), WK8(IJM0M0)
     *,          IWK(IJM0M0), IPR(IJM0M0),  JPR(IJM0M0),  KPR(IJM0M0)
     *,          CVV(2),   GLM00(2)
      DIMENSION  SGX(LM), S1A(LM,4)
      data ipx/158/
      SAVE
C
C
      FACDLM = (2.-FACLAM)*DLM

      JL  = JB - JA + 1
      IJL = IMP3*JL
      IJA = (JA-1)*IMP3
      IJB = IJA + 1
      IJE = IJA + IJL
C
C     CALCULATE SOUTH & NORTH POLE VALUES
C
      CALL VPOLE (IM, DLM, TWOPI, TWOOIM, PI
     *,       VST, CLM(1,2), SLM(1,2), CVV(1), GLM00(1))
C     CALL UPOLE (IM, DLM, TWOPI, TWOOIM, PI
C    *,       UST(I1+1,L), CLM(1,1), SLM(1,1), TEM1    , TEM2      )
C    *,       VST(I1+1,L), CLM(1,1), SLM(1,1), CVV(L,1), GLM00(L,1))
C     CVV(L,1)   = (CVV(L,1)   + TEM1) * 0.5
C     GLM00(L,1) = (GLM00(L,1) + TEM2) * 0.5
C
      CALL VPOLE (IM, DLM, TWOPI, TWOOIM, ZERO
     *,      VST(IJM2M0+1), CLM(1,2), SLM(1,2), CVV(2), GLM00(2))
C     CALL UPOLE (IM, DLM, TWOPI, TWOOIM, ZERO
C    *,      UST(IJM2M0+1,L), CLM(1,1), SLM(1,1), TEM1    , TEM2      )
C    *,      VST(IJM2M0+1,L), CLM(1,1), SLM(1,1), CVV(L,2), GLM00(L,2))
C     CVV(L,2)   = (CVV(L,2)   + TEM1) * 0.5
C     GLM00(L,2) = (GLM00(L,2) + TEM2) * 0.5

      CVV(:) = CVV(:) * DTOA
C
C
C  FIRST MAKE SURE DLM < LAM* < (TWOPI+DLM) FOR THE V AND FI POINTS
C  AND THAT 1.5*DLM < LAM* < (TWOPI+1.5*DLM) FOR THE U POINTS
C
      DO 30 I=ijb,ije
      WK1(I) = (BLMI(I)-FACDLM) * (ONE/TWOPI)
      IWK(I) =  WK1(I)
      IPR(I) =  IWK(I)
   30 CONTINUE
      DO 35 I=ijb,ije
      IF (WK1(I) .LT. 0.0) IWK(I) = IPR(I) - 1
   35 CONTINUE
      DO 40 I=ijb,ije
      TEM0   = (WK1(I) - IWK(I))*IM + TWO
      IPR(I) = TEM0
      TEM    = TEM0 - IPR(I)
      c5a(i) = tem
C
      TEM1   =  TEM * TEM
      TEM2   = (2.0 - TEM) * 0.5
C
      C1A(I) =  TEM2 * (1.0 - TEM1)
      C2A(I) = (TEM1 + TEM) * TEM2
      C3A(I) = (TEM1 - TEM) * TEM2 * THIRD
      C4A(I) = (TEM1 - 1.0) * TEM  * SIXTH
   40 CONTINUE
      DO 45 I=ijb,ije
      IF (IPR(I) .EQ. 1)  IPR(I) = IMP1
   45 CONTINUE
C
      DO 47 I=ijb,ije
      TEM0   = (PHII(I) + PIO2)*RDPH + FACPHI
      JPR(I) = TEM0

      TEM  = TEM0 - JPR(I)
      c5b(i) = tem
C
      TEM1 =  TEM * TEM
      TEM2 = (2.0 - TEM) * 0.5
C
      C1B(I) =  TEM2 * (1.0 - TEM1)
      C2B(I) = (TEM1 + TEM) * TEM2
      C3B(I) = (TEM1 - TEM) * TEM2 * THIRD
      C4B(I) = (TEM1 - 1.0) * TEM * SIXTH
   47 CONTINUE
C
C     CALCULATE POLE VALUES
C
      LP = KPR(1)
      LM1 = LP - 1
      LP1 = LP + 1
      LP2 = LP + 2

      IF (LP .LT. 2) THEN
         IF (LP .EQ. 0) LP=1
         LM1 = LP
      ENDIF
      IF (LP .GT. LMM2) THEN
         IF (LP .EQ. LM) LP1 = LM
         LP2 = LP1
      ENDIF

C
      CALL FIPOLE(W1I(1,LM1), WK1, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W1I(1,LP) , WK2, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W1I(1,LP1), WK3, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W1I(1,LP2), WK4, CVV(1), GLM00(1), .TRUE.)
C
      CALL FIPOLE(W2I, W2O, CVV(1), GLM00(1), .TRUE.)
C
      CALL FIPOLE(W3I(1,LM1), WK5, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W3I(1,LP) , WK6, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W3I(1,LP1), WK7, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W3I(1,LP2), WK8, CVV(1), GLM00(1), .TRUE.)
C
      IJP1 = IJM1M0 + 1
      LP = KPR(IJP1)
      LM1 = LP - 1
      LP1 = LP + 1
      LP2 = LP + 2
      IF (LP .LT. 2) THEN
         IF (LP .EQ. 0) LP=1
         LM1 = LP
      ENDIF
      IF (LP .GT. LMM2) THEN
         IF (LP .EQ. LM) LP1 = LM
         LP2 = LP1
      ENDIF
C
      CALL FIPOLE(W1I(1,LM1), WK1, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W1I(1,LP) , WK2, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W1I(1,LP1), WK3, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W1I(1,LP2), WK4, CVV(2), GLM00(2), .FALSE.)
C
      CALL FIPOLE(W2I, W2O, CVV(2), GLM00(2), .FALSE.)
C
      CALL FIPOLE(W3I(1,LM1), WK5, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W3I(1,LP) , WK6, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W3I(1,LP1), WK7, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W3I(1,LP2), WK8, CVV(2), GLM00(2), .FALSE.)
C
C     INTERIOR
C
      DO 70 I=ijb,ije
      IP = IPR(I) + (JPR(I)+1)*IMP3
C
      LP = KPR(I)
      LM1 = LP - 1
      LP1 = LP + 1
      LP2 = LP + 2
      IF (LP .LT. 2) THEN
         IF (LP .EQ. 0) LP=1
         LM1 = LP
      ENDIF
      IF (LP .GT. LMM2) THEN
         IF (LP .EQ. LM) LP1 = LM
         LP2 = LP1
      ENDIF
C
      WK4(I) = 1.0 - C5B(I)
C
      TEM1   = C5B(I)*W1I(IP+I1,   LM1)   + WK4(I)*W1I(IP,LM1)
      TEM2   = C5B(I)*W1I(IP+1+I1, LM1)   + WK4(I)*W1I(IP+1,LM1)
C
      WK1(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
C
      TEM1   = C1B(I)*W1I(IP,     LP)    + C2B(I)*W1I(IP+I1  ,LP)
     *       + C3B(I)*W1I(IP-I1,  LP)    + C4B(I)*W1I(IP+I2  ,LP)
      TEM2   = C1B(I)*W1I(IP+1,   LP)    + C2B(I)*W1I(IP+1+I1,LP)
     *       + C3B(I)*W1I(IP+1-I1,LP)    + C4B(I)*W1I(IP+1+I2,LP)
      TEM3   = C1B(I)*W1I(IP-1,   LP)    + C2B(I)*W1I(IP-1+I1,LP)
     *       + C3B(I)*W1I(IP-1-I1,LP)    + C4B(I)*W1I(IP-1+I2,LP)
      TEM4   = C1B(I)*W1I(IP+2   ,LP)    + C2B(I)*W1I(IP+2+I1,LP)
     *       + C3B(I)*W1I(IP+2-I1,LP)    + C4B(I)*W1I(IP+2+I2,LP)
C
      WK2(I)   = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
      TEM1   = C1B(I)*W1I(IP,     LP1)   + C2B(I)*W1I(IP+I1  ,LP1)
     *       + C3B(I)*W1I(IP-I1,  LP1)   + C4B(I)*W1I(IP+I2  ,LP1)
      TEM2   = C1B(I)*W1I(IP+1,   LP1)   + C2B(I)*W1I(IP+1+I1,LP1)
     *       + C3B(I)*W1I(IP+1-I1,LP1)   + C4B(I)*W1I(IP+1+I2,LP1)
      TEM3   = C1B(I)*W1I(IP-1,   LP1)   + C2B(I)*W1I(IP-1+I1,LP1)
     *       + C3B(I)*W1I(IP-1-I1,LP1)   + C4B(I)*W1I(IP-1+I2,LP1)
      TEM4   = C1B(I)*W1I(IP+2   ,LP1)   + C2B(I)*W1I(IP+2+I1,LP1)
     *       + C3B(I)*W1I(IP+2-I1,LP1)   + C4B(I)*W1I(IP+2+I2,LP1)
C
      WK3(I  ) = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
      TEM1   = C5B(I)*W1I(IP+I1,   LP2)   + WK4(I)*W1I(IP,LP2)
      TEM2   = C5B(I)*W1I(IP+1+I1, LP2)   + WK4(I)*W1I(IP+1,LP2)
C
      WK4(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
C
      TEM1   = C1B(I)*W2I(IP     )    + C2B(I)*W2I(IP+I1  )
     *       + C3B(I)*W2I(IP-I1  )    + C4B(I)*W2I(IP+I2  )
      TEM2   = C1B(I)*W2I(IP+1   )    + C2B(I)*W2I(IP+1+I1)
     *       + C3B(I)*W2I(IP+1-I1)    + C4B(I)*W2I(IP+1+I2)
      TEM3   = C1B(I)*W2I(IP-1   )    + C2B(I)*W2I(IP-1+I1)
     *       + C3B(I)*W2I(IP-1-I1)    + C4B(I)*W2I(IP-1+I2)
      TEM4   = C1B(I)*W2I(IP+2   )    + C2B(I)*W2I(IP+2+I1)
     *       + C3B(I)*W2I(IP+2-I1)    + C4B(I)*W2I(IP+2+I2)
C
      W2O(I) = C1A(I) * TEM1 + C2A(I) * TEM2
     *       + C3A(I) * TEM3 + C4A(I) * TEM4
C
C
      WK8(I) = 1.0 - C5B(I)
C
      TEM1   = C5B(I)*W3I(IP+I1,   LM1)   + WK8(I)*W3I(IP,LM1)
      TEM2   = C5B(I)*W3I(IP+1+I1, LM1)   + WK8(I)*W3I(IP+1,LM1)
C
      WK5(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
C
      TEM1   = C1B(I)*W3I(IP,     LP)    + C2B(I)*W3I(IP+I1  ,LP)
     *       + C3B(I)*W3I(IP-I1,  LP)    + C4B(I)*W3I(IP+I2  ,LP)
      TEM2   = C1B(I)*W3I(IP+1,   LP)    + C2B(I)*W3I(IP+1+I1,LP)
     *       + C3B(I)*W3I(IP+1-I1,LP)    + C4B(I)*W3I(IP+1+I2,LP)
      TEM3   = C1B(I)*W3I(IP-1,   LP)    + C2B(I)*W3I(IP-1+I1,LP)
     *       + C3B(I)*W3I(IP-1-I1,LP)    + C4B(I)*W3I(IP-1+I2,LP)
      TEM4   = C1B(I)*W3I(IP+2   ,LP)    + C2B(I)*W3I(IP+2+I1,LP)
     *       + C3B(I)*W3I(IP+2-I1,LP)    + C4B(I)*W3I(IP+2+I2,LP)
C
      WK6(I)   = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
      TEM1   = C1B(I)*W3I(IP,     LP1)   + C2B(I)*W3I(IP+I1  ,LP1)
     *       + C3B(I)*W3I(IP-I1,  LP1)   + C4B(I)*W3I(IP+I2  ,LP1)
      TEM2   = C1B(I)*W3I(IP+1,   LP1)   + C2B(I)*W3I(IP+1+I1,LP1)
     *       + C3B(I)*W3I(IP+1-I1,LP1)   + C4B(I)*W3I(IP+1+I2,LP1)
      TEM3   = C1B(I)*W3I(IP-1,   LP1)   + C2B(I)*W3I(IP-1+I1,LP1)
     *       + C3B(I)*W3I(IP-1-I1,LP1)   + C4B(I)*W3I(IP-1+I2,LP1)
      TEM4   = C1B(I)*W3I(IP+2   ,LP1)   + C2B(I)*W3I(IP+2+I1,LP1)
     *       + C3B(I)*W3I(IP+2-I1,LP1)   + C4B(I)*W3I(IP+2+I2,LP1)
C
      WK7(I  ) = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
      TEM1   = C5B(I)*W3I(IP+I1,   LP2)   + WK8(I)*W3I(IP,LP2)
      TEM2   = C5B(I)*W3I(IP+1+I1, LP2)   + WK8(I)*W3I(IP+1,LP2)
C
      WK8(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
   70 CONTINUE
C
C
      DO 90 I=1,IJM0M0
      LV = KPR(I)
      IF (LV .GE. 2 .AND. LV .LE. LMM2) THEN
         TEM  = SGI(I)
         TEMA = (TEM-SGB(LV+1))*(TEM-SGB(LV+2))
         TEMB = (TEM-SGB(LV-1))*(TEM-SGB(LV))
c
c        c1a(i) = tema * (sgi(i)-sgb(lv  )) * s1a(lv,1)
c        c2a(i) = tema * (sgi(i)-sgb(lv-1)) * s1a(lv,2)
c        c3a(i) = temb * (sgi(i)-sgb(lv+2)) * s1a(lv,3)
c        c4a(i) = temb * (sgi(i)-sgb(lv+1)) * s1a(lv,4)
c
c        w1o(i) = c1a(i) * wk1(i) + c2a(i) * wk2(i)
c    *          + c3a(i) * wk3(i) + c4a(i) * wk4(i)
c
c        w3o(i) = c1a(i) * wk5(i) + c2a(i) * wk6(i)
c    *          + c3a(i) * wk7(i) + c4a(i) * wk8(i)
C
         W1O(I) = TEMA * ((TEM-SGB(LV))   * S1A(LV,1) * WK1(I)
     *                  + (TEM-SGB(LV-1)) * S1A(LV,2) * WK2(I))
     *          + TEMB * ((TEM-SGB(LV+2)) * S1A(LV,3) * WK3(I)
     *                  + (TEM-SGB(LV+1)) * S1A(LV,4) * WK4(I))
C
         W3O(I) = TEMA * ((TEM-SGB(LV))   * S1A(LV,1) * WK5(I)
     *                  + (TEM-SGB(LV-1)) * S1A(LV,2) * WK6(I))
     *          + TEMB * ((TEM-SGB(LV+2)) * S1A(LV,3) * WK7(I)
     *                  + (TEM-SGB(LV+1)) * S1A(LV,4) * WK8(I))
C
      ELSEif (lv .eq. 1 .or. lv .eq. lmm1) then
         TEM = (SGI(I) - SGB(LV)) * SGX(LV)
         W1O(I) = (1.0 - TEM) * WK2(I) + TEM * WK3(I)
         W3O(I) = (1.0 - TEM) * WK6(I) + TEM * WK7(I)
      else
         W1O(I) = WK2(I)
         W3O(I) = WK6(I)
      ENDIF
   90 CONTINUE
C
  500 CONTINUE
C
      RETURN
      END
      SUBROUTINE BIC1nu(W1I, W1O, PHII, BLMI, SGI
     *,               IWK, IPR, JPR, KPR
     *,               C1A, C2A, C3A, C4A, C1B, C2B, C3B, C4B
     *,               WK1, WK2, WK3, WK4
     *,               SGX, SGB, S1A, FACLAM,FACPHI,JA,JB)
C
C
C  USES A BICUBIC SCHEME TO INTERPOLATE THE FIELD WIN TO THE POINT
C  (PHIIN,BLMIN) AND PUTS THE INTERPOLATED VALUE IN WOUT.
C  ON A C-GRID, FACLAM AND FACPHI MUST BE SET AS FOLLOWS DEPENDING ON
C  WHETHER THE POINTS BEING USED ARE THE U V OR FI POINTS.
C  U POINTS: FACLAM = 0.5; FACPHI = 1.0
C  V POINTS: FACLAM = 1.0; FACPHI = 0.5
C  H POINTS: FACLAM = 1.0; FACPHI = 1.0
C
#include "slgcm.slgprm"
C
      DIMENSION   W1I(IJP3M0,LM), W1O(IJM0M0)
     *,         PHII(IJM0M0), BLMI(IJM0M0), SGI(IJM0M0)
     *,          SGB(LM)
C
      DIMENSION  C1A(IJM0M0), C2A(IJM0M0), C3A(IJM0M0), C4A(IJM0M0)
      DIMENSION  C1B(IJM0M0), C2B(IJM0M0), C3B(IJM0M0), C4B(IJM0M0)
      dimension  c5a(ijm0m0), c5b(ijm0m0)
      DIMENSION  WK1(IJM0M0), WK2(IJM0M0), WK3(IJM0M0), WK4(IJM0M0)
     *,          IWK(IJM0M0), IPR(IJM0M0),  JPR(IJM0M0),  KPR(IJM0M0)
C
      DIMENSION  SGX(LM), S1A(LM,4)
C
C
      FACDLM = (2.-FACLAM)*DLM

      JL  = JB - JA + 1
      IJL = IMP3*JL
      IJA = (JA-1)*IMP3
      IJB = IJA + 1
      IJE = IJA + IJL
C
C
C  FIRST MAKE SURE DLM < LAM* < (TWOPI+DLM) FOR THE V AND FI POINTS
C  AND THAT 1.5*DLM < LAM* < (TWOPI+1.5*DLM) FOR THE U POINTS
C
      DO 30 I=IJB,IJE
      WK1(I) = (BLMI(I)-FACDLM) * (ONE/TWOPI)
      IWK(I) =  WK1(I)
      IPR(I) =  IWK(I)
   30 CONTINUE
      DO 35 I=ijb,ije
      IF (WK1(I) .LT. 0.0) IWK(I) = IPR(I) - 1
   35 CONTINUE
      DO 40 I=ijb,ije
      TEM0   = (WK1(I) - IWK(I))*IM + TWO
      IPR(I) = TEM0
      TEM    = TEM0 - IPR(I)
      c5a(i) = tem
C
      TEM1   =  TEM * TEM
      TEM2   = (2.0 - TEM) * 0.5
C
      C1A(I) =  TEM2 * (1.0 - TEM1)
      C2A(I) = (TEM1 + TEM) * TEM2
      C3A(I) = (TEM1 - TEM) * TEM2 * THIRD
      C4A(I) = (TEM1 - 1.0) * TEM  * SIXTH
   40 CONTINUE
      DO 45 I=ijb,ije
      IF (IPR(I) .EQ. 1)  IPR(I) = IMP1
   45 CONTINUE
C
      DO 47 I=ijb,ije
      TEM0   = (PHII(I) + PIO2)*RDPH + FACPHI
      JPR(I) = TEM0

      TEM  = TEM0 - JPR(I)
      c5b(i) = tem
C
      TEM1 =  TEM * TEM
      TEM2 = (2.0 - TEM) * 0.5
C
      C1B(I) =  TEM2 * (1.0 - TEM1)
      C2B(I) = (TEM1 + TEM) * TEM2
      C3B(I) = (TEM1 - TEM) * TEM2 * THIRD
      C4B(I) = (TEM1 - 1.0) * TEM * SIXTH
   47 CONTINUE
C
C
      DO 70 I=ijb,ije
      IP = IPR(I) + (JPR(I)+1)*IMP3
C
      LP = KPR(I)
      LM1 = LP - 1
      LP1 = LP + 1
      LP2 = LP + 2
      IF (LP .LT. 2) THEN
         IF (LP .EQ. 0) LP=1
         LM1 = LP
      ENDIF
      IF (LP .GT. LMM2) THEN
         IF (LP .EQ. LM) LP1 = LM
         LP2 = LP1
      ENDIF
C
C
      WK4(I) = 1.0 - C5B(I)
      TEM1   = WK4(I)*W1I(IP,   LM1) + C5B(I)*W1I(IP+I1  ,LM1)
      TEM2   = WK4(I)*W1I(IP+1, LM1) + C5B(I)*W1I(IP+1+I1,LM1)
C
      WK1(I) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
C
      TEM1   = C1B(I)*W1I(IP,     LP)    + C2B(I)*W1I(IP+I1  ,LP)
     *       + C3B(I)*W1I(IP-I1,  LP)    + C4B(I)*W1I(IP+I2  ,LP)
      TEM2   = C1B(I)*W1I(IP+1,   LP)    + C2B(I)*W1I(IP+1+I1,LP)
     *       + C3B(I)*W1I(IP+1-I1,LP)    + C4B(I)*W1I(IP+1+I2,LP)
      TEM3   = C1B(I)*W1I(IP-1,   LP)    + C2B(I)*W1I(IP-1+I1,LP)
     *       + C3B(I)*W1I(IP-1-I1,LP)    + C4B(I)*W1I(IP-1+I2,LP)
      TEM4   = C1B(I)*W1I(IP+2   ,LP)    + C2B(I)*W1I(IP+2+I1,LP)
     *       + C3B(I)*W1I(IP+2-I1,LP)    + C4B(I)*W1I(IP+2+I2,LP)
C
      WK2(I) = C1A(I) * TEM1 + C2A(I) * TEM2
     *       + C3A(I) * TEM3 + C4A(I) * TEM4
C
      TEM1   = C1B(I)*W1I(IP,     LP1)   + C2B(I)*W1I(IP+I1  ,LP1)
     *       + C3B(I)*W1I(IP-I1,  LP1)   + C4B(I)*W1I(IP+I2  ,LP1)
      TEM2   = C1B(I)*W1I(IP+1,   LP1)   + C2B(I)*W1I(IP+1+I1,LP1)
     *       + C3B(I)*W1I(IP+1-I1,LP1)   + C4B(I)*W1I(IP+1+I2,LP1)
      TEM3   = C1B(I)*W1I(IP-1,   LP1)   + C2B(I)*W1I(IP-1+I1,LP1)
     *       + C3B(I)*W1I(IP-1-I1,LP1)   + C4B(I)*W1I(IP-1+I2,LP1)
      TEM4   = C1B(I)*W1I(IP+2   ,LP1)   + C2B(I)*W1I(IP+2+I1,LP1)
     *       + C3B(I)*W1I(IP+2-I1,LP1)   + C4B(I)*W1I(IP+2+I2,LP1)
C
      WK3(I) = C1A(I) * TEM1 + C2A(I) * TEM2
     *       + C3A(I) * TEM3 + C4A(I) * TEM4
C
C
      TEM1   = WK4(I)*W1I(IP,   LP2) + C5B(I)*W1I(IP+I1  ,LP2)
      TEM2   = WK4(I)*W1I(IP+1, LP2) + C5B(I)*W1I(IP+1+I1,LP2)
C
      WK4(I) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
C
   70 CONTINUE
C
C
      DO 90 I=IJB,IJE
      LV = KPR(I)
      IF (LV .GE. 2 .AND. LV .LE. LMM2) THEN
         TEM  = SGI(I)
         TEMA = (TEM-SGB(LV+1))*(TEM-SGB(LV+2))
         TEMB = (TEM-SGB(LV-1))*(TEM-SGB(LV))
C
         TEM1 = (TEM-SGB(LV))   * TEMA * S1A(LV,1)
         TEM2 = (TEM-SGB(LV-1)) * TEMA * S1A(LV,2)
         TEM3 = (TEM-SGB(LV+2)) * TEMB * S1A(LV,3)
         TEM4 = (TEM-SGB(LV+1)) * TEMB * S1A(LV,4)
C
C        TEM1 = (TEM-SGB(LV))  *(TEM-SGB(LV+1))*(TEM-SGB(LV+2))*S1A(LV)
C        TEM2 = (TEM-SGB(LV-1))*(TEM-SGB(LV+1))*(TEM-SGB(LV+2))*S1B(LV)
C        TEM3 = (TEM-SGB(LV-1))*(TEM-SGB(LV))  *(TEM-SGB(LV+2))*S1C(LV)
C        TEM4 = (TEM-SGB(LV-1))*(TEM-SGB(LV))  *(TEM-SGB(LV+1))*S1D(LV)
C
         W1O(I) = TEM1*WK1(I) + TEM2*WK2(I) + TEM3*WK3(I) +TEM4*WK4(I)
      ELSEif (lv .eq. 1 .or. lv .eq. lmm1) then
         TEM = (SGI(I) - SGB(LV)) * SGX(LV)
         W1O(I) = (1.0 - TEM) * WK2(I) + TEM * WK3(I)
      else
         W1O(I) = WK2(I)
      ENDIF
   90 CONTINUE
      CALL WRAP(IJB,IJE,W1O)
  100 CONTINUE
C
      RETURN
      END
      SUBROUTINE BICPn1(VST, W1I, W2I, W3I, W1O, W2O, W3O
     *,                 PHII, BLMI,SGI
     *,                 IWK, IPR, JPR, KPR
     *,                 C1A, C2A, C3A, C4A, C1B, C2B, C3B, C4B
     *,                 WK1, WK2, WK3, WK4, WK5, WK6, WK7, WK8
     *,                 SGX, DTOA, SGB, S1A, FACLAM,FACPHI,JA,JB,UST)
C
C
C  USES A BICUBIC SCHEME TO INTERPOLATE THE FIELD W1I TO THE POINT
C  (PHII,BLMI) AND PUTS THE INTERPOLATED VALUE IN W1O.
C  ON A C-GRID, FACLAM AND FACPHI MUST BE SET AS FOLLOWS DEPENDING ON
C  WHETHER THE POINTS BEING USED ARE THE U V OR FI POINTS.
C  U POINTS: FACLAM = 0.5; FACPHI = 1.0
C  V POINTS: FACLAM = 1.0; FACPHI = 0.5
C  H POINTS: FACLAM = 1.0; FACPHI = 1.0
C
#include "slgcm.slgprm"
#include "slgcm.slgclm"
C
      DIMENSION  W1I(IJP3M0,LM),  W2I(IJP3M0), W3I(IJP3M0,LM)
     *,          W1O(IJM0M0),     W2O(IJM0M0), W3O(IJM0M0)
C
      DIMENSION  PHII(IJM0M0),    BLMI(IJM0M0), SGI(IJM0M0)
     *,          VST(IJM0M0),     SGB(LM),      UST(IJM0M0)
C
      DIMENSION  C1A(IJM0M0), C2A(IJM0M0), C3A(IJM0M0), C4A(IJM0M0)
      DIMENSION  C1B(IJM0M0), C2B(IJM0M0), C3B(IJM0M0), C4B(IJM0M0)
     *,          c5a(ijm0m0), c5b(ijm0m0)
      DIMENSION  WK1(IJM0M0), WK2(IJM0M0), WK3(IJM0M0), WK4(IJM0M0)
     *,          WK5(IJM0M0), WK6(IJM0M0), WK7(IJM0M0), WK8(IJM0M0)
     *,          IWK(IJM0M0), IPR(IJM0M0),  JPR(IJM0M0),  KPR(IJM0M0)
     *,          CVV(2),   GLM00(2)
      DIMENSION  SGX(LM), S1A(LM,4)
      data ipx/13200/
      SAVE
C
C
      FACDLM = (2.-FACLAM)*DLM

c
      jlin = 6
      ijs = (jlin+2) * imp3
      ijn = ijp3m0 - (jlin+1) * imp3
c
      JL  = JB - JA + 1
      IJL = IMP3*JL
      IJA = (JA-1)*IMP3
      IJB = IJA + 1
      IJE = IJA + IJL
C
C     CALCULATE SOUTH & NORTH POLE VALUES
C
      CALL VPOLE (IM, DLM, TWOPI, TWOOIM, PI
     *,       VST, CLM(1,2), SLM(1,2), CVV(1), GLM00(1))
C     CALL UPOLE (IM, DLM, TWOPI, TWOOIM, PI
C    *,       UST(I1+1,L), CLM(1,1), SLM(1,1), TEM1    , TEM2      )
C    *,       VST(I1+1,L), CLM(1,1), SLM(1,1), CVV(L,1), GLM00(L,1))
C     CVV(L,1)   = (CVV(L,1)   + TEM1) * 0.5
C     GLM00(L,1) = (GLM00(L,1) + TEM2) * 0.5
C
      CALL VPOLE (IM, DLM, TWOPI, TWOOIM, ZERO
     *,      VST(IJM2M0+1), CLM(1,2), SLM(1,2), CVV(2), GLM00(2))
C     CALL UPOLE (IM, DLM, TWOPI, TWOOIM, ZERO
C    *,      UST(IJM2M0+1,L), CLM(1,1), SLM(1,1), TEM1    , TEM2      )
C    *,      VST(IJM2M0+1,L), CLM(1,1), SLM(1,1), CVV(L,2), GLM00(L,2))
C     CVV(L,2)   = (CVV(L,2)   + TEM1) * 0.5
C     GLM00(L,2) = (GLM00(L,2) + TEM2) * 0.5

      CVV(:) = CVV(:) * DTOA
C
C
C  FIRST MAKE SURE DLM < LAM* < (TWOPI+DLM) FOR THE V AND FI POINTS
C  AND THAT 1.5*DLM < LAM* < (TWOPI+1.5*DLM) FOR THE U POINTS
C
      DO 30 I=ijb,ije
      WK1(I) = (BLMI(I)-FACDLM) * (ONE/TWOPI)
      IWK(I) =  WK1(I)
      IPR(I) =  IWK(I)
   30 CONTINUE
      DO 35 I=ijb,ije
      IF (WK1(I) .LT. 0.0) IWK(I) = IPR(I) - 1
   35 CONTINUE
      DO 40 I=ijb,ije
      TEM0   = (WK1(I) - IWK(I))*IM + TWO
      IPR(I) = TEM0
      TEM    = TEM0 - IPR(I)
      c5a(i) = tem
C
      TEM1   =  TEM * TEM
      TEM2   = (2.0 - TEM) * 0.5
C
      C1A(I) =  TEM2 * (1.0 - TEM1)
      C2A(I) = (TEM1 + TEM) * TEM2
      C3A(I) = (TEM1 - TEM) * TEM2 * THIRD
      C4A(I) = (TEM1 - 1.0) * TEM  * SIXTH
   40 CONTINUE
      DO 45 I=ijb,ije
      IF (IPR(I) .EQ. 1)  IPR(I) = IMP1
   45 CONTINUE
C
      DO 47 I=ijb,ije
      TEM0   = (PHII(I) + PIO2)*RDPH + FACPHI
      JPR(I) = TEM0

      TEM  = TEM0 - JPR(I)
      c5b(i) = tem
C
      TEM1 =  TEM * TEM
      TEM2 = (2.0 - TEM) * 0.5
C
      C1B(I) =  TEM2 * (1.0 - TEM1)
      C2B(I) = (TEM1 + TEM) * TEM2
      C3B(I) = (TEM1 - TEM) * TEM2 * THIRD
      C4B(I) = (TEM1 - 1.0) * TEM * SIXTH
   47 CONTINUE
C
C     CALCULATE POLE VALUES
C
      LP = KPR(1)
      LM1 = LP - 1
      LP1 = LP + 1
      LP2 = LP + 2

      IF (LP .LT. 2) THEN
         IF (LP .EQ. 0) LP=1
         LM1 = LP
      ENDIF
      IF (LP .GT. LMM2) THEN
         IF (LP .EQ. LM) LP1 = LM
         LP2 = LP1
      ENDIF

C
      CALL FIPOLE(W1I(1,LM1), WK1, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W1I(1,LP) , WK2, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W1I(1,LP1), WK3, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W1I(1,LP2), WK4, CVV(1), GLM00(1), .TRUE.)
C
      CALL FIPOLE(W2I, W2O, CVV(1), GLM00(1), .TRUE.)
C
      CALL FIPOLE(W3I(1,LM1), WK5, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W3I(1,LP) , WK6, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W3I(1,LP1), WK7, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W3I(1,LP2), WK8, CVV(1), GLM00(1), .TRUE.)
C
      IJP1 = IJM1M0 + 1
      LP = KPR(IJP1)
      LM1 = LP - 1
      LP1 = LP + 1
      LP2 = LP + 2
      IF (LP .LT. 2) THEN
         IF (LP .EQ. 0) LP=1
         LM1 = LP
      ENDIF
      IF (LP .GT. LMM2) THEN
         IF (LP .EQ. LM) LP1 = LM
         LP2 = LP1
      ENDIF
C
      CALL FIPOLE(W1I(1,LM1), WK1, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W1I(1,LP) , WK2, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W1I(1,LP1), WK3, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W1I(1,LP2), WK4, CVV(2), GLM00(2), .FALSE.)
C
      CALL FIPOLE(W2I, W2O, CVV(2), GLM00(2), .FALSE.)
C
      CALL FIPOLE(W3I(1,LM1), WK5, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W3I(1,LP) , WK6, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W3I(1,LP1), WK7, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W3I(1,LP2), WK8, CVV(2), GLM00(2), .FALSE.)
C
C     INTERIOR
C
      DO 70 I=ijb,ije
      IP = IPR(I) + (JPR(I)+1)*IMP3
C
      LP = KPR(I)
      LM1 = LP - 1
      LP1 = LP + 1
      LP2 = LP + 2
      IF (LP .LT. 2) THEN
         IF (LP .EQ. 0) LP=1
         LM1 = LP
      ENDIF
      IF (LP .GT. LMM2) THEN
         IF (LP .EQ. LM) LP1 = LM
         LP2 = LP1
      ENDIF
C
      WK4(I) = 1.0 - C5B(I)
C
      TEM1   = C5B(I)*W1I(IP+I1,   LM1)   + WK4(I)*W1I(IP,LM1)
      TEM2   = C5B(I)*W1I(IP+1+I1, LM1)   + WK4(I)*W1I(IP+1,LM1)
C
      WK1(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
      if (ip .gt. ijs .and. ip .lt. ijn) then
C
      TEM1   = C1B(I)*W1I(IP,     LP)    + C2B(I)*W1I(IP+I1  ,LP)
     *       + C3B(I)*W1I(IP-I1,  LP)    + C4B(I)*W1I(IP+I2  ,LP)
      TEM2   = C1B(I)*W1I(IP+1,   LP)    + C2B(I)*W1I(IP+1+I1,LP)
     *       + C3B(I)*W1I(IP+1-I1,LP)    + C4B(I)*W1I(IP+1+I2,LP)
      TEM3   = C1B(I)*W1I(IP-1,   LP)    + C2B(I)*W1I(IP-1+I1,LP)
     *       + C3B(I)*W1I(IP-1-I1,LP)    + C4B(I)*W1I(IP-1+I2,LP)
      TEM4   = C1B(I)*W1I(IP+2   ,LP)    + C2B(I)*W1I(IP+2+I1,LP)
     *       + C3B(I)*W1I(IP+2-I1,LP)    + C4B(I)*W1I(IP+2+I2,LP)
C
      WK2(I)   = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
      TEM1   = C1B(I)*W1I(IP,     LP1)   + C2B(I)*W1I(IP+I1  ,LP1)
     *       + C3B(I)*W1I(IP-I1,  LP1)   + C4B(I)*W1I(IP+I2  ,LP1)
      TEM2   = C1B(I)*W1I(IP+1,   LP1)   + C2B(I)*W1I(IP+1+I1,LP1)
     *       + C3B(I)*W1I(IP+1-I1,LP1)   + C4B(I)*W1I(IP+1+I2,LP1)
      TEM3   = C1B(I)*W1I(IP-1,   LP1)   + C2B(I)*W1I(IP-1+I1,LP1)
     *       + C3B(I)*W1I(IP-1-I1,LP1)   + C4B(I)*W1I(IP-1+I2,LP1)
      TEM4   = C1B(I)*W1I(IP+2   ,LP1)   + C2B(I)*W1I(IP+2+I1,LP1)
     *       + C3B(I)*W1I(IP+2-I1,LP1)   + C4B(I)*W1I(IP+2+I2,LP1)
C
      WK3(I  ) = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
      else
C
      TEM1   = C1B(I)*W1I(IP,     LP)   + C2B(I)*W1I(IP+I1  ,LP)
     *       + C3B(I)*W1I(IP-I1,  LP)   + C4B(I)*W1I(IP+I2  ,LP)
      TEM2   = C1B(I)*W1I(IP+1,   LP)   + C2B(I)*W1I(IP+1+I1,LP)
     *       + C3B(I)*W1I(IP+1-I1,LP)   + C4B(I)*W1I(IP+1+I2,LP)
C
      WK2(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
      TEM1   = C1B(I)*W1I(IP,     LP1)   + C2B(I)*W1I(IP+I1  ,LP1)
     *       + C3B(I)*W1I(IP-I1,  LP1)   + C4B(I)*W1I(IP+I2  ,LP1)
      TEM2   = C1B(I)*W1I(IP+1,   LP1)   + C2B(I)*W1I(IP+1+I1,LP1)
     *       + C3B(I)*W1I(IP+1-I1,LP1)   + C4B(I)*W1I(IP+1+I2,LP1)
C
      WK3(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
      endif
C
      TEM1   = C5B(I)*W1I(IP+I1,   LP2)   + WK4(I)*W1I(IP,LP2)
      TEM2   = C5B(I)*W1I(IP+1+I1, LP2)   + WK4(I)*W1I(IP+1,LP2)
C
      WK4(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
C
      TEM1   = C1B(I)*W2I(IP     )    + C2B(I)*W2I(IP+I1  )
     *       + C3B(I)*W2I(IP-I1  )    + C4B(I)*W2I(IP+I2  )
      TEM2   = C1B(I)*W2I(IP+1   )    + C2B(I)*W2I(IP+1+I1)
     *       + C3B(I)*W2I(IP+1-I1)    + C4B(I)*W2I(IP+1+I2)
      TEM3   = C1B(I)*W2I(IP-1   )    + C2B(I)*W2I(IP-1+I1)
     *       + C3B(I)*W2I(IP-1-I1)    + C4B(I)*W2I(IP-1+I2)
      TEM4   = C1B(I)*W2I(IP+2   )    + C2B(I)*W2I(IP+2+I1)
     *       + C3B(I)*W2I(IP+2-I1)    + C4B(I)*W2I(IP+2+I2)
C
      if (ip .gt. ijs .and. ip .lt. ijn) then

      W2O(I) = C1A(I) * TEM1 + C2A(I) * TEM2
     *       + C3A(I) * TEM3 + C4A(I) * TEM4
      else

      W2O(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1

      endif
C
C
      WK8(I) = 1.0 - C5B(I)
C
      TEM1   = C5B(I)*W3I(IP+I1,   LM1)   + WK8(I)*W3I(IP,LM1)
      TEM2   = C5B(I)*W3I(IP+1+I1, LM1)   + WK8(I)*W3I(IP+1,LM1)
C
      WK5(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
C
      if (ip .gt. ijs .and. ip .lt. ijn) then
C
      TEM1   = C1B(I)*W3I(IP,     LP)    + C2B(I)*W3I(IP+I1  ,LP)
     *       + C3B(I)*W3I(IP-I1,  LP)    + C4B(I)*W3I(IP+I2  ,LP)
      TEM2   = C1B(I)*W3I(IP+1,   LP)    + C2B(I)*W3I(IP+1+I1,LP)
     *       + C3B(I)*W3I(IP+1-I1,LP)    + C4B(I)*W3I(IP+1+I2,LP)
      TEM3   = C1B(I)*W3I(IP-1,   LP)    + C2B(I)*W3I(IP-1+I1,LP)
     *       + C3B(I)*W3I(IP-1-I1,LP)    + C4B(I)*W3I(IP-1+I2,LP)
      TEM4   = C1B(I)*W3I(IP+2   ,LP)    + C2B(I)*W3I(IP+2+I1,LP)
     *       + C3B(I)*W3I(IP+2-I1,LP)    + C4B(I)*W3I(IP+2+I2,LP)
C
      WK6(I)   = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
      TEM1   = C1B(I)*W3I(IP,     LP1)   + C2B(I)*W3I(IP+I1  ,LP1)
     *       + C3B(I)*W3I(IP-I1,  LP1)   + C4B(I)*W3I(IP+I2  ,LP1)
      TEM2   = C1B(I)*W3I(IP+1,   LP1)   + C2B(I)*W3I(IP+1+I1,LP1)
     *       + C3B(I)*W3I(IP+1-I1,LP1)   + C4B(I)*W3I(IP+1+I2,LP1)
      TEM3   = C1B(I)*W3I(IP-1,   LP1)   + C2B(I)*W3I(IP-1+I1,LP1)
     *       + C3B(I)*W3I(IP-1-I1,LP1)   + C4B(I)*W3I(IP-1+I2,LP1)
      TEM4   = C1B(I)*W3I(IP+2   ,LP1)   + C2B(I)*W3I(IP+2+I1,LP1)
     *       + C3B(I)*W3I(IP+2-I1,LP1)   + C4B(I)*W3I(IP+2+I2,LP1)
C
      WK7(I  ) = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C 
      else
C 
      TEM1   = C1B(I)*W3I(IP,     LP)   + C2B(I)*W3I(IP+I1  ,LP)
     *       + C3B(I)*W3I(IP-I1,  LP)   + C4B(I)*W3I(IP+I2  ,LP)
      TEM2   = C1B(I)*W3I(IP+1,   LP)   + C2B(I)*W3I(IP+1+I1,LP)
     *       + C3B(I)*W3I(IP+1-I1,LP)   + C4B(I)*W3I(IP+1+I2,LP)
C    
      WK6(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C    
      TEM1   = C1B(I)*W3I(IP,     LP1)   + C2B(I)*W3I(IP+I1  ,LP1)
     *       + C3B(I)*W3I(IP-I1,  LP1)   + C4B(I)*W3I(IP+I2  ,LP1)
      TEM2   = C1B(I)*W3I(IP+1,   LP1)   + C2B(I)*W3I(IP+1+I1,LP1)
     *       + C3B(I)*W3I(IP+1-I1,LP1)   + C4B(I)*W3I(IP+1+I2,LP1)
C 
      WK7(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
      endif

C
      TEM1   = C5B(I)*W3I(IP+I1,   LP2)   + WK8(I)*W3I(IP,LP2)
      TEM2   = C5B(I)*W3I(IP+1+I1, LP2)   + WK8(I)*W3I(IP+1,LP2)
C
      WK8(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
   70 CONTINUE
C
C
      DO 90 I=1,IJM0M0
      LV = KPR(I)
      IF (LV .GE. 2 .AND. LV .LE. LMM2) THEN
         TEM  = SGI(I)
         TEMA = (TEM-SGB(LV+1))*(TEM-SGB(LV+2))
         TEMB = (TEM-SGB(LV-1))*(TEM-SGB(LV))
C
         W1O(I) = TEMA * ((TEM-SGB(LV))   * S1A(LV,1) * WK1(I)
     *                  + (TEM-SGB(LV-1)) * S1A(LV,2) * WK2(I))
     *          + TEMB * ((TEM-SGB(LV+2)) * S1A(LV,3) * WK3(I)
     *                  + (TEM-SGB(LV+1)) * S1A(LV,4) * WK4(I))
C
         W3O(I) = TEMA * ((TEM-SGB(LV))   * S1A(LV,1) * WK5(I)
     *                  + (TEM-SGB(LV-1)) * S1A(LV,2) * WK6(I))
     *          + TEMB * ((TEM-SGB(LV+2)) * S1A(LV,3) * WK7(I)
     *                  + (TEM-SGB(LV+1)) * S1A(LV,4) * WK8(I))
C
      ELSEif (lv .eq. 1 .or. lv .eq. lmm1) then
         TEM = (SGI(I) - SGB(LV)) * SGX(LV)
         W1O(I) = (1.0 - TEM) * WK2(I) + TEM * WK3(I)
         W3O(I) = (1.0 - TEM) * WK6(I) + TEM * WK7(I)
      else
         W1O(I) = WK2(I)
         W3O(I) = WK6(I)
      ENDIF
   90 CONTINUE
C
  500 CONTINUE
C
      RETURN
      END
      SUBROUTINE BIC1n1(W1I, W1O, PHII, BLMI, SGI
     *,               IWK, IPR, JPR, KPR
     *,               C1A, C2A, C3A, C4A, C1B, C2B, C3B, C4B
     *,               WK1, WK2, WK3, WK4
     *,               SGX, SGB, S1A, FACLAM,FACPHI,JA,JB)
C
C
C  USES A BICUBIC SCHEME TO INTERPOLATE THE FIELD WIN TO THE POINT
C  (PHIIN,BLMIN) AND PUTS THE INTERPOLATED VALUE IN WOUT.
C  ON A C-GRID, FACLAM AND FACPHI MUST BE SET AS FOLLOWS DEPENDING ON
C  WHETHER THE POINTS BEING USED ARE THE U V OR FI POINTS.
C  U POINTS: FACLAM = 0.5; FACPHI = 1.0
C  V POINTS: FACLAM = 1.0; FACPHI = 0.5
C  H POINTS: FACLAM = 1.0; FACPHI = 1.0
C
#include "slgcm.slgprm"
C
      DIMENSION   W1I(IJP3M0,LM), W1O(IJM0M0)
     *,         PHII(IJM0M0), BLMI(IJM0M0), SGI(IJM0M0)
     *,          SGB(LM)
C
      DIMENSION  C1A(IJM0M0), C2A(IJM0M0), C3A(IJM0M0), C4A(IJM0M0)
      DIMENSION  C1B(IJM0M0), C2B(IJM0M0), C3B(IJM0M0), C4B(IJM0M0)
      dimension  c5a(ijm0m0), c5b(ijm0m0)
      DIMENSION  WK1(IJM0M0), WK2(IJM0M0), WK3(IJM0M0), WK4(IJM0M0)
     *,          IWK(IJM0M0), IPR(IJM0M0),  JPR(IJM0M0),  KPR(IJM0M0)
C
      DIMENSION  SGX(LM), S1A(LM,4)
C
C
      FACDLM = (2.-FACLAM)*DLM
C
      jlin = 6
      ijs = (jlin+2) * imp3
      ijn = ijp3m0 - (jlin+1) * imp3
C
      JL  = JB - JA + 1
      IJL = IMP3*JL
      IJA = (JA-1)*IMP3
      IJB = IJA + 1
      IJE = IJA + IJL
C
C
C
C  FIRST MAKE SURE DLM < LAM* < (TWOPI+DLM) FOR THE V AND FI POINTS
C  AND THAT 1.5*DLM < LAM* < (TWOPI+1.5*DLM) FOR THE U POINTS
C
      DO 30 I=IJB,IJE
      WK1(I) = (BLMI(I)-FACDLM) * (ONE/TWOPI)
      IWK(I) =  WK1(I)
      IPR(I) =  IWK(I)
   30 CONTINUE
      DO 35 I=ijb,ije
      IF (WK1(I) .LT. 0.0) IWK(I) = IPR(I) - 1
   35 CONTINUE
      DO 40 I=ijb,ije
      TEM0   = (WK1(I) - IWK(I))*IM + TWO
      IPR(I) = TEM0
      TEM    = TEM0 - IPR(I)
      c5a(i) = tem
C
      TEM1   =  TEM * TEM
      TEM2   = (2.0 - TEM) * 0.5
C
      C1A(I) =  TEM2 * (1.0 - TEM1)
      C2A(I) = (TEM1 + TEM) * TEM2
      C3A(I) = (TEM1 - TEM) * TEM2 * THIRD
      C4A(I) = (TEM1 - 1.0) * TEM  * SIXTH
   40 CONTINUE
      DO 45 I=ijb,ije
      IF (IPR(I) .EQ. 1)  IPR(I) = IMP1
   45 CONTINUE
C
      DO 47 I=ijb,ije
      TEM0   = (PHII(I) + PIO2)*RDPH + FACPHI
      JPR(I) = TEM0

      TEM  = TEM0 - JPR(I)
      c5b(i) = tem
C
      TEM1 =  TEM * TEM
      TEM2 = (2.0 - TEM) * 0.5
C
      C1B(I) =  TEM2 * (1.0 - TEM1)
      C2B(I) = (TEM1 + TEM) * TEM2
      C3B(I) = (TEM1 - TEM) * TEM2 * THIRD
      C4B(I) = (TEM1 - 1.0) * TEM * SIXTH
   47 CONTINUE
C
C
      DO 70 I=ijb,ije
      IP = IPR(I) + (JPR(I)+1)*IMP3
C
      LP = KPR(I)
      LM1 = LP - 1
      LP1 = LP + 1
      LP2 = LP + 2
      IF (LP .LT. 2) THEN
         IF (LP .EQ. 0) LP=1
         LM1 = LP
      ENDIF
      IF (LP .GT. LMM2) THEN
         IF (LP .EQ. LM) LP1 = LM
         LP2 = LP1
      ENDIF
C
C
      WK4(I) = 1.0 - C5B(I)
      TEM1   = WK4(I)*W1I(IP,   LM1) + C5B(I)*W1I(IP+I1  ,LM1)
      TEM2   = WK4(I)*W1I(IP+1, LM1) + C5B(I)*W1I(IP+1+I1,LM1)
C
      WK1(I) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
      if (ip .gt. ijs .and. ip .lt. ijn) then
C
      TEM1   = C1B(I)*W1I(IP,     LP)    + C2B(I)*W1I(IP+I1  ,LP)
     *       + C3B(I)*W1I(IP-I1,  LP)    + C4B(I)*W1I(IP+I2  ,LP)
      TEM2   = C1B(I)*W1I(IP+1,   LP)    + C2B(I)*W1I(IP+1+I1,LP)
     *       + C3B(I)*W1I(IP+1-I1,LP)    + C4B(I)*W1I(IP+1+I2,LP)
      TEM3   = C1B(I)*W1I(IP-1,   LP)    + C2B(I)*W1I(IP-1+I1,LP)
     *       + C3B(I)*W1I(IP-1-I1,LP)    + C4B(I)*W1I(IP-1+I2,LP)
      TEM4   = C1B(I)*W1I(IP+2   ,LP)    + C2B(I)*W1I(IP+2+I1,LP)
     *       + C3B(I)*W1I(IP+2-I1,LP)    + C4B(I)*W1I(IP+2+I2,LP)
C
      WK2(I) = C1A(I) * TEM1 + C2A(I) * TEM2
     *       + C3A(I) * TEM3 + C4A(I) * TEM4
C
      TEM1   = C1B(I)*W1I(IP,     LP1)   + C2B(I)*W1I(IP+I1  ,LP1)
     *       + C3B(I)*W1I(IP-I1,  LP1)   + C4B(I)*W1I(IP+I2  ,LP1)
      TEM2   = C1B(I)*W1I(IP+1,   LP1)   + C2B(I)*W1I(IP+1+I1,LP1)
     *       + C3B(I)*W1I(IP+1-I1,LP1)   + C4B(I)*W1I(IP+1+I2,LP1)
      TEM3   = C1B(I)*W1I(IP-1,   LP1)   + C2B(I)*W1I(IP-1+I1,LP1)
     *       + C3B(I)*W1I(IP-1-I1,LP1)   + C4B(I)*W1I(IP-1+I2,LP1)
      TEM4   = C1B(I)*W1I(IP+2   ,LP1)   + C2B(I)*W1I(IP+2+I1,LP1)
     *       + C3B(I)*W1I(IP+2-I1,LP1)   + C4B(I)*W1I(IP+2+I2,LP1)
C
      WK3(I) = C1A(I) * TEM1 + C2A(I) * TEM2
     *       + C3A(I) * TEM3 + C4A(I) * TEM4
C
      else
C
      TEM1   = C1B(I)*W1I(IP,     LP)   + C2B(I)*W1I(IP+I1  ,LP)
     *       + C3B(I)*W1I(IP-I1,  LP)   + C4B(I)*W1I(IP+I2  ,LP)
      TEM2   = C1B(I)*W1I(IP+1,   LP)   + C2B(I)*W1I(IP+1+I1,LP)
     *       + C3B(I)*W1I(IP+1-I1,LP)   + C4B(I)*W1I(IP+1+I2,LP)
C
      WK2(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
      TEM1   = C1B(I)*W1I(IP,     LP1)   + C2B(I)*W1I(IP+I1  ,LP1)
     *       + C3B(I)*W1I(IP-I1,  LP1)   + C4B(I)*W1I(IP+I2  ,LP1)
      TEM2   = C1B(I)*W1I(IP+1,   LP1)   + C2B(I)*W1I(IP+1+I1,LP1)
     *       + C3B(I)*W1I(IP+1-I1,LP1)   + C4B(I)*W1I(IP+1+I2,LP1)
C
      WK3(I  ) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
      endif
C
C
      TEM1   = WK4(I)*W1I(IP,   LP2) + C5B(I)*W1I(IP+I1  ,LP2)
      TEM2   = WK4(I)*W1I(IP+1, LP2) + C5B(I)*W1I(IP+1+I1,LP2)
C
      WK4(I) = C5A(I) * TEM2 + (1.0 - C5A(I)) * TEM1
C
C
   70 CONTINUE
C
C
      DO 90 I=IJB,IJE
      LV = KPR(I)
      IF (LV .GE. 2 .AND. LV .LE. LMM2) THEN
         TEM  = SGI(I)
         TEMA = (TEM-SGB(LV+1))*(TEM-SGB(LV+2))
         TEMB = (TEM-SGB(LV-1))*(TEM-SGB(LV))
C
         TEM1 = (TEM-SGB(LV))   * TEMA * S1A(LV,1)
         TEM2 = (TEM-SGB(LV-1)) * TEMA * S1A(LV,2)
         TEM3 = (TEM-SGB(LV+2)) * TEMB * S1A(LV,3)
         TEM4 = (TEM-SGB(LV+1)) * TEMB * S1A(LV,4)
C
C        TEM1 = (TEM-SGB(LV))  *(TEM-SGB(LV+1))*(TEM-SGB(LV+2))*S1A(LV)
C        TEM2 = (TEM-SGB(LV-1))*(TEM-SGB(LV+1))*(TEM-SGB(LV+2))*S1B(LV)
C        TEM3 = (TEM-SGB(LV-1))*(TEM-SGB(LV))  *(TEM-SGB(LV+2))*S1C(LV)
C        TEM4 = (TEM-SGB(LV-1))*(TEM-SGB(LV))  *(TEM-SGB(LV+1))*S1D(LV)
C
         W1O(I) = TEM1*WK1(I) + TEM2*WK2(I) + TEM3*WK3(I) +TEM4*WK4(I)
      ELSEif (lv .eq. 1 .or. lv .eq. lmm1) then
         TEM = (SGI(I) - SGB(LV)) * SGX(LV)
         W1O(I) = (1.0 - TEM) * WK2(I) + TEM * WK3(I)
      else
         W1O(I) = WK2(I)
      ENDIF
   90 CONTINUE
      CALL WRAP(IJB,IJE,W1O)
  100 CONTINUE
C
      RETURN
      END
      SUBROUTINE CHOPUV(FLD, IND)
C
C*********************************************************************
C  FILL IN THE 'GHOST' VALUES AT U - POINTS J=1 AND J=JM
C*********************************************************************
C
#include "slgcm.slgprm"
#include "slgcm.slgclm"
C
      DIMENSION FLD(IJM0M0), TXS(IM), TXN(IM)
      DIMENSION ANN(3), ANS(3), BNN(3), BNS(3)
C
C
      if (ind .eq. 0) then
         isf = 0
         inf = ijm2m0
      else
         isf = i1
         inf = ijm2m0
      endif
      tems = 0.0
      temn = 0.0
c     do 5 i=1,im
c     tems = tems + fld(i+isf)
c     temn = temn + fld(i+inf)
c  5  continue
c     tems = tems / im
c     temn = temn / im
      txs(:) = tems
      txn(:) = temn
C
      do 30 n=1,1
      ann(n) = 0.0
      ans(n) = 0.0
      bnn(n) = 0.0
      bns(n) = 0.0
      nn     = 2 * n - ind
      DO 10 I=1,IM
      ann(n) = ann(n) + fld(inf+i) * clm(i,nn)
      bnn(n) = bnn(n) + fld(inf+i) * slm(i,nn)
C
      ans(n) = ans(n) + fld(i+isf)     * clm(i,nn)
      bns(n) = bns(n) + fld(i+isf)     * slm(i,nn)
 10   CONTINUE
C
      ann(n) = ann(n) * twooim
      bnn(n) = bnn(n) * twooim
      ans(n) = ans(n) * twooim
      bns(n) = bns(n) * twooim
      do 20 i=1,im
      txs(i) = txs(i) + ans(n) * clm(i,nn) + bns(n) * slm(i,nn)
      txn(i) = txn(i) + ann(n) * clm(i,nn) + bnn(n) * slm(i,nn)
  20  continue
  30  continue
C
      do 40 i=1,im
      fld(i+isf) = txs(i)
      fld(i+inf) = txn(i)
   40 continue
C
      fld(IMP1+isf) = fld(isf+1)
      fld(IMP2+isf) = fld(isf+2)
      fld(IMP3+isf) = fld(isf+3)
c
      fld(IMP1+inf) = fld(1+inf)
      fld(IMP2+inf) = fld(2+inf)
      fld(IMP3+inf) = fld(3+inf)
C

      RETURN
      END
      SUBROUTINE CHOPT(FLD, JIND, NWAV)
C
C*********************************************************************
C  FILL IN THE 'GHOST' VALUES AT U - POINTS J=1 AND J=JM
C*********************************************************************
C
#include "slgcm.slgprm"
#include "slgcm.slgclm"
C
      DIMENSION FLD(IMP3,JM), TXS(IM), TXN(IM)
      DIMENSION ANN(3), ANS(3), BNN(3), BNS(3)
C
C
      js = 1  + jind
      jn = jm - jind
      tems = 0.0
      temn = 0.0
      do 5 i=1,im
      tems = tems + fld(i,js)
      temn = temn + fld(i,jn)
   5  continue
      tems = tems / im
      temn = temn / im
      txs(:) = tems
      txn(:) = temn
C
      do 30 n=1,nwav
      ann(n) = 0.0
      ans(n) = 0.0
      bnn(n) = 0.0
      bns(n) = 0.0
      nn     = 2 * n 
      DO 10 I=1,IM
      ann(n) = ann(n) + fld(i,jn) * clm(i,nn)
      bnn(n) = bnn(n) + fld(i,jn) * slm(i,nn)
C
      ans(n) = ans(n) + fld(i,js) * clm(i,nn)
      bns(n) = bns(n) + fld(i,js) * slm(i,nn)
 10   CONTINUE
C
      ann(n) = ann(n) * twooim
      bnn(n) = bnn(n) * twooim
      ans(n) = ans(n) * twooim
      bns(n) = bns(n) * twooim
      do 20 i=1,im
      txs(i) = txs(i) + ans(n) * clm(i,nn) + bns(n) * slm(i,nn)
      txn(i) = txn(i) + ann(n) * clm(i,nn) + bnn(n) * slm(i,nn)
  20  continue
  30  continue
C
      do 40 i=1,im
      fld(i,js) = txs(i)
      fld(i,jn) = txn(i)
   40 continue
C
      fld(IMP1,js) = fld(1,js)
      fld(IMP2,js) = fld(2,js)
      fld(IMP3,js) = fld(3,js)
c
      fld(IMP1,jn) = fld(1,jn)
      fld(IMP2,jn) = fld(2,jn)
      fld(IMP3,jn) = fld(3,jn)
C

      RETURN
      END
      SUBROUTINE BICPn2(VST, W1I, W2I, W3I, W1O, W2O, W3O
     *,                 PHII, BLMI,SGI
     *,                 IWK, IPR, JPR, KPR
     *,                 C1A, C2A, C3A, C4A, C1B, C2B, C3B, C4B
     *,                 WK1, WK2, WK3, WK4, WK5, WK6, WK7, WK8
     *,                 SGX, DTOA, SGB, S1A, FACLAM,FACPHI,JA,JB,UST)
C
C
C  USES A BICUBIC SCHEME TO INTERPOLATE THE FIELD W1I TO THE POINT
C  (PHII,BLMI) AND PUTS THE INTERPOLATED VALUE IN W1O.
C  ON A C-GRID, FACLAM AND FACPHI MUST BE SET AS FOLLOWS DEPENDING ON
C  WHETHER THE POINTS BEING USED ARE THE U V OR FI POINTS.
C  U POINTS: FACLAM = 0.5; FACPHI = 1.0
C  V POINTS: FACLAM = 1.0; FACPHI = 0.5
C  H POINTS: FACLAM = 1.0; FACPHI = 1.0
C
#include "slgcm.slgprm"
#include "slgcm.slgclm"
C
      DIMENSION  W1I(IJP3M0,LM),  W2I(IJP3M0), W3I(IJP3M0,LM)
     *,          W1O(IJM0M0),     W2O(IJM0M0), W3O(IJM0M0)
C
      DIMENSION  PHII(IJM0M0),    BLMI(IJM0M0), SGI(IJM0M0)
     *,          VST(IJM0M0),     SGB(LM),      UST(IJM0M0)
C
      DIMENSION  C1A(IJM0M0), C2A(IJM0M0), C3A(IJM0M0), C4A(IJM0M0)
      DIMENSION  C1B(IJM0M0), C2B(IJM0M0), C3B(IJM0M0), C4B(IJM0M0)
      DIMENSION  WK1(IJM0M0), WK2(IJM0M0), WK3(IJM0M0), WK4(IJM0M0)
     *,          WK5(IJM0M0), WK6(IJM0M0), WK7(IJM0M0), WK8(IJM0M0)
     *,          IWK(IJM0M0), IPR(IJM0M0),  JPR(IJM0M0),  KPR(IJM0M0)
     *,          CVV(2),   GLM00(2)
      DIMENSION  SGX(LM), S1A(LM,4)
      data ipx/13200/
      SAVE
C
C
      FACDLM = (2.-FACLAM)*DLM

      JL  = JB - JA + 1
      IJL = IMP3*JL
      IJA = (JA-1)*IMP3
      IJB = IJA + 1
      IJE = IJA + IJL
C
C     CALCULATE SOUTH & NORTH POLE VALUES
C
      CALL VPOLE (IM, DLM, TWOPI, TWOOIM, PI
     *,       VST, CLM(1,2), SLM(1,2), CVV(1), GLM00(1))
C     CALL UPOLE (IM, DLM, TWOPI, TWOOIM, PI
C    *,       UST(I1+1,L), CLM(1,1), SLM(1,1), TEM1    , TEM2      )
C    *,       VST(I1+1,L), CLM(1,1), SLM(1,1), CVV(L,1), GLM00(L,1))
C     CVV(L,1)   = (CVV(L,1)   + TEM1) * 0.5
C     GLM00(L,1) = (GLM00(L,1) + TEM2) * 0.5
C
      CALL VPOLE (IM, DLM, TWOPI, TWOOIM, ZERO
     *,      VST(IJM2M0+1), CLM(1,2), SLM(1,2), CVV(2), GLM00(2))
C     CALL UPOLE (IM, DLM, TWOPI, TWOOIM, ZERO
C    *,      UST(IJM2M0+1,L), CLM(1,1), SLM(1,1), TEM1    , TEM2      )
C    *,      VST(IJM2M0+1,L), CLM(1,1), SLM(1,1), CVV(L,2), GLM00(L,2))
C     CVV(L,2)   = (CVV(L,2)   + TEM1) * 0.5
C     GLM00(L,2) = (GLM00(L,2) + TEM2) * 0.5

      CVV(:) = CVV(:) * DTOA
C
C
C  FIRST MAKE SURE DLM < LAM* < (TWOPI+DLM) FOR THE V AND FI POINTS
C  AND THAT 1.5*DLM < LAM* < (TWOPI+1.5*DLM) FOR THE U POINTS
C
      DO 30 I=ijb,ije
      WK1(I) = (BLMI(I)-FACDLM) * (ONE/TWOPI)
      IWK(I) =  WK1(I)
      IPR(I) =  IWK(I)
   30 CONTINUE
      DO 35 I=ijb,ije
      IF (WK1(I) .LT. 0.0) IWK(I) = IPR(I) - 1
   35 CONTINUE
      DO 40 I=ijb,ije
      TEM0   = (WK1(I) - IWK(I))*IM + TWO
      IPR(I) = TEM0
      TEM    = TEM0 - IPR(I)
C
      TEM1   =  TEM * TEM
      TEM2   = (2.0 - TEM) * 0.5
C
      C1A(I) =  TEM2 * (1.0 - TEM1)
      C2A(I) = (TEM1 + TEM) * TEM2
      C3A(I) = (TEM1 - TEM) * TEM2 * THIRD
      C4A(I) = (TEM1 - 1.0) * TEM  * SIXTH
   40 CONTINUE
      DO 45 I=ijb,ije
      IF (IPR(I) .EQ. 1)  IPR(I) = IMP1
   45 CONTINUE
C
      DO 47 I=ijb,ije
      TEM0   = (PHII(I) + PIO2)*RDPH + FACPHI
      JPR(I) = TEM0

      TEM  = TEM0 - JPR(I)
C
      TEM1 =  TEM * TEM
      TEM2 = (2.0 - TEM) * 0.5
C
      C1B(I) =  TEM2 * (1.0 - TEM1)
      C2B(I) = (TEM1 + TEM) * TEM2
      C3B(I) = (TEM1 - TEM) * TEM2 * THIRD
      C4B(I) = (TEM1 - 1.0) * TEM * SIXTH
   47 CONTINUE
C
C     CALCULATE POLE VALUES
C
      LP = KPR(1)
      LM1 = LP - 1
      LP1 = LP + 1
      LP2 = LP + 2

      IF (LP .LT. 2) THEN
         IF (LP .EQ. 0) LP=1
         LM1 = LP
      ENDIF
      IF (LP .GT. LMM2) THEN
         IF (LP .EQ. LM) LP1 = LM
         LP2 = LP1
      ENDIF

C
      CALL FIPOLE(W1I(1,LM1), WK1, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W1I(1,LP) , WK2, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W1I(1,LP1), WK3, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W1I(1,LP2), WK4, CVV(1), GLM00(1), .TRUE.)
C
      CALL FIPOLE(W2I, W2O, CVV(1), GLM00(1), .TRUE.)
C
      CALL FIPOLE(W3I(1,LM1), WK5, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W3I(1,LP) , WK6, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W3I(1,LP1), WK7, CVV(1), GLM00(1), .TRUE.)
      CALL FIPOLE(W3I(1,LP2), WK8, CVV(1), GLM00(1), .TRUE.)
C
      IJP1 = IJM1M0 + 1
      LP = KPR(IJP1)
      LM1 = LP - 1
      LP1 = LP + 1
      LP2 = LP + 2
      IF (LP .LT. 2) THEN
         IF (LP .EQ. 0) LP=1
         LM1 = LP
      ENDIF
      IF (LP .GT. LMM2) THEN
         IF (LP .EQ. LM) LP1 = LM
         LP2 = LP1
      ENDIF
C
      CALL FIPOLE(W1I(1,LM1), WK1, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W1I(1,LP) , WK2, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W1I(1,LP1), WK3, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W1I(1,LP2), WK4, CVV(2), GLM00(2), .FALSE.)
C
      CALL FIPOLE(W2I, W2O, CVV(2), GLM00(2), .FALSE.)
C
      CALL FIPOLE(W3I(1,LM1), WK5, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W3I(1,LP) , WK6, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W3I(1,LP1), WK7, CVV(2), GLM00(2), .FALSE.)
      CALL FIPOLE(W3I(1,LP2), WK8, CVV(2), GLM00(2), .FALSE.)
C
C     INTERIOR
C
      DO 70 I=ijb,ije
      IP = IPR(I) + (JPR(I)+1)*IMP3
C
      LP = KPR(I)
      LM1 = LP - 1
      LP1 = LP + 1
      LP2 = LP + 2
      IF (LP .LT. 2) THEN
         IF (LP .EQ. 0) LP=1
         LM1 = LP
      ENDIF
      IF (LP .GT. LMM2) THEN
         IF (LP .EQ. LM) LP1 = LM
         LP2 = LP1
      ENDIF
C
C
      TEM1   = C1B(I)*W1I(IP,     LM1)    + C2B(I)*W1I(IP+I1  ,LM1)
     *       + C3B(I)*W1I(IP-I1,  LM1)    + C4B(I)*W1I(IP+I2  ,LM1)
      TEM2   = C1B(I)*W1I(IP+1,   LM1)    + C2B(I)*W1I(IP+1+I1,LM1)
     *       + C3B(I)*W1I(IP+1-I1,LM1)    + C4B(I)*W1I(IP+1+I2,LM1)
      TEM3   = C1B(I)*W1I(IP-1,   LM1)    + C2B(I)*W1I(IP-1+I1,LM1)
     *       + C3B(I)*W1I(IP-1-I1,LM1)    + C4B(I)*W1I(IP-1+I2,LM1)
      TEM4   = C1B(I)*W1I(IP+2   ,LM1)    + C2B(I)*W1I(IP+2+I1,LM1)
     *       + C3B(I)*W1I(IP+2-I1,LM1)    + C4B(I)*W1I(IP+2+I2,LM1)
C
      WK1(I)   = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
C
      TEM1   = C1B(I)*W1I(IP,     LP)    + C2B(I)*W1I(IP+I1  ,LP)
     *       + C3B(I)*W1I(IP-I1,  LP)    + C4B(I)*W1I(IP+I2  ,LP)
      TEM2   = C1B(I)*W1I(IP+1,   LP)    + C2B(I)*W1I(IP+1+I1,LP)
     *       + C3B(I)*W1I(IP+1-I1,LP)    + C4B(I)*W1I(IP+1+I2,LP)
      TEM3   = C1B(I)*W1I(IP-1,   LP)    + C2B(I)*W1I(IP-1+I1,LP)
     *       + C3B(I)*W1I(IP-1-I1,LP)    + C4B(I)*W1I(IP-1+I2,LP)
      TEM4   = C1B(I)*W1I(IP+2   ,LP)    + C2B(I)*W1I(IP+2+I1,LP)
     *       + C3B(I)*W1I(IP+2-I1,LP)    + C4B(I)*W1I(IP+2+I2,LP)
C
      WK2(I)   = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
      TEM1   = C1B(I)*W1I(IP,     LP1)   + C2B(I)*W1I(IP+I1  ,LP1)
     *       + C3B(I)*W1I(IP-I1,  LP1)   + C4B(I)*W1I(IP+I2  ,LP1)
      TEM2   = C1B(I)*W1I(IP+1,   LP1)   + C2B(I)*W1I(IP+1+I1,LP1)
     *       + C3B(I)*W1I(IP+1-I1,LP1)   + C4B(I)*W1I(IP+1+I2,LP1)
      TEM3   = C1B(I)*W1I(IP-1,   LP1)   + C2B(I)*W1I(IP-1+I1,LP1)
     *       + C3B(I)*W1I(IP-1-I1,LP1)   + C4B(I)*W1I(IP-1+I2,LP1)
      TEM4   = C1B(I)*W1I(IP+2   ,LP1)   + C2B(I)*W1I(IP+2+I1,LP1)
     *       + C3B(I)*W1I(IP+2-I1,LP1)   + C4B(I)*W1I(IP+2+I2,LP1)
C
      WK3(I  ) = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
      TEM1   = C1B(I)*W1I(IP,     LP2)   + C2B(I)*W1I(IP+I1  ,LP2)
     *       + C3B(I)*W1I(IP-I1,  LP2)   + C4B(I)*W1I(IP+I2  ,LP2)
      TEM2   = C1B(I)*W1I(IP+1,   LP2)   + C2B(I)*W1I(IP+1+I1,LP2)
     *       + C3B(I)*W1I(IP+1-I1,LP2)   + C4B(I)*W1I(IP+1+I2,LP2)
      TEM3   = C1B(I)*W1I(IP-1,   LP2)   + C2B(I)*W1I(IP-1+I1,LP2)
     *       + C3B(I)*W1I(IP-1-I1,LP2)   + C4B(I)*W1I(IP-1+I2,LP2)
      TEM4   = C1B(I)*W1I(IP+2   ,LP2)   + C2B(I)*W1I(IP+2+I1,LP2)
     *       + C3B(I)*W1I(IP+2-I1,LP2)   + C4B(I)*W1I(IP+2+I2,LP2)
C
      WK4(I  ) = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
C
      TEM1   = C1B(I)*W2I(IP     )    + C2B(I)*W2I(IP+I1  )
     *       + C3B(I)*W2I(IP-I1  )    + C4B(I)*W2I(IP+I2  )
      TEM2   = C1B(I)*W2I(IP+1   )    + C2B(I)*W2I(IP+1+I1)
     *       + C3B(I)*W2I(IP+1-I1)    + C4B(I)*W2I(IP+1+I2)
      TEM3   = C1B(I)*W2I(IP-1   )    + C2B(I)*W2I(IP-1+I1)
     *       + C3B(I)*W2I(IP-1-I1)    + C4B(I)*W2I(IP-1+I2)
      TEM4   = C1B(I)*W2I(IP+2   )    + C2B(I)*W2I(IP+2+I1)
     *       + C3B(I)*W2I(IP+2-I1)    + C4B(I)*W2I(IP+2+I2)
C
      W2O(I) = C1A(I) * TEM1 + C2A(I) * TEM2
     *       + C3A(I) * TEM3 + C4A(I) * TEM4
C
C
      TEM1   = C1B(I)*W3I(IP,     LM1)    + C2B(I)*W3I(IP+I1  ,LM1)
     *       + C3B(I)*W3I(IP-I1,  LM1)    + C4B(I)*W3I(IP+I2  ,LM1)
      TEM2   = C1B(I)*W3I(IP+1,   LM1)    + C2B(I)*W3I(IP+1+I1,LM1)
     *       + C3B(I)*W3I(IP+1-I1,LM1)    + C4B(I)*W3I(IP+1+I2,LM1)
      TEM3   = C1B(I)*W3I(IP-1,   LM1)    + C2B(I)*W3I(IP-1+I1,LM1)
     *       + C3B(I)*W3I(IP-1-I1,LM1)    + C4B(I)*W3I(IP-1+I2,LM1)
      TEM4   = C1B(I)*W3I(IP+2   ,LM1)    + C2B(I)*W3I(IP+2+I1,LM1)
     *       + C3B(I)*W3I(IP+2-I1,LM1)    + C4B(I)*W3I(IP+2+I2,LM1)
C
      WK5(I)   = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
C
C
      TEM1   = C1B(I)*W3I(IP,     LP)    + C2B(I)*W3I(IP+I1  ,LP)
     *       + C3B(I)*W3I(IP-I1,  LP)    + C4B(I)*W3I(IP+I2  ,LP)
      TEM2   = C1B(I)*W3I(IP+1,   LP)    + C2B(I)*W3I(IP+1+I1,LP)
     *       + C3B(I)*W3I(IP+1-I1,LP)    + C4B(I)*W3I(IP+1+I2,LP)
      TEM3   = C1B(I)*W3I(IP-1,   LP)    + C2B(I)*W3I(IP-1+I1,LP)
     *       + C3B(I)*W3I(IP-1-I1,LP)    + C4B(I)*W3I(IP-1+I2,LP)
      TEM4   = C1B(I)*W3I(IP+2   ,LP)    + C2B(I)*W3I(IP+2+I1,LP)
     *       + C3B(I)*W3I(IP+2-I1,LP)    + C4B(I)*W3I(IP+2+I2,LP)
C
      WK6(I)   = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
      TEM1   = C1B(I)*W3I(IP,     LP1)   + C2B(I)*W3I(IP+I1  ,LP1)
     *       + C3B(I)*W3I(IP-I1,  LP1)   + C4B(I)*W3I(IP+I2  ,LP1)
      TEM2   = C1B(I)*W3I(IP+1,   LP1)   + C2B(I)*W3I(IP+1+I1,LP1)
     *       + C3B(I)*W3I(IP+1-I1,LP1)   + C4B(I)*W3I(IP+1+I2,LP1)
      TEM3   = C1B(I)*W3I(IP-1,   LP1)   + C2B(I)*W3I(IP-1+I1,LP1)
     *       + C3B(I)*W3I(IP-1-I1,LP1)   + C4B(I)*W3I(IP-1+I2,LP1)
      TEM4   = C1B(I)*W3I(IP+2   ,LP1)   + C2B(I)*W3I(IP+2+I1,LP1)
     *       + C3B(I)*W3I(IP+2-I1,LP1)   + C4B(I)*W3I(IP+2+I2,LP1)
C
      WK7(I  ) = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
      TEM1   = C1B(I)*W3I(IP,     LP2)   + C2B(I)*W3I(IP+I1  ,LP2)
     *       + C3B(I)*W3I(IP-I1,  LP2)   + C4B(I)*W3I(IP+I2  ,LP2)
      TEM2   = C1B(I)*W3I(IP+1,   LP2)   + C2B(I)*W3I(IP+1+I1,LP2)
     *       + C3B(I)*W3I(IP+1-I1,LP2)   + C4B(I)*W3I(IP+1+I2,LP2)
      TEM3   = C1B(I)*W3I(IP-1,   LP2)   + C2B(I)*W3I(IP-1+I1,LP2)
     *       + C3B(I)*W3I(IP-1-I1,LP2)   + C4B(I)*W3I(IP-1+I2,LP2)
      TEM4   = C1B(I)*W3I(IP+2   ,LP2)   + C2B(I)*W3I(IP+2+I1,LP2)
     *       + C3B(I)*W3I(IP+2-I1,LP2)   + C4B(I)*W3I(IP+2+I2,LP2)
C
      WK8(I  ) = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4

C
   70 CONTINUE
C
C
      DO 90 I=1,IJM0M0
      LV = KPR(I)
      IF (LV .GE. 2 .AND. LV .LE. LMM2) THEN
 
c     LVB = KPR(I)
c     IF (LVB .GE. 2 .AND. LVB .LE. LMM2) THEN
c        LV = LVB 
c     ELSEIF (LVB .EQ. 1) THEN 
c        LV = 2 
c     ELSE 
c        LV = LMM2 
c     ENDIF 
c     IF (LVB .GE. 1 .AND. LVB .LE. LMM1) THEN

         TEM  = SGI(I)
         TEMA = (TEM-SGB(LV+1))*(TEM-SGB(LV+2))
         TEMB = (TEM-SGB(LV-1))*(TEM-SGB(LV))
C
         W1O(I) = TEMA * ((TEM-SGB(LV))   * S1A(LV,1) * WK1(I)
     *                  + (TEM-SGB(LV-1)) * S1A(LV,2) * WK2(I))
     *          + TEMB * ((TEM-SGB(LV+2)) * S1A(LV,3) * WK3(I)
     *                  + (TEM-SGB(LV+1)) * S1A(LV,4) * WK4(I))
C
         W3O(I) = TEMA * ((TEM-SGB(LV))   * S1A(LV,1) * WK5(I)
     *                  + (TEM-SGB(LV-1)) * S1A(LV,2) * WK6(I))
     *          + TEMB * ((TEM-SGB(LV+2)) * S1A(LV,3) * WK7(I)
     *                  + (TEM-SGB(LV+1)) * S1A(LV,4) * WK8(I))
C
      ELSEif (lv .eq. 1 .or. lv .eq. lmm1) then
         TEM = (SGI(I) - SGB(LV)) * SGX(LV)
         W1O(I) = (1.0 - TEM) * WK2(I) + TEM * WK3(I)
         W3O(I) = (1.0 - TEM) * WK6(I) + TEM * WK7(I)
      else
         W1O(I) = WK2(I)
         W3O(I) = WK6(I)
      ENDIF
   90 CONTINUE
C
  500 CONTINUE
C
      RETURN
      END
      SUBROUTINE BIC1n2(W1I, W1O, PHII, BLMI, SGI
     *,               IWK, IPR, JPR, KPR
     *,               C1A, C2A, C3A, C4A, C1B, C2B, C3B, C4B
     *,               WK1, WK2, WK3, WK4
     *,               SGX, SGB, S1A, FACLAM,FACPHI,JA,JB)
C
C
C  USES A BICUBIC SCHEME TO INTERPOLATE THE FIELD WIN TO THE POINT
C  (PHIIN,BLMIN) AND PUTS THE INTERPOLATED VALUE IN WOUT.
C  ON A C-GRID, FACLAM AND FACPHI MUST BE SET AS FOLLOWS DEPENDING ON
C  WHETHER THE POINTS BEING USED ARE THE U V OR FI POINTS.
C  U POINTS: FACLAM = 0.5; FACPHI = 1.0
C  V POINTS: FACLAM = 1.0; FACPHI = 0.5
C  H POINTS: FACLAM = 1.0; FACPHI = 1.0
C
#include "slgcm.slgprm"
C
      DIMENSION   W1I(IJP3M0,LM), W1O(IJM0M0)
     *,         PHII(IJM0M0), BLMI(IJM0M0), SGI(IJM0M0)
     *,          SGB(LM)
C
      DIMENSION  C1A(IJM0M0), C2A(IJM0M0), C3A(IJM0M0), C4A(IJM0M0)
      DIMENSION  C1B(IJM0M0), C2B(IJM0M0), C3B(IJM0M0), C4B(IJM0M0)
      dimension  c5a(ijm0m0), c5b(ijm0m0)
      DIMENSION  WK1(IJM0M0), WK2(IJM0M0), WK3(IJM0M0), WK4(IJM0M0)
     *,          IWK(IJM0M0), IPR(IJM0M0),  JPR(IJM0M0),  KPR(IJM0M0)
C
      DIMENSION  SGX(LM), S1A(LM,4)
C
C
      FACDLM = (2.-FACLAM)*DLM

      JL  = JB - JA + 1
      IJL = IMP3*JL
      IJA = (JA-1)*IMP3
      IJB = IJA + 1
      IJE = IJA + IJL
C
C
C
C  FIRST MAKE SURE DLM < LAM* < (TWOPI+DLM) FOR THE V AND FI POINTS
C  AND THAT 1.5*DLM < LAM* < (TWOPI+1.5*DLM) FOR THE U POINTS
C
      DO 30 I=IJB,IJE
      WK1(I) = (BLMI(I)-FACDLM) * (ONE/TWOPI)
      IWK(I) =  WK1(I)
      IPR(I) =  IWK(I)
   30 CONTINUE
      DO 35 I=ijb,ije
      IF (WK1(I) .LT. 0.0) IWK(I) = IPR(I) - 1
   35 CONTINUE
      DO 40 I=ijb,ije
      TEM0   = (WK1(I) - IWK(I))*IM + TWO
      IPR(I) = TEM0
      TEM    = TEM0 - IPR(I)
C
      TEM1   =  TEM * TEM
      TEM2   = (2.0 - TEM) * 0.5
C
      C1A(I) =  TEM2 * (1.0 - TEM1)
      C2A(I) = (TEM1 + TEM) * TEM2
      C3A(I) = (TEM1 - TEM) * TEM2 * THIRD
      C4A(I) = (TEM1 - 1.0) * TEM  * SIXTH
   40 CONTINUE
      DO 45 I=ijb,ije
      IF (IPR(I) .EQ. 1)  IPR(I) = IMP1
   45 CONTINUE
C
      DO 47 I=ijb,ije
      TEM0   = (PHII(I) + PIO2)*RDPH + FACPHI
      JPR(I) = TEM0

      TEM  = TEM0 - JPR(I)
C
      TEM1 =  TEM * TEM
      TEM2 = (2.0 - TEM) * 0.5
C
      C1B(I) =  TEM2 * (1.0 - TEM1)
      C2B(I) = (TEM1 + TEM) * TEM2
      C3B(I) = (TEM1 - TEM) * TEM2 * THIRD
      C4B(I) = (TEM1 - 1.0) * TEM * SIXTH
   47 CONTINUE
C
C
      DO 70 I=ijb,ije
      IP = IPR(I) + (JPR(I)+1)*IMP3
C
      LP = KPR(I)
      LM1 = LP - 1
      LP1 = LP + 1
      LP2 = LP + 2
      IF (LP .LT. 2) THEN
         IF (LP .EQ. 0) LP=1
         LM1 = LP
      ENDIF
      IF (LP .GT. LMM2) THEN
         IF (LP .EQ. LM) LP1 = LM
         LP2 = LP1
      ENDIF
C
C
      TEM1   = C1B(I)*W1I(IP,     LM1)    + C2B(I)*W1I(IP+I1  ,LM1)
     *       + C3B(I)*W1I(IP-I1,  LM1)    + C4B(I)*W1I(IP+I2  ,LM1)
      TEM2   = C1B(I)*W1I(IP+1,   LM1)    + C2B(I)*W1I(IP+1+I1,LM1)
     *       + C3B(I)*W1I(IP+1-I1,LM1)    + C4B(I)*W1I(IP+1+I2,LM1)
      TEM3   = C1B(I)*W1I(IP-1,   LM1)    + C2B(I)*W1I(IP-1+I1,LM1)
     *       + C3B(I)*W1I(IP-1-I1,LM1)    + C4B(I)*W1I(IP-1+I2,LM1)
      TEM4   = C1B(I)*W1I(IP+2   ,LM1)    + C2B(I)*W1I(IP+2+I1,LM1)
     *       + C3B(I)*W1I(IP+2-I1,LM1)    + C4B(I)*W1I(IP+2+I2,LM1)
C
      WK1(I)   = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
      TEM1   = C1B(I)*W1I(IP,     LP)    + C2B(I)*W1I(IP+I1  ,LP)
     *       + C3B(I)*W1I(IP-I1,  LP)    + C4B(I)*W1I(IP+I2  ,LP)
      TEM2   = C1B(I)*W1I(IP+1,   LP)    + C2B(I)*W1I(IP+1+I1,LP)
     *       + C3B(I)*W1I(IP+1-I1,LP)    + C4B(I)*W1I(IP+1+I2,LP)
      TEM3   = C1B(I)*W1I(IP-1,   LP)    + C2B(I)*W1I(IP-1+I1,LP)
     *       + C3B(I)*W1I(IP-1-I1,LP)    + C4B(I)*W1I(IP-1+I2,LP)
      TEM4   = C1B(I)*W1I(IP+2   ,LP)    + C2B(I)*W1I(IP+2+I1,LP)
     *       + C3B(I)*W1I(IP+2-I1,LP)    + C4B(I)*W1I(IP+2+I2,LP)
C
      WK2(I) = C1A(I) * TEM1 + C2A(I) * TEM2
     *       + C3A(I) * TEM3 + C4A(I) * TEM4
C
      TEM1   = C1B(I)*W1I(IP,     LP1)   + C2B(I)*W1I(IP+I1  ,LP1)
     *       + C3B(I)*W1I(IP-I1,  LP1)   + C4B(I)*W1I(IP+I2  ,LP1)
      TEM2   = C1B(I)*W1I(IP+1,   LP1)   + C2B(I)*W1I(IP+1+I1,LP1)
     *       + C3B(I)*W1I(IP+1-I1,LP1)   + C4B(I)*W1I(IP+1+I2,LP1)
      TEM3   = C1B(I)*W1I(IP-1,   LP1)   + C2B(I)*W1I(IP-1+I1,LP1)
     *       + C3B(I)*W1I(IP-1-I1,LP1)   + C4B(I)*W1I(IP-1+I2,LP1)
      TEM4   = C1B(I)*W1I(IP+2   ,LP1)   + C2B(I)*W1I(IP+2+I1,LP1)
     *       + C3B(I)*W1I(IP+2-I1,LP1)   + C4B(I)*W1I(IP+2+I2,LP1)
C
      WK3(I) = C1A(I) * TEM1 + C2A(I) * TEM2
     *       + C3A(I) * TEM3 + C4A(I) * TEM4
C
C
      TEM1   = C1B(I)*W1I(IP,     LP2)   + C2B(I)*W1I(IP+I1  ,LP2)
     *       + C3B(I)*W1I(IP-I1,  LP2)   + C4B(I)*W1I(IP+I2  ,LP2)
      TEM2   = C1B(I)*W1I(IP+1,   LP2)   + C2B(I)*W1I(IP+1+I1,LP2)
     *       + C3B(I)*W1I(IP+1-I1,LP2)   + C4B(I)*W1I(IP+1+I2,LP2)
      TEM3   = C1B(I)*W1I(IP-1,   LP2)   + C2B(I)*W1I(IP-1+I1,LP2)
     *       + C3B(I)*W1I(IP-1-I1,LP2)   + C4B(I)*W1I(IP-1+I2,LP2)
      TEM4   = C1B(I)*W1I(IP+2   ,LP2)   + C2B(I)*W1I(IP+2+I1,LP2)
     *       + C3B(I)*W1I(IP+2-I1,LP2)   + C4B(I)*W1I(IP+2+I2,LP2)
C
      WK4(I  ) = C1A(I) * TEM1 + C2A(I) * TEM2
     *         + C3A(I) * TEM3 + C4A(I) * TEM4
C
C
   70 CONTINUE
C
C
      DO 90 I=IJB,IJE
      LV = KPR(I)
      IF (LV .GE. 2 .AND. LV .LE. LMM2) THEN

c     LVB = KPR(I)
c     IF (LVB .GE. 2 .AND. LVB .LE. LMM2) THEN
c        LV = LVB
c     ELSEIF (LVB .EQ. 1) THEN
c        LV = 2
c     ELSE
c        LV = LMM2
c     ENDIF
c     IF (LVB .GE. 1 .AND. LVB .LE. LMM1) THEN

         TEM  = SGI(I)
         TEMA = (TEM-SGB(LV+1))*(TEM-SGB(LV+2))
         TEMB = (TEM-SGB(LV-1))*(TEM-SGB(LV))
C
         TEM1 = (TEM-SGB(LV))   * TEMA * S1A(LV,1)
         TEM2 = (TEM-SGB(LV-1)) * TEMA * S1A(LV,2)
         TEM3 = (TEM-SGB(LV+2)) * TEMB * S1A(LV,3)
         TEM4 = (TEM-SGB(LV+1)) * TEMB * S1A(LV,4)
C
C        TEM1 = (TEM-SGB(LV))  *(TEM-SGB(LV+1))*(TEM-SGB(LV+2))*S1A(LV)
C        TEM2 = (TEM-SGB(LV-1))*(TEM-SGB(LV+1))*(TEM-SGB(LV+2))*S1B(LV)
C        TEM3 = (TEM-SGB(LV-1))*(TEM-SGB(LV))  *(TEM-SGB(LV+2))*S1C(LV)
C        TEM4 = (TEM-SGB(LV-1))*(TEM-SGB(LV))  *(TEM-SGB(LV+1))*S1D(LV)
C
         W1O(I) = TEM1*WK1(I) + TEM2*WK2(I) + TEM3*WK3(I) +TEM4*WK4(I)
      ELSEif (lv .eq. 1 .or. lv .eq. lmm1) then
         TEM = (SGI(I) - SGB(LV)) * SGX(LV)
         W1O(I) = (1.0 - TEM) * WK2(I) + TEM * WK3(I)
      else
         W1O(I) = WK2(I)
      ENDIF
   90 CONTINUE
      CALL WRAP(IJB,IJE,W1O)
  100 CONTINUE
C
      RETURN
      END

