! The kind of code only a scientist could love.
! TODO: The cape routine needs work to remove the GOTOs

!======================================================================
!
! !IROUTINE: TVIRTUAL -- Calculate virtual temperature (K)
!
! !DESCRIPTION:
!
!   This function returns a single value of virtual temperature in
!   K, given temperature in K and mixing ratio in kg/kg.  For an
!   array of virtual temperatures, use subroutine VIRTUAL_TEMP.
!
! !INPUT:
!    RATMIX - water vapor mixing ratio (kg/kg)
!    TEMP   - temperature (K)
!
! !OUTPUT:
!    TV     - Virtual temperature (K)
!

REAL(KIND=8) FUNCTION tvirtual(temp,ratmix)
    IMPLICIT NONE
    REAL(KIND=8),INTENT(IN) :: temp,ratmix
    REAL(KIND=8),PARAMETER :: EPS = .622D0

    tvirtual = temp*(EPS + ratmix)/(EPS*(1.D0 + ratmix))
    RETURN
END FUNCTION tvirtual

REAL(KIND=8) FUNCTION tonpsadiabat(thte,prs,PSADITHTE,PSADIPRS,PSADITMK,GAMMA,&
                                    throw_exception)
    IMPLICIT NONE
    EXTERNAL throw_exception
    REAL(KIND=8),INTENT(IN) :: thte
    REAL(KIND=8),INTENT(IN) :: prs
    REAL(KIND=8),DIMENSION(150),INTENT(IN) :: PSADITHTE
    REAL(KIND=8),DIMENSION(150),INTENT(IN) :: PSADIPRS
    REAL(KIND=8),DIMENSION(150,150),INTENT(IN) :: PSADITMK
    REAL(KIND=8),INTENT(IN) :: GAMMA

    REAL(KIND=8) :: fracjt
    REAL(KIND=8) :: fracjt2
    REAL(KIND=8) :: fracip
    REAL(KIND=8) :: fracip2

    INTEGER :: ip, ipch, jt, jtch

    !   This function gives the temperature (in K) on a moist adiabat
    !   (specified by thte in K) given pressure in hPa.  It uses a
    !   lookup table, with data that was generated by the Bolton (1980)
    !   formula for theta_e.

    !     First check if pressure is less than min pressure in lookup table.
    !     If it is, assume parcel is so dry that the given theta-e value can
    !     be interpretted as theta, and get temperature from the simple dry
    !     theta formula.

    IF (prs.LE.PSADIPRS(150)) THEN
      tonpsadiabat = thte * (prs/1000.D0)**GAMMA
      RETURN
    END IF

    !   Otherwise, look for the given thte/prs point in the lookup table.

    jt = -1
    DO jtch = 1,150 - 1
      IF (thte .GE. PSADITHTE(jtch) .AND. thte .LT. PSADITHTE(jtch+1)) THEN
          jt = jtch
          EXIT
          !GO TO 213
      END IF
    END DO

!    JT = -1
!213 CONTINUE
    ip = -1
    DO ipch = 1,150 - 1
      IF (prs.LE.PSADIPRS(ipch) .AND. prs.GT.PSADIPRS(ipch+1)) THEN
          ip = ipch
          EXIT
          !GO TO 215
      END IF
    END DO

!    IP = -1
!215 CONTINUE
    IF (jt.EQ.-1 .OR. ip.EQ.-1) THEN
        ! Need an exception here
        CALL throw_exception('capecalc3d: ','Outside of lookup table bounds. prs,thte=',prs,thte)
        !STOP
    END IF
    fracjt = (thte-PSADITHTE(jt)) / (PSADITHTE(jt+1)-PSADITHTE(jt))
    fracjt2 = 1.D0 - fracjt
    fracip = (PSADIPRS(ip)-prs) / (PSADIPRS(ip)-PSADIPRS(ip+1))
    fracip2 = 1.D0 - fracip
    IF (PSADITMK(ip,jt) .GT. 1D9 .OR. PSADITMK(ip+1,jt) .GT. 1D9 .OR. &
        PSADITMK(ip,jt+1) .GT. 1D9 .OR. PSADITMK(ip+1,jt+1) .GT. 1D9) THEN
        CALL throw_exception('capecalc3d: ','Tried to access missing temperature in lookup table.',&
        'Prs and Thte probably unreasonable. prs,thte=',prs,thte)
      !STOP
    END IF
    tonpsadiabat = fracip2*fracjt2*PSADITMK(ip,jt)+fracip*fracjt2*PSADITMK(ip+1,jt)+&
            fracip2*fracjt*PSADITMK(ip,jt+1)+fracip*fracjt*PSADITMK(ip+1,jt+1)

    RETURN
END FUNCTION tonpsadiabat


!     Historically, this routine calculated the pressure at full sigma
!     levels when RIP was specifically designed for MM4/MM5 output.
!     With the new generalized RIP (Feb '02), this routine is still
!     intended to calculate a set of pressure levels that bound the
!     layers represented by the vertical grid points, although no such
!     layer boundaries are assumed to be defined.  The routine simply
!     uses the midpoint between the pressures of the vertical grid
!     points as the bounding levels.  The array only contains mkzh
!     levels, so the pressure of the top of the uppermost layer is
!     actually excluded.  The kth value of pf is the lower bounding
!     pressure for the layer represented by kth data level.  At the
!     lower bounding level of the lowest model layer, it uses the
!     surface pressure, unless the data set is pressure-level data, in
!     which case it assumes the lower bounding pressure level is as far
!     below the lowest vertical level as the upper bounding pressure
!     level is above.
SUBROUTINE dpfcalc(prs,sfp,pf,miy,mjx,mkzh,ter_follow)

    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: prs
    REAL(KIND=8),DIMENSION(miy,mjx),INTENT(IN) :: sfp
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(OUT) :: pf
    INTEGER,INTENT(IN) :: ter_follow,miy,mjx,mkzh

    INTEGER :: i,j,k

    !  do j=1,mjx-1  Artifact of MM5
    DO j = 1,mjx
    !  do i=1,miy-1  staggered grid
      DO i = 1,miy
          DO k = 1,mkzh
              IF (k .EQ. mkzh) THEN
    !  terrain-following data
                  IF (ter_follow .EQ. 1) THEN
                      pf(i,j,k) = sfp(i,j)
    !  pressure-level data
                  ELSE
                      pf(i,j,k) = .5D0 * (3.D0*prs(i,j,k) - prs(i,j,k-1))
                  END IF
              ELSE
                  pf(i,j,k) = .5D0 * (prs(i,j,k+1) + prs(i,j,k))
              END IF
          END DO
      END DO
    END DO

    RETURN
END SUBROUTINE dpfcalc

!======================================================================
!
! !IROUTINE: capecalc3d -- Calculate CAPE and CIN
!
! !DESCRIPTION:
!
!   If i3dflag=1, this routine calculates CAPE and CIN (in m**2/s**2,
!   or J/kg) for every grid point in the entire 3D domain (treating
!   each grid point as a parcel).  If i3dflag=0, then it
!   calculates CAPE and CIN only for the parcel with max theta-e in
!   the column, (i.e. something akin to Colman's MCAPE).  By "parcel",
!   we mean a 500-m deep parcel, with actual temperature and moisture
!   averaged over that depth.
!
!   In the case of i3dflag=0,
!   CAPE and CIN are 2D fields that are placed in the k=mkzh slabs of
!   the cape and cin arrays.  Also, if i3dflag=0, LCL and LFC heights
!   are put in the k=mkzh-1 and k=mkzh-2 slabs of the cin array.
!


! Important!  The z-indexes must be arranged so that mkzh (max z-index) is the
! surface pressure.  So, pressure must be ordered in ascending order before
! calling this routine.  Other variables must be ordered the same (p,tk,q,z).

! Also, be advised that missing data values are not checked during the computation.
! Also also, Pressure must be hPa
!
SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
            PSADITHTE,PSADIPRS,PSADITMK,cmsg,i3dflag,ter_follow,&
            throw_exception,miy,mjx,mkzh)

    IMPLICIT NONE
    EXTERNAL throw_exception
    INTEGER,INTENT(IN) :: miy,mjx,mkzh,i3dflag,ter_follow
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: prs
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: tmk
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: qvp
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: ght
    REAL(KIND=8),DIMENSION(miy,mjx),INTENT(IN) :: ter
    REAL(KIND=8),DIMENSION(miy,mjx),INTENT(IN) ::sfp
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(INOUT) :: cape
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(INOUT) :: cin
    REAL(KIND=8),DIMENSION(150),INTENT(IN) :: PSADITHTE,PSADIPRS
    REAL(KIND=8),DIMENSION(150,150),INTENT(IN) :: PSADITMK
    REAL(KIND=8),INTENT(IN) :: cmsg

    ! local variables
    INTEGER :: i,j,k,ilcl,kel,kk,klcl,klev,klfc,kmax,kpar,kpar1,kpar2
    REAL(KIND=8) :: davg,ethmax,q,t,p,e,eth,tlcl,zlcl
    REAL(KIND=8) :: pavg,tvirtual,p1,p2,pp1,pp2,th,totthe,totqvp,totprs
    REAL(KIND=8) :: cpm,deltap,ethpari,gammam,ghtpari,qvppari,prspari,tmkpari
    REAL(KIND=8) :: facden,fac1,fac2,qvplift,tmklift,tvenv,tvlift,ghtlift
    REAL(KIND=8) :: eslift,tmkenv,qvpenv,tonpsadiabat
    REAL(KIND=8) :: benamin,dz,pup,pdn
    REAL(KIND=8),DIMENSION(150) :: buoy,zrel,benaccum
    REAL(KIND=8),DIMENSION(miy,mjx,mkzh) :: prsf

    ! constants
    INTEGER,PARAMETER :: IUP = 6
    REAL(KIND=8),PARAMETER :: CELKEL = 273.15d0
    REAL(KIND=8),PARAMETER :: GRAV = 9.81d0
    ! hpa
    REAL(KIND=8),PARAMETER :: EZERO = 6.112d0
    REAL(KIND=8),PARAMETER :: ESLCON1 = 17.67d0
    REAL(KIND=8),PARAMETER :: ESLCON2 = 29.65d0
    REAL(KIND=8),PARAMETER :: EPS = 0.622d0
    ! j/k/kg
    REAL(KIND=8),PARAMETER :: RGAS = 287.04d0
    !  j/k/kg  note: not using bolton's value of 1005.7
    REAL(KIND=8),PARAMETER :: CP = 1004.d0
    REAL(KIND=8),PARAMETER :: GAMMA = RGAS/CP
    !  cp_moist=cp*(1.+cpmd*qvp)
    REAL(KIND=8),PARAMETER :: CPMD = .887d0
    !  rgas_moist=rgas*(1.+rgasmd*qvp)
    REAL(KIND=8),PARAMETER :: RGASMD = .608d0
    !  gamma_moist=gamma*(1.+gammamd*qvp)
    REAL(KIND=8),PARAMETER :: GAMMAMD = RGASMD - CPMD
    REAL(KIND=8),PARAMETER :: TLCLC1 = 2840.d0
    REAL(KIND=8),PARAMETER :: TLCLC2 = 3.5d0
    REAL(KIND=8),PARAMETER :: TLCLC3 = 4.805d0
    REAL(KIND=8),PARAMETER :: TLCLC4 = 55.d0
    !  k
    REAL(KIND=8),PARAMETER :: THTECON1 = 3376.d0
    REAL(KIND=8),PARAMETER :: THTECON2 = 2.54d0
    REAL(KIND=8),PARAMETER :: THTECON3 = .81d0

    ! To get rid of compiler warnings
    tmkpari = 0
    qvppari = 0
    klev = 0
    klcl = 0


    ! the comments were taken from a mark stoelinga email, 23 apr 2007,
    ! in response to a user getting the "outside of lookup table bounds"
    ! error message.

    ! tmkpari  - initial temperature of parcel, k
    !    values of 300 okay. (not sure how much from this you can stray.)

    ! prspari - initial pressure of parcel, hpa
    !    values of 980 okay. (not sure how much from this you can stray.)

    ! thtecon1, thtecon2, thtecon3
    !     these are all constants, the first in k and the other two have
    !     no units.  values of 3376, 2.54, and 0.81 were stated as being
    !     okay.

    ! tlcl - the temperature at the parcel's lifted condensation level, k
    !        should be a reasonable atmospheric temperature around 250-300 k
    !        (398 is "way too high")

    ! qvppari - the initial water vapor mixing ratio of the parcel,
    !           kg/kg (should range from 0.000 to 0.025)
    !

    !  calculated the pressure at full sigma levels (a set of pressure
    !  levels that bound the layers represented by the vertical grid points)

    CALL dpfcalc(prs,sfp,prsf,miy,mjx,mkzh,ter_follow)

    !  before looping, set lookup table for getting temperature on
    !  a pseudoadiabat.

    !call dlookup_table(psadithte,psadiprs,psaditmk,psafile)

    !   do j=1,mjx-1
    DO j = 1,mjx
    !   do i=1,miy-1
      DO i = 1,miy
          cape(i,j,1) = 0.d0
          cin(i,j,1) = 0.d0

          IF (i3dflag .EQ. 1) THEN
              kpar1 = 2
              kpar2 = mkzh
          ELSE
              ! find parcel with max theta-e in lowest 3 km agl.
              ethmax = -1.d0
              DO k = mkzh,1,-1
                  IF (ght(i,j,k)-ter(i,j) .LT. 3000.d0) THEN
                      q = MAX(qvp(i,j,k), 1.d-15)
                      t = tmk(i,j,k)
                      p = prs(i,j,k)
                      e = q*p/(EPS + q)
                      tlcl = TLCLC1/ (LOG(t**TLCLC2/e)-TLCLC3) + TLCLC4
                      eth = t * (1000.d0/p)**(GAMMA*(1.d0 + GAMMAMD*q))*&
                              EXP((THTECON1/tlcl - THTECON2)*q*(1.d0 + THTECON3*q))
                      IF (eth .GT. ethmax) THEN
                          klev = k
                          ethmax = eth
                      END IF
                  END IF
              END DO
              kpar1 = klev
              kpar2 = klev

              ! Establish average properties of that parcel
              ! (over depth of approximately davg meters)

              ! davg=.1
              davg = 500.d0
              pavg = davg*prs(i,j,kpar1)*&
                  GRAV/(RGAS*tvirtual(tmk(i,j,kpar1), qvp(i,j,kpar1)))
              p2 = MIN(prs(i,j,kpar1)+.5d0*pavg, prsf(i,j,mkzh))
              p1 = p2 - pavg
              totthe = 0.d0
              totqvp = 0.d0
              totprs = 0.d0
              DO k = mkzh,2,-1
                  IF (prsf(i,j,k) .LE. p1) GOTO 35
                  IF (prsf(i,j,k-1) .GE. p2) GOTO 34
                  p = prs(i,j,k)
                  pup = prsf(i,j,k)
                  pdn = prsf(i,j,k-1)
                  q = MAX(qvp(i,j,k),1.d-15)
                  th = tmk(i,j,k)*(1000.d0/p)**(GAMMA*(1.d0 + GAMMAMD*q))
                  pp1 = MAX(p1,pdn)
                  pp2 = MIN(p2,pup)
                  IF (pp2 .GT. pp1) THEN
                      deltap = pp2 - pp1
                      totqvp = totqvp + q*deltap
                      totthe = totthe + th*deltap
                      totprs = totprs + deltap
                  END IF
    34        CONTINUE
              END DO
    35        CONTINUE
              qvppari = totqvp/totprs
              tmkpari = (totthe/totprs)*&
                  (prs(i,j,kpar1)/1000.d0)**(GAMMA*(1.d0+GAMMAMD*qvp(i,j,kpar1)))
          END IF

          DO kpar = kpar1,kpar2

              ! Calculate temperature and moisture properties of parcel
              ! (note, qvppari and tmkpari already calculated above for 2d case.)

              IF (i3dflag .EQ. 1) THEN
                  qvppari = qvp(i,j,kpar)
                  tmkpari = tmk(i,j,kpar)
              END IF
              prspari = prs(i,j,kpar)
              ghtpari = ght(i,j,kpar)
              gammam = GAMMA * (1.d0 + GAMMAMD*qvppari)
              cpm = CP * (1.d0 + CPMD*qvppari)

              e = MAX(1.d-20,qvppari*prspari/(EPS + qvppari))
              tlcl = TLCLC1/(LOG(tmkpari**TLCLC2/e) - TLCLC3) + TLCLC4
              ethpari = tmkpari*(1000.d0/prspari)**(GAMMA*(1.d0 + GAMMAMD*qvppari))*&
                  EXP((THTECON1/tlcl - THTECON2)*qvppari*(1.d0 + THTECON3*qvppari))
              zlcl = ghtpari + (tmkpari - tlcl)/(GRAV/cpm)

              ! Calculate buoyancy and relative height of lifted parcel at
              ! all levels, and store in bottom up arrays.  add a level at the lcl,
              ! and at all points where buoyancy is zero.
              !
              ! For arrays that go bottom to top
              kk = 0
              ilcl = 0

              IF (ghtpari .GE. zlcl) THEN
                  ! Initial parcel already saturated or supersaturated.
                  ilcl = 2
                  klcl = 1
              END IF

              DO k = kpar,1,-1
                  ! For arrays that go bottom to top
    33            kk = kk + 1

                  ! Model level is below lcl
                  IF (ght(i,j,k).lt.zlcl) THEN
                      qvplift = qvppari
                      tmklift = tmkpari - GRAV/cpm*(ght(i,j,k) - ghtpari)
                      tvenv = tvirtual(tmk(i,j,k), qvp(i,j,k))
                      tvlift = tvirtual(tmklift, qvplift)
                      ghtlift = ght(i,j,k)
                  ELSE IF (ght(i,j,k) .GE. zlcl .AND. ilcl .EQ. 0) THEN
                      ! This model level and previous model level straddle the lcl,
                      ! so first create a new level in the bottom-up array, at the lcl.
                      tmklift = tlcl
                      qvplift = qvppari
                      facden = ght(i,j,k) - ght(i,j,k+1)
                      fac1 = (zlcl-ght(i,j,k+1))/facden
                      fac2 = (ght(i,j,k)-zlcl)/facden
                      tmkenv = tmk(i,j,k+1)*fac2 + tmk(i,j,k)*fac1
                      qvpenv = qvp(i,j,k+1)*fac2 + qvp(i,j,k)*fac1
                      tvenv = tvirtual(tmkenv, qvpenv)
                      tvlift = tvirtual(tmklift, qvplift)
                      ghtlift = zlcl
                      ilcl = 1
                  ELSE
                      tmklift = tonpsadiabat(ethpari, prs(i,j,k), PSADITHTE, PSADIPRS,&
                                             PSADITMK, GAMMA, throw_exception)
                      eslift = EZERO*EXP(ESLCON1*(tmklift - CELKEL)/(tmklift - ESLCON2))
                      qvplift = EPS*eslift/(prs(i,j,k) - eslift)
                      tvenv = tvirtual(tmk(i,j,k), qvp(i,j,k))
                      tvlift = tvirtual(tmklift, qvplift)
                      ghtlift = ght(i,j,k)
                  END IF
                  !  Buoyancy
                  buoy(kk) = GRAV*(tvlift - tvenv)/tvenv
                  zrel(kk) = ghtlift - ghtpari

                  IF ((kk .GT. 1) .AND. (buoy(kk)*buoy(kk-1) .LT. 0.0d0)) THEN
                      ! Parcel ascent curve crosses sounding curve, so create a new level
                      ! in the bottom-up array at the crossing.
                      kk = kk + 1
                      buoy(kk) = buoy(kk-1)
                      zrel(kk) = zrel(kk-1)
                      buoy(kk-1) = 0.d0
                      zrel(kk-1) = zrel(kk-2) + buoy(kk-2)/&
                          (buoy(kk-2) - buoy(kk))*(zrel(kk) - zrel(kk-2))
                  END IF

                  IF (ilcl .EQ. 1) THEN
                      klcl = kk
                      ilcl = 2
                      GOTO 33
                  END IF
              END DO

              kmax = kk
              IF (kmax .GT. 150) THEN
                  ! Need an exception here
                  CALL throw_exception('capecalc3d: kmax got too big. kmax=',kmax)
                  !STOP
              END IF

              ! If no lcl was found, set klcl to kmax.  it is probably not really
              ! at kmax, but this will make the rest of the routine behave
              ! properly.
              IF (ilcl .EQ. 0) klcl=kmax

              ! Get the accumulated buoyant energy from the parcel's starting
              ! point, at all levels up to the top level.
              benaccum(1) = 0.0d0
              benamin = 9d9
              DO k = 2,kmax
                  dz = zrel(k) - zrel(k-1)
                  benaccum(k) = benaccum(k-1) + .5d0*dz*(buoy(k-1) + buoy(k))
                  IF (benaccum(k) .LT. benamin) THEN
                      benamin = benaccum(k)
                  END IF
              END DO

              ! Determine equilibrium level (el), which we define as the highest
              ! level of non-negative buoyancy above the lcl. note, this may be
              ! the top level if the parcel is still buoyant there.

              DO k = kmax,klcl,-1
                  IF (buoy(k) .GE. 0.d0) THEN
                      ! k of equilibrium level
                      kel = k
                      GOTO 50
                  END IF
              END DO

              ! If we got through that loop, then there is no non-negative
              ! buoyancy above the lcl in the sounding.  in these situations,
              ! both cape and cin will be set to -0.1 j/kg. (see below about
              ! missing values in v6.1.0). also, where cape is
              ! non-zero, cape and cin will be set to a minimum of +0.1 j/kg, so
              ! that the zero contour in either the cin or cape fields will
              ! circumscribe regions of non-zero cape.

              ! In v6.1.0 of ncl, we added a _fillvalue attribute to the return
              ! value of this function. at that time we decided to change -0.1
              ! to a more appropriate missing value, which is passed into this
              ! routine as cmsg.

              ! cape(i,j,kpar) = -0.1d0
              ! cin(i,j,kpar) = -0.1d0
              cape(i,j,kpar) = cmsg
              cin(i,j,kpar)  = cmsg
              klfc = kmax

              GOTO 102

    50        CONTINUE

              !   If there is an equilibrium level, then cape is positive.  we'll
              !   define the level of free convection (lfc) as the point below the
              !   el, but at or above the lcl, where accumulated buoyant energy is a
              !   minimum.  the net positive area (accumulated buoyant energy) from
              !   the lfc up to the el will be defined as the cape, and the net
              !   negative area (negative of accumulated buoyant energy) from the
              !   parcel starting point to the lfc will be defined as the convective
              !   inhibition (cin).

              !   First get the lfc according to the above definition.
              benamin = 9d9
              klfc = kmax
              DO k = klcl,kel
                  IF (benaccum(k) .LT. benamin) THEN
                      benamin = benaccum(k)
                      klfc = k
                  END IF
              END DO

              ! Now we can assign values to cape and cin

              cape(i,j,kpar) = MAX(benaccum(kel)-benamin, 0.1d0)
              cin(i,j,kpar) = MAX(-benamin, 0.1d0)

              ! cin is uninteresting when cape is small (< 100 j/kg), so set
              ! cin to -0.1 (see note about missing values in v6.1.0) in
              ! that case.

              ! In v6.1.0 of ncl, we added a _fillvalue attribute to the return
              ! value of this function. at that time we decided to change -0.1
              ! to a more appropriate missing value, which is passed into this
              ! routine as cmsg.

              ! IF (cape(i,j,kpar).lt.100.d0) cin(i,j,kpar) = -0.1d0
              IF (cape(i,j,kpar) .LT. 100.d0) cin(i,j,kpar) = cmsg
    102   CONTINUE

          END DO

          IF (i3dflag .EQ. 0) THEN
              cape(i,j,mkzh) = cape(i,j,kpar1)
              cin(i,j,mkzh) = cin(i,j,kpar1)
    !  meters agl
              cin(i,j,mkzh-1) = zrel(klcl) + ghtpari - ter(i,j)
    !  meters agl
              cin(i,j,mkzh-2) = zrel(klfc) + ghtpari - ter(i,j)

          ENDIF
      END DO
    END DO

    RETURN
END SUBROUTINE f_computecape
