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

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

typo in file headers removed

  • Property svn:keywords set to Id
File size: 2.7 KB
RevLine 
[95]1 SUBROUTINE init_ocean
2
3!------------------------------------------------------------------------------!
[484]4! Current revisions:
[95]5! -----------------
[392]6!
[95]7!
8! Former revisions:
9! ------------------
[96]10! $Id: init_ocean.f90 484 2010-02-05 07:36:54Z raasch $
[95]11!
[392]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!
[139]16! 124 2007-10-19 15:47:46Z raasch
17! Bugfix: Initial density rho is calculated
18!
[98]19! 97 2007-06-21 08:23:15Z raasch
20! Initial revision
[95]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
[336]30    USE pegrid
[95]31    USE grid_variables
32    USE indices
33
34    IMPLICIT NONE
35
[336]36    INTEGER ::  k, n
[95]37
[388]38    REAL    ::  sa_l, pt_l
[95]39
[336]40    REAL, DIMENSION(nzb:nzt+1) ::  rho_init
41
[95]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)
[96]50!-- and the reference density (used later in buoyancy term)
[388]51!-- First step: Calculate pressure using reference density
[95]52    hyp(nzt+1) = surface_pressure * 100.0
53
[97]54    hyp(nzt)      = hyp(nzt+1) + rho_surface * g * 0.5 * dzu(nzt+1)
[336]55    rho_init(nzt) = rho_surface
[95]56
[366]57    DO  k = nzt-1, 1, -1
[336]58       hyp(k) = hyp(k+1) + rho_surface * g * dzu(k)
59    ENDDO
[366]60    hyp(0) = hyp(1) + rho_surface * g * dzu(1)
[95]61
[388]62!
63!-- Second step: Iteratively calculate in situ density (based on presssure)
64!-- and pressure (based on in situ density)
[336]65    DO  n = 1, 5
[95]66
[336]67       rho_reference = rho_surface * 0.5 * dzu(nzt+1)
[95]68
[336]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
[95]87    ENDDO
88
[97]89!
90!-- Calculate the reference potential density
91    prho_reference = 0.0
92    DO  k = 0, nzt
[96]93
[97]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) * &
[336]98                        eqn_state_seawater_func( 0.0, pt_l, sa_l )
[97]99
100    ENDDO
101
102    prho_reference = prho_reference / ( zu(nzt) - zu(nzb) )
103
[124]104!
[388]105!-- Calculate the 3d array of initial in situ and potential density,
106!-- based on the initial temperature and salinity profile
[124]107    CALL eqn_state_seawater
[97]108
[124]109
[95]110 END SUBROUTINE init_ocean
Note: See TracBrowser for help on using the repository browser.