Changeset 57 for palm/trunk/SOURCE
- Timestamp:
- Mar 9, 2007 12:05:41 PM (17 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 19 edited
-
CURRENT_MODIFICATIONS (modified) (5 diffs)
-
advec_particles.f90 (modified) (13 diffs)
-
buoyancy.f90 (modified) (3 diffs)
-
check_parameters.f90 (modified) (2 diffs)
-
diffusion_e.f90 (modified) (3 diffs)
-
diffusion_u.f90 (modified) (5 diffs)
-
diffusion_v.f90 (modified) (5 diffs)
-
diffusion_w.f90 (modified) (5 diffs)
-
diffusivities.f90 (modified) (2 diffs)
-
header.f90 (modified) (2 diffs)
-
init_particles.f90 (modified) (4 diffs)
-
modules.f90 (modified) (4 diffs)
-
palm.f90 (modified) (4 diffs)
-
parin.f90 (modified) (2 diffs)
-
production_e.f90 (modified) (3 diffs)
-
prognostic_equations.f90 (modified) (10 diffs)
-
read_var_list.f90 (modified) (2 diffs)
-
user_interface.f90 (modified) (2 diffs)
-
write_var_list.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/CURRENT_MODIFICATIONS
r54 r57 1 1 New: 2 2 --- 3 4 particle reflection from vertical walls implemented, particle SGS model adjusted to walls 3 5 4 6 Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall. 5 7 6 8 new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files 9 10 new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used). 11 12 new user interface user_advec_particles 7 13 8 14 new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays … … 12 18 samples added to the user interface which show how to add user-define time series quantities. 13 19 14 Makefile, check_open, check_parameters, header, init_3d_model, modules, netcdf, parin, user_interface 20 unit 9 opened for debug output (file DEBUG_<pe#>) 21 22 Makefile, advec_particles, buoyancy, check_open, check_parameters, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, user_interface, write_var_list 15 23 16 24 New: wall_fluxes … … 19 27 Changed: 20 28 ------- 29 +age_m in particle_type 21 30 22 31 Move call of user_actions( 'after_integration' ) below increment of times … … 27 36 Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc. 28 37 29 check_parameters, data_output_ts, flow_statistics, init_3d_model, modules, parin, time_integration 38 __vtk directives removed from main program. 39 40 check_parameters, data_output_ts, flow_statistics, init_particles, init_3d_model, modules, palm, parin, time_integration 30 41 31 42 … … 35 46 Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model 36 47 48 Bugfix in sample for reading user defined data from restart file (user_init) 49 37 50 in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in 38 51 # case of .mod files 39 52 40 53 Makefile 41 init_1d_model 54 init_1d_model, user_interface -
palm/trunk/SOURCE/advec_particles.f90
r39 r57 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Particle reflection at vertical walls (by Jin Zhang), regarding vertical 7 ! walls in the SGS model, 8 ! + user interface user_advec_particles 6 9 ! TEST: PRINT statements on unit 9 (commented out) 7 10 ! … … 37 40 ! sub-timesteps (if sgs velocities are included), +3d-arrays de_dx/y/z, 38 41 ! izuf renamed iran, output of particle time series 39 !40 ! Revision 1.28 2006/02/23 09:41:21 raasch41 ! Only a fraction of the particles may have tails, therefore output for unit42 ! 85 and 90 has changed,43 ! sorting of particles in gridbox order after each timestep, optimization of44 ! the droplet growth by collision,45 ! 0.5 replaced by 0.4999999999 in setting the random fluctuations,46 ! psl, psr, pdx, etc. are now 1d arrays, i.e. they depend on the particle group,47 ! improve particle release at PE boundaries, nt_anz renamed48 ! current_timestep_number, idum (particle_type) renamed tail_id,49 ! error number argument for handle_netcdf_error,50 !51 ! Revision 1.27 2005/10/20 14:04:15 raasch52 ! droplet growth by collision completed (still has to be optimized for speed),53 ! new routine collision_efficiency (see end of this file),54 ! 2*r replaced by r in the exponential terms (momentum equation for particles),55 ! number of particles really used is additionally output on the netcdf particle56 ! data file,57 ! test print statement removed58 !59 ! Revision 1.26 2005/06/26 19:50:14 raasch60 ! Cloud droplet features added, radius is used instead of diameter,61 ! particle type initial data assignment extended,62 ! output on unit 80 for non-parallel runs63 !64 ! Revision 1.25 2005/05/18 14:53:11 raasch65 ! Extensions for NetCDF output66 !67 ! Revision 1.24 2005/04/23 08:24:58 raasch68 ! Revised calculation of output time counters regarding a possible decrease of69 ! the output time interval in case of restart runs70 !71 ! Revision 1.23 2005/03/26 14:20:05 raasch72 ! calculation of vertical particle velocity (with inertia) corrected, exp_arg73 ! had a wrong sign74 !75 ! Revision 1.22 2004/04/30 07:48:57 raasch76 ! Particle type initial data assignment extended77 !78 ! Revision 1.21 2004/01/28 15:00:48 raasch79 ! Output of particle infos in subroutine allocate_prt_memory on demand only80 !81 ! Revision 1.20 2003/10/29 08:38:52 raasch82 ! Module random_function_mod is used, modifications for new particle group83 ! feature, version number is written on binary file84 !85 ! Revision 1.19 2003/03/16 09:19:32 raasch86 ! Two underscores (_) are placed in front of all define-strings87 !88 ! Revision 1.18 2003/03/14 13:36:31 raasch89 ! Error in reflection boundary condition removed: particle speed must also90 ! be inverted91 !92 ! Revision 1.17 2003/03/04 11:23:20 raasch93 ! Error in particle inertia part removed (exp_arg must not contain the timestep)94 ! Particle velocities are also stored in particles in case of zero density ratio95 !96 ! Revision 1.16 2002/09/12 12:55:29 raasch97 ! Transport of particles with inertia implemented, write density_ratio to98 ! restart file99 !100 ! Revision 1.15 2002/06/11 12:11:56 raasch101 ! Declaration of u_int_l, u_int_u, v_int_l, v_int_u added102 !103 ! Revision 1.14 2002/05/02 18:46:34 18:46:34 raasch (Siegfried Raasch)104 ! Horizontal velocity components for particle advection are now interpolated105 ! between horizontal grid levels106 !107 ! Revision 1.13 2002/04/24 18:46:11 18:46:11 raasch (Siegfried Raasch)108 ! Temporary particle-arrays trlp, trnp, trrp and trsp are now allocatable arrays109 !110 ! Revision 1.12 2002/04/16 07:44:07 raasch111 ! Additional boundary conditions added, particle data for later analysis can be112 ! written on file113 !114 ! Revision 1.11 2001/11/09 13:38:20 raasch115 ! Colour information for particle tails is stored in array116 ! particle_tail_coordinates, adding information to the particle tails117 ! moved after call of user_particle_attributes118 !119 ! Revision 1.10 2001/08/21 08:19:48 raasch120 ! Storage of particle tails, precision of vertical advection improved,121 ! particle attributes (e.g. size and color) can be defined by the user122 !123 ! Revision 1.9 2001/07/12 12:01:27 raasch124 ! Random start positions of initial particles possible125 !126 ! Revision 1.8 2001/03/29 17:27:10 raasch127 ! Translation of remaining German identifiers (variables, subroutines, etc.)128 !129 ! Revision 1.7 2001/02/13 10:24:00 raasch130 ! Particle advection possible in case of galilei transformation131 !132 ! Revision 1.6 2001/01/25 06:51:04 raasch133 ! Module test_variables removed, writing of particle informations is optional134 ! now135 !136 ! Revision 1.5 2001/01/02 17:15:29 raasch137 ! Unit 80 is immediately closed after usage every time. Particle positions138 ! are written on unit 90 instead of 81 (subroutine write_particles).139 !140 ! Revision 1.4 2000/12/28 12:43:54 raasch141 ! Packing of particles-array is done explicitly (not by PACK-function, which142 ! seems to have problems when a large number of PEs is used)143 ! Missing translations included,144 ! code is used only optionally (cpp-directives are added)145 !146 ! Revision 1.3 2000/01/20 08:51:37 letzel147 ! All comments translated into English148 !149 ! Revision 1.2 1999/11/25 16:41:48 raasch150 ! Indexfehler im nichtparallelen Teil korrigiert151 42 ! 152 43 ! Revision 1.1 1999/11/25 16:16:06 raasch … … 693 584 DO j = nys, nyn 694 585 DO k = nzb, nzt+1 695 de_dx(k,j,i) = sgs_wfu_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx 696 de_dy(k,j,i) = sgs_wfv_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy 586 587 IF ( k <= nzb_s_inner(j,i-1) .AND. & 588 k > nzb_s_inner(j,i) .AND. & 589 k > nzb_s_inner(j,i+1) ) THEN 590 de_dx(k,j,i) = 2.0 * sgs_wfu_part * & 591 ( e(k,j,i+1) - e(k,j,i) ) * ddx 592 ELSEIF ( k > nzb_s_inner(j,i-1) .AND. & 593 k > nzb_s_inner(j,i) .AND. & 594 k <= nzb_s_inner(j,i+1) ) THEN 595 de_dx(k,j,i) = 2.0 * sgs_wfu_part * & 596 ( e(k,j,i) - e(k,j,i-1) ) * ddx 597 ELSEIF ( k < nzb_s_inner(j,i) .AND. k < nzb_s_inner(j,i+1) ) & 598 THEN 599 de_dx(k,j,i) = 0.0 600 ELSEIF ( k < nzb_s_inner(j,i-1) .AND. k < nzb_s_inner(j,i) ) & 601 THEN 602 de_dx(k,j,i) = 0.0 603 ELSE 604 de_dx(k,j,i) = sgs_wfu_part * & 605 ( e(k,j,i+1) - e(k,j,i-1) ) * ddx 606 ENDIF 607 608 IF ( k <= nzb_s_inner(j-1,i) .AND. & 609 k > nzb_s_inner(j,i) .AND. & 610 k > nzb_s_inner(j+1,i) ) THEN 611 de_dy(k,j,i) = 2.0 * sgs_wfv_part * & 612 ( e(k,j+1,i) - e(k,j,i) ) * ddy 613 ELSEIF ( k > nzb_s_inner(j-1,i) .AND. & 614 k > nzb_s_inner(j,i) .AND. & 615 k <= nzb_s_inner(j+1,i) ) THEN 616 de_dy(k,j,i) = 2.0 * sgs_wfv_part * & 617 ( e(k,j,i) - e(k,j-1,i) ) * ddy 618 ELSEIF ( k < nzb_s_inner(j,i) .AND. k < nzb_s_inner(j+1,i) ) & 619 THEN 620 de_dy(k,j,i) = 0.0 621 ELSEIF ( k < nzb_s_inner(j-1,i) .AND. k < nzb_s_inner(j,i) ) & 622 THEN 623 de_dy(k,j,i) = 0.0 624 ELSE 625 de_dy(k,j,i) = sgs_wfv_part * & 626 ( e(k,j+1,i) - e(k,j-1,i) ) * ddy 627 ENDIF 628 697 629 ENDDO 698 630 ENDDO … … 703 635 DO i = nxl, nxr 704 636 DO j = nys, nyn 705 DO k = nzb+2, nzt-1 637 638 DO k = nzb_s_inner(j,i)+2, nzt-1 706 639 de_dz(k,j,i) = 2.0 * sgs_wfw_part * & 707 640 ( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1)-zu(k-1) ) 708 641 ENDDO 709 de_dz(nzb,j,i) = 0.0 710 de_dz(nzb+1,j,i) = 2.0 * sgs_wfw_part * & 711 ( e(nzb+2,j,i) - e(nzb+1,j,i) ) / & 712 ( zu(nzb+2) - zu(nzb+1) ) 642 643 k = nzb_s_inner(j,i) 644 de_dz(nzb:k,j,i) = 0.0 645 de_dz(k+1,j,i) = 2.0 * sgs_wfw_part * ( e(k+2,j,i) - e(k+1,j,i) ) & 646 / ( zu(k+2) - zu(k+1) ) 713 647 de_dz(nzt,j,i) = 0.0 714 648 de_dz(nzt+1,j,i) = 0.0 715 649 ENDDO 716 ENDDO 650 ENDDO 717 651 718 652 ! … … 967 901 k = ( particles(n)%z + 0.5 * dz ) / dz ! only exact if eq.dist 968 902 969 x = particles(n)%x - i * dx 970 y = particles(n)%y - j * dy 971 aa = x**2 + y**2 972 bb = ( dx - x )**2 + y**2 973 cc = x**2 + ( dy - y )**2 974 dd = ( dx - x )**2 + ( dy - y )**2 975 gg = aa + bb + cc + dd 976 977 e_int_l = ( ( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) & 978 + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) & 979 ) / ( 3.0 * gg ) 980 981 IF ( k+1 == nzt+1 ) THEN 982 e_int = e_int_l 903 IF ( topography == 'flat' ) THEN 904 905 x = particles(n)%x - i * dx 906 y = particles(n)%y - j * dy 907 aa = x**2 + y**2 908 bb = ( dx - x )**2 + y**2 909 cc = x**2 + ( dy - y )**2 910 dd = ( dx - x )**2 + ( dy - y )**2 911 gg = aa + bb + cc + dd 912 913 e_int_l = ( ( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) & 914 + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) & 915 ) / ( 3.0 * gg ) 916 917 IF ( k+1 == nzt+1 ) THEN 918 e_int = e_int_l 919 ELSE 920 e_int_u = ( ( gg - aa ) * e(k+1,j,i) + & 921 ( gg - bb ) * e(k+1,j,i+1) + & 922 ( gg - cc ) * e(k+1,j+1,i) + & 923 ( gg - dd ) * e(k+1,j+1,i+1) & 924 ) / ( 3.0 * gg ) 925 e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * & 926 ( e_int_u - e_int_l ) 927 ENDIF 928 929 ! 930 !-- Interpolate the TKE gradient along x (adopt incides i,j,k and 931 !-- all position variables from above (TKE)) 932 de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + & 933 ( gg - bb ) * de_dx(k,j,i+1) + & 934 ( gg - cc ) * de_dx(k,j+1,i) + & 935 ( gg - dd ) * de_dx(k,j+1,i+1) & 936 ) / ( 3.0 * gg ) 937 938 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 939 de_dx_int = de_dx_int_l 940 ELSE 941 de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + & 942 ( gg - bb ) * de_dx(k+1,j,i+1) + & 943 ( gg - cc ) * de_dx(k+1,j+1,i) + & 944 ( gg - dd ) * de_dx(k+1,j+1,i+1) & 945 ) / ( 3.0 * gg ) 946 de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * & 947 ( de_dx_int_u - de_dx_int_l ) 948 ENDIF 949 950 ! 951 !-- Interpolate the TKE gradient along y 952 de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + & 953 ( gg - bb ) * de_dy(k,j,i+1) + & 954 ( gg - cc ) * de_dy(k,j+1,i) + & 955 ( gg - dd ) * de_dy(k,j+1,i+1) & 956 ) / ( 3.0 * gg ) 957 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 958 de_dy_int = de_dy_int_l 959 ELSE 960 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + & 961 ( gg - bb ) * de_dy(k+1,j,i+1) + & 962 ( gg - cc ) * de_dy(k+1,j+1,i) + & 963 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 964 ) / ( 3.0 * gg ) 965 de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * & 966 ( de_dy_int_u - de_dy_int_l ) 967 ENDIF 968 969 ! 970 !-- Interpolate the TKE gradient along z 971 IF ( particles(n)%z < 0.5 * dz ) THEN 972 de_dz_int = 0.0 973 ELSE 974 de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + & 975 ( gg - bb ) * de_dz(k,j,i+1) + & 976 ( gg - cc ) * de_dz(k,j+1,i) + & 977 ( gg - dd ) * de_dz(k,j+1,i+1) & 978 ) / ( 3.0 * gg ) 979 980 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 981 de_dz_int = de_dz_int_l 982 ELSE 983 de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + & 984 ( gg - bb ) * de_dz(k+1,j,i+1) + & 985 ( gg - cc ) * de_dz(k+1,j+1,i) + & 986 ( gg - dd ) * de_dz(k+1,j+1,i+1) & 987 ) / ( 3.0 * gg ) 988 de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * & 989 ( de_dz_int_u - de_dz_int_l ) 990 ENDIF 991 ENDIF 992 993 ! 994 !-- Interpolate the dissipation of TKE 995 diss_int_l = ( ( gg - aa ) * diss(k,j,i) + & 996 ( gg - bb ) * diss(k,j,i+1) + & 997 ( gg - cc ) * diss(k,j+1,i) + & 998 ( gg - dd ) * diss(k,j+1,i+1) & 999 ) / ( 3.0 * gg ) 1000 1001 IF ( k+1 == nzt+1 ) THEN 1002 diss_int = diss_int_l 1003 ELSE 1004 diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + & 1005 ( gg - bb ) * diss(k+1,j,i+1) + & 1006 ( gg - cc ) * diss(k+1,j+1,i) + & 1007 ( gg - dd ) * diss(k+1,j+1,i+1) & 1008 ) / ( 3.0 * gg ) 1009 diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * & 1010 ( diss_int_u - diss_int_l ) 1011 ENDIF 1012 983 1013 ELSE 984 e_int_u = ( ( gg - aa ) * e(k+1,j,i) + & 985 ( gg - bb ) * e(k+1,j,i+1) + & 986 ( gg - cc ) * e(k+1,j+1,i) + & 987 ( gg - dd ) * e(k+1,j+1,i+1) & 988 ) / ( 3.0 * gg ) 989 e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * & 990 ( e_int_u - e_int_l ) 991 ENDIF 992 993 ! 994 !-- Interpolate the TKE gradient along x (adopt incides i,j,k and all 995 !-- position variables from above (TKE)) 996 de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + & 997 ( gg - bb ) * de_dx(k,j,i+1) + & 998 ( gg - cc ) * de_dx(k,j+1,i) + & 999 ( gg - dd ) * de_dx(k,j+1,i+1) & 1000 ) / ( 3.0 * gg ) 1001 1002 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1003 de_dx_int = de_dx_int_l 1004 ELSE 1005 de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + & 1006 ( gg - bb ) * de_dx(k+1,j,i+1) + & 1007 ( gg - cc ) * de_dx(k+1,j+1,i) + & 1008 ( gg - dd ) * de_dx(k+1,j+1,i+1) & 1009 ) / ( 3.0 * gg ) 1010 de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * & 1011 ( de_dx_int_u - de_dx_int_l ) 1012 ENDIF 1013 1014 ! 1015 !-- Interpolate the TKE gradient along y 1016 de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + & 1017 ( gg - bb ) * de_dy(k,j,i+1) + & 1018 ( gg - cc ) * de_dy(k,j+1,i) + & 1019 ( gg - dd ) * de_dy(k,j+1,i+1) & 1020 ) / ( 3.0 * gg ) 1021 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1022 de_dy_int = de_dy_int_l 1023 ELSE 1024 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + & 1025 ( gg - bb ) * de_dy(k+1,j,i+1) + & 1026 ( gg - cc ) * de_dy(k+1,j+1,i) + & 1027 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 1028 ) / ( 3.0 * gg ) 1029 de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * & 1030 ( de_dy_int_u - de_dy_int_l ) 1031 ENDIF 1032 1033 ! 1034 !-- Interpolate the TKE gradient along z 1035 de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + & 1036 ( gg - bb ) * de_dz(k,j,i+1) + & 1037 ( gg - cc ) * de_dz(k,j+1,i) + & 1038 ( gg - dd ) * de_dz(k,j+1,i+1) & 1039 ) / ( 3.0 * gg ) 1040 1041 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1042 de_dz_int = de_dz_int_l 1043 ELSE 1044 de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + & 1045 ( gg - bb ) * de_dz(k+1,j,i+1) + & 1046 ( gg - cc ) * de_dz(k+1,j+1,i) + & 1047 ( gg - dd ) * de_dz(k+1,j+1,i+1) & 1048 ) / ( 3.0 * gg ) 1049 de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * & 1050 ( de_dz_int_u - de_dz_int_l ) 1051 ENDIF 1052 1053 ! 1054 !-- Interpolate the dissipation of TKE 1055 diss_int_l = ( ( gg - aa ) * diss(k,j,i) + & 1056 ( gg - bb ) * diss(k,j,i+1) + & 1057 ( gg - cc ) * diss(k,j+1,i) + & 1058 ( gg - dd ) * diss(k,j+1,i+1) & 1059 ) / ( 3.0 * gg ) 1060 1061 IF ( k+1 == nzt+1 ) THEN 1062 diss_int = diss_int_l 1063 ELSE 1064 diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + & 1065 ( gg - bb ) * diss(k+1,j,i+1) + & 1066 ( gg - cc ) * diss(k+1,j+1,i) + & 1067 ( gg - dd ) * diss(k+1,j+1,i+1) & 1068 ) / ( 3.0 * gg ) 1069 diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * & 1070 ( diss_int_u - diss_int_l ) 1071 ENDIF 1014 1015 ! 1016 !-- In case that there are buildings it has to be determined 1017 !-- how many of the gridpoints defining the particle box are 1018 !-- situated within a building 1019 !-- gp_outside_of_building(1): i,j,k 1020 !-- gp_outside_of_building(2): i,j+1,k 1021 !-- gp_outside_of_building(3): i,j,k+1 1022 !-- gp_outside_of_building(4): i,j+1,k+1 1023 !-- gp_outside_of_building(5): i+1,j,k 1024 !-- gp_outside_of_building(6): i+1,j+1,k 1025 !-- gp_outside_of_building(7): i+1,j,k+1 1026 !-- gp_outside_of_building(8): i+1,j+1,k+1 1027 1028 gp_outside_of_building = 0 1029 location = 0.0 1030 num_gp = 0 1031 1032 IF ( k > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN 1033 num_gp = num_gp + 1 1034 gp_outside_of_building(1) = 1 1035 location(num_gp,1) = i * dx 1036 location(num_gp,2) = j * dx 1037 location(num_gp,3) = k * dz - 0.5 * dz 1038 ei(num_gp) = e(k,j,i) 1039 dissi(num_gp) = diss(k,j,i) 1040 de_dxi(num_gp) = de_dx(k,j,i) 1041 de_dyi(num_gp) = de_dy(k,j,i) 1042 de_dzi(num_gp) = de_dz(k,j,i) 1043 ENDIF 1044 1045 IF ( k > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) & 1046 THEN 1047 num_gp = num_gp + 1 1048 gp_outside_of_building(2) = 1 1049 location(num_gp,1) = i * dx 1050 location(num_gp,2) = (j+1) * dx 1051 location(num_gp,3) = k * dz - 0.5 * dz 1052 ei(num_gp) = e(k,j+1,i) 1053 dissi(num_gp) = diss(k,j+1,i) 1054 de_dxi(num_gp) = de_dx(k,j+1,i) 1055 de_dyi(num_gp) = de_dy(k,j+1,i) 1056 de_dzi(num_gp) = de_dz(k,j+1,i) 1057 ENDIF 1058 1059 IF ( k+1 > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN 1060 num_gp = num_gp + 1 1061 gp_outside_of_building(3) = 1 1062 location(num_gp,1) = i * dx 1063 location(num_gp,2) = j * dy 1064 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1065 ei(num_gp) = e(k+1,j,i) 1066 dissi(num_gp) = diss(k+1,j,i) 1067 de_dxi(num_gp) = de_dx(k+1,j,i) 1068 de_dyi(num_gp) = de_dy(k+1,j,i) 1069 de_dzi(num_gp) = de_dz(k+1,j,i) 1070 ENDIF 1071 1072 IF ( k+1 > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) & 1073 THEN 1074 num_gp = num_gp + 1 1075 gp_outside_of_building(4) = 1 1076 location(num_gp,1) = i * dx 1077 location(num_gp,2) = (j+1) * dy 1078 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1079 ei(num_gp) = e(k+1,j+1,i) 1080 dissi(num_gp) = diss(k+1,j+1,i) 1081 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1082 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1083 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1084 ENDIF 1085 1086 IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) & 1087 THEN 1088 num_gp = num_gp + 1 1089 gp_outside_of_building(5) = 1 1090 location(num_gp,1) = (i+1) * dx 1091 location(num_gp,2) = j * dy 1092 location(num_gp,3) = k * dz - 0.5 * dz 1093 ei(num_gp) = e(k,j,i+1) 1094 dissi(num_gp) = diss(k,j,i+1) 1095 de_dxi(num_gp) = de_dx(k,j,i+1) 1096 de_dyi(num_gp) = de_dy(k,j,i+1) 1097 de_dzi(num_gp) = de_dz(k,j,i+1) 1098 ENDIF 1099 1100 IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) & 1101 THEN 1102 num_gp = num_gp + 1 1103 gp_outside_of_building(6) = 1 1104 location(num_gp,1) = (i+1) * dx 1105 location(num_gp,2) = (j+1) * dy 1106 location(num_gp,3) = k * dz - 0.5 * dz 1107 ei(num_gp) = e(k,j+1,i+1) 1108 dissi(num_gp) = diss(k,j+1,i+1) 1109 de_dxi(num_gp) = de_dx(k,j+1,i+1) 1110 de_dyi(num_gp) = de_dy(k,j+1,i+1) 1111 de_dzi(num_gp) = de_dz(k,j+1,i+1) 1112 ENDIF 1113 1114 IF ( k+1 > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) & 1115 THEN 1116 num_gp = num_gp + 1 1117 gp_outside_of_building(7) = 1 1118 location(num_gp,1) = (i+1) * dx 1119 location(num_gp,2) = j * dy 1120 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1121 ei(num_gp) = e(k+1,j,i+1) 1122 dissi(num_gp) = diss(k+1,j,i+1) 1123 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1124 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1125 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1126 ENDIF 1127 1128 IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)& 1129 THEN 1130 num_gp = num_gp + 1 1131 gp_outside_of_building(8) = 1 1132 location(num_gp,1) = (i+1) * dx 1133 location(num_gp,2) = (j+1) * dy 1134 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1135 ei(num_gp) = e(k+1,j+1,i+1) 1136 dissi(num_gp) = diss(k+1,j+1,i+1) 1137 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1138 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 1139 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1140 ENDIF 1141 1142 ! 1143 !-- If all gridpoints are situated outside of a building, then the 1144 !-- ordinary interpolation scheme can be used. 1145 IF ( num_gp == 8 ) THEN 1146 1147 x = particles(n)%x - i * dx 1148 y = particles(n)%y - j * dy 1149 aa = x**2 + y**2 1150 bb = ( dx - x )**2 + y**2 1151 cc = x**2 + ( dy - y )**2 1152 dd = ( dx - x )**2 + ( dy - y )**2 1153 gg = aa + bb + cc + dd 1154 1155 e_int_l = (( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) & 1156 + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1)& 1157 ) / ( 3.0 * gg ) 1158 1159 IF ( k+1 == nzt+1 ) THEN 1160 e_int = e_int_l 1161 ELSE 1162 e_int_u = ( ( gg - aa ) * e(k+1,j,i) + & 1163 ( gg - bb ) * e(k+1,j,i+1) + & 1164 ( gg - cc ) * e(k+1,j+1,i) + & 1165 ( gg - dd ) * e(k+1,j+1,i+1) & 1166 ) / ( 3.0 * gg ) 1167 e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * & 1168 ( e_int_u - e_int_l ) 1169 ENDIF 1170 1171 ! 1172 !-- Interpolate the TKE gradient along x (adopt incides i,j,k 1173 !-- and all position variables from above (TKE)) 1174 de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + & 1175 ( gg - bb ) * de_dx(k,j,i+1) + & 1176 ( gg - cc ) * de_dx(k,j+1,i) + & 1177 ( gg - dd ) * de_dx(k,j+1,i+1) & 1178 ) / ( 3.0 * gg ) 1179 1180 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1181 de_dx_int = de_dx_int_l 1182 ELSE 1183 de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + & 1184 ( gg - bb ) * de_dx(k+1,j,i+1) + & 1185 ( gg - cc ) * de_dx(k+1,j+1,i) + & 1186 ( gg - dd ) * de_dx(k+1,j+1,i+1) & 1187 ) / ( 3.0 * gg ) 1188 de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / & 1189 dz * ( de_dx_int_u - de_dx_int_l ) 1190 ENDIF 1191 1192 ! 1193 !-- Interpolate the TKE gradient along y 1194 de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + & 1195 ( gg - bb ) * de_dy(k,j,i+1) + & 1196 ( gg - cc ) * de_dy(k,j+1,i) + & 1197 ( gg - dd ) * de_dy(k,j+1,i+1) & 1198 ) / ( 3.0 * gg ) 1199 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1200 de_dy_int = de_dy_int_l 1201 ELSE 1202 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + & 1203 ( gg - bb ) * de_dy(k+1,j,i+1) + & 1204 ( gg - cc ) * de_dy(k+1,j+1,i) + & 1205 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 1206 ) / ( 3.0 * gg ) 1207 de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / & 1208 dz * ( de_dy_int_u - de_dy_int_l ) 1209 ENDIF 1210 1211 ! 1212 !-- Interpolate the TKE gradient along z 1213 IF ( particles(n)%z < 0.5 * dz ) THEN 1214 de_dz_int = 0.0 1215 ELSE 1216 de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + & 1217 ( gg - bb ) * de_dz(k,j,i+1) + & 1218 ( gg - cc ) * de_dz(k,j+1,i) + & 1219 ( gg - dd ) * de_dz(k,j+1,i+1) & 1220 ) / ( 3.0 * gg ) 1221 1222 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1223 de_dz_int = de_dz_int_l 1224 ELSE 1225 de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + & 1226 ( gg - bb ) * de_dz(k+1,j,i+1) + & 1227 ( gg - cc ) * de_dz(k+1,j+1,i) + & 1228 ( gg - dd ) * de_dz(k+1,j+1,i+1) & 1229 ) / ( 3.0 * gg ) 1230 de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) /& 1231 dz * ( de_dz_int_u - de_dz_int_l ) 1232 ENDIF 1233 ENDIF 1234 1235 ! 1236 !-- Interpolate the dissipation of TKE 1237 diss_int_l = ( ( gg - aa ) * diss(k,j,i) + & 1238 ( gg - bb ) * diss(k,j,i+1) + & 1239 ( gg - cc ) * diss(k,j+1,i) + & 1240 ( gg - dd ) * diss(k,j+1,i+1) & 1241 ) / ( 3.0 * gg ) 1242 1243 IF ( k+1 == nzt+1 ) THEN 1244 diss_int = diss_int_l 1245 ELSE 1246 diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + & 1247 ( gg - bb ) * diss(k+1,j,i+1) + & 1248 ( gg - cc ) * diss(k+1,j+1,i) + & 1249 ( gg - dd ) * diss(k+1,j+1,i+1) & 1250 ) / ( 3.0 * gg ) 1251 diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz *& 1252 ( diss_int_u - diss_int_l ) 1253 ENDIF 1254 1255 ELSE 1256 1257 ! 1258 !-- If wall between gridpoint 1 and gridpoint 5, then 1259 !-- Neumann boundary condition has to be applied 1260 IF ( gp_outside_of_building(1) == 1 .AND. & 1261 gp_outside_of_building(5) == 0 ) THEN 1262 num_gp = num_gp + 1 1263 location(num_gp,1) = i * dx + 0.5 * dx 1264 location(num_gp,2) = j * dy 1265 location(num_gp,3) = k * dz - 0.5 * dz 1266 ei(num_gp) = e(k,j,i) 1267 dissi(num_gp) = diss(k,j,i) 1268 de_dxi(num_gp) = 0.0 1269 de_dyi(num_gp) = de_dy(k,j,i) 1270 de_dzi(num_gp) = de_dz(k,j,i) 1271 ENDIF 1272 1273 IF ( gp_outside_of_building(5) == 1 .AND. & 1274 gp_outside_of_building(1) == 0 ) THEN 1275 num_gp = num_gp + 1 1276 location(num_gp,1) = i * dx + 0.5 * dx 1277 location(num_gp,2) = j * dy 1278 location(num_gp,3) = k * dz - 0.5 * dz 1279 ei(num_gp) = e(k,j,i+1) 1280 dissi(num_gp) = diss(k,j,i+1) 1281 de_dxi(num_gp) = 0.0 1282 de_dyi(num_gp) = de_dy(k,j,i+1) 1283 de_dzi(num_gp) = de_dz(k,j,i+1) 1284 ENDIF 1285 1286 ! 1287 !-- If wall between gridpoint 5 and gridpoint 6, then 1288 !-- then Neumann boundary condition has to be applied 1289 IF ( gp_outside_of_building(5) == 1 .AND. & 1290 gp_outside_of_building(6) == 0 ) THEN 1291 num_gp = num_gp + 1 1292 location(num_gp,1) = (i+1) * dx 1293 location(num_gp,2) = j * dy + 0.5 * dy 1294 location(num_gp,3) = k * dz - 0.5 * dz 1295 ei(num_gp) = e(k,j,i+1) 1296 dissi(num_gp) = diss(k,j,i+1) 1297 de_dxi(num_gp) = de_dx(k,j,i+1) 1298 de_dyi(num_gp) = 0.0 1299 de_dzi(num_gp) = de_dz(k,j,i+1) 1300 ENDIF 1301 1302 IF ( gp_outside_of_building(6) == 1 .AND. & 1303 gp_outside_of_building(5) == 0 ) THEN 1304 num_gp = num_gp + 1 1305 location(num_gp,1) = (i+1) * dx 1306 location(num_gp,2) = j * dy + 0.5 * dy 1307 location(num_gp,3) = k * dz - 0.5 * dz 1308 ei(num_gp) = e(k,j+1,i+1) 1309 dissi(num_gp) = diss(k,j+1,i+1) 1310 de_dxi(num_gp) = de_dx(k,j+1,i+1) 1311 de_dyi(num_gp) = 0.0 1312 de_dzi(num_gp) = de_dz(k,j+1,i+1) 1313 ENDIF 1314 1315 ! 1316 !-- If wall between gridpoint 2 and gridpoint 6, then 1317 !-- Neumann boundary condition has to be applied 1318 IF ( gp_outside_of_building(2) == 1 .AND. & 1319 gp_outside_of_building(6) == 0 ) THEN 1320 num_gp = num_gp + 1 1321 location(num_gp,1) = i * dx + 0.5 * dx 1322 location(num_gp,2) = (j+1) * dy 1323 location(num_gp,3) = k * dz - 0.5 * dz 1324 ei(num_gp) = e(k,j+1,i) 1325 dissi(num_gp) = diss(k,j+1,i) 1326 de_dxi(num_gp) = 0.0 1327 de_dyi(num_gp) = de_dy(k,j+1,i) 1328 de_dzi(num_gp) = de_dz(k,j+1,i) 1329 ENDIF 1330 1331 IF ( gp_outside_of_building(6) == 1 .AND. & 1332 gp_outside_of_building(2) == 0 ) THEN 1333 num_gp = num_gp + 1 1334 location(num_gp,1) = i * dx + 0.5 * dx 1335 location(num_gp,2) = (j+1) * dy 1336 location(num_gp,3) = k * dz - 0.5 * dz 1337 ei(num_gp) = e(k,j+1,i+1) 1338 dissi(num_gp) = diss(k,j+1,i+1) 1339 de_dxi(num_gp) = 0.0 1340 de_dyi(num_gp) = de_dy(k,j+1,i+1) 1341 de_dzi(num_gp) = de_dz(k,j+1,i+1) 1342 ENDIF 1343 1344 ! 1345 !-- If wall between gridpoint 1 and gridpoint 2, then 1346 !-- Neumann boundary condition has to be applied 1347 IF ( gp_outside_of_building(1) == 1 .AND. & 1348 gp_outside_of_building(2) == 0 ) THEN 1349 num_gp = num_gp + 1 1350 location(num_gp,1) = i * dx 1351 location(num_gp,2) = j * dy + 0.5 * dy 1352 location(num_gp,3) = k * dz - 0.5 * dz 1353 ei(num_gp) = e(k,j,i) 1354 dissi(num_gp) = diss(k,j,i) 1355 de_dxi(num_gp) = de_dx(k,j,i) 1356 de_dyi(num_gp) = 0.0 1357 de_dzi(num_gp) = de_dz(k,j,i) 1358 ENDIF 1359 1360 IF ( gp_outside_of_building(2) == 1 .AND. & 1361 gp_outside_of_building(1) == 0 ) THEN 1362 num_gp = num_gp + 1 1363 location(num_gp,1) = i * dx 1364 location(num_gp,2) = j * dy + 0.5 * dy 1365 location(num_gp,3) = k * dz - 0.5 * dz 1366 ei(num_gp) = e(k,j+1,i) 1367 dissi(num_gp) = diss(k,j+1,i) 1368 de_dxi(num_gp) = de_dx(k,j+1,i) 1369 de_dyi(num_gp) = 0.0 1370 de_dzi(num_gp) = de_dz(k,j+1,i) 1371 ENDIF 1372 1373 ! 1374 !-- If wall between gridpoint 3 and gridpoint 7, then 1375 !-- Neumann boundary condition has to be applied 1376 IF ( gp_outside_of_building(3) == 1 .AND. & 1377 gp_outside_of_building(7) == 0 ) THEN 1378 num_gp = num_gp + 1 1379 location(num_gp,1) = i * dx + 0.5 * dx 1380 location(num_gp,2) = j * dy 1381 location(num_gp,3) = k * dz + 0.5 * dz 1382 ei(num_gp) = e(k+1,j,i) 1383 dissi(num_gp) = diss(k+1,j,i) 1384 de_dxi(num_gp) = 0.0 1385 de_dyi(num_gp) = de_dy(k+1,j,i) 1386 de_dzi(num_gp) = de_dz(k+1,j,i) 1387 ENDIF 1388 1389 IF ( gp_outside_of_building(7) == 1 .AND. & 1390 gp_outside_of_building(3) == 0 ) THEN 1391 num_gp = num_gp + 1 1392 location(num_gp,1) = i * dx + 0.5 * dx 1393 location(num_gp,2) = j * dy 1394 location(num_gp,3) = k * dz + 0.5 * dz 1395 ei(num_gp) = e(k+1,j,i+1) 1396 dissi(num_gp) = diss(k+1,j,i+1) 1397 de_dxi(num_gp) = 0.0 1398 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1399 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1400 ENDIF 1401 1402 ! 1403 !-- If wall between gridpoint 7 and gridpoint 8, then 1404 !-- Neumann boundary condition has to be applied 1405 IF ( gp_outside_of_building(7) == 1 .AND. & 1406 gp_outside_of_building(8) == 0 ) THEN 1407 num_gp = num_gp + 1 1408 location(num_gp,1) = (i+1) * dx 1409 location(num_gp,2) = j * dy + 0.5 * dy 1410 location(num_gp,3) = k * dz + 0.5 * dz 1411 ei(num_gp) = e(k+1,j,i+1) 1412 dissi(num_gp) = diss(k+1,j,i+1) 1413 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1414 de_dyi(num_gp) = 0.0 1415 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1416 ENDIF 1417 1418 IF ( gp_outside_of_building(8) == 1 .AND. & 1419 gp_outside_of_building(7) == 0 ) THEN 1420 num_gp = num_gp + 1 1421 location(num_gp,1) = (i+1) * dx 1422 location(num_gp,2) = j * dy + 0.5 * dy 1423 location(num_gp,3) = k * dz + 0.5 * dz 1424 ei(num_gp) = e(k+1,j+1,i+1) 1425 dissi(num_gp) = diss(k+1,j+1,i+1) 1426 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1427 de_dyi(num_gp) = 0.0 1428 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1429 ENDIF 1430 1431 ! 1432 !-- If wall between gridpoint 4 and gridpoint 8, then 1433 !-- Neumann boundary condition has to be applied 1434 IF ( gp_outside_of_building(4) == 1 .AND. & 1435 gp_outside_of_building(8) == 0 ) THEN 1436 num_gp = num_gp + 1 1437 location(num_gp,1) = i * dx + 0.5 * dx 1438 location(num_gp,2) = (j+1) * dy 1439 location(num_gp,3) = k * dz + 0.5 * dz 1440 ei(num_gp) = e(k+1,j+1,i) 1441 dissi(num_gp) = diss(k+1,j+1,i) 1442 de_dxi(num_gp) = 0.0 1443 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1444 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1445 ENDIF 1446 1447 IF ( gp_outside_of_building(8) == 1 .AND. & 1448 gp_outside_of_building(4) == 0 ) THEN 1449 num_gp = num_gp + 1 1450 location(num_gp,1) = i * dx + 0.5 * dx 1451 location(num_gp,2) = (j+1) * dy 1452 location(num_gp,3) = k * dz + 0.5 * dz 1453 ei(num_gp) = e(k+1,j+1,i+1) 1454 dissi(num_gp) = diss(k+1,j+1,i+1) 1455 de_dxi(num_gp) = 0.0 1456 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 1457 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1458 ENDIF 1459 1460 ! 1461 !-- If wall between gridpoint 3 and gridpoint 4, then 1462 !-- Neumann boundary condition has to be applied 1463 IF ( gp_outside_of_building(3) == 1 .AND. & 1464 gp_outside_of_building(4) == 0 ) THEN 1465 num_gp = num_gp + 1 1466 location(num_gp,1) = i * dx 1467 location(num_gp,2) = j * dy + 0.5 * dy 1468 location(num_gp,3) = k * dz + 0.5 * dz 1469 ei(num_gp) = e(k+1,j,i) 1470 dissi(num_gp) = diss(k+1,j,i) 1471 de_dxi(num_gp) = de_dx(k+1,j,i) 1472 de_dyi(num_gp) = 0.0 1473 de_dzi(num_gp) = de_dz(k+1,j,i) 1474 ENDIF 1475 1476 IF ( gp_outside_of_building(4) == 1 .AND. & 1477 gp_outside_of_building(3) == 0 ) THEN 1478 num_gp = num_gp + 1 1479 location(num_gp,1) = i * dx 1480 location(num_gp,2) = j * dy + 0.5 * dy 1481 location(num_gp,3) = k * dz + 0.5 * dz 1482 ei(num_gp) = e(k+1,j+1,i) 1483 dissi(num_gp) = diss(k+1,j+1,i) 1484 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1485 de_dyi(num_gp) = 0.0 1486 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1487 ENDIF 1488 1489 ! 1490 !-- If wall between gridpoint 1 and gridpoint 3, then 1491 !-- Neumann boundary condition has to be applied 1492 !-- (only one case as only building beneath is possible) 1493 IF ( gp_outside_of_building(1) == 1 .AND. & 1494 gp_outside_of_building(3) == 0 ) THEN 1495 num_gp = num_gp + 1 1496 location(num_gp,1) = i * dx 1497 location(num_gp,2) = j * dy 1498 location(num_gp,3) = k * dz 1499 ei(num_gp) = e(k+1,j,i) 1500 dissi(num_gp) = diss(k+1,j,i) 1501 de_dxi(num_gp) = de_dx(k+1,j,i) 1502 de_dyi(num_gp) = de_dy(k+1,j,i) 1503 de_dzi(num_gp) = 0.0 1504 ENDIF 1505 1506 ! 1507 !-- If wall between gridpoint 5 and gridpoint 7, then 1508 !-- Neumann boundary condition has to be applied 1509 !-- (only one case as only building beneath is possible) 1510 IF ( gp_outside_of_building(5) == 1 .AND. & 1511 gp_outside_of_building(7) == 0 ) THEN 1512 num_gp = num_gp + 1 1513 location(num_gp,1) = (i+1) * dx 1514 location(num_gp,2) = j * dy 1515 location(num_gp,3) = k * dz 1516 ei(num_gp) = e(k+1,j,i+1) 1517 dissi(num_gp) = diss(k+1,j,i+1) 1518 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1519 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1520 de_dzi(num_gp) = 0.0 1521 ENDIF 1522 1523 ! 1524 !-- If wall between gridpoint 2 and gridpoint 4, then 1525 !-- Neumann boundary condition has to be applied 1526 !-- (only one case as only building beneath is possible) 1527 IF ( gp_outside_of_building(2) == 1 .AND. & 1528 gp_outside_of_building(4) == 0 ) THEN 1529 num_gp = num_gp + 1 1530 location(num_gp,1) = i * dx 1531 location(num_gp,2) = (j+1) * dy 1532 location(num_gp,3) = k * dz 1533 ei(num_gp) = e(k+1,j+1,i) 1534 dissi(num_gp) = diss(k+1,j+1,i) 1535 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1536 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1537 de_dzi(num_gp) = 0.0 1538 ENDIF 1539 1540 ! 1541 !-- If wall between gridpoint 6 and gridpoint 8, then 1542 !-- Neumann boundary condition has to be applied 1543 !-- (only one case as only building beneath is possible) 1544 IF ( gp_outside_of_building(6) == 1 .AND. & 1545 gp_outside_of_building(8) == 0 ) THEN 1546 num_gp = num_gp + 1 1547 location(num_gp,1) = (i+1) * dx 1548 location(num_gp,2) = (j+1) * dy 1549 location(num_gp,3) = k * dz 1550 ei(num_gp) = e(k+1,j+1,i+1) 1551 dissi(num_gp) = diss(k+1,j+1,i+1) 1552 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1553 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 1554 de_dzi(num_gp) = 0.0 1555 ENDIF 1556 1557 ! 1558 !-- Carry out the interpolation 1559 IF ( num_gp == 1 ) THEN 1560 ! 1561 !-- If only one of the gridpoints is situated outside of the 1562 !-- building, it follows that the values at the particle 1563 !-- location are the same as the gridpoint values 1564 e_int = ei(num_gp) 1565 1566 ELSE IF ( num_gp > 1 ) THEN 1567 1568 d_sum = 0.0 1569 ! 1570 !-- Evaluation of the distances between the gridpoints 1571 !-- contributing to the interpolated values, and the particle 1572 !-- location 1573 DO agp = 1, num_gp 1574 d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2 & 1575 + ( particles(n)%y-location(agp,2) )**2 & 1576 + ( particles(n)%z-location(agp,3) )**2 1577 d_sum = d_sum + d_gp_pl(agp) 1578 ENDDO 1579 1580 ! 1581 !-- Finally the interpolation can be carried out 1582 e_int = 0.0 1583 diss_int = 0.0 1584 de_dx_int = 0.0 1585 de_dy_int = 0.0 1586 de_dz_int = 0.0 1587 DO agp = 1, num_gp 1588 e_int = e_int + ( d_sum - d_gp_pl(agp) ) * & 1589 ei(agp) / ( (num_gp-1) * d_sum ) 1590 diss_int = diss_int + ( d_sum - d_gp_pl(agp) ) * & 1591 dissi(agp) / ( (num_gp-1) * d_sum ) 1592 de_dx_int = de_dx_int + ( d_sum - d_gp_pl(agp) ) * & 1593 de_dxi(agp) / ( (num_gp-1) * d_sum ) 1594 de_dy_int = de_dy_int + ( d_sum - d_gp_pl(agp) ) * & 1595 de_dyi(agp) / ( (num_gp-1) * d_sum ) 1596 de_dz_int = de_dz_int + ( d_sum - d_gp_pl(agp) ) * & 1597 de_dzi(agp) / ( (num_gp-1) * d_sum ) 1598 ENDDO 1599 1600 ENDIF 1601 1602 ENDIF 1603 1604 ENDIF 1072 1605 1073 1606 ! … … 1148 1681 1149 1682 ELSE 1683 1684 ! 1685 !-- Restriction of the size of the new timestep: compared to the 1686 !-- previous timestep the increase must not exceed 200% 1687 1688 dt_particle_old = particles(n)%age - particles(n)%age_m 1689 IF ( dt_particle > 2.0 * dt_particle_old ) THEN 1690 dt_particle = 2.0 * dt_particle_old 1691 ENDIF 1692 1150 1693 ! 1151 1694 !-- For old particles the SGS components are correlated with the 1152 1695 !-- values from the previous timestep. Random numbers have also to 1153 1696 !-- be limited (see above). 1697 !-- As negative values for the subgrid TKE are not allowed, the 1698 !-- change of the subgrid TKE with time cannot be smaller than 1699 !-- -e_int/dt_particle. This value is used as a lower boundary 1700 !-- value for the change of TKE 1701 1702 de_dt_min = - e_int / dt_particle 1703 1704 de_dt = ( e_int - particles(n)%e_m ) / dt_particle_old 1705 1706 IF ( de_dt < de_dt_min ) THEN 1707 de_dt = de_dt_min 1708 ENDIF 1709 1154 1710 particles(n)%speed_x_sgs = particles(n)%speed_x_sgs - & 1155 1711 fs_int * c_0 * diss_int * particles(n)%speed_x_sgs * & 1156 1712 dt_particle / ( 4.0 * sgs_wfu_part * e_int ) + & 1157 ( ( 2.0 * sgs_wfu_part * ( e_int - particles(n)%e_m ) /&1158 dt_particle ) * particles(n)%speed_x_sgs /&1713 ( 2.0 * sgs_wfu_part * de_dt * & 1714 particles(n)%speed_x_sgs / & 1159 1715 ( 2.0 * sgs_wfu_part * e_int ) + de_dx_int & 1160 1716 ) * dt_particle / 2.0 + & … … 1166 1722 fs_int * c_0 * diss_int * particles(n)%speed_y_sgs * & 1167 1723 dt_particle / ( 4.0 * sgs_wfv_part * e_int ) + & 1168 ( ( 2.0 * sgs_wfv_part * ( e_int - particles(n)%e_m ) /&1169 dt_particle ) * particles(n)%speed_y_sgs /&1724 ( 2.0 * sgs_wfv_part * de_dt * & 1725 particles(n)%speed_y_sgs / & 1170 1726 ( 2.0 * sgs_wfv_part * e_int ) + de_dy_int & 1171 1727 ) * dt_particle / 2.0 + & … … 1177 1733 fs_int * c_0 * diss_int * particles(n)%speed_z_sgs * & 1178 1734 dt_particle / ( 4.0 * sgs_wfw_part * e_int ) + & 1179 ( ( 2.0 * sgs_wfw_part * ( e_int - particles(n)%e_m ) /&1180 dt_particle ) * particles(n)%speed_z_sgs /&1735 ( 2.0 * sgs_wfw_part * de_dt * & 1736 particles(n)%speed_z_sgs / & 1181 1737 ( 2.0 * sgs_wfw_part * e_int ) + de_dz_int & 1182 1738 ) * dt_particle / 2.0 + & … … 1203 1759 1204 1760 ENDIF 1761 1762 ! 1763 !-- Remember the old age of the particle ( needed to prevent that a 1764 !-- particle crosses several PEs during one timestep and for the 1765 !-- evaluation of the subgrid particle velocity fluctuations ) 1766 particles(n)%age_m = particles(n)%age 1767 1205 1768 1206 1769 ! … … 1270 1833 1271 1834 ENDDO ! advection loop 1835 1836 1837 ! 1838 !-- Particle reflection from walls 1839 CALL cpu_log( log_point_s(48), 'advec_part_refle', 'start' ) 1840 1841 DO n = 1, number_of_particles 1842 1843 !------------------------------------------------------- 1844 i2 = (particles(n)%x + 0.5*dx) * ddx 1845 j2 = (particles(n)%y + 0.5*dy) * ddy 1846 k2 = particles(n)%z / dz + 1 1847 1848 particles_x = particles(n)%x 1849 particles_y = particles(n)%y 1850 particles_z = particles(n)%z 1851 1852 1853 IF (k2<=nzb_s_inner(j2,i2).and.nzb_s_inner(j2,i2)/=0) THEN 1854 1855 old_positions_x = particles(n)%x - particles(n)%speed_x * dt_particle 1856 old_positions_y = particles(n)%y - particles(n)%speed_y * dt_particle 1857 old_positions_z = particles(n)%z - particles(n)%speed_z * dt_particle 1858 i1= (old_positions_x +0.5*dx) * ddx 1859 j1 =(old_positions_y +0.5*dy) * ddy 1860 k1 = old_positions_z / dz 1861 1862 1863 !-- CASE 1 1864 1865 IF ( particles(n)%x > old_positions_x .and. & 1866 particles(n)%y > old_positions_y )THEN 1867 1868 1869 1870 1871 t_index = 1 1872 1873 Do i = i1, i2 1874 xline = i*dx+0.5*dx 1875 1876 t(t_index) =(xline-old_positions_x)/(particles(n)%x-old_positions_x) 1877 t_index = t_index + 1 1878 ENDDO 1879 1880 Do j = j1, j2 1881 yline = j*dy+0.5*dy 1882 1883 t(t_index) =(yline-old_positions_y)/(particles(n)%y-old_positions_y) 1884 t_index = t_index + 1 1885 ENDDO 1886 1887 IF ( particles(n)%z < old_positions_z ) THEN 1888 Do k = k1, k2 , -1 1889 zline = k*dz 1890 t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z) 1891 t_index = t_index + 1 1892 ENDDO 1893 ENDIF 1894 t_index_number = t_index-1 1895 !-- sorting t 1896 inc = 1 1897 jr = 1 1898 DO WHILE ( inc <= t_index_number ) 1899 inc = 3 * inc + 1 1900 ENDDO 1901 1902 DO WHILE ( inc > 1 ) 1903 inc = inc / 3 1904 DO ir = inc+1, t_index_number 1905 tmp_t = t(ir) 1906 jr = ir 1907 1908 DO WHILE ( t(jr-inc) > tmp_t ) 1909 t(jr) = t(jr-inc) 1910 jr = jr - inc 1911 IF ( jr <= inc ) THEN 1912 EXIT 1913 ENDIF 1914 ENDDO 1915 t(jr) = tmp_t 1916 1917 ENDDO 1918 ENDDO 1919 1920 1921 1922 1923 Do t_index = 1, t_index_number 1924 1925 positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x) 1926 positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y) 1927 positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z) 1928 1929 1930 1931 i3 = (positions_x + 0.5*dx) * ddx 1932 j3 = (positions_y + 0.5*dy) * ddy 1933 k3 = positions_z / dz 1934 1935 i5 = positions_x * ddx 1936 j5 = positions_y * ddy 1937 k5 = positions_z / dz 1938 1939 1940 1941 1942 1943 IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN 1944 1945 IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. & 1946 k3 == nzb_s_inner(j5,i3) ) THEN 1947 1948 particles(n)%z = 2*positions_z - particles_z 1949 particles(n)%speed_z = - particles(n)%speed_z 1950 1951 IF ( use_sgs_for_particles ) THEN 1952 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 1953 ENDIF 1954 GOTO 999 1955 ENDIF 1956 ENDIF 1957 1958 IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN 1959 1960 IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. & 1961 k3 == nzb_s_inner(j3,i5) ) THEN 1962 1963 particles(n)%z = 2*positions_z - particles_z 1964 particles(n)%speed_z = - particles(n)%speed_z 1965 1966 IF ( use_sgs_for_particles ) THEN 1967 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 1968 ENDIF 1969 1970 GOTO 999 1971 ENDIF 1972 ENDIF 1973 1974 IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN 1975 1976 IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. & 1977 k3 == nzb_s_inner(j3,i3) ) THEN 1978 1979 particles(n)%z = 2*positions_z - particles_z 1980 particles(n)%speed_z = - particles(n)%speed_z 1981 1982 IF ( use_sgs_for_particles ) THEN 1983 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 1984 ENDIF 1985 1986 GOTO 999 1987 ENDIF 1988 1989 IF ( positions_y == j3*dy-0.5*dy .and. & 1990 positions_z < nzb_s_inner(j3,i3)*dz) THEN 1991 1992 particles(n)%y = 2*positions_y - particles_y 1993 particles(n)%speed_y = - particles(n)%speed_y 1994 1995 IF ( use_sgs_for_particles ) THEN 1996 particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs 1997 ENDIF 1998 1999 GOTO 999 2000 ENDIF 2001 2002 2003 IF ( positions_x == i3*dx-0.5*dx.and. & 2004 positions_z < nzb_s_inner(j3,i3)*dz ) THEN 2005 2006 particles(n)%x = 2*positions_x - particles_x 2007 particles(n)%speed_x = - particles(n)%speed_x 2008 2009 IF ( use_sgs_for_particles ) THEN 2010 particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs 2011 GOTO 999 2012 ENDIF 2013 2014 ENDIF 2015 2016 ENDIF 2017 ENDDO 2018 ENDIF 2019 2020 !-- CASE 2 2021 2022 IF ( particles(n)%x > old_positions_x .and. & 2023 particles(n)%y < old_positions_y )THEN 2024 2025 2026 t_index = 1 2027 2028 Do i = i1, i2 2029 xline = i*dx+0.5*dx 2030 2031 t(t_index) =(xline-old_positions_x)/(particles(n)%x-old_positions_x) 2032 t_index = t_index + 1 2033 ENDDO 2034 2035 Do j = j1, j2,-1 2036 yline = j*dy-0.5*dy 2037 2038 t(t_index) =(old_positions_y-yline)/(old_positions_y-particles(n)%y) 2039 t_index = t_index + 1 2040 ENDDO 2041 2042 IF ( particles(n)%z < old_positions_z ) THEN 2043 Do k = k1, k2 , -1 2044 zline = k*dz 2045 t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z) 2046 t_index = t_index + 1 2047 ENDDO 2048 ENDIF 2049 t_index_number = t_index-1 2050 2051 !-- sorting t 2052 inc = 1 2053 jr = 1 2054 DO WHILE ( inc <= t_index_number ) 2055 inc = 3 * inc + 1 2056 ENDDO 2057 2058 DO WHILE ( inc > 1 ) 2059 inc = inc / 3 2060 DO ir = inc+1, t_index_number 2061 tmp_t = t(ir) 2062 jr = ir 2063 2064 DO WHILE ( t(jr-inc) > tmp_t ) 2065 t(jr) = t(jr-inc) 2066 jr = jr - inc 2067 IF ( jr <= inc ) THEN 2068 EXIT 2069 ENDIF 2070 ENDDO 2071 t(jr) = tmp_t 2072 2073 ENDDO 2074 ENDDO 2075 2076 2077 2078 2079 2080 2081 Do t_index = 1, t_index_number 2082 2083 positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x) 2084 positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y) 2085 positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z) 2086 2087 2088 2089 i3 = (positions_x + 0.5*dx) * ddx 2090 j3 = (positions_y + 0.5*dy) * ddy 2091 k3 = positions_z / dz 2092 2093 i5 = positions_x * ddx 2094 j5 = positions_y * ddy 2095 k5 = positions_z / dz 2096 2097 2098 2099 IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN 2100 2101 IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. & 2102 k3 == nzb_s_inner(j3,i5) ) THEN 2103 2104 particles(n)%z = 2*positions_z - particles_z 2105 particles(n)%speed_z = - particles(n)%speed_z 2106 2107 IF ( use_sgs_for_particles ) THEN 2108 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2109 ENDIF 2110 2111 GOTO 999 2112 ENDIF 2113 ENDIF 2114 2115 IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN 2116 2117 IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. & 2118 k3 == nzb_s_inner(j3,i3) ) THEN 2119 2120 particles(n)%z = 2*positions_z - particles_z 2121 particles(n)%speed_z = - particles(n)%speed_z 2122 2123 IF ( use_sgs_for_particles ) THEN 2124 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2125 ENDIF 2126 2127 GOTO 999 2128 ENDIF 2129 2130 2131 IF ( positions_x == i3*dx-0.5*dx.and.& 2132 positions_z < nzb_s_inner(j3,i3)*dz ) THEN 2133 2134 particles(n)%x = 2*positions_x - particles_x 2135 particles(n)%speed_x = - particles(n)%speed_x 2136 2137 IF ( use_sgs_for_particles ) THEN 2138 particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs 2139 ENDIF 2140 2141 GOTO 999 2142 ENDIF 2143 2144 ENDIF 2145 2146 2147 IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN 2148 2149 IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. & 2150 k3 == nzb_s_inner(j5,i3) ) THEN 2151 2152 particles(n)%z = 2*positions_z - particles_z 2153 particles(n)%speed_z = - particles(n)%speed_z 2154 IF ( use_sgs_for_particles ) THEN 2155 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2156 ENDIF 2157 2158 GOTO 999 2159 ENDIF 2160 2161 IF ( positions_y == j5 * dy + 0.5*dy.and.& 2162 positions_z < nzb_s_inner(j5,i3)*dz ) THEN 2163 2164 particles(n)%y = 2*positions_y - particles_y 2165 particles(n)%speed_y = - particles(n)%speed_y 2166 2167 IF ( use_sgs_for_particles ) THEN 2168 particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs 2169 ENDIF 2170 GOTO 999 2171 ENDIF 2172 2173 ENDIF 2174 ENDDO 2175 ENDIF 2176 2177 2178 !-- CASE 3 2179 2180 2181 2182 IF ( particles(n)%x < old_positions_x .and. & 2183 particles(n)%y > old_positions_y )THEN 2184 2185 2186 t_index = 1 2187 2188 Do i = i1, i2,-1 2189 xline = i*dx-0.5*dx 2190 2191 t(t_index) =(old_positions_x-xline)/(old_positions_x-particles(n)%x) 2192 t_index = t_index + 1 2193 ENDDO 2194 2195 Do j = j1, j2 2196 yline = j*dy+0.5*dy 2197 2198 t(t_index) =(yline-old_positions_y)/(particles(n)%y-old_positions_y) 2199 t_index = t_index + 1 2200 ENDDO 2201 2202 IF ( particles(n)%z < old_positions_z ) THEN 2203 Do k = k1, k2 , -1 2204 zline = k*dz 2205 t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z) 2206 t_index = t_index + 1 2207 ENDDO 2208 ENDIF 2209 2210 t_index_number = t_index-1 2211 2212 !-- sorting t 2213 inc = 1 2214 jr = 1 2215 DO WHILE ( inc <= t_index_number ) 2216 inc = 3 * inc + 1 2217 ENDDO 2218 2219 DO WHILE ( inc > 1 ) 2220 inc = inc / 3 2221 DO ir = inc+1, t_index_number 2222 tmp_t = t(ir) 2223 jr = ir 2224 2225 DO WHILE ( t(jr-inc) > tmp_t ) 2226 t(jr) = t(jr-inc) 2227 jr = jr - inc 2228 IF ( jr <= inc ) THEN 2229 EXIT 2230 ENDIF 2231 ENDDO 2232 t(jr) = tmp_t 2233 2234 ENDDO 2235 ENDDO 2236 2237 2238 2239 Do t_index = 1, t_index_number 2240 2241 positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x) 2242 positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y) 2243 positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z) 2244 2245 2246 2247 i3 = (positions_x + 0.5*dx) * ddx 2248 j3 = (positions_y + 0.5*dy) * ddy 2249 k3 = positions_z / dz 2250 2251 i5 = positions_x * ddx 2252 j5 = positions_y * ddy 2253 k5 = positions_z / dz 2254 2255 2256 2257 IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN 2258 2259 IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. & 2260 k3 == nzb_s_inner(j5,i3) ) THEN 2261 2262 particles(n)%z = 2*positions_z - particles_z 2263 particles(n)%speed_z = - particles(n)%speed_z 2264 2265 IF ( use_sgs_for_particles ) THEN 2266 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2267 ENDIF 2268 2269 GOTO 999 2270 ENDIF 2271 ENDIF 2272 2273 2274 2275 IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN 2276 2277 2278 IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. & 2279 k3 == nzb_s_inner(j3,i3) ) THEN 2280 2281 particles(n)%z = 2*positions_z - particles_z 2282 particles(n)%speed_z = - particles(n)%speed_z 2283 2284 IF ( use_sgs_for_particles ) THEN 2285 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2286 ENDIF 2287 2288 GOTO 999 2289 ENDIF 2290 2291 IF ( positions_y == j3*dy-0.5*dy .and. & 2292 positions_z < nzb_s_inner(j3,i3)*dz) THEN 2293 2294 particles(n)%y = 2*positions_y - particles_y 2295 particles(n)%speed_y = - particles(n)%speed_y 2296 2297 IF ( use_sgs_for_particles ) THEN 2298 particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs 2299 ENDIF 2300 2301 GOTO 999 2302 ENDIF 2303 2304 ENDIF 2305 2306 2307 IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN 2308 2309 IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. & 2310 k3 == nzb_s_inner(j3,i5) ) THEN 2311 2312 particles(n)%z = 2*positions_z - particles_z 2313 particles(n)%speed_z = - particles(n)%speed_z 2314 2315 IF ( use_sgs_for_particles ) THEN 2316 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2317 ENDIF 2318 2319 2320 GOTO 999 2321 ENDIF 2322 2323 IF ( positions_x == i5 * dx + 0.5*dx .and. & 2324 positions_z < nzb_s_inner(j3,i5)*dz ) THEN 2325 2326 particles(n)%x= 2*positions_x - particles_x 2327 particles(n)%speed_x = - particles(n)%speed_x 2328 2329 IF ( use_sgs_for_particles ) THEN 2330 particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs 2331 ENDIF 2332 2333 GOTO 999 2334 ENDIF 2335 2336 ENDIF 2337 ENDDO 2338 ENDIF 2339 2340 2341 !-- CASE 4 2342 2343 2344 IF ( particles(n)%x < old_positions_x .and. & 2345 particles(n)%y < old_positions_y )THEN 2346 2347 2348 t_index = 1 2349 2350 Do i = i1, i2,-1 2351 xline = i*dx-0.5*dx 2352 2353 t(t_index) =(old_positions_x-xline)/(old_positions_x-particles(n)%x) 2354 t_index = t_index + 1 2355 2356 2357 2358 ENDDO 2359 2360 Do j = j1, j2,-1 2361 yline = j*dy-0.5*dy 2362 2363 t(t_index) =(old_positions_y-yline)/(old_positions_y-particles(n)%y) 2364 t_index = t_index + 1 2365 2366 2367 ENDDO 2368 2369 IF ( particles(n)%z < old_positions_z ) THEN 2370 Do k = k1, k2 , -1 2371 zline = k*dz 2372 t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z) 2373 t_index = t_index + 1 2374 2375 2376 ENDDO 2377 ENDIF 2378 2379 t_index_number = t_index-1 2380 2381 !-- sorting t 2382 inc = 1 2383 jr = 1 2384 DO WHILE ( inc <= t_index_number ) 2385 inc = 3 * inc + 1 2386 ENDDO 2387 2388 DO WHILE ( inc > 1 ) 2389 inc = inc / 3 2390 DO ir = inc+1, t_index_number 2391 tmp_t = t(ir) 2392 jr = ir 2393 2394 DO WHILE ( t(jr-inc) > tmp_t ) 2395 t(jr) = t(jr-inc) 2396 jr = jr - inc 2397 IF ( jr <= inc ) THEN 2398 EXIT 2399 ENDIF 2400 ENDDO 2401 t(jr) = tmp_t 2402 2403 ENDDO 2404 ENDDO 2405 2406 2407 Do t_index = 1, t_index_number 2408 2409 positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x) 2410 positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y) 2411 positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z) 2412 2413 2414 2415 i3 = (positions_x + 0.5*dx) * ddx 2416 j3 = (positions_y + 0.5*dy) * ddy 2417 k3 = positions_z / dz 2418 2419 i5 = positions_x * ddx 2420 j5 = positions_y * ddy 2421 k5 = positions_z / dz 2422 2423 2424 2425 2426 IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN 2427 2428 IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. & 2429 k3 == nzb_s_inner(j3,i3) ) THEN 2430 2431 particles(n)%z = 2*positions_z - particles_z 2432 particles(n)%speed_z = - particles(n)%speed_z 2433 2434 IF ( use_sgs_for_particles ) THEN 2435 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2436 ENDIF 2437 2438 2439 GOTO 999 2440 ENDIF 2441 ENDIF 2442 2443 2444 IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN 2445 2446 IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. & 2447 k3 == nzb_s_inner(j3,i5) ) THEN 2448 2449 particles(n)%z = 2*positions_z - particles_z 2450 particles(n)%speed_z = - particles(n)%speed_z 2451 2452 IF ( use_sgs_for_particles ) THEN 2453 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2454 ENDIF 2455 2456 GOTO 999 2457 ENDIF 2458 2459 IF ( positions_x == i5 * dx + 0.5*dx .and. nzb_s_inner(j3,i5)/=0.and.& 2460 positions_z < nzb_s_inner(j3,i5)*dz) THEN 2461 2462 particles(n)%x= 2*positions_x - particles_x 2463 particles(n)%speed_x = - particles(n)%speed_x 2464 2465 IF ( use_sgs_for_particles ) THEN 2466 particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs 2467 ENDIF 2468 GOTO 999 2469 ENDIF 2470 2471 ENDIF 2472 2473 IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0) THEN 2474 2475 2476 IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. & 2477 k5 == nzb_s_inner(j5,i3) ) THEN 2478 2479 particles(n)%z = 2*positions_z - particles_z 2480 particles(n)%speed_z = - particles(n)%speed_z 2481 2482 IF ( use_sgs_for_particles ) THEN 2483 particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs 2484 ENDIF 2485 2486 GOTO 999 2487 ENDIF 2488 2489 IF ( positions_y == j5 * dy + 0.5*dy .and. nzb_s_inner(j5,i3)/=0.and.& 2490 positions_z < nzb_s_inner(j5,i3)*dz ) THEN 2491 2492 particles(n)%y= 2*positions_y - particles_y 2493 particles(n)%speed_y = - particles(n)%speed_y 2494 2495 IF ( use_sgs_for_particles ) THEN 2496 particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs 2497 ENDIF 2498 2499 GOTO 999 2500 ENDIF 2501 2502 2503 ENDIF 2504 ENDDO 2505 ENDIF 2506 2507 999 CONTINUE 2508 2509 ENDIF 2510 2511 !--------------------------------------------------- 2512 2513 ENDDO ! particle reflection from walls 2514 2515 CALL cpu_log( log_point_s(48), 'advec_part_refle', 'stop' ) 2516 2517 ! 2518 !-- User-defined actions after the evaluation of the new particle position 2519 CALL user_advec_particles 1272 2520 1273 2521 ! … … 1500 2748 trlp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1501 2749 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1502 0 , 0, 0, 0 )2750 0.0, 0, 0, 0, 0 ) 1503 2751 trrp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1504 2752 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1505 0 , 0, 0, 0 )2753 0.0, 0, 0, 0, 0 ) 1506 2754 1507 2755 IF ( use_particle_tails ) THEN … … 1933 3181 trsp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1934 3182 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1935 0 , 0, 0, 0 )3183 0.0, 0, 0, 0, 0 ) 1936 3184 trnp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1937 3185 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 1938 0 , 0, 0, 0 )3186 0.0, 0, 0, 0, 0 ) 1939 3187 1940 3188 IF ( use_particle_tails ) THEN … … 2460 3708 2461 3709 nn = particles(n)%tail_id 3710 3711 ! 3712 !-- Stop if particles have moved further than the length of one 3713 !-- PE subdomain 3714 IF ( ABS(particles(n)%speed_x) > & 3715 ((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m) .OR. & 3716 ABS(particles(n)%speed_y) > & 3717 ((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) ) THEN 3718 3719 PRINT*, '+++ advec_particles: particle too fast. n = ', n 3720 #if defined( __parallel ) 3721 CALL MPI_ABORT( comm2d, 9999, ierr ) 3722 #else 3723 CALL local_stop 3724 #endif 3725 ENDIF 2462 3726 2463 3727 IF ( particles(n)%age > particle_maximum_age .AND. & -
palm/trunk/SOURCE/buoyancy.f90
r4 r57 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Reference temperature pt_reference can be used. 7 7 ! 8 8 ! Former revisions: … … 59 59 ! 60 60 !-- Normal case: horizontal surface 61 DO i = nxl, nxr 62 DO j = nys, nyn 63 DO k = nzb_s_inner(j,i)+1, nzt-1 64 tend(k,j,i) = tend(k,j,i) + g * 0.5 * ( & 61 IF ( use_pt_reference ) THEN 62 DO i = nxl, nxr 63 DO j = nys, nyn 64 DO k = nzb_s_inner(j,i)+1, nzt-1 65 tend(k,j,i) = tend(k,j,i) + g * 0.5 * ( & 66 ( theta(k,j,i) - hom(k,1,pr,0) ) / pt_reference + & 67 ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / pt_reference & 68 ) 69 ENDDO 70 ENDDO 71 ENDDO 72 ELSE 73 DO i = nxl, nxr 74 DO j = nys, nyn 75 DO k = nzb_s_inner(j,i)+1, nzt-1 76 tend(k,j,i) = tend(k,j,i) + g * 0.5 * ( & 65 77 ( theta(k,j,i) - hom(k,1,pr,0) ) / hom(k,1,pr,0) + & 66 78 ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) & 67 ) 68 ENDDO 69 ENDDO 70 ENDDO 79 ) 80 ENDDO 81 ENDDO 82 ENDDO 83 ENDIF 71 84 72 85 ELSE … … 136 149 ! 137 150 !-- Normal case: horizontal surface 138 DO k = nzb_s_inner(j,i)+1, nzt-1 139 tend(k,j,i) = tend(k,j,i) + g * 0.5 * ( & 151 IF ( use_pt_reference ) THEN 152 DO k = nzb_s_inner(j,i)+1, nzt-1 153 tend(k,j,i) = tend(k,j,i) + g * 0.5 * ( & 154 ( theta(k,j,i) - hom(k,1,pr,0) ) / pt_reference + & 155 ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / pt_reference & 156 ) 157 ENDDO 158 ELSE 159 DO k = nzb_s_inner(j,i)+1, nzt-1 160 tend(k,j,i) = tend(k,j,i) + g * 0.5 * ( & 140 161 ( theta(k,j,i) - hom(k,1,pr,0) ) / hom(k,1,pr,0) + & 141 162 ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) & 142 ) 143 ENDDO 163 ) 164 ENDDO 165 ENDIF 144 166 145 167 ELSE -
palm/trunk/SOURCE/check_parameters.f90
r51 r57 5 5 ! ----------------- 6 6 ! "by_user" allowed as initializing action, -data_output_ts, 7 ! leapfrog with non-flat topography not allowed any more 7 ! leapfrog with non-flat topography not allowed any more, pt_reference is 8 ! checked 8 9 ! 9 10 ! Former revisions: … … 574 575 575 576 ! 577 !-- Reference temperature to be used in buoyancy terms 578 IF ( pt_reference /= 9999999.9 ) use_pt_reference = .TRUE. 579 580 ! 576 581 !-- In case of a given slope, compute the relevant quantities 577 582 IF ( alpha_surface /= 0.0 ) THEN -
palm/trunk/SOURCE/diffusion_e.f90
r39 r57 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Reference temperature pt_reference can be used in buoyancy term 7 7 ! 8 8 ! Former revisions: … … 82 82 dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k) 83 83 IF ( dpt_dz > 0.0 ) THEN 84 l_stable = 0.76 * SQRT( e(k,j,i) ) / & 85 SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5 84 IF ( use_pt_reference ) THEN 85 l_stable = 0.76 * SQRT( e(k,j,i) ) / & 86 SQRT( g / pt_reference * dpt_dz ) + 1E-5 87 ELSE 88 l_stable = 0.76 * SQRT( e(k,j,i) ) / & 89 SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5 90 ENDIF 86 91 ELSE 87 92 l_stable = l_grid(k) … … 195 200 dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k) 196 201 IF ( dpt_dz > 0.0 ) THEN 197 l_stable = 0.76 * SQRT( e(k,j,i) ) / & 198 SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5 202 IF ( use_pt_reference ) THEN 203 l_stable = 0.76 * SQRT( e(k,j,i) ) / & 204 SQRT( g / pt_reference * dpt_dz ) + 1E-5 205 ELSE 206 l_stable = 0.76 * SQRT( e(k,j,i) ) / & 207 SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5 208 ENDIF 199 209 ELSE 200 210 l_stable = l_grid(k) -
palm/trunk/SOURCE/diffusion_u.f90
r56 r57 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes 6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes, 7 ! z0 removed from argument list 7 8 ! 8 9 ! Former revisions: … … 50 51 ! Call for all grid points 51 52 !------------------------------------------------------------------------------! 52 SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w , z0)53 SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w ) 53 54 54 55 USE control_parameters … … 61 62 REAL :: kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp 62 63 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1) 63 REAL :: z0(nys-1:nyn+1,nxl-1:nxr+1)64 64 REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) 65 65 REAL, DIMENSION(:,:), POINTER :: usws … … 207 207 !------------------------------------------------------------------------------! 208 208 SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, & 209 v, w , z0)209 v, w ) 210 210 211 211 USE control_parameters … … 218 218 REAL :: kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp 219 219 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1) 220 REAL :: z0(nys-1:nyn+1,nxl-1:nxr+1)221 220 REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) 222 221 REAL, DIMENSION(nzb:nzt+1) :: usvs -
palm/trunk/SOURCE/diffusion_v.f90
r56 r57 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes 6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes, 7 ! z0 removed from argument list 7 8 ! 8 9 ! Former revisions: … … 48 49 ! Call for all grid points 49 50 !------------------------------------------------------------------------------! 50 SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w , z0)51 SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w ) 51 52 52 53 USE control_parameters … … 59 60 REAL :: kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp 60 61 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1) 61 REAL :: z0(nys-1:nyn+1,nxl-1:nxr+1)62 62 REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) 63 63 REAL, DIMENSION(:,:), POINTER :: vsws … … 205 205 !------------------------------------------------------------------------------! 206 206 SUBROUTINE diffusion_v_ij( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, & 207 vsws, w , z0)207 vsws, w ) 208 208 209 209 USE control_parameters … … 216 216 REAL :: kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp 217 217 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1) 218 REAL :: z0(nys-1:nyn+1,nxl-1:nxr+1)219 218 REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) 220 219 REAL, DIMENSION(nzb:nzt+1) :: vsus -
palm/trunk/SOURCE/diffusion_w.f90
r56 r57 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes 6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes, 7 ! z0 removed from argument list 7 8 ! 8 9 ! Former revisions: … … 46 47 !------------------------------------------------------------------------------! 47 48 SUBROUTINE diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, & 48 w , z0)49 w ) 49 50 50 51 USE control_parameters … … 59 60 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), & 60 61 km_damp_y(nys-1:nyn+1) 61 REAL :: z0(nys-1:nyn+1,nxl-1:nxr+1)62 62 REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) 63 63 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w … … 194 194 !------------------------------------------------------------------------------! 195 195 SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, & 196 tend, u, v, w , z0)196 tend, u, v, w ) 197 197 198 198 USE control_parameters … … 207 207 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), & 208 208 km_damp_y(nys-1:nyn+1) 209 REAL :: z0(nys-1:nyn+1,nxl-1:nxr+1)210 209 REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) 211 210 REAL, DIMENSION(nzb:nzt+1) :: wsus, wsvs -
palm/trunk/SOURCE/diffusivities.f90
r4 r57 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Reference temperature pt_reference can be used in buoyancy term 7 7 ! 8 8 ! Former revisions: … … 90 90 dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k) 91 91 IF ( dpt_dz > 0.0 ) THEN 92 l_stable = 0.76 * sqrt_e(k) / & 93 SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5 92 IF ( use_pt_reference ) THEN 93 l_stable = 0.76 * sqrt_e(k) / & 94 SQRT( g / pt_reference * dpt_dz ) + 1E-5 95 ELSE 96 l_stable = 0.76 * sqrt_e(k) / & 97 SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5 98 ENDIF 94 99 ELSE 95 100 l_stable = l_grid(k) -
palm/trunk/SOURCE/header.f90
r46 r57 861 861 !-- Other quantities 862 862 WRITE ( io, 411 ) g 863 IF ( use_pt_reference ) WRITE ( io, 412 ) pt_reference 863 864 864 865 ! 865 866 !-- Cloud physics parameters 866 867 IF ( cloud_physics ) THEN 867 WRITE ( io, 41 2)868 WRITE ( io, 41 3) surface_pressure, r_d, rho_surface, cp, l_v868 WRITE ( io, 415 ) 869 WRITE ( io, 416 ) surface_pressure, r_d, rho_surface, cp, l_v 869 870 ENDIF 870 871 … … 1326 1327 ' f* = ',F9.6,' 1/s') 1327 1328 411 FORMAT (/' Gravity : g = ',F4.1,' m/s**2') 1328 412 FORMAT (/' Cloud physics parameters:'/ & 1329 412 FORMAT (/' Reference temperature in buoyancy terms: ',F8.4,' K') 1330 415 FORMAT (/' Cloud physics parameters:'/ & 1329 1331 ' ------------------------'/) 1330 41 3FORMAT (' Surface pressure : p_0 = ',F7.2,' hPa'/ &1332 416 FORMAT (' Surface pressure : p_0 = ',F7.2,' hPa'/ & 1331 1333 ' Gas constant : R = ',F5.1,' J/(kg K)'/ & 1332 1334 ' Density of air : rho_0 = ',F5.3,' kg/m**3'/ & -
palm/trunk/SOURCE/init_particles.f90
r39 r57 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! displacements for mpi_particle_type changed, age_m initialized 7 7 ! 8 8 ! Former revisions: … … 61 61 blocklengths(1) = 18; blocklengths(2) = 4; blocklengths(3) = 1 62 62 #if defined( __t3eb ) 63 displacements(1) = 0; displacements(2) = 1 44; displacements(3) = 17663 displacements(1) = 0; displacements(2) = 152; displacements(3) = 184 64 64 #else 65 displacements(1) = 0; displacements(2) = 1 44; displacements(3) = 16065 displacements(1) = 0; displacements(2) = 152; displacements(3) = 168 66 66 #endif 67 67 types(1) = MPI_REAL … … 192 192 particles = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 193 193 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 194 0 , 0, 0, 0 )194 0.0, 0, 0, 0, 0 ) 195 195 particle_groups = particle_groups_type( 0.0, 0.0, 0.0, 0.0 ) 196 196 … … 267 267 particles(n)%z = pos_z 268 268 particles(n)%age = 0.0 269 particles(n)%age_m = 0.0 269 270 particles(n)%dt_sum = 0.0 270 271 particles(n)%dvrp_psize = dvrp_psize -
palm/trunk/SOURCE/modules.f90
r51 r57 7 7 ! +array rif_wall 8 8 ! +netcdf_64bit_3d, zu_s_inner, zw_w_inner, id_var_zusi_*, id_var_zwwi_*, 9 ! ts_value, u_nzb_p1_for_vfc, v_nzb_p1_for_vfc 9 ! ts_value, u_nzb_p1_for_vfc, v_nzb_p1_for_vfc, pt_reference, use_pt_reference, 10 ! +age_m in particle_type 10 11 ! -data_output_ts, dots_n 11 12 ! arrays dots_label and dots_unit now dimensioned with dots_max … … 296 297 sloping_surface = .FALSE., stop_dt = .FALSE., & 297 298 terminate_run = .FALSE., use_prior_plot1d_parameters = .FALSE.,& 298 use_ surface_fluxes = .FALSE., use_top_fluxes = .FALSE., &299 use_ ug_for_galilei_tr = .TRUE., &299 use_pt_reference = .FALSE., use_surface_fluxes = .FALSE., & 300 use_top_fluxes = .FALSE., use_ug_for_galilei_tr = .TRUE., & 300 301 use_upstream_for_tke = .FALSE., wall_adjustment = .TRUE. 301 302 … … 330 331 overshoot_limit_u = 0.0, overshoot_limit_v = 0.0, & 331 332 overshoot_limit_w = 0.0, particle_maximum_age = 9999999.9, & 332 phi = 55.0, prandtl_number = 1.0, pt_slope_offset = 0.0, & 333 pt_surface = 300.0, pt_surface_initial_change = 0.0, & 334 q_surface = 0.0, q_surface_initial_change = 0.0, & 335 rayleigh_damping_factor = -1.0, rayleigh_damping_height = -1.0, & 336 residual_limit = 1.0E-4, restart_time = 9999999.9, rif_max = 1.0, & 337 rif_min = -5.0, roughness_length = 0.1, simulated_time = 0.0, & 333 phi = 55.0, prandtl_number = 1.0, pt_reference = 9999999.9, & 334 pt_slope_offset = 0.0, pt_surface = 300.0, & 335 pt_surface_initial_change = 0.0, q_surface = 0.0, & 336 q_surface_initial_change = 0.0, rayleigh_damping_factor = -1.0, & 337 rayleigh_damping_height = -1.0, residual_limit = 1.0E-4, & 338 restart_time = 9999999.9, rif_max = 1.0, rif_min = -5.0, & 339 roughness_length = 0.1, simulated_time = 0.0, & 338 340 simulated_time_at_begin, sin_alpha_surface, & 339 341 skip_time_data_output = 0.0, skip_time_data_output_av = 9999999.9,& … … 828 830 TYPE particle_type 829 831 SEQUENCE 830 REAL :: age, dt_sum, dvrp_psize, e_m, origin_x, origin_y, origin_z, &831 radius, speed_x, speed_x_sgs, speed_y, speed_y_sgs, &832 speed_ z, speed_z_sgs, weight_factor, x, y, z832 REAL :: age, age_m, dt_sum, dvrp_psize, e_m, origin_x, origin_y, & 833 origin_z, radius, speed_x, speed_x_sgs, speed_y, & 834 speed_y_sgs, speed_z, speed_z_sgs, weight_factor, x, y, z 833 835 INTEGER :: color, group, tailpoints, tail_id 834 836 END TYPE particle_type -
palm/trunk/SOURCE/palm.f90
r4 r57 1 #if defined( __vtk )2 SUBROUTINE palm3 #else4 1 PROGRAM palm 5 #endif6 2 7 3 !------------------------------------------------------------------------------! 8 4 ! Actual revisions: 9 5 ! ----------------- 10 ! TEST: open(9)6 ! __vtk directives removed, open unit 9 for debug output 11 7 ! 12 8 ! Former revisions: … … 99 95 100 96 ! 97 !-- Open a file for debug output 98 OPEN( 9, FILE='DEBUG'//myid_char, FORM='FORMATTED' ) 99 100 ! 101 101 !-- Generate grid parameters 102 102 CALL init_grid … … 106 106 CALL check_parameters 107 107 108 OPEN( 9, FILE='DEBUG'//myid_char, FORM='FORMATTED' )109 108 ! 110 109 !-- Initialize all necessary variables … … 169 168 #endif 170 169 171 #if defined( __vtk )172 END SUBROUTINE palm173 #else174 170 END PROGRAM palm 175 #endif176 171 177 -
palm/trunk/SOURCE/parin.f90
r48 r57 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! +netcdf_64bit_3d in d3par, -data_output_ts6 ! +netcdf_64bit_3d in d3par, +pt_reference in inipar, -data_output_ts 7 7 ! 8 8 ! Former revisions: … … 65 65 overshoot_limit_pt, overshoot_limit_u, & 66 66 overshoot_limit_v, overshoot_limit_w, passive_scalar, & 67 phi, prandtl_layer, precipitation, pt_surface, & 68 pt_surface_initial_change, pt_vertical_gradient, & 69 pt_vertical_gradient_level, q_surface, & 70 q_surface_initial_change, q_vertical_gradient, & 71 q_vertical_gradient_level, radiation, random_generator, & 72 random_heatflux, rif_max, rif_min, roughness_length, & 73 scalar_advec, statistic_regions, surface_heatflux, & 74 surface_pressure, surface_scalarflux, surface_waterflux,& 75 s_surface, s_surface_initial_change, & 76 s_vertical_gradient, s_vertical_gradient_level, & 77 top_heatflux, timestep_scheme, topography, ug_surface, & 67 phi, prandtl_layer, precipitation, pt_reference, & 68 pt_surface, pt_surface_initial_change, & 69 pt_vertical_gradient, pt_vertical_gradient_level, & 70 q_surface, q_surface_initial_change, & 71 q_vertical_gradient, q_vertical_gradient_level, & 72 radiation, random_generator, random_heatflux, rif_max, & 73 rif_min, roughness_length, scalar_advec, & 74 statistic_regions, surface_heatflux, surface_pressure, & 75 surface_scalarflux, surface_waterflux, s_surface, & 76 s_surface_initial_change, s_vertical_gradient, & 77 s_vertical_gradient_level, top_heatflux, & 78 timestep_scheme, topography, ug_surface, & 78 79 ug_vertical_gradient, ug_vertical_gradient_level, & 79 80 ups_limit_e, ups_limit_pt, ups_limit_u, ups_limit_v, & -
palm/trunk/SOURCE/production_e.f90
r56 r57 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes_e 6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes_e, 7 ! reference temperature pt_reference can be used in buoyancy term 7 8 ! 8 9 ! Former revisions: … … 359 360 IF ( .NOT. moisture ) THEN 360 361 361 DO j = nys, nyn 362 363 DO k = nzb_diff_s_inner(j,i), nzt_diff 364 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * & 365 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 362 IF ( use_pt_reference ) THEN 363 364 DO j = nys, nyn 365 DO k = nzb_diff_s_inner(j,i), nzt_diff 366 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g/pt_reference * & 367 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 368 ENDDO 369 370 IF ( use_surface_fluxes ) THEN 371 k = nzb_diff_s_inner(j,i)-1 372 tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i) 373 ENDIF 374 375 IF ( use_top_fluxes ) THEN 376 k = nzt 377 tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i) 378 ENDIF 366 379 ENDDO 367 380 368 IF ( use_surface_fluxes ) THEN 369 k = nzb_diff_s_inner(j,i)-1 370 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i) 371 ENDIF 372 373 IF ( use_top_fluxes ) THEN 374 k = nzt 375 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i) 376 ENDIF 377 378 ENDDO 381 ELSE 382 383 DO j = nys, nyn 384 DO k = nzb_diff_s_inner(j,i), nzt_diff 385 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * & 386 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 387 ENDDO 388 389 IF ( use_surface_fluxes ) THEN 390 k = nzb_diff_s_inner(j,i)-1 391 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i) 392 ENDIF 393 394 IF ( use_top_fluxes ) THEN 395 k = nzt 396 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i) 397 ENDIF 398 ENDDO 399 400 ENDIF 379 401 380 402 ELSE … … 736 758 IF ( .NOT. moisture ) THEN 737 759 738 DO k = nzb_diff_s_inner(j,i), nzt_diff 739 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * & 740 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 741 ENDDO 742 743 IF ( use_surface_fluxes ) THEN 744 k = nzb_diff_s_inner(j,i)-1 745 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i) 746 ENDIF 747 748 IF ( use_top_fluxes ) THEN 749 k = nzt 750 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i) 751 ENDIF 760 IF ( use_pt_reference ) THEN 761 762 DO k = nzb_diff_s_inner(j,i), nzt_diff 763 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * & 764 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 765 ENDDO 766 767 IF ( use_surface_fluxes ) THEN 768 k = nzb_diff_s_inner(j,i)-1 769 tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i) 770 ENDIF 771 772 IF ( use_top_fluxes ) THEN 773 k = nzt 774 tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i) 775 ENDIF 776 777 ELSE 778 779 DO k = nzb_diff_s_inner(j,i), nzt_diff 780 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * & 781 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 782 ENDDO 783 784 IF ( use_surface_fluxes ) THEN 785 k = nzb_diff_s_inner(j,i)-1 786 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i) 787 ENDIF 788 789 IF ( use_top_fluxes ) THEN 790 k = nzt 791 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i) 792 ENDIF 793 794 ENDIF 752 795 753 796 ELSE -
palm/trunk/SOURCE/prognostic_equations.f90
r39 r57 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! z0 removed from arguments in calls of diffusion_u/v/w 7 7 ! 8 8 ! Former revisions: … … 130 130 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 131 131 CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, u_m, & 132 usws_m, v_m, w_m , z0)132 usws_m, v_m, w_m ) 133 133 ELSE 134 134 CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, & 135 v, w , z0)135 v, w ) 136 136 ENDIF 137 137 CALL coriolis( i, j, 1 ) … … 198 198 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 199 199 CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, u_m, & 200 v_m, vsws_m, w_m , z0)200 v_m, vsws_m, w_m ) 201 201 ELSE 202 202 CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, & 203 vsws, w , z0)203 vsws, w ) 204 204 ENDIF 205 205 CALL coriolis( i, j, 2 ) … … 265 265 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 266 266 CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, km_damp_y, & 267 tend, u_m, v_m, w_m , z0)267 tend, u_m, v_m, w_m ) 268 268 ELSE 269 269 CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, & 270 tend, u, v, w , z0)270 tend, u, v, w ) 271 271 ENDIF 272 272 CALL coriolis( i, j, 3 ) … … 668 668 THEN 669 669 CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, & 670 u_m, usws_m, v_m, w_m , z0)670 u_m, usws_m, v_m, w_m ) 671 671 ELSE 672 672 CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, & 673 usws, v, w , z0)673 usws, v, w ) 674 674 ENDIF 675 675 CALL coriolis( i, j, 1 ) … … 718 718 THEN 719 719 CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, & 720 u_m, v_m, vsws_m, w_m , z0)720 u_m, v_m, vsws_m, w_m ) 721 721 ELSE 722 722 CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, & 723 vsws, w , z0)723 vsws, w ) 724 724 ENDIF 725 725 CALL coriolis( i, j, 2 ) … … 767 767 THEN 768 768 CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, & 769 km_damp_y, tend, u_m, v_m, w_m , z0)769 km_damp_y, tend, u_m, v_m, w_m ) 770 770 ELSE 771 771 CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, & 772 tend, u, v, w , z0)772 tend, u, v, w ) 773 773 ENDIF 774 774 CALL coriolis( i, j, 3 ) … … 1041 1041 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 1042 1042 CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, v_m, & 1043 w_m , z0)1043 w_m ) 1044 1044 ELSE 1045 CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w , z0)1045 CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w ) 1046 1046 ENDIF 1047 1047 CALL coriolis( 1 ) … … 1113 1113 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 1114 1114 CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, & 1115 w_m , z0)1115 w_m ) 1116 1116 ELSE 1117 CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w , z0)1117 CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w ) 1118 1118 ENDIF 1119 1119 CALL coriolis( 2 ) … … 1184 1184 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 1185 1185 CALL diffusion_w( ddzu, ddzw, km_m, km_damp_x, km_damp_y, tend, u_m, & 1186 v_m, w_m , z0)1186 v_m, w_m ) 1187 1187 ELSE 1188 CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w, & 1189 z0 ) 1188 CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w ) 1190 1189 ENDIF 1191 1190 CALL coriolis( 3 ) -
palm/trunk/SOURCE/read_var_list.f90
r39 r57 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! +pt_reference 7 7 ! 8 8 ! Former revisions: … … 252 252 CASE ( 'pt_init' ) 253 253 READ ( 13 ) pt_init 254 CASE ( 'pt_reference' ) 255 READ ( 13 ) pt_reference 254 256 CASE ( 'pt_surface' ) 255 257 READ ( 13 ) pt_surface -
palm/trunk/SOURCE/user_interface.f90
r48 r57 8 8 ! routine user_statistics now has one argument (sr), 9 9 ! sample for generating time series quantities added 10 ! Bugfix in sample for reading user defined data from restart file (user_init) 10 11 ! 11 12 ! Former revisions: … … 174 175 ! 175 176 ! END SELECT 177 ! 178 ! READ ( 13 ) field_chr 179 ! 176 180 ! ENDDO 177 181 ! ENDIF -
palm/trunk/SOURCE/write_var_list.f90
r39 r57 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! +pt_refrence 7 7 ! 8 8 ! Former revisions: … … 218 218 WRITE ( 14 ) 'pt_init ' 219 219 WRITE ( 14 ) pt_init 220 WRITE ( 14 ) 'pt_reference ' 221 WRITE ( 14 ) pt_reference 220 222 WRITE ( 14 ) 'pt_surface ' 221 223 WRITE ( 14 ) pt_surface
Note: See TracChangeset
for help on using the changeset viewer.
