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

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

preliminary update of bugfixes and extensions for non-cyclic BCs

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