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

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

several small bugfixes; some more dvrp changes

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