Ignore:
Timestamp:
Sep 23, 2009 9:40:33 AM (12 years ago)
Author:
raasch
Message:

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/eqn_state_seawater.f90

    r336 r388  
    44! Actual revisions:
    55! -----------------
    6 ! First constant in array den also defined as type double
     6! Potential density is additionally calculated in eqn_state_seawater,
     7! first constant in array den also defined as type double
    78!
    89! Former revisions:
     
    6970       INTEGER ::  i, j, k
    7071
    71        REAL ::  p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
     72       REAL ::  pden, pnom, p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
    7273
    7374       DO  i = nxl, nxr
     
    7677!
    7778!--             Pressure is needed in dbar
    78 !                p1 = hyp(0) * 1E-4
    79 !                p1 = 0.0
    8079                p1 = hyp(k) * 1E-4
    8180                p2 = p1 * p1
     
    9392                sa2  = sa1 * sa1
    9493
    95                 rho(k,j,i) =                                                   &
    96                         ( nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
    97                           nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
    98                           nom(7)*sa2       + nom(8)*p1      + nom(9)*p1*pt2  + &
    99                           nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2   &
    100                         ) /                                                    &
    101                         ( den(1)           + den(2)*pt1     + den(3)*pt2     + &
    102                           den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
    103                           den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
    104                           den(10)*sa15*pt2 + den(11)*p1     + den(12)*p2*pt3 + &
    105                           den(13)*p3*pt1                                       &
    106                         )
     94                pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
     95                       nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
     96                       nom(7)*sa2
     97
     98                pden = den(1)           + den(2)*pt1     + den(3)*pt2     + &
     99                       den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
     100                       den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
     101                       den(10)*sa15*pt2
     102
     103!
     104!--             Potential density (without pressure terms)
     105                prho(k,j,i) = pnom / pden
     106
     107                pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  + &
     108                       nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
     109
     110                pden = pden +             den(11)*p1     + den(12)*p2*pt3 + &
     111                       den(13)*p3*pt1
     112
     113!
     114!--             In-situ density
     115                rho(k,j,i) = pnom / pden
    107116
    108117             ENDDO
    109118!
    110119!--          Neumann conditions are assumed at bottom and top boundary
    111              rho(nzt+1,j,i)            = rho(nzt,j,i)
    112              rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
     120             prho(nzt+1,j,i)            = prho(nzt,j,i)
     121             prho(nzb_s_inner(j,i),j,i) = prho(nzb_s_inner(j,i)+1,j,i)
     122             rho(nzt+1,j,i)             = rho(nzt,j,i)
     123             rho(nzb_s_inner(j,i),j,i)  = rho(nzb_s_inner(j,i)+1,j,i)
     124
    113125          ENDDO
    114126       ENDDO
     
    129141       INTEGER ::  i, j, k
    130142
    131        REAL ::  p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
     143       REAL ::  pden, pnom, p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
    132144
    133145       DO  k = nzb_s_inner(j,i)+1, nzt
    134146!
    135147!--       Pressure is needed in dbar
    136 !          p1 = hyp(0) * 1E-4
    137 !          p1 = 0.0
    138148          p1 = hyp(k) * 1E-4
    139149          p2 = p1 * p1
     
    151161          sa2  = sa1 * sa1
    152162
    153           rho(k,j,i) = ( nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
    154                          nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
    155                          nom(7)*sa2       + nom(8)*p1      + nom(9)*p1*pt2  + &
    156                          nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2   &
    157                        ) /                                                    &
    158                        ( den(1)           + den(2)*pt1     + den(3)*pt2     + &
    159                          den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
    160                          den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
    161                          den(10)*sa15*pt2 + den(11)*p1     + den(12)*p2*pt3 + &
    162                          den(13)*p3*pt1                                       &
    163                        )
     163          pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
     164                 nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
     165                 nom(7)*sa2
     166
     167          pden = den(1)           + den(2)*pt1     + den(3)*pt2     + &
     168                 den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
     169                 den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
     170                 den(10)*sa15*pt2
     171
     172!
     173!--       Potential density (without pressure terms)
     174          prho(k,j,i) = pnom / pden
     175
     176          pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  + &
     177                 nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
     178          pden = pden +             den(11)*p1     + den(12)*p2*pt3 + &
     179                 den(13)*p3*pt1
     180
     181!
     182!--       In-situ density
     183          rho(k,j,i) = pnom / pden
     184
     185
    164186       ENDDO
     187
    165188!
    166189!--    Neumann conditions are assumed at bottom and top boundary
    167        rho(nzt+1,j,i)            = rho(nzt,j,i)
    168        rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
     190       prho(nzt+1,j,i)            = prho(nzt,j,i)
     191       prho(nzb_s_inner(j,i),j,i) = prho(nzb_s_inner(j,i)+1,j,i)
     192       rho(nzt+1,j,i)             = rho(nzt,j,i)
     193       rho(nzb_s_inner(j,i),j,i)  = rho(nzb_s_inner(j,i)+1,j,i)
    169194
    170195    END SUBROUTINE eqn_state_seawater_ij
Note: See TracChangeset for help on using the changeset viewer.