Changeset 978 for palm/trunk/SOURCE/diffusion_w.f90
- Timestamp:
- Aug 9, 2012 8:28:32 AM (12 years ago)
- Location:
- palm/trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk
- Property svn:mergeinfo changed
/palm/branches/fricke (added) merged: 942-944,967-968,971-972,977
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE
- Property svn:mergeinfo changed
/palm/branches/fricke/SOURCE (added) merged: 967-968,971-972,977
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE/diffusion_w.f90
r668 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! outflow damping layer removed 7 ! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp 8 ! kmym_y/_z and kmyp_y/_z change to kmym and kmyp 6 9 ! 7 10 ! Former revisions: … … 54 57 ! Call for all grid points 55 58 !------------------------------------------------------------------------------! 56 SUBROUTINE diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, & 57 w ) 59 SUBROUTINE diffusion_w( ddzu, ddzw, km, tend, u, v, w ) 58 60 59 61 USE control_parameters … … 64 66 65 67 INTEGER :: i, j, k 66 REAL :: kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, & 67 kmyp_z 68 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg), & 69 km_damp_y(nysg:nyng) 68 REAL :: kmxm, kmxp, kmym, kmyp 69 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 70 70 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 71 71 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w … … 88 88 ! 89 89 !-- Interpolate eddy diffusivities on staggered gridpoints 90 kmxp_x = 0.25 * & 91 ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 92 kmxm_x = 0.25 * & 93 ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 94 kmxp_z = kmxp_x 95 kmxm_z = kmxm_x 96 kmyp_y = 0.25 * & 97 ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 98 kmym_y = 0.25 * & 99 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 100 kmyp_z = kmyp_y 101 kmym_z = kmym_y 102 ! 103 !-- Increase diffusion at the outflow boundary in case of 104 !-- non-cyclic lateral boundaries. Damping is only needed for 105 !-- velocity components parallel to the outflow boundary in 106 !-- the direction normal to the outflow boundary. 107 IF ( .NOT. bc_lr_cyc ) THEN 108 kmxp_x = MAX( kmxp_x, km_damp_x(i) ) 109 kmxm_x = MAX( kmxm_x, km_damp_x(i) ) 110 ENDIF 111 IF ( .NOT. bc_ns_cyc ) THEN 112 kmyp_y = MAX( kmyp_y, km_damp_y(j) ) 113 kmym_y = MAX( kmym_y, km_damp_y(j) ) 114 ENDIF 90 kmxp = 0.25 * & 91 ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 92 kmxm = 0.25 * & 93 ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 94 kmyp = 0.25 * & 95 ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 96 kmym = 0.25 * & 97 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 115 98 116 99 tend(k,j,i) = tend(k,j,i) & 117 & + ( kmxp _x * ( w(k,j,i+1) - w(k,j,i) ) * ddx&118 & + kmxp _z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)&119 & - kmxm _x * ( w(k,j,i) - w(k,j,i-1) ) * ddx&120 & - kmxm _z * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)&100 & + ( kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx & 101 & + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 102 & - kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx & 103 & - kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 121 104 & ) * ddx & 122 & + ( kmyp _y * ( w(k,j+1,i) - w(k,j,i) ) * ddy&123 & + kmyp _z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)&124 & - kmym _y * ( w(k,j,i) - w(k,j-1,i) ) * ddy&125 & - kmym _z * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)&105 & + ( kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy & 106 & + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 107 & - kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 108 & - kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 126 109 & ) * ddy & 127 110 & + 2.0 * ( & … … 138 121 ! 139 122 !-- Interpolate eddy diffusivities on staggered gridpoints 140 kmxp_x = 0.25 * & 141 ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 142 kmxm_x = 0.25 * & 143 ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 144 kmxp_z = kmxp_x 145 kmxm_z = kmxm_x 146 kmyp_y = 0.25 * & 147 ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 148 kmym_y = 0.25 * & 149 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 150 kmyp_z = kmyp_y 151 kmym_z = kmym_y 152 ! 153 !-- Increase diffusion at the outflow boundary in case of 154 !-- non-cyclic lateral boundaries. Damping is only needed for 155 !-- velocity components parallel to the outflow boundary in 156 !-- the direction normal to the outflow boundary. 157 IF ( .NOT. bc_lr_cyc ) THEN 158 kmxp_x = MAX( kmxp_x, km_damp_x(i) ) 159 kmxm_x = MAX( kmxm_x, km_damp_x(i) ) 160 ENDIF 161 IF ( .NOT. bc_ns_cyc ) THEN 162 kmyp_y = MAX( kmyp_y, km_damp_y(j) ) 163 kmym_y = MAX( kmym_y, km_damp_y(j) ) 164 ENDIF 123 kmxp = 0.25 * & 124 ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 125 kmxm = 0.25 * & 126 ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 127 kmyp = 0.25 * & 128 ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 129 kmym = 0.25 * & 130 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 165 131 166 132 tend(k,j,i) = tend(k,j,i) & 167 133 + ( fwxp(j,i) * ( & 168 kmxp _x * ( w(k,j,i+1) - w(k,j,i) ) * ddx&169 + kmxp _z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)&134 kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx & 135 + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 170 136 ) & 171 137 - fwxm(j,i) * ( & 172 kmxm _x * ( w(k,j,i) - w(k,j,i-1) ) * ddx&173 + kmxm _z * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)&138 kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx & 139 + kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 174 140 ) & 175 141 + wall_w_x(j,i) * wsus(k,j,i) & 176 142 ) * ddx & 177 143 + ( fwyp(j,i) * ( & 178 kmyp _y * ( w(k,j+1,i) - w(k,j,i) ) * ddy&179 + kmyp _z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)&144 kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy & 145 + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 180 146 ) & 181 147 - fwym(j,i) * ( & 182 kmym _y * ( w(k,j,i) - w(k,j-1,i) ) * ddy&183 + kmym _z * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)&148 kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 149 + kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 184 150 ) & 185 151 + wall_w_y(j,i) * wsvs(k,j,i) & … … 201 167 ! Call for grid point i,j 202 168 !------------------------------------------------------------------------------! 203 SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, & 204 tend, u, v, w ) 169 SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, tend, u, v, w ) 205 170 206 171 USE control_parameters … … 211 176 212 177 INTEGER :: i, j, k 213 REAL :: kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, & 214 kmyp_z 215 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg), & 216 km_damp_y(nysg:nyng) 178 REAL :: kmxm, kmxp, kmym, kmyp 179 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 217 180 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 218 181 REAL, DIMENSION(nzb:nzt+1) :: wsus, wsvs … … 223 186 ! 224 187 !-- Interpolate eddy diffusivities on staggered gridpoints 225 kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 226 kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 227 kmxp_z = kmxp_x 228 kmxm_z = kmxm_x 229 kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 230 kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 231 kmyp_z = kmyp_y 232 kmym_z = kmym_y 233 ! 234 !-- Increase diffusion at the outflow boundary in case of non-cyclic 235 !-- lateral boundaries. Damping is only needed for velocity components 236 !-- parallel to the outflow boundary in the direction normal to the 237 !-- outflow boundary. 238 IF ( .NOT. bc_lr_cyc ) THEN 239 kmxp_x = MAX( kmxp_x, km_damp_x(i) ) 240 kmxm_x = MAX( kmxm_x, km_damp_x(i) ) 241 ENDIF 242 IF ( .NOT. bc_ns_cyc ) THEN 243 kmyp_y = MAX( kmyp_y, km_damp_y(j) ) 244 kmym_y = MAX( kmym_y, km_damp_y(j) ) 245 ENDIF 188 kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 189 kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 190 kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 191 kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 246 192 247 193 tend(k,j,i) = tend(k,j,i) & 248 & + ( kmxp _x * ( w(k,j,i+1) - w(k,j,i) ) * ddx&249 & + kmxp _z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)&250 & - kmxm _x * ( w(k,j,i) - w(k,j,i-1) ) * ddx&251 & - kmxm _z * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)&194 & + ( kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx & 195 & + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 196 & - kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx & 197 & - kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 252 198 & ) * ddx & 253 & + ( kmyp _y * ( w(k,j+1,i) - w(k,j,i) ) * ddy&254 & + kmyp _z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)&255 & - kmym _y * ( w(k,j,i) - w(k,j-1,i) ) * ddy&256 & - kmym _z * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)&199 & + ( kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy & 200 & + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 201 & - kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 202 & - kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 257 203 & ) * ddy & 258 204 & + 2.0 * ( & … … 285 231 ! 286 232 !-- Interpolate eddy diffusivities on staggered gridpoints 287 kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 288 kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 289 kmxp_z = kmxp_x 290 kmxm_z = kmxm_x 291 kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 292 kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 293 kmyp_z = kmyp_y 294 kmym_z = kmym_y 295 ! 296 !-- Increase diffusion at the outflow boundary in case of 297 !-- non-cyclic lateral boundaries. Damping is only needed for 298 !-- velocity components parallel to the outflow boundary in 299 !-- the direction normal to the outflow boundary. 300 IF ( .NOT. bc_lr_cyc ) THEN 301 kmxp_x = MAX( kmxp_x, km_damp_x(i) ) 302 kmxm_x = MAX( kmxm_x, km_damp_x(i) ) 303 ENDIF 304 IF ( .NOT. bc_ns_cyc ) THEN 305 kmyp_y = MAX( kmyp_y, km_damp_y(j) ) 306 kmym_y = MAX( kmym_y, km_damp_y(j) ) 307 ENDIF 233 kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 234 kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 235 kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 236 kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 308 237 309 238 tend(k,j,i) = tend(k,j,i) & 310 239 + ( fwxp(j,i) * ( & 311 kmxp _x * ( w(k,j,i+1) - w(k,j,i) ) * ddx&312 + kmxp _z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)&240 kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx & 241 + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 313 242 ) & 314 243 - fwxm(j,i) * ( & 315 kmxm _x * ( w(k,j,i) - w(k,j,i-1) ) * ddx&316 + kmxm _z * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)&244 kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx & 245 + kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 317 246 ) & 318 247 + wall_w_x(j,i) * wsus(k) & 319 248 ) * ddx & 320 249 + ( fwyp(j,i) * ( & 321 kmyp _y * ( w(k,j+1,i) - w(k,j,i) ) * ddy&322 + kmyp _z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)&250 kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy & 251 + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 323 252 ) & 324 253 - fwym(j,i) * ( & 325 kmym _y * ( w(k,j,i) - w(k,j-1,i) ) * ddy&326 + kmym _z * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)&254 kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 255 + kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 327 256 ) & 328 257 + wall_w_y(j,i) * wsvs(k) &
Note: See TracChangeset
for help on using the changeset viewer.
