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

Last change on this file since 1179 was 1179, checked in by raasch, 8 years ago

New:
---
Initial profiles can be used as reference state in the buoyancy term. New parameter
reference_state introduced. Calculation and handling of reference state in buoyancy term revised.
binary version for restart files changed from 3.9 to 3.9a (no downward compatibility!),
initial profile for rho added to hom (id=77)

Errors:


small bugfix for background communication (time_integration)

  • Property svn:keywords set to Id
File size: 3.9 KB
Line 
1 SUBROUTINE init_ocean
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22! Initial density profile is stored in array hom
23!
24! Former revisions:
25! ------------------
26! $Id: init_ocean.f90 1179 2013-06-14 05:57:58Z raasch $
27!
28! 1036 2012-10-22 13:43:42Z raasch
29! code put under GPL (PALM 3.9)
30!
31! 388 2009-09-23 09:40:33Z raasch
32! Bugfix: Initial profiles of hydrostatic pressure and density are calculated
33! iteratively. First calculation of hyp(0) changed.
34!
35! 124 2007-10-19 15:47:46Z raasch
36! Bugfix: Initial density rho is calculated
37!
38! 97 2007-06-21 08:23:15Z raasch
39! Initial revision
40!
41! Description:
42! ------------
43! Initialization of quantities needed for the ocean version
44!------------------------------------------------------------------------------!
45
46    USE arrays_3d
47    USE control_parameters
48    USE eqn_state_seawater_mod
49    USE grid_variables
50    USE indices
51    USE pegrid
52    USE statistics
53
54    IMPLICIT NONE
55
56    INTEGER ::  k, n
57
58    REAL    ::  sa_l, pt_l
59
60    REAL, DIMENSION(nzb:nzt+1) ::  rho_init
61
62    ALLOCATE( hyp(nzb:nzt+1) )
63
64!
65!-- Set water density near the ocean surface
66    rho_surface = 1027.62
67
68!
69!-- Calculate initial vertical profile of hydrostatic pressure (in Pa)
70!-- and the reference density (used later in buoyancy term)
71!-- First step: Calculate pressure using reference density
72    hyp(nzt+1) = surface_pressure * 100.0
73
74    hyp(nzt)      = hyp(nzt+1) + rho_surface * g * 0.5 * dzu(nzt+1)
75    rho_init(nzt) = rho_surface
76
77    DO  k = nzt-1, 1, -1
78       hyp(k) = hyp(k+1) + rho_surface * g * dzu(k)
79    ENDDO
80    hyp(0) = hyp(1) + rho_surface * g * dzu(1)
81
82!
83!-- Second step: Iteratively calculate in situ density (based on presssure)
84!-- and pressure (based on in situ density)
85    DO  n = 1, 5
86
87       rho_reference = rho_surface * 0.5 * dzu(nzt+1)
88
89       DO  k = nzt-1, 0, -1
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          rho_init(k) = eqn_state_seawater_func( hyp(k), pt_l, sa_l )
95
96          rho_reference = rho_reference + rho_init(k) * dzu(k+1)
97
98       ENDDO
99
100       rho_reference = rho_reference / ( zw(nzt) - zu(nzb) )
101
102       DO  k = nzt-1, 0, -1
103          hyp(k) = hyp(k+1) + g * 0.5 * ( rho_init(k) + rho_init(k+1 ) ) * &
104                              dzu(k+1)
105       ENDDO
106
107    ENDDO
108
109!
110!-- Calculate the reference potential density
111    prho_reference = 0.0
112    DO  k = 0, nzt
113
114       sa_l = 0.5 * ( sa_init(k) + sa_init(k+1) )
115       pt_l = 0.5 * ( pt_init(k) + pt_init(k+1) )
116
117       prho_reference = prho_reference + dzu(k+1) * &
118                        eqn_state_seawater_func( 0.0, pt_l, sa_l )
119
120    ENDDO
121
122    prho_reference = prho_reference / ( zu(nzt) - zu(nzb) )
123
124!
125!-- Calculate the 3d array of initial in situ and potential density,
126!-- based on the initial temperature and salinity profile
127    CALL eqn_state_seawater
128
129!
130!-- Store initial density profile
131    hom(:,1,77,:)  = SPREAD( rho_init(:), 2, statistic_regions+1 )
132
133!
134!-- Set the reference state to be used in the buoyancy terms
135    IF ( use_single_reference_value )  THEN
136       ref_state(:) = prho_reference
137    ELSE
138       ref_state(:) = rho_init(:)
139    ENDIF
140
141
142 END SUBROUTINE init_ocean
Note: See TracBrowser for help on using the repository browser.