Changeset 2037 for palm/trunk/SOURCE/init_3d_model.f90
- Timestamp:
- Oct 26, 2016 11:15:40 AM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/init_3d_model.f90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/init_3d_model.f90
r2032 r2037 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Anelastic approximation implemented 23 23 ! 24 24 ! Former revisions: … … 315 315 USE arrays_3d 316 316 317 USE cloud_parameters, & 318 ONLY: cp, l_v, r_d 319 317 320 USE constants, & 318 321 ONLY: pi … … 324 327 325 328 USE grid_variables, & 326 ONLY: dx, dy 329 ONLY: dx, dy, ddx2_mg, ddy2_mg 327 330 328 331 USE indices … … 393 396 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: ngp_2dh_outer_l !< 394 397 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: ngp_2dh_s_inner_l !< 398 399 REAL(wp) :: t_surface !< air temperature at the surface 400 401 REAL(wp), DIMENSION(:), ALLOCATABLE :: p_hydrostatic !< hydrostatic pressure 402 403 INTEGER(iwp) :: l !< loop variable 404 INTEGER(iwp) :: nzt_l !< index of top PE boundary for multigrid level 405 REAL(wp) :: dx_l !< grid spacing along x on different multigrid level 406 REAL(wp) :: dy_l !< grid spacing along y on different multigrid level 395 407 396 408 REAL(wp), DIMENSION(1:3) :: volume_flow_area_l !< … … 647 659 648 660 ! 661 !-- Allocation of anelastic and Boussinesq approximation specific arrays 662 ALLOCATE( p_hydrostatic(nzb:nzt+1) ) 663 ALLOCATE( rho_air(nzb:nzt+1) ) 664 ALLOCATE( rho_air_zw(nzb:nzt+1) ) 665 ALLOCATE( drho_air(nzb:nzt+1) ) 666 ALLOCATE( drho_air_zw(nzb:nzt+1) ) 667 668 ! 669 !-- Density profile calculation for anelastic approximation 670 IF ( TRIM( approximation ) == 'anelastic' ) THEN 671 t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**( r_d / cp ) 672 DO k = nzb, nzt+1 673 p_hydrostatic(k) = surface_pressure * 100.0_wp * & 674 ( 1 - ( g * zu(k) ) / ( cp * t_surface ) & 675 )**( cp / r_d ) 676 rho_air(k) = ( p_hydrostatic(k) * & 677 ( 100000.0_wp / p_hydrostatic(k) & 678 )**( r_d / cp ) & 679 ) / ( r_d * pt_init(k) ) 680 ENDDO 681 DO k = nzb, nzt 682 rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) ) 683 ENDDO 684 rho_air_zw(nzt+1) = rho_air_zw(nzt) & 685 + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt) ) 686 ELSE 687 rho_air = 1.0_wp 688 rho_air_zw = 1.0_wp 689 ENDIF 690 691 !-- compute the inverse density array in order to avoid expencive divisions 692 drho_air = 1.0_wp / rho_air 693 drho_air_zw = 1.0_wp / rho_air_zw 694 695 ! 696 !-- Allocation of flux conversion arrays 697 ALLOCATE( heatflux_input_conversion(nzb:nzt+1) ) 698 ALLOCATE( waterflux_input_conversion(nzb:nzt+1) ) 699 ALLOCATE( momentumflux_input_conversion(nzb:nzt+1) ) 700 ALLOCATE( heatflux_output_conversion(nzb:nzt+1) ) 701 ALLOCATE( waterflux_output_conversion(nzb:nzt+1) ) 702 ALLOCATE( momentumflux_output_conversion(nzb:nzt+1) ) 703 704 ! 705 !-- calculate flux conversion factors according to approximation and in-/output mode 706 DO k = nzb, nzt+1 707 708 IF ( TRIM( flux_input_mode ) == 'kinematic' ) THEN 709 heatflux_input_conversion(k) = rho_air_zw(k) 710 waterflux_input_conversion(k) = rho_air_zw(k) 711 momentumflux_input_conversion(k) = rho_air_zw(k) 712 ELSEIF ( TRIM( flux_input_mode ) == 'dynamic' ) THEN 713 heatflux_input_conversion(k) = 1.0_wp / cp 714 waterflux_input_conversion(k) = 1.0_wp / l_v 715 momentumflux_input_conversion(k) = 1.0_wp 716 ENDIF 717 718 IF ( TRIM( flux_output_mode ) == 'kinematic' ) THEN 719 heatflux_output_conversion(k) = drho_air_zw(k) 720 waterflux_output_conversion(k) = drho_air_zw(k) 721 momentumflux_output_conversion(k) = drho_air_zw(k) 722 ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN 723 heatflux_output_conversion(k) = cp 724 waterflux_output_conversion(k) = l_v 725 momentumflux_output_conversion(k) = 1.0_wp 726 ENDIF 727 728 IF ( .NOT. humidity ) THEN 729 waterflux_input_conversion(k) = 1.0_wp 730 waterflux_output_conversion(k) = 1.0_wp 731 ENDIF 732 733 ENDDO 734 735 ! 736 !-- In case of multigrid method, compute grid lengths and grid factors for the 737 !-- grid levels with respective density on each grid 738 IF ( psolver(1:9) == 'multigrid' ) THEN 739 740 ALLOCATE( ddx2_mg(maximum_grid_level) ) 741 ALLOCATE( ddy2_mg(maximum_grid_level) ) 742 ALLOCATE( dzu_mg(nzb+1:nzt+1,maximum_grid_level) ) 743 ALLOCATE( dzw_mg(nzb+1:nzt+1,maximum_grid_level) ) 744 ALLOCATE( f1_mg(nzb+1:nzt,maximum_grid_level) ) 745 ALLOCATE( f2_mg(nzb+1:nzt,maximum_grid_level) ) 746 ALLOCATE( f3_mg(nzb+1:nzt,maximum_grid_level) ) 747 ALLOCATE( rho_air_mg(nzb:nzt+1,maximum_grid_level) ) 748 ALLOCATE( rho_air_zw_mg(nzb:nzt+1,maximum_grid_level) ) 749 750 dzu_mg(:,maximum_grid_level) = dzu 751 rho_air_mg(:,maximum_grid_level) = rho_air 752 ! 753 !-- Next line to ensure an equally spaced grid. 754 dzu_mg(1,maximum_grid_level) = dzu(2) 755 rho_air_mg(nzb,maximum_grid_level) = rho_air(nzb) + & 756 (rho_air(nzb) - rho_air(nzb+1)) 757 758 dzw_mg(:,maximum_grid_level) = dzw 759 rho_air_zw_mg(:,maximum_grid_level) = rho_air_zw 760 nzt_l = nzt 761 DO l = maximum_grid_level-1, 1, -1 762 dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1) 763 dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1) 764 rho_air_mg(nzb,l) = rho_air_mg(nzb,l+1) + (rho_air_mg(nzb,l+1) - rho_air_mg(nzb+1,l+1)) 765 rho_air_zw_mg(nzb,l) = rho_air_zw_mg(nzb,l+1) + (rho_air_zw_mg(nzb,l+1) - rho_air_zw_mg(nzb+1,l+1)) 766 rho_air_mg(nzb+1,l) = rho_air_mg(nzb+1,l+1) 767 rho_air_zw_mg(nzb+1,l) = rho_air_zw_mg(nzb+1,l+1) 768 nzt_l = nzt_l / 2 769 DO k = 2, nzt_l+1 770 dzu_mg(k,l) = dzu_mg(2*k-2,l+1) + dzu_mg(2*k-1,l+1) 771 dzw_mg(k,l) = dzw_mg(2*k-2,l+1) + dzw_mg(2*k-1,l+1) 772 rho_air_mg(k,l) = rho_air_mg(2*k-1,l+1) 773 rho_air_zw_mg(k,l) = rho_air_zw_mg(2*k-1,l+1) 774 ENDDO 775 ENDDO 776 777 nzt_l = nzt 778 dx_l = dx 779 dy_l = dy 780 DO l = maximum_grid_level, 1, -1 781 ddx2_mg(l) = 1.0_wp / dx_l**2 782 ddy2_mg(l) = 1.0_wp / dy_l**2 783 DO k = nzb+1, nzt_l 784 f2_mg(k,l) = rho_air_zw_mg(k,l) / ( dzu_mg(k+1,l) * dzw_mg(k,l) ) 785 f3_mg(k,l) = rho_air_zw_mg(k-1,l) / ( dzu_mg(k,l) * dzw_mg(k,l) ) 786 f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) & 787 * rho_air_mg(k,l) + f2_mg(k,l) + f3_mg(k,l) 788 ENDDO 789 nzt_l = nzt_l / 2 790 dx_l = dx_l * 2.0_wp 791 dy_l = dy_l * 2.0_wp 792 ENDDO 793 794 ENDIF 795 796 ! 649 797 !-- 3D-array for storing the dissipation, needed for calculating the sgs 650 798 !-- particle velocities … … 883 1031 vsws = 0.0_wp 884 1032 ENDIF 885 uswst = top_momentumflux_u 886 vswst = top_momentumflux_v 1033 uswst = top_momentumflux_u * momentumflux_input_conversion(nzt+1) 1034 vswst = top_momentumflux_v * momentumflux_input_conversion(nzt+1) 887 1035 888 1036 ! … … 1043 1191 us = 1E-30_wp 1044 1192 usws = 0.0_wp 1045 uswst = top_momentumflux_u 1193 uswst = top_momentumflux_u * momentumflux_input_conversion(nzt+1) 1046 1194 vsws = 0.0_wp 1047 vswst = top_momentumflux_v 1195 vswst = top_momentumflux_v * momentumflux_input_conversion(nzt+1) 1048 1196 IF ( humidity ) qs = 0.0_wp 1049 1197 IF ( passive_scalar ) ss = 0.0_wp … … 1165 1313 CALL disturb_heatflux 1166 1314 ELSE 1167 shf = surface_heatflux 1315 shf = surface_heatflux * heatflux_input_conversion(nzb) 1168 1316 ! 1169 1317 !-- Initialize shf with data from external file LSF_DATA … … 1178 1326 DO j = nysg, nyng 1179 1327 IF ( nzb_s_inner(j,i) /= 0 ) THEN 1180 shf(j,i) = wall_heatflux(0) 1328 shf(j,i) = wall_heatflux(0) & 1329 * heatflux_input_conversion(nzb_s_inner(j,i)) 1181 1330 ENDIF 1182 1331 ENDDO … … 1194 1343 ENDIF 1195 1344 IF ( constant_waterflux ) THEN 1196 qsws = surface_waterflux 1345 qsws = surface_waterflux * waterflux_input_conversion(nzb) 1197 1346 ! 1198 1347 !-- Over topography surface_waterflux is replaced by … … 1203 1352 DO j = nysg, nyng 1204 1353 IF ( nzb_s_inner(j,i) /= 0 ) THEN 1205 qsws(j,i) = wall_qflux(0) 1354 qsws(j,i) = wall_qflux(0) & 1355 * waterflux_input_conversion(nzb_s_inner(j,i)) 1206 1356 ENDIF 1207 1357 ENDDO … … 1241 1391 ! 1242 1392 !-- Prescribe to heat flux 1243 IF ( constant_top_heatflux ) tswst = top_heatflux 1393 IF ( constant_top_heatflux ) tswst = top_heatflux & 1394 * heatflux_input_conversion(nzt+1) 1244 1395 ! 1245 1396 !-- Prescribe zero latent flux at the top
Note: See TracChangeset
for help on using the changeset viewer.
