SUBROUTINE lpm_set_attributes !--------------------------------------------------------------------------------! ! This file is part of PALM. ! ! PALM is free software: you can redistribute it and/or modify it under the terms ! of the GNU General Public License as published by the Free Software Foundation, ! either version 3 of the License, or (at your option) any later version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 1997-2014 Leibniz Universitaet Hannover !--------------------------------------------------------------------------------! ! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: lpm_set_attributes.f90 1321 2014-03-20 09:40:40Z raasch $ ! ! 1320 2014-03-20 08:40:49Z raasch ! ONLY-attribute added to USE-statements, ! kind-parameters added to all INTEGER and REAL declaration statements, ! kinds are defined in new module kinds, ! revision history before 2012 removed, ! comment fields (!:) to be used for variable explanations added to ! all variable declaration statements ! ! 1318 2014-03-17 13:35:16Z raasch ! module interfaces removed ! ! 1036 2012-10-22 13:43:42Z raasch ! code put under GPL (PALM 3.9) ! ! 849 2012-03-15 10:35:09Z raasch ! routine renamed: set_particle_attributes -> lpm_set_attributes ! ! 828 2012-02-21 12:00:36Z raasch ! particle feature color renamed class ! ! 271 2009-03-26 00:47:14Z raasch ! Initial version ! ! Description: ! ------------ ! This routine sets certain particle attributes depending on the values that ! other PALM variables have at the current particle position. !------------------------------------------------------------------------------! USE arrays_3d, & ONLY: pt, u, v, w, zu, zw USE control_parameters, & ONLY: atmos_ocean_sign, u_gtrans, v_gtrans, dz USE cpulog, & ONLY: cpu_log, log_point_s USE dvrp_variables, & ONLY: color_interval, dvrp_colortable_entries_prt, dvrpsize_interval, & particle_color, particle_dvrpsize USE grid_variables, & ONLY: ddx, dx, ddy, dy USE indices, & ONLY: ngp_2dh, nxl, nxr, nyn, nys, nzb, nzt USE kinds USE particle_attributes, & ONLY: number_of_particles, offset_ocean_nzt, particles USE pegrid USE statistics, & ONLY: sums, sums_l IMPLICIT NONE INTEGER(iwp) :: i !: INTEGER(iwp) :: j !: INTEGER(iwp) :: k !: INTEGER(iwp) :: n !: REAL(wp) :: aa !: REAL(wp) :: absuv !: REAL(wp) :: bb !: REAL(wp) :: cc !: REAL(wp) :: dd !: REAL(wp) :: gg !: REAL(wp) :: height !: REAL(wp) :: pt_int !: REAL(wp) :: pt_int_l !: REAL(wp) :: pt_int_u !: REAL(wp) :: u_int !: REAL(wp) :: u_int_l !: REAL(wp) :: u_int_u !: REAL(wp) :: v_int !: REAL(wp) :: v_int_l !: REAL(wp) :: v_int_u !: REAL(wp) :: w_int !: REAL(wp) :: w_int_l !: REAL(wp) :: w_int_u !: REAL(wp) :: x !: REAL(wp) :: y !: CALL cpu_log( log_point_s(49), 'lpm_set_attributes', 'start' ) ! !-- Set particle color IF ( particle_color == 'absuv' ) THEN ! !-- Set particle color depending on the absolute value of the horizontal !-- velocity DO n = 1, number_of_particles ! !-- Interpolate u velocity-component, determine left, front, bottom !-- index of u-array i = ( particles(n)%x + 0.5 * dx ) * ddx j = particles(n)%y * ddy k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz & + offset_ocean_nzt ! only exact if equidistant ! !-- Interpolation of the velocity components in the xy-plane x = particles(n)%x + ( 0.5 - i ) * dx y = particles(n)%y - j * dy aa = x**2 + y**2 bb = ( dx - x )**2 + y**2 cc = x**2 + ( dy - y )**2 dd = ( dx - x )**2 + ( dy - y )**2 gg = aa + bb + cc + dd u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) & + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) & ) / ( 3.0 * gg ) - u_gtrans IF ( k+1 == nzt+1 ) THEN u_int = u_int_l ELSE u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) & + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) & ) / ( 3.0 * gg ) - u_gtrans u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * & ( u_int_u - u_int_l ) ENDIF ! !-- Same procedure for interpolation of the v velocity-component (adopt !-- index k from u velocity-component) i = particles(n)%x * ddx j = ( particles(n)%y + 0.5 * dy ) * ddy x = particles(n)%x - i * dx y = particles(n)%y + ( 0.5 - j ) * dy aa = x**2 + y**2 bb = ( dx - x )**2 + y**2 cc = x**2 + ( dy - y )**2 dd = ( dx - x )**2 + ( dy - y )**2 gg = aa + bb + cc + dd v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) & + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) & ) / ( 3.0 * gg ) - v_gtrans IF ( k+1 == nzt+1 ) THEN v_int = v_int_l ELSE v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) & + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) & ) / ( 3.0 * gg ) - v_gtrans v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * & ( v_int_u - v_int_l ) ENDIF absuv = SQRT( u_int**2 + v_int**2 ) ! !-- Limit values by the given interval and normalize to interval [0,1] absuv = MIN( absuv, color_interval(2) ) absuv = MAX( absuv, color_interval(1) ) absuv = ( absuv - color_interval(1) ) / & ( color_interval(2) - color_interval(1) ) ! !-- Number of available colors is defined in init_dvrp particles(n)%class = 1 + absuv * ( dvrp_colortable_entries_prt - 1 ) ENDDO ELSEIF ( particle_color == 'pt*' ) THEN ! !-- Set particle color depending on the resolved scale temperature !-- fluctuation. !-- First, calculate the horizontal average of the potential temperature !-- (This is also done in flow_statistics, but flow_statistics is called !-- after this routine.) sums_l(:,4,0) = 0.0 DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 sums_l(k,4,0) = sums_l(k,4,0) + pt(k,j,i) ENDDO ENDDO ENDDO #if defined( __parallel ) IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, & MPI_REAL, MPI_SUM, comm2d, ierr ) #else sums(:,4) = sums_l(:,4,0) #endif sums(:,4) = sums(:,4) / ngp_2dh(0) DO n = 1, number_of_particles ! !-- Interpolate temperature to the current particle position i = particles(n)%x * ddx j = particles(n)%y * ddy k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz & + offset_ocean_nzt ! only exact if equidistant x = particles(n)%x - i * dx y = particles(n)%y - j * dy aa = x**2 + y**2 bb = ( dx - x )**2 + y**2 cc = x**2 + ( dy - y )**2 dd = ( dx - x )**2 + ( dy - y )**2 gg = aa + bb + cc + dd pt_int_l = ( ( gg - aa ) * pt(k,j,i) + ( gg - bb ) * pt(k,j,i+1) & + ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) & ) / ( 3.0 * gg ) - sums(k,4) pt_int_u = ( ( gg-aa ) * pt(k+1,j,i) + ( gg-bb ) * pt(k+1,j,i+1) & + ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) & ) / ( 3.0 * gg ) - sums(k,4) pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * & ( pt_int_u - pt_int_l ) ! !-- Limit values by the given interval and normalize to interval [0,1] pt_int = MIN( pt_int, color_interval(2) ) pt_int = MAX( pt_int, color_interval(1) ) pt_int = ( pt_int - color_interval(1) ) / & ( color_interval(2) - color_interval(1) ) ! !-- Number of available colors is defined in init_dvrp particles(n)%class = 1 + pt_int * ( dvrp_colortable_entries_prt - 1 ) ENDDO ELSEIF ( particle_color == 'z' ) THEN ! !-- Set particle color depending on the height above the bottom !-- boundary (z=0) DO n = 1, number_of_particles height = particles(n)%z ! !-- Limit values by the given interval and normalize to interval [0,1] height = MIN( height, color_interval(2) ) height = MAX( height, color_interval(1) ) height = ( height - color_interval(1) ) / & ( color_interval(2) - color_interval(1) ) ! !-- Number of available colors is defined in init_dvrp particles(n)%class = 1 + height * ( dvrp_colortable_entries_prt - 1 ) ENDDO ENDIF ! !-- Set particle size for dvrp graphics IF ( particle_dvrpsize == 'absw' ) THEN DO n = 1, number_of_particles ! !-- Interpolate w-component to the current particle position i = particles(n)%x * ddx j = particles(n)%y * ddy k = particles(n)%z / dz x = particles(n)%x - i * dx y = particles(n)%y - j * dy aa = x**2 + y**2 bb = ( dx - x )**2 + y**2 cc = x**2 + ( dy - y )**2 dd = ( dx - x )**2 + ( dy - y )**2 gg = aa + bb + cc + dd w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) & + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) & ) / ( 3.0 * gg ) IF ( k+1 == nzt+1 ) THEN w_int = w_int_l ELSE w_int_u = ( ( gg-aa ) * w(k+1,j,i) + ( gg-bb ) * w(k+1,j,i+1) & + ( gg-cc ) * w(k+1,j+1,i) + ( gg-dd ) * w(k+1,j+1,i+1) & ) / ( 3.0 * gg ) w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * & ( w_int_u - w_int_l ) ENDIF ! !-- Limit values by the given interval and normalize to interval [0,1] w_int = ABS( w_int ) w_int = MIN( w_int, dvrpsize_interval(2) ) w_int = MAX( w_int, dvrpsize_interval(1) ) w_int = ( w_int - dvrpsize_interval(1) ) / & ( dvrpsize_interval(2) - dvrpsize_interval(1) ) particles(n)%dvrp_psize = ( 0.25 + w_int * 0.6 ) * dx ENDDO ENDIF CALL cpu_log( log_point_s(49), 'lpm_set_attributes', 'stop' ) END SUBROUTINE lpm_set_attributes