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

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

New:
---
ocean version including prognostic equation for salinity and equation of state for seawater. Routine buoyancy can be used with both temperature and density.
+ inipar-parameters bc_sa_t, bottom_salinityflux, ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, top_salinityflux

advec_s_bc, average_3d_data, boundary_conds, buoyancy, check_parameters, data_output_2d, data_output_3d, diffusion_e, flow_statistics, header, init_grid, init_3d_model, modules, netcdf, parin, production_e, prognostic_equations, read_var_list, sum_up_3d_data, swap_timelevel, time_integration, user_interface, write_var_list, write_3d_binary

New:
eqn_state_seawater, init_ocean

Changed:


inipar-parameter use_pt_reference renamed use_reference

hydro_press renamed hyp, routine calc_mean_pt_profile renamed calc_mean_profile

format adjustments for the ocean version (run_control)

advec_particles, buoyancy, calc_liquid_water_content, check_parameters, diffusion_e, diffusivities, header, init_cloud_physics, modules, production_e, prognostic_equations, run_control

Errors:


Bugfix: height above topography instead of height above level k=0 is used for calculating the mixing length (diffusion_e and diffusivities).

Bugfix: error in boundary condition for TKE removed (advec_s_bc)

advec_s_bc, diffusion_e, prognostic_equations

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