MODULE production_e_mod !------------------------------------------------------------------------------! ! Actual revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Log: production_e.f90,v $ ! Revision 1.21 2006/04/26 12:45:35 raasch ! OpenMP parallelization of production_e_init ! ! Revision 1.20 2006/02/23 12:51:25 raasch ! nzb_2d and nzb_diff_2d replaced by nzb_s_outer and nzb_diff_s_outer ! respectively, momentum flux calculation at vertical walls using profile ! functions, loop rearrangement for better vectorization ! ! Revision 1.19 2005/06/29 08:15:33 steinfeld ! Error in calculating u_0 and v_0 removed (multiplication instead of ! division by 2dz) ! ! Revision 1.18 2004/04/30 12:45:54 raasch ! dptdz eliminated ! ! Revision 1.17 2004/01/30 10:36:00 raasch ! Velocity gradients at the surface limited (see u_0, v_0), ! scalar lower k index nzb_diff replaced by 2d-array nzb_diff_2d ! ! Revision 1.16 2003/03/12 16:40:22 raasch ! Full code replaced in the call for all gridpoints instead of calling the ! _ij version (required by NEC, because otherwise no vectorization) ! ! Revision 1.15 2002/12/19 16:16:34 raasch ! Calculation of deformation tensor re-designed (most of the finite ! differences are now formed over two grid spacings, errors in derivatives ! along y removed) ! ! Revision 1.14 2002/09/12 13:09:44 raasch ! Error in calculating v_0 removed ! ! Revision 1.13 2002/06/11 13:19:04 raasch ! Former subroutine changed to a module which allows to be called for all grid ! points of a single vertical column with index i,j or for all grid points by ! using function overloading. ! Calculation of u_0 and v_0 moved to the new subroutine production_e_init ! due to the global communication necessary. These arrays are allocated only ! once during the first call. ! ! Revision 1.12 2001/08/21 09:59:47 raasch ! Special treatment at k=1 generally if surface fluxes are prescribed (not only ! in case of a Prandtl layer). In these cases, the boundary conditions are ! also included in the shear production term. ! ! Revision 1.11 2001/03/30 07:48:54 raasch ! Calculation of shear production and buoyancy production is split into ! two loops, arguments k1 and k2 are eliminated and now used as scalars, ! Translation of remaining German identifiers (variables, subroutines, etc.) ! ! Revision 1.10 2001/01/22 07:56:19 raasch ! Module test_variables removed ! ! Revision 1.9 2001/01/02 17:34:12 raasch ! -dpt_dz_d, dpt_dz_u ! ! Revision 1.8 2000/07/03 13:00:24 raasch ! Arguments k1 and k2 declared as pointers, ! all comments tranlated into English ! ! Revision 1.7 2000/04/18 08:13:00 schroeter ! Revision 1.5 rueckgaengig gemacht, Temperaturgrandient im ! Auftriebtrem wird wieder durch zentrale Differenzen gebildet ! ! Revision 1.6 2000/04/13 14:24:02 schroeter ! condsidering the influence of humidity to TKE-production ! ! Revision 1.5 99/02/17 09:29:18 09:29:18 raasch (Siegfried Raasch) ! Vertikalscherung des Windes exakter formuliert ! Temperaturgradient im Auftriebsterm enger gefasst, im stabilen Fall wird ! jetzt das Minimum vom oberen und unteren Differenzenquotienten verwendet ! ! Revision 1.4 1998/07/06 12:31:52 raasch ! + USE test_variables ! ! Revision 1.3 1998/04/15 11:23:49 raasch ! Im Fall einer Prandtl-Schicht wird die Energieproduktion bei nzb+1 ! jetzt immer mittels shf und sonst ueber dpt/dz berechnet (bisher nur bei ! vorgegebenem Waermestrom) ! ! Revision 1.2 1998/02/19 07:10:51 raasch ! vorgegebener Waermestrom wird gegebenenfalls im Auftriebsterm bei ! k=nzb+1 direkt angegeben ! ! Revision 1.1 1997/09/19 07:45:35 raasch ! Initial revision ! ! ! Description: ! ------------ ! Production terms (shear + buoyancy) of the TKE !------------------------------------------------------------------------------! PRIVATE PUBLIC production_e, production_e_init LOGICAL, SAVE :: first_call = .TRUE. REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: u_0, v_0 INTERFACE production_e MODULE PROCEDURE production_e MODULE PROCEDURE production_e_ij END INTERFACE production_e INTERFACE production_e_init MODULE PROCEDURE production_e_init END INTERFACE production_e_init CONTAINS !------------------------------------------------------------------------------! ! Call for all grid points !------------------------------------------------------------------------------! SUBROUTINE production_e USE arrays_3d USE cloud_parameters USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, j, k REAL :: def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, & k1, k2, theta, temp, usvs, vsus, wsus, wsvs ! !-- Calculate TKE production by shear DO i = nxl, nxr DO j = nys, nyn DO k = nzb_diff_s_outer(j,i), nzt-1 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) IF ( def < 0.0 ) def = 0.0 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def ENDDO ENDDO IF ( use_surface_fluxes ) THEN ! !-- Position neben Gebaeudewand !-- 2 - Wird immer ausgefuehrt. 'Boden und Wand: !-- u_0,v_0 und Wall functions' DO j = nys, nyn IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) & THEN k = nzb_diff_s_inner(j,i) - 1 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx IF ( wall_e_y(j,i) /= 0.0 ) THEN usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) & / LOG( 0.5 * dy / z0(j,i) ) usvs = usvs * ABS( usvs ) dudy = wall_e_y(j,i) * usvs / km(k,j,i) ELSE dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy ENDIF dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) IF ( wall_e_x(j,i) /= 0.0 ) THEN vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) & / LOG( 0.5 * dx / z0(j,i)) vsus = vsus * ABS( vsus ) dvdx = wall_e_x(j,i) * vsus / km(k,j,i) ELSE dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx ENDIF dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) IF ( wall_e_x(j,i) /= 0.0 ) THEN wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) & / LOG( 0.5 * dx / z0(j,i)) wsus = wsus * ABS( wsus ) dwdx = wall_e_x(j,i) * wsus / km(k,j,i) ELSE dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx ENDIF IF ( wall_e_y(j,i) /= 0.0 ) THEN wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) & / LOG( 0.5 * dy / z0(j,i)) wsvs = wsvs * ABS( wsvs ) dwdy = wall_e_y(j,i) * wsvs / km(k,j,i) ELSE dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy ENDIF dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) IF ( def < 0.0 ) def = 0.0 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def ENDIF ENDDO ! !-- 3 - Wird nur ausgefuehrt, wenn mindestens ein Niveau !-- zwischen 2 und 4 liegt, d.h. ab einer Gebaeudemindest- !-- hoehe von 2 dz. 'Nur Wand: Wall functions' DO j = nys, nyn IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) & THEN DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx IF ( wall_e_y(j,i) /= 0.0 ) THEN usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) & / LOG( 0.5 * dy / z0(j,i) ) usvs = usvs * ABS( usvs ) dudy = wall_e_y(j,i) * usvs / km(k,j,i) ELSE dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy ENDIF dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) IF ( wall_e_x(j,i) /= 0.0 ) THEN vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) & / LOG( 0.5 * dx / z0(j,i)) vsus = vsus * ABS( vsus ) dvdx = wall_e_x(j,i) * vsus / km(k,j,i) ELSE dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx ENDIF dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) IF ( wall_e_x(j,i) /= 0.0 ) THEN wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) & / LOG( 0.5 * dx / z0(j,i)) wsus = wsus * ABS( wsus ) dwdx = wall_e_x(j,i) * wsus / km(k,j,i) ELSE dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx ENDIF IF ( wall_e_y(j,i) /= 0.0 ) THEN wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) & / LOG( 0.5 * dy / z0(j,i)) wsvs = wsvs * ABS( wsvs ) dwdy = wall_e_y(j,i) * wsvs / km(k,j,i) ELSE dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy ENDIF dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) IF ( def < 0.0 ) def = 0.0 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def ENDDO ENDIF ENDDO ! !-- 4 - Wird immer ausgefuehrt. !-- 'Sonderfall: Freie Atmosphaere' (wie bei 0) DO j = nys, nyn IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) & THEN k = nzb_diff_s_outer(j,i)-1 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) IF ( def < 0.0 ) def = 0.0 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def ENDIF ENDDO ! !-- Position ohne angrenzende Gebaeudewand !-- 1 - Wird immer ausgefuehrt. 'Nur Boden: u_0,v_0' DO j = nys, nyn IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) & THEN k = nzb_diff_s_inner(j,i)-1 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) IF ( def < 0.0 ) def = 0.0 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def ENDIF ENDDO ENDIF ! !-- Calculate TKE production by buoyancy IF ( .NOT. moisture ) THEN DO j = nys, nyn DO k = nzb_diff_s_inner(j,i), nzt-1 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * & ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) ENDDO IF ( use_surface_fluxes ) THEN k = nzb_diff_s_inner(j,i)-1 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i) ENDIF ENDDO ELSE DO j = nys, nyn DO k = nzb_diff_s_inner(j,i), nzt-1 IF ( .NOT. cloud_physics ) THEN k1 = 1.0 + 0.61 * q(k,j,i) k2 = 0.61 * pt(k,j,i) ELSE IF ( ql(k,j,i) == 0.0 ) THEN k1 = 1.0 + 0.61 * q(k,j,i) k2 = 0.61 * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0 - q(k,j,i) + 1.61 * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0 + 0.622 * l_d_r / temp ) ) / & ( 1.0 + 0.622 * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) ENDIF ENDIF tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & ) * dd2zu(k) ENDDO ENDDO IF ( use_surface_fluxes ) THEN DO j = nys, nyn k = nzb_diff_s_inner(j,i)-1 IF ( .NOT. cloud_physics ) THEN k1 = 1.0 + 0.61 * q(k,j,i) k2 = 0.61 * pt(k,j,i) ELSE IF ( ql(k,j,i) == 0.0 ) THEN k1 = 1.0 + 0.61 * q(k,j,i) k2 = 0.61 * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0 - q(k,j,i) + 1.61 * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0 + 0.622 * l_d_r / temp ) ) / & ( 1.0 + 0.622 * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) ENDIF ENDIF tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & ( k1* shf(j,i) + k2 * qsws(j,i) ) ENDDO ENDIF ENDIF ENDDO END SUBROUTINE production_e !------------------------------------------------------------------------------! ! Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE production_e_ij( i, j ) USE arrays_3d USE cloud_parameters USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, j, k REAL :: def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, & k1, k2, theta, temp, usvs, vsus, wsus,wsvs ! !-- Calculate TKE production by shear DO k = nzb_diff_s_outer(j,i), nzt-1 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) & + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 & + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) IF ( def < 0.0 ) def = 0.0 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def ENDDO IF ( use_surface_fluxes ) THEN IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) THEN ! !-- Position neben Gebaeudewand !-- 2 - Wird immer ausgefuehrt. 'Boden und Wand: !-- u_0,v_0 und Wall functions' k = nzb_diff_s_inner(j,i)-1 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx IF ( wall_e_y(j,i) /= 0.0 ) THEN usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) & / LOG( 0.5 * dy / z0(j,i) ) usvs = usvs * ABS( usvs ) dudy = wall_e_y(j,i) * usvs / km(k,j,i) ELSE dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy ENDIF dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) IF ( wall_e_x(j,i) /= 0.0 ) THEN vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) & / LOG( 0.5 * dx / z0(j,i)) vsus = vsus * ABS( vsus ) dvdx = wall_e_x(j,i) * vsus / km(k,j,i) ELSE dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx ENDIF dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) IF ( wall_e_x(j,i) /= 0.0 ) THEN wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) & / LOG( 0.5 * dx / z0(j,i)) wsus = wsus * ABS( wsus ) dwdx = wall_e_x(j,i) * wsus / km(k,j,i) ELSE dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx ENDIF IF ( wall_e_y(j,i) /= 0.0 ) THEN wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) & / LOG( 0.5 * dy / z0(j,i)) wsvs = wsvs * ABS( wsvs ) dwdy = wall_e_y(j,i) * wsvs / km(k,j,i) ELSE dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy ENDIF dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) IF ( def < 0.0 ) def = 0.0 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def ! !-- 3 - Wird nur ausgefuehrt, wenn mindestens ein Niveau !-- zwischen 2 und 4 liegt, d.h. ab einer Gebaeudemindest- !-- hoehe von 2 dz. 'Nur Wand: Wall functions' DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx IF ( wall_e_y(j,i) /= 0.0 ) THEN usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) & / LOG( 0.5 * dy / z0(j,i) ) usvs = usvs * ABS( usvs ) dudy = wall_e_y(j,i) * usvs / km(k,j,i) ELSE dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy ENDIF dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) IF ( wall_e_x(j,i) /= 0.0 ) THEN vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) & / LOG( 0.5 * dx / z0(j,i)) vsus = vsus * ABS( vsus ) dvdx = wall_e_x(j,i) * vsus / km(k,j,i) ELSE dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx ENDIF dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) IF ( wall_e_x(j,i) /= 0.0 ) THEN wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) & / LOG( 0.5 * dx / z0(j,i)) wsus = wsus * ABS( wsus ) dwdx = wall_e_x(j,i) * wsus / km(k,j,i) ELSE dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx ENDIF IF ( wall_e_y(j,i) /= 0.0 ) THEN wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) & / LOG( 0.5 * dy / z0(j,i)) wsvs = wsvs * ABS( wsvs ) dwdy = wall_e_y(j,i) * wsvs / km(k,j,i) ELSE dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy ENDIF dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) IF ( def < 0.0 ) def = 0.0 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def ENDDO ! !-- 4 - Wird immer ausgefuehrt. !-- 'Sonderfall: Freie Atmosphaere' (wie bei 0) k = nzb_diff_s_outer(j,i)-1 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) IF ( def < 0.0 ) def = 0.0 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def ELSE ! !-- Position ohne angrenzende Gebaeudewand !-- 1 - Wird immer ausgefuehrt. 'Nur Boden: u_0,v_0' k = nzb_diff_s_inner(j,i)-1 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) & + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 & + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) IF ( def < 0.0 ) def = 0.0 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def ENDIF ENDIF ! !-- Calculate TKE production by buoyancy IF ( .NOT. moisture ) THEN DO k = nzb_diff_s_inner(j,i), nzt-1 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * & ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) ENDDO IF ( use_surface_fluxes ) THEN k = nzb_diff_s_inner(j,i)-1 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i) ENDIF ELSE DO k = nzb_diff_s_inner(j,i), nzt-1 IF ( .NOT. cloud_physics ) THEN k1 = 1.0 + 0.61 * q(k,j,i) k2 = 0.61 * pt(k,j,i) ELSE IF ( ql(k,j,i) == 0.0 ) THEN k1 = 1.0 + 0.61 * q(k,j,i) k2 = 0.61 * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0 - q(k,j,i) + 1.61 * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0 + 0.622 * l_d_r / temp ) ) / & ( 1.0 + 0.622 * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) ENDIF ENDIF tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & ) * dd2zu(k) ENDDO IF ( use_surface_fluxes ) THEN k = nzb_diff_s_inner(j,i)-1 IF ( .NOT. cloud_physics ) THEN k1 = 1.0 + 0.61 * q(k,j,i) k2 = 0.61 * pt(k,j,i) ELSE IF ( ql(k,j,i) == 0.0 ) THEN k1 = 1.0 + 0.61 * q(k,j,i) k2 = 0.61 * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0 - q(k,j,i) + 1.61 * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0 + 0.622 * l_d_r / temp ) ) / & ( 1.0 + 0.622 * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) ENDIF ENDIF tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & ( k1* shf(j,i) + k2 * qsws(j,i) ) ENDIF ENDIF END SUBROUTINE production_e_ij SUBROUTINE production_e_init USE arrays_3d USE control_parameters USE grid_variables USE indices IMPLICIT NONE INTEGER :: i, j, ku, kv IF ( use_surface_fluxes ) THEN IF ( first_call ) THEN ALLOCATE( u_0(nys-1:nyn+1,nxl-1:nxr+1), & v_0(nys-1:nyn+1,nxl-1:nxr+1) ) first_call = .FALSE. ENDIF ! !-- Calculate a virtual velocity at the surface in a way that the !-- vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the !-- Prandtl law (-w'u'/km). This gradient is used in the TKE shear !-- production term at k=1 (see production_e_ij). !-- The velocity gradient has to be limited in case of too small km !-- (otherwise the timestep may be significantly reduced by large !-- surface winds). !-- WARNING: the exact analytical solution would require the determination !-- of the eddy diffusivity by km = u* * kappa * zp / phi_m. !$OMP PARALLEL DO PRIVATE( ku, kv ) DO i = nxl, nxr DO j = nys, nyn ku = nzb_u_inner(j,i)+1 kv = nzb_v_inner(j,i)+1 u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / & ( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) + & 1.0E-20 ) ! ( us(j,i) * kappa * zu(1) ) v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / & ( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) + & 1.0E-20 ) ! ( us(j,i) * kappa * zu(1) ) IF ( ABS( u(ku+1,j,i) - u_0(j,i) ) > & ABS( u(ku+1,j,i) - u(ku-1,j,i) ) ) u_0(j,i) = u(ku-1,j,i) IF ( ABS( v(kv+1,j,i) - v_0(j,i) ) > & ABS( v(kv+1,j,i) - v(kv-1,j,i) ) ) v_0(j,i) = v(kv-1,j,i) ENDDO ENDDO CALL exchange_horiz_2d( u_0 ) CALL exchange_horiz_2d( v_0 ) ENDIF END SUBROUTINE production_e_init END MODULE production_e_mod