var/convertor/wave2grid_kma/PRESUB.inc

   1       SUBROUTINE CRH2SH( WV, T, PS, A, B, IJMAX, LMAX, LARHM,
   2      \                   QMIN, ITERMX )
   3
   4       DIMENSION
   5      1        WV(IJMAX,LMAX), T(IJMAX,LMAX), PS(IJMAX), A(LMAX), B(LMAX)
   6
   7 C     INPUT  WV : RELATIVE HUMIDITY( 0-1 )   ( L=<LARHM ) : CHANGED
   8 C               : SPECIFIC HUMIDITY( KG/KG ) ( L >LARHM ) : UNCHANGED
   9 C            T  : VIRTUAL TEMPERATURE(K)                  : CHANGED
  10 C            PS : SURFACE PRESSURE(HPA)
  11 C     OUTPUT WV : SPECIFIC HUMIDITY( KG/KG ) ( FOR ALL LEVELS )
  12 C            T  : REAL TEMPERATURE(K)
  13
  14       EE = 1.D0/(0.608D0*t_kelvin)
  15       FF = 1.D0/(0.608D0*35.9D0)
  16       GG = 6.11D0*0.622D0
  17       CC = t_kelvin/35.9D0*7.5D0*LOG(10.D0)
  18
  19       DO 100 IJ = 1, IJMAX
  20
  21         DO 110 L = 1, LARHM
  22
  23           VT = T (IJ,L)
  24           RH = WV(IJ,L)
  25           IF( L.LE.LMAX-1 ) THEN
  26             PHF1 = A(L  )+B(L  )*PS(IJ)
  27             PHF2 = A(L+1)+B(L+1)*PS(IJ)
  28             PL   = EXP( ( PHF1*LOG(PHF1)-PHF2*LOG(PHF2) )
  29      \                    /( PHF1-PHF2 )-1.0 )
  30           ELSE
  31             PL   = 0.5*( A(L)+B(L)*PS(IJ) )
  32           END IF
  33
  34           AA = GG*RH/PL
  35           BB = EE*(VT-t_kelvin)
  36           DD = FF*(VT-35.9)
  37           HH = CC*(BB-DD)
  38
  39 C         ŒJ‚è•Ô‚µŒvŽZŠJŽn : QQ = INITIAL VALUE OF SPECIFIC HUMIDITY
  40
  41           QQ = AA*EXP(CC*BB/DD)
  42
  43           DO 120 ITER = 1, ITERMX
  44             QV  = AA*EXP(CC*(QQ-BB)/(QQ-DD))
  45             FQ  = QV-QQ/(1.0+0.608*QQ)
  46             DFQ = HH*QV/((QQ-DD)**2)-1.0/(1.0+0.608*QQ)**2
  47             QQ  = QQ-FQ/DFQ
  48             IF( QQ.LT.QMIN ) QQ = QMIN
  49   120     CONTINUE
  50
  51 C         ŒJ‚è•Ô‚µŒvŽZ�I—¹�¨ŽÀ‰·“x‚ÌŒvŽZ
  52
  53           TT = VT/(1.0+0.608*QQ)
  54           T (IJ,L) = TT
  55           WV(IJ,L) = QQ
  56
  57   110   CONTINUE
  58
  59         IF( LARHM.LT.LMAX ) THEN
  60           DO 130 L = LARHM+1, LMAX
  61             T(IJ,L) = T(IJ,L)/(1.0+0.608*WV(IJ,L))
  62   130     CONTINUE
  63         ENDIF
  64
  65   100 CONTINUE
  66 C
  67       RETURN
  68       END SUBROUTINE CRH2SH
  69 C**********************************************************************
  70       SUBROUTINE CRH2SHA
  71      I(IJMAX,KMAX,PS,A,B,GRAV,GASR,TLAPS,QCONS,QMIN,KST,ITERMX,
  72 *-->
  73      I IDX, LARHM,
  74 *<--
  75      O WV,T)
  76 C
  77 C***  CALCULATE SPECIFIC HUMIDITY AND REAL TEMPERATURE
  78 C   FROM RELATIVE HUMIDITY ,REAL TEMPERATURE AND VIRTUAL TEMPERATURE
  79 C   CORRECTION
  80 C     SPECIFIC HUMIDITY = CONSTANT IN STRATOSPHERE
  81 C
  82 C*** (ARRAYS)    (INPUT) ******************* (OUTPUT) ******************
  83 C     WV         RELATIVE HUMIDITY(NOT %)    SPECIFIC HUMIDITY(KG/KG)
  84 C     T          VIRTUAL TEMPERATURE(K)      REAL TEMPERATURE(K)
  85 C     PS         SURFACE PRESSURE (MB)       (UNCHANGED)
  86 C     A,B        A+B*PS DEFINES INTER-LAYER LEVEL
  87 C***********************************************************************
  88       DIMENSION WV(IJMAX,KMAX),T(IJMAX,KMAX),PS(IJMAX),A(KMAX),B(KMAX)
  89       DATA Q90/2.7E-6/   ! SAME AS IN WV300M
  90       COMMON/CTETEN/TABLE(25000)
  91       COMMON/DTETEN/DTABLE(25000)
  92       REAL*8 TABLE,DTABLE,X
  93       LOGICAL FIRST/.TRUE./
  94 C
  95       IF(FIRST) THEN
  96           GBYR=GRAV/GASR
  97       ENDIF
  98 C
  99       DO 100 K=1,KMAX
 100       DO 100 I=  1,IJMAX
 101       VT=T  (I,K)
 102       IF( K.LT.LARHM ) THEN
 103       RH=WV (I,K)
 104       ENDIF
 105 C
 106       IF(K.LE.KMAX-1) THEN
 107       PHF1=A(K  )+B(K  )*PS(I)
 108       PHF2=A(K+1)+B(K+1)*PS(I)
 109       PL  =EXP( (PHF1*LOG(PHF1)-PHF2*LOG(PHF2))/(PHF1-PHF2)-1.0 )
 110       ELSE
 111       PL  =0.5*(A(K)+B(K)*PS(I))
 112       END IF
 113 C
 114 C
 115           T0 = VT
 116        YI = (T0-123.2D0)*100.0D0
 117        IY = YI
 118        IY = MAX( IY, 1 )                                     ! 96/01/31
 119        IY = MIN( IY, 24999 )                                 ! 96/01/31
 120        X = YI - IY
 121        QSAT = ((1.0D0-X)*TABLE(IY)+X*TABLE(IY+1))/PL
 122 CMM     IF((IDX.EQ.1).AND.(PL.LE.90.)) THEN
 123 CMM       QQ=MIN(Q90,QSAT*0.9)
 124 CMM     ELSE
 125       IF( K.LT.LARHM ) THEN
 126        DQSAT = ((1.0D0-X)*DTABLE(IY)+X*DTABLE(IY+1))/PL
 127        FT = VT-T0-0.608*T0*QSAT*RH
 128        DFT = -1.0-0.608*QSAT*RH-0.608*T0*DQSAT*RH
 129        T0 = T0 - FT/DFT
 130 C
 131        YI = (T0-123.2D0)*100.0D0
 132        IY = YI
 133        IY = MAX( IY, 1 )                                    ! 96/01/31
 134        IY = MIN( IY, 24999 )                                ! 96/01/31
 135        X = YI - IY
 136        QSAT = ((1.0D0-X)*TABLE(IY)+X*TABLE(IY+1))/PL
 137        DQSAT = ((1.0D0-X)*DTABLE(IY)+X*DTABLE(IY+1))/PL
 138        FT = VT-T0-0.608*T0*QSAT*RH
 139        DFT = -1.0-0.608*QSAT*RH-0.608*T0*DQSAT*RH
 140        T0 = T0 - FT/DFT
 141 C
 142        YI = (T0-123.2D0)*100.0D0
 143        IY = YI
 144        IY = MAX( IY, 1 )                                    ! 96/01/31
 145        IY = MIN( IY, 24999 )                                 ! 96/01/31
 146        X = YI - IY
 147        QQ = ((1.0D0-X)*TABLE(IY)+X*TABLE(IY+1))/PL*RH
 148       ELSE
 149        QQ = WV(I,K)
 150       ENDIF
 151 C
 152       IF((IDX.EQ.1).AND.(PL.LE.90.)) THEN
 153         QQ=MIN(QQ,QSAT*0.9)
 154       ENDIF
 155 C
 156 C +++ CALCULATE REAL TEMPERATURE +++
 157 C
 158           TT=VT/(1.+0.608*QQ)
 159 C
 160       T (I,K)=TT
 161       WV(I,K)=QQ
 162 C
 163   100 CONTINUE
 164 C
 165       RETURN
 166       END SUBROUTINE CRH2SHA
 167 C ---------------------------------------------------------------------
 168       SUBROUTINE SPLDIF( ZOUT, POUT, LMXOUT, ZIN, PIN, LMXIN )
 169
 170       DIMENSION ZOUT(LMXOUT), POUT(LMXOUT), ZIN(LMXIN), PIN(LMXIN)
 171       DIMENSION SM(40), H(40), AL(40), AM(40), AP(40), C(40)
 172
 173 C     INPUT / ZIN (L), PIN (L), LMXIN  : INPUT.DATA, PRES(LOG), NUMBER
 174 C     OUTPUT/ ZOUT(L), POUT(L), LMXOUT : OUTPUT-VAL, PRES(LOG), NUMBER
 175
 176       LM1 = LMXIN-1
 177       GR = -gravity/gas_constant
 178
 179       DO 110 L = 2, LMXIN
 180         H(L) = PIN(L)-PIN(L-1)
 181   110 CONTINUE
 182
 183       DO 120 L = 2, LM1
 184         AL(L) = 0.5*H(L+1)/(H(L)+H(L+1))
 185         AM(L) = 0.5-AL(L)
 186   120 CONTINUE
 187
 188 C     �I’[�ðŒ�( LAPSE RATE IS CONSTANT )
 189 C     SM(1) = SM(2) ; SM(LMXIN-1) = SM(LMXIN) : SECOND DERIVATIVE
 190
 191       AL(1)     = -1.0
 192       AM(LMXIN) = -1.0
 193       AL(LMXIN) =  0.0
 194
 195       DO 130 L = 2, LMXIN
 196         AP(L) = 1.0/(1.0-AL(L-1)*AM(L))
 197         AL(L) = AL(L)*AP(L)
 198   130 CONTINUE
 199
 200       C(1)     = 0.0
 201       C(LMXIN) = 0.0
 202       DO 160 L = 2, LM1
 203         C(L) = 3.0*((ZIN(L+1)-ZIN(L))/H(L+1) - (ZIN(L)-ZIN(L-1))/H(L))
 204      \            /(H(L)+H(L+1))
 205   160 CONTINUE
 206
 207 C     FORWARD SUBSTITUTION
 208
 209       DO 200 L = 2, LMXIN
 210         C(L) = (C(L)-C(L-1)*AM(L))*AP(L)
 211   200 CONTINUE
 212       SM(LMXIN) = C(LMXIN)
 213
 214 C     BACKWARD SUBSTUTUTION
 215
 216       DO 220 K = 1, LM1
 217         L = LMXIN-K
 218         SM(L) = C(L)-AL(L)*SM(L+1)
 219   220 CONTINUE
 220
 221 C     INTERPOLATION
 222
 223       LB = 2
 224
 225       DO 500 LOUT = 1, LMXOUT
 226
 227         X = POUT(LOUT)
 228         DO 300 L = LB, LMXIN
 229           IF( X.GE.PIN(L) ) GO TO 310
 230   300   CONTINUE
 231         L = LMXIN
 232   310   LB = L
 233
 234 C       ‚RŽŸƒXƒvƒ‰ƒCƒ“ŠÖ�”‚Ì”÷•ª
 235
 236         ZOUT(LOUT) = SM(L-1)*(  -(PIN(L)-X)**2 /(2.0*H(L))+H(L)/6.0 )
 237      \             + SM(L)  *( (X-PIN(L-1))**2 /(2.0*H(L))-H(L)/6.0 )
 238      \             + ( ZIN(L)-ZIN(L-1) )/H(L)
 239         ZOUT(LOUT) = ZOUT(LOUT)*GR
 240
 241   500 CONTINUE
 242
 243       RETURN
 244       END SUBROUTINE SPLDIF
 245 C ---------------------------------------------------------------------
 246       SUBROUTINE SPLDIF3( ZOUT, POUT, LMXOUT, ZIN, PIN, LMXIN, IJMAX,
 247      W SM, H, AL, AM, AP, C )
 248       DIMENSION ZOUT(IJMAX,LMXOUT), POUT(IJMAX,LMXOUT),
 249      1          ZIN (IJMAX,LMXIN),  PIN (IJMAX,LMXIN)
 250 CMM   DIMENSION ZOUT(LMXOUT), POUT(LMXOUT), ZIN(LMXIN), PIN(LMXIN)
 251       DIMENSION SM(IJMAX,LMXIN), H(IJMAX,LMXIN), AL(IJMAX,LMXIN),
 252      1          AM(IJMAX,LMXIN), AP(IJMAX,LMXIN), C(IJMAX,LMXIN)
 253 CMM   DIMENSION SM(40), H(40), AL(40), AM(40), AP(40), C(40)
 254
 255 C     INPUT / ZIN (L), PIN (L), LMXIN  : INPUT.DATA, PRES(LOG), NUMBER
 256 C     OUTPUT/ ZOUT(L), POUT(L), LMXOUT : OUTPUT-VAL, PRES(LOG), NUMBER
 257
 258       LM1 = LMXIN-1
 259       GR = -gravity/gas_constant
 260
 261       DO 110 L = 2, LMXIN
 262       DO 110 I = 1, IJMAX
 263         H(I,L) = PIN(I,L)-PIN(I,L-1)
 264   110 CONTINUE
 265
 266       DO 120 L = 2, LM1
 267       DO 120 I = 1, IJMAX
 268         AL(I,L) = 0.5*H(I,L+1)/(H(I,L)+H(I,L+1))
 269         AM(I,L) = 0.5-AL(I,L)
 270   120 CONTINUE
 271
 272 C     �I’[�ðŒ�( LAPSE RATE IS CONSTANT )
 273 C     SM(1) = SM(2) ; SM(LMXIN-1) = SM(LMXIN) : SECOND DERIVATIVE
 274       DO 125 I = 1, IJMAX
 275         AL(I,1)     = -1.0
 276         AM(I,LMXIN) = -1.0
 277         AL(I,LMXIN) =  0.0
 278   125 CONTINUE
 279
 280       DO 130 L = 2, LMXIN
 281       DO 130 I = 1, IJMAX
 282 C
 283         AP(I,L) = 1.0/(1.0-AL(I,L-1)*AM(I,L))
 284         AL(I,L) = AL(I,L)*AP(I,L)
 285   130 CONTINUE
 286
 287       DO 155 I = 1, IJMAX
 288       C(I,1)     = 0.0
 289       C(I,LMXIN) = 0.0
 290   155 CONTINUE
 291 C
 292       DO 160 L = 2, LM1
 293       DO 160 I = 1, IJMAX
 294         C(I,L) = 3.0*((ZIN(I,L+1)-ZIN(I,L))/H(I,L+1)
 295      1         - (ZIN(I,L)-ZIN(I,L-1))/H(I,L))
 296      2            /(H(I,L)+H(I,L+1))
 297   160 CONTINUE
 298
 299 C     FORWARD SUBSTITUTION
 300
 301       DO 200 L = 2, LMXIN
 302       DO 200 I = 1, IJMAX
 303         C(I,L) = (C(I,L)-C(I,L-1)*AM(I,L))*AP(I,L)
 304   200 CONTINUE
 305       DO 205 I = 1, IJMAX
 306       SM(I,LMXIN) = C(I,LMXIN)
 307   205 CONTINUE
 308
 309 C     BACKWARD SUBSTUTUTION
 310
 311       DO 220 K = 1, LM1
 312       DO 220 I = 1, IJMAX
 313         L = LMXIN-K
 314         SM(I,L) = C(I,L)-AL(I,L)*SM(I,L+1)
 315   220 CONTINUE
 316
 317 C     INTERPOLATION
 318
 319 C     LB = 2
 320
 321       DO 500 LOUT = 1, LMXOUT
 322       DO 500 I = 1, IJMAX
 323
 324         X = POUT(I,LOUT)
 325         L = LOUT ! FOR ONLY PIN.EQ.POUT
 326         IF( L.LT.2 ) L=2
 327 CM      DO 300 L = LB, LMXIN
 328 CM        IF( X.GE.PIN(L) ) GO TO 310
 329 CM300   CONTINUE
 330 CM      L = LMXIN
 331 CM310   LB = L
 332
 333         ZOUT(I,LOUT) = SM(I,L-1)*(  -(PIN(I,L)-X)**2
 334      1             /(2.0*H(I,L))+H(I,L)/6.0 )
 335      2             + SM(I,L)  *( (X-PIN(I,L-1))**2
 336      3             /(2.0*H(I,L))-H(I,L)/6.0 )
 337      \             + ( ZIN(I,L)-ZIN(I,L-1) )/H(I,L)
 338         ZOUT(I,LOUT) = ZOUT(I,LOUT)*GR
 339
 340   500 CONTINUE
 341
 342       RETURN
 343       END SUBROUTINE SPLDIF3
 344 C======================================================================
 345       SUBROUTINE SETWHT
 346      1(IMAX,JMAX,DP,IP,DCOSCL,DGW,COCLT)
 347       DIMENSION DP(4,IMAX,JMAX)
 348       INTEGER*2 IP(2,IMAX,JMAX)
 349       REAL*8  DCOSCL(JMAX),DGW(JMAX),COCLT(JMAX)
 350 C
 351 Crizvi     Already defined in module_wave2grid_kma
 352 C      JMAXHF=JMAX/2
 353       DELX=360.0/FLOAT(IMAX)
 354       CALL GAUSS(DCOSCL,DGW,JMAX)
 355 C
 356       R2D=180.0/pi
 357       DO 100 J=1,JMAXHF
 358       COCLT(       J)=      R2D*ACOS(DCOSCL(J))
 359       COCLT(JMAX+1-J)=180.0-R2D*ACOS(DCOSCL(J))
 360   100 CONTINUE
 361 C     WRITE(6,*) 'CHECK OF COLATTITUDE',COCLT
 362 C
 363       DO 120 J=1,JMAX
 364       PX=COCLT(J)+1.5
 365       LAT=PX
 366       DLAT=PX-FLOAT(LAT)
 367       DO 120 I=1,IMAX
 368       QX=DELX*FLOAT(I-1)+1.5
 369       LON=QX
 370       DLON=QX-FLOAT(LON)
 371          IP(1,I,J)=LON
 372          IP(2,I,J)=LAT
 373          DP(1,I,J)=(1.0-DLAT)*(1.0-DLON)
 374          DP(2,I,J)=(1.0-DLAT)*DLON
 375          DP(3,I,J)=DLAT*DLON
 376          DP(4,I,J)=DLAT*(1.0-DLON)
 377   120 CONTINUE
 378 C     DO 1000 J=1,JMAX
 379 C     WRITE(6,1010) J,IP(1,1,J),IP(2,1,J),(DP(I,1,J),I=1,4)
 380 C1000 CONTINUE
 381 C     DO 1100 J=1,JMAX
 382 C     WRITE(6,1010) J,IP(1,IMAX,J),IP(2,IMAX,J),(DP(I,IMAX,J),I=1,4)
 383 C1100 CONTINUE
 384 C     DO 1200 J=1,IMAX
 385 C     WRITE(6,1020) J,IP(1,J,1),IP(2,J,1),(DP(I,J,1),I=1,4)
 386 C1200 CONTINUE
 387 C     DO 1300 J=1,IMAX
 388 C     WRITE(6,1020) J,IP(1,J,JMAX),IP(2,J,JMAX),(DP(I,J,JMAX),I=1,4)
 389 C1300 CONTINUE
 390 C1010 FORMAT(1H ,'J=',I3,5X,2I5,5X,4F10.3)
 391 C1020 FORMAT(1H ,'I=',I3,5X,2I5,5X,4F10.3)
 392 C
 393       RETURN
 394       END SUBROUTINE SETWHT
 395       SUBROUTINE INTERP
 396      1(WORK,DATA,IMAX,JMAX,DP,IP)
 397 C
 398       DIMENSION WORK(362,182),DATA(IMAX,JMAX),DP(4,IMAX,JMAX)
 399       INTEGER*2 IP(2,IMAX,JMAX)
 400 C
 401       DO 100 J=1,JMAX
 402       DO 100 I=1,IMAX
 403        II=IP(1,I,J)
 404        JJ=IP(2,I,J)
 405        DATA(I,J)=DP(1,I,J)*WORK(II,JJ  )+DP(2,I,J)*WORK(II+1,JJ  )
 406      *          +DP(4,I,J)*WORK(II,JJ+1)+DP(3,I,J)*WORK(II+1,JJ+1)
 407   100 CONTINUE
 408 C
 409       RETURN
 410       END SUBROUTINE INTERP