      SUBROUTINE GEOSTR(U, V, T, SGD, DIV, GTP, LGP, DT, TX1, TX2)
Status: R

C
#include "slgcm.slgprm"
#include "slgcm.slgcon"
#include "slgcm.slgclm"
C
      DIMENSION   U(IJM0M0,LM),   V(IJM0M0,LM),   T(IJM0M0,LM)
      DIMENSION SGD(IJM0M0,LM), DIV(IJM0M0,LM), GTP(IJM0M0,LM)
      DIMENSION TX1(IJM0M0),    TX2(IJM0M0)
      DIMENSION LGP(IJM0M0)
      REAL LGP
      INTEGER DT
      DATA U0/30.0/, DTDZ/0.006/
C
      CALL CONSL(DT)
C
      DO 20 J=2,JMM1
      II = (J-1)*IMP3
      DO 10 I=1,IM
      TX1(I+II) = - U0 * SLM(I,1)*SPH(J,2)*(4.0*CPH(J,2)*CPH(J,2)-1.0)
   10 CONTINUE
   20 CONTINUE
      DO 40 J=1,JMM1
      II = (J-1)*IMP3
      DO 30 I=1,IM
      TX2(I+II) = U0 * CLM(I,2)*SPH(J,1)*SPH(J,1)
   30 CONTINUE
   40 CONTINUE
      DO 47 J=1,JM
      II = (J-1)*IMP3
      DO 45 I=1,IM
      DIV(I+II,1) = TWOM*AE*U0*SPH(J,2)**3*CPH(J,2)*SLM(I,2)
   45 CONTINUE
   47 CONTINUE
C
      CALL WRAP(I1,IJM1M0,TX1)
      CALL FILPU(TX1)
      CALL WRAP(1,IJM1M0,TX2)
      CALL FILPV(TX2)
      CALL WRAP(1,IJM0M0,DIV(1,1))
C
      DO 70 I=1,IJM0M0
      LGP(I) = ALOG(100000.0)
   70 CONTINUE
C
      TEM = 0.0
      DO 100 L=LM,1,-1
C     T0 = 270.0 * EXP(RGAS*DTDZ/GRAV)*SGB(L)
      TEM = TEM + RGAS * TBR(L) * DSL(L)
      DO 50 I=1,IJM0M0
      U(I,L)   = TX1(I)
      V(I,L)   = TX2(I)
      T(I,L)   = TBR(L)
      GTP(I,L) = DIV(I,1) + RGAS * TBR(L) * LGP(I) + TEM
      SGD(I,L) = ZERO
      DIV(I,L) = ZERO
   50 CONTINUE
  100 CONTINUE
C
      RETURN
      END
C
      SUBROUTINE GEOPOL(U, V, T, Q, SGD, DIV, GTP, LGP, PIB, TSF)
C
#include "slgcm.slgprm"
#include "slgcm.slgcon"
#include "slgcm.slgclm"
C
      DIMENSION   U(IJM0M0,LM),   V(IJM0M0,LM),   T(IJM0M0,LM)
     *,           Q(IJM0M0,LM),   PIB(IJM0M0),    TSF(IJM0M0)
      DIMENSION SGD(IJM0M0,LM), DIV(IJM0M0,LM), GTP(IJM0M0,LM)
      DIMENSION LGP(IJM0M0)
      REAL LGP
      DATA V0B/00.0/, DVDSG/100.0/
C
C PHI FIELD
C
      PIKAPI  = 1.0 / (984.7 ** RKAP)
      LGP (:) = 0.0
      PIB(:)  = 984.7
      TSF(:)  = TBR(LM)
      Q(:,:)  = 0.010
      PBR     = ALOG(98470.0)
      DO 10 J=1,JM
      II = (J-1)*IMP3
      TEM = TWOM*AE*V0B*SPH(J,2)*SPH(J,2)*SPH(J,2)*CPH(J,2)
      DO 5 I=1,IM
      PHIS(I+II) = TEM * SLM(I,2)
    5 CONTINUE
   10 CONTINUE
      CALL WRAP(1,IJM0M0,PHIS)
      GTP(:,2) = PHIS(:)
C
      U(:,:) = 0.0
      V(:,:) = 0.0
      T(:,:) = 0.0
C
      DO 60 L=LM,1,-1
C
      V0 = V0B + (1.0 - SGB(L)) * DVDSG
      print *,' l=',l,' V0=',V0
C
C U FIELD
C
      DO 20 J=2,JMM1
      II = (J-1)*IMP3
      TEM = V0 * SPH(J,2) * (4.0*CPH(J,2)*CPH(J,2) - 1.0)
C     TEM1 = TWOM * AE * CPH(J,2) * 0.5
C     TEM2 = TWOM * AE * V0 * SPH(J,2) * CPH(J,2)
C    *            * (4.0*CPH(J,2)*CPH(J,2) - 1.0)
C     PRINT *,' TEM2=',TEM2,' AE=',AE,' V0=',V0,' SPH=',SPH(J,2)
C    *,     ' CPH=',CPH(J,2)
      DO 15 I=1,IM
C     TEM3 =  TEM1*TEM1 - TEM2 * SLM(I,1)
C     IF (TEM3 .LE. 0.0) THEN
C        PRINT *,' I=',I,' J=',J,' TEM3=',TEM3,' TEM1=',TEM1,' TEM2='
C    *,   TEM2,' SLM=',SLM(I,1)
C     ENDIF

C     TEM3 =  SQRT(TEM3)
C     U1   = -TEM1 + TEM3
C     U2   = -TEM1 - TEM3
C     U(I+II,1) = U1
      U(I+II,L) = - TEM * SLM(I,1)
   15 CONTINUE
   20 CONTINUE
      CALL WRAP(1,IJM0M0,U(1,L))
C
C V FIELD
C
      DO 30 J=1,JMM1
      II = (J-1)*IMP3
      TEM = V0 * SPH(J,1) * SPH(J,1)
C     TEM1 = TWOM * AE * CPH(J,1) * 0.5
C     TEM2 = TWOM * AE * V0 * SPH(J,1) * CPH(J,1)
C    *            * (4.0*CPH(J,1)*CPH(J,1) - 1.0)
C     TEM4 = TWOM * AE * V0 * SPH(J,1) * SPH(J,1) * CPH(J,1)
      DO 25 I=1,IM
C     TEM3 =  TEM1*TEM1 - TEM2 * SLM(I,2)
C     IF (TEM3 .LE. 0.0) THEN
C        PRINT *,' I=',I,' J=',J,' TEM3=',TEM3,' TEM1=',TEM1,' TEM2='
C    *,   TEM2,' SLM2=',SLM(I,2)
C     ENDIF
C     TEM3 =  SQRT(TEM3)
C     U1   = -TEM1 + TEM3
C     U2   = -TEM1 - TEM3
C     V(I+II,1) = TEM4 * CLM(I,2) / (TWOM*AE*CPH(J,1) + U1)
      V(I+II,L) = TEM * CLM(I,2)
   25 CONTINUE
   30 CONTINUE
      CALL WRAP(1,IJM0M0,V(1,L))
C
      V00 = DSG(L) * DVDSG / (RGAS*DSL(L))
      DO 35 J=1,JM
      II = (J-1)*IMP3
      TEM = TWOM*AE*V00*SPH(J,2)*SPH(J,2)*SPH(J,2)*CPH(J,2)
c*    TEM = TWOM*AE*V0*SPH(J,2)*SPH(J,2)*SPH(J,2)*CPH(J,2)
      DO 35 I=1,IM
c     GTP(I+II,1) = TEM * SLM(I,2)
      T(I+II,L) = (TBR(L) + TEM * SLM(I,2)) * (PIKAPI/PRH(L))
   35 CONTINUE
      CALL WRAP (1,IJM0M0, T(1,L))
c     CALL WRAP (1,IJM0M0, GTP(1,1))
c     DO 37 I=1,IJM0M0
c     GTP(I,1) = 2.0 * GTP(I,1) - GTP(I,2)
c     T(I,L) = TBR(L) + (GTP(I,1) - GTP(I,2)) / (RGAS*DSL(L))
c     T(I,L) = (GTP(I,1) - GTP(I,2)) / (RGAS*DSL(L)*TBR(L))
c     GTP(I,2) = GTP(I,1)
c  37 CONTINUE
   60 CONTINUE
C
C     DO 40 L=2,LM
C     U(:,L) = U(:,1)
C     V(:,L) = V(:,1)
C  40 CONTINUE
C
      TEMS = 0.0
      DO 100 L=LM,1,-1
      TEM = RGAS * TBR(L) * DSL(L)
      TEMA = TEMS + TEM * 0.5
      TEMS = TEMS + TEM
C
C V, T, G, SIGDOT, DIV FIELDS
C
      IF (MOD(LM-L,2) .EQ. 0) THEN
         TEM   =  (2.0 / (RGAS*DSL(L)*TBR(L)))
      ELSE
         TEM   = - (2.0 / (RGAS*DSL(L)*TBR(L)))
      ENDIF

      DO 50 I=1,IJM0M0
C     T(I,L)   = 0.0
C     T(I,L)   = PHIS(I) * TEM
C     GTP(I,L) = PHIS(I) + RGAS * TBR(L) * LGP(I) + TEMA
      GTP(I,L) = 0.0
      SGD(I,L) = 0.0
      DIV(I,L) = 0.0
   50 CONTINUE
  100 CONTINUE
c     PHIS(:) = 0.0
C
      RETURN
      END
      SUBROUTINE GEOC(U, V, T, SGD, DIV, GTP, LGP, DT)
C
#include "slgcm.slgprm"
#include "slgcm.slgcon"
C
      DIMENSION   U(IJM0M0,LM),   V(IJM0M0,LM),   T(IJM0M0,LM)
      DIMENSION SGD(IJM0M0,LM), DIV(IJM0M0,LM), GTP(IJM0M0,LM)
      DIMENSION LGP(IJM0M0)
      REAL LGP
      INTEGER DT
      DATA U0/10.0/
C
C PRESET ARRAYS
C
      CALL CONSL(DT)
C
C PHI FIELD
C
      DO 47 J=1,JM
      II = (J-1)*IMP3
      TEM = -(TWOM*AE+U0)*U0*SPH(J,2)*SPH(J,2)*0.5
      DO 45 I=1,IMP3
      PHIS(I+II) = 0.0
      LGP(I+II)  = TEM + 100000.0
   45 CONTINUE
   47 CONTINUE
C
      DO 70 I=1,IJM0M0
      LGP(I) = ALOG(LGP(I))
   70 CONTINUE
C
      TEMS = 0.0
      DO 100 L=LM,1,-1
      TEM = RGAS * TBR(L) * DSL(L)
      TEMA = TEMS + TEM * 0.5
      TEMS = TEMS + TEM
C
C U FIELD
C
      DO 40 J=1,JM
      II = (J-1) * IMP3
      DO 40 I=1,IMP3
      U(II+I,L) = U0 * CPH(J,2)
   40 CONTINUE
C
C V, T, G, SIGDOT, DIV FIELDS
C
      DO 50 I=1,IJM0M0
      V(I,L)   = 0.0
      T(I,L)   = 0.0
      GTP(I,L) = PHIS(I) + RGAS * TBR(L) * LGP(I) + TEMA
      SGD(I,L) = ZERO
      DIV(I,L) = ZERO
   50 CONTINUE
  100 CONTINUE
C
      RETURN
      END
      SUBROUTINE GEOCN(U, V, T, SGD, DIV, GTP, LGP, DT)
C
#include "slgcm.slgprm"
#include "slgcm.slgcon"
C
      DIMENSION   U(IJM0M0,LM),   V(IJM0M0,LM),   T(IJM0M0,LM)
      DIMENSION SGD(IJM0M0,LM), DIV(IJM0M0,LM), GTP(IJM0M0,LM)
      DIMENSION LGP(IJM0M0), U0(50), WK1(IMP3,JM)
      REAL LGP
      INTEGER DT
      DATA DUDSG/50.0/
C
C PRESET ARRAYS
C
      CALL CONSL(DT)
C
C PHI FIELD
C
      DO 10 L=1,LMP1
      U0(L) = DUDSG * (SIG(LMP1) - SIG(L))
   10 CONTINUE
C
C     READ (KBC) ((WK1(I,J),I=1,IMP1),J=1,JM)
C     CALL WRAP(1,IJM0M0,WK1(1,1))
      DO 12 I=1,IJM0M0
C     PHIS(I) = WK1(I,1)
      PHIS(I) = 0.0
   12 CONTINUE
C
      DO 20 L=1,LM
      DO 47 J=1,JM
      II = (J-1)*IMP3
      GTP(II+1,L) = -(TWOM*AE+U0(L))*U0(L)*SPH(J,2)*SPH(J,2)*0.5
   47 CONTINUE
   20 CONTINUE
      DO 25 L=1,LM
      IF (L .LT. LM) THEN
         DO 30 J=1,JM
         II = (J-1)*IMP3
         TEMA = (GTP(II+1,L+1) - GTP(II+1,L)) / (DSL(L)*RGAS*TBR(L))
         TEMB = (GTP(II+1,L+1) + GTP(II+1,L)) * 0.5
C
         TEMC = 0.0
         IF ((J.GT.10 .AND. J.LT.20) .OR. (J.GT.30 .AND. J.LT.40))
     *      TEMC = 0.02
C
         DO 35 I=1,IMP3
         IF (I .LE. IMB2+1) THEN
            T(I+II,L)   = - TEMA  + TEMC * I/FLOAT(IMB2+1)
         ELSE
            T(I+II,L)   = - TEMA  + TEMC * (IMP2-I)/FLOAT(IMB2+1)
         ENDIF
         GTP(I+II,L) = TEMB
   35    CONTINUE
   30    CONTINUE
         CALL WRAP(1,IJM0M0,T(1,L))
      ELSE
         DO 38 J=1,JM
         II = (J-1)*IMP3
         TEMA =  GTP(II+1,L) / (DSL(L)*RGAS*TBR(L))
         TEMB =  GTP(II+1,L) * 0.5
         DO 36 I=1,IMP3
         T(I+II,L)   = TEMA
         GTP(I+II,L) = TEMB
   36    CONTINUE
   38    CONTINUE
      ENDIF
   25 CONTINUE
C
      DO 70 I=1,IJM0M0
      LGP(I) = ALOG(100000.0)
   70 CONTINUE
C
      TEMS = 0.0
      DO 100 L=LM,1,-1
      TEM = RGAS * TBR(L) * DSL(L)
      TEMA = TEMS + TEM * 0.5
      TEMS = TEMS + TEM
C
C U FIELD
C
      DO 40 J=1,JM
      II = (J-1) * IMP3
      DO 40 I=1,IMP3
      U(II+I,L) = 0.5 * (U0(L)+U0(L+1)) * CPH(J,2)
   40 CONTINUE
C
C V, T, G, SIGDOT, DIV FIELDS
C
      DO 50 I=1,IJM0M0
      V(I,L)   = 0.0
      GTP(I,L) = PHIS(I) + RGAS * TBR(L) * LGP(I) + TEMA + GTP(I,L)
      SGD(I,L) = ZERO
      DIV(I,L) = ZERO
   50 CONTINUE
  100 CONTINUE
C
      RETURN
      END
      SUBROUTINE RHWTZ(U, V, T, SGD, DIV, GTP, LGP, DT)
C
#include "slgcm.slgprm"
#include "slgcm.slgcon"
C
      PARAMETER (OMG=1.0, RW=1.0, KH=00,  FRQ=0.0)
      DIMENSION   U(IJM0M0,LM),   V(IJM0M0,LM),   T(IJM0M0,LM)
      DIMENSION SGD(IJM0M0,LM), DIV(IJM0M0,LM), GTP(IJM0M0,LM)
      DIMENSION LGP(IJM0M0), U0(50), WK1(IMP3,JM)
      REAL LGP
      INTEGER DT

C
C PRESET ARRAYS
C
      CALL CONSL(DT)
C
C PHI FIELD
C
C
      READ (KBC) ((WK1(I,J),I=1,IMP1),J=1,JM)
      CALL WRAP(1,IJM0M0,WK1(1,1))
      DO 12 I=1,IJM0M0
C     PHIS(I) = WK1(I,1)
      PHIS(I) = 0.0
   12 CONTINUE
C
      DO 20 L=1,LM
      DO 47 J=1,JM
      II = (J-1)*IMP3
      GTP(II+1,L) = -(TWOM*AE+U0(L))*U0(L)*SPH(J,2)*SPH(J,2)*0.5
   47 CONTINUE
   20 CONTINUE
      DO 25 L=1,LM
      IF (L .LT. LM) THEN
         DO 30 J=1,JM
         II = (J-1)*IMP3
         TEMA = (GTP(II+1,L+1) - GTP(II+1,L)) / (DSL(L)*RGAS*TBR(L))
         TEMB = (GTP(II+1,L+1) + GTP(II+1,L)) * 0.5
C
         TEMC = 0.0
         IF ((J.GT.10 .AND. J.LT.20) .OR. (J.GT.30 .AND. J.LT.40))
     *      TEMC = 0.005
C
         DO 35 I=1,IMP3
         IF (I .LE. IMB2+1) THEN
            T(I+II,L)   = - TEMA  + TEMC * I/FLOAT(IMB2+1)
         ELSE
            T(I+II,L)   = - TEMA  + TEMC * (IMP2-I)/FLOAT(IMB2+1)
         ENDIF
         GTP(I+II,L) = TEMB
   35    CONTINUE
   30    CONTINUE
         CALL WRAP(1,IJM0M0,T(1,L))
      ELSE
         DO 38 J=1,JM
         II = (J-1)*IMP3
         TEMA =  GTP(II+1,L) / (DSL(L)*RGAS*TBR(L))
         TEMB =  GTP(II+1,L) * 0.5
         DO 36 I=1,IMP3
         T(I+II,L)   = TEMA
         GTP(I+II,L) = TEMB
   36    CONTINUE
   38    CONTINUE
      ENDIF
   25 CONTINUE
C
      DO 70 I=1,IJM0M0
      LGP(I) = ALOG(100000.0)
   70 CONTINUE
C
      TEMS = 0.0
      DO 100 L=LM,1,-1
      TEM = RGAS * TBR(L) * DSL(L)
      TEMA = TEMS + TEM * 0.5
      TEMS = TEMS + TEM
C
C U FIELD
C
      DO 40 J=1,JM
      II = (J-1) * IMP3
      DO 40 I=1,IMP3
      U(II+I,L) = 0.5 * (U0(L)+U0(L+1)) * CPH(J,2)
   40 CONTINUE
C
C V, T, G, SIGDOT, DIV FIELDS
C
      DO 50 I=1,IJM0M0
      V(I,L)   = 0.0
      GTP(I,L) = PHIS(I) + RGAS * TBR(L) * LGP(I) + TEMA + GTP(I,L)
      SGD(I,L) = ZERO
      DIV(I,L) = ZERO
   50 CONTINUE
  100 CONTINUE
C
      RETURN
      END
      SUBROUTINE GEOPO1(U, V, T, SGD, DIV, GTP, LGP, DT)
C
#include "slgcm.slgprm"
#include "slgcm.slgcon"
#include "slgcm.slgclm"
C
      DIMENSION   U(IJM0M0,LM),   V(IJM0M0,LM),   T(IJM0M0,LM)
      DIMENSION SGD(IJM0M0,LM), DIV(IJM0M0,LM), GTP(IJM0M0,LM)
      DIMENSION LGP(IJM0M0)
      REAL LGP
      INTEGER DT
      DATA V0B/100.0/, DVDSG/00.0/
C
C PRESET ARRAYS
C
      CALL CONSL(DT)
C
C PHI FIELD
C
      cdrg(:) = 0.0
      LGP (:) = ALOG(100000.0)
      CALL GLMNP2(IMP3, JM, LGP     , DX(1,2), DY, ARE2, DX(1,1), PBR)
      LGP (:) = LGP(:) - PBR
      DO 10 J=1,JM
      II = (J-1)*IMP3
      TEM = TWOM*AE*V0B*SPH(J,2)*SPH(J,2)*SPH(J,2)*CPH(J,2)
      DO 5 I=1,IM
      PHIS(I+II) = TEM * SLM(I,2)
    5 CONTINUE
   10 CONTINUE
      CALL WRAP(1,IJM0M0,PHIS)
      GTP(:,2) = PHIS(:)
C
C U FIELD
C
      DO 60 L=LM,1,-1
      V0 = V0B + (1.0 - SGB(L)) * DVDSG
      DO 20 J=2,JMM1
      II = (J-1)*IMP3
      TEM = V0 * SPH(J,2) * (4.0*CPH(J,2)*CPH(J,2) - 1.0)
C     TEM1 = TWOM * AE * CPH(J,2) * 0.5
C     TEM2 = TWOM * AE * V0 * SPH(J,2) * CPH(J,2)
C    *            * (4.0*CPH(J,2)*CPH(J,2) - 1.0)
      DO 15 I=1,IM
C     TEM3 =  TEM1*TEM1 - TEM2 * SLM(I,1)

C     TEM3 =  SQRT(TEM3)
C     U1   = -TEM1 + TEM3
C     U2   = -TEM1 - TEM3
C     U(I+II,1) = U1
      U(I+II,L) = - TEM * SLM(I,1)
   15 CONTINUE
   20 CONTINUE
      CALL WRAP(1,IJM0M0,U(1,1))
C
C V FIELD
C
      DO 30 J=1,JMM1
      II = (J-1)*IMP3
      TEM = V0 * SPH(J,1) * SPH(J,1)
C     TEM1 = TWOM * AE * CPH(J,1) * 0.5
C     TEM2 = TWOM * AE * V0 * SPH(J,1) * CPH(J,1)
C    *            * (4.0*CPH(J,1)*CPH(J,1) - 1.0)
C     TEM4 = TWOM * AE * V0 * SPH(J,1) * SPH(J,1) * CPH(J,1)
      DO 25 I=1,IM
C     TEM3 =  TEM1*TEM1 - TEM2 * SLM(I,2)
C     TEM3 =  SQRT(TEM3)
C     U1   = -TEM1 + TEM3
C     U2   = -TEM1 - TEM3
C     V(I+II,1) = TEM4 * CLM(I,2) / (TWOM*AE*CPH(J,1) + U1)
      V(I+II,L) = TEM * CLM(I,2)
   25 CONTINUE
   30 CONTINUE
      CALL WRAP(1,IJM0M0,V(1,1))
C
      DO 35 J=1,JM
      II = (J-1)*IMP3
      TEM = TWOM*AE*V0*SPH(J,2)*SPH(J,2)*SPH(J,2)*CPH(J,2)
      DO 35 I=1,IM
      GTP(I+II,1) = TEM * SLM(I,2)
   35 CONTINUE
      CALL WRAP (1,IJM0M0, GTP(1,1))
      DO 37 I=1,IJM0M0
      GTP(I,1) = 2.0 * GTP(I,1) - GTP(I,2)
      T(I,L) = (GTP(I,1) - GTP(I,2)) / (RGAS*DSL(L)*TBR(L))
      GTP(I,2) = GTP(I,1)
   37 CONTINUE
   60 CONTINUE
C
C     DO 40 L=2,LM
C     U(:,L) = U(:,1)
C     V(:,L) = V(:,1)
C  40 CONTINUE
C
      TEMS = 0.0
      DO 100 L=LM,1,-1
      TEM = RGAS * TBR(L) * DSL(L)
      TEMA = TEMS + TEM * 0.5
      TEMS = TEMS + TEM
C
C V, T, G, SIGDOT, DIV FIELDS
C
      IF (MOD(LM-L,2) .EQ. 0) THEN
         TEM   =  (2.0 / (RGAS*DSL(L)*TBR(L)))
      ELSE
         TEM   = - (2.0 / (RGAS*DSL(L)*TBR(L)))
      ENDIF

      DO 50 I=1,IJM0M0
C     T(I,L)   = 0.0
C     T(I,L)   = PHIS(I) * TEM
C     GTP(I,L) = PHIS(I) + RGAS * TBR(L) * LGP(I) + TEMA
      GTP(I,L) = 0.0
      SGD(I,L) = 0.0
      DIV(I,L) = 0.0
   50 CONTINUE
  100 CONTINUE
C     PHIS(:) = 0.0
C
      RETURN
      END

