var/convertor/wave2grid_kma/pvchkdv.F

   1 CSHCO PROGRAM CPVFULL
   2       SUBROUTINE CPVFULL !SHCN
   3       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
   4       PARAMETER( IMAX=1,JMAX=1,KMAX=3 )
   5       DIMENSION GPV(IMAX,JMAX,KMAX),GPS(IMAX,JMAX),GPS9(IMAX,JMAX),
   6      .          A(KMAX),B(KMAX),GPS8(IMAX,JMAX),GPS1(IMAX,JMAX),
   7      .          GPV9(IMAX,JMAX,KMAX),GPV8(IMAX,JMAX,KMAX)
   8       DIMENSION GPVC9(IMAX,JMAX,KMAX)
   9 C : eta-level coefficients
  10       A(1) = 0.01
  11       A(2) = 0.021
  12       A(3) = 0.03
  13       B(1) = 1.0
  14       B(2) = 0.95
  15       B(3) = 0.8
  16 C : base state
  17       GPS9(1,1)= 1013.
  18 C : input variable
  19       GPS(1,1) = GPS9(1,1)*0.01
  20 C : original log P on full model levels
  21       JSTA=1
  22       JFIN=1
  23       CALL PVFULL
  24      I  (GPS9,IMAX,JMAX,KMAX,A,B,
  25      I   JSTA,JFIN,
  26      O   GPV9)
  27 C     WRITE(6,*) ' GPV9=',(DEXP(GPV9(1,1,K)),K=1,KMAX)
  28       CALL PVFULLC
  29      I  (GPS9,IMAX,JMAX,KMAX,A,B,
  30      I   JSTA,JFIN,
  31      O   GPVC9)
  32 C : tangent code
  33       CALL TPVFULL
  34      I  (GPS,IMAX,JMAX,KMAX,
  35      I   GPVC9,
  36      I   JSTA,JFIN,
  37      O   GPV)
  38 C : input variable
  39 C     DO 8000 N=1,10
  40       DO 8000 N=-4,8
  41         ALFA = 10.D0**(-DFLOAT(N)/2.D0)
  42         GPS8(1,1) = GPS9(1,1)+GPS(1,1)*ALFA
  43 C : original code
  44       CALL PVFULL
  45      I  (GPS8,IMAX,JMAX,KMAX,A,B,
  46      I   JSTA,JFIN,
  47      O   GPV8)
  48 C : incremental check
  49         DOT1 = 0.D0
  50         DOT2 = 0.D0
  51         DOT  = 0.D0
  52           DOT1 = DOT1 + (GPS8(1,1)-GPS9(1,1))**2
  53           DOT2 = DOT2 +  GPS(1,1)**2
  54           DOT  = DOT  + (GPS8(1,1)-GPS9(1,1))*GPS(1,1)
  55 C         write(6,*) ' GPS DOT1,DOT2,DOT=',DOT1,DOT2,DOT
  56         DO 8100 K=1,KMAX
  57           DOT1V=        (GPV8(1,1,K)-GPV9(1,1,K))**2.D0
  58      . *1.D5**2
  59           DOT2V=         GPV(1,1,K)**2
  60      . *1.D5**2
  61           DOTV =        (GPV8(1,1,K)-GPV9(1,1,K))*GPV(1,1,K)
  62      . *1.D5**2
  63           DOT1 = DOT1 + DOT1V
  64           DOT2 = DOT2 + DOT2V
  65           DOT  = DOT  + DOTV
  66 C         write(6,*) ' GPV DOT1,DOT2,DOT=',DOT1V,DOT2V,DOTV
  67  8100   CONTINUE
  68           WRITE(6,*) ' GSI DEV=',N,DOT/DSQRT(DOT1*DOT2),ALFA*GPS(1,1)
  69  8000 CONTINUE
  70 C : left-hand side calculation (inner product of tagent code output)
  71       RLEFT = 0.
  72       RLEFT = RLEFT + GPS(1,1)*GPS(1,1)
  73       DO 3100 K=1,KMAX
  74       RLEFT = RLEFT + GPV(1,1,K)*GPV(1,1,K)
  75  3100 CONTINUE
  76       WRITE(6,*) ' RLEFT=',RLEFT
  77 C : adjoint code
  78       GPS1(1,1) = GPS(1,1)
  79       CALL APVFULL
  80      I  (GPS1,IMAX,JMAX,KMAX,
  81      I   GPVC9,
  82      I   JSTA,JFIN,
  83      O   GPV)
  84       WRITE(6,*)
  85       WRITE(6,*)
  86       WRITE(6,*) ' AGPV=',GPV
  87       WRITE(6,*) ' GPS9=',GPS9
  88       WRITE(6,*) ' AGPS=',GPS1
  89 C : right-hand side calculation
  90       RIGHT= 0.
  91       RIGHT = RIGHT + GPS(1,1)*GPS1(1,1)
  92       WRITE(6,*) ' LEFT,RIGHT,DEV=',RLEFT,RIGHT,RLEFT-RIGHT
  93       STOP
  94       END
  95       SUBROUTINE PVFULL
  96      I  (GPS,IMAX,JMAX,KMAX,A,B,
  97      I   JSTA,JFIN,
  98      O   GPV)
  99 C**********************************************************************
 100 C full-level (L) pressure (PV) (HPA) calcualtion (log P)
 101 C                                            2000.01.19 Y.TAKEUCHI
 102 C<INPUT>
 103 C GPS(IMAX,JMAX): surface pressure (hPa)
 104 C<OUTPUT>
 105 C GPV(IMAX,JMAX,KMAX): full-level log P (NOUNIT)
 106 C**********************************************************************
 107       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 108       DIMENSION GPS(IMAX,JMAX),GPV(IMAX,JMAX,KMAX)
 109       REAL*8    A(50), B(50)
 110 C******************** PROCEDURE ***************************************
 111 C : half-level (L+1/2,L-1/2) pressure (HPA) calculation
 112 C : parallelization, modified for vectorization
 113 C     IJFACT = 8
 114 C     IJFACT = 1
 115 CSHCO*POPTION PARALLEL
 116 CSHCO*POPTION TLOCAL(I,J,K,PU,PD)
 117 CSHCO*POPTION INDEP(GPV)
 118       DO 160 J=JSTA,JFIN
 119 C     DO 160 J1= 1,JMAX/IJFACT
 120 CSHCO*VOPTION INDEP
 121 *vdir novector !SHCN
 122       DO 161 I=1,IMAX
 123 C     DO 161 I1= 1,IMAX*IJFACT
 124 C       I = MOD(I1-1,IMAX)+1
 125 C       J = (J1-1)*IJFACT+(I1-1)/IMAX+1
 126 C       write(6,*) ' I,J,I1,J1=',I,J,I1,J1
 127 *vdir novector  !SHCN
 128       DO 165 K = 1,KMAX-1
 129 C : half-level (L+1/2,L-1/2) pressure (HPA) calculation
 130           PU  = A(K+1) + B(K+1)*GPS(I,J)
 131           PD  = A(K  ) + B(K  )*GPS(I,J)
 132 c         if (pd.eq.0.0.or.pu.eq.0.0) write(999,*) i,j,pd,pu  !SHCN
 133 C : full-level (L) pressure (HPA) calculation
 134           GPV(I,J,K) = ( PD*DLOG(PD)-PU*DLOG(PU) )/(PD-PU) -1.D0
 135 C       if(I.EQ.1.AND.J.EQ.1) write(6,*) ' GPV',K,PU,PD,GPV(I,J,K)
 136 C         GPV(I,J,K) = GPV(I,J,K)*1000.D0
 137   165 CONTINUE
 138           GPV(I,J,KMAX) = DLOG((A(KMAX)+B(KMAX)*GPS(I,J))/2.D0)
 139 C         GPV(I,J,KMAX) = GPV(I,J,KMAX)*1000.D0
 140 C       if(I.EQ.1.AND.J.EQ.1) write(6,*) ' GPV',KMAX,PU,PD,GPV(I,J,KMAX)
 141   161 CONTINUE
 142   160 CONTINUE
 143       RETURN
 144       END
 145       SUBROUTINE TPVFULL
 146      I  (GPS,IMAX,JMAX,KMAX,
 147      I   GPVC9,
 148      I   JSTA,JFIN,
 149      O   GPV)
 150 C**********************************************************************
 151 C : PV: full-level (L) pressure (HPA) calculation  GPV DELETE
 152 C : tangent code
 153 C : log P
 154 C                                            2000.01.19 Y.TAKEUCHI
 155 C<INPUT>
 156 C GPS(IMAX,JMAX): surface pressure increment (hPa)
 157 C GPVC9(IMAX,JMAX,KMAX): coefficient
 158 C<OUTPUT>
 159 C GPV(IMAX,JMAX,KMAX): full-level log P increment (NOUNIT)
 160 C**********************************************************************
 161       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 162       DIMENSION GPS(IMAX,JMAX),GPVC9(IMAX,JMAX,KMAX)
 163       DIMENSION GPV(IMAX,JMAX,KMAX)
 164 C******************** PROCEDURE ***************************************
 165 C : parallelization, modified for vectorization
 166 C     IJFACT = 8
 167 C     IJFACT = 1
 168 *POPTION PARALLEL
 169 *POPTION TLOCAL(I,J,K)
 170 *POPTION INDEP(GPVC9)
 171       DO 960 J= JSTA,JFIN
 172 C     DO 960 J1= 1,JMAX/IJFACT
 173 *VOPTION INDEP
 174       DO 961 I = 1,IMAX
 175 C     DO 961 I1= 1,IMAX*IJFACT
 176 C       I = MOD(I1-1,IMAX)+1
 177 C       J = (J1-1)*IJFACT+(I1-1)/IMAX+1
 178
 179 C     DO 960 J = 1,JMAX
 180 C     DO 961 I = 1,IMAX
 181
 182           GPV(I,J,KMAX) = GPS(I,J) * GPVC9(I,J,KMAX)
 183   961 CONTINUE
 184       DO 965 K = 1,KMAX-1
 185       DO 966 I = 1,IMAX
 186           GPV(I,J,K) = GPS(I,J)*GPVC9(I,J,K)
 187   966 CONTINUE
 188   965 CONTINUE
 189   960 CONTINUE
 190       RETURN
 191       END
 192       SUBROUTINE APVFULL
 193      I  (GPS,IMAX,JMAX,KMAX,
 194      I   GPVC9,
 195      I   JSTA,JFIN,
 196      O   GPV)
 197 C**********************************************************************
 198 C : PV: full-level (L) pressure (HPA) calculation  GPV DELETE
 199 C : verification of the adjoint code TDVAR.FORT(PVFULL)
 200 C : log P
 201 C                                            2000.01.19 Y.TAKEUCHI
 202 C**********************************************************************
 203       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 204       DIMENSION GPS(IMAX,JMAX),GPV(IMAX,JMAX,KMAX),
 205      .          GPVC9(IMAX,JMAX,KMAX)
 206 C******************** PROCEDURE ***************************************
 207 C2000.09.11
 208 *POPTION PARALLEL,PRIND((J,1))
 209 *POPTION TLOCAL(I,J,K)
 210       DO 960 J = JSTA,JFIN
 211       DO 961 I = 1,IMAX
 212           GPS(I,J) = GPS(I,J)+GPV(I,J,KMAX)*GPVC9(I,J,KMAX)
 213   961 CONTINUE
 214       DO 965 K = 1,KMAX-1
 215       DO 966 I = 1,IMAX
 216           GPS(I,J) = GPS(I,J) +GPV(I,J,K)*GPVC9(I,J,K)
 217   966 CONTINUE
 218   965 CONTINUE
 219   960 CONTINUE
 220       RETURN
 221       END
 222       SUBROUTINE PVFULLC
 223      I  (GPS9,IMAX,JMAX,KMAX,A,B,
 224      I   JSTA,JFIN,
 225      O   GPVC9)
 226 C**********************************************************************
 227 C : coefficients for full-level pressure calculation
 228 C : log P
 229 C                                            2000.01.19 Y.TAKEUCHI
 230 C<INPUT>
 231 C GPS9(IMAX,JMAX): base state of surface pressure (hPa)
 232 C<OUTPUT>
 233 C GPVC9(IMAX,JMAX,KMAX):
 234 C**********************************************************************
 235       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 236       DIMENSION GPS9(IMAX,JMAX),GPVC9(IMAX,JMAX,KMAX)
 237       REAL*8    A(50), B(50)
 238 C******************** PROCEDURE ***************************************
 239 C : half-level (L+1/2,L-1/2) pressure (HPA) calculation
 240
 241 C : parallelization, modified for vectorization
 242 C     IJFACT = 8
 243 C     IJFACT = 1
 244 *POPTION PARALLEL
 245 *POPTION TLOCAL(I,J,K,PU9,PD9)
 246 *POPTION INDEP(GPVC9)
 247       DO 960 J=JSTA,JFIN
 248 C     DO 960 J1= 1,JMAX/IJFACT
 249 *VOPTION INDEP
 250       DO 961 I=1,IMAX
 251 C     DO 961 I1= 1,IMAX*IJFACT
 252 C       I = MOD(I1-1,IMAX)+1
 253 C       J = (J1-1)*IJFACT+(I1-1)/IMAX+1
 254
 255 C     DO 960 J = 1,JMAX
 256 C     DO 961 I = 1,IMAX
 257
 258           GPVC9(I,J,KMAX) =
 259      .       B(KMAX)/(A(KMAX)+B(KMAX)*GPS9(I,J))
 260       DO 965 K = 1,KMAX-1
 261           PU9 = A(K+1) + B(K+1)*GPS9(I,J)
 262           PD9 = A(K  ) + B(K  )*GPS9(I,J)
 263
 264           GPVC9(I,J,K) = B(K  )*
 265      .     1.D0/(PD9-PU9)**2*(-PU9*(DLOG(PD9)-DLOG(PU9))+PD9-PU9)
 266      .                 + B(K+1)*
 267      .     1.D0/(PD9-PU9)**2*(PD9*(DLOG(PD9)-DLOG(PU9))-PD9+PU9)
 268   965 CONTINUE
 269   961 CONTINUE
 270   960 CONTINUE
 271       RETURN
 272       END