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

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

speed optomizations +bugfix in init_ocean

  • Property svn:keywords set to Id
File size: 2.8 KB
RevLine 
[95]1 SUBROUTINE init_ocean
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
[366]6! Bugfix: First calculation of hyp(0) changed
[95]7!
8! Former revisions:
9! ------------------
[96]10! $Id: init_ocean.f90 366 2009-08-25 08:06:27Z raasch $
[95]11!
[139]12! 124 2007-10-19 15:47:46Z raasch
13! Bugfix: Initial density rho is calculated
14!
[98]15! 97 2007-06-21 08:23:15Z raasch
16! Initial revision
[95]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
[336]26    USE pegrid
[95]27    USE grid_variables
28    USE indices
29
30    IMPLICIT NONE
31
[336]32    INTEGER ::  k, n
[95]33
34    REAL    ::  sa_l, pt_l, rho_l
35
[336]36    REAL, DIMENSION(nzb:nzt+1) ::  rho_init
37
[95]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)
[96]46!-- and the reference density (used later in buoyancy term)
[95]47    hyp(nzt+1) = surface_pressure * 100.0
48
[97]49    hyp(nzt)      = hyp(nzt+1) + rho_surface * g * 0.5 * dzu(nzt+1)
[336]50    rho_init(nzt) = rho_surface
[95]51
[366]52    DO  k = nzt-1, 1, -1
[336]53       hyp(k) = hyp(k+1) + rho_surface * g * dzu(k)
54    ENDDO
[366]55    hyp(0) = hyp(1) + rho_surface * g * dzu(1)
[95]56
[336]57    IF ( myid == 0 )  THEN
58       print*,'hydro pres using rho_surface'
59       DO  k = nzt+1, 0, -1
60          print*, 'k = ', k, ' hyp = ', hyp(k)
61       ENDDO
62       print*, ' '
63    ENDIF
[95]64
[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
87       IF ( myid == 0 )  THEN
88          print*,'hydro pres / rho  n = ', n
89          DO  k = nzt+1, 0, -1
90             print*, 'k = ', k, ' hyp = ', hyp(k), ' rho = ', rho_init(k)
91          ENDDO
92          print*, ' '
93       ENDIF
94
[95]95    ENDDO
96
[97]97!
98!-- Calculate the reference potential density
99    prho_reference = 0.0
100    DO  k = 0, nzt
[96]101
[97]102       sa_l = 0.5 * ( sa_init(k) + sa_init(k+1) )
103       pt_l = 0.5 * ( pt_init(k) + pt_init(k+1) )
104
105       prho_reference = prho_reference + dzu(k+1) * &
[336]106                        eqn_state_seawater_func( 0.0, pt_l, sa_l )
[97]107
108    ENDDO
109
110    prho_reference = prho_reference / ( zu(nzt) - zu(nzb) )
111
[124]112!
113!-- Calculate the initial potential density, based on the initial
114!-- temperature and salinity profile
115    CALL eqn_state_seawater
[97]116
[124]117
[95]118 END SUBROUTINE init_ocean
Note: See TracBrowser for help on using the repository browser.