Changeset 1353 for palm/trunk/SOURCE/init_1d_model.f90
- Timestamp:
- Apr 8, 2014 3:21:23 PM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/init_1d_model.f90 (modified) (27 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/init_1d_model.f90
r1347 r1353 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants provided with KIND-attribute 23 23 ! 24 24 ! Former revisions: … … 113 113 kh1d = km_constant / prandtl_number 114 114 ELSE 115 e1d = 0.0 ; e1d_p = 0.0116 kh1d = 0.0 ; km1d = 0.0117 rif1d = 0.0 115 e1d = 0.0_wp; e1d_p = 0.0_wp 116 kh1d = 0.0_wp; km1d = 0.0_wp 117 rif1d = 0.0_wp 118 118 ! 119 119 !-- Compute the mixing length 120 l_black(nzb) = 0.0 120 l_black(nzb) = 0.0_wp 121 121 122 122 IF ( TRIM( mixing_length_1d ) == 'blackadar' ) THEN 123 123 ! 124 124 !-- Blackadar mixing length 125 IF ( f /= 0.0 ) THEN126 lambda = 2.7E-4 * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /&127 ABS( f ) + 1E-10 125 IF ( f /= 0.0_wp ) THEN 126 lambda = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) / & 127 ABS( f ) + 1E-10_wp 128 128 ELSE 129 lambda = 30.0 129 lambda = 30.0_wp 130 130 ENDIF 131 131 132 132 DO k = nzb+1, nzt+1 133 l_black(k) = kappa * zu(k) / ( 1.0 + kappa * zu(k) / lambda )133 l_black(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / lambda ) 134 134 ENDDO 135 135 … … 152 152 !-- value in order to avoid too small time steps caused by the diffusion limit 153 153 !-- in the initial phase of a run (at k=1, dz/2 occurs in the limiting formula!) 154 u1d(0:1) = 0.1 155 u1d_p(0:1) = 0.1 156 v1d(0:1) = 0.1 157 v1d_p(0:1) = 0.1 154 u1d(0:1) = 0.1_wp 155 u1d_p(0:1) = 0.1_wp 156 v1d(0:1) = 0.1_wp 157 v1d_p(0:1) = 0.1_wp 158 158 159 159 ! 160 160 !-- For u*, theta* and the momentum fluxes plausible values are set 161 161 IF ( prandtl_layer ) THEN 162 us1d = 0.1 ! without initial friction the flow would not change162 us1d = 0.1_wp ! without initial friction the flow would not change 163 163 ELSE 164 e1d(nzb+1) = 1.0 165 km1d(nzb+1) = 1.0 166 us1d = 0.0 164 e1d(nzb+1) = 1.0_wp 165 km1d(nzb+1) = 1.0_wp 166 us1d = 0.0_wp 167 167 ENDIF 168 ts1d = 0.0 169 usws1d = 0.0 170 vsws1d = 0.0 168 ts1d = 0.0_wp 169 usws1d = 0.0_wp 170 vsws1d = 0.0_wp 171 171 z01d = roughness_length 172 172 z0h1d = z0h_factor * z01d 173 IF ( humidity .OR. passive_scalar ) qs1d = 0.0 173 IF ( humidity .OR. passive_scalar ) qs1d = 0.0_wp 174 174 175 175 ! 176 176 !-- Tendencies must be preset in order to avoid runtime errors within the 177 177 !-- first Runge-Kutta step 178 te_em = 0.0 179 te_um = 0.0 180 te_vm = 0.0 178 te_em = 0.0_wp 179 te_um = 0.0_wp 180 te_vm = 0.0_wp 181 181 182 182 ! … … 270 270 DO k = nzb_diff, nzt 271 271 272 kmzm = 0.5 * ( km1d(k-1) + km1d(k) )273 kmzp = 0.5 * ( km1d(k) + km1d(k+1) )272 kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) ) 273 kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) ) 274 274 ! 275 275 !-- u-component … … 289 289 ! 290 290 !-- TKE 291 kmzm = 0.5 * ( km1d(k-1) + km1d(k) )292 kmzp = 0.5 * ( km1d(k) + km1d(k+1) )291 kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) ) 292 kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) ) 293 293 IF ( .NOT. humidity ) THEN 294 294 pt_0 = pt_init(k) 295 295 flux = ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k) 296 296 ELSE 297 pt_0 = pt_init(k) * ( 1.0 + 0.61* q_init(k) )298 flux = ( ( pt_init(k+1) - pt_init(k-1) ) + &299 0.61 * pt_init(k) * ( q_init(k+1) - q_init(k-1) )&300 ) * dd2zu(k)297 pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) ) 298 flux = ( ( pt_init(k+1) - pt_init(k-1) ) + & 299 0.61_wp * pt_init(k) * & 300 ( q_init(k+1) - q_init(k-1) ) ) * dd2zu(k) 301 301 ENDIF 302 302 … … 304 304 ! 305 305 !-- According to Detering, c_e=0.064 306 dissipation = 0.064 * e1d(k) * SQRT( e1d(k) ) / l1d(k)306 dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d(k) 307 307 ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN 308 dissipation = ( 0.19 + 0.74 * l1d(k) / l_grid(k) )&308 dissipation = ( 0.19_wp + 0.74_wp * l1d(k) / l_grid(k) ) & 309 309 * e1d(k) * SQRT( e1d(k) ) / l1d(k) 310 310 ENDIF … … 329 329 330 330 k = nzb+1 331 kmzm = 0.5 * ( km1d(k-1) + km1d(k) )332 kmzp = 0.5 * ( km1d(k) + km1d(k+1) )331 kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) ) 332 kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) ) 333 333 IF ( .NOT. humidity ) THEN 334 334 pt_0 = pt_init(k) 335 335 flux = ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k) 336 336 ELSE 337 pt_0 = pt_init(k) * ( 1.0 + 0.61* q_init(k) )338 flux = ( ( pt_init(k+1) - pt_init(k-1) ) + &339 0.61 * pt_init(k) * ( q_init(k+1) - q_init(k-1) )&337 pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) ) 338 flux = ( ( pt_init(k+1) - pt_init(k-1) ) + & 339 0.61_wp * pt_init(k) * ( q_init(k+1) - q_init(k-1) ) & 340 340 ) * dd2zu(k) 341 341 ENDIF … … 344 344 ! 345 345 !-- According to Detering, c_e=0.064 346 dissipation = 0.064 * e1d(k) * SQRT( e1d(k) ) / l1d(k)346 dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d(k) 347 347 ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN 348 dissipation = ( 0.19 + 0.74 * l1d(k) / l_grid(k) )&348 dissipation = ( 0.19_wp + 0.74_wp * l1d(k) / l_grid(k) ) & 349 349 * e1d(k) * SQRT( e1d(k) ) / l1d(k) 350 350 ENDIF … … 354 354 te_u(k) = f * ( v1d(k) - vg(k) ) + ( & 355 355 kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) + usws1d & 356 ) * 2.0 * ddzw(k)356 ) * 2.0_wp * ddzw(k) 357 357 ! 358 358 !-- v-component 359 359 te_v(k) = -f * ( u1d(k) - ug(k) ) + ( & 360 360 kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) + vsws1d & 361 ) * 2.0 * ddzw(k)361 ) * 2.0_wp * ddzw(k) 362 362 ! 363 363 !-- TKE … … 394 394 !-- integration due to numerical inaccuracies. In such cases the TKE 395 395 !-- value is reduced to 10 percent of its old value. 396 WHERE ( e1d_p < 0.0 ) e1d_p = 0.1* e1d396 WHERE ( e1d_p < 0.0_wp ) e1d_p = 0.1_wp * e1d 397 397 ENDIF 398 398 … … 417 417 418 418 DO k = nzb+1, nzt 419 te_um(k) = -9.5625 * te_u(k) + 5.3125* te_um(k)420 te_vm(k) = -9.5625 * te_v(k) + 5.3125* te_vm(k)419 te_um(k) = -9.5625_wp * te_u(k) + 5.3125_wp * te_um(k) 420 te_vm(k) = -9.5625_wp * te_v(k) + 5.3125_wp * te_vm(k) 421 421 ENDDO 422 422 423 423 IF ( .NOT. constant_diffusion ) THEN 424 424 DO k = nzb+1, nzt 425 te_em(k) = -9.5625 * te_e(k) + 5.3125* te_em(k)425 te_em(k) = -9.5625_wp * te_e(k) + 5.3125_wp * te_em(k) 426 426 ENDDO 427 427 ENDIF … … 440 440 ! v1d_p(nzb) = -v1d_p(nzb+1) 441 441 442 u1d_p(nzb) = 0.0 443 v1d_p(nzb) = 0.0 442 u1d_p(nzb) = 0.0_wp 443 v1d_p(nzb) = 0.0_wp 444 444 445 445 ! … … 460 460 ! 461 461 !-- Compute theta* using Rif numbers of the previous time step 462 IF ( rif1d(1) >= 0.0 ) THEN462 IF ( rif1d(1) >= 0.0_wp ) THEN 463 463 ! 464 464 !-- Stable stratification 465 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / &466 ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * &467 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) &465 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & 466 ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * rif1d(nzb+1) * & 467 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 468 468 ) 469 469 ELSE 470 470 ! 471 471 !-- Unstable stratification 472 a = SQRT( 1.0 - 16.0 * rif1d(nzb+1) ) 473 b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z0h1d ) 472 a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) ) 473 b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / & 474 zu(nzb+1) * z0h1d ) 474 475 ! 475 476 !-- In the borderline case the formula for stable stratification 476 477 !-- must be applied, because otherwise a zero division would 477 478 !-- occur in the argument of the logarithm. 478 IF ( a == 0.0 .OR. b == 0.0) THEN479 IF ( a == 0.0_wp .OR. b == 0.0_wp ) THEN 479 480 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & 480 ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * & 481 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 481 ( LOG( zu(nzb+1) / z0h1d ) + & 482 5.0_wp * rif1d(nzb+1) * & 483 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 482 484 ) 483 485 ELSE 484 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & 485 LOG( (a-1.0) / (a+1.0) * (b+1.0) / (b-1.0) ) 486 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & 487 LOG( (a-1.0_wp) / (a+1.0_wp) * & 488 (b+1.0_wp) / (b-1.0_wp) ) 486 489 ENDIF 487 490 ENDIF … … 500 503 flux = ts1d 501 504 ELSE 502 pt_0 = pt_init(nzb+1) * ( 1.0 + 0.61* q_init(nzb+1) )503 flux = ts1d + 0.61 * pt_init(k) * qs1d505 pt_0 = pt_init(nzb+1) * ( 1.0_wp + 0.61_wp * q_init(nzb+1) ) 506 flux = ts1d + 0.61_wp * pt_init(k) * qs1d 504 507 ENDIF 505 508 rif1d(nzb+1) = zu(nzb+1) * kappa * g * flux / & 506 ( pt_0 * ( us1d**2 + 1E-30 ) )509 ( pt_0 * ( us1d**2 + 1E-30_wp ) ) 507 510 ENDIF 508 511 … … 512 515 flux = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k) 513 516 ELSE 514 pt_0 = pt_init(k) * ( 1.0 + 0.61* q_init(k) )517 pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) ) 515 518 flux = ( ( pt_init(k+1) - pt_init(k-1) ) & 516 + 0.61 * pt_init(k) * ( q_init(k+1) - q_init(k-1) )& 519 + 0.61_wp * pt_init(k) & 520 * ( q_init(k+1) - q_init(k-1) ) & 517 521 ) * dd2zu(k) 518 522 ENDIF 519 IF ( rif1d(k) >= 0.0 ) THEN520 rif1d(k) = g / pt_0 * flux / &521 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &522 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &523 + 1E-30 523 IF ( rif1d(k) >= 0.0_wp ) THEN 524 rif1d(k) = g / pt_0 * flux / & 525 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & 526 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & 527 + 1E-30_wp & 524 528 ) 525 529 ELSE 526 rif1d(k) = g / pt_0 * flux / &527 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &528 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &529 + 1E-30 530 ) * ( 1.0 - 16.0 * rif1d(k) )**0.25530 rif1d(k) = g / pt_0 * flux / & 531 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & 532 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & 533 + 1E-30_wp & 534 ) * ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp 531 535 ENDIF 532 536 ENDDO … … 543 547 uv_total = SQRT( u1d(nzb+1)**2 + v1d(nzb+1)**2 ) 544 548 545 IF ( rif1d(nzb+1) >= 0.0 ) THEN549 IF ( rif1d(nzb+1) >= 0.0_wp ) THEN 546 550 ! 547 551 !-- Stable stratification 548 552 us1d = kappa * uv_total / ( & 549 LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) *&553 LOG( zu(nzb+1) / z01d ) + 5.0_wp * rif1d(nzb+1) * & 550 554 ( zu(nzb+1) - z01d ) / zu(nzb+1) & 551 555 ) … … 553 557 ! 554 558 !-- Unstable stratification 555 a = 1.0 / SQRT( SQRT( 1.0 - 16.0* rif1d(nzb+1) ) )556 b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1)&557 * z01d ) )559 a = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) ) ) 560 b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / & 561 zu(nzb+1) * z01d ) ) 558 562 ! 559 563 !-- In the borderline case the formula for stable stratification 560 564 !-- must be applied, because otherwise a zero division would 561 565 !-- occur in the argument of the logarithm. 562 IF ( a == 1.0 .OR. b == 1.0) THEN563 us1d = kappa * uv_total / ( &564 LOG( zu(nzb+1) / z01d ) + &565 5.0 * rif1d(nzb+1) * ( zu(nzb+1) - z01d ) /&566 IF ( a == 1.0_wp .OR. b == 1.0_wp ) THEN 567 us1d = kappa * uv_total / ( & 568 LOG( zu(nzb+1) / z01d ) + & 569 5.0_wp * rif1d(nzb+1) * ( zu(nzb+1) - z01d ) / & 566 570 zu(nzb+1) ) 567 571 ELSE 568 572 us1d = kappa * uv_total / ( & 569 LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) +& 570 2.0 * ( ATAN( b ) - ATAN( a ) ) & 573 LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) / & 574 (1.0_wp+a) ) + & 575 2.0_wp * ( ATAN( b ) - ATAN( a ) ) & 571 576 ) 572 577 ENDIF … … 584 589 !-- compatibility with the 3D model. 585 590 IF ( ibc_e_b == 2 ) THEN 586 e1d(nzb+1) = ( us1d / 0.1 )**2587 ! e1d(nzb+1) = ( us1d / 0.4 )**2 !not used so far, see also588 !prandtl_fluxes591 e1d(nzb+1) = ( us1d / 0.1_wp )**2 592 ! e1d(nzb+1) = ( us1d / 0.4_wp )**2 !not used so far, see also 593 !prandtl_fluxes 589 594 ENDIF 590 595 e1d(nzb) = e1d(nzb+1) … … 593 598 ! 594 599 !-- Compute q* 595 IF ( rif1d(1) >= 0.0 ) THEN600 IF ( rif1d(1) >= 0.0_wp ) THEN 596 601 ! 597 602 !-- Stable stratification 598 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / &599 ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * &600 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) &603 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & 604 ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * rif1d(nzb+1) * & 605 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 601 606 ) 602 607 ELSE 603 608 ! 604 609 !-- Unstable stratification 605 a = SQRT( 1.0 - 16.0 * rif1d(nzb+1) ) 606 b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z0h1d ) 610 a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) ) 611 b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / & 612 zu(nzb+1) * z0h1d ) 607 613 ! 608 614 !-- In the borderline case the formula for stable stratification 609 615 !-- must be applied, because otherwise a zero division would 610 616 !-- occur in the argument of the logarithm. 611 IF ( a == 1.0 .OR. b == 1.0) THEN617 IF ( a == 1.0_wp .OR. b == 1.0_wp ) THEN 612 618 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & 613 ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * & 614 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 619 ( LOG( zu(nzb+1) / z0h1d ) + & 620 5.0_wp * rif1d(nzb+1) * & 621 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 615 622 ) 616 623 ELSE 617 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & 618 LOG( (a-1.0) / (a+1.0) * (b+1.0) / (b-1.0) ) 624 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & 625 LOG( (a-1.0_wp) / (a+1.0_wp) * & 626 (b+1.0_wp) / (b-1.0_wp) ) 619 627 ENDIF 620 628 ENDIF 621 629 ELSE 622 qs1d = 0.0 630 qs1d = 0.0_wp 623 631 ENDIF 624 632 … … 629 637 IF ( mixing_length_1d == 'blackadar' ) THEN 630 638 DO k = nzb+1, nzt 631 IF ( rif1d(k) >= 0.0 ) THEN632 l1d(k) = l_black(k) / ( 1.0 + 5.0* rif1d(k) )639 IF ( rif1d(k) >= 0.0_wp ) THEN 640 l1d(k) = l_black(k) / ( 1.0_wp + 5.0_wp * rif1d(k) ) 633 641 ELSE 634 l1d(k) = l_black(k) * ( 1.0 - 16.0 * rif1d(k) )**0.25 642 l1d(k) = l_black(k) * & 643 ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp 635 644 ENDIF 636 645 l1d(k) = l_black(k) … … 640 649 DO k = nzb+1, nzt 641 650 dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k) 642 IF ( dpt_dz > 0.0 ) THEN643 l_stable = 0.76 * SQRT( e1d(k) ) /&644 SQRT( g / pt_init(k) * dpt_dz ) + 1E-5 651 IF ( dpt_dz > 0.0_wp ) THEN 652 l_stable = 0.76_wp * SQRT( e1d(k) ) / & 653 SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp 645 654 ELSE 646 655 l_stable = l_grid(k) … … 657 666 !-- already been taken account of via the TKE (cf. also Diss.). 658 667 IF ( prandtl_layer ) THEN 659 IF ( rif1d(nzb+1) >= 0.0 ) THEN660 km1d(nzb+1) = us1d * kappa * zu(nzb+1) / &661 ( 1.0 + 5.0* rif1d(nzb+1) )662 ELSE 663 km1d(nzb+1) = us1d * kappa * zu(nzb+1) * &664 ( 1.0 - 16.0 * rif1d(nzb+1) )**0.25668 IF ( rif1d(nzb+1) >= 0.0_wp ) THEN 669 km1d(nzb+1) = us1d * kappa * zu(nzb+1) / & 670 ( 1.0_wp + 5.0_wp * rif1d(nzb+1) ) 671 ELSE 672 km1d(nzb+1) = us1d * kappa * zu(nzb+1) * & 673 ( 1.0_wp - 16.0_wp * rif1d(nzb+1) )**0.25_wp 665 674 ENDIF 666 675 ENDIF 667 676 DO k = nzb_diff, nzt 668 677 ! km1d(k) = 0.4 * SQRT( e1d(k) ) !changed: adjustment to 3D-model 669 km1d(k) = 0.1 * SQRT( e1d(k) )670 IF ( rif1d(k) >= 0.0 ) THEN678 km1d(k) = 0.1_wp * SQRT( e1d(k) ) 679 IF ( rif1d(k) >= 0.0_wp ) THEN 671 680 km1d(k) = km1d(k) * l1d(k) 672 681 ELSE … … 678 687 !-- Add damping layer 679 688 DO k = damp_level_ind_1d+1, nzt+1 680 km1d(k) = 1.1 * km1d(k-1)689 km1d(k) = 1.1_wp * km1d(k-1) 681 690 km1d(k) = MIN( km1d(k), 10.0_wp ) 682 691 ENDDO … … 686 695 !-- kh = phim / phih * km 687 696 DO k = nzb+1, nzt 688 IF ( rif1d(k) >= 0.0 ) THEN697 IF ( rif1d(k) >= 0.0_wp ) THEN 689 698 kh1d(k) = km1d(k) 690 699 ELSE 691 kh1d(k) = km1d(k) * ( 1.0 - 16.0 * rif1d(k) )**0.25700 kh1d(k) = km1d(k) * ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp 692 701 ENDIF 693 702 ENDDO … … 778 787 !-- Compute control quantities 779 788 !-- grid level nzp is excluded due to mirror boundary condition 780 umax = 0.0 ; vmax = 0.0; energy = 0.0789 umax = 0.0_wp; vmax = 0.0_wp; energy = 0.0_wp 781 790 DO k = nzb+1, nzt+1 782 791 umax = MAX( ABS( umax ), ABS( u1d(k) ) ) 783 792 vmax = MAX( ABS( vmax ), ABS( v1d(k) ) ) 784 energy = energy + 0.5 * ( u1d(k)**2 + v1d(k)**2 )793 energy = energy + 0.5_wp * ( u1d(k)**2 + v1d(k)**2 ) 785 794 ENDDO 786 795 energy = energy / REAL( nzt - nzb + 1, KIND=wp ) 787 796 788 797 uv_total = SQRT( u1d(nzb+1)**2 + v1d(nzb+1)**2 ) 789 IF ( ABS( v1d(nzb+1) ) .LT. 1.0E-5 ) THEN798 IF ( ABS( v1d(nzb+1) ) .LT. 1.0E-5_wp ) THEN 790 799 alpha = ACOS( SIGN( 1.0_wp , u1d(nzb+1) ) ) 791 800 ELSE 792 801 alpha = ACOS( u1d(nzb+1) / uv_total ) 793 IF ( v1d(nzb+1) <= 0.0 ) alpha = 2.0* pi - alpha802 IF ( v1d(nzb+1) <= 0.0_wp ) alpha = 2.0_wp * pi - alpha 794 803 ENDIF 795 alpha = alpha / ( 2.0 * pi ) * 360.0804 alpha = alpha / ( 2.0_wp * pi ) * 360.0_wp 796 805 797 806 WRITE ( 15, 101 ) current_timestep_number_1d, simulated_time_chr, & … … 853 862 !-- Compute the currently feasible time step according to the diffusion 854 863 !-- criterion. At nzb+1 the half grid length is used. 855 fac = 0.35 864 fac = 0.35_wp 856 865 dt_diff = dt_max_1d 857 866 DO k = nzb+2, nzt 858 value = fac * dzu(k) * dzu(k) / ( km1d(k) + 1E-20 )867 value = fac * dzu(k) * dzu(k) / ( km1d(k) + 1E-20_wp ) 859 868 dt_diff = MIN( value, dt_diff ) 860 869 ENDDO 861 value = fac * zu(nzb+1) * zu(nzb+1) / ( km1d(nzb+1) + 1E-20 )870 value = fac * zu(nzb+1) * zu(nzb+1) / ( km1d(nzb+1) + 1E-20_wp ) 862 871 dt_1d = MIN( value, dt_diff ) 863 872 864 873 ! 865 874 !-- Set flag when the time step becomes too small 866 IF ( dt_1d < ( 0.00001 * dt_max_1d ) ) THEN875 IF ( dt_1d < ( 0.00001_wp * dt_max_1d ) ) THEN 867 876 stop_dt_1d = .TRUE. 868 877 … … 876 885 !-- A more or less simple new time step value is obtained taking only the 877 886 !-- first two significant digits 878 div = 1000.0 887 div = 1000.0_wp 879 888 DO WHILE ( dt_1d < div ) 880 div = div / 10.0 889 div = div / 10.0_wp 881 890 ENDDO 882 dt_1d = NINT( dt_1d * 100.0 / div ) * div / 100.0891 dt_1d = NINT( dt_1d * 100.0_wp / div ) * div / 100.0_wp 883 892 884 893 old_dt_1d = dt_1d
Note: See TracChangeset
for help on using the changeset viewer.
