- Timestamp:
- Mar 14, 2014 6:25:17 PM (11 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
lpm_advec.f90 (modified) (3 diffs)
-
lpm_init.f90 (modified) (3 diffs)
-
modules.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_advec.f90
r1310 r1314 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Vertical logarithmic interpolation of horizontal particle speed for particles 23 ! between roughness height and first vertical grid level. 23 24 ! 24 25 ! Former revisions: … … 65 66 e_mean_int, fs_int, lagr_timescale, random_gauss, vv_int 66 67 68 REAL :: height_int, height_p, log_z_z0_int, us_int, z_p, d_z_p_z0 69 67 70 REAL :: location(1:30,1:3) 68 71 REAL, DIMENSION(1:30) :: de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei 69 72 73 ! 74 !-- Determine height of Prandtl layer and distance between Prandtl-layer 75 !-- height and horizontal mean roughness height, which are required for 76 !-- vertical logarithmic interpolation of horizontal particle speeds 77 !-- (for particles below first vertical grid level). 78 z_p = zu(nzb+1) - zw(nzb) 79 d_z_p_z0 = 1.0 / ( z_p - z0_av_global ) 70 80 71 81 DO n = 1, number_of_particles … … 75 85 !-- reached 76 86 IF ( ( dt_3d - particles(n)%dt_sum ) < 1E-8 ) CYCLE 77 78 ! 79 !-- Interpolate u velocity-component, determine left, front, bottom 80 !-- index of u-array 81 i = ( particles(n)%x + 0.5 * dx ) * ddx 82 j = particles(n)%y * ddy 87 ! 88 !-- Determine bottom index 83 89 k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz & 84 + offset_ocean_nzt ! only exact if equidistant 85 86 ! 87 !-- Interpolation of the velocity components in the xy-plane 88 x = particles(n)%x + ( 0.5 - i ) * dx 89 y = particles(n)%y - j * dy 90 aa = x**2 + y**2 91 bb = ( dx - x )**2 + y**2 92 cc = x**2 + ( dy - y )**2 93 dd = ( dx - x )**2 + ( dy - y )**2 94 gg = aa + bb + cc + dd 95 96 u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) & 97 + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) & 98 ) / ( 3.0 * gg ) - u_gtrans 99 IF ( k+1 == nzt+1 ) THEN 100 u_int = u_int_l 90 + offset_ocean_nzt ! only exact if equidistant 91 ! 92 !-- Interpolation of the u velocity component onto particle position. 93 !-- Particles are interpolation bi-linearly in the horizontal and a 94 !-- linearly in the vertical. An exception is made for particles below 95 !-- the first vertical grid level in case of a prandtl layer. In this 96 !-- case the horizontal particle velocity components are determined using 97 !-- Monin-Obukhov relations (if branch). 98 !-- First, check if particle is located below first vertical grid level 99 !-- (Prandtl-layer height) 100 IF ( prandtl_layer .AND. particles(n)%z < z_p ) THEN 101 ! 102 !-- Resolved-scale horizontal particle velocity is zero below z0. 103 IF ( particles(n)%z < z0_av_global ) THEN 104 105 u_int = 0.0 106 107 ELSE 108 ! 109 !-- Determine the sublayer. Further used as index. 110 height_p = ( particles(n)%z - z0_av_global ) & 111 * REAL( number_of_sublayers ) & 112 * d_z_p_z0 113 ! 114 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 115 !-- interpolate linearly between precalculated logarithm. 116 log_z_z0_int = log_z_z0(INT(height_p)) & 117 + ( height_p - INT(height_p) ) & 118 * ( log_z_z0(INT(height_p)+1) & 119 - log_z_z0(INT(height_p)) & 120 ) 121 ! 122 !-- Neutral solution is applied for all situations, e.g. also for 123 !-- unstable and stable situations. Even though this is not exact 124 !-- this saves a lot of CPU time since several calls of intrinsic 125 !-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified 126 !-- as sensitivity studies revealed no significant effect of 127 !-- using the neutral solution also for un/stable situations. 128 !-- Calculated left and bottom index on u grid. 129 i = ( particles(n)%x + 0.5 * dx ) * ddx 130 j = particles(n)%y * ddy 131 132 us_int = 0.5 * ( us(j,i) + us(j,i-1) ) 133 134 u_int = -usws(j,i) / ( us_int * kappa + 1E-10 ) & 135 * log_z_z0_int 136 137 ENDIF 138 ! 139 !-- Particle above the first grid level. Bi-linear interpolation in the 140 !-- horizontal and linear interpolation in the vertical direction. 101 141 ELSE 102 u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) & 103 + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) & 142 ! 143 !-- Interpolate u velocity-component, determine left, front, bottom 144 !-- index of u-array. Adopt k index from above 145 i = ( particles(n)%x + 0.5 * dx ) * ddx 146 j = particles(n)%y * ddy 147 ! 148 !-- Interpolation of the velocity components in the xy-plane 149 x = particles(n)%x + ( 0.5 - i ) * dx 150 y = particles(n)%y - j * dy 151 aa = x**2 + y**2 152 bb = ( dx - x )**2 + y**2 153 cc = x**2 + ( dy - y )**2 154 dd = ( dx - x )**2 + ( dy - y )**2 155 gg = aa + bb + cc + dd 156 157 u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) & 158 + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1)& 104 159 ) / ( 3.0 * gg ) - u_gtrans 105 u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * & 160 161 IF ( k+1 == nzt+1 ) THEN 162 163 u_int = u_int_l 164 165 ELSE 166 167 u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) & 168 + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) & 169 ) / ( 3.0 * gg ) - u_gtrans 170 171 u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * & 106 172 ( u_int_u - u_int_l ) 173 174 ENDIF 175 107 176 ENDIF 108 177 109 178 ! 110 !-- Same procedure for interpolation of the v velocity-component (adopt 111 !-- index k from u velocity-component) 112 i = particles(n)%x * ddx 113 j = ( particles(n)%y + 0.5 * dy ) * ddy 114 115 x = particles(n)%x - i * dx 116 y = particles(n)%y + ( 0.5 - j ) * dy 117 aa = x**2 + y**2 118 bb = ( dx - x )**2 + y**2 119 cc = x**2 + ( dy - y )**2 120 dd = ( dx - x )**2 + ( dy - y )**2 121 gg = aa + bb + cc + dd 122 123 v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) & 124 + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) & 125 ) / ( 3.0 * gg ) - v_gtrans 126 IF ( k+1 == nzt+1 ) THEN 127 v_int = v_int_l 179 !-- Same procedure for interpolation of the v velocity-component. 180 IF ( prandtl_layer .AND. particles(n)%z < z_p ) THEN 181 ! 182 !-- Resolved-scale horizontal particle velocity is zero below z0. 183 IF ( particles(n)%z < z0_av_global ) THEN 184 185 v_int = 0.0 186 187 ELSE 188 ! 189 !-- Neutral solution is applied for all situations, e.g. also for 190 !-- unstable and stable situations. Even though this is not exact 191 !-- this saves a lot of CPU time since several calls of intrinsic 192 !-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified 193 !-- as sensitivity studies revealed no significant effect of 194 !-- using the neutral solution also for un/stable situations. 195 !-- Calculated left and bottom index on v grid. 196 i = particles(n)%x * ddx 197 j = ( particles(n)%y + 0.5 * dy ) * ddy 198 199 us_int = 0.5 * ( us(j,i) + us(j-1,i) ) 200 201 v_int = -vsws(j,i) / ( us_int * kappa + 1E-10 ) & 202 * log_z_z0_int 203 204 ENDIF 205 ! 206 !-- Particle above the first grid level. Bi-linear interpolation in the 207 !-- horizontal and linear interpolation in the vertical direction. 128 208 ELSE 129 v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) & 130 + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) & 209 i = particles(n)%x * ddx 210 j = ( particles(n)%y + 0.5 * dy ) * ddy 211 x = particles(n)%x - i * dx 212 y = particles(n)%y + ( 0.5 - j ) * dy 213 aa = x**2 + y**2 214 bb = ( dx - x )**2 + y**2 215 cc = x**2 + ( dy - y )**2 216 dd = ( dx - x )**2 + ( dy - y )**2 217 gg = aa + bb + cc + dd 218 219 v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) & 220 + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1)& 131 221 ) / ( 3.0 * gg ) - v_gtrans 132 v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * & 222 IF ( k+1 == nzt+1 ) THEN 223 v_int = v_int_l 224 ELSE 225 v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) & 226 + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) & 227 ) / ( 3.0 * gg ) - v_gtrans 228 v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * & 133 229 ( v_int_u - v_int_l ) 230 ENDIF 231 134 232 ENDIF 135 233 136 234 ! 137 !-- Same procedure for interpolation of the w velocity-component (adopt 138 !-- index i from v velocity-component) 235 !-- Same procedure for interpolation of the w velocity-component 139 236 IF ( vertical_particle_advection(particles(n)%group) ) THEN 237 i = particles(n)%x * ddx 140 238 j = particles(n)%y * ddy 141 239 k = particles(n)%z / dz + offset_ocean_nzt_m1 -
palm/trunk/SOURCE/lpm_init.f90
r1310 r1314 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Vertical logarithmic interpolation of horizontal particle speed for particles 23 ! between roughness height and first vertical grid level. 23 24 ! 24 25 ! Former revisions: … … 110 111 IMPLICIT NONE 111 112 112 INTEGER :: i, j, n, nn113 INTEGER :: i, j, k, n, nn 113 114 #if defined( __parallel ) 114 115 INTEGER, DIMENSION(3) :: blocklengths, displacements, types 115 116 #endif 116 117 LOGICAL :: uniform_particles_l 117 REAL :: pos_x, pos_y, pos_z 118 REAL :: height_int, height_p, pos_x, pos_y, pos_z, z_p, & 119 z0_av_local = 0.0 118 120 119 121 … … 184 186 de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 185 187 de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 188 ENDIF 189 190 ! 191 !-- Allocate array required for logarithmic vertical interpolation of 192 !-- horizontal particle velocities between the surface and the first vertical 193 !-- grid level. In order to avoid repeated CPU cost-intensive CALLS of 194 !-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for 195 !-- several heights. Splitting into 20 sublayers turned out to be sufficient. 196 !-- To obtain exact height levels of particles, linear interpolation is applied 197 !-- (see lpm_advec.f90). 198 IF ( prandtl_layer ) THEN 199 200 ALLOCATE ( log_z_z0(0:number_of_sublayers) ) 201 z_p = zu(nzb+1) - zw(nzb) 202 203 ! 204 !-- Calculate horizontal mean value of z0 used for logartihmic 205 !-- interpolation. Note: this is not exact for heterogeneous z0. 206 !-- However, sensitivity studies showed that the effect is 207 !-- negligible. 208 z0_av_local = SUM( z0(nys:nyn,nxl:nxr) ) 209 z0_av_global = 0.0 210 211 CALL MPI_ALLREDUCE(z0_av_local, z0_av_global, 1, MPI_REAL, MPI_SUM, & 212 comm2d, ierr ) 213 214 z0_av_global = z0_av_global / ( ( ny + 1 ) * ( nx + 1 ) ) 215 ! 216 !-- Horizontal wind speed is zero below and at z0 217 log_z_z0(0) = 0.0 218 ! 219 !-- Calculate vertical depth of the sublayers 220 height_int = ( z_p - z0_av_global ) / REAL( number_of_sublayers ) 221 ! 222 !-- Precalculate LOG(z/z0) 223 height_p = 0.0 224 DO k = 1, number_of_sublayers 225 226 height_p = height_p + height_int 227 log_z_z0(k) = LOG( height_p / z0_av_global ) 228 229 ENDDO 230 231 186 232 ENDIF 187 233 -
palm/trunk/SOURCE/modules.f90
r1310 r1314 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! + log_z_z0, number_of_sublayers, z0_av_global 23 23 ! 24 24 ! Former revisions: … … 1428 1428 maximum_number_of_tailpoints = 100, & 1429 1429 maximum_number_of_tails = 0, & 1430 number_of_sublayers = 20, & 1430 1431 number_of_initial_particles = 0, number_of_particles = 0, & 1431 1432 number_of_particle_groups = 1, number_of_tails = 0, & … … 1465 1466 sgs_wfv_part = 0.3333333, sgs_wfw_part = 0.3333333, & 1466 1467 time_prel = 0.0, time_sort_particles = 0.0, & 1467 time_write_particle_data = 0.0 1468 time_write_particle_data = 0.0, z0_av_global 1468 1469 1469 1470 REAL, DIMENSION(max_number_of_particle_groups) :: & … … 1472 1473 psn = 9999999.9, psr = 9999999.9, pss = 9999999.9, & 1473 1474 pst = 9999999.9, radius = 9999999.9 1475 1476 REAL, DIMENSION(:), ALLOCATABLE :: log_z_z0 1474 1477 1475 1478 REAL, DIMENSION(:,:,:), ALLOCATABLE :: particle_tail_coordinates
Note: See TracChangeset
for help on using the changeset viewer.
