Changeset 3668 for palm/trunk/SOURCE
- Timestamp:
- Jan 14, 2019 12:49:24 PM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 7 edited
-
check_parameters.f90 (modified) (2 diffs)
-
land_surface_model_mod.f90 (modified) (4 diffs)
-
modules.f90 (modified) (2 diffs)
-
parin.f90 (modified) (3 diffs)
-
read_restart_data_mod.f90 (modified) (3 diffs)
-
surface_layer_fluxes_mod.f90 (modified) (9 diffs)
-
write_restart_data_mod.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/check_parameters.f90
r3655 r3668 25 25 ! ----------------- 26 26 ! $Id$ 27 ! most_method removed 28 ! 29 ! 3655 2019-01-07 16:51:22Z knoop 27 30 ! Formatting 28 31 ! … … 3790 3793 ENDIF 3791 3794 3792 !3793 !-- Check for valid setting of most_method3794 IF ( TRIM( most_method ) /= 'circular' .AND. &3795 TRIM( most_method ) /= 'newton' .AND. &3796 TRIM( most_method ) /= 'lookup' ) THEN3797 message_string = 'most_method = "' // TRIM( most_method ) // &3798 '" is unknown'3799 CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )3800 ENDIF3801 3795 3802 3796 ! -
palm/trunk/SOURCE/land_surface_model_mod.f90
r3655 r3668 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Removed most_method 28 ! 29 ! 3655 2019-01-07 16:51:22Z knoop 27 30 ! nopointer option removed 28 31 ! … … 1324 1327 1325 1328 USE control_parameters, & 1326 ONLY: bc_pt_b, bc_q_b, constant_flux_layer, message_string, & 1327 most_method 1329 ONLY: bc_pt_b, bc_q_b, constant_flux_layer, message_string 1328 1330 1329 1331 … … 1446 1448 1447 1449 IF ( TRIM( surface_type ) == 'water' ) THEN 1448 1449 IF ( TRIM( most_method ) == 'lookup' ) THEN1450 WRITE( message_string, * ) 'surface_type = ', surface_type, &1451 ' is not allowed in combination with ', &1452 'most_method = ', most_method1453 CALL message( 'lsm_check_parameters', 'PA0414', 1, 2, 0, 6, 0 )1454 ENDIF1455 1450 1456 1451 IF ( water_type == 0 ) THEN … … 1531 1526 1532 1527 IF ( TRIM( surface_type ) == 'netcdf' ) THEN 1533 IF ( ANY( water_type_f%var /= water_type_f%fill ) .AND. &1534 TRIM( most_method ) == 'lookup' ) THEN1535 WRITE( message_string, * ) 'water-surfaces are not allowed in ' // &1536 'combination with most_method = ', &1537 TRIM( most_method )1538 CALL message( 'lsm_check_parameters', 'PA0999', 2, 2, 0, 6, 0 )1539 ENDIF1540 1528 ! 1541 1529 !-- MS: Some problme here, after calling message everythings stucks at -
palm/trunk/SOURCE/modules.f90
r3648 r3668 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Removed most_method 28 ! 29 ! 3648 2019-01-02 16:35:46Z suehring 27 30 ! -surface_data_output +surface_output 28 31 ! … … 1090 1093 CHARACTER (LEN=8) :: coupling_char = '' !< appended to filenames in coupled or nested runs ('_O': ocean PE, 1091 1094 !< '_NV': vertically nested atmosphere PE, '_N##': PE of nested domain ## 1092 CHARACTER (LEN=8) :: most_method = 'newton' !< namelist parameter1093 1095 CHARACTER (LEN=10) :: run_date = ' ' !< date of simulation run 1094 1096 CHARACTER (LEN=8) :: run_time = ' ' !< time of simulation run -
palm/trunk/SOURCE/parin.f90
r3649 r3668 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Removed most_method 28 ! 29 ! 3649 2019-01-02 16:52:21Z suehring 27 30 ! Delete debug-print statements 28 31 ! … … 573 576 loop_optimization, lsf_exception, masking_method, mg_cycles, & 574 577 mg_switch_to_pe0_level, mixing_length_1d, momentum_advec, & 575 most_method, &576 578 netcdf_precision, neutral, ngsrb, & 577 579 nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega, & … … 644 646 loop_optimization, lsf_exception, masking_method, mg_cycles, & 645 647 mg_switch_to_pe0_level, mixing_length_1d, momentum_advec, & 646 most_method, &647 648 netcdf_precision, neutral, ngsrb, & 648 649 nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega, & -
palm/trunk/SOURCE/read_restart_data_mod.f90
r3655 r3668 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Removed most_method and increased binary version 28 ! 29 ! 3655 2019-01-07 16:51:22Z knoop 27 30 ! Implementation of the PALM module interface 28 31 ! … … 203 206 READ ( 13 ) version_on_file 204 207 205 binary_version_global = '4. 7'208 binary_version_global = '4.8' 206 209 IF ( TRIM( version_on_file ) /= TRIM( binary_version_global ) ) THEN 207 210 WRITE( message_string, * ) 'version mismatch concerning ', & … … 488 491 CASE ( 'momentum_advec' ) 489 492 READ ( 13 ) momentum_advec 490 CASE ( 'most_method' )491 READ ( 13 ) most_method492 493 CASE ( 'netcdf_precision' ) 493 494 READ ( 13 ) netcdf_precision -
palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
r3655 r3668 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Removed methods "circular" and "lookup" 29 ! 30 ! 3655 2019-01-07 16:51:22Z knoop 28 31 ! OpenACC port for SPEC 29 32 ! … … 228 231 ! ------------ 229 232 !> Diagnostic computation of vertical fluxes in the constant flux layer from the 230 !> values of the variables at grid point k=1. Three different methods are 231 !> available: 232 !> 1) the "old" version (most_method = 'circular') which is fast, but inaccurate 233 !> 2) a Newton iteration method (most_method = 'newton'), which is accurate, but 234 !> slower 235 !> 3) a method using a lookup table which is fast and accurate. Note, however, 236 !> that this method cannot be used in case of roughness heterogeneity 233 !> values of the variables at grid point k=1 based on Newton iteration 237 234 !> 238 235 !> @todo (re)move large_scale_forcing actions … … 261 258 intermediate_timestep_count, intermediate_timestep_count_max, & 262 259 land_surface, large_scale_forcing, lsf_surf, message_string, & 263 most_method, neutral, passive_scalar, pt_surface, q_surface,&260 neutral, passive_scalar, pt_surface, q_surface, & 264 261 run_coupled, surface_pressure, simulated_time, terminate_run, & 265 262 time_since_reference_point, urban_surface, & … … 400 397 !-- First, calculate the new Obukhov length, then new friction velocity, 401 398 !-- followed by the new scaling parameters (th*, q*, etc.), and the new 402 !-- surface fluxes if required. The old routine ("circular") requires a 403 !-- different order of calls as the scaling parameters from the previous time 404 !-- steps are used to calculate the Obukhov length 405 406 ! 407 !-- Depending on setting of most_method use the "old" routine 408 !-- Note, each routine is called for different surface types. 399 !-- surface fluxes if required. Note, each routine is called for different surface types. 409 400 !-- First call for default-type horizontal surfaces, for natural- and 410 401 !-- urban-type surfaces. Note, at this place only upward-facing horizontal 411 !-- surfaces are treted. 412 IF ( most_method == 'circular' ) THEN 413 ! 414 !-- Default-type upward-facing horizontal surfaces 415 IF ( surf_def_h(0)%ns >= 1 ) THEN 416 surf => surf_def_h(0) 417 CALL calc_scaling_parameters 418 CALL calc_uvw_abs 419 IF ( .NOT. neutral ) CALL calc_ol 420 CALL calc_us 421 CALL calc_surface_fluxes 422 423 IF ( do_output_at_2m ) THEN 424 CALL calc_pt_near_surface ( '2m' ) 425 ENDIF 426 ENDIF 427 ! 428 !-- Natural-type horizontal surfaces 429 IF ( surf_lsm_h%ns >= 1 ) THEN 430 surf => surf_lsm_h 431 CALL calc_scaling_parameters 432 CALL calc_uvw_abs 433 IF ( .NOT. neutral ) CALL calc_ol 434 CALL calc_us 435 CALL calc_surface_fluxes 436 IF ( do_output_at_2m ) THEN 437 CALL calc_pt_near_surface ( '2m' ) 438 ENDIF 439 ENDIF 440 ! 441 !-- Urban-type horizontal surfaces 442 IF ( surf_usm_h%ns >= 1 ) THEN 443 surf => surf_usm_h 444 CALL calc_scaling_parameters 445 CALL calc_uvw_abs 446 IF ( .NOT. neutral ) CALL calc_ol 447 CALL calc_us 448 CALL calc_surface_fluxes 449 IF ( do_output_at_2m ) THEN 450 CALL calc_pt_near_surface ( '2m' ) 451 ENDIF 452 IF ( indoor_model ) THEN 453 CALL calc_pt_near_surface ( '10cm' ) 454 ENDIF 455 ENDIF 456 ! 457 !-- Use either Newton iteration or a lookup table for the bulk Richardson 458 !-- number to calculate the Obukhov length 459 ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' ) THEN 460 ! 461 !-- Default-type upward-facing horizontal surfaces 462 IF ( surf_def_h(0)%ns >= 1 ) THEN 463 surf => surf_def_h(0) 464 CALL calc_uvw_abs 465 IF ( .NOT. neutral ) CALL calc_ol 466 CALL calc_us 467 CALL calc_scaling_parameters 468 CALL calc_surface_fluxes 469 IF ( do_output_at_2m ) THEN 470 CALL calc_pt_near_surface ( '2m' ) 471 ENDIF 472 ENDIF 473 ! 474 !-- Natural-type horizontal surfaces 475 IF ( surf_lsm_h%ns >= 1 ) THEN 476 surf => surf_lsm_h 477 CALL calc_uvw_abs 478 IF ( .NOT. neutral ) CALL calc_ol 479 CALL calc_us 480 CALL calc_scaling_parameters 481 CALL calc_surface_fluxes 482 IF ( do_output_at_2m ) THEN 483 CALL calc_pt_near_surface ( '2m' ) 484 ENDIF 485 ENDIF 486 ! 487 !-- Urban-type horizontal surfaces 488 IF ( surf_usm_h%ns >= 1 ) THEN 489 surf => surf_usm_h 490 CALL calc_uvw_abs 491 IF ( .NOT. neutral ) CALL calc_ol 492 CALL calc_us 493 CALL calc_scaling_parameters 494 CALL calc_surface_fluxes 495 IF ( do_output_at_2m ) THEN 496 CALL calc_pt_near_surface ( '2m' ) 497 ENDIF 498 IF ( indoor_model ) THEN 499 CALL calc_pt_near_surface ( '10cm' ) 500 ENDIF 501 ENDIF 502 402 !-- surfaces are treated. 403 404 ! 405 !-- Default-type upward-facing horizontal surfaces 406 IF ( surf_def_h(0)%ns >= 1 ) THEN 407 surf => surf_def_h(0) 408 CALL calc_uvw_abs 409 IF ( .NOT. neutral ) CALL calc_ol 410 CALL calc_us 411 CALL calc_scaling_parameters 412 CALL calc_surface_fluxes 413 IF ( do_output_at_2m ) THEN 414 CALL calc_pt_near_surface ( '2m' ) 415 ENDIF 503 416 ENDIF 417 ! 418 !-- Natural-type horizontal surfaces 419 IF ( surf_lsm_h%ns >= 1 ) THEN 420 surf => surf_lsm_h 421 CALL calc_uvw_abs 422 IF ( .NOT. neutral ) CALL calc_ol 423 CALL calc_us 424 CALL calc_scaling_parameters 425 CALL calc_surface_fluxes 426 IF ( do_output_at_2m ) THEN 427 CALL calc_pt_near_surface ( '2m' ) 428 ENDIF 429 ENDIF 430 ! 431 !-- Urban-type horizontal surfaces 432 IF ( surf_usm_h%ns >= 1 ) THEN 433 surf => surf_usm_h 434 CALL calc_uvw_abs 435 IF ( .NOT. neutral ) CALL calc_ol 436 CALL calc_us 437 CALL calc_scaling_parameters 438 CALL calc_surface_fluxes 439 IF ( do_output_at_2m ) THEN 440 CALL calc_pt_near_surface ( '2m' ) 441 ENDIF 442 IF ( indoor_model ) THEN 443 CALL calc_pt_near_surface ( '10cm' ) 444 ENDIF 445 ENDIF 446 504 447 ! 505 448 !-- Treat downward-facing horizontal surfaces. Note, so far, these are … … 549 492 !-- flux in land-surface model in case of aero_resist_kray is not true. 550 493 IF ( .NOT. aero_resist_kray ) THEN 551 IF ( most_method == 'circular' ) THEN 552 DO l = 0, 1 553 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 554 surf => surf_lsm_v(l) 555 ! 556 !-- Compute scaling parameters 557 CALL calc_scaling_parameters 558 ! 559 !-- Compute surface-parallel velocity 560 CALL calc_uvw_abs_v_ugrid 561 ! 562 !-- Compute Obukhov length 563 IF ( .NOT. neutral ) CALL calc_ol 564 ! 565 !-- Compute respective friction velocity on staggered grid 566 CALL calc_us 567 ! 568 !-- Compute respective surface fluxes for momentum and TKE 569 CALL calc_surface_fluxes 570 ENDIF 571 ENDDO 572 ELSE 573 DO l = 0, 1 574 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 575 surf => surf_lsm_v(l) 576 ! 577 !-- Compute surface-parallel velocity 578 CALL calc_uvw_abs_v_ugrid 579 ! 580 !-- Compute Obukhov length 581 IF ( .NOT. neutral ) CALL calc_ol 582 ! 583 !-- Compute respective friction velocity on staggered grid 584 CALL calc_us 585 ! 586 !-- Compute scaling parameters 587 CALL calc_scaling_parameters 588 ! 589 !-- Compute respective surface fluxes for momentum and TKE 590 CALL calc_surface_fluxes 591 ENDIF 592 ENDDO 593 ENDIF 494 DO l = 0, 1 495 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 496 surf => surf_lsm_v(l) 497 ! 498 !-- Compute surface-parallel velocity 499 CALL calc_uvw_abs_v_ugrid 500 ! 501 !-- Compute Obukhov length 502 IF ( .NOT. neutral ) CALL calc_ol 503 ! 504 !-- Compute respective friction velocity on staggered grid 505 CALL calc_us 506 ! 507 !-- Compute scaling parameters 508 CALL calc_scaling_parameters 509 ! 510 !-- Compute respective surface fluxes for momentum and TKE 511 CALL calc_surface_fluxes 512 ENDIF 513 ENDDO 594 514 ! 595 515 !-- No ts is required, so scaling parameters and Obukhov length do not need … … 652 572 !-- flux in land-surface model in case of aero_resist_kray is not true. 653 573 IF ( .NOT. aero_resist_kray ) THEN 654 IF ( most_method == 'circular' ) THEN 655 DO l = 2, 3 656 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 657 surf => surf_lsm_v(l) 658 ! 659 !-- Compute scaling parameters 660 CALL calc_scaling_parameters 661 ! 662 !-- Compute surface-parallel velocity 663 CALL calc_uvw_abs_v_vgrid 664 ! 665 !-- Compute Obukhov length 666 IF ( .NOT. neutral ) CALL calc_ol 667 ! 668 !-- Compute respective friction velocity on staggered grid 669 CALL calc_us 670 ! 671 !-- Compute respective surface fluxes for momentum and TKE 672 CALL calc_surface_fluxes 673 ENDIF 674 ENDDO 675 ELSE 676 DO l = 2, 3 677 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 678 surf => surf_lsm_v(l) 679 ! 680 !-- Compute surface-parallel velocity 681 CALL calc_uvw_abs_v_vgrid 682 ! 683 !-- Compute Obukhov length 684 IF ( .NOT. neutral ) CALL calc_ol 685 ! 686 !-- Compute respective friction velocity on staggered grid 687 CALL calc_us 688 ! 689 !-- Compute scaling parameters 690 CALL calc_scaling_parameters 691 ! 692 !-- Compute respective surface fluxes for momentum and TKE 693 CALL calc_surface_fluxes 694 ENDIF 695 ENDDO 696 ENDIF 574 DO l = 2, 3 575 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 576 surf => surf_lsm_v(l) 577 ! 578 !-- Compute surface-parallel velocity 579 CALL calc_uvw_abs_v_vgrid 580 ! 581 !-- Compute Obukhov length 582 IF ( .NOT. neutral ) CALL calc_ol 583 ! 584 !-- Compute respective friction velocity on staggered grid 585 CALL calc_us 586 ! 587 !-- Compute scaling parameters 588 CALL calc_scaling_parameters 589 ! 590 !-- Compute respective surface fluxes for momentum and TKE 591 CALL calc_surface_fluxes 592 ENDIF 593 ENDDO 697 594 ELSE 698 595 DO l = 2, 3 … … 841 738 IF ( surf_lsm_h%ns >= 1 ) surf_lsm_h%ol = 1.0E10_wp 842 739 IF ( surf_usm_h%ns >= 1 ) surf_usm_h%ol = 1.0E10_wp 843 ENDIF844 845 IF ( most_method == 'lookup' ) THEN846 847 !848 !-- Check for roughness heterogeneity. In that case terminate run and849 !-- inform user. Check for both, natural and non-natural walls.850 IF ( surf_def_h(0)%ns >= 1 ) THEN851 IF ( MINVAL( surf_def_h(0)%z0h ) /= MAXVAL( surf_def_h(0)%z0h ) .OR. &852 MINVAL( surf_def_h(0)%z0 ) /= MAXVAL( surf_def_h(0)%z0 ) ) THEN853 terminate_run_l = .TRUE.854 ENDIF855 ENDIF856 IF ( surf_lsm_h%ns >= 1 ) THEN857 IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h ) .OR. &858 MINVAL( surf_lsm_h%z0 ) /= MAXVAL( surf_lsm_h%z0 ) ) THEN859 terminate_run_l = .TRUE.860 ENDIF861 ENDIF862 IF ( surf_usm_h%ns >= 1 ) THEN863 IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_usm_h%z0h ) .OR. &864 MINVAL( surf_usm_h%z0 ) /= MAXVAL( surf_usm_h%z0 ) ) THEN865 terminate_run_l = .TRUE.866 ENDIF867 ENDIF868 !869 !-- Check roughness homogeneity between differt surface types.870 IF ( surf_lsm_h%ns >= 1 .AND. surf_def_h(0)%ns >= 1 ) THEN871 IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h ) .OR. &872 MINVAL( surf_lsm_h%z0 ) /= MAXVAL( surf_def_h(0)%z0 ) ) THEN873 terminate_run_l = .TRUE.874 ENDIF875 ENDIF876 IF ( surf_usm_h%ns >= 1 .AND. surf_def_h(0)%ns >= 1 ) THEN877 IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h ) .OR. &878 MINVAL( surf_usm_h%z0 ) /= MAXVAL( surf_def_h(0)%z0 ) ) THEN879 terminate_run_l = .TRUE.880 ENDIF881 ENDIF882 IF ( surf_usm_h%ns >= 1 .AND. surf_lsm_h%ns >= 1 ) THEN883 IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h ) .OR. &884 MINVAL( surf_usm_h%z0 ) /= MAXVAL( surf_lsm_h%z0 ) ) THEN885 terminate_run_l = .TRUE.886 ENDIF887 ENDIF888 889 890 #if defined( __parallel )891 !892 !-- Make a logical OR for all processes. Force termiation of model if result893 !-- is TRUE894 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )895 CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL, &896 MPI_LOR, comm2d, ierr )897 #else898 terminate_run = terminate_run_l899 #endif900 901 IF ( terminate_run ) THEN902 message_string = 'most_method = "lookup" cannot be used in ' // &903 'combination with a prescribed roughness ' // &904 'heterogeneity'905 CALL message( 'surface_layer_fluxes', 'PA0116', 1, 2, 0, 6, 0 )906 ENDIF907 908 ALLOCATE( zeta_tmp(0:num_steps-1) )909 910 !911 !-- Use the lowest possible value for z_mo912 k = nzb913 z_mo = zu(k+1) - zw(k)914 915 !916 !-- Calculate z/L range from zeta_stretch to zeta_max using 90% of the917 !-- available steps (num_steps). The calculation is done with negative918 !-- values of zeta in order to simplify the stretching in the free919 !-- convection limit for the remaining 10% of steps.920 zeta_tmp(0) = - zeta_max921 num_steps_n = ( num_steps * 9 / 10 ) - 1922 zeta_step = (zeta_max - zeta_stretch) / REAL(num_steps_n)923 924 DO li = 1, num_steps_n925 zeta_tmp(li) = zeta_tmp(li-1) + zeta_step926 ENDDO927 928 !929 !-- Calculate stretching factor for the free convection range930 DO WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )931 regr_old = regr932 regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr ) &933 )**( 10.0_wp / REAL(num_steps) )934 ENDDO935 936 !937 !-- Calculate z/L range from zeta_min to zeta_stretch938 DO li = num_steps_n+1, num_steps-1939 zeta_tmp(li) = zeta_tmp(li-1) + zeta_step940 zeta_step = zeta_step * regr941 ENDDO942 943 !944 !-- Save roughness lengths to temporary variables945 IF ( surf_def_h(0)%ns >= 1 ) THEN946 z0h_min = surf_def_h(0)%z0h(1)947 z0_min = surf_def_h(0)%z0(1)948 ELSEIF ( surf_lsm_h%ns >= 1 ) THEN949 z0h_min = surf_lsm_h%z0h(1)950 z0_min = surf_lsm_h%z0(1)951 ELSEIF ( surf_usm_h%ns >= 1 ) THEN952 z0h_min = surf_usm_h%z0h(1)953 z0_min = surf_usm_h%z0(1)954 ENDIF955 !956 !-- Calculate lookup table for the Richardson number versus Obukhov length957 !-- The Richardson number (rib) is defined depending on the choice of958 !-- boundary conditions for temperature959 IF ( ibc_pt_b == 1 ) THEN960 DO li = 0, num_steps-1961 ol_tab(li) = - z_mo / zeta_tmp(num_steps-1-li)962 rib_tab(li) = z_mo / ol_tab(li) / ( LOG( z_mo / z0_min ) &963 - psi_m( z_mo / ol_tab(li) ) &964 + psi_m( z0_min / ol_tab(li) ) &965 )**3966 ENDDO967 ELSE968 DO li = 0, num_steps-1969 ol_tab(li) = - z_mo / zeta_tmp(num_steps-1-li)970 rib_tab(li) = z_mo / ol_tab(li) * ( LOG( z_mo / z0h_min ) &971 - psi_h( z_mo / ol_tab(li) ) &972 + psi_h( z0h_min / ol_tab(li) ) &973 ) &974 / ( LOG( z_mo / z0_min ) &975 - psi_m( z_mo / ol_tab(li) ) &976 + psi_m( z0_min / ol_tab(li) ) &977 )**2978 ENDDO979 ENDIF980 981 !982 !-- Determine minimum values of rib in the lookup table. Set upper limit983 !-- to critical Richardson number (0.25)984 rib_min = MINVAL(rib_tab)985 rib_max = 0.25 !MAXVAL(rib_tab)986 987 DEALLOCATE( zeta_tmp )988 740 ENDIF 989 741 … … 1256 1008 ol_u !< Upper bound of L for Newton iteration 1257 1009 1258 IF ( TRIM( most_method ) /= 'circular' ) THEN 1259 ! 1260 !-- Evaluate bulk Richardson number (calculation depends on 1261 !-- definition based on setting of boundary conditions 1262 IF ( ibc_pt_b /= 1 ) THEN 1263 IF ( humidity ) THEN 1264 !$OMP PARALLEL DO PRIVATE( z_mo ) 1265 DO m = 1, surf%ns 1266 1267 z_mo = surf%z_mo(m) 1268 1269 surf%rib(m) = g * z_mo & 1270 * ( surf%vpt1(m) - surf%vpt_surface(m) ) & 1271 / ( surf%uvw_abs(m)**2 * surf%vpt1(m) & 1272 + 1.0E-20_wp ) 1273 ENDDO 1274 ELSE 1275 !$OMP PARALLEL DO PRIVATE( z_mo ) 1276 DO m = 1, surf%ns 1277 1278 z_mo = surf%z_mo(m) 1279 1280 surf%rib(m) = g * z_mo & 1281 * ( surf%pt1(m) - surf%pt_surface(m) ) & 1282 / ( surf%uvw_abs(m)**2 * surf%pt1(m) & 1283 + 1.0E-20_wp ) 1284 ENDDO 1285 ENDIF 1010 ! 1011 !-- Evaluate bulk Richardson number (calculation depends on 1012 !-- definition based on setting of boundary conditions 1013 IF ( ibc_pt_b /= 1 ) THEN 1014 IF ( humidity ) THEN 1015 !$OMP PARALLEL DO PRIVATE( z_mo ) 1016 DO m = 1, surf%ns 1017 1018 z_mo = surf%z_mo(m) 1019 1020 surf%rib(m) = g * z_mo & 1021 * ( surf%vpt1(m) - surf%vpt_surface(m) ) & 1022 / ( surf%uvw_abs(m)**2 * surf%vpt1(m) & 1023 + 1.0E-20_wp ) 1024 ENDDO 1286 1025 ELSE 1287 IF ( humidity ) THEN 1288 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 1289 DO m = 1, surf%ns 1290 1291 k = surf%k(m) 1292 1293 z_mo = surf%z_mo(m) 1294 1295 surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp & 1026 !$OMP PARALLEL DO PRIVATE( z_mo ) 1027 DO m = 1, surf%ns 1028 1029 z_mo = surf%z_mo(m) 1030 1031 surf%rib(m) = g * z_mo & 1032 * ( surf%pt1(m) - surf%pt_surface(m) ) & 1033 / ( surf%uvw_abs(m)**2 * surf%pt1(m) + 1.0E-20_wp ) 1034 ENDDO 1035 ENDIF 1036 ELSE 1037 IF ( humidity ) THEN 1038 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 1039 DO m = 1, surf%ns 1040 1041 k = surf%k(m) 1042 1043 z_mo = surf%z_mo(m) 1044 1045 surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp & 1296 1046 * surf%qv1(m) ) * surf%shf(m) + 0.61_wp & 1297 1047 * surf%pt1(m) * surf%qsws(m) ) * & … … 1299 1049 ( surf%uvw_abs(m)**3 * surf%vpt1(m) * kappa**2 & 1300 1050 + 1.0E-20_wp ) 1301 ENDDO 1051 ENDDO 1052 ELSE 1053 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 1054 !$ACC PARALLEL LOOP PRIVATE(k, z_mo) & 1055 !$ACC PRESENT(surf, drho_air_zw) 1056 DO m = 1, surf%ns 1057 1058 k = surf%k(m) 1059 1060 z_mo = surf%z_mo(m) 1061 1062 surf%rib(m) = - g * z_mo * surf%shf(m) * & 1063 drho_air_zw(k-1) / & 1064 ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2 & 1065 + 1.0E-20_wp ) 1066 ENDDO 1067 ENDIF 1068 ENDIF 1069 1070 ! 1071 !-- Calculate the Obukhov length using Newton iteration 1072 !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) & 1073 !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) & 1074 !$ACC PRESENT(surf) 1075 DO m = 1, surf%ns 1076 1077 i = surf%i(m) 1078 j = surf%j(m) 1079 1080 z_mo = surf%z_mo(m) 1081 1082 ! 1083 !-- Store current value in case the Newton iteration fails 1084 ol_old = surf%ol(m) 1085 1086 ! 1087 !-- Ensure that the bulk Richardson number and the Obukhov 1088 !-- length have the same sign 1089 IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. & 1090 ABS( surf%ol(m) ) == ol_max ) THEN 1091 IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp 1092 IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp 1093 ENDIF 1094 ! 1095 !-- Iteration to find Obukhov length 1096 iter = 0 1097 DO 1098 iter = iter + 1 1099 ! 1100 !-- In case of divergence, use the value of the previous time step 1101 IF ( iter > 1000 ) THEN 1102 surf%ol(m) = ol_old 1103 EXIT 1104 ENDIF 1105 1106 ol_m = surf%ol(m) 1107 ol_l = ol_m - 0.001_wp * ol_m 1108 ol_u = ol_m + 0.001_wp * ol_m 1109 1110 1111 IF ( ibc_pt_b /= 1 ) THEN 1112 ! 1113 !-- Calculate f = Ri - [...]/[...]^2 = 0 1114 f = surf%rib(m) - ( z_mo / ol_m ) * ( & 1115 LOG( z_mo / surf%z0h(m) ) & 1116 - psi_h( z_mo / ol_m ) & 1117 + psi_h( surf%z0h(m) / & 1118 ol_m ) & 1119 ) & 1120 / ( LOG( z_mo / surf%z0(m) ) & 1121 - psi_m( z_mo / ol_m ) & 1122 + psi_m( surf%z0(m) / ol_m ) & 1123 )**2 1124 1125 ! 1126 !-- Calculate df/dL 1127 f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / & 1128 surf%z0h(m) ) & 1129 - psi_h( z_mo / ol_u ) & 1130 + psi_h( surf%z0h(m) / ol_u ) & 1131 ) & 1132 / ( LOG( z_mo / surf%z0(m) ) & 1133 - psi_m( z_mo / ol_u ) & 1134 + psi_m( surf%z0(m) / ol_u ) & 1135 )**2 & 1136 + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) ) & 1137 - psi_h( z_mo / ol_l ) & 1138 + psi_h( surf%z0h(m) / ol_l ) & 1139 ) & 1140 / ( LOG( z_mo / surf%z0(m) ) & 1141 - psi_m( z_mo / ol_l ) & 1142 + psi_m( surf%z0(m) / ol_l ) & 1143 )**2 & 1144 ) / ( ol_u - ol_l ) 1302 1145 ELSE 1303 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 1304 !$ACC PARALLEL LOOP PRIVATE(k, z_mo) & 1305 !$ACC PRESENT(surf, drho_air_zw) 1306 DO m = 1, surf%ns 1307 1308 k = surf%k(m) 1309 1310 z_mo = surf%z_mo(m) 1311 1312 surf%rib(m) = - g * z_mo * surf%shf(m) * & 1313 drho_air_zw(k-1) / & 1314 ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2 & 1315 + 1.0E-20_wp ) 1316 ENDDO 1146 ! 1147 !-- Calculate f = Ri - 1 /[...]^3 = 0 1148 f = surf%rib(m) - ( z_mo / ol_m ) / & 1149 ( LOG( z_mo / surf%z0(m) ) & 1150 - psi_m( z_mo / ol_m ) & 1151 + psi_m( surf%z0(m) / ol_m ) & 1152 )**3 1153 1154 ! 1155 !-- Calculate df/dL 1156 f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) ) & 1157 - psi_m( z_mo / ol_u ) & 1158 + psi_m( surf%z0(m) / ol_u ) & 1159 )**3 & 1160 + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) ) & 1161 - psi_m( z_mo / ol_l ) & 1162 + psi_m( surf%z0(m) / ol_l ) & 1163 )**3 & 1164 ) / ( ol_u - ol_l ) 1317 1165 ENDIF 1318 ENDIF 1319 1320 ENDIF 1321 1322 1323 ! 1324 !-- Calculate the Obukhov length either using a Newton iteration 1325 !-- method, via a lookup table, or using the old circular way 1326 IF ( TRIM( most_method ) == 'newton' ) THEN 1327 1328 !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) & 1329 !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) & 1330 !$ACC PRESENT(surf) 1331 DO m = 1, surf%ns 1332 1333 i = surf%i(m) 1334 j = surf%j(m) 1335 1336 z_mo = surf%z_mo(m) 1337 1338 ! 1339 !-- Store current value in case the Newton iteration fails 1340 ol_old = surf%ol(m) 1166 ! 1167 !-- Calculate new L 1168 surf%ol(m) = ol_m - f / f_d_ol 1341 1169 1342 1170 ! 1343 1171 !-- Ensure that the bulk Richardson number and the Obukhov 1344 !-- length have the same sign 1345 IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. & 1346 ABS( surf%ol(m) ) == ol_max ) THEN 1347 IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp 1348 IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp 1172 !-- length have the same sign and ensure convergence. 1173 IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp 1174 ! 1175 !-- If unrealistic value occurs, set L to the maximum 1176 !-- value that is allowed 1177 IF ( ABS( surf%ol(m) ) > ol_max ) THEN 1178 surf%ol(m) = ol_max 1179 EXIT 1349 1180 ENDIF 1350 1181 ! 1351 !-- Iteration to find Obukhov length 1352 iter = 0 1353 DO 1354 iter = iter + 1 1355 ! 1356 !-- In case of divergence, use the value of the previous time step 1357 IF ( iter > 1000 ) THEN 1358 surf%ol(m) = ol_old 1359 EXIT 1360 ENDIF 1361 1362 ol_m = surf%ol(m) 1363 ol_l = ol_m - 0.001_wp * ol_m 1364 ol_u = ol_m + 0.001_wp * ol_m 1365 1366 1367 IF ( ibc_pt_b /= 1 ) THEN 1368 ! 1369 !-- Calculate f = Ri - [...]/[...]^2 = 0 1370 f = surf%rib(m) - ( z_mo / ol_m ) * ( & 1371 LOG( z_mo / surf%z0h(m) ) & 1372 - psi_h( z_mo / ol_m ) & 1373 + psi_h( surf%z0h(m) / & 1374 ol_m ) & 1375 ) & 1376 / ( LOG( z_mo / surf%z0(m) ) & 1377 - psi_m( z_mo / ol_m ) & 1378 + psi_m( surf%z0(m) / & 1379 ol_m ) & 1380 )**2 1381 1382 ! 1383 !-- Calculate df/dL 1384 f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / & 1385 surf%z0h(m) ) & 1386 - psi_h( z_mo / ol_u ) & 1387 + psi_h( surf%z0h(m) / ol_u ) & 1388 ) & 1389 / ( LOG( z_mo / surf%z0(m) ) & 1390 - psi_m( z_mo / ol_u ) & 1391 + psi_m( surf%z0(m) / ol_u ) & 1392 )**2 & 1393 + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) ) & 1394 - psi_h( z_mo / ol_l ) & 1395 + psi_h( surf%z0h(m) / ol_l ) & 1396 ) & 1397 / ( LOG( z_mo / surf%z0(m) ) & 1398 - psi_m( z_mo / ol_l ) & 1399 + psi_m( surf%z0(m) / ol_l ) & 1400 )**2 & 1401 ) / ( ol_u - ol_l ) 1402 ELSE 1403 ! 1404 !-- Calculate f = Ri - 1 /[...]^3 = 0 1405 f = surf%rib(m) - ( z_mo / ol_m ) / & 1406 ( LOG( z_mo / surf%z0(m) ) & 1407 - psi_m( z_mo / ol_m ) & 1408 + psi_m( surf%z0(m) / ol_m ) & 1409 )**3 1410 1411 ! 1412 !-- Calculate df/dL 1413 f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / & 1414 surf%z0(m) ) & 1415 - psi_m( z_mo / ol_u ) & 1416 + psi_m( surf%z0(m) / ol_u ) & 1417 )**3 & 1418 + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) ) & 1419 - psi_m( z_mo / ol_l ) & 1420 + psi_m( surf%z0(m) / ol_l ) & 1421 )**3 & 1422 ) / ( ol_u - ol_l ) 1423 ENDIF 1424 ! 1425 !-- Calculate new L 1426 surf%ol(m) = ol_m - f / f_d_ol 1427 1428 ! 1429 !-- Ensure that the bulk Richardson number and the Obukhov 1430 !-- length have the same sign and ensure convergence. 1431 IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp 1432 ! 1433 !-- If unrealistic value occurs, set L to the maximum 1434 !-- value that is allowed 1435 IF ( ABS( surf%ol(m) ) > ol_max ) THEN 1436 surf%ol(m) = ol_max 1437 EXIT 1438 ENDIF 1439 ! 1440 !-- Check for convergence 1441 IF ( ABS( ( surf%ol(m) - ol_m ) / & 1442 surf%ol(m) ) < 1.0E-4_wp ) THEN 1443 EXIT 1444 ELSE 1445 CYCLE 1446 ENDIF 1447 1448 ENDDO 1182 !-- Check for convergence 1183 IF ( ABS( ( surf%ol(m) - ol_m ) / surf%ol(m) ) < 1.0E-4_wp ) THEN 1184 EXIT 1185 ELSE 1186 CYCLE 1187 ENDIF 1188 1449 1189 ENDDO 1450 1451 ELSEIF ( TRIM( most_method ) == 'lookup' ) THEN 1452 1453 !$OMP PARALLEL DO PRIVATE( i, j, z_mo, li ) FIRSTPRIVATE( li_bnd ) LASTPRIVATE( li_bnd ) 1454 DO m = 1, surf%ns 1455 1456 i = surf%i(m) 1457 j = surf%j(m) 1458 ! 1459 !-- If the bulk Richardson number is outside the range of the lookup 1460 !-- table, set it to the exceeding threshold value 1461 IF ( surf%rib(m) < rib_min ) surf%rib(m) = rib_min 1462 IF ( surf%rib(m) > rib_max ) surf%rib(m) = rib_max 1463 1464 ! 1465 !-- Find the correct index bounds for linear interpolation. As the 1466 !-- Richardson number will not differ very much from time step to 1467 !-- time step , use the index from the last step and search in the 1468 !-- correct direction 1469 li = li_bnd 1470 IF ( rib_tab(li) - surf%rib(m) > 0.0_wp ) THEN 1471 DO WHILE ( rib_tab(li-1) - surf%rib(m) > 0.0_wp .AND. li > 0 ) 1472 li = li-1 1473 ENDDO 1474 ELSE 1475 DO WHILE ( rib_tab(li) - surf%rib(m) < 0.0_wp & 1476 .AND. li < num_steps-1 ) 1477 li = li+1 1478 ENDDO 1479 ENDIF 1480 li_bnd = li 1481 1482 ! 1483 !-- Linear interpolation to find the correct value of z/L 1484 surf%ol(m) = ( ol_tab(li-1) + ( ol_tab(li) - ol_tab(li-1) ) & 1485 / ( rib_tab(li) - rib_tab(li-1) ) & 1486 * ( surf%rib(m) - rib_tab(li-1) ) ) 1487 1488 ENDDO 1489 1490 ELSEIF ( TRIM( most_method ) == 'circular' ) THEN 1491 1492 IF ( .NOT. humidity ) THEN 1493 !$OMP PARALLEL DO PRIVATE( z_mo ) 1494 DO m = 1, surf%ns 1495 1496 z_mo = surf%z_mo(m) 1497 1498 surf%ol(m) = ( surf%pt1(m) * surf%us(m)**2 ) / & 1499 ( kappa * g * & 1500 surf%ts(m) + 1E-30_wp ) 1501 ! 1502 !-- Limit the value range of the Obukhov length. 1503 !-- This is necessary for very small velocities (u,v --> 1), because 1504 !-- the absolute value of ol can then become very small, which in 1505 !-- consequence would result in very large shear stresses and very 1506 !-- small momentum fluxes (both are generally unrealistic). 1507 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min ) & 1508 surf%ol(m) = z_mo / zeta_min 1509 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max ) & 1510 surf%ol(m) = z_mo / zeta_max 1511 1512 ENDDO 1513 ELSEIF ( bulk_cloud_model .OR. cloud_droplets ) THEN 1514 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1515 DO m = 1, surf%ns 1516 1517 i = surf%i(m) 1518 j = surf%j(m) 1519 k = surf%k(m) 1520 1521 z_mo = surf%z_mo(m) 1522 1523 1524 surf%ol(m) = ( surf%vpt1(m) * surf%us(m)**2 ) / & 1525 ( kappa * g * ( surf%ts(m) + & 1526 0.61_wp * surf%pt1(m) * surf%us(m) & 1527 + 0.61_wp * surf%qv1(m) * surf%ts(m) - & 1528 surf%ts(m) * ql(k,j,i) ) + 1E-30_wp ) 1529 ! 1530 !-- Limit the value range of the Obukhov length. 1531 !-- This is necessary for very small velocities (u,v --> 1), because 1532 !-- the absolute value of ol can then become very small, which in 1533 !-- consequence would result in very large shear stresses and very 1534 !-- small momentum fluxes (both are generally unrealistic). 1535 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min ) & 1536 surf%ol(m) = z_mo / zeta_min 1537 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max ) & 1538 surf%ol(m) = z_mo / zeta_max 1539 1540 ENDDO 1541 ELSE 1542 1543 !$OMP PARALLEL DO PRIVATE( z_mo ) 1544 DO m = 1, surf%ns 1545 1546 z_mo = surf%z_mo(m) 1547 1548 surf%ol(m) = ( surf%vpt1(m) * surf%us(m)**2 ) / & 1549 ( kappa * g * ( surf%ts(m) + 0.61_wp * surf%pt1(m) * & 1550 surf%qs(m) + 0.61_wp * surf%qv1(m) * & 1551 surf%ts(m) ) + 1E-30_wp ) 1552 1553 ! 1554 !-- Limit the value range of the Obukhov length. 1555 !-- This is necessary for very small velocities (u,v --> 1), because 1556 !-- the absolute value of ol can then become very small, which in 1557 !-- consequence would result in very large shear stresses and very 1558 !-- small momentum fluxes (both are generally unrealistic). 1559 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min ) & 1560 surf%ol(m) = z_mo / zeta_min 1561 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max ) & 1562 surf%ol(m) = z_mo / zeta_max 1563 1564 ENDDO 1565 1566 ENDIF 1567 1568 ENDIF 1190 ENDDO 1569 1191 1570 1192 END SUBROUTINE calc_ol -
palm/trunk/SOURCE/write_restart_data_mod.f90
r3655 r3668 163 163 164 164 165 binary_version_global = '4. 7'165 binary_version_global = '4.8' 166 166 167 167 CALL wrd_write_string( 'binary_version_global' ) … … 466 466 WRITE ( 14 ) momentum_advec 467 467 468 CALL wrd_write_string( 'most_method' )469 WRITE ( 14 ) most_method470 471 468 CALL wrd_write_string( 'netcdf_precision' ) 472 469 WRITE ( 14 ) netcdf_precision
Note: See TracChangeset
for help on using the changeset viewer.
