ungrib/src/rrpr.F

   1 subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_format, prefix)
   2 !                                                                             !
   3 ! In case you are wondering, RRPR stands for "Read, ReProcess, and wRite"     !
   4 !                                                                             !
   5 !*****************************************************************************!
   6 !                                                                             !
   7
   8   use filelist
   9   use gridinfo
  10   use storage_module
  11   use table
  12   use module_debug
  13   use misc_definitions_module
  14   use stringutil
  15
  16   implicit none
  17
  18 !------------------------------------------------------------------------------
  19 ! Arguments:
  20
  21 ! HSTART:  Starting date of times to process
  22   character (LEN=19) :: hstart
  23
  24 ! NTIMES:  Number of time periods to process
  25   integer :: ntimes
  26
  27 ! INTERVAL:  Time inteval (seconds) of time periods to process.
  28   integer :: interval
  29
  30 ! NLVL:  The number of levels in the stored data.
  31   integer :: nlvl
  32
  33 ! MAXLVL: The parameterized maximum number of levels to allow.
  34   integer :: maxlvl
  35
  36 ! PLVL:  Array of pressure levels (Pa) in the dataset
  37   real , dimension(maxlvl) :: plvl
  38
  39 ! DEBUG_LEVEL:  Integer level of debug printing (from namelist)
  40   integer :: debug_level
  41
  42 !------------------------------------------------------------------------------
  43
  44   character (LEN=25) :: units
  45   character (LEN=46) :: Desc
  46   real, allocatable, dimension(:,:) :: scr2d, tmp2d
  47   real, pointer, dimension(:,:) :: ptr2d
  48
  49   integer :: k, kk, mm, n, ierr, ifv
  50   integer :: iunit=13
  51
  52   character(LEN=19) :: hdate, hend
  53   character(LEN=24) :: hdate_output
  54   character(LEN=3)  :: out_format
  55   character(LEN=MAX_FILENAME_LEN)  :: prefix
  56   real :: xfcst, level
  57   character(LEN=9) :: field
  58
  59   integer :: ntime, idts
  60
  61 ! DATELEN:  length of date strings to use for our output file names.
  62   integer :: datelen
  63
  64 ! Decide the length of date strings to use for output file names.
  65 ! DATELEN is 13 for hours, 16 for minutes, and 19 for seconds.
  66
  67   if (mod(interval,3600) == 0) then
  68      datelen = 13
  69   else if (mod(interval, 60) == 0) then
  70      datelen = 16
  71   else
  72      datelen = 19
  73   endif
  74
  75   if ( debug_level .gt. 100 ) then
  76     call mprintf(.true.,DEBUG,"Begin rrpr")
  77     call mprintf(.true.,DEBUG,"nfiles = %i , ntimes = %i )",i1=nfiles,i2=ntimes)
  78     do n = 1, nfiles
  79       call mprintf(.true.,DEBUG,"filedates(%i) = %s",i1=n,s1=filedates(n))
  80     enddo
  81   endif
  82
  83 ! Compute the ending time:
  84
  85   call geth_newdate(hend, hstart, interval*ntimes)
  86
  87   call clear_storage
  88
  89 ! We want to do something for each of the requested times:
  90   TIMELOOP : do ntime = 1, ntimes
  91      idts = (ntime-1) * interval
  92      call geth_newdate(hdate, hstart, idts)
  93      call mprintf(.true.,DEBUG, &
  94      "RRPR: hstart = %s , hdate = %s , idts = %i",s1=hstart,s2=hdate,i1=idts)
  95
  96 ! Loop over the output file dates, and do stuff if the file date matches
  97 ! the requested time we are working on now.
  98
  99      FILELOOP : do n = 1, nfiles
 100        if ( debug_level .gt. 100 ) then
 101          call mprintf(.true.,DEBUG, &
 102             "hstart = %s , hend = %s",s1=hstart,s2=hend)
 103          call mprintf(.true.,DEBUG, &
 104             "filedates(n) = %s",s1=filedates(n))
 105          call mprintf(.true.,DEBUG, &
 106             "filedates(n) = %s",s1=filedates(n)(1:datelen))
 107        end if
 108        if (filedates(n)(1:datelen).ne.hdate(1:datelen)) cycle FILELOOP
 109        if (debug_level .gt. 50 ) then
 110          call mprintf(.true.,INFORM, &
 111             "RRPR Processing : %s",s1=filedates(n)(1:datelen))
 112        endif
 113        open(iunit, file=trim(get_path(prefix))//'PFILE:'//filedates(n)(1:datelen), &
 114           form='unformatted',status='old')
 115
 116 ! Read the file:
 117
 118      rdloop: do
 119         read (iunit, iostat=ierr) ifv
 120         if (ierr.ne.0) exit rdloop
 121         if ( ifv .eq. 5) then     ! WPS
 122           read (iunit) hdate_output, xfcst, map%source, field, units, Desc, &
 123                level, map%nx, map%ny, map%igrid
 124           hdate = hdate_output(1:19)
 125           select case (map%igrid)
 126           case (0, 4)
 127              read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%r_earth
 128           case (3)
 129            read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 130                 map%truelat1, map%truelat2, map%r_earth
 131           case (5)
 132              read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 133                 map%truelat1, map%r_earth
 134           case (1)
 135            read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
 136                 map%truelat1, map%r_earth
 137           case default
 138              call mprintf(.true.,ERROR, &
 139                 "Unrecognized map%%igrid: %i in RRPR 1",i1=map%igrid)
 140           end select
 141           read (iunit) map%grid_wind
 142
 143         else if ( ifv .eq. 4 ) then          ! SI
 144           read (iunit) hdate_output, xfcst, map%source, field, units, desc, level, &
 145                 map%nx, map%ny, map%igrid
 146           hdate = hdate_output(1:19)
 147           select case (map%igrid)
 148           case (0, 4)
 149              read(iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx
 150           case (3)
 151              read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 152                 map%lov, map%truelat1, map%truelat2
 153           case (5)
 154              read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 155                 map%lov, map%truelat1
 156           case default
 157              call mprintf(.true.,ERROR, &
 158                 "Unrecognized map%%igrid: %i in RRPR 2",i1=map%igrid)
 159           end select
 160
 161         else if ( ifv .eq. 3 ) then          ! MM5
 162           read(iunit) hdate_output, xfcst, field, units, desc, level,&
 163                 map%nx, map%ny, map%igrid
 164           hdate = hdate_output(1:19)
 165           select case (map%igrid)
 166           case (3)      ! lamcon
 167             read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 168                     map%truelat1, map%truelat2
 169            case (5)      ! Polar Stereographic
 170               read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 171                    map%truelat1
 172            case (0, 4)      ! lat/lon
 173               read (iunit) map%lat1, map%lon1, map%dy, map%dx
 174            case (1)      ! Mercator
 175               read (iunit) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
 176            case default
 177              call mprintf(.true.,ERROR, &
 178                 "Unrecognized map%%igrid: %i in RRPR 3",i1=map%igrid)
 179            end select
 180         else
 181            call mprintf(.true.,ERROR, &
 182               "unknown out_format, ifv = %i",i1=ifv)
 183         endif
 184
 185         allocate(ptr2d(map%nx,map%ny))
 186         read (iunit) ptr2d
 187         call refw_storage(nint(level), field, ptr2d, map%nx, map%ny)
 188         nullify (ptr2d)
 189      enddo rdloop
 190
 191    write (0,*) 'Name of source model =>',map%source
 192 !
 193 ! We have reached the end of file, so time to close it.
 194 !
 195      close(iunit)
 196      if (debug_level .gt. 100 ) call print_storage
 197 !
 198 ! By now the file has been read completely.  Now, see if we need to fill in
 199 ! missing fields:
 200 !
 201
 202 ! Retrieve the number of levels in storage:
 203 !
 204      call get_plvls(plvl, maxlvl, nlvl)
 205 !
 206 ! Fill the surface level (code 200100) from higher 200100s, as necessary
 207 !
 208         do k = 1, nlvl
 209            if ((plvl(k).gt.200100) .and. (plvl(k).lt.200200)) then
 210            ! We found a level between 200100 and 200200, now find the field
 211            ! corresponding to that level.
 212               MLOOP : do mm = 1, maxvar
 213                  if (is_there(nint(plvl(k)), namvar(mm))) then
 214                     INLOOP : do kk = 200101, nint(plvl(k))
 215                        if (is_there(kk, namvar(mm))) then
 216                           if ( debug_level .gt. 100 ) then
 217                             call mprintf(.true.,DEBUG, &
 218                "Copying %s at level %i to level 200100.",s1=namvar(mm),i1=kk)
 219                           end if
 220                           call get_dims(kk, namvar(mm))
 221                           allocate(scr2d(map%nx,map%ny))
 222                           call get_storage &
 223                                (kk, namvar(mm), scr2d, map%nx, map%ny)
 224                           call put_storage &
 225                                (200100,namvar(mm), scr2d,map%nx,map%ny)
 226                           deallocate(scr2d)
 227                           EXIT INLOOP
 228                        endif
 229                     enddo INLOOP
 230                  endif
 231               enddo MLOOP
 232            endif
 233         enddo
 234
 235 !
 236 ! If upper-air U is missing, see if we can interpolate from surrounding levels.
 237 ! This is a simple vertical interpolation, linear in pressure.
 238 ! Currently, this simply fills in one missing level between two present levels.
 239 !
 240
 241         do k = 2, nlvl-1, 1
 242            if (plvl(k-1) .lt. 200000.) then
 243               if ( (.not. is_there(nint(plvl(k)),'UU')) .and. &
 244                    ( is_there(nint(plvl(k-1)), 'UU')) .and.&
 245                    ( is_there(nint(plvl(k+1)), 'UU')) ) then
 246                  call get_dims(nint(plvl(k+1)), 'UU')
 247                  call vntrp(plvl, maxlvl, k, "UU      ", map%nx, map%ny)
 248               endif
 249            endif
 250         enddo
 251
 252 !
 253 ! If upper-air V is missing, see if we can interpolate from surrounding levels.
 254 ! This is a simple vertical interpolation, linear in pressure.
 255 ! Currently, this simply fills in one missing level between two present levels.
 256 !
 257
 258         do k = 2, nlvl-1, 1
 259            if (plvl(k-1) .lt. 200000.) then
 260               if ( (.not. is_there(nint(plvl(k)),'VV')) .and. &
 261                    ( is_there(nint(plvl(k-1)), 'VV')) .and.&
 262                    ( is_there(nint(plvl(k+1)), 'VV')) ) then
 263                  call get_dims(nint(plvl(k+1)), 'VV')
 264                  call vntrp(plvl, maxlvl, k, "VV      ", map%nx, map%ny)
 265               endif
 266            endif
 267         enddo
 268
 269 !
 270 ! If upper-air SPECHUMD is missing, see if we can compute SPECHUMD from QVAPOR:
 271 !--- Tanya's change for initializing WRF with RUC
 272
 273         do k = 1, nlvl
 274            if (plvl(k).lt.200000.) then
 275               if (.not. is_there(nint(plvl(k)), 'SPECHUMD').and. &
 276                    is_there(nint(plvl(k)), 'QV')) then
 277                  call get_dims(nint(plvl(k)), 'QV')
 278                  call compute_spechumd_qvapor(map%nx, map%ny, plvl(k))
 279               endif
 280            endif
 281         enddo
 282
 283 !--- Tanya's change for initializing WRF with RUC
 284 !   This allows for the ingestion for RUC isentropic data
 285 !
 286         do k = 1, nlvl
 287            if (plvl(k).lt.200000.) then
 288               if (.not. is_there(nint(plvl(k)), 'TT').and. &
 289                    is_there(nint(plvl(k)), 'VPTMP').and. &
 290                    is_there(nint(plvl(k)), 'SPECHUMD')) then
 291                  call get_dims(nint(plvl(k)), 'VPTMP')
 292                  call compute_t_vptmp(map%nx, map%ny, plvl(k))
 293               endif
 294            endif
 295         enddo
 296 !!!
 297 !
 298 ! If upper-air T is missing, see if we can interpolate from surrounding levels.
 299 ! This is a simple vertical interpolation, linear in pressure.
 300 ! Currently, this simply fills in one missing level between two present levels.
 301 !
 302
 303         do k = 2, nlvl-1, 1
 304            if (plvl(k-1) .lt. 200000.) then
 305               if ( (.not. is_there(nint(plvl(k)),'TT')) .and. &
 306                    ( is_there(nint(plvl(k-1)), 'TT')) .and.&
 307                    ( is_there(nint(plvl(k+1)), 'TT')) ) then
 308                  call get_dims(nint(plvl(k+1)), 'TT')
 309                  call vntrp(plvl, maxlvl, k, "TT      ", map%nx, map%ny)
 310               endif
 311            endif
 312         enddo
 313
 314 !
 315 ! Check to see if we need to fill HGT from GEOPT.
 316 !
 317         do k = 1, nlvl
 318            if (plvl(k).lt.200000.) then
 319               if (.not. is_there(nint(plvl(k)), 'HGT').and. &
 320                    is_there(nint(plvl(k)), 'GEOPT')) then
 321                  call get_dims(nint(plvl(k)), 'GEOPT')
 322                  allocate(scr2d(map%nx,map%ny))
 323                  call get_storage(nint(plvl(k)), 'GEOPT', scr2d, map%nx, map%ny)
 324                  scr2d = scr2d / 9.81
 325                  call put_storage(nint(plvl(k)), 'HGT',   scr2d, map%nx, map%ny)
 326                  call mprintf(.true.,DEBUG, &
 327                    "RRPR:   Computing GHT from GEOPT ")
 328                  deallocate(scr2d)
 329               endif
 330            endif
 331         enddo
 332
 333 !
 334 ! If this is GFS data, we might have data at the level of max wind speed,
 335 ! or the level of the tropopause.  If so, we want to replicate the pressures
 336 ! at those levels (new names).  The replicated names are to allow the
 337 ! metgrid program to interpolate the 2d pressure array with both a nearest
 338 ! neighbor AND a 4-pt technique.  Those two pressures are used in ARW real
 339 ! for vertical interpolation of the trop and max wind level data.
 340 !
 341
 342         if (index(map%source,'NCEP GFS') .ne. 0 ) then
 343           call mprintf(.true.,DEBUG, &
 344              "RRPR:   Replicating GFS pressures for max wind and trop")
 345             if ( is_there(200100,'PMAXW  ') .or. &
 346                  is_there(200100,'PTROP  ') ) then
 347               call gfs_trop_maxw_pressures (map%nx, map%ny)
 348             endif
 349         endif
 350
 351 ! Repair GFS and ECMWF pressure-level RH
 352         if (index(map%source,'NCEP GFS') .ne. 0 .or.  &
 353             index(map%source,'NCEP CDAS CFSV2') .ne. 0 .or.  &
 354             index(map%source,'ECMWF') .ne. 0 ) then
 355           call mprintf(.true.,DEBUG, &
 356              "RRPR:   Adjusting RH values ")
 357           do k = 1, nlvl
 358             if ( is_there(nint(plvl(k)),'RH') .and. &
 359                  is_there(nint(plvl(k)),'TT') ) then
 360               call fix_gfs_rh (map%nx, map%ny, plvl(k))
 361             endif
 362           enddo
 363         endif
 364
 365 ! If upper-air RH is missing, see if we can compute RH from Specific Humidity:
 366
 367         do k = 1, nlvl
 368            if (plvl(k).lt.200000.) then
 369               if (.not. is_there(nint(plvl(k)), 'RH') .and. &
 370                    is_there(nint(plvl(k)), 'TT') .and. &
 371                    is_there(nint(plvl(k)), 'SPECHUMD')) then
 372                  call get_dims(nint(plvl(k)), 'TT')
 373                  call compute_rh_spechumd_upa(map%nx, map%ny, plvl(k))
 374               endif
 375            endif
 376         enddo
 377
 378 ! If upper-air RH is missing, see if we can compute RH from Vapor Pressure:
 379 !   (Thanks to Bob Hart of PSU ESSC -- 1999-05-27.)
 380
 381         do k = 1, nlvl
 382            if (plvl(k).lt.200000.) then
 383               if (.not. is_there(nint(plvl(k)),'RH').and. &
 384                    is_there(nint(plvl(k)), 'TT') .and. &
 385                    is_there(nint(plvl(k)),'VAPP')) then
 386                  call get_dims(nint(plvl(k)),'TT')
 387                  call compute_rh_vapp_upa(map%nx, map%ny, plvl(k))
 388               endif
 389            endif
 390         enddo
 391
 392 ! If upper-air RH is missing, see if we can compute RH from Dewpoint Depression:
 393
 394         do k = 1, nlvl
 395            if (plvl(k).lt.200000.) then
 396               if (.not. is_there(nint(plvl(k)),'RH').and. &
 397                    is_there(nint(plvl(k)), 'TT') .and. &
 398                    is_there(nint(plvl(k)),'DEPR')) then
 399                  call get_dims(nint(plvl(k)),'TT')
 400                  call compute_rh_depr(map%nx, map%ny, plvl(k))
 401               endif
 402            endif
 403         enddo
 404 !
 405 ! If upper-air RH is missing, see if we can interpolate from surrounding levels.
 406 ! This is a simple vertical interpolation, linear in pressure.
 407 ! Currently, this simply fills in one missing level between two present levels.
 408 ! May expand this in the future to fill in additional levels.  May also expand
 409 ! this in the future to vertically interpolate other variables.
 410 !
 411
 412         do k = 2, nlvl-1, 1
 413            if (plvl(k-1) .lt. 200000.) then
 414               if ( (.not. is_there(nint(plvl(k)),'RH')) .and. &
 415                    ( is_there(nint(plvl(k-1)), 'RH')) .and.&
 416                    ( is_there(nint(plvl(k+1)), 'RH')) ) then
 417                  call get_dims(nint(plvl(k+1)), 'RH')
 418                  call vntrp(plvl, maxlvl, k, "RH      ", map%nx, map%ny)
 419               endif
 420            endif
 421         enddo
 422
 423 !
 424 ! Check to see if we need to fill RH above 300 mb:
 425 !
 426         if (is_there(30000, 'RH')) then
 427            call get_dims(30000, 'RH')
 428            allocate(scr2d(map%nx,map%ny))
 429
 430            do k = 1, nlvl
 431 !   Set missing RH to 5% between 300 and 70 hPa. Set RH to 0 above 70 hPa.
 432 !   The stratospheric RH will be adjusted further in real.
 433               if (plvl(k).le.7000.) then
 434                 scr2d = 0.
 435               else if (plvl(k).lt.30000.) then
 436                 scr2d = 5.
 437               endif
 438               if (plvl(k).lt.30000. .and. plvl(k) .gt. 10. ) then
 439               ! levels higher than .1 mb are special - do not fill
 440                  if (.not. is_there(nint(plvl(k)), 'RH')) then
 441                     call put_storage(nint(plvl(k)),'RH',scr2d,map%nx,map%ny)
 442                     call mprintf(.true.,DEBUG, &
 443                  "RRPR:   RH missing at %i hPa, inserting synthetic RH ",i1=nint(plvl(k)/100.))
 444                  endif
 445               endif
 446            enddo
 447            deallocate(scr2d)
 448         endif
 449 !
 450 ! If surface RH is missing, see if we can compute RH from Specific Humidity
 451 ! or Dewpoint or Dewpoint depression:
 452 !
 453         if (.not. is_there (200100, 'RH')) then
 454            if (is_there(200100, 'TT').and. &
 455                 is_there(200100, 'PSFC'    )   .and. &
 456                 is_there(200100, 'SPECHUMD')) then
 457               call get_dims(200100, 'TT')
 458               call compute_rh_spechumd(map%nx, map%ny)
 459               call mprintf(.true.,DEBUG, &
 460                 "RRPR:   SURFACE RH is computed")
 461            elseif (is_there(200100, 'TT'       ).and. &
 462                 is_there(200100, 'DEWPT')) then
 463               call get_dims(200100, 'TT')
 464               call compute_rh_dewpt(map%nx, map%ny)
 465            elseif (is_there(200100, 'TT').and. &
 466                 is_there(200100, 'DEPR')) then
 467               call get_dims(200100, 'TT')
 468               call compute_rh_depr(map%nx, map%ny, 200100.)
 469            endif
 470         endif
 471
 472 !
 473 ! If surface SNOW is missing, see if we can compute SNOW from SNOWRUC
 474 ! (From Wei Wang, 2007 June 21, modified 12/28/2007)
 475 !
 476         if (.not. is_there(200100, 'SNOW') .and. &
 477              is_there(200100, 'SNOWRUC')) then
 478            call get_dims(200100, 'SNOWRUC')
 479            allocate(scr2d(map%nx,map%ny))
 480            call get_storage(200100, 'SNOWRUC', scr2d, map%nx, map%ny)
 481            scr2d = scr2d * 1000.
 482            call put_storage(200100, 'SNOW',   scr2d, map%nx, map%ny)
 483            deallocate(scr2d)
 484         endif
 485
 486 ! compute snow water equivalent (SNOW) for NCEP RUC  models
 487 ! As of Sept. 14  2011
 488         if ( index(map%source,'NCEP RUC Model') .ne. 0) then
 489           if (is_there(200100, 'SNOWH') .and. .not. is_there(200100, 'SNOW')) then
 490           call get_dims(200100, 'SNOWH')
 491           allocate(scr2d(map%nx,map%ny))
 492           call get_storage(200100, 'SNOWH', scr2d, map%nx, map%ny)
 493           call mprintf(.true.,DEBUG, &
 494              "RRPR:   Computing SNOWH from SNOW")
 495             if (is_there(200100, 'RHOSN')) then        ! If we have snow density, use it to compute snowh
 496               call get_dims(200100, 'RHOSN')
 497               allocate(tmp2d(map%nx,map%ny))
 498               call get_storage(200100, 'RHOSN', tmp2d, map%nx, map%ny)
 499               scr2d = scr2d * tmp2d
 500               deallocate(tmp2d)
 501             else
 502               scr2d = scr2d * 200.0          ! Assume 200:1 ratio
 503             endif
 504           call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
 505           deallocate(scr2d)
 506           endif
 507         endif
 508
 509 ! Add residual soil moisture to SOILM* if initialized from the GSD RUC model or from NCEP RUC
 510       if (index(map%source,'NOAA GSD') .ne. 0 .or.    &
 511           index(map%source,'NCEP RUC Model') .ne. 0) then
 512             if ( .not. is_there(200100, 'SOILM000') .and.&
 513                        is_there(200100, 'SM000ruc') ) then
 514            call get_dims(200100, 'SM000ruc')
 515              print *,'Adjust RUC soil moisture'
 516           call mprintf(.true.,DEBUG, &
 517              "RRPR:   Adjusting RUC soil moisture ")
 518               call fix_ruc_soilm (map%nx, map%ny)
 519             endif
 520        endif
 521
 522 !
 523 ! Check to see if we need to fill SOILHGT from SOILGEO.
 524 ! (From Wei Wang, 2007 June 21)
 525 !
 526         if (.not. is_there(200100, 'SOILHGT') .and. &
 527              is_there(200100, 'SOILGEO')) then
 528            call get_dims(200100, 'SOILGEO')
 529            allocate(scr2d(map%nx,map%ny))
 530            call get_storage(200100, 'SOILGEO', scr2d, map%nx, map%ny)
 531            scr2d = scr2d / 9.81
 532            call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny)
 533            call mprintf(.true.,DEBUG, &
 534              "RRPR:   Computing SOILGHT from SOILGEO ")
 535            deallocate(scr2d)
 536         endif
 537
 538 ! For hybrid-level input, soilgeo is in level 1 (e.g. ERA40)
 539         if (.not. is_there(200100, 'SOILHGT') .and. &
 540              is_there(1, 'SOILGEO')) then
 541            call get_dims(1, 'SOILGEO')
 542            allocate(scr2d(map%nx,map%ny))
 543            call get_storage(1, 'SOILGEO', scr2d, map%nx, map%ny)
 544            scr2d = scr2d / 9.81
 545            call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny)
 546            deallocate(scr2d)
 547         endif
 548
 549 ! For NCEP RR (using the same ID as for RUC) native-level input,
 550 ! may need to move PSFC from level 1 to 2001.
 551 ! From TGS 8 Sept. 2011
 552         if ( index(map%source,'NCEP RUC Model') .ne. 0) then
 553         if (.not. is_there(200100, 'PSFC') .and. &
 554              is_there(1, 'PRESSURE')) then
 555     print *,'Process PSFC for NCEP RR'
 556            call get_dims(1, 'PRESSURE')
 557            allocate(scr2d(map%nx,map%ny))
 558            call get_storage(1, 'PRESSURE', scr2d, map%nx, map%ny)
 559            call put_storage(200100, 'PSFC', scr2d, map%nx, map%ny)
 560            deallocate(scr2d)
 561         endif
 562         endif
 563
 564 ! For ECMWF hybrid-level input, may need to move psfc from level 1 to 2001.
 565         if ( index(map%source,'ECMWF') .ne. 0) then
 566         if (.not. is_there(200100, 'PSFC') .and. &
 567              is_there(1, 'PSFCH')) then
 568            call get_dims(1, 'PSFCH')
 569            allocate(scr2d(map%nx,map%ny))
 570            call get_storage(1, 'PSFCH', scr2d, map%nx, map%ny)
 571            call put_storage(200100, 'PSFC', scr2d, map%nx, map%ny)
 572            deallocate(scr2d)
 573         endif
 574         endif
 575
 576 ! ECMWF snow depth in meters of water equivalent (Table 128). Convert to kg/m2
 577 !
 578         if (is_there(200100, 'SNOW_EC')) then
 579            call get_dims(200100, 'SNOW_EC')
 580            allocate(scr2d(map%nx,map%ny))
 581            call get_storage(200100, 'SNOW_EC', scr2d, map%nx, map%ny)
 582            scr2d = scr2d * 1000.
 583            call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
 584            deallocate(scr2d)
 585         endif
 586
 587 ! Convert the ECMWF LANDSEA mask from a fraction to a flag
 588
 589         if ( index(map%source,'ECMWF') .ne. 0) then
 590         if (is_there(200100, 'LANDSEA')) then
 591            call get_dims(200100, 'LANDSEA')
 592            call make_zero_or_one(map%nx, map%ny, 'LANDSEA')
 593         endif
 594         endif
 595
 596 ! NCEP GFS weasd is one-half of the NAM value. Increase it for use in WRF.
 597 ! The GFS-based reanalyses values should be OK as is.
 598         if ((index(map%source,'NCEP GFS') .ne. 0 .or. &
 599             index(map%source,'NCEP GEFS') .ne. 0) .and. &
 600             is_there(200100, 'SNOW')) then
 601            call mprintf(.true.,DEBUG, &
 602               "RRPR:   Recomputing SNOW for NCEP GFS")
 603            call get_dims(200100, 'SNOW')
 604            allocate(scr2d(map%nx,map%ny))
 605            call get_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
 606            scr2d = scr2d * 2.
 607            call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
 608            deallocate(scr2d)
 609         endif
 610
 611 ! compute physical snow depth (SNOWH) for various models
 612 ! As of March 2011, this is done here instead of real because we have model
 613 ! source information.
 614         if (is_there(200100, 'SNOW') .and. .not. is_there(200100, 'SNOWH')) then
 615           call get_dims(200100, 'SNOW')
 616           allocate(scr2d(map%nx,map%ny))
 617           call get_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
 618           call mprintf(.true.,DEBUG, &
 619              "RRPR:   Computing SNOWH from SNOW")
 620           if ( index(map%source,'NCEP ') .ne. 0) then
 621             scr2d = scr2d * 0.005          ! Assume 200:1 ratio as used at NCEP and in NOAH
 622           else if (index(map%source,'ECMWF') .ne. 0) then
 623             if (is_there(200100, 'SNOW_DEN')) then        ! If we have snow density, use it to compute snowh
 624               call get_dims(200100, 'SNOW_DEN')
 625               allocate(tmp2d(map%nx,map%ny))
 626               call get_storage(200100, 'SNOW_DEN', tmp2d, map%nx, map%ny)
 627               scr2d = scr2d / tmp2d
 628               deallocate(tmp2d)
 629             else
 630               scr2d = scr2d * 0.004     ! otherwise, assume a density of 250 mm/m (i.e. 250:1 ratio).
 631             endif
 632           else                       ! Other models
 633             scr2d = scr2d * 0.005    ! Use real's default method (200:1)
 634           endif
 635           call put_storage(200100, 'SNOWH', scr2d, map%nx, map%ny)
 636           deallocate(scr2d)
 637         endif
 638
 639 ! As of March 2011, SEAICE can be a flag or a fraction. It will be converted
 640 ! to the appropriate values in real depending on whether or not the polar mods are used.
 641
 642 !! If we've got a SEAICE field, make sure that it is all Zeros and Ones:
 643
 644 !       if (is_there(200100, 'SEAICE')) then
 645 !          call get_dims(200100, 'SEAICE')
 646 !          call make_zero_or_one(map%nx, map%ny, 'SEAICE')
 647 !       endif
 648
 649 ! If we've got an ICEMASK field, re-flag it for output to met_em and real:
 650 !     Field  | GRIB In  |  Out
 651 !    -------------------------
 652 !    water   |    0     |  0
 653 !    land    |   -1     |  1
 654 !    ice     |    1     |  0
 655
 656         if (is_there(200100, 'ICEMASK')) then
 657            call get_dims(200100, 'ICEMASK')
 658            call re_flag_ice_mask(map%nx, map%ny)
 659         endif
 660
 661 ! If we have an ICEFRAC field, convert from % to fraction
 662         if (is_there(200100, 'ICEFRAC')) then
 663            call get_dims(200100, 'ICEFRAC')
 664            allocate(scr2d(map%nx,map%ny))
 665            call get_storage(200100, 'ICEFRAC', scr2d, map%nx, map%ny)
 666            scr2d = scr2d / 100.
 667            call put_storage(200100, 'ICEFRAC', scr2d, map%nx, map%ny)
 668            deallocate(scr2d)
 669         endif
 670
 671
 672         call mprintf(.true.,INFORM, &
 673            "RRPR: hdate = %s ",s1=hdate)
 674         call output(hdate, nlvl, maxlvl, plvl, interval, 2, out_format, prefix, debug_level)
 675         call clear_storage
 676         exit FILELOOP
 677      enddo FILELOOP
 678    enddo TIMELOOP
 679 end subroutine rrpr
 680
 681 subroutine make_zero_or_one(ix, jx, infield)
 682 ! Make sure the input field (SEAICE or LANDSEA) is zero or one.
 683   use storage_module
 684   implicit none
 685   integer :: ix, jx
 686   real, dimension(ix,jx) :: seaice
 687   character(len=*) :: infield
 688
 689   call get_storage(200100, infield, seaice, ix, jx)
 690   where(seaice > 0.5)
 691      seaice = 1.0
 692   elsewhere
 693      seaice = 0.0
 694   end where
 695   call put_storage(200100, infield, seaice, ix, jx)
 696 end subroutine make_zero_or_one
 697
 698 subroutine re_flag_ice_mask(ix, jx)
 699 !
 700 ! Change land points from -1 to 1
 701 ! Change ice  points from  1 to 0
 702 ! Water       points stay    at 0
 703 !
 704   use storage_module
 705   implicit none
 706   integer :: ix, jx
 707   real, dimension(ix,jx) :: iceflag
 708
 709   call get_storage(200100, 'ICEMASK',iceflag, ix, jx)
 710   where(iceflag > 0.5)     ! Ice points, set to water value
 711      iceflag = 0.0
 712   end where
 713   where(iceflag < -0.5)    ! Land points
 714      iceflag = 1.0
 715   end where
 716   call put_storage(200100, 'ICEMASK',iceflag, ix, jx)
 717 end subroutine re_flag_ice_mask
 718
 719 subroutine compute_spechumd_qvapor(ix, jx, plvl)
 720 ! Compute specific humidity from water vapor mixing ratio.
 721   use storage_module
 722   implicit none
 723   integer :: ix, jx
 724   real :: plvl
 725   real, dimension(ix,jx) :: QVAPOR, SPECHUMD
 726
 727   call get_storage(nint(plvl), 'QV', QVAPOR, ix, jx)
 728
 729   SPECHUMD = QVAPOR/(1.+QVAPOR)
 730
 731   call put_storage(nint(plvl), 'SPECHUMD', spechumd, ix, jx)
 732  if(nint(plvl).eq.1) then
 733   call put_storage(200100,'SPECHUMD', spechumd, ix, jx)
 734  endif
 735
 736 end subroutine compute_spechumd_qvapor
 737
 738 subroutine compute_t_vptmp(ix, jx, plvl)
 739 ! Compute temperature from virtual potential temperature
 740   use storage_module
 741   implicit none
 742   integer :: ix, jx
 743   real :: plvl
 744   real, dimension(ix,jx) :: T, VPTMP, P, Q
 745
 746   real, parameter :: rovcp=0.28571
 747
 748   call get_storage(nint(plvl), 'VPTMP',  VPTMP, ix, jx)
 749   IF (nint(plvl) .LT. 200) THEN
 750     call get_storage(nint(plvl), 'PRESSURE',   P, ix, jx)
 751   ELSE
 752     p = plvl
 753   ENDIF
 754   call get_storage(nint(plvl), 'SPECHUMD',   Q, ix, jx)
 755
 756    t=vptmp * (p*1.e-5)**rovcp * (1./(1.+0.6078*Q))
 757
 758   call put_storage(nint(plvl), 'TT', t, ix, jx)
 759        if(nint(plvl).eq.1) then
 760   call put_storage(200100, 'PSFC', p, ix, jx)
 761        endif
 762
 763 end subroutine compute_t_vptmp
 764
 765
 766 subroutine compute_rh_spechumd(ix, jx)
 767 ! Compute relative humidity from specific humidity.
 768   use storage_module
 769   implicit none
 770   integer :: ix, jx
 771   real, dimension(ix,jx) :: T, P, RH, Q
 772
 773   real, parameter :: svp1=611.2
 774   real, parameter :: svp2=17.67
 775   real, parameter :: svp3=29.65
 776   real, parameter :: svpt0=273.15
 777   real, parameter :: eps = 0.622
 778
 779   call get_storage(200100, 'TT',        T, ix, jx)
 780   call get_storage(200100, 'PSFC',     P, ix, jx)
 781   call get_storage(200100, 'SPECHUMD', Q, ix, jx)
 782
 783   rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
 784
 785   call put_storage(200100, 'RH', rh, ix, jx)
 786
 787 end subroutine compute_rh_spechumd
 788
 789 subroutine compute_rh_spechumd_upa(ix, jx, plvl)
 790 ! Compute relative humidity from specific humidity.
 791   use storage_module
 792   implicit none
 793   integer :: ix, jx
 794   real :: plvl
 795   real, dimension(ix,jx) :: T, P, RH, Q
 796
 797   real, parameter :: svp1=611.2
 798   real, parameter :: svp2=17.67
 799   real, parameter :: svp3=29.65
 800   real, parameter :: svpt0=273.15
 801   real, parameter :: eps = 0.622
 802
 803   IF ( nint(plvl).LT. 200) THEN
 804     if (is_there(nint(plvl), 'PRESSURE')) then
 805       call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
 806     else
 807       return     ! if we don't have pressure on model levels, return
 808     endif
 809   ELSE
 810     P = plvl
 811   ENDIF
 812   call get_storage(nint(plvl), 'TT',        T, ix, jx)
 813   call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx)
 814   Q=MAX(1.E-10,Q)
 815
 816   rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
 817
 818   call put_storage(nint(plvl), 'RH', rh, ix, jx)
 819
 820 end subroutine compute_rh_spechumd_upa
 821
 822 subroutine compute_rh_vapp_upa(ix, jx, plvl)
 823 ! Compute relative humidity from vapor pressure.
 824 ! Thanks to Bob Hart of PSU ESSC -- 1999-05-27.
 825   use storage_module
 826   implicit none
 827   integer :: ix, jx
 828   real :: plvl
 829   real, dimension(ix,jx) :: P, ES
 830   real, pointer, dimension(:,:) :: T, E, RH
 831
 832   real, parameter :: svp1=611.2
 833   real, parameter :: svp2=17.67
 834   real, parameter :: svp3=29.65
 835   real, parameter :: svpt0=273.15
 836
 837   allocate(RH(ix,jx))
 838
 839   IF ( nint(plvl).LT. 200) THEN
 840     if (is_there(nint(plvl), 'PRESSURE')) then
 841       call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
 842     else
 843       return     ! if we don't have pressure on model levels, return
 844     endif
 845   ELSE
 846     P = plvl
 847   ENDIF
 848   call refr_storage(nint(plvl), 'TT',    T, ix, jx)
 849   call refr_storage(nint(plvl), 'VAPP', E, ix, jx)
 850
 851   ES=svp1*exp(svp2*(T-svpt0)/(T-svp3))
 852   rh=min(1.E2*(P-ES)*E/((P-E)*ES), 1.E2)
 853
 854   call refw_storage(nint(plvl), 'RH', rh, ix, jx)
 855
 856   nullify(T,E)
 857
 858 end subroutine compute_rh_vapp_upa
 859
 860 subroutine compute_rh_depr(ix, jx, plvl)
 861 ! Compute relative humidity from Dewpoint Depression
 862   use storage_module
 863   implicit none
 864   integer :: ix, jx
 865   real :: plvl
 866   real, dimension(ix,jx) :: t, depr, rh
 867
 868   real, parameter :: Xlv = 2.5e6
 869   real, parameter :: Rv = 461.5
 870
 871   integer :: i, j
 872
 873   call get_storage(nint(plvl), 'TT', T,  ix, jx)
 874   call get_storage(nint(plvl), 'DEPR', DEPR, ix, jx)
 875
 876   where(DEPR < 100.)
 877      rh = exp(Xlv/Rv*(1./T - 1./(T-depr))) * 1.E2
 878   elsewhere
 879      rh = 0.0
 880   endwhere
 881
 882   call put_storage(nint(plvl),'RH      ', rh, ix, jx)
 883
 884 end subroutine compute_rh_depr
 885
 886 subroutine compute_rh_dewpt(ix,jx)
 887 ! Compute relative humidity from Dewpoint
 888   use storage_module
 889   implicit none
 890   integer :: ix, jx
 891   real, dimension(ix,jx) :: t, dp, rh
 892
 893   real, parameter :: Xlv = 2.5e6
 894   real, parameter :: Rv = 461.5
 895
 896   call get_storage(200100, 'TT      ', T,  ix, jx)
 897   call get_storage(200100, 'DEWPT   ', DP, ix, jx)
 898
 899   rh = exp(Xlv/Rv*(1./T - 1./dp)) * 1.E2
 900
 901   call put_storage(200100,'RH      ', rh, ix, jx)
 902
 903 end subroutine compute_rh_dewpt
 904
 905 subroutine gfs_trop_maxw_pressures(ix,jx)
 906 ! These are duplicate pressure values from the GFS, for
 907 ! the level of max wind speed and for the trop level.
 908 ! The duplicates are saved with a different name, so that
 909 ! the metgrid program can horizontally interpolate them to
 910 ! the model domain with a nearest neighbor method.
 911   use storage_module
 912   implicit none
 913   integer :: ix, jx
 914   real, dimension(ix,jx) :: pmaxw, pmaxwnn, ptrop, ptropnn
 915
 916   if ( is_there(200100, 'PMAXW   ') ) then
 917      call get_storage(200100, 'PMAXW   ', pmaxw  , ix, jx)
 918      pmaxwnn = pmaxw
 919      call put_storage(200100, 'PMAXWNN ', pmaxwnn, ix, jx)
 920   end if
 921
 922   if ( is_there(200100, 'PTROP   ') ) then
 923      call get_storage(200100, 'PTROP   ', ptrop  , ix, jx)
 924      ptropnn = ptrop
 925      call put_storage(200100, 'PTROPNN ', ptropnn, ix, jx)
 926   end if
 927
 928 end subroutine gfs_trop_maxw_pressures
 929
 930 subroutine vntrp(plvl, maxlvl, k, name, ix, jx)
 931   use storage_module
 932   implicit none
 933   integer :: ix, jx, k, maxlvl
 934   real, dimension(maxlvl) :: plvl
 935   character(len=8) :: name
 936   real, dimension(ix,jx) :: a, b, c
 937   real :: frc
 938
 939   write(*,'("Interpolating to fill in ", A, " at level ", I8)') trim(name), nint(plvl(k))
 940
 941   call  get_storage(nint(plvl(k-1)), name, a, ix, jx)
 942   call  get_storage(nint(plvl(k+1)), name, c, ix, jx)
 943
 944   frc = (plvl(k) - plvl(k+1)) / ( plvl(k-1)-plvl(k+1))
 945
 946   b = (1.-frc)*c + frc*a
 947 !KWM  b = 0.5 * (a + c)
 948   call  put_storage(nint(plvl(k)), name, b, ix, jx)
 949
 950 end subroutine vntrp
 951
 952 subroutine fix_gfs_rh (ix, jx, plvl)
 953 ! This routine replaces GFS RH (wrt ice) with RH wrt liquid (which is what is assumed in real.exe).
 954   use storage_module
 955   implicit none
 956   integer :: ix, jx, i, j
 957   real :: plvl, eis, ews, r
 958   real, allocatable, dimension(:,:) :: rh, tt
 959
 960   allocate(rh(ix,jx))
 961   allocate(tt(ix,jx))
 962   call get_storage(nint(plvl), 'RH', rh, ix, jx)
 963   call get_storage(nint(plvl), 'TT', tt, ix, jx)
 964   do j = 1, jx
 965   do i = 1, ix
 966     if ( tt(i,j) .le. 273.15 ) then
 967       ! Murphy and Koop 2005 ice saturation vapor pressure.
 968       ! eis and ews in hPA, tt is in K
 969       eis = .01 * exp (9.550426 - (5723.265 / tt(i,j)) + (3.53068 * alog(tt(i,j))) &
 970          - (0.00728332 * tt(i,j)))
 971       ! Bolton 1980 liquid saturation vapor pressure. For water saturation, most
 972       ! formulae are very similar from 0 to -20, so we don't need a more exact formula.
 973
 974       ews = 6.112 * exp(17.67 * (tt(i,j)-273.15) / ((tt(i,j)-273.15)+243.5))
 975       if ( tt(i,j) .gt. 253.15 ) then
 976         ! A linear approximation to the GFS blending region ( -20 > T < 0 )
 977         r = ((273.15 - tt(i,j)) / 20.)
 978         r = (r * eis) + ((1-r)*ews)
 979       else
 980         r = eis
 981       endif
 982       rh(i,j) = rh(i,j) * (r / ews)
 983     endif
 984   enddo
 985   enddo
 986   call put_storage(nint(plvl), 'RH', rh, ix, jx)
 987   deallocate (rh)
 988   deallocate (tt)
 989 end subroutine fix_gfs_rh
 990
 991
 992 subroutine fix_ruc_soilm (ix, jx)
 993 ! This routine adds residual soil moisture if initialized fron RUC
 994   use storage_module
 995   implicit none
 996   integer :: ix, jx, i, j
 997   REAL , DIMENSION(100) :: lqmi
 998   real, allocatable, dimension(:,:) :: soilm000, soilm005, soilm020, &
 999                      soilm040, soilm160, soilm300,soilcat
1000   allocate(soilm000(ix,jx))
1001   allocate(soilm005(ix,jx))
1002   allocate(soilm020(ix,jx))
1003   allocate(soilm040(ix,jx))
1004   allocate(soilm160(ix,jx))
1005   allocate(soilm300(ix,jx))
1006   allocate(soilcat(ix,jx))
1007   call get_storage(200100, 'SM000ruc', soilm000, ix, jx)
1008   call get_storage(200100, 'SM005ruc', soilm005, ix, jx)
1009   call get_storage(200100, 'SM020ruc', soilm020, ix, jx)
1010   call get_storage(200100, 'SM040ruc', soilm040, ix, jx)
1011   call get_storage(200100, 'SM160ruc', soilm160, ix, jx)
1012   call get_storage(200100, 'SM300ruc', soilm300, ix, jx)
1013
1014   call get_storage(200100, 'SOILCAT', soilcat, ix, jx)
1015
1016       lqmi(1:16) = &
1017       (/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10,     &
1018         0.089, 0.095, 0.10,  0.070, 0.068, 0.078, 0.0,      &
1019         0.004, 0.065 /)
1020
1021   do j = 1, jx
1022   do i = 1, ix
1023
1024          SOILM000(i,j)=SOILM000(i,j) + lqmi(nint(soilcat(i,j)))
1025          SOILM005(i,j)=SOILM005(i,j) + lqmi(nint(soilcat(i,j)))
1026          SOILM020(i,j)=SOILM020(i,j) + lqmi(nint(soilcat(i,j)))
1027          SOILM040(i,j)=SOILM040(i,j) + lqmi(nint(soilcat(i,j)))
1028          SOILM160(i,j)=SOILM160(i,j) + lqmi(nint(soilcat(i,j)))
1029          SOILM300(i,j)=SOILM300(i,j) + lqmi(nint(soilcat(i,j)))
1030   enddo
1031   enddo
1032   call put_storage(200100, 'SOILM000', soilm000, ix, jx)
1033   call put_storage(200100, 'SOILM005', soilm005, ix, jx)
1034   call put_storage(200100, 'SOILM020', soilm020, ix, jx)
1035   call put_storage(200100, 'SOILM040', soilm040, ix, jx)
1036   call put_storage(200100, 'SOILM160', soilm160, ix, jx)
1037   call put_storage(200100, 'SOILM300', soilm300, ix, jx)
1038
1039  print *,'fix_ruc_soilm is done!'
1040
1041   deallocate(soilm000)
1042   deallocate(soilm005)
1043   deallocate(soilm020)
1044   deallocate(soilm040)
1045   deallocate(soilm160)
1046   deallocate(soilm300)
1047   deallocate(soilcat)
1048
1049 end subroutine fix_ruc_soilm
1050
1051