Changeset 1314 for palm/trunk/SOURCE/lpm_advec.f90
 Timestamp:
 Mar 14, 2014 6:25:17 PM (8 years ago)
 File:

 1 edited
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 Prandtllayer 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 ) < 1E8 ) CYCLE 77 78 ! 79 ! Interpolate u velocitycomponent, determine left, front, bottom 80 ! index of uarray 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 xyplane 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 bilinearly 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 ! MoninObukhov relations (if branch). 98 ! First, check if particle is located below first vertical grid level 99 ! (Prandtllayer height) 100 IF ( prandtl_layer .AND. particles(n)%z < z_p ) THEN 101 ! 102 ! Resolvedscale 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,i1) ) 133 134 u_int = usws(j,i) / ( us_int * kappa + 1E10 ) & 135 * log_z_z0_int 136 137 ENDIF 138 ! 139 ! Particle above the first grid level. Bilinear interpolation in the 140 ! horizontal and linear interpolation in the vertical direction. 101 141 ELSE 102 u_int_u = ( ( ggaa ) * u(k+1,j,i) + ( ggbb ) * u(k+1,j,i+1) & 103 + ( ggcc ) * u(k+1,j+1,i) + ( ggdd ) * u(k+1,j+1,i+1) & 142 ! 143 ! Interpolate u velocitycomponent, determine left, front, bottom 144 ! index of uarray. 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 xyplane 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 = ( ( ggaa ) * u(k+1,j,i) + ( ggbb ) * u(k+1,j,i+1) & 168 + ( ggcc ) * u(k+1,j+1,i) + ( ggdd ) * 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 velocitycomponent (adopt 111 ! index k from u velocitycomponent) 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 velocitycomponent. 180 IF ( prandtl_layer .AND. particles(n)%z < z_p ) THEN 181 ! 182 ! Resolvedscale 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(j1,i) ) 200 201 v_int = vsws(j,i) / ( us_int * kappa + 1E10 ) & 202 * log_z_z0_int 203 204 ENDIF 205 ! 206 ! Particle above the first grid level. Bilinear interpolation in the 207 ! horizontal and linear interpolation in the vertical direction. 128 208 ELSE 129 v_int_u = ( ( ggaa ) * v(k+1,j,i) + ( ggbb ) * v(k+1,j,i+1) & 130 + ( ggcc ) * v(k+1,j+1,i) + ( ggdd ) * 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 = ( ( ggaa ) * v(k+1,j,i) + ( ggbb ) * v(k+1,j,i+1) & 226 + ( ggcc ) * v(k+1,j+1,i) + ( ggdd ) * 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 velocitycomponent (adopt 138 ! index i from v velocitycomponent) 235 ! Same procedure for interpolation of the w velocitycomponent 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
Note: See TracChangeset
for help on using the changeset viewer.