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