      PROGRAM CCMI_PCALC
C     PROGRAM CCMI_PCALC (TEMP,  SRFP,  XOUT1,  INPUT,  OUTPUT,)
C    1         TAPE1=TEMP,TAPE2=LNSP,TAPE11=XOUT1, TAPE12=XOUT2, TAPE5=INPUT,
C              TAPE6=OUTPUT)
C     ------------------------------------------------------------------
C
C     Aug 13/15 - D. Plummer - derived from gridcalc_eta
C     Aug 11/11 - D. Plummer (DERIVED FROM VSINTH_TROP)
C
C        - CALCULATES THE ATMOSPHERIC PRESSURE AT THE LAYER CENTRE OR LAYER
C          FACE
C
CINPUT FILE...
C
C      TEMP = AT LEAST ONE TIMESTEP OF A 3-D FIELD USED TO DERIVE THE
C               VERTICAL PLACEMENT OF THE LAYERS
C      SRFP = INPUT SERIES OF SURFACE PRESSURE IN MB
C                -- CURRENTLY ASSUMING PS IS LABELLED AS ' EXP'
C
COUTPUT FILE...
C
C      XOUT1 = PRESSURE AT LAYER CENTER OR LAYER FACE DEPENDING ON
C               OPTION SPECIFIED
C
CINPUT PARAMETERS...
C                                                                               J5
C      LEVTYP = FULL FOR MOMENTUM VARIABLE, AND                                 J5
C               HALF FOR THERMODYNAMIC ONE.                                     J5
C      LAY    = DEFINES THE POSITION OF LAYER INTERFACES IN RELATION            J5
C               TO LAYER CENTRES (SEE BASCAL).                                  J5
C               (ZERO DEFAULTS TO THE FORMER STAGGERING CONVENTION).            J5
C      ICOORD = 4H ETA/4H SIG FOR ETA/SIGMA LEVELS OF INPUT FIELDS.             J5
C      SIGTOP = VALUE OF SIGMA AT TOP OF DOMAIN FOR VERTICAL INTEGRAL.          J5
C               IF .LT.0., THEN INTERNALLY DEFINED BASED ON LEVTYP              J5
C               FOR UPWARD COMPATIBILITY.                                       J5
C      PTOIT  = PRESSURE (PA) OF THE RIGID LID OF THE MODEL.                    J5
C                                                                               J5
C      NOTE - LAY AND LEVTYP DEFINE THE TYPE OF LEVELLING FOR THE VARIABLE.     J5
C                                                                               J5
CEXAMPLE OF INPUT CARD...                                                       J5
C                                                                               J5
C*ACCGRID. HALF        1.    0  SIG  -1.        0.                              J5
C------------------------------------------------------------------------------
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
C
      PARAMETER (JPSIZ=65341, JPSIZB=130682, MXLEV=100)
C
      LOGICAL OK,SPEC
      INTEGER LEV(MXLEV)
C
      REAL ETA(MXLEV), ETAB(MXLEV)
      REAL ABM(MXLEV), BBM (MXLEV)
      REAL ABB(MXLEV), BBB (MXLEV)
C
      REAL SIGB(JPSIZ)
      REAL PSMBLN(JPSIZ), PRESS(JPSIZ)
C
      COMMON /ICOM/IBUF(8),IDAT(JPSIZB)
      DATA MAXX/JPSIZB/, MAXLEV/MXLEV/
C---------------------------------------------------------------------
      NFIL = 6
      CALL JCLPNT(NFIL,1,2,11,12,5,6)
      REWIND 1
      REWIND 2

C     * READ CONTROL DIRECTIVES.

      READ(5,5010,END=991) LEVTYP,LAY,ICOORD,IPOS,SIGTOP,PTOIT
      IF(ICOORD.EQ.NC4TO8("    ")) ICOORD=NC4TO8(" SIG")
      IF(ICOORD.EQ.NC4TO8(" SIG")) THEN
        PTOIT=MAX(PTOIT,0.00)
      ELSE
        PTOIT=MAX(PTOIT, 0.100E-09)
      ENDIF
C
      ICALC=0
      IF(IPOS.EQ.NC4TO8("MIDP")) ICALC=1
      IF(IPOS.EQ.NC4TO8("FACE")) ICALC=2
      IF(ICALC.EQ.0) CALL                          XIT('PCALC',-1)

      WRITE(6,6005) LEVTYP,LAY,ICOORD,IPOS,SIGTOP,PTOIT

C     * GET ETA VALUES OF LAYER CENTRES FROM XIN.

      CALL FILEV (LEV,ILEV,IBUF,1)
      IF(ILEV.LT.1 .OR. ILEV.GT.MAXLEV)  CALL      XIT('PCALC',-2)
      WRITE(6,6007) ILEV,(LEV(L),L=1,ILEV)

      CALL LVDCODE(ETA,LEV,ILEV)
      DO 150 L=1,ILEV
         ETA(L)=ETA(L)*.001
  150 CONTINUE

C     * DETERMINE THE FIELD SIZE.

      NWDS=IBUF(5)*IBUF(6)

C     * EVALUATE LAYER INTERFACES FROM LEVTYP AND LAY.
      IF    (LEVTYP.EQ.NC4TO8("FULL")) THEN
         CALL BASCAL (ETAB,IBUF, ETA,ETA,ILEV,LAY)
      ELSEIF(LEVTYP.EQ.NC4TO8("HALF")) THEN
         CALL BASCAL (IBUF,ETAB, ETA,ETA,ILEV,LAY)
      ELSE
         CALL                                   XIT('GRIDCALC',-3)
      ENDIF

C     * EVALUATE THE PARAMETERS A AND B OF THE VERTICAL DISCRETIZATION
C     * PERFORMED FOR BOTH THE LAYER MIDPOINTS AND LOWER FACE

      CALL COORDAB (ABM,BBM, ILEV,ETA ,ICOORD,PTOIT)
      CALL COORDAB (ABB,BBB, ILEV,ETAB,ICOORD,PTOIT)

C---------------------------------------------------------------------

      NSETS=0
  200 CONTINUE

C     * READ NEXT PS (HPA) NEEDED TO EVALUATE THE SIGMA LEVELS.

c     CALL GETFLD2 (2,PSMBLN,-1,-1,NC4TO8(" EXP"),1,IBUF,MAXX,OK)
      CALL GETFLD2 (2,PSMBLN,-1,-1,-1,1,IBUF,MAXX,OK)
      IF(NSETS.EQ.0) WRITE(6,6025) IBUF
      IF(.NOT.OK) THEN
         IF(NSETS.GT.0) THEN
            WRITE(6,6025) IBUF
            WRITE(6,6020) NSETS
            CALL                                XIT('GRIDCALC',0)
         ELSE
            CALL                                XIT('GRIDCALC',-4)
         ENDIF
      ENDIF

      ITIM=IBUF(2)
      DO 300 I=1,NWDS
         PRESS(I)=100.0*PSMBLN(I)
  300 CONTINUE

C     * CALCULATE THE GEOMETRIC THICKNESS AND MASS THICKNESS OF EACH LAYER

      IGH = 1
      LOW = 2

      IF(SIGTOP.LT.0.E0)THEN
         IF(LEVTYP.EQ.NC4TO8("FULL"))THEN
            DO 440 I=1,NWDS
               SIGB(I)=0.0E0
  440       CONTINUE
         ELSE
            DO 450 I=1,NWDS
               SIGB(I)=PTOIT/PRESS(I)
  450       CONTINUE
         ENDIF
      ELSE
         DO 460 I=1,NWDS
            SIGB(I)=SIGTOP
  460    CONTINUE
      ENDIF
C
C     * CALCULATE PRESSURE AT THE MID-POINT (ICALC=1) OR 
C     * THE TOP FACE OF EACH LAYER
C
      IF (ICALC.EQ.1) THEN

        DO L=1,ILEV
           CALL NIVCAL  (SIGB, ABM(L),BBM(L),PRESS,1,NWDS,NWDS)

           DO I=1,NWDS
             SIGB(I) = SIGB(I)*PRESS(I)
           ENDDO

           IBUF(3)=NC4TO8("PMID")
           IBUF(4)=LEV(L)
           CALL PUTFLD2 (11,SIGB,IBUF,MAXX)
        ENDDO

      ELSE

        DO I=1,NWDS
          SIGB(I) = SIGB(I)*PRESS(I)
        ENDDO

        IBUF(3)=NC4TO8("PFCE")
        IBUF(4)=LEV(1)
        CALL PUTFLD2 (11,SIGB,IBUF,MAXX)

        DO L=1,ILEV-1
           CALL NIVCAL  (SIGB, ABB(L),BBB(L),PRESS,1,NWDS,NWDS)

           DO I=1,NWDS
             SIGB(I) = SIGB(I)*PRESS(I)
           ENDDO

           IBUF(3)=NC4TO8("PFCE")
           IBUF(4)=LEV(L+1)
           CALL PUTFLD2 (11,SIGB,IBUF,MAXX)
        ENDDO

      ENDIF

      NSETS=NSETS+1

      GO TO 200

C     * E.O.F. ON INPUT.

  991 CALL                                      XIT('GRIDCALC',-6)
C---------------------------------------------------------------------
 5010 FORMAT(11X,A4,I5,1X,A4,1X,A4,E5.0,E10.0)
 6005 FORMAT('  LEVTYP = ',A4,', LAY = ',0P,I5, 
     1       ', COORD=',1X,A4,', IPOS =',1X,A4,' SIGTOP =',F15.10,
     2       ', P.LID (PA)=',E10.3)
 6007 FORMAT(' ILEV,ETA =',I5,2X,20I5/(18X,20I5))
 6020 FORMAT('0',I6,' SETS WERE PROCESSED')
 6025 FORMAT(' ',A4,I10,2X,A4,5I10)
      END
