var/convertor/average_be/module_readwrf.F

   1 module readwrf_module
   2
   3   real*4, parameter :: sclht = 287.04 * 256.0 / 9.81
   4   real*4, parameter :: eps = 0.622
   5   real*4, parameter :: ezero = 6.112
   6   real*4, parameter :: eslcon1 = 17.67
   7   real*4, parameter :: eslcon2 = 29.65
   8   real*4, parameter :: gamma = 287.04/1004.
   9   real*4, parameter :: gammamd = 0.608-0.887
  10   real*4, parameter :: pi = 3.1415926           ! Value used in WRF.
  11   real*4, parameter :: earth_radius = 6370.0    ! Value used in WRF.
  12   real*4, parameter :: deg_to_rad = pi/180.0
  13   real*4, parameter :: rad_to_deg = 1.0/deg_to_rad
  14
  15   real*4 :: ptop
  16   real*4, allocatable, dimension(:) :: znw, znu, znfac
  17   integer, allocatable, dimension(:,:) :: landmask
  18   real*4, allocatable, dimension(:,:) :: wrf_psfc, wrf_t2, wrf_q2, wrf_u10, wrf_v10, wrf_rh2, wrf_xlat, wrf_xlong
  19   real*4, allocatable, dimension(:,:,:) :: wrf_t, wrf_u, wrf_v, wrf_p, wrf_pb, wrf_ph, wrf_phb, wrf_ght, wrf_q, wrf_rh, wrf_tdp
  20
  21   integer :: kdim_stag
  22
  23   contains
  24     function readwrf(filename, west_east_dim, south_north_dim, dx, dy, cen_lat, cen_lon, &
  25                  stand_lon, true1, true2, mapproj, idim, jdim, kdim, idim_stag, jdim_stag, &
  26                  ratio, miycors, mjxcors)
  27
  28     implicit none
  29
  30 #include "netcdf.inc"
  31
  32     ! Arguments
  33     integer :: mapproj, west_east_dim, south_north_dim, idim, jdim, kdim, idim_stag, jdim_stag
  34     real*4 :: dx, dy, cen_lat, cen_lon, stand_lon, true1, true2, ratio, miycors, mjxcors
  35     character (len=*) :: filename
  36
  37     ! Local variables
  38     integer :: readwrf
  39     integer :: status, ncid, n_dimensions, n_variables, n_attributes, unltd_dim_id
  40     integer :: xtype, ndims, n_atts, i, j, xtype_att, l, att_num, dimid
  41     integer :: ii, jj, kk
  42     integer, dimension(NF_MAX_VAR_DIMS) :: dimids
  43     real*4 :: q, e, es, elog, gammam
  44     character (len=NF_MAX_NAME) :: varname, attname, dimname
  45
  46     ! Set for compatibility with
  47     ratio = 1.0
  48     miycors = 1.0
  49     mjxcors = 1.0
  50
  51     status = nf_open(trim(filename), 0, ncid)
  52
  53     if (status /= NF_NOERR) then
  54       write(6,*) 'Error :: readwrf() : Could not open file '//trim(filename)
  55       readwrf = 0
  56       return
  57     end if
  58
  59     status = nf_inq(ncid, n_dimensions, n_variables, n_attributes, unltd_dim_id)
  60
  61     if (status == NF_NOERR) then
  62
  63     do i=1,n_attributes
  64       status = nf_inq_attname(ncid, NF_GLOBAL, i, attname)
  65       if (status == NF_NOERR) then
  66
  67         status = nf_inq_att(ncid, NF_GLOBAL, attname, xtype, l)
  68         if (status == NF_NOERR) then
  69           if (index(trim(attname),'WEST-EAST_GRID_DIMENSION') /= 0) then
  70             status = nf_get_att_int(ncid, NF_GLOBAL, attname, west_east_dim)
  71             if (status /= NF_NOERR) write(6,*) 'Error getting west-east dimension'
  72           else if (index(trim(attname),'SOUTH-NORTH_GRID_DIMENSION') /= 0) then
  73             status = nf_get_att_int(ncid, NF_GLOBAL, attname, south_north_dim)
  74             if (status /= NF_NOERR) write(6,*) 'Error getting south-north dimension'
  75           else if (index(trim(attname),'MAP_PROJ') /= 0) then
  76             status = nf_get_att_int(ncid, NF_GLOBAL, attname, mapproj)
  77             if (status /= NF_NOERR) write(6,*) 'Error getting map_projection'
  78           else if (index(trim(attname),'DX') /= 0) then
  79             status = nf_get_att_real(ncid, NF_GLOBAL, attname, dx)
  80             if (status /= NF_NOERR) write(6,*) 'Error getting dx'
  81 !            dx = dx / 1000.
  82           else if (index(trim(attname),'DY') /= 0) then
  83             status = nf_get_att_real(ncid, NF_GLOBAL, attname, dy)
  84             if (status /= NF_NOERR) write(6,*) 'Error getting dy'
  85 !            dy = dy / 1000.
  86           else if (index(trim(attname),'CEN_LAT') /= 0 .and. len_trim(attname) == 7) then
  87             status = nf_get_att_real(ncid, NF_GLOBAL, attname, cen_lat)
  88             if (status /= NF_NOERR) write(6,*) 'Error getting cen_lat'
  89           else if (index(trim(attname),'CEN_LON') /= 0) then
  90             status = nf_get_att_real(ncid, NF_GLOBAL, attname, cen_lon)
  91             if (status /= NF_NOERR) write(6,*) 'Error getting cen_lon'
  92           else if (index(trim(attname),'STAND_LON') /= 0) then
  93             status = nf_get_att_real(ncid, NF_GLOBAL, attname, stand_lon)
  94             if (status /= NF_NOERR) write(6,*) 'Error getting stand_lon'
  95           else if (index(trim(attname),'TRUELAT1') /= 0) then
  96             status = nf_get_att_real(ncid, NF_GLOBAL, attname, true1)
  97             if (status /= NF_NOERR) write(6,*) 'Error getting true1'
  98           else if (index(trim(attname),'TRUELAT2') /= 0) then
  99             status = nf_get_att_real(ncid, NF_GLOBAL, attname, true2)
 100             if (status /= NF_NOERR) write(6,*) 'Error getting true2'
 101           end if
 102         end if
 103       end if
 104     end do
 105
 106     do i=1,n_variables
 107       status = nf_inq_var(ncid, i, varname, xtype, ndims, dimids, n_atts)
 108
 109       if (status == NF_NOERR) then
 110         if (index(trim(varname),'P_TOP') /= 0) then !{
 111
 112          status = nf_get_var_real(ncid, i, ptop)
 113
 114         else if (index(trim(varname),'LANDMASK') /= 0) then !{
 115
 116           do j=1, ndims
 117             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 118             if (status == NF_NOERR) then
 119               if (index(trim(dimname),'west_east') /= 0) then
 120                 idim = l
 121               else if (index(trim(dimname),'west-east') /= 0) then
 122                 idim = l
 123               else if (index(trim(dimname),'south_north') /= 0) then
 124                 jdim = l
 125               else if (index(trim(dimname),'south-north') /= 0) then
 126                 jdim = l
 127               end if
 128             end if
 129           end do
 130
 131          allocate(landmask(idim, jdim))
 132          status = nf_get_var_real(ncid, i, landmask)
 133
 134         else if (index(trim(varname),'XLAT') /= 0 .and. len_trim(varname) == 4) then !{
 135
 136           do j=1, ndims
 137             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 138             if (status == NF_NOERR) then
 139               if (index(trim(dimname),'west_east') /= 0) then
 140                 idim = l
 141               else if (index(trim(dimname),'west-east') /= 0) then
 142                 idim = l
 143               else if (index(trim(dimname),'south_north') /= 0) then
 144                 jdim = l
 145               else if (index(trim(dimname),'south-north') /= 0) then
 146                 jdim = l
 147               end if
 148             end if
 149           end do
 150
 151          allocate(wrf_xlat(idim, jdim))
 152          status = nf_get_var_real(ncid, i, wrf_xlat)
 153
 154         else if (index(trim(varname),'XLONG') /= 0 .and. len_trim(varname) == 5) then !{
 155
 156           do j=1, ndims
 157             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 158             if (status == NF_NOERR) then
 159               if (index(trim(dimname),'west_east') /= 0) then
 160                 idim = l
 161               else if (index(trim(dimname),'west-east') /= 0) then
 162                 idim = l
 163               else if (index(trim(dimname),'south_north') /= 0) then
 164                 jdim = l
 165               else if (index(trim(dimname),'south-north') /= 0) then
 166                 jdim = l
 167               end if
 168             end if
 169           end do
 170
 171          allocate(wrf_xlong(idim, jdim))
 172          status = nf_get_var_real(ncid, i, wrf_xlong)
 173
 174         else if (index(trim(varname),'T') /= 0 .and. len_trim(varname) == 1) then !{
 175
 176           do j=1, ndims
 177             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 178             if (status == NF_NOERR) then
 179               if (index(trim(dimname),'west_east') /= 0) then
 180                 idim = l
 181               else if (index(trim(dimname),'west-east') /= 0) then
 182                 idim = l
 183               else if (index(trim(dimname),'south_north') /= 0) then
 184                 jdim = l
 185               else if (index(trim(dimname),'south-north') /= 0) then
 186                 jdim = l
 187               else if (index(trim(dimname),'bottom_top') /= 0) then
 188                 kdim = l
 189               else if (index(trim(dimname),'bottom-top') /= 0) then
 190                 kdim = l
 191               end if
 192             end if
 193           end do
 194
 195          allocate(wrf_t(idim, jdim, kdim))
 196          status = nf_get_var_real(ncid, i, wrf_t)
 197
 198         else if (index(trim(varname),'U') /= 0 .and. len_trim(varname) == 1) then !{
 199
 200           do j=1, ndims
 201             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 202             if (status == NF_NOERR) then
 203               if (index(trim(dimname),'west_east_stag') /= 0) then
 204                 idim_stag = l
 205               else if (index(trim(dimname),'west-east_stag') /= 0) then
 206                 idim_stag = l
 207               else if (index(trim(dimname),'south_north') /= 0) then
 208                 jdim = l
 209               else if (index(trim(dimname),'south-north') /= 0) then
 210                 jdim = l
 211               else if (index(trim(dimname),'bottom_top') /= 0) then
 212                 kdim = l
 213               else if (index(trim(dimname),'bottom-top') /= 0) then
 214                 kdim = l
 215               end if
 216             end if
 217           end do
 218
 219          allocate(wrf_u(idim_stag, jdim, kdim))
 220          status = nf_get_var_real(ncid, i, wrf_u)
 221
 222         else if (index(trim(varname),'V') /= 0 .and. len_trim(varname) == 1) then !{
 223
 224           do j=1, ndims
 225             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 226             if (status == NF_NOERR) then
 227               if (index(trim(dimname),'west_east') /= 0) then
 228                 idim = l
 229               else if (index(trim(dimname),'west-east') /= 0) then
 230                 idim = l
 231               else if (index(trim(dimname),'south_north_stag') /= 0) then
 232                 jdim_stag = l
 233               else if (index(trim(dimname),'south-north_stag') /= 0) then
 234                 jdim_stag = l
 235               else if (index(trim(dimname),'bottom_top') /= 0) then
 236                 kdim = l
 237               else if (index(trim(dimname),'bottom-top') /= 0) then
 238                 kdim = l
 239               end if
 240             end if
 241           end do
 242
 243          allocate(wrf_v(idim, jdim_stag, kdim))
 244          status = nf_get_var_real(ncid, i, wrf_v)
 245
 246         else if (index(trim(varname),'PSFC') /= 0 .and. len_trim(varname) == 4) then !{
 247
 248           do j=1, ndims
 249             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 250             if (status == NF_NOERR) then
 251               if (index(trim(dimname),'west_east') /= 0) then
 252                 idim = l
 253               else if (index(trim(dimname),'west-east') /= 0) then
 254                 idim = l
 255               else if (index(trim(dimname),'south_north') /= 0) then
 256                 jdim = l
 257               else if (index(trim(dimname),'south-north') /= 0) then
 258                 jdim = l
 259               end if
 260             end if
 261           end do
 262
 263          allocate(wrf_psfc(idim, jdim))
 264          status = nf_get_var_real(ncid, i, wrf_psfc)
 265
 266         else if (index(trim(varname),'T2') /= 0 .and. len_trim(varname) == 2) then !{
 267
 268           do j=1, ndims
 269             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 270             if (status == NF_NOERR) then
 271               if (index(trim(dimname),'west_east') /= 0) then
 272                 idim = l
 273               else if (index(trim(dimname),'west-east') /= 0) then
 274                 idim = l
 275               else if (index(trim(dimname),'south_north') /= 0) then
 276                 jdim = l
 277               else if (index(trim(dimname),'south-north') /= 0) then
 278                 jdim = l
 279               end if
 280             end if
 281           end do
 282
 283          allocate(wrf_t2(idim, jdim))
 284          status = nf_get_var_real(ncid, i, wrf_t2)
 285
 286         else if (index(trim(varname),'Q2') /= 0 .and. len_trim(varname) == 2) then !{
 287
 288           do j=1, ndims
 289             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 290             if (status == NF_NOERR) then
 291               if (index(trim(dimname),'west_east') /= 0) then
 292                 idim = l
 293               else if (index(trim(dimname),'west-east') /= 0) then
 294                 idim = l
 295               else if (index(trim(dimname),'south_north') /= 0) then
 296                 jdim = l
 297               else if (index(trim(dimname),'south-north') /= 0) then
 298                 jdim = l
 299               end if
 300             end if
 301           end do
 302
 303          allocate(wrf_q2(idim, jdim))
 304          status = nf_get_var_real(ncid, i, wrf_q2)
 305
 306         else if (index(trim(varname),'U10') /= 0 .and. len_trim(varname) == 3) then !{
 307
 308           do j=1, ndims
 309             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 310             if (status == NF_NOERR) then
 311               if (index(trim(dimname),'west_east') /= 0) then
 312                 idim = l
 313               else if (index(trim(dimname),'west-east') /= 0) then
 314                 idim = l
 315               else if (index(trim(dimname),'south_north') /= 0) then
 316                 jdim = l
 317               else if (index(trim(dimname),'south-north') /= 0) then
 318                 jdim = l
 319               end if
 320             end if
 321           end do
 322
 323          allocate(wrf_u10(idim, jdim))
 324          status = nf_get_var_real(ncid, i, wrf_u10)
 325
 326         else if (index(trim(varname),'V10') /= 0 .and. len_trim(varname) == 3) then !{
 327
 328           do j=1, ndims
 329             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 330             if (status == NF_NOERR) then
 331               if (index(trim(dimname),'west_east') /= 0) then
 332                 idim = l
 333               else if (index(trim(dimname),'west-east') /= 0) then
 334                 idim = l
 335               else if (index(trim(dimname),'south_north') /= 0) then
 336                 jdim = l
 337               else if (index(trim(dimname),'south-north') /= 0) then
 338                 jdim = l
 339               end if
 340             end if
 341           end do
 342
 343          allocate(wrf_v10(idim, jdim))
 344          status = nf_get_var_real(ncid, i, wrf_v10)
 345
 346         else if (index(trim(varname),'QVAPOR') /= 0 .and. len_trim(varname) == 6) then !{
 347
 348           do j=1, ndims
 349             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 350             if (status == NF_NOERR) then
 351               if (index(trim(dimname),'west_east') /= 0) then
 352                 idim = l
 353               else if (index(trim(dimname),'west-east') /= 0) then
 354                 idim = l
 355               else if (index(trim(dimname),'south_north') /= 0) then
 356                 jdim = l
 357               else if (index(trim(dimname),'south-north') /= 0) then
 358                 jdim = l
 359               else if (index(trim(dimname),'bottom_top') /= 0) then
 360                 kdim = l
 361               else if (index(trim(dimname),'bottom-top') /= 0) then
 362                 kdim = l
 363               end if
 364             end if
 365           end do
 366
 367          allocate(wrf_q(idim, jdim, kdim))
 368          status = nf_get_var_real(ncid, i, wrf_q)
 369
 370         else if (index(trim(varname),'PB') /= 0 .and. len_trim(varname) == 2) then !{
 371
 372           do j=1, ndims
 373             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 374             if (status == NF_NOERR) then
 375               if (index(trim(dimname),'west_east') /= 0) then
 376                 idim = l
 377               else if (index(trim(dimname),'west-east') /= 0) then
 378                 idim = l
 379               else if (index(trim(dimname),'south_north') /= 0) then
 380                 jdim = l
 381               else if (index(trim(dimname),'south-north') /= 0) then
 382                 jdim = l
 383               else if (index(trim(dimname),'bottom_top') /= 0) then
 384                 kdim = l
 385               else if (index(trim(dimname),'bottom-top') /= 0) then
 386                 kdim = l
 387               end if
 388             end if
 389           end do
 390
 391          allocate(wrf_pb(idim, jdim, kdim))
 392          status = nf_get_var_real(ncid, i, wrf_pb)
 393
 394         else if (index(trim(varname),'P') /= 0 .and. len_trim(varname) == 1) then !{
 395
 396           do j=1, ndims
 397             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 398             if (status == NF_NOERR) then
 399               if (index(trim(dimname),'west_east') /= 0) then
 400                 idim = l
 401               else if (index(trim(dimname),'west-east') /= 0) then
 402                 idim = l
 403               else if (index(trim(dimname),'south_north') /= 0) then
 404                 jdim = l
 405               else if (index(trim(dimname),'south-north') /= 0) then
 406                 jdim = l
 407               else if (index(trim(dimname),'bottom_top') /= 0) then
 408                 kdim = l
 409               else if (index(trim(dimname),'bottom-top') /= 0) then
 410                 kdim = l
 411               end if
 412             end if
 413           end do
 414
 415          allocate(wrf_p(idim, jdim, kdim))
 416          status = nf_get_var_real(ncid, i, wrf_p)
 417
 418         else if (index(trim(varname),'PHB') /= 0 .and. len_trim(varname) == 3) then !{
 419
 420           do j=1, ndims
 421             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 422             if (status == NF_NOERR) then
 423               if (index(trim(dimname),'west_east') /= 0) then
 424                 idim = l
 425               else if (index(trim(dimname),'west-east') /= 0) then
 426                 idim = l
 427               else if (index(trim(dimname),'south_north') /= 0) then
 428                 jdim = l
 429               else if (index(trim(dimname),'south-north') /= 0) then
 430                 jdim = l
 431               else if (index(trim(dimname),'bottom_top_stag') /= 0) then
 432                 kdim_stag = l
 433               else if (index(trim(dimname),'bottom-top_stag') /= 0) then
 434                 kdim_stag = l
 435               end if
 436             end if
 437           end do
 438
 439          allocate(wrf_phb(idim, jdim, kdim_stag))
 440          status = nf_get_var_real(ncid, i, wrf_phb)
 441
 442         else if (index(trim(varname),'PH') /= 0 .and. len_trim(varname) == 2) then !{
 443
 444           do j=1, ndims
 445             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 446             if (status == NF_NOERR) then
 447               if (index(trim(dimname),'west_east') /= 0) then
 448                 idim = l
 449               else if (index(trim(dimname),'west-east') /= 0) then
 450                 idim = l
 451               else if (index(trim(dimname),'south_north') /= 0) then
 452                 jdim = l
 453               else if (index(trim(dimname),'south-north') /= 0) then
 454                 jdim = l
 455               else if (index(trim(dimname),'bottom_top_stag') /= 0) then
 456                 kdim_stag = l
 457               else if (index(trim(dimname),'bottom-top_stag') /= 0) then
 458                 kdim_stag = l
 459               end if
 460             end if
 461           end do
 462
 463          allocate(wrf_ph(idim, jdim, kdim_stag))
 464          status = nf_get_var_real(ncid, i, wrf_ph)
 465
 466         else if (index(trim(varname),'ZNU') /= 0 .and. len_trim(varname) == 3) then !{
 467
 468           do j=1, ndims
 469             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 470             if (status == NF_NOERR) then
 471               if (index(trim(dimname),'bottom_top') /= 0) then
 472                 kdim = l
 473               else if (index(trim(dimname),'bottom-top') /= 0) then
 474                 kdim = l
 475               end if
 476             end if
 477           end do
 478
 479          allocate(znu(kdim))
 480          status = nf_get_var_real(ncid, i, znu)
 481
 482         else if (index(trim(varname),'ZNW') /= 0 .and. len_trim(varname) == 3) then !{
 483
 484           do j=1, ndims
 485             status = nf_inq_dim(ncid, dimids(j), dimname, l)
 486             if (status == NF_NOERR) then
 487               if (index(trim(dimname),'bottom-top_stag') /= 0) then
 488                 kdim_stag = l
 489               else if (index(trim(dimname),'bottom_top_stag') /= 0) then
 490                 kdim_stag = l
 491               end if
 492             end if
 493           end do
 494
 495          allocate(znw(kdim_stag))
 496          status = nf_get_var_real(ncid, i, znw)
 497
 498         end if !}
 499
 500       end if
 501     end do
 502
 503     end if
 504
 505     readwrf = 1
 506
 507     allocate(znfac(kdim))
 508
 509     do kk=1,kdim
 510        znfac(kk)=(znw(kk)-znu(kk))/(znw(kk)-znw(kk+1))
 511     enddo
 512
 513     do kk=1,kdim_stag
 514       do jj=1,jdim
 515         do ii=1,idim
 516           wrf_ph(ii,jj,kk)=exp(-(wrf_ph(ii,jj,kk)+wrf_phb(ii,jj,kk))/(9.81*sclht))
 517         end do
 518       end do
 519     end do
 520
 521     allocate(wrf_ght(idim, jdim, kdim))
 522
 523     do kk=1,kdim
 524       do jj=1,jdim
 525         do ii=1,idim
 526           wrf_ght(ii,jj,kk)=znfac(kk)*wrf_ph(ii,jj,kk+1)+(1.-znfac(kk))*wrf_ph(ii,jj,kk)
 527           wrf_ght(ii,jj,kk)=-sclht*log(wrf_ght(ii,jj,kk))
 528         end do
 529       end do
 530     end do
 531
 532 !    do kk=1,kdim
 533 !      do jj=1,jdim
 534 !        do ii=1,idim
 535 !           wrf_p(ii,jj,kk)=wrf_p(ii,jj,kk) + wrf_pb(ii,jj,kk)
 536 !        end do
 537 !      end do
 538 !    end do
 539
 540 !    do kk=1,kdim
 541 !      do jj=1,jdim
 542 !        do ii=1,idim
 543 !         gammam=gamma*(1.+gammamd*.001*wrf_q(ii,jj,kk))
 544 !         wrf_t(ii,jj,kk)=(wrf_t(ii,jj,kk)+300.)*(wrf_p(ii,jj,kk)/100000.)**gammam
 545 !        end do
 546 !      end do
 547 !    end do
 548
 549     allocate(wrf_rh(idim, jdim, kdim))
 550
 551 !    do kk=1,kdim
 552 !      do jj=1,jdim
 553 !        do ii=1,idim
 554 !           e = wrf_q(ii,jj,kk)*(wrf_p(ii,jj,kk)/100.)/(eps+wrf_q(ii,jj,kk))
 555 !           es = ezero * exp( eslcon1*(wrf_t(ii,jj,kk)-273.15)/(wrf_t(ii,jj,kk)-eslcon2) )
 556 !           wrf_rh(ii,jj,kk)=100.*(e*((wrf_p(ii,jj,kk)/100.)-es))/(es*((wrf_p(ii,jj,kk)/100.)-e))
 557 !        end do
 558 !      end do
 559 !    end do
 560
 561     allocate(wrf_tdp(idim, jdim, kdim))
 562
 563 !    do kk=1,kdim
 564 !      do jj=1,jdim
 565 !        do ii=1,idim
 566 !           q=max(wrf_q(ii,jj,kk),1.e-15)
 567 !           e=q*(wrf_p(ii,jj,kk)/100.)/(eps+q)
 568 !           elog=alog(e/ezero)
 569 !           wrf_tdp(ii,jj,kk)=(eslcon2*elog-eslcon1*273.15)/(elog-eslcon1)
 570 !        end do
 571 !      end do
 572 !    end do
 573
 574     allocate(wrf_rh2(idim, jdim))
 575
 576     do jj=1,jdim
 577       do ii=1,idim
 578          e = wrf_q2(ii,jj)*(wrf_psfc(ii,jj)/100.)/(eps+wrf_q2(ii,jj))
 579          es = ezero * exp( eslcon1*(wrf_t2(ii,jj)-273.15)/(wrf_t2(ii,jj)-eslcon2) )
 580          wrf_rh2(ii,jj)=100.*(e*((wrf_psfc(ii,jj)/100.)-es))/(es*((wrf_psfc(ii,jj)/100.)-e))
 581       end do
 582     end do
 583
 584     end function readwrf
 585
 586 end module readwrf_module