source: palm/trunk/SOURCE/init_ocean.f90 @ 388

Last change on this file since 388 was 388, checked in by raasch, 12 years ago

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

  • Property svn:keywords set to Id
File size: 2.7 KB
Line 
1 SUBROUTINE init_ocean
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Bugfix: Initial profiles of hydrostatic pressure and density are calculated
7! iteratively. First calculation of hyp(0) changed.
8!
9! Former revisions:
10! ------------------
11! $Id: init_ocean.f90 388 2009-09-23 09:40:33Z raasch $
12!
13! 124 2007-10-19 15:47:46Z raasch
14! Bugfix: Initial density rho is calculated
15!
16! 97 2007-06-21 08:23:15Z raasch
17! Initial revision
18!
19! Description:
20! ------------
21! Initialization of quantities needed for the ocean version
22!------------------------------------------------------------------------------!
23
24    USE arrays_3d
25    USE control_parameters
26    USE eqn_state_seawater_mod
27    USE pegrid
28    USE grid_variables
29    USE indices
30
31    IMPLICIT NONE
32
33    INTEGER ::  k, n
34
35    REAL    ::  sa_l, pt_l
36
37    REAL, DIMENSION(nzb:nzt+1) ::  rho_init
38
39    ALLOCATE( hyp(nzb:nzt+1) )
40
41!
42!-- Set water density near the ocean surface
43    rho_surface = 1027.62
44
45!
46!-- Calculate initial vertical profile of hydrostatic pressure (in Pa)
47!-- and the reference density (used later in buoyancy term)
48!-- First step: Calculate pressure using reference density
49    hyp(nzt+1) = surface_pressure * 100.0
50
51    hyp(nzt)      = hyp(nzt+1) + rho_surface * g * 0.5 * dzu(nzt+1)
52    rho_init(nzt) = rho_surface
53
54    DO  k = nzt-1, 1, -1
55       hyp(k) = hyp(k+1) + rho_surface * g * dzu(k)
56    ENDDO
57    hyp(0) = hyp(1) + rho_surface * g * dzu(1)
58
59!
60!-- Second step: Iteratively calculate in situ density (based on presssure)
61!-- and pressure (based on in situ density)
62    DO  n = 1, 5
63
64       rho_reference = rho_surface * 0.5 * dzu(nzt+1)
65
66       DO  k = nzt-1, 0, -1
67
68          sa_l = 0.5 * ( sa_init(k) + sa_init(k+1) )
69          pt_l = 0.5 * ( pt_init(k) + pt_init(k+1) )
70
71          rho_init(k) = eqn_state_seawater_func( hyp(k), pt_l, sa_l )
72
73          rho_reference = rho_reference + rho_init(k) * dzu(k+1)
74
75       ENDDO
76
77       rho_reference = rho_reference / ( zw(nzt) - zu(nzb) )
78
79       DO  k = nzt-1, 0, -1
80          hyp(k) = hyp(k+1) + g * 0.5 * ( rho_init(k) + rho_init(k+1 ) ) * &
81                              dzu(k+1)
82       ENDDO
83
84    ENDDO
85
86!
87!-- Calculate the reference potential density
88    prho_reference = 0.0
89    DO  k = 0, nzt
90
91       sa_l = 0.5 * ( sa_init(k) + sa_init(k+1) )
92       pt_l = 0.5 * ( pt_init(k) + pt_init(k+1) )
93
94       prho_reference = prho_reference + dzu(k+1) * &
95                        eqn_state_seawater_func( 0.0, pt_l, sa_l )
96
97    ENDDO
98
99    prho_reference = prho_reference / ( zu(nzt) - zu(nzb) )
100
101!
102!-- Calculate the 3d array of initial in situ and potential density,
103!-- based on the initial temperature and salinity profile
104    CALL eqn_state_seawater
105
106
107 END SUBROUTINE init_ocean
Note: See TracBrowser for help on using the repository browser.