Changeset 1001 for palm/trunk/SOURCE/init_1d_model.f90
- Timestamp:
- Sep 13, 2012 2:08:46 PM (12 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/init_1d_model.f90 (modified) (19 diffs)
Legend:
- Unmodified
- Added
- Removed
-
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
Note: See TracChangeset
for help on using the changeset viewer.
