source: palm/trunk/SOURCE/buoyancy.f90 @ 132

Last change on this file since 132 was 132, checked in by letzel, 16 years ago

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars and velocity components) and
ngp_2dh (staggered velocity components and their products), respectively.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

  • Property svn:keywords set to Id
File size: 9.7 KB
Line 
1 MODULE buoyancy_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Vertical scalar profiles now based on nzb_s_inner and ngp_2dh_s_inner.
7!
8! Former revisions:
9! -----------------
10! $Id: buoyancy.f90 132 2007-11-20 09:46:11Z letzel $
11!
12! 106 2007-08-16 14:30:26Z raasch
13! i loop for u-component (sloping surface) is starting from nxlu (needed for
14! non-cyclic boundary conditions)
15!
16! 97 2007-06-21 08:23:15Z raasch
17! Routine reneralized to be used with temperature AND density:
18! argument theta renamed var, new argument var_reference,
19! use_pt_reference renamed use_reference,
20! calc_mean_pt_profile renamed calc_mean_profile
21!
22! 57 2007-03-09 12:05:41Z raasch
23! Reference temperature pt_reference can be used.
24!
25! RCS Log replace by Id keyword, revision history cleaned up
26!
27! Revision 1.19  2006/04/26 12:09:56  raasch
28! OpenMP optimization (one dimension added to sums_l)
29!
30! Revision 1.1  1997/08/29 08:56:48  raasch
31! Initial revision
32!
33!
34! Description:
35! ------------
36! Buoyancy term of the third component of the equation of motion.
37! WARNING: humidity is not regarded when using a sloping surface!
38!------------------------------------------------------------------------------!
39
40    PRIVATE
41    PUBLIC buoyancy, calc_mean_profile
42
43    INTERFACE buoyancy
44       MODULE PROCEDURE buoyancy
45       MODULE PROCEDURE buoyancy_ij
46    END INTERFACE buoyancy
47
48    INTERFACE calc_mean_profile
49       MODULE PROCEDURE calc_mean_profile
50    END INTERFACE calc_mean_profile
51
52 CONTAINS
53
54
55!------------------------------------------------------------------------------!
56! Call for all grid points
57!------------------------------------------------------------------------------!
58    SUBROUTINE buoyancy( var, var_reference, wind_component, pr )
59
60       USE arrays_3d
61       USE control_parameters
62       USE indices
63       USE pegrid
64       USE statistics
65
66       IMPLICIT NONE
67
68       INTEGER ::  i, j, k, pr, wind_component
69       REAL    ::  var_reference
70       REAL, DIMENSION(:,:,:), POINTER  ::  var
71
72
73       IF ( .NOT. sloping_surface )  THEN
74!
75!--       Normal case: horizontal surface
76          IF ( use_reference )  THEN
77             DO  i = nxl, nxr
78                DO  j = nys, nyn
79                   DO  k = nzb_s_inner(j,i)+1, nzt-1
80                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
81                                                            (                  &
82                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
83                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
84                                                            )
85                   ENDDO
86                ENDDO
87             ENDDO
88          ELSE
89             DO  i = nxl, nxr
90                DO  j = nys, nyn
91                   DO  k = nzb_s_inner(j,i)+1, nzt-1
92                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
93                                                            (                  &
94                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
95                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
96                                                            )
97                   ENDDO
98                ENDDO
99             ENDDO
100          ENDIF
101
102       ELSE
103!
104!--       Buoyancy term for a surface with a slope in x-direction. The equations
105!--       for both the u and w velocity-component contain proportionate terms.
106!--       Temperature field at time t=0 serves as environmental temperature.
107!--       Reference temperature (pt_surface) is the one at the lower left corner
108!--       of the total domain.
109          IF ( wind_component == 1 )  THEN
110
111             DO  i = nxlu, nxr
112                DO  j = nys, nyn
113                   DO  k = nzb_s_inner(j,i)+1, nzt-1
114                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
115                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
116                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
117                                 ) / pt_surface
118                   ENDDO
119                ENDDO
120             ENDDO
121
122          ELSEIF ( wind_component == 3 )  THEN
123
124             DO  i = nxl, nxr
125                DO  j = nys, nyn
126                   DO  k = nzb_s_inner(j,i)+1, nzt-1
127                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
128                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
129                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
130                                 ) / pt_surface
131                   ENDDO
132                ENDDO
133            ENDDO
134
135          ELSE
136
137             IF ( myid == 0 )  PRINT*, '+++ buoyancy: no term for component "',&
138                                       wind_component,'"'
139             CALL local_stop
140
141          ENDIF
142
143       ENDIF
144
145    END SUBROUTINE buoyancy
146
147
148!------------------------------------------------------------------------------!
149! Call for grid point i,j
150!------------------------------------------------------------------------------!
151    SUBROUTINE buoyancy_ij( i, j, var, var_reference, wind_component, pr )
152
153       USE arrays_3d
154       USE control_parameters
155       USE indices
156       USE pegrid
157       USE statistics
158
159       IMPLICIT NONE
160
161       INTEGER ::  i, j, k, pr, wind_component
162       REAL    ::  var_reference
163       REAL, DIMENSION(:,:,:), POINTER  ::  var
164
165
166       IF ( .NOT. sloping_surface )  THEN
167!
168!--       Normal case: horizontal surface
169          IF ( use_reference )  THEN
170             DO  k = nzb_s_inner(j,i)+1, nzt-1
171                 tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (   &
172                         ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
173                         ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
174                                                                          )
175             ENDDO
176          ELSE
177             DO  k = nzb_s_inner(j,i)+1, nzt-1
178                 tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (    &
179                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
180                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
181                                                                          )
182             ENDDO
183          ENDIF
184
185       ELSE
186!
187!--       Buoyancy term for a surface with a slope in x-direction. The equations
188!--       for both the u and w velocity-component contain proportionate terms.
189!--       Temperature field at time t=0 serves as environmental temperature.
190!--       Reference temperature (pt_surface) is the one at the lower left corner
191!--       of the total domain.
192          IF ( wind_component == 1 )  THEN
193
194             DO  k = nzb_s_inner(j,i)+1, nzt-1
195                tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *            &
196                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
197                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
198                                 ) / pt_surface
199             ENDDO
200
201          ELSEIF ( wind_component == 3 )  THEN
202
203             DO  k = nzb_s_inner(j,i)+1, nzt-1
204                tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *            &
205                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
206                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
207                                 ) / pt_surface
208             ENDDO
209
210          ELSE
211
212             IF ( myid == 0 )  PRINT*, '+++ buoyancy: no term for component "',&
213                                       wind_component,'"'
214             CALL local_stop
215
216          ENDIF
217
218       ENDIF
219
220    END SUBROUTINE buoyancy_ij
221
222
223    SUBROUTINE calc_mean_profile( var, pr )
224
225!------------------------------------------------------------------------------!
226! Description:
227! ------------
228! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
229! of potential temperature and 44 in case of virtual potential temperature).
230!------------------------------------------------------------------------------!
231
232       USE control_parameters
233       USE indices
234       USE pegrid
235       USE statistics
236
237       IMPLICIT NONE
238
239       INTEGER ::  i, j, k, omp_get_thread_num, pr, tn
240       REAL, DIMENSION(:,:,:), POINTER  ::  var
241
242!
243!--    Computation of the horizontally averaged profile of variable var, unless
244!--    already done by the relevant call from flow_statistics. The calculation
245!--    is done only for the first respective intermediate timestep in order to
246!--    spare communication time and to produce identical model results with jobs
247!--    which are calling flow_statistics at different time intervals.
248       IF ( .NOT. flow_statistics_called  .AND. &
249            intermediate_timestep_count == 1 )  THEN
250
251!
252!--       Horizontal average of variable var
253          tn           =   0  ! Default thread number in case of one thread
254          !$OMP PARALLEL PRIVATE( i, j, k, tn )
255!$        tn = omp_get_thread_num()
256          sums_l(:,pr,tn) = 0.0
257          !$OMP DO
258          DO  i = nxl, nxr
259             DO  j =  nys, nyn
260                DO  k = nzb_s_inner(j,i), nzt+1
261                   sums_l(k,pr,tn) = sums_l(k,pr,tn) + var(k,j,i)
262                ENDDO
263             ENDDO
264          ENDDO
265          !$OMP END PARALLEL
266
267          DO  i = 1, threads_per_task-1
268             sums_l(:,pr,0) = sums_l(:,pr,0) + sums_l(:,pr,i)
269          ENDDO
270
271#if defined( __parallel )
272
273          CALL MPI_ALLREDUCE( sums_l(nzb,pr,0), sums(nzb,pr), nzt+2-nzb, &
274                              MPI_REAL, MPI_SUM, comm2d, ierr )
275
276#else
277
278          sums(:,pr) = sums_l(:,pr,0)
279
280#endif
281
282          hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_s_inner(:,0)
283
284       ENDIF
285
286    END SUBROUTINE calc_mean_profile
287
288 END MODULE buoyancy_mod
Note: See TracBrowser for help on using the repository browser.