subroutine gwdc::gwdc_run	(	integer, intent(in)	im,
		integer, intent(in)	ix,
		integer, intent(in)	km,
		integer, intent(in)	lat,
		real(kind=kind_phys), dimension(ix,km), intent(in)	u1,
		real(kind=kind_phys), dimension(ix,km), intent(in)	v1,
		real(kind=kind_phys), dimension(ix,km), intent(in)	t1,
		real(kind=kind_phys), dimension(ix,km), intent(in)	q1,
		real(kind=kind_phys), intent(in)	deltim,
		real(kind=kind_phys), dimension(ix,km), intent(in)	pmid1,
		real(kind=kind_phys), dimension(ix,km+1), intent(in)	pint1,
		real(kind=kind_phys), dimension(ix,km), intent(in)	dpmid1,
		real(kind=kind_phys), dimension(im), intent(in)	qmax,
		integer, dimension(im), intent(in)	ktop,
		integer, dimension(im), intent(in)	kbot,
		integer, dimension(im), intent(in)	kcnv,
		real(kind=kind_phys), dimension(im), intent(in)	cldf,
		real(kind=kind_phys), intent(in)	grav,
		real(kind=kind_phys), intent(in)	cp,
		real(kind=kind_phys), intent(in)	rd,
		real(kind=kind_phys), intent(in)	fv,
		real(kind=kind_phys), intent(in)	pi,
		real(kind=kind_phys), dimension(im), intent(in)	dlength,
		logical, intent(in)	lprnt,
		integer, intent(in)	ipr,
		real(kind=kind_phys), intent(in)	fhour,
		real(kind=kind_phys), dimension(ix,km), intent(out)	utgwc,
		real(kind=kind_phys), dimension(ix,km), intent(out)	vtgwc,
		real(kind=kind_phys), dimension(im), intent(out)	tauctx,
		real(kind=kind_phys), dimension(im), intent(out)	taucty,
		character(len=*), intent(out)	errmsg,
		integer, intent(out)	errflg
	)

• Determine if deep convection occurs and activate convection-induced GWD scheme.
• Create local arrays with reversed vertical indices and Initialize local variables.
  • The top interface temperature, density, and Brunt-Vaisala frequencies ( $N$) are calculated assuming an isothermal atmosphere above the top mid level.
  • The bottom interface temperature, density, and Brunt-Vaisala frequencies ( $N$) are calculated assuming an isothermal atmosphere below the bottom mid level.
  • The interface level temperature, density, and Brunt-Vaisala frequencies ( $N$) are calculated based on linear interpolation of temperature in ln(P).
  • The mid-level Brunt-Vaisala frequencies ( $N$) are calculated based on interpolated interface temperatures.
• Calculate the cloud top wind components, speed and direction.
• Calculate the basic state wind profiles projected in the direction of the cloud top wind at mid level and interface level.
• Calculate the local Richardson number:

  $$Ri=N^2/\eta^2$$

  where $\eta$ is the vertical shear ( $dU/dz$).
• Calculate the gravity wave stress at the interface level cloud top.
  • Wave stress at cloud top is calculated when the atmosphere is dynamically stable at the cloud top
  • The cloud top wave stress and nonlinear parameter are calculated using density, temperature, and wind that are defined at mid level just below the interface level in which cloud top wave stress is defined. The parameter $\mu$ is the nonlinearity factor of thermally induced internal gravity waves defined by eq.(17) in Chun and Baik (1998) [18] :

    $$\mu=\frac{gQ_{0}a_{1}}{c_{p}T_{0}NU^{2}}$$

    where $Q_{0}$ is the maximum deep convective heating rate in a horizontal grid point calculated from cumulus parameterization. $a_{1}$ is the half-width of the forcing function. $g$ is gravity. $c_{p}$ is specific heat at constant pressure. $T_{0}$ is the layer mean temperature (T1). As eqs.(18) and (19) Chun and Baik (1998) [18], the zonal momentum flux is given by

    $$\tau_{x}=-[\rho U^{3}/(N\triangle x)]G(\mu)$$

    where

    $$G(\mu)=c_{1}c_2^2 \mu^{2}$$

    wher $\rho$ is the local density. The tunable parameter $c_1$ is related to the horizontal structure of thermal forcing. The tunable parameter $c_2$ is related to the basic-state wind and stability and the bottom and top heights of thermal forcing. If the atmosphere is dynamically unstable at the cloud top, the convective GWD calculation is skipped at that grid point.
  • The stress is capped at tauctmax = - 5 $n/m^2$ in order to prevent numerical instability.
    • The stress is capped at $\tau_{x} = - 5n/m^2$ in order to prevent numerical instability.
  • If the atmosphere is dynamically unstable at the cloud top, GWDC calculation in current horizontal grid is skipped.
  • If mean wind at the cloud top is less than zero, GWDC calculation in current horizontal grid is skipped.
• Calculate the minimum Richardson number including both the basic-state condition and wave effects.

  $$Ri_{min}\approx\frac{Ri(1-\mu|c_{2}|)}{(1+\mu Ri^{1/2}|c_{2}|)^{2}}$$

• Calculate the gravity wave stress profile using the wave saturation hypothesis of Lindzen (1981) [67].
  • When $Ri_{min}$ is set to 1/4 based on Lindzen's (1981) [67] saturation hypothesis, the nonlinearity factor for wave saturation can be derived by

    $$\mu_{s}=\frac{1}{|c_{2}|}[2\sqrt{2+\frac{1}{\sqrt{Ri}}}-(2+\frac{1}{\sqrt{Ri}})]$$

    Then the saturation zonal momentum flux is given by

    $$\tau_{s}=-[\rho U^{3}/(N\triangle x)]c_{1}c_2^2\mu_s^2$$

  • If the minimum $R_{i}$ at interface cloud top is less than or equal to 1/4, the convective GWD calculation is skipped at that grid point.
  • As an upper boundary condition, upward propagation of gravity wave energy is permitted.
• Calculate wind tendency in direction to the wind vector,zonal wind tendency and meridional wind tendency above the cloud top level due to convectively generated gravity waves.
• Convert back local convective GWD tendency arrays to GFS model vertical indices.