WPS/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
  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   write(0,*) 'Begin rrpr'
  68
  69   if (mod(interval,3600) == 0) then
  70      datelen = 13
  71   else if (mod(interval, 60) == 0) then
  72      datelen = 16
  73   else
  74      datelen = 19
  75   endif
  76
  77   if ( debug_level .gt. 100 ) then
  78     call mprintf(.true.,DEBUG,"Begin rrpr")
  79     call mprintf(.true.,DEBUG,"nfiles = %i , ntimes = %i )",i1=nfiles,i2=ntimes)
  80     do n = 1, nfiles
  81       call mprintf(.true.,DEBUG,"filedates(%i) = %s",i1=n,s1=filedates(n))
  82     enddo
  83   endif
  84
  85 ! Compute the ending time:
  86
  87   call geth_newdate(hend, hstart, interval*ntimes)
  88
  89   call clear_storage
  90
  91 ! We want to do something for each of the requested times:
  92   TIMELOOP : do ntime = 1, ntimes
  93      idts = (ntime-1) * interval
  94      call geth_newdate(hdate, hstart, idts)
  95      call mprintf(.true.,DEBUG, &
  96      "RRPR: hstart = %s , hdate = %s , idts = %i",s1=hstart,s2=hdate,i1=idts)
  97
  98 ! Loop over the output file dates, and do stuff if the file date matches
  99 ! the requested time we are working on now.
 100
 101      FILELOOP : do n = 1, nfiles
 102        if ( debug_level .gt. 100 ) then
 103          call mprintf(.true.,DEBUG, &
 104             "hstart = %s , hend = %s",s1=hstart,s2=hend)
 105          call mprintf(.true.,DEBUG, &
 106             "filedates(n) = %s",s1=filedates(n))
 107          call mprintf(.true.,DEBUG, &
 108             "filedates(n) = %s",s1=filedates(n)(1:datelen))
 109        end if
 110        if (filedates(n)(1:datelen).ne.hdate(1:datelen)) cycle FILELOOP
 111        if (debug_level .gt. 50 ) then
 112          call mprintf(.true.,INFORM, &
 113             "RRPR Processing : %s",s1=filedates(n)(1:datelen))
 114        endif
 115        open(iunit, file=trim(get_path(prefix))//'PFILE:'//filedates(n)(1:datelen), &
 116           form='unformatted',status='old')
 117
 118 ! Read the file:
 119
 120      rdloop: do
 121         read (iunit, iostat=ierr) ifv
 122         if (ierr.ne.0) exit rdloop
 123         if ( ifv .eq. 5) then     ! WPS
 124           read (iunit) hdate_output, xfcst, map%source, field, units, Desc, &
 125                level, map%nx, map%ny, map%igrid
 126           hdate = hdate_output(1:19)
 127           select case (map%igrid)
 128           case (0, 4)
 129              read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%r_earth
 130           case (3)
 131            read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 132                 map%truelat1, map%truelat2, map%r_earth
 133           case (5)
 134              read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 135                 map%truelat1, map%r_earth
 136           case (1)
 137            read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
 138                 map%truelat1, map%r_earth
 139           case default
 140              call mprintf(.true.,ERROR, &
 141                 "Unrecognized map%%igrid: %i in RRPR 1",i1=map%igrid)
 142           end select
 143           read (iunit) map%grid_wind
 144
 145         else if ( ifv .eq. 4 ) then          ! SI
 146           read (iunit) hdate_output, xfcst, map%source, field, units, desc, level, &
 147                 map%nx, map%ny, map%igrid
 148           hdate = hdate_output(1:19)
 149           select case (map%igrid)
 150           case (0, 4)
 151              read(iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx
 152           case (3)
 153              read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 154                 map%lov, map%truelat1, map%truelat2
 155           case (5)
 156              read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
 157                 map%lov, map%truelat1
 158           case default
 159              call mprintf(.true.,ERROR, &
 160                 "Unrecognized map%%igrid: %i in RRPR 2",i1=map%igrid)
 161           end select
 162
 163         else if ( ifv .eq. 3 ) then          ! MM5
 164           read(iunit) hdate_output, xfcst, field, units, desc, level,&
 165                 map%nx, map%ny, map%igrid
 166           hdate = hdate_output(1:19)
 167           select case (map%igrid)
 168           case (3)      ! lamcon
 169             read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 170                     map%truelat1, map%truelat2
 171            case (5)      ! Polar Stereographic
 172               read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
 173                    map%truelat1
 174            case (0, 4)      ! lat/lon
 175               read (iunit) map%lat1, map%lon1, map%dy, map%dx
 176            case (1)      ! Mercator
 177               read (iunit) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
 178            case default
 179              call mprintf(.true.,ERROR, &
 180                 "Unrecognized map%%igrid: %i in RRPR 3",i1=map%igrid)
 181            end select
 182         else
 183            call mprintf(.true.,ERROR, &
 184               "unknown out_format, ifv = %i",i1=ifv)
 185         endif
 186
 187         allocate(ptr2d(map%nx,map%ny))
 188         read (iunit) ptr2d
 189         call refw_storage(nint(level), field, ptr2d, map%nx, map%ny)
 190         nullify (ptr2d)
 191      enddo rdloop
 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 upper-air RH is missing, see if we can compute RH from Specific Humidity:
 333
 334         do k = 1, nlvl
 335            if (plvl(k).lt.200000.) then
 336               if (.not. is_there(nint(plvl(k)), 'RH').and. &
 337                    is_there(nint(plvl(k)), 'SPECHUMD')) then
 338                  call get_dims(nint(plvl(k)), 'TT')
 339                  call compute_rh_spechumd_upa(map%nx, map%ny, plvl(k))
 340               endif
 341            endif
 342         enddo
 343
 344 ! If upper-air RH is missing, see if we can compute RH from Vapor Pressure:
 345 !   (Thanks to Bob Hart of PSU ESSC -- 1999-05-27.)
 346
 347         do k = 1, nlvl
 348            if (plvl(k).lt.200000.) then
 349               if (.not. is_there(nint(plvl(k)),'RH').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)),'DEPR')) then
 363                  call get_dims(nint(plvl(k)),'TT')
 364                  call compute_rh_depr(map%nx, map%ny, plvl(k))
 365               endif
 366            endif
 367         enddo
 368 !
 369 ! If upper-air RH is missing, see if we can interpolate from surrounding levels.
 370 ! This is a simple vertical interpolation, linear in pressure.
 371 ! Currently, this simply fills in one missing level between two present levels.
 372 ! May expand this in the future to fill in additional levels.  May also expand
 373 ! this in the future to vertically interpolate other variables.
 374 !
 375
 376         do k = 2, nlvl-1, 1
 377            if (plvl(k-1) .lt. 200000.) then
 378               if ( (.not. is_there(nint(plvl(k)),'RH')) .and. &
 379                    ( is_there(nint(plvl(k-1)), 'RH')) .and.&
 380                    ( is_there(nint(plvl(k+1)), 'RH')) ) then
 381                  call get_dims(nint(plvl(k+1)), 'RH')
 382                  call vntrp(plvl, maxlvl, k, "RH      ", map%nx, map%ny)
 383               endif
 384            endif
 385         enddo
 386
 387 ! Repair GFS and ECMWF pressure-level RH
 388         if (index(map%source,'NCEP GFS') .ne. 0 .or.  &
 389             index(map%source,'ECMWF') .ne. 0 ) then
 390           call mprintf(.true.,DEBUG, &
 391              "RRPR:   Adjusting GFS/ECMWF RH values ")
 392           do k = 1, nlvl
 393             if ( is_there(nint(plvl(k)),'RH') .and. &
 394                  is_there(nint(plvl(k)),'TT') ) then
 395               call fix_gfs_rh (map%nx, map%ny, plvl(k))
 396             endif
 397           enddo
 398         endif
 399 !
 400 ! Check to see if we need to fill RH above 300 mb:
 401 !
 402         if (is_there(30000, 'RH')) then
 403            call get_dims(30000, 'RH')
 404            allocate(scr2d(map%nx,map%ny))
 405
 406            do k = 1, nlvl
 407 !   Set missing RH to 5% between 300 and 70 hPa. Set RH to 0 above 70 hPa.
 408 !   The stratospheric RH will be adjusted further in real.
 409               if (plvl(k).le.7000.) then
 410                 scr2d = 0.
 411               else if (plvl(k).lt.30000.) then
 412                 scr2d = 5.
 413               endif
 414               if (plvl(k).lt.30000. .and. plvl(k) .gt. 10. ) then
 415               ! levels higher than .1 mb are special - do not fill
 416                  if (.not. is_there(nint(plvl(k)), 'RH')) then
 417                     call put_storage(nint(plvl(k)),'RH',scr2d,map%nx,map%ny)
 418                     call mprintf(.true.,DEBUG, &
 419                  "RRPR:   RH missing at %i hPa, inserting synthetic RH ",i1=nint(plvl(k)/100.))
 420                  endif
 421               endif
 422            enddo
 423            deallocate(scr2d)
 424         endif
 425 !
 426 ! If surface RH is missing, see if we can compute RH from Specific Humidity
 427 ! or Dewpoint or Dewpoint depression:
 428 !
 429         if (.not. is_there (200100, 'RH')) then
 430            if (is_there(200100, 'TT').and. &
 431                 is_there(200100, 'PSFC'    )   .and. &
 432                 is_there(200100, 'SPECHUMD')) then
 433               call get_dims(200100, 'TT')
 434               call compute_rh_spechumd(map%nx, map%ny)
 435               call mprintf(.true.,DEBUG, &
 436                 "RRPR:   SURFACE RH is computed")
 437            elseif (is_there(200100, 'TT'       ).and. &
 438                 is_there(200100, 'DEWPT')) then
 439               call get_dims(200100, 'TT')
 440               call compute_rh_dewpt(map%nx, map%ny)
 441            elseif (is_there(200100, 'TT').and. &
 442                 is_there(200100, 'DEPR')) then
 443               call get_dims(200100, 'TT')
 444               call compute_rh_depr(map%nx, map%ny, 200100.)
 445            endif
 446         endif
 447
 448 !
 449 ! If surface SNOW is missing, see if we can compute SNOW from SNOWRUC
 450 ! (From Wei Wang, 2007 June 21, modified 12/28/2007)
 451 !
 452         if (.not. is_there(200100, 'SNOW') .and. &
 453              is_there(200100, 'SNOWRUC')) then
 454            call get_dims(200100, 'SNOWRUC')
 455            allocate(scr2d(map%nx,map%ny))
 456            call get_storage(200100, 'SNOWRUC', scr2d, map%nx, map%ny)
 457            scr2d = scr2d * 1000.
 458            call put_storage(200100, 'SNOW',   scr2d, map%nx, map%ny)
 459            deallocate(scr2d)
 460         endif
 461 !
 462 ! Check to see if we need to fill SOILHGT from SOILGEO.
 463 ! (From Wei Wang, 2007 June 21)
 464 !
 465         if (.not. is_there(200100, 'SOILHGT') .and. &
 466              is_there(200100, 'SOILGEO')) then
 467            call get_dims(200100, 'SOILGEO')
 468            allocate(scr2d(map%nx,map%ny))
 469            call get_storage(200100, 'SOILGEO', scr2d, map%nx, map%ny)
 470            scr2d = scr2d / 9.81
 471            call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny)
 472            deallocate(scr2d)
 473         endif
 474
 475 ! For hybrid-level input, soilgeo is in level 1 (e.g. ERA40)
 476         if (.not. is_there(200100, 'SOILHGT') .and. &
 477              is_there(1, 'SOILGEO')) then
 478            call get_dims(1, 'SOILGEO')
 479            allocate(scr2d(map%nx,map%ny))
 480            call get_storage(1, 'SOILGEO', scr2d, map%nx, map%ny)
 481            scr2d = scr2d / 9.81
 482            call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny)
 483            deallocate(scr2d)
 484         endif
 485
 486 ! For ECMWF hybrid-level input, may need to move psfc from level 1 to 2001.
 487         if ( index(map%source,'ECMWF') .ne. 0) then
 488         if (.not. is_there(200100, 'PSFC') .and. &
 489              is_there(1, 'PSFCH')) then
 490            call get_dims(1, 'PSFCH')
 491            allocate(scr2d(map%nx,map%ny))
 492            call get_storage(1, 'PSFCH', scr2d, map%nx, map%ny)
 493            call put_storage(200100, 'PSFC', scr2d, map%nx, map%ny)
 494            deallocate(scr2d)
 495         endif
 496         endif
 497
 498 ! ECMWF snow depth is in meters (Table 128). Convert to kg/m2
 499 !
 500         if (is_there(200100, 'SNOW_EC')) then
 501            call get_dims(200100, 'SNOW_EC')
 502            allocate(scr2d(map%nx,map%ny))
 503            call get_storage(200100, 'SNOW_EC', scr2d, map%nx, map%ny)
 504            scr2d = scr2d * 1000.
 505            call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
 506            deallocate(scr2d)
 507         endif
 508
 509 ! If we've got a SEAICE field, make sure that it is all Zeros and Ones:
 510
 511         if (is_there(200100, 'SEAICE')) then
 512            call get_dims(200100, 'SEAICE')
 513            call make_zero_or_one(map%nx, map%ny)
 514         endif
 515
 516 ! If we've got an ICEMASK field, re-flag it for output to met_em and real:
 517 !     Field  | GRIB In  |  Out
 518 !    -------------------------
 519 !    water   |    0     |  0
 520 !    land    |   -1     |  1
 521 !    ice     |    1     |  0
 522
 523         if (is_there(200100, 'ICEMASK')) then
 524            call get_dims(200100, 'ICEMASK')
 525            call re_flag_ice_mask(map%nx, map%ny)
 526         endif
 527
 528 ! If we have an ICEFRAC field, convert from % to fraction
 529         if (is_there(200100, 'ICEFRAC')) then
 530            call get_dims(200100, 'ICEFRAC')
 531            allocate(scr2d(map%nx,map%ny))
 532            call get_storage(200100, 'ICEFRAC', scr2d, map%nx, map%ny)
 533            scr2d = scr2d / 100.
 534            call put_storage(200100, 'ICEFRAC', scr2d, map%nx, map%ny)
 535            deallocate(scr2d)
 536         endif
 537
 538
 539         call mprintf(.true.,INFORM, &
 540            "RRPR: hdate = %s ",s1=hdate)
 541         call output(hdate, nlvl, maxlvl, plvl, interval, 2, out_format, prefix, debug_level)
 542         call clear_storage
 543         exit FILELOOP
 544      enddo FILELOOP
 545    enddo TIMELOOP
 546 end subroutine rrpr
 547
 548 subroutine make_zero_or_one(ix, jx)
 549 ! Make sure the SEAICE field is zero or one.
 550   use storage_module
 551   implicit none
 552   integer :: ix, jx
 553   real, dimension(ix,jx) :: seaice
 554
 555   call get_storage(200100, 'SEAICE',seaice, ix, jx)
 556   where(seaice > 0.5)
 557      seaice = 1.0
 558   elsewhere
 559      seaice = 0.0
 560   end where
 561   call put_storage(200100, 'SEAICE',seaice, ix, jx)
 562 end subroutine make_zero_or_one
 563
 564 subroutine re_flag_ice_mask(ix, jx)
 565 !
 566 ! Change land points from -1 to 1
 567 ! Change ice  points from  1 to 0
 568 ! Water       points stay    at 0
 569 !
 570   use storage_module
 571   implicit none
 572   integer :: ix, jx
 573   real, dimension(ix,jx) :: iceflag
 574
 575   call get_storage(200100, 'ICEMASK',iceflag, ix, jx)
 576   where(iceflag > 0.5)     ! Ice points, set to water value
 577      iceflag = 0.0
 578   end where
 579   where(iceflag < -0.5)    ! Land points
 580      iceflag = 1.0
 581   end where
 582   call put_storage(200100, 'ICEMASK',iceflag, ix, jx)
 583 end subroutine re_flag_ice_mask
 584
 585 subroutine compute_spechumd_qvapor(ix, jx, plvl)
 586 ! Compute specific humidity from water vapor mixing ratio.
 587   use storage_module
 588   implicit none
 589   integer :: ix, jx
 590   real :: plvl
 591   real, dimension(ix,jx) :: QVAPOR, SPECHUMD
 592
 593   call get_storage(nint(plvl), 'QV', QVAPOR, ix, jx)
 594
 595   SPECHUMD = QVAPOR/(1.+QVAPOR)
 596
 597   call put_storage(nint(plvl), 'SPECHUMD', spechumd, ix, jx)
 598  if(nint(plvl).eq.1) then
 599   call put_storage(200100,'SPECHUMD', spechumd, ix, jx)
 600  endif
 601
 602 end subroutine compute_spechumd_qvapor
 603
 604 subroutine compute_t_vptmp(ix, jx, plvl)
 605 ! Compute temperature from virtual potential temperature
 606   use storage_module
 607   implicit none
 608   integer :: ix, jx
 609   real :: plvl
 610   real, dimension(ix,jx) :: T, VPTMP, P, Q
 611
 612   real, parameter :: rovcp=0.28571
 613
 614   call get_storage(nint(plvl), 'VPTMP',  VPTMP, ix, jx)
 615   IF (nint(plvl) .LT. 200) THEN
 616     call get_storage(nint(plvl), 'PRESSURE',   P, ix, jx)
 617   ELSE
 618     p = plvl
 619   ENDIF
 620   call get_storage(nint(plvl), 'SPECHUMD',   Q, ix, jx)
 621
 622    t=vptmp * (p*1.e-5)**rovcp * (1./(1.+0.6078*Q))
 623
 624   call put_storage(nint(plvl), 'TT', t, ix, jx)
 625        if(nint(plvl).eq.1) then
 626   call put_storage(200100, 'PSFC', p, ix, jx)
 627        endif
 628
 629 end subroutine compute_t_vptmp
 630
 631
 632 subroutine compute_rh_spechumd(ix, jx)
 633 ! Compute relative humidity from specific humidity.
 634   use storage_module
 635   implicit none
 636   integer :: ix, jx
 637   real, dimension(ix,jx) :: T, P, RH, Q
 638
 639   real, parameter :: svp1=611.2
 640   real, parameter :: svp2=17.67
 641   real, parameter :: svp3=29.65
 642   real, parameter :: svpt0=273.15
 643   real, parameter :: eps = 0.622
 644
 645   call get_storage(200100, 'TT',        T, ix, jx)
 646   call get_storage(200100, 'PSFC',     P, ix, jx)
 647   call get_storage(200100, 'SPECHUMD', Q, ix, jx)
 648
 649   rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
 650
 651   call put_storage(200100, 'RH', rh, ix, jx)
 652
 653 end subroutine compute_rh_spechumd
 654
 655 subroutine compute_rh_spechumd_upa(ix, jx, plvl)
 656 ! Compute relative humidity from specific humidity.
 657   use storage_module
 658   implicit none
 659   integer :: ix, jx
 660   real :: plvl
 661   real, dimension(ix,jx) :: T, P, RH, Q
 662
 663   real, parameter :: svp1=611.2
 664   real, parameter :: svp2=17.67
 665   real, parameter :: svp3=29.65
 666   real, parameter :: svpt0=273.15
 667   real, parameter :: eps = 0.622
 668
 669   IF ( nint(plvl).LT. 200) THEN
 670     call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
 671   ELSE
 672     P = plvl
 673   ENDIF
 674   call get_storage(nint(plvl), 'TT',        T, ix, jx)
 675   call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx)
 676   Q=MAX(1.E-10,Q)
 677
 678   rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
 679
 680   call put_storage(nint(plvl), 'RH', rh, ix, jx)
 681
 682 end subroutine compute_rh_spechumd_upa
 683
 684 subroutine compute_rh_vapp_upa(ix, jx, plvl)
 685 ! Compute relative humidity from vapor pressure.
 686 ! Thanks to Bob Hart of PSU ESSC -- 1999-05-27.
 687   use storage_module
 688   implicit none
 689   integer :: ix, jx
 690   real :: plvl
 691   real, dimension(ix,jx) :: P, ES
 692   real, pointer, dimension(:,:) :: T, E, RH
 693
 694   real, parameter :: svp1=611.2
 695   real, parameter :: svp2=17.67
 696   real, parameter :: svp3=29.65
 697   real, parameter :: svpt0=273.15
 698
 699   allocate(RH(ix,jx))
 700
 701   P = plvl
 702   call refr_storage(nint(plvl), 'TT',    T, ix, jx)
 703   call refr_storage(nint(plvl), 'VAPP', E, ix, jx)
 704
 705   ES=svp1*exp(svp2*(T-svpt0)/(T-svp3))
 706   rh=min(1.E2*(P-ES)*E/((P-E)*ES), 1.E2)
 707
 708   call refw_storage(nint(plvl), 'RH', rh, ix, jx)
 709
 710   nullify(T,E)
 711
 712 end subroutine compute_rh_vapp_upa
 713
 714 subroutine compute_rh_depr(ix, jx, plvl)
 715 ! Compute relative humidity from Dewpoint Depression
 716   use storage_module
 717   implicit none
 718   integer :: ix, jx
 719   real :: plvl
 720   real, dimension(ix,jx) :: t, depr, rh
 721
 722   real, parameter :: Xlv = 2.5e6
 723   real, parameter :: Rv = 461.5
 724
 725   integer :: i, j
 726
 727   call get_storage(nint(plvl), 'TT', T,  ix, jx)
 728   call get_storage(nint(plvl), 'DEPR', DEPR, ix, jx)
 729
 730   where(DEPR < 100.)
 731      rh = exp(Xlv/Rv*(1./T - 1./(T-depr))) * 1.E2
 732   elsewhere
 733      rh = 0.0
 734   endwhere
 735
 736   call put_storage(nint(plvl),'RH      ', rh, ix, jx)
 737
 738 end subroutine compute_rh_depr
 739
 740 subroutine compute_rh_dewpt(ix,jx)
 741 ! Compute relative humidity from Dewpoint
 742   use storage_module
 743   implicit none
 744   integer :: ix, jx
 745   real, dimension(ix,jx) :: t, dp, rh
 746
 747   real, parameter :: Xlv = 2.5e6
 748   real, parameter :: Rv = 461.5
 749
 750   call get_storage(200100, 'TT      ', T,  ix, jx)
 751   call get_storage(200100, 'DEWPT   ', DP, ix, jx)
 752
 753   rh = exp(Xlv/Rv*(1./T - 1./dp)) * 1.E2
 754
 755   call put_storage(200100,'RH      ', rh, ix, jx)
 756
 757 end subroutine compute_rh_dewpt
 758
 759 subroutine vntrp(plvl, maxlvl, k, name, ix, jx)
 760   use storage_module
 761   implicit none
 762   integer :: ix, jx, k, maxlvl
 763   real, dimension(maxlvl) :: plvl
 764   character(len=8) :: name
 765   real, dimension(ix,jx) :: a, b, c
 766   real :: frc
 767
 768   write(*,'("Interpolating to fill in ", A, " at level ", I8)') trim(name), nint(plvl(k))
 769
 770   call  get_storage(nint(plvl(k-1)), name, a, ix, jx)
 771   call  get_storage(nint(plvl(k+1)), name, c, ix, jx)
 772
 773   frc = (plvl(k) - plvl(k+1)) / ( plvl(k-1)-plvl(k+1))
 774
 775   b = (1.-frc)*a + frc*c
 776 !KWM  b = 0.5 * (a + c)
 777   call  put_storage(nint(plvl(k)), name, b, ix, jx)
 778
 779 end subroutine vntrp
 780
 781 subroutine fix_gfs_rh (ix, jx, plvl)
 782 ! This routine replaces GFS RH (wrt ice) with RH wrt liquid (which is what is assumed in real.exe).
 783   use storage_module
 784   implicit none
 785   integer :: ix, jx, i, j
 786   real :: plvl, eis, ews, r
 787   real, allocatable, dimension(:,:) :: rh, tt
 788
 789   allocate(rh(ix,jx))
 790   allocate(tt(ix,jx))
 791   call get_storage(nint(plvl), 'RH', rh, ix, jx)
 792   call get_storage(nint(plvl), 'TT', tt, ix, jx)
 793   do j = 1, jx
 794   do i = 1, ix
 795     if ( tt(i,j) .le. 273.15 ) then
 796       ! Murphy and Koop 2005 ice saturation vapor pressure.
 797       ! eis and ews in hPA, tt is in K
 798       eis = .01 * exp (9.550426 - (5723.265 / tt(i,j)) + (3.53068 * alog(tt(i,j))) &
 799          - (0.00728332 * tt(i,j)))
 800       ! Bolton 1980 liquid saturation vapor pressure. For water saturation, most
 801       ! formulae are very similar from 0 to -20, so we don't need a more exact formula.
 802
 803       ews = 6.112 * exp(17.67 * (tt(i,j)-273.15) / ((tt(i,j)-273.15)+243.5))
 804       if ( tt(i,j) .gt. 253.15 ) then
 805         ! A linear approximation to the GFS blending region ( -20 > T < 0 )
 806         r = ((273.15 - tt(i,j)) / 20.)
 807         r = (r * eis) + ((1-r)*ews)
 808       else
 809         r = eis
 810       endif
 811       rh(i,j) = rh(i,j) * (r / ews)
 812     endif
 813   enddo
 814   enddo
 815   call put_storage(nint(plvl), 'RH', rh, ix, jx)
 816   deallocate (rh)
 817   deallocate (tt)
 818 end subroutine fix_gfs_rh