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

Last change on this file since 515 was 515, checked in by raasch, 12 years ago

simple installation method documented and updated; bugfix for PGI compiler in buoyancy

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