revised implementation of equation of state

src/gotm/diagnostics.F90

@@ -89,7 +89,7 @@ subroutine do_diagnostics(nlev)
 !
 ! !USES:
    use airsea_driver,only: sst,tx,ty
-   use density,      only: rho0, rho, calculate_density
+   use density,      only: rho0
    use meanflow,     only: gravity,cp,drag
    use meanflow,     only: h,u,v,s,t,NN,SS,buoy,rad
    use turbulence,   only: turb_method

src/gotm/gotm.F90

@@ -44,7 +44,7 @@ module gotm
    use diagnostics
    use settings
 
-   use density, only: init_density,calculate_density,do_density
+   use density, only: init_density,do_density
    use density, only: density_method,T0,S0,p0,rho0,alpha0,beta0
    use density, only: rho, rho_p
    use meanflow
@@ -700,7 +700,7 @@ subroutine initialize_gotm()
    ! Call stratification to make sure density has sensible value.
    ! This is needed to ensure the initial density is saved correctly, and also for FABM.
    call shear(nlev,cnpar)
-   call do_density(nlev,S,T,-zi)
+   call do_density(nlev,S,T,-z,-zi)
    buoy(1:) = -gravity*(rho_p(1:)-rho0)/rho0
    call stratification(nlev)
 
@@ -909,7 +909,7 @@ subroutine integrate_gotm()
 !GSW
 !  update shear and stratification
    call shear(nlev,cnpar)
-   call do_density(nlev,S,T,-zi)
+   call do_density(nlev,S,T,-z,-zi)
    buoy(1:nlev) = -gravity*(rho_p(1:nlev)-rho0)/rho0
    call stratification(nlev)
 

src/gotm/register_all_variables.F90

@@ -391,13 +391,9 @@ subroutine register_density_variables(nlev)
    LEVEL2 'register_density_variables()'
    call fm%register('rho', 'kg/m3', 'density (in-situ)', standard_name='??', dimensions=(/id_dim_z/), data1d=rho(1:nlev),category='density')
    call fm%register('rho_p', 'kg/m3', 'density (potential)', standard_name='??', dimensions=(/id_dim_z/), data1d=rho_p(1:nlev),category='density')
-#if 0
-   call fm%register('alpha', '1/K', 'thermal expansion coefficient', standard_name='??', dimensions=(/id_dim_z/), data1d=alpha(1:nlev),category='density',output_level=output_level_debug)
-   call fm%register('beta', 'kg/g', 'saline contraction coefficient', standard_name='??', dimensions=(/id_dim_z/), data1d=beta(1:nlev),category='density',output_level=output_level_debug)
-#else
-   call fm%register('alpha', '1/K', 'thermal expansion coefficient', standard_name='??', dimensions=(/id_dim_z/), data1d=alpha(1:nlev),category='density')
-   call fm%register('beta', 'kg/g', 'saline contraction coefficient', standard_name='??', dimensions=(/id_dim_z/), data1d=beta(1:nlev),category='density')
-#endif
+   call fm%register('alpha', '1/K', 'thermal expansion coefficient', standard_name='??', dimensions=(/id_dim_zi/), data1d=alpha(0:nlev),category='density',output_level=output_level_debug)
+   call fm%register('beta', 'kg/g', 'saline contraction coefficient', standard_name='??', dimensions=(/id_dim_zi/), data1d=beta(0:nlev),category='density',output_level=output_level_debug)
+
    end subroutine register_density_variables
 !EOC
 

src/meanflow/convectiveadjustment.F90

@@ -28,7 +28,7 @@ subroutine convectiveadjustment(nlev,num,nuh,const_num,const_nuh, &
 ! the eddy viscosity $\nu_t$ and the eddy diffusivity $\nu'_t$, respectively.
 !
 ! !USES:
-   use density, only: calculate_density
+   use density, only:  get_rho_p
    use meanflow, only: h,t,s,buoy,NN
 !
    IMPLICIT NONE
@@ -67,8 +67,8 @@ subroutine convectiveadjustment(nlev,num,nuh,const_num,const_nuh, &
       hint=h(nlev)
       i=nlev
 111   i=i-1
-      buoyupp=calculate_density(Sint,Tint,hint,g)
-      buoylow=calculate_density(S(i),T(i),hint,g)
+      buoyupp =-g*( get_rho_p(Sint,Tint,hint) - rho_0 )/rho_0
+      buoylow =-g*( get_rho_p(S(i),T(i),hint) - rho_0 )/rho_0
       if (buoyupp.lt.buoylow) then     ! instable stratification
          NN(i)=0.
          Tint=(Tint*hint+T(i)*h(i))/(hint+h(i))
@@ -84,7 +84,7 @@ subroutine convectiveadjustment(nlev,num,nuh,const_num,const_nuh, &
       do i=ii,nlev
          T(i)=Tint
          S(i)=Sint
-         buoy(i)=calculate_density(Sint,Tint,zero,g)
+         buoy(i)=-g*( get_rho_p(Sint,Tint,zero) - rho_0 )/rho_0
       end do
    else   ! if (buoy_method.eq.2)
       buoyint=buoy(nlev)

src/meanflow/internal_pressure.F90

@@ -152,7 +152,7 @@ subroutine internal_pressure(nlev)
 ! into the plume.
 !
 ! !USES:
-   use density,       only: calculate_density
+   use density,       only: get_rho,rho0
    use meanflow,      only: T,S
    use meanflow,      only: gravity,h
    use meanflow,      only: buoy
@@ -196,18 +196,18 @@ subroutine internal_pressure(nlev)
          z=z+0.5*h(i)
 
 !        buoyancy gradient in x direction
-         dSS=dx*dsdx_input%data(i)
-         dTT=dx*dtdx_input%data(i)
-         Bl=calculate_density(S(i),T(i),z,gravity)
-         Br=calculate_density(S(i)+dSS,T(i)+dTT,z,gravity)
-         dxB(i)=(Br-Bl)/dx
+         dSS    = dx*dsdx_input%data(i)
+         dTT    = dx*dtdx_input%data(i)
+         Bl     = -gravity*(get_rho(S(i)    ,T(i)    ,z) - rho0)/rho0
+         Br     = -gravity*(get_rho(S(i)+dSS,T(i)+dTT,z) - rho0)/rho0
+         dxB(i) = (Br-Bl)/dx
 
 !        buoyancy gradient in y direction
-         dSS=dy*dsdy_input%data(i)
-         dTT=dy*dtdy_input%data(i)
-         Bl=calculate_density(S(i),T(i),z,gravity)
-         Br=calculate_density(S(i)+dSS,T(i)+dTT,z,gravity)
-         dyB(i)=(Br-Bl)/dy
+         dSS    = dy*dsdy_input%data(i)
+         dTT    = dy*dtdy_input%data(i)
+         Bl     = -gravity*(get_rho(S(i)     ,T(i)   ,z) - rho0)/rho0
+         Br     = -gravity*(get_rho(S(i)+dSS,T(i)+dTT,z) - rho0)/rho0
+         dyB(i) = (Br-Bl)/dy
 
          z=z+0.5*h(i)
       end do

src/observations/const_NNS.F90

@@ -12,7 +12,7 @@ subroutine const_NNS(nlev,z,zi,S_top,T_const,NN,gravity,S)
 ! {\tt v1}
 !
 ! !USES:
-   use density, only: get_beta, beta
+   use density, only: get_beta
    IMPLICIT NONE
 !
 ! !INPUT PARAMETERS:
@@ -29,17 +29,28 @@ subroutine const_NNS(nlev,z,zi,S_top,T_const,NN,gravity,S)
 !  Original author(s): Lars Umlauf
 !
 ! !LOCAL VARIABLES:
-   integer :: n
+   integer                               :: i
+    REALTYPE                             :: lbeta   
 !EOP
 !-----------------------------------------------------------------------
 !BOC
-   S(nlev) = S_top ! must be done here
-   call get_beta(nlev,gravity,T_const,NN,z,zi,S)
+    S(nlev) = S_top ! must be done here
 
-   do n=nlev-1,1,-1
-      S(n) = S(n+1) + _ONE_/(gravity*beta(n))*NN*(z(n+1)-z(n)) 
+    do i=nlev-1,1,-1
+      ! estimate interface beta based on S above the interface
+      lbeta = get_beta(S(i+1),T_const,-zi(i))
+
+      ! use this to estimate S below the interface
+      S(i)  = S(i+1) + _ONE_/(gravity*lbeta)*NN*(z(i+1)-z(i))
+
+      ! compute improved interface beta
+      lbeta  = get_beta(0.5*(S(i+1)+S(i)),T_const,-zi(i))
+
+      ! compute final salinity profile
+      S(i)   = S(i+1) + _ONE_/(gravity*lbeta)*NN*(z(i+1)-z(i))
    end do
-   end subroutine const_NNS
+
+   end subroutine const_NNS   
 !EOC
 
 !-----------------------------------------------------------------------

src/observations/const_NNT.F90

@@ -12,7 +12,7 @@ subroutine const_NNT(nlev,z,zi,T_top,S_const,NN,gravity,T)
 ! {\tt v1}
 !
 ! !USES:
-   use density, only: get_alpha, alpha
+   use density, only: get_alpha
    IMPLICIT NONE
 !
 ! !INPUT PARAMETERS:
@@ -29,16 +29,28 @@ subroutine const_NNT(nlev,z,zi,T_top,S_const,NN,gravity,T)
 !  Original author(s): Lars Umlauf
 !
 ! !LOCAL VARIABLES:
-   integer :: n
+   integer                             :: i
+   REALTYPE                           :: lalpha
+
 !EOP
 !-----------------------------------------------------------------------
 !BOC
-   T(nlev) = T_top ! must be done here
-   call get_alpha(nlev,gravity,S_const,NN,z,zi,T)
+    T(nlev) = T_top ! must be done here
 
-   do n=nlev-1,1,-1
-      T(n) = T(n+1) - _ONE_/(gravity*alpha(n))*NN*(z(n+1)-z(n))
+    do i=nlev-1,1,-1
+      ! estimate interface alpha based on T above the interface
+      lalpha = get_alpha(S_const,T(i+1),-zi(i))
+
+      ! use this to estimate T below the interface
+      T(i)   = T(i+1) - _ONE_/(gravity*lalpha)*NN*(z(i+1)-z(i))
+
+      ! compute improved interface alpha
+      lalpha = get_alpha(S_const,0.5*(T(i+1)+T(i)),-zi(i))
+
+      ! compute final temperature profile
+      T(i)   = T(i+1) - _ONE_/(gravity*lalpha)*NN*(z(i+1)-z(i))
    end do
+
    end subroutine const_NNT
 !EOC