pkg/fizhi/fizhi_utils.F

a456aa407c Andr*0001 #include "FIZHI_OPTIONS.h"
d8b3764985 Andr*0002       function minval (q,im)
                0003       implicit none
                0004       integer  im, i
a456aa407c Andr*0005       _RL q(im), minval
d8b3764985 Andr*0006       minval = 1.e15
                0007       do i=1,im
                0008       if( q(i).lt.minval ) minval = q(i)
                0009       enddo
                0010       return
                0011       end
                0012       FUNCTION ERRF (ARG)
                0013 C***********************************************************************
                0014 C  FUNCTION ERRF
                0015 C  PURPOSE
                0016 C     COMPUTES ERROR FUNCTION OF ARGUMENT
                0017 C  USAGE
                0018 C     CALLED BY TRBFLX
                0019 C  DESCRIPTION OF PARAMETERS
                0020 C     ARG   -  INPUTED ARGUMENT
                0021 C  REMARKS:
                0022 C        USED TO COMPUTE FRACTIONAL CLOUD COVER AND LIQUID WATER CONTENT
                0023 C          FROM TURBULENCE STATISTICS
                0024 C **********************************************************************
28e90e66e5 Andr*0025       implicit none
a456aa407c Andr*0026       _RL arg,errf
0889f02121 Jean*0027
a456aa407c Andr*0028       _RL aa1,aa2,aa3,aa4,aa5,pp,x2,x3,x4,x5
d8b3764985 Andr*0029       PARAMETER ( AA1  =   0.254829592 )
                0030       PARAMETER ( AA2  =  -0.284496736 )
                0031       PARAMETER ( AA3  =   1.421413741 )
                0032       PARAMETER ( AA4  =  -1.453152027 )
                0033       PARAMETER ( AA5  =   1.061405429 )
                0034       PARAMETER ( PP   =   0.3275911   )
                0035       PARAMETER ( X2   =   AA2 / AA1   )
                0036       PARAMETER ( X3   =   AA3 / AA2   )
                0037       PARAMETER ( X4   =   AA5 / AA3   )
                0038       PARAMETER ( X5   =   AA5 / AA4   )
0889f02121 Jean*0039
a456aa407c Andr*0040       _RL aarg,tt
0889f02121 Jean*0041
d8b3764985 Andr*0042       ERRF = 1.
                0043       AARG=ABS(ARG)
0889f02121 Jean*0044
d8b3764985 Andr*0045       IF ( AARG .LT. 4.0 )  THEN
                0046        TT = 1./(1.+PP*AARG)
                0047        ERRF = 1. -
                0048      1  (AA1*TT*(1.+X2*TT*(1.+X3*TT*(1.+X4*TT*(1.+X5*TT)))))
                0049      2   * EXP(-AARG*AARG)
                0050       ENDIF
0889f02121 Jean*0051
d8b3764985 Andr*0052       IF ( ARG .LT. 0.0 )  ERRF = -ERRF
0889f02121 Jean*0053
d8b3764985 Andr*0054       RETURN
                0055       END
                0056
0889f02121 Jean*0057       SUBROUTINE STRIP(A,B,IA,IB,L,K)
28e90e66e5 Andr*0058       implicit none
                0059       integer ia,ib,L,K
0889f02121 Jean*0060       _RL A(IA,L), B(IB,L)
                0061
471f7f004b Andr*0062       INTEGER OFFSET,Lena,i,j
0889f02121 Jean*0063
                0064       OFFSET = IB*(K-1)
                0065       Lena    = MIN(IB,IA-OFFSET)
                0066       OFFSET = OFFSET+1
                0067
                0068       IF(Lena.EQ.IB) THEN
                0069         DO 100 J=1,L
                0070         DO 100 I=1,Lena
                0071         B(I,J) = A(I+OFFSET-1,J)
                0072 100     CONTINUE
                0073       ELSE
                0074         DO 200 J=1,L
                0075         DO 300 I=1,Lena
                0076         B(I,J) = A(I+OFFSET-1,J)
                0077 300     CONTINUE
                0078         DO 400 I=1,IB-Lena
                0079         B(Lena+I,J) = A(Lena+OFFSET-1,J)
                0080 400     CONTINUE
                0081 200     CONTINUE
                0082       ENDIF
                0083
                0084       RETURN
                0085       END
                0086       SUBROUTINE STRIPINT(A,B,IA,IB,L,K)
d416560405 Andr*0087       implicit none
                0088       integer ia,ib,L,K
0889f02121 Jean*0089       INTEGER A(IA,L), B(IB,L)
                0090
d416560405 Andr*0091       INTEGER OFFSET,Lena,i,j
0889f02121 Jean*0092
                0093       OFFSET = IB*(K-1)
                0094       Lena    = MIN(IB,IA-OFFSET)
                0095       OFFSET = OFFSET+1
                0096
                0097       IF(Lena.EQ.IB) THEN
                0098         DO 100 J=1,L
                0099         DO 100 I=1,Lena
                0100         B(I,J) = A(I+OFFSET-1,J)
                0101 100     CONTINUE
                0102       ELSE
                0103         DO 200 J=1,L
                0104         DO 300 I=1,Lena
                0105         B(I,J) = A(I+OFFSET-1,J)
                0106 300     CONTINUE
                0107         DO 400 I=1,IB-Lena
                0108         B(Lena+I,J) = A(Lena+OFFSET-1,J)
                0109 400     CONTINUE
                0110 200     CONTINUE
                0111       ENDIF
                0112
                0113       RETURN
                0114       END
                0115       SUBROUTINE PASTE(B,A,IB,IA,L,K)
28e90e66e5 Andr*0116       implicit none
                0117       integer ia,ib,L,K
0889f02121 Jean*0118       _RL A(IA,L), B(IB,L)
28e90e66e5 Andr*0119
471f7f004b Andr*0120       INTEGER OFFSET,Lena,i,j
0889f02121 Jean*0121
                0122       OFFSET = IB*(K-1)
                0123       Lena    = MIN(IB,IA-OFFSET)
                0124       OFFSET = OFFSET+1
                0125
                0126       DO 100 J=1,L
                0127       DO 100 I=1,Lena
                0128       A(I+OFFSET-1,J) = B(I,J)
                0129 100   CONTINUE
                0130
                0131       RETURN
                0132       END
                0133       SUBROUTINE PSTBMP(B,A,IB,IA,L,K)
28e90e66e5 Andr*0134       implicit none
                0135       integer ia,ib,L,K
0889f02121 Jean*0136       _RL A(IA,L), B(IB,L)
28e90e66e5 Andr*0137
471f7f004b Andr*0138       INTEGER OFFSET,Lena,i,j
0889f02121 Jean*0139
                0140       OFFSET = IB*(K-1)
                0141       Lena    = MIN(IB,IA-OFFSET)
                0142       OFFSET = OFFSET+1
                0143
                0144       DO 100 J=1,L
                0145       DO 100 I=1,Lena
                0146       A(I+OFFSET-1,J) = A(I+OFFSET-1,J) + B(I,J)
                0147 100   CONTINUE
                0148 C
                0149       RETURN
                0150       END
d8b3764985 Andr*0151       SUBROUTINE STRINT(A,B,IA,IB,L,K)
28e90e66e5 Andr*0152       implicit none
                0153       integer ia,ib,L,K
                0154       INTEGER A(IA,L), B(IB,L)
                0155
471f7f004b Andr*0156       INTEGER OFFSET,Lena,i,j
d8b3764985 Andr*0157
                0158       OFFSET = IB*(K-1)
471f7f004b Andr*0159       Lena = MIN(IB,IA-OFFSET)
d8b3764985 Andr*0160       OFFSET = OFFSET+1
                0161
471f7f004b Andr*0162       IF(Lena.EQ.IB) THEN
d8b3764985 Andr*0163         DO 100 J=1,L
471f7f004b Andr*0164         DO 100 I=1,Lena
d8b3764985 Andr*0165         B(I,J) = A(I+OFFSET-1,J)
                0166 100     CONTINUE
                0167       ELSE
                0168         DO 200 J=1,L
471f7f004b Andr*0169         DO 300 I=1,Lena
d8b3764985 Andr*0170         B(I,J) = A(I+OFFSET-1,J)
                0171 300     CONTINUE
471f7f004b Andr*0172         DO 400 I=1,IB-Lena
                0173         B(Lena+I,J) = A(Lena+OFFSET-1,J)
d8b3764985 Andr*0174 400     CONTINUE
                0175 200     CONTINUE
                0176       ENDIF
                0177
                0178       RETURN
                0179       END
                0180       SUBROUTINE QSAT (TT,P,Q,DQDT,LDQDT)
                0181 C***********************************************************************
                0182 C
                0183 C  PURPOSE:
                0184 C  ========
                0185 C    Compute Saturation Specific Humidity
                0186 C
                0187 C  INPUT:
                0188 C  ======
                0189 C    TT ......... Temperature (Kelvin)
                0190 C    P .......... Pressure (mb)
                0191 C    LDQDT ...... Logical Flag to compute QSAT Derivative
                0192 C
                0193 C  OUTPUT:
                0194 C  =======
                0195 C    Q .......... Saturation Specific Humidity
                0196 C    DQDT ....... Saturation Specific Humidity Derivative wrt Temperature
                0197 C
                0198 C
                0199 C***********************************************************************
                0200
                0201       IMPLICIT NONE
a456aa407c Andr*0202       _RL TT, P, Q, DQDT
d8b3764985 Andr*0203       LOGICAL LDQDT
28e90e66e5 Andr*0204
a456aa407c Andr*0205       _RL AIRMW, H2OMW
0889f02121 Jean*0206
                0207       PARAMETER ( AIRMW  = 28.97      )
                0208       PARAMETER ( H2OMW  = 18.01      )
d8b3764985 Andr*0209
a456aa407c Andr*0210       _RL ESFAC, ERFAC
d8b3764985 Andr*0211       PARAMETER ( ESFAC = H2OMW/AIRMW       )
                0212       PARAMETER ( ERFAC = (1.0-ESFAC)/ESFAC )
                0213
a456aa407c Andr*0214       _RL aw0, aw1, aw2, aw3, aw4, aw5, aw6
                0215       _RL bw0, bw1, bw2, bw3, bw4, bw5, bw6
                0216       _RL ai0, ai1, ai2, ai3, ai4, ai5, ai6
                0217       _RL bi0, bi1, bi2, bi3, bi4, bi5, bi6
d8b3764985 Andr*0218
a456aa407c Andr*0219       _RL d0, d1, d2, d3, d4, d5, d6
                0220       _RL e0, e1, e2, e3, e4, e5, e6
                0221       _RL f0, f1, f2, f3, f4, f5, f6
                0222       _RL g0, g1, g2, g3, g4, g5, g6
d8b3764985 Andr*0223
                0224 c ********************************************************
                0225 c ***  Polynomial Coefficients WRT Water (Lowe, 1977) ****
                0226 c ***              (Valid +50 C to -50 C)             ****
                0227 c ********************************************************
                0228
                0229       parameter ( aw0 =  6.107799961e+00 * esfac )
                0230       parameter ( aw1 =  4.436518521e-01 * esfac )
                0231       parameter ( aw2 =  1.428945805e-02 * esfac )
                0232       parameter ( aw3 =  2.650648471e-04 * esfac )
                0233       parameter ( aw4 =  3.031240396e-06 * esfac )
                0234       parameter ( aw5 =  2.034080948e-08 * esfac )
                0235       parameter ( aw6 =  6.136820929e-11 * esfac )
                0236
                0237       parameter ( bw0 = +4.438099984e-01 * esfac )
                0238       parameter ( bw1 = +2.857002636e-02 * esfac )
                0239       parameter ( bw2 = +7.938054040e-04 * esfac )
                0240       parameter ( bw3 = +1.215215065e-05 * esfac )
                0241       parameter ( bw4 = +1.036561403e-07 * esfac )
                0242       parameter ( bw5 = +3.532421810e-10 * esfac )
                0243       parameter ( bw6 = -7.090244804e-13 * esfac )
                0244
                0245 c ********************************************************
                0246 c ***   Polynomial Coefficients WRT Ice  (Lowe, 1977) ****
                0247 c ***              (Valid  +0 C to -50 C)             ****
                0248 c ********************************************************
                0249
                0250       parameter ( ai0 = +6.109177956e+00 * esfac )
                0251       parameter ( ai1 = +5.034698970e-01 * esfac )
                0252       parameter ( ai2 = +1.886013408e-02 * esfac )
                0253       parameter ( ai3 = +4.176223716e-04 * esfac )
                0254       parameter ( ai4 = +5.824720280e-06 * esfac )
                0255       parameter ( ai5 = +4.838803174e-08 * esfac )
                0256       parameter ( ai6 = +1.838826904e-10 * esfac )
                0257
                0258       parameter ( bi0 = +5.030305237e-01 * esfac )
                0259       parameter ( bi1 = +3.773255020e-02 * esfac )
                0260       parameter ( bi2 = +1.267995369e-03 * esfac )
                0261       parameter ( bi3 = +2.477563108e-05 * esfac )
                0262       parameter ( bi4 = +3.005693132e-07 * esfac )
                0263       parameter ( bi5 = +2.158542548e-09 * esfac )
                0264       parameter ( bi6 = +7.131097725e-12 * esfac )
                0265
                0266 c ********************************************************
                0267 c ***         Polynomial Coefficients WRT Ice         ****
                0268 c ***   Starr and Cox (1985) (Valid -40 C to -70 C)   ****
                0269 c ********************************************************
                0270
                0271       parameter ( d0 = 0.535098336e+01 * esfac )
                0272       parameter ( d1 = 0.401390832e+00 * esfac )
                0273       parameter ( d2 = 0.129690326e-01 * esfac )
                0274       parameter ( d3 = 0.230325039e-03 * esfac )
                0275       parameter ( d4 = 0.236279781e-05 * esfac )
                0276       parameter ( d5 = 0.132243858e-07 * esfac )
                0277       parameter ( d6 = 0.314296723e-10 * esfac )
                0278
                0279       parameter ( e0 = 0.469290530e+00 * esfac )
                0280       parameter ( e1 = 0.333092511e-01 * esfac )
                0281       parameter ( e2 = 0.102164528e-02 * esfac )
                0282       parameter ( e3 = 0.172979242e-04 * esfac )
                0283       parameter ( e4 = 0.170017544e-06 * esfac )
                0284       parameter ( e5 = 0.916466531e-09 * esfac )
                0285       parameter ( e6 = 0.210844486e-11 * esfac )
                0286
                0287 c ********************************************************
                0288 c ***         Polynomial Coefficients WRT Ice         ****
                0289 c ***   Starr and Cox (1985) (Valid -65 C to -95 C)   ****
                0290 c ********************************************************
                0291
                0292       parameter ( f0 = 0.298152339e+01 * esfac )
                0293       parameter ( f1 = 0.191372282e+00 * esfac )
                0294       parameter ( f2 = 0.517609116e-02 * esfac )
                0295       parameter ( f3 = 0.754129933e-04 * esfac )
                0296       parameter ( f4 = 0.623439266e-06 * esfac )
                0297       parameter ( f5 = 0.276961083e-08 * esfac )
                0298       parameter ( f6 = 0.516000335e-11 * esfac )
                0299
                0300       parameter ( g0 = 0.312654072e+00 * esfac )
                0301       parameter ( g1 = 0.195789002e-01 * esfac )
                0302       parameter ( g2 = 0.517837908e-03 * esfac )
                0303       parameter ( g3 = 0.739410547e-05 * esfac )
                0304       parameter ( g4 = 0.600331350e-07 * esfac )
                0305       parameter ( g5 = 0.262430726e-09 * esfac )
                0306       parameter ( g6 = 0.481960676e-12 * esfac )
                0307
a456aa407c Andr*0308       _RL TMAX, TICE
d8b3764985 Andr*0309       PARAMETER ( TMAX=323.15, TICE=273.16)
0889f02121 Jean*0310
a456aa407c Andr*0311       _RL T, D, W, QX, DQX
d8b3764985 Andr*0312       T = MIN(TT,TMAX) - TICE
                0313       DQX = 0.
                0314       QX  = 0.
                0315
                0316 c Fitting for temperatures above 0 degrees centigrade
                0317 c ---------------------------------------------------
                0318       if(t.gt.0.) then
                0319        qx = aw0+T*(aw1+T*(aw2+T*(aw3+T*(aw4+T*(aw5+T*aw6)))))
                0320       if (ldqdt)  then
                0321       dqx = bw0+T*(bw1+T*(bw2+T*(bw3+T*(bw4+T*(bw5+T*bw6)))))
                0322       endif
                0323       endif
                0324
                0325 c Fitting for temperatures between 0 and -40
                0326 c ------------------------------------------
                0327       if( t.le.0. .and. t.gt.-40.0 ) then
                0328         w = (40.0 + t)/40.0
                0329        qx =     w *(aw0+T*(aw1+T*(aw2+T*(aw3+T*(aw4+T*(aw5+T*aw6))))))
                0330      .    + (1.-w)*(ai0+T*(ai1+T*(ai2+T*(ai3+T*(ai4+T*(ai5+T*ai6))))))
                0331       if (ldqdt)  then
                0332       dqx =     w *(bw0+T*(bw1+T*(bw2+T*(bw3+T*(bw4+T*(bw5+T*bw6))))))
                0333      .    + (1.-w)*(bi0+T*(bi1+T*(bi2+T*(bi3+T*(bi4+T*(bi5+T*bi6))))))
                0334       endif
                0335       endif
                0336
                0337 c Fitting for temperatures between -40 and -70
                0338 c --------------------------------------------
                0339       if( t.le.-40.0 .and. t.ge.-70.0 ) then
                0340        qx = d0+T*(d1+T*(d2+T*(d3+T*(d4+T*(d5+T*d6)))))
                0341       if (ldqdt) then
                0342       dqx = e0+T*(e1+T*(e2+T*(e3+T*(e4+T*(e5+T*e6)))))
                0343       endif
                0344       endif
                0345
                0346 c Fitting for temperatures less than -70
                0347 c --------------------------------------
                0348       if(t.lt.-70.0) then
                0349        qx = f0+t*(f1+t*(f2+t*(f3+t*(f4+t*(f5+t*f6)))))
                0350       if (ldqdt) then
                0351       dqx = g0+t*(g1+t*(g2+t*(g3+t*(g4+t*(g5+t*g6)))))
                0352       endif
                0353       endif
                0354
                0355 c Compute Saturation Specific Humidity
                0356 c ------------------------------------
                0357       D = (P-ERFAC*QX)
                0358       IF(D.LT.0.) THEN
                0359        Q = 1.0
                0360        IF (LDQDT)  DQDT = 0.
                0361       ELSE
                0362        D = 1.0 / D
32362b8fa8 Cons*0363        Q = MIN(QX * D,1.0 _d 0)
d8b3764985 Andr*0364        IF (LDQDT)  DQDT = (1.0 + ERFAC*Q) * D * DQX
                0365       ENDIF
                0366       RETURN
                0367       END
                0368       subroutine vqsat (tt,p,q,dqdt,ldqdt,n)
                0369       implicit none
                0370       integer  i,n
                0371       logical  ldqdt
a456aa407c Andr*0372       _RL tt(n), p(n), q(n), dqdt(n)
06d4643e1f Jean*0373 #ifdef FIZHI_CRAY
4e1c053948 Jean*0374 #ifdef FIZHI_F77_COMPIL
d8b3764985 Andr*0375 cfpp$ expand (qsat)
                0376 #endif
                0377 #endif
                0378       do i=1,n
                0379       call qsat ( tt(i),p(i),q(i),dqdt(i),ldqdt )
                0380       enddo
                0381       return
                0382       end
                0383
                0384       subroutine stripit(a,b,irun,ia,ib,l,k)
                0385       implicit none
                0386       integer ia,ib,irun,l,k
a456aa407c Andr*0387       _RL a(ia,l), b(ib,l)
471f7f004b Andr*0388       integer   i,j,Lena,offset
d8b3764985 Andr*0389
                0390       offset = ib*(k-1)
471f7f004b Andr*0391       Lena    = min(ib,irun-offset)
d8b3764985 Andr*0392       offset = offset+1
                0393
471f7f004b Andr*0394       if(Lena.eq.ib) then
d8b3764985 Andr*0395         do 100 j=1,l
471f7f004b Andr*0396         do 100 i=1,Lena
d8b3764985 Andr*0397         b(i,j) = a(i+offset-1,j)
                0398 100     continue
                0399       else
                0400         do 200 j=1,l
471f7f004b Andr*0401         do 300 i=1,Lena
d8b3764985 Andr*0402         b(i,j) = a(i+offset-1,j)
                0403 300     continue
471f7f004b Andr*0404         do 400 i=1,ib-Lena
                0405         b(Lena+i,j) = a(Lena+offset-1,j)
d8b3764985 Andr*0406 400     continue
                0407 200     continue
                0408       endif
                0409       return
                0410       end
                0411
                0412       subroutine stripitint(a,b,irun,ia,ib,l,k)
                0413       implicit none
                0414       integer ia,ib,irun,l,k,a(ia,l),b(ib,l)
471f7f004b Andr*0415       integer i,j,Lena,offset
d8b3764985 Andr*0416
                0417       offset = ib*(k-1)
471f7f004b Andr*0418       Lena = min(ib,irun-offset)
d8b3764985 Andr*0419       offset = offset+1
                0420
471f7f004b Andr*0421       if(Lena.eq.ib) then
d8b3764985 Andr*0422         do 100 j=1,l
471f7f004b Andr*0423         do 100 i=1,Lena
d8b3764985 Andr*0424         b(i,j) = a(i+offset-1,j)
                0425 100     continue
                0426       else
                0427         do 200 j=1,l
471f7f004b Andr*0428         do 300 i=1,Lena
d8b3764985 Andr*0429         b(i,j) = a(i+offset-1,j)
                0430 300     continue
471f7f004b Andr*0431         do 400 i=1,ib-Lena
                0432         b(Lena+i,j) = a(Lena+offset-1,j)
d8b3764985 Andr*0433 400     continue
                0434 200     continue
                0435       endif
                0436       return
                0437       end
                0438
                0439       subroutine pastit(b,a,ib,ia,irun,L,k)
                0440       implicit none
471f7f004b Andr*0441       integer ib,ia,L,k,irun,Lena,offset
d8b3764985 Andr*0442       integer i,j
a456aa407c Andr*0443       _RL a(ia,l), b(ib,l)
d8b3764985 Andr*0444
                0445       offset = ib*(k-1)
471f7f004b Andr*0446       Lena    = min(ib,irun-offset)
d8b3764985 Andr*0447       offset = offset+1
                0448
                0449       do 100 j=1,L
471f7f004b Andr*0450       do 100 i=1,Lena
d8b3764985 Andr*0451       a(i+offset-1,j) = b(i,j)
                0452 100   continue
                0453       return
                0454       end
                0455
                0456       subroutine pstbitint(b,a,ib,ia,irun,l,k)
                0457       implicit none
471f7f004b Andr*0458       integer ib,ia,L,k,irun,Lena,offset
a456aa407c Andr*0459       _RL a(ia,l)
28e90e66e5 Andr*0460       integer b(ib,l)
d8b3764985 Andr*0461       integer i,j
                0462
                0463       offset = ib*(k-1)
471f7f004b Andr*0464       Lena    = min(ib,irun-offset)
d8b3764985 Andr*0465       offset = offset+1
0889f02121 Jean*0466
d8b3764985 Andr*0467       do 100 j=1,L
471f7f004b Andr*0468       do 100 i=1,Lena
d8b3764985 Andr*0469       a(i+offset-1,j) = a(i+offset-1,j) + float(b(i,j))
                0470 100   continue
0889f02121 Jean*0471       return
                0472       end
d8b3764985 Andr*0473
                0474       subroutine pstbmpit(b,a,ib,ia,irun,l,k)
                0475       implicit none
471f7f004b Andr*0476       integer ib,ia,L,k,irun,Lena,offset
a456aa407c Andr*0477       _RL a(ia,l), b(ib,l)
d8b3764985 Andr*0478       integer i,j
                0479
                0480       offset = ib*(k-1)
471f7f004b Andr*0481       Lena    = min(ib,irun-offset)
d8b3764985 Andr*0482       offset = offset+1
0889f02121 Jean*0483
d8b3764985 Andr*0484       do 100 j=1,L
471f7f004b Andr*0485       do 100 i=1,Lena
d8b3764985 Andr*0486       a(i+offset-1,j) = a(i+offset-1,j) + b(i,j)
                0487 100   continue
0889f02121 Jean*0488       return
                0489       end
d8b3764985 Andr*0490
471f7f004b Andr*0491       subroutine strip2tile(a,indx,b,irun,ia,ib,levs,npeice)
d8b3764985 Andr*0492 c-----------------------------------------------------------------------
                0493 c subroutine strip2tile - extract one processors worth of grid points
                0494 c                         from a grid space array to a stripped tile
                0495 c                         space array
0889f02121 Jean*0496 c
d8b3764985 Andr*0497 c  input:
                0498 c      a      - array to be stripped FROM [ia,levs]
471f7f004b Andr*0499 c      indx   - array of horizontal indeces of grid points to convert to
d8b3764985 Andr*0500 c               tile space
                0501 c      irun   - number of points in array a that need to be stripped
                0502 c      ia     - inner of dimension of source array
                0503 c      ib     - inner dimension of target array AND the number of points
                0504 c               in the target array to be filled
                0505 c      levs   - number of vertical levels AND outer dimension of arrays
                0506 c      npeice - the current strip number to be filled
                0507 c output:
                0508 c      b      - array to be filled, ie, one processors field [ib,levs]
                0509 c-----------------------------------------------------------------------
                0510       implicit none
                0511       integer ia,ib,irun,levs,npeice
a456aa407c Andr*0512       _RL a(ia,levs), b(ib,levs)
471f7f004b Andr*0513       integer indx(irun)
                0514       integer i,k,Lena,offset
d8b3764985 Andr*0515
                0516       offset = ib*(npeice-1)
471f7f004b Andr*0517       Lena    = min(ib,irun-offset)
d8b3764985 Andr*0518       offset = offset+1
                0519
471f7f004b Andr*0520       if(Lena.eq.ib) then
d8b3764985 Andr*0521         do 100 k=1,levs
471f7f004b Andr*0522         do 100 i=1,Lena
                0523         b(i,k) = a(indx(i+offset-1),k)
d8b3764985 Andr*0524 100     continue
                0525       else
                0526         do 200 k=1,levs
471f7f004b Andr*0527         do 300 i=1,Lena
                0528         b(i,k) = a(indx(i+offset-1),k)
d8b3764985 Andr*0529 300     continue
471f7f004b Andr*0530         do 400 i=1,ib-Lena
                0531         b(Lena+i,k) = a(indx(Lena+offset-1),k)
d8b3764985 Andr*0532 400     continue
                0533 200     continue
                0534       endif
                0535       return
                0536       end
                0537
471f7f004b Andr*0538       subroutine paste2grd_old(b,indx,chfr,ib,numpts,a,ia,levs,npeice)
d8b3764985 Andr*0539 c-----------------------------------------------------------------------
                0540 c subroutine paste2grd - paste one processors worth of grid points
                0541 c                        from a stripped tile array to a grid
                0542 c                        space array
0889f02121 Jean*0543 c
d8b3764985 Andr*0544 c  input:
                0545 c      b      - array to be pasted back into grid space [ib,levs]
471f7f004b Andr*0546 c      indx   - array of horizontal indeces of grid points to convert to
d8b3764985 Andr*0547 c               tile space[numpts]
                0548 c      chfr   - fractional area covered by the tile [ib]
0889f02121 Jean*0549 c      ib     - inner dimension of source array AND number of points in
d8b3764985 Andr*0550 c               array a that need to be pasted
                0551 c      numpts - total number of points which were stripped
                0552 c      ia     - inner of dimension of target array
                0553 c      levs   - number of vertical levels AND outer dimension of arrays
                0554 c      npeice - the current strip number to be filled
                0555 c output:
                0556 c      a      - grid space array to be filled [ia,levs]
                0557 c
                0558 c IMPORTANT NOTE:
                0559 c
                0560 c This routine will result in roundoff differences if called from
                0561 c within a parallel region.
                0562 c-----------------------------------------------------------------------
                0563
                0564       implicit none
                0565       integer ia,ib,levs,numpts,npeice
471f7f004b Andr*0566       integer indx(numpts)
a456aa407c Andr*0567       _RL a(ia,levs), b(ib,levs), chfr(ib)
d8b3764985 Andr*0568
471f7f004b Andr*0569       integer i,L,offset,Lena
d8b3764985 Andr*0570
                0571       offset = ib*(npeice-1)
471f7f004b Andr*0572       Lena    = min(ib,numpts-offset)
d8b3764985 Andr*0573       offset = offset+1
                0574
                0575       do L = 1,levs
471f7f004b Andr*0576       do i=1,Lena
                0577        a(indx(i+offset-1),L) = a(indx(i+offset-1),L) + b(i,L)*chfr(i)
d8b3764985 Andr*0578       enddo
                0579       enddo
                0580       return
                0581       end
471f7f004b Andr*0582       subroutine paste2grd (b,indx,chfr,ib,numpts,a,ia,levs,npeice,
28e90e66e5 Andr*0583      .                                                            check)
d8b3764985 Andr*0584 c-----------------------------------------------------------------------
                0585 c subroutine paste2grd - paste one processors worth of grid points
                0586 c                        from a stripped tile array to a grid
                0587 c                        space array
0889f02121 Jean*0588 c
d8b3764985 Andr*0589 c  input:
                0590 c      b      - array to be pasted back into grid space [ib,levs]
471f7f004b Andr*0591 c      indx   - array of horizontal indeces of grid points to convert to
d8b3764985 Andr*0592 c               tile space[numpts]
                0593 c      chfr   - fractional area covered by the tile [ib]
0889f02121 Jean*0594 c      ib     - inner dimension of source array AND number of points in
d8b3764985 Andr*0595 c               array a that need to be pasted
                0596 c      numpts - total number of points which were stripped
                0597 c      ia     - inner of dimension of target array
                0598 c      levs   - number of vertical levels AND outer dimension of arrays
                0599 c      npeice - the current strip number to be filled
                0600 c      check  - logical to check for undefined values
                0601 c output:
                0602 c      a      - grid space array to be filled [ia,levs]
                0603 c
                0604 c IMPORTANT NOTE:
                0605 c
                0606 c This routine will result in roundoff differences if called from
                0607 c within a parallel region.
                0608 c-----------------------------------------------------------------------
                0609
                0610       implicit none
                0611       integer ia,ib,levs,numpts,npeice
471f7f004b Andr*0612       integer indx(numpts)
a456aa407c Andr*0613       _RL a(ia,levs), b(ib,levs), chfr(ib)
d8b3764985 Andr*0614       logical check
                0615
471f7f004b Andr*0616       integer i,L,offset,Lena
a456aa407c Andr*0617       _RL    undef,getcon
d8b3764985 Andr*0618
                0619       offset = ib*(npeice-1)
471f7f004b Andr*0620       Lena    = min(ib,numpts-offset)
d8b3764985 Andr*0621       offset = offset+1
                0622
                0623       if( check ) then
                0624       undef = getcon('UNDEF')
28e90e66e5 Andr*0625        do L= 1,levs
471f7f004b Andr*0626        do i= 1,Lena
                0627         if( a(indx(i+offset-1),L).eq.undef .or. b(i,L).eq.undef ) then
                0628          a(indx(i+offset-1),L) = undef
28e90e66e5 Andr*0629         else
471f7f004b Andr*0630          a(indx(i+offset-1),L)=a(indx(i+offset-1),L) + b(i,L)*chfr(i)
28e90e66e5 Andr*0631         endif
                0632        enddo
                0633        enddo
d8b3764985 Andr*0634       else
28e90e66e5 Andr*0635        do L= 1,levs
471f7f004b Andr*0636        do i= 1,Lena
                0637         a(indx(i+offset-1),L)=a(indx(i+offset-1),L) + b(i,L)*chfr(i)
28e90e66e5 Andr*0638        enddo
                0639        enddo
d8b3764985 Andr*0640       endif
                0641
                0642       return
                0643       end
                0644       SUBROUTINE GRD2MSC(A,IM,JM,IGRD,B,MXCHPS,NCHP)
                0645
                0646       implicit none
                0647       integer im,jm,mxchps,nchp
                0648       integer igrd(mxchps)
0889f02121 Jean*0649 c     _RL A(IM,JM), B(MXCHPS)
                0650       _RL A(IM*JM), B(MXCHPS)
d8b3764985 Andr*0651
                0652       integer i
                0653
                0654       IF(NCHP.GE.0) THEN
                0655         DO I=1,NCHP
0889f02121 Jean*0656 c        B(I) = A(IGRD(I),1)
                0657          B(I) = A(IGRD(I))
d8b3764985 Andr*0658         ENDDO
                0659       ELSE
                0660         PRINT *, 'ERROR IN GRD2MSC'
                0661       ENDIF
                0662
                0663       RETURN
                0664       END
                0665
                0666       SUBROUTINE MSC2GRD(IGRD,CHFR,B,MXCHPS,NCHP,FRACG,A,IM,JM)
                0667
                0668       implicit none
a456aa407c Andr*0669       _RL zero,one
d8b3764985 Andr*0670       parameter ( one = 1.)
                0671       parameter (zero = 0.)
                0672       integer im,jm,mxchps,nchp
                0673       integer igrd(mxchps)
0889f02121 Jean*0674 c     _RL A(IM,JM), B(MXCHPS), CHFR(MXCHPS), FRACG(IM,JM)
                0675       _RL A(IM*JM), B(MXCHPS), CHFR(MXCHPS), FRACG(IM*JM)
d8b3764985 Andr*0676
0889f02121 Jean*0677 c     _RL VT1(IM,JM)
                0678       _RL VT1(IM*JM)
d8b3764985 Andr*0679       integer i
                0680
                0681       IF(NCHP.GE.0) THEN
                0682         DO I=1,IM*JM
0889f02121 Jean*0683 c        VT1(I,1) = ZERO
                0684          VT1(I) = ZERO
d8b3764985 Andr*0685         ENDDO
                0686
                0687         DO I=1,NCHP
0889f02121 Jean*0688 c        VT1(IGRD(I),1) = VT1(IGRD(I),1) + B(I)*CHFR(I)
                0689          VT1(IGRD(I)) = VT1(IGRD(I)) + B(I)*CHFR(I)
d8b3764985 Andr*0690         ENDDO
                0691
                0692         DO I=1,IM*JM
0889f02121 Jean*0693 c        A(I,1) = A(I,1)*(ONE-FRACG(I,1)) + VT1(I,1)
                0694          A(I) = A(I)*(ONE-FRACG(I)) + VT1(I)
d8b3764985 Andr*0695         ENDDO
                0696       ELSE
                0697         PRINT *, 'ERROR IN MSC2GRD'
                0698       ENDIF
                0699
                0700       RETURN
                0701       END
                0702
                0703       subroutine chpprm(nymd,nhms,mxchps,nchp,chlt,ityp,alai,
                0704      1       agrn,zoch,z2ch,cdrc,cdsc,sqsc,ufac,rsl1,rsl2,rdcs)
                0705
                0706        implicit none
                0707
                0708        integer nymd,nhms,nchp,mxchps,ityp(mxchps)
a456aa407c Andr*0709        _RL chlt(mxchps)
                0710        _RL alai(mxchps),agrn(mxchps)
                0711        _RL zoch(mxchps), z2ch(mxchps), cdrc(mxchps), cdsc(mxchps)
                0712        _RL sqsc(mxchps), ufac(mxchps), rsl1(mxchps), rsl2(mxchps)
                0713        _RL rdcs(mxchps)
d8b3764985 Andr*0714
                0715 C*********************************************************************
                0716 C********************* SUBROUTINE CHPPRM  ****************************
                0717 C**********************  14 JUNE 1991   ******************************
                0718 C*********************************************************************
                0719
                0720       integer ntyps
                0721       parameter (ntyps=10)
                0722
a456aa407c Andr*0723       _RL pblzet
d8b3764985 Andr*0724       parameter (pblzet = 50.)
                0725       integer k1,k2,nymd1,nhms1,nymd2,nhms2,i
a456aa407c Andr*0726       _RL getcon,vkrm,rootl,vroot,dum1,dum2,alphaf
                0727       _RL facm,facp
                0728       _RL scat,d
d8b3764985 Andr*0729
a456aa407c Andr*0730       _RL
d8b3764985 Andr*0731      &     vgdd(12, ntyps),    vgz0(12, ntyps),
                0732      &     vgrd(12, ntyps),    vgrt(12, ntyps),
                0733
                0734      &     vgrf11(ntyps),         vgrf12(ntyps),
                0735      &     vgtr11(ntyps),         vgtr12(ntyps),
                0736      &     vgroca(ntyps),         vgrotd(ntyps),
                0737      &     vgrdrs(ntyps),         vgz2  (ntyps)
                0738
                0739       data vgz0  /
0889f02121 Jean*0740      1   2.6530, 2.6530, 2.6530, 2.6530, 2.6530, 2.6530, 2.6530,
                0741      1   2.6530, 2.6530, 2.6530, 2.6530, 2.6530,
d8b3764985 Andr*0742      2   0.5200, 0.5200, 0.6660, 0.9100, 1.0310, 1.0440, 1.0420,
0889f02121 Jean*0743      2   1.0370, 1.0360, 0.9170, 0.6660, 0.5200,
d8b3764985 Andr*0744      3   1.1120, 1.1030, 1.0880, 1.0820, 1.0760, 1.0680, 1.0730,
0889f02121 Jean*0745      3   1.0790, 1.0820, 1.0880, 1.1030, 1.1120,
d8b3764985 Andr*0746      4   0.0777, 0.0778, 0.0778, 0.0779, 0.0778, 0.0771, 0.0759,
0889f02121 Jean*0747      4   0.0766, 0.0778, 0.0779, 0.0778, 0.0778,
d8b3764985 Andr*0748      5   0.2450, 0.2450, 0.2270, 0.2000, 0.2000, 0.2000, 0.2000,
0889f02121 Jean*0749      5   0.267,  0.292,  0.280,  0.258,  0.2450,
d8b3764985 Andr*0750      6   0.0752, 0.0752, 0.0752, 0.0752, 0.0752, 0.0757, 0.0777,
0889f02121 Jean*0751      6   0.0778, 0.0774, 0.0752, 0.0752, 0.0752,
                0752      7   0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
                0753      7   0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
                0754      8   0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
                0755      8   0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
                0756      9   0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
d8b3764985 Andr*0757      9   0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
0889f02121 Jean*0758      1   0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
d8b3764985 Andr*0759      1   0.0112, 0.0112, 0.0112, 0.0112, 0.0112
                0760      &  /
                0761
                0762       data vgrd    /
                0763      1   285.87, 285.87, 285.87, 285.87, 285.87, 285.87, 285.87,
0889f02121 Jean*0764      1   285.87, 285.87, 285.87, 285.87, 285.87,
d8b3764985 Andr*0765      2   211.32, 211.32, 218.78, 243.40, 294.87, 345.90, 355.18,
0889f02121 Jean*0766      2   341.84, 307.22, 244.84, 218.78, 211.32,
                0767      3   565.41, 587.05, 623.46, 638.13, 652.86, 675.04, 660.24,
                0768      3   645.49, 638.13, 623.46, 587.05, 565.41,
d8b3764985 Andr*0769      4    24.43,  24.63,  24.80,  24.96,  25.72,  27.74,  30.06,
0889f02121 Jean*0770      4    28.86,  25.90,  25.11,  24.80,  24.63,
d8b3764985 Andr*0771      5   103.60, 103.60, 102.35, 100.72, 100.72, 100.72, 100.72,
0889f02121 Jean*0772      5   105.30, 107.94, 106.59, 104.49, 103.60,
d8b3764985 Andr*0773      6    22.86,  22.86,  22.86,  22.86,  22.86,  23.01,  24.36,
0889f02121 Jean*0774      6    24.69,  24.04,  22.86,  22.86,  22.86,
d8b3764985 Andr*0775      7    23.76,  23.76,  23.76,  23.76,  23.76,  23.76,  23.76,
0889f02121 Jean*0776      7    23.76,  23.76,  23.76,  23.76,  23.76,
d8b3764985 Andr*0777      8    23.76,  23.76,  23.76,  23.76,  23.76,  23.76,  23.76,
0889f02121 Jean*0778      8    23.76,  23.76,  23.76,  23.76,  23.76,
d8b3764985 Andr*0779      9    23.76,  23.76,  23.76,  23.76,  23.76,  23.76,  23.76,
0889f02121 Jean*0780      9    23.76,  23.76,  23.76,  23.76,  23.76,
d8b3764985 Andr*0781      1    23.76,  23.76,  23.76,  23.76,  23.76,  23.76,  23.76,
                0782      1    23.76,  23.76,  23.76,  23.76,  23.76
                0783      &  /
                0784
                0785       data vgrt  /
                0786      1   19737.8, 19737.8, 19737.8, 19737.8, 19737.8, 19737.8, 19737.8,
                0787      1   19737.8, 19737.8, 19737.8, 19737.8, 19737.8,
                0788      2    5010.0,  5010.0,  5270.0,  6200.0,  8000.0,  9700.0,  9500.0,
                0789      2    8400.0,  6250.0,  5270.0,  5010.0,  5010.0,
0889f02121 Jean*0790      3    9000.0,  9200.0,  9533.3,  9666.7,  9800.0,  9866.7,  9733.3,
d8b3764985 Andr*0791      3    9666.7,  9533.3,  9200.0,  9000.0,  9000.0,
                0792      4    5500.0,  5625.0,  5750.0,  5875.0,  6625.0,  8750.0,  9375.0,
                0793      4    6875.0,  6000.0,  5750.0,  5625.0,  5500.0,
                0794      5    6500.0,  6000.0,  5500.0,  5500.0,  5500.0,  5500.0,  5500.0,
                0795      5    7500.0,  8500.0,  7000.0,  6500.0,  6500.0,
                0796      6   10625.0, 10625.0, 10625.0, 10625.0, 10625.0, 11250.0, 18750.0,
                0797      6   17500.0, 10625.0, 10625.0, 10625.0, 10625.0,
                0798      7       1.0,     1.0,     1.0,     1.0,     1.0,     1.0,     1.0,
                0799      7       1.0,     1.0,     1.0,     1.0,     1.0,
                0800      8       1.0,     1.0,     1.0,     1.0,     1.0,     1.0,     1.0,
                0801      8       1.0,     1.0,     1.0,     1.0,     1.0,
                0802      9       1.0,     1.0,     1.0,     1.0,     1.0,     1.0,     1.0,
                0803      9       1.0,     1.0,     1.0,     1.0,     1.0,
                0804      1       1.0,     1.0,     1.0,     1.0,     1.0,     1.0,     1.0,
                0805      1       1.0,     1.0,     1.0,     1.0,     1.0
                0806      &  /
                0807
                0808       data vgdd  /
                0809      1   27.37,   27.37,   27.37,   27.37,   27.37,   27.37,  27.37,
0889f02121 Jean*0810      1   27.37,   27.37,   27.37,   27.37,   27.37,
d8b3764985 Andr*0811      2   13.66,   13.66,   14.62,   15.70,   16.33,   16.62,   16.66,
0889f02121 Jean*0812      2   16.60,   16.41,   15.73,   14.62,   13.66,
d8b3764985 Andr*0813      3   13.76,   13.80,   13.86,   13.88,   13.90,   13.93,   13.91,
0889f02121 Jean*0814      3   13.89,   13.88,   13.86,   13.80,   13.76,
d8b3764985 Andr*0815      4    0.218,   0.227,   0.233,   0.239,   0.260,   0.299,   0.325,
0889f02121 Jean*0816      4    0.313,   0.265,   0.244,   0.233,   0.227,
d8b3764985 Andr*0817      5    2.813,   2.813,   2.662,   2.391,   2.391,   2.391,   2.391,
0889f02121 Jean*0818      5    2.975,   3.138,   3.062,   2.907,   2.813,
d8b3764985 Andr*0819      6    0.10629, 0.10629, 0.10629, 0.10629, 0.10629, 0.12299, 0.21521,
0889f02121 Jean*0820      6    0.22897, 0.19961, 0.10629, 0.10629, 0.10629,
d8b3764985 Andr*0821      7    0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,
0889f02121 Jean*0822      7    0.0001,  0.0001,  0.0001,  0.0001,  0.0001,
d8b3764985 Andr*0823      8    0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,
0889f02121 Jean*0824      8    0.0001,  0.0001,  0.0001,  0.0001,  0.0001,
d8b3764985 Andr*0825      9    0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,
                0826      9    0.0001,  0.0001,  0.0001,  0.0001,  0.0001,
                0827      1    0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,
0889f02121 Jean*0828      1    0.0001,  0.0001,  0.0001,  0.0001,  0.0001
d8b3764985 Andr*0829      &  /
                0830
                0831       data vgrf11 /0.10,0.10,0.07,0.105,0.10,0.10,.001,.001,.001,.001/
                0832
                0833       data vgrf12 /0.16,0.16,0.16,0.360,0.16,0.16,.001,.001,.001,.001/
                0834
                0835       data vgtr11 /0.05,0.05,0.05,0.070,0.05,0.05,.001,.001,.001,.001/
                0836
                0837       data vgtr12 /.001,.001,.001, .220,.001,.001,.001,.001,.001,.001/
                0838
                0839       data vgroca /
                0840      &  0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6,
0889f02121 Jean*0841      &   .1E-6, .1E-6, .1E-6, .1E-6  /
d8b3764985 Andr*0842
0889f02121 Jean*0843       data vgrotd /1.00,1.00,0.50,0.50,0.50,0.20,0.10,0.10,0.10,0.10/
d8b3764985 Andr*0844
                0845       data vgrdrs  /
                0846      &   0.75E13, 0.75E13, 0.75E13, 0.40E13, 0.75E13, 0.75E13,
                0847      &   0.10E13, 0.10E13, 0.10E13, 0.10E13  /
                0848
                0849       data vgz2 /35.0, 20.0, 17.0, 0.6, 5.0, 0.6, 0.1, 0.1, 0.1, 0.1/
                0850
                0851       vkrm = GETCON('VON KARMAN')
                0852
                0853       call time_bound  ( nymd,nhms, nymd1,nhms1, nymd2,nhms2, k1,k2 )
28e90e66e5 Andr*0854       call interp_time ( nymd,nhms, nymd1,nhms1, nymd2,nhms2, facm,facp)
d8b3764985 Andr*0855
                0856       do i=1,nchp
                0857
                0858       zoch(i) = vgz0(k2,ityp(i))*facp + vgz0(k1,ityp(i))*facm
                0859       rdcs(i) = vgrd(k2,ityp(i))*facp + vgrd(k1,ityp(i))*facm
                0860
                0861       rootl   = vgrt(k2,ityp(i))*facp + vgrt(k1,ityp(i))*facm
                0862
                0863       vroot  = rootl * vgroca(ityp (i))
                0864       dum1   = log (vroot / (1. - vroot))
                0865       dum2   = 1. / (8. * 3.14159 * rootl)
                0866       alphaf = dum2 * (vroot - 3. -2. * dum1)
                0867
                0868       rsl1(i) = vgrdrs (ityp (i)) / (rootl * vgrotd (ityp (i)))
                0869       rsl2(i) = alphaf / vgrotd (ityp (i))
                0870
                0871       scat =       agrn(i) *(vgtr11(ityp(i)) + vgrf11(ityp(i)))
                0872      &     + (1. - agrn(i))*(vgtr12(ityp(i)) + vgrf12(ityp(i)))
                0873       sqsc(i) = sqrt (1. - scat)
                0874
                0875       d       = vgdd(k2,ityp(i))*facp + vgdd(k1,ityp(i))*facm
0889f02121 Jean*0876       ufac(i) = log( (vgz2(ityp(i)) - d) / zoch(i) )
d8b3764985 Andr*0877      *        / log(  pblzet             / zoch(i) )
                0878
                0879       z2ch(i) = vgz2(ityp (i))
                0880
                0881       cdsc(i) = pblzet/zoch(i)+1.
                0882       cdrc(i) = vkrm/log(cdsc(i))
                0883       cdrc(i) = cdrc(i)*cdrc(i)
                0884       cdsc(i) = sqrt(cdsc(i))
                0885       cdsc(i) = cdrc(i)*cdsc(i)
                0886
                0887       enddo
                0888
                0889       return
                0890       end
2a3ae9099b Andr*0891
082e38725b Andr*0892       subroutine pkappa (im,jm,lm,ple,pkle,pkz)
2a3ae9099b Andr*0893 C***********************************************************************
                0894 C  Purpose
                0895 C     Calculate Phillips P**Kappa
                0896 C
                0897 C  Arguments
                0898 C     PLE .... edge-level pressure
                0899 C     PKLE ... edge-level pressure**kappa
                0900 C     IM ..... longitude  dimension
                0901 C     JM ..... latitude   dimension
                0902 C     LM ..... vertical   dimension
                0903 C     PKZ .... mid-level pressure**kappa
                0904 C***********************************************************************
                0905       implicit  none
                0906
                0907       integer   im,jm,lm
a456aa407c Andr*0908       _RL ple(im,jm,lm+1)
                0909       _RL pkle(im,jm,lm+1)
                0910       _RL pkz(im,jm,lm)
2a3ae9099b Andr*0911
a456aa407c Andr*0912       _RL akap1,getcon
2a3ae9099b Andr*0913       integer i,j,L
0889f02121 Jean*0914
2a3ae9099b Andr*0915       akap1 = 1.0 + getcon('KAPPA')
0889f02121 Jean*0916
2a3ae9099b Andr*0917       do L = 1,lm
                0918       do j = 1,jm
                0919       do i = 1,im
                0920        pkz(i,j,L) = ( ple (i,j,l+1)*pkle(i,j,l+1)
                0921      .             - ple (i,j,l)*pkle(i,j,l) )
                0922      .           /  ( akap1* (ple (i,j,l+1)-ple (i,j,l)) )
                0923       enddo
                0924       enddo
                0925       enddo
                0926
                0927       return
                0928       end