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

Last change on this file since 1 was 1, checked in by raasch, 15 years ago

Initial repository layout and content

File size: 10.2 KB
Line 
1 MODULE buoyancy_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: buoyancy.f90,v $
11! Revision 1.19  2006/04/26 12:09:56  raasch
12! OpenMP optimization (one dimension added to sums_l)
13!
14! Revision 1.18  2006/02/23 09:55:33  raasch
15! nzb_2d replaced by nzb_s_inner/outer, nanz_2dh(0) replaced by
16! ngp_2dh_outer(:,0)
17!
18! Revision 1.17  2005/12/06 15:29:16  raasch
19! Mean pt profile is calculated only in case of the first intermediate timestep
20!
21! Revision 1.16  2004/01/30 10:14:31  raasch
22! Scalar lower k index nzb replaced by 2d-array nzb_2d
23!
24! Revision 1.15  2003/10/29 08:41:47  raasch
25! Horizontal mean temperature is now taken from array hom instead of array sums
26!
27! Revision 1.14  2003/03/16 09:27:22  raasch
28! Two underscores (_) are placed in front of all define-strings
29!
30! Revision 1.13  2003/03/12 16:20:47  raasch
31! Full code replaced in the call for all gridpoints instead of calling the
32! _ij version (required by NEC, because otherwise no vectorization)
33!
34! Revision 1.12  2002/12/19 13:47:38  raasch
35! STOP statement replaced by call of subroutine local_stop
36!
37! Revision 1.11  2002/06/11 12:33:07  raasch
38! Former subroutine changed to a module which allows to be called for all grid
39! points of a single vertical column with index i,j or for all grid points by
40! using function overloading.
41! Calculation of the horizontally averaged temperature profile moved to new
42! subroutine calc_mean_pt_profile, which is part of this module.
43!
44! Revision 1.10  2001/03/30 06:54:56  raasch
45! Translation of remaining German identifiers (variables, subroutines, etc.)
46!
47! Revision 1.9  2001/01/22 05:27:03  raasch
48! Module test_variables removed
49!
50! Revision 1.8  2000/07/03 12:55:56  raasch
51! Argument theta declared as pointer
52!
53! Revision 1.7  2000/04/27 07:04:19  raasch
54! in case of a sloping surface the temperature difference is computed relative
55! to the initial (2D) temperature field
56!
57! Revision 1.6  2000/04/13 14:27:34  schroeter
58! considering the influence of humidity to the buoyancy term
59! (only for sloping_surface = false),
60!
61! Revision 1.5  2000/01/21  16:17:23  16:17:23  letzel (Marcus Letzel)
62! All comments translated into English
63!
64! Revision 1.4  1998/09/22 17:16:22  raasch
65! Auftriebsterm erweitert fuer Rechnungen mit geneigter Ebene
66!
67! Revision 1.3  1998/07/06 12:07:27  raasch
68! + USE test_variables
69!
70! Revision 1.2  1998/03/30 11:33:29  raasch
71! nanz_2dh in nanz_2dh(0) geaendert
72!
73! Revision 1.1  1997/08/29 08:56:48  raasch
74! Initial revision
75!
76!
77! Description:
78! ------------
79! Buoyancy term of the third component of the equation of motion.
80! WARNING: humidity is not regarded when using a sloping surface!
81!------------------------------------------------------------------------------!
82
83    PRIVATE
84    PUBLIC buoyancy, calc_mean_pt_profile
85
86    INTERFACE buoyancy
87       MODULE PROCEDURE buoyancy
88       MODULE PROCEDURE buoyancy_ij
89    END INTERFACE buoyancy
90
91    INTERFACE calc_mean_pt_profile
92       MODULE PROCEDURE calc_mean_pt_profile
93    END INTERFACE calc_mean_pt_profile
94
95 CONTAINS
96
97
98!------------------------------------------------------------------------------!
99! Call for all grid points
100!------------------------------------------------------------------------------!
101    SUBROUTINE buoyancy( theta, wind_component, pr )
102
103       USE arrays_3d
104       USE control_parameters
105       USE indices
106       USE pegrid
107       USE statistics
108
109       IMPLICIT NONE
110
111       INTEGER ::  i, j, k, pr, wind_component
112       REAL, DIMENSION(:,:,:), POINTER  ::  theta
113
114
115       IF ( .NOT. sloping_surface )  THEN
116!
117!--       Normal case: horizontal surface
118          DO  i = nxl, nxr
119             DO  j = nys, nyn
120                DO  k = nzb_s_inner(j,i)+1, nzt-1
121                    tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                    &
122                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
123                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
124                                                          )
125                ENDDO
126             ENDDO
127          ENDDO
128
129       ELSE
130!
131!--       Buoyancy term for a surface with a slope in x-direction. The equations
132!--       for both the u and w velocity-component contain proportionate terms.
133!--       Temperature field at time t=0 serves as environmental temperature.
134!--       Reference temperature (pt_surface) is the one at the lower left corner
135!--       of the total domain.
136          IF ( wind_component == 1 )  THEN
137
138             DO  i = nxl, nxr
139                DO  j = nys, nyn
140                   DO  k = nzb_s_inner(j,i)+1, nzt-1
141                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
142                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
143                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
144                                 ) / pt_surface
145                   ENDDO
146                ENDDO
147             ENDDO
148
149          ELSEIF ( wind_component == 3 )  THEN
150
151             DO  i = nxl, nxr
152                DO  j = nys, nyn
153                   DO  k = nzb_s_inner(j,i)+1, nzt-1
154                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
155                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
156                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
157                                 ) / pt_surface
158                   ENDDO
159                ENDDO
160            ENDDO
161
162          ELSE
163
164             IF ( myid == 0 )  PRINT*, '+++ buoyancy: no term for component "',&
165                                       wind_component,'"'
166             CALL local_stop
167
168          ENDIF
169
170       ENDIF
171
172    END SUBROUTINE buoyancy
173
174
175!------------------------------------------------------------------------------!
176! Call for grid point i,j
177!------------------------------------------------------------------------------!
178    SUBROUTINE buoyancy_ij( i, j, theta, wind_component, pr )
179
180       USE arrays_3d
181       USE control_parameters
182       USE indices
183       USE pegrid
184       USE statistics
185
186       IMPLICIT NONE
187
188       INTEGER ::  i, j, k, pr, wind_component
189       REAL, DIMENSION(:,:,:), POINTER  ::  theta
190
191
192       IF ( .NOT. sloping_surface )  THEN
193!
194!--       Normal case: horizontal surface
195          DO  k = nzb_s_inner(j,i)+1, nzt-1
196              tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                          &
197                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
198                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
199                                                    )
200          ENDDO
201
202       ELSE
203!
204!--       Buoyancy term for a surface with a slope in x-direction. The equations
205!--       for both the u and w velocity-component contain proportionate terms.
206!--       Temperature field at time t=0 serves as environmental temperature.
207!--       Reference temperature (pt_surface) is the one at the lower left corner
208!--       of the total domain.
209          IF ( wind_component == 1 )  THEN
210
211             DO  k = nzb_s_inner(j,i)+1, nzt-1
212                tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *            &
213                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
214                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
215                                 ) / pt_surface
216             ENDDO
217
218          ELSEIF ( wind_component == 3 )  THEN
219
220             DO  k = nzb_s_inner(j,i)+1, nzt-1
221                tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *            &
222                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
223                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
224                                 ) / pt_surface
225             ENDDO
226
227          ELSE
228
229             IF ( myid == 0 )  PRINT*, '+++ buoyancy: no term for component "',&
230                                       wind_component,'"'
231             CALL local_stop
232
233          ENDIF
234
235       ENDIF
236
237    END SUBROUTINE buoyancy_ij
238
239
240    SUBROUTINE calc_mean_pt_profile( theta, pr )
241
242!------------------------------------------------------------------------------!
243! Description:
244! ------------
245! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
246! of potential temperature and 44 in case of virtual potential temperature).
247!------------------------------------------------------------------------------!
248
249       USE control_parameters
250       USE indices
251       USE pegrid
252       USE statistics
253
254       IMPLICIT NONE
255
256       INTEGER ::  i, j, k, omp_get_thread_num, pr, tn
257       REAL, DIMENSION(:,:,:), POINTER  ::  theta
258
259!
260!--    Computation of the horizontally averaged temperature profile, unless
261!--    already done by the relevant call from flow_statistics. The calculation
262!--    is done only for the first respective intermediate timestep in order to
263!--    spare communication time and to produce identical model results with jobs
264!--    which are calling flow_statistics at different time intervals.
265!--    Although this calculation is not required for model runs with a slope,
266!--    it is nevertheless also computed.
267       IF ( .NOT. flow_statistics_called  .AND. &
268            intermediate_timestep_count == 1 )  THEN
269
270!
271!--       Horizontal average of the potential temperature
272          tn           =   0  ! Default thread number in case of one thread
273          !$OMP PARALLEL PRIVATE( i, j, k, tn )
274!$        tn = omp_get_thread_num()
275          sums_l(:,pr,tn) = 0.0
276          !$OMP DO
277          DO  i = nxl, nxr
278             DO  j =  nys, nyn
279                DO  k = nzb_s_outer(j,i), nzt+1
280                   sums_l(k,pr,tn) = sums_l(k,pr,tn) + theta(k,j,i)
281                ENDDO
282             ENDDO
283          ENDDO
284          !$OMP END PARALLEL
285
286          DO  i = 1, threads_per_task-1
287             sums_l(:,pr,0) = sums_l(:,pr,0) + sums_l(:,pr,i)
288          ENDDO
289
290#if defined( __parallel )
291
292          CALL MPI_ALLREDUCE( sums_l(nzb,pr,0), sums(nzb,pr), nzt+2-nzb, &
293                              MPI_REAL, MPI_SUM, comm2d, ierr )
294
295#else
296
297          sums(:,pr) = sums_l(:,pr,0)
298
299#endif
300
301          hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_outer(:,0)
302
303       ENDIF
304
305    END SUBROUTINE calc_mean_pt_profile
306
307 END MODULE buoyancy_mod
Note: See TracBrowser for help on using the repository browser.