.. _sub_phys_pkg_bulk_force: BULK_FORCE: Bulk Formula Package --------------------------------- author: Stephanie Dutkiewicz Instead of forcing the model with heat and fresh water flux data, this package calculates these fluxes using the changing sea surface temperature. We need to read in some atmospheric data: **air temperature, air humidity, down shortwave radiation, down longwave radiation, precipitation, wind speed**. The current setup also reads in **wind stress**, but this can be changed so that the stresses are calculated from the wind speed. The current setup requires that there is the thermodynamic-seaice package (*pkg/thsice*, also refered below as seaice) is also used. It would be useful though to have it also setup to run with some very simple parametrization of the sea ice. The heat and fresh water fluxes are calculated in *bulkf\_forcing.F* called from *forward\_step.F*. These fluxes are used over open water, fluxes over seaice are recalculated in the sea-ice package. Before the call to *bulkf\_forcing.F* we call *bulkf\_fields\_load.F* to find the current atmospheric conditions. The only other changes to the model code come from the initializing and writing diagnostics of these fluxes. subroutine BULKF_FIELDS_LOAD ++++++++++++++++++++++++++++ Here we find the atmospheric data needed for the bulk formula calculations. These are read in at periodic intervals and values are interpolated to the current time. The data file names come from **data.blk**. The values that can be read in are: air temperature, air humidity, precipitation, down solar radiation, down long wave radiation, zonal and meridional wind speeds, total wind speed, net heat flux, net freshwater forcing, cloud cover, snow fall, zonal and meridional wind stresses, and SST and SSS used for relaxation terms. Not all these files are necessary or used. For instance cloud cover and snow fall are not used in the current bulk formula calculation. If total wind speed is not supplied, wind speed is calculate from the zonal and meridional components. If wind stresses are not read in, then the stresses are calculated from the wind speed. Net heat flux and net freshwater can be read in and used over open ocean instead of the bulk formula calculations (but over seaice the bulkf formula is always used). This is “hardwired” into *bulkf\_forcing* and the “ch” in the variable names suggests that this is “cheating”. SST and SSS need to be read in if there is any relaxation used. subroutine BULKF_FORCING ++++++++++++++++++++++++ In *bulkf\_forcing.F*, we calculate heat and fresh water fluxes (and wind stress, if necessary) for each grid cell. First we determine if the grid cell is open water or seaice and this information is carried by **iceornot**. There is a provision here for a different designation if there is snow cover (but currently this does not make any difference). We then call *bulkf\_formula\_lanl.F* which provides values for: up long wave radiation, latent and sensible heat fluxes, the derivative of these three with respect to surface temperature, wind stress, evaporation. Net long wave radiation is calculated from the combination of the down long wave read in and the up long wave calculated. We then find the albedo of the surface - with a call to *sfc\_albedo* if there is sea-ice (see the seaice package for information on the subroutine). If the grid cell is open ocean the albedo is set as 0.1. Note that this is a parameter that can be used to tune the results. The net short wave radiation is then the down shortwave radiation minus the amount reflected. If the wind stress needed to be calculated in *bulkf\_formula\_lanl.F*, it was calculated to grid cell center points, so in *bulkf\_forcing.F* we regrid to **u** and **v** points. We let the model know if it has read in stresses or calculated stresses by the switch **readwindstress** which is can be set in data.blk, and defaults to **.TRUE.**. We then calculate **Qnet** and **EmPmR** that will be used as the fluxes over the open ocean. There is a provision for using runoff. If we are “cheating” and using observed fluxes over the open ocean, then there is a provision here to use read in **Qnet** and **EmPmR**. The final call is to calculate averages of the terms found in this subroutine. subroutine BULKF_FORMULA_LANL +++++++++++++++++++++++++++++ This is the main program of the package where the heat fluxes and freshwater fluxes over ice and open water are calculated. Note that this subroutine is also called from the seaice package during the iterations to find the ice surface temperature. Latent heat (:math:`L`) used in this subroutine depends on the state of the surface: vaporization for open water, fusion and vaporization for ice surfaces. Air temperature is converted from Celsius to Kelvin. If there is no wind speed (:math:`u_s`) given, then the wind speed is calculated from the zonal and meridional components. We calculate the virtual temperature: .. math:: T_o = T_{air} (1+\gamma q_{air}) where :math:`T_{air}` is the air temperature at :math:`h_T`, :math:`q_{air}` is humidity at :math:`h_q` and :math:`\gamma` is a constant. The saturated vapor pressure is calculate (QQ ref): .. math:: q_{sat} = \frac{a}{p_o} e^{L (b-\frac{c}{T_{srf}})} where :math:`a,b,c` are constants, :math:`T_{srf}` is surface temperature and :math:`p_o` is the surface pressure. The two values crucial for the bulk formula calculations are the difference between air at sea surface and sea surface temperature: .. math:: \Delta T = T_{air} - T_{srf} +\alpha h_T where :math:`\alpha` is adiabatic lapse rate and :math:`h_T` is the height where the air temperature was taken; and the difference between the air humidity and the saturated humidity .. math:: \Delta q = q_{air} - q_{sat}. We then calculate the turbulent exchange coefficients following Bryan et al (1996) and the numerical scheme of Hunke and Lipscombe (1998). We estimate initial values for the exchange coefficients, :math:`c_u`, :math:`c_T` and :math:`c_q` as .. math:: \frac{\kappa}{ln(z_{ref}/z_{rou})} where :math:`\kappa` is the Von Karman constant, :math:`z_{ref}` is a reference height and :math:`z_{rou}` is a roughness length scale which could be a function of type of surface, but is here set as a constant. Turbulent scales are: .. math:: \begin{aligned} u^* & = & c_u u_s \nonumber\\ T^* & = & c_T \Delta T \nonumber\\ q^* & = & c_q \Delta q \nonumber\end{aligned} We find the “integrated flux profile” for momentum and stability if there are stable QQ conditions (:math:`\Upsilon>0`) : .. math:: \psi_m = \psi_s = -5 \Upsilon and for unstable QQ conditions (:math:`\Upsilon<0`): .. math:: \begin{aligned} \psi_m & = & 2 ln(0.5(1+\chi)) + ln(0.5(1+\chi^2)) - 2 \tan^{-1} \chi + \pi/2 \nonumber \\ \psi_s & = & 2 ln(0.5(1+\chi^2)) \nonumber\end{aligned} where .. math:: \Upsilon = \frac{\kappa g z_{ref}}{u^{*2}} (\frac{T^*}{T_o} + \frac{q^*}{1/\gamma + q_a}) and :math:`\chi=(1-16\Upsilon)^{1/2}`. The coefficients are updated through 5 iterations as: .. math:: \begin{aligned} c_u & = & \frac {\hat{c_u}}{1+\hat{c_u}(\lambda - \psi_m)/\kappa} \nonumber \\ c_T & = & \frac {\hat{c_T}}{1+\hat{c_T}(\lambda - \psi_s)/\kappa} \nonumber \\ c_q & = & c'_T\end{aligned} where :math:`\lambda =ln(h_T/z_{ref})`. We can then find the bulk formula heat fluxes: Sensible heat flux: .. math:: Q_s=\rho_{air} c_{p_{air}} u_s c_u c_T \Delta T Latent heat flux: .. math:: Q_l=\rho_{air} L u_s c_u c_q \Delta q Up long wave radiation .. math:: Q_{lw}^{up}=\epsilon \sigma T_{srf}^4 where :math:`\epsilon` is emissivity (which can be different for open ocean, ice and snow), :math:`\sigma` is Stefan-Boltzman constant. We calculate the derivatives of the three above functions with respect to surface temperature .. math:: \begin{aligned} \frac{dQ_s}{d_T} & = & \rho_{air} c_{p_{air}} u_s c_u c_T \nonumber \\ \frac{dQ_l}{d_T} & = & \frac{\rho_{air} L^2 u_s c_u c_q c}{T_{srf}^2} \nonumber \\ \frac{dQ_{]lw}^{up}}{d_T} & = & 4 \epsilon \sigma t_{srf}^3 \nonumber\end{aligned} And total derivative :math:`\frac{dQ_o}{dT}= \frac{dQ_s}{dT} + \frac{dQ_l}{dT} + \frac{dQ_{lw}^{up}}{dT}`. If we do not read in the wind stress, it is calculated here. Initializing subroutines ++++++++++++++++++++++++ :code:`bulkf_init.F`: Set bulkf variables to zero. :code:`bulkf_readparms.F`: Reads **data.blk** Diagnostic subroutines ++++++++++++++++++++++ :code:`bulkf_ave.F`: Keeps track of means of the bulkf variables :code:`bulkf_diags.F`: Finds averages and writes out diagnostics Common Blocks +++++++++++++ :code:`BULKF.h`: BULKF Variables, data file names, and logicals **readwindstress** and **readsurface** :code:`BULKF_DIAGS.h`: matrices for diagnostics: averages of fields from *bulkf_diags.F* :code:`BULKF_ICE_CONSTANTS.h`: all the parameters needed by the ice model and in the bulkf formula calculations. Input file DATA.ICE +++++++++++++++++++ We read in the file names of atmospheric data used in the bulk formula calculations. Here we can also set the logicals: **readwindstress** if we read in the wind stress rather than calculate it from the wind speed; and **readsurface** to read in the surface temperature and salinity if these will be used as part of a relaxing term. Important Notes +++++++++++++++ #. heat fluxes have different signs in the ocean and ice models. #. **StartIceModel** must be changed in **data.ice**: 1 (if starting from no ice), 0 (if using pickup.ic file). References ++++++++++ Bryan F.O., B.G Kauffman, W.G. Large, P.R. Gent, 1996: The NCAR CSM flux coupler. Technical note TN-425+STR, NCAR. Hunke, E.C and W.H. Lipscomb, circa 2001: CICE: the Los Alamos Sea Ice Model Documentation and Software User’s Manual. LACC-98-16v.2. (note: this documentation is no longer available as CICE has progressed to a very different version 3) Experiments and tutorials that use bulk\_force ++++++++++++++++++++++++++++++++++++++++++++++ - Global ocean experiment in global\_ocean.cs32x15 verification directory, input from input.thsice directory.