gpxyb.F90

Go to the documentation of this file.

SUBROUTINE gpxyb(KPROMA,KSTART,KPROF,KFLEV,PVDELB,PVC,&
 & pres,pdelp,prdelp,plnpr,palph,prtgr,&
 & prpres,prpp)  

!**** *GPXYB* - Computes auxillary arrays

!     Purpose.
!     --------
!           COMPUTES AUXILLARY ARRAYS RELATED TO THE HYBRID COORDINATE

!**   Interface.
!     ----------
!        *CALL* *GPXYB(..)

!        Explicit arguments :
!        --------------------
!     KPROMA : dimensioning
!     KSTART : start of work
!     KPROF  : depth of work
!     KFLEV     : vert. dimensioning

!     PVDELB(KPROMA,0:KFLEV) : related to vert. coordinate        (input)
!     PVC   (KPROMA,0:KFLEV) :  "     "      "     "    "         (input)
!     PRES (KPROMA,0:KFLEV)  : HALF LEVEL PRESSURE                (input)
!     PDELP (KPROMA,KFLEV)   : PRESSURE DIFFERENCE ACROSS LAYERS  (output)
!     PRDELP(KPROMA,KFLEV)   : THEIR INVERSE                      (output)
!     PLNPR (KPROMA,KFLEV)   : LOGARITHM OF RATIO OF PRESSURE     (output)
!     PALPH (KPROMA,KFLEV)   : COEFFICIENTS OF THE HYDROSTATICS   (output)
!     PRTGR (KPROMA,KFLEV)   : FOR PRES. GRAD. TERM AND ENE. CONV.(output)
!                              ((rssavnabla prehyd/prehyd)_[layer]
!                              = prtgr_[layer] * (rssavnabla prehyds))
!     PRPRES(KPROMA,KFLEV)   : inverse of HALF LEVEL PRESSURE     (output)
!     PRPP  (KPROMA,KFLEV)   : inverse of PRES(J)*PRES(J-1)       (output)

!        Implicit arguments :  None.
!        --------------------

!     Method.
!     -------
!        See documentation

!     Externals.      None.
!     ----------

!     Reference.
!     ----------
!        ECMWF Research Department documentation of the IFS

!     Author.
!     -------
!        Mats Hamrud and Philippe Courtier  *ECMWF*

!     Modifications.
!     --------------
!        Original : 88-02-04
!        Modified : 94-10-11 by Radmila Bubnova: correction in the case
!                            of the other approximation of d (ln p).
!        Modified : 99-06-04 Optimisation   D.SALMOND
!        Modified : 02-03-08 K. YESSAD: consistent discretisations of
!                    "alpha" (PALPH) and "prtgr" (PRTGR)
!                    for finite element vertical discretisation
!                    to allow model to run with MF-physics and DDH.
!        Modified : 03-07-07 J. Hague:  Replace divides with reciprocal
!        M.Hamrud      01-Oct-2003 CY28 Cleaning
!        Modified : 15-Feb-2005 by K. YESSAD: ZTOPPRES becomes TOPPRES
!     ------------------------------------------------------------------

USE parkind1  ,ONLY : jpim     ,jprb
USE yomhook   ,ONLY : lhook,   dr_hook

USE yomdyn   , ONLY : ndlnpr   ,rhydr0
USE yomcst   , ONLY : rd       ,rcvd
USE yomgem   , ONLY : vdela    ,vaf      ,vbf      ,toppres
USE yomcver  , ONLY : lvertfe

!     ------------------------------------------------------------------

IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN)    :: KPROMA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KFLEV 
INTEGER(KIND=JPIM),INTENT(IN)    :: KSTART 
INTEGER(KIND=JPIM),INTENT(IN)    :: KPROF 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PVDELB(kflev) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PVC(kflev) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PRES(kproma,0:kflev) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PDELP(kproma,kflev) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PRDELP(kproma,kflev) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PLNPR(kproma,kflev) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PALPH(kproma,kflev) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PRTGR(kproma,kflev) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PRPRES(kproma,kflev) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PRPP(kproma,kflev) 

!     ------------------------------------------------------------------

INTEGER(KIND=JPIM) :: IFIRST, JLEV, JLON, JJ, JTEMP, JM

REAL(KIND=JPRB) :: ZPRESF
REAL(KIND=JPRB) :: ZRPRES(kproma,2)
REAL(KIND=JPRB) :: ZPRESFD
REAL(KIND=JPRB) :: ZHOOK_HANDLE

!     ------------------------------------------------------------------

#include "abor1.intfb.h"

!     ------------------------------------------------------------------

IF (lhook) CALL dr_hook('GPXYB',0,zhook_handle)

!     ------------------------------------------------------------------

!*       0.    Level to begin normal computations
!              ----------------------------------

! This is introduced to allow the use of GPXYB without the implicit
! assumption that the top level input for pressure is 0 hPa. This
! is used in the surface observation operators where you do not want
! to compute geopotential at all model levels.
! The first block if is for economy (no do loop start up) and the second
! for safety.
!print *,'GPXYB: NDLNPR RHYDR0=',NDLNPR,RHYDR0
toppres=0.1  !!!!! A REVOIR (MPL) 29042010 passe de 0 a 0.1 comme ARPEGE
IF(pres(kstart,0) <= toppres)THEN
  ifirst=2
ELSE
  ifirst=1
  DO jlon=kstart,kprof
    IF(pres(jlon,0) <= toppres)then
      ifirst=2
      EXIT
    ENDIF
  ENDDO
ENDIF
!     ------------------------------------------------------------------

!*       1.    COMPUTES EVERYTHING.
!              --------------------

!print *,'NDLNPR LVERTFE',NDLNPR,LVERTFE
IF(ndlnpr == 0) THEN

  IF(lvertfe) THEN
    DO jlev=1,kflev
      DO jlon=kstart,kprof
        pdelp(jlon,jlev)=vdela(jlev) + pvdelb(jlev)*pres(jlon,kflev)
        prdelp(jlon,jlev)=1.0_jprb/pdelp(jlon,jlev)
        zpresf =vaf(jlev) + vbf(jlev)*pres(jlon,kflev)
        zpresfd=1.0_jprb/zpresf
        plnpr(jlon,jlev)=pdelp(jlon,jlev)*zpresfd
        ! * PRTGR needed for DDH and option LVERCOR=T.
        !   for finite element vertical discretisation,
        !   "prtgr_[layer]" is simply B_[layer]/prehyd_[layer]
        prtgr(jlon,jlev)=vbf(jlev)*zpresfd
        ! * PALPH needed for MF physics:
        palph(jlon,jlev)=(pres(jlon,jlev)-zpresf)*zpresfd
      ENDDO
    ENDDO
  ELSE
    jj=1
    jm=2
    DO jlon=kstart,kprof
      zrpres(jlon,jm)=1.0_jprb/pres(jlon,ifirst-1)
    ENDDO
    DO jlev=ifirst,kflev
      DO jlon=kstart,kprof
        zrpres(jlon,jj)=1.0_jprb/pres(jlon,jlev)
        pdelp(jlon,jlev)=pres(jlon,jlev)-pres(jlon,jlev-1)
        prdelp(jlon,jlev)=1.0_jprb/pdelp(jlon,jlev)
        plnpr(jlon,jlev)=log(pres(jlon,jlev)*zrpres(jlon,jm))
        prpres(jlon,jlev)=zrpres(jlon,jj)
        palph(jlon,jlev)=1.0_jprb-pres(jlon,jlev-1)*prdelp(jlon,jlev)&
         & *plnpr(jlon,jlev)  
        prpp(jlon,jlev)=zrpres(jlon,jj)*zrpres(jlon,jm)
        prtgr(jlon,jlev)=prdelp(jlon,jlev)&
         & *(pvdelb(jlev)+pvc(jlev)*plnpr(jlon,jlev)*prdelp(jlon,&
         & jlev))  
!       print *,'GPXYB JLEV JLON JJ PRES ZPRES PDELP ', JLEV,JLON,JJ,PRES(JLON,JLEV),ZRPRES(JLON,JJ),PDELP(JLON,JLEV)
!       print *,'GPXYB JLEV JLON JM PRDELP PLNPR ', JLEV,JLON,JM,PRDELP(JLON,JLEV),PLNPR (JLON,JLEV)
!       print *,'GPXYB JLEV JLON JJ PRPRES PALPH ', JLEV,JLON,JJ,PRPRES(JLON,JLEV),PALPH (JLON,JLEV)
!       print *,'GPXYB JLEV JLON PRPP PRTGR PVDELB PVC ', JLEV,JLON,PRPP  (JLON,JLEV),PRTGR (JLON,JLEV),PVDELB(JLEV),PVC(JLEV)
      ENDDO
      jtemp=jm
      jm=jj
      jj=jtemp
    ENDDO
    DO jlev=1,ifirst-1
      DO jlon=kstart,kprof
        pdelp(jlon,jlev)=pres(jlon,jlev)-pres(jlon,jlev-1)
        prdelp(jlon,jlev)=1.0_jprb/pdelp(jlon,jlev)
        plnpr(jlon,jlev)=log(pres(jlon,1)/toppres)
        prpres(jlon,jlev)=1.0_jprb/pres(jlon,1)
        palph(jlon,jlev)=rhydr0
        prpp(jlon,jlev)=1.0_jprb/(pres(jlon,1)*toppres)
        prtgr(jlon,jlev)=prdelp(jlon,jlev)*pvdelb(jlev)
      ENDDO
    ENDDO
  ENDIF

ELSEIF(ndlnpr == 1) THEN
  IF(lvertfe) THEN
    CALL abor1(' LVERTFE=.T. NOT COMPATIBLE WITH NDLNPR == 1')
  ENDIF

  DO jlev=ifirst,kflev
    DO jlon=kstart,kprof
      pdelp(jlon,jlev)=pres(jlon,jlev)-pres(jlon,jlev-1)
      prdelp(jlon,jlev)=1.0_jprb/pdelp(jlon,jlev)
      prpp(jlon,jlev)=1.0_jprb/(pres(jlon,jlev)*pres(jlon,jlev-1))
      plnpr(jlon,jlev)=pdelp(jlon,jlev)*sqrt(prpp(jlon,jlev))
      palph(jlon,jlev)=1.0_jprb-pres(jlon,jlev-1)*prdelp(jlon,jlev)&
       & *plnpr(jlon,jlev)  
      prtgr(jlon,jlev)=prdelp(jlon,jlev)&
       & *(pvdelb(jlev)+pvc(jlev)*plnpr(jlon,jlev)*prdelp(jlon,&
       & jlev))  
      prpres(jlon,jlev)=1.0_jprb/pres(jlon,jlev)
    ENDDO
  ENDDO

  DO jlev=1,ifirst-1
    DO jlon=kstart,kprof
      pdelp(jlon,jlev)=pres(jlon,jlev)
      prdelp(jlon,jlev)=1.0_jprb/pdelp(jlon,jlev)
      plnpr(jlon,jlev)=2.0_jprb+rcvd/rd
      palph(jlon,jlev)=1.0_jprb
      prtgr(jlon,jlev)=prdelp(jlon,jlev)*pvdelb(jlev)
      prpres(jlon,jlev)=1.0_jprb/pres(jlon,1)
      prpp(jlon,jlev)=(plnpr(jlon,jlev)*prdelp(jlon,jlev))**2
    ENDDO
  ENDDO

ENDIF

! (PLNPR(JLON,1) AND PRPP(JLON,1) ARE A PRIORI NOT USED LATER)

!     ------------------------------------------------------------------

IF (lhook) CALL dr_hook('GPXYB',1,zhook_handle)
END SUBROUTINE gpxyb