Changeset 388 for palm/trunk/SOURCE/eqn_state_seawater.f90

Timestamp:

Sep 23, 2009 9:40:33 AM (15 years ago)

Author:

raasch

Message:

in-situ AND potential density are calculated and used in the ocean version

File:

• palm/trunk/SOURCE/eqn_state_seawater.f90 (modified) (6 diffs)

palm/trunk/SOURCE/eqn_state_seawater.f90

@@ -4 +4 @@
 ! Actual revisions:
 ! -----------------
-! First constant in array den also defined as type double
+! Potential density is additionally calculated in eqn_state_seawater,
+! first constant in array den also defined as type double
 !
 ! Former revisions:
@@ -69 +70 @@
        INTEGER ::  i, j, k
 
-       REAL ::  p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
+       REAL ::  pden, pnom, p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
 
        DO  i = nxl, nxr
@@ -76 +77 @@
 !
 !--             Pressure is needed in dbar
-!                p1 = hyp(0) * 1E-4
-!                p1 = 0.0
                 p1 = hyp(k) * 1E-4
                 p2 = p1 * p1
@@ -93 +92 @@
                 sa2  = sa1 * sa1
 
-                rho(k,j,i) =                                                   &
-                        ( nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
-                          nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
-                          nom(7)*sa2       + nom(8)*p1      + nom(9)*p1*pt2  + &
-                          nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2   &
-                        ) /                                                    &
-                        ( den(1)           + den(2)*pt1     + den(3)*pt2     + &
-                          den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
-                          den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
-                          den(10)*sa15*pt2 + den(11)*p1     + den(12)*p2*pt3 + &
-                          den(13)*p3*pt1                                       &
-                        )
+                pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
+                       nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
+                       nom(7)*sa2
+
+                pden = den(1)           + den(2)*pt1     + den(3)*pt2     + &
+                       den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
+                       den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
+                       den(10)*sa15*pt2
+
+!
+!--             Potential density (without pressure terms)
+                prho(k,j,i) = pnom / pden
+
+                pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  + &
+                       nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
+
+                pden = pden +             den(11)*p1     + den(12)*p2*pt3 + &
+                       den(13)*p3*pt1
+
+!
+!--             In-situ density
+                rho(k,j,i) = pnom / pden
 
              ENDDO
 !
 !--          Neumann conditions are assumed at bottom and top boundary
-             rho(nzt+1,j,i)            = rho(nzt,j,i)
-             rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
+             prho(nzt+1,j,i)            = prho(nzt,j,i)
+             prho(nzb_s_inner(j,i),j,i) = prho(nzb_s_inner(j,i)+1,j,i)
+             rho(nzt+1,j,i)             = rho(nzt,j,i)
+             rho(nzb_s_inner(j,i),j,i)  = rho(nzb_s_inner(j,i)+1,j,i)
+
           ENDDO
        ENDDO
@@ -129 +141 @@
        INTEGER ::  i, j, k
 
-       REAL ::  p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
+       REAL ::  pden, pnom, p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
 
        DO  k = nzb_s_inner(j,i)+1, nzt
 !
 !--       Pressure is needed in dbar
-!          p1 = hyp(0) * 1E-4
-!          p1 = 0.0
           p1 = hyp(k) * 1E-4
           p2 = p1 * p1
@@ -151 +161 @@
           sa2  = sa1 * sa1
 
-          rho(k,j,i) = ( nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
-                         nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
-                         nom(7)*sa2       + nom(8)*p1      + nom(9)*p1*pt2  + &
-                         nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2   &
-                       ) /                                                    &
-                       ( den(1)           + den(2)*pt1     + den(3)*pt2     + &
-                         den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
-                         den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
-                         den(10)*sa15*pt2 + den(11)*p1     + den(12)*p2*pt3 + &
-                         den(13)*p3*pt1                                       &
-                       )
+          pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
+                 nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
+                 nom(7)*sa2
+
+          pden = den(1)           + den(2)*pt1     + den(3)*pt2     + &
+                 den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
+                 den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
+                 den(10)*sa15*pt2
+
+!
+!--       Potential density (without pressure terms)
+          prho(k,j,i) = pnom / pden
+
+          pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  + &
+                 nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
+          pden = pden +             den(11)*p1     + den(12)*p2*pt3 + &
+                 den(13)*p3*pt1
+
+!
+!--       In-situ density
+          rho(k,j,i) = pnom / pden
+
+
        ENDDO
+
 !
 !--    Neumann conditions are assumed at bottom and top boundary
-       rho(nzt+1,j,i)            = rho(nzt,j,i)
-       rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
+       prho(nzt+1,j,i)            = prho(nzt,j,i)
+       prho(nzb_s_inner(j,i),j,i) = prho(nzb_s_inner(j,i)+1,j,i)
+       rho(nzt+1,j,i)             = rho(nzt,j,i)
+       rho(nzb_s_inner(j,i),j,i)  = rho(nzb_s_inner(j,i)+1,j,i)
 
     END SUBROUTINE eqn_state_seawater_ij