Home Contact Us Site Map  
 
       
    next up previous contents
Next: 6.9 Generic Advection/Diffusion Up: 6. Physical Parameterization and Previous: 6.7 Sea Ice Package:   Contents


6.8 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 ($ 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 ($ u_s$) given, then the wind speed is calculated from the zonal and meridional components.

We calculate the virtual temperature:

$\displaystyle T_o = T_{air} (1+\gamma q_{air}) $

where $ T_{air}$ is the air temperature at $ h_T$, $ q_{air}$ is humidity at $ h_q$ and $ \gamma$ is a constant.

The saturated vapor pressure is calculate (QQ ref):

$\displaystyle q_{sat} = \frac{a}{p_o} e^{L (b-\frac{c}{T_{srf}})} $

where $ a,b,c$ are constants, $ T_{srf}$ is surface temperature and $ 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:

$\displaystyle \Delta T = T_{air} - T_{srf} +\alpha h_T $

where $ \alpha $ is adiabatic lapse rate and $ h_T$ is the height where the air temperature was taken; and the difference between the air humidity and the saturated humidity

$\displaystyle \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, $ c_u$, $ c_T$ and $ c_q$ as

$\displaystyle \frac{\kappa}{ln(z_{ref}/z_{rou})} $

where $ \kappa$ is the Von Karman constant, $ z_{ref}$ is a reference height and $ 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:
$\displaystyle u^*$ $\displaystyle =$ $\displaystyle c_u u_s$  
$\displaystyle T^*$ $\displaystyle =$ $\displaystyle c_T \Delta T$  
$\displaystyle q^*$ $\displaystyle =$ $\displaystyle c_q \Delta q$  

We find the "integrated flux profile" for momentum and stability if there are stable QQ conditions ( $ \Upsilon>0$) :

$\displaystyle \psi_m = \psi_s = -5 \Upsilon $

and for unstable QQ conditions ( $ \Upsilon<0$):
$\displaystyle \psi_m$ $\displaystyle =$ $\displaystyle 2 ln(0.5(1+\chi)) + ln(0.5(1+\chi^2)) - 2 \tan^{-1} \chi + \pi/2$  
$\displaystyle \psi_s$ $\displaystyle =$ $\displaystyle 2 ln(0.5(1+\chi^2))$  

where

$\displaystyle \Upsilon = \frac{\kappa g z_{ref}}{u^{*2}} (\frac{T^*}{T_o} + \frac{q^*}{1/\gamma + q_a}) $

and $ \chi=(1-16\Upsilon)^{1/2}$.

The coefficients are updated through 5 iterations as:

$\displaystyle c_u$ $\displaystyle =$ $\displaystyle \frac {\hat{c_u}}{1+\hat{c_u}(\lambda - \psi_m)/\kappa}$  
$\displaystyle c_T$ $\displaystyle =$ $\displaystyle \frac {\hat{c_T}}{1+\hat{c_T}(\lambda - \psi_s)/\kappa}$  
$\displaystyle c_q$ $\displaystyle =$ $\displaystyle c'_T$ (6.23)

where $ \lambda =ln(h_T/z_{ref})$.

We can then find the bulk formula heat fluxes:


Sensible heat flux:

$\displaystyle Q_s=\rho_{air} c_{p_{air}} u_s c_u c_T \Delta T $


Latent heat flux:

$\displaystyle Q_l=\rho_{air} L u_s c_u c_q \Delta q $


Up long wave radiation

$\displaystyle Q_{lw}^{up}=\epsilon \sigma T_{srf}^4 $

where $ \epsilon$ is emissivity (which can be different for open ocean, ice and snow), $ \sigma $ is Stefan-Boltzman constant.

We calculate the derivatives of the three above functions with respect to surface temperature

$\displaystyle \frac{dQ_s}{d_T}$ $\displaystyle =$ $\displaystyle \rho_{air} c_{p_{air}} u_s c_u c_T$  
$\displaystyle \frac{dQ_l}{d_T}$ $\displaystyle =$ $\displaystyle \frac{\rho_{air} L^2 u_s c_u c_q c}{T_{srf}^2}$  
$\displaystyle \frac{dQ_{]lw}^{up}}{d_T}$ $\displaystyle =$ $\displaystyle 4 \epsilon \sigma t_{srf}^3$  

And total derivative $ \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

bulkf_init.F: Set bulkf variables to zero.

bulkf_readparms.F: Reads data.blk




Diagnostic subroutines

bulkf_ave.F: Keeps track of means of the bulkf variables

bulkf_diags.F: Finds averages and writes out diagnostics




Common Blocks

BULKF.h: BULKF Variables, data file names, and logicals readwindstress and readsurface

BULKF_DIAGS.h: matrices for diagnostics: averages of fields from bulkf_diags.F

BULKF_ICE_CONSTANTS.h: all the parameters need 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

1) heat fluxes have different signs in the ocean and ice models.

2) 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)

next up previous contents
Next: 6.9 Generic Advection/Diffusion Up: 6. Physical Parameterization and Previous: 6.7 Sea Ice Package:   Contents
mitgcm-support@dev.mitgcm.org
Copyright © 2002 Massachusetts Institute of Technology