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

Last change on this file since 96 was 96, checked in by raasch, 14 years ago

more preliminary uncomplete changes for ocean version

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