      SUBROUTINE TRIGML4(NLATH,ILAT,SR,WR,CR,RADR,WOSQ,RADBU,RADBD,
     1                   SRX,WRX,CRX,RADRX,WOSQX,RADBXU,RADBXD,PI)
C
C     * MAR 22/19 - J. COLE    ADD CODE TO COMPUTE THE LATITUDE BOUNDARIES.
C     * JUN 12/06 - M. LAZARE. NEW VERSION FOR GCM15F:
C     *                        - PASSES IN PI AND REMOVES
C     *                          "PARAMS" COMMON BLOCK.
C     * JUN 30/03 - M. LAZARE. PREVIOUS VERSION TRIGML3:
C     *                        LIKE PREVIOUS VERSION TRIGML2, EXCEPT
C     *                        RE-ORDERED CHAINING TO ENSURE OPTIMAL
C     *                        LOAD-BALANCING FOR ALL PHYSICS
C     *                        PROCESSES.
C     *                        ALSO, ONLY GENERATES FIELDS FOR A GIVEN
C     *                        NODE, BASED ON "MPINFO".
C
C     * THIS ROUTINE IS BASED ON TRIGL, EXCEPT THAT THE TRIGNOMETRIC
C     * ARRAYS ARE RE-ORDERED IN CHAINED S-N PAIRS, BASED ON THE SYMMETRY
C     * PROPERTIES OF THE PARTICULAR TRIGNOMETRIC FUNCTION, FOR USE
C     * IN A MULTIPLE-LATITUDE MODEL SIMULATION. LATITUDES ARE ALTERNATED
C     * BETWEEN POLES AND EQUATOR TO ENSURE AVERAGE OPTIMAL LOAD-BALANCING
C     * FOR ALL PHYSICS PROCESSES.
C
C     * THE ORDERING IS THUS: (1,ILAT,ILATH,ILATH+1,2,ILAT-1,...).

C     * THE ROUTINE GAUSSG FILLS ONLY THE N HEM ORDERED N TO S AND IS
C     * CALLED IMMEDIATELY PRIOR TO THIS ROUTINE. BE AWARE THAT ITS
C     * OUTPUT FOR "RADR" DEFINES CO-LATITUDE; THIS ROUTINE CONVERTS
C     * THAT TO LATITUDE.
C
C     *      SR=SIN(LAT),  CR=COS(LAT),  RADR=LATITUDE IN RADIANS.
C     *      WR = GAUSSIAN WEIGHTS,  WOSQ = WR/(SR**2).
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
C
      REAL*8 SRX(1),WRX(1),CRX(1),RADRX(1),WOSQX(1)
      REAL*8 SR (1),WR (1),CR (1),RADR (1),WOSQ (1)
      REAL*8 RADBU(1), RADBD(1), RADBXU(1), RADBXD(1)
C
      INTEGER*4 MYNODE
C
      COMMON /MPINFO/ MYNODE
C--------------------------------------------------------------------
      PIH=PI/2.
      NLAT=NLATH*2
      IF(MOD(NLATH,2).NE.0)       CALL XIT('TRIGML4',-1)
      NLATQ=NLATH/2
C
C     * MPI HOOK.
C
      JSTART=MYNODE*ILAT+1
      JEND  =MYNODE*ILAT+ILAT
C
C     * FILL OUPUT ARRAYS BASED ON SYMMETRY ARGUEMENTS.
C     * CR,WR,WOSQ ARE SYMMETRIC ABOUT THE EQUATOR.
C     * SR AND RADR ARE ANTISYMMETRIC.
C     * NOTE THAT RADR IS CONVERTED FROM CO-LATITUDE TO LATITUDE HERE!
C     * NOTE ALSO THAT ONLY THE VALUES RELEVANT TO THE NODE ARE OBTAINED!
C     * "J" IS THE COUNTER INDEX FOR THE WHOLE SET OF GAUSSIAN LATITUDES
C     * WHILE "NJ" IS FOR THE PARTICULAR NODE ONLY.
C
      NJ=0
      J=0
      DO 150 JQ=1,NLATQ
        J0      = NLATH-JQ+1
C
        J=J+1
        IF(J.GE.JSTART .AND. J.LE.JEND) THEN
          NJ      = NJ+1
          SR(NJ)  =-SRX(JQ)
          CR(NJ)  = CRX(JQ)
          WR(NJ)  = WRX(JQ)
          WOSQ(NJ)= WOSQX(JQ)
          RADR(NJ)= RADRX(JQ)-PIH
          RADBU(NJ) = -RADBXD(JQ)
          RADBD(NJ) = -RADBXU(JQ)
        ENDIF
C
        J=J+1
        IF(J.GE.JSTART .AND. J.LE.JEND) THEN
          NJ      = NJ+1
          SR(NJ)  = SRX(JQ)
          CR(NJ)  = CRX(JQ)
          WR(NJ)  = WRX(JQ)
          WOSQ(NJ)= WOSQX(JQ)
          RADR(NJ)= PIH-RADRX(JQ)
          RADBU(NJ) = RADBXU(JQ)
          RADBD(NJ) = RADBXD(JQ)
        ENDIF
C
        J=J+1
        IF(J.GE.JSTART .AND. J.LE.JEND) THEN
          NJ      = NJ+1
          SR(NJ)  =-SRX(J0)
          CR(NJ)  = CRX(J0)
          WR(NJ)  = WRX(J0)
          WOSQ(NJ)= WOSQX(J0)
          RADR(NJ)= RADRX(J0)-PIH
          RADBU(NJ) = -RADBXD(J0)
          RADBD(NJ) = -RADBXU(J0)
        ENDIF
C
        J=J+1
        IF(J.GE.JSTART .AND. J.LE.JEND) THEN
          NJ      = NJ+1
          SR(NJ)  = SRX(J0)
          CR(NJ)  = CRX(J0)
          WR(NJ)  = WRX(J0)
          WOSQ(NJ)= WOSQX(J0)
          RADR(NJ)= PIH-RADRX(J0)
          RADBU(NJ) = RADBXU(J0)
          RADBD(NJ) = RADBXD(J0)
        ENDIF
  150 CONTINUE
C
      RETURN
      END
