Changeset 53 for palm/trunk/SOURCE/production_e.f90
- Timestamp:
- Mar 7, 2007 12:33:47 PM (17 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/production_e.f90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/production_e.f90
r39 r53 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes 7 7 ! 8 8 ! Former revisions: … … 67 67 68 68 REAL :: def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, & 69 k1, k2, theta, temp, usvs, vsus, wsus, wsvs 69 k1, k2, theta, temp 70 71 REAL, DIMENSION(nzb:nzt+1) :: usvs, vsus, wsus, wsvs 70 72 71 73 … … 119 121 k = nzb_diff_s_inner(j,i) - 1 120 122 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 123 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & 124 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 125 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 126 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & 127 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 128 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 129 121 130 IF ( wall_e_y(j,i) /= 0.0 ) THEN 122 usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) & 123 / LOG( 0.5 * dy / z0(j,i) ) 124 usvs = usvs * ABS( usvs ) 125 dudy = wall_e_y(j,i) * usvs / km(k,j,i) 131 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 132 usvs, 1.0, 0.0, 0.0, 0.0 ) 133 134 dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i) 135 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 136 wsvs, 0.0, 0.0, 1.0, 0.0 ) 137 dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i) 126 138 ELSE 127 139 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & 128 140 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 141 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & 142 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 129 143 ENDIF 130 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &131 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)132 144 133 145 IF ( wall_e_x(j,i) /= 0.0 ) THEN 134 vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) & 135 / LOG( 0.5 * dx / z0(j,i)) 136 vsus = vsus * ABS( vsus ) 137 dvdx = wall_e_x(j,i) * vsus / km(k,j,i) 146 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 147 vsus, 0.0, 1.0, 0.0, 0.0 ) 148 dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i) 149 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 150 wsus, 0.0, 0.0, 0.0, 1.0 ) 151 dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i) 138 152 ELSE 139 153 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & 140 154 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 141 ENDIF142 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy143 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &144 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)145 146 IF ( wall_e_x(j,i) /= 0.0 ) THEN147 wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &148 / LOG( 0.5 * dx / z0(j,i))149 wsus = wsus * ABS( wsus )150 dwdx = wall_e_x(j,i) * wsus / km(k,j,i)151 ELSE152 155 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & 153 156 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 154 157 ENDIF 155 IF ( wall_e_y(j,i) /= 0.0 ) THEN156 wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &157 / LOG( 0.5 * dy / z0(j,i))158 wsvs = wsvs * ABS( wsvs )159 dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)160 ELSE161 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &162 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy163 ENDIF164 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)165 158 166 159 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & … … 172 165 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 173 166 174 ENDIF 175 176 ENDDO 177 178 ! 179 !-- 3 - Wird nur ausgefuehrt, wenn mindestens ein Niveau 180 !-- zwischen 2 und 4 liegt, d.h. ab einer Gebaeudemindest- 181 !-- hoehe von 2 dz. 'Nur Wand: Wall functions' 182 DO j = nys, nyn 183 184 IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) & 185 THEN 167 168 ! 169 !-- 3 - Wird nur ausgefuehrt, wenn mindestens ein Niveau 170 !-- zwischen 2 und 4 liegt, d.h. ab einer Gebaeudemindest- 171 !-- hoehe von 2 dz. Die Wall fluxes sind in diesem Fall 172 !-- schon unter (2) berechnet worden. 'Nur Wand: Wall functions' 186 173 187 174 DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2 188 175 189 176 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 177 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & 178 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 179 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 180 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & 181 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 182 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 183 190 184 IF ( wall_e_y(j,i) /= 0.0 ) THEN 191 usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) & 192 / LOG( 0.5 * dy / z0(j,i) ) 193 usvs = usvs * ABS( usvs ) 194 dudy = wall_e_y(j,i) * usvs / km(k,j,i) 185 dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i) 186 dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i) 195 187 ELSE 196 188 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & 197 189 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 190 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & 191 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 198 192 ENDIF 199 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &200 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)201 193 202 194 IF ( wall_e_x(j,i) /= 0.0 ) THEN 203 vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) & 204 / LOG( 0.5 * dx / z0(j,i)) 205 vsus = vsus * ABS( vsus ) 206 dvdx = wall_e_x(j,i) * vsus / km(k,j,i) 195 dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i) 196 dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i) 207 197 ELSE 208 198 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & 209 199 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 210 ENDIF211 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy212 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &213 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)214 215 IF ( wall_e_x(j,i) /= 0.0 ) THEN216 wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &217 / LOG( 0.5 * dx / z0(j,i))218 wsus = wsus * ABS( wsus )219 dwdx = wall_e_x(j,i) * wsus / km(k,j,i)220 ELSE221 200 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & 222 201 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 223 202 ENDIF 224 IF ( wall_e_y(j,i) /= 0.0 ) THEN225 wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &226 / LOG( 0.5 * dy / z0(j,i))227 wsvs = wsvs * ABS( wsvs )228 dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)229 ELSE230 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &231 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy232 ENDIF233 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)234 203 235 204 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & … … 423 392 DO j = nys, nyn 424 393 425 k = nzb_diff_s_inner(j,i) 394 k = nzb_diff_s_inner(j,i)-1 426 395 427 396 IF ( .NOT. cloud_physics ) THEN … … 505 474 506 475 REAL :: def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, & 507 k1, k2, theta, temp, usvs, vsus, wsus,wsvs 476 k1, k2, theta, temp 477 478 REAL, DIMENSION(nzb:nzt+1) :: usvs, vsus, wsus,wsvs 508 479 509 480 ! … … 549 520 550 521 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 522 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & 523 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 524 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 525 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & 526 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 527 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 528 551 529 IF ( wall_e_y(j,i) /= 0.0 ) THEN 552 usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) & 553 / LOG( 0.5 * dy / z0(j,i) ) 554 usvs = usvs * ABS( usvs ) 555 dudy = wall_e_y(j,i) * usvs / km(k,j,i) 530 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 531 usvs, 1.0, 0.0, 0.0, 0.0 ) 532 dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i) 533 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 534 wsvs, 0.0, 0.0, 1.0, 0.0 ) 535 dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i) 556 536 ELSE 557 537 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & 558 538 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 539 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & 540 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 559 541 ENDIF 560 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &561 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)562 542 563 543 IF ( wall_e_x(j,i) /= 0.0 ) THEN 564 vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) & 565 / LOG( 0.5 * dx / z0(j,i)) 566 vsus = vsus * ABS( vsus ) 567 dvdx = wall_e_x(j,i) * vsus / km(k,j,i) 544 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 545 vsus, 0.0, 1.0, 0.0, 0.0 ) 546 dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i) 547 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 548 wsus, 0.0, 0.0, 0.0, 1.0 ) 549 dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i) 568 550 ELSE 569 551 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & 570 552 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 571 ENDIF572 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy573 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &574 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)575 576 IF ( wall_e_x(j,i) /= 0.0 ) THEN577 wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &578 / LOG( 0.5 * dx / z0(j,i))579 wsus = wsus * ABS( wsus )580 dwdx = wall_e_x(j,i) * wsus / km(k,j,i)581 ELSE582 553 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & 583 554 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 584 555 ENDIF 585 IF ( wall_e_y(j,i) /= 0.0 ) THEN586 wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &587 / LOG( 0.5 * dy / z0(j,i))588 wsvs = wsvs * ABS( wsvs )589 dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)590 ELSE591 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &592 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy593 ENDIF594 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)595 556 596 557 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & … … 605 566 !-- 3 - Wird nur ausgefuehrt, wenn mindestens ein Niveau 606 567 !-- zwischen 2 und 4 liegt, d.h. ab einer Gebaeudemindest- 607 !-- hoehe von 2 dz. 'Nur Wand: Wall functions' 568 !-- hoehe von 2 dz. Di Wall fluxes sind in diesem Fall bereits vorab 569 !-- unter (2) berechnet worden. 'Nur Wand: Wall functions' 608 570 DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2 609 571 610 572 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 573 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & 574 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 575 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 576 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & 577 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 578 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 579 611 580 IF ( wall_e_y(j,i) /= 0.0 ) THEN 612 usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) & 613 / LOG( 0.5 * dy / z0(j,i) ) 614 usvs = usvs * ABS( usvs ) 615 dudy = wall_e_y(j,i) * usvs / km(k,j,i) 581 dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i) 582 dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i) 616 583 ELSE 617 584 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & 618 585 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 619 ENDIF620 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &621 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)586 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & 587 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 588 ENDIF 622 589 623 590 IF ( wall_e_x(j,i) /= 0.0 ) THEN 624 vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) & 625 / LOG( 0.5 * dx / z0(j,i)) 626 vsus = vsus * ABS( vsus ) 627 dvdx = wall_e_x(j,i) * vsus / km(k,j,i) 591 dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i) 592 dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i) 628 593 ELSE 629 594 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & 630 595 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 631 ENDIF632 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy633 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &634 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)635 636 IF ( wall_e_x(j,i) /= 0.0 ) THEN637 wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &638 / LOG( 0.5 * dx / z0(j,i))639 wsus = wsus * ABS( wsus )640 dwdx = wall_e_x(j,i) * wsus / km(k,j,i)641 ELSE642 596 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & 643 597 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 644 598 ENDIF 645 IF ( wall_e_y(j,i) /= 0.0 ) THEN646 wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &647 / LOG( 0.5 * dy / z0(j,i))648 wsvs = wsvs * ABS( wsvs )649 dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)650 ELSE651 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &652 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy653 ENDIF654 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)655 599 656 600 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &
Note: See TracChangeset
for help on using the changeset viewer.
