source: LMDZ6/trunk/libf/phylmdiso/rrtm/gpxyb.F90 @ 3927

Last change on this file since 3927 was 3927, checked in by Laurent Fairhead, 3 years ago

Initial import of the physics wih isotopes from Camille Risi
CR
File size: 8.5 KB
Line 
1SUBROUTINE GPXYB(KPROMA,KSTART,KPROF,KFLEV,PVDELB,PVC,&
2 & PRES,PDELP,PRDELP,PLNPR,PALPH,PRTGR,&
3 & PRPRES,PRPP) 
4
5!**** *GPXYB* - Computes auxillary arrays
6
7!     Purpose.
8!     --------
9!           COMPUTES AUXILLARY ARRAYS RELATED TO THE HYBRID COORDINATE
10
11!**   Interface.
12!     ----------
13!        *CALL* *GPXYB(..)
14
15!        Explicit arguments :
16!        --------------------
17!     KPROMA : dimensioning
18!     KSTART : start of work
19!     KPROF  : depth of work
20!     KFLEV     : vert. dimensioning
21
22!     PVDELB(KPROMA,0:KFLEV) : related to vert. coordinate        (input)
23!     PVC   (KPROMA,0:KFLEV) :  "     "      "     "    "         (input)
24!     PRES (KPROMA,0:KFLEV)  : HALF LEVEL PRESSURE                (input)
25!     PDELP (KPROMA,KFLEV)   : PRESSURE DIFFERENCE ACROSS LAYERS  (output)
26!     PRDELP(KPROMA,KFLEV)   : THEIR INVERSE                      (output)
27!     PLNPR (KPROMA,KFLEV)   : LOGARITHM OF RATIO OF PRESSURE     (output)
28!     PALPH (KPROMA,KFLEV)   : COEFFICIENTS OF THE HYDROSTATICS   (output)
29!     PRTGR (KPROMA,KFLEV)   : FOR PRES. GRAD. TERM AND ENE. CONV.(output)
30!                              ((rssavnabla prehyd/prehyd)_[layer]
31!                              = prtgr_[layer] * (rssavnabla prehyds))
32!     PRPRES(KPROMA,KFLEV)   : inverse of HALF LEVEL PRESSURE     (output)
33!     PRPP  (KPROMA,KFLEV)   : inverse of PRES(J)*PRES(J-1)       (output)
34
35!        Implicit arguments :  None.
36!        --------------------
37
38!     Method.
39!     -------
40!        See documentation
41
42!     Externals.      None.
43!     ----------
44
45!     Reference.
46!     ----------
47!        ECMWF Research Department documentation of the IFS
48
49!     Author.
50!     -------
51!        Mats Hamrud and Philippe Courtier  *ECMWF*
52
53!     Modifications.
54!     --------------
55!        Original : 88-02-04
56!        Modified : 94-10-11 by Radmila Bubnova: correction in the case
57!                            of the other approximation of d (ln p).
58!        Modified : 99-06-04 Optimisation   D.SALMOND
59!        Modified : 02-03-08 K. YESSAD: consistent discretisations of
60!                    "alpha" (PALPH) and "prtgr" (PRTGR)
61!                    for finite element vertical discretisation
62!                    to allow model to run with MF-physics and DDH.
63!        Modified : 03-07-07 J. Hague:  Replace divides with reciprocal
64!        M.Hamrud      01-Oct-2003 CY28 Cleaning
65!        Modified : 15-Feb-2005 by K. YESSAD: ZTOPPRES becomes TOPPRES
66!     ------------------------------------------------------------------
67
68USE PARKIND1  ,ONLY : JPIM     ,JPRB
69USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
70
71USE YOMDYN   , ONLY : NDLNPR   ,RHYDR0
72USE YOMCST   , ONLY : RD       ,RCVD
73USE YOMGEM   , ONLY : VDELA    ,VAF      ,VBF      ,TOPPRES
74USE YOMCVER  , ONLY : LVERTFE
75
76!     ------------------------------------------------------------------
77
78IMPLICIT NONE
79
80INTEGER(KIND=JPIM),INTENT(IN)    :: KPROMA
81INTEGER(KIND=JPIM),INTENT(IN)    :: KFLEV
82INTEGER(KIND=JPIM),INTENT(IN)    :: KSTART
83INTEGER(KIND=JPIM),INTENT(IN)    :: KPROF
84REAL(KIND=JPRB)   ,INTENT(IN)    :: PVDELB(KFLEV)
85REAL(KIND=JPRB)   ,INTENT(IN)    :: PVC(KFLEV)
86REAL(KIND=JPRB)   ,INTENT(IN)    :: PRES(KPROMA,0:KFLEV)
87REAL(KIND=JPRB)   ,INTENT(INOUT) :: PDELP(KPROMA,KFLEV)
88REAL(KIND=JPRB)   ,INTENT(INOUT) :: PRDELP(KPROMA,KFLEV)
89REAL(KIND=JPRB)   ,INTENT(INOUT) :: PLNPR(KPROMA,KFLEV)
90REAL(KIND=JPRB)   ,INTENT(OUT)   :: PALPH(KPROMA,KFLEV)
91REAL(KIND=JPRB)   ,INTENT(OUT)   :: PRTGR(KPROMA,KFLEV)
92REAL(KIND=JPRB)   ,INTENT(OUT)   :: PRPRES(KPROMA,KFLEV)
93REAL(KIND=JPRB)   ,INTENT(INOUT) :: PRPP(KPROMA,KFLEV)
94
95!     ------------------------------------------------------------------
96
97INTEGER(KIND=JPIM) :: IFIRST, JLEV, JLON, JJ, JTEMP, JM
98
99REAL(KIND=JPRB) :: ZPRESF
100REAL(KIND=JPRB) :: ZRPRES(KPROMA,2)
101REAL(KIND=JPRB) :: ZPRESFD
102REAL(KIND=JPRB) :: ZHOOK_HANDLE
103
104!     ------------------------------------------------------------------
105
106#include "abor1.intfb.h"
107
108!     ------------------------------------------------------------------
109
110IF (LHOOK) CALL DR_HOOK('GPXYB',0,ZHOOK_HANDLE)
111
112!     ------------------------------------------------------------------
113
114!*       0.    Level to begin normal computations
115!              ----------------------------------
116
117! This is introduced to allow the use of GPXYB without the implicit
118! assumption that the top level input for pressure is 0 hPa. This
119! is used in the surface observation operators where you do not want
120! to compute geopotential at all model levels.
121! The first block if is for economy (no do loop start up) and the second
122! for safety.
123!print *,'GPXYB: NDLNPR RHYDR0=',NDLNPR,RHYDR0
124TOPPRES=0.1  !!!!! A REVOIR (MPL) 29042010 passe de 0 a 0.1 comme ARPEGE
125IF(PRES(KSTART,0) <= TOPPRES)THEN
126  IFIRST=2
127ELSE
128  IFIRST=1
129  DO JLON=KSTART,KPROF
130    IF(PRES(JLON,0) <= TOPPRES)then
131      IFIRST=2
132      EXIT
133    ENDIF
134  ENDDO
135ENDIF
136!     ------------------------------------------------------------------
137
138!*       1.    COMPUTES EVERYTHING.
139!              --------------------
140
141!print *,'NDLNPR LVERTFE',NDLNPR,LVERTFE
142IF(NDLNPR == 0) THEN
143
144  IF(LVERTFE) THEN
145    DO JLEV=1,KFLEV
146      DO JLON=KSTART,KPROF
147        PDELP(JLON,JLEV)=VDELA(JLEV) + PVDELB(JLEV)*PRES(JLON,KFLEV)
148        PRDELP(JLON,JLEV)=1.0_JPRB/PDELP(JLON,JLEV)
149        ZPRESF =VAF(JLEV) + VBF(JLEV)*PRES(JLON,KFLEV)
150        ZPRESFD=1.0_JPRB/ZPRESF
151        PLNPR(JLON,JLEV)=PDELP(JLON,JLEV)*ZPRESFD
152        ! * PRTGR needed for DDH and option LVERCOR=T.
153        !   for finite element vertical discretisation,
154        !   "prtgr_[layer]" is simply B_[layer]/prehyd_[layer]
155        PRTGR (JLON,JLEV)=VBF(JLEV)*ZPRESFD
156        ! * PALPH needed for MF physics:
157        PALPH(JLON,JLEV)=(PRES(JLON,JLEV)-ZPRESF)*ZPRESFD
158      ENDDO
159    ENDDO
160  ELSE
161    JJ=1
162    JM=2
163    DO JLON=KSTART,KPROF
164      ZRPRES(JLON,JM)=1.0_JPRB/PRES(JLON,IFIRST-1)
165    ENDDO
166    DO JLEV=IFIRST,KFLEV
167      DO JLON=KSTART,KPROF
168        ZRPRES(JLON,JJ)=1.0_JPRB/PRES(JLON,JLEV)
169        PDELP (JLON,JLEV)=PRES(JLON,JLEV)-PRES(JLON,JLEV-1)
170        PRDELP(JLON,JLEV)=1.0_JPRB/PDELP(JLON,JLEV)
171        PLNPR (JLON,JLEV)=LOG(PRES(JLON,JLEV)*ZRPRES(JLON,JM))
172        PRPRES(JLON,JLEV)=ZRPRES(JLON,JJ)
173        PALPH (JLON,JLEV)=1.0_JPRB-PRES(JLON,JLEV-1)*PRDELP(JLON,JLEV)&
174         & *PLNPR(JLON,JLEV) 
175        PRPP  (JLON,JLEV)=ZRPRES(JLON,JJ)*ZRPRES(JLON,JM)
176        PRTGR (JLON,JLEV)=PRDELP(JLON,JLEV)&
177         & *(PVDELB(JLEV)+PVC(JLEV)*PLNPR(JLON,JLEV)*PRDELP(JLON,&
178         & JLEV)) 
179!       print *,'GPXYB JLEV JLON JJ PRES ZPRES PDELP ', JLEV,JLON,JJ,PRES(JLON,JLEV),ZRPRES(JLON,JJ),PDELP(JLON,JLEV)
180!       print *,'GPXYB JLEV JLON JM PRDELP PLNPR ', JLEV,JLON,JM,PRDELP(JLON,JLEV),PLNPR (JLON,JLEV)
181!       print *,'GPXYB JLEV JLON JJ PRPRES PALPH ', JLEV,JLON,JJ,PRPRES(JLON,JLEV),PALPH (JLON,JLEV)
182!       print *,'GPXYB JLEV JLON PRPP PRTGR PVDELB PVC ', JLEV,JLON,PRPP  (JLON,JLEV),PRTGR (JLON,JLEV),PVDELB(JLEV),PVC(JLEV)
183      ENDDO
184      JTEMP=JM
185      JM=JJ
186      JJ=JTEMP
187    ENDDO
188    DO JLEV=1,IFIRST-1
189      DO JLON=KSTART,KPROF
190        PDELP (JLON,JLEV)=PRES(JLON,JLEV)-PRES(JLON,JLEV-1)
191        PRDELP(JLON,JLEV)=1.0_JPRB/PDELP(JLON,JLEV)
192        PLNPR (JLON,JLEV)=LOG(PRES(JLON,1)/TOPPRES)
193        PRPRES(JLON,JLEV)=1.0_JPRB/PRES(JLON,1)
194        PALPH (JLON,JLEV)=RHYDR0
195        PRPP  (JLON,JLEV)=1.0_JPRB/(PRES(JLON,1)*TOPPRES)
196        PRTGR (JLON,JLEV)=PRDELP(JLON,JLEV)*PVDELB(JLEV)
197      ENDDO
198    ENDDO
199  ENDIF
200
201ELSEIF(NDLNPR == 1) THEN
202  IF(LVERTFE) THEN
203    CALL ABOR1(' LVERTFE=.T. NOT COMPATIBLE WITH NDLNPR == 1')
204  ENDIF
205
206  DO JLEV=IFIRST,KFLEV
207    DO JLON=KSTART,KPROF
208      PDELP (JLON,JLEV)=PRES(JLON,JLEV)-PRES(JLON,JLEV-1)
209      PRDELP(JLON,JLEV)=1.0_JPRB/PDELP(JLON,JLEV)
210      PRPP  (JLON,JLEV)=1.0_JPRB/(PRES(JLON,JLEV)*PRES(JLON,JLEV-1))
211      PLNPR (JLON,JLEV)=PDELP(JLON,JLEV)*SQRT(PRPP(JLON,JLEV))
212      PALPH (JLON,JLEV)=1.0_JPRB-PRES(JLON,JLEV-1)*PRDELP(JLON,JLEV)&
213       & *PLNPR(JLON,JLEV) 
214      PRTGR (JLON,JLEV)=PRDELP(JLON,JLEV)&
215       & *(PVDELB(JLEV)+PVC(JLEV)*PLNPR(JLON,JLEV)*PRDELP(JLON,&
216       & JLEV)) 
217      PRPRES(JLON,JLEV)=1.0_JPRB/PRES(JLON,JLEV)
218    ENDDO
219  ENDDO
220
221  DO JLEV=1,IFIRST-1
222    DO JLON=KSTART,KPROF
223      PDELP (JLON,JLEV)=PRES(JLON,JLEV)
224      PRDELP(JLON,JLEV)=1.0_JPRB/PDELP(JLON,JLEV)
225      PLNPR (JLON,JLEV)=2.0_JPRB+RCVD/RD
226      PALPH (JLON,JLEV)=1.0_JPRB
227      PRTGR (JLON,JLEV)=PRDELP(JLON,JLEV)*PVDELB(JLEV)
228      PRPRES(JLON,JLEV)=1.0_JPRB/PRES(JLON,1)
229      PRPP  (JLON,JLEV)=(PLNPR(JLON,JLEV)*PRDELP(JLON,JLEV))**2
230    ENDDO
231  ENDDO
232
233ENDIF
234
235! (PLNPR(JLON,1) AND PRPP(JLON,1) ARE A PRIORI NOT USED LATER)
236
237!     ------------------------------------------------------------------
238
239IF (LHOOK) CALL DR_HOOK('GPXYB',1,ZHOOK_HANDLE)
240END SUBROUTINE GPXYB
Note: See TracBrowser for help on using the repository browser.