Changeset 4586 for palm/trunk/SOURCE
- Timestamp:
- Jul 1, 2020 4:16:43 PM (4 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 4 edited
-
header.f90 (modified) (25 diffs)
-
model_1d_mod.f90 (modified) (21 diffs)
-
surface_mod.f90 (modified) (3 diffs)
-
turbulence_closure_mod.f90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/header.f90
r4573 r4586 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 23 ! 22 ! 23 ! 24 24 ! Former revisions: 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Renamed rif to Ri (gradient Richardson number, 1D model) 28 ! and zeta (= z_mo / ol, stability parameter, 3D model) 29 ! 30 ! 4573 2020-06-24 13:08:47Z oliver.maas 27 31 ! added statement for pt_surface_heating_rate 28 ! 32 ! 29 33 ! 4536 2020-05-17 17:24:13Z raasch 30 34 ! output of restart data format added 31 ! 35 ! 32 36 ! 4473 2020-03-25 21:04:07Z gronemeier 33 37 ! revised message if wall_adjustment is used 34 ! 38 ! 35 39 ! 4444 2020-03-05 15:59:50Z raasch 36 40 ! bugfix: cpp-directives for serial mode added 37 ! 41 ! 38 42 ! 4360 2020-01-07 11:25:50Z suehring 39 43 ! Bugfix, character length too short, caused crash on NEC. 40 ! 44 ! 41 45 ! 4309 2019-11-26 18:49:59Z suehring 42 46 ! replaced recycling_yshift by y_shift 43 ! 47 ! 44 48 ! 4301 2019-11-22 12:09:09Z oliver.maas 45 49 ! 46 50 ! 4297 2019-11-21 10:37:50Z oliver.maas 47 51 ! Adjusted format for simulated time and related quantities 48 ! 52 ! 49 53 ! 4297 2019-11-21 10:37:50Z oliver.maas 50 54 ! adjusted message to the changed parameter recycling_yshift 51 ! 55 ! 52 56 ! 4227 2019-09-10 18:04:34Z gronemeier 53 57 ! implement new palm_date_time_mod 54 ! 58 ! 55 59 ! 4223 2019-09-10 09:20:47Z gronemeier 56 60 ! Write information about rotation angle 57 ! 61 ! 58 62 ! 4182 2019-08-22 15:20:23Z scharf 59 63 ! Corrected "Former revisions" section 60 ! 64 ! 61 65 ! 4168 2019-08-16 13:50:17Z suehring 62 66 ! Replace function get_topography_top_index by topo_top_ind 63 ! 67 ! 64 68 ! 4069 2019-07-01 14:05:51Z Giersch 65 ! Masked output running index mid has been introduced as a local variable to 66 ! avoid runtime error (Loop variable has been modified) in time_integration 67 ! 69 ! Masked output running index mid has been introduced as a local variable to 70 ! avoid runtime error (Loop variable has been modified) in time_integration 71 ! 68 72 ! 4023 2019-06-12 13:20:01Z maronga 69 73 ! Renamed "coupling start time" to "spinup time" 70 ! 74 ! 71 75 ! 4017 2019-06-06 12:16:46Z schwenkel 72 76 ! unused variable removed 73 ! 77 ! 74 78 ! 3655 2019-01-07 16:51:22Z knoop 75 79 ! Implementation of the PALM module interface … … 82 86 ! ------------ 83 87 !> Writing a header with all important information about the current run. 84 !> This subroutine is called three times, two times at the beginning 88 !> This subroutine is called three times, two times at the beginning 85 89 !> (writing information on files RUN_CONTROL and HEADER) and one time at the 86 90 !> end of the run, then writing additional information about CPU-usage on file … … 88 92 !-----------------------------------------------------------------------------! 89 93 SUBROUTINE header 90 94 91 95 92 96 USE arrays_3d, & … … 148 152 IMPLICIT NONE 149 153 150 154 151 155 CHARACTER (LEN=2) :: do2d_mode !< mode of 2D data output (xy, xz, yz) 152 156 153 157 CHARACTER (LEN=5) :: section_chr !< string indicating grid information where to output 2D slices 154 158 155 159 CHARACTER (LEN=10) :: host_chr !< string for hostname 156 160 157 161 CHARACTER (LEN=16) :: begin_chr !< string indication start time for the data output 158 162 CHARACTER (LEN=16) :: coor_chr !< dummy string 159 163 160 164 CHARACTER (LEN=26) :: ver_rev !< string for run identification 161 165 162 166 #if defined( __parallel ) 163 CHARACTER (LEN=32) :: cpl_name !< name of child domain (nesting mode only) 167 CHARACTER (LEN=32) :: cpl_name !< name of child domain (nesting mode only) 164 168 #endif 165 169 166 170 CHARACTER (LEN=40) :: output_format !< netcdf format 167 171 168 172 CHARACTER (LEN=70) :: char1 !< dummy varialbe used for various strings 169 173 CHARACTER (LEN=70) :: char2 !< string containing informating about the advected distance in case of Galilei transformation … … 175 179 CHARACTER (LEN=70) :: do3d_chr !< string indicating 3D output variables 176 180 CHARACTER (LEN=70) :: domask_chr !< string indicating masked output variables 177 CHARACTER (LEN=70) :: run_classification !< string classifying type of run, e.g. nested, coupled, etc. 178 181 CHARACTER (LEN=70) :: run_classification !< string classifying type of run, e.g. nested, coupled, etc. 182 179 183 CHARACTER (LEN=85) :: r_upper !< string indicating model top boundary condition for various quantities 180 184 CHARACTER (LEN=85) :: r_lower !< string indicating bottom boundary condition for various quantities 181 185 182 186 CHARACTER (LEN=86) :: coordinates !< string indicating height coordinates for profile-prescribed variables 183 187 CHARACTER (LEN=86) :: gradients !< string indicating gradients of profile-prescribed variables between the prescribed height coordinates … … 218 222 INTEGER(iwp) :: npe_total !< number of total PEs in a coupler (parent + child) 219 223 #endif 220 224 221 225 222 226 REAL(wp) :: cpuseconds_per_simulated_second !< CPU time (in s) per simulated second … … 506 510 IF ( .NOT. ocean_mode ) THEN 507 511 WRITE ( io, 250 ) dx, dy 508 512 509 513 DO i = 1, number_stretch_level_start+1 510 514 WRITE ( io, 253 ) i, dz(i) 511 515 ENDDO 512 516 513 517 WRITE( io, 251 ) (nx+1)*dx, (ny+1)*dy, zu(nzt+1) 514 518 515 519 IF ( ANY( dz_stretch_level_start_index < nzt+1 ) ) THEN 516 520 WRITE( io, '(A)', advance='no') ' Vertical stretching starts at height:' … … 535 539 ENDDO 536 540 ENDIF 537 541 538 542 ELSE 539 543 WRITE ( io, 250 ) dx, dy … … 541 545 WRITE ( io, 253 ) i, dz(i) 542 546 ENDDO 543 547 544 548 WRITE ( io, 251 ) (nx+1)*dx, (ny+1)*dy, zu(0) 545 549 546 550 IF ( ANY( dz_stretch_level_start_index > 0 ) ) THEN 547 551 WRITE( io, '(A)', advance='no') ' Vertical stretching starts at height:' … … 603 607 ENDDO 604 608 605 609 606 610 IF ( .NOT. large_scale_forcing ) THEN 607 611 WRITE ( io, 426 ) TRIM( coordinates ), TRIM( temperatures ), & … … 620 624 i = 1 621 625 DO WHILE ( ug_vertical_gradient_level_ind(i) /= -9999 ) 622 626 623 627 WRITE (coor_chr,'(F6.2,1X)') ug(ug_vertical_gradient_level_ind(i)) 624 628 ugcomponent = TRIM( ugcomponent ) // ' ' // TRIM( coor_chr ) … … 672 676 i = i + 1 673 677 ENDIF 674 678 675 679 ENDDO 676 680 … … 782 786 !-- Complex terrain 783 787 IF ( complex_terrain ) THEN 784 WRITE( io, 280 ) 788 WRITE( io, 280 ) 785 789 IF ( turbulent_inflow ) THEN 786 790 WRITE( io, 281 ) zu(topo_top_ind(0,0,0)) … … 846 850 r_upper = 'e(nzt+1) = e(nzt) = e(nzt-1)' 847 851 848 WRITE ( io, 301 ) 'e', r_lower, r_upper 852 WRITE ( io, 301 ) 'e', r_lower, r_upper 849 853 850 854 ENDIF … … 953 957 WRITE ( io, 317 ) bc_lr, bc_ns 954 958 IF ( .NOT. bc_lr_cyc .OR. .NOT. bc_ns_cyc ) THEN 955 WRITE ( io, 318 ) use_cmax, pt_damping_width, pt_damping_factor 959 WRITE ( io, 318 ) use_cmax, pt_damping_width, pt_damping_factor 956 960 IF ( turbulent_inflow ) THEN 957 961 IF ( y_shift == 0 ) THEN … … 1009 1013 TRIM( gradients ), TRIM( slices ) 1010 1014 ELSE 1011 WRITE ( io, 428 ) 1015 WRITE ( io, 428 ) 1012 1016 ENDIF 1013 1017 … … 1022 1026 i = 1 1023 1027 DO WHILE ( q_vertical_gradient_level_ind(i) /= -9999 ) 1024 1028 1025 1029 WRITE (coor_chr,'(E8.1,4X)') q_init(q_vertical_gradient_level_ind(i)) 1026 1030 temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr ) … … 1028 1032 WRITE (coor_chr,'(E8.1,4X)') q_vertical_gradient(i) 1029 1033 gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr ) 1030 1034 1031 1035 WRITE (coor_chr,'(I8,4X)') q_vertical_gradient_level_ind(i) 1032 1036 slices = TRIM( slices ) // ' ' // TRIM( coor_chr ) 1033 1037 1034 1038 WRITE (coor_chr,'(F8.1,4X)') q_vertical_gradient_level(i) 1035 1039 coordinates = TRIM( coordinates ) // ' ' // TRIM( coor_chr ) … … 1058 1062 i = 1 1059 1063 DO WHILE ( s_vertical_gradient_level_ind(i) /= -9999 ) 1060 1064 1061 1065 WRITE (coor_chr,'(E8.1,4X)') s_init(s_vertical_gradient_level_ind(i)) 1062 1066 temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr ) … … 1064 1068 WRITE (coor_chr,'(E8.1,4X)') s_vertical_gradient(i) 1065 1069 gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr ) 1066 1070 1067 1071 WRITE (coor_chr,'(I8,4X)') s_vertical_gradient_level_ind(i) 1068 1072 slices = TRIM( slices ) // ' ' // TRIM( coor_chr ) 1069 1073 1070 1074 WRITE (coor_chr,'(F8.1,4X)') s_vertical_gradient_level(i) 1071 1075 coordinates = TRIM( coordinates ) // ' ' // TRIM( coor_chr ) … … 1081 1085 WRITE ( io, 422 ) TRIM( coordinates ), TRIM( temperatures ), & 1082 1086 TRIM( gradients ), TRIM( slices ) 1083 ENDIF 1087 ENDIF 1084 1088 1085 1089 ! … … 1586 1590 ENDIF 1587 1591 IF ( humidity .AND. q_surface_initial_change /= 0.0_wp ) THEN 1588 WRITE ( io, 477 ) q_surface_initial_change 1592 WRITE ( io, 477 ) q_surface_initial_change 1589 1593 ENDIF 1590 1594 IF ( passive_scalar .AND. q_surface_initial_change /= 0.0_wp ) THEN 1591 WRITE ( io, 478 ) q_surface_initial_change 1595 WRITE ( io, 478 ) q_surface_initial_change 1592 1596 ENDIF 1593 1597 … … 1727 1731 ' Tunnel width: ', F6.2 ) 1728 1732 274 FORMAT ( ' Tunnel in ', A, ' direction.' / & 1729 ' Tunnel height: ', F6.2, / & 1733 ' Tunnel height: ', F6.2, / & 1730 1734 ' Tunnel-wall depth: ', F6.2 / & 1731 1735 ' Tunnel width: ', F6.2, / & … … 1756 1760 ' z_mo = ',F6.2,' m z0 =',F7.4,' m z0h =',F8.5,& 1757 1761 ' m kappa =',F5.2/ & 1758 ' Rif value range: ',F8.2,' <= rif<=',F6.2)1762 ' zeta value range: ',F8.2,' <= zeta <=',F6.2) 1759 1763 306 FORMAT (' Predefined constant heatflux: ',F9.6,' K m/s') 1760 1764 307 FORMAT (' Heatflux has a random normal distribution') … … 1762 1766 309 FORMAT (' Predefined constant salinityflux: ',F9.6,' psu m/s') 1763 1767 310 FORMAT (//' 1D-Model:'// & 1764 ' Ri f value range: ',F6.2,' <= rif<=',F6.2)1768 ' Ri value range: ',F6.2,' <= Ri <=',F6.2) 1765 1769 311 FORMAT (' Predefined constant humidity flux: ',E10.3,' kg/kg m/s') 1766 1770 312 FORMAT (' Predefined surface humidity') -
palm/trunk/SOURCE/model_1d_mod.f90
r4449 r4586 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 23 ! 22 ! 23 ! 24 24 ! Former revisions: 25 25 ! ----------------- 26 26 ! $Id$ 27 ! renamed Richardson flux number into gradient Richardson number 28 ! 29 ! 4449 2020-03-09 14:43:16Z suehring 27 30 ! Set intermediate_timestep_count back to zero after 1D-model integration. 28 31 ! This is required e.g. for initial calls of calc_mean_profile. … … 120 123 REAL(wp), DIMENSION(:), ALLOCATABLE :: l1d_init !< initial mixing length (1d-model) 121 124 REAL(wp), DIMENSION(:), ALLOCATABLE :: l1d_diss !< mixing length for dissipation (1d-model) 122 REAL(wp), DIMENSION(:), ALLOCATABLE :: ri f1d !< Richardson fluxnumber (1d-model)125 REAL(wp), DIMENSION(:), ALLOCATABLE :: ri1d !< gradient Richardson number (1d-model) 123 126 REAL(wp), DIMENSION(:), ALLOCATABLE :: te_diss !< tendency of diss (1d-model) 124 127 REAL(wp), DIMENSION(:), ALLOCATABLE :: te_dissm !< weighted tendency of diss for previous sub-timestep (1d-model) … … 174 177 !-- Public variables 175 178 PUBLIC damp_level_1d, damp_level_ind_1d, diss1d, dt_pr_1d, & 176 dt_run_control_1d, e1d, end_time_1d, kh1d, km1d, l1d, ri f1d, u1d,&179 dt_run_control_1d, e1d, end_time_1d, kh1d, km1d, l1d, ri1d, u1d, & 177 180 us1d, usws1d, v1d, vsws1d 178 181 … … 196 199 e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1), & 197 200 km1d(nzb:nzt+1), l1d(nzb:nzt+1), l1d_init(nzb:nzt+1), & 198 l1d_diss(nzb:nzt+1), ri f1d(nzb:nzt+1), te_diss(nzb:nzt+1),&201 l1d_diss(nzb:nzt+1), ri1d(nzb:nzt+1), te_diss(nzb:nzt+1), & 199 202 te_dissm(nzb:nzt+1), te_e(nzb:nzt+1), & 200 203 te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1), & … … 211 214 e1d = 0.0_wp; e1d_p = 0.0_wp 212 215 kh1d = 0.0_wp; km1d = 0.0_wp 213 ri f1d = 0.0_wp216 ri1d = 0.0_wp 214 217 ! 215 218 !-- Compute the mixing length … … 414 417 ! 415 418 !-- dissipation rate 416 IF ( ri f1d(k) >= 0.0_wp ) THEN419 IF ( ri1d(k) >= 0.0_wp ) THEN 417 420 alpha_buoyancy = 1.0_wp - l1d(k) / lambda 418 421 ELSE … … 604 607 IF ( constant_flux_layer ) THEN 605 608 ! 606 !-- Compute theta* using Ri fnumbers of the previous time step607 IF ( ri f1d(nzb+1) >= 0.0_wp ) THEN609 !-- Compute theta* using Ri numbers of the previous time step 610 IF ( ri1d(nzb+1) >= 0.0_wp ) THEN 608 611 ! 609 612 !-- Stable stratification 610 613 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & 611 ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri f1d(nzb+1) *&614 ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) * & 612 615 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 613 616 ) … … 615 618 ! 616 619 !-- Unstable stratification 617 a = SQRT( 1.0_wp - 16.0_wp * ri f1d(nzb+1) )618 b = SQRT( 1.0_wp - 16.0_wp * ri f1d(nzb+1) /&620 a = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) ) 621 b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / & 619 622 zu(nzb+1) * z0h1d ) 620 623 … … 630 633 !> 2018-04-23, gronemeier 631 634 ! 632 !-- Compute the Richardson-fluxnumbers,635 !-- Compute the gradient Richardson numbers, 633 636 !-- first at the top of the constant-flux layer using u* of the 634 637 !-- previous time step (+1E-30, if u* = 0), then in the remaining area. 635 !-- There the rif-numbers of the previous time step are used.638 !-- There, the Ri numbers of the previous time step are used. 636 639 637 640 IF ( constant_flux_layer ) THEN … … 643 646 flux = ts1d + 0.61_wp * pt_init(k) * qs1d 644 647 ENDIF 645 ri f1d(nzb+1) = zu(nzb+1) * kappa * g * flux /&648 ri1d(nzb+1) = zu(nzb+1) * kappa * g * flux / & 646 649 ( pt_0 * ( us1d**2 + 1E-30_wp ) ) 647 650 ENDIF … … 659 662 ) * dd2zu(k) 660 663 ENDIF 661 IF ( ri f1d(k) >= 0.0_wp ) THEN662 ri f1d(k) = g / pt_0 * flux /&664 IF ( ri1d(k) >= 0.0_wp ) THEN 665 ri1d(k) = g / pt_0 * flux / & 663 666 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & 664 667 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & … … 666 669 ) 667 670 ELSE 668 ri f1d(k) = g / pt_0 * flux /&671 ri1d(k) = g / pt_0 * flux / & 669 672 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & 670 673 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & 671 674 + 1E-30_wp & 672 ) * ( 1.0_wp - 16.0_wp * ri f1d(k) )**0.25_wp675 ) * ( 1.0_wp - 16.0_wp * ri1d(k) )**0.25_wp 673 676 ENDIF 674 677 ENDDO 675 678 ! 676 !-- Richardson -numbers must remain restricted to a realistic value679 !-- Richardson numbers must remain restricted to a realistic value 677 680 !-- range. It is exceeded excessively for very small velocities 678 681 !-- (u,v --> 0). 679 WHERE ( ri f1d < -5.0_wp ) rif1d = -5.0_wp680 WHERE ( ri f1d > 1.0_wp ) rif1d = 1.0_wp682 WHERE ( ri1d < -5.0_wp ) ri1d = -5.0_wp 683 WHERE ( ri1d > 1.0_wp ) ri1d = 1.0_wp 681 684 682 685 ! … … 685 688 uv_total = SQRT( u1d(nzb+1)**2 + v1d(nzb+1)**2 ) 686 689 687 IF ( ri f1d(nzb+1) >= 0.0_wp ) THEN690 IF ( ri1d(nzb+1) >= 0.0_wp ) THEN 688 691 ! 689 692 !-- Stable stratification 690 693 us1d = kappa * uv_total / ( & 691 LOG( zu(nzb+1) / z01d ) + 5.0_wp * ri f1d(nzb+1) *&694 LOG( zu(nzb+1) / z01d ) + 5.0_wp * ri1d(nzb+1) * & 692 695 ( zu(nzb+1) - z01d ) / zu(nzb+1) & 693 696 ) … … 695 698 ! 696 699 !-- Unstable stratification 697 a = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri f1d(nzb+1) ) )698 b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri f1d(nzb+1) /&700 a = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) ) ) 701 b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / & 699 702 zu(nzb+1) * z01d ) ) 700 703 us1d = kappa * uv_total / ( & … … 728 731 ! 729 732 !-- Compute q* 730 IF ( ri f1d(nzb+1) >= 0.0_wp ) THEN733 IF ( ri1d(nzb+1) >= 0.0_wp ) THEN 731 734 ! 732 735 !-- Stable stratification 733 736 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & 734 ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri f1d(nzb+1) *&737 ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) * & 735 738 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 736 739 ) … … 738 741 ! 739 742 !-- Unstable stratification 740 a = SQRT( 1.0_wp - 16.0_wp * ri f1d(nzb+1) )741 b = SQRT( 1.0_wp - 16.0_wp * ri f1d(nzb+1) /&743 a = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) ) 744 b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / & 742 745 zu(nzb+1) * z0h1d ) 743 746 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & … … 759 762 IF ( mixing_length_1d == 'blackadar' ) THEN 760 763 DO k = nzb+1, nzt 761 IF ( ri f1d(k) >= 0.0_wp ) THEN762 l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * ri f1d(k) )764 IF ( ri1d(k) >= 0.0_wp ) THEN 765 l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * ri1d(k) ) 763 766 l1d_diss(k) = l1d(k) 764 767 ELSE 765 768 l1d(k) = l1d_init(k) 766 769 l1d_diss(k) = l1d_init(k) * & 767 SQRT( 1.0_wp - 16.0_wp * ri f1d(k) )770 SQRT( 1.0_wp - 16.0_wp * ri1d(k) ) 768 771 ENDIF 769 772 ENDDO … … 793 796 !-- Prandtl-Kolmogorov, respectively 794 797 IF ( constant_flux_layer ) THEN 795 IF ( ri f1d(nzb+1) >= 0.0_wp ) THEN798 IF ( ri1d(nzb+1) >= 0.0_wp ) THEN 796 799 km1d(nzb+1) = us1d * kappa * zu(nzb+1) / & 797 ( 1.0_wp + 5.0_wp * ri f1d(nzb+1) )800 ( 1.0_wp + 5.0_wp * ri1d(nzb+1) ) 798 801 ELSE 799 802 km1d(nzb+1) = us1d * kappa * zu(nzb+1) * & 800 ( 1.0_wp - 16.0_wp * ri f1d(nzb+1) )**0.25_wp803 ( 1.0_wp - 16.0_wp * ri1d(nzb+1) )**0.25_wp 801 804 ENDIF 802 805 ENDIF … … 823 826 !-- kh = phim / phih * km 824 827 DO k = nzb+1, nzt 825 IF ( ri f1d(k) >= 0.0_wp ) THEN828 IF ( ri1d(k) >= 0.0_wp ) THEN 826 829 kh1d(k) = km1d(k) 827 830 ELSE 828 kh1d(k) = km1d(k) * ( 1.0_wp - 16.0_wp * ri f1d(k) )**0.25_wp831 kh1d(k) = km1d(k) * ( 1.0_wp - 16.0_wp * ri1d(k) )**0.25_wp 829 832 ENDIF 830 833 ENDDO … … 1026 1029 WRITE ( 17, 101 ) 1027 1030 DO k = nzt+1, nzb, -1 1028 WRITE ( 17, 103) k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), &1029 ri f1d(k), km1d(k), kh1d(k), l1d(k), diss1d(k)1031 WRITE ( 17, 103) k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), & 1032 ri1d(k), km1d(k), kh1d(k), l1d(k), diss1d(k) 1030 1033 ENDDO 1031 1034 WRITE ( 17, 101 ) … … 1045 1048 101 FORMAT ('#',111('-')) 1046 1049 102 FORMAT ('# k zu u v pt e ', & 1047 ' rifKm Kh l diss ')1050 'Ri Km Kh l diss ') 1048 1051 103 FORMAT (1X,I4,1X,F7.1,9(1X,E10.3)) 1049 1052 -
palm/trunk/SOURCE/surface_mod.f90
r4559 r4586 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 23 ! 22 ! 23 ! 24 24 ! Former revisions: 25 25 ! ----------------- 26 26 ! $Id$ 27 ! renamed Richardson flux number into gradient Richardson number (1D model) 28 ! 29 ! 4559 2020-06-11 08:51:48Z raasch 27 30 ! File re-formatted to follow the PALM coding standard 28 31 ! … … 166 169 167 170 USE model_1d_mod, & 168 ONLY: ri f1d,&171 ONLY: ri1d, & 169 172 us1d, & 170 173 usws1d, & … … 2535 2538 IF ( .NOT. constant_diffusion ) THEN 2536 2539 IF ( constant_flux_layer ) THEN 2537 surf%ol(num_h) = surf%z_mo(num_h) / ( ri f1d(nzb+1) + 1.0E-20_wp )2540 surf%ol(num_h) = surf%z_mo(num_h) / ( ri1d(nzb+1) + 1.0E-20_wp ) 2538 2541 surf%us(num_h) = us1d 2539 2542 surf%usws(num_h) = usws1d -
palm/trunk/SOURCE/turbulence_closure_mod.f90
r4510 r4586 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 23 ! 22 ! 23 ! 24 24 ! Former revisions: 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Renamed rif to richardson_number (Richardson Flux Number to Gradient Richardson Number) 28 ! 29 ! 4510 2020-04-29 14:19:18Z raasch 27 30 ! file re-formatted to follow the PALM coding standard 28 31 ! … … 4010 4013 4011 4014 REAL(wp), DIMENSION(nzb+1:nzt) :: ml_stratification !< mixing length according to stratification 4012 REAL(wp), DIMENSION(nzb+1:nzt) :: ri f !< Richardson fluxnumber4015 REAL(wp), DIMENSION(nzb+1:nzt) :: richardson_number !< gradient Richardson number 4013 4016 4014 4017 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) :: var !< temperature … … 4111 4114 4112 4115 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) & 4113 !$ACC PRIVATE(ml_stratification, duv2_dz2, ri f, dvar_dz) &4116 !$ACC PRIVATE(ml_stratification, duv2_dz2, richardson_number, dvar_dz) & 4114 4117 !$ACC PRESENT(diss, e, u, v, var, wall_flags_total_0) & 4115 4118 !$ACC PRESENT(dd2zu, ml_blackadar, delta) … … 4117 4120 DO j = nys, nyn 4118 4121 ! 4119 !-- Calculate Richardson -fluxnumber4122 !-- Calculate Richardson number 4120 4123 IF ( use_single_reference_value ) THEN 4121 4124 … … 4128 4131 + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 + 1E-30_wp 4129 4132 4130 rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ), 1.0_wp ) 4133 richardson_number(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, & 4134 -5.0_wp ), & 4135 1.0_wp ) 4131 4136 ENDDO 4132 4137 ELSE … … 4143 4148 + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 + 1E-30_wp 4144 4149 4145 rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ), 1.0_wp ) 4150 richardson_number(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, & 4151 -5.0_wp ), & 4152 1.0_wp ) 4146 4153 ENDDO 4147 4154 ENDIF … … 4150 4157 !$ACC LOOP PRIVATE(k) 4151 4158 DO k = nzb+1, nzt 4152 IF ( rif(k) >= 0.0_wp ) THEN 4153 ml_stratification(k) = ml_blackadar(k) / ( 1.0_wp + 5.0_wp * rif(k) ) 4159 IF ( richardson_number(k) >= 0.0_wp ) THEN 4160 ml_stratification(k) = ml_blackadar(k) & 4161 / ( 1.0_wp + 5.0_wp * richardson_number(k) ) 4154 4162 ELSE 4155 ml_stratification(k) = ml_blackadar(k) * SQRT( 1.0_wp - 16.0_wp * rif(k) ) 4163 ml_stratification(k) = ml_blackadar(k) & 4164 * SQRT( 1.0_wp - 16.0_wp * richardson_number(k) ) 4156 4165 ENDIF 4157 4166 ENDDO … … 4274 4283 4275 4284 REAL(wp), DIMENSION(nzb+1:nzt) :: ml_stratification !< mixing length according to stratification 4276 REAL(wp), DIMENSION(nzb+1:nzt) :: ri f !< Richardson fluxnumber4285 REAL(wp), DIMENSION(nzb+1:nzt) :: richardson_number !< gradient Richardson number 4277 4286 4278 4287 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) :: var !< temperature … … 4348 4357 4349 4358 ! 4350 !-- Calculate Richardson -fluxnumber4359 !-- Calculate Richardson number 4351 4360 IF ( use_single_reference_value ) THEN 4352 4361 DO k = nzb+1, nzt … … 4356 4365 + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 + 1E-30_wp 4357 4366 4358 rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ), 1.0_wp ) 4367 richardson_number(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ), & 4368 1.0_wp ) 4359 4369 ENDDO 4360 4370 ELSE … … 4365 4375 + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 + 1E-30_wp 4366 4376 4367 rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ), 1.0_wp ) 4377 richardson_number(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ), & 4378 1.0_wp ) 4368 4379 ENDDO 4369 4380 ENDIF … … 4371 4382 !-- Calculate diabatic mixing length using Dyer-profile functions 4372 4383 DO k = nzb+1, nzt 4373 IF ( ri f(k) >= 0.0_wp ) THEN4374 ml_stratification(k) = ml_blackadar(k) / ( 1.0_wp + 5.0_wp * ri f(k) )4384 IF ( richardson_number(k) >= 0.0_wp ) THEN 4385 ml_stratification(k) = ml_blackadar(k) / ( 1.0_wp + 5.0_wp * richardson_number(k) ) 4375 4386 ELSE 4376 ml_stratification(k) = ml_blackadar(k) * SQRT( 1.0_wp - 16.0_wp * rif(k) ) 4387 ml_stratification(k) = ml_blackadar(k) & 4388 * SQRT( 1.0_wp - 16.0_wp * richardson_number(k) ) 4377 4389 ENDIF 4378 4390 … … 4772 4784 REAL(wp) :: var_reference !< reference temperature 4773 4785 4774 !DIR$ ATTRIBUTES ALIGN:64:: ml_local_profile, ml_stratification, ri f4786 !DIR$ ATTRIBUTES ALIGN:64:: ml_local_profile, ml_stratification, richardson_number 4775 4787 REAL(wp), DIMENSION(nzb+1:nzt) :: ml_local_profile !< mixing length (all heights) 4776 4788 REAL(wp), DIMENSION(nzb+1:nzt) :: ml_stratification !< mixing length according to stratification 4777 REAL(wp), DIMENSION(nzb+1:nzt) :: ri f !< Richardson fluxnumber4789 REAL(wp), DIMENSION(nzb+1:nzt) :: richardson_number !< gradient Richardson number 4778 4790 4779 4791 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) :: var !< temperature … … 4937 4949 !$OMP DO 4938 4950 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) & 4939 !$ACC PRIVATE(dvar_dz, duv2_dz2, ml_stratification, ml_local_profile, ri f) &4951 !$ACC PRIVATE(dvar_dz, duv2_dz2, ml_stratification, ml_local_profile, richardson_number) & 4940 4952 !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, u, v, delta, ml_blackadar) & 4941 4953 !$ACC PRESENT(kh, km, sums_l_l, rmask) … … 4943 4955 DO j = nysg, nyng 4944 4956 ! 4945 !-- Calculate Richardson -fluxnumber4957 !-- Calculate Richardson number 4946 4958 IF ( use_single_reference_value ) THEN 4947 4959 !$ACC LOOP PRIVATE(k) … … 4953 4965 + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 + 1E-30_wp 4954 4966 4955 rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ), 1.0_wp ) 4967 richardson_number(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, & 4968 -5.0_wp ), & 4969 1.0_wp ) 4956 4970 ENDDO 4957 4971 ELSE … … 4964 4978 + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 + 1E-30_wp 4965 4979 4966 rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ), 1.0_wp ) 4980 richardson_number(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, & 4981 -5.0_wp ), & 4982 1.0_wp ) 4967 4983 ENDDO 4968 4984 ENDIF … … 4972 4988 !$ACC LOOP PRIVATE(k) 4973 4989 DO k = nzb+1, nzt 4974 IF ( rif(k) >= 0.0_wp ) THEN 4975 ml_stratification(k) = ml_blackadar(k) / ( 1.0_wp + 5.0_wp * rif(k) ) 4990 IF ( richardson_number(k) >= 0.0_wp ) THEN 4991 ml_stratification(k) = ml_blackadar(k) & 4992 / ( 1.0_wp + 5.0_wp * richardson_number(k) ) 4976 4993 ELSE 4977 4994 ml_stratification(k) = ml_blackadar(k)
Note: See TracChangeset
for help on using the changeset viewer.
