[1682] | 1 | !> @file disturb_field.f90 |
---|
[4583] | 2 | !--------------------------------------------------------------------------------------------------! |
---|
[2696] | 3 | ! This file is part of the PALM model system. |
---|
[1036] | 4 | ! |
---|
[4583] | 5 | ! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General |
---|
| 6 | ! Public License as published by the Free Software Foundation, either version 3 of the License, or |
---|
| 7 | ! (at your option) any later version. |
---|
[1036] | 8 | ! |
---|
[4583] | 9 | ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the |
---|
| 10 | ! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General |
---|
| 11 | ! Public License for more details. |
---|
[1036] | 12 | ! |
---|
[4583] | 13 | ! You should have received a copy of the GNU General Public License along with PALM. If not, see |
---|
| 14 | ! <http://www.gnu.org/licenses/>. |
---|
[1036] | 15 | ! |
---|
[4360] | 16 | ! Copyright 1997-2020 Leibniz Universitaet Hannover |
---|
[4583] | 17 | !--------------------------------------------------------------------------------------------------! |
---|
[1036] | 18 | ! |
---|
[484] | 19 | ! Current revisions: |
---|
[1318] | 20 | ! ------------------ |
---|
[1354] | 21 | ! |
---|
[2233] | 22 | ! |
---|
[1321] | 23 | ! Former revisions: |
---|
| 24 | ! ----------------- |
---|
| 25 | ! $Id: disturb_field.f90 4583 2020-06-29 12:36:47Z oliver.maas $ |
---|
[4583] | 26 | ! file re-formatted to follow the PALM coding standard |
---|
| 27 | ! |
---|
| 28 | ! 4457 2020-03-11 14:20:43Z raasch |
---|
[4457] | 29 | ! use statement for exchange horiz added |
---|
[4583] | 30 | ! |
---|
[4457] | 31 | ! 4360 2020-01-07 11:25:50Z suehring |
---|
[4583] | 32 | ! Introduction of wall_flags_total_0, which currently sets bits based on static topography |
---|
| 33 | ! information used in wall_flags_static_0 |
---|
| 34 | ! |
---|
[4346] | 35 | ! 4329 2019-12-10 15:46:36Z motisi |
---|
[4329] | 36 | ! Renamed wall_flags_0 to wall_flags_static_0 |
---|
[4583] | 37 | ! |
---|
[4329] | 38 | ! 4237 2019-09-25 11:33:42Z knoop |
---|
[4237] | 39 | ! Added missing OpenMP directives |
---|
[4583] | 40 | ! |
---|
[4237] | 41 | ! 4182 2019-08-22 15:20:23Z scharf |
---|
[2716] | 42 | ! Corrected "Former revisions" section |
---|
[4583] | 43 | ! |
---|
[4182] | 44 | ! 3849 2019-04-01 16:35:16Z knoop |
---|
| 45 | ! Corrected "Former revisions" section |
---|
[2716] | 46 | ! |
---|
[4182] | 47 | ! Revision 1.1 1998/02/04 15:40:45 raasch |
---|
| 48 | ! Initial revision |
---|
| 49 | ! |
---|
| 50 | ! |
---|
[1] | 51 | ! Description: |
---|
| 52 | ! ------------ |
---|
[1682] | 53 | !> Imposing a random perturbation on a 3D-array. |
---|
[4583] | 54 | !> On parallel computers, the random number generator is as well called for all gridpoints of the |
---|
| 55 | !> total domain to ensure, regardless of the number of PEs used, that the elements of the array have |
---|
| 56 | !> the same values in the same order in every case. The perturbation range is steered by dist_range. |
---|
| 57 | !--------------------------------------------------------------------------------------------------! |
---|
[2232] | 58 | SUBROUTINE disturb_field( var_char, dist1, field ) |
---|
[1] | 59 | |
---|
[4583] | 60 | |
---|
| 61 | USE control_parameters, & |
---|
| 62 | ONLY: dist_nxl, dist_nxr, dist_nyn, dist_nys, dist_range, disturbance_amplitude, & |
---|
| 63 | disturbance_created, disturbance_level_ind_b, disturbance_level_ind_t, iran, & |
---|
[1320] | 64 | random_generator, topography |
---|
[4583] | 65 | |
---|
| 66 | USE cpulog, & |
---|
[1320] | 67 | ONLY: cpu_log, log_point |
---|
[4583] | 68 | |
---|
| 69 | USE exchange_horiz_mod, & |
---|
[4457] | 70 | ONLY: exchange_horiz |
---|
| 71 | |
---|
[4583] | 72 | USE indices, & |
---|
| 73 | ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, nzt, & |
---|
| 74 | wall_flags_total_0 |
---|
| 75 | |
---|
[1320] | 76 | USE kinds |
---|
[4583] | 77 | |
---|
| 78 | USE random_function_mod, & |
---|
[1320] | 79 | ONLY: random_function |
---|
[1] | 80 | |
---|
[4583] | 81 | USE random_generator_parallel, & |
---|
| 82 | ONLY: random_number_parallel, random_seed_parallel, random_dummy, seq_random_array |
---|
| 83 | |
---|
[1] | 84 | IMPLICIT NONE |
---|
| 85 | |
---|
[2232] | 86 | CHARACTER (LEN = *) :: var_char !< flag to distinguish betwenn u- and v-component |
---|
[1] | 87 | |
---|
[4583] | 88 | INTEGER(iwp) :: flag_nr !< number of respective flag for u- or v-grid |
---|
[2232] | 89 | INTEGER(iwp) :: i !< index variable |
---|
| 90 | INTEGER(iwp) :: j !< index variable |
---|
| 91 | INTEGER(iwp) :: k !< index variable |
---|
| 92 | |
---|
[1682] | 93 | REAL(wp) :: randomnumber !< |
---|
[4583] | 94 | |
---|
[1682] | 95 | REAL(wp) :: dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !< |
---|
| 96 | REAL(wp) :: field(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !< |
---|
[4583] | 97 | |
---|
[1682] | 98 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dist2 !< |
---|
[1] | 99 | |
---|
| 100 | |
---|
| 101 | CALL cpu_log( log_point(20), 'disturb_field', 'start' ) |
---|
| 102 | ! |
---|
[2232] | 103 | !-- Set flag number, 20 for u-grid, 21 for v-grid, required to mask topography |
---|
| 104 | flag_nr = MERGE( 20, 21, TRIM(var_char) == 'u' ) |
---|
| 105 | ! |
---|
[4583] | 106 | !-- Create an additional temporary array and initialize the arrays needed to store the disturbance |
---|
[667] | 107 | ALLOCATE( dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
[3849] | 108 | !$ACC DATA CREATE(dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) |
---|
| 109 | |
---|
| 110 | ! |
---|
| 111 | !-- dist1 is initialized on the host (see below) and then updated on the device. |
---|
[1353] | 112 | dist1 = 0.0_wp |
---|
[3849] | 113 | !$ACC KERNELS PRESENT(dist2) |
---|
[1353] | 114 | dist2 = 0.0_wp |
---|
[3849] | 115 | !$ACC END KERNELS |
---|
[1] | 116 | |
---|
| 117 | ! |
---|
| 118 | !-- Create the random perturbation and store it on temporary array |
---|
| 119 | IF ( random_generator == 'numerical-recipes' ) THEN |
---|
| 120 | DO i = dist_nxl(dist_range), dist_nxr(dist_range) |
---|
| 121 | DO j = dist_nys(dist_range), dist_nyn(dist_range) |
---|
| 122 | DO k = disturbance_level_ind_b, disturbance_level_ind_t |
---|
[4583] | 123 | randomnumber = 3.0_wp * disturbance_amplitude * ( random_function( iran ) - 0.5_wp ) |
---|
| 124 | IF ( nxl <= i .AND. nxr >= i .AND. nys <= j .AND. nyn >= j ) THEN |
---|
[1] | 125 | dist1(k,j,i) = randomnumber |
---|
| 126 | ENDIF |
---|
| 127 | ENDDO |
---|
| 128 | ENDDO |
---|
| 129 | ENDDO |
---|
[1400] | 130 | ELSEIF ( random_generator == 'random-parallel' ) THEN |
---|
| 131 | DO i = dist_nxl(dist_range), dist_nxr(dist_range) |
---|
| 132 | DO j = dist_nys(dist_range), dist_nyn(dist_range) |
---|
| 133 | CALL random_seed_parallel( put=seq_random_array(:, j, i) ) |
---|
| 134 | DO k = disturbance_level_ind_b, disturbance_level_ind_t |
---|
| 135 | CALL random_number_parallel( random_dummy ) |
---|
[4583] | 136 | randomnumber = 3.0_wp * disturbance_amplitude * ( random_dummy - 0.5_wp ) |
---|
| 137 | IF ( nxl <= i .AND. nxr >= i .AND. nys <= j .AND. nyn >= j ) THEN |
---|
[1400] | 138 | dist1(k,j,i) = randomnumber |
---|
| 139 | ENDIF |
---|
| 140 | ENDDO |
---|
| 141 | CALL random_seed_parallel( get=seq_random_array(:, j, i) ) |
---|
| 142 | ENDDO |
---|
| 143 | ENDDO |
---|
[1] | 144 | ELSEIF ( random_generator == 'system-specific' ) THEN |
---|
| 145 | DO i = dist_nxl(dist_range), dist_nxr(dist_range) |
---|
| 146 | DO j = dist_nys(dist_range), dist_nyn(dist_range) |
---|
| 147 | DO k = disturbance_level_ind_b, disturbance_level_ind_t |
---|
| 148 | CALL RANDOM_NUMBER( randomnumber ) |
---|
[4583] | 149 | randomnumber = 3.0_wp * disturbance_amplitude * ( randomnumber - 0.5_wp ) |
---|
| 150 | IF ( nxl <= i .AND. nxr >= i .AND. nys <= j .AND. nyn >= j ) THEN |
---|
[1] | 151 | dist1(k,j,i) = randomnumber |
---|
| 152 | ENDIF |
---|
| 153 | ENDDO |
---|
| 154 | ENDDO |
---|
| 155 | ENDDO |
---|
| 156 | |
---|
| 157 | ENDIF |
---|
| 158 | |
---|
| 159 | ! |
---|
[3849] | 160 | !-- Update dist1 on the device, this is expected by exchange_horiz! |
---|
| 161 | !$ACC UPDATE DEVICE(dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) |
---|
| 162 | |
---|
| 163 | ! |
---|
[1] | 164 | !-- Exchange of ghost points for the random perturbation |
---|
[420] | 165 | |
---|
[667] | 166 | CALL exchange_horiz( dist1, nbgp ) |
---|
[1] | 167 | ! |
---|
| 168 | !-- Applying the Shuman filter in order to smooth the perturbations. |
---|
[4583] | 169 | !-- Neighboured grid points in all three directions are used for the filter operation. |
---|
| 170 | !-- Loop has been splitted to make runs reproducible on HLRN systems using compiler option -O3 |
---|
[3849] | 171 | !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) PRESENT(dist1, dist2) |
---|
[4237] | 172 | !$OMP PARALLEL DO PRIVATE(i, j, k) |
---|
[420] | 173 | DO i = nxl, nxr |
---|
| 174 | DO j = nys, nyn |
---|
[1] | 175 | DO k = disturbance_level_ind_b-1, disturbance_level_ind_t+1 |
---|
[4583] | 176 | dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1) & |
---|
| 177 | + dist1(k,j+1,i) + dist1(k+1,j,i) & |
---|
[1322] | 178 | ) / 12.0_wp |
---|
[420] | 179 | ENDDO |
---|
| 180 | DO k = disturbance_level_ind_b-1, disturbance_level_ind_t+1 |
---|
[4583] | 181 | dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i) & |
---|
| 182 | + 6.0_wp * dist1(k,j,i) & |
---|
| 183 | ) / 12.0_wp |
---|
[1] | 184 | ENDDO |
---|
[420] | 185 | ENDDO |
---|
| 186 | ENDDO |
---|
[1] | 187 | |
---|
| 188 | ! |
---|
| 189 | !-- Exchange of ghost points for the filtered perturbation. |
---|
| 190 | !-- Afterwards, filter operation and exchange of ghost points are repeated. |
---|
[667] | 191 | CALL exchange_horiz( dist2, nbgp ) |
---|
[420] | 192 | |
---|
[3849] | 193 | !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) PRESENT(dist1, dist2) |
---|
[4237] | 194 | !$OMP PARALLEL DO PRIVATE(i, j, k) |
---|
[1] | 195 | DO i = nxl, nxr |
---|
| 196 | DO j = nys, nyn |
---|
| 197 | DO k = disturbance_level_ind_b-2, disturbance_level_ind_t+2 |
---|
[4583] | 198 | dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i) & |
---|
| 199 | + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i) & |
---|
| 200 | + 6.0_wp * dist2(k,j,i) & |
---|
[1322] | 201 | ) / 12.0_wp |
---|
[1] | 202 | ENDDO |
---|
| 203 | ENDDO |
---|
| 204 | ENDDO |
---|
[420] | 205 | |
---|
[667] | 206 | CALL exchange_horiz( dist1, nbgp ) |
---|
[1] | 207 | |
---|
| 208 | ! |
---|
[4583] | 209 | !-- Remove perturbations below topography (including one gridpoint above it in order to allow for |
---|
| 210 | !-- larger timesteps at the beginning of the simulation (diffusion criterion)) |
---|
[1] | 211 | IF ( TRIM( topography ) /= 'flat' ) THEN |
---|
[667] | 212 | DO i = nxlg, nxrg |
---|
| 213 | DO j = nysg, nyng |
---|
[2232] | 214 | DO k = nzb, nzb_max |
---|
[4583] | 215 | dist1(k,j,i) = MERGE( dist1(k,j,i), 0.0_wp, & |
---|
| 216 | BTEST( wall_flags_total_0(k,j,i), flag_nr ) & |
---|
[2232] | 217 | ) |
---|
| 218 | ENDDO |
---|
[1] | 219 | ENDDO |
---|
| 220 | ENDDO |
---|
| 221 | ENDIF |
---|
| 222 | |
---|
| 223 | ! |
---|
| 224 | !-- Random perturbation is added to the array to be disturbed. |
---|
[3849] | 225 | !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) PRESENT(field, dist1) |
---|
[4237] | 226 | !$OMP PARALLEL DO PRIVATE(i, j, k) |
---|
[667] | 227 | DO i = nxlg, nxrg |
---|
| 228 | DO j = nysg, nyng |
---|
[1] | 229 | DO k = disturbance_level_ind_b-2, disturbance_level_ind_t+2 |
---|
| 230 | field(k,j,i) = field(k,j,i) + dist1(k,j,i) |
---|
| 231 | ENDDO |
---|
| 232 | ENDDO |
---|
| 233 | ENDDO |
---|
| 234 | |
---|
[3849] | 235 | !$ACC END DATA |
---|
| 236 | |
---|
[1] | 237 | ! |
---|
| 238 | !-- Deallocate the temporary array |
---|
| 239 | DEALLOCATE( dist2 ) |
---|
| 240 | |
---|
| 241 | ! |
---|
| 242 | !-- Set a flag, which indicates that a random perturbation is imposed |
---|
| 243 | disturbance_created = .TRUE. |
---|
| 244 | |
---|
| 245 | |
---|
| 246 | CALL cpu_log( log_point(20), 'disturb_field', 'stop' ) |
---|
| 247 | |
---|
| 248 | |
---|
| 249 | END SUBROUTINE disturb_field |
---|