subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_format, prefix)
! !
! In case you are wondering, RRPR stands for "Read, ReProcess, and wRite" !
! !
!*****************************************************************************!
! !
! Recent changes: !
! !
! 2004-10-29: !
! - Sync rrpr.F WRFSI v2.0.1 with MM5 v3.5 !
! Added DATELEN: length of date strings to use for !
! our output file names. !
! Added MM5 interpolation from surrounding levels !
! if upper-air U or V are missing !
! Added test to see if we've got a SEAICE field, !
! make sure that it is all Zeros and Ones: !
! !
! 2002-05-16: !
! - Handle the Mercator projection. !
! This change also required changes to output.F, rd_grib.F, !
! datint.F, gribcode.F !
! !
! 2002-02-13: !
! - Added vertical interpolation in pressure in case of missing !
! U, V, T (the check for RH was already there) !
! !
! 2001-02-14: !
! - Allow file names to have date stamps out to minutes or !
! seconds, if the user requests a time interval (in seconds) !
! that is not evenly divisible into hours or minutes. !
! INTERVAL is checked for divisibility into 3600 (for hours) !
! or 60 (for minutes). The local variable DATELEN is set !
! to be the number of characters to use in our character !
! dates. Valid values for DATELEN are 13 (for hours), !
! 16 (for minutes), and 19 (for seconds). !
! !
! This change also requires changes to pregrid_grib.F, !
! output.F, datint.F, file_delete.F !
! !
! - Do processing not just if the requested date matches one we !
! want, but if the requested date falls between the startdate !
! and the enddate. !
! !
!*****************************************************************************!
use filelist
use gridinfo
use storage_module
use table
use module_debug
use misc_definitions_module
use stringutil
implicit none
!------------------------------------------------------------------------------
! Arguments:
! HSTART: Starting date of times to process
character (LEN=19) :: hstart
! NTIMES: Number of time periods to process
integer :: ntimes
! INTERVAL: Time inteval (seconds) of time periods to process.
integer :: interval
! NLVL: The number of levels in the stored data.
integer :: nlvl
! MAXLVL: The parameterized maximum number of levels to allow.
integer :: maxlvl
! PLVL: Array of pressure levels (Pa) in the dataset
real , dimension(maxlvl) :: plvl
! DEBUG_LEVEL: Integer level of debug printing (from namelist)
integer :: debug_level
!------------------------------------------------------------------------------
character (LEN=25) :: units
character (LEN=46) :: Desc
real, allocatable, dimension(:,:) :: scr2d
real, pointer, dimension(:,:) :: ptr2d
integer :: k, kk, mm, n, ierr, ifv
integer :: iunit=13
character(LEN=19) :: hdate, hend
character(LEN=24) :: hdate_output
character(LEN=3) :: out_format
character(LEN=MAX_FILENAME_LEN) :: prefix
real :: xfcst, level
character(LEN=9) :: field
integer :: ntime, idts
! DATELEN: length of date strings to use for our output file names.
integer :: datelen
! Decide the length of date strings to use for output file names.
! DATELEN is 13 for hours, 16 for minutes, and 19 for seconds.
write(0,*) 'Begin rrpr'
if (mod(interval,3600) == 0) then
datelen = 13
else if (mod(interval, 60) == 0) then
datelen = 16
else
datelen = 19
endif
if ( debug_level .gt. 100 ) then
call mprintf(.true.,DEBUG,"Begin rrpr")
call mprintf(.true.,DEBUG,"nfiles = %i , ntimes = %i )",i1=nfiles,i2=ntimes)
do n = 1, nfiles
call mprintf(.true.,DEBUG,"filedates(%i) = %s",i1=n,s1=filedates(n))
enddo
endif
! Compute the ending time:
call geth_newdate(hend, hstart, interval*ntimes)
call clear_storage
! We want to do something for each of the requested times:
TIMELOOP : do ntime = 1, ntimes
idts = (ntime-1) * interval
call geth_newdate(hdate, hstart, idts)
call mprintf(.true.,DEBUG, &
"RRPR: hstart = %s , hdate = %s , idts = %i",s1=hstart,s2=hdate,i1=idts)
! Loop over the output file dates, and do stuff if the file date matches
! the requested time we are working on now.
FILELOOP : do n = 1, nfiles
if ( debug_level .gt. 100 ) then
call mprintf(.true.,DEBUG, &
"hstart = %s , hend = %s",s1=hstart,s2=hend)
call mprintf(.true.,DEBUG, &
"filedates(n) = %s",s1=filedates(n))
call mprintf(.true.,DEBUG, &
"filedates(n) = %s",s1=filedates(n)(1:datelen))
end if
if (filedates(n)(1:datelen).ne.hdate(1:datelen)) cycle FILELOOP
if (debug_level .gt. 50 ) then
call mprintf(.true.,INFORM, &
"RRPR Processing : %s",s1=filedates(n)(1:datelen))
endif
open(iunit, file=trim(get_path(prefix))//'PFILE:'//filedates(n)(1:datelen), &
form='unformatted',status='old')
! Read the file:
rdloop: do
read (iunit, iostat=ierr) ifv
if (ierr.ne.0) exit rdloop
if ( ifv .eq. 5) then ! WPS
read (iunit) hdate_output, xfcst, map%source, field, units, Desc, &
level, map%nx, map%ny, map%igrid
hdate = hdate_output(1:19)
select case (map%igrid)
case (0, 4)
read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%r_earth
case (3)
read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, &
map%truelat1, map%truelat2, map%r_earth
case (5)
read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, &
map%truelat1, map%r_earth
case (1)
read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
map%truelat1, map%r_earth
case default
call mprintf(.true.,ERROR, &
"Unrecognized map%%igrid: %i in RRPR 1",i1=map%igrid)
end select
read (iunit) map%grid_wind
else if ( ifv .eq. 4 ) then ! SI
read (iunit) hdate_output, xfcst, map%source, field, units, desc, level, &
map%nx, map%ny, map%igrid
hdate = hdate_output(1:19)
select case (map%igrid)
case (0, 4)
read(iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx
case (3)
read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
map%lov, map%truelat1, map%truelat2
case (5)
read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
map%lov, map%truelat1
case default
call mprintf(.true.,ERROR, &
"Unrecognized map%%igrid: %i in RRPR 2",i1=map%igrid)
end select
else if ( ifv .eq. 3 ) then ! MM5
read(iunit) hdate_output, xfcst, field, units, desc, level,&
map%nx, map%ny, map%igrid
hdate = hdate_output(1:19)
select case (map%igrid)
case (3) ! lamcon
read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
map%truelat1, map%truelat2
case (5) ! Polar Stereographic
read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
map%truelat1
case (0, 4) ! lat/lon
read (iunit) map%lat1, map%lon1, map%dy, map%dx
case (1) ! Mercator
read (iunit) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
case default
call mprintf(.true.,ERROR, &
"Unrecognized map%%igrid: %i in RRPR 3",i1=map%igrid)
end select
else
call mprintf(.true.,ERROR, &
"unknown out_format, ifv = %i",i1=ifv)
endif
allocate(ptr2d(map%nx,map%ny))
read (iunit) ptr2d
call refw_storage(nint(level), field, ptr2d, map%nx, map%ny)
nullify (ptr2d)
enddo rdloop
!
! We have reached the end of file, so time to close it.
!
close(iunit)
if (debug_level .gt. 100 ) call print_storage
!
! By now the file has been read completely. Now, see if we need to fill in
! missing fields:
!
! Retrieve the number of levels in storage:
!
call get_plvls(plvl, maxlvl, nlvl)
!
! Fill the surface level (code 200100) from higher 200100s, as necessary
!
do k = 1, nlvl
if ((plvl(k).gt.200100) .and. (plvl(k).lt.200200)) then
! We found a level between 200100 and 200200, now find the field
! corresponding to that level.
MLOOP : do mm = 1, maxvar
if (is_there(nint(plvl(k)), namvar(mm))) then
INLOOP : do kk = 200101, nint(plvl(k))
if (is_there(kk, namvar(mm))) then
if ( debug_level .gt. 100 ) then
call mprintf(.true.,DEBUG, &
"Copying %s at level %i to level 200100.",s1=namvar(mm),i1=kk)
end if
call get_dims(kk, namvar(mm))
allocate(scr2d(map%nx,map%ny))
call get_storage &
(kk, namvar(mm), scr2d, map%nx, map%ny)
call put_storage &
(200100,namvar(mm), scr2d,map%nx,map%ny)
deallocate(scr2d)
EXIT INLOOP
endif
enddo INLOOP
endif
enddo MLOOP
endif
enddo
!
! If upper-air U is missing, see if we can interpolate from surrounding levels.
! This is a simple vertical interpolation, linear in pressure.
! Currently, this simply fills in one missing level between two present levels.
!
do k = 2, nlvl-1, 1
if (plvl(k-1) .lt. 200000.) then
if ( (.not. is_there(nint(plvl(k)),'UU')) .and. &
( is_there(nint(plvl(k-1)), 'UU')) .and.&
( is_there(nint(plvl(k+1)), 'UU')) ) then
call get_dims(nint(plvl(k+1)), 'UU')
call vntrp(plvl, maxlvl, k, "UU ", map%nx, map%ny)
endif
endif
enddo
!
! If upper-air V is missing, see if we can interpolate from surrounding levels.
! This is a simple vertical interpolation, linear in pressure.
! Currently, this simply fills in one missing level between two present levels.
!
do k = 2, nlvl-1, 1
if (plvl(k-1) .lt. 200000.) then
if ( (.not. is_there(nint(plvl(k)),'VV')) .and. &
( is_there(nint(plvl(k-1)), 'VV')) .and.&
( is_there(nint(plvl(k+1)), 'VV')) ) then
call get_dims(nint(plvl(k+1)), 'VV')
call vntrp(plvl, maxlvl, k, "VV ", map%nx, map%ny)
endif
endif
enddo
!
! If upper-air SPECHUMD is missing, see if we can compute SPECHUMD from QVAPOR:
!--- Tanya's change for initializing WRF with RUC
do k = 1, nlvl
if (plvl(k).lt.200000.) then
if (.not. is_there(nint(plvl(k)), 'SPECHUMD').and. &
is_there(nint(plvl(k)), 'QV')) then
call get_dims(nint(plvl(k)), 'QV')
call compute_spechumd_qvapor(map%nx, map%ny, plvl(k))
endif
endif
enddo
!--- Tanya's change for initializing WRF with RUC
! This allows for the ingestion for RUC isentropic data
!
do k = 1, nlvl
if (plvl(k).lt.200000.) then
if (.not. is_there(nint(plvl(k)), 'TT').and. &
is_there(nint(plvl(k)), 'VPTMP').and. &
is_there(nint(plvl(k)), 'SPECHUMD')) then
call get_dims(nint(plvl(k)), 'VPTMP')
call compute_t_vptmp(map%nx, map%ny, plvl(k))
endif
endif
enddo
!!!
!
! If upper-air T is missing, see if we can interpolate from surrounding levels.
! This is a simple vertical interpolation, linear in pressure.
! Currently, this simply fills in one missing level between two present levels.
!
do k = 2, nlvl-1, 1
if (plvl(k-1) .lt. 200000.) then
if ( (.not. is_there(nint(plvl(k)),'TT')) .and. &
( is_there(nint(plvl(k-1)), 'TT')) .and.&
( is_there(nint(plvl(k+1)), 'TT')) ) then
call get_dims(nint(plvl(k+1)), 'TT')
call vntrp(plvl, maxlvl, k, "TT ", map%nx, map%ny)
endif
endif
enddo
!
! Check to see if we need to fill HGT from GEOPT.
!
do k = 1, nlvl
if (plvl(k).lt.200000.) then
if (.not. is_there(nint(plvl(k)), 'HGT').and. &
is_there(nint(plvl(k)), 'GEOPT')) then
call get_dims(nint(plvl(k)), 'GEOPT')
allocate(scr2d(map%nx,map%ny))
call get_storage(nint(plvl(k)), 'GEOPT', scr2d, map%nx, map%ny)
scr2d = scr2d / 9.81
call put_storage(nint(plvl(k)), 'HGT', scr2d, map%nx, map%ny)
deallocate(scr2d)
endif
endif
enddo
!!! GFS cloud fix
do k = 1, nlvl
if (plvl(k).lt.200000.) then
if (.not. is_there(nint(10000.),'SNMR').and..not.is_there(nint(plvl(k)),'CLWMR').and. &
is_there(nint(10000.), 'CLWMR')) then
call get_dims(nint(10000.), 'CWLMR')
allocate(scr2d(map%nx,map%ny))
call get_storage(nint(10000.), 'CLWMR', scr2d, map%nx, map%ny)
scr2d = scr2d / 10.
print*, 'put fake CLWMR that is 10% of the 100 hPa level'
write(0,*) 'put fake CLWMR that is 10% of the 100 hPa level'
call mprintf(.true.,DEBUG,"Fake 10% write of 100 hPa CLWMR")
call put_storage(nint(plvl(k)), 'CLWMR', scr2d, map%nx, map%ny)
deallocate(scr2d)
endif
endif
enddo
! If upper-air RH is missing, see if we can compute RH from Specific Humidity:
do k = 1, nlvl
if (plvl(k).lt.200000.) then
if (.not. is_there(nint(plvl(k)), 'RH').and. &
is_there(nint(plvl(k)), 'SPECHUMD')) then
call get_dims(nint(plvl(k)), 'TT')
call compute_rh_spechumd_upa(map%nx, map%ny, plvl(k))
endif
endif
enddo
! If upper-air RH is missing, see if we can compute RH from Vapor Pressure:
! (Thanks to Bob Hart of PSU ESSC -- 1999-05-27.)
do k = 1, nlvl
if (plvl(k).lt.200000.) then
if (.not. is_there(nint(plvl(k)),'RH').and. &
is_there(nint(plvl(k)),'VAPP')) then
call get_dims(nint(plvl(k)),'TT')
call compute_rh_vapp_upa(map%nx, map%ny, plvl(k))
endif
endif
enddo
! If upper-air RH is missing, see if we can compute RH from Dewpoint Depression:
do k = 1, nlvl
if (plvl(k).lt.200000.) then
if (.not. is_there(nint(plvl(k)),'RH').and. &
is_there(nint(plvl(k)),'DEPR')) then
call get_dims(nint(plvl(k)),'TT')
call compute_rh_depr(map%nx, map%ny, plvl(k))
endif
endif
enddo
!
! If upper-air RH is missing, see if we can interpolate from surrounding levels.
! This is a simple vertical interpolation, linear in pressure.
! Currently, this simply fills in one missing level between two present levels.
! May expand this in the future to fill in additional levels. May also expand
! this in the future to vertically interpolate other variables.
!
do k = 2, nlvl-1, 1
if (plvl(k-1) .lt. 200000.) then
if ( (.not. is_there(nint(plvl(k)),'RH')) .and. &
( is_there(nint(plvl(k-1)), 'RH')) .and.&
( is_there(nint(plvl(k+1)), 'RH')) ) then
call get_dims(nint(plvl(k+1)), 'RH')
call vntrp(plvl, maxlvl, k, "RH ", map%nx, map%ny)
endif
endif
enddo
! Repair GFS RH
if (index(map%source,'NCEP GFS') .ne. 0) then
call mprintf(.true.,DEBUG, &
"RRPR: Adjusting GFS RH values ")
do k = 1, nlvl
if ( is_there(nint(plvl(k)),'RH') .and. &
is_there(nint(plvl(k)),'TT') ) then
call fix_gfs_rh (map%nx, map%ny, plvl(k))
endif
enddo
endif
!
! Check to see if we need to fill RH above 300 mb:
!
if (is_there(30000, 'RH')) then
call get_dims(30000, 'RH')
allocate(scr2d(map%nx,map%ny))
do k = 1, nlvl
! Set missing RH to 10% between 300 and 70 hPa. Use P/1000 above 70 hPa.
! The stratospheric RH will be adjusted further in real.
if (plvl(k).le.7000.) then
scr2d = plvl(k) / 1000.
else if (plvl(k).lt.30000.) then
scr2d = 10.
endif
if (plvl(k).lt.30000.) then
if (.not. is_there(nint(plvl(k)), 'RH')) then
call put_storage(nint(plvl(k)),'RH',scr2d,map%nx,map%ny)
endif
endif
enddo
deallocate(scr2d)
endif
!
! If surface RH is missing, see if we can compute RH from Specific Humidity
! or Dewpoint or Dewpoint depression:
!
if (.not. is_there (200100, 'RH')) then
if (is_there(200100, 'TT').and. &
is_there(200100, 'PSFC' ) .and. &
is_there(200100, 'SPECHUMD')) then
call get_dims(200100, 'TT')
call compute_rh_spechumd(map%nx, map%ny)
call mprintf(.true.,DEBUG, &
"RRPR: SURFACE RH is computed")
elseif (is_there(200100, 'TT' ).and. &
is_there(200100, 'DEWPT')) then
call get_dims(200100, 'TT')
call compute_rh_dewpt(map%nx, map%ny)
elseif (is_there(200100, 'TT').and. &
is_there(200100, 'DEPR')) then
call get_dims(200100, 'TT')
call compute_rh_depr(map%nx, map%ny, 200100.)
endif
endif
!
! If surface SNOW is missing, see if we can compute SNOW and SNOWH from SNOWRUC
! (From Wei Wang, 2007 June 21, modified 10/20/2007)
!
if (.not. is_there(200100, 'SNOW') .and. &
is_there(200100, 'SNOWRUC')) then
call get_dims(200100, 'SNOWRUC')
allocate(scr2d(map%nx,map%ny))
call get_storage(200100, 'SNOWRUC', scr2d, map%nx, map%ny)
scr2d = scr2d / 20.
call put_storage(200100, 'SNOWH', scr2d, map%nx, map%ny)
scr2d = scr2d / 6.
call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
deallocate(scr2d)
endif
!
! Check to see if we need to fill SOILHGT from SOILGEO.
! (From Wei Wang, 2007 June 21)
!
if (.not. is_there(200100, 'SOILHGT') .and. &
is_there(200100, 'SOILGEO')) then
call get_dims(200100, 'SOILGEO')
allocate(scr2d(map%nx,map%ny))
call get_storage(200100, 'SOILGEO', scr2d, map%nx, map%ny)
scr2d = scr2d / 9.81
call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny)
deallocate(scr2d)
endif
! If we've got a SEAICE field, make sure that it is all Zeros and Ones:
if (is_there(200100, 'SEAICE')) then
call get_dims(200100, 'SEAICE')
call make_zero_or_one(map%nx, map%ny)
endif
call mprintf(.true.,INFORM, &
"RRPR: hdate = %s ",s1=hdate)
call output(hdate, nlvl, maxlvl, plvl, interval, 2, out_format, prefix, debug_level)
call clear_storage
exit FILELOOP
enddo FILELOOP
enddo TIMELOOP
end subroutine rrpr
subroutine make_zero_or_one(ix, jx)
! Make sure the SEAICE field is zero or one.
use storage_module
implicit none
integer :: ix, jx
real, dimension(ix,jx) :: seaice
call get_storage(200100, 'SEAICE',seaice, ix, jx)
where(seaice > 0.5)
seaice = 1.0
elsewhere
seaice = 0.0
end where
call put_storage(200100, 'SEAICE',seaice, ix, jx)
end subroutine make_zero_or_one
subroutine compute_spechumd_qvapor(ix, jx, plvl)
! Compute specific humidity from water vapor mixing ratio.
use storage_module
implicit none
integer :: ix, jx
real :: plvl
real :: lat1, lon1, dx, dy
real, dimension(ix,jx) :: QVAPOR, SPECHUMD
real startlat, startlon, deltalat, deltalon
call get_storage(nint(plvl), 'QV', QVAPOR, ix, jx)
SPECHUMD = QVAPOR/(1.+QVAPOR)
call put_storage(nint(plvl), 'SPECHUMD', spechumd, ix, jx)
if(nint(plvl).eq.1) then
call put_storage(200100,'SPECHUMD', spechumd, ix, jx)
endif
end subroutine compute_spechumd_qvapor
subroutine compute_t_vptmp(ix, jx, plvl)
! Compute relative humidity from specific humidity.
use storage_module
implicit none
integer :: ix, jx
real :: plvl
real :: lat1, lon1, dx, dy
real, dimension(ix,jx) :: T, VPTMP, P, Q
real, parameter :: rovcp=0.28571
real startlat, startlon, deltalat, deltalon
call get_storage(nint(plvl), 'VPTMP', VPTMP, ix, jx)
IF (nint(plvl) .LT. 200) THEN
call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
ELSE
p = plvl
ENDIF
call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx)
t=vptmp * (p*1.e-5)**rovcp * (1./(1.+0.6078*Q))
call put_storage(nint(plvl), 'TT', t, ix, jx)
if(nint(plvl).eq.1) then
call put_storage(200100, 'PSFC', p, ix, jx)
endif
end subroutine compute_t_vptmp
subroutine compute_rh_spechumd(ix, jx)
! Compute relative humidity from specific humidity.
use storage_module
implicit none
integer :: ix, jx
real :: lat1, lon1, dx, dy
real, dimension(ix,jx) :: T, P, RH, Q
real, parameter :: svp1=611.2
real, parameter :: svp2=17.67
real, parameter :: svp3=29.65
real, parameter :: svpt0=273.15
real, parameter :: eps = 0.622
real startlat, startlon, deltalat, deltalon
call get_storage(200100, 'TT', T, ix, jx)
call get_storage(200100, 'PSFC', P, ix, jx)
call get_storage(200100, 'SPECHUMD', Q, ix, jx)
rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
call put_storage(200100, 'RH', rh, ix, jx)
end subroutine compute_rh_spechumd
subroutine compute_rh_spechumd_upa(ix, jx, plvl)
! Compute relative humidity from specific humidity.
use storage_module
implicit none
integer :: ix, jx
real :: plvl
real :: lat1, lon1, dx, dy
real, dimension(ix,jx) :: T, P, RH, Q
real, parameter :: svp1=611.2
real, parameter :: svp2=17.67
real, parameter :: svp3=29.65
real, parameter :: svpt0=273.15
real, parameter :: eps = 0.622
real startlat, startlon, deltalat, deltalon
IF ( nint(plvl).LT. 200) THEN
call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
ELSE
P = plvl
ENDIF
call get_storage(nint(plvl), 'TT', T, ix, jx)
call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx)
rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
call put_storage(nint(plvl), 'RH', rh, ix, jx)
end subroutine compute_rh_spechumd_upa
subroutine compute_rh_vapp_upa(ix, jx, plvl)
! Compute relative humidity from vapor pressure.
! Thanks to Bob Hart of PSU ESSC -- 1999-05-27.
use storage_module
implicit none
integer :: ix, jx
real :: plvl
real :: lat1, lon1, dx, dy
real, dimension(ix,jx) :: P, ES
real, pointer, dimension(:,:) :: T, E, RH
real, parameter :: svp1=611.2
real, parameter :: svp2=17.67
real, parameter :: svp3=29.65
real, parameter :: svpt0=273.15
real, parameter :: eps = 0.622
real :: startlat, startlon, deltalat, deltalon
allocate(RH(ix,jx))
P = plvl
call refr_storage(nint(plvl), 'TT', T, ix, jx)
call refr_storage(nint(plvl), 'VAPP', E, ix, jx)
ES=svp1*exp(svp2*(T-svpt0)/(T-svp3))
rh=min(1.E2*(P-ES)*E/((P-E)*ES), 1.E2)
call refw_storage(nint(plvl), 'RH', rh, ix, jx)
nullify(T,E)
end subroutine compute_rh_vapp_upa
subroutine compute_rh_depr(ix, jx, plvl)
! Compute relative humidity from Dewpoint Depression
use storage_module
implicit none
integer :: ix, jx
real :: plvl
real :: lat1, lon1, dx, dy
real, dimension(ix,jx) :: t, depr, rh
real, parameter :: Xlv = 2.5e6
real, parameter :: Rv = 461.5
integer :: i, j
call get_storage(nint(plvl), 'TT', T, ix, jx)
call get_storage(nint(plvl), 'DEPR', DEPR, ix, jx)
where(DEPR < 100.)
rh = exp(Xlv/Rv*(1./T - 1./(T-depr))) * 1.E2
elsewhere
rh = 0.0
endwhere
call put_storage(nint(plvl),'RH ', rh, ix, jx)
end subroutine compute_rh_depr
subroutine compute_rh_dewpt(ix,jx)
! Compute relative humidity from Dewpoint
use storage_module
implicit none
integer :: ix, jx
real :: lat1, lon1, dx, dy
real, dimension(ix,jx) :: t, dp, rh
real, parameter :: Xlv = 2.5e6
real, parameter :: Rv = 461.5
call get_storage(200100, 'TT ', T, ix, jx)
call get_storage(200100, 'DEWPT ', DP, ix, jx)
rh = exp(Xlv/Rv*(1./T - 1./dp)) * 1.E2
call put_storage(200100,'RH ', rh, ix, jx)
end subroutine compute_rh_dewpt
subroutine vntrp(plvl, maxlvl, k, name, ix, jx)
use storage_module
implicit none
integer :: ix, jx, k, maxlvl
real, dimension(maxlvl) :: plvl
character(len=8) :: name
real, dimension(ix,jx) :: a, b, c
real :: frc
write(*,'("Interpolating to fill in ", A, " at level ", I8)') trim(name), nint(plvl(k))
call get_storage(nint(plvl(k-1)), name, a, ix, jx)
call get_storage(nint(plvl(k+1)), name, c, ix, jx)
frc = (plvl(k) - plvl(k+1)) / ( plvl(k-1)-plvl(k+1))
b = (1.-frc)*a + frc*c
!KWM b = 0.5 * (a + c)
call put_storage(nint(plvl(k)), name, b, ix, jx)
end subroutine vntrp
subroutine fix_gfs_rh (ix, jx, plvl)
! This routine replaces GFS RH (wrt ice) with RH wrt liquid (which is what is assumed in real.exe).
use storage_module
implicit none
integer :: ix, jx, i, j
real :: plvl, eis, ews, r
real, allocatable, dimension(:,:) :: rh, tt
allocate(rh(ix,jx))
allocate(tt(ix,jx))
call get_storage(nint(plvl), 'RH', rh, ix, jx)
call get_storage(nint(plvl), 'TT', tt, ix, jx)
do j = 1, jx
do i = 1, ix
if ( tt(i,j) .le. 273.15 ) then
! Murphy and Koop 2005 ice saturation vapor pressure.
! eis and ews in hPA, tt is in K
eis = .01 * exp (9.550426 - (5723.265 / tt(i,j)) + (3.53068 * alog(tt(i,j))) &
- (0.00728332 * tt(i,j)))
! Bolton 1980 liquid saturation vapor pressure. For water saturation, most
! formulae are very similar from 0 to -20, so we don't need a more exact formula.
ews = 6.112 * exp(17.67 * (tt(i,j)-273.15) / ((tt(i,j)-273.15)+243.5))
if ( tt(i,j) .gt. 253.15 ) then
! A linear approximation to the GFS blending region ( -20 > T < 0 )
r = ((273.15 - tt(i,j)) / 20.)
r = (r * eis) + ((1-r)*ews)
else
r = eis
endif
rh(i,j) = rh(i,j) * (r / ews)
endif
enddo
enddo
call put_storage(nint(plvl), 'RH', rh, ix, jx)
deallocate (rh)
deallocate (tt)
end subroutine fix_gfs_rh