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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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