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

Last change on this file since 42 was 4, checked in by raasch, 17 years ago

Id keyword set as property for all *.f90 files

  • Property svn:keywords set to Id
File size: 8.0 KB
Line 
1 MODULE buoyancy_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: buoyancy.f90 4 2007-02-13 11:33:16Z raasch $
11! RCS Log replace by Id keyword, revision history cleaned up
12!
13! Revision 1.19  2006/04/26 12:09:56  raasch
14! OpenMP optimization (one dimension added to sums_l)
15!
16! Revision 1.1  1997/08/29 08:56:48  raasch
17! Initial revision
18!
19!
20! Description:
21! ------------
22! Buoyancy term of the third component of the equation of motion.
23! WARNING: humidity is not regarded when using a sloping surface!
24!------------------------------------------------------------------------------!
25
26    PRIVATE
27    PUBLIC buoyancy, calc_mean_pt_profile
28
29    INTERFACE buoyancy
30       MODULE PROCEDURE buoyancy
31       MODULE PROCEDURE buoyancy_ij
32    END INTERFACE buoyancy
33
34    INTERFACE calc_mean_pt_profile
35       MODULE PROCEDURE calc_mean_pt_profile
36    END INTERFACE calc_mean_pt_profile
37
38 CONTAINS
39
40
41!------------------------------------------------------------------------------!
42! Call for all grid points
43!------------------------------------------------------------------------------!
44    SUBROUTINE buoyancy( theta, wind_component, pr )
45
46       USE arrays_3d
47       USE control_parameters
48       USE indices
49       USE pegrid
50       USE statistics
51
52       IMPLICIT NONE
53
54       INTEGER ::  i, j, k, pr, wind_component
55       REAL, DIMENSION(:,:,:), POINTER  ::  theta
56
57
58       IF ( .NOT. sloping_surface )  THEN
59!
60!--       Normal case: horizontal surface
61          DO  i = nxl, nxr
62             DO  j = nys, nyn
63                DO  k = nzb_s_inner(j,i)+1, nzt-1
64                    tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                    &
65                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
66                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
67                                                          )
68                ENDDO
69             ENDDO
70          ENDDO
71
72       ELSE
73!
74!--       Buoyancy term for a surface with a slope in x-direction. The equations
75!--       for both the u and w velocity-component contain proportionate terms.
76!--       Temperature field at time t=0 serves as environmental temperature.
77!--       Reference temperature (pt_surface) is the one at the lower left corner
78!--       of the total domain.
79          IF ( wind_component == 1 )  THEN
80
81             DO  i = nxl, nxr
82                DO  j = nys, nyn
83                   DO  k = nzb_s_inner(j,i)+1, nzt-1
84                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
85                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
86                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
87                                 ) / pt_surface
88                   ENDDO
89                ENDDO
90             ENDDO
91
92          ELSEIF ( wind_component == 3 )  THEN
93
94             DO  i = nxl, nxr
95                DO  j = nys, nyn
96                   DO  k = nzb_s_inner(j,i)+1, nzt-1
97                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
98                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
99                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
100                                 ) / pt_surface
101                   ENDDO
102                ENDDO
103            ENDDO
104
105          ELSE
106
107             IF ( myid == 0 )  PRINT*, '+++ buoyancy: no term for component "',&
108                                       wind_component,'"'
109             CALL local_stop
110
111          ENDIF
112
113       ENDIF
114
115    END SUBROUTINE buoyancy
116
117
118!------------------------------------------------------------------------------!
119! Call for grid point i,j
120!------------------------------------------------------------------------------!
121    SUBROUTINE buoyancy_ij( i, j, theta, wind_component, pr )
122
123       USE arrays_3d
124       USE control_parameters
125       USE indices
126       USE pegrid
127       USE statistics
128
129       IMPLICIT NONE
130
131       INTEGER ::  i, j, k, pr, wind_component
132       REAL, DIMENSION(:,:,:), POINTER  ::  theta
133
134
135       IF ( .NOT. sloping_surface )  THEN
136!
137!--       Normal case: horizontal surface
138          DO  k = nzb_s_inner(j,i)+1, nzt-1
139              tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                          &
140                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
141                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
142                                                    )
143          ENDDO
144
145       ELSE
146!
147!--       Buoyancy term for a surface with a slope in x-direction. The equations
148!--       for both the u and w velocity-component contain proportionate terms.
149!--       Temperature field at time t=0 serves as environmental temperature.
150!--       Reference temperature (pt_surface) is the one at the lower left corner
151!--       of the total domain.
152          IF ( wind_component == 1 )  THEN
153
154             DO  k = nzb_s_inner(j,i)+1, nzt-1
155                tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *            &
156                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
157                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
158                                 ) / pt_surface
159             ENDDO
160
161          ELSEIF ( wind_component == 3 )  THEN
162
163             DO  k = nzb_s_inner(j,i)+1, nzt-1
164                tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *            &
165                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
166                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
167                                 ) / pt_surface
168             ENDDO
169
170          ELSE
171
172             IF ( myid == 0 )  PRINT*, '+++ buoyancy: no term for component "',&
173                                       wind_component,'"'
174             CALL local_stop
175
176          ENDIF
177
178       ENDIF
179
180    END SUBROUTINE buoyancy_ij
181
182
183    SUBROUTINE calc_mean_pt_profile( theta, pr )
184
185!------------------------------------------------------------------------------!
186! Description:
187! ------------
188! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
189! of potential temperature and 44 in case of virtual potential temperature).
190!------------------------------------------------------------------------------!
191
192       USE control_parameters
193       USE indices
194       USE pegrid
195       USE statistics
196
197       IMPLICIT NONE
198
199       INTEGER ::  i, j, k, omp_get_thread_num, pr, tn
200       REAL, DIMENSION(:,:,:), POINTER  ::  theta
201
202!
203!--    Computation of the horizontally averaged temperature profile, unless
204!--    already done by the relevant call from flow_statistics. The calculation
205!--    is done only for the first respective intermediate timestep in order to
206!--    spare communication time and to produce identical model results with jobs
207!--    which are calling flow_statistics at different time intervals.
208!--    Although this calculation is not required for model runs with a slope,
209!--    it is nevertheless also computed.
210       IF ( .NOT. flow_statistics_called  .AND. &
211            intermediate_timestep_count == 1 )  THEN
212
213!
214!--       Horizontal average of the potential temperature
215          tn           =   0  ! Default thread number in case of one thread
216          !$OMP PARALLEL PRIVATE( i, j, k, tn )
217!$        tn = omp_get_thread_num()
218          sums_l(:,pr,tn) = 0.0
219          !$OMP DO
220          DO  i = nxl, nxr
221             DO  j =  nys, nyn
222                DO  k = nzb_s_outer(j,i), nzt+1
223                   sums_l(k,pr,tn) = sums_l(k,pr,tn) + theta(k,j,i)
224                ENDDO
225             ENDDO
226          ENDDO
227          !$OMP END PARALLEL
228
229          DO  i = 1, threads_per_task-1
230             sums_l(:,pr,0) = sums_l(:,pr,0) + sums_l(:,pr,i)
231          ENDDO
232
233#if defined( __parallel )
234
235          CALL MPI_ALLREDUCE( sums_l(nzb,pr,0), sums(nzb,pr), nzt+2-nzb, &
236                              MPI_REAL, MPI_SUM, comm2d, ierr )
237
238#else
239
240          sums(:,pr) = sums_l(:,pr,0)
241
242#endif
243
244          hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_outer(:,0)
245
246       ENDIF
247
248    END SUBROUTINE calc_mean_pt_profile
249
250 END MODULE buoyancy_mod
Note: See TracBrowser for help on using the repository browser.