stratification.F90 6.8 KB
 dumoda01 committed Jan 19, 2011 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 !$Id: stratification.F90,v 1.7 2005-07-18 08:54:33 lars Exp$ #include"cppdefs.h" !----------------------------------------------------------------------- !BOP ! ! !ROUTINE: Calculation of the stratification\label{sec:stratification} ! ! !INTERFACE: subroutine stratification(nlev,buoy_method,dt,cnpar,nub,gamB) ! ! !DESCRIPTION: ! This routine computes the mean potential density, $\mean{\rho}$, the mean ! potential buoyancy, $B$, defined in \eq{DefBuoyancy}, and the mean buoyancy ! frequency, ! ! \label{DefBuoyancyFrequency} ! N^2 = - \dfrac{g}{\rho_0} \partder{\rho}{z} = \partder{B}{z} ! \comma ! ! which is based on potential density or buoyancy such that for $N^2=0$, the entropy ! is constant in the whole water column and mixing does not work against buoyancy ! forces. If GOTM used as a turbulence library in your own three-dimensional model, ! you have to insure that the $N^2$ computed by you, and passed to the turbulence ! routines in GOTM, is consistent with the concept of potential density and your ! equation of state. ! ! The mean potential density is evaluated from the equation of state, \eq{DefEOS}, ! according to ! ! \label{DefPotentialDensity} ! \mean{\rho} = \hat{\rho} (\Theta,S,P_R) ! \comma ! ! where $\Theta$ denotes the mean potential temperature, $S$ the mean salinity ! and $P_R$ the mean reference pressure. The buoyancy frequency defined in ! \eq{DefBuoyancyFrequency} can be decomposed into contributions due to ! potential temperature and salinity stratification, ! ! \label{NDecompostionA} ! N^2 = N_\Theta^2 + N_S^2 ! \comma ! ! where we introduced the quantities ! ! \label{NNT} ! N_\Theta^2 = - \dfrac{g}{\rho_0} \partder{\rho}{z} \Big|_{S} ! = g \alpha(\Theta,S,P_R) \partder{\Theta}{z} ! \comma ! ! with the thermal expansion coefficient defined in \eq{eosAlpha}, and ! ! \label{NNS} ! N_S^2 = - \dfrac{g}{\rho_0} \partder{\rho}{z} \Big|_{\Theta} ! = - g \beta(\Theta,S,P_R) \partder{S}{z} ! \comma ! ! with the saline contraction coefficient defined in \eq{eosBeta}. It is important ! to note that in the actual code the reference pressure, $P_R$, has been replaced by ! the (approximate) hydrostatic pressure. Only if this dependence is replaced by ! the constant reference pressure at the surface in the equation of state, ! see \sect{sec:eqstate}, the model is truely based on potential temperature and density. ! Otherwise, the model is based on \emph{in-situ} quantities. ! ! Alternatively to the procedure outlined above, depending on the values of the ! parameter {\tt buoy\_method}, the buoyancy may be calculated by means of the ! transport equation \eq{bEq}. This equation then replaces the computation of $\Theta$ ! and $S$ and is only recommended for idealized studies. ! ! !USES: use meanflow, only: h,S,T,buoy,rho use meanflow, only: NN,NNT,NNS use meanflow, only: gravity,rho_0 use eqstate, only: eqstate1 IMPLICIT NONE ! ! !INPUT PARAMETERS: ! number of vertical layers integer, intent(in) :: nlev ! method to compute buoyancy integer, intent(in) :: buoy_method ! time step (s) REALTYPE, intent(in) :: dt ! numerical "implicitness" parameter REALTYPE, intent(in) :: cnpar ! diffusivity of buoyancy (m^2/s) REALTYPE, intent(in) :: nub(0:nlev) ! non-local buoyancy flux (m^2/s^3) REALTYPE, intent(in) :: gamb(0:nlev) ! ! !OUTPUT PARAMETERS: ! ! !REVISION HISTORY: ! Original author(s): Karsten Bolding & Hans Burchard ! ! $Log: stratification.F90,v$ ! Revision 1.7 2005-07-18 08:54:33 lars ! changed docu for html compliance ! ! Revision 1.6 2005/06/27 13:44:07 kbk ! modified + removed traling blanks ! ! Revision 1.5 2003/03/28 09:20:35 kbk ! added new copyright to files ! ! Revision 1.4 2003/03/28 08:56:56 kbk ! removed tabs ! ! Revision 1.3 2003/03/10 08:50:07 gotm ! Improved documentation and cleaned up code ! ! Revision 1.2 2001/11/18 11:50:37 gotm ! Cleaned ! ! Revision 1.1.1.1 2001/02/12 15:55:57 gotm ! initial import into CVS ! !EOP ! ! !LOCAL VARIABLES: integer :: i REALTYPE :: buoyp,buoym,Sface,Tface REALTYPE :: zCenter,zFace,dz REALTYPE :: pFace,pCenter integer, parameter :: USEEQSTATE=1 REALTYPE, parameter :: oneTenth=_ONE_/10. ! !----------------------------------------------------------------------- !BOC if (buoy_method .EQ. USEEQSTATE) then ! initialize parameters for uppermost grid box zFace = _ZERO_ zCenter = 0.5*h(nlev) pCenter = oneTenth*zCenter buoy(nlev) = eqstate1(S(nlev),T(nlev),pCenter,gravity,rho_0) rho(nlev) = rho_0 - rho_0/gravity*buoy(nlev) do i=nlev-1,1,-1 ! compute distances between centers dz = 0.5*(h(i)+h(i+1)) ! compute local depths zFace = zFace + h(i+1) zCenter = zCenter + dz ! compute approximate local pressure pFace = oneTenth*zFace pCenter = oneTenth*zCenter ! interpolate T and S from centers to faces Sface = ( S(i+1)*h(i) + S(i)*h(i+1) ) / ( h(i+1) + h(i) ) Tface = ( T(i+1)*h(i) + T(i)*h(i+1) ) / ( h(i+1) + h(i) ) ! T contribution to buoyancy frequency buoyp = eqstate1(Sface,T(i+1),pFace,gravity,rho_0) buoym = eqstate1(Sface,T(i ),pFace,gravity,rho_0) NNT(i) = (buoyp-buoym)/dz ! S contribution to buoyancy frequency buoyp = eqstate1(S(i+1),Tface,pFace,gravity,rho_0) buoym = eqstate1(S(i ),Tface,pFace,gravity,rho_0) NNS(i) = (buoyp-buoym)/dz ! total buoyancy frequency is the sum NN(i) = NNT(i) + NNS(i) ! compute buoyancy and density buoy(i) = eqstate1(S(i),T(i),pCenter,gravity,rho_0) rho(i) = rho_0 - rho_0/gravity*buoy(i) end do else ! for some idealized cases, compute buoyancy from ! prognostic equation call buoyancy(nlev,dt,cnpar,nub,gamb) ! compute density and buoyancy frequency rho(nlev) = rho_0 - rho_0/gravity*buoy(nlev) do i=nlev-1,1,-1 dz = 0.5*(h(i)+h(i+1)) NN(i) = (buoy(i+1)-buoy(i))/dz rho(i) = rho_0 - rho_0/gravity*buoy(i) end do end if ! update boundary values NN(0) = _ZERO_ NN(nlev) = _ZERO_ return end subroutine stratification !EOC !----------------------------------------------------------------------- ! Copyright by the GOTM-team under the GNU Public License - www.gnu.org !-----------------------------------------------------------------------