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 ! We have reached the end of file, so time to close it.
 192 !
 193      close(iunit)
 194      if (debug_level .gt. 100 ) call print_storage
 195 !
 196 ! By now the file has been read completely.  Now, see if we need to fill in
 197 ! missing fields:
 198 !
 199
 200 ! Retrieve the number of levels in storage:
 201 !
 202      call get_plvls(plvl, maxlvl, nlvl)
 203 !
 204 ! Fill the surface level (code 200100) from higher 200100s, as necessary
 205 !
 206         do k = 1, nlvl
 207            if ((plvl(k).gt.200100) .and. (plvl(k).lt.200200)) then
 208            ! We found a level between 200100 and 200200, now find the field
 209            ! corresponding to that level.
 210               MLOOP : do mm = 1, maxvar
 211                  if (is_there(nint(plvl(k)), namvar(mm))) then
 212                     INLOOP : do kk = 200101, nint(plvl(k))
 213                        if (is_there(kk, namvar(mm))) then
 214                           if ( debug_level .gt. 100 ) then
 215                             call mprintf(.true.,DEBUG, &
 216                "Copying %s at level %i to level 200100.",s1=namvar(mm),i1=kk)
 217                           end if
 218                           call get_dims(kk, namvar(mm))
 219                           allocate(scr2d(map%nx,map%ny))
 220                           call get_storage &
 221                                (kk, namvar(mm), scr2d, map%nx, map%ny)
 222                           call put_storage &
 223                                (200100,namvar(mm), scr2d,map%nx,map%ny)
 224                           deallocate(scr2d)
 225                           EXIT INLOOP
 226                        endif
 227                     enddo INLOOP
 228                  endif
 229               enddo MLOOP
 230            endif
 231         enddo
 232
 233 !
 234 ! If upper-air U is missing, see if we can interpolate from surrounding levels.
 235 ! This is a simple vertical interpolation, linear in pressure.
 236 ! Currently, this simply fills in one missing level between two present levels.
 237 !
 238
 239         do k = 2, nlvl-1, 1
 240            if (plvl(k-1) .lt. 200000.) then
 241               if ( (.not. is_there(nint(plvl(k)),'UU')) .and. &
 242                    ( is_there(nint(plvl(k-1)), 'UU')) .and.&
 243                    ( is_there(nint(plvl(k+1)), 'UU')) ) then
 244                  call get_dims(nint(plvl(k+1)), 'UU')
 245                  call vntrp(plvl, maxlvl, k, "UU      ", map%nx, map%ny)
 246               endif
 247            endif
 248         enddo
 249
 250 !
 251 ! If upper-air V is missing, see if we can interpolate from surrounding levels.
 252 ! This is a simple vertical interpolation, linear in pressure.
 253 ! Currently, this simply fills in one missing level between two present levels.
 254 !
 255
 256         do k = 2, nlvl-1, 1
 257            if (plvl(k-1) .lt. 200000.) then
 258               if ( (.not. is_there(nint(plvl(k)),'VV')) .and. &
 259                    ( is_there(nint(plvl(k-1)), 'VV')) .and.&
 260                    ( is_there(nint(plvl(k+1)), 'VV')) ) then
 261                  call get_dims(nint(plvl(k+1)), 'VV')
 262                  call vntrp(plvl, maxlvl, k, "VV      ", map%nx, map%ny)
 263               endif
 264            endif
 265         enddo
 266
 267 !
 268 ! If upper-air SPECHUMD is missing, see if we can compute SPECHUMD from QVAPOR:
 269 !--- Tanya's change for initializing WRF with RUC
 270
 271         do k = 1, nlvl
 272            if (plvl(k).lt.200000.) then
 273               if (.not. is_there(nint(plvl(k)), 'SPECHUMD').and. &
 274                    is_there(nint(plvl(k)), 'QV')) then
 275                  call get_dims(nint(plvl(k)), 'QV')
 276                  call compute_spechumd_qvapor(map%nx, map%ny, plvl(k))
 277               endif
 278            endif
 279         enddo
 280
 281 !--- Tanya's change for initializing WRF with RUC
 282 !   This allows for the ingestion for RUC isentropic data
 283 !
 284         do k = 1, nlvl
 285            if (plvl(k).lt.200000.) then
 286               if (.not. is_there(nint(plvl(k)), 'TT').and. &
 287                    is_there(nint(plvl(k)), 'VPTMP').and. &
 288                    is_there(nint(plvl(k)), 'SPECHUMD')) then
 289                  call get_dims(nint(plvl(k)), 'VPTMP')
 290                  call compute_t_vptmp(map%nx, map%ny, plvl(k))
 291               endif
 292            endif
 293         enddo
 294 !!!
 295 !
 296 ! If upper-air T is missing, see if we can interpolate from surrounding levels.
 297 ! This is a simple vertical interpolation, linear in pressure.
 298 ! Currently, this simply fills in one missing level between two present levels.
 299 !
 300
 301         do k = 2, nlvl-1, 1
 302            if (plvl(k-1) .lt. 200000.) then
 303               if ( (.not. is_there(nint(plvl(k)),'TT')) .and. &
 304                    ( is_there(nint(plvl(k-1)), 'TT')) .and.&
 305                    ( is_there(nint(plvl(k+1)), 'TT')) ) then
 306                  call get_dims(nint(plvl(k+1)), 'TT')
 307                  call vntrp(plvl, maxlvl, k, "TT      ", map%nx, map%ny)
 308               endif
 309            endif
 310         enddo
 311
 312 !
 313 ! Check to see if we need to fill HGT from GEOPT.
 314 !
 315         do k = 1, nlvl
 316            if (plvl(k).lt.200000.) then
 317               if (.not. is_there(nint(plvl(k)), 'HGT').and. &
 318                    is_there(nint(plvl(k)), 'GEOPT')) then
 319                  call get_dims(nint(plvl(k)), 'GEOPT')
 320                  allocate(scr2d(map%nx,map%ny))
 321                  call get_storage(nint(plvl(k)), 'GEOPT', scr2d, map%nx, map%ny)
 322                  scr2d = scr2d / 9.81
 323                  call put_storage(nint(plvl(k)), 'HGT',   scr2d, map%nx, map%ny)
 324                  deallocate(scr2d)
 325               endif
 326            endif
 327         enddo
 328
 329
 330 ! If upper-air RH is missing, see if we can compute RH from Specific Humidity:
 331
 332         do k = 1, nlvl
 333            if (plvl(k).lt.200000.) then
 334               if (.not. is_there(nint(plvl(k)), 'RH') .and. &
 335                    is_there(nint(plvl(k)), 'TT') .and. &
 336                    is_there(nint(plvl(k)), 'SPECHUMD')) then
 337                  call get_dims(nint(plvl(k)), 'TT')
 338                  call compute_rh_spechumd_upa(map%nx, map%ny, plvl(k))
 339               endif
 340            endif
 341         enddo
 342
 343 ! If upper-air RH is missing, see if we can compute RH from Vapor Pressure:
 344 !   (Thanks to Bob Hart of PSU ESSC -- 1999-05-27.)
 345
 346         do k = 1, nlvl
 347            if (plvl(k).lt.200000.) then
 348               if (.not. is_there(nint(plvl(k)),'RH').and. &
 349                    is_there(nint(plvl(k)), 'TT') .and. &
 350                    is_there(nint(plvl(k)),'VAPP')) then
 351                  call get_dims(nint(plvl(k)),'TT')
 352                  call compute_rh_vapp_upa(map%nx, map%ny, plvl(k))
 353               endif
 354            endif
 355         enddo
 356
 357 ! If upper-air RH is missing, see if we can compute RH from Dewpoint Depression:
 358
 359         do k = 1, nlvl
 360            if (plvl(k).lt.200000.) then
 361               if (.not. is_there(nint(plvl(k)),'RH').and. &
 362                    is_there(nint(plvl(k)), 'TT') .and. &
 363                    is_there(nint(plvl(k)),'DEPR')) then
 364                  call get_dims(nint(plvl(k)),'TT')
 365                  call compute_rh_depr(map%nx, map%ny, plvl(k))
 366               endif
 367            endif
 368         enddo
 369 !
 370 ! If upper-air RH is missing, see if we can interpolate from surrounding levels.
 371 ! This is a simple vertical interpolation, linear in pressure.
 372 ! Currently, this simply fills in one missing level between two present levels.
 373 ! May expand this in the future to fill in additional levels.  May also expand
 374 ! this in the future to vertically interpolate other variables.
 375 !
 376
 377         do k = 2, nlvl-1, 1
 378            if (plvl(k-1) .lt. 200000.) then
 379               if ( (.not. is_there(nint(plvl(k)),'RH')) .and. &
 380                    ( is_there(nint(plvl(k-1)), 'RH')) .and.&
 381                    ( is_there(nint(plvl(k+1)), 'RH')) ) then
 382                  call get_dims(nint(plvl(k+1)), 'RH')
 383                  call vntrp(plvl, maxlvl, k, "RH      ", map%nx, map%ny)
 384               endif
 385            endif
 386         enddo
 387
 388 ! Repair GFS and ECMWF pressure-level RH
 389         if (index(map%source,'NCEP GFS') .ne. 0 .or.  &
 390             index(map%source,'ECMWF') .ne. 0 ) then
 391           call mprintf(.true.,DEBUG, &
 392              "RRPR:   Adjusting GFS/ECMWF RH values ")
 393           do k = 1, nlvl
 394             if ( is_there(nint(plvl(k)),'RH') .and. &
 395                  is_there(nint(plvl(k)),'TT') ) then
 396               call fix_gfs_rh (map%nx, map%ny, plvl(k))
 397             endif
 398           enddo
 399         endif
 400 !
 401 ! Check to see if we need to fill RH above 300 mb:
 402 !
 403         if (is_there(30000, 'RH')) then
 404            call get_dims(30000, 'RH')
 405            allocate(scr2d(map%nx,map%ny))
 406
 407            do k = 1, nlvl
 408 !   Set missing RH to 5% between 300 and 70 hPa. Set RH to 0 above 70 hPa.
 409 !   The stratospheric RH will be adjusted further in real.
 410               if (plvl(k).le.7000.) then
 411                 scr2d = 0.
 412               else if (plvl(k).lt.30000.) then
 413                 scr2d = 5.
 414               endif
 415               if (plvl(k).lt.30000. .and. plvl(k) .gt. 10. ) then
 416               ! levels higher than .1 mb are special - do not fill
 417                  if (.not. is_there(nint(plvl(k)), 'RH')) then
 418                     call put_storage(nint(plvl(k)),'RH',scr2d,map%nx,map%ny)
 419                     call mprintf(.true.,DEBUG, &
 420                  "RRPR:   RH missing at %i hPa, inserting synthetic RH ",i1=nint(plvl(k)/100.))
 421                  endif
 422               endif
 423            enddo
 424            deallocate(scr2d)
 425         endif
 426 !
 427 ! If surface RH is missing, see if we can compute RH from Specific Humidity
 428 ! or Dewpoint or Dewpoint depression:
 429 !
 430         if (.not. is_there (200100, 'RH')) then
 431            if (is_there(200100, 'TT').and. &
 432                 is_there(200100, 'PSFC'    )   .and. &
 433                 is_there(200100, 'SPECHUMD')) then
 434               call get_dims(200100, 'TT')
 435               call compute_rh_spechumd(map%nx, map%ny)
 436               call mprintf(.true.,DEBUG, &
 437                 "RRPR:   SURFACE RH is computed")
 438            elseif (is_there(200100, 'TT'       ).and. &
 439                 is_there(200100, 'DEWPT')) then
 440               call get_dims(200100, 'TT')
 441               call compute_rh_dewpt(map%nx, map%ny)
 442            elseif (is_there(200100, 'TT').and. &
 443                 is_there(200100, 'DEPR')) then
 444               call get_dims(200100, 'TT')
 445               call compute_rh_depr(map%nx, map%ny, 200100.)
 446            endif
 447         endif
 448
 449 !
 450 ! If surface SNOW is missing, see if we can compute SNOW from SNOWRUC
 451 ! (From Wei Wang, 2007 June 21, modified 12/28/2007)
 452 !
 453         if (.not. is_there(200100, 'SNOW') .and. &
 454              is_there(200100, 'SNOWRUC')) then
 455            call get_dims(200100, 'SNOWRUC')
 456            allocate(scr2d(map%nx,map%ny))
 457            call get_storage(200100, 'SNOWRUC', scr2d, map%nx, map%ny)
 458            scr2d = scr2d * 1000.
 459            call put_storage(200100, 'SNOW',   scr2d, map%nx, map%ny)
 460            deallocate(scr2d)
 461         endif
 462 !
 463 ! Check to see if we need to fill SOILHGT from SOILGEO.
 464 ! (From Wei Wang, 2007 June 21)
 465 !
 466         if (.not. is_there(200100, 'SOILHGT') .and. &
 467              is_there(200100, 'SOILGEO')) then
 468            call get_dims(200100, 'SOILGEO')
 469            allocate(scr2d(map%nx,map%ny))
 470            call get_storage(200100, 'SOILGEO', scr2d, map%nx, map%ny)
 471            scr2d = scr2d / 9.81
 472            call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny)
 473            deallocate(scr2d)
 474         endif
 475
 476 ! For hybrid-level input, soilgeo is in level 1 (e.g. ERA40)
 477         if (.not. is_there(200100, 'SOILHGT') .and. &
 478              is_there(1, 'SOILGEO')) then
 479            call get_dims(1, 'SOILGEO')
 480            allocate(scr2d(map%nx,map%ny))
 481            call get_storage(1, 'SOILGEO', scr2d, map%nx, map%ny)
 482            scr2d = scr2d / 9.81
 483            call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny)
 484            deallocate(scr2d)
 485         endif
 486
 487 ! For ECMWF hybrid-level input, may need to move psfc from level 1 to 2001.
 488         if ( index(map%source,'ECMWF') .ne. 0) then
 489         if (.not. is_there(200100, 'PSFC') .and. &
 490              is_there(1, 'PSFCH')) then
 491            call get_dims(1, 'PSFCH')
 492            allocate(scr2d(map%nx,map%ny))
 493            call get_storage(1, 'PSFCH', scr2d, map%nx, map%ny)
 494            call put_storage(200100, 'PSFC', scr2d, map%nx, map%ny)
 495            deallocate(scr2d)
 496         endif
 497         endif
 498
 499 ! ECMWF snow depth in meters of water equivalent (Table 128). Convert to kg/m2
 500 !
 501         if (is_there(200100, 'SNOW_EC')) then
 502            call get_dims(200100, 'SNOW_EC')
 503            allocate(scr2d(map%nx,map%ny))
 504            call get_storage(200100, 'SNOW_EC', scr2d, map%nx, map%ny)
 505            scr2d = scr2d * 1000.
 506            call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
 507            deallocate(scr2d)
 508         endif
 509
 510 ! Convert the ECMWF LANDSEA mask from a fraction to a flag
 511
 512         if ( index(map%source,'ECMWF') .ne. 0) then
 513         if (is_there(200100, 'LANDSEA')) then
 514            call get_dims(200100, 'LANDSEA')
 515            call make_zero_or_one(map%nx, map%ny, 'LANDSEA')
 516         endif
 517         endif
 518
 519 ! NCEP GFS weasd is one-half of the NAM value. Increase it for use in WRF.
 520 ! The GFS-based reanalyses values should be OK as is.
 521         if ((index(map%source,'NCEP GFS') .ne. 0 .or. &
 522             index(map%source,'NCEP GEFS') .ne. 0) .and. &
 523             is_there(200100, 'SNOW')) then
 524            call mprintf(.true.,DEBUG, &
 525               "RRPR:   Recomputing SNOW for NCEP GFS")
 526            call get_dims(200100, 'SNOW')
 527            allocate(scr2d(map%nx,map%ny))
 528            call get_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
 529            scr2d = scr2d * 2.
 530            call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
 531            deallocate(scr2d)
 532         endif
 533
 534 ! compute physical snow depth (SNOWH) for various models
 535 ! As of March 2011, this is done here instead of real because we have model
 536 ! source information.
 537         if (is_there(200100, 'SNOW') .and. .not. is_there(200100, 'SNOWH')) then
 538           call get_dims(200100, 'SNOW')
 539           allocate(scr2d(map%nx,map%ny))
 540           call get_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
 541           call mprintf(.true.,DEBUG, &
 542              "RRPR:   Computing SNOWH from SNOW")
 543           if ( index(map%source,'NCEP ') .ne. 0) then
 544             scr2d = scr2d * 0.005          ! Assume 200:1 ratio as used at NCEP and in NOAH
 545           else if (index(map%source,'ECMWF') .ne. 0) then
 546             if (is_there(200100, 'SNOW_DEN')) then        ! If we have snow density, use it to compute snowh
 547               call get_dims(200100, 'SNOW_DEN')
 548               allocate(tmp2d(map%nx,map%ny))
 549               call get_storage(200100, 'SNOW_DEN', tmp2d, map%nx, map%ny)
 550               scr2d = scr2d / tmp2d
 551               deallocate(tmp2d)
 552             else
 553               scr2d = scr2d * 0.004     ! otherwise, assume a density of 250 mm/m (i.e. 250:1 ratio).
 554             endif
 555           else                       ! Other models
 556             scr2d = scr2d * 0.005    ! Use real's default method (200:1)
 557           endif
 558           call put_storage(200100, 'SNOWH', scr2d, map%nx, map%ny)
 559           deallocate(scr2d)
 560         endif
 561
 562 ! As of March 2011, SEAICE can be a flag or a fraction. It will be converted
 563 ! to the appropriate values in real depending on whether or not the polar mods are used.
 564
 565 !! If we've got a SEAICE field, make sure that it is all Zeros and Ones:
 566
 567 !       if (is_there(200100, 'SEAICE')) then
 568 !          call get_dims(200100, 'SEAICE')
 569 !          call make_zero_or_one(map%nx, map%ny, 'SEAICE')
 570 !       endif
 571
 572 ! If we've got an ICEMASK field, re-flag it for output to met_em and real:
 573 !     Field  | GRIB In  |  Out
 574 !    -------------------------
 575 !    water   |    0     |  0
 576 !    land    |   -1     |  1
 577 !    ice     |    1     |  0
 578
 579         if (is_there(200100, 'ICEMASK')) then
 580            call get_dims(200100, 'ICEMASK')
 581            call re_flag_ice_mask(map%nx, map%ny)
 582         endif
 583
 584 ! If we have an ICEFRAC field, convert from % to fraction
 585         if (is_there(200100, 'ICEFRAC')) then
 586            call get_dims(200100, 'ICEFRAC')
 587            allocate(scr2d(map%nx,map%ny))
 588            call get_storage(200100, 'ICEFRAC', scr2d, map%nx, map%ny)
 589            scr2d = scr2d / 100.
 590            call put_storage(200100, 'ICEFRAC', scr2d, map%nx, map%ny)
 591            deallocate(scr2d)
 592         endif
 593
 594
 595         call mprintf(.true.,INFORM, &
 596            "RRPR: hdate = %s ",s1=hdate)
 597         call output(hdate, nlvl, maxlvl, plvl, interval, 2, out_format, prefix, debug_level)
 598         call clear_storage
 599         exit FILELOOP
 600      enddo FILELOOP
 601    enddo TIMELOOP
 602 end subroutine rrpr
 603
 604 subroutine make_zero_or_one(ix, jx, infield)
 605 ! Make sure the input field (SEAICE or LANDSEA) is zero or one.
 606   use storage_module
 607   implicit none
 608   integer :: ix, jx
 609   real, dimension(ix,jx) :: seaice
 610   character(len=*) :: infield
 611
 612   call get_storage(200100, infield, seaice, ix, jx)
 613   where(seaice > 0.5)
 614      seaice = 1.0
 615   elsewhere
 616      seaice = 0.0
 617   end where
 618   call put_storage(200100, infield, seaice, ix, jx)
 619 end subroutine make_zero_or_one
 620
 621 subroutine re_flag_ice_mask(ix, jx)
 622 !
 623 ! Change land points from -1 to 1
 624 ! Change ice  points from  1 to 0
 625 ! Water       points stay    at 0
 626 !
 627   use storage_module
 628   implicit none
 629   integer :: ix, jx
 630   real, dimension(ix,jx) :: iceflag
 631
 632   call get_storage(200100, 'ICEMASK',iceflag, ix, jx)
 633   where(iceflag > 0.5)     ! Ice points, set to water value
 634      iceflag = 0.0
 635   end where
 636   where(iceflag < -0.5)    ! Land points
 637      iceflag = 1.0
 638   end where
 639   call put_storage(200100, 'ICEMASK',iceflag, ix, jx)
 640 end subroutine re_flag_ice_mask
 641
 642 subroutine compute_spechumd_qvapor(ix, jx, plvl)
 643 ! Compute specific humidity from water vapor mixing ratio.
 644   use storage_module
 645   implicit none
 646   integer :: ix, jx
 647   real :: plvl
 648   real, dimension(ix,jx) :: QVAPOR, SPECHUMD
 649
 650   call get_storage(nint(plvl), 'QV', QVAPOR, ix, jx)
 651
 652   SPECHUMD = QVAPOR/(1.+QVAPOR)
 653
 654   call put_storage(nint(plvl), 'SPECHUMD', spechumd, ix, jx)
 655  if(nint(plvl).eq.1) then
 656   call put_storage(200100,'SPECHUMD', spechumd, ix, jx)
 657  endif
 658
 659 end subroutine compute_spechumd_qvapor
 660
 661 subroutine compute_t_vptmp(ix, jx, plvl)
 662 ! Compute temperature from virtual potential temperature
 663   use storage_module
 664   implicit none
 665   integer :: ix, jx
 666   real :: plvl
 667   real, dimension(ix,jx) :: T, VPTMP, P, Q
 668
 669   real, parameter :: rovcp=0.28571
 670
 671   call get_storage(nint(plvl), 'VPTMP',  VPTMP, ix, jx)
 672   IF (nint(plvl) .LT. 200) THEN
 673     call get_storage(nint(plvl), 'PRESSURE',   P, ix, jx)
 674   ELSE
 675     p = plvl
 676   ENDIF
 677   call get_storage(nint(plvl), 'SPECHUMD',   Q, ix, jx)
 678
 679    t=vptmp * (p*1.e-5)**rovcp * (1./(1.+0.6078*Q))
 680
 681   call put_storage(nint(plvl), 'TT', t, ix, jx)
 682        if(nint(plvl).eq.1) then
 683   call put_storage(200100, 'PSFC', p, ix, jx)
 684        endif
 685
 686 end subroutine compute_t_vptmp
 687
 688
 689 subroutine compute_rh_spechumd(ix, jx)
 690 ! Compute relative humidity from specific humidity.
 691   use storage_module
 692   implicit none
 693   integer :: ix, jx
 694   real, dimension(ix,jx) :: T, P, RH, Q
 695
 696   real, parameter :: svp1=611.2
 697   real, parameter :: svp2=17.67
 698   real, parameter :: svp3=29.65
 699   real, parameter :: svpt0=273.15
 700   real, parameter :: eps = 0.622
 701
 702   call get_storage(200100, 'TT',        T, ix, jx)
 703   call get_storage(200100, 'PSFC',     P, ix, jx)
 704   call get_storage(200100, 'SPECHUMD', Q, ix, jx)
 705
 706   rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
 707
 708   call put_storage(200100, 'RH', rh, ix, jx)
 709
 710 end subroutine compute_rh_spechumd
 711
 712 subroutine compute_rh_spechumd_upa(ix, jx, plvl)
 713 ! Compute relative humidity from specific humidity.
 714   use storage_module
 715   implicit none
 716   integer :: ix, jx
 717   real :: plvl
 718   real, dimension(ix,jx) :: T, P, RH, Q
 719
 720   real, parameter :: svp1=611.2
 721   real, parameter :: svp2=17.67
 722   real, parameter :: svp3=29.65
 723   real, parameter :: svpt0=273.15
 724   real, parameter :: eps = 0.622
 725
 726   IF ( nint(plvl).LT. 200) THEN
 727     if (is_there(nint(plvl), 'PRESSURE')) then
 728       call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
 729     else
 730       return     ! if we don't have pressure on model levels, return
 731     endif
 732   ELSE
 733     P = plvl
 734   ENDIF
 735   call get_storage(nint(plvl), 'TT',        T, ix, jx)
 736   call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx)
 737   Q=MAX(1.E-10,Q)
 738
 739   rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
 740
 741   call put_storage(nint(plvl), 'RH', rh, ix, jx)
 742
 743 end subroutine compute_rh_spechumd_upa
 744
 745 subroutine compute_rh_vapp_upa(ix, jx, plvl)
 746 ! Compute relative humidity from vapor pressure.
 747 ! Thanks to Bob Hart of PSU ESSC -- 1999-05-27.
 748   use storage_module
 749   implicit none
 750   integer :: ix, jx
 751   real :: plvl
 752   real, dimension(ix,jx) :: P, ES
 753   real, pointer, dimension(:,:) :: T, E, RH
 754
 755   real, parameter :: svp1=611.2
 756   real, parameter :: svp2=17.67
 757   real, parameter :: svp3=29.65
 758   real, parameter :: svpt0=273.15
 759
 760   allocate(RH(ix,jx))
 761
 762   IF ( nint(plvl).LT. 200) THEN
 763     if (is_there(nint(plvl), 'PRESSURE')) then
 764       call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
 765     else
 766       return     ! if we don't have pressure on model levels, return
 767     endif
 768   ELSE
 769     P = plvl
 770   ENDIF
 771   call refr_storage(nint(plvl), 'TT',    T, ix, jx)
 772   call refr_storage(nint(plvl), 'VAPP', E, ix, jx)
 773
 774   ES=svp1*exp(svp2*(T-svpt0)/(T-svp3))
 775   rh=min(1.E2*(P-ES)*E/((P-E)*ES), 1.E2)
 776
 777   call refw_storage(nint(plvl), 'RH', rh, ix, jx)
 778
 779   nullify(T,E)
 780
 781 end subroutine compute_rh_vapp_upa
 782
 783 subroutine compute_rh_depr(ix, jx, plvl)
 784 ! Compute relative humidity from Dewpoint Depression
 785   use storage_module
 786   implicit none
 787   integer :: ix, jx
 788   real :: plvl
 789   real, dimension(ix,jx) :: t, depr, rh
 790
 791   real, parameter :: Xlv = 2.5e6
 792   real, parameter :: Rv = 461.5
 793
 794   integer :: i, j
 795
 796   call get_storage(nint(plvl), 'TT', T,  ix, jx)
 797   call get_storage(nint(plvl), 'DEPR', DEPR, ix, jx)
 798
 799   where(DEPR < 100.)
 800      rh = exp(Xlv/Rv*(1./T - 1./(T-depr))) * 1.E2
 801   elsewhere
 802      rh = 0.0
 803   endwhere
 804
 805   call put_storage(nint(plvl),'RH      ', rh, ix, jx)
 806
 807 end subroutine compute_rh_depr
 808
 809 subroutine compute_rh_dewpt(ix,jx)
 810 ! Compute relative humidity from Dewpoint
 811   use storage_module
 812   implicit none
 813   integer :: ix, jx
 814   real, dimension(ix,jx) :: t, dp, rh
 815
 816   real, parameter :: Xlv = 2.5e6
 817   real, parameter :: Rv = 461.5
 818
 819   call get_storage(200100, 'TT      ', T,  ix, jx)
 820   call get_storage(200100, 'DEWPT   ', DP, ix, jx)
 821
 822   rh = exp(Xlv/Rv*(1./T - 1./dp)) * 1.E2
 823
 824   call put_storage(200100,'RH      ', rh, ix, jx)
 825
 826 end subroutine compute_rh_dewpt
 827
 828 subroutine vntrp(plvl, maxlvl, k, name, ix, jx)
 829   use storage_module
 830   implicit none
 831   integer :: ix, jx, k, maxlvl
 832   real, dimension(maxlvl) :: plvl
 833   character(len=8) :: name
 834   real, dimension(ix,jx) :: a, b, c
 835   real :: frc
 836
 837   write(*,'("Interpolating to fill in ", A, " at level ", I8)') trim(name), nint(plvl(k))
 838
 839   call  get_storage(nint(plvl(k-1)), name, a, ix, jx)
 840   call  get_storage(nint(plvl(k+1)), name, c, ix, jx)
 841
 842   frc = (plvl(k) - plvl(k+1)) / ( plvl(k-1)-plvl(k+1))
 843
 844   b = (1.-frc)*a + frc*c
 845 !KWM  b = 0.5 * (a + c)
 846   call  put_storage(nint(plvl(k)), name, b, ix, jx)
 847
 848 end subroutine vntrp
 849
 850 subroutine fix_gfs_rh (ix, jx, plvl)
 851 ! This routine replaces GFS RH (wrt ice) with RH wrt liquid (which is what is assumed in real.exe).
 852   use storage_module
 853   implicit none
 854   integer :: ix, jx, i, j
 855   real :: plvl, eis, ews, r
 856   real, allocatable, dimension(:,:) :: rh, tt
 857
 858   allocate(rh(ix,jx))
 859   allocate(tt(ix,jx))
 860   call get_storage(nint(plvl), 'RH', rh, ix, jx)
 861   call get_storage(nint(plvl), 'TT', tt, ix, jx)
 862   do j = 1, jx
 863   do i = 1, ix
 864     if ( tt(i,j) .le. 273.15 ) then
 865       ! Murphy and Koop 2005 ice saturation vapor pressure.
 866       ! eis and ews in hPA, tt is in K
 867       eis = .01 * exp (9.550426 - (5723.265 / tt(i,j)) + (3.53068 * alog(tt(i,j))) &
 868          - (0.00728332 * tt(i,j)))
 869       ! Bolton 1980 liquid saturation vapor pressure. For water saturation, most
 870       ! formulae are very similar from 0 to -20, so we don't need a more exact formula.
 871
 872       ews = 6.112 * exp(17.67 * (tt(i,j)-273.15) / ((tt(i,j)-273.15)+243.5))
 873       if ( tt(i,j) .gt. 253.15 ) then
 874         ! A linear approximation to the GFS blending region ( -20 > T < 0 )
 875         r = ((273.15 - tt(i,j)) / 20.)
 876         r = (r * eis) + ((1-r)*ews)
 877       else
 878         r = eis
 879       endif
 880       rh(i,j) = rh(i,j) * (r / ews)
 881     endif
 882   enddo
 883   enddo
 884   call put_storage(nint(plvl), 'RH', rh, ix, jx)
 885   deallocate (rh)
 886   deallocate (tt)
 887 end subroutine fix_gfs_rh