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

Last change on this file since 1241 was 1241, checked in by heinze, 8 years ago

Nudging and large scale forcing from external file implemented

  • Property svn:keywords set to Id
File size: 14.4 KB
Line 
1 MODULE buoyancy_mod
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Currrent revisions:
21! ------------------
22! Generalize calc_mean_profile for wider use: use additional steering
23! character loc
24!
25!
26! Former revisions:
27! -----------------
28! $Id: buoyancy.f90 1241 2013-10-30 11:36:58Z heinze $
29!
30! 1179 2013-06-14 05:57:58Z raasch
31! steering of reference state revised (var_reference and pr removed from
32! parameter list), use_reference renamed use_single_reference_value
33!
34! 1171 2013-05-30 11:27:45Z raasch
35! openacc statements added to use_reference-case in accelerator version
36!
37! 1153 2013-05-10 14:33:08Z raasch
38! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
39!
40! 1128 2013-04-12 06:19:32Z raasch
41! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
42! j_north
43!
44! 1036 2012-10-22 13:43:42Z raasch
45! code put under GPL (PALM 3.9)
46!
47! 1015 2012-09-27 09:23:24Z raasch
48! accelerator version (*_acc) added
49!
50! 1010 2012-09-20 07:59:54Z raasch
51! cpp switch __nopointer added for pointer free version
52!
53! 622 2010-12-10 08:08:13Z raasch
54! optional barriers included in order to speed up collective operations
55!
56! 515 2010-03-18 02:30:38Z raasch
57! PGI-compiler creates SIGFPE in routine buoyancy, if opt>1 is used! Therefore,
58! opt=1 is forced by compiler-directive.
59!
60! 247 2009-02-27 14:01:30Z heinze
61! Output of messages replaced by message handling routine
62!
63! 132 2007-11-20 09:46:11Z letzel
64! Vertical scalar profiles now based on nzb_s_inner and ngp_2dh_s_inner.
65!
66! 106 2007-08-16 14:30:26Z raasch
67! i loop for u-component (sloping surface) is starting from nxlu (needed for
68! non-cyclic boundary conditions)
69!
70! 97 2007-06-21 08:23:15Z raasch
71! Routine reneralized to be used with temperature AND density:
72! argument theta renamed var, new argument var_reference,
73! use_pt_reference renamed use_reference,
74! calc_mean_pt_profile renamed calc_mean_profile
75!
76! 57 2007-03-09 12:05:41Z raasch
77! Reference temperature pt_reference can be used.
78!
79! RCS Log replace by Id keyword, revision history cleaned up
80!
81! Revision 1.19  2006/04/26 12:09:56  raasch
82! OpenMP optimization (one dimension added to sums_l)
83!
84! Revision 1.1  1997/08/29 08:56:48  raasch
85! Initial revision
86!
87!
88! Description:
89! ------------
90! Buoyancy term of the third component of the equation of motion.
91! WARNING: humidity is not regarded when using a sloping surface!
92!------------------------------------------------------------------------------!
93
94    PRIVATE
95    PUBLIC buoyancy, buoyancy_acc, calc_mean_profile
96
97    INTERFACE buoyancy
98       MODULE PROCEDURE buoyancy
99       MODULE PROCEDURE buoyancy_ij
100    END INTERFACE buoyancy
101
102    INTERFACE buoyancy_acc
103       MODULE PROCEDURE buoyancy_acc
104    END INTERFACE buoyancy_acc
105
106    INTERFACE calc_mean_profile
107       MODULE PROCEDURE calc_mean_profile
108    END INTERFACE calc_mean_profile
109
110 CONTAINS
111
112
113!------------------------------------------------------------------------------!
114! Call for all grid points
115!------------------------------------------------------------------------------!
116    SUBROUTINE buoyancy( var, wind_component )
117
118       USE arrays_3d
119       USE control_parameters
120       USE indices
121       USE pegrid
122
123       IMPLICIT NONE
124
125       INTEGER ::  i, j, k, wind_component
126#if defined( __nopointer )
127       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
128#else
129       REAL, DIMENSION(:,:,:), POINTER ::  var
130#endif
131
132
133       IF ( .NOT. sloping_surface )  THEN
134!
135!--       Normal case: horizontal surface
136          DO  i = nxl, nxr
137             DO  j = nys, nyn
138                DO  k = nzb_s_inner(j,i)+1, nzt-1
139                   tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (  &
140                          ( var(k,j,i)   - ref_state(k) )   / ref_state(k) +   &
141                          ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)   &
142                                                                            )
143                ENDDO
144             ENDDO
145          ENDDO
146
147       ELSE
148!
149!--       Buoyancy term for a surface with a slope in x-direction. The equations
150!--       for both the u and w velocity-component contain proportionate terms.
151!--       Temperature field at time t=0 serves as environmental temperature.
152!--       Reference temperature (pt_surface) is the one at the lower left corner
153!--       of the total domain.
154          IF ( wind_component == 1 )  THEN
155
156             DO  i = nxlu, nxr
157                DO  j = nys, nyn
158                   DO  k = nzb_s_inner(j,i)+1, nzt-1
159                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
160                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
161                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
162                                 ) / pt_surface
163                   ENDDO
164                ENDDO
165             ENDDO
166
167          ELSEIF ( wind_component == 3 )  THEN
168
169             DO  i = nxl, nxr
170                DO  j = nys, nyn
171                   DO  k = nzb_s_inner(j,i)+1, nzt-1
172                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
173                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
174                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
175                                 ) / pt_surface
176                   ENDDO
177                ENDDO
178            ENDDO
179
180          ELSE
181             
182             WRITE( message_string, * ) 'no term for component "',&
183                                       wind_component,'"'
184             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
185
186          ENDIF
187
188       ENDIF
189
190    END SUBROUTINE buoyancy
191
192
193!------------------------------------------------------------------------------!
194! Call for all grid points - accelerator version
195!------------------------------------------------------------------------------!
196    SUBROUTINE buoyancy_acc( var, wind_component )
197
198       USE arrays_3d
199       USE control_parameters
200       USE indices
201       USE pegrid
202
203       IMPLICIT NONE
204
205       INTEGER ::  i, j, k, wind_component
206#if defined( __nopointer )
207       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
208#else
209       REAL, DIMENSION(:,:,:), POINTER ::  var
210#endif
211
212
213       IF ( .NOT. sloping_surface )  THEN
214!
215!--       Normal case: horizontal surface
216          !$acc kernels present( nzb_s_inner, ref_state, tend, var )
217          !$acc loop
218          DO  i = i_left, i_right
219             DO  j = j_south, j_north
220                !$acc loop independent vector(32)
221                DO  k = nzb_s_inner(j,i)+1, nzt-1
222                   tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (  &
223                          ( var(k,j,i)   - ref_state(k) )   / ref_state(k) +   &
224                          ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)   &
225                                                                            )
226                ENDDO
227             ENDDO
228          ENDDO
229          !$acc end kernels
230
231       ELSE
232!
233!--       Buoyancy term for a surface with a slope in x-direction. The equations
234!--       for both the u and w velocity-component contain proportionate terms.
235!--       Temperature field at time t=0 serves as environmental temperature.
236!--       Reference temperature (pt_surface) is the one at the lower left corner
237!--       of the total domain.
238          IF ( wind_component == 1 )  THEN
239
240             DO  i = nxlu, nxr
241                DO  j = nys, nyn
242                   DO  k = nzb_s_inner(j,i)+1, nzt-1
243                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
244                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
245                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
246                                 ) / pt_surface
247                   ENDDO
248                ENDDO
249             ENDDO
250
251          ELSEIF ( wind_component == 3 )  THEN
252
253             DO  i = nxl, nxr
254                DO  j = nys, nyn
255                   DO  k = nzb_s_inner(j,i)+1, nzt-1
256                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
257                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
258                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
259                                 ) / pt_surface
260                   ENDDO
261                ENDDO
262            ENDDO
263
264          ELSE
265
266             WRITE( message_string, * ) 'no term for component "',&
267                                       wind_component,'"'
268             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
269
270          ENDIF
271
272       ENDIF
273
274    END SUBROUTINE buoyancy_acc
275
276
277!------------------------------------------------------------------------------!
278! Call for grid point i,j
279! ATTENTION: PGI-compiler creates SIGFPE if opt>1 is used! Therefore, opt=1 is
280!            forced by compiler-directive.
281!------------------------------------------------------------------------------!
282!pgi$r opt=1
283    SUBROUTINE buoyancy_ij( i, j, var, wind_component )
284
285       USE arrays_3d
286       USE control_parameters
287       USE indices
288       USE pegrid
289
290       IMPLICIT NONE
291
292       INTEGER ::  i, j, k, pr, wind_component
293#if defined( __nopointer )
294       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
295#else
296       REAL, DIMENSION(:,:,:), POINTER ::  var
297#endif
298
299
300       IF ( .NOT. sloping_surface )  THEN
301!
302!--       Normal case: horizontal surface
303          DO  k = nzb_s_inner(j,i)+1, nzt-1
304              tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (     &
305                        ( var(k,j,i)   - ref_state(k)   ) / ref_state(k)   + &
306                        ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)   &
307                                                                       )
308          ENDDO
309
310       ELSE
311!
312!--       Buoyancy term for a surface with a slope in x-direction. The equations
313!--       for both the u and w velocity-component contain proportionate terms.
314!--       Temperature field at time t=0 serves as environmental temperature.
315!--       Reference temperature (pt_surface) is the one at the lower left corner
316!--       of the total domain.
317          IF ( wind_component == 1 )  THEN
318
319             DO  k = nzb_s_inner(j,i)+1, nzt-1
320                tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *            &
321                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
322                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
323                                 ) / pt_surface
324             ENDDO
325
326          ELSEIF ( wind_component == 3 )  THEN
327
328             DO  k = nzb_s_inner(j,i)+1, nzt-1
329                tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *            &
330                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
331                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
332                                 ) / pt_surface
333             ENDDO
334
335          ELSE
336
337             WRITE( message_string, * ) 'no term for component "',&
338                                       wind_component,'"'
339             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
340
341          ENDIF
342
343       ENDIF
344
345    END SUBROUTINE buoyancy_ij
346
347
348    SUBROUTINE calc_mean_profile( var, pr, loc )
349
350!------------------------------------------------------------------------------!
351! Description:
352! ------------
353! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
354! of potential temperature, 44 in case of virtual potential temperature, and 64
355! in case of density (ocean runs)).
356!------------------------------------------------------------------------------!
357
358       USE arrays_3d,  ONLY: ref_state
359       USE control_parameters
360       USE indices
361       USE pegrid
362       USE statistics
363
364       IMPLICIT NONE
365
366       INTEGER ::  i, j, k, omp_get_thread_num, pr, tn
367       CHARACTER (LEN=*) ::  loc
368#if defined( __nopointer )
369       REAL, DIMENSION(:,:,:) ::  var
370#else
371       REAL, DIMENSION(:,:,:), POINTER ::  var
372#endif
373
374!
375!--    Computation of the horizontally averaged profile of variable var, unless
376!--    already done by the relevant call from flow_statistics. The calculation
377!--    is done only for the first respective intermediate timestep in order to
378!--    spare communication time and to produce identical model results with jobs
379!--    which are calling flow_statistics at different time intervals.
380       IF ( .NOT. flow_statistics_called  .AND. &
381            intermediate_timestep_count == 1 )  THEN
382
383!
384!--       Horizontal average of variable var
385          tn           =   0  ! Default thread number in case of one thread
386          !$OMP PARALLEL PRIVATE( i, j, k, tn )
387!$        tn = omp_get_thread_num()
388          sums_l(:,pr,tn) = 0.0
389          !$OMP DO
390          DO  i = nxl, nxr
391             DO  j =  nys, nyn
392                DO  k = nzb_s_inner(j,i), nzt+1
393                   sums_l(k,pr,tn) = sums_l(k,pr,tn) + var(k,j,i)
394                ENDDO
395             ENDDO
396          ENDDO
397          !$OMP END PARALLEL
398
399          DO  i = 1, threads_per_task-1
400             sums_l(:,pr,0) = sums_l(:,pr,0) + sums_l(:,pr,i)
401          ENDDO
402
403#if defined( __parallel )
404
405          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
406          CALL MPI_ALLREDUCE( sums_l(nzb,pr,0), sums(nzb,pr), nzt+2-nzb, &
407                              MPI_REAL, MPI_SUM, comm2d, ierr )
408
409#else
410
411          sums(:,pr) = sums_l(:,pr,0)
412
413#endif
414
415          hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_s_inner(:,0)
416
417       ENDIF
418
419       SELECT CASE ( loc )
420
421          CASE ( 'time_int' )
422
423             ref_state(:)  = hom(:,1,pr,0)   ! this is used in the buoyancy term
424
425
426          CASE ( 'nudging' )
427             !nothing to be done
428
429
430          CASE DEFAULT
431             message_string = 'unknown location "' // loc // '"'
432             CALL message( 'calc_mean_profile', 'PA0379', 1, 2, 0, 6, 0 )
433
434       END SELECT
435
436
437
438    END SUBROUTINE calc_mean_profile
439
440 END MODULE buoyancy_mod
Note: See TracBrowser for help on using the repository browser.