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