source: palm/tags/release-3.3/SOURCE/sor.f90 @ 2238

Last change on this file since 2238 was 77, checked in by raasch, 17 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:


General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:


Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

  • Property svn:keywords set to Id
File size: 5.7 KB
Line 
1 SUBROUTINE sor( d, ddzu, ddzw, p )
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: sor.f90 77 2007-03-29 04:26:56Z suehring $
11!
12! 75 2007-03-22 09:54:05Z raasch
13! 2nd+3rd argument removed from exchange horiz
14!
15! RCS Log replace by Id keyword, revision history cleaned up
16!
17! Revision 1.9  2005/03/26 21:02:23  raasch
18! Implementation of non-cyclic (Neumann) horizontal boundary conditions,
19! dx2,dy2 replaced by ddx2,ddy2
20!
21! Revision 1.1  1997/08/11 06:25:56  raasch
22! Initial revision
23!
24!
25! Description:
26! ------------
27! Solve the Poisson-equation with the SOR-Red/Black-scheme.
28!------------------------------------------------------------------------------!
29
30    USE grid_variables
31    USE indices
32    USE pegrid
33    USE control_parameters
34
35    IMPLICIT NONE
36
37    INTEGER ::  i, j, k, n, nxl1, nxl2, nys1, nys2
38    REAL    ::  ddzu(1:nz+1), ddzw(1:nz)
39    REAL    ::  d(nzb+1:nzt,nys:nyn,nxl:nxr),         &
40                p(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
41    REAL, DIMENSION(:), ALLOCATABLE ::  f1, f2, f3
42
43    ALLOCATE( f1(1:nz), f2(1:nz), f3(1:nz) )
44
45!
46!-- Compute pre-factors.
47    DO  k = 1, nz
48         f2(k) = ddzu(k+1) * ddzw(k)
49         f3(k) = ddzu(k)   * ddzw(k)
50         f1(k) = 2.0 * ( ddx2 + ddy2 ) + f2(k) + f3(k)
51    ENDDO
52
53!
54!-- Limits for RED- and BLACK-part.
55    IF ( MOD( nxl , 2 ) == 0 )  THEN
56       nxl1 = nxl
57       nxl2 = nxl + 1
58    ELSE
59       nxl1 = nxl + 1
60       nxl2 = nxl
61    ENDIF
62    IF ( MOD( nys , 2 ) == 0 )  THEN
63       nys1 = nys
64       nys2 = nys + 1
65    ELSE
66       nys1 = nys + 1
67       nys2 = nys
68    ENDIF
69
70    DO  n = 1, n_sor
71
72!
73!--    RED-part
74       DO  i = nxl1, nxr, 2
75          DO  j = nys2, nyn, 2
76             DO  k = nzb+1, nzt
77                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
78                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
79                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
80                               f2(k) * p(k+1,j,i)                 +    &
81                               f3(k) * p(k-1,j,i)                 -    &
82                               d(k,j,i)                           -    &
83                               f1(k) * p(k,j,i)           )
84             ENDDO
85          ENDDO
86       ENDDO
87
88       DO  i = nxl2, nxr, 2
89          DO  j = nys1, nyn, 2
90             DO  k = nzb+1, nzt
91                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
92                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
93                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
94                               f2(k) * p(k+1,j,i)                 +    &
95                               f3(k) * p(k-1,j,i)                 -    &
96                               d(k,j,i)                           -    &
97                               f1(k) * p(k,j,i)           )
98             ENDDO
99          ENDDO
100       ENDDO
101
102!
103!--    Exchange of boundary values for p.
104       CALL exchange_horiz( p )
105
106!
107!--    Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
108       IF ( bc_lr /= 'cyclic' )  THEN
109          IF ( inflow_l .OR. outflow_l )  p(:,:,nxl-1) = p(:,:,nxl)
110          IF ( inflow_r .OR. outflow_r )  p(:,:,nxr+1) = p(:,:,nxr)
111       ENDIF
112       IF ( bc_ns /= 'cyclic' )  THEN
113          IF ( inflow_n .OR. outflow_n )  p(:,nyn+1,:) = p(:,nyn,:)
114          IF ( inflow_s .OR. outflow_s )  p(:,nys-1,:) = p(:,nys,:)
115       ENDIF
116
117!
118!--    BLACK-part
119       DO  i = nxl1, nxr, 2
120          DO  j = nys1, nyn, 2
121             DO  k = nzb+1, nzt
122                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
123                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
124                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
125                               f2(k) * p(k+1,j,i)                 +    &
126                               f3(k) * p(k-1,j,i)                 -    &
127                               d(k,j,i)                           -    &
128                               f1(k) * p(k,j,i)           )
129             ENDDO
130          ENDDO
131       ENDDO
132
133       DO  i = nxl2, nxr, 2
134          DO  j = nys2, nyn, 2
135             DO  k = nzb+1, nzt
136                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
137                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
138                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
139                               f2(k) * p(k+1,j,i)                 +    &
140                               f3(k) * p(k-1,j,i)                 -    &
141                               d(k,j,i)                           -    &
142                               f1(k) * p(k,j,i)           )
143             ENDDO
144          ENDDO
145       ENDDO
146
147!
148!--    Exchange of boundary values for p.
149       CALL exchange_horiz( p )
150
151!
152!--    Boundary conditions top/bottom.
153!--    Bottom boundary
154       IF ( ibc_p_b == 1 )  THEN
155!
156!--       Neumann
157          p(nzb,:,:) = p(nzb+1,:,:)
158       ELSE
159!
160!--       Dirichlet
161          p(nzb,:,:) = 0.0
162       ENDIF
163
164!
165!--    Top boundary
166       IF ( ibc_p_t == 1 )  THEN
167!
168!--       Neumann
169          p(nzt+1,:,:) = p(nzt,:,:)
170       ELSE
171!
172!--       Dirichlet
173          p(nzt+1,:,:) = 0.0
174       ENDIF
175
176!
177!--    Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
178       IF ( bc_lr /= 'cyclic' )  THEN
179          IF ( inflow_l .OR. outflow_l )  p(:,:,nxl-1) = p(:,:,nxl)
180          IF ( inflow_r .OR. outflow_r )  p(:,:,nxr+1) = p(:,:,nxr)
181       ENDIF
182       IF ( bc_ns /= 'cyclic' )  THEN
183          IF ( inflow_n .OR. outflow_n )  p(:,nyn+1,:) = p(:,nyn,:)
184          IF ( inflow_s .OR. outflow_s )  p(:,nys-1,:) = p(:,nys,:)
185       ENDIF
186
187    ENDDO
188
189    DEALLOCATE( f1, f2, f3 )
190
191 END SUBROUTINE sor
Note: See TracBrowser for help on using the repository browser.