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

Last change on this file since 550 was 484, checked in by raasch, 14 years ago

typo in file headers removed

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