source: palm/trunk/SOURCE/subsidence_mod.f90 @ 3198

Last change on this file since 3198 was 3045, checked in by Giersch, 6 years ago

Code adjusted according to coding standards, renamed namelists, error messages revised until PA0347, output CASE 108 disabled

  • Property svn:keywords set to Id
File size: 15.0 KB
Line 
1!> @file subsidence_mod.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-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: subsidence_mod.f90 3045 2018-05-28 07:55:41Z sward $
27! Error message revised
28!
29! 2718 2018-01-02 08:49:38Z maronga
30! Corrected "Former revisions" section
31!
32! 2696 2017-12-14 17:12:51Z kanani
33! Change in file header (GPL part)
34!
35! 2233 2017-05-30 18:08:54Z suehring
36!
37! 2232 2017-05-30 17:47:52Z suehring
38! Adjustments to new topography and surface concept
39!
40! 2000 2016-08-20 18:09:15Z knoop
41! Forced header and separation lines into 80 columns
42!
43! 1862 2016-04-14 09:07:42Z hoffmann
44! Bugfix: In case of vector machine optimized model runs, sums_ls_l should not
45! be addressed if large_scale_forcing is false because sums_ls_l is not
46! allocated.
47!
48! 1850 2016-04-08 13:29:27Z maronga
49! Module renamed
50!
51! 1729 2015-11-20 11:01:48Z gronemeier
52! Bugfix: shifting of initial profiles only once each time step
53!
54! 1682 2015-10-07 23:56:08Z knoop
55! Code annotations made doxygen readable
56!
57! 1489 2014-10-30 08:09:12Z raasch
58! bugfix: sums_ls_l can only be used if large_scale_forcing is switched on
59!
60! 1382 2014-04-30 12:15:41Z boeske
61! Changed the weighting factor that is used in the summation of subsidence
62! tendencies for profile data output from weight_pres to weight_substep
63! added Neumann boundary conditions for profile data output of subsidence terms
64! at nzt+1
65!
66! 1380 2014-04-28 12:40:45Z heinze
67! Shifting only necessary in case of scalar Rayleigh damping
68!
69! 1365 2014-04-22 15:03:56Z boeske
70! Summation of subsidence tendencies for data output added
71! +ls_index, sums_ls_l, tmp_tend
72!
73! 1353 2014-04-08 15:21:23Z heinze
74! REAL constants provided with KIND-attribute
75!
76! 1320 2014-03-20 08:40:49Z raasch
77! ONLY-attribute added to USE-statements,
78! kind-parameters added to all INTEGER and REAL declaration statements,
79! kinds are defined in new module kinds,
80! old module precision_kind is removed,
81! revision history before 2012 removed,
82! comment fields (!:) to be used for variable explanations added to
83! all variable declaration statements
84!
85! 1036 2012-10-22 13:43:42Z raasch
86! code put under GPL (PALM 3.9)
87!
88! Revision 3.7 2009-12-11 14:15:58Z heinze
89! Initial revision
90!
91! Description:
92! ------------
93!> Impact of large-scale subsidence or ascent as tendency term for use
94!> in the prognostic equation of potential temperature. This enables the
95!> construction of a constant boundary layer height z_i with time.
96!-----------------------------------------------------------------------------!
97 MODULE subsidence_mod
98 
99
100
101    IMPLICIT NONE
102
103    PRIVATE
104    PUBLIC  init_w_subsidence, subsidence
105
106    INTERFACE init_w_subsidence
107       MODULE PROCEDURE init_w_subsidence
108    END INTERFACE init_w_subsidence
109
110    INTERFACE subsidence
111       MODULE PROCEDURE subsidence
112       MODULE PROCEDURE subsidence_ij
113    END INTERFACE subsidence
114
115 CONTAINS
116
117!------------------------------------------------------------------------------!
118! Description:
119! ------------
120!> @todo Missing subroutine description.
121!------------------------------------------------------------------------------!
122    SUBROUTINE init_w_subsidence 
123
124       USE arrays_3d,                                                          &
125           ONLY:  dzu, w_subs, zu
126
127       USE control_parameters,                                                 &
128           ONLY:  message_string, ocean, subs_vertical_gradient,               &
129                  subs_vertical_gradient_level, subs_vertical_gradient_level_i
130
131       USE indices,                                                            &
132           ONLY:  nzb, nzt
133
134       USE kinds
135
136       IMPLICIT NONE
137
138       INTEGER(iwp) ::  i !<
139       INTEGER(iwp) ::  k !<
140
141       REAL(wp)     ::  gradient   !<
142       REAL(wp)     ::  ws_surface !<
143
144       IF ( .NOT. ALLOCATED( w_subs ) )  THEN
145          ALLOCATE( w_subs(nzb:nzt+1) )
146          w_subs = 0.0_wp
147       ENDIF
148
149       IF ( ocean )  THEN
150          message_string = 'Applying large scale vertical motion is not ' //   &
151                           'allowed for ocean runs'
152          CALL message( 'init_w_subsidence', 'PA0324', 2, 2, 0, 6, 0 )
153       ENDIF
154
155!
156!--   Compute the profile of the subsidence/ascent velocity
157!--   using the given gradients
158      i = 1
159      gradient = 0.0_wp
160      ws_surface = 0.0_wp
161     
162
163      subs_vertical_gradient_level_i(1) = 0
164      DO  k = 1, nzt+1
165         IF ( i < 11 )  THEN
166            IF ( subs_vertical_gradient_level(i) < zu(k)  .AND. &
167                 subs_vertical_gradient_level(i) >= 0.0_wp )  THEN
168               gradient = subs_vertical_gradient(i) / 100.0_wp
169               subs_vertical_gradient_level_i(i) = k - 1
170               i = i + 1
171            ENDIF
172         ENDIF
173         IF ( gradient /= 0.0_wp )  THEN
174            IF ( k /= 1 )  THEN
175               w_subs(k) = w_subs(k-1) + dzu(k) * gradient
176            ELSE
177               w_subs(k) = ws_surface   + 0.5_wp * dzu(k) * gradient
178            ENDIF
179         ELSE
180            w_subs(k) = w_subs(k-1)
181         ENDIF
182      ENDDO
183
184!
185!--   In case of no given gradients for the subsidence/ascent velocity,
186!--   choose zero gradient
187      IF ( subs_vertical_gradient_level(1) == -9999999.9_wp )  THEN
188         subs_vertical_gradient_level(1) = 0.0_wp
189      ENDIF
190
191    END SUBROUTINE init_w_subsidence
192
193
194!------------------------------------------------------------------------------!
195! Description:
196! ------------
197!> @todo Missing subroutine description.
198!------------------------------------------------------------------------------!
199    SUBROUTINE subsidence( tendency, var, var_init, ls_index ) 
200
201       USE arrays_3d,                                                          &
202           ONLY:  ddzu, w_subs
203
204       USE control_parameters,                                                 &
205           ONLY:  dt_3d, intermediate_timestep_count, large_scale_forcing,     &
206                  scalar_rayleigh_damping
207
208       USE indices,                                                            &
209           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,        &
210                  wall_flags_0
211
212       USE kinds
213
214       USE statistics,                                                         &
215           ONLY:  sums_ls_l, weight_substep
216
217       IMPLICIT NONE
218 
219       INTEGER(iwp) ::  i !<
220       INTEGER(iwp) ::  j !<
221       INTEGER(iwp) ::  k !<
222       INTEGER(iwp) ::  ls_index !<
223
224       REAL(wp)     ::  tmp_tend !<
225       REAL(wp)     ::  tmp_grad !<
226   
227       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var      !<
228       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  tendency !<
229       REAL(wp), DIMENSION(nzb:nzt+1) ::  var_init                     !<
230       REAL(wp), DIMENSION(nzb:nzt+1) ::  var_mod                      !<
231
232       var_mod = var_init
233
234!
235!--    Influence of w_subsidence on the current tendency term
236       DO  i = nxl, nxr
237          DO  j = nys, nyn
238
239             DO  k = nzb+1, nzt 
240                IF ( w_subs(k) < 0.0_wp )  THEN    ! large-scale subsidence
241                   tmp_tend = - w_subs(k) *                                    &
242                              ( var(k+1,j,i) - var(k,j,i) ) * ddzu(k+1) *      &
243                                        MERGE( 1.0_wp, 0.0_wp,                 &
244                                               BTEST( wall_flags_0(k,j,i), 0 ) )
245                ELSE                               ! large-scale ascent
246                   tmp_tend = - w_subs(k) *                                    &
247                              ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k) *        &
248                                        MERGE( 1.0_wp, 0.0_wp,                 &
249                                               BTEST( wall_flags_0(k,j,i), 0 ) )
250                ENDIF
251
252                tendency(k,j,i) = tendency(k,j,i) + tmp_tend
253
254                IF ( large_scale_forcing )  THEN
255                   sums_ls_l(k,ls_index) = sums_ls_l(k,ls_index) + tmp_tend    &
256                                 * weight_substep(intermediate_timestep_count) &
257                                 * MERGE( 1.0_wp, 0.0_wp,                      &
258                                          BTEST( wall_flags_0(k,j,i), 0 ) )
259                ENDIF
260             ENDDO
261
262             IF ( large_scale_forcing )  THEN
263                sums_ls_l(nzt+1,ls_index) = sums_ls_l(nzt,ls_index)
264             ENDIF
265
266          ENDDO
267       ENDDO
268
269!
270!--    Shifting of the initial profile is especially necessary with Rayleigh
271!--    damping switched on
272       IF ( scalar_rayleigh_damping .AND.                                      &
273            intermediate_timestep_count == 1 )  THEN
274          DO  k = nzb, nzt
275             IF ( w_subs(k) < 0.0_wp )  THEN      ! large-scale subsidence
276                var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
277                                  ( var_init(k+1) - var_init(k) ) * ddzu(k+1)
278             ENDIF
279          ENDDO
280!
281!--      At the upper boundary, the initial profile is shifted with aid of
282!--      the gradient tmp_grad. (This is ok if the gradients are linear.)
283         IF ( w_subs(nzt) < 0.0_wp )  THEN
284            tmp_grad = ( var_init(nzt+1) - var_init(nzt) ) * ddzu(nzt+1)
285            var_mod(nzt+1) = var_init(nzt+1) -  &
286                                 dt_3d * w_subs(nzt+1) * tmp_grad
287         ENDIF
288       
289
290         DO  k = nzt+1, nzb+1, -1
291            IF ( w_subs(k) >= 0.0_wp )  THEN  ! large-scale ascent
292               var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
293                                  ( var_init(k) - var_init(k-1) ) * ddzu(k) 
294            ENDIF
295         ENDDO
296!
297!--      At the lower boundary shifting is not necessary because the
298!--      subsidence velocity w_subs(nzb) vanishes.
299         IF ( w_subs(nzb+1) >= 0.0_wp )  THEN
300            var_mod(nzb) = var_init(nzb)
301         ENDIF
302
303         var_init = var_mod
304      ENDIF
305
306
307 END SUBROUTINE subsidence
308
309!------------------------------------------------------------------------------!
310! Description:
311! ------------
312!> @todo Missing subroutine description.
313!------------------------------------------------------------------------------!
314 SUBROUTINE subsidence_ij( i, j, tendency, var, var_init, ls_index ) 
315
316       USE arrays_3d,                                                          &
317           ONLY:  ddzu, w_subs
318
319       USE control_parameters,                                                 &
320           ONLY:  dt_3d, intermediate_timestep_count, large_scale_forcing,     &
321                  scalar_rayleigh_damping
322
323       USE indices,                                                            &
324           ONLY:  nxl, nxlg, nxrg, nyng, nys, nysg, nzb, nzt, wall_flags_0
325
326       USE kinds
327
328       USE statistics,                                                         &
329           ONLY:  sums_ls_l, weight_substep
330
331       IMPLICIT NONE
332 
333       INTEGER(iwp) ::  i !<
334       INTEGER(iwp) ::  j !<
335       INTEGER(iwp) ::  k !<
336       INTEGER(iwp) ::  ls_index !<
337
338       REAL(wp)     ::  tmp_tend !<
339       REAL(wp)     ::  tmp_grad !<
340   
341       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var      !<
342       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  tendency !<
343       REAL(wp), DIMENSION(nzb:nzt+1) ::  var_init                     !<
344       REAL(wp), DIMENSION(nzb:nzt+1) ::  var_mod                      !<
345
346       var_mod = var_init
347
348!
349!--    Influence of w_subsidence on the current tendency term
350       DO  k = nzb+1, nzt 
351          IF ( w_subs(k) < 0.0_wp )  THEN      ! large-scale subsidence
352             tmp_tend = - w_subs(k) * ( var(k+1,j,i) - var(k,j,i) )            &
353                                    * ddzu(k+1)                                &
354                                    * MERGE( 1.0_wp, 0.0_wp,                   &
355                                             BTEST( wall_flags_0(k,j,i), 0 ) )
356          ELSE                                 ! large-scale ascent
357             tmp_tend = - w_subs(k) * ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k)  &
358                                    * MERGE( 1.0_wp, 0.0_wp,                   &
359                                             BTEST( wall_flags_0(k,j,i), 0 ) )
360          ENDIF
361
362          tendency(k,j,i) = tendency(k,j,i) + tmp_tend
363
364          IF ( large_scale_forcing )  THEN
365             sums_ls_l(k,ls_index) = sums_ls_l(k,ls_index) + tmp_tend          &
366                                  * weight_substep(intermediate_timestep_count)&
367                                  * MERGE( 1.0_wp, 0.0_wp,                     &
368                                           BTEST( wall_flags_0(k,j,i), 0 ) )
369          ENDIF
370       ENDDO
371
372       IF ( large_scale_forcing )  THEN
373          sums_ls_l(nzt+1,ls_index) = sums_ls_l(nzt,ls_index)
374       ENDIF
375
376!
377!--    Shifting of the initial profile is especially necessary with Rayleigh
378!--    damping switched on
379       IF ( scalar_rayleigh_damping .AND.                                      &
380            intermediate_timestep_count == 1 )  THEN
381          IF ( i == nxl .AND. j == nys )  THEN ! shifting only once per PE
382
383             DO  k = nzb, nzt
384                IF ( w_subs(k) < 0.0_wp )  THEN      ! large-scale subsidence
385                   var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
386                                     ( var_init(k+1) - var_init(k) ) * ddzu(k+1)
387                ENDIF
388             ENDDO
389!
390!--          At the upper boundary, the initial profile is shifted with aid of
391!--          the gradient tmp_grad. (This is ok if the gradients are linear.)
392             IF ( w_subs(nzt) < 0.0_wp )  THEN
393                tmp_grad = ( var_init(nzt+1) - var_init(nzt) ) * ddzu(nzt+1)
394                var_mod(nzt+1) = var_init(nzt+1) -  &
395                                     dt_3d * w_subs(nzt+1) * tmp_grad
396             ENDIF
397       
398
399             DO  k = nzt+1, nzb+1, -1
400                IF ( w_subs(k) >= 0.0_wp )  THEN  ! large-scale ascent
401                   var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
402                                      ( var_init(k) - var_init(k-1) ) * ddzu(k)
403                ENDIF
404             ENDDO
405!
406!--          At the lower boundary shifting is not necessary because the
407!--          subsidence velocity w_subs(nzb) vanishes.
408             IF ( w_subs(nzb+1) >= 0.0_wp )  THEN
409                var_mod(nzb) = var_init(nzb)
410             ENDIF
411
412             var_init = var_mod 
413
414          ENDIF
415       ENDIF
416
417 END SUBROUTINE subsidence_ij
418
419
420 END MODULE subsidence_mod
Note: See TracBrowser for help on using the repository browser.