source: palm/trunk/SOURCE/init_rankine.f90 @ 250

Last change on this file since 250 was 110, checked in by raasch, 16 years ago

New:
---
Allows runs for a coupled atmosphere-ocean LES,
coupling frequency is controlled by new d3par-parameter dt_coupling,
the coupling mode (atmosphere_to_ocean or ocean_to_atmosphere) for the
respective processes is read from environment variable coupling_mode,
which is set by the mpiexec-command,
communication between the two models is done using the intercommunicator
comm_inter,
local files opened by the ocean model get the additional suffic "_O".
Assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean.

A momentum flux can be set as top boundary condition using the new
inipar parameter top_momentumflux_u|v.

Non-cyclic boundary conditions can be used along all horizontal directions.

Quantities w*p* and w"e can be output as vertical profiles.

Initial profiles are reset to constant profiles in case that initializing_actions /= 'set_constant_profiles'. (init_rankine)

Optionally calculate km and kh from initial TKE e_init.

Changed:


Remaining variables iran changed to iran_part (advec_particles, init_particles).

In case that the presure solver is not called for every Runge-Kutta substep
(call_psolver_at_all_substeps = .F.), it is called after the first substep
instead of the last. In that case, random perturbations are also added to the
velocity field after the first substep.

Initialization of km,kh = 0.00001 for ocean = .T. (for ocean = .F. it remains 0.01).

Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T. (instead of cloud_physics = .T.).

Errors:


Bugs from code parts for non-cyclic boundary conditions are removed: loops for
u and v are starting from index nxlu, nysv, respectively. The radiation boundary
condition is used for every Runge-Kutta substep. Velocity phase speeds for
the radiation boundary conditions are calculated for the first Runge-Kutta
substep only and reused for the further substeps. New arrays c_u, c_v, and c_w
are defined for this purpose. Several index errors are removed from the
radiation boundary condition code parts. Upper bounds for calculating
u_0 and v_0 (in production_e) are nxr+1 and nyn+1 because otherwise these
values are not available in case of non-cyclic boundary conditions.

+dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface)

Bugfix: wrong sign removed from the buoyancy production term in the case use_reference = .T. (production_e)

Bugfix: Error message concerning output of particle concentration (pc) modified (check_parameters).

Bugfix: Rayleigh damping for ocean fixed.

  • Property svn:keywords set to Id
File size: 4.0 KB
Line 
1 SUBROUTINE init_rankine
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: init_rankine.f90 110 2007-10-05 05:13:14Z letzel $
11!
12! 107 2007-08-17 13:54:45Z raasch
13! Initial profiles are reset to constant profiles
14!
15! 75 2007-03-22 09:54:05Z raasch
16! uxrp, vynp eliminated, 2nd+3rd argument removed from exchange horiz
17!
18! RCS Log replace by Id keyword, revision history cleaned up
19!
20! Revision 1.11  2005/03/26 20:38:49  raasch
21! Arguments for non-cyclic boundary conditions added to argument list of
22! routine exchange_horiz
23!
24! Revision 1.1  1997/08/11 06:18:43  raasch
25! Initial revision
26!
27!
28! Description:
29! ------------
30! Initialize a (nondivergent) Rankine eddy with a vertical axis in order to test
31! the advection terms and the pressure solver.
32!------------------------------------------------------------------------------!
33
34    USE arrays_3d
35    USE constants
36    USE grid_variables
37    USE indices
38    USE control_parameters
39
40    IMPLICIT NONE
41
42    INTEGER ::  i, ic, j, jc, k, kc1, kc2
43    REAL    ::  alpha, betrag, radius, rc, uw, vw, x, y
44
45!
46!-- Default: eddy radius rc, eddy strength z,
47!--          position of eddy centre: ic, jc, kc1, kc2
48    rc  =  4.0 * dx
49    ic  =  ( nx+1 ) / 2
50    jc  =  ic
51    kc1 = nzb
52    kc2 = nzt+1
53
54!
55!-- Reset initial profiles to constant profiles
56    IF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )  THEN
57       DO  i = nxl-1, nxr+1
58          DO  j = nys-1, nyn+1
59             pt(:,j,i) = pt_init
60             u(:,j,i)  = u_init
61             v(:,j,i)  = v_init
62          ENDDO
63       ENDDO
64    ENDIF
65
66!
67!-- Compute the u-component.
68    DO  i = nxl, nxr
69       DO  j = nys, nyn
70          x = ( i - ic - 0.5 ) * dx
71          y = ( j - jc ) * dy
72          radius = SQRT( x**2 + y**2 )
73          IF ( radius <= 2.0 * rc )  THEN
74             betrag = radius / ( 2.0 * rc ) * 0.08
75          ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
76             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
77          ELSE
78             betrag = 0.0
79          ENDIF
80
81          IF ( x == 0.0 )  THEN
82             IF ( y > 0.0 )  THEN
83                alpha = pi / 2.0
84             ELSEIF ( y < 0.0 )  THEN
85                alpha = 3.0 * pi / 2.0
86             ENDIF
87          ELSE
88             IF ( x < 0.0 )  THEN
89                alpha = ATAN( y / x ) + pi
90             ELSE
91                IF ( y < 0.0 )  THEN
92                   alpha = ATAN( y / x ) + 2.0 * pi
93                ELSE
94                   alpha = ATAN( y / x )
95                ENDIF
96             ENDIF
97          ENDIF
98
99          uw = -SIN( alpha ) * betrag
100
101          DO  k = kc1, kc2
102             u(k,j,i) = u(k,j,i) + uw
103          ENDDO
104       ENDDO
105    ENDDO
106
107!
108!-- Compute the v-component.
109    DO  i = nxl, nxr
110       DO  j = nys, nyn
111          x = ( i - ic ) * dx
112          y = ( j - jc - 0.5) * dy
113          radius = SQRT( x**2 + y**2 )
114          IF ( radius <= 2.0 * rc )  THEN
115             betrag = radius / ( 2.0 * rc ) * 0.08
116          ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
117             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
118          ELSE
119             betrag = 0.0
120          ENDIF
121
122          IF ( x == 0.0 )  THEN
123             IF ( y > 0.0 )  THEN
124                alpha = pi / 2.0
125             ELSEIF ( y < 0.0 )  THEN
126                alpha = 3.0 * pi / 2.0
127             ENDIF
128          ELSE
129             IF ( x < 0.0 )  THEN
130                alpha = ATAN( y / x ) + pi
131             ELSE
132                IF ( y < 0.0 )  THEN
133                   alpha = ATAN( y / x ) + 2.0 * pi
134                ELSE
135                   alpha = ATAN( y / x )
136                ENDIF
137             ENDIF
138          ENDIF
139
140          vw = COS( alpha ) * betrag
141
142          DO  k = kc1, kc2
143             v(k,j,i) = v(k,j,i) + vw
144          ENDDO
145       ENDDO
146    ENDDO
147
148!
149!-- Exchange of boundary values for the velocities.
150    CALL exchange_horiz( u )
151    CALL exchange_horiz( v )
152!
153!-- Make velocity field nondivergent.
154    n_sor = nsor_ini
155    CALL pres
156    n_sor = nsor
157
158 END SUBROUTINE init_rankine
Note: See TracBrowser for help on using the repository browser.