From 12c6370379becea4922781800bff4c187e34a90d Mon Sep 17 00:00:00 2001 From: Jim Bresch Date: Sat, 19 Mar 2011 00:49:24 +0000 Subject: [PATCH] Add vtables for ERA-interim files. Update g1print for ECMWF table128 Update rd_grib1 for ECMWF earth radius. It's unclear what they actually use, so code it to the standard. Update rrpr so RH is computed on model levels only if we have all the required variables. Allow SEAICE to be a fraction or flag. Calculate SNOWH in rrpr when we have model-dependent information. SNOWH must be in the Vtable for this to occur, otherwise, real.exe's default behavior is used (correct only for NAM). git-svn-id: https://svn-wrf-wps.cgd.ucar.edu/trunk@583 86b71a92-4018-0410-97f8-d555beccfc3a --- ungrib/Variable_Tables/Vtable.ERA-interim.ml | 42 +++++++++++++ ungrib/Variable_Tables/Vtable.ERA-interim.pl | 45 ++++++++++++++ ungrib/src/g1print.F | 49 +++++++++++++++- ungrib/src/rd_grib1.F | 14 ++++- ungrib/src/rrpr.F | 88 +++++++++++++++++++++++----- 5 files changed, 218 insertions(+), 20 deletions(-) create mode 100644 ungrib/Variable_Tables/Vtable.ERA-interim.ml create mode 100644 ungrib/Variable_Tables/Vtable.ERA-interim.pl diff --git a/ungrib/Variable_Tables/Vtable.ERA-interim.ml b/ungrib/Variable_Tables/Vtable.ERA-interim.ml new file mode 100644 index 0000000..5f077b3 --- /dev/null +++ b/ungrib/Variable_Tables/Vtable.ERA-interim.ml @@ -0,0 +1,42 @@ +GRIB | Level| Level| Level| metgrid | metgrid | metgrid | +Code | Code | 1 | 2 | Name | Units | Description | +-----+------+------+------+----------+----------+------------------------------------------+ + 130 | 109 | * | | TT | K | Temperature | + 131 | 109 | * | | UU | m s-1 | U | + 132 | 109 | * | | VV | m s-1 | V | + 133 | 109 | * | | SPECHUMD | kg kg-1 | Specific humidity | + 165 | 1 | 0 | | UU | m s-1 | U | At 10 m + 166 | 1 | 0 | | VV | m s-1 | V | At 10 m + 167 | 1 | 0 | | TT | K | Temperature | At 2 m + 168 | 1 | 0 | | DEWPT | K | | At 2 m + | 1 | 0 | | RH | % | Relative Humidity at 2 m | At 2 m + 172 | 1 | 0 | | LANDSEA | 0/1 Flag | Land/Sea flag | + 129 | 1 | 0 | | SOILGEO | m2 s-2 | | + | 1 | 0 | | SOILHGT | m | Terrain field of source analysis | + 134 | 1 | 0 | | PSFC | Pa | Surface Pressure | + 151 | 1 | 0 | | PMSL | Pa | Sea-level Pressure | + 235 | 1 | 0 | | SKINTEMP | K | Sea-Surface Temperature | + 31 | 1 | 0 | | SEAICE | fraction | Sea-Ice-Fraction | + 34 | 1 | 0 | | SST | K | Sea-Surface Temperature | + 33 | 1 | 0 | | SNOW_DEN | kg m-3 | | + 141 | 1 | 0 | | SNOW_EC | m | | + | 1 | 0 | | SNOW | kg m-2 |Water Equivalent of Accumulated Snow Depth| + | 1 | 0 | | SNOWH | m | Physical Snow Depth | + 139 | 112 | 0 | 7 | ST000007 | K | T of 0-7 cm ground layer | + 170 | 112 | 7 | 28 | ST007028 | K | T of 7-28 cm ground layer | + 183 | 112 | 28 | 100 | ST028100 | K | T of 28-100 cm ground layer | + 236 | 112 | 100 | 255 | ST100255 | K | T of 100-255 cm ground layer | + 39 | 112 | 0 | 7 | SM000007 | m3 m-3 | Soil moisture of 0-7 cm ground layer | + 40 | 112 | 7 | 28 | SM007028 | m3 m-3 | Soil moisture of 7-28 cm ground layer | + 41 | 112 | 28 | 100 | SM028100 | m3 m-3 | Soil moisture of 28-100 cm ground layer | + 42 | 112 | 100 | 255 | SM100255 | m3 m-3 | Soil moisture of 100-255 cm ground layer | +-----+------+------+------+----------+----------+------------------------------------------+ +# +# For use with ERA-interim model-level output. +# +# Grib codes are from Table 128 +# http://www.ecmwf.int/services/archive/d/parameters/order=grib_parameter/table=128/ +# +# snow depth is converted to the proper units in rrpr.F +# +# For ERA-interim data at NCAR, use the ml (sc and uv) and sfc sc files. diff --git a/ungrib/Variable_Tables/Vtable.ERA-interim.pl b/ungrib/Variable_Tables/Vtable.ERA-interim.pl new file mode 100644 index 0000000..9efd93b --- /dev/null +++ b/ungrib/Variable_Tables/Vtable.ERA-interim.pl @@ -0,0 +1,45 @@ +GRIB | Level| Level| Level| metgrid | metgrid | metgrid | +Code | Code | 1 | 2 | Name | Units | Description | +-----+------+------+------+----------+----------+------------------------------------------+ + 129 | 100 | * | | GEOPT | m2 s-2 | | + | 100 | * | | HGT | m | Height | + 130 | 100 | * | | TT | K | Temperature | + 131 | 100 | * | | UU | m s-1 | U | + 132 | 100 | * | | VV | m s-1 | V | + 157 | 100 | * | | RH | % | Relative Humidity | + 165 | 1 | 0 | | UU | m s-1 | U | At 10 m + 166 | 1 | 0 | | VV | m s-1 | V | At 10 m + 167 | 1 | 0 | | TT | K | Temperature | At 2 m + 168 | 1 | 0 | | DEWPT | K | | At 2 m + | 1 | 0 | | RH | % | Relative Humidity | At 2 m + 172 | 1 | 0 | | LANDSEA | 0/1 Flag | Land/Sea flag | + 129 | 1 | 0 | | SOILGEO | m2 s-2 | | + | 1 | 0 | | SOILHGT | m | Terrain field of source analysis | + 134 | 1 | 0 | | PSFC | Pa | Surface Pressure | + 151 | 1 | 0 | | PMSL | Pa | Sea-level Pressure | + 235 | 1 | 0 | | SKINTEMP | K | Sea-Surface Temperature | + 31 | 1 | 0 | | SEAICE | fraction | Sea-Ice Fraction | + 34 | 1 | 0 | | SST | K | Sea-Surface Temperature | + 33 | 1 | 0 | | SNOW_DEN | kg m-3 | | + 141 | 1 | 0 | | SNOW_EC | m | | + | 1 | 0 | | SNOW | kg m-2 |Water Equivalent of Accumulated Snow Depth| + | 1 | 0 | | SNOWH | m | Physical Snow Depth | + 139 | 112 | 0 | 7 | ST000007 | K | T of 0-7 cm ground layer | + 170 | 112 | 7 | 28 | ST007028 | K | T of 7-28 cm ground layer | + 183 | 112 | 28 | 100 | ST028100 | K | T of 28-100 cm ground layer | + 236 | 112 | 100 | 255 | ST100255 | K | T of 100-255 cm ground layer | + 39 | 112 | 0 | 7 | SM000007 | m3 m-3 | Soil moisture of 0-7 cm ground layer | + 40 | 112 | 7 | 28 | SM007028 | m3 m-3 | Soil moisture of 7-28 cm ground layer | + 41 | 112 | 28 | 100 | SM028100 | m3 m-3 | Soil moisture of 28-100 cm ground layer | + 42 | 112 | 100 | 255 | SM100255 | m3 m-3 | Soil moisture of 100-255 cm ground layer | +-----+------+------+------+----------+----------+------------------------------------------+ +# +# For use with ERA-interim pressure-level output. +# +# Grib codes are from Table 128 +# http://www.ecmwf.int/services/archive/d/parameters/order=grib_parameter/table=128/ +# +# snow depth is converted to the proper units in rrpr.F +# +# For ERA-interim data at NCAR, use the pl (sc and uv) and sfc sc files. + diff --git a/ungrib/src/g1print.F b/ungrib/src/g1print.F index 3c419bb..9e34d09 100644 --- a/ungrib/src/g1print.F +++ b/ungrib/src/g1print.F @@ -350,15 +350,26 @@ end subroutine fieldname afwa(235) = 'SSRUN' ! ECMWF -! from http://www.ecmwf.int/products/data/technical/GRIB_tables/table_128.html +! from http://www.ecmwf.int/services/archive/d/parameters/order=grib_parameter/table=128/ do i = 1, 254 ecmwf(i) = ' ' enddo + ecmwf(1) = 'STRF ' + ecmwf(2) = 'VPOT ' + ecmwf(3) = 'PT ' + ecmwf(4) = 'EQPT ' + ecmwf(5) = 'SEPT ' + ecmwf(8) = 'SRO ' + ecmwf(9) = 'SSRO ' + ecmwf(10) = 'WS ' + ecmwf(26) = 'CL ' ecmwf(27) = 'CVL ' ecmwf(28) = 'CVH ' ecmwf(29) = 'TVL ' ecmwf(30) = 'TVH ' ecmwf(31) = 'CI ' + ecmwf(32) = 'ASN ' + ecmwf(33) = 'RSN ' ecmwf(34) = 'SSTK ' ecmwf(35) = 'ISTL1' ecmwf(36) = 'ISTL2' @@ -368,9 +379,22 @@ end subroutine fieldname ecmwf(40) = 'SWVL2' ecmwf(41) = 'SWVL3' ecmwf(42) = 'SWVL4' + ecmwf(43) = 'SLT ' ecmwf(44) = 'ES ' ecmwf(45) = 'SMLT ' ecmwf(60) = 'PV ' + ecmwf(74) = 'SDFOR' + ecmwf(75) = 'CRWC ' + ecmwf(76) = 'CSWC ' + ecmwf(77) = 'ETADT' + ecmwf(78) = 'TCLW ' + ecmwf(79) = 'TCIW ' + ecmwf(121) = 'MX2T6' + ecmwf(122) = 'MN2T6' + ecmwf(123) = '10FG6' + ecmwf(124) = 'EMIS ' + ecmwf(127) = 'AT ' + ecmwf(128) = 'BV ' ecmwf(129) = 'Z ' ecmwf(130) = 'T ' ecmwf(131) = 'U ' @@ -390,9 +414,15 @@ end subroutine fieldname ecmwf(145) = 'BLD ' ecmwf(146) = 'SSHF ' ecmwf(147) = 'SLHF ' + ecmwf(148) = 'CHNK ' + ecmwf(149) = 'SNR ' + ecmwf(150) = 'TNR ' ecmwf(151) = 'MSL ' ecmwf(152) = 'LNSP ' + ecmwf(153) = 'SWHR ' + ecmwf(154) = 'LWHR ' ecmwf(155) = 'D ' + ecmwf(156) = 'GH ' ecmwf(157) = 'R ' ecmwf(159) = 'BLH ' ecmwf(160) = 'SDOR ' @@ -425,8 +455,15 @@ end subroutine fieldname ecmwf(187) = 'MCC ' ecmwf(188) = 'HCC ' ecmwf(189) = 'SUND ' + ecmwf(194) = 'BTMP ' ecmwf(195) = 'LGWS ' ecmwf(196) = 'MGWS ' + ecmwf(197) = 'GWD ' + ecmwf(198) = 'SRC ' + ecmwf(199) = 'VEG ' + ecmwf(200) = 'VSO ' + ecmwf(201) = 'MX2T ' + ecmwf(202) = 'MN2T ' ecmwf(203) = 'O3 ' ecmwf(204) = 'PAW ' ecmwf(205) = 'RO ' @@ -434,8 +471,12 @@ end subroutine fieldname ecmwf(207) = '10SI ' ecmwf(208) = 'TSRC ' ecmwf(209) = 'TTRC ' - ecmwf(210) = 'STRC ' + ecmwf(210) = 'SSRC ' ecmwf(211) = 'STRC ' + ecmwf(212) = 'TISR ' + ecmwf(213) = 'VIMD ' + ecmwf(214) = 'DHR ' + ecmwf(227) = 'CRNH ' ecmwf(229) = 'IEWS ' ecmwf(230) = 'INSS ' ecmwf(231) = 'ISHF ' @@ -444,7 +485,11 @@ end subroutine fieldname ecmwf(234) = 'LSRH ' ecmwf(235) = 'SKT ' ecmwf(236) = 'STL4 ' + ecmwf(237) = 'SWL4 ' ecmwf(238) = 'TSN ' + ecmwf(239) = 'CSF ' + ecmwf(240) = 'LSF ' + ecmwf(248) = 'CC ' end subroutine init_tables block data ptables diff --git a/ungrib/src/rd_grib1.F b/ungrib/src/rd_grib1.F index 98e678f..5c1bf7d 100644 --- a/ungrib/src/rd_grib1.F +++ b/ungrib/src/rd_grib1.F @@ -368,8 +368,15 @@ SUBROUTINE rd_grib1(IUNIT, gribflnm, level, field, hdate, & map%startloc = 'SWCORNER' map%grid_wind = .true. -! NCEP's grib1 messages say they use 6367.47, but they really use 6371.229. Hardcode it. - map%r_earth = 6371.229 +! NCEP's grib1 messages (in GDS Octet 17, the Resolution and Component Flags) +! all have '0' for the earth radius flag which the documentation (written by NCEP) +! says is 6367.47, but they really use 6371.229. Hardcode it. +! It's not clear what ECMWF uses. One place says 6367.47 and another 6371.229. + if ( index(map%source,'NCEP') .ne. 0 ) then + map%r_earth = 6371.229 + else + map%r_earth = 6367.47 + endif if (ksec2(4).eq.0) then ! Lat/Lon grid map%igrid = 0 @@ -455,7 +462,8 @@ SUBROUTINE rd_grib1(IUNIT, gribflnm, level, field, hdate, & if (tmp8(5:5) .eq. '0') map%grid_wind = .false. else - print*, 'Unknown ksec2(4): ', ksec2(4) + print*, 'Unknown Data Representation Type, ksec2(4)= ', ksec2(4) + stop 'rd_grib1' endif 111 format(' igrid : ', i3, /, & diff --git a/ungrib/src/rrpr.F b/ungrib/src/rrpr.F index d7fa542..e03eac7 100644 --- a/ungrib/src/rrpr.F +++ b/ungrib/src/rrpr.F @@ -43,7 +43,7 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f character (LEN=25) :: units character (LEN=46) :: Desc - real, allocatable, dimension(:,:) :: scr2d + real, allocatable, dimension(:,:) :: scr2d, tmp2d real, pointer, dimension(:,:) :: ptr2d integer :: k, kk, mm, n, ierr, ifv @@ -64,8 +64,6 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f ! 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 @@ -333,7 +331,8 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f do k = 1, nlvl if (plvl(k).lt.200000.) then - if (.not. is_there(nint(plvl(k)), 'RH').and. & + if (.not. is_there(nint(plvl(k)), 'RH') .and. & + is_there(nint(plvl(k)), 'TT') .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)) @@ -347,6 +346,7 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f do k = 1, nlvl if (plvl(k).lt.200000.) then if (.not. is_there(nint(plvl(k)),'RH').and. & + is_there(nint(plvl(k)), 'TT') .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)) @@ -359,6 +359,7 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f do k = 1, nlvl if (plvl(k).lt.200000.) then if (.not. is_there(nint(plvl(k)),'RH').and. & + is_there(nint(plvl(k)), 'TT') .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)) @@ -495,7 +496,7 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f endif endif -! ECMWF snow depth is in meters (Table 128). Convert to kg/m2 +! ECMWF snow depth in meters of water equivalent (Table 128). Convert to kg/m2 ! if (is_there(200100, 'SNOW_EC')) then call get_dims(200100, 'SNOW_EC') @@ -506,13 +507,57 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f deallocate(scr2d) endif -! If we've got a SEAICE field, make sure that it is all Zeros and Ones: +! Convert the ECMWF LANDSEA mask from a fraction to a flag - if (is_there(200100, 'SEAICE')) then - call get_dims(200100, 'SEAICE') - call make_zero_or_one(map%nx, map%ny) + if ( index(map%source,'ECMWF') .ne. 0) then + if (is_there(200100, 'LANDSEA')) then + call get_dims(200100, 'LANDSEA') + call make_zero_or_one(map%nx, map%ny, 'LANDSEA') + endif endif +! compute physical snow depth (SNOWH) for various models +! As of March 2011, this is done here instead of real. + if (is_there(200100, 'SNOW') .and. .not. is_there(200100, 'SNOWH')) then + call get_dims(200100, 'SNOW') + allocate(scr2d(map%nx,map%ny)) + call get_storage(200100, 'SNOW', scr2d, map%nx, map%ny) + if ( index(map%source,'NCEP ') .ne. 0) then + if ( index(map%source,'NMM ') .ne. 0 .and. & + index(map%source,'NAM ') .ne. 0 .and. & + index(map%source,'RTMA') .ne. 0 .and. & + index(map%source,'RUC ') .ne. 0 ) then ! it's a GFS variant + scr2d = scr2d * 0.01 ! NCEP uses 100:1 ratio for GFS + else + scr2d = scr2d * 0.005 ! NCEP uses 200:1 ratio for NAM + endif + else if (index(map%source,'ECMWF') .ne. 0) then + if (is_there(200100, 'SNOW_DEN')) then ! If we have snow density, use it to compute snowh + call get_dims(200100, 'SNOW_DEN') + allocate(tmp2d(map%nx,map%ny)) + call get_storage(200100, 'SNOW_DEN', tmp2d, map%nx, map%ny) + scr2d = scr2d / tmp2d + deallocate(tmp2d) + else + scr2d = scr2d * 0.004 ! otherwise, assume a density of 250 mm/m (i.e. 250:1 ratio). + endif + else ! Other models + scr2d = scr2d * 0.005 ! Use real's default method (200:1) + endif + call put_storage(200100, 'SNOWH', scr2d, map%nx, map%ny) + deallocate(scr2d) + endif + +! As of March 2011, SEAICE can be a flag or a fraction. It will be converted +! to the appropriate values in real depending on whether or not the polar mods are used. + +!! 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, 'SEAICE') +! endif + ! If we've got an ICEMASK field, re-flag it for output to met_em and real: ! Field | GRIB In | Out ! ------------------------- @@ -545,20 +590,21 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f enddo TIMELOOP end subroutine rrpr -subroutine make_zero_or_one(ix, jx) -! Make sure the SEAICE field is zero or one. +subroutine make_zero_or_one(ix, jx, infield) +! Make sure the input field (SEAICE or LANDSEA) is zero or one. use storage_module implicit none integer :: ix, jx real, dimension(ix,jx) :: seaice + character(len=*) :: infield - call get_storage(200100, 'SEAICE',seaice, ix, jx) + call get_storage(200100, infield, seaice, ix, jx) where(seaice > 0.5) seaice = 1.0 elsewhere seaice = 0.0 end where - call put_storage(200100, 'SEAICE',seaice, ix, jx) + call put_storage(200100, infield, seaice, ix, jx) end subroutine make_zero_or_one subroutine re_flag_ice_mask(ix, jx) @@ -667,7 +713,11 @@ subroutine compute_rh_spechumd_upa(ix, jx, plvl) real, parameter :: eps = 0.622 IF ( nint(plvl).LT. 200) THEN - call get_storage(nint(plvl), 'PRESSURE', P, ix, jx) + if (is_there(nint(plvl(k)), 'PRESSURE')) then + call get_storage(nint(plvl), 'PRESSURE', P, ix, jx) + else + return ! if we don't have pressure on model levels, return + endif ELSE P = plvl ENDIF @@ -698,7 +748,15 @@ subroutine compute_rh_vapp_upa(ix, jx, plvl) allocate(RH(ix,jx)) - P = plvl + IF ( nint(plvl).LT. 200) THEN + if (is_there(nint(plvl(k)), 'PRESSURE')) then + call get_storage(nint(plvl), 'PRESSURE', P, ix, jx) + else + return ! if we don't have pressure on model levels, return + endif + ELSE + P = plvl + ENDIF call refr_storage(nint(plvl), 'TT', T, ix, jx) call refr_storage(nint(plvl), 'VAPP', E, ix, jx) -- 2.11.4.GIT