Changeset 1001
- Timestamp:
- Sep 13, 2012 2:08:46 PM (12 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 8 deleted
- 24 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r873 r1001 4 4 # Current revisions: 5 5 # ------------------ 6 # 6 # -asselin_filter, advec_s|u|v|w_ups, spline_x|y|z 7 7 # 8 8 # Former revisions: … … 86 86 PROG = palm 87 87 88 RCS = advec_s_bc.f90 advec_s_pw.f90 advec_s_up.f90 \ 89 advec_ws.f90 advec_s_ups.f90 advec_u_pw.f90 advec_u_up.f90 \ 90 advec_u_ups.f90 advec_v_pw.f90 advec_v_up.f90 advec_v_ups.f90 \ 91 advec_w_pw.f90 advec_w_up.f90 advec_w_ups.f90 asselin_filter.f90 \ 92 average_3d_data.f90 boundary_conds.f90 buoyancy.f90 \ 93 calc_liquid_water_content.f90 calc_precipitation.f90 \ 88 RCS = advec_s_bc.f90 advec_s_pw.f90 advec_s_up.f90 advec_ws.f90 \ 89 advec_u_pw.f90 advec_u_up.f90 advec_v_pw.f90 advec_v_up.f90 \ 90 advec_w_pw.f90 advec_w_up.f90 average_3d_data.f90 boundary_conds.f90 \ 91 buoyancy.f90 calc_liquid_water_content.f90 calc_precipitation.f90 \ 94 92 calc_radiation.f90 calc_spectra.f90 check_for_restart.f90 \ 95 93 check_open.f90 check_parameters.f90 close_file.f90 compute_vpt.f90 \ … … 123 121 production_e.f90 prognostic_equations.f90 random_function.f90 \ 124 122 random_gauss.f90 read_3d_binary.f90 read_var_list.f90 run_control.f90 \ 125 set_slicer_attributes_dvrp.f90 singleton.f90 sor.f90 spline_x.f90\126 s pline_y.f90 spline_z.f90 subsidence.f90 sum_up_3d_data.f90 \123 set_slicer_attributes_dvrp.f90 singleton.f90 sor.f90 \ 124 subsidence.f90 sum_up_3d_data.f90 \ 127 125 surface_coupler.f90 swap_timelevel.f90 temperton_fft.f90 \ 128 126 time_integration.f90 time_to_string.f90 timestep.f90 \ … … 141 139 write_3d_binary.f90 write_compressed.f90 write_var_list.f90 142 140 143 OBJS = advec_s_bc.o advec_s_pw.o advec_s_up.o advec_s_ups.o advec_u_pw.o \ 144 advec_u_up.o advec_u_ups.o advec_ws.o advec_v_pw.o advec_v_up.o \ 145 advec_v_ups.o advec_w_pw.o advec_w_up.o advec_w_ups.o asselin_filter.o \ 141 OBJS = advec_s_bc.o advec_s_pw.o advec_s_up.o advec_u_pw.o advec_u_up.o \ 142 advec_ws.o advec_v_pw.o advec_v_up.o advec_w_pw.o advec_w_up.o \ 146 143 average_3d_data.o boundary_conds.o buoyancy.o \ 147 144 calc_liquid_water_content.o calc_precipitation.o calc_radiation.o \ … … 173 170 production_e.o prognostic_equations.o random_function.o random_gauss.o \ 174 171 read_3d_binary.o read_var_list.o run_control.o \ 175 set_slicer_attributes_dvrp.o singleton.o sor.o spline_x.o spline_y.o\176 s pline_z.o subsidence.o sum_up_3d_data.o surface_coupler.o \172 set_slicer_attributes_dvrp.o singleton.o sor.o \ 173 subsidence.o sum_up_3d_data.o surface_coupler.o \ 177 174 swap_timelevel.o temperton_fft.o time_integration.o time_to_string.o \ 178 175 timestep.o timestep_scheme_steering.o transpose.o \ … … 219 216 advec_s_pw.o: modules.o 220 217 advec_s_up.o: modules.o 221 advec_s_ups.o: modules.o222 218 advec_u_pw.o: modules.o 223 219 advec_u_up.o: modules.o 224 advec_u_ups.o: modules.o225 220 advec_v_pw.o: modules.o 226 221 advec_v_up.o: modules.o 227 advec_v_ups.o: modules.o228 222 advec_ws.o: modules.o 229 223 advec_w_pw.o: modules.o 230 224 advec_w_up.o: modules.o 231 advec_w_ups.o: modules.o232 asselin_filter.o: modules.o233 225 average_3d_data.o: modules.o 234 226 boundary_conds.o: modules.o … … 340 332 singleton.o: singleton.f90 341 333 sor.o: modules.o 342 spline_x.o: modules.o343 spline_y.o: modules.o344 spline_z.o: modules.o345 334 subsidence.o: modules.o 346 335 sum_up_3d_data.o: modules.o -
palm/trunk/SOURCE/check_parameters.f90
r997 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! all actions concerning leapfrog- and upstream-spline-scheme removed 7 7 ! 8 8 ! Former revisions: … … 12 12 ! 996 2012-09-07 10:41:47Z raasch 13 13 ! little reformatting 14 ! 14 15 15 ! 978 2012-08-09 08:28:32Z fricke 16 16 ! setting of bc_lr/ns_dirneu/neudir … … 492 492 THEN 493 493 WRITE( action, '(A,A)' ) 'momentum_advec = ', momentum_advec 494 ENDIF495 IF ( timestep_scheme(1:8) == 'leapfrog' ) THEN496 WRITE( action, '(A,A)' ) 'timestep_scheme = ', timestep_scheme497 494 ENDIF 498 495 IF ( psolver == 'sor' ) THEN … … 562 559 563 560 action = ' ' 564 IF ( timestep_scheme(1:8) == 'leapfrog' ) THEN565 WRITE( action, '(A,A)' ) 'timestep_scheme = ', timestep_scheme566 ENDIF567 IF ( momentum_advec == 'ups-scheme' ) THEN568 WRITE( action, '(A,A)' ) 'momentum_advec = ', momentum_advec569 ENDIF570 561 IF ( action /= ' ' ) THEN 571 562 message_string = 'ocean = .T. does not allow ' // TRIM( action ) … … 653 644 IF ( scalar_advec == 'ws-scheme' ) ws_scheme_sca = .TRUE. 654 645 655 IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' .AND.&656 momentum_advec /= 'ups-scheme' )THEN646 IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' ) & 647 THEN 657 648 message_string = 'unknown advection scheme: momentum_advec = "' // & 658 649 TRIM( momentum_advec ) // '"' 659 650 CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 ) 660 651 ENDIF 661 IF ((( momentum_advec == 'ups-scheme' .OR. scalar_advec == 'ups-scheme' )& 662 .AND. timestep_scheme /= 'euler' ) .OR. (( momentum_advec == 'ws-scheme'& 663 .OR. scalar_advec == 'ws-scheme') .AND. (timestep_scheme == 'euler' .OR. & 664 timestep_scheme == 'leapfrog+euler' .OR. timestep_scheme == 'leapfrog' & 665 .OR. timestep_scheme == 'runge-kutta-2'))) THEN 652 IF ( ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' ) & 653 .AND. ( timestep_scheme == 'euler' .OR. & 654 timestep_scheme == 'runge-kutta-2' ) ) & 655 THEN 666 656 message_string = 'momentum_advec or scalar_advec = "' & 667 657 // TRIM( momentum_advec ) // '" is not allowed with timestep_scheme = "' // & … … 670 660 ENDIF 671 661 IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme' .AND. & 672 scalar_advec /= 'bc-scheme' .AND. scalar_advec /= 'ups-scheme' ) THEN 662 scalar_advec /= 'bc-scheme' ) & 663 THEN 673 664 message_string = 'unknown advection scheme: scalar_advec = "' // & 674 665 TRIM( scalar_advec ) // '"' … … 683 674 ENDIF 684 675 685 IF ( use_upstream_for_tke .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN686 message_string = 'use_upstream_for_tke = .TRUE. not allowed with ' // &687 'timestep_scheme = "' // TRIM( timestep_scheme ) // '"'688 CALL message( 'check_parameters', 'PA0026', 1, 2, 0, 6, 0 )689 ENDIF690 691 676 IF ( use_sgs_for_particles .AND. curvature_solution_effects ) THEN 692 677 message_string = 'use_sgs_for_particles = .TRUE. not allowed with ' // & … … 701 686 CASE ( 'euler' ) 702 687 intermediate_timestep_count_max = 1 703 asselin_filter_factor = 0.0704 705 CASE ( 'leapfrog', 'leapfrog+euler' )706 intermediate_timestep_count_max = 1707 688 708 689 CASE ( 'runge-kutta-2' ) 709 690 intermediate_timestep_count_max = 2 710 asselin_filter_factor = 0.0711 691 712 692 CASE ( 'runge-kutta-3' ) 713 693 intermediate_timestep_count_max = 3 714 asselin_filter_factor = 0.0715 694 716 695 CASE DEFAULT … … 720 699 721 700 END SELECT 722 723 IF ( scalar_advec == 'ups-scheme' .AND. timestep_scheme(1:5) == 'runge' )&724 THEN725 message_string = 'scalar advection scheme "' // TRIM( scalar_advec ) // &726 '" & does not work with timestep_scheme "' // &727 TRIM( timestep_scheme ) // '"'728 CALL message( 'check_parameters', 'PA0028', 1, 2, 0, 6, 0 )729 ENDIF730 701 731 702 IF ( (momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme') & … … 831 802 ENDIF 832 803 833 IF ( humidity .AND. scalar_advec == 'ups-scheme' ) THEN834 message_string = 'UPS-scheme is not implemented for humidity = .TRUE.'835 CALL message( 'check_parameters', 'PA0037', 1, 2, 0, 6, 0 )836 ENDIF837 838 804 IF ( passive_scalar .AND. humidity ) THEN 839 805 message_string = 'humidity = .TRUE. and passive_scalar = .TRUE. ' // & 840 806 'is not allowed simultaneously' 841 807 CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 ) 842 ENDIF843 844 IF ( passive_scalar .AND. scalar_advec == 'ups-scheme' ) THEN845 message_string = 'UPS-scheme is not implemented for passive_scalar' // &846 ' = .TRUE.'847 CALL message( 'check_parameters', 'PA0039', 1, 2, 0, 6, 0 )848 808 ENDIF 849 809 … … 1319 1279 IF ( cfl_factor <= 0.0 .OR. cfl_factor > 1.0 ) THEN 1320 1280 IF ( cfl_factor == -1.0 ) THEN 1321 IF ( momentum_advec == 'ups-scheme' .OR. & 1322 scalar_advec == 'ups-scheme' ) THEN 1281 IF ( timestep_scheme == 'runge-kutta-2' ) THEN 1323 1282 cfl_factor = 0.8 1283 ELSEIF ( timestep_scheme == 'runge-kutta-3' ) THEN 1284 cfl_factor = 0.9 1324 1285 ELSE 1325 IF ( timestep_scheme == 'runge-kutta-2' ) THEN 1326 cfl_factor = 0.8 1327 ELSEIF ( timestep_scheme == 'runge-kutta-3' ) THEN 1328 cfl_factor = 0.9 1329 ELSE 1330 cfl_factor = 0.1 1331 ENDIF 1286 cfl_factor = 0.9 1332 1287 ENDIF 1333 1288 ELSE … … 1717 1672 !-- Compute and check, respectively, the Rayleigh Damping parameter 1718 1673 IF ( rayleigh_damping_factor == -1.0 ) THEN 1719 IF ( momentum_advec == 'ups-scheme' ) THEN 1720 rayleigh_damping_factor = 0.01 1721 ELSE 1722 rayleigh_damping_factor = 0.0 1723 ENDIF 1674 rayleigh_damping_factor = 0.0 1724 1675 ELSE 1725 1676 IF ( rayleigh_damping_factor < 0.0 .OR. rayleigh_damping_factor > 1.0 ) & … … 1753 1704 ENDIF 1754 1705 ENDIF 1755 ENDIF1756 1757 !1758 !-- Check limiters for Upstream-Spline scheme1759 IF ( overshoot_limit_u < 0.0 .OR. overshoot_limit_v < 0.0 .OR. &1760 overshoot_limit_w < 0.0 .OR. overshoot_limit_pt < 0.0 .OR. &1761 overshoot_limit_e < 0.0 ) THEN1762 message_string = 'overshoot_limit_... < 0.0 is not allowed'1763 CALL message( 'check_parameters', 'PA0080', 1, 2, 0, 6, 0 )1764 ENDIF1765 IF ( ups_limit_u < 0.0 .OR. ups_limit_v < 0.0 .OR. ups_limit_w < 0.0 .OR. &1766 ups_limit_pt < 0.0 .OR. ups_limit_e < 0.0 ) THEN1767 message_string = 'ups_limit_... < 0.0 is not allowed'1768 CALL message( 'check_parameters', 'PA0081', 1, 2, 0, 6, 0 )1769 1706 ENDIF 1770 1707 -
palm/trunk/SOURCE/diffusion_e.f90
r826 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! most arrays comunicated by module instead of parameter list 7 7 ! 8 8 ! Former revisions: … … 64 64 ! Call for all grid points 65 65 !------------------------------------------------------------------------------! 66 SUBROUTINE diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, var, &67 var_reference, rif, tend, zu, zw ) 68 66 SUBROUTINE diffusion_e( var, var_reference ) 67 68 USE arrays_3d 69 69 USE control_parameters 70 70 USE grid_variables … … 75 75 76 76 INTEGER :: i, j, k 77 REAL :: dvar_dz, l_stable, phi_m, var_reference 78 REAL :: ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), & 79 l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1) 80 REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: diss, tend 81 REAL, DIMENSION(:,:), POINTER :: rif 82 REAL, DIMENSION(:,:,:), POINTER :: e, km, var 77 REAL :: dvar_dz, l_stable, phi_m, var_reference 78 79 REAL, DIMENSION(:,:,:), POINTER :: var 83 80 REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation, l, ll 84 81 … … 285 282 ! Call for grid point i,j 286 283 !------------------------------------------------------------------------------! 287 SUBROUTINE diffusion_e_ij( i, j, ddzu, dd2zu, ddzw, diss, e, km, l_grid, &288 var, var_reference, rif, tend, zu, zw ) 289 284 SUBROUTINE diffusion_e_ij( i, j, var, var_reference ) 285 286 USE arrays_3d 290 287 USE control_parameters 291 288 USE grid_variables … … 295 292 IMPLICIT NONE 296 293 297 INTEGER :: i, j, k 298 REAL :: dvar_dz, l_stable, phi_m, var_reference 299 REAL :: ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), & 300 l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1) 301 REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: diss, tend 302 REAL, DIMENSION(:,:), POINTER :: rif 303 REAL, DIMENSION(:,:,:), POINTER :: e, km, var 304 REAL, DIMENSION(nzb+1:nzt) :: dissipation, l, ll 294 INTEGER :: i, j, k 295 REAL :: dvar_dz, l_stable, phi_m, var_reference 296 297 REAL, DIMENSION(:,:,:), POINTER :: var 298 REAL, DIMENSION(nzb+1:nzt) :: dissipation, l, ll 305 299 306 300 -
palm/trunk/SOURCE/diffusion_s.f90
r668 r1001 3 3 !------------------------------------------------------------------------------! 4 4 ! Current revisions: 5 ! ----------------- 5 ! ------------------ 6 ! some arrays comunicated by module instead of parameter list 6 7 ! 7 8 ! Former revisions: … … 54 55 ! Call for all grid points 55 56 !------------------------------------------------------------------------------! 56 SUBROUTINE diffusion_s( ddzu, ddzw, kh, s, s_flux_b, s_flux_t, &57 wall_s_flux, tend ) 58 57 SUBROUTINE diffusion_s( s, s_flux_b, s_flux_t, wall_s_flux ) 58 59 USE arrays_3d 59 60 USE control_parameters 60 61 USE grid_variables … … 65 66 INTEGER :: i, j, k 66 67 REAL :: vertical_gridspace 67 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1)68 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)69 68 REAL :: wall_s_flux(0:4) 70 REAL, DIMENSION( :,:), POINTER:: s_flux_b, s_flux_t71 REAL, DIMENSION(:,:,:), POINTER :: kh,s69 REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: s_flux_b, s_flux_t 70 REAL, DIMENSION(:,:,:), POINTER :: s 72 71 73 72 DO i = nxl, nxr … … 166 165 ! Call for grid point i,j 167 166 !------------------------------------------------------------------------------! 168 SUBROUTINE diffusion_s_ij( i, j, ddzu, ddzw, kh, s, s_flux_b, s_flux_t, &169 wall_s_flux, tend ) 170 167 SUBROUTINE diffusion_s_ij( i, j, s, s_flux_b, s_flux_t, wall_s_flux ) 168 169 USE arrays_3d 171 170 USE control_parameters 172 171 USE grid_variables … … 177 176 INTEGER :: i, j, k 178 177 REAL :: vertical_gridspace 179 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1)180 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)181 178 REAL :: wall_s_flux(0:4) 182 REAL, DIMENSION( :,:), POINTER:: s_flux_b, s_flux_t183 REAL, DIMENSION(:,:,:), POINTER :: kh,s179 REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: s_flux_b, s_flux_t 180 REAL, DIMENSION(:,:,:), POINTER :: s 184 181 185 182 ! -
palm/trunk/SOURCE/diffusion_u.f90
r979 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! arrays comunicated by module instead of parameter list 7 7 ! 8 8 ! Former revisions: … … 68 68 ! Call for all grid points 69 69 !------------------------------------------------------------------------------! 70 SUBROUTINE diffusion_u( ddzu, ddzw, km, tend, u, usws, uswst, v, w ) 71 70 SUBROUTINE diffusion_u 71 72 USE arrays_3d 72 73 USE control_parameters 73 74 USE grid_variables … … 78 79 INTEGER :: i, j, k 79 80 REAL :: kmym, kmyp, kmzm, kmzp 80 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 81 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 82 REAL, DIMENSION(:,:), POINTER :: usws, uswst 83 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w 81 84 82 REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs 85 83 … … 93 91 94 92 DO i = nxlu, nxr 95 DO j = nys, nyn93 DO j = nys, nyn 96 94 ! 97 95 !-- Compute horizontal diffusion … … 221 219 ! Call for grid point i,j 222 220 !------------------------------------------------------------------------------! 223 SUBROUTINE diffusion_u_ij( i, j , ddzu, ddzw, km, tend, u, usws, &224 uswst, v, w ) 225 221 SUBROUTINE diffusion_u_ij( i, j ) 222 223 USE arrays_3d 226 224 USE control_parameters 227 225 USE grid_variables … … 232 230 INTEGER :: i, j, k 233 231 REAL :: kmym, kmyp, kmzm, kmzp 234 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 235 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 236 REAL, DIMENSION(nzb:nzt+1) :: usvs 237 REAL, DIMENSION(:,:), POINTER :: usws, uswst 238 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w 232 233 REAL, DIMENSION(nzb:nzt+1) :: usvs 239 234 240 235 ! -
palm/trunk/SOURCE/diffusion_v.f90
r979 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! arrays comunicated by module instead of parameter list 7 7 ! 8 8 ! Former revisions: … … 66 66 ! Call for all grid points 67 67 !------------------------------------------------------------------------------! 68 SUBROUTINE diffusion_v( ddzu, ddzw, km, tend, u, v, vsws, vswst, w ) 69 68 SUBROUTINE diffusion_v 69 70 USE arrays_3d 70 71 USE control_parameters 71 72 USE grid_variables … … 76 77 INTEGER :: i, j, k 77 78 REAL :: kmxm, kmxp, kmzm, kmzp 78 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 79 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 80 REAL, DIMENSION(:,:), POINTER :: vsws, vswst 81 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w 79 82 80 REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: vsus 83 81 … … 219 217 ! Call for grid point i,j 220 218 !------------------------------------------------------------------------------! 221 SUBROUTINE diffusion_v_ij( i, j , ddzu, ddzw, km, tend, u, v, &222 vsws, vswst, w ) 223 219 SUBROUTINE diffusion_v_ij( i, j ) 220 221 USE arrays_3d 224 222 USE control_parameters 225 223 USE grid_variables … … 230 228 INTEGER :: i, j, k 231 229 REAL :: kmxm, kmxp, kmzm, kmzp 232 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 233 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 234 REAL, DIMENSION(nzb:nzt+1) :: vsus 235 REAL, DIMENSION(:,:), POINTER :: vsws, vswst 236 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w 230 231 REAL, DIMENSION(nzb:nzt+1) :: vsus 237 232 238 233 ! -
palm/trunk/SOURCE/diffusion_w.f90
r979 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! arrays comunicated by module instead of parameter list 7 7 ! 8 8 ! Former revisions: … … 60 60 ! Call for all grid points 61 61 !------------------------------------------------------------------------------! 62 SUBROUTINE diffusion_w( ddzu, ddzw, km, tend, u, v, w ) 63 62 SUBROUTINE diffusion_w 63 64 USE arrays_3d 64 65 USE control_parameters 65 66 USE grid_variables … … 70 71 INTEGER :: i, j, k 71 72 REAL :: kmxm, kmxp, kmym, kmyp 72 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 73 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 74 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w 73 75 74 REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsus, wsvs 76 75 … … 170 169 ! Call for grid point i,j 171 170 !------------------------------------------------------------------------------! 172 SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, tend, u, v, w ) 173 171 SUBROUTINE diffusion_w_ij( i, j ) 172 173 USE arrays_3d 174 174 USE control_parameters 175 175 USE grid_variables … … 180 180 INTEGER :: i, j, k 181 181 REAL :: kmxm, kmxp, kmym, kmyp 182 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 183 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 184 REAL, DIMENSION(nzb:nzt+1) :: wsus, wsvs 185 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w 182 183 REAL, DIMENSION(nzb:nzt+1) :: wsus, wsvs 186 184 187 185 -
palm/trunk/SOURCE/header.f90
r979 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! all actions concerning leapfrog- and upstream-spline-scheme removed 7 7 ! 8 8 ! Former revisions: … … 296 296 ELSEIF (momentum_advec == 'ws-scheme' ) THEN 297 297 WRITE ( io, 503 ) 298 ELSEIF (momentum_advec == 'ups-scheme' ) THEN299 WRITE ( io, 114 )300 IF ( cut_spline_overshoot ) WRITE ( io, 124 )301 IF ( overshoot_limit_u /= 0.0 .OR. overshoot_limit_v /= 0.0 .OR. &302 overshoot_limit_w /= 0.0 ) THEN303 WRITE ( io, 127 ) overshoot_limit_u, overshoot_limit_v, &304 overshoot_limit_w305 ENDIF306 IF ( ups_limit_u /= 0.0 .OR. ups_limit_v /= 0.0 .OR. &307 ups_limit_w /= 0.0 ) &308 THEN309 WRITE ( io, 125 ) ups_limit_u, ups_limit_v, ups_limit_w310 ENDIF311 IF ( long_filter_factor /= 0.0 ) WRITE ( io, 115 ) long_filter_factor312 298 ENDIF 313 299 IF ( scalar_advec == 'pw-scheme' ) THEN … … 315 301 ELSEIF ( scalar_advec == 'ws-scheme' ) THEN 316 302 WRITE ( io, 504 ) 317 ELSEIF ( scalar_advec == 'ups-scheme' ) THEN318 WRITE ( io, 117 )319 IF ( cut_spline_overshoot ) WRITE ( io, 124 )320 IF ( overshoot_limit_e /= 0.0 .OR. overshoot_limit_pt /= 0.0 ) THEN321 WRITE ( io, 128 ) overshoot_limit_e, overshoot_limit_pt322 ENDIF323 IF ( ups_limit_e /= 0.0 .OR. ups_limit_pt /= 0.0 ) THEN324 WRITE ( io, 126 ) ups_limit_e, ups_limit_pt325 ENDIF326 303 ELSE 327 304 WRITE ( io, 118 ) … … 344 321 advected_distance_x/1000.0, advected_distance_y/1000.0 345 322 ENDIF 346 IF ( timestep_scheme == 'leapfrog' ) THEN 347 WRITE ( io, 120 ) 348 ELSEIF ( timestep_scheme == 'leapfrog+euler' ) THEN 349 WRITE ( io, 121 ) 350 ELSE 351 WRITE ( io, 122 ) timestep_scheme 352 ENDIF 323 WRITE ( io, 122 ) timestep_scheme 353 324 IF ( use_upstream_for_tke ) WRITE ( io, 143 ) 354 325 IF ( rayleigh_damping_factor /= 0.0 ) THEN … … 1646 1617 113 FORMAT (' --> Momentum advection via Piascek-Williams-Scheme (Form C3)', & 1647 1618 ' or Upstream') 1648 114 FORMAT (' --> Momentum advection via Upstream-Spline-Scheme')1649 115 FORMAT (' Tendencies are smoothed via Long-Filter with factor ',F5.3)1650 1619 116 FORMAT (' --> Scalar advection via Piascek-Williams-Scheme (Form C3)', & 1651 1620 ' or Upstream') 1652 117 FORMAT (' --> Scalar advection via Upstream-Spline-Scheme')1653 1621 118 FORMAT (' --> Scalar advection via Bott-Chlond-Scheme') 1654 1622 119 FORMAT (' --> Galilei-Transform applied to horizontal advection', & 1655 1623 ' Translation velocity = ',A/ & 1656 1624 ' distance advected ',A,': ',F8.3,' km(x) ',F8.3,' km(y)') 1657 120 FORMAT (' --> Time differencing scheme: leapfrog only (no euler in case', &1658 ' of timestep changes)')1659 121 FORMAT (' --> Time differencing scheme: leapfrog + euler in case of', &1660 ' timestep changes')1661 1625 122 FORMAT (' --> Time differencing scheme: ',A) 1662 1626 123 FORMAT (' --> Rayleigh-Damping active, starts ',A,' z = ',F8.2,' m'/ & 1663 1627 ' maximum damping coefficient: ',F5.3, ' 1/s') 1664 124 FORMAT (' Spline-overshoots are being suppressed')1665 125 FORMAT (' Upstream-Scheme is used if Upstream-differences fall short', &1666 ' of'/ &1667 ' delta_u = ',F6.4,4X,'delta_v = ',F6.4,4X,'delta_w = ',F6.4)1668 126 FORMAT (' Upstream-Scheme is used if Upstream-differences fall short', &1669 ' of'/ &1670 ' delta_e = ',F6.4,4X,'delta_pt = ',F6.4)1671 127 FORMAT (' The following absolute overshoot differences are tolerated:'/&1672 ' delta_u = ',F6.4,4X,'delta_v = ',F6.4,4X,'delta_w = ',F6.4)1673 128 FORMAT (' The following absolute overshoot differences are tolerated:'/&1674 ' delta_e = ',F6.4,4X,'delta_pt = ',F6.4)1675 1628 129 FORMAT (' --> Additional prognostic equation for the specific humidity') 1676 1629 130 FORMAT (' --> Additional prognostic equation for the total water content') -
palm/trunk/SOURCE/init_1d_model.f90
r997 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! all actions concerning leapfrog scheme removed 7 7 ! 8 8 ! Former revisions: … … 68 68 ! 69 69 !-- Allocate required 1D-arrays 70 ALLOCATE( e1d(nzb:nzt+1), e1d_ m(nzb:nzt+1), e1d_p(nzb:nzt+1), &71 kh1d(nzb:nzt+1), k h1d_m(nzb:nzt+1), km1d(nzb:nzt+1), &72 km1d_m(nzb:nzt+1),l_black(nzb:nzt+1), l1d(nzb:nzt+1), &73 l1d_m(nzb:nzt+1),rif1d(nzb:nzt+1), te_e(nzb:nzt+1), &70 ALLOCATE( e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), & 71 kh1d(nzb:nzt+1), km1d(nzb:nzt+1), & 72 l_black(nzb:nzt+1), l1d(nzb:nzt+1), & 73 rif1d(nzb:nzt+1), te_e(nzb:nzt+1), & 74 74 te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1), & 75 75 te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1), & 76 u1d_ m(nzb:nzt+1), u1d_p(nzb:nzt+1),v1d(nzb:nzt+1), &77 v1d_ m(nzb:nzt+1), v1d_p(nzb:nzt+1) )76 u1d_p(nzb:nzt+1), v1d(nzb:nzt+1), & 77 v1d_p(nzb:nzt+1) ) 78 78 79 79 ! 80 80 !-- Initialize arrays 81 81 IF ( constant_diffusion ) THEN 82 km1d = km_constant 83 km1d_m = km_constant 84 kh1d = km_constant / prandtl_number 85 kh1d_m = km_constant / prandtl_number 82 km1d = km_constant 83 kh1d = km_constant / prandtl_number 86 84 ELSE 87 e1d = 0.0; e1d_ m = 0.0; e1d_p = 0.088 kh1d = 0.0; k h1d_m = 0.0; km1d = 0.0; km1d_m= 0.085 e1d = 0.0; e1d_p = 0.0 86 kh1d = 0.0; km1d = 0.0 89 87 rif1d = 0.0 90 88 ! … … 123 121 ENDIF 124 122 l1d = l_black 125 l1d_m = l_black126 123 u1d = u_init 127 u1d_m = u_init128 124 u1d_p = u_init 129 125 v1d = v_init 130 v1d_m = v_init131 126 v1d_p = v_init 132 127 … … 136 131 !-- in the initial phase of a run (at k=1, dz/2 occurs in the limiting formula!) 137 132 u1d(0:1) = 0.1 138 u1d_m(0:1) = 0.1139 133 u1d_p(0:1) = 0.1 140 134 v1d(0:1) = 0.1 141 v1d_m(0:1) = 0.1142 135 v1d_p(0:1) = 0.1 143 136 … … 152 145 ENDIF 153 146 ts1d = 0.0 154 usws1d = 0.0 ; usws1d_m = 0.0155 vsws1d = 0.0 ; vsws1d_m = 0.0147 usws1d = 0.0 148 vsws1d = 0.0 156 149 z01d = roughness_length 157 150 z0h1d = z0h_factor * z01d … … 226 219 DO k = nzb_diff, nzt 227 220 228 kmzm = 0.5 * ( km1d _m(k-1) + km1d_m(k) )229 kmzp = 0.5 * ( km1d _m(k) + km1d_m(k+1) )221 kmzm = 0.5 * ( km1d(k-1) + km1d(k) ) 222 kmzp = 0.5 * ( km1d(k) + km1d(k+1) ) 230 223 ! 231 224 !-- u-component 232 225 te_u(k) = f * ( v1d(k) - vg(k) ) + ( & 233 kmzp * ( u1d _m(k+1) - u1d_m(k) ) * ddzu(k+1) &234 - kmzm * ( u1d _m(k) - u1d_m(k-1) ) * ddzu(k) &235 ) * ddzw(k)226 kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) & 227 - kmzm * ( u1d(k) - u1d(k-1) ) * ddzu(k) & 228 ) * ddzw(k) 236 229 ! 237 230 !-- v-component 238 te_v(k) = -f * ( u1d(k) - ug(k) ) + ( &239 kmzp * ( v1d _m(k+1) - v1d_m(k) ) * ddzu(k+1) &240 - kmzm * ( v1d _m(k) - v1d_m(k-1) ) * ddzu(k) &241 ) * ddzw(k)231 te_v(k) = -f * ( u1d(k) - ug(k) ) + ( & 232 kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) & 233 - kmzm * ( v1d(k) - v1d(k-1) ) * ddzu(k) & 234 ) * ddzw(k) 242 235 ENDDO 243 236 IF ( .NOT. constant_diffusion ) THEN … … 245 238 ! 246 239 !-- TKE 247 kmzm = 0.5 * ( km1d _m(k-1) + km1d_m(k) )248 kmzp = 0.5 * ( km1d _m(k) + km1d_m(k+1) )240 kmzm = 0.5 * ( km1d(k-1) + km1d(k) ) 241 kmzp = 0.5 * ( km1d(k) + km1d(k+1) ) 249 242 IF ( .NOT. humidity ) THEN 250 243 pt_0 = pt_init(k) … … 260 253 ! 261 254 !-- According to Detering, c_e=0.064 262 dissipation = 0.064 * e1d _m(k) * SQRT( e1d_m(k) ) / l1d_m(k)255 dissipation = 0.064 * e1d(k) * SQRT( e1d(k) ) / l1d(k) 263 256 ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN 264 dissipation = ( 0.19 + 0.74 * l1d _m(k) / l_grid(k) ) &265 * e1d _m(k) * SQRT( e1d_m(k) ) / l1d_m(k)257 dissipation = ( 0.19 + 0.74 * l1d(k) / l_grid(k) ) & 258 * e1d(k) * SQRT( e1d(k) ) / l1d(k) 266 259 ENDIF 267 260 … … 271 264 - g / pt_0 * kh1d(k) * flux & 272 265 + ( & 273 kmzp * ( e1d_m(k+1) - e1d_m(k) ) * ddzu(k+1) &274 - kmzm * ( e1d_m(k) - e1d_m(k-1) ) * ddzu(k) &266 kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1) & 267 - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k) & 275 268 ) * ddzw(k) & 276 - dissipation269 - dissipation 277 270 ENDDO 278 271 ENDIF … … 285 278 286 279 k = nzb+1 287 kmzm = 0.5 * ( km1d _m(k-1) + km1d_m(k) )288 kmzp = 0.5 * ( km1d _m(k) + km1d_m(k+1) )280 kmzm = 0.5 * ( km1d(k-1) + km1d(k) ) 281 kmzp = 0.5 * ( km1d(k) + km1d(k+1) ) 289 282 IF ( .NOT. humidity ) THEN 290 283 pt_0 = pt_init(k) … … 300 293 ! 301 294 !-- According to Detering, c_e=0.064 302 dissipation = 0.064 * e1d _m(k) * SQRT( e1d_m(k) ) / l1d_m(k)295 dissipation = 0.064 * e1d(k) * SQRT( e1d(k) ) / l1d(k) 303 296 ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN 304 dissipation = ( 0.19 + 0.74 * l1d _m(k) / l_grid(k) ) &305 * e1d _m(k) * SQRT( e1d_m(k) ) / l1d_m(k)297 dissipation = ( 0.19 + 0.74 * l1d(k) / l_grid(k) ) & 298 * e1d(k) * SQRT( e1d(k) ) / l1d(k) 306 299 ENDIF 307 300 308 301 ! 309 302 !-- u-component 310 te_u(k) = f * ( v1d(k) - vg(k) ) + ( &311 kmzp * ( u1d _m(k+1) - u1d_m(k) ) * ddzu(k+1) + usws1d_m&312 ) * 2.0 * ddzw(k)303 te_u(k) = f * ( v1d(k) - vg(k) ) + ( & 304 kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) + usws1d & 305 ) * 2.0 * ddzw(k) 313 306 ! 314 307 !-- v-component 315 te_v(k) = -f * ( u1d(k) - ug(k) ) + ( &316 kmzp * ( v1d _m(k+1) - v1d_m(k) ) * ddzu(k+1) + vsws1d_m&317 ) * 2.0 * ddzw(k)308 te_v(k) = -f * ( u1d(k) - ug(k) ) + ( & 309 kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) + vsws1d & 310 ) * 2.0 * ddzw(k) 318 311 ! 319 312 !-- TKE … … 323 316 - g / pt_0 * kh1d(k) * flux & 324 317 + ( & 325 kmzp * ( e1d_m(k+1) - e1d_m(k) ) * ddzu(k+1) &326 - kmzm * ( e1d_m(k) - e1d_m(k-1) ) * ddzu(k) &318 kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1) & 319 - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k) & 327 320 ) * ddzw(k) & 328 - dissipation321 - dissipation 329 322 ENDIF 330 323 … … 333 326 DO k = nzb+1, nzt 334 327 335 u1d_p(k) = ( 1. - tsc(1) ) * u1d_m(k) + & 336 tsc(1) * u1d(k) + dt_1d * ( tsc(2) * te_u(k) + & 337 tsc(3) * te_um(k) ) 338 v1d_p(k) = ( 1. - tsc(1) ) * v1d_m(k) + & 339 tsc(1) * v1d(k) + dt_1d * ( tsc(2) * te_v(k) + & 340 tsc(3) * te_vm(k) ) 328 u1d_p(k) = u1d(k) + dt_1d * ( tsc(2) * te_u(k) + & 329 tsc(3) * te_um(k) ) 330 v1d_p(k) = v1d(k) + dt_1d * ( tsc(2) * te_v(k) + & 331 tsc(3) * te_vm(k) ) 341 332 342 333 ENDDO … … 344 335 DO k = nzb+1, nzt 345 336 346 e1d_p(k) = ( 1. - tsc(1) ) * e1d_m(k) + & 347 tsc(1) * e1d(k) + dt_1d * ( tsc(2) * te_e(k) + & 348 tsc(3) * te_em(k) ) 337 e1d_p(k) = e1d(k) + dt_1d * ( tsc(2) * te_e(k) + & 338 tsc(3) * te_em(k) ) 349 339 350 340 ENDDO … … 403 393 404 394 ! 405 !-- If necessary, apply the time filter406 IF ( asselin_filter_factor /= 0.0 .AND. &407 timestep_scheme(1:5) /= 'runge' ) THEN408 409 u1d = u1d + asselin_filter_factor * ( u1d_p - 2.0 * u1d + u1d_m )410 v1d = v1d + asselin_filter_factor * ( v1d_p - 2.0 * v1d + v1d_m )411 412 IF ( .NOT. constant_diffusion ) THEN413 e1d = e1d + asselin_filter_factor * &414 ( e1d_p - 2.0 * e1d + e1d_m )415 ENDIF416 417 ENDIF418 419 !420 395 !-- Swap the time levels in preparation for the next time step. 421 IF ( timestep_scheme(1:4) == 'leap' ) THEN422 u1d_m = u1d423 v1d_m = v1d424 IF ( .NOT. constant_diffusion ) THEN425 e1d_m = e1d426 kh1d_m = kh1d ! The old diffusion quantities are required for427 km1d_m = km1d ! explicit diffusion in the leap-frog scheme.428 l1d_m = l1d429 IF ( prandtl_layer ) THEN430 usws1d_m = usws1d431 vsws1d_m = vsws1d432 ENDIF433 ENDIF434 ENDIF435 396 u1d = u1d_p 436 397 v1d = v1d_p … … 696 657 ENDIF ! .NOT. constant_diffusion 697 658 698 !699 !-- The Runge-Kutta scheme needs the recent diffusion quantities700 IF ( timestep_scheme(1:5) == 'runge' ) THEN701 u1d_m = u1d702 v1d_m = v1d703 IF ( .NOT. constant_diffusion ) THEN704 e1d_m = e1d705 kh1d_m = kh1d706 km1d_m = km1d707 l1d_m = l1d708 IF ( prandtl_layer ) THEN709 usws1d_m = usws1d710 vsws1d_m = vsws1d711 ENDIF712 ENDIF713 ENDIF714 715 716 659 ENDDO ! intermediate step loop 717 660 … … 836 779 837 780 INTEGER :: k 838 REAL :: div, dt_diff, fac, percent_change,value781 REAL :: div, dt_diff, fac, value 839 782 840 783 … … 842 785 !-- Compute the currently feasible time step according to the diffusion 843 786 !-- criterion. At nzb+1 the half grid length is used. 844 IF ( timestep_scheme(1:4) == 'leap' ) THEN 845 fac = 0.25 846 ELSE 847 fac = 0.35 848 ENDIF 787 fac = 0.35 849 788 dt_diff = dt_max_1d 850 789 DO k = nzb+2, nzt … … 866 805 ENDIF 867 806 868 IF ( timestep_scheme(1:4) == 'leap' ) THEN 869 870 ! 871 !-- The current time step will only be changed if the new time step exceeds 872 !-- its previous value by 5 % or falls 2 % below. After a time step 873 !-- reduction at least 30 iterations must be done with this value before a 874 !-- new reduction will be allowed again. 875 !-- The control parameters for application of Euler- or leap-frog schemes are 876 !-- set accordingly. 877 percent_change = dt_1d / old_dt_1d - 1.0 878 IF ( percent_change > 0.05 .OR. percent_change < -0.02 ) THEN 879 880 ! 881 !-- Each time step increase is by at most 2 % 882 IF ( percent_change > 0.0 .AND. simulated_time_1d /= 0.0 ) THEN 883 dt_1d = 1.02 * old_dt_1d 884 ENDIF 885 886 ! 887 !-- A more or less simple new time step value is obtained taking only the 888 !-- first two significant digits 889 div = 1000.0 890 DO WHILE ( dt_1d < div ) 891 div = div / 10.0 892 ENDDO 893 dt_1d = NINT( dt_1d * 100.0 / div ) * div / 100.0 894 895 ! 896 !-- Now the time step can be changed. 897 IF ( percent_change < 0.0 ) THEN 898 ! 899 !-- Time step reduction 900 old_dt_1d = dt_1d 901 last_dt_change_1d = current_timestep_number_1d 902 ELSE 903 ! 904 !-- Time step will only be increased if at least 30 iterations have 905 !-- been done since the previous time step change, and of course at 906 !-- simulation start, respectively. 907 IF ( current_timestep_number_1d >= last_dt_change_1d + 30 .OR. & 908 simulated_time_1d == 0.0 ) THEN 909 old_dt_1d = dt_1d 910 last_dt_change_1d = current_timestep_number_1d 911 ELSE 912 dt_1d = old_dt_1d 913 ENDIF 914 ENDIF 915 ELSE 916 ! 917 !-- No time step change since the difference is too small 918 dt_1d = old_dt_1d 919 ENDIF 920 921 ELSE ! Runge-Kutta 922 923 !-- A more or less simple new time step value is obtained taking only the 924 !-- first two significant digits 925 div = 1000.0 926 DO WHILE ( dt_1d < div ) 927 div = div / 10.0 928 ENDDO 929 dt_1d = NINT( dt_1d * 100.0 / div ) * div / 100.0 930 931 old_dt_1d = dt_1d 932 last_dt_change_1d = current_timestep_number_1d 933 934 ENDIF 807 ! 808 !-- A more or less simple new time step value is obtained taking only the 809 !-- first two significant digits 810 div = 1000.0 811 DO WHILE ( dt_1d < div ) 812 div = div / 10.0 813 ENDDO 814 dt_1d = NINT( dt_1d * 100.0 / div ) * div / 100.0 815 816 old_dt_1d = dt_1d 817 935 818 936 819 END SUBROUTINE timestep_1d -
palm/trunk/SOURCE/init_3d_model.f90
r997 r1001 7 7 ! Current revisions: 8 8 ! ------------------ 9 ! 9 ! all actions concerning leapfrog scheme removed 10 10 ! 11 11 ! Former revisions: … … 225 225 ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) ) 226 226 227 ALLOCATE( rif_1(nysg:nyng,nxlg:nxrg), shf_1(nysg:nyng,nxlg:nxrg), & 228 ts(nysg:nyng,nxlg:nxrg), tswst_1(nysg:nyng,nxlg:nxrg), & 229 us(nysg:nyng,nxlg:nxrg), usws_1(nysg:nyng,nxlg:nxrg), & 230 uswst_1(nysg:nyng,nxlg:nxrg), & 231 vsws_1(nysg:nyng,nxlg:nxrg), & 232 vswst_1(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg), & 227 ALLOCATE( rif(nysg:nyng,nxlg:nxrg), shf(nysg:nyng,nxlg:nxrg), & 228 ts(nysg:nyng,nxlg:nxrg), tswst(nysg:nyng,nxlg:nxrg), & 229 us(nysg:nyng,nxlg:nxrg), usws(nysg:nyng,nxlg:nxrg), & 230 uswst(nysg:nyng,nxlg:nxrg), vsws(nysg:nyng,nxlg:nxrg), & 231 vswst(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg), & 233 232 z0h(nysg:nyng,nxlg:nxrg) ) 234 235 IF ( timestep_scheme(1:5) /= 'runge' ) THEN236 !237 !-- Leapfrog scheme needs two timelevels of diffusion quantities238 ALLOCATE( rif_2(nysg:nyng,nxlg:nxrg), &239 shf_2(nysg:nyng,nxlg:nxrg), &240 tswst_2(nysg:nyng,nxlg:nxrg), &241 usws_2(nysg:nyng,nxlg:nxrg), &242 uswst_2(nysg:nyng,nxlg:nxrg), &243 vswst_2(nysg:nyng,nxlg:nxrg), &244 vsws_2(nysg:nyng,nxlg:nxrg) )245 ENDIF246 233 247 234 ALLOCATE( d(nzb+1:nzta,nys:nyna,nxl:nxra), & … … 249 236 e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 250 237 e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 251 kh _1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),&252 km _1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),&238 kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 239 km(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 253 240 p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 254 241 pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & … … 278 265 ENDIF 279 266 280 IF ( timestep_scheme(1:5) /= 'runge' ) THEN281 ALLOCATE( kh_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &282 km_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )283 ENDIF284 285 267 IF ( humidity .OR. passive_scalar ) THEN 286 268 ! 287 269 !-- 2D-humidity/scalar arrays 288 ALLOCATE ( qs(nysg:nyng,nxlg:nxrg), & 289 qsws_1(nysg:nyng,nxlg:nxrg), & 290 qswst_1(nysg:nyng,nxlg:nxrg) ) 291 292 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 293 ALLOCATE( qsws_2(nysg:nyng,nxlg:nxrg), & 294 qswst_2(nysg:nyng,nxlg:nxrg) ) 295 ENDIF 270 ALLOCATE ( qs(nysg:nyng,nxlg:nxrg), & 271 qsws(nysg:nyng,nxlg:nxrg), & 272 qswst(nysg:nyng,nxlg:nxrg) ) 273 296 274 ! 297 275 !-- 3D-humidity/scalar arrays … … 304 282 IF ( humidity ) THEN 305 283 ALLOCATE( vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 306 307 IF ( timestep_scheme(1:5) /= 'runge' ) THEN308 ALLOCATE( vpt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )309 ENDIF310 284 311 285 IF ( cloud_physics ) THEN … … 334 308 335 309 IF ( ocean ) THEN 336 ALLOCATE( saswsb _1(nysg:nyng,nxlg:nxrg), &337 saswst _1(nysg:nyng,nxlg:nxrg) )310 ALLOCATE( saswsb(nysg:nyng,nxlg:nxrg), & 311 saswst(nysg:nyng,nxlg:nxrg) ) 338 312 ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 339 313 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & … … 434 408 ! 435 409 !-- Initial assignment of the pointers 436 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 437 438 rif_m => rif_1; rif => rif_2 439 shf_m => shf_1; shf => shf_2 440 tswst_m => tswst_1; tswst => tswst_2 441 usws_m => usws_1; usws => usws_2 442 uswst_m => uswst_1; uswst => uswst_2 443 vsws_m => vsws_1; vsws => vsws_2 444 vswst_m => vswst_1; vswst => vswst_2 445 e_m => e_1; e => e_2; e_p => e_3; te_m => e_3 446 kh_m => kh_1; kh => kh_2 447 km_m => km_1; km => km_2 448 pt_m => pt_1; pt => pt_2; pt_p => pt_3; tpt_m => pt_3 449 u_m => u_1; u => u_2; u_p => u_3; tu_m => u_3 450 v_m => v_1; v => v_2; v_p => v_3; tv_m => v_3 451 w_m => w_1; w => w_2; w_p => w_3; tw_m => w_3 452 453 IF ( humidity .OR. passive_scalar ) THEN 454 qsws_m => qsws_1; qsws => qsws_2 455 qswst_m => qswst_1; qswst => qswst_2 456 q_m => q_1; q => q_2; q_p => q_3; tq_m => q_3 457 IF ( humidity ) vpt_m => vpt_1; vpt => vpt_2 458 IF ( cloud_physics ) ql => ql_1 459 IF ( cloud_droplets ) THEN 460 ql => ql_1 461 ql_c => ql_2 462 ENDIF 463 ENDIF 464 465 ELSE 466 467 rif => rif_1 468 shf => shf_1 469 tswst => tswst_1 470 usws => usws_1 471 uswst => uswst_1 472 vsws => vsws_1 473 vswst => vswst_1 474 e => e_1; e_p => e_2; te_m => e_3; e_m => e_3 475 kh => kh_1 476 km => km_1 477 pt => pt_1; pt_p => pt_2; tpt_m => pt_3; pt_m => pt_3 478 u => u_1; u_p => u_2; tu_m => u_3; u_m => u_3 479 v => v_1; v_p => v_2; tv_m => v_3; v_m => v_3 480 w => w_1; w_p => w_2; tw_m => w_3; w_m => w_3 481 482 IF ( humidity .OR. passive_scalar ) THEN 483 qsws => qsws_1 484 qswst => qswst_1 485 q => q_1; q_p => q_2; tq_m => q_3; q_m => q_3 486 IF ( humidity ) vpt => vpt_1 487 IF ( cloud_physics ) ql => ql_1 488 IF ( cloud_droplets ) THEN 489 ql => ql_1 490 ql_c => ql_2 491 ENDIF 492 ENDIF 493 494 IF ( ocean ) THEN 495 saswsb => saswsb_1 496 saswst => saswst_1 497 sa => sa_1; sa_p => sa_2; tsa_m => sa_3 498 ENDIF 499 500 ENDIF 501 410 e => e_1; e_p => e_2; te_m => e_3 411 pt => pt_1; pt_p => pt_2; tpt_m => pt_3 412 u => u_1; u_p => u_2; tu_m => u_3 413 v => v_1; v_p => v_2; tv_m => v_3 414 w => w_1; w_p => w_2; tw_m => w_3 415 416 IF ( humidity .OR. passive_scalar ) THEN 417 q => q_1; q_p => q_2; tq_m => q_3 418 IF ( humidity ) vpt => vpt_1 419 IF ( cloud_physics ) ql => ql_1 420 IF ( cloud_droplets ) THEN 421 ql => ql_1 422 ql_c => ql_2 423 ENDIF 424 ENDIF 425 426 IF ( ocean ) THEN 427 sa => sa_1; sa_p => sa_2; tsa_m => sa_3 428 ENDIF 429 502 430 ! 503 431 !-- Allocate arrays containing the RK coefficient for calculation of … … 785 713 ENDIF 786 714 ENDIF 787 IF ( ASSOCIATED( shf_m ) ) shf_m = shf788 715 ENDIF 789 716 … … 806 733 ENDDO 807 734 ENDIF 808 IF ( ASSOCIATED( qsws_m ) ) qsws_m = qsws809 735 ENDIF 810 736 ENDIF … … 822 748 !-- Heat flux is prescribed 823 749 tswst = top_heatflux 824 IF ( ASSOCIATED( tswst_m ) ) tswst_m = tswst 825 826 IF ( humidity .OR. passive_scalar ) THEN 827 qswst = 0.0 828 IF ( ASSOCIATED( qswst_m ) ) qswst_m = qswst 829 ENDIF 750 751 IF ( humidity .OR. passive_scalar ) qswst = 0.0 830 752 831 753 IF ( ocean ) THEN … … 839 761 IF ( coupling_mode == 'ocean_to_atmosphere' ) THEN 840 762 tswst = 0.0 841 IF ( ASSOCIATED( tswst_m ) ) tswst_m = tswst842 763 ENDIF 843 764 … … 859 780 !-- value in the course of the first few time steps. 860 781 shf = 0.0 861 IF ( ASSOCIATED( shf_m ) ) shf_m = 0.0862 782 ENDIF 863 783 864 784 IF ( humidity .OR. passive_scalar ) THEN 865 IF ( .NOT. constant_waterflux ) THEN 866 qsws = 0.0 867 IF ( ASSOCIATED( qsws_m ) ) qsws_m = 0.0 868 ENDIF 785 IF ( .NOT. constant_waterflux ) qsws = 0.0 869 786 ENDIF 870 787 … … 923 840 ! 924 841 !-- Initialize old and new time levels. 925 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 926 e_m = e; pt_m = pt; u_m = u; v_m = v; w_m = w; kh_m = kh; km_m = km 927 ELSE 928 te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0 929 ENDIF 842 te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0 930 843 e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w 931 844 932 845 IF ( humidity .OR. passive_scalar ) THEN 933 IF ( ASSOCIATED( q_m ) ) q_m = q 934 IF ( timestep_scheme(1:5) == 'runge' ) tq_m = 0.0 846 tq_m = 0.0 935 847 q_p = q 936 IF ( humidity .AND. ASSOCIATED( vpt_m ) ) vpt_m = vpt937 848 ENDIF 938 849 … … 1080 991 !-- Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at present, 1081 992 !-- maybe revise later. 1082 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1083 DO i = nxlg, nxrg 1084 DO j = nysg, nyng 1085 u (nzb:nzb_u_inner(j,i),j,i) = 0.0 1086 v (nzb:nzb_v_inner(j,i),j,i) = 0.0 1087 w (nzb:nzb_w_inner(j,i),j,i) = 0.0 1088 e (nzb:nzb_w_inner(j,i),j,i) = 0.0 1089 u_m(nzb:nzb_u_inner(j,i),j,i) = 0.0 1090 v_m(nzb:nzb_v_inner(j,i),j,i) = 0.0 1091 w_m(nzb:nzb_w_inner(j,i),j,i) = 0.0 1092 e_m(nzb:nzb_w_inner(j,i),j,i) = 0.0 1093 tu_m(nzb:nzb_u_inner(j,i),j,i) = 0.0 1094 tv_m(nzb:nzb_v_inner(j,i),j,i) = 0.0 1095 tw_m(nzb:nzb_w_inner(j,i),j,i) = 0.0 1096 te_m(nzb:nzb_w_inner(j,i),j,i) = 0.0 1097 tpt_m(nzb:nzb_w_inner(j,i),j,i) = 0.0 1098 ENDDO 1099 ENDDO 1100 ELSE 1101 DO i = nxlg, nxrg 1102 DO j = nysg, nyng 1103 u (nzb:nzb_u_inner(j,i),j,i) = 0.0 1104 v (nzb:nzb_v_inner(j,i),j,i) = 0.0 1105 w (nzb:nzb_w_inner(j,i),j,i) = 0.0 1106 e (nzb:nzb_w_inner(j,i),j,i) = 0.0 1107 u_m(nzb:nzb_u_inner(j,i),j,i) = 0.0 1108 v_m(nzb:nzb_v_inner(j,i),j,i) = 0.0 1109 w_m(nzb:nzb_w_inner(j,i),j,i) = 0.0 1110 e_m(nzb:nzb_w_inner(j,i),j,i) = 0.0 1111 u_p(nzb:nzb_u_inner(j,i),j,i) = 0.0 1112 v_p(nzb:nzb_v_inner(j,i),j,i) = 0.0 1113 w_p(nzb:nzb_w_inner(j,i),j,i) = 0.0 1114 e_p(nzb:nzb_w_inner(j,i),j,i) = 0.0 1115 ENDDO 1116 ENDDO 1117 ENDIF 993 DO i = nxlg, nxrg 994 DO j = nysg, nyng 995 u (nzb:nzb_u_inner(j,i),j,i) = 0.0 996 v (nzb:nzb_v_inner(j,i),j,i) = 0.0 997 w (nzb:nzb_w_inner(j,i),j,i) = 0.0 998 e (nzb:nzb_w_inner(j,i),j,i) = 0.0 999 tu_m(nzb:nzb_u_inner(j,i),j,i) = 0.0 1000 tv_m(nzb:nzb_v_inner(j,i),j,i) = 0.0 1001 tw_m(nzb:nzb_w_inner(j,i),j,i) = 0.0 1002 te_m(nzb:nzb_w_inner(j,i),j,i) = 0.0 1003 tpt_m(nzb:nzb_w_inner(j,i),j,i) = 0.0 1004 ENDDO 1005 ENDDO 1118 1006 1119 1007 ENDIF … … 1135 1023 !-- have to be predefined here because they are used (but multiplied with 0) 1136 1024 !-- there before they are set. 1137 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1138 te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0 1139 IF ( humidity .OR. passive_scalar ) tq_m = 0.0 1140 IF ( ocean ) tsa_m = 0.0 1141 ENDIF 1025 te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0 1026 IF ( humidity .OR. passive_scalar ) tq_m = 0.0 1027 IF ( ocean ) tsa_m = 0.0 1142 1028 1143 1029 ELSE … … 1424 1310 shf(:,:) = canopy_heat_flux(0,:,:) 1425 1311 1426 IF ( ASSOCIATED( shf_m ) ) shf_m = shf1427 1428 1312 ENDIF 1429 1313 … … 1477 1361 weight_pres(2) = 1./2. 1478 1362 1479 ELSE ! for Euler- and leapfrog-method1363 ELSE ! for Euler-method 1480 1364 1481 1365 weight_substep(1) = 1.0 -
palm/trunk/SOURCE/init_advec.f90
r667 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! all actions concerning upstream-spline-method removed 7 7 ! 8 8 ! Former revisions: … … 80 80 ENDIF 81 81 82 IF ( momentum_advec == 'ups-scheme' .OR. scalar_advec == 'ups-scheme' ) &83 THEN84 85 !86 !-- Provide the constant parameters for the Upstream-Spline advection scheme.87 !-- In x- und y-direction the Sherman-Morrison formula is applied88 !-- (cf. Press et al, 1986 (Numerical Recipes)).89 !90 !-- Allocate nonlocal arrays91 ALLOCATE( spl_z_x(0:nx), spl_z_y(0:ny), spl_tri_x(5,0:nx), &92 spl_tri_y(5,0:ny), spl_tri_zu(5,nzb:nzt+1), &93 spl_tri_zw(5,nzb:nzt+1) )94 95 !96 !-- Provide diagonal elements of the tridiagonal matrices for all97 !-- directions98 spl_tri_x(1,:) = 2.099 spl_tri_y(1,:) = 2.0100 spl_tri_zu(1,:) = 2.0101 spl_tri_zw(1,:) = 2.0102 103 !104 !-- Elements of the cyclic tridiagonal matrix105 !-- (same for all horizontal directions)106 spl_alpha = 0.5107 spl_beta = 0.5108 109 !110 !-- Sub- and superdiagonal elements, x-direction111 spl_tri_x(2,0:nx) = 0.5112 spl_tri_x(3,0:nx) = 0.5113 114 !115 !-- mMdify the diagonal elements (Sherman-Morrison)116 spl_gamma_x = -spl_tri_x(1,0)117 spl_tri_x(1,0) = spl_tri_x(1,0) - spl_gamma_x118 spl_tri_x(1,nx) = spl_tri_x(1,nx) - spl_alpha * spl_beta / spl_gamma_x119 120 !121 !-- Split the tridiagonal matrix for Thomas algorithm122 spl_tri_x(4,0) = spl_tri_x(1,0)123 DO i = 1, nx124 spl_tri_x(5,i) = spl_tri_x(2,i) / spl_tri_x(4,i-1)125 spl_tri_x(4,i) = spl_tri_x(1,i) - spl_tri_x(5,i) * spl_tri_x(3,i-1)126 ENDDO127 128 !129 !-- Allocate arrays required locally130 ALLOCATE( temp(0:nx), spl_u(0:nx) )131 132 !133 !-- Provide "corrective vector", x-direction134 spl_u(0) = spl_gamma_x135 spl_u(1:nx-1) = 0.0136 spl_u(nx) = spl_alpha137 138 !139 !-- Solve the system of equations for the corrective vector140 !-- (Sherman-Morrison)141 temp(0) = spl_u(0)142 DO i = 1, nx143 temp(i) = spl_u(i) - spl_tri_x(5,i) * temp(i-1)144 ENDDO145 spl_z_x(nx) = temp(nx) / spl_tri_x(4,nx)146 DO i = nx-1, 0, -1147 spl_z_x(i) = ( temp(i) - spl_tri_x(3,i) * spl_z_x(i+1) ) / &148 spl_tri_x(4,i)149 ENDDO150 151 !152 !-- Deallocate local arrays, for they are allocated in a different way for153 !-- operations in y-direction154 DEALLOCATE( temp, spl_u )155 156 !157 !-- Provide sub- and superdiagonal elements, y-direction158 spl_tri_y(2,0:ny) = 0.5159 spl_tri_y(3,0:ny) = 0.5160 161 !162 !-- Modify the diagonal elements (Sherman-Morrison)163 spl_gamma_y = -spl_tri_y(1,0)164 spl_tri_y(1,0) = spl_tri_y(1,0) - spl_gamma_y165 spl_tri_y(1,ny) = spl_tri_y(1,ny) - spl_alpha * spl_beta / spl_gamma_y166 167 !168 !-- Split the tridiagonal matrix for Thomas algorithm169 spl_tri_y(4,0) = spl_tri_y(1,0)170 DO j = 1, ny171 spl_tri_y(5,j) = spl_tri_y(2,j) / spl_tri_y(4,j-1)172 spl_tri_y(4,j) = spl_tri_y(1,j) - spl_tri_y(5,j) * spl_tri_y(3,j-1)173 ENDDO174 175 !176 !-- Allocate arrays required locally177 ALLOCATE( temp(0:ny), spl_u(0:ny) )178 179 !180 !-- Provide "corrective vector", y-direction181 spl_u(0) = spl_gamma_y182 spl_u(1:ny-1) = 0.0183 spl_u(ny) = spl_alpha184 185 !186 !-- Solve the system of equations for the corrective vector187 !-- (Sherman-Morrison)188 temp = 0.0189 spl_z_y = 0.0190 temp(0) = spl_u(0)191 DO j = 1, ny192 temp(j) = spl_u(j) - spl_tri_y(5,j) * temp(j-1)193 ENDDO194 spl_z_y(ny) = temp(ny) / spl_tri_y(4,ny)195 DO j = ny-1, 0, -1196 spl_z_y(j) = ( temp(j) - spl_tri_y(3,j) * spl_z_y(j+1) ) / &197 spl_tri_y(4,j)198 ENDDO199 200 !201 !-- deallocate local arrays, for they are no longer required202 DEALLOCATE( temp, spl_u )203 204 !205 !-- provide sub- and superdiagonal elements, z-direction206 spl_tri_zu(2,nzb) = 0.0207 spl_tri_zu(2,nzt+1) = 1.0208 spl_tri_zw(2,nzb) = 0.0209 spl_tri_zw(2,nzt+1) = 1.0210 211 spl_tri_zu(3,nzb) = 1.0212 spl_tri_zu(3,nzt+1) = 0.0213 spl_tri_zw(3,nzb) = 1.0214 spl_tri_zw(3,nzt+1) = 0.0215 216 DO k = nzb+1, nzt217 spl_tri_zu(2,k) = dzu(k) / ( dzu(k) + dzu(k+1) )218 spl_tri_zw(2,k) = dzw(k) / ( dzw(k) + dzw(k+1) )219 spl_tri_zu(3,k) = 1.0 - spl_tri_zu(2,k)220 spl_tri_zw(3,k) = 1.0 - spl_tri_zw(2,k)221 ENDDO222 223 spl_tri_zu(4,nzb) = spl_tri_zu(1,nzb)224 spl_tri_zw(4,nzb) = spl_tri_zw(1,nzb)225 DO k = nzb+1, nzt+1226 spl_tri_zu(5,k) = spl_tri_zu(2,k) / spl_tri_zu(4,k-1)227 spl_tri_zw(5,k) = spl_tri_zw(2,k) / spl_tri_zw(4,k-1)228 spl_tri_zu(4,k) = spl_tri_zu(1,k) - spl_tri_zu(5,k) * &229 spl_tri_zu(3,k-1)230 spl_tri_zw(4,k) = spl_tri_zw(1,k) - spl_tri_zw(5,k) * &231 spl_tri_zw(3,k-1)232 ENDDO233 234 ENDIF235 236 82 END SUBROUTINE init_advec -
palm/trunk/SOURCE/init_pegrid.f90
r979 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! all actions concerning upstream-spline-method removed 7 7 ! 8 8 ! Former revisions: … … 384 384 ! 385 385 !-- 1. transposition z --> x 386 !-- This transposition is not neccessary in case of a 1d-decomposition along x, 387 !-- except that the uptream-spline method is switched on 388 IF ( pdims(2) /= 1 .OR. momentum_advec == 'ups-scheme' .OR. & 389 scalar_advec == 'ups-scheme' ) THEN 390 391 IF ( pdims(2) == 1 .AND. ( momentum_advec == 'ups-scheme' .OR. & 392 scalar_advec == 'ups-scheme' ) ) THEN 393 message_string = '1d-decomposition along x ' // & 394 'chosen but nz restrictions may occur' // & 395 '& since ups-scheme is activated' 396 CALL message( 'init_pegrid', 'PA0229', 0, 1, 0, 6, 0 ) 397 ENDIF 386 !-- This transposition is not neccessary in case of a 1d-decomposition along x 387 IF ( pdims(2) /= 1 ) THEN 388 398 389 nys_x = nys 399 390 nyn_xa = nyna … … 448 439 !-- 3. transposition y --> z (ELSE: x --> y in case of 1D-decomposition 449 440 !-- along x) 450 IF ( pdims(2) /= 1 .OR. momentum_advec == 'ups-scheme' .OR. & 451 scalar_advec == 'ups-scheme' ) THEN 441 IF ( pdims(2) /= 1 ) THEN 452 442 ! 453 443 !-- y --> z -
palm/trunk/SOURCE/modules.f90
r997 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! -asselin_filter_factor, cut_spline_overshoot, dt_changed, last_dt_change, 7 ! last_dt_change_1d, long_filter_factor, overshoot_limit_*, ups_limit_* 8 ! several pointer/target arrays converted to normal ones 7 9 ! 8 10 ! Former revisions: … … 343 345 diss_s_v, diss_s_w, dzu_mg, dzw_mg, flux_s_e, flux_s_pt, flux_s_q, & 344 346 flux_s_sa, flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg, & 345 mean_inflow_profiles, pt_slope_ref, qs, qswst_remote, total_2d_a, & 346 total_2d_o, ts, us, z0, z0h 347 348 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET :: & 349 qsws_1, qsws_2, qswst_1, qswst_2, rif_1, rif_2, saswsb_1, saswst_1, & 350 shf_1, shf_2, tswst_1, tswst_2, usws_1, usws_2, uswst_1, uswst_2, & 351 vsws_1, vsws_2, vswst_1, vswst_2 352 353 REAL, DIMENSION(:,:), POINTER :: & 354 qsws, qsws_m, qswst, qswst_m, rif, rif_m, saswsb, saswst, shf, & 355 shf_m, tswst, tswst_m, usws, uswst, usws_m, uswst_m, vsws, vswst, & 356 vsws_m, vswst_m 347 mean_inflow_profiles, pt_slope_ref, qs, qsws, qswst, qswst_remote, & 348 rif, saswsb, saswst, shf, total_2d_a, total_2d_o, ts, tswst, us, & 349 usws, uswst, vsws, vswst, z0, z0h 357 350 358 351 REAL, DIMENSION(:,:,:), ALLOCATABLE :: & … … 360 353 diss_l_pt, diss_l_q, diss_l_sa, diss_l_u, diss_l_v, diss_l_w, & 361 354 flux_l_e, flux_l_pt, flux_l_q, flux_l_sa, flux_l_u, flux_l_v, & 362 flux_l_w, lad_s, lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls,&363 tend, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l,&364 w_m_n, w_m_r, w_m_s355 flux_l_w, kh, km, lad_s, lad_u, lad_v, lad_w, lai, l_wall, p_loc, & 356 sec, sls, tend, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, & 357 v_m_s, w_m_l, w_m_n, w_m_r, w_m_s 365 358 366 359 REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & 367 ql_v, ql_vp 368 369 REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & 370 e_1, e_2, e_3, kh_1, kh_2, km_1, km_2, p, prho_1, pt_1, pt_2, pt_3, & 371 q_1, q_2, q_3, ql_1, ql_2, rho_1, sa_1, sa_2, sa_3, u_1, u_2, u_3, & 372 v_1, v_2, v_3, vpt_1, vpt_2, w_1, w_2, w_3 360 e_1, e_2, e_3, p, prho_1, pt_1, pt_2, pt_3, q_1, q_2, q_3, ql_v, & 361 ql_vp, ql_1, ql_2, rho_1, sa_1, sa_2, sa_3, u_1, u_2, u_3, & 362 v_1, v_2, v_3, vpt_1, w_1, w_2, w_3 373 363 374 364 REAL, DIMENSION(:,:,:), POINTER :: & 375 e, e_m, e_p, kh, kh_m, km, km_m, prho, pt, pt_m, pt_p, q, q_m, q_p, & 376 ql, ql_c, rho, sa, sa_p, te_m, tpt_m, tq_m, tsa_m, tu_m, tv_m, tw_m, & 377 u, u_m, u_p, v, v_m, v_p, vpt, vpt_m, w, w_m, w_p 365 e, e_p, prho, pt, pt_p, q, q_p, ql, ql_c, rho, sa, sa_p, te_m, & 366 tpt_m, tq_m, tsa_m, tu_m, tv_m, tw_m, u, u_p, v, v_p, vpt, w, w_p 378 367 379 368 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: rif_wall … … 555 544 intermediate_timestep_count, intermediate_timestep_count_max, & 556 545 io_group = 0, io_blocks = 1, iran = -1234567, & 557 last_dt_change = 0,masks = 0, maximum_grid_level, &546 masks = 0, maximum_grid_level, & 558 547 maximum_parallel_io_streams = -1, max_pr_user = 0, & 559 548 mgcycles = 0, mg_cycles = -1, mg_switch_to_pe0_level = 0, mid, & … … 605 594 constant_top_salinityflux = .TRUE., & 606 595 constant_waterflux = .TRUE., create_disturbances = .TRUE., & 607 cut_spline_overshoot = .TRUE., &608 596 data_output_2d_on_each_pe = .TRUE., & 609 597 dissipation_control = .FALSE., disturbance_created = .FALSE., & 610 598 do2d_at_begin = .FALSE., do3d_at_begin = .FALSE., & 611 599 do3d_compress = .FALSE., do_sum = .FALSE., & 612 dp_external = .FALSE., dp_smooth = .FALSE., & 613 dt_changed = .FALSE., dt_fixed = .FALSE., & 600 dp_external = .FALSE., dp_smooth = .FALSE., dt_fixed = .FALSE., & 614 601 dt_3d_reached, dt_3d_reached_l, exchange_mg = .FALSE., & 615 602 first_call_lpm = .TRUE., & … … 640 627 641 628 REAL :: advected_distance_x = 0.0, advected_distance_y = 0.0, & 642 alpha_surface = 0.0, asselin_filter_factor = 0.1, & 643 atmos_ocean_sign = 1.0, & 629 alpha_surface = 0.0, atmos_ocean_sign = 1.0, & 644 630 averaging_interval = 0.0, averaging_interval_pr = 9999999.9, & 645 631 averaging_interval_sp = 9999999.9, bc_pt_t_val, bc_q_t_val, & … … 671 657 f = 0.0, fs = 0.0, g = 9.81, inflow_damping_height = 9999999.9, & 672 658 inflow_damping_width = 9999999.9, kappa = 0.4, km_constant = -1.0,& 673 lad_surface = 0.0, & 674 leaf_surface_concentration = 0.0, long_filter_factor = 0.0, & 659 lad_surface = 0.0, leaf_surface_concentration = 0.0, & 675 660 mask_scale_x = 1.0, mask_scale_y = 1.0, mask_scale_z = 1.0, & 676 661 maximum_cpu_time_allowed = 0.0, & 677 662 molecular_viscosity = 1.461E-5, & 678 663 old_dt = 1.0E-10, omega = 7.29212E-5, omega_sor = 1.8, & 679 overshoot_limit_e = 0.0, overshoot_limit_pt = 0.0, & 680 overshoot_limit_u = 0.0, overshoot_limit_v = 0.0, & 681 overshoot_limit_w = 0.0, particle_maximum_age = 9999999.9, & 664 particle_maximum_age = 9999999.9, & 682 665 phi = 55.0, prandtl_number = 1.0, & 683 666 precipitation_amount_interval = 9999999.9, prho_reference, & … … 710 693 top_momentumflux_v = 9999999.9, top_salinityflux = 9999999.9, & 711 694 ug_surface = 0.0, u_bulk = 0.0, u_gtrans = 0.0, & 712 ups_limit_e = 0.0, ups_limit_pt = 0.0, ups_limit_u = 0.0, & 713 ups_limit_v = 0.0, ups_limit_w = 0.0, vg_surface = 0.0, & 695 vg_surface = 0.0, & 714 696 v_bulk = 0.0, v_gtrans = 0.0, wall_adjustment_factor = 1.8, & 715 697 z_max_do2d = -1.0, z0h_factor = 1.0 … … 730 712 skip_time_domask(max_masks) = 9999999.9, threshold(20) = 0.0, & 731 713 time_domask(max_masks) = 0.0, & 732 tsc(10) = (/ 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /), &714 tsc(10) = (/ 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /), & 733 715 u_profile(100) = 9999999.9, uv_heights(100) = 9999999.9, & 734 716 v_profile(100) = 9999999.9, & … … 1060 1042 !------------------------------------------------------------------------------! 1061 1043 1062 INTEGER :: current_timestep_number_1d = 0, damp_level_ind_1d, & 1063 last_dt_change_1d = 0 1044 INTEGER :: current_timestep_number_1d = 0, damp_level_ind_1d 1064 1045 1065 1046 LOGICAL :: run_control_header_1d = .FALSE., stop_dt_1d = .FALSE. … … 1069 1050 end_time_1d = 864000.0, old_dt_1d = 1.0E-10, & 1070 1051 qs1d, simulated_time_1d = 0.0, time_pr_1d = 0.0, & 1071 time_run_control_1d = 0.0, ts1d, us1d, usws1d, usws1d_m, & 1072 vsws1d, vsws1d_m, z01d, z0h1d 1073 1074 1075 REAL, DIMENSION(:), ALLOCATABLE :: e1d, e1d_m, e1d_p, kh1d, kh1d_m, km1d, & 1076 km1d_m, l_black, l1d, l1d_m, rif1d, & 1077 te_e, te_em, te_u, te_um, te_v, te_vm, & 1078 u1d, u1d_m, u1d_p, v1d, v1d_m, v1d_p 1052 time_run_control_1d = 0.0, ts1d, us1d, usws1d, & 1053 vsws1d, z01d, z0h1d 1054 1055 1056 REAL, DIMENSION(:), ALLOCATABLE :: e1d, e1d_p, kh1d, km1d, l_black, l1d, & 1057 rif1d, te_e, te_em, te_u, te_um, te_v, & 1058 te_vm, u1d, u1d_p, v1d, v1d_p 1079 1059 1080 1060 SAVE -
palm/trunk/SOURCE/parin.f90
r997 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! -cut_spline_overshoot, long_filter_factor, overshoot_limit_*, ups_limit_* 7 7 ! 8 8 ! Former revisions: … … 168 168 collective_wait, conserve_volume_flow, conserve_volume_flow_mode, & 169 169 coupling_start_time, cthf, curvature_solution_effects, & 170 cut_spline_overshoot, &171 170 cycle_mg, damp_level_1d, dissipation_1d, & !dissipation_control, & 172 171 dp_external, dp_level_b, dp_smooth, dpdxy, drag_coefficient, & … … 178 177 initializing_actions, km_constant, lad_surface, & 179 178 lad_vertical_gradient, lad_vertical_gradient_level, & 180 leaf_surface_concentration, long_filter_factor,&179 leaf_surface_concentration, & 181 180 loop_optimization, masking_method, mg_cycles, & 182 181 mg_switch_to_pe0_level, mixing_length_1d, momentum_advec, & 183 182 netcdf_precision, neutral, ngsrb, nsor, & 184 183 nsor_ini, nx, ny, nz, ocean, omega, omega_sor, & 185 overshoot_limit_e, overshoot_limit_pt, &186 overshoot_limit_u, overshoot_limit_v, overshoot_limit_w, &187 184 passive_scalar, pch_index, phi, plant_canopy, prandtl_layer, & 188 185 prandtl_number, precipitation, psolver, pt_damping_factor, & … … 204 201 top_momentumflux_u, top_momentumflux_v, top_salinityflux, & 205 202 turbulent_inflow, ug_surface, ug_vertical_gradient, & 206 ug_vertical_gradient_level, ups_limit_e, ups_limit_pt, & 207 ups_limit_u, ups_limit_v, ups_limit_w, use_surface_fluxes, & 203 ug_vertical_gradient_level, use_surface_fluxes, & 208 204 use_top_fluxes, use_ug_for_galilei_tr, use_upstream_for_tke, & 209 205 uv_heights, u_bulk, u_profile, vg_surface, vg_vertical_gradient, & -
palm/trunk/SOURCE/prognostic_equations.f90
r979 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! all actions concerning leapfrog- and upstream-spline-scheme removed 7 7 ! 8 8 ! Former revisions: … … 163 163 CHARACTER (LEN=9) :: time_to_string 164 164 INTEGER :: i, i_omp_start, j, k, tn = 0 165 REAL :: s at, sbt165 REAL :: sbt 166 166 167 167 ! … … 178 178 CALL cpu_log( log_point(5), 'u-equation', 'start' ) 179 179 180 !181 !-- u-tendency terms with communication182 IF ( momentum_advec == 'ups-scheme' ) THEN183 tend = 0.0184 CALL advec_u_ups185 ENDIF186 187 !188 !-- u-tendency terms with no communication189 180 i_omp_start = nxlu 190 181 DO i = nxlu, nxr … … 192 183 ! 193 184 !-- Tendency terms 194 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN195 tend(:,j,i) = 0.0185 tend(:,j,i) = 0.0 186 IF ( timestep_scheme(1:5) == 'runge' ) THEN 196 187 IF ( ws_scheme_mom ) THEN 197 188 CALL advec_u_ws( i, j, i_omp_start, tn ) … … 199 190 CALL advec_u_pw( i, j ) 200 191 ENDIF 201 202 192 ELSE 203 IF ( momentum_advec /= 'ups-scheme' ) THEN 204 tend(:,j,i) = 0.0 205 CALL advec_u_up( i, j ) 206 ENDIF 207 ENDIF 208 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 209 CALL diffusion_u( i, j, ddzu, ddzw, km_m, tend, u_m, usws_m, & 210 uswst_m, v_m, w_m ) 211 ELSE 212 CALL diffusion_u( i, j, ddzu, ddzw, km, tend, u, usws, uswst, & 213 v, w ) 214 ENDIF 193 CALL advec_u_up( i, j ) 194 ENDIF 195 CALL diffusion_u( i, j ) 215 196 CALL coriolis( i, j, 1 ) 216 197 IF ( sloping_surface .AND. .NOT. neutral ) THEN … … 235 216 !-- Prognostic equation for u-velocity component 236 217 DO k = nzb_u_inner(j,i)+1, nzt 237 u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + & 238 dt_3d * ( & 239 tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) & 240 ) - & 241 tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) ) 218 u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 219 tsc(3) * tu_m(k,j,i) ) & 220 - tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) ) 242 221 ENDDO 243 222 … … 266 245 CALL cpu_log( log_point(6), 'v-equation', 'start' ) 267 246 268 !269 !-- v-tendency terms with communication270 IF ( momentum_advec == 'ups-scheme' ) THEN271 tend = 0.0272 CALL advec_v_ups273 ENDIF274 275 !276 !-- v-tendency terms with no communication277 247 i_omp_start = nxl 278 248 DO i = nxl, nxr … … 280 250 ! 281 251 !-- Tendency terms 282 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN283 tend(:,j,i) = 0.0252 tend(:,j,i) = 0.0 253 IF ( timestep_scheme(1:5) == 'runge' ) THEN 284 254 IF ( ws_scheme_mom ) THEN 285 255 CALL advec_v_ws( i, j, i_omp_start, tn ) … … 289 259 290 260 ELSE 291 IF ( momentum_advec /= 'ups-scheme' ) THEN 292 tend(:,j,i) = 0.0 293 CALL advec_v_up( i, j ) 294 ENDIF 295 ENDIF 296 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 297 CALL diffusion_v( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, & 298 vsws_m, vswst_m, w_m ) 299 ELSE 300 CALL diffusion_v( i, j, ddzu, ddzw, km, tend, u, v, vsws, & 301 vswst, w ) 302 ENDIF 261 CALL advec_v_up( i, j ) 262 ENDIF 263 CALL diffusion_v( i, j ) 303 264 CALL coriolis( i, j, 2 ) 304 265 … … 320 281 !-- Prognostic equation for v-velocity component 321 282 DO k = nzb_v_inner(j,i)+1, nzt 322 v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + & 323 dt_3d * ( & 324 tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) & 325 ) - & 326 tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) ) 283 v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 284 tsc(3) * tv_m(k,j,i) ) & 285 - tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) ) 327 286 ENDDO 328 287 … … 351 310 CALL cpu_log( log_point(7), 'w-equation', 'start' ) 352 311 353 !354 !-- w-tendency terms with communication355 IF ( momentum_advec == 'ups-scheme' ) THEN356 tend = 0.0357 CALL advec_w_ups358 ENDIF359 360 !361 !-- w-tendency terms with no communication362 312 DO i = nxl, nxr 363 313 DO j = nys, nyn 364 314 ! 365 315 !-- Tendency terms 366 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN367 tend(:,j,i) = 0.0316 tend(:,j,i) = 0.0 317 IF ( timestep_scheme(1:5) == 'runge' ) THEN 368 318 IF ( ws_scheme_mom ) THEN 369 319 CALL advec_w_ws( i, j, i_omp_start, tn ) … … 373 323 374 324 ELSE 375 IF ( momentum_advec /= 'ups-scheme' ) THEN 376 tend(:,j,i) = 0.0 377 CALL advec_w_up( i, j ) 378 ENDIF 379 ENDIF 380 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 381 CALL diffusion_w( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, w_m ) 382 ELSE 383 CALL diffusion_w( i, j, ddzu, ddzw, km, tend, u, v, w ) 384 ENDIF 325 CALL advec_w_up( i, j ) 326 ENDIF 327 CALL diffusion_w( i, j ) 385 328 CALL coriolis( i, j, 3 ) 386 329 … … 406 349 !-- Prognostic equation for w-velocity component 407 350 DO k = nzb_w_inner(j,i)+1, nzt-1 408 w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + & 409 dt_3d * ( & 410 tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) & 411 ) - & 412 tsc(5) * rdf(k) * w(k,j,i) 351 w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 352 tsc(3) * tw_m(k,j,i) ) & 353 - tsc(5) * rdf(k) * w(k,j,i) 413 354 ENDDO 414 355 … … 439 380 CALL cpu_log( log_point(13), 'pt-equation', 'start' ) 440 381 441 ! 442 !-- pt-tendency terms with communication 443 sat = tsc(1) 382 ! 383 !-- pt-tendency terms with communication 444 384 sbt = tsc(2) 445 385 IF ( scalar_advec == 'bc-scheme' ) THEN 446 386 447 387 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 448 ! 449 !-- Bott-Chlond scheme always uses Euler time step when leapfrog is 450 !-- switched on. Thus: 451 sat = 1.0 388 ! 389 !-- Bott-Chlond scheme always uses Euler time step. Thus: 452 390 sbt = 1.0 453 391 ENDIF 454 392 tend = 0.0 455 393 CALL advec_s_bc( pt, 'pt' ) 456 ELSE 457 IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN 458 tend = 0.0 459 CALL advec_s_ups( pt, 'pt' ) 460 ENDIF 461 ENDIF 462 463 ! 464 !-- pt-tendency terms with no communication 394 395 ENDIF 396 397 ! 398 !-- pt-tendency terms with no communication 465 399 DO i = nxl, nxr 466 400 DO j = nys, nyn 467 ! 468 !-- Tendency terms 469 IF ( scalar_advec == 'bc-scheme' ) THEN 470 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, & 471 wall_heatflux, tend ) 472 ELSE 473 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 474 THEN 475 tend(:,j,i) = 0.0 401 ! 402 !-- Tendency terms 403 IF ( scalar_advec /= 'bc-scheme' ) THEN 404 tend(:,j,i) = 0.0 405 IF ( timestep_scheme(1:5) == 'runge' ) THEN 476 406 IF ( ws_scheme_sca ) THEN 477 407 CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, & … … 482 412 ENDIF 483 413 ELSE 484 IF ( scalar_advec /= 'ups-scheme' ) THEN 485 tend(:,j,i) = 0.0 486 CALL advec_s_up( i, j, pt ) 487 ENDIF 488 ENDIF 489 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 490 THEN 491 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, & 492 tswst_m, wall_heatflux, tend ) 493 ELSE 494 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, & 495 wall_heatflux, tend ) 496 ENDIF 497 ENDIF 498 499 ! 500 !-- If required compute heating/cooling due to long wave radiation 501 !-- processes 414 CALL advec_s_up( i, j, pt ) 415 ENDIF 416 ENDIF 417 418 CALL diffusion_s( i, j, pt, shf, tswst, wall_heatflux ) 419 420 ! 421 !-- If required compute heating/cooling due to long wave radiation 422 !-- processes 502 423 IF ( radiation ) THEN 503 424 CALL calc_radiation( i, j ) 504 425 ENDIF 505 426 506 507 !--If required compute impact of latent heat due to precipitation427 ! 428 !-- If required compute impact of latent heat due to precipitation 508 429 IF ( precipitation ) THEN 509 430 CALL impact_of_latent_heat( i, j ) 510 431 ENDIF 511 432 512 513 !--Consideration of heat sources within the plant canopy433 ! 434 !-- Consideration of heat sources within the plant canopy 514 435 IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN 515 436 CALL plant_canopy_model( i, j, 4 ) 516 437 ENDIF 517 438 518 519 !--If required compute influence of large-scale subsidence/ascent439 ! 440 !-- If required compute influence of large-scale subsidence/ascent 520 441 IF ( large_scale_subsidence ) THEN 521 442 CALL subsidence( i, j, tend, pt, pt_init ) … … 524 445 CALL user_actions( i, j, 'pt-tendency' ) 525 446 526 527 !--Prognostic equation for potential temperature447 ! 448 !-- Prognostic equation for potential temperature 528 449 DO k = nzb_s_inner(j,i)+1, nzt 529 pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + & 530 dt_3d * ( & 531 sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 532 ) - & 533 tsc(5) * ( pt(k,j,i) - pt_init(k) ) * & 534 ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) ) 535 ENDDO 536 537 ! 538 !-- Calculate tendencies for the next Runge-Kutta step 450 pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 451 tsc(3) * tpt_m(k,j,i) ) & 452 - tsc(5) * ( pt(k,j,i) - pt_init(k) ) *& 453 ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) ) 454 ENDDO 455 456 ! 457 !-- Calculate tendencies for the next Runge-Kutta step 539 458 IF ( timestep_scheme(1:5) == 'runge' ) THEN 540 459 IF ( intermediate_timestep_count == 1 ) THEN … … 566 485 ! 567 486 !-- sa-tendency terms with communication 568 sat = tsc(1)569 487 sbt = tsc(2) 570 488 IF ( scalar_advec == 'bc-scheme' ) THEN … … 572 490 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 573 491 ! 574 !-- Bott-Chlond scheme always uses Euler time step when leapfrog is 575 !-- switched on. Thus: 576 sat = 1.0 492 !-- Bott-Chlond scheme always uses Euler time step. Thus: 577 493 sbt = 1.0 578 494 ENDIF 579 495 tend = 0.0 580 496 CALL advec_s_bc( sa, 'sa' ) 581 ELSE 582 IF ( tsc(2) /= 2.0 ) THEN 583 IF ( scalar_advec == 'ups-scheme' ) THEN 584 tend = 0.0 585 CALL advec_s_ups( sa, 'sa' ) 586 ENDIF 587 ENDIF 497 588 498 ENDIF 589 499 … … 594 504 ! 595 505 !-- Tendency-terms 596 IF ( scalar_advec == 'bc-scheme' ) THEN 597 CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, & 598 wall_salinityflux, tend ) 599 ELSE 600 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 601 tend(:,j,i) = 0.0 506 IF ( scalar_advec /= 'bc-scheme' ) THEN 507 tend(:,j,i) = 0.0 508 IF ( timestep_scheme(1:5) == 'runge' ) THEN 602 509 IF ( ws_scheme_sca ) THEN 603 510 CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa, & … … 608 515 609 516 ELSE 610 IF ( scalar_advec /= 'ups-scheme' ) THEN 611 tend(:,j,i) = 0.0 612 CALL advec_s_up( i, j, sa ) 613 ENDIF 614 ENDIF 615 CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, & 616 wall_salinityflux, tend ) 617 ENDIF 517 CALL advec_s_up( i, j, sa ) 518 ENDIF 519 ENDIF 520 521 CALL diffusion_s( i, j, sa, saswsb, saswst, wall_salinityflux ) 618 522 619 523 CALL user_actions( i, j, 'sa-tendency' ) … … 622 526 !-- Prognostic equation for salinity 623 527 DO k = nzb_s_inner(j,i)+1, nzt 624 sa_p(k,j,i) = sat * sa(k,j,i) + & 625 dt_3d * ( & 626 sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) & 627 ) - & 628 tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) ) 528 sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 529 tsc(3) * tsa_m(k,j,i) ) & 530 - tsc(5) * rdf_sc(k) * & 531 ( sa(k,j,i) - sa_init(k) ) 629 532 IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i) 630 533 ENDDO … … 665 568 ! 666 569 !-- Scalar/q-tendency terms with communication 667 sat = tsc(1)668 570 sbt = tsc(2) 669 571 IF ( scalar_advec == 'bc-scheme' ) THEN … … 671 573 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 672 574 ! 673 !-- Bott-Chlond scheme always uses Euler time step when leapfrog is 674 !-- switched on. Thus: 675 sat = 1.0 575 !-- Bott-Chlond scheme always uses Euler time step. Thus: 676 576 sbt = 1.0 677 577 ENDIF 678 578 tend = 0.0 679 579 CALL advec_s_bc( q, 'q' ) 680 ELSE 681 IF ( tsc(2) /= 2.0 ) THEN 682 IF ( scalar_advec == 'ups-scheme' ) THEN 683 tend = 0.0 684 CALL advec_s_ups( q, 'q' ) 685 ENDIF 686 ENDIF 580 687 581 ENDIF 688 582 … … 693 587 ! 694 588 !-- Tendency-terms 695 IF ( scalar_advec == 'bc-scheme' ) THEN 696 CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, & 697 wall_qflux, tend ) 698 ELSE 699 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 700 tend(:,j,i) = 0.0 589 IF ( scalar_advec /= 'bc-scheme' ) THEN 590 tend(:,j,i) = 0.0 591 IF ( timestep_scheme(1:5) == 'runge' ) THEN 701 592 IF ( ws_scheme_sca ) THEN 702 593 CALL advec_s_ws( i, j, q, 'q', flux_s_q, & … … 706 597 ENDIF 707 598 ELSE 708 IF ( scalar_advec /= 'ups-scheme' ) THEN 709 tend(:,j,i) = 0.0 710 CALL advec_s_up( i, j, q ) 711 ENDIF 712 ENDIF 713 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )& 714 THEN 715 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, & 716 qswst_m, wall_qflux, tend ) 717 ELSE 718 CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, & 719 wall_qflux, tend ) 720 ENDIF 721 ENDIF 599 CALL advec_s_up( i, j, q ) 600 ENDIF 601 ENDIF 602 603 CALL diffusion_s( i, j, q, qsws, qswst, wall_qflux ) 722 604 723 605 ! … … 743 625 !-- Prognostic equation for total water content / scalar 744 626 DO k = nzb_s_inner(j,i)+1, nzt 745 q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + & 746 dt_3d * ( & 747 sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) & 748 ) - & 749 tsc(5) * rdf_sc(k) * ( q(k,j,i) - q_init(k) ) 627 q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 628 tsc(3) * tq_m(k,j,i) ) & 629 - tsc(5) * rdf_sc(k) * & 630 ( q(k,j,i) - q_init(k) ) 750 631 IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i) 751 632 ENDDO … … 784 665 CALL production_e_init 785 666 786 sat = tsc(1)787 667 sbt = tsc(2) 788 668 IF ( .NOT. use_upstream_for_tke ) THEN … … 791 671 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 792 672 ! 793 !-- Bott-Chlond scheme always uses Euler time step when leapfrog is 794 !-- switched on. Thus: 795 sat = 1.0 673 !-- Bott-Chlond scheme always uses Euler time step. Thus: 796 674 sbt = 1.0 797 675 ENDIF 798 676 tend = 0.0 799 677 CALL advec_s_bc( e, 'e' ) 800 ELSE801 IF ( tsc(2) /= 2.0 ) THEN802 IF ( scalar_advec == 'ups-scheme' ) THEN803 tend = 0.0804 CALL advec_s_ups( e, 'e' )805 ENDIF806 ENDIF807 678 ENDIF 808 679 ENDIF … … 814 685 ! 815 686 !-- Tendency-terms 816 IF ( scalar_advec == 'bc-scheme' .AND. & 817 .NOT. use_upstream_for_tke ) THEN 818 IF ( .NOT. humidity ) THEN 819 IF ( ocean ) THEN 820 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, & 821 l_grid, prho, prho_reference, rif, & 822 tend, zu, zw ) 823 ELSE 824 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, & 825 l_grid, pt, pt_reference, rif, tend, & 826 zu, zw ) 827 ENDIF 828 ELSE 829 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, & 830 l_grid, vpt, pt_reference, rif, tend, zu, & 831 zw ) 832 ENDIF 833 ELSE 687 IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN 834 688 IF ( use_upstream_for_tke ) THEN 835 689 tend(:,j,i) = 0.0 836 690 CALL advec_s_up( i, j, e ) 837 691 ELSE 838 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 839 THEN 840 tend(:,j,i) = 0.0 692 tend(:,j,i) = 0.0 693 IF ( timestep_scheme(1:5) == 'runge' ) THEN 841 694 IF ( ws_scheme_sca ) THEN 842 695 CALL advec_s_ws( i, j, e, 'e', flux_s_e, & … … 846 699 ENDIF 847 700 ELSE 848 IF ( scalar_advec /= 'ups-scheme' ) THEN 849 tend(:,j,i) = 0.0 850 CALL advec_s_up( i, j, e ) 851 ENDIF 701 CALL advec_s_up( i, j, e ) 852 702 ENDIF 853 703 ENDIF 854 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )& 855 THEN 856 IF ( .NOT. humidity ) THEN 857 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, & 858 km_m, l_grid, pt_m, pt_reference, & 859 rif_m, tend, zu, zw ) 860 ELSE 861 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, & 862 km_m, l_grid, vpt_m, pt_reference, & 863 rif_m, tend, zu, zw ) 864 ENDIF 704 ENDIF 705 706 IF ( .NOT. humidity ) THEN 707 IF ( ocean ) THEN 708 CALL diffusion_e( i, j, prho, prho_reference ) 865 709 ELSE 866 IF ( .NOT. humidity ) THEN 867 IF ( ocean ) THEN 868 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, & 869 km, l_grid, prho, prho_reference, & 870 rif, tend, zu, zw ) 871 ELSE 872 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, & 873 km, l_grid, pt, pt_reference, rif, & 874 tend, zu, zw ) 875 ENDIF 876 ELSE 877 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, & 878 l_grid, vpt, pt_reference, rif, tend, & 879 zu, zw ) 880 ENDIF 881 ENDIF 882 ENDIF 710 CALL diffusion_e( i, j, pt, pt_reference ) 711 ENDIF 712 ELSE 713 CALL diffusion_e( i, j, vpt, pt_reference ) 714 ENDIF 715 883 716 CALL production_e( i, j ) 884 717 … … 895 728 !-- value is reduced by 90%. 896 729 DO k = nzb_s_inner(j,i)+1, nzt 897 e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + & 898 dt_3d * ( & 899 sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) & 900 ) 730 e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 731 tsc(3) * te_m(k,j,i) ) 901 732 IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i) 902 733 ENDDO … … 986 817 987 818 tend(:,j,i) = 0.0 988 IF ( t sc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN819 IF ( timestep_scheme(1:5) == 'runge' ) THEN 989 820 IF ( ws_scheme_mom ) THEN 990 821 IF ( ( inflow_l .OR. outflow_l ) .AND. i_omp_start == nxl ) THEN 991 ! CALL local_diss( i, j, u) ! dissipation control992 822 CALL advec_u_ws( i, j, i_omp_start + 1, tn ) 993 823 ELSE … … 1000 830 CALL advec_u_up( i, j ) 1001 831 ENDIF 1002 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 1003 THEN 1004 CALL diffusion_u( i, j, ddzu, ddzw, km_m, tend, u_m, & 1005 usws_m, uswst_m, v_m, w_m ) 1006 ELSE 1007 CALL diffusion_u( i, j, ddzu, ddzw, km, tend, u, usws, & 1008 uswst, v, w ) 1009 ENDIF 832 CALL diffusion_u( i, j ) 1010 833 CALL coriolis( i, j, 1 ) 1011 834 IF ( sloping_surface .AND. .NOT. neutral ) THEN … … 1030 853 !-- Prognostic equation for u-velocity component 1031 854 DO k = nzb_u_inner(j,i)+1, nzt 1032 u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + & 1033 dt_3d * ( & 1034 tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) & 1035 ) - & 1036 tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) ) 855 u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 856 tsc(3) * tu_m(k,j,i) ) & 857 - tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) ) 1037 858 ENDDO 1038 859 … … 1059 880 1060 881 tend(:,j,i) = 0.0 1061 IF ( t sc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN882 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1062 883 IF ( ws_scheme_mom ) THEN 1063 ! CALL local_diss( i, j, v)1064 884 CALL advec_v_ws( i, j, i_omp_start, tn ) 1065 885 ELSE … … 1069 889 CALL advec_v_up( i, j ) 1070 890 ENDIF 1071 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 1072 THEN 1073 CALL diffusion_v( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, & 1074 vsws_m, vswst_m, w_m ) 1075 ELSE 1076 CALL diffusion_v( i, j, ddzu, ddzw, km, tend, u, v, vsws, & 1077 vswst, w ) 1078 ENDIF 891 CALL diffusion_v( i, j ) 1079 892 CALL coriolis( i, j, 2 ) 1080 893 … … 1096 909 !-- Prognostic equation for v-velocity component 1097 910 DO k = nzb_v_inner(j,i)+1, nzt 1098 v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + & 1099 dt_3d * ( & 1100 tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) & 1101 ) - & 1102 tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) ) 911 v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 912 tsc(3) * tv_m(k,j,i) ) & 913 - tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) ) 1103 914 ENDDO 1104 915 … … 1123 934 !-- Tendency terms for w-velocity component 1124 935 tend(:,j,i) = 0.0 1125 IF ( t sc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN936 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1126 937 IF ( ws_scheme_mom ) THEN 1127 ! CALL local_diss( i, j, w)1128 938 CALL advec_w_ws( i, j, i_omp_start, tn ) 1129 939 ELSE … … 1133 943 CALL advec_w_up( i, j ) 1134 944 ENDIF 1135 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 1136 THEN 1137 CALL diffusion_w( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, & 1138 w_m ) 1139 ELSE 1140 CALL diffusion_w( i, j, ddzu, ddzw, km, tend, u, v, w ) 1141 ENDIF 945 CALL diffusion_w( i, j ) 1142 946 CALL coriolis( i, j, 3 ) 1143 947 … … 1163 967 !-- Prognostic equation for w-velocity component 1164 968 DO k = nzb_w_inner(j,i)+1, nzt-1 1165 w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + & 1166 dt_3d * ( & 1167 tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) & 1168 ) - & 1169 tsc(5) * rdf(k) * w(k,j,i) 969 w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 970 tsc(3) * tw_m(k,j,i) ) & 971 - tsc(5) * rdf(k) * w(k,j,i) 1170 972 ENDDO 1171 973 … … 1191 993 !-- Tendency terms for potential temperature 1192 994 tend(:,j,i) = 0.0 1193 IF ( t sc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN995 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1194 996 IF ( ws_scheme_sca ) THEN 1195 997 CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt, & … … 1201 1003 CALL advec_s_up( i, j, pt ) 1202 1004 ENDIF 1203 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 1204 THEN 1205 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, & 1206 tswst_m, wall_heatflux, tend ) 1207 ELSE 1208 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, & 1209 wall_heatflux, tend ) 1210 ENDIF 1005 CALL diffusion_s( i, j, pt, shf, tswst, wall_heatflux ) 1211 1006 1212 1007 ! … … 1241 1036 !-- Prognostic equation for potential temperature 1242 1037 DO k = nzb_s_inner(j,i)+1, nzt 1243 pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +& 1244 dt_3d * ( & 1245 tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 1246 ) - & 1247 tsc(5) * ( pt(k,j,i) - pt_init(k) ) * & 1248 ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) ) 1038 pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 1039 tsc(3) * tpt_m(k,j,i) ) & 1040 - tsc(5) * ( pt(k,j,i) - pt_init(k) ) *& 1041 ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) ) 1249 1042 ENDDO 1250 1043 … … 1274 1067 !-- Tendency-terms for salinity 1275 1068 tend(:,j,i) = 0.0 1276 IF ( t sc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &1069 IF ( timestep_scheme(1:5) == 'runge' ) & 1277 1070 THEN 1278 1071 IF ( ws_scheme_sca ) THEN 1279 ! CALL local_diss( i, j, sa )1280 1072 CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa, & 1281 1073 diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn ) … … 1286 1078 CALL advec_s_up( i, j, sa ) 1287 1079 ENDIF 1288 CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, & 1289 wall_salinityflux, tend ) 1080 CALL diffusion_s( i, j, sa, saswsb, saswst, wall_salinityflux ) 1290 1081 1291 1082 CALL user_actions( i, j, 'sa-tendency' ) … … 1294 1085 !-- Prognostic equation for salinity 1295 1086 DO k = nzb_s_inner(j,i)+1, nzt 1296 sa_p(k,j,i) = tsc(1) * sa(k,j,i) + & 1297 dt_3d * ( & 1298 tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) & 1299 ) - & 1300 tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) ) 1087 sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 1088 tsc(3) * tsa_m(k,j,i) ) & 1089 - tsc(5) * rdf_sc(k) * & 1090 ( sa(k,j,i) - sa_init(k) ) 1301 1091 IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i) 1302 1092 ENDDO … … 1332 1122 !-- Tendency-terms for total water content / scalar 1333 1123 tend(:,j,i) = 0.0 1334 IF ( t sc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &1124 IF ( timestep_scheme(1:5) == 'runge' ) & 1335 1125 THEN 1336 1126 IF ( ws_scheme_sca ) THEN 1337 ! CALL local_diss( i, j, q )1338 1127 CALL advec_s_ws( i, j, q, 'q', flux_s_q, & 1339 1128 diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn ) … … 1344 1133 CALL advec_s_up( i, j, q ) 1345 1134 ENDIF 1346 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )& 1347 THEN 1348 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, & 1349 qswst_m, wall_qflux, tend ) 1350 ELSE 1351 CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, & 1352 wall_qflux, tend ) 1353 ENDIF 1135 CALL diffusion_s( i, j, q, qsws, qswst, wall_qflux ) 1354 1136 1355 1137 ! … … 1374 1156 !-- Prognostic equation for total water content / scalar 1375 1157 DO k = nzb_s_inner(j,i)+1, nzt 1376 q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +& 1377 dt_3d * ( & 1378 tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) & 1379 ) - & 1380 tsc(5) * rdf_sc(k) * ( q(k,j,i) - q_init(k) ) 1158 q_p(k,j,i) = q(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 1159 tsc(3) * tq_m(k,j,i) ) & 1160 - tsc(5) * rdf_sc(k) * & 1161 ( q(k,j,i) - q_init(k) ) 1381 1162 IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i) 1382 1163 ENDDO … … 1408 1189 !-- Tendency-terms for TKE 1409 1190 tend(:,j,i) = 0.0 1410 IF ( ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' )&1191 IF ( timestep_scheme(1:5) == 'runge' & 1411 1192 .AND. .NOT. use_upstream_for_tke ) THEN 1412 1193 IF ( ws_scheme_sca ) THEN 1413 ! CALL local_diss( i, j, e ) 1414 CALL advec_s_ws( i, j, e, 'e', flux_s_e, & 1415 diss_s_e, flux_l_e, diss_l_e , i_omp_start, tn ) 1194 CALL advec_s_ws( i, j, e, 'e', flux_s_e, diss_s_e, & 1195 flux_l_e, diss_l_e , i_omp_start, tn ) 1416 1196 ELSE 1417 1197 CALL advec_s_pw( i, j, e ) … … 1420 1200 CALL advec_s_up( i, j, e ) 1421 1201 ENDIF 1422 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )& 1423 THEN 1424 IF ( .NOT. humidity ) THEN 1425 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, & 1426 km_m, l_grid, pt_m, pt_reference, & 1427 rif_m, tend, zu, zw ) 1202 IF ( .NOT. humidity ) THEN 1203 IF ( ocean ) THEN 1204 CALL diffusion_e( i, j, prho, prho_reference ) 1428 1205 ELSE 1429 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, & 1430 km_m, l_grid, vpt_m, pt_reference, & 1431 rif_m, tend, zu, zw ) 1206 CALL diffusion_e( i, j, pt, pt_reference ) 1432 1207 ENDIF 1433 1208 ELSE 1434 IF ( .NOT. humidity ) THEN 1435 IF ( ocean ) THEN 1436 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, & 1437 km, l_grid, prho, prho_reference, & 1438 rif, tend, zu, zw ) 1439 ELSE 1440 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, & 1441 km, l_grid, pt, pt_reference, rif, & 1442 tend, zu, zw ) 1443 ENDIF 1444 ELSE 1445 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, & 1446 l_grid, vpt, pt_reference, rif, tend, & 1447 zu, zw ) 1448 ENDIF 1209 CALL diffusion_e( i, j, vpt, pt_reference ) 1449 1210 ENDIF 1450 1211 CALL production_e( i, j ) … … 1462 1223 !-- TKE value is reduced by 90%. 1463 1224 DO k = nzb_s_inner(j,i)+1, nzt 1464 e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +& 1465 dt_3d * ( & 1466 tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) & 1467 ) 1225 e_p(k,j,i) = e(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 1226 tsc(3) * te_m(k,j,i) ) 1468 1227 IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i) 1469 1228 ENDDO … … 1507 1266 CHARACTER (LEN=9) :: time_to_string 1508 1267 INTEGER :: i, j, k 1509 REAL :: s at, sbt1268 REAL :: sbt 1510 1269 1511 1270 ! … … 1522 1281 CALL cpu_log( log_point(5), 'u-equation', 'start' ) 1523 1282 1524 ! 1525 !-- u-tendency terms with communication 1526 IF ( momentum_advec == 'ups-scheme' ) THEN 1527 tend = 0.0 1528 CALL advec_u_ups 1529 ENDIF 1530 1531 ! 1532 !-- u-tendency terms with no communication 1533 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 1534 tend = 0.0 1283 tend = 0.0 1284 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1535 1285 IF ( ws_scheme_mom ) THEN 1536 1286 CALL advec_u_ws … … 1539 1289 ENDIF 1540 1290 ELSE 1541 IF ( momentum_advec /= 'ups-scheme' ) THEN 1542 tend = 0.0 1543 CALL advec_u_up 1544 ENDIF 1291 CALL advec_u_up 1545 1292 ENDIF 1546 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 1547 CALL diffusion_u( ddzu, ddzw, km_m, tend, u_m, usws_m, uswst_m, & 1548 v_m, w_m ) 1549 ELSE 1550 CALL diffusion_u( ddzu, ddzw, km, tend, u, usws, uswst, v, w ) 1551 ENDIF 1293 CALL diffusion_u 1552 1294 CALL coriolis( 1 ) 1553 1295 IF ( sloping_surface .AND. .NOT. neutral ) THEN … … 1578 1320 DO j = nys, nyn 1579 1321 DO k = nzb_u_inner(j,i)+1, nzt 1580 u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + & 1581 dt_3d * ( & 1582 tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) & 1583 ) - & 1584 tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) ) 1322 u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 1323 tsc(3) * tu_m(k,j,i) ) & 1324 - tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) ) 1585 1325 ENDDO 1586 1326 ENDDO … … 1616 1356 CALL cpu_log( log_point(6), 'v-equation', 'start' ) 1617 1357 1618 ! 1619 !-- v-tendency terms with communication 1620 IF ( momentum_advec == 'ups-scheme' ) THEN 1621 tend = 0.0 1622 CALL advec_v_ups 1623 ENDIF 1624 1625 ! 1626 !-- v-tendency terms with no communication 1627 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 1628 tend = 0.0 1358 tend = 0.0 1359 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1629 1360 IF ( ws_scheme_mom ) THEN 1630 1361 CALL advec_v_ws … … 1633 1364 END IF 1634 1365 ELSE 1635 IF ( momentum_advec /= 'ups-scheme' ) THEN 1636 tend = 0.0 1637 CALL advec_v_up 1638 ENDIF 1366 CALL advec_v_up 1639 1367 ENDIF 1640 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 1641 CALL diffusion_v( ddzu, ddzw, km_m, tend, u_m, v_m, vsws_m, vswst_m, & 1642 w_m ) 1643 ELSE 1644 CALL diffusion_v( ddzu, ddzw, km, tend, u, v, vsws, vswst, w ) 1645 ENDIF 1368 CALL diffusion_v 1646 1369 CALL coriolis( 2 ) 1647 1370 … … 1669 1392 DO j = nysv, nyn 1670 1393 DO k = nzb_v_inner(j,i)+1, nzt 1671 v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + & 1672 dt_3d * ( & 1673 tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) & 1674 ) - & 1675 tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) ) 1394 v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 1395 tsc(3) * tv_m(k,j,i) ) & 1396 - tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) ) 1676 1397 ENDDO 1677 1398 ENDDO … … 1707 1428 CALL cpu_log( log_point(7), 'w-equation', 'start' ) 1708 1429 1709 ! 1710 !-- w-tendency terms with communication 1711 IF ( momentum_advec == 'ups-scheme' ) THEN 1712 tend = 0.0 1713 CALL advec_w_ups 1714 ENDIF 1715 1716 ! 1717 !-- w-tendency terms with no communication 1718 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 1719 tend = 0.0 1430 tend = 0.0 1431 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1720 1432 IF ( ws_scheme_mom ) THEN 1721 1433 CALL advec_w_ws … … 1724 1436 ENDIF 1725 1437 ELSE 1726 IF ( momentum_advec /= 'ups-scheme' ) THEN 1727 tend = 0.0 1728 CALL advec_w_up 1729 ENDIF 1438 CALL advec_w_up 1730 1439 ENDIF 1731 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 1732 CALL diffusion_w( ddzu, ddzw, km_m, tend, u_m, v_m, w_m ) 1733 ELSE 1734 CALL diffusion_w( ddzu, ddzw, km, tend, u, v, w ) 1735 ENDIF 1440 CALL diffusion_w 1736 1441 CALL coriolis( 3 ) 1737 1442 … … 1759 1464 DO j = nys, nyn 1760 1465 DO k = nzb_w_inner(j,i)+1, nzt-1 1761 w_p(k,j,i) = ( 1-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + & 1762 dt_3d * ( & 1763 tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) & 1764 ) - & 1765 tsc(5) * rdf(k) * w(k,j,i) 1466 w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + & 1467 tsc(3) * tw_m(k,j,i) ) & 1468 - tsc(5) * rdf(k) * w(k,j,i) 1766 1469 ENDDO 1767 1470 ENDDO … … 1802 1505 ! 1803 1506 !-- pt-tendency terms with communication 1804 sat = tsc(1)1805 1507 sbt = tsc(2) 1806 1508 IF ( scalar_advec == 'bc-scheme' ) THEN … … 1808 1510 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 1809 1511 ! 1810 !-- Bott-Chlond scheme always uses Euler time step when leapfrog is 1811 !-- switched on. Thus: 1812 sat = 1.0 1512 !-- Bott-Chlond scheme always uses Euler time step. Thus: 1813 1513 sbt = 1.0 1814 1514 ENDIF 1815 1515 tend = 0.0 1816 1516 CALL advec_s_bc( pt, 'pt' ) 1817 ELSE 1818 IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN 1819 tend = 0.0 1820 CALL advec_s_ups( pt, 'pt' ) 1821 ENDIF 1517 1822 1518 ENDIF 1823 1519 1824 1520 ! 1825 1521 !-- pt-tendency terms with no communication 1826 IF ( scalar_advec == 'bc-scheme' ) THEN 1827 CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, & 1828 tend ) 1829 ELSE 1830 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 1831 tend = 0.0 1522 IF ( scalar_advec /= 'bc-scheme' ) THEN 1523 tend = 0.0 1524 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1832 1525 IF ( ws_scheme_sca ) THEN 1833 1526 CALL advec_s_ws( pt, 'pt' ) … … 1836 1529 ENDIF 1837 1530 ELSE 1838 IF ( scalar_advec /= 'ups-scheme' ) THEN 1839 tend = 0.0 1840 CALL advec_s_up( pt ) 1841 ENDIF 1842 ENDIF 1843 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 1844 CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, & 1845 wall_heatflux, tend ) 1846 ELSE 1847 CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, & 1848 tend ) 1849 ENDIF 1850 ENDIF 1531 CALL advec_s_up( pt ) 1532 ENDIF 1533 ENDIF 1534 1535 CALL diffusion_s( pt, shf, tswst, wall_heatflux ) 1851 1536 1852 1537 ! … … 1881 1566 DO j = nys, nyn 1882 1567 DO k = nzb_s_inner(j,i)+1, nzt 1883 pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + & 1884 dt_3d * ( & 1885 sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i)& 1886 ) - & 1887 tsc(5) * ( pt(k,j,i) - pt_init(k) ) * & 1888 ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) ) 1568 pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 1569 tsc(3) * tpt_m(k,j,i) ) & 1570 - tsc(5) * ( pt(k,j,i) - pt_init(k) ) *& 1571 ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) ) 1889 1572 ENDDO 1890 1573 ENDDO … … 1927 1610 ! 1928 1611 !-- sa-tendency terms with communication 1929 sat = tsc(1)1930 1612 sbt = tsc(2) 1931 1613 IF ( scalar_advec == 'bc-scheme' ) THEN … … 1933 1615 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 1934 1616 ! 1935 !-- Bott-Chlond scheme always uses Euler time step when leapfrog is 1936 !-- switched on. Thus: 1937 sat = 1.0 1617 !-- Bott-Chlond scheme always uses Euler time step. Thus: 1938 1618 sbt = 1.0 1939 1619 ENDIF 1940 1620 tend = 0.0 1941 1621 CALL advec_s_bc( sa, 'sa' ) 1942 ELSE 1943 IF ( tsc(2) /= 2.0 ) THEN 1944 IF ( scalar_advec == 'ups-scheme' ) THEN 1945 tend = 0.0 1946 CALL advec_s_ups( sa, 'sa' ) 1947 ENDIF 1948 ENDIF 1622 1949 1623 ENDIF 1950 1624 1951 1625 ! 1952 1626 !-- sa-tendency terms with no communication 1953 IF ( scalar_advec == 'bc-scheme' ) THEN 1954 CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, & 1955 wall_salinityflux, tend ) 1956 ELSE 1957 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 1958 tend = 0.0 1627 IF ( scalar_advec /= 'bc-scheme' ) THEN 1628 tend = 0.0 1629 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1959 1630 IF ( ws_scheme_sca ) THEN 1960 1631 CALL advec_s_ws( sa, 'sa' ) … … 1963 1634 ENDIF 1964 1635 ELSE 1965 IF ( scalar_advec /= 'ups-scheme' ) THEN 1966 tend = 0.0 1967 CALL advec_s_up( sa ) 1968 ENDIF 1969 ENDIF 1970 CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, & 1971 wall_salinityflux, tend ) 1972 ENDIF 1636 CALL advec_s_up( sa ) 1637 ENDIF 1638 ENDIF 1639 1640 CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux ) 1973 1641 1974 1642 CALL user_actions( 'sa-tendency' ) … … 1979 1647 DO j = nys, nyn 1980 1648 DO k = nzb_s_inner(j,i)+1, nzt 1981 sa_p(k,j,i) = sat * sa(k,j,i) + & 1982 dt_3d * ( & 1983 sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) & 1984 ) - & 1985 tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) ) 1649 sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 1650 tsc(3) * tsa_m(k,j,i) ) & 1651 - tsc(5) * rdf_sc(k) * & 1652 ( sa(k,j,i) - sa_init(k) ) 1986 1653 IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i) 1987 1654 ENDDO … … 2031 1698 ! 2032 1699 !-- Scalar/q-tendency terms with communication 2033 sat = tsc(1)2034 1700 sbt = tsc(2) 2035 1701 IF ( scalar_advec == 'bc-scheme' ) THEN … … 2037 1703 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 2038 1704 ! 2039 !-- Bott-Chlond scheme always uses Euler time step when leapfrog is 2040 !-- switched on. Thus: 2041 sat = 1.0 1705 !-- Bott-Chlond scheme always uses Euler time step. Thus: 2042 1706 sbt = 1.0 2043 1707 ENDIF 2044 1708 tend = 0.0 2045 1709 CALL advec_s_bc( q, 'q' ) 2046 ELSE 2047 IF ( tsc(2) /= 2.0 ) THEN 2048 IF ( scalar_advec == 'ups-scheme' ) THEN 2049 tend = 0.0 2050 CALL advec_s_ups( q, 'q' ) 2051 ENDIF 2052 ENDIF 1710 2053 1711 ENDIF 2054 1712 2055 1713 ! 2056 1714 !-- Scalar/q-tendency terms with no communication 2057 IF ( scalar_advec == 'bc-scheme' ) THEN 2058 CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, wall_qflux, tend ) 2059 ELSE 2060 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 2061 tend = 0.0 1715 IF ( scalar_advec /= 'bc-scheme' ) THEN 1716 tend = 0.0 1717 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2062 1718 IF ( ws_scheme_sca ) THEN 2063 1719 CALL advec_s_ws( q, 'q' ) … … 2066 1722 ENDIF 2067 1723 ELSE 2068 IF ( scalar_advec /= 'ups-scheme' ) THEN 2069 tend = 0.0 2070 CALL advec_s_up( q ) 2071 ENDIF 2072 ENDIF 2073 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 2074 CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, qswst_m, & 2075 wall_qflux, tend ) 2076 ELSE 2077 CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, & 2078 wall_qflux, tend ) 2079 ENDIF 2080 ENDIF 1724 CALL advec_s_up( q ) 1725 ENDIF 1726 ENDIF 1727 1728 CALL diffusion_s( q, qsws, qswst, wall_qflux ) 2081 1729 2082 1730 ! … … 2104 1752 DO j = nys, nyn 2105 1753 DO k = nzb_s_inner(j,i)+1, nzt 2106 q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + & 2107 dt_3d * ( & 2108 sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) & 2109 ) - & 2110 tsc(5) * rdf_sc(k) * ( q(k,j,i) - q_init(k) ) 1754 q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 1755 tsc(3) * tq_m(k,j,i) ) & 1756 - tsc(5) * rdf_sc(k) * & 1757 ( q(k,j,i) - q_init(k) ) 2111 1758 IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i) 2112 1759 ENDDO … … 2152 1799 CALL production_e_init 2153 1800 2154 sat = tsc(1)2155 1801 sbt = tsc(2) 2156 1802 IF ( .NOT. use_upstream_for_tke ) THEN … … 2159 1805 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 2160 1806 ! 2161 !-- Bott-Chlond scheme always uses Euler time step when leapfrog is 2162 !-- switched on. Thus: 2163 sat = 1.0 1807 !-- Bott-Chlond scheme always uses Euler time step. Thus: 2164 1808 sbt = 1.0 2165 1809 ENDIF 2166 1810 tend = 0.0 2167 1811 CALL advec_s_bc( e, 'e' ) 2168 ELSE 2169 IF ( tsc(2) /= 2.0 ) THEN 2170 IF ( scalar_advec == 'ups-scheme' ) THEN 2171 tend = 0.0 2172 CALL advec_s_ups( e, 'e' ) 2173 ENDIF 2174 ENDIF 1812 2175 1813 ENDIF 2176 1814 ENDIF … … 2178 1816 ! 2179 1817 !-- TKE-tendency terms with no communication 2180 IF ( scalar_advec == 'bc-scheme' .AND. .NOT. use_upstream_for_tke ) & 2181 THEN 2182 IF ( .NOT. humidity ) THEN 2183 IF ( ocean ) THEN 2184 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, & 2185 prho, prho_reference, rif, tend, zu, zw ) 2186 ELSE 2187 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, & 2188 pt_reference, rif, tend, zu, zw ) 2189 ENDIF 2190 ELSE 2191 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, & 2192 pt_reference, rif, tend, zu, zw ) 2193 ENDIF 2194 ELSE 1818 IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN 2195 1819 IF ( use_upstream_for_tke ) THEN 2196 1820 tend = 0.0 2197 1821 CALL advec_s_up( e ) 2198 1822 ELSE 2199 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN2200 tend = 0.01823 tend = 0.0 1824 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2201 1825 IF ( ws_scheme_sca ) THEN 2202 1826 CALL advec_s_ws( e, 'e' ) … … 2205 1829 ENDIF 2206 1830 ELSE 2207 IF ( scalar_advec /= 'ups-scheme' ) THEN 2208 tend = 0.0 2209 CALL advec_s_up( e ) 2210 ENDIF 2211 ENDIF 2212 ENDIF 2213 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 2214 IF ( .NOT. humidity ) THEN 2215 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, & 2216 pt_m, pt_reference, rif_m, tend, zu, zw ) 2217 ELSE 2218 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, & 2219 vpt_m, pt_reference, rif_m, tend, zu, zw ) 2220 ENDIF 1831 CALL advec_s_up( e ) 1832 ENDIF 1833 ENDIF 1834 ENDIF 1835 1836 IF ( .NOT. humidity ) THEN 1837 IF ( ocean ) THEN 1838 CALL diffusion_e( prho, prho_reference ) 2221 1839 ELSE 2222 IF ( .NOT. humidity ) THEN 2223 IF ( ocean ) THEN 2224 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, & 2225 prho, prho_reference, rif, tend, zu, zw ) 2226 ELSE 2227 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, & 2228 pt, pt_reference, rif, tend, zu, zw ) 2229 ENDIF 2230 ELSE 2231 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, & 2232 pt_reference, rif, tend, zu, zw ) 2233 ENDIF 2234 ENDIF 2235 ENDIF 1840 CALL diffusion_e( pt, pt_reference ) 1841 ENDIF 1842 ELSE 1843 CALL diffusion_e( vpt, pt_reference ) 1844 ENDIF 1845 2236 1846 CALL production_e 2237 1847 … … 2249 1859 DO j = nys, nyn 2250 1860 DO k = nzb_s_inner(j,i)+1, nzt 2251 e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + & 2252 dt_3d * ( & 2253 sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) & 2254 ) 1861 e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 1862 tsc(3) * te_m(k,j,i) ) 2255 1863 IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i) 2256 1864 ENDDO -
palm/trunk/SOURCE/read_3d_binary.f90
r979 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! all actions concerning leapfrog scheme removed 7 7 ! 8 8 ! Former revisions: … … 376 376 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 377 377 378 CASE ( 'e_m' )379 IF ( k == 1 ) READ ( 13 ) tmp_3d380 e_m(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &381 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)382 383 378 CASE ( 'iran' ) ! matching random numbers is still unresolved 384 379 ! issue … … 390 385 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 391 386 392 CASE ( 'kh_m' )393 IF ( k == 1 ) READ ( 13 ) tmp_3d394 kh_m(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &395 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)396 397 387 CASE ( 'km' ) 398 388 IF ( k == 1 ) READ ( 13 ) tmp_3d 399 389 km(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 400 390 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 401 402 CASE ( 'km_m' )403 IF ( k == 1 ) READ ( 13 ) tmp_3d404 km_m(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &405 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)406 391 407 392 CASE ( 'lpt_av' ) … … 476 461 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 477 462 478 CASE ( 'pt_m' )479 IF ( k == 1 ) READ ( 13 ) tmp_3d480 pt_m(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &481 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)482 483 463 CASE ( 'q' ) 484 464 IF ( k == 1 ) READ ( 13 ) tmp_3d … … 494 474 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 495 475 496 CASE ( 'q_m' )497 IF ( k == 1 ) READ ( 13 ) tmp_3d498 q_m(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &499 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)500 501 476 CASE ( 'ql' ) 502 477 IF ( k == 1 ) READ ( 13 ) tmp_3d … … 544 519 IF ( k == 1 ) READ ( 13 ) tmp_2d 545 520 qsws(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 546 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)547 548 CASE ( 'qsws_m' )549 IF ( k == 1 ) READ ( 13 ) tmp_2d550 qsws_m(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &551 521 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 552 522 … … 564 534 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 565 535 566 CASE ( 'qswst_m' )567 IF ( k == 1 ) READ ( 13 ) tmp_2d568 qswst_m(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &569 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)570 571 536 CASE ( 'qv_av' ) 572 537 IF ( .NOT. ALLOCATED( qv_av ) ) THEN … … 594 559 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 595 560 596 CASE ( 'rif_m' )597 IF ( k == 1 ) READ ( 13 ) tmp_2d598 rif_m(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &599 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)600 601 561 CASE ( 'rif_wall' ) 602 562 IF ( k == 1 ) THEN … … 644 604 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 645 605 646 CASE ( 'shf_m' )647 IF ( k == 1 ) READ ( 13 ) tmp_2d648 shf_m(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &649 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)650 606 CASE ( 'shf_av' ) 651 607 IF ( .NOT. ALLOCATED( shf_av ) ) THEN … … 655 611 shf_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 656 612 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 613 657 614 CASE ( 'spectrum_x' ) 658 615 IF ( k == 1 ) THEN … … 703 660 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 704 661 705 CASE ( 'tswst_m' )706 IF ( k == 1 ) READ ( 13 ) tmp_2d707 tswst_m(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &708 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)709 710 662 CASE ( 'u' ) 711 663 IF ( k == 1 ) READ ( 13 ) tmp_3d … … 721 673 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 722 674 723 CASE ( 'u_m' )724 IF ( k == 1 ) READ ( 13 ) tmp_3d725 u_m(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &726 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)727 728 675 CASE ( 'u_m_l' ) 729 676 IF ( k == 1 ) THEN … … 781 728 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 782 729 783 CASE ( 'usws_m' )784 IF ( k == 1 ) READ ( 13 ) tmp_2d785 usws_m(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &786 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)787 788 CASE ( 'uswst_m' )789 IF ( k == 1 ) READ ( 13 ) tmp_2d790 uswst_m(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &791 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)792 793 730 CASE ( 'us_av' ) 794 731 IF ( .NOT. ALLOCATED( us_av ) ) THEN … … 812 749 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 813 750 814 CASE ( 'v_m' )815 IF ( k == 1 ) READ ( 13 ) tmp_3d816 v_m(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &817 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)818 819 751 CASE ( 'v_m_l' ) 820 752 IF ( k == 1 ) THEN … … 870 802 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 871 803 872 CASE ( 'vpt_m' )873 IF ( k == 1 ) READ ( 13 ) tmp_3d874 vpt_m(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &875 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)876 877 804 CASE ( 'vsws' ) 878 805 IF ( k == 1 ) READ ( 13 ) tmp_2d … … 884 811 vswst(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 885 812 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 886 887 CASE ( 'vsws_m' )888 IF ( k == 1 ) READ ( 13 ) tmp_2d889 vsws_m(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &890 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)891 892 CASE ( 'vswst_m' )893 IF ( k == 1 ) READ ( 13 ) tmp_2d894 vswst_m(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &895 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)896 813 897 814 CASE ( 'w' ) … … 907 824 w_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 908 825 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 909 910 CASE ( 'w_m' )911 IF ( k == 1 ) READ ( 13 ) tmp_3d912 w_m(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &913 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)914 826 915 827 CASE ( 'w_m_l' ) -
palm/trunk/SOURCE/read_var_list.f90
r979 r1001 4 4 ! Current revisions: 5 5 ! ------------------ 6 ! 6 ! -cut_spline_overshoot, dt_fixed, last_dt_change, long_filter_factor, 7 ! overshoot_limit_*, ups_limit_* 7 8 ! 8 9 ! Former revisions: … … 313 314 CASE ( 'curvature_solution_effects' ) 314 315 READ ( 13 ) curvature_solution_effects 315 CASE ( 'cut_spline_overshoot' )316 READ ( 13 ) cut_spline_overshoot317 316 CASE ( 'cycle_mg' ) 318 317 READ ( 13 ) cycle_mg … … 331 330 CASE ( 'drag_coefficient' ) 332 331 READ ( 13 ) drag_coefficient 333 CASE ( 'dt_fixed' )334 READ ( 13 ) ldum ! restart files created before rev 333335 ! contained dt_fixed by mistake; it is still336 ! read here in order to allow usage of these337 ! older restart files; can be removed in a338 ! later version339 332 CASE ( 'dt_pr_1d' ) 340 333 READ ( 13 ) dt_pr_1d … … 400 393 CASE ( 'lad_vertical_gradient_level_in' ) 401 394 READ ( 13 ) lad_vertical_gradient_level_ind 402 CASE ( 'last_dt_change' )403 READ ( 13 ) last_dt_change404 395 CASE ( 'large_scale_subsidence' ) 405 396 READ ( 13 ) large_scale_subsidence 406 397 CASE ( 'leaf_surface_concentration' ) 407 398 READ ( 13 ) leaf_surface_concentration 408 CASE ( 'long_filter_factor' )409 READ ( 13 ) long_filter_factor410 399 CASE ( 'loop_optimization' ) 411 400 READ ( 13 ) loop_optimization … … 451 440 CASE ( 'output_for_t0' ) 452 441 READ (13) output_for_t0 453 CASE ( 'overshoot_limit_e' )454 READ ( 13 ) overshoot_limit_e455 CASE ( 'overshoot_limit_pt' )456 READ ( 13 ) overshoot_limit_pt457 CASE ( 'overshoot_limit_u' )458 READ ( 13 ) overshoot_limit_u459 CASE ( 'overshoot_limit_v' )460 READ ( 13 ) overshoot_limit_v461 CASE ( 'overshoot_limit_w' )462 READ ( 13 ) overshoot_limit_w463 442 CASE ( 'passive_scalar' ) 464 443 READ ( 13 ) passive_scalar … … 637 616 CASE ( 'ug_vertical_gradient_level_ind' ) 638 617 READ ( 13 ) ug_vertical_gradient_level_ind 639 CASE ( 'ups_limit_e' )640 READ ( 13 ) ups_limit_e641 CASE ( 'ups_limit_pt' )642 READ ( 13 ) ups_limit_pt643 CASE ( 'ups_limit_u' )644 READ ( 13 ) ups_limit_u645 CASE ( 'ups_limit_v' )646 READ ( 13 ) ups_limit_v647 CASE ( 'ups_limit_w' )648 READ ( 13 ) ups_limit_w649 618 CASE ( 'use_surface_fluxes' ) 650 619 READ ( 13 ) use_surface_fluxes -
palm/trunk/SOURCE/run_control.f90
r484 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! all actions concerning leapfrog scheme removed 7 7 ! 8 8 ! Former revisions: … … 43 43 IMPLICIT NONE 44 44 45 CHARACTER (LEN=1) :: change_chr,disturb_chr45 CHARACTER (LEN=1) :: disturb_chr 46 46 47 47 ! … … 70 70 71 71 ! 72 !-- Output the the beginning of the run receives no information about an73 !-- Euler-timestep74 IF ( dt_changed .AND. simulated_time /= 0.0 .AND. &75 timestep_scheme(1:5) /= 'runge' ) THEN76 IF ( timestep_scheme == 'leapfrog' ) THEN77 change_chr = 'L'78 ELSE79 change_chr = 'E'80 ENDIF81 ELSE82 change_chr = ' '83 ENDIF84 !85 72 !-- If required, set disturbance flag 86 73 IF ( disturbance_created ) THEN … … 91 78 WRITE ( 15, 101 ) runnr, current_timestep_number, simulated_time_chr, & 92 79 simulated_time-INT( simulated_time ), dt_3d, & 93 timestep_reason, change_chr, u_max, disturb_chr,&80 timestep_reason, u_max, disturb_chr, & 94 81 v_max, disturb_chr, w_max, hom(nzb,1,pr_palm,0), & 95 82 hom(nzb+8,1,pr_palm,0), hom(nzb+3,1,pr_palm,0), & … … 124 111 &'----------------------------------------------------------------', & 125 112 &'-----') 126 101 FORMAT (I3,1X,I6,1X,A8,F3.2,1X,F8.4,A1, A1,F8.4,A1,F8.4,A1,F8.4,2X,F5.3,2X, &113 101 FORMAT (I3,1X,I6,1X,A8,F3.2,1X,F8.4,A1,1X,F8.4,A1,F8.4,A1,F8.4,2X,F5.3,2X, & 127 114 F4.2, & 128 115 2X,F6.3,2X,F6.0,1X,4(E10.3,1X),3(3(I4),1X),F8.3,1X,F8.3,5X,I3) -
palm/trunk/SOURCE/swap_timelevel.f90
r484 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! all actions concerning leapfrog scheme removed 7 7 ! 8 8 ! Former revisions: … … 51 51 52 52 ! 53 !-- Swap of 3-level variables 54 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 55 56 SELECT CASE ( MOD( timestep_count, 3 ) ) 57 58 CASE ( 0 ) 59 60 u_m => u_1; u => u_2; u_p => u_3 61 v_m => v_1; v => v_2; v_p => v_3 62 w_m => w_1; w => w_2; w_p => w_3 63 pt_m => pt_1; pt => pt_2; pt_p => pt_3 64 IF ( .NOT. constant_diffusion ) THEN 65 e_m => e_1; e => e_2; e_p => e_3 66 ENDIF 67 IF ( humidity .OR. passive_scalar ) THEN 68 q_m => q_1; q => q_2; q_p => q_3 69 ENDIF 70 71 CASE ( 1 ) 72 73 u_m => u_2; u => u_3; u_p => u_1 74 v_m => v_2; v => v_3; v_p => v_1 75 w_m => w_2; w => w_3; w_p => w_1 76 pt_m => pt_2; pt => pt_3; pt_p => pt_1 77 IF ( .NOT. constant_diffusion ) THEN 78 e_m => e_2; e => e_3; e_p => e_1 79 ENDIF 80 IF ( humidity .OR. passive_scalar ) THEN 81 q_m => q_2; q => q_3; q_p => q_1 82 ENDIF 83 84 CASE ( 2 ) 85 86 u_m => u_3; u => u_1; u_p => u_2 87 v_m => v_3; v => v_1; v_p => v_2 88 w_m => w_3; w => w_1; w_p => w_2 89 pt_m => pt_3; pt => pt_1; pt_p => pt_2 90 IF ( .NOT. constant_diffusion ) THEN 91 e_m => e_3; e => e_1; e_p => e_2 92 ENDIF 93 IF ( humidity .OR. passive_scalar ) THEN 94 q_m => q_3; q => q_1; q_p => q_2 95 ENDIF 96 97 END SELECT 98 99 ENDIF 100 101 ! 102 !-- Swap of 2-level variables 53 !-- Swap of variables 103 54 SELECT CASE ( MOD( timestep_count, 2 ) ) 104 55 105 56 CASE ( 0 ) 106 57 107 IF ( timestep_scheme(1:5) == 'runge' ) THEN 58 u => u_1; u_p => u_2 59 v => v_1; v_p => v_2 60 w => w_1; w_p => w_2 61 pt => pt_1; pt_p => pt_2 62 IF ( .NOT. constant_diffusion ) THEN 63 e => e_1; e_p => e_2 64 ENDIF 65 IF ( ocean ) THEN 66 sa => sa_1; sa_p => sa_2 67 ENDIF 68 IF ( humidity .OR. passive_scalar ) THEN 69 q => q_1; q_p => q_2 70 ENDIF 108 71 109 u => u_1; u_p => u_2110 v => v_1; v_p => v_2111 w => w_1; w_p => w_2112 pt => pt_1; pt_p => pt_2113 IF ( .NOT. constant_diffusion ) THEN114 e => e_1; e_p => e_2115 ENDIF116 IF ( ocean ) THEN117 sa => sa_1; sa_p => sa_2118 ENDIF119 IF ( humidity .OR. passive_scalar ) THEN120 q => q_1; q_p => q_2121 ENDIF122 123 ELSE124 !125 !-- Old timelevels are needed for explicit diffusion within leapfrog126 IF ( .NOT. constant_diffusion ) THEN127 kh_m => kh_1; kh => kh_2128 km_m => km_1; km => km_2129 IF ( use_surface_fluxes ) THEN130 usws_m => usws_1; usws => usws_2131 vsws_m => vsws_1; vsws => vsws_2132 shf_m => shf_1; shf => shf_2133 IF ( humidity .OR. passive_scalar ) THEN134 qsws_m => qsws_1; qsws => qsws_2135 ENDIF136 ENDIF137 IF ( prandtl_layer ) THEN138 rif_m => rif_1; rif => rif_2139 ENDIF140 IF ( use_top_fluxes ) THEN141 uswst_m => uswst_1; uswst => uswst_2142 vswst_m => vswst_1; vswst => vswst_2143 tswst_m => tswst_1; tswst => tswst_2144 IF ( humidity .OR. passive_scalar ) THEN145 qswst_m => qswst_1; qswst => qswst_2146 ENDIF147 ENDIF148 ENDIF149 150 IF ( humidity ) THEN151 vpt_m => vpt_1; vpt => vpt_2152 ENDIF153 154 ENDIF155 72 156 73 CASE ( 1 ) 157 74 158 IF ( timestep_scheme(1:5) == 'runge' ) THEN 75 u => u_2; u_p => u_1 76 v => v_2; v_p => v_1 77 w => w_2; w_p => w_1 78 pt => pt_2; pt_p => pt_1 79 IF ( .NOT. constant_diffusion ) THEN 80 e => e_2; e_p => e_1 81 ENDIF 82 IF ( ocean ) THEN 83 sa => sa_2; sa_p => sa_1 84 ENDIF 85 IF ( humidity .OR. passive_scalar ) THEN 86 q => q_2; q_p => q_1 87 ENDIF 159 88 160 u => u_2; u_p => u_1161 v => v_2; v_p => v_1162 w => w_2; w_p => w_1163 pt => pt_2; pt_p => pt_1164 IF ( .NOT. constant_diffusion ) THEN165 e => e_2; e_p => e_1166 ENDIF167 IF ( ocean ) THEN168 sa => sa_2; sa_p => sa_1169 ENDIF170 IF ( humidity .OR. passive_scalar ) THEN171 q => q_2; q_p => q_1172 ENDIF173 174 ELSE175 176 IF ( .NOT. constant_diffusion ) THEN177 kh_m => kh_2; kh => kh_1178 km_m => km_2; km => km_1179 IF ( use_surface_fluxes ) THEN180 usws_m => usws_2; usws => usws_1181 vsws_m => vsws_2; vsws => vsws_1182 shf_m => shf_2; shf => shf_1183 IF ( humidity .OR. passive_scalar ) THEN184 qsws_m => qsws_2; qsws => qsws_1185 ENDIF186 ENDIF187 IF ( prandtl_layer ) THEN188 rif_m => rif_2; rif => rif_1189 ENDIF190 IF ( use_top_fluxes ) THEN191 uswst_m => uswst_2; uswst => uswst_1192 vswst_m => vswst_2; vswst => vswst_1193 tswst_m => tswst_2; tswst => tswst_1194 IF ( humidity .OR. passive_scalar ) THEN195 qswst_m => qswst_2; qswst => qswst_1196 ENDIF197 ENDIF198 ENDIF199 200 IF ( humidity ) THEN201 vpt_m => vpt_2; vpt => vpt_1202 ENDIF203 204 ENDIF205 89 206 90 END SELECT -
palm/trunk/SOURCE/time_integration.f90
r850 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! all actions concerning leapfrog- and upstream-spline-scheme removed 7 7 ! 8 8 ! Former revisions: … … 174 174 CALL prognostic_equations_vector 175 175 ELSE 176 IF ( momentum_advec == 'ups-scheme' .OR. & 177 scalar_advec == 'ups-scheme' .OR. & 178 scalar_advec == 'bc-scheme' ) & 179 THEN 176 IF ( scalar_advec == 'bc-scheme' ) THEN 180 177 CALL prognostic_equations_noopt 181 178 ELSE … … 227 224 228 225 CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' ) 229 230 !231 !-- Apply time filter in case of leap-frog timestep232 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN233 CALL asselin_filter234 ENDIF235 226 236 227 ! … … 459 450 ! 460 451 !-- Computation and output of run control parameters. 461 !-- This is also done whenever the time step has changed or perturbations 462 !-- have been imposed 452 !-- This is also done whenever perturbations have been imposed 463 453 IF ( time_run_control >= dt_run_control .OR. & 464 ( dt_changed .AND. timestep_scheme(1:5) /= 'runge' ) .OR. & 465 disturbance_created ) & 454 timestep_scheme(1:5) /= 'runge' .OR. disturbance_created ) & 466 455 THEN 467 456 CALL run_control -
palm/trunk/SOURCE/timestep.f90
r979 r1001 4 4 ! Current revisions: 5 5 ! ------------------ 6 ! 6 ! all actions concerning leapfrog scheme removed 7 7 ! 8 8 ! Former revisions: … … 68 68 dt_plant_canopy_l, & 69 69 dt_plant_canopy_u, dt_plant_canopy_v, dt_plant_canopy_w, & 70 dt_u, dt_v, dt_w, lad_max, percent_change,&70 dt_u, dt_v, dt_w, lad_max, & 71 71 u_gtrans_l, u_max_cfl, vabs_max, value, v_gtrans_l, v_max_cfl 72 72 … … 277 277 ! 278 278 !-- The time step is the minimum of the 3-4 components and the diffusion time 279 !-- step minus a reduction to be on the safe side. Factor 0.5 is necessary 280 !-- since the leap-frog scheme always progresses by 2 * delta t. 281 !-- The user has to set the cfl_factor small enough to ensure that the 282 !-- divergences do not become too large. 279 !-- step minus a reduction (cfl_factor) to be on the safe side. 283 280 !-- The time step must not exceed the maximum allowed value. 284 IF ( timestep_scheme(1:5) == 'runge' ) THEN 285 dt_3d = cfl_factor * MIN( dt_diff, dt_plant_canopy, dt_u, dt_v, dt_w ) 286 ELSE 287 dt_3d = cfl_factor * 0.5 * & 288 MIN( dt_diff, dt_plant_canopy, dt_u, dt_v, dt_w ) 289 ENDIF 281 dt_3d = cfl_factor * MIN( dt_diff, dt_plant_canopy, dt_u, dt_v, dt_w ) 290 282 dt_3d = MIN( dt_3d, dt_max ) 291 283 … … 339 331 340 332 ! 341 !-- Ensure a smooth value (two significant digits) of the timestep. For 342 !-- other schemes than Runge-Kutta, the following restrictions appear: 343 !-- The current timestep is only then changed, if the change relative to 344 !-- its previous value exceeds +5 % or -2 %. In case of a timestep 345 !-- reduction, at least 30 iterations have to be performed before a timestep 346 !-- enlargement is permitted again. 347 percent_change = dt_3d / old_dt - 1.0 348 IF ( percent_change > 0.05 .OR. percent_change < -0.02 .OR. & 349 timestep_scheme(1:5) == 'runge' ) THEN 350 351 ! 352 !-- Time step enlargement by no more than 2 %. 353 IF ( percent_change > 0.0 .AND. simulated_time /= 0.0 .AND. & 354 timestep_scheme(1:5) /= 'runge' ) THEN 355 dt_3d = 1.02 * old_dt 356 ENDIF 357 358 ! 359 !-- A relatively smooth value of the time step is ensured by taking 360 !-- only the first two significant digits. 361 div = 1000.0 362 DO WHILE ( dt_3d < div ) 363 div = div / 10.0 364 ENDDO 365 dt_3d = NINT( dt_3d * 100.0 / div ) * div / 100.0 366 367 ! 368 !-- Now the time step can be adjusted. 369 IF ( percent_change < 0.0 .OR. timestep_scheme(1:5) == 'runge' ) & 370 THEN 371 ! 372 !-- Time step reduction. 373 old_dt = dt_3d 374 dt_changed = .TRUE. 375 ELSE 376 ! 377 !-- For other timestep schemes , the time step is only enlarged 378 !-- after at least 30 iterations since the previous time step 379 !-- change or, of course, after model initialization. 380 IF ( current_timestep_number >= last_dt_change + 30 .OR. & 381 simulated_time == 0.0 ) THEN 382 old_dt = dt_3d 383 dt_changed = .TRUE. 384 ELSE 385 dt_3d = old_dt 386 dt_changed = .FALSE. 387 ENDIF 388 389 ENDIF 390 ELSE 391 ! 392 !-- No time step change since the difference is too small. 393 dt_3d = old_dt 394 dt_changed = .FALSE. 395 ENDIF 396 397 IF ( dt_changed ) last_dt_change = current_timestep_number 333 !-- Ensure a smooth value (two significant digits) of the timestep. 334 div = 1000.0 335 DO WHILE ( dt_3d < div ) 336 div = div / 10.0 337 ENDDO 338 dt_3d = NINT( dt_3d * 100.0 / div ) * div / 100.0 339 340 ! 341 !-- Adjust the time step 342 old_dt = dt_3d 398 343 399 344 ENDIF 345 400 346 CALL cpu_log( log_point(12), 'calculate_timestep', 'stop' ) 401 347 -
palm/trunk/SOURCE/timestep_scheme_steering.f90
r674 r1001 3 3 !------------------------------------------------------------------------------! 4 4 ! Current revisions: 5 ! ----------------- 5 ! ------------------ 6 ! all actions concerning leapfrog scheme removed 6 7 ! 7 8 ! Former revisions: … … 52 53 ENDIF 53 54 54 ELSE 55 56 IF ( .NOT. dt_fixed ) THEN 55 ELSEIF ( timestep_scheme == 'euler' ) THEN 57 56 ! 58 !-- Leapfrog and Euler schemes 59 !-- Determine whether after the time step adjustment the Euler- or the 60 !-- leapfrog scheme will be applied. The very first time step must always 61 !-- be an Euler step. 62 IF ( dt_changed ) THEN 63 IF ( timestep_scheme == 'leapfrog+euler' .OR. & 64 timestep_scheme == 'euler' .OR. simulated_time == 0.0 ) THEN 65 tsc(1:5) = (/ 1.0, 1.0, 0.0, 0.0, 1.0 /) 66 ELSE 67 tsc(1:5) = (/ 0.0, 2.0, 0.0, 0.0, 2.0 /) 68 ENDIF 69 ELSE 70 ! 71 !-- No time step change, hence continue with the scheme set by the 72 !-- user. 73 IF ( timestep_scheme == 'euler' ) THEN 74 tsc(1:5) = (/ 1.0, 1.0, 0.0, 0.0, 1.0 /) 75 ELSE 76 tsc(1:5) = (/ 0.0, 2.0, 0.0, 0.0, 2.0 /) 77 ENDIF 78 ENDIF 79 80 ELSE 81 82 ! 83 !-- Fixed time step: 84 ! 85 !-- In any case, the very first time step must always be an Euler step. 86 timestep_reason = 'F' 87 IF ( simulated_time == 0.0 ) THEN 88 dt_changed = .TRUE. 89 tsc(1:5) = (/ 1.0, 1.0, 0.0, 0.0, 1.0 /) 90 ELSE 91 dt_changed = .FALSE. 92 IF ( timestep_scheme == 'euler' ) THEN 93 tsc(1:5) = (/ 1.0, 1.0, 0.0, 0.0, 1.0 /) 94 ELSE 95 tsc(1:5) = (/ 0.0, 2.0, 0.0, 0.0, 2.0 /) 96 ENDIF 97 ENDIF 98 99 ENDIF 57 !-- Euler scheme 58 tsc(1:5) = (/ 1.0, 1.0, 0.0, 0.0, 1.0 /) 100 59 101 60 ENDIF -
palm/trunk/SOURCE/write_3d_binary.f90
r979 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! all actions concerning leapfrog scheme removed 7 7 ! 8 8 ! Former revisions: … … 104 104 WRITE ( 14 ) 'e_av '; WRITE ( 14 ) e_av 105 105 ENDIF 106 WRITE ( 14 ) 'e_m '; WRITE ( 14 ) e_m107 106 WRITE ( 14 ) 'iran '; WRITE ( 14 ) iran, iran_part 108 107 WRITE ( 14 ) 'kh '; WRITE ( 14 ) kh 109 IF ( timestep_scheme(1:5) /= 'runge' ) THEN110 WRITE ( 14 ) 'kh_m '; WRITE ( 14 ) kh_m111 ENDIF112 108 WRITE ( 14 ) 'km '; WRITE ( 14 ) km 113 IF ( timestep_scheme(1:5) /= 'runge' ) THEN114 WRITE ( 14 ) 'km_m '; WRITE ( 14 ) km_m115 ENDIF116 109 IF ( ALLOCATED( lpt_av ) ) THEN 117 110 WRITE ( 14 ) 'lpt_av '; WRITE ( 14 ) lpt_av … … 141 134 WRITE ( 14 ) 'pt_av '; WRITE ( 14 ) pt_av 142 135 ENDIF 143 WRITE ( 14 ) 'pt_m '; WRITE ( 14 ) pt_m144 136 IF ( humidity .OR. passive_scalar ) THEN 145 137 WRITE ( 14 ) 'q '; WRITE ( 14 ) q … … 147 139 WRITE ( 14 ) 'q_av '; WRITE ( 14 ) q_av 148 140 ENDIF 149 WRITE ( 14 ) 'q_m '; WRITE ( 14 ) q_m150 141 IF ( cloud_physics .OR. cloud_droplets ) THEN 151 142 WRITE ( 14 ) 'ql '; WRITE ( 14 ) ql … … 156 147 WRITE ( 14 ) 'qs '; WRITE ( 14 ) qs 157 148 WRITE ( 14 ) 'qsws '; WRITE ( 14 ) qsws 158 IF ( timestep_scheme(1:5) /= 'runge' ) THEN159 WRITE ( 14 ) 'qsws_m '; WRITE ( 14 ) qsws_m160 ENDIF161 149 IF ( ALLOCATED( qsws_av ) ) THEN 162 150 WRITE ( 14 ) 'qsws_av '; WRITE ( 14 ) qsws_av 163 151 ENDIF 164 152 WRITE ( 14 ) 'qswst '; WRITE ( 14 ) qswst 165 IF ( timestep_scheme(1:5) /= 'runge' ) THEN166 WRITE ( 14 ) 'qswst_m '; WRITE ( 14 ) qswst_m167 ENDIF168 153 ENDIF 169 154 IF ( ocean ) THEN … … 193 178 WRITE ( 14 ) random_iy 194 179 WRITE ( 14 ) 'rif '; WRITE ( 14 ) rif 195 IF ( timestep_scheme(1:5) /= 'runge' ) THEN196 WRITE ( 14 ) 'rif_m '; WRITE ( 14 ) rif_m197 ENDIF198 180 IF ( topography /= 'flat' ) THEN 199 181 WRITE ( 14 ) 'rif_wall '; WRITE ( 14 ) rif_wall … … 203 185 ENDIF 204 186 WRITE ( 14 ) 'shf '; WRITE ( 14 ) shf 205 IF ( timestep_scheme(1:5) /= 'runge' ) THEN206 WRITE ( 14 ) 'shf_m '; WRITE ( 14 ) shf_m207 ENDIF208 187 IF ( ALLOCATED( shf_av ) ) THEN 209 188 WRITE ( 14 ) 'shf_av '; WRITE ( 14 ) shf_av … … 218 197 ENDIF 219 198 WRITE ( 14 ) 'tswst '; WRITE ( 14 ) tswst 220 IF ( timestep_scheme(1:5) /= 'runge' ) THEN221 WRITE ( 14 ) 'tswst_m '; WRITE ( 14 ) tswst_m222 ENDIF223 199 WRITE ( 14 ) 'u '; WRITE ( 14 ) u 224 200 IF ( ALLOCATED( u_av ) ) THEN 225 201 WRITE ( 14 ) 'u_av '; WRITE ( 14 ) u_av 226 202 ENDIF 227 WRITE ( 14 ) 'u_m '; WRITE ( 14 ) u_m228 203 IF ( ALLOCATED( u_m_l ) ) THEN 229 204 WRITE ( 14 ) 'u_m_l '; WRITE ( 14 ) u_m_l … … 241 216 WRITE ( 14 ) 'usws '; WRITE ( 14 ) usws 242 217 WRITE ( 14 ) 'uswst '; WRITE ( 14 ) uswst 243 IF ( timestep_scheme(1:5) /= 'runge' ) THEN244 WRITE ( 14 ) 'usws_m '; WRITE ( 14 ) usws_m245 WRITE ( 14 ) 'uswst_m '; WRITE ( 14 ) uswst_m246 ENDIF247 218 IF ( ALLOCATED( us_av ) ) THEN 248 219 WRITE ( 14 ) 'us_av '; WRITE ( 14 ) us_av … … 252 223 WRITE ( 14 ) 'v_av '; WRITE ( 14 ) v_av 253 224 ENDIF 254 WRITE ( 14 ) 'v_m '; WRITE ( 14 ) v_m255 225 IF ( ALLOCATED( v_m_l ) ) THEN 256 226 WRITE ( 14 ) 'v_m_l '; WRITE ( 14 ) v_m_l … … 270 240 WRITE ( 14 ) 'vpt_av '; WRITE ( 14 ) vpt_av 271 241 ENDIF 272 IF ( timestep_scheme(1:5) /= 'runge' ) THEN273 WRITE ( 14 ) 'vpt_m '; WRITE ( 14 ) vpt_m274 ENDIF275 242 ENDIF 276 243 WRITE ( 14 ) 'vsws '; WRITE ( 14 ) vsws 277 244 WRITE ( 14 ) 'vswst '; WRITE ( 14 ) vswst 278 IF ( timestep_scheme(1:5) /= 'runge' ) THEN279 WRITE ( 14 ) 'vsws_m '; WRITE ( 14 ) vsws_m280 WRITE ( 14 ) 'vswst_m '; WRITE ( 14 ) vswst_m281 ENDIF282 245 WRITE ( 14 ) 'w '; WRITE ( 14 ) w 283 246 IF ( ALLOCATED( w_av ) ) THEN 284 247 WRITE ( 14 ) 'w_av '; WRITE ( 14 ) w_av 285 248 ENDIF 286 WRITE ( 14 ) 'w_m '; WRITE ( 14 ) w_m287 249 IF ( ALLOCATED( w_m_l ) ) THEN 288 250 WRITE ( 14 ) 'w_m_l '; WRITE ( 14 ) w_m_l -
palm/trunk/SOURCE/write_var_list.f90
r979 r1001 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! -cut_spline_overshoot, last_dt_change, long_filter_factor, overshoot_limit_*, 7 ! ups_limit_* 7 8 ! 8 9 ! Former revisions: … … 235 236 WRITE ( 14 ) 'cthf ' 236 237 WRITE ( 14 ) cthf 237 WRITE ( 14 ) 'cut_spline_overshoot '238 WRITE ( 14 ) cut_spline_overshoot239 238 WRITE ( 14 ) 'cycle_mg ' 240 239 WRITE ( 14 ) cycle_mg … … 315 314 WRITE ( 14 ) 'lad_vertical_gradient_level_in' 316 315 WRITE ( 14 ) lad_vertical_gradient_level_ind 317 WRITE ( 14 ) 'last_dt_change '318 WRITE ( 14 ) last_dt_change319 316 WRITE ( 14 ) 'large_scale_subsidence ' 320 317 WRITE ( 14 ) large_scale_subsidence 321 318 WRITE ( 14 ) 'leaf_surface_concentration ' 322 319 WRITE ( 14 ) leaf_surface_concentration 323 WRITE ( 14 ) 'long_filter_factor '324 WRITE ( 14 ) long_filter_factor325 320 WRITE ( 14 ) 'loop_optimization ' 326 321 WRITE ( 14 ) loop_optimization … … 363 358 WRITE ( 14 ) 'output_for_t0 ' 364 359 WRITE ( 14 ) output_for_t0 365 WRITE ( 14 ) 'overshoot_limit_e '366 WRITE ( 14 ) overshoot_limit_e367 WRITE ( 14 ) 'overshoot_limit_pt '368 WRITE ( 14 ) overshoot_limit_pt369 WRITE ( 14 ) 'overshoot_limit_u '370 WRITE ( 14 ) overshoot_limit_u371 WRITE ( 14 ) 'overshoot_limit_v '372 WRITE ( 14 ) overshoot_limit_v373 WRITE ( 14 ) 'overshoot_limit_w '374 WRITE ( 14 ) overshoot_limit_w375 360 WRITE ( 14 ) 'passive_scalar ' 376 361 WRITE ( 14 ) passive_scalar … … 549 534 WRITE ( 14 ) 'ug_vertical_gradient_level_ind' 550 535 WRITE ( 14 ) ug_vertical_gradient_level_ind 551 WRITE ( 14 ) 'ups_limit_e '552 WRITE ( 14 ) ups_limit_e553 WRITE ( 14 ) 'ups_limit_pt '554 WRITE ( 14 ) ups_limit_pt555 WRITE ( 14 ) 'ups_limit_u '556 WRITE ( 14 ) ups_limit_u557 WRITE ( 14 ) 'ups_limit_v '558 WRITE ( 14 ) ups_limit_v559 WRITE ( 14 ) 'ups_limit_w '560 WRITE ( 14 ) ups_limit_w561 536 WRITE ( 14 ) 'use_surface_fluxes ' 562 537 WRITE ( 14 ) use_surface_fluxes
Note: See TracChangeset
for help on using the changeset viewer.