copyright (C) 2001  MSC-RPN COMM  %%%MC2%%%
***s/r hydprt1 -- Computes hydrostatic pressure on thermodymic levels
*

      subroutine hydprt1 ( prt,psm,ngf,nkf,tp,ht,hm,sbxy,odx,ody,area, 1,10
     $                                    lminx,lmaxx,lminy,lmaxy,lnk )
      implicit none
*
      integer ngf,nkf,lminx,lmaxx,lminy,lmaxy,lnk
      real prt(ngf,nkf), psm(ngf), area(ngf),
     $  tp(lminx:lmaxx,lminy:lmaxy,*), ht(lminx:lmaxx,lminy:lmaxy,*),
     $  hm(lminx:lmaxx,lminy:lmaxy,0:lnk), sbxy(lminx:lmaxx,lminy:lmaxy)
      real*8 odx(*),ody(*)
*
*AUTHOR   Michel Desgagne           Mar   1996
*
*ARGUMENTS
*    NAMES     I/O  TYPE  A/S        DESCRIPTION
*
*    prt       I/O    R    A    hydrostatic presure on thermo. lvls
*    psm       I/O    R    A    surface pressure (t-dt)
*    area       O     R    A    cell area (m**2)
*    ngf        I     I    S    physics horizontal dimension
*    nkf        I     I    S    physics vertical   dimension
*    tp         I     R    A    temperature   (t)
*    nk         I     I    S    vertical dimension
*
*IMPLICIT
#include "lcldim.cdk"
#include "yomdyn1.cdk"
#include "consdyn_8.cdk"
#include "levels.cdk"
c#include "dynmem.cdk"
*
*MODULES
*
**
      integer i,j,k
      real p00,gamma,w1(minx:maxx,miny:maxy)
*
**----------------------------------------------------------------------
*
      p00  = 100000.
      gamma= 6.5e-3
*
!$omp parallel
!$omp do
      do j=1,ldnj
         do i=1,ldni
            prt ((j-1)*ldni+i,nkf-1)= prt ((j-1)*ldni+i,nkf  )*
     $                 exp(-grav_8*(ht(i,j,1)-hm(i,j,0))/
     $             (rgasd_8*(tp(i,j,1)+gamma*(hm(i,j,1)-hm(i,j,0))/2.)))
            w1(i,j)= prt ((j-1)*ldni+i,nkf-1)*
     $                 exp(-grav_8*(hm(i,j,1)-ht(i,j,1))/
     $             (rgasd_8*(5.*tp(i,j,1)/6.+tp(i,j,2)/6.)))
            prt ((j-1)*ldni+i,nkf-2)= w1(i,j)*
     $                 exp(-grav_8*(ht(i,j,2)-hm(i,j,1))/
     $             (rgasd_8*(2.*tp(i,j,2)/3.+tp(i,j,1)/3.)))
            w1(i,j)= w1(i,j)*exp(-grav_8*(hm(i,j,2)-hm(i,j,1))
     $                         /(rgasd_8*tp(i,j,2)))
            area((j-1)*ldni+i) = 1./(odx(1)*ody(j)) / sbxy(i,j)
         end do
*
         do k=3,lnk
         do i=1,ldni
            prt((j-1)*ldni+i,nkf-k)=w1(i,j)*
     $             exp(-grav_8*(ht(i,j,k)-hm(i,j,k-1))/
     $             (rgasd_8*(0.75*tp(i,j,k)+0.25*tp(i,j,k-1))))
            w1(i,j)=w1(i,j)*exp(-grav_8*(hm(i,j,k)-hm(i,j,k-1))/
     $             (rgasd_8*tp(i,j,k)))
         end do
         end do
*
      end do
!$omp end do
!$omp end parallel
*
      do k=1,nkf
      do i=ldni*ldnj+1,ngf
         prt(i,k) = prt(ldni*ldnj,k)
      end do
      end do
*
      do i=ldni*ldnj+1,ngf
         psm (i) = psm (ldni*ldnj)
         area(i) = area(ldni*ldnj)
      end do
*
*----------------------------------------------------------------------
      return
      end