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

Last change on this file since 4867 was 4828, checked in by Giersch, 4 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

  • Property svn:keywords set to Id
File size: 8.3 KB
Line 
1!> @file buoyancy.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16!
17! Copyright 1997-2021 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: buoyancy.f90 4828 2021-01-05 11:21:41Z moh.hefny $
27! file re-formatted to follow the PALM coding standard
28!
29! 4360 2020-01-07 11:25:50Z suehring
30! Corrected "Former revisions" section
31!
32! 3725 2019-02-07 10:11:02Z raasch
33! unused variables removed
34!
35! 3655 2019-01-07 16:51:22Z knoop
36! nopointer option removed
37!
38! Revision 1.1  1997/08/29 08:56:48  raasch
39! Initial revision
40!
41!
42! Description:
43! ------------
44!> Buoyancy term of the third component of the equation of motion.
45!> @attention Humidity is not regarded when using a sloping surface!
46!--------------------------------------------------------------------------------------------------!
47 MODULE buoyancy_mod
48
49    USE basic_constants_and_equations_mod,                                                         &
50        ONLY:  g
51
52    USE arrays_3d,                                                                                 &
53        ONLY:  pt, pt_slope_ref, ref_state, tend
54
55    USE control_parameters,                                                                        &
56        ONLY:  atmos_ocean_sign, cos_alpha_surface, message_string, pt_surface, sin_alpha_surface, &
57               sloping_surface
58
59    USE kinds
60
61    PRIVATE
62    PUBLIC buoyancy
63
64    INTERFACE buoyancy
65       MODULE PROCEDURE buoyancy
66       MODULE PROCEDURE buoyancy_ij
67    END INTERFACE buoyancy
68
69 CONTAINS
70
71
72!--------------------------------------------------------------------------------------------------!
73! Description:
74! ------------
75!> Call for all grid points
76!--------------------------------------------------------------------------------------------------!
77    SUBROUTINE buoyancy( var, wind_component )
78
79       USE indices,                                                                                &
80           ONLY:  nxl, nxlu, nxr, nyn, nys, nzb, nzt
81
82
83       IMPLICIT NONE
84
85       INTEGER(iwp) ::  i              !<
86       INTEGER(iwp) ::  j              !<
87       INTEGER(iwp) ::  k              !<
88       INTEGER(iwp) ::  wind_component !<
89
90       REAL(wp), DIMENSION(:,:,:), POINTER ::  var
91
92
93       IF ( .NOT. sloping_surface )  THEN
94!
95!--       Normal case: horizontal surface
96          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
97          !$ACC PRESENT(var) &
98          !$ACC PRESENT(ref_state) &
99          !$ACC PRESENT(tend)
100          DO  i = nxl, nxr
101             DO  j = nys, nyn
102                DO  k = nzb+1, nzt-1
103                   tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp *                     &
104                                 (                                                                 &
105                                    ( var(k,j,i)   - ref_state(k) )   / ref_state(k) +             &
106                                    ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)             &
107                                 )
108                ENDDO
109             ENDDO
110          ENDDO
111
112       ELSE
113!
114!--       Buoyancy term for a surface with a slope in x-direction. The equations for both the u and
115!         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 of the total
118!--       domain.
119          IF ( wind_component == 1 )  THEN
120
121             DO  i = nxlu, nxr
122                DO  j = nys, nyn
123                   DO  k = nzb+1, nzt-1
124                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *                          &
125                                    0.5_wp * ( ( 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+1, nzt-1
137                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *                          &
138                                    0.5_wp * ( ( 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 "', wind_component, '"'
148             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
149
150          ENDIF
151
152       ENDIF
153
154    END SUBROUTINE buoyancy
155
156
157!--------------------------------------------------------------------------------------------------!
158! Description:
159! ------------
160!> Call for grid point i,j
161!> @attention PGI-compiler creates SIGFPE if opt>1 is used! Therefore, opt=1 is forced by
162!>            compiler-directive.
163!--------------------------------------------------------------------------------------------------!
164!pgi$r opt=1
165    SUBROUTINE buoyancy_ij( i, j, var, wind_component )
166
167
168       USE indices,                                                                                &
169           ONLY:  nzb, nzt
170
171       IMPLICIT NONE
172
173       INTEGER(iwp) ::  i              !<
174       INTEGER(iwp) ::  j              !<
175       INTEGER(iwp) ::  k              !<
176       INTEGER(iwp) ::  wind_component !<
177
178       REAL(wp), DIMENSION(:,:,:), POINTER ::  var
179
180
181       IF ( .NOT. sloping_surface )  THEN
182!
183!--       Normal case: horizontal surface
184          DO  k = nzb+1, nzt-1
185              tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * (                        &
186                            ( var(k,j,i)   - ref_state(k)   ) / ref_state(k)   +                   &
187                            ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)                     &
188                                                                          )
189          ENDDO
190
191       ELSE
192!
193!--       Buoyancy term for a surface with a slope in x-direction. The equations for both the u and
194!--       w velocity-component contain proportionate terms.
195!--       Temperature field at time t=0 serves as environmental temperature.
196!--       Reference temperature (pt_surface) is the one at the lower left corner of the total
197!--       domain.
198          IF ( wind_component == 1 )  THEN
199
200             DO  k = nzb+1, nzt-1
201                tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *                                &
202                              0.5_wp * ( ( pt(k,j,i-1)         + pt(k,j,i)         )               &
203                                       - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) )               &
204                                       ) / pt_surface
205             ENDDO
206
207          ELSEIF ( wind_component == 3 )  THEN
208
209             DO  k = nzb+1, nzt-1
210                tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *                                &
211                              0.5_wp * ( ( pt(k,j,i)         + pt(k+1,j,i)         )               &
212                                       - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) )               &
213                                       ) / pt_surface
214             ENDDO
215
216          ELSE
217
218             WRITE( message_string, * ) 'no term for component "', wind_component, '"'
219             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
220
221          ENDIF
222
223       ENDIF
224
225    END SUBROUTINE buoyancy_ij
226
227 END MODULE buoyancy_mod
Note: See TracBrowser for help on using the repository browser.