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

Last change on this file since 856 was 623, checked in by raasch, 14 years ago

last commit documented

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