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

Last change on this file since 622 was 622, checked in by raasch, 11 years ago

New:
---

Optional barriers included in order to speed up collective operations
MPI_ALLTOALL and MPI_ALLREDUCE. This feature is controlled with new initial
parameter collective_wait. Default is .FALSE, but .TRUE. on SGI-type
systems. (advec_particles, advec_s_bc, buoyancy, check_for_restart,
cpu_statistics, data_output_2d, data_output_ptseries, flow_statistics,
global_min_max, inflow_turbulence, init_3d_model, init_particles, init_pegrid,
init_slope, parin, pres, poismg, set_particle_attributes, timestep,
read_var_list, user_statistics, write_compressed, write_var_list)

Adjustments for Kyushu Univ. (lcrte, ibmku). Concerning hybrid
(MPI/openMP) runs, the number of openMP threads per MPI tasks can now
be given as an argument to mrun-option -O. (mbuild, mrun, subjob)

Changed:


Initialization of the module command changed for SGI-ICE/lcsgi (mbuild, subjob)

Errors:


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