Impact of Subgrid Variation of Water Vapor on Longwave Radiation in a General Circulation Model

Most general circulation models compute radiation fluxes by assuming that water vapor is uniform within individual grid layers, which leads to an underestimation of satellite‐observed longwave (LW) cloud radiative forcing (LWCF). To fix this problem, we calculated water vapor content separately for clear and cloudy portions and used them to compute LW radiation. The impacts of this modification were examined by comparing two global simulations with and without the modification (NEW and OLD, respectively). Global‐annual mean LWCF from NEW was 1.8 W m−2 higher than that of OLD, thus remedying a long‐standing negative bias of LWCF. This improvement is a combined result of more clear‐sky and less all‐sky upward LW flux at the top of the atmosphere than OLD. Large increases in LWCF and clear‐sky LW flux occurred in the tropical deep convection and midlatitude storm track regions where upper‐ and middle‐level clouds are abundant. Although only the LW radiation scheme was modified, global‐annual mean shortwave cloud radiative forcing also increased, particularly in the vicinity of the eastern subtropical marine stratocumulus decks through radiative feedback processes. With this improved treatment, it may be possible to tune general circulation models in a more flexible and physical way without introducing compensating errors.


Introduction
Mainly designed for simulating global climate over a long period, general circulation models (GCMs) are required to maintain a global radiation balance at the top of the atmosphere (TOA). This radiation balance is often achieved by adjusting low-level cloud fraction and associated shortwave (SW) radiation. However, if there are systematic biases in the simulated longwave (LW) radiation, this strategy introduces a compensating error in addition to the other existing errors, which may distort physical processes sensitive to individual SW and LW radiation, exerting undesirable effects on the simulated climate. Thus, the success of this tuning strategy depends on the accuracy of GCM-simulated LW radiation. It was reported that many GCMs, which participated in Phase 3 of the Coupled Model Intercomparison Project (CMIP3), systematically underestimate the observed downward LW radiation (DLR) at the surface (Wild, 2008). These biases were reduced from the CMIP3 to CMIP5 models (Li et al., 2013), but the majority of CMIP5 models still suffer from the positive biases in outgoing LW radiation (OLR) at TOA and the negative biases in LW cloud radiative forcing (LWCF) at TOA on global-annual mean (Calisto et al., 2014;Dolinar et al., 2015).
To correct the biases in the GCM-simulated LW radiation, many studies have focused on parameterizing LW scattering effects of cloud condensates (Chou et al., 1999;Edwards & Slingo, 1996;Ritter & Geleyn, 1992). Stephens et al. (2001) showed that OLR calculated by a nonscattering LW radiation model is approximately 8 W m −2 larger than that from a scattering one. Joseph and Min (2003) also reported that the model neglecting multiple scattering of LW radiation by thin cirrus clouds overestimates the observed OLR by 6-8 W m −2 . Costa and Shine (2006) estimated that the impact of LW scattering by cloud condensates on global OLR is about −3 W m −2 . Schmidt et al. (2006) also reported that multiple scattering of LW radiation by cloud condensates decreases the global OLR by about 1.5 W m −2 and increases the global mean DLR by about 0.4 W m −2 . Based on the analysis of observation data, Kuo et al. (2017) reported that the LW scattering by cloud condensates decreases OLR by 2.6 W m −2 and increases DLR by 1.2 W m −2 . More recently, Jin et al. (2019) reported that when the LW scattering by cloud condensates is included in a GCM, the global-annual mean OLR decreases by 2.7 W m −2 and DLR increases by 1.6 W m −2 .
Although not much investigated by previous studies, another important source for the biases in the GCM-simulated LW radiation is the neglect of the water vapor differences between the cloud and clear portions. Most GCMs compute radiation fluxes by assuming that water vapor is uniform within individual grid layers (Bentsen et al., 2013;Donner et al., 2011;Ji et al., 2014;Martin et al., 2006;Park et al., 2014;Sekiguchi & Nakajima, 2008;Stevens et al., 2013;Voldoire et al., 2013;Wu et al., 2014). For example, the radiation scheme in the Community Atmosphere Model version 5 (CAM5) generates a set of radiative subcolumns by assuming that all thermodynamic variables (e.g., water vapor, temperature, and concentrations of gases and aerosols) except the mass and number concentration of cloud condensates are horizontally uniform within each grid layer (Park et al., 2014). In reality, however, water vapor in the cloud portion is higher than that in the clear portion, such that satellite-observed LWCF and clear-sky upward LW radiation at the TOA are likely to be larger than those from GCM simulations. In spite of this apparent conceptual flaw, the GCM modeling community has been conservative in fixing this problem, presumably due in part to computational cost, complex structure of the existing radiation schemes, or the notion that this is not an important factor. However, from simple off-line GCM simulations, Park (2018) showed that accounting for subgrid fluctuations of water vapors between cloud and clear portions in the LW radiation scheme substantially increases global mean LWCF from 22.6 to 26.9 W m −2 . This change in LWCF is comparable in magnitude to the aforementioned biases of OLR in association with the neglect of LW scattering by cloud condensates. This indicates that systematic negative biases of global mean LWCF in many GCMs as well as the biases in other LW radiation components may be due to the neglect of subgrid fluctuation of water vapor as well as the neglect of the LW scattering effect by cloud condensates.
In this paper, we attempt to address this issue by incorporating bulk subgrid variations of water vapor between clear and cloudy portions into the LW radiation scheme. From the given grid mean water vapor and cloud fraction in each grid layer, we calculated water vapor content in the clear and cloudy portions separately and used them to compute clear-sky LW radiation, all-sky LW radiation, and LWCF. The responses of the simulated climate to our treatment were examined. We note that treating water vapor variability within each grid cell has one impact on diagnostic calculations of cloud radiative effect and a separate opportunity to change the simulation via radiative heating rates and surface fluxes. Park (2018) examined the former effect. Our study is designed to test the latter as well as the former effects.
The structure of this paper is as follows. Section 2.1 explains how the water vapor content in the clear and cloud portions was computed within individual GCM grid layers. Section 2.2 provides a description on the Seoul National University Atmosphere Model version 0 with a Unified Convection Scheme (SAM0-UNICON), an atmospheric GCM used in our study, and how our new treatment of subgrid variations of water vapor was implemented into the LW radiation scheme. Section 3 compares various aspects of the simulated climates in two global simulations with and without the modification. A summary and conclusions are provided in section 4.

Computation of Subgrid Water Vapor in the Clear and Cloud Portions
From a given set of grid mean water vapor (q v ), grid mean condensates of cloud liquid (q l ), and cloud ice (q i ), and cloud fraction (0 ≤ A ≤ 1) in each grid layer, we need to compute the bulk water vapor content in the clear (q v,clr ) and cloud (q v,cld ) portions. First, we guess that the air in the cloud portion is saturated, so q v,cld ≈ q s . Here, q s is the saturation specific humidity computed as where q s,l and q s,i are the saturation specific humidity over water and ice, respectively, and we set q s ≥q v . The conservation constraint of the grid mean water vapor can be written asq v = A · q s + (1 − A) · q v,clr . From this, we computed q v,clr as which ensures that q v,clr ≤q v since q s ≥q v . To prevent the onset of unreasonable values, we imposed a lower bound (c ·q v ) and an upper bound (q v ) on q v,clr , where the bounding constant, c, was externally specified in the range of 0 ≤ c ≤ 1. The lower bound value of c ·q v is necessary since all cloud fractions in the SAM0-UNICON (e.g., ice stratus fraction, cumulus fraction, and detrained cumulus fraction) except the liquid stratus fraction are defined in a bulk way without an explicit connection to subgrid variations of water vapor (Park et al., 2014), such that q v,clr from equation (2) can be negative when A is large. The sensitivity of simulated LW radiations to c has been discussed later.
The upper bound value ofq v is necessary since q v,clr should be smaller thanq v whenever cloud (the water vapor within the cloud portion is larger than the grid mean water vapor) exists in the grid layer. Then, we derived q v,cld using the conservation constraint of the grid mean water vapor, which ensures that q cld ≥q v sinceq v ≥ q clr . Due to the use of the conservation constraint, q v,cld can be different from q s . To prevent the onset of unreasonable values, we imposed a lower bound (q v ) and an upper bound The computation of q v,clr and q v,cld ensured that q v,clr ≤q v ≤ q v,cld and conserved the grid mean water vapor, clr for all cases of 0 ≤ A ≤ 1. The global-annual mean frequencies for q v,clr to hit the lower bound (c ·q v ) and the upper bound (q v ) in equation (3) are 0.026 and 0.014, respectively. In addition, the global-annual mean frequency for q v,cld to hit the lower bound (q v ) in equation (5) is 0.004. This indicates that the changes made by our implementation are mostly introduced by the physical formulation (e.g., equations (2) and (4)) rather than the limiters (e.g., equations (3) and (5)).

Implementation of q v,clr and q v,cld Into the LW Radiation Scheme
Our study used the SAM0-UNICON  as an atmospheric GCM. The SAM0-UNICON, one of the international GCMs participating in Phase 6 of the Coupled Model Intercomparison Project (CMIP6, Eyring et al., 2016), is based on the Community Atmosphere Model version 5 (CAM5) (Park et al., 2014). However, CAM5's shallow (Park & Bretherton, 2009) and deep convection schemes (Zhang & McFarlane, 1995) are replaced by the UNICON (Park, 2014a(Park, , 2014b, and the treatment of the convective detrainment process is modified . To handle vertical cloud overlap consistently in various physics parameterization schemes, Park (2017Park ( , 2018 developed an interactive vertical overlap parameterization for cumulus and stratus clouds and implemented it into the convection, stratus microphysics, radiation, aerosol wet deposition, and activation schemes of SAM0-UNICON . From now on, SAM0-UNICON with the interactive vertical overlap parameterization will be referred to as the default SAM. In this study, we modified the LW radiation scheme of the default SAM. The bulk subgrid variations in water vapor content between clear and cloudy portions were taken into account while computing clear-sky LW radiation, all-sky LW radiation, and LWCF. In contrast to CAM5 and SAM0-UNICON, which merged cumulus and stratus into a single cloud in the grid layer, the default SAM treats the radiative effects of individual cumulus, stratus, and stratiform snow, each of which has its own fractional area, optical property, and vertical overlap structure. More specifically, by assuming the independent precipitation approximation (IPA, Park, 2014aPark, , 2017Park, , 2018) that convective (stratiform) precipitation flux falling into the stratus (cumulus) does not generate the stratiform (convective) precipitation by accretion, the integrated vertical overlap parameterization constructs a total of 112 subcolumns with multiple cloud types and multiple precipitation types based on the vertical overlap probabilities between the precipitation areas at the top interface of an individual grid layer and the cloud areas at the layer midpoint . The details on constructing a set of stochastic subcolumns with the Monte Carlo Independent Column In reality, the temperature within cumulus may be different from that within stratus in the same grid layer. In addition, compared to stratus, cumulus with strong updraft vertical velocity is more likely to be supersaturated. However, our default SAM cannot handle these differences in water vapors between cumulus and stratus in a fully comprehensive manner. Thus, we simply assumed that cumulus and stratus have the same q v,cld , which also saves computation time. For each subcolumn, cloud fraction in each grid layer was either 0 or 1. To compute clear-sky LW radiation, we used q v,clr in all layers. To compute all-sky LW radiation, however, we used q v,cld in the cloud layer where the cloud fraction was 1 and q v,clr in the clear layer where the cloud fraction was 0. This is similar to the way how satellites observe clear-sky and all-sky radiations. To conserveq v , the water vapor in stratiform snow, falling into the clear portion, was set to q v,clr . In nature, due to the evaporation of precipitation, the water vapor content in the stratiform snow area may be larger than that in the clear portion. However, to be consistent with the formulation described in the previous section that computesq v as the area-weighted averages of cloud and clear-sky water vapors, we assume that the water vapor within the stratiform snow area is the same as water vapor in the cloud (clear) portion, if the stratiform snow is overlapped with the cloud (clear) portion. LWCF was calculated as the difference between all-sky LW radiation and clear-sky LW radiation at TOA. We did not modify the SW radiation scheme because water vapor is nearly transparent to SW radiation. Ideally, the other physics parameterizations as well as the LW radiation scheme should consistently use the same q v,clr and q v,cld computed in equations (2) and (4). For example, the water vapor of environmental air entrained into a convective updraft plume rising into the clear portion should be q v,clr rather thanq v . In addition, the water vapor in the environmental air into which precipitation falls and evaporates should be q v,clr . While the latter is taken into account, the former has not been treated yet. Imposing a full interprocess consistency on the treatment of q v,clr and q v,cld in various physics, parameterizations is a future research subject. The default SAM with the aforementioned treatment of bulk subgrid variations of water vapor between clear and cloudy portions in the LW radiation scheme has been referred to as a modified SAM.  (3)) and then used c = 0 in the AMIP simulation to examine the maximum impact of the bulk subgrid variations of water vapor. The AMIP simulations with the default and modified SAM with c = 0 have been referred to as OLD and NEW, respectively. We analyzed the changes in LW and SW radiative fluxes and other related climates. To reduce the impacts of random noise, we ran four ensemble simulations for each OLD and NEW and then averaged them. Figure 1 shows the variations in annual-global mean clear-sky upward LW radiative flux at TOA (FLUTC), all-sky upward LW radiative flux at TOA (FLUT), and LWCF at TOA as a function of the bounding constant, c (equation 3). The simulation with c = 1 was identical to the default simulation that neglects subgrid variations of water vapor. As c decreased, the FLUTC increased, since the transparency of the clear portion to LW radiation increased (q v,clr <q v ) and a higher amount of LW radiation emitted from the surface reached TOA. In addition, the opacity of the cloudy portion to LW radiation increased (q v,cld >q v ). Interestingly, although there were some fluctuations, the resulting FLUT decreased as c decreased, presumably as a result of nonlinear dependency of LW emissivity to the amount of water vapor and the complex feedback processes involving the changes in cloud fraction. Similar to FLUTC, LWCF (=FLUTC−FLUT) increased almost monotonically as c decreased. The differences of global-annual mean FLUTC, FLUT, and LWCF between c = 0 and c = 1 were approximately 1.3, −0.5, and 1.8 (W m −2 ), respectively, which are large enough to consider an alternative tuning strategy to obtain the global radiation balance at TOA in a GCM. Since the biases in the simulated FLUTC, FLUT, and LWCF can be caused by many uncertainties in other physics parameterizations, we cannot determine what is the correct value of c. Henceforth, we compared the simulations with c = 0 (NEW) and c = 1 (OLD), which will provide information on the maximum impacts of subgrid variations of water vapor between clear and cloudy portions on the LW radiation scheme. Figure 2 shows the annual zonal mean cross sections of grid mean water vapor (q v ) and the difference between cloudy and clear portion water vapor (Δq v = q v,cld − q v,clr ) obtained from NEW simulation. The magnitude of Δq v is comparable to that ofq v , implying that the changes in LW radiation between NEW and OLD may be significant. Sinceq v is large in the tropical lower troposphere, Δq v was also found to be large in the tropical lower troposphere. In the tropical and subtropical regions where cloud fraction A is relatively small, Δq v tends to increase as A decreased. According to equations (2) and (3), if q s −q v ≈q v − q v,clr , Δq v increased as A decreased when A < 0.5, which explains the bimodal shape of Δq v in the tropical and subtropical regions. An additional aspect that should be noted is a nonnegligible Δq v in the Arctic area, which was about 1/3 ofq v . The absorptivity and emissivity ( ) of a single cloud layer is = 1 − exp(− ) and the optical thickness is =k Δz, where is air density, Δz is cloud layer thickness, andk is the absorption cross section. We can further decomposek into the contributions of cloud liquid mass (q l ), cloud ice mass (q i ), and water vapor content (q v ) ask =k l ·q l +k i ·q i +k v ·q v , wherek l ,k i , andk v are the absorption cross sections of in-cloud liquid, ice, and water vapor, respectively. Since Arctic clouds have very lowq l andq i (Park et al., 2014), the use of q v,cld instead ofq v forq v is likely to increase LW radiation emitted from clouds in a nonlinear manner.  , 2009) observations and the differences between the simulations and the observations. The observed FLUTC was high in the regions where underlying surface temperature was high (e.g., tropics) and the atmospheric water vapor content was low (e.g., deserts and subtropical subsidence regions). However, the FLUTC was low where underlying surface temperature was low (e.g., polar regions) and the atmospheric water vapor content was large (e.g., tropical deep convection and midlatitude storm track regions). In almost all regions, NEW simulated a higher FLUTC than OLD, since q v,clr was smaller thanq v (equation (2)), and a higher portion of the LW radiation emitted from surface can reach TOA. The increase in FLUTC was more pronounced in the tropical deep convection and midlatitude storm track regions where upper-and middle-level clouds are abundant (see Figure 2), since according to equation (2), the difference betweenq v and q v,clr increases as A increases. OLD underestimated the observed global mean FLUTC by 4.2 W m −2 (Figure 3c), which reduced to 2.8 W m −2 under NEW (Figure 3d). The rmse as well as the global mean bias of FLUTC reduced, mainly due to the reduction of the negative biases in the midlatitude oceans. However, the positive biases in the tropical regions were further degraded and the midlatitude oceans still had negative biases. These biases are closely associated with the biases in column-integrated water vapor, PREH2O (hatched lines in Figure 3d). The biases in PREH2O could be reduced by reducing the production rate of convective precipitation in the tropics (i.e., weakening the parameterized convective activity) and enhancing the condensation rate or production rate of stratiform precipitation in the midlatitude regions. Figure 4 shows the annual mean all-sky upward LW radiative flux at TOA (FLUT) obtained from the CERES-EBAF observations and the differences between the simulations and the observations. The observed FLUT was high in the regions where underlying surface temperature was high (e.g., tropics) and no overlying clouds exist (e.g., deserts and subtropical subsidence regions) but was small where underlying surface temperature was low (e.g., polar regions) and middle-level clouds (e.g., midlatitude storm track region) or upper-level clouds (tropical deep convection region) exist. Global mean FLUT simulated by NEW was about 0.5 W m −2 lower than that of OLD. In contrast to the ΔFLUTC shown in Figure 3b, ΔFLUT showed complex geographical variations due to the changes in cloud fraction through complex feedback processes (see the hatched lines in Figure 4b) as well as the nonlinear dependency of cloud emissivity on water vapor content. The rmse in NEW was similar to that in OLD. The difference between Figures 4d and 3d indicates that substantial portions of the biases in FLUT were associated with the biases in the simulated clouds. Overall, NEW underestimated the observed global mean FLUT by 0.3 W m −2 , which can be remedied by reducing the amount or increasing cloud top temperature of upper-level clouds in the tropical continental regions. Figure 5 shows the annual mean LWCF at TOA obtained from the CERES-EBAF observations and the differences between the simulations and the observations. The observed LWCF was mostly positive and high in the tropical deep convection and midlatitude storm track regions where upper-and middle-level clouds are abundant (Figure 5a). Globally, NEW simulated stronger LWCF than OLD with the largest increases in the tropical deep convection and midlatitude storm track regions, which seems to be mainly due to positive ΔFLUTC (Figure 3b) and then ΔFLUT (Figure 4b). OLD simulated significantly weak global mean LWCF, approximately 4.3 W m −2 lower than the observations, which was substantially remedied by NEW that had a global mean bias of −2.5 W m −2 . In the high-latitude regions, Δq v = q v,cld − q v,clr is much smaller than in the tropics (Figure 2b) but ΔLWCF is not small. Presumably, this is becauseq v in the high-latitude region is so small that even a small Δq v exerts a relatively large impact on the simulated LWCF. Along with the global mean bias, the rmse is also improved from 7.0 W m −2 in OLD to 6.4 W m −2 in NEW. However, NEW still had positive biases in the tropical continents and far eastern equatorial Pacific Ocean and negative biases in the subtropical and midlatitude regions. It seems to be necessary to weaken convective activity in the tropical continental region in order to reduce the amount of upper-level cirrus clouds and increase cirrus-top temperature. In regions where strong negative anomalies were centered (e.g., midlatitude North and South Pacific, Atlantic, and Indian Oceans), more middle-or upper-level clouds need to be generated, which may be achieved by reducing the critical relative humidity used for computing stratus cloud fraction in the middle-and upper-troposphere above 700 hPa (Park et al., 2014). Figure 6 shows the annual mean SWCF at TOA obtained from the CERES-EBAF observations and the differences between the simulations and the observations. Even though only LW radiation scheme was modified, significant changes occurred in the simulated SW radiation through complex feedback processes. The global-annual mean SWCF decreased by 0.9 (W m −2 ) (i.e., enhanced cooling) from OLD to NEW, mainly due to enhanced SWCF in the eastern subtropical regions in which marine stratocumulus is abundant. Enhanced SWCF in these regions is associated with enhanced low-level cloud fraction (Figure 6b). We found a decrease, from OLD to NEW, in downward LW flux in the lower troposphere at 850-900 hPa above the planetary boundary layer (not shown). It is speculated that LW cooling rate at the stratocumulus top increased due to the reduced downward LW flux from above marine stratocumulus, turbulent moisture transport from the sea surface to overlying stratocumulus deck at the planetary boundary layer top was enhanced, and low-level cloud fraction (CLDLOW) and SWCF increased. With more CLDLOW, grid mean LW cooling rate at the stratocumulus top was further enhanced. It seems that this radiative positive feedback was responsible for the increase in CLDLOW and SWCF from OLD to NEW in the subtropical stratocumulus decks. Several systematic biases still exist in the simulated SWCF in NEW, which could be addressed by decreasing (increasing) CLDLOW in the upstream (downstream) of the subtropical marine stratocumulus decks, enhancing (reducing) convective activity over the tropical oceans (continents), and increasing CLDLOW or reducing sea ice fraction in the southern hemispheric circumpolar region. Overall, global-annual mean net downward SW and SW+LW radiations at TOA from NEW were smaller than those of OLD by 0.9 and 0.44 (W m −2 ), respectively (not shown).

Summary and Conclusions
Most GCMs compute radiation fluxes by assuming that water vapor is uniform within individual grid layers. This assumption leads to underestimation of satellite-observed LWCF since water vapor content in cloudy portion is higher than that in clear portion. To address this issue, we calculated the water vapor content in the clear (q v,clr ) and cloud (q v,cld ) portions separately within individual grid layers, constructed a set of radiative subcolumns by using q v,clr and q v,cld instead of the grid mean water vapor (q v ), and then computed LW radiation using these new subcolumns. The impacts of this treatment on the global climate were examined using the SAM0-UNICON. A set of AMIP simulations were conducted with the modified and default codes, which are referred to as NEW and OLD, respectively.
The differences in global-annual mean clear-sky (FLUTC) and all-sky (FLUT) upward LW fluxes and LWCF at TOA between NEW and OLD were 1.3, −0.5, and 1.8 (W m −2 ), respectively, which are large enough to influence a tuning strategy. The magnitude of Δq v = q v,cld −q v,clr was comparable to that ofq v and was higher in the tropical lower troposphere, extending into the Arctic area. As expected, NEW simulated more FLUTC than OLD, particularly in the tropical deep convection and midlatitude storm track regions where upper-and middle-level clouds are abundant. The global mean bias of FLUTC significantly reduced from −4.2 W m −2 in OLD to −2.8 W m −2 in NEW. The remaining positive (negative) biases of FLUTC in the tropical (midlatitude) regions were closely associated with the biases in the column-integrated water vapor. In contrast to ΔFLUTC, ΔFLUT between NEW and OLD showed complex geographical variations in association with the changes in cloud fraction through complex feedback processes. Overall, NEW underestimated the observed global mean FLUT by 0.3 W m −2 . Similar to FLUTC, NEW simulated stronger LWCF than OLD with the largest increases in the tropical deep convection and midlatitude storm track regions, which resulted in substantial reduction of a long-standing global mean negative bias (from −4.3 to −2.5 W m −2 ) and rmse error of LWCF.
Even though only the LW radiation scheme was modified, global-annual mean SWCF also increased significantly (i.e., more cooling), particularly in the vicinity of the eastern subtropical marine stratocumulus decks. This could be due to the positive radiative feedback triggered by the reduced downward LW flux falling into stratocumulus from above and an associated increase in LW cooling rate at the top of stratocumulus and enhanced turbulent moisture transport from the sea surface to the overlying stratocumulus deck. Global-annual mean net downward SW and SW + LW radiations at TOA from NEW were 0.9 and 0.44 (W m −2 ) lower than those of OLD, respectively.
Our modification of an LW radiation scheme to take subgrid variation of water vapor between clear and cloudy portions into account helped overcome a long-standing negative bias of global mean LWCF in a GCM. This allows us to diagnose the source of the remaining biases (particularly, the biases that persist with the same sign between NEW and OLD) in a more accurate way, reducing the risk of introducing a compensating error. For example, the positive (negative) biases of FLUTC in the tropical (midlatitude) regions are now clearly associated with the biases in the column-integrated water vapor, which could be potentially addressed by reducing (enhancing) the production rate of convective (stratiform) precipitation. The positive (negative) biases of LWCF over the tropical continents (subtropical and midlatitude regions) and similar biases in FLUT with an opposite sign provide a clear indication that our GCM simulates a very strong convection over the tropical continents and too few middle-and upper-level clouds in the midlatitude regions. The easiest way to address these biases associated with convection is to enhance the fractional entrainment rate for convective updraft plumes. A more desirable but more complicated approach is to impose stochastic characteristics on convective updraft plumes, such that a strong top heavy heating in current UNICON (as implied by the positive biases of LWCF over the tropical continents) can be alleviated by allowing multiple convective updraft plumes with different top heights to coexist in the same environment (Shin & Park, 2020). The production rate of convective precipitation can be computed more accurately with a double moment cumulus microphysics scheme interacting with aerosols instead of a current single-moment scheme. The production of stratiform precipitation in the midlatitude regions can be enhanced by decreasing the critical diameter of ice crystals at which stratus ice condensate starts to be converted into snow (Park et al., 2014). Too few middle-and upper-level clouds in the midlatitude regions and associated negative biases of LWCF imply a need to retune or revise cloud macrophysics scheme for ice stratus. Accounting for LW scattering by cloud droplets may also help to reduce the biases in LWCF and FLUT (Jin et al., 2019;Kuo et al., 2017;Schmidt et al., 2006). Our study has focused on the biases in the LW radiation. However, we note that accounting for subgrid variability of cloud condensate within cloud portion can exert significant effects on SW radiation (Barker, 1996;Barker et al., 1996;Cahalan et al., 1994;Pincus et al., 1999). In future, we plan to implement these neglected processes and retune our GCM in a comprehensive manner.