Changeset 19 for palm/trunk/SOURCE
- Timestamp:
- Feb 23, 2007 4:53:48 AM (18 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 23 edited
-
advec_s_pw.f90 (modified) (3 diffs)
-
boundary_conds.f90 (modified) (4 diffs)
-
calc_liquid_water_content.f90 (modified) (2 diffs)
-
calc_precipitation.f90 (modified) (3 diffs)
-
check_parameters.f90 (modified) (6 diffs)
-
diffusion_e.f90 (modified) (8 diffs)
-
diffusion_s.f90 (modified) (11 diffs)
-
flow_statistics.f90 (modified) (3 diffs)
-
header.f90 (modified) (4 diffs)
-
impact_of_latent_heat.f90 (modified) (3 diffs)
-
init_3d_model.f90 (modified) (8 diffs)
-
init_grid.f90 (modified) (3 diffs)
-
init_pt_anomaly.f90 (modified) (2 diffs)
-
modules.f90 (modified) (9 diffs)
-
parin.f90 (modified) (2 diffs)
-
production_e.f90 (modified) (12 diffs)
-
prognostic_equations.f90 (modified) (33 diffs)
-
read_3d_binary.f90 (modified) (3 diffs)
-
read_var_list.f90 (modified) (2 diffs)
-
spline_z.f90 (modified) (2 diffs)
-
swap_timelevel.f90 (modified) (3 diffs)
-
write_3d_binary.f90 (modified) (3 diffs)
-
write_var_list.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_s_pw.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Calculation extended for gridpoint nzt 7 7 ! 8 8 ! Former revisions: … … 58 58 DO i = nxl, nxr 59 59 DO j = nys, nyn 60 DO k = nzb_s_inner(j,i)+1, nzt -160 DO k = nzb_s_inner(j,i)+1, nzt 61 61 tend(k,j,i) = tend(k,j,i) & 62 62 -0.5 * ( ( u(k,j,i+1) - u_gtrans ) * ( sk(k,j,i+1) - sk(k,j,i) ) & … … 93 93 94 94 95 DO k = nzb_s_inner(j,i)+1, nzt -195 DO k = nzb_s_inner(j,i)+1, nzt 96 96 tend(k,j,i) = tend(k,j,i) & 97 97 -0.5 * ( ( u(k,j,i+1) - u_gtrans ) * ( sk(k,j,i+1) - sk(k,j,i) ) & -
palm/trunk/SOURCE/boundary_conds.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these 7 ! values are now calculated by the prognostic equation, 8 ! Dirichlet and zero gradient condition for pt established at top boundary 7 9 ! 8 10 ! Former revisions: … … 100 102 ! 101 103 !-- Temperature at top boundary 102 IF ( ibc_pt_t == 1 ) THEN 103 pt(nzt,:,:) = pt(nzt-1,:,:) + bc_pt_t_val * dzu(nzt) 104 IF ( ibc_pt_t == 0 ) THEN 105 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 106 pt(nzt+1,:,:) = pt_m(nzt+1,:,:) 107 ELSE 108 pt(nzt+1,:,:) = pt_p(nzt+1,:,:) ! pt_m not used for Runge-Kutta 109 ENDIF 110 ELSEIF ( ibc_pt_t == 1 ) THEN 111 pt(nzt+1,:,:) = pt(nzt,:,:) 112 ELSEIF ( ibc_pt_t == 2 ) THEN 104 113 pt(nzt+1,:,:) = pt(nzt,:,:) + bc_pt_t_val * dzu(nzt+1) 105 114 ENDIF … … 114 123 ENDDO 115 124 ENDDO 116 e(nzt,:,:) = e(nzt-1,:,:)117 125 e(nzt+1,:,:) = e(nzt,:,:) 118 126 ENDIF … … 147 155 ! 148 156 !-- Top boundary 149 q(nzt,:,:) = q(nzt-1,:,:) + bc_q_t_val * dzu(nzt)150 157 q(nzt+1,:,:) = q(nzt,:,:) + bc_q_t_val * dzu(nzt+1) 151 158 ENDIF -
palm/trunk/SOURCE/calc_liquid_water_content.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Old comment removed 7 7 ! 8 8 ! Former revisions: … … 79 79 ENDDO 80 80 81 !82 !-- Setting boundary values of ql83 ! CALL exchange_horiz( ql, 0, 0 )84 85 81 END SUBROUTINE calc_liquid_water_content -
palm/trunk/SOURCE/calc_precipitation.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Calculation extended for gridpoint nzt 7 7 ! 8 8 ! Former revisions: … … 54 54 DO i = nxl, nxr 55 55 DO j = nys, nyn 56 DO k = nzb_2d(j,i)+1, nzt -156 DO k = nzb_2d(j,i)+1, nzt 57 57 58 58 IF ( ql(k,j,i) > ql_crit ) THEN … … 87 87 88 88 89 DO k = nzb_2d(j,i)+1, nzt -189 DO k = nzb_2d(j,i)+1, nzt 90 90 91 91 IF ( ql(k,j,i) > ql_crit ) THEN -
palm/trunk/SOURCE/check_parameters.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Temperature and humidity gradients at top are now calculated for nzt+1, 7 ! top_heatflux and respective boundary condition bc_pt_t is checked 7 8 ! 8 9 ! Former revisions: … … 483 484 !-- Store temperature gradient at the top boundary for possile Neumann 484 485 !-- boundary condition 485 bc_pt_t_val = ( pt_init(nzt ) - pt_init(nzt-1) ) / dzu(nzt)486 bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1) 486 487 487 488 ! … … 539 540 !-- Store humidity gradient at the top boundary for possile Neumann 540 541 !-- boundary condition 541 bc_q_t_val = ( q_init(nzt ) - q_init(nzt-1) ) / dzu(nzt)542 bc_q_t_val = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1) 542 543 543 544 ENDIF … … 795 796 ELSEIF ( bc_pt_t == 'neumann' ) THEN 796 797 ibc_pt_t = 1 798 ELSEIF ( bc_pt_t == 'initial_gradient' ) THEN 799 ibc_pt_t = 2 797 800 ELSE 798 801 IF ( myid == 0 ) THEN … … 804 807 805 808 IF ( surface_heatflux == 9999999.9 ) constant_heatflux = .FALSE. 809 IF ( top_heatflux == 9999999.9 ) THEN 810 constant_top_heatflux = .FALSE. 811 ELSE 812 use_top_fluxes = .TRUE. ! because this is currently the only choice 813 ENDIF 806 814 807 815 ! … … 823 831 PRINT*, ' allowed with pt_surface_initial_change (/=0) = ', & 824 832 pt_surface_initial_change 833 ENDIF 834 CALL local_stop 835 ENDIF 836 837 ! 838 !-- A given temperature at the top implies Dirichlet boundary condition for 839 !-- temperature. In this case specification of a constant heat flux is 840 !-- forbidden. 841 IF ( ibc_pt_t == 0 .AND. constant_top_heatflux .AND. & 842 top_heatflux /= 0.0 ) THEN 843 IF ( myid == 0 ) THEN 844 PRINT*, '+++ check_parameters:' 845 PRINT*, ' boundary_condition: bc_pt_t = ', bc_pt_t 846 PRINT*, ' is not allowed with constant_top_heatflux = .TRUE.' 825 847 ENDIF 826 848 CALL local_stop -
palm/trunk/SOURCE/diffusion_e.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Calculation extended for gridpoint nzt 7 7 ! 8 8 ! Former revisions: … … 56 56 REAL, DIMENSION(:,:), POINTER :: rif 57 57 REAL, DIMENSION(:,:,:), POINTER :: e, km, theta 58 REAL, DIMENSION(nzb+1:nzt -1,nys:nyn) :: dissipation, l, ll58 REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation, l, ll 59 59 60 60 … … 72 72 ENDIF 73 73 74 DO k = nzb_s_inner(j,i)+1, nzt -174 DO k = nzb_s_inner(j,i)+1, nzt 75 75 ! 76 76 !-- Calculate the mixing length (for dissipation) … … 102 102 !-- Calculate the tendency terms 103 103 DO j = nys, nyn 104 DO k = nzb_s_inner(j,i)+1, nzt -1104 DO k = nzb_s_inner(j,i)+1, nzt 105 105 106 106 dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * & … … 130 130 IF ( use_sgs_for_particles ) THEN 131 131 DO j = nys, nyn 132 DO k = nzb_s_inner(j,i)+1, nzt -1132 DO k = nzb_s_inner(j,i)+1, nzt 133 133 diss(k,j,i) = dissipation(k,j) 134 134 ENDDO … … 171 171 REAL, DIMENSION(:,:), POINTER :: rif 172 172 REAL, DIMENSION(:,:,:), POINTER :: e, km, theta 173 REAL, DIMENSION(nzb+1:nzt -1) :: dissipation, l, ll173 REAL, DIMENSION(nzb+1:nzt) :: dissipation, l, ll 174 174 175 175 … … 187 187 ! 188 188 !-- Calculate the mixing length (for dissipation) 189 DO k = nzb_s_inner(j,i)+1, nzt -1189 DO k = nzb_s_inner(j,i)+1, nzt 190 190 dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k) 191 191 IF ( dpt_dz > 0.0 ) THEN … … 234 234 !-- Store dissipation if needed for calculating the sgs particle velocities 235 235 IF ( use_sgs_for_particles ) THEN 236 DO k = nzb_s_inner(j,i)+1, nzt -1236 DO k = nzb_s_inner(j,i)+1, nzt 237 237 diss(k,j,i) = dissipation(k) 238 238 ENDDO -
palm/trunk/SOURCE/diffusion_s.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Calculation extended for gridpoint nzt, fluxes can be given at top, 7 ! +s_flux_t in parameter list, s_flux renamed s_flux_b 7 8 ! 8 9 ! Former revisions: … … 39 40 ! Call for all grid points 40 41 !------------------------------------------------------------------------------! 41 SUBROUTINE diffusion_s( ddzu, ddzw, kh, s, s_flux , tend )42 SUBROUTINE diffusion_s( ddzu, ddzw, kh, s, s_flux_b, s_flux_t, tend ) 42 43 43 44 USE control_parameters … … 51 52 REAL :: ddzu(1:nzt+1), ddzw(1:nzt) 52 53 REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) 53 REAL, DIMENSION(:,:), POINTER :: s_flux 54 REAL, DIMENSION(:,:), POINTER :: s_flux_b, s_flux_t 54 55 REAL, DIMENSION(:,:,:), POINTER :: kh, s 55 56 … … 58 59 ! 59 60 !-- Compute horizontal diffusion 60 DO k = nzb_s_outer(j,i)+1, nzt -161 DO k = nzb_s_outer(j,i)+1, nzt 61 62 62 63 tend(k,j,i) = tend(k,j,i) & … … 97 98 ! 98 99 !-- Compute vertical diffusion. In case that surface fluxes have been 99 !-- presribed or computed, index k starts at nzb+2. 100 DO k = nzb_diff_s_inner(j,i), nzt-1 100 !-- prescribed or computed at bottom and/or top, index k starts/ends at 101 !-- nzb+2 or nzt-1, respectively. 102 DO k = nzb_diff_s_inner(j,i), nzt_diff 101 103 102 104 tend(k,j,i) = tend(k,j,i) & … … 108 110 109 111 ! 110 !-- Vertical diffusion at the first computational gridpoint in &112 !-- Vertical diffusion at the first computational gridpoint along 111 113 !-- z-direction 112 114 IF ( use_surface_fluxes ) THEN … … 118 120 * ( s(k+1,j,i)-s(k,j,i) ) & 119 121 * ddzu(k+1) & 120 + s_flux(j,i) & 122 + s_flux_b(j,i) & 123 ) * ddzw(k) 124 125 ENDIF 126 127 ! 128 !-- Vertical diffusion at the last computational gridpoint along 129 !-- z-direction 130 IF ( use_top_fluxes ) THEN 131 132 k = nzt 133 134 tend(k,j,i) = tend(k,j,i) & 135 + ( - s_flux_t(j,i) & 136 - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) ) & 137 * ( s(k,j,i)-s(k-1,j,i) ) & 138 * ddzu(k) & 121 139 ) * ddzw(k) 122 140 … … 132 150 ! Call for grid point i,j 133 151 !------------------------------------------------------------------------------! 134 SUBROUTINE diffusion_s_ij( i, j, ddzu, ddzw, kh, s, s_flux, tend ) 152 SUBROUTINE diffusion_s_ij( i, j, ddzu, ddzw, kh, s, s_flux_b, s_flux_t, & 153 tend ) 135 154 136 155 USE control_parameters … … 144 163 REAL :: ddzu(1:nzt+1), ddzw(1:nzt) 145 164 REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) 146 REAL, DIMENSION(:,:), POINTER :: s_flux 165 REAL, DIMENSION(:,:), POINTER :: s_flux_b, s_flux_t 147 166 REAL, DIMENSION(:,:,:), POINTER :: kh, s 148 167 149 168 ! 150 169 !-- Compute horizontal diffusion 151 DO k = nzb_s_outer(j,i)+1, nzt -1170 DO k = nzb_s_outer(j,i)+1, nzt 152 171 153 172 tend(k,j,i) = tend(k,j,i) & … … 188 207 ! 189 208 !-- Compute vertical diffusion. In case that surface fluxes have been 190 !-- presribed or computed, index k starts at nzb+2. 191 DO k = nzb_diff_s_inner(j,i), nzt-1 209 !-- prescribed or computed at bottom and/or top, index k starts/ends at 210 !-- nzb+2 or nzt-1, respectively. 211 DO k = nzb_diff_s_inner(j,i), nzt_diff 192 212 193 213 tend(k,j,i) = tend(k,j,i) & … … 199 219 200 220 ! 201 !-- Vertical diffusion at the first computational gridpoint inz-direction221 !-- Vertical diffusion at the first computational gridpoint along z-direction 202 222 IF ( use_surface_fluxes ) THEN 203 223 204 224 k = nzb_s_inner(j,i)+1 205 225 206 tend(k,j,i) = tend(k,j,i) & 207 + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) ) & 208 * ( s(k+1,j,i)-s(k,j,i) ) & 209 * ddzu(k+1) & 210 + s_flux(j,i) & 211 ) * ddzw(k) 226 tend(k,j,i) = tend(k,j,i) + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) ) & 227 * ( s(k+1,j,i)-s(k,j,i) ) & 228 * ddzu(k+1) & 229 + s_flux_b(j,i) & 230 ) * ddzw(k) 212 231 213 232 ENDIF 214 233 234 ! 235 !-- Vertical diffusion at the last computational gridpoint along z-direction 236 IF ( use_top_fluxes ) THEN 237 238 k = nzt 239 240 tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(j,i) & 241 - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) ) & 242 * ( s(k,j,i)-s(k-1,j,i) ) & 243 * ddzu(k) & 244 ) * ddzw(k) 245 246 ENDIF 247 215 248 END SUBROUTINE diffusion_s_ij 216 249 -
palm/trunk/SOURCE/flow_statistics.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! fluxes at top modified (tswst, qswst) 7 7 ! 8 8 ! Former revisions: … … 347 347 + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 348 348 ) * rmask(j,i,sr) 349 ENDDO 350 351 DO k = nzb_diff_s_outer(j,i)-1, nzt_diff 349 352 ! 350 353 !-- Heat flux w"pt" … … 417 420 sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + & 418 421 qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 422 ENDIF 423 ENDIF 424 425 ! 426 !-- Subgridscale fluxes at the top surface 427 IF ( use_top_fluxes ) THEN 428 sums_l(nzt,16,tn) = sums_l(nzt,16,tn) + & 429 tswst(j,i) * rmask(j,i,sr) ! w"pt" 430 sums_l(nzt,58,tn) = sums_l(nzt,58,tn) + & 431 0.0 * rmask(j,i,sr) ! u"pt" 432 sums_l(nzt,61,tn) = sums_l(nzt,61,tn) + & 433 0.0 * rmask(j,i,sr) ! v"pt" 434 IF ( moisture ) THEN 435 sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + & 436 qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 437 IF ( cloud_physics ) THEN 438 sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( & 439 ( 1.0 + 0.61 * q(nzt,j,i) ) * & 440 tswst(j,i) + 0.61 * pt(nzt,j,i) * & 441 qsws(j,i) & 442 ) 443 ! 444 !-- Formula does not work if ql(nzb) /= 0.0 445 sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + & ! w"q" (w"qv") 446 qswst(j,i) * rmask(j,i,sr) 447 ENDIF 448 ENDIF 449 IF ( passive_scalar ) THEN 450 sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + & 451 qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 419 452 ENDIF 420 453 ENDIF -
palm/trunk/SOURCE/header.f90
r4 r19 345 345 ENDIF 346 346 IF ( ibc_pt_t == 0 ) THEN 347 roben = TRIM( roben ) // ' pt(nzt) = pt_top' 348 ELSE 349 roben = TRIM( roben ) // ' pt(nzt) = pt(nzt-1) + dpt/dz' 347 roben = TRIM( roben ) // ' pt(nzt+1) = pt_top' 348 ELSEIF( ibc_pt_t == 1 ) THEN 349 roben = TRIM( roben ) // ' pt(nzt+1) = pt(nzt)' 350 ELSEIF( ibc_pt_t == 2 ) THEN 351 roben = TRIM( roben ) // ' pt(nzt+1) = pt(nzt) + dpt/dz_ini' 350 352 ENDIF 351 353 … … 404 406 IF ( passive_scalar .AND. constant_waterflux ) THEN 405 407 WRITE ( io, 313 ) surface_waterflux 408 ENDIF 409 ENDIF 410 411 IF ( use_top_fluxes ) THEN 412 WRITE ( io, 304 ) 413 IF ( constant_top_heatflux ) THEN 414 WRITE ( io, 306 ) top_heatflux 415 ENDIF 416 IF ( moisture .OR. passive_scalar ) THEN 417 WRITE ( io, 315 ) 406 418 ENDIF 407 419 ENDIF … … 1220 1232 ' B. bound.: ',A/ & 1221 1233 ' T. bound.: ',A) 1222 303 FORMAT (/' Surface fluxes are used in diffusion terms at k=1') 1223 305 FORMAT (//' Prandtl-Layer between surface and first computational ', & 1224 'u,v-level:'// & 1234 303 FORMAT (/' Bottom surface fluxes are used in diffusion terms at k=1') 1235 304 FORMAT (/' Top surface fluxes are used in diffusion terms at k=nzt') 1236 305 FORMAT (//' Prandtl-Layer between bottom surface and first ', & 1237 'computational u,v-level:'// & 1225 1238 ' zp = ',F6.2,' m z0 = ',F6.4,' m kappa = ',F4.2/ & 1226 1239 ' Rif value range: ',F6.2,' <= rif <=',F6.2) … … 1234 1247 313 FORMAT (' Predefined constant scalar flux: ',E10.3,' kg/(m**2 s)') 1235 1248 314 FORMAT (' Predefined scalar value at the surface') 1249 315 FORMAT (' Humidity / scalar flux at top surface is 0.0') 1236 1250 317 FORMAT (//' Lateral boundaries:'/ & 1237 1251 ' left/right: ',A/ & -
palm/trunk/SOURCE/impact_of_latent_heat.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Calculation extended for gridpoint nzt 7 7 ! 8 8 ! Former revisions: … … 53 53 DO i = nxl, nxr 54 54 DO j = nys, nyn 55 DO k = nzb_2d(j,i)+1, nzt -155 DO k = nzb_2d(j,i)+1, nzt 56 56 57 57 IF ( ql(k,j,i) > ql_crit ) THEN … … 87 87 88 88 89 DO k = nzb_2d(j,i)+1, nzt -189 DO k = nzb_2d(j,i)+1, nzt 90 90 91 91 IF ( ql(k,j,i) > ql_crit ) THEN -
palm/trunk/SOURCE/init_3d_model.f90
r4 r19 7 7 ! Actual revisions: 8 8 ! ----------------- 9 ! 9 ! +handling of top fluxes 10 10 ! 11 11 ! Former revisions: … … 77 77 ALLOCATE( km_damp_x(nxl-1:nxr+1), km_damp_y(nys-1:nyn+1) ) 78 78 79 ALLOCATE( rif_1(nys-1:nyn+1,nxl-1:nxr+1), shf_1(nys-1:nyn+1,nxl-1:nxr+1), 80 ts(nys-1:nyn+1,nxl-1:nxr+1), us(nys-1:nyn+1,nxl-1:nxr+1),&81 us ws_1(nys-1:nyn+1,nxl-1:nxr+1), vsws_1(nys-1:nyn+1,nxl-1:nxr+1),&82 z0(nys-1:nyn+1,nxl-1:nxr+1) )79 ALLOCATE( rif_1(nys-1:nyn+1,nxl-1:nxr+1), shf_1(nys-1:nyn+1,nxl-1:nxr+1), & 80 ts(nys-1:nyn+1,nxl-1:nxr+1), tswst_1(nys-1:nyn+1,nxl-1:nxr+1), & 81 us(nys-1:nyn+1,nxl-1:nxr+1), usws_1(nys-1:nyn+1,nxl-1:nxr+1), & 82 vsws_1(nys-1:nyn+1,nxl-1:nxr+1), z0(nys-1:nyn+1,nxl-1:nxr+1) ) 83 83 84 84 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 85 85 ! 86 86 !-- Leapfrog scheme needs two timelevels of diffusion quantities 87 ALLOCATE( rif_2(nys-1:nyn+1,nxl-1:nxr+1), & 88 shf_2(nys-1:nyn+1,nxl-1:nxr+1), & 89 usws_2(nys-1:nyn+1,nxl-1:nxr+1), & 87 ALLOCATE( rif_2(nys-1:nyn+1,nxl-1:nxr+1), & 88 shf_2(nys-1:nyn+1,nxl-1:nxr+1), & 89 tswst_2(nys-1:nyn+1,nxl-1:nxr+1), & 90 usws_2(nys-1:nyn+1,nxl-1:nxr+1), & 90 91 vsws_2(nys-1:nyn+1,nxl-1:nxr+1) ) 91 92 ENDIF … … 121 122 !-- 2D-moisture/scalar arrays 122 123 ALLOCATE ( qs(nys-1:nyn+1,nxl-1:nxr+1), & 123 qsws_1(nys-1:nyn+1,nxl-1:nxr+1) ) 124 qsws_1(nys-1:nyn+1,nxl-1:nxr+1), & 125 qswst_1(nys-1:nyn+1,nxl-1:nxr+1) ) 124 126 125 127 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 126 ALLOCATE( qsws_2(nys-1:nyn+1,nxl-1:nxr+1) ) 128 ALLOCATE( qsws_2(nys-1:nyn+1,nxl-1:nxr+1), & 129 qswst_2(nys-1:nyn+1,nxl-1:nxr+1) ) 127 130 ENDIF 128 131 ! … … 177 180 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 178 181 179 rif_m => rif_1; rif => rif_2 180 shf_m => shf_1; shf => shf_2 181 usws_m => usws_1; usws => usws_2 182 vsws_m => vsws_1; vsws => vsws_2 182 rif_m => rif_1; rif => rif_2 183 shf_m => shf_1; shf => shf_2 184 tswst_m => tswst_1; tswst => tswst_2 185 usws_m => usws_1; usws => usws_2 186 vsws_m => vsws_1; vsws => vsws_2 183 187 e_m => e_1; e => e_2; e_p => e_3; te_m => e_3 184 188 kh_m => kh_1; kh => kh_2 … … 190 194 191 195 IF ( moisture .OR. passive_scalar ) THEN 192 qsws_m => qsws_1; qsws => qsws_2 196 qsws_m => qsws_1; qsws => qsws_2 197 qswst_m => qswst_1; qswst => qswst_2 193 198 q_m => q_1; q => q_2; q_p => q_3; tq_m => q_3 194 199 IF ( moisture ) vpt_m => vpt_1; vpt => vpt_2 … … 202 207 ELSE 203 208 204 rif => rif_1 205 shf => shf_1 206 usws => usws_1 207 vsws => vsws_1 208 e => e_1; e_p => e_2; te_m => e_3; e_m => e_3 209 kh => kh_1 210 km => km_1 211 pt => pt_1; pt_p => pt_2; tpt_m => pt_3; pt_m => pt_3 212 u => u_1; u_p => u_2; tu_m => u_3; u_m => u_3 213 v => v_1; v_p => v_2; tv_m => v_3; v_m => v_3 214 w => w_1; w_p => w_2; tw_m => w_3; w_m => w_3 209 rif => rif_1 210 shf => shf_1 211 tswst => tswst_1 212 usws => usws_1 213 vsws => vsws_1 214 e => e_1; e_p => e_2; te_m => e_3; e_m => e_3 215 kh => kh_1 216 km => km_1 217 pt => pt_1; pt_p => pt_2; tpt_m => pt_3; pt_m => pt_3 218 u => u_1; u_p => u_2; tu_m => u_3; u_m => u_3 219 v => v_1; v_p => v_2; tv_m => v_3; v_m => v_3 220 w => w_1; w_p => w_2; tw_m => w_3; w_m => w_3 215 221 216 222 IF ( moisture .OR. passive_scalar ) THEN 217 223 qsws => qsws_1 224 qswst => qswst_1 218 225 q => q_1; q_p => q_2; tq_m => q_3; q_m => q_3 219 226 IF ( moisture ) vpt => vpt_1 … … 453 460 454 461 ! 455 !-- Initialize surface fluxes462 !-- Initialize fluxes at bottom surface 456 463 IF ( use_surface_fluxes ) THEN 457 464 … … 486 493 ENDIF 487 494 ENDIF 495 496 ENDIF 497 498 ! 499 !-- Initialize fluxes at top surface 500 !-- Currently, only the heatflux can be prescribed. The latent flux is 501 !-- zeri in this case! 502 IF ( use_top_fluxes ) THEN 503 504 IF ( constant_top_heatflux ) THEN 505 ! 506 !-- Heat flux is prescribed 507 tswst = top_heatflux 508 IF ( ASSOCIATED( tswst_m ) ) tswst_m = tswst 509 510 IF ( moisture .OR. passive_scalar ) THEN 511 qswst = 0.0 512 IF ( ASSOCIATED( qswst_m ) ) qswst_m = qswst 513 ENDIF 514 ENDIF 488 515 489 516 ENDIF -
palm/trunk/SOURCE/init_grid.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Setting of nzt_diff 7 7 ! 8 8 ! Former revisions: … … 202 202 203 203 ! 204 !-- Define vertical gridpoint from which on the usual finite difference204 !-- Define vertical gridpoint from (or to) which on the usual finite difference 205 205 !-- form (which does not use surface fluxes) is applied 206 206 IF ( prandtl_layer .OR. use_surface_fluxes ) THEN … … 208 208 ELSE 209 209 nzb_diff = nzb + 1 210 ENDIF 211 IF ( use_top_fluxes ) THEN 212 nzt_diff = nzt - 1 213 ELSE 214 nzt_diff = nzt 210 215 ENDIF 211 216 -
palm/trunk/SOURCE/init_pt_anomaly.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Calculation extended for gridpoint nzt 7 7 ! 8 8 ! Former revisions: … … 47 47 DO i = nxl, nxr 48 48 DO j = nys, nyn 49 DO k = nzb+1, nzt -149 DO k = nzb+1, nzt 50 50 x = ( i - ic ) * dx 51 51 y = ( j - jc ) * dy -
palm/trunk/SOURCE/modules.f90
r4 r19 5 5 ! Actual revisions: 6 6 ! ----------------- 7 ! 7 ! +constant_top_heatflux, top_heatflux, use_top_fluxes, +arrays for top fluxes, 8 ! +nzt_diff, default of bc_pt_t renamed "initial_gradient" 9 ! Bugfix: p is not a pointer 8 10 ! 9 11 ! Former revisions: … … 75 77 76 78 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET :: & 77 qsws_1, qsws_2, rif_1, rif_2, shf_1, shf_2, usws_1, usws_2,&78 vsws_1, vsws_279 qsws_1, qsws_2, qswst_1, qswst_2, rif_1, rif_2, shf_1, shf_2, & 80 tswst_1, tswst_2, usws_1, usws_2, vsws_1, vsws_2 79 81 80 82 REAL, DIMENSION(:,:), POINTER :: & 81 qsws, qsws_m, rif, rif_m, shf, shf_m, usws, usws_m, vsws, vsws_m 83 qsws, qsws_m, qswst, qswst_m, rif, rif_m, shf, shf_m, tswst, & 84 tswst_m, usws, usws_m, vsws, vsws_m 82 85 83 86 REAL, DIMENSION(:,:,:), ALLOCATABLE :: & … … 88 91 89 92 REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & 90 e_1, e_2, e_3, kh_1, kh_2, km_1, km_2, p t_1, pt_2, pt_3, q_1, q_2,&91 q_ 3, ql_1, ql_2, u_1, u_2, u_3, v_1, v_2, v_3, vpt_1, vpt_2, w_1, &92 w_ 2, w_393 e_1, e_2, e_3, kh_1, kh_2, km_1, km_2, p, pt_1, pt_2, pt_3, q_1, & 94 q_2, q_3, ql_1, ql_2, u_1, u_2, u_3, v_1, v_2, v_3, vpt_1, vpt_2, & 95 w_1, w_2, w_3 93 96 94 97 REAL, DIMENSION(:,:,:), POINTER :: & 95 e, e_m, e_p, kh, kh_m, km, km_m, p , pt, pt_m, pt_p, q, q_m, q_p, ql,&98 e, e_m, e_p, kh, kh_m, km, km_m, pt, pt_m, pt_p, q, q_m, q_p, ql, & 96 99 ql_c, te_m, tpt_m, tq_m, tu_m, tv_m, tw_m, u, u_m, u_p, v, v_m, v_p, & 97 100 vpt, vpt_m, w, w_m, w_p … … 201 204 scalar_advec = 'pw-scheme' 202 205 CHARACTER (LEN=20) :: bc_e_b = 'neumann', bc_lr = 'cyclic', & 203 bc_ns = 'cyclic', &204 bc_p_ b = 'neumann', bc_p_t= 'dirichlet', &205 bc_pt_ b = 'dirichlet', bc_pt_t = 'neumann', &206 bc_ns = 'cyclic', bc_p_b = 'neumann', & 207 bc_p_t = 'dirichlet', bc_pt_b = 'dirichlet', & 208 bc_pt_t = 'initial_gradient', & 206 209 bc_q_b = 'dirichlet', bc_q_t = 'neumann', & 207 210 bc_s_b = 'dirichlet', bc_s_t = 'neumann', & … … 263 266 cloud_droplets = .FALSE., cloud_physics = .FALSE., & 264 267 conserve_volume_flow = .FALSE., constant_diffusion = .FALSE., & 265 constant_heatflux = .TRUE., constant_waterflux = .TRUE., & 266 create_disturbances = .TRUE., cut_spline_overshoot = .TRUE., & 268 constant_heatflux = .TRUE., constant_top_heatflux = .TRUE., & 269 constant_waterflux = .TRUE., create_disturbances = .TRUE., & 270 cut_spline_overshoot = .TRUE., & 267 271 data_output_2d_on_each_pe = .TRUE., do2d_at_begin = .FALSE., & 268 272 do3d_at_begin = .FALSE., do3d_compress = .FALSE., & … … 270 274 disturbance_created = .FALSE., & 271 275 first_call_advec_particles = .TRUE., & 272 force_print_header = .FALSE., galilei_transformation = .FALSE., 276 force_print_header = .FALSE., galilei_transformation = .FALSE.,& 273 277 inflow_l = .FALSE., inflow_n = .FALSE., inflow_r = .FALSE., & 274 278 inflow_s = .FALSE., iso2d_output = .FALSE., & … … 281 285 random_heatflux = .FALSE., run_control_header = .FALSE., & 282 286 sloping_surface = .FALSE., stop_dt = .FALSE., & 283 terminate_run = .FALSE., use_prior_plot1d_parameters = .FALSE., & 284 use_surface_fluxes = .FALSE., use_ug_for_galilei_tr = .TRUE., & 287 terminate_run = .FALSE., use_prior_plot1d_parameters = .FALSE.,& 288 use_surface_fluxes = .FALSE., use_top_fluxes = .FALSE., & 289 use_ug_for_galilei_tr = .TRUE., & 285 290 use_upstream_for_tke = .FALSE., wall_adjustment = .TRUE. 286 291 … … 336 341 time_do_sla = 0.0, time_dvrp = 0.0, time_prel = 0.0, & 337 342 time_restart = 9999999.9, time_run_control = 0.0, & 338 ug_surface = 0.0, u_gtrans = 0.0, ups_limit_e= 0.0, &339 ups_limit_ pt = 0.0, ups_limit_u = 0.0, ups_limit_v= 0.0, &340 ups_limit_ w = 0.0, vg_surface = 0.0, v_gtrans= 0.0, &341 wall_adjustment_factor = 1.8, z_max_do1d = -1.0, &343 top_heatflux = 9999999.9, ug_surface = 0.0, u_gtrans = 0.0, & 344 ups_limit_e = 0.0, ups_limit_pt = 0.0, ups_limit_u = 0.0, & 345 ups_limit_v = 0.0, ups_limit_w = 0.0, vg_surface = 0.0, & 346 v_gtrans = 0.0, wall_adjustment_factor = 1.8, z_max_do1d = -1.0, & 342 347 z_max_do1d_normalized = -1.0, z_max_do2d = -1.0 343 348 … … 500 505 INTEGER :: ngp_sums, nnx, nx = 0, nxa, nxl, nxr, nxra, nny, ny = 0, nya, & 501 506 nyn, nyna, nys, nnz, nz = 0, nza, nzb, nzb_diff, nzt, nzta, & 502 uxrp = 0, vynp = 0507 nzt_diff, uxrp = 0, vynp = 0 503 508 504 509 INTEGER, DIMENSION(:), ALLOCATABLE :: & -
palm/trunk/SOURCE/parin.f90
r7 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! +top_heatflux in inipar 7 7 ! 8 8 ! Former revisions: … … 71 71 s_surface, s_surface_initial_change, & 72 72 s_vertical_gradient, s_vertical_gradient_level, & 73 t imestep_scheme, topography, ug_surface, &73 top_heatflux, timestep_scheme, topography, ug_surface, & 74 74 ug_vertical_gradient, ug_vertical_gradient_level, & 75 75 ups_limit_e, ups_limit_pt, ups_limit_u, ups_limit_v, & -
palm/trunk/SOURCE/production_e.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Calculation extended for gridpoint nzt, extended for given temperature / 7 ! humidity fluxes at the top 7 8 ! 8 9 ! Former revisions: … … 67 68 68 69 DO j = nys, nyn 69 DO k = nzb_diff_s_outer(j,i), nzt -170 DO k = nzb_diff_s_outer(j,i), nzt 70 71 71 72 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx … … 327 328 DO j = nys, nyn 328 329 329 DO k = nzb_diff_s_inner(j,i), nzt -1330 DO k = nzb_diff_s_inner(j,i), nzt_diff 330 331 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * & 331 332 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 332 333 ENDDO 334 333 335 IF ( use_surface_fluxes ) THEN 334 336 k = nzb_diff_s_inner(j,i)-1 … … 336 338 ENDIF 337 339 340 IF ( use_top_fluxes ) THEN 341 k = nzt 342 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i) 343 ENDIF 344 338 345 ENDDO 339 346 … … 342 349 DO j = nys, nyn 343 350 344 DO k = nzb_diff_s_inner(j,i), nzt -1351 DO k = nzb_diff_s_inner(j,i), nzt_diff 345 352 346 353 IF ( .NOT. cloud_physics ) THEN … … 375 382 DO j = nys, nyn 376 383 377 k = nzb_diff_s_inner(j,i) -1384 k = nzb_diff_s_inner(j,i) 378 385 379 386 IF ( .NOT. cloud_physics ) THEN … … 402 409 ENDIF 403 410 411 IF ( use_top_fluxes ) THEN 412 413 DO j = nys, nyn 414 415 k = nzt 416 417 IF ( .NOT. cloud_physics ) THEN 418 k1 = 1.0 + 0.61 * q(k,j,i) 419 k2 = 0.61 * pt(k,j,i) 420 ELSE 421 IF ( ql(k,j,i) == 0.0 ) THEN 422 k1 = 1.0 + 0.61 * q(k,j,i) 423 k2 = 0.61 * pt(k,j,i) 424 ELSE 425 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 426 temp = theta * t_d_pt(k) 427 k1 = ( 1.0 - q(k,j,i) + 1.61 * & 428 ( q(k,j,i) - ql(k,j,i) ) * & 429 ( 1.0 + 0.622 * l_d_r / temp ) ) / & 430 ( 1.0 + 0.622 * l_d_r * l_d_cp * & 431 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 432 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 433 ENDIF 434 ENDIF 435 436 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 437 ( k1* tswst(j,i) + k2 * qswst(j,i) ) 438 ENDDO 439 440 ENDIF 441 404 442 ENDIF 405 443 … … 430 468 ! 431 469 !-- Calculate TKE production by shear 432 DO k = nzb_diff_s_outer(j,i), nzt -1470 DO k = nzb_diff_s_outer(j,i), nzt 433 471 434 472 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx … … 657 695 IF ( .NOT. moisture ) THEN 658 696 659 DO k = nzb_diff_s_inner(j,i), nzt -1697 DO k = nzb_diff_s_inner(j,i), nzt_diff 660 698 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * & 661 699 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 662 700 ENDDO 701 663 702 IF ( use_surface_fluxes ) THEN 664 703 k = nzb_diff_s_inner(j,i)-1 … … 666 705 ENDIF 667 706 707 IF ( use_top_fluxes ) THEN 708 k = nzt 709 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i) 710 ENDIF 711 668 712 ELSE 669 713 670 DO k = nzb_diff_s_inner(j,i), nzt -1714 DO k = nzb_diff_s_inner(j,i), nzt_diff 671 715 672 716 IF ( .NOT. cloud_physics ) THEN … … 694 738 ) * dd2zu(k) 695 739 ENDDO 740 696 741 IF ( use_surface_fluxes ) THEN 697 742 k = nzb_diff_s_inner(j,i)-1 … … 720 765 ENDIF 721 766 767 IF ( use_top_fluxes ) THEN 768 k = nzt 769 770 IF ( .NOT. cloud_physics ) THEN 771 k1 = 1.0 + 0.61 * q(k,j,i) 772 k2 = 0.61 * pt(k,j,i) 773 ELSE 774 IF ( ql(k,j,i) == 0.0 ) THEN 775 k1 = 1.0 + 0.61 * q(k,j,i) 776 k2 = 0.61 * pt(k,j,i) 777 ELSE 778 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 779 temp = theta * t_d_pt(k) 780 k1 = ( 1.0 - q(k,j,i) + 1.61 * & 781 ( q(k,j,i) - ql(k,j,i) ) * & 782 ( 1.0 + 0.622 * l_d_r / temp ) ) / & 783 ( 1.0 + 0.622 * l_d_r * l_d_cp * & 784 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 785 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 786 ENDIF 787 ENDIF 788 789 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 790 ( k1* tswst(j,i) + k2 * qswst(j,i) ) 791 ENDIF 792 722 793 ENDIF 723 794 -
palm/trunk/SOURCE/prognostic_equations.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Calculation of e, q, and pt extended for gridpoint nzt, 7 ! handling of given temperature/humidity/scalar fluxes at top surface 7 8 ! 8 9 ! Former revisions: … … 22 23 ! Description: 23 24 ! ------------ 24 ! Solving the prognostic equations and advecting particles.25 ! Solving the prognostic equations. 25 26 !------------------------------------------------------------------------------! 26 27 … … 333 334 !-- Tendency terms 334 335 IF ( scalar_advec == 'bc-scheme' ) THEN 335 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, t end )336 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend ) 336 337 ELSE 337 338 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN … … 346 347 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 347 348 THEN 348 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, tend ) 349 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, & 350 tswst_m, tend ) 349 351 ELSE 350 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, t end )352 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend ) 351 353 ENDIF 352 354 ENDIF … … 368 370 ! 369 371 !-- Prognostic equation for potential temperature 370 DO k = nzb_s_inner(j,i)+1, nzt -1372 DO k = nzb_s_inner(j,i)+1, nzt 371 373 pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + & 372 374 dt_3d * ( & … … 380 382 IF ( timestep_scheme(1:5) == 'runge' ) THEN 381 383 IF ( intermediate_timestep_count == 1 ) THEN 382 DO k = nzb_s_inner(j,i)+1, nzt -1384 DO k = nzb_s_inner(j,i)+1, nzt 383 385 tpt_m(k,j,i) = tend(k,j,i) 384 386 ENDDO 385 387 ELSEIF ( intermediate_timestep_count < & 386 388 intermediate_timestep_count_max ) THEN 387 DO k = nzb_s_inner(j,i)+1, nzt -1389 DO k = nzb_s_inner(j,i)+1, nzt 388 390 tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i) 389 391 ENDDO … … 429 431 !-- Tendency-terms 430 432 IF ( scalar_advec == 'bc-scheme' ) THEN 431 CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, tend )433 CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, tend ) 432 434 ELSE 433 435 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN … … 443 445 THEN 444 446 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, & 447 qswst_m, tend ) 448 ELSE 449 CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, & 445 450 tend ) 446 ELSE447 CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, tend )448 451 ENDIF 449 452 ENDIF … … 459 462 ! 460 463 !-- Prognostic equation for total water content / scalar 461 DO k = nzb_s_inner(j,i)+1, nzt -1464 DO k = nzb_s_inner(j,i)+1, nzt 462 465 q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + & 463 466 dt_3d * ( & … … 471 474 IF ( timestep_scheme(1:5) == 'runge' ) THEN 472 475 IF ( intermediate_timestep_count == 1 ) THEN 473 DO k = nzb_s_inner(j,i)+1, nzt -1476 DO k = nzb_s_inner(j,i)+1, nzt 474 477 tq_m(k,j,i) = tend(k,j,i) 475 478 ENDDO 476 479 ELSEIF ( intermediate_timestep_count < & 477 480 intermediate_timestep_count_max ) THEN 478 DO k = nzb_s_inner(j,i)+1, nzt -1481 DO k = nzb_s_inner(j,i)+1, nzt 479 482 tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i) 480 483 ENDDO … … 577 580 !-- reasons in the course of the integration. In such cases the old TKE 578 581 !-- value is reduced by 90%. 579 DO k = nzb_s_inner(j,i)+1, nzt -1582 DO k = nzb_s_inner(j,i)+1, nzt 580 583 e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + & 581 584 dt_3d * ( & … … 589 592 IF ( timestep_scheme(1:5) == 'runge' ) THEN 590 593 IF ( intermediate_timestep_count == 1 ) THEN 591 DO k = nzb_s_inner(j,i)+1, nzt -1594 DO k = nzb_s_inner(j,i)+1, nzt 592 595 te_m(k,j,i) = tend(k,j,i) 593 596 ENDDO 594 597 ELSEIF ( intermediate_timestep_count < & 595 598 intermediate_timestep_count_max ) THEN 596 DO k = nzb_s_inner(j,i)+1, nzt -1599 DO k = nzb_s_inner(j,i)+1, nzt 597 600 te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i) 598 601 ENDDO … … 809 812 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 810 813 THEN 811 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, tend ) 814 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, & 815 tswst_m, tend ) 812 816 ELSE 813 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, t end )817 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend ) 814 818 ENDIF 815 819 … … 830 834 ! 831 835 !-- Prognostic equation for potential temperature 832 DO k = nzb_s_inner(j,i)+1, nzt -1836 DO k = nzb_s_inner(j,i)+1, nzt 833 837 pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +& 834 838 dt_3d * ( & … … 842 846 IF ( timestep_scheme(1:5) == 'runge' ) THEN 843 847 IF ( intermediate_timestep_count == 1 ) THEN 844 DO k = nzb_s_inner(j,i)+1, nzt -1848 DO k = nzb_s_inner(j,i)+1, nzt 845 849 tpt_m(k,j,i) = tend(k,j,i) 846 850 ENDDO 847 851 ELSEIF ( intermediate_timestep_count < & 848 852 intermediate_timestep_count_max ) THEN 849 DO k = nzb_s_inner(j,i)+1, nzt -1853 DO k = nzb_s_inner(j,i)+1, nzt 850 854 tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + & 851 855 5.3125 * tpt_m(k,j,i) … … 870 874 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )& 871 875 THEN 872 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, tend ) 876 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, & 877 qswst_m, tend ) 873 878 ELSE 874 CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, tend ) 879 CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, & 880 tend ) 875 881 ENDIF 876 882 … … 885 891 ! 886 892 !-- Prognostic equation for total water content / scalar 887 DO k = nzb_s_inner(j,i)+1, nzt -1893 DO k = nzb_s_inner(j,i)+1, nzt 888 894 q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +& 889 895 dt_3d * ( & … … 897 903 IF ( timestep_scheme(1:5) == 'runge' ) THEN 898 904 IF ( intermediate_timestep_count == 1 ) THEN 899 DO k = nzb_s_inner(j,i)+1, nzt -1905 DO k = nzb_s_inner(j,i)+1, nzt 900 906 tq_m(k,j,i) = tend(k,j,i) 901 907 ENDDO 902 908 ELSEIF ( intermediate_timestep_count < & 903 909 intermediate_timestep_count_max ) THEN 904 DO k = nzb_s_inner(j,i)+1, nzt -1910 DO k = nzb_s_inner(j,i)+1, nzt 905 911 tq_m(k,j,i) = -9.5625 * tend(k,j,i) + & 906 912 5.3125 * tq_m(k,j,i) … … 951 957 !-- reasons in the course of the integration. In such cases the old 952 958 !-- TKE value is reduced by 90%. 953 DO k = nzb_s_inner(j,i)+1, nzt -1959 DO k = nzb_s_inner(j,i)+1, nzt 954 960 e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +& 955 961 dt_3d * ( & … … 963 969 IF ( timestep_scheme(1:5) == 'runge' ) THEN 964 970 IF ( intermediate_timestep_count == 1 ) THEN 965 DO k = nzb_s_inner(j,i)+1, nzt -1971 DO k = nzb_s_inner(j,i)+1, nzt 966 972 te_m(k,j,i) = tend(k,j,i) 967 973 ENDDO 968 974 ELSEIF ( intermediate_timestep_count < & 969 975 intermediate_timestep_count_max ) THEN 970 DO k = nzb_s_inner(j,i)+1, nzt -1976 DO k = nzb_s_inner(j,i)+1, nzt 971 977 te_m(k,j,i) = -9.5625 * tend(k,j,i) + & 972 978 5.3125 * te_m(k,j,i) … … 1252 1258 !-- pt-tendency terms with no communication 1253 1259 IF ( scalar_advec == 'bc-scheme' ) THEN 1254 CALL diffusion_s( ddzu, ddzw, kh, pt, shf, t end )1260 CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, tend ) 1255 1261 ELSE 1256 1262 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN … … 1264 1270 ENDIF 1265 1271 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 1266 CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, t end )1272 CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, tend ) 1267 1273 ELSE 1268 CALL diffusion_s( ddzu, ddzw, kh, pt, shf, t end )1274 CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, tend ) 1269 1275 ENDIF 1270 1276 ENDIF … … 1288 1294 DO i = nxl, nxr 1289 1295 DO j = nys, nyn 1290 DO k = nzb_s_inner(j,i)+1, nzt -11296 DO k = nzb_s_inner(j,i)+1, nzt 1291 1297 pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + & 1292 1298 dt_3d * ( & … … 1304 1310 DO i = nxl, nxr 1305 1311 DO j = nys, nyn 1306 DO k = nzb_s_inner(j,i)+1, nzt -11312 DO k = nzb_s_inner(j,i)+1, nzt 1307 1313 tpt_m(k,j,i) = tend(k,j,i) 1308 1314 ENDDO … … 1313 1319 DO i = nxl, nxr 1314 1320 DO j = nys, nyn 1315 DO k = nzb_s_inner(j,i)+1, nzt -11321 DO k = nzb_s_inner(j,i)+1, nzt 1316 1322 tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i) 1317 1323 ENDDO … … 1352 1358 !-- Scalar/q-tendency terms with no communication 1353 1359 IF ( scalar_advec == 'bc-scheme' ) THEN 1354 CALL diffusion_s( ddzu, ddzw, kh, q, qsws, tend )1360 CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, tend ) 1355 1361 ELSE 1356 1362 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN … … 1364 1370 ENDIF 1365 1371 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 1366 CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, tend )1372 CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, qswst_m, tend ) 1367 1373 ELSE 1368 CALL diffusion_s( ddzu, ddzw, kh, q, qsws, tend )1374 CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, tend ) 1369 1375 ENDIF 1370 1376 ENDIF … … 1382 1388 DO i = nxl, nxr 1383 1389 DO j = nys, nyn 1384 DO k = nzb_s_inner(j,i)+1, nzt -11390 DO k = nzb_s_inner(j,i)+1, nzt 1385 1391 q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + & 1386 1392 dt_3d * ( & … … 1398 1404 DO i = nxl, nxr 1399 1405 DO j = nys, nyn 1400 DO k = nzb_s_inner(j,i)+1, nzt -11406 DO k = nzb_s_inner(j,i)+1, nzt 1401 1407 tq_m(k,j,i) = tend(k,j,i) 1402 1408 ENDDO … … 1407 1413 DO i = nxl, nxr 1408 1414 DO j = nys, nyn 1409 DO k = nzb_s_inner(j,i)+1, nzt -11415 DO k = nzb_s_inner(j,i)+1, nzt 1410 1416 tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i) 1411 1417 ENDDO … … 1503 1509 DO i = nxl, nxr 1504 1510 DO j = nys, nyn 1505 DO k = nzb_s_inner(j,i)+1, nzt -11511 DO k = nzb_s_inner(j,i)+1, nzt 1506 1512 e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + & 1507 1513 dt_3d * ( & … … 1519 1525 DO i = nxl, nxr 1520 1526 DO j = nys, nyn 1521 DO k = nzb_s_inner(j,i)+1, nzt -11527 DO k = nzb_s_inner(j,i)+1, nzt 1522 1528 te_m(k,j,i) = tend(k,j,i) 1523 1529 ENDDO … … 1528 1534 DO i = nxl, nxr 1529 1535 DO j = nys, nyn 1530 DO k = nzb_s_inner(j,i)+1, nzt -11536 DO k = nzb_s_inner(j,i)+1, nzt 1531 1537 te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i) 1532 1538 ENDDO -
palm/trunk/SOURCE/read_3d_binary.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! +qswst, qswst_m, tswst, tswst_m 7 7 ! 8 8 ! Former revisions: … … 249 249 CASE ( 'qsws_m' ) 250 250 READ ( 13 ) qsws_m 251 CASE ( 'qswst' ) 252 READ ( 13 ) qswst 253 CASE ( 'qswst_m' ) 254 READ ( 13 ) qswst_m 251 255 CASE ( 'qv_av' ) 252 256 ALLOCATE( qv_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) … … 266 270 CASE ( 'shf_m' ) 267 271 READ ( 13 ) shf_m 272 CASE ( 'tswst' ) 273 READ ( 13 ) tswst 274 CASE ( 'tswst_m' ) 275 READ ( 13 ) tswst_m 268 276 CASE ( 'spectrum_x' ) 269 277 READ ( 13 ) spectrum_x -
palm/trunk/SOURCE/read_var_list.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! +top_heatflux 7 7 ! 8 8 ! Former revisions: … … 340 340 CASE ( 'topography' ) 341 341 READ ( 13 ) topography 342 CASE ( 'top_heatflux' ) 343 READ ( 13 ) top_heatflux 342 344 CASE ( 'tsc' ) 343 345 READ ( 13 ) tsc -
palm/trunk/SOURCE/spline_z.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Boundary condition for pt at top adjusted 7 7 ! 8 8 ! Former revisions: … … 156 156 ! 157 157 !-- Top boundary for temperature 158 IF ( ibc_pt_t == 1 ) THEN 159 vad(nzt,:) = vad(nzt-1,:) + bc_pt_t_val * dz_spline(nzt) 158 IF ( ibc_pt_t == 0 ) THEN 159 vad(nzt+1,:) = pt(nzt+1,nys:nyn,i) 160 ELSEIF ( ibc_pt_t == 1 ) THEN 161 vad(nzt+1,:) = vad(nzt,:) 162 ELSEIF ( ibc_pt_t == 2 ) THEN 160 163 vad(nzt+1,:) = vad(nzt,:) + bc_pt_t_val * dz_spline(nzt+1) 161 ELSE162 vad(nzt,:) = pt(nzt,nys:nyn,i)163 vad(nzt+1,:) = pt(nzt+1,nys:nyn,i)164 164 ENDIF 165 165 -
palm/trunk/SOURCE/swap_timelevel.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Swaping of top fluxes 7 7 ! 8 8 ! Former revisions: … … 122 122 rif_m => rif_1; rif => rif_2 123 123 ENDIF 124 IF ( use_top_fluxes ) THEN 125 tswst_m => tswst_1; tswst => tswst_2 126 IF ( moisture .OR. passive_scalar ) THEN 127 qswst_m => qswst_1; qswst => qswst_2 128 ENDIF 129 ENDIF 124 130 ENDIF 125 131 … … 161 167 rif_m => rif_2; rif => rif_1 162 168 ENDIF 169 IF ( use_top_fluxes ) THEN 170 tswst_m => tswst_2; tswst => tswst_1 171 IF ( moisture .OR. passive_scalar ) THEN 172 qswst_m => qswst_2; qswst => qswst_1 173 ENDIF 174 ENDIF 163 175 ENDIF 164 176 -
palm/trunk/SOURCE/write_3d_binary.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! +qswst, qswst_m, tswst, tswst_m 7 7 ! 8 8 ! Former revisions: … … 111 111 WRITE ( 14 ) 'qsws_m '; WRITE ( 14 ) qsws_m 112 112 ENDIF 113 WRITE ( 14 ) 'qswst '; WRITE ( 14 ) qswst 114 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 115 WRITE ( 14 ) 'qswst_m '; WRITE ( 14 ) qswst_m 116 ENDIF 113 117 ENDIF 114 118 IF ( ALLOCATED( ql_c_av ) ) THEN … … 140 144 IF ( ALLOCATED( ts_av ) ) THEN 141 145 WRITE ( 14 ) 'ts_av '; WRITE ( 14 ) ts_av 146 ENDIF 147 WRITE ( 14 ) 'tswst '; WRITE ( 14 ) tswst 148 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 149 WRITE ( 14 ) 'tswst_m '; WRITE ( 14 ) tswst_m 142 150 ENDIF 143 151 WRITE ( 14 ) 'u '; WRITE ( 14 ) u -
palm/trunk/SOURCE/write_var_list.f90
r4 r19 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! +top_heatflux 7 7 ! 8 8 ! Former revisions: … … 306 306 WRITE ( 14 ) 'topography ' 307 307 WRITE ( 14 ) topography 308 WRITE ( 14 ) 'top_heatflux ' 309 WRITE ( 14 ) top_heatflux 308 310 WRITE ( 14 ) 'tsc ' 309 311 WRITE ( 14 ) tsc
Note: See TracChangeset
for help on using the changeset viewer.
