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

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

updating comments and rc-file

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