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

Last change on this file since 3523 was 3302, checked in by raasch, 6 years ago

Craik-Leibovich force and wave breaking effects added to ocean mode, header output for ocean mode

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