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

Last change on this file since 4555 was 4360, checked in by suehring, 5 years ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

  • Property svn:keywords set to Id
File size: 14.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!
[4360]17! Copyright 1997-2020 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 4360 2020-01-07 11:25:50Z moh.hefny $
[4346]27! Introduction of wall_flags_total_0, which currently sets bits based on static
28! topography information used in wall_flags_static_0
29!
30! 4329 2019-12-10 15:46:36Z motisi
[4329]31! Renamed wall_flags_0 to wall_flags_static_0
32!
33! 4182 2019-08-22 15:20:23Z scharf
[4182]34! Corrected "Former revisions" section
35!
36! 3655 2019-01-07 16:51:22Z knoop
[3554]37! add subroutine and variable description
[1383]38!
[4182]39! Revision 3.7 2009-12-11 14:15:58Z heinze
40! Initial revision
41!
[411]42! Description:
43! ------------
[1682]44!> Impact of large-scale subsidence or ascent as tendency term for use
45!> in the prognostic equation of potential temperature. This enables the
46!> construction of a constant boundary layer height z_i with time.
[411]47!-----------------------------------------------------------------------------!
[1682]48 MODULE subsidence_mod
49 
[411]50
51
52    IMPLICIT NONE
53
54    PRIVATE
55    PUBLIC  init_w_subsidence, subsidence
56
57    INTERFACE init_w_subsidence
58       MODULE PROCEDURE init_w_subsidence
59    END INTERFACE init_w_subsidence
60
61    INTERFACE subsidence
62       MODULE PROCEDURE subsidence
63       MODULE PROCEDURE subsidence_ij
64    END INTERFACE subsidence
65
66 CONTAINS
67
[1682]68!------------------------------------------------------------------------------!
69! Description:
70! ------------
[3554]71!> Initialize vertical subsidence velocity w_subs.
[1682]72!------------------------------------------------------------------------------!
[411]73    SUBROUTINE init_w_subsidence 
74
[1320]75       USE arrays_3d,                                                          &
76           ONLY:  dzu, w_subs, zu
[411]77
[1320]78       USE control_parameters,                                                 &
[3294]79           ONLY:  message_string, ocean_mode, subs_vertical_gradient,          &
[1320]80                  subs_vertical_gradient_level, subs_vertical_gradient_level_i
81
82       USE indices,                                                            &
83           ONLY:  nzb, nzt
84
85       USE kinds
86
[411]87       IMPLICIT NONE
88
[3554]89       INTEGER(iwp) ::  i !< loop index
90       INTEGER(iwp) ::  k !< loop index
[411]91
[3554]92       REAL(wp)     ::  gradient   !< vertical gradient of subsidence velocity
93       REAL(wp)     ::  ws_surface !< subsidence velocity at the surface
[1320]94
[1862]95       IF ( .NOT. ALLOCATED( w_subs ) )  THEN
[411]96          ALLOCATE( w_subs(nzb:nzt+1) )
[1353]97          w_subs = 0.0_wp
[411]98       ENDIF
99
[3294]100       IF ( ocean_mode )  THEN
[3302]101          message_string = 'applying large scale vertical motion is not ' //   &
[3294]102                           'allowed for ocean mode'
[411]103          CALL message( 'init_w_subsidence', 'PA0324', 2, 2, 0, 6, 0 )
104       ENDIF
105
106!
107!--   Compute the profile of the subsidence/ascent velocity
108!--   using the given gradients
109      i = 1
[1353]110      gradient = 0.0_wp
111      ws_surface = 0.0_wp
[411]112     
113
[580]114      subs_vertical_gradient_level_i(1) = 0
[411]115      DO  k = 1, nzt+1
[1365]116         IF ( i < 11 )  THEN
[580]117            IF ( subs_vertical_gradient_level(i) < zu(k)  .AND. &
[1353]118                 subs_vertical_gradient_level(i) >= 0.0_wp )  THEN
119               gradient = subs_vertical_gradient(i) / 100.0_wp
[580]120               subs_vertical_gradient_level_i(i) = k - 1
[411]121               i = i + 1
122            ENDIF
123         ENDIF
[1353]124         IF ( gradient /= 0.0_wp )  THEN
[411]125            IF ( k /= 1 )  THEN
126               w_subs(k) = w_subs(k-1) + dzu(k) * gradient
127            ELSE
[1353]128               w_subs(k) = ws_surface   + 0.5_wp * dzu(k) * gradient
[411]129            ENDIF
130         ELSE
131            w_subs(k) = w_subs(k-1)
132         ENDIF
133      ENDDO
134
135!
136!--   In case of no given gradients for the subsidence/ascent velocity,
137!--   choose zero gradient
[1353]138      IF ( subs_vertical_gradient_level(1) == -9999999.9_wp )  THEN
139         subs_vertical_gradient_level(1) = 0.0_wp
[411]140      ENDIF
141
142    END SUBROUTINE init_w_subsidence
143
144
[1682]145!------------------------------------------------------------------------------!
146! Description:
147! ------------
[3554]148!> Add effect of large-scale subsidence to variable.
[1682]149!------------------------------------------------------------------------------!
[1365]150    SUBROUTINE subsidence( tendency, var, var_init, ls_index ) 
[411]151
[1320]152       USE arrays_3d,                                                          &
153           ONLY:  ddzu, w_subs
[411]154
[1320]155       USE control_parameters,                                                 &
[1380]156           ONLY:  dt_3d, intermediate_timestep_count, large_scale_forcing,     &
157                  scalar_rayleigh_damping
[1320]158
159       USE indices,                                                            &
[2232]160           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,        &
[4346]161                  wall_flags_total_0
[1320]162
163       USE kinds
164
[1365]165       USE statistics,                                                         &
[1382]166           ONLY:  sums_ls_l, weight_substep
[1365]167
[411]168       IMPLICIT NONE
169 
[3554]170       INTEGER(iwp) ::  i        !< loop index
171       INTEGER(iwp) ::  j        !< loop index
172       INTEGER(iwp) ::  k        !< loop index
173       INTEGER(iwp) ::  ls_index !< index of large-scale subsidence in sums_ls_l
[411]174
[3554]175       REAL(wp)     ::  tmp_tend !< temporary tendency
176       REAL(wp)     ::  tmp_grad !< temporary gradient
[411]177   
[3554]178       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var      !< variable where to add subsidence
179       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  tendency !< tendency of var
180       REAL(wp), DIMENSION(nzb:nzt+1) ::  var_init                     !< initialization profile of var
181       REAL(wp), DIMENSION(nzb:nzt+1) ::  var_mod                      !< modified profile of var
[411]182
[464]183       var_mod = var_init
[411]184
185!
186!--    Influence of w_subsidence on the current tendency term
187       DO  i = nxl, nxr
188          DO  j = nys, nyn
[1382]189
[2232]190             DO  k = nzb+1, nzt 
[1365]191                IF ( w_subs(k) < 0.0_wp )  THEN    ! large-scale subsidence
192                   tmp_tend = - w_subs(k) *                                    &
[2232]193                              ( var(k+1,j,i) - var(k,j,i) ) * ddzu(k+1) *      &
194                                        MERGE( 1.0_wp, 0.0_wp,                 &
[4346]195                                        BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1365]196                ELSE                               ! large-scale ascent
197                   tmp_tend = - w_subs(k) *                                    &
[2232]198                              ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k) *        &
199                                        MERGE( 1.0_wp, 0.0_wp,                 &
[4346]200                                        BTEST( wall_flags_total_0(k,j,i), 0 ) )
[411]201                ENDIF
[1365]202
203                tendency(k,j,i) = tendency(k,j,i) + tmp_tend
204
205                IF ( large_scale_forcing )  THEN
206                   sums_ls_l(k,ls_index) = sums_ls_l(k,ls_index) + tmp_tend    &
[2232]207                                 * weight_substep(intermediate_timestep_count) &
208                                 * MERGE( 1.0_wp, 0.0_wp,                      &
[4346]209                                   BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1365]210                ENDIF
[411]211             ENDDO
[1382]212
[1862]213             IF ( large_scale_forcing )  THEN
214                sums_ls_l(nzt+1,ls_index) = sums_ls_l(nzt,ls_index)
215             ENDIF
[1382]216
[411]217          ENDDO
218       ENDDO
219
220!
221!--    Shifting of the initial profile is especially necessary with Rayleigh
222!--    damping switched on
[1729]223       IF ( scalar_rayleigh_damping .AND.                                      &
224            intermediate_timestep_count == 1 )  THEN
[1380]225          DO  k = nzb, nzt
226             IF ( w_subs(k) < 0.0_wp )  THEN      ! large-scale subsidence
227                var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
228                                  ( var_init(k+1) - var_init(k) ) * ddzu(k+1)
229             ENDIF
230          ENDDO
[411]231!
[1380]232!--      At the upper boundary, the initial profile is shifted with aid of
233!--      the gradient tmp_grad. (This is ok if the gradients are linear.)
234         IF ( w_subs(nzt) < 0.0_wp )  THEN
235            tmp_grad = ( var_init(nzt+1) - var_init(nzt) ) * ddzu(nzt+1)
236            var_mod(nzt+1) = var_init(nzt+1) -  &
237                                 dt_3d * w_subs(nzt+1) * tmp_grad
238         ENDIF
[411]239       
240
[1380]241         DO  k = nzt+1, nzb+1, -1
242            IF ( w_subs(k) >= 0.0_wp )  THEN  ! large-scale ascent
243               var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
244                                  ( var_init(k) - var_init(k-1) ) * ddzu(k) 
245            ENDIF
246         ENDDO
247!
248!--      At the lower boundary shifting is not necessary because the
249!--      subsidence velocity w_subs(nzb) vanishes.
250         IF ( w_subs(nzb+1) >= 0.0_wp )  THEN
251            var_mod(nzb) = var_init(nzb)
[411]252         ENDIF
[1380]253
254         var_init = var_mod
[411]255      ENDIF
256
257
258 END SUBROUTINE subsidence
259
[1682]260!------------------------------------------------------------------------------!
261! Description:
262! ------------
[3554]263!> Add effect of large-scale subsidence to variable.
[1682]264!------------------------------------------------------------------------------!
[1365]265 SUBROUTINE subsidence_ij( i, j, tendency, var, var_init, ls_index ) 
[411]266
[1320]267       USE arrays_3d,                                                          &
268           ONLY:  ddzu, w_subs
[411]269
[1320]270       USE control_parameters,                                                 &
[1380]271           ONLY:  dt_3d, intermediate_timestep_count, large_scale_forcing,     &
272                  scalar_rayleigh_damping
[1320]273
274       USE indices,                                                            &
[4346]275           ONLY:  nxl, nxlg, nxrg, nyng, nys, nysg, nzb, nzt,                  &
276                  wall_flags_total_0
[1320]277
278       USE kinds
279
[1365]280       USE statistics,                                                         &
[1382]281           ONLY:  sums_ls_l, weight_substep
[1365]282
[411]283       IMPLICIT NONE
284 
[3554]285       INTEGER(iwp) ::  i        !< loop variable
286       INTEGER(iwp) ::  j        !< loop variable
287       INTEGER(iwp) ::  k        !< loop variable
288       INTEGER(iwp) ::  ls_index !< index of large-scale subsidence in sums_ls_l
[411]289
[3554]290       REAL(wp)     ::  tmp_tend !< temporary tendency
291       REAL(wp)     ::  tmp_grad !< temporary gradient
[411]292   
[3554]293       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var      !< variable where to add subsidence
294       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  tendency !< tendency of var
295       REAL(wp), DIMENSION(nzb:nzt+1) ::  var_init                     !< initialization profile of var
296       REAL(wp), DIMENSION(nzb:nzt+1) ::  var_mod                      !< modified profile of var
[411]297
[464]298       var_mod = var_init
[411]299
300!
301!--    Influence of w_subsidence on the current tendency term
[2232]302       DO  k = nzb+1, nzt 
[1365]303          IF ( w_subs(k) < 0.0_wp )  THEN      ! large-scale subsidence
[2232]304             tmp_tend = - w_subs(k) * ( var(k+1,j,i) - var(k,j,i) )            &
305                                    * ddzu(k+1)                                &
306                                    * MERGE( 1.0_wp, 0.0_wp,                   &
[4346]307                                      BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1365]308          ELSE                                 ! large-scale ascent
[2232]309             tmp_tend = - w_subs(k) * ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k)  &
310                                    * MERGE( 1.0_wp, 0.0_wp,                   &
[4346]311                                      BTEST( wall_flags_total_0(k,j,i), 0 ) )
[411]312          ENDIF
[1365]313
314          tendency(k,j,i) = tendency(k,j,i) + tmp_tend
315
316          IF ( large_scale_forcing )  THEN
317             sums_ls_l(k,ls_index) = sums_ls_l(k,ls_index) + tmp_tend          &
[2232]318                                  * weight_substep(intermediate_timestep_count)&
319                                  * MERGE( 1.0_wp, 0.0_wp,                     &
[4346]320                                    BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1365]321          ENDIF
[411]322       ENDDO
323
[1489]324       IF ( large_scale_forcing )  THEN
325          sums_ls_l(nzt+1,ls_index) = sums_ls_l(nzt,ls_index)
326       ENDIF
[411]327
328!
329!--    Shifting of the initial profile is especially necessary with Rayleigh
330!--    damping switched on
[1729]331       IF ( scalar_rayleigh_damping .AND.                                      &
332            intermediate_timestep_count == 1 )  THEN
[1380]333          IF ( i == nxl .AND. j == nys )  THEN ! shifting only once per PE
[411]334
[1380]335             DO  k = nzb, nzt
336                IF ( w_subs(k) < 0.0_wp )  THEN      ! large-scale subsidence
337                   var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
338                                     ( var_init(k+1) - var_init(k) ) * ddzu(k+1)
339                ENDIF
340             ENDDO
341!
342!--          At the upper boundary, the initial profile is shifted with aid of
343!--          the gradient tmp_grad. (This is ok if the gradients are linear.)
344             IF ( w_subs(nzt) < 0.0_wp )  THEN
345                tmp_grad = ( var_init(nzt+1) - var_init(nzt) ) * ddzu(nzt+1)
346                var_mod(nzt+1) = var_init(nzt+1) -  &
347                                     dt_3d * w_subs(nzt+1) * tmp_grad
[411]348             ENDIF
349       
350
[1380]351             DO  k = nzt+1, nzb+1, -1
352                IF ( w_subs(k) >= 0.0_wp )  THEN  ! large-scale ascent
353                   var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
354                                      ( var_init(k) - var_init(k-1) ) * ddzu(k)
355                ENDIF
356             ENDDO
357!
358!--          At the lower boundary shifting is not necessary because the
359!--          subsidence velocity w_subs(nzb) vanishes.
360             IF ( w_subs(nzb+1) >= 0.0_wp )  THEN
361                var_mod(nzb) = var_init(nzb)
[411]362             ENDIF
363
[1380]364             var_init = var_mod 
[411]365
[1380]366          ENDIF
[411]367       ENDIF
368
369 END SUBROUTINE subsidence_ij
370
371
372 END MODULE subsidence_mod
Note: See TracBrowser for help on using the repository browser.