Changeset 1340
- Timestamp:
- Mar 25, 2014 7:45:13 PM (11 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 10 edited
-
diffusion_e.f90 (modified) (11 diffs)
-
diffusion_s.f90 (modified) (14 diffs)
-
diffusion_u.f90 (modified) (17 diffs)
-
diffusion_v.f90 (modified) (18 diffs)
-
diffusion_w.f90 (modified) (8 diffs)
-
diffusivities.f90 (modified) (4 diffs)
-
init_3d_model.f90 (modified) (55 diffs)
-
package_parin.f90 (modified) (2 diffs)
-
plant_canopy_model.f90 (modified) (13 diffs)
-
prandtl_fluxes.f90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/diffusion_e.f90
r1321 r1340 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 148 148 dvar_dz = atmos_ocean_sign * & 149 149 ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 150 IF ( dvar_dz > 0.0 ) THEN151 l_stable = 0.76 * SQRT( e(k,j,i) ) / &152 SQRT( g / var_reference * dvar_dz ) + 1E-5150 IF ( dvar_dz > 0.0_wp ) THEN 151 l_stable = 0.76_wp * SQRT( e(k,j,i) ) / & 152 SQRT( g / var_reference * dvar_dz ) + 1E-5_wp 153 153 ELSE 154 154 l_stable = l_grid(k) … … 176 176 DO k = nzb_s_inner(j,i)+1, nzt 177 177 178 dissipation(k,j) = ( 0.19 + 0.74* l(k,j) / ll(k,j) ) * &178 dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * & 179 179 e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j) 180 180 … … 220 220 dvar_dz = atmos_ocean_sign * & 221 221 ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 222 IF ( dvar_dz > 0.0 ) THEN223 l_stable = 0.76 * SQRT( e(k,j,i) ) / &224 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5222 IF ( dvar_dz > 0.0_wp ) THEN 223 l_stable = 0.76_wp * SQRT( e(k,j,i) ) / & 224 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp 225 225 ELSE 226 226 l_stable = l_grid(k) … … 248 248 DO k = nzb_s_inner(j,i)+1, nzt 249 249 250 dissipation(k,j) = ( 0.19 + 0.74* l(k,j) / ll(k,j) ) * &251 e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)250 dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * & 251 e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j) 252 252 253 253 tend(k,j,i) = tend(k,j,i) & … … 356 356 dvar_dz = atmos_ocean_sign * & 357 357 ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 358 IF ( dvar_dz > 0.0 ) THEN359 l_stable = 0.76 * SQRT( e(k,j,i) ) / &360 SQRT( g / var_reference * dvar_dz ) + 1E-5358 IF ( dvar_dz > 0.0_wp ) THEN 359 l_stable = 0.76_wp * SQRT( e(k,j,i) ) / & 360 SQRT( g / var_reference * dvar_dz ) + 1E-5_wp 361 361 ELSE 362 362 l_stable = l_grid(k) … … 377 377 ! 378 378 !-- Calculate the tendency terms 379 dissipation = ( 0.19 + 0.74* l / ll ) * &380 e(k,j,i) * SQRT( e(k,j,i) ) / l379 dissipation = ( 0.19_wp + 0.74_wp * l / ll ) * & 380 e(k,j,i) * SQRT( e(k,j,i) ) / l 381 381 382 382 tend(k,j,i) = tend(k,j,i) & … … 423 423 dvar_dz = atmos_ocean_sign * & 424 424 ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 425 IF ( dvar_dz > 0.0 ) THEN426 l_stable = 0.76 * SQRT( e(k,j,i) ) / &427 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5425 IF ( dvar_dz > 0.0_wp ) THEN 426 l_stable = 0.76_wp * SQRT( e(k,j,i) ) / & 427 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp 428 428 ELSE 429 429 l_stable = l_grid(k) … … 444 444 ! 445 445 !-- Calculate the tendency terms 446 dissipation = ( 0.19 + 0.74* l / ll ) * &447 e(k,j,i) * SQRT( e(k,j,i) ) / l446 dissipation = ( 0.19_wp + 0.74_wp * l / ll ) * & 447 e(k,j,i) * SQRT( e(k,j,i) ) / l 448 448 449 449 tend(k,j,i) = tend(k,j,i) & … … 541 541 dvar_dz = atmos_ocean_sign * & 542 542 ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 543 IF ( dvar_dz > 0.0 ) THEN543 IF ( dvar_dz > 0.0_wp ) THEN 544 544 IF ( use_single_reference_value ) THEN 545 l_stable = 0.76 * SQRT( e(k,j,i) ) / &546 SQRT( g / var_reference * dvar_dz ) + 1E-5545 l_stable = 0.76_wp * SQRT( e(k,j,i) ) / & 546 SQRT( g / var_reference * dvar_dz ) + 1E-5_wp 547 547 ELSE 548 l_stable = 0.76 * SQRT( e(k,j,i) ) / &549 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5548 l_stable = 0.76_wp * SQRT( e(k,j,i) ) / & 549 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp 550 550 ENDIF 551 551 ELSE … … 566 566 ! 567 567 !-- Calculate the tendency term 568 dissipation(k) = ( 0.19 + 0.74* l(k) / ll(k) ) * e(k,j,i) * &569 SQRT( e(k,j,i) ) / l(k)568 dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) * & 569 SQRT( e(k,j,i) ) / l(k) 570 570 571 571 tend(k,j,i) = tend(k,j,i) & -
palm/trunk/SOURCE/diffusion_s.f90
r1321 r1340 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 120 120 121 121 tend(k,j,i) = tend(k,j,i) & 122 + 0.5 * (&122 + 0.5_wp * ( & 123 123 ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) & 124 124 - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) & 125 ) * ddx2&126 + 0.5 * (&125 ) * ddx2 & 126 + 0.5_wp * ( & 127 127 ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) & 128 128 - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) & 129 ) * ddy2129 ) * ddy2 130 130 ENDDO 131 131 132 132 ! 133 133 !-- Apply prescribed horizontal wall heatflux where necessary 134 IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0) ) &134 IF ( ( wall_w_x(j,i) .NE. 0.0_wp ) .OR. ( wall_w_y(j,i) .NE. 0.0_wp ) ) & 135 135 THEN 136 136 DO k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i) 137 137 138 138 tend(k,j,i) = tend(k,j,i) & 139 + ( fwxp(j,i) * 0.5 *&140 ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) & 141 + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)&142 -fwxm(j,i) * 0.5 *&139 + ( fwxp(j,i) * 0.5_wp * & 140 ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) & 141 + ( 1.0_wp - fwxp(j,i) ) * wall_s_flux(1) & 142 -fwxm(j,i) * 0.5_wp * & 143 143 ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) & 144 + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)&144 + ( 1.0_wp - fwxm(j,i) ) * wall_s_flux(2) & 145 145 ) * ddx2 & 146 + ( fwyp(j,i) * 0.5 *&147 ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) & 148 + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)&149 -fwym(j,i) * 0.5 *&146 + ( fwyp(j,i) * 0.5_wp * & 147 ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) & 148 + ( 1.0_wp - fwyp(j,i) ) * wall_s_flux(3) & 149 -fwym(j,i) * 0.5_wp * & 150 150 ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) & 151 + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)&151 + ( 1.0_wp - fwym(j,i) ) * wall_s_flux(4) & 152 152 ) * ddy2 153 153 ENDDO … … 161 161 162 162 tend(k,j,i) = tend(k,j,i) & 163 + 0.5 * (&163 + 0.5_wp * ( & 164 164 ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) & 165 165 - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) & 166 ) * ddzw(k)166 ) * ddzw(k) 167 167 ENDDO 168 168 … … 175 175 176 176 tend(k,j,i) = tend(k,j,i) & 177 + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )&178 * ( s(k+1,j,i)-s(k,j,i) )&179 * ddzu(k+1)&177 + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) ) & 178 * ( s(k+1,j,i)-s(k,j,i) ) & 179 * ddzu(k+1) & 180 180 + s_flux_b(j,i) & 181 181 ) * ddzw(k) … … 192 192 tend(k,j,i) = tend(k,j,i) & 193 193 + ( - s_flux_t(j,i) & 194 - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )&195 * ( s(k,j,i)-s(k-1,j,i) )&196 * ddzu(k)&194 - 0.5_wp * ( kh(k-1,j,i)+kh(k,j,i) )& 195 * ( s(k,j,i)-s(k-1,j,i) ) & 196 * ddzu(k) & 197 197 ) * ddzw(k) 198 198 … … 251 251 252 252 tend(k,j,i) = tend(k,j,i) & 253 + 0.5 * (&253 + 0.5_wp * ( & 254 254 ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) & 255 255 - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) & 256 ) * ddx2&257 + 0.5 * (&256 ) * ddx2 & 257 + 0.5_wp * ( & 258 258 ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) & 259 259 - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) & 260 ) * ddy2260 ) * ddy2 261 261 ENDIF 262 262 ENDDO … … 266 266 DO k = 1, nzt 267 267 IF ( k > nzb_s_inner(j,i) .AND. k <= nzb_s_outer(j,i) .AND. & 268 ( wall_w_x(j,i) /= 0.0 .OR. wall_w_y(j,i) /= 0.0) ) &268 ( wall_w_x(j,i) /= 0.0_wp .OR. wall_w_y(j,i) /= 0.0_wp ) ) & 269 269 THEN 270 270 tend(k,j,i) = tend(k,j,i) & 271 + ( fwxp(j,i) * 0.5 *&272 ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) & 273 + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)&274 -fwxm(j,i) * 0.5 *&271 + ( fwxp(j,i) * 0.5_wp * & 272 ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) & 273 + ( 1.0_wp - fwxp(j,i) ) * wall_s_flux(1) & 274 -fwxm(j,i) * 0.5_wp * & 275 275 ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) & 276 + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)&276 + ( 1.0_wp - fwxm(j,i) ) * wall_s_flux(2) & 277 277 ) * ddx2 & 278 + ( fwyp(j,i) * 0.5 *&279 ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) & 280 + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)&281 -fwym(j,i) * 0.5 *&278 + ( fwyp(j,i) * 0.5_wp * & 279 ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) & 280 + ( 1.0_wp - fwyp(j,i) ) * wall_s_flux(3) & 281 -fwym(j,i) * 0.5_wp * & 282 282 ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) & 283 + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)&283 + ( 1.0_wp - fwym(j,i) ) * wall_s_flux(4) & 284 284 ) * ddy2 285 285 ENDIF … … 293 293 IF ( k >= nzb_diff_s_inner(j,i) ) THEN 294 294 tend(k,j,i) = tend(k,j,i) & 295 + 0.5 * (&295 + 0.5_wp * ( & 296 296 ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) & 297 297 - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) & 298 ) * ddzw(k)298 ) * ddzw(k) 299 299 ENDIF 300 300 ENDDO … … 306 306 IF ( use_surface_fluxes .AND. k == nzb_s_inner(j,i)+1 ) THEN 307 307 tend(k,j,i) = tend(k,j,i) & 308 + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )&309 * ( s(k+1,j,i)-s(k,j,i) )&310 * ddzu(k+1)&308 + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) )& 309 * ( s(k+1,j,i)-s(k,j,i) ) & 310 * ddzu(k+1) & 311 311 + s_flux_b(j,i) & 312 312 ) * ddzw(k) … … 319 319 tend(k,j,i) = tend(k,j,i) & 320 320 + ( - s_flux_t(j,i) & 321 - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )&322 * ( s(k,j,i)-s(k-1,j,i) ) &323 * ddzu(k) &321 - 0.5_wp * ( kh(k-1,j,i)+kh(k,j,i) )& 322 * ( s(k,j,i)-s(k-1,j,i) ) & 323 * ddzu(k) & 324 324 ) * ddzw(k) 325 325 ENDIF … … 372 372 373 373 tend(k,j,i) = tend(k,j,i) & 374 + 0.5 * (&374 + 0.5_wp * ( & 375 375 ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) & 376 376 - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) & 377 ) * ddx2&378 + 0.5 * (&377 ) * ddx2 & 378 + 0.5_wp * ( & 379 379 ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) & 380 380 - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) & 381 ) * ddy2381 ) * ddy2 382 382 ENDDO 383 383 384 384 ! 385 385 !-- Apply prescribed horizontal wall heatflux where necessary 386 IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0) ) &386 IF ( ( wall_w_x(j,i) .NE. 0.0_wp ) .OR. ( wall_w_y(j,i) .NE. 0.0_wp ) ) & 387 387 THEN 388 388 DO k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i) 389 389 390 390 tend(k,j,i) = tend(k,j,i) & 391 + ( fwxp(j,i) * 0.5 *&392 ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) & 393 + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)&394 -fwxm(j,i) * 0.5 *&391 + ( fwxp(j,i) * 0.5_wp * & 392 ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) & 393 + ( 1.0_wp - fwxp(j,i) ) * wall_s_flux(1) & 394 -fwxm(j,i) * 0.5_wp * & 395 395 ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) & 396 + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)&396 + ( 1.0_wp - fwxm(j,i) ) * wall_s_flux(2) & 397 397 ) * ddx2 & 398 + ( fwyp(j,i) * 0.5 *&399 ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) & 400 + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)&401 -fwym(j,i) * 0.5 *&398 + ( fwyp(j,i) * 0.5_wp * & 399 ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) & 400 + ( 1.0_wp - fwyp(j,i) ) * wall_s_flux(3) & 401 -fwym(j,i) * 0.5_wp * & 402 402 ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) & 403 + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)&403 + ( 1.0_wp - fwym(j,i) ) * wall_s_flux(4) & 404 404 ) * ddy2 405 405 ENDDO … … 413 413 414 414 tend(k,j,i) = tend(k,j,i) & 415 + 0.5 * (&415 + 0.5_wp * ( & 416 416 ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) & 417 417 - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) & 418 ) * ddzw(k)418 ) * ddzw(k) 419 419 ENDDO 420 420 … … 425 425 k = nzb_s_inner(j,i)+1 426 426 427 tend(k,j,i) = tend(k,j,i) + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )&428 * ( s(k+1,j,i)-s(k,j,i) )&429 * ddzu(k+1)&427 tend(k,j,i) = tend(k,j,i) + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) ) & 428 * ( s(k+1,j,i)-s(k,j,i) ) & 429 * ddzu(k+1) & 430 430 + s_flux_b(j,i) & 431 431 ) * ddzw(k) … … 440 440 441 441 tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(j,i) & 442 - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )&443 * ( s(k,j,i)-s(k-1,j,i) )&444 * ddzu(k)&442 - 0.5_wp * ( kh(k-1,j,i)+kh(k,j,i) ) & 443 * ( s(k,j,i)-s(k-1,j,i) ) & 444 * ddzu(k) & 445 445 ) * ddzw(k) 446 446 -
palm/trunk/SOURCE/diffusion_u.f90
r1321 r1340 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 131 131 ! 132 132 !-- Interpolate eddy diffusivities on staggered gridpoints 133 kmyp = 0.25 *&133 kmyp = 0.25_wp * & 134 134 ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 135 kmym = 0.25 *&135 kmym = 0.25_wp * & 136 136 ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 137 137 138 138 tend(k,j,i) = tend(k,j,i) & 139 & + 2.0 * (&139 & + 2.0_wp * ( & 140 140 & km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) & 141 141 & - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) & 142 & ) * ddx2&142 & ) * ddx2 & 143 143 & + ( kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy & 144 144 & + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx & … … 150 150 ! 151 151 !-- Wall functions at the north and south walls, respectively 152 IF ( wall_u(j,i) /= 0.0 ) THEN152 IF ( wall_u(j,i) /= 0.0_wp ) THEN 153 153 154 154 DO k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i) 155 kmyp = 0.25 *&155 kmyp = 0.25_wp * & 156 156 ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 157 kmym = 0.25 *&157 kmym = 0.25_wp * & 158 158 ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 159 159 160 160 tend(k,j,i) = tend(k,j,i) & 161 + 2.0 * (&161 + 2.0_wp * ( & 162 162 km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) & 163 163 - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) & 164 ) * ddx2&164 ) * ddx2 & 165 165 + ( fyp(j,i) * ( & 166 166 kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy & … … 182 182 ! 183 183 !-- Interpolate eddy diffusivities on staggered gridpoints 184 kmzp = 0.25 *&184 kmzp = 0.25_wp * & 185 185 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) ) 186 kmzm = 0.25 *&186 kmzm = 0.25_wp * & 187 187 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) ) 188 188 … … 211 211 ! 212 212 !-- Interpolate eddy diffusivities on staggered gridpoints 213 kmzp = 0.25 *&213 kmzp = 0.25_wp * & 214 214 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) ) 215 kmzm = 0.25 *&215 kmzm = 0.25_wp * & 216 216 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) ) 217 217 … … 231 231 ! 232 232 !-- Interpolate eddy diffusivities on staggered gridpoints 233 kmzp = 0.25 *&233 kmzp = 0.25_wp * & 234 234 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) ) 235 kmzm = 0.25 *&235 kmzm = 0.25_wp * & 236 236 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) ) 237 237 … … 303 303 ! 304 304 !-- Interpolate eddy diffusivities on staggered gridpoints 305 kmyp = 0.25 *&305 kmyp = 0.25_wp * & 306 306 ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 307 kmym = 0.25 *&307 kmym = 0.25_wp * & 308 308 ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 309 309 310 310 tend(k,j,i) = tend(k,j,i) & 311 & + 2.0 * (&311 & + 2.0_wp * ( & 312 312 & km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) & 313 313 & - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) & 314 & ) * ddx2&314 & ) * ddx2 & 315 315 & + ( kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy & 316 316 & + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx & … … 325 325 DO k = 1, nzt 326 326 IF( k > nzb_u_inner(j,i) .AND. k <= nzb_u_outer(j,i) .AND. & 327 wall_u(j,i) /= 0.0 ) THEN328 329 kmyp = 0.25 *&327 wall_u(j,i) /= 0.0_wp ) THEN 328 329 kmyp = 0.25_wp * & 330 330 ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 331 kmym = 0.25 *&331 kmym = 0.25_wp * & 332 332 ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 333 333 334 334 tend(k,j,i) = tend(k,j,i) & 335 + 2.0 * (&335 + 2.0_wp * ( & 336 336 km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) & 337 337 - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) & 338 ) * ddx2&338 ) * ddx2 & 339 339 + ( fyp(j,i) * ( & 340 340 kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy & … … 357 357 ! 358 358 !-- Interpolate eddy diffusivities on staggered gridpoints 359 kmzp = 0.25 *&359 kmzp = 0.25_wp * & 360 360 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) ) 361 kmzm = 0.25 *&361 kmzm = 0.25_wp * & 362 362 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) ) 363 363 … … 394 394 ! 395 395 !-- Interpolate eddy diffusivities on staggered gridpoints 396 kmzp = 0.25 *&396 kmzp = 0.25_wp * & 397 397 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) ) 398 kmzm = 0.25 *&398 kmzm = 0.25_wp * & 399 399 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) ) 400 400 … … 422 422 ! 423 423 !-- Interpolate eddy diffusivities on staggered gridpoints 424 kmzp = 0.25 *&424 kmzp = 0.25_wp * & 425 425 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) ) 426 kmzm = 0.25 *&426 kmzm = 0.25_wp * & 427 427 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) ) 428 428 … … 478 478 ! 479 479 !-- Interpolate eddy diffusivities on staggered gridpoints 480 kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )481 kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )480 kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 481 kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 482 482 483 483 tend(k,j,i) = tend(k,j,i) & 484 & + 2.0 * (&484 & + 2.0_wp * ( & 485 485 & km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) & 486 486 & - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) & 487 & ) * ddx2&487 & ) * ddx2 & 488 488 & + ( kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy & 489 489 & + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx & … … 495 495 ! 496 496 !-- Wall functions at the north and south walls, respectively 497 IF ( wall_u(j,i) .NE. 0.0 ) THEN497 IF ( wall_u(j,i) .NE. 0.0_wp ) THEN 498 498 499 499 ! … … 503 503 504 504 DO k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i) 505 kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )506 kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )505 kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 506 kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 507 507 508 508 tend(k,j,i) = tend(k,j,i) & 509 + 2.0 * (&509 + 2.0_wp * ( & 510 510 km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) & 511 511 - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) & 512 ) * ddx2&512 ) * ddx2 & 513 513 + ( fyp(j,i) * ( & 514 514 kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy & … … 530 530 ! 531 531 !-- Interpolate eddy diffusivities on staggered gridpoints 532 kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )533 kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )532 kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) ) 533 kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) ) 534 534 535 535 tend(k,j,i) = tend(k,j,i) & … … 557 557 ! 558 558 !-- Interpolate eddy diffusivities on staggered gridpoints 559 kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )560 kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )559 kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) ) 560 kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) ) 561 561 562 562 tend(k,j,i) = tend(k,j,i) & … … 575 575 ! 576 576 !-- Interpolate eddy diffusivities on staggered gridpoints 577 kmzp = 0.25 *&577 kmzp = 0.25_wp * & 578 578 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) ) 579 kmzm = 0.25 *&579 kmzm = 0.25_wp * & 580 580 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) ) 581 581 -
palm/trunk/SOURCE/diffusion_v.f90
r1321 r1340 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 129 129 ! 130 130 !-- Interpolate eddy diffusivities on staggered gridpoints 131 kmxp = 0.25 * &131 kmxp = 0.25_wp * & 132 132 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 133 kmxm = 0.25 * &133 kmxm = 0.25_wp * & 134 134 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 135 135 … … 140 140 & - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 141 141 & ) * ddx & 142 & + 2.0 * (&142 & + 2.0_wp * ( & 143 143 & km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) & 144 144 & - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) & 145 & ) * ddy2145 & ) * ddy2 146 146 ENDDO 147 147 148 148 ! 149 149 !-- Wall functions at the left and right walls, respectively 150 IF ( wall_v(j,i) /= 0.0 ) THEN150 IF ( wall_v(j,i) /= 0.0_wp ) THEN 151 151 152 152 DO k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i) 153 kmxp = 0.25 *&153 kmxp = 0.25_wp * & 154 154 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 155 kmxm = 0.25 *&155 kmxm = 0.25_wp * & 156 156 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 157 157 158 158 tend(k,j,i) = tend(k,j,i) & 159 + 2.0 * (&159 + 2.0_wp * ( & 160 160 km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) & 161 161 - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) & 162 ) * ddy2&162 ) * ddy2 & 163 163 + ( fxp(j,i) * ( & 164 164 kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx & … … 180 180 ! 181 181 !-- Interpolate eddy diffusivities on staggered gridpoints 182 kmzp = 0.25 * &182 kmzp = 0.25_wp * & 183 183 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 184 kmzm = 0.25 * &184 kmzm = 0.25_wp * & 185 185 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 186 186 … … 209 209 ! 210 210 !-- Interpolate eddy diffusivities on staggered gridpoints 211 kmzp = 0.25 *&211 kmzp = 0.25_wp * & 212 212 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 213 kmzm = 0.25 *&213 kmzm = 0.25_wp * & 214 214 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 215 215 … … 229 229 ! 230 230 !-- Interpolate eddy diffusivities on staggered gridpoints 231 kmzp = 0.25 *&231 kmzp = 0.25_wp * & 232 232 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 233 kmzm = 0.25 *&233 kmzm = 0.25_wp * & 234 234 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 235 235 … … 301 301 ! 302 302 !-- Interpolate eddy diffusivities on staggered gridpoints 303 kmxp = 0.25 *&303 kmxp = 0.25_wp * & 304 304 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 305 kmxm = 0.25 *&305 kmxm = 0.25_wp * & 306 306 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 307 307 … … 312 312 & - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 313 313 & ) * ddx & 314 & + 2.0 * (&314 & + 2.0_wp * ( & 315 315 & km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) & 316 316 & - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) & 317 & ) * ddy2317 & ) * ddy2 318 318 ENDIF 319 319 ENDDO … … 323 323 DO k = 1, nzt 324 324 IF( k > nzb_v_inner(j,i) .AND. k <= nzb_v_outer(j,i) .AND. & 325 wall_v(j,i) /= 0.0 ) THEN326 327 kmxp = 0.25 *&325 wall_v(j,i) /= 0.0_wp ) THEN 326 327 kmxp = 0.25_wp * & 328 328 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 329 kmxm = 0.25 *&329 kmxm = 0.25_wp * & 330 330 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 331 331 332 332 tend(k,j,i) = tend(k,j,i) & 333 + 2.0 * (&333 + 2.0_wp * ( & 334 334 km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) & 335 335 - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) & 336 ) * ddy2&336 ) * ddy2 & 337 337 + ( fxp(j,i) * ( & 338 338 kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx & … … 355 355 ! 356 356 !-- Interpolate eddy diffusivities on staggered gridpoints 357 kmzp = 0.25 *&357 kmzp = 0.25_wp * & 358 358 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 359 kmzm = 0.25 *&359 kmzm = 0.25_wp * & 360 360 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 361 361 … … 392 392 ! 393 393 !-- Interpolate eddy diffusivities on staggered gridpoints 394 kmzp = 0.25 *&394 kmzp = 0.25_wp * & 395 395 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 396 kmzm = 0.25 *&396 kmzm = 0.25_wp * & 397 397 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 398 398 … … 420 420 ! 421 421 !-- Interpolate eddy diffusivities on staggered gridpoints 422 kmzp = 0.25 *&422 kmzp = 0.25_wp * & 423 423 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 424 kmzm = 0.25 *&424 kmzm = 0.25_wp * & 425 425 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 426 426 … … 476 476 ! 477 477 !-- Interpolate eddy diffusivities on staggered gridpoints 478 kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )479 kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )478 kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 479 kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 480 480 481 481 tend(k,j,i) = tend(k,j,i) & … … 485 485 & - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 486 486 & ) * ddx & 487 & + 2.0 * (&487 & + 2.0_wp * ( & 488 488 & km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) & 489 489 & - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) & 490 & ) * ddy2490 & ) * ddy2 491 491 ENDDO 492 492 493 493 ! 494 494 !-- Wall functions at the left and right walls, respectively 495 IF ( wall_v(j,i) /= 0.0 ) THEN495 IF ( wall_v(j,i) /= 0.0_wp ) THEN 496 496 497 497 ! … … 501 501 502 502 DO k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i) 503 kmxp = 0.25 *&503 kmxp = 0.25_wp * & 504 504 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 505 kmxm = 0.25 *&505 kmxm = 0.25_wp * & 506 506 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 507 507 508 508 tend(k,j,i) = tend(k,j,i) & 509 + 2.0 * (&509 + 2.0_wp * ( & 510 510 km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) & 511 511 - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) & 512 ) * ddy2&512 ) * ddy2 & 513 513 + ( fxp(j,i) * ( & 514 514 kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx & … … 530 530 ! 531 531 !-- Interpolate eddy diffusivities on staggered gridpoints 532 kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )533 kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )532 kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 533 kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 534 534 535 535 tend(k,j,i) = tend(k,j,i) & … … 557 557 ! 558 558 !-- Interpolate eddy diffusivities on staggered gridpoints 559 kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )560 kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )559 kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 560 kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 561 561 562 562 tend(k,j,i) = tend(k,j,i) & … … 575 575 ! 576 576 !-- Interpolate eddy diffusivities on staggered gridpoints 577 kmzp = 0.25 * &577 kmzp = 0.25_wp * & 578 578 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 579 kmzm = 0.25 * &579 kmzm = 0.25_wp * & 580 580 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 581 581 -
palm/trunk/SOURCE/diffusion_w.f90
r1323 r1340 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 158 158 & km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 159 159 & - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 160 & ) * ddzu(k+1)160 & ) * ddzu(k+1) 161 161 ENDDO 162 162 163 163 ! 164 164 !-- Wall functions at all vertical walls, where necessary 165 IF ( wall_w_x(j,i) /= 0.0 .OR. wall_w_y(j,i) /= 0.0) THEN165 IF ( wall_w_x(j,i) /= 0.0_wp .OR. wall_w_y(j,i) /= 0.0_wp ) THEN 166 166 167 167 DO k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i) … … 198 198 + wall_w_y(j,i) * wsvs(k,j,i) & 199 199 ) * ddy & 200 + 2.0 * (&200 + 2.0_wp * ( & 201 201 km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 202 202 - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 203 ) * ddzu(k+1)203 ) * ddzu(k+1) 204 204 ENDDO 205 205 ENDIF … … 288 288 & km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 289 289 & - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 290 & ) * ddzu(k+1)290 & ) * ddzu(k+1) 291 291 ENDIF 292 292 ENDDO … … 297 297 298 298 IF ( k > nzb_w_inner(j,i) .AND. k <= nzb_w_outer(j,i) .AND. & 299 wall_w_x(j,i) /= 0.0 .AND. wall_w_y(j,i) /= 0.0) THEN299 wall_w_x(j,i) /= 0.0_wp .AND. wall_w_y(j,i) /= 0.0_wp ) THEN 300 300 ! 301 301 !-- Interpolate eddy diffusivities on staggered gridpoints … … 330 330 + wall_w_y(j,i) * wsvs(k,j,i) & 331 331 ) * ddy & 332 + 2.0 * (&332 + 2.0_wp * ( & 333 333 km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 334 334 - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 335 ) * ddzu(k+1)335 ) * ddzu(k+1) 336 336 ENDIF 337 337 ENDDO … … 400 400 & km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 401 401 & - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 402 & ) * ddzu(k+1)402 & ) * ddzu(k+1) 403 403 ENDDO 404 404 405 405 ! 406 406 !-- Wall functions at all vertical walls, where necessary 407 IF ( wall_w_x(j,i) /= 0.0 .OR. wall_w_y(j,i) /= 0.0) THEN407 IF ( wall_w_x(j,i) /= 0.0_wp .OR. wall_w_y(j,i) /= 0.0_wp ) THEN 408 408 409 409 ! 410 410 !-- Calculate the horizontal momentum fluxes w'u' and/or w'v' 411 IF ( wall_w_x(j,i) /= 0.0 ) THEN411 IF ( wall_w_x(j,i) /= 0.0_wp ) THEN 412 412 CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i), & 413 413 wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp ) 414 414 ELSE 415 wsus = 0.0 415 wsus = 0.0_wp 416 416 ENDIF 417 417 418 IF ( wall_w_y(j,i) /= 0.0 ) THEN418 IF ( wall_w_y(j,i) /= 0.0_wp ) THEN 419 419 CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i), & 420 420 wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp ) 421 421 ELSE 422 wsvs = 0.0 422 wsvs = 0.0_wp 423 423 ENDIF 424 424 … … 452 452 + wall_w_y(j,i) * wsvs(k) & 453 453 ) * ddy & 454 + 2.0 * (&454 + 2.0_wp * ( & 455 455 km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 456 456 - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 457 ) * ddzu(k+1)457 ) * ddzu(k+1) 458 458 ENDDO 459 459 ENDIF -
palm/trunk/SOURCE/diffusivities.f90
r1323 r1340 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 99 99 ! 100 100 !-- Initialization for calculation of the mixing length profile 101 sums_l_l = 0.0 101 sums_l_l = 0.0_wp 102 102 103 103 ! … … 113 113 ! 114 114 !-- Introduce an optional minimum tke 115 IF ( e_min > 0.0 ) THEN115 IF ( e_min > 0.0_wp ) THEN 116 116 !$OMP DO 117 117 !$acc loop … … 142 142 dvar_dz = atmos_ocean_sign * & ! inverse effect of pt/rho gradient 143 143 ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 144 IF ( dvar_dz > 0.0 ) THEN144 IF ( dvar_dz > 0.0_wp ) THEN 145 145 IF ( use_single_reference_value ) THEN 146 146 l_stable = 0.76_wp * sqrt_e / & 147 SQRT( g / var_reference * dvar_dz ) + 1E-5_wp147 SQRT( g / var_reference * dvar_dz ) + 1E-5_wp 148 148 ELSE 149 149 l_stable = 0.76_wp * sqrt_e / & 150 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp150 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp 151 151 ENDIF 152 152 ELSE -
palm/trunk/SOURCE/init_3d_model.f90
r1323 r1340 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 441 441 IF ( humidity_remote ) THEN 442 442 ALLOCATE( qswst_remote(nysg:nyng,nxlg:nxrg)) 443 qswst_remote = 0.0 443 qswst_remote = 0.0_wp 444 444 ENDIF 445 445 ENDIF … … 453 453 ENDIF 454 454 455 IF ( dt_dosp /= 9999999.9 ) THEN455 IF ( dt_dosp /= 9999999.9_wp ) THEN 456 456 ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ), & 457 457 spectrum_y( 1:ny/2, 1:10, 1:10 ) ) 458 spectrum_x = 0.0 459 spectrum_y = 0.0 458 spectrum_x = 0.0_wp 459 spectrum_y = 0.0_wp 460 460 ENDIF 461 461 … … 463 463 !-- 1D-array for large scale subsidence velocity 464 464 ALLOCATE ( w_subs(nzb:nzt+1) ) 465 w_subs = 0.0 465 w_subs = 0.0_wp 466 466 467 467 … … 480 480 ENDIF 481 481 482 IF ( cthf /= 0.0 ) THEN482 IF ( cthf /= 0.0_wp ) THEN 483 483 ALLOCATE ( lai(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 484 484 canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) … … 491 491 IF ( topography /= 'flat' ) THEN 492 492 ALLOCATE( rif_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg,1:4) ) 493 rif_wall = 0.0 493 rif_wall = 0.0_wp 494 494 ENDIF 495 495 … … 578 578 ALLOCATE( weight_substep(1:intermediate_timestep_count_max), & 579 579 weight_pres(1:intermediate_timestep_count_max) ) 580 weight_substep = 1.0 581 weight_pres = 1.0 580 weight_substep = 1.0_wp 581 weight_pres = 1.0_wp 582 582 intermediate_timestep_count = 1 ! needed when simulated_time = 0.0 583 583 … … 617 617 DO i = nxlg, nxrg 618 618 DO j = nysg, nyng 619 qr(:,j,i) = 0.0 620 nr(:,j,i) = 0.0 619 qr(:,j,i) = 0.0_wp 620 nr(:,j,i) = 0.0_wp 621 621 ENDDO 622 622 ENDDO … … 640 640 IF ( prandtl_layer ) THEN 641 641 rif = rif1d(nzb+1) 642 ts = 0.0 ! could actually be computed more accurately in the643 ! 1D model. Update when opportunity arises.642 ts = 0.0_wp ! could actually be computed more accurately in the 643 ! 1D model. Update when opportunity arises. 644 644 us = us1d 645 645 usws = usws1d 646 646 vsws = vsws1d 647 647 ELSE 648 ts = 0.0 ! must be set, because used in649 rif = 0.0 ! flowste650 us = 0.0 651 usws = 0.0 652 vsws = 0.0 648 ts = 0.0_wp ! must be set, because used in 649 rif = 0.0_wp ! flowste 650 us = 0.0_wp 651 usws = 0.0_wp 652 vsws = 0.0_wp 653 653 ENDIF 654 654 655 655 ELSE 656 e = 0.0 ! must be set, because used in657 rif = 0.0 ! flowste658 ts = 0.0 659 us = 0.0 660 usws = 0.0 661 vsws = 0.0 656 e = 0.0_wp ! must be set, because used in 657 rif = 0.0_wp ! flowste 658 ts = 0.0_wp 659 us = 0.0_wp 660 usws = 0.0_wp 661 vsws = 0.0_wp 662 662 ENDIF 663 663 uswst = top_momentumflux_u … … 669 669 !-- Update when opportunity arises! 670 670 IF ( humidity .OR. passive_scalar ) THEN 671 qs = 0.0 671 qs = 0.0_wp 672 672 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 673 673 precipitation ) THEN 674 qrs = 0.0 675 nrs = 0.0 674 qrs = 0.0_wp 675 nrs = 0.0_wp 676 676 ENDIF 677 677 ENDIF … … 682 682 DO i = nxl-1, nxr+1 683 683 DO j = nys-1, nyn+1 684 u(nzb:nzb_u_inner(j,i),j,i) = 0.0 685 v(nzb:nzb_v_inner(j,i),j,i) = 0.0 684 u(nzb:nzb_u_inner(j,i),j,i) = 0.0_wp 685 v(nzb:nzb_v_inner(j,i),j,i) = 0.0_wp 686 686 ENDDO 687 687 ENDDO … … 745 745 DO i = nxlg, nxrg 746 746 DO j = nysg, nyng 747 u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0 748 v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0 747 u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0_wp 748 v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0_wp 749 749 ENDDO 750 750 ENDDO … … 764 764 DO i = nxlg, nxrg 765 765 DO j = nysg, nyng 766 qr(:,j,i) = 0.0 767 nr(:,j,i) = 0.0 766 qr(:,j,i) = 0.0_wp 767 nr(:,j,i) = 0.0_wp 768 768 ENDDO 769 769 ENDDO … … 784 784 km = km_constant 785 785 kh = km / prandtl_number 786 e = 0.0 787 ELSEIF ( e_init > 0.0 ) THEN786 e = 0.0_wp 787 ELSEIF ( e_init > 0.0_wp ) THEN 788 788 DO k = nzb+1, nzt 789 km(k,:,:) = 0.1 * l_grid(k) * SQRT( e_init )789 km(k,:,:) = 0.1_wp * l_grid(k) * SQRT( e_init ) 790 790 ENDDO 791 791 km(nzb,:,:) = km(nzb+1,:,:) … … 795 795 ELSE 796 796 IF ( .NOT. ocean ) THEN 797 kh = 0.01 ! there must exist an initial diffusion, because798 km = 0.01 ! otherwise no TKE would be produced by the797 kh = 0.01_wp ! there must exist an initial diffusion, because 798 km = 0.01_wp ! otherwise no TKE would be produced by the 799 799 ! production terms, as long as not yet 800 800 ! e = (u*/cm)**2 at k=nzb+1 801 801 ELSE 802 kh = 0.00001 803 km = 0.00001 804 ENDIF 805 e = 0.0 806 ENDIF 807 rif = 0.0 808 ts = 0.0 809 us = 0.0 810 usws = 0.0 802 kh = 0.00001_wp 803 km = 0.00001_wp 804 ENDIF 805 e = 0.0_wp 806 ENDIF 807 rif = 0.0_wp 808 ts = 0.0_wp 809 us = 0.0_wp 810 usws = 0.0_wp 811 811 uswst = top_momentumflux_u 812 vsws = 0.0 812 vsws = 0.0_wp 813 813 vswst = top_momentumflux_v 814 IF ( humidity .OR. passive_scalar ) qs = 0.0 814 IF ( humidity .OR. passive_scalar ) qs = 0.0_wp 815 815 816 816 ! … … 829 829 !-- Bottom boundary 830 830 IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN 831 u(nzb,:,:) = 0.0 832 v(nzb,:,:) = 0.0 831 u(nzb,:,:) = 0.0_wp 832 v(nzb,:,:) = 0.0_wp 833 833 ENDIF 834 834 … … 836 836 !-- Apply channel flow boundary condition 837 837 IF ( TRIM( bc_uv_t ) == 'dirichlet_0' ) THEN 838 u(nzt+1,:,:) = 0.0 839 v(nzt+1,:,:) = 0.0 838 u(nzt+1,:,:) = 0.0_wp 839 v(nzt+1,:,:) = 0.0_wp 840 840 ENDIF 841 841 842 842 ! 843 843 !-- Calculate virtual potential temperature 844 IF ( humidity ) vpt = pt * ( 1.0 + 0.61* q )844 IF ( humidity ) vpt = pt * ( 1.0_wp + 0.61_wp * q ) 845 845 846 846 ! … … 849 849 hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 ) 850 850 IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2) THEN 851 hom(nzb,1,5,:) = 0.0 852 hom(nzb,1,6,:) = 0.0 851 hom(nzb,1,5,:) = 0.0_wp 852 hom(nzb,1,6,:) = 0.0_wp 853 853 ENDIF 854 854 hom(:,1,7,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 ) … … 919 919 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 920 920 precipitation ) THEN 921 qrsws = 0.0 922 nrsws = 0.0 921 qrsws = 0.0_wp 922 nrsws = 0.0_wp 923 923 ENDIF 924 924 IF ( constant_waterflux ) THEN … … 954 954 955 955 IF ( humidity .OR. passive_scalar ) THEN 956 qswst = 0.0 956 qswst = 0.0_wp 957 957 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 958 958 precipitation ) THEN 959 nrswst = 0.0 960 qrswst = 0.0 959 nrswst = 0.0_wp 960 qrswst = 0.0_wp 961 961 ENDIF 962 962 ENDIF … … 971 971 !-- Initialization in case of a coupled model run 972 972 IF ( coupling_mode == 'ocean_to_atmosphere' ) THEN 973 tswst = 0.0 973 tswst = 0.0_wp 974 974 ENDIF 975 975 … … 990 990 !-- to be zero before the first time step. It approaches its correct 991 991 !-- value in the course of the first few time steps. 992 shf = 0.0 992 shf = 0.0_wp 993 993 ENDIF 994 994 995 995 IF ( humidity .OR. passive_scalar ) THEN 996 IF ( .NOT. constant_waterflux ) qsws = 0.0 996 IF ( .NOT. constant_waterflux ) qsws = 0.0_wp 997 997 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 998 998 precipitation ) THEN 999 qrsws = 0.0 1000 nrsws = 0.0 999 qrsws = 0.0_wp 1000 nrsws = 0.0_wp 1001 1001 ENDIF 1002 1002 ENDIF … … 1023 1023 ! 1024 1024 !-- For the moment, vertical velocity is zero 1025 w = 0.0 1025 w = 0.0_wp 1026 1026 1027 1027 ! 1028 1028 !-- Initialize array sums (must be defined in first call of pres) 1029 sums = 0.0 1029 sums = 0.0_wp 1030 1030 1031 1031 ! 1032 1032 !-- In case of iterative solvers, p must get an initial value 1033 IF ( psolver == 'multigrid' .OR. psolver == 'sor' ) p = 0.0 1033 IF ( psolver == 'multigrid' .OR. psolver == 'sor' ) p = 0.0_wp 1034 1034 1035 1035 ! … … 1037 1037 !-- are zero at beginning of the simulation 1038 1038 IF ( cloud_physics ) THEN 1039 ql = 0.0 1040 IF ( precipitation ) precipitation_amount = 0.0 1039 ql = 0.0_wp 1040 IF ( precipitation ) precipitation_amount = 0.0_wp 1041 1041 IF ( icloud_scheme == 0 ) THEN 1042 qc = 0.0 1042 qc = 0.0_wp 1043 1043 nc_1d = nc_const 1044 1044 ENDIF … … 1058 1058 ! 1059 1059 !-- If required, change the surface temperature at the start of the 3D run 1060 IF ( pt_surface_initial_change /= 0.0 ) THEN1060 IF ( pt_surface_initial_change /= 0.0_wp ) THEN 1061 1061 pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change 1062 1062 ENDIF … … 1066 1066 !-- run 1067 1067 IF ( ( humidity .OR. passive_scalar ) .AND. & 1068 q_surface_initial_change /= 0.0 ) THEN1068 q_surface_initial_change /= 0.0_wp ) THEN 1069 1069 q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change 1070 1070 ENDIF … … 1075 1075 ! 1076 1076 !-- Initialize old and new time levels. 1077 te_m = 0.0 ; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.01077 te_m = 0.0_wp; tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp 1078 1078 e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w 1079 1079 1080 1080 IF ( humidity .OR. passive_scalar ) THEN 1081 tq_m = 0.0 1081 tq_m = 0.0_wp 1082 1082 q_p = q 1083 1083 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 1084 1084 precipitation ) THEN 1085 tqr_m = 0.0 1085 tqr_m = 0.0_wp 1086 1086 qr_p = qr 1087 tnr_m = 0.0 1087 tnr_m = 0.0_wp 1088 1088 nr_p = nr 1089 1089 ENDIF … … 1091 1091 1092 1092 IF ( ocean ) THEN 1093 tsa_m = 0.0 1093 tsa_m = 0.0_wp 1094 1094 sa_p = sa 1095 1095 ENDIF … … 1171 1171 u(k,j,nxlg:-1) = mean_inflow_profiles(k,1) 1172 1172 v(k,j,nxlg:-1) = mean_inflow_profiles(k,2) 1173 w(k,j,nxlg:-1) = 0.0 1173 w(k,j,nxlg:-1) = 0.0_wp 1174 1174 pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4) 1175 1175 e(k,j,nxlg:-1) = mean_inflow_profiles(k,5) … … 1183 1183 !-- vertically because otherwise the turbulent inflow layer will grow 1184 1184 !-- in time. 1185 IF ( inflow_damping_height == 9999999.9 ) THEN1185 IF ( inflow_damping_height == 9999999.9_wp ) THEN 1186 1186 ! 1187 1187 !-- Default: use the inversion height calculated by the prerun; if 1188 1188 !-- this is zero, inflow_damping_height must be explicitly 1189 1189 !-- specified. 1190 IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0 ) THEN1190 IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0_wp ) THEN 1191 1191 inflow_damping_height = hom_sum(nzb+6,pr_palm,0) 1192 1192 ELSE … … 1199 1199 ENDIF 1200 1200 1201 IF ( inflow_damping_width == 9999999.9 ) THEN1201 IF ( inflow_damping_width == 9999999.9_wp ) THEN 1202 1202 ! 1203 1203 !-- Default for the transition range: one tenth of the undamped 1204 1204 !-- layer 1205 inflow_damping_width = 0.1 * inflow_damping_height1205 inflow_damping_width = 0.1_wp * inflow_damping_height 1206 1206 1207 1207 ENDIF … … 1212 1212 1213 1213 IF ( zu(k) <= inflow_damping_height ) THEN 1214 inflow_damping_factor(k) = 1.0 1214 inflow_damping_factor(k) = 1.0_wp 1215 1215 ELSEIF ( zu(k) <= ( inflow_damping_height + inflow_damping_width ) ) THEN 1216 inflow_damping_factor(k) = 1.0 -&1216 inflow_damping_factor(k) = 1.0_wp - & 1217 1217 ( zu(k) - inflow_damping_height ) / & 1218 1218 inflow_damping_width 1219 1219 ELSE 1220 inflow_damping_factor(k) = 0.0 1220 inflow_damping_factor(k) = 0.0_wp 1221 1221 ENDIF 1222 1222 … … 1235 1235 DO i = nxlg, nxrg 1236 1236 DO j = nysg, nyng 1237 u (nzb:nzb_u_inner(j,i),j,i) = 0.0 1238 v (nzb:nzb_v_inner(j,i),j,i) = 0.0 1239 w (nzb:nzb_w_inner(j,i),j,i) = 0.0 1240 e (nzb:nzb_w_inner(j,i),j,i) = 0.0 1241 tu_m(nzb:nzb_u_inner(j,i),j,i) = 0.0 1242 tv_m(nzb:nzb_v_inner(j,i),j,i) = 0.0 1243 tw_m(nzb:nzb_w_inner(j,i),j,i) = 0.0 1244 te_m(nzb:nzb_w_inner(j,i),j,i) = 0.0 1245 tpt_m(nzb:nzb_w_inner(j,i),j,i) = 0.0 1237 u (nzb:nzb_u_inner(j,i),j,i) = 0.0_wp 1238 v (nzb:nzb_v_inner(j,i),j,i) = 0.0_wp 1239 w (nzb:nzb_w_inner(j,i),j,i) = 0.0_wp 1240 e (nzb:nzb_w_inner(j,i),j,i) = 0.0_wp 1241 tu_m(nzb:nzb_u_inner(j,i),j,i) = 0.0_wp 1242 tv_m(nzb:nzb_v_inner(j,i),j,i) = 0.0_wp 1243 tw_m(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp 1244 te_m(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp 1245 tpt_m(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp 1246 1246 ENDDO 1247 1247 ENDDO … … 1272 1272 !-- have to be predefined here because they are used (but multiplied with 0) 1273 1273 !-- there before they are set. 1274 te_m = 0.0 ; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.01274 te_m = 0.0_wp; tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp 1275 1275 IF ( humidity .OR. passive_scalar ) THEN 1276 tq_m = 0.0 1276 tq_m = 0.0_wp 1277 1277 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 1278 1278 precipitation ) THEN 1279 tqr_m = 0.0 1280 tnr_m = 0.0 1281 ENDIF 1282 ENDIF 1283 IF ( ocean ) tsa_m = 0.0 1279 tqr_m = 0.0_wp 1280 tnr_m = 0.0_wp 1281 ENDIF 1282 ENDIF 1283 IF ( ocean ) tsa_m = 0.0_wp 1284 1284 1285 1285 ELSE … … 1323 1323 IF ( use_prescribed_profile_data ) THEN 1324 1324 1325 volume_flow_initial_l = 0.0 1326 volume_flow_area_l = 0.0 1325 volume_flow_initial_l = 0.0_wp 1326 volume_flow_area_l = 0.0_wp 1327 1327 1328 1328 IF ( nxr == nx ) THEN … … 1359 1359 ELSEIF ( TRIM( initializing_actions ) == 'cyclic_fill' ) THEN 1360 1360 1361 volume_flow_initial_l = 0.0 1362 volume_flow_area_l = 0.0 1361 volume_flow_initial_l = 0.0_wp 1362 volume_flow_area_l = 0.0_wp 1363 1363 1364 1364 IF ( nxr == nx ) THEN … … 1395 1395 ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 1396 1396 1397 volume_flow_initial_l = 0.0 1398 volume_flow_area_l = 0.0 1397 volume_flow_initial_l = 0.0_wp 1398 volume_flow_area_l = 0.0_wp 1399 1399 1400 1400 IF ( nxr == nx ) THEN … … 1513 1513 DO j = nys, nyn 1514 1514 DO k = nzb, nzt+1 1515 IF ( lad_s(k,j,i) > 0.0 ) THEN1515 IF ( lad_s(k,j,i) > 0.0_wp ) THEN 1516 1516 lad_u(k,j,i) = lad_s(k,j,i) 1517 1517 lad_u(k,j,i+1) = lad_s(k,j,i) … … 1521 1521 ENDDO 1522 1522 DO k = nzb, nzt 1523 lad_w(k,j,i) = 0.5 * ( lad_s(k+1,j,i) + lad_s(k,j,i) )1523 lad_w(k,j,i) = 0.5_wp * ( lad_s(k+1,j,i) + lad_s(k,j,i) ) 1524 1524 ENDDO 1525 1525 ENDDO 1526 1526 ENDDO 1527 1527 1528 lad_w(pch_index,:,:) = 0.0 1528 lad_w(pch_index,:,:) = 0.0_wp 1529 1529 lad_w(nzt+1,:,:) = lad_w(nzt,:,:) 1530 1530 … … 1535 1535 ! 1536 1536 !-- Initialisation of the canopy heat source distribution 1537 IF ( cthf /= 0.0 ) THEN1537 IF ( cthf /= 0.0_wp ) THEN 1538 1538 ! 1539 1539 !-- Piecewise evaluation of the leaf area index by 1540 1540 !-- integration of the leaf area density 1541 lai(:,:,:) = 0.0 1541 lai(:,:,:) = 0.0_wp 1542 1542 DO i = nxlg, nxrg 1543 1543 DO j = nysg, nyng 1544 1544 DO k = pch_index-1, 0, -1 1545 1545 lai(k,j,i) = lai(k+1,j,i) + & 1546 ( 0.5 * ( lad_w(k+1,j,i) +&1546 ( 0.5_wp * ( lad_w(k+1,j,i) + & 1547 1547 lad_s(k+1,j,i) ) * & 1548 1548 ( zw(k+1) - zu(k+1) ) ) + & 1549 ( 0.5 * ( lad_w(k,j,i) +&1549 ( 0.5_wp * ( lad_w(k,j,i) + & 1550 1550 lad_s(k+1,j,i) ) * & 1551 1551 ( zu(k+1) - zw(k) ) ) … … 1561 1561 DO k = 0, pch_index 1562 1562 canopy_heat_flux(k,j,i) = cthf * & 1563 exp( -0.6 * lai(k,j,i) )1563 exp( -0.6_wp * lai(k,j,i) ) 1564 1564 ENDDO 1565 1565 ENDDO … … 1577 1577 ! 1578 1578 !-- If required, initialize dvrp-software 1579 IF ( dt_dvrp /= 9999999.9 ) CALL init_dvrp1579 IF ( dt_dvrp /= 9999999.9_wp ) CALL init_dvrp 1580 1580 1581 1581 IF ( ocean ) THEN … … 1624 1624 ELSE ! for Euler-method 1625 1625 1626 weight_substep(1) = 1.0 1627 weight_pres(1) = 1.0 1626 weight_substep(1) = 1.0_wp 1627 weight_pres(1) = 1.0_wp 1628 1628 1629 1629 ENDIF … … 1631 1631 ! 1632 1632 !-- Initialize Rayleigh damping factors 1633 rdf = 0.0 1634 rdf_sc = 0.0 1635 IF ( rayleigh_damping_factor /= 0.0 ) THEN1633 rdf = 0.0_wp 1634 rdf_sc = 0.0_wp 1635 IF ( rayleigh_damping_factor /= 0.0_wp ) THEN 1636 1636 IF ( .NOT. ocean ) THEN 1637 1637 DO k = nzb+1, nzt 1638 1638 IF ( zu(k) >= rayleigh_damping_height ) THEN 1639 1639 rdf(k) = rayleigh_damping_factor * & 1640 ( SIN( pi * 0.5 * ( zu(k) - rayleigh_damping_height )&1641 / ( zu(nzt) - rayleigh_damping_height ) )&1640 ( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height ) & 1641 / ( zu(nzt) - rayleigh_damping_height ) )& 1642 1642 )**2 1643 1643 ENDIF … … 1647 1647 IF ( zu(k) <= rayleigh_damping_height ) THEN 1648 1648 rdf(k) = rayleigh_damping_factor * & 1649 ( SIN( pi * 0.5 * ( rayleigh_damping_height - zu(k) )&1650 / ( rayleigh_damping_height - zu(nzb+1)))&1649 ( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) ) & 1650 / ( rayleigh_damping_height - zu(nzb+1)))& 1651 1651 )**2 1652 1652 ENDIF … … 1659 1659 !-- Initialize the starting level and the vertical smoothing factor used for 1660 1660 !-- the external pressure gradient 1661 dp_smooth_factor = 1.0 1661 dp_smooth_factor = 1.0_wp 1662 1662 IF ( dp_external ) THEN 1663 1663 ! … … 1670 1670 ENDIF 1671 1671 IF ( dp_smooth ) THEN 1672 dp_smooth_factor(:dp_level_ind_b) = 0.0 1672 dp_smooth_factor(:dp_level_ind_b) = 0.0_wp 1673 1673 DO k = dp_level_ind_b+1, nzt 1674 dp_smooth_factor(k) = 0.5 * ( 1.0 + SIN( pi *&1675 ( REAL( k - dp_level_ind_b, KIND=wp ) /&1676 REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5) ) )1674 dp_smooth_factor(k) = 0.5_wp * ( 1.0_wp + SIN( pi * & 1675 ( REAL( k - dp_level_ind_b, KIND=wp ) / & 1676 REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5_wp ) ) ) 1677 1677 ENDDO 1678 1678 ENDIF … … 1683 1683 !-- non-cyclic lateral boundaries. The damping zone has the maximum value 1684 1684 !-- at the inflow boundary and decreases to zero at pt_damping_width. 1685 ptdf_x = 0.0 1686 ptdf_y = 0.0 1685 ptdf_x = 0.0_wp 1686 ptdf_y = 0.0_wp 1687 1687 IF ( bc_lr_dirrad ) THEN 1688 1688 DO i = nxl, nxr 1689 1689 IF ( ( i * dx ) < pt_damping_width ) THEN 1690 ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5 * &1691 REAL( pt_damping_width - i * dx, KIND=wp ) / (&1692 REAL( pt_damping_width, KIND=wp ) ) ) )**21690 ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp * & 1691 REAL( pt_damping_width - i * dx, KIND=wp ) / ( & 1692 REAL( pt_damping_width, KIND=wp ) ) ) )**2 1693 1693 ENDIF 1694 1694 ENDDO … … 1697 1697 IF ( ( i * dx ) > ( nx * dx - pt_damping_width ) ) THEN 1698 1698 ptdf_x(i) = pt_damping_factor * & 1699 SIN( pi * 0.5 *&1700 ( ( i - nx ) * dx + pt_damping_width ) /&1701 REAL( pt_damping_width, KIND=wp ) )**21699 SIN( pi * 0.5_wp * & 1700 ( ( i - nx ) * dx + pt_damping_width ) / & 1701 REAL( pt_damping_width, KIND=wp ) )**2 1702 1702 ENDIF 1703 1703 ENDDO … … 1706 1706 IF ( ( j * dy ) > ( ny * dy - pt_damping_width ) ) THEN 1707 1707 ptdf_y(j) = pt_damping_factor * & 1708 SIN( pi * 0.5 *&1709 ( ( j - ny ) * dy + pt_damping_width ) /&1710 REAL( pt_damping_width, KIND=wp ) )**21708 SIN( pi * 0.5_wp * & 1709 ( ( j - ny ) * dy + pt_damping_width ) / & 1710 REAL( pt_damping_width, KIND=wp ) )**2 1711 1711 ENDIF 1712 1712 ENDDO … … 1715 1715 IF ( ( j * dy ) < pt_damping_width ) THEN 1716 1716 ptdf_y(j) = pt_damping_factor * & 1717 SIN( pi * 0.5 *&1718 ( pt_damping_width - j * dy ) /&1719 REAL( pt_damping_width, KIND=wp ) )**21717 SIN( pi * 0.5_wp * & 1718 ( pt_damping_width - j * dy ) / & 1719 REAL( pt_damping_width, KIND=wp ) )**2 1720 1720 ENDIF 1721 1721 ENDDO … … 1727 1727 !-- are called from flow_statistics (or - depending on the chosen model run - 1728 1728 !-- are never initialized) 1729 sums_divnew_l = 0.0 1730 sums_divold_l = 0.0 1731 sums_l_l = 0.0 1732 sums_up_fraction_l = 0.0 1733 sums_wsts_bc_l = 0.0 1729 sums_divnew_l = 0.0_wp 1730 sums_divold_l = 0.0_wp 1731 sums_l_l = 0.0_wp 1732 sums_up_fraction_l = 0.0_wp 1733 sums_wsts_bc_l = 0.0_wp 1734 1734 1735 1735 ! … … 1737 1737 !-- Ghost points are excluded because counting values at the ghost boundaries 1738 1738 !-- would bias the statistics 1739 rmask = 1.0 1740 rmask(:,nxlg:nxl-1,:) = 0.0 ; rmask(:,nxr+1:nxrg,:) = 0.01741 rmask(nysg:nys-1,:,:) = 0.0 ; rmask(nyn+1:nyng,:,:) = 0.01739 rmask = 1.0_wp 1740 rmask(:,nxlg:nxl-1,:) = 0.0_wp; rmask(:,nxr+1:nxrg,:) = 0.0_wp 1741 rmask(nysg:nys-1,:,:) = 0.0_wp; rmask(nyn+1:nyng,:,:) = 0.0_wp 1742 1742 1743 1743 ! … … 1769 1769 ngp_2dh_l = 0 1770 1770 ngp_2dh = 0 1771 ngp_3d_inner_l = 0.0 1771 ngp_3d_inner_l = 0.0_wp 1772 1772 ngp_3d_inner = 0 1773 1773 ngp_3d = 0 … … 1777 1777 DO i = nxl, nxr 1778 1778 DO j = nys, nyn 1779 IF ( rmask(j,i,sr) == 1.0 ) THEN1779 IF ( rmask(j,i,sr) == 1.0_wp ) THEN 1780 1780 ! 1781 1781 !-- All xy-grid points -
palm/trunk/SOURCE/package_parin.f90
r1325 r1340 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! REAL constants defined as wp-kinds 22 23 ! 23 24 ! Former revisions: … … 196 197 !-- Default setting of dt_dosp here (instead of check_parameters), because its 197 198 !-- current value is needed in init_pegrid 198 IF ( dt_dosp == 9999999.9 ) dt_dosp = dt_data_output199 IF ( dt_dosp == 9999999.9_wp ) dt_dosp = dt_data_output 199 200 200 201 30 CONTINUE -
palm/trunk/SOURCE/plant_canopy_model.f90
r1321 r1340 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 100 100 v(k,j+1,i) + & 101 101 v(k,j+1,i-1) ) & 102 / 4.0 )**2+ &102 / 4.0_wp )**2 + & 103 103 ( ( w(k-1,j,i-1) + & 104 104 w(k-1,j,i) + & 105 105 w(k,j,i-1) + & 106 106 w(k,j,i) ) & 107 / 4.0 )**2 )* &107 / 4.0_wp )**2 ) * & 108 108 u(k,j,i) 109 109 ENDDO … … 123 123 u(k,j,i) + & 124 124 u(k,j,i+1) ) & 125 / 4.0 )**2+ &125 / 4.0_wp )**2 + & 126 126 v(k,j,i)**2 + & 127 127 ( ( w(k-1,j-1,i) + & … … 129 129 w(k,j-1,i) + & 130 130 w(k,j,i) ) & 131 / 4.0 )**2 )* &131 / 4.0_wp )**2 ) * & 132 132 v(k,j,i) 133 133 ENDDO … … 147 147 u(k+1,j,i) + & 148 148 u(k+1,j,i+1) ) & 149 / 4.0 )**2+ &149 / 4.0_wp )**2 + & 150 150 ( ( v(k,j,i) + & 151 151 v(k,j+1,i) + & 152 152 v(k+1,j,i) + & 153 153 v(k+1,j+1,i) ) & 154 / 4.0 )**2+ &154 / 4.0_wp )**2 + & 155 155 w(k,j,i)**2 ) * & 156 156 w(k,j,i) … … 183 183 SQRT( ( ( u(k,j,i) + & 184 184 u(k,j,i+1) ) & 185 / 2.0 )**2+ &185 / 2.0_wp )**2 + & 186 186 ( ( v(k,j,i) + & 187 187 v(k,j+1,i) ) & 188 / 2.0 )**2+ &188 / 2.0_wp )**2 + & 189 189 ( ( w(k-1,j,i) + & 190 190 w(k,j,i) ) & 191 / 2.0 )**2 )* &191 / 2.0_wp )**2 ) * & 192 192 ( q(k,j,i) - sls(k,j,i) ) 193 193 ENDDO … … 202 202 DO k = nzb_s_inner(j,i)+1, pch_index 203 203 tend(k,j,i) = tend(k,j,i) - & 204 2.0 * cdc(k,j,i) * lad_s(k,j,i) * &204 2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * & 205 205 SQRT( ( ( u(k,j,i) + & 206 206 u(k,j,i+1) ) & 207 / 2.0 )**2+ &207 / 2.0_wp )**2 + & 208 208 ( ( v(k,j,i) + & 209 209 v(k,j+1,i) ) & 210 / 2.0 )**2+ &210 / 2.0_wp )**2 + & 211 211 ( ( w(k,j,i) + & 212 212 w(k+1,j,i) ) & 213 / 2.0 )**2 )* &213 / 2.0_wp )**2 ) * & 214 214 e(k,j,i) 215 215 ENDDO … … 267 267 v(k,j+1,i) + & 268 268 v(k,j+1,i-1) ) & 269 / 4.0 )**2+ &269 / 4.0_wp )**2 + & 270 270 ( ( w(k-1,j,i-1) + & 271 271 w(k-1,j,i) + & 272 272 w(k,j,i-1) + & 273 273 w(k,j,i) ) & 274 / 4.0 )**2 ) *&274 / 4.0_wp )**2 ) * & 275 275 u(k,j,i) 276 276 ENDDO … … 286 286 u(k,j,i) + & 287 287 u(k,j,i+1) ) & 288 / 4.0 )**2+ &288 / 4.0_wp )**2 + & 289 289 v(k,j,i)**2 + & 290 290 ( ( w(k-1,j-1,i) + & … … 292 292 w(k,j-1,i) + & 293 293 w(k,j,i) ) & 294 / 4.0 )**2 ) *&294 / 4.0_wp )**2 ) * & 295 295 v(k,j,i) 296 296 ENDDO … … 306 306 u(k+1,j,i) + & 307 307 u(k+1,j,i+1) ) & 308 / 4.0 )**2+ &308 / 4.0_wp )**2 + & 309 309 ( ( v(k,j,i) + & 310 310 v(k,j+1,i) + & 311 311 v(k+1,j,i) + & 312 312 v(k+1,j+1,i) ) & 313 / 4.0 )**2+ &313 / 4.0_wp )**2 + & 314 314 w(k,j,i)**2 ) * & 315 315 w(k,j,i) … … 336 336 SQRT( ( ( u(k,j,i) + & 337 337 u(k,j,i+1) ) & 338 / 2.0 )**2+ &338 / 2.0_wp )**2 + & 339 339 ( ( v(k,j,i) + & 340 340 v(k,j+1,i) ) & 341 / 2.0 )**2+ &341 / 2.0_wp )**2 + & 342 342 ( ( w(k-1,j,i) + & 343 343 w(k,j,i) ) & 344 / 2.0 )**2 )* &344 / 2.0_wp )**2 ) * & 345 345 ( q(k,j,i) - sls(k,j,i) ) 346 346 ENDDO … … 351 351 DO k = nzb_s_inner(j,i)+1, pch_index 352 352 tend(k,j,i) = tend(k,j,i) - & 353 2.0 * cdc(k,j,i) * lad_s(k,j,i) * &353 2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * & 354 354 SQRT( ( ( u(k,j,i) + & 355 355 u(k,j,i+1) ) & 356 / 2.0 )**2+ &356 / 2.0_wp )**2 + & 357 357 ( ( v(k,j,i) + & 358 358 v(k,j+1,i) ) & 359 / 2.0 )**2+ &359 / 2.0_wp )**2 + & 360 360 ( ( w(k,j,i) + & 361 361 w(k+1,j,i) ) & 362 / 2.0 )**2 )* &362 / 2.0_wp )**2 ) * & 363 363 e(k,j,i) 364 364 ENDDO -
palm/trunk/SOURCE/prandtl_fluxes.f90
r1321 r1340 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 105 105 DO i = nxlg, nxrg 106 106 DO j = nysg, nyng 107 ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30 )107 ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30_wp ) 108 108 ! 109 109 !-- ts must be limited, because otherwise overflow may occur in case of 110 110 !-- us=0 when computing rif further below 111 IF ( ts(j,i) < -1.05E5 ) ts(j,i) = -1.0E5112 IF ( ts(j,i) > 1.0E5 ) ts(j,i) = 1.0E5111 IF ( ts(j,i) < -1.05E5_wp ) ts(j,i) = -1.0E5_wp 112 IF ( ts(j,i) > 1.0E5_wp ) ts(j,i) = 1.0E5_wp 113 113 ENDDO 114 114 ENDDO … … 131 131 z_p = zu(k+1) - zw(k) 132 132 133 IF ( rif(j,i) >= 0.0 ) THEN133 IF ( rif(j,i) >= 0.0_wp ) THEN 134 134 ! 135 135 !-- Stable stratification 136 136 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( & 137 137 LOG( z_p / z0h(j,i) ) + & 138 5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &138 5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p & 139 139 ) 140 140 ELSE 141 141 ! 142 142 !-- Unstable stratification 143 a = SQRT( 1.0 - 16.0* rif(j,i) )144 b = SQRT( 1.0 - 16.0* rif(j,i) * z0h(j,i) / z_p )143 a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) 144 b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p ) 145 145 146 146 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( & 147 147 LOG( z_p / z0h(j,i) ) - & 148 2.0 * LOG( ( 1.0 + a ) / ( 1.0+ b ) ) )148 2.0_wp * LOG( ( 1.0_wp + a ) / ( 1.0_wp + b ) ) ) 149 149 ENDIF 150 150 … … 163 163 z_p = zu(k+1) - zw(k) 164 164 rif(j,i) = z_p * kappa * g * ts(j,i) / & 165 ( pt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )165 ( pt(k+1,j,i) * ( us(j,i)**2 + 1E-30_wp ) ) 166 166 ! 167 167 !-- Limit the value range of the Richardson numbers. … … 182 182 z_p = zu(k+1) - zw(k) 183 183 rif(j,i) = z_p * kappa * g * & 184 ( ts(j,i) + 0.61 * pt(k+1,j,i) * qs(j,i) ) / &185 ( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )184 ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i) ) / & 185 ( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30_wp ) ) 186 186 ! 187 187 !-- Limit the value range of the Richardson numbers. … … 209 209 !-- Compute the absolute value of the horizontal velocity 210 210 !-- (relative to the surface) 211 uv_total = SQRT( ( 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) &212 - u(k,j,i) - u(k,j,i+1) ) )**2 + &213 ( 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) &214 - v(k,j,i) - v(k,j+1,i) ) )**2 )215 216 217 IF ( rif(j,i) >= 0.0 ) THEN211 uv_total = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) & 212 - u(k,j,i) - u(k,j,i+1) ) )**2 + & 213 ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) & 214 - v(k,j,i) - v(k,j+1,i) ) )**2 ) 215 216 217 IF ( rif(j,i) >= 0.0_wp ) THEN 218 218 ! 219 219 !-- Stable stratification 220 220 us(j,i) = kappa * uv_total / ( & 221 221 LOG( z_p / z0(j,i) ) + & 222 5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &222 5.0_wp * rif(j,i) * ( z_p - z0(j,i) ) / z_p & 223 223 ) 224 224 ELSE 225 225 ! 226 226 !-- Unstable stratification 227 a = SQRT( SQRT( 1.0 - 16.0* rif(j,i) ) )228 b = SQRT( SQRT( 1.0 - 16.0* rif(j,i) / z_p * z0(j,i) ) )227 a = SQRT( SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) ) 228 b = SQRT( SQRT( 1.0_wp - 16.0_wp * rif(j,i) / z_p * z0(j,i) ) ) 229 229 230 230 us(j,i) = kappa * uv_total / ( & 231 231 LOG( z_p / z0(j,i) ) - & 232 LOG( ( 1.0 + a )**2 * ( 1.0+ a**2 ) / ( &233 ( 1.0 + b )**2 * ( 1.0+ b**2 ) ) ) + &234 2.0* ( ATAN( a ) - ATAN( b ) ) &232 LOG( ( 1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / ( & 233 ( 1.0_wp + b )**2 * ( 1.0_wp + b**2 ) ) ) + & 234 2.0_wp * ( ATAN( a ) - ATAN( b ) ) & 235 235 ) 236 236 ENDIF … … 258 258 ! 259 259 !-- Compute Richardson-flux number for this point 260 rifm = 0.5 * ( rif(j,i-1) + rif(j,i) )261 IF ( rifm >= 0.0 ) THEN260 rifm = 0.5_wp * ( rif(j,i-1) + rif(j,i) ) 261 IF ( rifm >= 0.0_wp ) THEN 262 262 ! 263 263 !-- Stable stratification 264 264 usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )/ ( & 265 265 LOG( z_p / z0(j,i) ) + & 266 5.0 * rifm * ( z_p - z0(j,i) ) / z_p &267 )266 5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p & 267 ) 268 268 ELSE 269 269 ! 270 270 !-- Unstable stratification 271 a = SQRT( SQRT( 1.0 - 16.0* rifm ) )272 b = SQRT( SQRT( 1.0 - 16.0* rifm / z_p * z0(j,i) ) )271 a = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm ) ) 272 b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) ) 273 273 274 274 usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) ) / ( & 275 275 LOG( z_p / z0(j,i) ) - & 276 LOG( (1.0 + a )**2 * ( 1.0+ a**2 ) / ( &277 (1.0 + b )**2 * ( 1.0+ b**2 ) ) ) + &278 2.0* ( ATAN( a ) - ATAN( b ) ) &276 LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / ( & 277 (1.0_wp + b )**2 * ( 1.0_wp + b**2 ) ) ) + & 278 2.0_wp * ( ATAN( a ) - ATAN( b ) ) & 279 279 ) 280 280 ENDIF 281 usws(j,i) = -usws(j,i) * 0.5 * ( us(j,i-1) + us(j,i) )281 usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) ) 282 282 ENDDO 283 283 ENDDO … … 296 296 ! 297 297 !-- Compute Richardson-flux number for this point 298 rifm = 0.5 * ( rif(j-1,i) + rif(j,i) )299 IF ( rifm >= 0.0 ) THEN298 rifm = 0.5_wp * ( rif(j-1,i) + rif(j,i) ) 299 IF ( rifm >= 0.0_wp ) THEN 300 300 ! 301 301 !-- Stable stratification 302 302 vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / ( & 303 303 LOG( z_p / z0(j,i) ) + & 304 5.0 * rifm * ( z_p - z0(j,i) ) / z_p &305 )304 5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p & 305 ) 306 306 ELSE 307 307 ! 308 308 !-- Unstable stratification 309 a = SQRT( SQRT( 1.0 - 16.0* rifm ) )310 b = SQRT( SQRT( 1.0 - 16.0* rifm / z_p * z0(j,i) ) )309 a = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm ) ) 310 b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) ) 311 311 312 312 vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / ( & 313 313 LOG( z_p / z0(j,i) ) - & 314 LOG( (1.0 + a )**2 * ( 1.0+ a**2 ) / ( &315 (1.0 + b )**2 * ( 1.0+ b**2 ) ) ) + &316 2.0* ( ATAN( a ) - ATAN( b ) ) &314 LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / ( & 315 (1.0_wp + b )**2 * ( 1.0_wp + b**2 ) ) ) + & 316 2.0_wp * ( ATAN( a ) - ATAN( b ) ) & 317 317 ) 318 318 ENDIF 319 vsws(j,i) = -vsws(j,i) * 0.5 * ( us(j-1,i) + us(j,i) )319 vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j-1,i) + us(j,i) ) 320 320 ENDDO 321 321 ENDDO … … 331 331 DO i = nxlg, nxrg 332 332 DO j = nysg, nyng 333 qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30 )333 qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp ) 334 334 ENDDO 335 335 ENDDO … … 355 355 !-- in case of precursor runs) 356 356 IF ( coupled_run ) THEN 357 e_q = 6.1 * &358 EXP( 0.07 * ( MIN(pt(0,j,i),pt(1,j,i)) - 273.15) )359 q(k,j,i) = 0.622 * e_q / ( surface_pressure - e_q )357 e_q = 6.1_wp * & 358 EXP( 0.07_wp * ( MIN(pt(0,j,i),pt(1,j,i)) - 273.15_wp ) ) 359 q(k,j,i) = 0.622_wp * e_q / ( surface_pressure - e_q ) 360 360 ENDIF 361 IF ( rif(j,i) >= 0.0 ) THEN361 IF ( rif(j,i) >= 0.0_wp ) THEN 362 362 ! 363 363 !-- Stable stratification 364 364 qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( & 365 365 LOG( z_p / z0h(j,i) ) + & 366 5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &366 5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p & 367 367 ) 368 368 ELSE 369 369 ! 370 370 !-- Unstable stratification 371 a = SQRT( 1.0 - 16.0* rif(j,i) )372 b = SQRT( 1.0 - 16.0* rif(j,i) * z0h(j,i) / z_p )371 a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) 372 b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p ) 373 373 374 374 qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( & 375 375 LOG( z_p / z0h(j,i) ) - & 376 2.0 * LOG( (1.0 + a ) / ( 1.0+ b ) ) )376 2.0_wp * LOG( (1.0_wp + a ) / ( 1.0_wp + b ) ) ) 377 377 ENDIF 378 378 … … 429 429 !$acc loop independent 430 430 DO j = nysg, nyng 431 e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1 )**2431 e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1_wp )**2 432 432 ! 433 433 !-- As a test: cm = 0.4 434 ! e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4 )**2434 ! e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4_wp )**2 435 435 e(nzb_s_inner(j,i),j,i) = e(nzb_s_inner(j,i)+1,j,i) 436 436 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
