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

Last change on this file since 4381 was 4360, checked in by suehring, 5 years ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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