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

Last change on this file since 336 was 247, checked in by heinze, 15 years ago

Output of messages replaced by message handling routin

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