source: palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90 @ 3994

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

new module for diagnostic output quantities added + output of turbulence intensity

  • Property svn:keywords set to Id
File size: 11.6 KB
RevLine 
[3994]1!> @file diagnostic_output_quantities_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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: diagnostic_output_quantities_mod.f90 3994 2019-05-22 18:08:09Z suehring $
27! Initial revision
28!
29!
30! @author Farah Kanani-Suehring
31!
32! Description:
33! ------------
34!> ...
35!------------------------------------------------------------------------------!
36 MODULE diagnostic_output_quantities_mod
37 
38
39!     USE arrays_3d,                                                             &
40!         ONLY:  dzw, e, heatflux_output_conversion, nc, nr, p, pt,              &
41!                precipitation_amount, prr, q, qc, ql, ql_c, ql_v, qr, s, tend,  &
42!                u, v, vpt, w, zu, zw, waterflux_output_conversion, hyrho, d_exner
43!
44!     USE averaging
45!
46!     USE basic_constants_and_equations_mod,                                     &
47!         ONLY:  c_p, lv_d_cp, l_v
48!
49!     USE bulk_cloud_model_mod,                                                  &
50!         ONLY:  bulk_cloud_model
51!
52    USE control_parameters,                                                    &
53        ONLY:  current_timestep_number, varnamelength
54!
55!     USE cpulog,                                                                &
56!         ONLY:  cpu_log, log_point
57!
58!     USE indices,                                                               &
59!         ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,       &
60!                nzb, nzt, wall_flags_0
61!
62    USE kinds
63!
64!     USE land_surface_model_mod,                                                &
65!         ONLY:  zs
66!
67!     USE module_interface,                                                      &
68!         ONLY:  module_interface_data_output_2d
69!
70! #if defined( __netcdf )
71!     USE NETCDF
72! #endif
73!
74!     USE netcdf_interface,                                                      &
75!         ONLY:  fill_value, id_set_xy, id_set_xz, id_set_yz, id_var_do2d,       &
76!                id_var_time_xy, id_var_time_xz, id_var_time_yz, nc_stat,        &
77!                netcdf_data_format, netcdf_handle_error
78!
79!     USE particle_attributes,                                                   &
80!         ONLY:  grid_particles, number_of_particles, particle_advection_start,  &
81!                particles, prt_count
82!     
83!     USE pegrid
84!
85!     USE surface_mod,                                                           &
86!         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h,           &
87!                surf_lsm_h, surf_usm_h
88!
89!     USE turbulence_closure_mod,                                                &
90!         ONLY:  tcm_data_output_2d
91
92
93    IMPLICIT NONE
94
95    CHARACTER(LEN=varnamelength), DIMENSION(500) ::  do_all = ' '
96
97    INTEGER(iwp) ::  timestep_number_at_prev_calc = 0  !< ...at previous diagnostic output calculation
98
99    LOGICAL ::  initialized_diagnostic_output_quantities = .FALSE.
100    LOGICAL ::  prepared_diagnostic_output_quantities = .FALSE.
101
102    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti     !< rotation(u,v,w) aka turbulence intensity
103    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti_av  !< avg. rotation(u,v,w) aka turbulence intensity
104
105
106    SAVE
107
108    PRIVATE
109
110!
111!-- Public variables
112    PUBLIC do_all,                                                                                 &
113           initialized_diagnostic_output_quantities,                                               &
114           prepared_diagnostic_output_quantities,                                                  &
115           ti, ti_av,                                                                              &
116           timestep_number_at_prev_calc
117!
118!-- Public routines
119    PUBLIC diagnostic_output_quantities_calculate,                                                 &
120           diagnostic_output_quantities_init,                                                      &
121           diagnostic_output_quantities_prepare
122
123
124    INTERFACE diagnostic_output_quantities_init
125       MODULE PROCEDURE diagnostic_output_quantities_init
126    END INTERFACE diagnostic_output_quantities_init
127
128    INTERFACE diagnostic_output_quantities_calculate
129       MODULE PROCEDURE diagnostic_output_quantities_calculate
130    END INTERFACE diagnostic_output_quantities_calculate
131
132    INTERFACE diagnostic_output_quantities_prepare
133       MODULE PROCEDURE diagnostic_output_quantities_prepare
134    END INTERFACE diagnostic_output_quantities_prepare
135
136
137 CONTAINS
138
139!------------------------------------------------------------------------------!
140! Description:
141! ------------
142!> ...
143!------------------------------------------------------------------------------!
144 SUBROUTINE diagnostic_output_quantities_init
145
146
147    USE indices,                                                                                   &
148        ONLY:  nxl, nxr, nyn, nys, nzb, nzt
149
150    IMPLICIT NONE
151
152
153    initialized_diagnostic_output_quantities = .FALSE.
154
155    IF ( .NOT. ALLOCATED( ti ) )  THEN
156       ALLOCATE( ti(nzb:nzt+1,nys:nyn,nxl:nxr) )
157       ti = 0.0_wp
158    ENDIF
159
160    initialized_diagnostic_output_quantities = .TRUE.
161
162 END SUBROUTINE diagnostic_output_quantities_init
163
164
165!--------------------------------------------------------------------------------------------------!
166! Description:
167! ------------
168!> ...
169!--------------------------------------------------------------------------------------------------!
170 SUBROUTINE diagnostic_output_quantities_calculate
171
172
173    USE arrays_3d,                                                                                 &
174        ONLY:  ddzu, u, v, w
175
176    USE grid_variables,                                                                            &
177        ONLY:  ddx, ddy
178
179    USE indices,                                                                                   &
180        ONLY:  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
181
182    IMPLICIT NONE
183
184    INTEGER(iwp) ::  i, j, k    !< grid loop index in x-, y-, z-direction
185    INTEGER(iwp) ::  ivar_all   !< loop index over all 2d/3d/mask output quantities
186
187
188!     CALL cpu_log( log_point(41), 'calculate_quantities', 'start' )
189!
190!-- Preparatory steps and initialization of output arrays
191    IF ( .NOT.  prepared_diagnostic_output_quantities )     CALL diagnostic_output_quantities_prepare
192    IF ( .NOT.  initialized_diagnostic_output_quantities )  CALL diagnostic_output_quantities_init
193!
194!-- Save timestep number to check in time_integration if diagnostic_output_quantities_calculate
195!-- has been called already, since the CALL occurs at two locations, but the calculations need to be
196!-- done only once per timestep.
197    timestep_number_at_prev_calc = current_timestep_number
198
199
200    ivar_all = 1
201
202    DO  WHILE ( do_all(ivar_all)(1:1) /= ' ' )
203
204       SELECT CASE ( TRIM( do_all(ivar_all) ) )
205!
206!--       Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point
207          CASE ( 'ti' )
208             DO  i = nxl, nxr
209                DO  j = nys, nyn
210                   DO  k = nzb+1, nzt
211
212                      ti(k,j,i) = 0.25_wp * SQRT(                                                  &
213                        (   (   w(k,j+1,i) + w(k-1,j+1,i)                                          &
214                              - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy                                  &
215                          - (   v(k+1,j,i) + v(k+1,j+1,i)                                          &
216                              - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2                         &
217                      + (   (   u(k+1,j,i) + u(k+1,j,i+1)                                          &
218                              - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k)                              &
219                          - (   w(k,j,i+1) + w(k-1,j,i+1)                                          &
220                              - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx     )**2                         &
221                      + (   (   v(k,j,i+1) + v(k,j+1,i+1)                                          &
222                              - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx                                  &
223                          - (   u(k,j+1,i) + u(k,j+1,i+1)                                          &
224                              - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy     )**2  )                      &
225                                          * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0) )
226                   ENDDO
227                ENDDO
228             ENDDO
229
230          END SELECT
231
232          ivar_all = ivar_all + 1
233    ENDDO
234
235!     CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' )
236
237 END SUBROUTINE diagnostic_output_quantities_calculate
238
239
240!------------------------------------------------------------------------------!
241! Description:
242! ------------
243!> ...
244!------------------------------------------------------------------------------!
245 SUBROUTINE diagnostic_output_quantities_prepare
246
247
248    USE control_parameters,                                                                        &
249        ONLY:  do2d, do3d, domask, masks, mid
250
251    IMPLICIT NONE
252
253    CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) ::  do2d_var = ' '  !<
254                                                          !< label array for 2d output quantities
255
256    INTEGER(iwp) ::  av         !< index defining type of output, av=0 instantaneous, av=1 averaged
257    INTEGER(iwp) ::  ivar       !< loop index
258    INTEGER(iwp) ::  ivar_all   !< loop index
259    INTEGER(iwp) ::  l          !< index for cutting string
260
261
262    prepared_diagnostic_output_quantities = .FALSE.
263
264    ivar     = 1
265    ivar_all = 1
266
267    DO  av = 0, 1
268!
269!--    Remove _xy, _xz, or _yz from string
270       l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
271       do2d_var(av,ivar) = do2d(av,ivar)(0:l-2)
272!
273!--    Gather 2d output quantity names, check for double occurrence of output quantity
274       DO  WHILE ( do2d_var(av,ivar)(1:1) /= ' ' )
275          IF ( .NOT.  ANY( do_all == do2d_var(av,ivar) ) )  THEN
276             do_all(ivar_all) = do2d_var(av,ivar)
277          ENDIF
278          ivar = ivar + 1
279          ivar_all = ivar_all + 1
280          l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
281          do2d_var(av,ivar) = do2d(av,ivar)(0:l-2)
282       ENDDO
283
284       ivar = 1
285!
286!--    Gather 3d output quantity names, check for double occurrence of output quantity
287       DO  WHILE ( do3d(av,ivar)(1:1) /= ' ' )
288          do_all(ivar_all) = do3d(av,ivar)
289
290          ivar = ivar + 1
291          ivar_all = ivar_all + 1
292       ENDDO
293
294       ivar = 1
295!
296!--    Gather masked output quantity names, check for double occurrence of output quantity
297       DO  mid = 1, masks
298          DO  WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
299             do_all(ivar_all) = domask(mid,av,ivar)
300
301             ivar = ivar + 1
302             ivar_all = ivar_all + 1
303          ENDDO
304          ivar = 1
305       ENDDO
306
307    ENDDO
308
309    prepared_diagnostic_output_quantities = .TRUE.
310
311 END SUBROUTINE diagnostic_output_quantities_prepare
312
313 END MODULE diagnostic_output_quantities_mod
Note: See TracBrowser for help on using the repository browser.