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

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

preliminary update of further changes, advec_particles is not running!

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