source: palm/trunk/SOURCE/prandtl_fluxes.f90 @ 173

Last change on this file since 173 was 110, checked in by raasch, 16 years ago

New:
---
Allows runs for a coupled atmosphere-ocean LES,
coupling frequency is controlled by new d3par-parameter dt_coupling,
the coupling mode (atmosphere_to_ocean or ocean_to_atmosphere) for the
respective processes is read from environment variable coupling_mode,
which is set by the mpiexec-command,
communication between the two models is done using the intercommunicator
comm_inter,
local files opened by the ocean model get the additional suffic "_O".
Assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean.

A momentum flux can be set as top boundary condition using the new
inipar parameter top_momentumflux_u|v.

Non-cyclic boundary conditions can be used along all horizontal directions.

Quantities w*p* and w"e can be output as vertical profiles.

Initial profiles are reset to constant profiles in case that initializing_actions /= 'set_constant_profiles'. (init_rankine)

Optionally calculate km and kh from initial TKE e_init.

Changed:


Remaining variables iran changed to iran_part (advec_particles, init_particles).

In case that the presure solver is not called for every Runge-Kutta substep
(call_psolver_at_all_substeps = .F.), it is called after the first substep
instead of the last. In that case, random perturbations are also added to the
velocity field after the first substep.

Initialization of km,kh = 0.00001 for ocean = .T. (for ocean = .F. it remains 0.01).

Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T. (instead of cloud_physics = .T.).

Errors:


Bugs from code parts for non-cyclic boundary conditions are removed: loops for
u and v are starting from index nxlu, nysv, respectively. The radiation boundary
condition is used for every Runge-Kutta substep. Velocity phase speeds for
the radiation boundary conditions are calculated for the first Runge-Kutta
substep only and reused for the further substeps. New arrays c_u, c_v, and c_w
are defined for this purpose. Several index errors are removed from the
radiation boundary condition code parts. Upper bounds for calculating
u_0 and v_0 (in production_e) are nxr+1 and nyn+1 because otherwise these
values are not available in case of non-cyclic boundary conditions.

+dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface)

Bugfix: wrong sign removed from the buoyancy production term in the case use_reference = .T. (production_e)

Bugfix: Error message concerning output of particle concentration (pc) modified (check_parameters).

Bugfix: Rayleigh damping for ocean fixed.

  • Property svn:keywords set to Id
File size: 14.2 KB
Line 
1 SUBROUTINE prandtl_fluxes
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: prandtl_fluxes.f90 110 2007-10-05 05:13:14Z raasch $
11!
12! 108 2007-08-24 15:10:38Z letzel
13! assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean
14!
15! 75 2007-03-22 09:54:05Z raasch
16! moisture renamed humidity
17!
18! RCS Log replace by Id keyword, revision history cleaned up
19!
20! Revision 1.19  2006/04/26 12:24:35  raasch
21! +OpenMP directives and optimization (array assignments replaced by DO loops)
22!
23! Revision 1.1  1998/01/23 10:06:06  raasch
24! Initial revision
25!
26!
27! Description:
28! ------------
29! Diagnostic computation of vertical fluxes in the Prandtl layer from the
30! values of the variables at grid point k=1
31!------------------------------------------------------------------------------!
32
33    USE arrays_3d
34    USE control_parameters
35    USE grid_variables
36    USE indices
37
38    IMPLICIT NONE
39
40    INTEGER ::  i, j, k
41    REAL    ::  a, b, e_q, rifm, uv_total, z_p
42
43!
44!-- Compute theta*
45    IF ( constant_heatflux )  THEN
46!
47!--    For a given heat flux in the Prandtl layer:
48!--    for u* use the value from the previous time step
49       !$OMP PARALLEL DO
50       DO  i = nxl-1, nxr+1
51          DO  j = nys-1, nyn+1
52             ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30 )
53!
54!--          ts must be limited, because otherwise overflow may occur in case of
55!--          us=0 when computing rif further below
56             IF ( ts(j,i) < -1.05E5 )  ts = -1.0E5
57             IF ( ts(j,i) >   1.0E5 )  ts =  1.0E5
58          ENDDO
59       ENDDO
60
61    ELSE
62!
63!--    For a given surface temperature:
64!--    (the Richardson number is still the one from the previous time step)
65       !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
66       DO  i = nxl-1, nxr+1
67          DO  j = nys-1, nyn+1
68
69             k   = nzb_s_inner(j,i)
70             z_p = zu(k+1) - zw(k)
71
72             IF ( rif(j,i) >= 0.0 )  THEN
73!
74!--             Stable stratification
75                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
76                                  LOG( z_p / z0(j,i) ) +                   &
77                                  5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
78                                                                )
79             ELSE
80!
81!--             Unstable stratification
82                a = SQRT( 1.0 - 16.0 * rif(j,i) )
83                b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) )
84!
85!--             If a borderline case occurs, the formula for stable
86!--             stratification must be used anyway, or else a zero division
87!--             would occur in the argument of the logarithm
88                IF ( a == 1.0  .OR.  b == 1.0 )  THEN
89                   ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
90                                     LOG( z_p / z0(j,i) ) +                   &
91                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
92                                                                   )
93                ELSE
94                   ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
95                                 LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) &
96                                                                   )
97                ENDIF
98             ENDIF
99
100          ENDDO
101       ENDDO
102    ENDIF
103
104!
105!-- Compute z_p/L (corresponds to the Richardson-flux number)
106    IF ( .NOT. humidity )  THEN
107       !$OMP PARALLEL DO PRIVATE( k, z_p )
108       DO  i = nxl-1, nxr+1
109          DO  j = nys-1, nyn+1
110             k   = nzb_s_inner(j,i)
111             z_p = zu(k+1) - zw(k)
112             rif(j,i) = z_p * kappa * g * ts(j,i) / &
113                        ( pt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )
114!
115!--          Limit the value range of the Richardson numbers.
116!--          This is necessary for very small velocities (u,v --> 0), because
117!--          the absolute value of rif can then become very large, which in
118!--          consequence would result in very large shear stresses and very
119!--          small momentum fluxes (both are generally unrealistic).
120             IF ( rif(j,i) < rif_min )  rif(j,i) = rif_min
121             IF ( rif(j,i) > rif_max )  rif(j,i) = rif_max
122          ENDDO
123       ENDDO
124    ELSE
125       !$OMP PARALLEL DO PRIVATE( k, z_p )
126       DO  i = nxl-1, nxr+1
127          DO  j = nys-1, nyn+1
128             k   = nzb_s_inner(j,i)
129             z_p = zu(k+1) - zw(k)
130             rif(j,i) = z_p * kappa * g *                            &
131                        ( ts(j,i) + 0.61 * pt(k+1,j,i) * qs(j,i) ) / &
132                        ( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )
133!
134!--          Limit the value range of the Richardson numbers.
135!--          This is necessary for very small velocities (u,v --> 0), because
136!--          the absolute value of rif can then become very large, which in
137!--          consequence would result in very large shear stresses and very
138!--          small momentum fluxes (both are generally unrealistic).
139             IF ( rif(j,i) < rif_min )  rif(j,i) = rif_min
140             IF ( rif(j,i) > rif_max )  rif(j,i) = rif_max
141          ENDDO
142       ENDDO       
143    ENDIF
144
145!
146!-- Compute u* at the scalars' grid points
147    !$OMP PARALLEL DO PRIVATE( a, b, k, uv_total, z_p )
148    DO  i = nxl, nxr
149       DO  j = nys, nyn
150
151          k   = nzb_s_inner(j,i)
152          z_p = zu(k+1) - zw(k)
153
154!
155!--       Compute the absolute value of the horizontal velocity
156          uv_total = SQRT( ( 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) ) )**2 + &
157                           ( 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) ) )**2   &
158                         )
159
160          IF ( rif(j,i) >= 0.0 )  THEN
161!
162!--          Stable stratification
163             us(j,i) = kappa * uv_total / (                                &
164                                  LOG( z_p / z0(j,i) ) +                   &
165                                  5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
166                                          )
167          ELSE
168!
169!--          Unstable stratification
170             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
171             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) )
172!
173!--          If a borderline case occurs, the formula for stable stratification
174!--          must be used anyway, or else a zero division would occur in the
175!--          argument of the logarithm.
176             IF ( a == 1.0  .OR.  b == 1.0 )  THEN
177                us(j,i) = kappa * uv_total / (                                &
178                                     LOG( z_p / z0(j,i) ) +                   &
179                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
180                                             )
181             ELSE
182                us(j,i) = kappa * uv_total / (                               &
183                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
184                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
185                                             )
186             ENDIF
187          ENDIF
188       ENDDO
189    ENDDO
190
191!
192!-- Compute u'w' for the total model domain.
193!-- First compute the corresponding component of u* and square it.
194    !$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
195    DO  i = nxl, nxr
196       DO  j = nys, nyn
197
198          k   = nzb_u_inner(j,i)
199          z_p = zu(k+1) - zw(k)
200
201!
202!--       Compute Richardson-flux number for this point
203          rifm = 0.5 * ( rif(j,i-1) + rif(j,i) )
204          IF ( rifm >= 0.0 )  THEN
205!
206!--          Stable stratification
207             usws(j,i) = kappa * u(k+1,j,i) / (                           &
208                                     LOG( z_p / z0(j,i) ) +               &
209                                     5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
210                                              )
211          ELSE
212!
213!--          Unstable stratification
214             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) )
215             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
216!
217!--          If a borderline case occurs, the formula for stable stratification
218!--          must be used anyway, or else a zero division would occur in the
219!--          argument of the logarithm.
220             IF ( a == 1.0  .OR.  B == 1.0 )  THEN
221                usws(j,i) = kappa * u(k+1,j,i) / (                           &
222                                        LOG( z_p / z0(j,i) ) +               &
223                                        5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
224                                                 )
225             ELSE
226                usws(j,i) = kappa * u(k+1,j,i) / (                           &
227                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
228                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
229                                                 )
230             ENDIF
231          ENDIF
232          usws(j,i) = -usws(j,i) * ABS( usws(j,i) )
233       ENDDO
234    ENDDO
235
236!
237!-- Compute v'w' for the total model domain.
238!-- First compute the corresponding component of u* and square it.
239    !$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
240    DO  i = nxl, nxr
241       DO  j = nys, nyn
242
243          k   = nzb_v_inner(j,i)
244          z_p = zu(k+1) - zw(k)
245
246!
247!--       Compute Richardson-flux number for this point
248          rifm = 0.5 * ( rif(j-1,i) + rif(j,i) )
249          IF ( rifm >= 0.0 )  THEN
250!
251!--          Stable stratification
252             vsws(j,i) = kappa * v(k+1,j,i) / (                           &
253                                     LOG( z_p / z0(j,i) ) +               &
254                                     5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
255                                              )
256          ELSE
257!
258!--          Unstable stratification
259             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) )
260             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
261!
262!--          If a borderline case occurs, the formula for stable stratification
263!--          must be used anyway, or else a zero division would occur in the
264!--          argument of the logarithm.
265             IF ( a == 1.0  .OR.  b == 1.0 )  THEN
266                vsws(j,i) = kappa * v(k+1,j,i) / (                           &
267                                        LOG( z_p / z0(j,i) ) +               &
268                                        5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
269                                                 )
270             ELSE
271                vsws(j,i) = kappa * v(k+1,j,i) / (                           &
272                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
273                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
274                                                 )
275             ENDIF
276          ENDIF
277          vsws(j,i) = -vsws(j,i) * ABS( vsws(j,i) )
278       ENDDO
279    ENDDO
280
281!
282!-- If required compute q*
283    IF ( humidity  .OR.  passive_scalar )  THEN
284       IF ( constant_waterflux )  THEN
285!
286!--       For a given water flux in the Prandtl layer:
287          !$OMP PARALLEL DO
288          DO  i = nxl-1, nxr+1
289             DO  j = nys-1, nyn+1
290                qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30 )
291             ENDDO
292          ENDDO
293         
294       ELSE         
295          !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
296          DO  i = nxl-1, nxr+1
297             DO  j = nys-1, nyn+1
298
299                k   = nzb_s_inner(j,i)
300                z_p = zu(k+1) - zw(k)
301
302!
303!--             assume saturation for atmosphere coupled to ocean
304                IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
305                   e_q = 6.1 * &
306                        EXP( 0.07 * ( MIN(pt(0,j,i),pt(1,j,i)) - 273.15 ) )
307                   q(k,j,i) = 0.622 * e_q / ( surface_pressure - e_q )
308                ENDIF
309                IF ( rif(j,i) >= 0.0 )  THEN
310!
311!--                Stable stratification
312                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (         &
313                                  LOG( z_p / z0(j,i) ) +                   &
314                                  5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
315                                                                 )
316                ELSE
317!
318!--                Unstable stratification
319                   a = SQRT( 1.0 - 16.0 * rif(j,i) )
320                   b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) )
321!
322!--                If a borderline case occurs, the formula for stable
323!--                stratification must be used anyway, or else a zero division
324!--                would occur in the argument of the logarithm.
325                   IF ( a == 1.0  .OR.  b == 1.0 )  THEN
326                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (         &
327                                     LOG( z_p / z0(j,i) ) +                   &
328                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
329                                                                    )
330                   ELSE
331                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (          &
332                                  LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) &
333                                                                    )
334                   ENDIF
335                ENDIF
336
337             ENDDO
338          ENDDO
339       ENDIF
340    ENDIF
341
342!
343!-- Exchange the boundaries for u* and the momentum fluxes (fluxes only for
344!-- completeness's sake).
345    CALL exchange_horiz_2d( us )
346    CALL exchange_horiz_2d( usws )
347    CALL exchange_horiz_2d( vsws )
348    IF ( humidity  .OR.  passive_scalar )  CALL exchange_horiz_2d( qsws )
349
350!
351!-- Compute the vertical kinematic heat flux
352    IF ( .NOT. constant_heatflux )  THEN
353       !$OMP PARALLEL DO
354       DO  i = nxl-1, nxr+1
355          DO  j = nys-1, nyn+1
356             shf(j,i) = -ts(j,i) * us(j,i)
357          ENDDO
358       ENDDO
359    ENDIF
360
361!
362!-- Compute the vertical water/scalar flux
363    IF ( .NOT. constant_heatflux .AND. ( humidity .OR. passive_scalar ) ) THEN
364       !$OMP PARALLEL DO
365       DO  i = nxl-1, nxr+1
366          DO  j = nys-1, nyn+1
367             qsws(j,i) = -qs(j,i) * us(j,i)
368          ENDDO
369       ENDDO
370    ENDIF
371
372!
373!-- Bottom boundary condition for the TKE
374    IF ( ibc_e_b == 2 )  THEN
375       !$OMP PARALLEL DO
376       DO  i = nxl-1, nxr+1
377          DO  j = nys-1, nyn+1
378             e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1 )**2
379!
380!--          As a test: cm = 0.4
381!            e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4 )**2
382             e(nzb_s_inner(j,i),j,i)   = e(nzb_s_inner(j,i)+1,j,i)
383          ENDDO
384       ENDDO
385    ENDIF
386
387
388 END SUBROUTINE prandtl_fluxes
Note: See TracBrowser for help on using the repository browser.