Skip Navigation Links www.nws.noaa.gov
NOAA logo - Click to go to the NOAA homepage National Weather Service NWS logo - Click to go to the NWS homepage
EMC Logo
Navigation Bar Left Cap Home News Organization
Navigation Bar End Cap

         MISSION / VISION    |    About EMC

EMC > GLOBAL BRANCH > GFS > DOCUMENTATION
GLOBAL FORECAST SYSTEM

The Global Forecast System (GFS) - Global Spectral Model (GSM)

(GFS Version 9.0.1)

Black (2003 info, no recent change) - Blue (change since 2003, previous implementations) - Red (new change as of latest implementation)

___________________________________________________________________________________________________________________

CURRENT MODEL DOCUMENTATION (as of 5/9/11):


Numerical/Computational Properties:

Dynamical/Physics Properties:
___________________________________________________________________________________________________________________

PREVIOUS MODEL DOCUMENTATION:


___________________________________________________________________________________________________________________


  • Numerical/Computational Properties


    • Horizontal Representation


      • Spectral (spherical harmonic basis functions) with transformation to a Gaussian grid for calculation of nonlinear quantities and physics.

    • Horizontal Resolution


      • Spectral triangular

      Forecast periodResolutionMaster & Flux file resolutionSigma file resolution
      Week one forecast, 0-192 hrs574 (T574)1760x880T574
      Week two forecast, 192-384 hrs190 (T190)576x288T190

    • Vertical Domain


      • The vertical domain is a sigma pressure hybrid coordinate system with 64 layers. The sigma pressure hybrid coordinate involves near-surface sigma levels, blended sigma/contant pressure levels in the mid-atmosphere, and constant pressure levels above. For more information see Office Note #461: The Implementation Of The Sigma Pressure Hybrid Coordinate Into The GFS and Office Note #462: The Derivation of the Sigma Pressure Hybrid Coordinate Semi-Lagrangian Model Equations for the GFS.

    • Vertical Representation


      • Sigma coordinate. Lorenz grid. Quadratic-conserving finite difference scheme by Arakawa and Mintz (1974).

    • Vertical Resolution


      • 64 unequally-spaced sigma levels. For a surface pressure of 1000 hPa, 15 levels are below 800 hPa, and 24 levels are above 100 hPa with the top level at .3 hPa.

    • Computer/Operating System


      • IBM RS/6000 SP (Class VIII) in an AIX environment.

    • Computational Performance


      • Operations: 30 nodes, ~8.5 minutes
      Parallel: 12 nodes, ~11 minutes

    • Initialization


      • Initialization is not necessary because the statistical spectral interpolation analysis scheme eliminates the unbalanced initial state.

    • Time Integration Scheme(s)


      • The main time integration is leapfrog for nonlinear advection terms, and semi-implicit for gravity waves and for zonal advection of vorticity and moisture. An Asselin (1972) time filter is used to reduce computational modes. The dynamics and physics are split. The physics are written in the form of an adjustment and executed in sequence. For physical processes, implicit integration with a special time filter (Kalnay and Kanamitsu, 1988) is used for vertical diffusion. In order to incorporate physical tendencies into the semi-implicit integration scheme, a special adjustment scheme is performed (Kanamitsu et al., 1991). The time step is 7.5 minutes for computation of dynamics and physics, except that the full calculation of longwave radiation is done once every 3 hours and shortwave radiation every hour (but with corrections made at every time step for diurnal variations in the shortwave fluxes and in the surface upward longwave flux).

    • Smoothing/Filling


      • Mean orographic heights on the Gaussian grid are used (see Orography). Negative atmospheric moisture values are not filled for moisture conservation, except for a temporary moisture filling that is applied in the radiation calculation.

  • Dynamical/Physical Properties


    • Atmospheric Dynamics


      • Primitive equations with vorticity, divergence, logarithm of surface pressure, specific humidity virtual temperature, and cloud condensate as dependent variables.

    • Horizontal Diffusion


      • Scale-selective, second-order horizontal diffusion after Leith (1971) is applied to vorticity, divergence, virtual temperature, and specific humidity and cloud condensate. The diffusion of temperature, specific humidity, and cloud condensate are performed on quasi-constant pressure surfaces (Kanamitsu et al. 1991).

    • Vertical Diffusion


      • The nonlocal planetary boundary layer (PBL) scheme in the NCEP GFS [so called MRF (Medium-range Forecast model) PBL model] was originally proposed by Troen and Mahrt (1986) and implemented by Hong and Pan (1996). The scheme is a first-order vertical diffusion scheme. There is a diagnostically determined PBL height that uses the bulk-Richardson approach to iteratively estimate a PBL height starting from the ground upward. Once the PBL height is determined, the profile of the coefficient of diffusivity is specified as a cubic function of the height. The actual values of the coefficients are determined by matching with the surface-layer fluxes. There is also a counter-gradient flux parameterization that is based on the fluxes at the surface and the convective velocity scale.

      In the recent update (Han and Pan, 2010), a stratocumulus top driven vertical diffusion scheme is incorporated into the MRF PBL model to increase vertical diffusion in the cloudy region of the lower troposphere. The stratocumulus top driven diffusion is further enhanced when the condition for cloud top entrainment instability is met. For the nighttime stable PBL, the revised scheme uses a local diffusivity rather than a surface layer stability based diffusivity profile used in the MRF PBL scheme. On the other hand, the background diffusivity for heat and moisture is exponentially decreased with height. To avoid too much erosion of stratocumulus along the costal area, the background diffusivity in the lower inversion layers is further reduced to 30% of that at the surface. Background diffusivity for momentum has been substantially increased to 3.0 m2s-1 everywhere, which helped reduce the wind forecast errors significantly.
       See also: Planetary Boundary Layer

    • Gravity-wave Drag Parameterization


      • Gravity-wave drag is simulated as described by Alpert et al. (1988). The parameterization includes determination of the momentum flux due to gravity waves at the surface, as well as at higher levels. The surface stress is a nonlinear function of the surface wind speed and the local Froude number, following Pierrehumbert (1987). Vertical variations in the momentum flux occur when the local Richardson number is less than 0.25 (the stress vanishes), or when wave breaking occurs (local Froude number becomes critical); in the latter case, the momentum flux is reduced according to the Lindzen (1981) wave saturation hypothesis. Modifications are made to avoid instability when the critical layer is near the surface, since the time scale for gravity-wave drag is shorter than the model time step (see also Time Integration Schemes and Orography). The treatment of the gravity-wave drag parameterization in the lower troposphere is improved by the use of the Kim and Arakawa (1995) enhancement. Included is a dependence of variance on wind direction relative to the mountain as well as subgrid statisical details of mountain distribution. This improves the prediction of lee cyclogenesis and the accompanying movement of cold outbreaks (Alpert et al, 199x).

      The gravity wave drag (GWD) code is modified to automatically scale mountain blocking and GWD stress with model resolution change (T382L64 -> T574L64). Compared to the T382L64 version of the GFS, the T574L64 GFS uses four times stronger mountain blocking and one half the strength of GWD.

    • Radiation


      • Longwave (LW) radiation:


        The longwave (LW) radiation in NCEP's operational GFS employs a Rapid Radiative Transfer Model (RRTM) developed at AER (Mlawer et al. 1997).
        • The parameterization scheme uses a correlated-k distribution method and a linear-in-tau transmittance table look-up to achieve high accuracy and efficiency.
        • The algorithm contains 140 unevenly distributed intervals (g-point) in 16 broad spectral bands.
        • In addition to the major atmospheric absorbing gases of ozone, water vapor, and carbon dioxide, the algorithm also includes various minor absorbing species such as methane, nitrous oxide, oxygen, and up to four types of halocarbons (CFCs).
        • In water vapor continuum absorption calculations, RRTM-LW employs an advanced CKD_2.4 scheme (Clough et al. 1992).
        • A maximum-random cloud overlapping method is used in the GFS application.
        • Cloud liquid/ice water path and effective radius for liquid water and ice are used for calculation of cloud-radiative properties. Hu and Stamnes' method (1993) is used to treat liquid water clouds, while Ebert and Curry's method (1992) is used for ice cloud.

      • Shortwave (SW) radiation:


        The shortwave (SW) radiation parameterization now employs the Rapid Radiative Transfer Model version 2 (RRTM2) developed at AER Inc. and updated/modified at NCEP.
        • A maximum-random cloud overlap scheme is used with this parameterization, consistent with the maximum-random overlap used in the current operational RRTM longwave (LW) radiation parameterization.
        • The aerosol SW single scattering albedo and asymmetry factor now reflect more recent data.

      • For both LW and SW:


        • Both LW and SW radiation parameterizations take into account the existence of atmospheric aerosols; both in the troposphere and stratosphere (capable of handling volcanic aerosols).

      • Additional upgrades/changes due to T574 implementation:


        • The frequency at which the LW radiation is computed was changed from three hours to one hour.
          • With this change and the change in cloud overlap for SW, both solar and infrared radiations now see the clouds consistently.
        • Additional upgrade to radiation include using more realistic time varying observed global mean CO2 (previous operational model used a constant value).
        • A new scheme was added to treat the dependence of direct-beam surface albedo on solar zenith angle over snow-free land surface.
        • Minimum value of specific humidity used in radation changed from 1.0e-5 to 1.0e-20.


    • Convection


      • Deep Convection


        • Penetrative convection is simulated following Pan and Wu (1994), which is based on Arakawa and Schubert(1974) as simplified by Grell (1993) and with a saturated downdraft. Convection occurs when the cloud work function (CWF) exceeds a certain threshold. Mass flux of the cloud is determined using a quasi-equilibrium assumption based on this threshold CWF. The CWF is a function of temperature and moisture in each air column of the model gridpoint. The temperature and moisture profiles are adjusted towards the equilibrium CWF within a specified time scale using the deduced mass flux. A major simplification of the original Arakawa-Shubert scheme is to consider only the deepest cloud and not the spectrum of clouds. The cloud model incorporates a downdraft mechanism as well as the evaporation of precipitation. Entrainment of the updraft and detrainment of the downdraft in the sub-cloud layers are included. Downdraft strength is based on the vertical wind shear through the cloud. The critical CWF is a function of the cloud base vertical motion. As the large-scale rising motion becomes strong, the CWF [similar to convective available potential energy (CAPE)] is allowed to approach zero (therefore approaching neutral stability).

        Mass fluxes induced in the updraft and the downdraft are allowed to transport momentum. The momentum exchange is calculated through the mass flux formulation in a manner similar to that for heat and moisture. The effect of the convection-induced pressure gradient force on cumulus momentum transport is parameterized in terms of mass flux and vertical wind shear (Han and Pan, 2006). As a result, the cumulus momentum exchange is reduced by about 55 % compared to the full exchange.

        The entrainment rate in cloud layers is dependent upon environmental humidity (Han and Pan, 2010). A drier environment increases the entrainment, suppressing the convection. The entrainment rate in sub-cloud layers is given as inversely proportional to height. The detrainment rate is assumed to be a constant in all layers and equal to the entrainment rate value at cloud base, which is O(10-4). The liquid water in the updraft layer is assumed to be detrained from the layers above the level of the minimum moist static energy into the grid-scale cloud water with conversion parameter of 0.002 m-1, which is same as the rain conversion parameter.

        Following Han and Pan (2010), the trigger condition is that a parcel lifted from the convection starting level without entrainment must reach its level of free convection within 120-180 hPa of ascent, proportional to the large-scale vertical velocity. This is intended to produce more convection in large-scale convergent regions but less convection in large-scale subsidence regions. Another important trigger mechanism is to include the effect of environmental humidity in the sub-cloud layer, taking into account convection inhibition due to existence of dry layers below cloud base. On the other hand, the cloud parcel might overshoot beyond the level of neutral buoyancy due to its inertia, eventually stopping its overshoot at cloud top. The CWF is used to model the overshoot. The overshoot of the cloud top is stopped at the height where a parcel lifted from the neutral buoyancy level with energy equal to 10% of the CWF would first have zero energy.

        Deep convection parameterization (SAS) modifications include:
        • Detraining cloud water from every updraft layer
        • Starting convection from the level of maximum moist static energy within PBL
        • Random cloud top is eliminated and only deepest cloud is considered
        • Cloud water is detrained from every cloud layer
        • Finite entrainment and detrainment rates for heat, moisture, and momentum are specified
        • Similar to shallow convection scheme,
          • entrainment rate is given to be inversely proportional to height in sub-cloud layers
          • detrainment rate is set to be a constant as entrainment rate at the cloud base.
        • Above cloud base, an organized entrainment is added, which is a function of environmental relative humidity

      • Shallow convection


        • The shallow convection (SC) scheme (Han and Pan, 2010) uses a bulk mass-flux parameterization, similar to the deep convection scheme. While the SC scheme is developed based on the deep convection scheme, many parameters such as cloud base mass flux, entrainment and detrainment, are modified to accommodate the SC. A cloud thickness criterion distinguishes shallow from deep convection. Deep convection is checked first: if the cloud is thicker than 150 hPa, deep convection is activated; otherwise the convection is treated as shallow. Cloud top in the SC is limited to P/Ps=0.7 (where P is the layer pressure and the subscript s represents the ground surface).

        The entrainment rate in the SC cloud layers is given to be inversely proportional to height and much smaller than that in the deep convection. The detrainment rate is assumed to be a constant in all layers and equal to the entrainment rate value at cloud base. In this way, the mass flux decreases with height above the cloud base, consistent with the large-eddy simulation studies. The liquid water is assumed to be detrained from every updraft layer. It is assumed that there exist no convective scale downdrafts. The precipitation process is allowed, but precipitation by the SC is quite small.

        Mass flux at cloud base is given as a function of the surface buoyancy flux. This differs from the deep convection scheme, which uses a quasi-equilibrium closure where the destabilization of an air column by the large-scale atmosphere is nearly balanced by the stabilization due to the cumulus.

        A new mass flux based shallow convection scheme was added to the model during the T574 implementation in 2010.
        • The shallow convection is based on the deep convection scheme (SAS), but with...
          • Cloud top limitied to 700 hPa.
          • Entrainment rate assumed to be inversely proportional to height.
          • Detrainment rate is set to a constant as the entrainment rate at the cloud base.
        • The cloud base mass flux is given as a function of convective boundary layer velocity scale.

    • Cloud Fraction


      • Following Xu and Randall (1996), the fractional cloud cover within grid box (σ) is given by

      where RH is the environmental relative humidity, ql the liquid water mixing ratio, qs the saturation specific humidity, k1, k2 and k3 the empirical coefficients. Following Xu and Randall, the values of k1=0.25, k2=100, and k3=0.49 are used.
      (See also Radiation)

    • Grid-scale Condensation and Precipitation


      • The prognostic cloud condensate has two sources, namely convective detrainment (see convection) and grid-scale condensation. The grid-scale condensation is based on Zhao and Carr (1997), which in turn is based on Sundqvist et al. (1989). The sinks of cloud condensate are grid-scale precipitation which is parameterized following Zhao and Carr (1997) for ice, and Sundqvist et al. (1989) for liquid water, and evaporation of the cloud condensate which also follows Zhao and Carr (1997).

      Evaporation of rain in the unsaturated layers below the level of condensation is also taken into account. All precipitation that penetrates the bottom atmospheric layer is allowed to fall to the surface (see also Snow Cover).

    • Planetary Boundary Layer


      • A new scheme based on the Troen and Mahrt (1986) paper was implemented on 25 October, 1995. The scheme is still a first-order vertical diffusion scheme. There is a diagnostically determined pbl height that uses the bulk-Richardson approach to iteratively estimate a pbl height starting from the ground upward. Once the pbl height is determined, the profile of the coefficient of diffusivity is specified as a cubic function of the pbl height. The actual values of the coefficients are determined by matching with the surface-layer fluxes. There is also a counter-gradient flux parameterization that is based on the fluxes at the surface and the convective velocity scale. (See Hong and Pan(1996) for a description of the scheme as well as a description of the convection scheme in the model).

      Boundary layer model significantly updated during T574 implementation to include:
      • stratocumulus top driven turbulence mixing
      • stratocumulus top driven diffusion when condition for cloud top entrainment instability is met
      • the use of local diffusion for the nighttime stable PBL

      The background vertical diffusion for heat and tracers in inversion layers below 2.5 km over ocean is reduced by 70% to decrease the erosion of stratocumulus along the coastal area.
      • A constant value of 3.0 m^2/s background vertical diffusion is added to momentum at all levels.


    • Orography


      • New orography data sets are constructed based on a United States Geological Survey (USGS) global digital elevation model (DEM) with a horizontal grid spacing of 30 arc seconds (approximately 1 km). Orography statistics including average height, mountain variance, maximum orography, land-sea-lake masks are directly derived from a 30-arc second DEM for a given resolution. See NCEP Office Note 424 (Hong, 1999) for more details. (see also Gravity-wave Drag).

      How the model handles orography/topography:

      The GFS is a global spectral model, terrain following, meaning that the topography is not actually used directly in the integration. Currently the model is T574 in operations and a corresponding physics grid of 1760x880 Gaussian grid points. The topography is originally from the USGS 30" elevations, the center of each physics grid Gaussian grid point, up to the next surrounding boxes, is used to find all the 30" grid point elevations in that box, which are averaged and used as the elevation value across each Gaussian grid box at the center point. Thus a 1760x880 grid is close to 0.2 degrees and the topography is calculated offline from model runs. This is done for any resolution that the model runs at although the present production is as stated above. Other fields like the variance of the elevation and other moments are also calculated from the USGS 30" grid and the interpolations/averaging are done to the Gaussian grid. A land sea mask is also created similarly with lakes added from the University Maryland vegetative index file which contains a global 30" lake mask and Antarctica elevations are from a RADSAT data (instead of USGS). There is also a spectral quadratic filter applied to the elevations so there is a smoothing for the highest spectral waves, starting at T384.

      The model takes the Gaussian grid and transforms the elevations to spectral space if it is the first time used, or if there is a change (like a resolution change) from previous runs and writes the new spectralized topography to the initial condition sigma-spectral file, otherwise the model expects the spectral topography to be present in its initial condition file and is used to set the surface pressure hydrostatically. Thus, normally the model starts from the initial condition sigma vertical coordinate spectral coefficient file containing the topography in spectral space coefficients (of spherical harmonics) for the first time step which may not change until the model resolution is changed. After each return to physical space, all the dependent variables as well as the surface pressure, (from the topography) are calculated on a reduced physical space Gaussian grid. The reduced grid reduces the physics grid used in all the physical calculations by a preset number of longitudes per latitude, with no reduction at the equator to a maximum reduction near each pole. The process is designed to deliver all the spectral information of the spectral T574 to a reduced grid that is tuned to receive the full information content from the model dynamics step. All physical grid data is returned to the original full Gaussian grid before writing out including the surface dynamic height and surface pressure.

    • Ocean


      • A daily OI sea surface temperature analysis that assimilates observations from past seven days is used (Reynolds and Smith, 1994, available here ). The sea surface temperature anomaly is damped with an e-folding time of 90 days during the course of the forecast.

    • Relative Humidity


      • Relative humidity in the GFS is computed the following way:

      1) at T < -20 C, w.r.t. ice
      2) at T > 0 C, w.r.t. water
      3) -20 < T < 0, mix phased.

    • Sea Ice


      • Sea-ice is obtained from the analysis by the Marine Modeling Branch, available daily. The sea ice is assumed to have a constant thickness of 3 meters, and the ocean temperature below the ice is specified to be 271.2 K. The surface temperature of sea ice is determined from an energy balance that includes the surface heat fluxes (see Surface Fluxes) and the heat capacity of the ice. Snow accumulation does not affect the albedo or the heat capacity of the ice.

    • Snow Cover


      • Snow cover is obtained from an analysis by NESDIS (the IMS system) and the Air Force, updated daily. When the snow cover analysis is not available, the predicted snow in the data assimilation is used. Precipitation falls as snow if the temperature at sigma=.85 is below 0 C. Snow mass is determined prognostically from a budget equation that accounts for accumulation and melting. Snow melt contributes to soil moisture, and sublimation of snow to surface evaporation. Snow cover affects the surface albedo and heat transfer/capacity of the soil, but not of sea ice. See also Sea Ice, Surface Characteristics, Surface Fluxes, and Land Surface Processes.

    • Surface Characteristics


      • Roughness lengths over oceans are determined from the surface wind stress after the method of Charnock (1955). Over sea ice the roughness is a uniform 0.01 cm. Roughness lengths over land are prescribed from data of Dorman and Sellers (1989) which include 12 vegetation types. Note that the surface roughness is not a function of orography. Over oceans the surface albedo depends on zenith angle. The albedo of sea ice is a function of surface skin temperature and nearby atmospheric temperature as well as snow cover (Grumbine, 1994), with values ranging from 0.65-0.8 for snow-covered sea ice and from 0.45-0.65 for bare sea ice. Albedoes for land surfaces are based on Matthews (1985) surface vegetation distribution (See Radiation).

      Longwave emissivity is prescribed to be unity (black body emission) for all surfaces. Soil type and Vegetation type data base from GCIP is used. Vegetation fraction monthly climatology based on NESDIS NDVI 5-year climatology is used.

    • Surface Fluxes


      • Surface solar absorption is determined from the surface albedos, and longwave emission from the Planck equation with emissivity of 1.0 (see Surface Characteristics).
      • The lowest model layer is assumed to be the surface layer (sigma=0.996) and the Monin-Obukhov similarity profile relationship is applied to obtain the surface stress and sensible and latent heat fluxes. The formulation was based on Miyakoda and Sirutis (1986) and has been modified by P. Long in the very stable and very unstable situations. A bulk aerodynamic formula is used to calculate the fluxes once the turbulent exchange coefficients have be obtained.
      • Roughness length over ocean is updated with a Charnock formula after surface stress has been obtained.
      • Thermal roughness over the ocean is based on a formulation derived from TOGA COARE(Zeng et al, 1998). A new thermal roughness scheme was added during the 9.1.0 upgrade. This new scheme gives the different values to the thermal roughness length and the momentum roughness length, thus yielding lower and more realistic canopy conductance compared to the old scheme. Greenness vegetation fraction is used for weighting on the computation of both thermal and momentum roughness length to achieve the turbulent formulation convergence from a fully-covered vegetation area to bare soil.
      • Land surface evaporation is comprised of three components: direct evaporation from the soil and from the canopy, and transpiration from the vegetation. The formulation follows Pan and Mahrt (1987).

    • Land Surface Processes


      • Soil temperature and soil volumetric water content are computed in two layers at depths 0.1 and 1.0 meters by a fully implicit time integration scheme (Pan and Mahrt, 1987). For sea ice, the layer depths were specified as 1.5 and 3 meters. Heat capacity, thermal and hydraulic diffusivity and hydraulic conductivity coefficients are strong functions of the soil moisture content. A climatological deep-soil temperature is specified at the third layer of 4 meters for soil and a constant value of 272 K is specified as the ice-water interface temperature for sea ice. The vegetation canopy is allowed to intercept precipitation and re-evaporation. Runoff from the surface and drainage from the bottom layer are also calculated.

    • Chemistry


      • Ozone is a prognostic variable that is updated in the analysis and transported in the model. The sources and sinks of ozone are computed using zonally averaged seasonally varying production and destruction rates provided by NASA/GSFC.