Changeset 1788 for palm/trunk/SOURCE
- Timestamp:
- Mar 10, 2016 11:01:04 AM (9 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 12 edited
-
check_parameters.f90 (modified) (203 diffs)
-
data_output_2d.f90 (modified) (3 diffs)
-
header.f90 (modified) (4 diffs)
-
init_3d_model.f90 (modified) (48 diffs)
-
land_surface_model.f90 (modified) (46 diffs)
-
modules.f90 (modified) (3 diffs)
-
package_parin.f90 (modified) (4 diffs)
-
radiation_model.f90 (modified) (9 diffs)
-
read_3d_binary.f90 (modified) (4 diffs)
-
sum_up_3d_data.f90 (modified) (5 diffs)
-
surface_layer_fluxes.f90 (modified) (20 diffs)
-
write_3d_binary.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/check_parameters.f90
r1787 r1788 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 22 ! 21 ! Added check for use of most_method = 'lookup' in combination with water 22 ! surface presribed in the land surface model. Added output of z0q. 23 ! Syntax layout improved. 24 ! 23 25 ! Former revisions: 24 26 ! ----------------- … … 397 399 398 400 IF ( dt_coupling == 9999999.9_wp ) THEN 399 message_string = 'dt_coupling is not set but required for coup' // &401 message_string = 'dt_coupling is not set but required for coup' // & 400 402 'ling mode "' // TRIM( coupling_mode ) // '"' 401 403 CALL message( 'check_parameters', 'PA0003', 1, 2, 0, 6, 0 ) … … 419 421 #if ! defined( __check ) 420 422 IF ( myid == 0 ) THEN 421 CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, &423 CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, & 422 424 ierr ) 423 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, &425 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, & 424 426 status, ierr ) 425 427 ENDIF … … 436 438 IF ( myid == 0 ) THEN 437 439 CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr ) 438 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, &440 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, & 439 441 status, ierr ) 440 442 ENDIF … … 451 453 CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, & 452 454 ierr ) 453 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, &455 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, & 454 456 status, ierr ) 455 457 ENDIF … … 464 466 #if ! defined( __check ) 465 467 IF ( myid == 0 ) THEN 466 CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, &468 CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, & 467 469 ierr ) 468 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, &470 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, & 469 471 status, ierr ) 470 472 ENDIF … … 483 485 CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, & 484 486 14, comm_inter, ierr ) 485 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, &487 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, & 486 488 status, ierr ) 487 489 ENDIF … … 499 501 IF ( myid == 0 ) THEN 500 502 CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr ) 501 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, &503 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, & 502 504 status, ierr ) 503 505 ENDIF … … 508 510 509 511 IF ( dx < remote ) THEN 510 WRITE( message_string, * ) 'coupling mode "', &511 TRIM( coupling_mode ), &512 WRITE( message_string, * ) 'coupling mode "', & 513 TRIM( coupling_mode ), & 512 514 '": dx in Atmosphere is not equal to or not larger then dx in ocean' 513 515 CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 ) … … 515 517 516 518 IF ( (nx_a+1)*dx /= (nx_o+1)*remote ) THEN 517 WRITE( message_string, * ) 'coupling mode "', &518 TRIM( coupling_mode ), &519 WRITE( message_string, * ) 'coupling mode "', & 520 TRIM( coupling_mode ), & 519 521 '": Domain size in x-direction is not equal in ocean and atmosphere' 520 522 CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 ) … … 526 528 IF ( myid == 0) THEN 527 529 CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr ) 528 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, &530 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, & 529 531 status, ierr ) 530 532 ENDIF … … 534 536 535 537 IF ( dy < remote ) THEN 536 WRITE( message_string, * ) 'coupling mode "', &537 TRIM( coupling_mode ), &538 WRITE( message_string, * ) 'coupling mode "', & 539 TRIM( coupling_mode ), & 538 540 '": dy in Atmosphere is not equal to or not larger then dy in ocean' 539 541 CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 ) … … 541 543 542 544 IF ( (ny_a+1)*dy /= (ny_o+1)*remote ) THEN 543 WRITE( message_string, * ) 'coupling mode "', &544 TRIM( coupling_mode ), &545 WRITE( message_string, * ) 'coupling mode "', & 546 TRIM( coupling_mode ), & 545 547 '": Domain size in y-direction is not equal in ocean and atmosphere' 546 548 CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 ) … … 548 550 549 551 IF ( MOD(nx_o+1,nx_a+1) /= 0 ) THEN 550 WRITE( message_string, * ) 'coupling mode "', &551 TRIM( coupling_mode ), &552 WRITE( message_string, * ) 'coupling mode "', & 553 TRIM( coupling_mode ), & 552 554 '": nx+1 in ocean is not divisible without remainder with nx+1 in', & 553 555 ' atmosphere' … … 556 558 557 559 IF ( MOD(ny_o+1,ny_a+1) /= 0 ) THEN 558 WRITE( message_string, * ) 'coupling mode "', &559 TRIM( coupling_mode ), &560 WRITE( message_string, * ) 'coupling mode "', & 561 TRIM( coupling_mode ), & 560 562 '": ny+1 in ocean is not divisible without remainder with ny+1 in', & 561 563 ' atmosphere' … … 565 567 ENDIF 566 568 #else 567 WRITE( message_string, * ) 'coupling requires PALM to be called with', &569 WRITE( message_string, * ) 'coupling requires PALM to be called with', & 568 570 ' ''mrun -K parallel''' 569 571 CALL message( 'check_parameters', 'PA0141', 1, 2, 0, 6, 0 ) … … 575 577 !-- Exchange via intercommunicator 576 578 IF ( coupling_mode == 'atmosphere_to_ocean' .AND. myid == 0 ) THEN 577 CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter, &579 CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter, & 578 580 ierr ) 579 581 ELSEIF ( coupling_mode == 'ocean_to_atmosphere' .AND. myid == 0) THEN 580 CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19, &582 CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19, & 581 583 comm_inter, status, ierr ) 582 584 ENDIF … … 633 635 634 636 CASE DEFAULT 635 message_string = 'illegal value given for loop_optimization: "' // &637 message_string = 'illegal value given for loop_optimization: "' // & 636 638 TRIM( loop_optimization ) // '"' 637 639 CALL message( 'check_parameters', 'PA0013', 1, 2, 0, 6, 0 ) … … 659 661 IF ( topography /= 'flat' ) THEN 660 662 action = ' ' 661 IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme' &663 IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme' & 662 664 .AND. scalar_advec /= 'ws-scheme-mono' ) THEN 663 665 WRITE( action, '(A,A)' ) 'scalar_advec = ', scalar_advec 664 666 ENDIF 665 IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' ) 667 IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' )& 666 668 THEN 667 669 WRITE( action, '(A,A)' ) 'momentum_advec = ', momentum_advec … … 686 688 ENDIF 687 689 IF ( action /= ' ' ) THEN 688 message_string = 'a non-flat topography does not allow ' // &690 message_string = 'a non-flat topography does not allow ' // & 689 691 TRIM( action ) 690 692 CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 ) … … 695 697 !-- is applicable. If this is not possible, abort. 696 698 IF ( TRIM( topography_grid_convention ) == ' ' ) THEN 697 IF ( TRIM( topography ) /= 'single_building' .AND. &698 TRIM( topography ) /= 'single_street_canyon' .AND. &699 IF ( TRIM( topography ) /= 'single_building' .AND. & 700 TRIM( topography ) /= 'single_street_canyon' .AND. & 699 701 TRIM( topography ) /= 'read_from_file' ) THEN 700 702 !-- The default value is not applicable here, because it is only valid 701 703 !-- for the two standard cases 'single_building' and 'read_from_file' 702 704 !-- defined in init_grid. 703 WRITE( message_string, * ) &704 'The value for "topography_grid_convention" ', &705 'is not set. Its default value is & only valid for ', &706 '"topography" = ''single_building'', ', &707 '''single_street_canyon'' & or ''read_from_file''.', &705 WRITE( message_string, * ) & 706 'The value for "topography_grid_convention" ', & 707 'is not set. Its default value is & only valid for ', & 708 '"topography" = ''single_building'', ', & 709 '''single_street_canyon'' & or ''read_from_file''.', & 708 710 ' & Choose ''cell_edge'' or ''cell_center''.' 709 711 CALL message( 'user_check_parameters', 'PA0239', 1, 2, 0, 6, 0 ) … … 711 713 !-- The default value is applicable here. 712 714 !-- Set convention according to topography. 713 IF ( TRIM( topography ) == 'single_building' .OR. &715 IF ( TRIM( topography ) == 'single_building' .OR. & 714 716 TRIM( topography ) == 'single_street_canyon' ) THEN 715 717 topography_grid_convention = 'cell_edge' … … 718 720 ENDIF 719 721 ENDIF 720 ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND. &722 ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND. & 721 723 TRIM( topography_grid_convention ) /= 'cell_center' ) THEN 722 WRITE( message_string, * ) &723 'The value for "topography_grid_convention" is ', &724 WRITE( message_string, * ) & 725 'The value for "topography_grid_convention" is ', & 724 726 'not recognized. & Choose ''cell_edge'' or ''cell_center''.' 725 727 CALL message( 'user_check_parameters', 'PA0240', 1, 2, 0, 6, 0 ) … … 738 740 ENDIF 739 741 740 ELSEIF ( TRIM( coupling_mode ) == 'uncoupled' .AND. &742 ELSEIF ( TRIM( coupling_mode ) == 'uncoupled' .AND. & 741 743 TRIM( coupling_char ) == '_O' ) THEN 742 744 … … 744 746 !-- Check whether an (uncoupled) atmospheric run has been declared as an 745 747 !-- ocean run (this setting is done via mrun-option -y) 746 747 message_string = 'ocean = .F. does not allow coupling_char = "' // & 748 message_string = 'ocean = .F. does not allow coupling_char = "' // & 748 749 TRIM( coupling_char ) // '" set by mrun-option "-y"' 749 750 CALL message( 'check_parameters', 'PA0317', 1, 2, 0, 6, 0 ) … … 777 778 gamma_mg = 1 778 779 ELSE 779 message_string = 'unknown multigrid cycle: cycle_mg = "' // &780 message_string = 'unknown multigrid cycle: cycle_mg = "' // & 780 781 TRIM( cycle_mg ) // '"' 781 782 CALL message( 'check_parameters', 'PA0020', 1, 2, 0, 6, 0 ) … … 783 784 ENDIF 784 785 785 IF ( fft_method /= 'singleton-algorithm' .AND. &786 fft_method /= 'temperton-algorithm' .AND. &787 fft_method /= 'fftw' .AND. &786 IF ( fft_method /= 'singleton-algorithm' .AND. & 787 fft_method /= 'temperton-algorithm' .AND. & 788 fft_method /= 'fftw' .AND. & 788 789 fft_method /= 'system-specific' ) THEN 789 message_string = 'unknown fft-algorithm: fft_method = "' // &790 message_string = 'unknown fft-algorithm: fft_method = "' // & 790 791 TRIM( fft_method ) // '"' 791 792 CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 ) 792 793 ENDIF 793 794 794 IF( momentum_advec == 'ws-scheme' .AND. &795 IF( momentum_advec == 'ws-scheme' .AND. & 795 796 .NOT. call_psolver_at_all_substeps ) THEN 796 message_string = 'psolver must be called at each RK3 substep when "'// &797 message_string = 'psolver must be called at each RK3 substep when "'// & 797 798 TRIM(momentum_advec) // ' "is used for momentum_advec' 798 799 CALL message( 'check_parameters', 'PA0344', 1, 2, 0, 6, 0 ) … … 802 803 IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' ) & 803 804 THEN 804 message_string = 'unknown advection scheme: momentum_advec = "' // &805 message_string = 'unknown advection scheme: momentum_advec = "' // & 805 806 TRIM( momentum_advec ) // '"' 806 807 CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 ) … … 811 812 timestep_scheme == 'runge-kutta-2' ) ) & 812 813 THEN 813 message_string = 'momentum_advec or scalar_advec = "' &814 message_string = 'momentum_advec or scalar_advec = "' & 814 815 // TRIM( momentum_advec ) // '" is not allowed with timestep_scheme = "' // & 815 816 TRIM( timestep_scheme ) // '"' … … 819 820 scalar_advec /= 'ws-scheme-mono' .AND. scalar_advec /= 'bc-scheme' ) & 820 821 THEN 821 message_string = 'unknown advection scheme: scalar_advec = "' // &822 message_string = 'unknown advection scheme: scalar_advec = "' // & 822 823 TRIM( scalar_advec ) // '"' 823 824 CALL message( 'check_parameters', 'PA0024', 1, 2, 0, 6, 0 ) … … 825 826 IF ( scalar_advec == 'bc-scheme' .AND. loop_optimization == 'cache' ) & 826 827 THEN 827 message_string = 'advection_scheme scalar_advec = "' &828 message_string = 'advection_scheme scalar_advec = "' & 828 829 // TRIM( scalar_advec ) // '" not implemented for & loop_optimization = "' // & 829 830 TRIM( loop_optimization ) // '"' … … 851 852 !-- Set LOGICAL switches to enhance performance 852 853 IF ( momentum_advec == 'ws-scheme' ) ws_scheme_mom = .TRUE. 853 IF ( scalar_advec == 'ws-scheme' .OR. &854 IF ( scalar_advec == 'ws-scheme' .OR. & 854 855 scalar_advec == 'ws-scheme-mono' ) ws_scheme_sca = .TRUE. 855 856 IF ( scalar_advec == 'ws-scheme-mono' ) monotonic_adjustment = .TRUE. … … 908 909 IF ( collision_kernel(6:9) == 'fast' ) use_kernel_tables = .TRUE. 909 910 910 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &911 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. & 911 912 TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN 912 913 ! 913 914 !-- No restart run: several initialising actions are possible 914 915 action = initializing_actions 915 DO WHILE ( TRIM( action ) /= '' )916 DO WHILE ( TRIM( action ) /= '' ) 916 917 position = INDEX( action, ' ' ) 917 918 SELECT CASE ( action(1:position-1) ) 918 919 919 CASE ( 'set_constant_profiles', 'set_1d-model_profiles', &920 CASE ( 'set_constant_profiles', 'set_1d-model_profiles', & 920 921 'by_user', 'initialize_vortex', 'initialize_ptanom' ) 921 922 action = action(position+1:) 922 923 923 924 CASE DEFAULT 924 message_string = 'initializing_action = "' // &925 message_string = 'initializing_action = "' // & 925 926 TRIM( action ) // '" unkown or not allowed' 926 927 CALL message( 'check_parameters', 'PA0030', 1, 2, 0, 6, 0 ) … … 930 931 ENDIF 931 932 932 IF ( TRIM( initializing_actions ) == 'initialize_vortex' .AND.&933 IF ( TRIM( initializing_actions ) == 'initialize_vortex' .AND. & 933 934 conserve_volume_flow ) THEN 934 message_string = 'initializing_actions = "initialize_vortex"' // &935 message_string = 'initializing_actions = "initialize_vortex"' // & 935 936 ' ist not allowed with conserve_volume_flow = .T.' 936 937 CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 ) … … 938 939 939 940 940 IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .AND. &941 IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .AND. & 941 942 INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN 942 message_string = 'initializing_actions = "set_constant_profiles"' // &943 ' and "set_1d-model_profiles" are not allowed ' // &943 message_string = 'initializing_actions = "set_constant_profiles"' // & 944 ' and "set_1d-model_profiles" are not allowed ' // & 944 945 'simultaneously' 945 946 CALL message( 'check_parameters', 'PA0031', 1, 2, 0, 6, 0 ) 946 947 ENDIF 947 948 948 IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .AND. &949 IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .AND. & 949 950 INDEX( initializing_actions, 'by_user' ) /= 0 ) THEN 950 message_string = 'initializing_actions = "set_constant_profiles"' // &951 message_string = 'initializing_actions = "set_constant_profiles"' // & 951 952 ' and "by_user" are not allowed simultaneously' 952 953 CALL message( 'check_parameters', 'PA0032', 1, 2, 0, 6, 0 ) 953 954 ENDIF 954 955 955 IF ( INDEX( initializing_actions, 'by_user' ) /= 0 .AND. &956 IF ( INDEX( initializing_actions, 'by_user' ) /= 0 .AND. & 956 957 INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN 957 message_string = 'initializing_actions = "by_user" and ' // &958 message_string = 'initializing_actions = "by_user" and ' // & 958 959 '"set_1d-model_profiles" are not allowed simultaneously' 959 960 CALL message( 'check_parameters', 'PA0033', 1, 2, 0, 6, 0 ) 960 961 ENDIF 961 962 962 IF ( cloud_physics .AND. .NOT. humidity ) THEN963 WRITE( message_string, * ) 'cloud_physics = ', cloud_physics, ' is ', &963 IF ( cloud_physics .AND. .NOT. humidity ) THEN 964 WRITE( message_string, * ) 'cloud_physics = ', cloud_physics, ' is ', & 964 965 'not allowed with humidity = ', humidity 965 966 CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 ) … … 967 968 968 969 IF ( precipitation .AND. .NOT. cloud_physics ) THEN 969 WRITE( message_string, * ) 'precipitation = ', precipitation, ' is ', &970 WRITE( message_string, * ) 'precipitation = ', precipitation, ' is ', & 970 971 'not allowed with cloud_physics = ', cloud_physics 971 972 CALL message( 'check_parameters', 'PA0035', 1, 2, 0, 6, 0 ) … … 973 974 974 975 IF ( humidity .AND. sloping_surface ) THEN 975 message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' // &976 message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' // & 976 977 'are not allowed simultaneously' 977 978 CALL message( 'check_parameters', 'PA0036', 1, 2, 0, 6, 0 ) … … 979 980 980 981 IF ( passive_scalar .AND. humidity ) THEN 981 message_string = 'humidity = .TRUE. and passive_scalar = .TRUE. ' // &982 message_string = 'humidity = .TRUE. and passive_scalar = .TRUE. ' // & 982 983 'is not allowed simultaneously' 983 984 CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 ) … … 1032 1033 !-- calculated from the temperature/humidity gradients in the land surface 1033 1034 !-- model 1034 IF ( bc_pt_b == 'neumann' .OR.bc_q_b == 'neumann' ) THEN1035 IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' ) THEN 1035 1036 message_string = 'lsm requires setting of'// & 1036 1037 'bc_pt_b = "dirichlet" and '// & … … 1039 1040 ENDIF 1040 1041 1041 IF ( .NOT.constant_flux_layer ) THEN1042 IF ( .NOT. constant_flux_layer ) THEN 1042 1043 message_string = 'lsm requires '// & 1043 1044 'constant_flux_layer = .T.' … … 1046 1047 1047 1048 IF ( topography /= 'flat' ) THEN 1048 message_string = 'lsm cannot be used ' // &1049 message_string = 'lsm cannot be used ' // & 1049 1050 'in combination with topography /= "flat"' 1050 1051 CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 ) 1051 1052 ENDIF 1052 1053 1054 IF ( ( veg_type == 14 .OR. veg_type == 15 ) .AND. & 1055 most_method == 'lookup' ) THEN 1056 WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ', & 1057 'allowed in combination with ', & 1058 'most_method = ', most_method 1059 CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 ) 1060 ENDIF 1061 1053 1062 IF ( veg_type == 0 ) THEN 1054 IF ( SUM( root_fraction) /= 1.0_wp) THEN1063 IF ( SUM( root_fraction ) /= 1.0_wp ) THEN 1055 1064 message_string = 'veg_type = 0 (user_defined)'// & 1056 1065 'requires setting of root_fraction(0:3)'// & … … 1059 1068 ENDIF 1060 1069 1061 IF ( min_canopy_resistance == 9999999.9_wp ) THEN1070 IF ( min_canopy_resistance == 9999999.9_wp ) THEN 1062 1071 message_string = 'veg_type = 0 (user defined)'// & 1063 1072 'requires setting of min_canopy_resistance'// & … … 1066 1075 ENDIF 1067 1076 1068 IF ( leaf_area_index == 9999999.9_wp ) THEN1077 IF ( leaf_area_index == 9999999.9_wp ) THEN 1069 1078 message_string = 'veg_type = 0 (user_defined)'// & 1070 1079 'requires setting of leaf_area_index'// & … … 1073 1082 ENDIF 1074 1083 1075 IF ( vegetation_coverage == 9999999.9_wp ) THEN1084 IF ( vegetation_coverage == 9999999.9_wp ) THEN 1076 1085 message_string = 'veg_type = 0 (user_defined)'// & 1077 1086 'requires setting of vegetation_coverage'// & … … 1087 1096 ENDIF 1088 1097 1089 IF ( lambda_surface_stable == 9999999.9_wp ) THEN1098 IF ( lambda_surface_stable == 9999999.9_wp ) THEN 1090 1099 message_string = 'veg_type = 0 (user_defined)'// & 1091 1100 'requires setting of lambda_surface_stable'// & … … 1094 1103 ENDIF 1095 1104 1096 IF ( lambda_surface_unstable == 9999999.9_wp ) THEN1105 IF ( lambda_surface_unstable == 9999999.9_wp ) THEN 1097 1106 message_string = 'veg_type = 0 (user_defined)'// & 1098 1107 'requires setting of lambda_surface_unstable'// & … … 1101 1110 ENDIF 1102 1111 1103 IF ( f_shortwave_incoming == 9999999.9_wp ) THEN1112 IF ( f_shortwave_incoming == 9999999.9_wp ) THEN 1104 1113 message_string = 'veg_type = 0 (user_defined)'// & 1105 1114 'requires setting of f_shortwave_incoming'// & … … 1108 1117 ENDIF 1109 1118 1110 IF ( z0_eb == 9999999.9_wp ) THEN1119 IF ( z0_eb == 9999999.9_wp ) THEN 1111 1120 message_string = 'veg_type = 0 (user_defined)'// & 1112 1121 'requires setting of z0_eb'// & … … 1115 1124 ENDIF 1116 1125 1117 IF ( z0h_eb == 9999999.9_wp ) THEN1126 IF ( z0h_eb == 9999999.9_wp ) THEN 1118 1127 message_string = 'veg_type = 0 (user_defined)'// & 1119 1128 'requires setting of z0h_eb'// & … … 1127 1136 IF ( soil_type == 0 ) THEN 1128 1137 1129 IF ( alpha_vangenuchten == 9999999.9_wp ) THEN1138 IF ( alpha_vangenuchten == 9999999.9_wp ) THEN 1130 1139 message_string = 'soil_type = 0 (user_defined)'// & 1131 1140 'requires setting of alpha_vangenuchten'// & … … 1134 1143 ENDIF 1135 1144 1136 IF ( l_vangenuchten == 9999999.9_wp ) THEN1145 IF ( l_vangenuchten == 9999999.9_wp ) THEN 1137 1146 message_string = 'soil_type = 0 (user_defined)'// & 1138 1147 'requires setting of l_vangenuchten'// & … … 1141 1150 ENDIF 1142 1151 1143 IF ( n_vangenuchten == 9999999.9_wp ) THEN1152 IF ( n_vangenuchten == 9999999.9_wp ) THEN 1144 1153 message_string = 'soil_type = 0 (user_defined)'// & 1145 1154 'requires setting of n_vangenuchten'// & … … 1148 1157 ENDIF 1149 1158 1150 IF ( hydraulic_conductivity == 9999999.9_wp ) THEN1159 IF ( hydraulic_conductivity == 9999999.9_wp ) THEN 1151 1160 message_string = 'soil_type = 0 (user_defined)'// & 1152 1161 'requires setting of hydraulic_conductivity'// & … … 1155 1164 ENDIF 1156 1165 1157 IF ( saturation_moisture == 9999999.9_wp ) THEN1166 IF ( saturation_moisture == 9999999.9_wp ) THEN 1158 1167 message_string = 'soil_type = 0 (user_defined)'// & 1159 1168 'requires setting of saturation_moisture'// & … … 1162 1171 ENDIF 1163 1172 1164 IF ( field_capacity == 9999999.9_wp ) THEN1173 IF ( field_capacity == 9999999.9_wp ) THEN 1165 1174 message_string = 'soil_type = 0 (user_defined)'// & 1166 1175 'requires setting of field_capacity'// & … … 1169 1178 ENDIF 1170 1179 1171 IF ( wilting_point == 9999999.9_wp ) THEN1180 IF ( wilting_point == 9999999.9_wp ) THEN 1172 1181 message_string = 'soil_type = 0 (user_defined)'// & 1173 1182 'requires setting of wilting_point'// & … … 1176 1185 ENDIF 1177 1186 1178 IF ( residual_moisture == 9999999.9_wp ) THEN1187 IF ( residual_moisture == 9999999.9_wp ) THEN 1179 1188 message_string = 'soil_type = 0 (user_defined)'// & 1180 1189 'requires setting of residual_moisture'// & … … 1185 1194 ENDIF 1186 1195 1187 IF ( .NOT.radiation ) THEN1196 IF ( .NOT. radiation ) THEN 1188 1197 message_string = 'lsm requires '// & 1189 1198 'radiation = .T.' … … 1194 1203 1195 1204 IF ( radiation ) THEN 1196 IF ( radiation_scheme /= 'constant' .AND.&1197 radiation_scheme /= 'clear-sky' .AND.&1205 IF ( radiation_scheme /= 'constant' .AND. & 1206 radiation_scheme /= 'clear-sky' .AND. & 1198 1207 radiation_scheme /= 'rrtmg' ) THEN 1199 1208 message_string = 'unknown radiation_scheme = '// & … … 1215 1224 1216 1225 ENDIF 1217 IF ( albedo_type == 0 .AND. albedo == 9999999.9_wp .AND.&1226 IF ( albedo_type == 0 .AND. albedo == 9999999.9_wp .AND. & 1218 1227 radiation_scheme == 'clear-sky') THEN 1219 1228 message_string = 'radiation_scheme = "clear-sky" in combination' // & … … 1222 1231 CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 ) 1223 1232 ENDIF 1224 IF ( albedo_type == 0 .AND. radiation_scheme == 'rrtmg' .AND.&1233 IF ( albedo_type == 0 .AND. radiation_scheme == 'rrtmg' .AND. & 1225 1234 ( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp& 1226 1235 .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp& … … 1235 1244 ENDIF 1236 1245 IF ( topography /= 'flat' ) THEN 1237 message_string = 'radiation scheme cannot be used ' // &1246 message_string = 'radiation scheme cannot be used ' // & 1238 1247 'in combination with topography /= "flat"' 1239 1248 CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 ) … … 1244 1253 loop_optimization == 'vector' ) & 1245 1254 .AND. cloud_physics .AND. icloud_scheme == 0 ) THEN 1246 message_string = 'cloud_scheme = seifert_beheng requires ' // &1255 message_string = 'cloud_scheme = seifert_beheng requires ' // & 1247 1256 'loop_optimization = "cache" or "vector"' 1248 1257 CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 ) … … 1270 1279 gradient = 0.0_wp 1271 1280 1272 IF ( .NOT.ocean ) THEN1281 IF ( .NOT. ocean ) THEN 1273 1282 1274 1283 ug_vertical_gradient_level_ind(1) = 0 1275 1284 ug(0) = ug_surface 1276 1285 DO k = 1, nzt+1 1277 IF ( i < 11 ) THEN1278 IF ( ug_vertical_gradient_level(i) < zu(k) .AND. &1286 IF ( i < 11 ) THEN 1287 IF ( ug_vertical_gradient_level(i) < zu(k) .AND. & 1279 1288 ug_vertical_gradient_level(i) >= 0.0_wp ) THEN 1280 1289 gradient = ug_vertical_gradient(i) / 100.0_wp … … 1299 1308 ug(nzt+1) = ug_surface 1300 1309 DO k = nzt, nzb, -1 1301 IF ( i < 11 ) THEN1302 IF ( ug_vertical_gradient_level(i) > zu(k) .AND. &1310 IF ( i < 11 ) THEN 1311 IF ( ug_vertical_gradient_level(i) > zu(k) .AND. & 1303 1312 ug_vertical_gradient_level(i) <= 0.0_wp ) THEN 1304 1313 gradient = ug_vertical_gradient(i) / 100.0_wp … … 1334 1343 gradient = 0.0_wp 1335 1344 1336 IF ( .NOT.ocean ) THEN1345 IF ( .NOT. ocean ) THEN 1337 1346 1338 1347 vg_vertical_gradient_level_ind(1) = 0 1339 1348 vg(0) = vg_surface 1340 1349 DO k = 1, nzt+1 1341 IF ( i < 11 ) THEN1342 IF ( vg_vertical_gradient_level(i) < zu(k) .AND. &1350 IF ( i < 11 ) THEN 1351 IF ( vg_vertical_gradient_level(i) < zu(k) .AND. & 1343 1352 vg_vertical_gradient_level(i) >= 0.0_wp ) THEN 1344 1353 gradient = vg_vertical_gradient(i) / 100.0_wp … … 1363 1372 vg(nzt+1) = vg_surface 1364 1373 DO k = nzt, nzb, -1 1365 IF ( i < 11 ) THEN1366 IF ( vg_vertical_gradient_level(i) > zu(k) .AND. &1374 IF ( i < 11 ) THEN 1375 IF ( vg_vertical_gradient_level(i) > zu(k) .AND. & 1367 1376 vg_vertical_gradient_level(i) <= 0.0_wp ) THEN 1368 1377 gradient = vg_vertical_gradient(i) / 100.0_wp … … 1415 1424 1416 1425 IF ( kk < 100 ) THEN 1417 DO WHILE ( uv_heights(kk+1) <= zu(k) )1426 DO WHILE ( uv_heights(kk+1) <= zu(k) ) 1418 1427 kk = kk + 1 1419 1428 IF ( kk == 100 ) EXIT … … 1421 1430 ENDIF 1422 1431 1423 IF ( kk < 100 .AND.uv_heights(kk+1) /= 9999999.9_wp ) THEN1432 IF ( kk < 100 .AND. uv_heights(kk+1) /= 9999999.9_wp ) THEN 1424 1433 u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) / & 1425 1434 ( uv_heights(kk+1) - uv_heights(kk) ) * & … … 1444 1453 ! 1445 1454 !-- Compute initial temperature profile using the given temperature gradients 1446 IF ( .NOT.neutral ) THEN1455 IF ( .NOT. neutral ) THEN 1447 1456 1448 1457 i = 1 1449 1458 gradient = 0.0_wp 1450 1459 1451 IF ( .NOT.ocean ) THEN1460 IF ( .NOT. ocean ) THEN 1452 1461 1453 1462 pt_vertical_gradient_level_ind(1) = 0 1454 1463 DO k = 1, nzt+1 1455 IF ( i < 11 ) THEN1456 IF ( pt_vertical_gradient_level(i) < zu(k) .AND. &1464 IF ( i < 11 ) THEN 1465 IF ( pt_vertical_gradient_level(i) < zu(k) .AND. & 1457 1466 pt_vertical_gradient_level(i) >= 0.0_wp ) THEN 1458 1467 gradient = pt_vertical_gradient(i) / 100.0_wp … … 1476 1485 pt_vertical_gradient_level_ind(1) = nzt+1 1477 1486 DO k = nzt, 0, -1 1478 IF ( i < 11 ) THEN1479 IF ( pt_vertical_gradient_level(i) > zu(k) .AND. &1487 IF ( i < 11 ) THEN 1488 IF ( pt_vertical_gradient_level(i) > zu(k) .AND. & 1480 1489 pt_vertical_gradient_level(i) <= 0.0_wp ) THEN 1481 1490 gradient = pt_vertical_gradient(i) / 100.0_wp … … 1533 1542 q_vertical_gradient_level_ind(1) = 0 1534 1543 DO k = 1, nzt+1 1535 IF ( i < 11 ) THEN1536 IF ( q_vertical_gradient_level(i) < zu(k) .AND. &1544 IF ( i < 11 ) THEN 1545 IF ( q_vertical_gradient_level(i) < zu(k) .AND. & 1537 1546 q_vertical_gradient_level(i) >= 0.0_wp ) THEN 1538 1547 gradient = q_vertical_gradient(i) / 100.0_wp … … 1579 1588 sa_vertical_gradient_level_ind(1) = nzt+1 1580 1589 DO k = nzt, 0, -1 1581 IF ( i < 11 ) THEN1582 IF ( sa_vertical_gradient_level(i) > zu(k) .AND. &1590 IF ( i < 11 ) THEN 1591 IF ( sa_vertical_gradient_level(i) > zu(k) .AND. & 1583 1592 sa_vertical_gradient_level(i) <= 0.0_wp ) THEN 1584 1593 gradient = sa_vertical_gradient(i) / 100.0_wp … … 1606 1615 ! 1607 1616 !-- Check if the control parameter use_subsidence_tendencies is used correctly 1608 IF ( use_subsidence_tendencies .AND. .NOT. large_scale_subsidence ) THEN1609 message_string = 'The usage of use_subsidence_tendencies ' // &1617 IF ( use_subsidence_tendencies .AND. .NOT. large_scale_subsidence ) THEN 1618 message_string = 'The usage of use_subsidence_tendencies ' // & 1610 1619 'requires large_scale_subsidence = .T..' 1611 1620 CALL message( 'check_parameters', 'PA0396', 1, 2, 0, 6, 0 ) 1612 1621 ELSEIF ( use_subsidence_tendencies .AND. .NOT. large_scale_forcing ) THEN 1613 message_string = 'The usage of use_subsidence_tendencies ' // &1622 message_string = 'The usage of use_subsidence_tendencies ' // & 1614 1623 'requires large_scale_forcing = .T..' 1615 1624 CALL message( 'check_parameters', 'PA0397', 1, 2, 0, 6, 0 ) … … 1619 1628 !-- Initialize large scale subsidence if required 1620 1629 If ( large_scale_subsidence ) THEN 1621 IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp .AND.&1622 .NOT. large_scale_forcing ) THEN1630 IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp .AND. & 1631 .NOT. large_scale_forcing ) THEN 1623 1632 CALL init_w_subsidence 1624 1633 ENDIF … … 1627 1636 !-- are read in from file LSF_DATA 1628 1637 1629 IF ( subs_vertical_gradient_level(1) == -9999999.9_wp .AND.&1630 .NOT.large_scale_forcing ) THEN1631 message_string = 'There is no default large scale vertical ' // &1632 'velocity profile set. Specify the subsidence ' // &1633 'velocity profile via subs_vertical_gradient and ' //&1634 ' subs_vertical_gradient_level.'1638 IF ( subs_vertical_gradient_level(1) == -9999999.9_wp .AND. & 1639 .NOT. large_scale_forcing ) THEN 1640 message_string = 'There is no default large scale vertical ' // & 1641 'velocity profile set. Specify the subsidence ' // & 1642 'velocity profile via subs_vertical_gradient ' // & 1643 'and subs_vertical_gradient_level.' 1635 1644 CALL message( 'check_parameters', 'PA0380', 1, 2, 0, 6, 0 ) 1636 1645 ENDIF 1637 1646 ELSE 1638 1647 IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp ) THEN 1639 message_string = 'Enable usage of large scale subsidence by ' // &1648 message_string = 'Enable usage of large scale subsidence by ' // & 1640 1649 'setting large_scale_subsidence = .T..' 1641 1650 CALL message( 'check_parameters', 'PA0381', 1, 2, 0, 6, 0 ) … … 1659 1668 vpt_reference = pt_reference * ( 1.0_wp + 0.61_wp * q_surface ) 1660 1669 ELSE 1661 message_string = 'illegal value for reference_state: "' // &1670 message_string = 'illegal value for reference_state: "' // & 1662 1671 TRIM( reference_state ) // '"' 1663 1672 CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 ) … … 1686 1695 IF ( alpha_surface /= 0.0_wp ) THEN 1687 1696 IF ( ABS( alpha_surface ) > 90.0_wp ) THEN 1688 WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface, &1697 WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface, & 1689 1698 ' ) must be < 90.0' 1690 1699 CALL message( 'check_parameters', 'PA0043', 1, 2, 0, 6, 0 ) … … 1716 1725 ENDIF 1717 1726 ELSE 1718 WRITE( message_string, * ) 'cfl_factor = ', cfl_factor, &1727 WRITE( message_string, * ) 'cfl_factor = ', cfl_factor, & 1719 1728 ' out of range & 0.0 < cfl_factor <= 1.0 is required' 1720 1729 CALL message( 'check_parameters', 'PA0045', 1, 2, 0, 6, 0 ) … … 1740 1749 !-- Set wind speed in the Galilei-transformed system 1741 1750 IF ( galilei_transformation ) THEN 1742 IF ( use_ug_for_galilei_tr .AND. &1743 ug_vertical_gradient_level(1) == 0.0_wp .AND. &1744 ug_vertical_gradient(1) == 0.0_wp .AND.&1745 vg_vertical_gradient_level(1) == 0.0_wp .AND. &1751 IF ( use_ug_for_galilei_tr .AND. & 1752 ug_vertical_gradient_level(1) == 0.0_wp .AND. & 1753 ug_vertical_gradient(1) == 0.0_wp .AND. & 1754 vg_vertical_gradient_level(1) == 0.0_wp .AND. & 1746 1755 vg_vertical_gradient(1) == 0.0_wp ) THEN 1747 1756 u_gtrans = ug_surface * 0.6_wp 1748 1757 v_gtrans = vg_surface * 0.6_wp 1749 ELSEIF ( use_ug_for_galilei_tr .AND. &1750 ( ug_vertical_gradient_level(1) /= 0.0_wp .OR. &1758 ELSEIF ( use_ug_for_galilei_tr .AND. & 1759 ( ug_vertical_gradient_level(1) /= 0.0_wp .OR. & 1751 1760 ug_vertical_gradient(1) /= 0.0_wp ) ) THEN 1752 message_string = 'baroclinity (ug) not allowed simultaneously' // &1761 message_string = 'baroclinity (ug) not allowed simultaneously' // & 1753 1762 ' with galilei transformation' 1754 1763 CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 ) 1755 ELSEIF ( use_ug_for_galilei_tr .AND. &1756 ( vg_vertical_gradient_level(1) /= 0.0_wp .OR. &1764 ELSEIF ( use_ug_for_galilei_tr .AND. & 1765 ( vg_vertical_gradient_level(1) /= 0.0_wp .OR. & 1757 1766 vg_vertical_gradient(1) /= 0.0_wp ) ) THEN 1758 message_string = 'baroclinity (vg) not allowed simultaneously' // &1767 message_string = 'baroclinity (vg) not allowed simultaneously' // & 1759 1768 ' with galilei transformation' 1760 1769 CALL message( 'check_parameters', 'PA0047', 1, 2, 0, 6, 0 ) 1761 1770 ELSE 1762 message_string = 'variable translation speed used for galilei-' // &1763 'transformation, which may cause & instabilities in stably ' // &1771 message_string = 'variable translation speed used for galilei-' // & 1772 'transformation, which may cause & instabilities in stably ' // & 1764 1773 'stratified regions' 1765 1774 CALL message( 'check_parameters', 'PA0048', 0, 1, 0, 6, 0 ) … … 1782 1791 IF ( bc_lr /= 'cyclic' .OR. bc_ns /= 'cyclic' ) THEN 1783 1792 IF ( psolver(1:9) /= 'multigrid' ) THEN 1784 message_string = 'non-cyclic lateral boundaries do not allow ' // &1793 message_string = 'non-cyclic lateral boundaries do not allow ' // & 1785 1794 'psolver = "' // TRIM( psolver ) // '"' 1786 1795 CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 ) 1787 1796 ENDIF 1788 IF ( momentum_advec /= 'pw-scheme' .AND. &1789 ( momentum_advec /= 'ws-scheme' .AND. &1790 momentum_advec /= 'ws-scheme-mono' ) &1797 IF ( momentum_advec /= 'pw-scheme' .AND. & 1798 ( momentum_advec /= 'ws-scheme' .AND. & 1799 momentum_advec /= 'ws-scheme-mono' ) & 1791 1800 ) THEN 1792 1801 1793 message_string = 'non-cyclic lateral boundaries do not allow ' // &1802 message_string = 'non-cyclic lateral boundaries do not allow ' // & 1794 1803 'momentum_advec = "' // TRIM( momentum_advec ) // '"' 1795 1804 CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 ) 1796 1805 ENDIF 1797 IF ( scalar_advec /= 'pw-scheme' .AND.&1798 ( scalar_advec /= 'ws-scheme' .AND.&1806 IF ( scalar_advec /= 'pw-scheme' .AND. & 1807 ( scalar_advec /= 'ws-scheme' .AND. & 1799 1808 scalar_advec /= 'ws-scheme-mono' ) & 1800 1809 ) THEN 1801 message_string = 'non-cyclic lateral boundaries do not allow ' // &1810 message_string = 'non-cyclic lateral boundaries do not allow ' // & 1802 1811 'scalar_advec = "' // TRIM( scalar_advec ) // '"' 1803 1812 CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 ) 1804 1813 ENDIF 1805 1814 IF ( galilei_transformation ) THEN 1806 message_string = 'non-cyclic lateral boundaries do not allow ' // &1815 message_string = 'non-cyclic lateral boundaries do not allow ' // & 1807 1816 'galilei_transformation = .T.' 1808 1817 CALL message( 'check_parameters', 'PA0054', 1, 2, 0, 6, 0 ) … … 1819 1828 bc_e_b = 'neumann' 1820 1829 ibc_e_b = 1 1821 message_string = 'boundary condition bc_e_b changed to "' // &1830 message_string = 'boundary condition bc_e_b changed to "' // & 1822 1831 TRIM( bc_e_b ) // '"' 1823 1832 CALL message( 'check_parameters', 'PA0057', 0, 1, 0, 6, 0 ) 1824 1833 ENDIF 1825 1834 ELSE 1826 message_string = 'unknown boundary condition: bc_e_b = "' // &1835 message_string = 'unknown boundary condition: bc_e_b = "' // & 1827 1836 TRIM( bc_e_b ) // '"' 1828 1837 CALL message( 'check_parameters', 'PA0058', 1, 2, 0, 6, 0 ) … … 1836 1845 ibc_p_b = 1 1837 1846 ELSE 1838 message_string = 'unknown boundary condition: bc_p_b = "' // &1847 message_string = 'unknown boundary condition: bc_p_b = "' // & 1839 1848 TRIM( bc_p_b ) // '"' 1840 1849 CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 ) … … 1847 1856 ibc_p_t = 1 1848 1857 ELSE 1849 message_string = 'unknown boundary condition: bc_p_t = "' // &1858 message_string = 'unknown boundary condition: bc_p_t = "' // & 1850 1859 TRIM( bc_p_t ) // '"' 1851 1860 CALL message( 'check_parameters', 'PA0061', 1, 2, 0, 6, 0 ) … … 1862 1871 ibc_pt_b = 1 1863 1872 ELSE 1864 message_string = 'unknown boundary condition: bc_pt_b = "' // &1873 message_string = 'unknown boundary condition: bc_pt_b = "' // & 1865 1874 TRIM( bc_pt_b ) // '"' 1866 1875 CALL message( 'check_parameters', 'PA0062', 1, 2, 0, 6, 0 ) … … 1877 1886 ibc_pt_t = 3 1878 1887 ELSE 1879 message_string = 'unknown boundary condition: bc_pt_t = "' // &1888 message_string = 'unknown boundary condition: bc_pt_t = "' // & 1880 1889 TRIM( bc_pt_t ) // '"' 1881 1890 CALL message( 'check_parameters', 'PA0063', 1, 2, 0, 6, 0 ) … … 1889 1898 ELSEIF ( ibc_pt_b == 1 ) THEN 1890 1899 constant_heatflux = .TRUE. 1891 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND.&1892 .NOT. land_surface ) THEN1900 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. & 1901 .NOT. land_surface ) THEN 1893 1902 surface_heatflux = shf_surf(1) 1894 1903 ELSE … … 1899 1908 ELSE 1900 1909 constant_heatflux = .TRUE. 1901 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND.&1902 large_scale_forcing .AND. .NOT.land_surface ) THEN1910 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. & 1911 large_scale_forcing .AND. .NOT. land_surface ) THEN 1903 1912 surface_heatflux = shf_surf(1) 1904 1913 ENDIF … … 1909 1918 IF ( neutral ) THEN 1910 1919 1911 IF ( surface_heatflux /= 0.0_wp .AND. surface_heatflux /= 9999999.9_wp ) & 1920 IF ( surface_heatflux /= 0.0_wp .AND. & 1921 surface_heatflux /= 9999999.9_wp ) THEN 1922 message_string = 'heatflux must not be set for pure neutral flow' 1923 CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 ) 1924 ENDIF 1925 1926 IF ( top_heatflux /= 0.0_wp .AND. top_heatflux /= 9999999.9_wp ) & 1912 1927 THEN 1913 1928 message_string = 'heatflux must not be set for pure neutral flow' … … 1915 1930 ENDIF 1916 1931 1917 IF ( top_heatflux /= 0.0_wp .AND. top_heatflux /= 9999999.9_wp ) & 1918 THEN 1919 message_string = 'heatflux must not be set for pure neutral flow' 1920 CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 ) 1921 ENDIF 1922 1923 ENDIF 1924 1925 IF ( top_momentumflux_u /= 9999999.9_wp .AND. & 1932 ENDIF 1933 1934 IF ( top_momentumflux_u /= 9999999.9_wp .AND. & 1926 1935 top_momentumflux_v /= 9999999.9_wp ) THEN 1927 1936 constant_top_momentumflux = .TRUE. 1928 ELSEIF ( .NOT. ( top_momentumflux_u == 9999999.9_wp .AND. &1937 ELSEIF ( .NOT. ( top_momentumflux_u == 9999999.9_wp .AND. & 1929 1938 top_momentumflux_v == 9999999.9_wp ) ) THEN 1930 message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' // &1939 message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' // & 1931 1940 'must be set' 1932 1941 CALL message( 'check_parameters', 'PA0064', 1, 2, 0, 6, 0 ) … … 1937 1946 !-- temperature. In this case specification of a constant heat flux is 1938 1947 !-- forbidden. 1939 IF ( ibc_pt_b == 0 .AND. constant_heatflux .AND.&1948 IF ( ibc_pt_b == 0 .AND. constant_heatflux .AND. & 1940 1949 surface_heatflux /= 0.0_wp ) THEN 1941 1950 message_string = 'boundary_condition: bc_pt_b = "' // TRIM( bc_pt_b ) //& … … 1944 1953 ENDIF 1945 1954 IF ( constant_heatflux .AND. pt_surface_initial_change /= 0.0_wp ) THEN 1946 WRITE ( message_string, * ) 'constant_heatflux = .TRUE. is not allo', &1947 'wed with pt_surface_initial_change (/=0) = ', &1955 WRITE ( message_string, * ) 'constant_heatflux = .TRUE. is not allo', & 1956 'wed with pt_surface_initial_change (/=0) = ', & 1948 1957 pt_surface_initial_change 1949 1958 CALL message( 'check_parameters', 'PA0066', 1, 2, 0, 6, 0 ) … … 1954 1963 !-- temperature. In this case specification of a constant heat flux is 1955 1964 !-- forbidden. 1956 IF ( ibc_pt_t == 0 .AND. constant_top_heatflux .AND.&1965 IF ( ibc_pt_t == 0 .AND. constant_top_heatflux .AND. & 1957 1966 top_heatflux /= 0.0_wp ) THEN 1958 1967 message_string = 'boundary_condition: bc_pt_t = "' // TRIM( bc_pt_t ) //& … … 1969 1978 ibc_sa_t = 1 1970 1979 ELSE 1971 message_string = 'unknown boundary condition: bc_sa_t = "' // &1980 message_string = 'unknown boundary condition: bc_sa_t = "' // & 1972 1981 TRIM( bc_sa_t ) // '"' 1973 1982 CALL message( 'check_parameters', 'PA0068', 1, 2, 0, 6, 0 ) … … 1975 1984 1976 1985 IF ( top_salinityflux == 9999999.9_wp ) constant_top_salinityflux = .FALSE. 1977 IF ( ibc_sa_t == 1 .AND. 1978 message_string = 'boundary condition: bc_sa_t = "' // &1979 TRIM( bc_sa_t ) // '" requires to set ' // &1986 IF ( ibc_sa_t == 1 .AND. top_salinityflux == 9999999.9_wp ) THEN 1987 message_string = 'boundary condition: bc_sa_t = "' // & 1988 TRIM( bc_sa_t ) // '" requires to set ' // & 1980 1989 'top_salinityflux' 1981 1990 CALL message( 'check_parameters', 'PA0069', 1, 2, 0, 6, 0 ) … … 1986 1995 !-- salinity. In this case specification of a constant salinity flux is 1987 1996 !-- forbidden. 1988 IF ( ibc_sa_t == 0 .AND. constant_top_salinityflux .AND.&1997 IF ( ibc_sa_t == 0 .AND. constant_top_salinityflux .AND. & 1989 1998 top_salinityflux /= 0.0_wp ) THEN 1990 message_string = 'boundary condition: bc_sa_t = "' // &1991 TRIM( bc_sa_t ) // '" is not allowed with ' // &1999 message_string = 'boundary condition: bc_sa_t = "' // & 2000 TRIM( bc_sa_t ) // '" is not allowed with ' // & 1992 2001 'constant_top_salinityflux = .TRUE.' 1993 2002 CALL message( 'check_parameters', 'PA0070', 1, 2, 0, 6, 0 ) … … 2010 2019 ibc_q_b = 1 2011 2020 ELSE 2012 message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // &2021 message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // & 2013 2022 '_b ="' // TRIM( bc_q_b ) // '"' 2014 2023 CALL message( 'check_parameters', 'PA0071', 1, 2, 0, 6, 0 ) … … 2021 2030 ibc_q_t = 3 2022 2031 ELSE 2023 message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // &2032 message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // & 2024 2033 '_t ="' // TRIM( bc_q_t ) // '"' 2025 2034 CALL message( 'check_parameters', 'PA0072', 1, 2, 0, 6, 0 ) … … 2040 2049 ELSE 2041 2050 constant_waterflux = .TRUE. 2042 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND.&2043 large_scale_forcing )THEN2051 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. & 2052 large_scale_forcing ) THEN 2044 2053 surface_waterflux = qsws_surf(1) 2045 2054 ENDIF … … 2058 2067 IF ( constant_waterflux .AND. q_surface_initial_change /= 0.0_wp ) THEN 2059 2068 WRITE( message_string, * ) 'a prescribed surface flux is not allo', & 2060 'wed with ', sq, '_surface_initial_change (/=0) = ', &2069 'wed with ', sq, '_surface_initial_change (/=0) = ', & 2061 2070 q_surface_initial_change 2062 2071 CALL message( 'check_parameters', 'PA0074', 1, 2, 0, 6, 0 ) … … 2071 2080 ibc_uv_b = 1 2072 2081 IF ( constant_flux_layer ) THEN 2073 message_string = 'boundary condition: bc_uv_b = "' // &2082 message_string = 'boundary condition: bc_uv_b = "' // & 2074 2083 TRIM( bc_uv_b ) // '" is not allowed with constant_flux_layer' & 2075 2084 // ' = .TRUE.' … … 2077 2086 ENDIF 2078 2087 ELSE 2079 message_string = 'unknown boundary condition: bc_uv_b = "' // &2088 message_string = 'unknown boundary condition: bc_uv_b = "' // & 2080 2089 TRIM( bc_uv_b ) // '"' 2081 2090 CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 ) … … 2106 2115 ibc_uv_t = 3 2107 2116 ELSE 2108 message_string = 'unknown boundary condition: bc_uv_t = "' // &2117 message_string = 'unknown boundary condition: bc_uv_t = "' // & 2109 2118 TRIM( bc_uv_t ) // '"' 2110 2119 CALL message( 'check_parameters', 'PA0077', 1, 2, 0, 6, 0 ) … … 2117 2126 rayleigh_damping_factor = 0.0_wp 2118 2127 ELSE 2119 IF ( rayleigh_damping_factor < 0.0_wp .OR. rayleigh_damping_factor > 1.0_wp )&2120 THEN2121 WRITE( message_string, * ) 'rayleigh_damping_factor = ', &2128 IF ( rayleigh_damping_factor < 0.0_wp .OR. & 2129 rayleigh_damping_factor > 1.0_wp ) THEN 2130 WRITE( message_string, * ) 'rayleigh_damping_factor = ', & 2122 2131 rayleigh_damping_factor, ' out of range [0.0,1.0]' 2123 2132 CALL message( 'check_parameters', 'PA0078', 1, 2, 0, 6, 0 ) … … 2126 2135 2127 2136 IF ( rayleigh_damping_height == -1.0_wp ) THEN 2128 IF ( .NOT.ocean ) THEN2137 IF ( .NOT. ocean ) THEN 2129 2138 rayleigh_damping_height = 0.66666666666_wp * zu(nzt) 2130 2139 ELSE … … 2132 2141 ENDIF 2133 2142 ELSE 2134 IF ( .NOT.ocean ) THEN2135 IF ( rayleigh_damping_height < 0.0_wp .OR. &2143 IF ( .NOT. ocean ) THEN 2144 IF ( rayleigh_damping_height < 0.0_wp .OR. & 2136 2145 rayleigh_damping_height > zu(nzt) ) THEN 2137 WRITE( message_string, * ) 'rayleigh_damping_height = ', &2146 WRITE( message_string, * ) 'rayleigh_damping_height = ', & 2138 2147 rayleigh_damping_height, ' out of range [0.0,', zu(nzt), ']' 2139 2148 CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 ) 2140 2149 ENDIF 2141 2150 ELSE 2142 IF ( rayleigh_damping_height > 0.0_wp .OR. &2151 IF ( rayleigh_damping_height > 0.0_wp .OR. & 2143 2152 rayleigh_damping_height < zu(nzb) ) THEN 2144 WRITE( message_string, * ) 'rayleigh_damping_height = ', &2153 WRITE( message_string, * ) 'rayleigh_damping_height = ', & 2145 2154 rayleigh_damping_height, ' out of range [0.0,', zu(nzb), ']' 2146 2155 CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 ) … … 2154 2163 !-- be opened (cf. check_open) 2155 2164 IF ( statistic_regions > 9 .OR. statistic_regions < 0 ) THEN 2156 WRITE ( message_string, * ) 'number of statistic_regions = ', &2165 WRITE ( message_string, * ) 'number of statistic_regions = ', & 2157 2166 statistic_regions+1, ' but only 10 regions are allowed' 2158 2167 CALL message( 'check_parameters', 'PA0082', 1, 2, 0, 6, 0 ) 2159 2168 ENDIF 2160 IF ( normalizing_region > statistic_regions .OR. &2169 IF ( normalizing_region > statistic_regions .OR. & 2161 2170 normalizing_region < 0) THEN 2162 WRITE ( message_string, * ) 'normalizing_region = ', &2171 WRITE ( message_string, * ) 'normalizing_region = ', & 2163 2172 normalizing_region, ' must be >= 0 and <= ',statistic_regions, & 2164 2173 ' (value of statistic_regions)' … … 2184 2193 ! 2185 2194 !-- Set the default skip time intervals for data output, if necessary 2186 IF ( skip_time_dopr == 9999999.9_wp ) &2195 IF ( skip_time_dopr == 9999999.9_wp ) & 2187 2196 skip_time_dopr = skip_time_data_output 2188 IF ( skip_time_dosp == 9999999.9_wp ) &2197 IF ( skip_time_dosp == 9999999.9_wp ) & 2189 2198 skip_time_dosp = skip_time_data_output 2190 IF ( skip_time_do2d_xy == 9999999.9_wp ) &2199 IF ( skip_time_do2d_xy == 9999999.9_wp ) & 2191 2200 skip_time_do2d_xy = skip_time_data_output 2192 IF ( skip_time_do2d_xz == 9999999.9_wp ) &2201 IF ( skip_time_do2d_xz == 9999999.9_wp ) & 2193 2202 skip_time_do2d_xz = skip_time_data_output 2194 IF ( skip_time_do2d_yz == 9999999.9_wp ) &2203 IF ( skip_time_do2d_yz == 9999999.9_wp ) & 2195 2204 skip_time_do2d_yz = skip_time_data_output 2196 IF ( skip_time_do3d == 9999999.9_wp ) &2205 IF ( skip_time_do3d == 9999999.9_wp ) & 2197 2206 skip_time_do3d = skip_time_data_output 2198 IF ( skip_time_data_output_av == 9999999.9_wp ) &2207 IF ( skip_time_data_output_av == 9999999.9_wp ) & 2199 2208 skip_time_data_output_av = skip_time_data_output 2200 2209 DO mid = 1, max_masks 2201 IF ( skip_time_domask(mid) == 9999999.9_wp ) &2210 IF ( skip_time_domask(mid) == 9999999.9_wp ) & 2202 2211 skip_time_domask(mid) = skip_time_data_output 2203 2212 ENDDO … … 2207 2216 !-- spectra) 2208 2217 IF ( averaging_interval > dt_data_output_av ) THEN 2209 WRITE( message_string, * ) 'averaging_interval = ', &2218 WRITE( message_string, * ) 'averaging_interval = ', & 2210 2219 averaging_interval, ' must be <= dt_data_output = ', dt_data_output 2211 2220 CALL message( 'check_parameters', 'PA0085', 1, 2, 0, 6, 0 ) … … 2217 2226 2218 2227 IF ( averaging_interval_pr > dt_dopr ) THEN 2219 WRITE( message_string, * ) 'averaging_interval_pr = ', &2228 WRITE( message_string, * ) 'averaging_interval_pr = ', & 2220 2229 averaging_interval_pr, ' must be <= dt_dopr = ', dt_dopr 2221 2230 CALL message( 'check_parameters', 'PA0086', 1, 2, 0, 6, 0 ) … … 2227 2236 2228 2237 IF ( averaging_interval_sp > dt_dosp ) THEN 2229 WRITE( message_string, * ) 'averaging_interval_sp = ', &2238 WRITE( message_string, * ) 'averaging_interval_sp = ', & 2230 2239 averaging_interval_sp, ' must be <= dt_dosp = ', dt_dosp 2231 2240 CALL message( 'check_parameters', 'PA0087', 1, 2, 0, 6, 0 ) … … 2252 2261 !-- Check the sample rate for averaging (first for 3d-data, then for profiles) 2253 2262 IF ( dt_averaging_input > averaging_interval ) THEN 2254 WRITE( message_string, * ) 'dt_averaging_input = ', &2255 dt_averaging_input, ' must be <= averaging_interval = ', &2263 WRITE( message_string, * ) 'dt_averaging_input = ', & 2264 dt_averaging_input, ' must be <= averaging_interval = ', & 2256 2265 averaging_interval 2257 2266 CALL message( 'check_parameters', 'PA0088', 1, 2, 0, 6, 0 ) … … 2259 2268 2260 2269 IF ( dt_averaging_input_pr > averaging_interval_pr ) THEN 2261 WRITE( message_string, * ) 'dt_averaging_input_pr = ', &2270 WRITE( message_string, * ) 'dt_averaging_input_pr = ', & 2262 2271 dt_averaging_input_pr, ' must be <= averaging_interval_pr = ', & 2263 2272 averaging_interval_pr … … 2272 2281 ELSE 2273 2282 IF ( precipitation_amount_interval > dt_do2d_xy ) THEN 2274 WRITE( message_string, * ) 'precipitation_amount_interval = ', &2275 precipitation_amount_interval, ' must not be larger than ', &2283 WRITE( message_string, * ) 'precipitation_amount_interval = ', & 2284 precipitation_amount_interval, ' must not be larger than ', & 2276 2285 'dt_do2d_xy = ', dt_do2d_xy 2277 2286 CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 ) … … 2532 2541 CASE ( 's', '#s' ) 2533 2542 IF ( .NOT. passive_scalar ) THEN 2534 message_string = 'data_output_pr = ' // &2543 message_string = 'data_output_pr = ' // & 2535 2544 TRIM( data_output_pr(i) ) // ' is not imp' // & 2536 2545 'lemented for passive_scalar = .FALSE.' … … 2573 2582 CASE ( 'lpt', '#lpt' ) 2574 2583 IF ( .NOT. cloud_physics ) THEN 2575 message_string = 'data_output_pr = ' // &2584 message_string = 'data_output_pr = ' // & 2576 2585 TRIM( data_output_pr(i) ) // ' is not imp' // & 2577 2586 'lemented for cloud_physics = .FALSE.' … … 2616 2625 2617 2626 CASE ( 'w"q"' ) 2618 IF ( .NOT.humidity ) THEN2619 message_string = 'data_output_pr = ' // &2627 IF ( .NOT. humidity ) THEN 2628 message_string = 'data_output_pr = ' // & 2620 2629 TRIM( data_output_pr(i) ) // ' is not imp' // & 2621 2630 'lemented for humidity = .FALSE.' … … 2628 2637 2629 2638 CASE ( 'w*q*' ) 2630 IF ( .NOT.humidity ) THEN2631 message_string = 'data_output_pr = ' // &2639 IF ( .NOT. humidity ) THEN 2640 message_string = 'data_output_pr = ' // & 2632 2641 TRIM( data_output_pr(i) ) // ' is not imp' // & 2633 2642 'lemented for humidity = .FALSE.' … … 2640 2649 2641 2650 CASE ( 'wq' ) 2642 IF ( .NOT.humidity ) THEN2643 message_string = 'data_output_pr = ' // &2651 IF ( .NOT. humidity ) THEN 2652 message_string = 'data_output_pr = ' // & 2644 2653 TRIM( data_output_pr(i) ) // ' is not imp' // & 2645 2654 'lemented for humidity = .FALSE.' … … 2652 2661 2653 2662 CASE ( 'w"s"' ) 2654 IF ( .NOT. passive_scalar )THEN2655 message_string = 'data_output_pr = ' // &2663 IF ( .NOT. passive_scalar ) THEN 2664 message_string = 'data_output_pr = ' // & 2656 2665 TRIM( data_output_pr(i) ) // ' is not imp' // & 2657 2666 'lemented for passive_scalar = .FALSE.' … … 2664 2673 2665 2674 CASE ( 'w*s*' ) 2666 IF ( .NOT. passive_scalar )THEN2667 message_string = 'data_output_pr = ' // &2675 IF ( .NOT. passive_scalar ) THEN 2676 message_string = 'data_output_pr = ' // & 2668 2677 TRIM( data_output_pr(i) ) // ' is not imp' // & 2669 2678 'lemented for passive_scalar = .FALSE.' … … 2676 2685 2677 2686 CASE ( 'ws' ) 2678 IF ( .NOT. passive_scalar )THEN2679 message_string = 'data_output_pr = ' // &2687 IF ( .NOT. passive_scalar ) THEN 2688 message_string = 'data_output_pr = ' // & 2680 2689 TRIM( data_output_pr(i) ) // ' is not imp' // & 2681 2690 'lemented for passive_scalar = .FALSE.' … … 2688 2697 2689 2698 CASE ( 'w"qv"' ) 2690 IF ( humidity .AND. .NOT. cloud_physics ) & 2691 THEN 2699 IF ( humidity .AND. .NOT. cloud_physics ) THEN 2692 2700 dopr_index(i) = 48 2693 2701 dopr_unit(i) = 'kg/kg m/s' 2694 2702 hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 ) 2695 ELSEIF ( humidity .AND. cloud_physics )THEN2703 ELSEIF ( humidity .AND. cloud_physics ) THEN 2696 2704 dopr_index(i) = 51 2697 2705 dopr_unit(i) = 'kg/kg m/s' 2698 2706 hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 ) 2699 2707 ELSE 2700 message_string = 'data_output_pr = ' // &2708 message_string = 'data_output_pr = ' // & 2701 2709 TRIM( data_output_pr(i) ) // ' is not imp' // & 2702 2710 'lemented for cloud_physics = .FALSE. an&' // & … … 2706 2714 2707 2715 CASE ( 'w*qv*' ) 2708 IF ( humidity .AND. .NOT. cloud_physics ) &2716 IF ( humidity .AND. .NOT. cloud_physics ) & 2709 2717 THEN 2710 2718 dopr_index(i) = 49 … … 2716 2724 hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 ) 2717 2725 ELSE 2718 message_string = 'data_output_pr = ' // &2726 message_string = 'data_output_pr = ' // & 2719 2727 TRIM( data_output_pr(i) ) // ' is not imp' // & 2720 2728 'lemented for cloud_physics = .FALSE. an&' // & … … 2724 2732 2725 2733 CASE ( 'wqv' ) 2726 IF ( humidity .AND. .NOT. cloud_physics ) & 2727 THEN 2734 IF ( humidity .AND. .NOT. cloud_physics ) THEN 2728 2735 dopr_index(i) = 50 2729 2736 dopr_unit(i) = 'kg/kg m/s' 2730 2737 hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 ) 2731 ELSEIF ( humidity .AND. cloud_physics )THEN2738 ELSEIF ( humidity .AND. cloud_physics ) THEN 2732 2739 dopr_index(i) = 53 2733 2740 dopr_unit(i) = 'kg/kg m/s' … … 2742 2749 2743 2750 CASE ( 'ql' ) 2744 IF ( .NOT. cloud_physics .AND. .NOT.cloud_droplets ) THEN2751 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 2745 2752 message_string = 'data_output_pr = ' // & 2746 2753 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 2801 2808 2802 2809 CASE ( 'rho' ) 2803 IF ( .NOT.ocean ) THEN2810 IF ( .NOT. ocean ) THEN 2804 2811 message_string = 'data_output_pr = ' // & 2805 2812 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 2819 2826 2820 2827 CASE ( 'w"sa"' ) 2821 IF ( .NOT.ocean ) THEN2828 IF ( .NOT. ocean ) THEN 2822 2829 message_string = 'data_output_pr = ' // & 2823 2830 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 2831 2838 2832 2839 CASE ( 'w*sa*' ) 2833 IF ( .NOT. ocean) THEN2840 IF ( .NOT. ocean ) THEN 2834 2841 message_string = 'data_output_pr = ' // & 2835 2842 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 2843 2850 2844 2851 CASE ( 'wsa' ) 2845 IF ( .NOT.ocean ) THEN2852 IF ( .NOT. ocean ) THEN 2846 2853 message_string = 'data_output_pr = ' // & 2847 2854 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 2865 2872 2866 2873 CASE ( 'q*2' ) 2867 IF ( .NOT.humidity ) THEN2874 IF ( .NOT. humidity ) THEN 2868 2875 message_string = 'data_output_pr = ' // & 2869 2876 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 2877 2884 2878 2885 CASE ( 'prho' ) 2879 IF ( .NOT.ocean ) THEN2886 IF ( .NOT. ocean ) THEN 2880 2887 message_string = 'data_output_pr = ' // & 2881 2888 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 2894 2901 2895 2902 CASE ( 'nr' ) 2896 IF ( .NOT.cloud_physics ) THEN2903 IF ( .NOT. cloud_physics ) THEN 2897 2904 message_string = 'data_output_pr = ' // & 2898 2905 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 2904 2911 'lemented for cloud_scheme /= seifert_beheng' 2905 2912 CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 ) 2906 ELSEIF ( .NOT.precipitation ) THEN2913 ELSEIF ( .NOT. precipitation ) THEN 2907 2914 message_string = 'data_output_pr = ' // & 2908 2915 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 2916 2923 2917 2924 CASE ( 'qr' ) 2918 IF ( .NOT.cloud_physics ) THEN2925 IF ( .NOT. cloud_physics ) THEN 2919 2926 message_string = 'data_output_pr = ' // & 2920 2927 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 2926 2933 'lemented for cloud_scheme /= seifert_beheng' 2927 2934 CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 ) 2928 ELSEIF ( .NOT.precipitation ) THEN2935 ELSEIF ( .NOT. precipitation ) THEN 2929 2936 message_string = 'data_output_pr = ' // & 2930 2937 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 2938 2945 2939 2946 CASE ( 'qc' ) 2940 IF ( .NOT.cloud_physics ) THEN2947 IF ( .NOT. cloud_physics ) THEN 2941 2948 message_string = 'data_output_pr = ' // & 2942 2949 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 2955 2962 2956 2963 CASE ( 'prr' ) 2957 IF ( .NOT.cloud_physics ) THEN2964 IF ( .NOT. cloud_physics ) THEN 2958 2965 message_string = 'data_output_pr = ' // & 2959 2966 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 2965 2972 'lemented for cloud_scheme /= seifert_beheng' 2966 2973 CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 ) 2967 ELSEIF ( .NOT.precipitation ) THEN2974 ELSEIF ( .NOT. precipitation ) THEN 2968 2975 message_string = 'data_output_pr = ' // & 2969 2976 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 2988 2995 2989 2996 CASE ( 'w_subs' ) 2990 IF ( .NOT.large_scale_subsidence ) THEN2997 IF ( .NOT. large_scale_subsidence ) THEN 2991 2998 message_string = 'data_output_pr = ' // & 2992 2999 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 3000 3007 3001 3008 CASE ( 'td_lsa_lpt' ) 3002 IF ( .NOT.large_scale_forcing ) THEN3009 IF ( .NOT. large_scale_forcing ) THEN 3003 3010 message_string = 'data_output_pr = ' // & 3004 3011 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 3012 3019 3013 3020 CASE ( 'td_lsa_q' ) 3014 IF ( .NOT.large_scale_forcing ) THEN3021 IF ( .NOT. large_scale_forcing ) THEN 3015 3022 message_string = 'data_output_pr = ' // & 3016 3023 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 3024 3031 3025 3032 CASE ( 'td_sub_lpt' ) 3026 IF ( .NOT.large_scale_forcing ) THEN3033 IF ( .NOT. large_scale_forcing ) THEN 3027 3034 message_string = 'data_output_pr = ' // & 3028 3035 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 3036 3043 3037 3044 CASE ( 'td_sub_q' ) 3038 IF ( .NOT.large_scale_forcing ) THEN3045 IF ( .NOT. large_scale_forcing ) THEN 3039 3046 message_string = 'data_output_pr = ' // & 3040 3047 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 3048 3055 3049 3056 CASE ( 'td_nud_lpt' ) 3050 IF ( .NOT.nudging ) THEN3057 IF ( .NOT. nudging ) THEN 3051 3058 message_string = 'data_output_pr = ' // & 3052 3059 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 3060 3067 3061 3068 CASE ( 'td_nud_q' ) 3062 IF ( .NOT.nudging ) THEN3069 IF ( .NOT. nudging ) THEN 3063 3070 message_string = 'data_output_pr = ' // & 3064 3071 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 3072 3079 3073 3080 CASE ( 'td_nud_u' ) 3074 IF ( .NOT.nudging ) THEN3081 IF ( .NOT. nudging ) THEN 3075 3082 message_string = 'data_output_pr = ' // & 3076 3083 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 3084 3091 3085 3092 CASE ( 'td_nud_v' ) 3086 IF ( .NOT.nudging ) THEN3093 IF ( .NOT. nudging ) THEN 3087 3094 message_string = 'data_output_pr = ' // & 3088 3095 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 3096 3103 3097 3104 CASE ( 't_soil', '#t_soil' ) 3098 IF ( .NOT.land_surface ) THEN3105 IF ( .NOT. land_surface ) THEN 3099 3106 message_string = 'data_output_pr = ' // & 3100 3107 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 3113 3120 3114 3121 CASE ( 'm_soil', '#m_soil' ) 3115 IF ( .NOT.land_surface ) THEN3122 IF ( .NOT. land_surface ) THEN 3116 3123 message_string = 'data_output_pr = ' // & 3117 3124 TRIM( data_output_pr(i) ) // ' is not imp' // & … … 3130 3137 3131 3138 CASE ( 'rad_net' ) 3132 IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' ) THEN 3139 IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )& 3140 THEN 3133 3141 message_string = 'data_output_pr = ' // & 3134 3142 TRIM( data_output_pr(i) ) // ' is not ava' // & … … 3143 3151 3144 3152 CASE ( 'rad_lw_in' ) 3145 IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' ) THEN 3153 IF ( ( .NOT. radiation) .OR. radiation_scheme == 'constant' ) & 3154 THEN 3146 3155 message_string = 'data_output_pr = ' // & 3147 3156 TRIM( data_output_pr(i) ) // ' is not ava' // & … … 3156 3165 3157 3166 CASE ( 'rad_lw_out' ) 3158 IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' ) THEN 3167 IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) & 3168 THEN 3159 3169 message_string = 'data_output_pr = ' // & 3160 3170 TRIM( data_output_pr(i) ) // ' is not ava' // & … … 3169 3179 3170 3180 CASE ( 'rad_sw_in' ) 3171 IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' ) THEN 3181 IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) & 3182 THEN 3172 3183 message_string = 'data_output_pr = ' // & 3173 3184 TRIM( data_output_pr(i) ) // ' is not ava' // & … … 3182 3193 3183 3194 CASE ( 'rad_sw_out') 3184 IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' ) THEN 3195 IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) & 3196 THEN 3185 3197 message_string = 'data_output_pr = ' // & 3186 3198 TRIM( data_output_pr(i) ) // ' is not ava' // & … … 3195 3207 3196 3208 CASE ( 'rad_lw_cs_hr' ) 3197 IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' ) THEN 3209 IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) & 3210 THEN 3198 3211 message_string = 'data_output_pr = ' // & 3199 3212 TRIM( data_output_pr(i) ) // ' is not ava' // & … … 3208 3221 3209 3222 CASE ( 'rad_lw_hr' ) 3210 IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' ) THEN 3223 IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) & 3224 THEN 3211 3225 message_string = 'data_output_pr = ' // & 3212 3226 TRIM( data_output_pr(i) ) // ' is not ava' // & … … 3221 3235 3222 3236 CASE ( 'rad_sw_cs_hr' ) 3223 IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' ) THEN 3237 IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) & 3238 THEN 3224 3239 message_string = 'data_output_pr = ' // & 3225 3240 TRIM( data_output_pr(i) ) // ' is not ava' // & … … 3234 3249 3235 3250 CASE ( 'rad_sw_hr' ) 3236 IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' ) THEN 3251 IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) & 3252 THEN 3237 3253 message_string = 'data_output_pr = ' // & 3238 3254 TRIM( data_output_pr(i) ) // ' is not ava' // & … … 3329 3345 3330 3346 CASE ( 'lpt' ) 3331 IF ( .NOT.cloud_physics ) THEN3347 IF ( .NOT. cloud_physics ) THEN 3332 3348 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3333 3349 'res cloud_physics = .TRUE.' … … 3337 3353 3338 3354 CASE ( 'm_soil' ) 3339 IF ( .NOT.land_surface ) THEN3355 IF ( .NOT. land_surface ) THEN 3340 3356 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3341 3357 'land_surface = .TRUE.' … … 3345 3361 3346 3362 CASE ( 'nr' ) 3347 IF ( .NOT.cloud_physics ) THEN3363 IF ( .NOT. cloud_physics ) THEN 3348 3364 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3349 3365 'res cloud_physics = .TRUE.' … … 3357 3373 3358 3374 CASE ( 'pc', 'pr' ) 3359 IF ( .NOT.particle_advection ) THEN3375 IF ( .NOT. particle_advection ) THEN 3360 3376 message_string = 'output of "' // TRIM( var ) // '" requir' // & 3361 3377 'es a "particles_par"-NAMELIST in the parameter file (PARIN)' … … 3366 3382 3367 3383 CASE ( 'prr' ) 3368 IF ( .NOT.cloud_physics ) THEN3384 IF ( .NOT. cloud_physics ) THEN 3369 3385 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3370 3386 'res cloud_physics = .TRUE.' … … 3374 3390 'res cloud_scheme = seifert_beheng' 3375 3391 CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 ) 3376 ELSEIF ( .NOT.precipitation ) THEN3392 ELSEIF ( .NOT. precipitation ) THEN 3377 3393 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3378 3394 'res precipitation = .TRUE.' … … 3382 3398 3383 3399 CASE ( 'q', 'vpt' ) 3384 IF ( .NOT.humidity ) THEN3400 IF ( .NOT. humidity ) THEN 3385 3401 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3386 3402 'res humidity = .TRUE.' … … 3391 3407 3392 3408 CASE ( 'qc' ) 3393 IF ( .NOT.cloud_physics ) THEN3409 IF ( .NOT. cloud_physics ) THEN 3394 3410 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3395 3411 'res cloud_physics = .TRUE.' … … 3403 3419 3404 3420 CASE ( 'ql' ) 3405 IF ( .NOT. ( cloud_physics .OR. cloud_droplets ) ) THEN3421 IF ( .NOT. ( cloud_physics .OR. cloud_droplets ) ) THEN 3406 3422 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3407 3423 'res cloud_physics = .TRUE. or cloud_droplets = .TRUE.' … … 3411 3427 3412 3428 CASE ( 'ql_c', 'ql_v', 'ql_vp' ) 3413 IF ( .NOT.cloud_droplets ) THEN3429 IF ( .NOT. cloud_droplets ) THEN 3414 3430 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3415 3431 'res cloud_droplets = .TRUE.' … … 3421 3437 3422 3438 CASE ( 'qr' ) 3423 IF ( .NOT.cloud_physics ) THEN3439 IF ( .NOT. cloud_physics ) THEN 3424 3440 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3425 3441 'res cloud_physics = .TRUE.' … … 3429 3445 'res cloud_scheme = seifert_beheng' 3430 3446 CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 ) 3431 ELSEIF ( .NOT.precipitation ) THEN3447 ELSEIF ( .NOT. precipitation ) THEN 3432 3448 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3433 3449 'res precipitation = .TRUE.' … … 3437 3453 3438 3454 CASE ( 'qv' ) 3439 IF ( .NOT.cloud_physics ) THEN3455 IF ( .NOT. cloud_physics ) THEN 3440 3456 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3441 3457 'res cloud_physics = .TRUE.' … … 3446 3462 CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_lw_cs_hr', 'rad_lw_hr', & 3447 3463 'rad_sw_in', 'rad_sw_out', 'rad_sw_cs_hr', 'rad_sw_hr' ) 3448 IF ( .NOT. radiation .OR.radiation_scheme /= 'rrtmg' ) THEN3464 IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' ) THEN 3449 3465 message_string = '"output of "' // TRIM( var ) // '" requi' // & 3450 3466 'res radiation = .TRUE. and ' // & … … 3455 3471 3456 3472 CASE ( 'rho' ) 3457 IF ( .NOT.ocean ) THEN3473 IF ( .NOT. ocean ) THEN 3458 3474 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3459 3475 'res ocean = .TRUE.' … … 3463 3479 3464 3480 CASE ( 's' ) 3465 IF ( .NOT.passive_scalar ) THEN3481 IF ( .NOT. passive_scalar ) THEN 3466 3482 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3467 3483 'res passive_scalar = .TRUE.' … … 3471 3487 3472 3488 CASE ( 'sa' ) 3473 IF ( .NOT.ocean ) THEN3489 IF ( .NOT. ocean ) THEN 3474 3490 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3475 3491 'res ocean = .TRUE.' … … 3479 3495 3480 3496 CASE ( 't_soil' ) 3481 IF ( .NOT.land_surface ) THEN3497 IF ( .NOT. land_surface ) THEN 3482 3498 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3483 3499 'land_surface = .TRUE.' … … 3491 3507 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*', 'rad_net*', & 3492 3508 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*', & 3493 'r_a*', 'r_s*', 'shf*', 'shf_eb*', 't*', 'u*', 'z0*', 'z0h*' ) 3509 'r_a*', 'r_s*', 'shf*', 'shf_eb*', 't*', 'u*', 'z0*', 'z0h*', & 3510 'z0q*' ) 3494 3511 IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN 3495 3512 message_string = 'illegal value for data_output: "' // & … … 3498 3515 CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 ) 3499 3516 ENDIF 3500 IF ( .NOT. radiation .OR.radiation_scheme /= "rrtmg" ) THEN3501 IF ( TRIM( var ) == 'rrtm_aldif*' .OR.&3502 TRIM( var ) == 'rrtm_aldir*' .OR.&3503 TRIM( var ) == 'rrtm_asdif*' .OR.&3517 IF ( .NOT. radiation .OR. radiation_scheme /= "rrtmg" ) THEN 3518 IF ( TRIM( var ) == 'rrtm_aldif*' .OR. & 3519 TRIM( var ) == 'rrtm_aldir*' .OR. & 3520 TRIM( var ) == 'rrtm_asdif*' .OR. & 3504 3521 TRIM( var ) == 'rrtm_asdir*' ) & 3505 3522 THEN … … 3511 3528 ENDIF 3512 3529 3513 IF ( TRIM( var ) == 'c_liq*' .AND. .NOT. land_surface ) THEN3530 IF ( TRIM( var ) == 'c_liq*' .AND. .NOT. land_surface ) THEN 3514 3531 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3515 3532 'res land_surface = .TRUE.' 3516 3533 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 3517 3534 ENDIF 3518 IF ( TRIM( var ) == 'c_soil*' .AND. .NOT. land_surface ) THEN3535 IF ( TRIM( var ) == 'c_soil*' .AND. .NOT. land_surface ) THEN 3519 3536 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3520 3537 'res land_surface = .TRUE.' … … 3526 3543 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 3527 3544 ENDIF 3528 IF ( TRIM( var ) == 'ghf_eb*' .AND. .NOT. land_surface ) THEN3545 IF ( TRIM( var ) == 'ghf_eb*' .AND. .NOT. land_surface ) THEN 3529 3546 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3530 3547 'res land_surface = .TRUE.' 3531 3548 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 3532 3549 ENDIF 3533 IF ( TRIM( var ) == 'lai*' .AND. .NOT. land_surface ) THEN3550 IF ( TRIM( var ) == 'lai*' .AND. .NOT. land_surface ) THEN 3534 3551 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3535 3552 'res land_surface = .TRUE.' … … 3541 3558 CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 ) 3542 3559 ENDIF 3543 IF ( TRIM( var ) == 'm_liq_eb*' .AND. .NOT. land_surface ) THEN3560 IF ( TRIM( var ) == 'm_liq_eb*' .AND. .NOT. land_surface ) THEN 3544 3561 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3545 3562 'res land_surface = .TRUE.' … … 3556 3573 CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 ) 3557 3574 ENDIF 3558 IF ( TRIM( var ) == 'prr*' .AND. .NOT. precipitation ) THEN3575 IF ( TRIM( var ) == 'prr*' .AND. .NOT. precipitation ) THEN 3559 3576 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3560 3577 'res precipitation = .TRUE.' 3561 3578 CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 ) 3562 3579 ENDIF 3563 IF ( TRIM( var ) == 'qsws*' .AND. .NOT. humidity ) THEN3580 IF ( TRIM( var ) == 'qsws*' .AND. .NOT. humidity ) THEN 3564 3581 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3565 3582 'res humidity = .TRUE.' 3566 3583 CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 ) 3567 3584 ENDIF 3568 IF ( TRIM( var ) == 'qsws_eb*' .AND. .NOT. land_surface ) THEN3585 IF ( TRIM( var ) == 'qsws_eb*' .AND. .NOT. land_surface ) THEN 3569 3586 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3570 3587 'res land_surface = .TRUE.' 3571 3588 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 3572 3589 ENDIF 3573 IF ( TRIM( var ) == 'qsws_liq_eb*' .AND. .NOT. land_surface ) &3590 IF ( TRIM( var ) == 'qsws_liq_eb*' .AND. .NOT. land_surface ) & 3574 3591 THEN 3575 3592 message_string = 'output of "' // TRIM( var ) // '" requi' // & … … 3577 3594 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 3578 3595 ENDIF 3579 IF ( TRIM( var ) == 'qsws_soil_eb*' .AND. .NOT. land_surface ) &3596 IF ( TRIM( var ) == 'qsws_soil_eb*' .AND. .NOT. land_surface ) & 3580 3597 THEN 3581 3598 message_string = 'output of "' // TRIM( var ) // '" requi' // & … … 3583 3600 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 3584 3601 ENDIF 3585 IF ( TRIM( var ) == 'qsws_veg_eb*' .AND. .NOT. land_surface ) &3602 IF ( TRIM( var ) == 'qsws_veg_eb*' .AND. .NOT. land_surface ) & 3586 3603 THEN 3587 3604 message_string = 'output of "' // TRIM( var ) // '" requi' // & … … 3589 3606 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 3590 3607 ENDIF 3591 IF ( TRIM( var ) == 'r_a*' .AND. .NOT. land_surface ) &3608 IF ( TRIM( var ) == 'r_a*' .AND. .NOT. land_surface ) & 3592 3609 THEN 3593 3610 message_string = 'output of "' // TRIM( var ) // '" requi' // & … … 3595 3612 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 3596 3613 ENDIF 3597 IF ( TRIM( var ) == 'r_s*' .AND. .NOT. land_surface ) &3614 IF ( TRIM( var ) == 'r_s*' .AND. .NOT. land_surface ) & 3598 3615 THEN 3599 3616 message_string = 'output of "' // TRIM( var ) // '" requi' // & … … 3757 3774 !-- Check mask conditions 3758 3775 DO mid = 1, max_masks 3759 IF ( data_output_masks(mid,1) /= ' ' .OR.&3776 IF ( data_output_masks(mid,1) /= ' ' .OR. & 3760 3777 data_output_masks_user(mid,1) /= ' ' ) THEN 3761 3778 masks = masks + 1 … … 3763 3780 ENDDO 3764 3781 3765 IF ( masks < 0 .OR.masks > max_masks ) THEN3782 IF ( masks < 0 .OR. masks > max_masks ) THEN 3766 3783 WRITE( message_string, * ) 'illegal value: masks must be >= 0 and ', & 3767 3784 '<= ', max_masks, ' (=max_masks)' … … 3908 3925 ! 3909 3926 !-- Check random generator 3910 IF ( (random_generator /= 'system-specific' .AND.&3911 random_generator /= 'random-parallel' ) .AND.&3927 IF ( (random_generator /= 'system-specific' .AND. & 3928 random_generator /= 'random-parallel' ) .AND. & 3912 3929 random_generator /= 'numerical-recipes' ) THEN 3913 3930 message_string = 'unknown random generator: random_generator = "' // & … … 3919 3936 !-- Determine upper and lower hight level indices for random perturbations 3920 3937 IF ( disturbance_level_b == -9999999.9_wp ) THEN 3921 IF ( ocean ) THEN3938 IF ( ocean ) THEN 3922 3939 disturbance_level_b = zu((nzt*2)/3) 3923 3940 disturbance_level_ind_b = ( nzt * 2 ) / 3 … … 4182 4199 ! 4183 4200 !-- Check pressure gradient conditions 4184 IF ( dp_external .AND.conserve_volume_flow ) THEN4201 IF ( dp_external .AND. conserve_volume_flow ) THEN 4185 4202 WRITE( message_string, * ) 'Both dp_external and conserve_volume_flo', & 4186 4203 'w are .TRUE. but one of them must be .FALSE.' … … 4188 4205 ENDIF 4189 4206 IF ( dp_external ) THEN 4190 IF ( dp_level_b < zu(nzb) .OR.dp_level_b > zu(nzt) ) THEN4207 IF ( dp_level_b < zu(nzb) .OR. dp_level_b > zu(nzt) ) THEN 4191 4208 WRITE( message_string, * ) 'dp_level_b = ', dp_level_b, ' is out ', & 4192 4209 ' of range' … … 4199 4216 ENDIF 4200 4217 ENDIF 4201 IF ( ANY( dpdxy /= 0.0_wp ) .AND. .NOT.dp_external ) THEN4218 IF ( ANY( dpdxy /= 0.0_wp ) .AND. .NOT. dp_external ) THEN 4202 4219 WRITE( message_string, * ) 'dpdxy is nonzero but dp_external is ', & 4203 4220 '.FALSE., i.e. the external pressure gradient & will not be applied' … … 4230 4247 ENDIF 4231 4248 ENDIF 4232 IF ( ( u_bulk /= 0.0_wp .OR. v_bulk /= 0.0_wp ) .AND.&4233 ( .NOT. conserve_volume_flow .OR.&4249 IF ( ( u_bulk /= 0.0_wp .OR. v_bulk /= 0.0_wp ) .AND. & 4250 ( .NOT. conserve_volume_flow .OR. & 4234 4251 TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) ) THEN 4235 4252 WRITE( message_string, * ) 'nonzero bulk velocity requires ', & … … 4270 4287 ! 4271 4288 !-- Check nudging and large scale forcing from external file 4272 IF ( nudging .AND. ( .NOT.large_scale_forcing ) ) THEN4289 IF ( nudging .AND. ( .NOT. large_scale_forcing ) ) THEN 4273 4290 message_string = 'Nudging requires large_scale_forcing = .T.. &'// & 4274 4291 'Surface fluxes and geostrophic wind should be &'// & … … 4277 4294 ENDIF 4278 4295 4279 IF ( large_scale_forcing .AND. ( bc_lr /= 'cyclic' .OR.&4296 IF ( large_scale_forcing .AND. ( bc_lr /= 'cyclic' .OR. & 4280 4297 bc_ns /= 'cyclic' ) ) THEN 4281 4298 message_string = 'Non-cyclic lateral boundaries do not allow for &' // & … … 4284 4301 ENDIF 4285 4302 4286 IF ( large_scale_forcing .AND. ( .NOT.humidity ) ) THEN4303 IF ( large_scale_forcing .AND. ( .NOT. humidity ) ) THEN 4287 4304 message_string = 'The usage of large scale forcing from external &'// & 4288 4305 'file LSF_DATA requires humidity = .T..' … … 4290 4307 ENDIF 4291 4308 4292 IF ( large_scale_forcing .AND.topography /= 'flat' ) THEN4309 IF ( large_scale_forcing .AND. topography /= 'flat' ) THEN 4293 4310 message_string = 'The usage of large scale forcing from external &'// & 4294 4311 'file LSF_DATA is not implemented for non-flat topography' … … 4296 4313 ENDIF 4297 4314 4298 IF ( large_scale_forcing .AND. ocean ) THEN4315 IF ( large_scale_forcing .AND. ocean ) THEN 4299 4316 message_string = 'The usage of large scale forcing from external &'// & 4300 4317 'file LSF_DATA is not implemented for ocean runs' … … 4320 4337 !-- Check for valid setting of most_method 4321 4338 IF ( TRIM( most_method ) /= 'circular' .AND. & 4322 TRIM( most_method ) /= 'newton' .AND.&4339 TRIM( most_method ) /= 'newton' .AND. & 4323 4340 TRIM( most_method ) /= 'lookup' ) THEN 4324 message_string = 'most_method = "' // TRIM( most_method ) // &4341 message_string = 'most_method = "' // TRIM( most_method ) // & 4325 4342 '" is unknown' 4326 4343 CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 ) … … 4352 4369 IF ( dt_do == 0.0_wp ) THEN 4353 4370 IF ( dt_fixed ) THEN 4354 WRITE( message_string, '(A,F9.4,A)' ) 'Output at every ' // &4355 'timestep is desired (' // dt_do_name // ' = 0.0).&'// &4356 'Setting the output interval to the fixed timestep '// &4371 WRITE( message_string, '(A,F9.4,A)' ) 'Output at every ' // & 4372 'timestep is desired (' // dt_do_name // ' = 0.0).&'// & 4373 'Setting the output interval to the fixed timestep '// & 4357 4374 'dt = ', dt, 's.' 4358 4375 CALL message( 'check_parameters', 'PA0060', 0, 0, 0, 6, 0 ) 4359 4376 dt_do = dt 4360 4377 ELSE 4361 message_string = dt_do_name // ' = 0.0 while using a ' // &4362 'variable timestep and parallel netCDF4 ' // &4378 message_string = dt_do_name // ' = 0.0 while using a ' // & 4379 'variable timestep and parallel netCDF4 ' // & 4363 4380 'is not allowed.' 4364 4381 CALL message( 'check_parameters', 'PA0081', 1, 2, 0, 6, 0 ) -
palm/trunk/SOURCE/data_output_2d.f90
r1784 r1788 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! Added output of z0q 21 22 ! 22 !23 23 ! Former revisions: 24 24 ! ----------------- … … 129 129 USE arrays_3d, & 130 130 ONLY: dzw, e, nr, ol, p, pt, q, qc, ql, ql_c, ql_v, ql_vp, qr, qsws, & 131 rho, sa, shf, tend, ts, u, us, v, vpt, w, z0, z0h, z u, zw131 rho, sa, shf, tend, ts, u, us, v, vpt, w, z0, z0h, z0q, zu, zw 132 132 133 133 USE averaging … … 1247 1247 DO j = nysg, nyng 1248 1248 local_pf(i,j,nzb+1) = z0h_av(j,i) 1249 ENDDO 1250 ENDDO 1251 ENDIF 1252 resorted = .TRUE. 1253 two_d = .TRUE. 1254 level_z(nzb+1) = zu(nzb+1) 1255 1256 CASE ( 'z0q*_xy' ) ! 2d-array 1257 IF ( av == 0 ) THEN 1258 DO i = nxlg, nxrg 1259 DO j = nysg, nyng 1260 local_pf(i,j,nzb+1) = z0q(j,i) 1261 ENDDO 1262 ENDDO 1263 ELSE 1264 DO i = nxlg, nxrg 1265 DO j = nysg, nyng 1266 local_pf(i,j,nzb+1) = z0q_av(j,i) 1249 1267 ENDDO 1250 1268 ENDDO -
palm/trunk/SOURCE/header.f90
r1787 r1788 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Parameter dewfall removed 22 22 ! 23 23 ! Former revisions: … … 248 248 249 249 USE land_surface_model_mod, & 250 ONLY: conserve_water_content, dewfall, land_surface, nzb_soil,&250 ONLY: conserve_water_content, land_surface, nzb_soil, & 251 251 nzt_soil, root_fraction, soil_moisture, soil_temperature, & 252 252 soil_type, soil_type_name, veg_type, veg_type_name, zs … … 964 964 ELSE 965 965 WRITE( io, 441 ) 966 ENDIF967 968 IF ( dewfall ) THEN969 WRITE( io, 442 )970 ELSE971 WRITE( io, 443 )972 966 ENDIF 973 967 … … 2341 2335 ' Root fraction: ',A,' '/ & 2342 2336 ' Gridpoint: ',A) 2343 440 FORMAT (/' --> Dewfall is allowed (default)') 2344 441 FORMAT (' --> Dewfall is inhibited') 2345 442 FORMAT (' --> Soil bottom is closed (water content is conserved, default)') 2346 443 FORMAT (' --> Soil bottom is open (water content is not conserved)') 2337 440 FORMAT (' --> Soil bottom is closed (water content is conserved, default)') 2338 441 FORMAT (' --> Soil bottom is open (water content is not conserved)') 2347 2339 444 FORMAT (//' Radiation model information:'/ & 2348 2340 ' ----------------------------'/) -
palm/trunk/SOURCE/init_3d_model.f90
r1784 r1788 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Added z0q. 22 ! Syntax layout improved. 22 23 ! 23 24 ! Former revisions: … … 320 321 ! 321 322 !-- Allocate arrays 322 ALLOCATE( mean_surface_level_height(0:statistic_regions), &323 mean_surface_level_height_l(0:statistic_regions), &324 ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), &325 ngp_3d(0:statistic_regions), &326 ngp_3d_inner(0:statistic_regions), &327 ngp_3d_inner_l(0:statistic_regions), &328 ngp_3d_inner_tmp(0:statistic_regions), &329 sums_divnew_l(0:statistic_regions), &323 ALLOCATE( mean_surface_level_height(0:statistic_regions), & 324 mean_surface_level_height_l(0:statistic_regions), & 325 ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), & 326 ngp_3d(0:statistic_regions), & 327 ngp_3d_inner(0:statistic_regions), & 328 ngp_3d_inner_l(0:statistic_regions), & 329 ngp_3d_inner_tmp(0:statistic_regions), & 330 sums_divnew_l(0:statistic_regions), & 330 331 sums_divold_l(0:statistic_regions) ) 331 332 ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt) ) 332 ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions), &333 ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions), &334 ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions), &335 ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions), &336 rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions), &337 sums(nzb:nzt+1,pr_palm+max_pr_user), &338 sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1), &339 sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), &340 sums_up_fraction_l(10,3,0:statistic_regions), &341 sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions), &333 ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions), & 334 ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions), & 335 ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions), & 336 ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions), & 337 rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions), & 338 sums(nzb:nzt+1,pr_palm+max_pr_user), & 339 sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1), & 340 sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), & 341 sums_up_fraction_l(10,3,0:statistic_regions), & 342 sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions), & 342 343 ts_value(dots_max,0:statistic_regions) ) 343 344 ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) ) 344 345 345 ALLOCATE( ol(nysg:nyng,nxlg:nxrg), shf(nysg:nyng,nxlg:nxrg), &346 ts(nysg:nyng,nxlg:nxrg), tswst(nysg:nyng,nxlg:nxrg), &347 us(nysg:nyng,nxlg:nxrg), usws(nysg:nyng,nxlg:nxrg), &348 uswst(nysg:nyng,nxlg:nxrg), vsws(nysg:nyng,nxlg:nxrg), &349 vswst(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg), &350 z0h(nysg:nyng,nxlg:nxrg) )351 352 ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr), &353 kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &354 km(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &355 p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &346 ALLOCATE( ol(nysg:nyng,nxlg:nxrg), shf(nysg:nyng,nxlg:nxrg), & 347 ts(nysg:nyng,nxlg:nxrg), tswst(nysg:nyng,nxlg:nxrg), & 348 us(nysg:nyng,nxlg:nxrg), usws(nysg:nyng,nxlg:nxrg), & 349 uswst(nysg:nyng,nxlg:nxrg), vsws(nysg:nyng,nxlg:nxrg), & 350 vswst(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg), & 351 z0h(nysg:nyng,nxlg:nxrg), z0q(nysg:nyng,nxlg:nxrg) ) 352 353 ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr), & 354 kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 355 km(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 356 p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 356 357 tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 357 358 358 359 #if defined( __nopointer ) 359 ALLOCATE( e(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &360 e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &361 pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &362 pt_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &363 u(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &364 u_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &365 v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &366 v_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &367 w(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &368 w_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &369 te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &370 tpt_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &371 tu_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &372 tv_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &360 ALLOCATE( e(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 361 e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 362 pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 363 pt_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 364 u(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 365 u_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 366 v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 367 v_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 368 w(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 369 w_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 370 te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 371 tpt_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 372 tu_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 373 tv_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 373 374 tw_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 374 375 #else 375 ALLOCATE( e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &376 e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &377 e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &378 pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &379 pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &380 u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &381 u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &382 u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &383 v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &384 v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &385 v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &386 w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &387 w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &376 ALLOCATE( e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 377 e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 378 e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 379 pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 380 pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 381 u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 382 u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 383 u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 384 v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 385 v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 386 v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 387 w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 388 w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 388 389 w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 389 IF ( .NOT.neutral ) THEN390 IF ( .NOT. neutral ) THEN 390 391 ALLOCATE( pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 391 392 ENDIF … … 412 413 ENDIF 413 414 414 IF ( humidity .OR. passive_scalar ) THEN415 IF ( humidity .OR. passive_scalar ) THEN 415 416 ! 416 417 !-- 2D-humidity/scalar arrays 417 ALLOCATE ( qs(nysg:nyng,nxlg:nxrg), &418 qsws(nysg:nyng,nxlg:nxrg), &418 ALLOCATE ( qs(nysg:nyng,nxlg:nxrg), & 419 qsws(nysg:nyng,nxlg:nxrg), & 419 420 qswst(nysg:nyng,nxlg:nxrg) ) 420 421 … … 422 423 !-- 3D-humidity/scalar arrays 423 424 #if defined( __nopointer ) 424 ALLOCATE( q(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &425 q_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &425 ALLOCATE( q(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 426 q_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 426 427 tq_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 427 428 #else 428 ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &429 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &429 ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 430 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 430 431 q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 431 432 #endif … … 440 441 #endif 441 442 442 IF ( cloud_physics ) THEN443 IF ( cloud_physics ) THEN 443 444 444 445 ! … … 451 452 ! 452 453 !-- Precipitation amount and rate (only needed if output is switched) 453 ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg), &454 ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg), & 454 455 precipitation_rate(nysg:nyng,nxlg:nxrg) ) 455 456 … … 457 458 ! 458 459 !-- 1D-arrays 459 ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1), &460 ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1), & 460 461 q_1d(nzb:nzt+1), qc_1d(nzb:nzt+1) ) 461 462 ! … … 474 475 ! 475 476 !-- 2D-rain water content and rain drop concentration arrays 476 ALLOCATE ( qrs(nysg:nyng,nxlg:nxrg), &477 qrsws(nysg:nyng,nxlg:nxrg), &478 qrswst(nysg:nyng,nxlg:nxrg), &479 nrs(nysg:nyng,nxlg:nxrg), &480 nrsws(nysg:nyng,nxlg:nxrg), &477 ALLOCATE ( qrs(nysg:nyng,nxlg:nxrg), & 478 qrsws(nysg:nyng,nxlg:nxrg), & 479 qrswst(nysg:nyng,nxlg:nxrg), & 480 nrs(nysg:nyng,nxlg:nxrg), & 481 nrsws(nysg:nyng,nxlg:nxrg), & 481 482 nrswst(nysg:nyng,nxlg:nxrg) ) 482 483 ! 483 484 !-- 3D-rain water content, rain drop concentration arrays 484 485 #if defined( __nopointer ) 485 ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &486 nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &487 qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &488 qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &489 tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &486 ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 487 nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 488 qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 489 qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 490 tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 490 491 tqr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 491 492 #else 492 ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &493 nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &494 nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &495 qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &496 qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &493 ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 494 nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 495 nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 496 qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 497 qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 497 498 qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 498 499 #endif … … 509 510 !-- Liquid water content, change in liquid water content 510 511 #if defined( __nopointer ) 511 ALLOCATE ( ql(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &512 ALLOCATE ( ql(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 512 513 ql_c(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 513 514 #else 514 ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &515 ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 515 516 ql_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 516 517 #endif 517 518 ! 518 519 !-- Real volume of particles (with weighting), volume of particles 519 ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &520 ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 520 521 ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 521 522 ENDIF … … 526 527 527 528 IF ( ocean ) THEN 528 ALLOCATE( saswsb(nysg:nyng,nxlg:nxrg), &529 ALLOCATE( saswsb(nysg:nyng,nxlg:nxrg), & 529 530 saswst(nysg:nyng,nxlg:nxrg) ) 530 531 #if defined( __nopointer ) 531 ALLOCATE( prho(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &532 rho(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &533 sa(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &534 sa_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &532 ALLOCATE( prho(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 533 rho(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 534 sa(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 535 sa_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 535 536 tsa_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 536 537 #else 537 ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &538 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &539 sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &540 sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &538 ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 539 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 540 sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 541 sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 541 542 sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 542 543 prho => prho_1 … … 553 554 !-- 3D-array for storing the dissipation, needed for calculating the sgs 554 555 !-- particle velocities 555 IF ( use_sgs_for_particles .OR. wang_kernel .OR. turbulence .OR. &556 IF ( use_sgs_for_particles .OR. wang_kernel .OR. turbulence .OR. & 556 557 num_acc_per_node > 0 ) THEN 557 558 ALLOCATE( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) … … 559 560 560 561 IF ( dt_dosp /= 9999999.9_wp ) THEN 561 ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ), &562 ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ), & 562 563 spectrum_y( 1:ny/2, 1:10, 1:10 ) ) 563 564 spectrum_x = 0.0_wp … … 595 596 !-- are needed for radiation boundary conditions 596 597 IF ( outflow_l ) THEN 597 ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2), &598 v_m_l(nzb:nzt+1,nysg:nyng,0:1), &598 ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2), & 599 v_m_l(nzb:nzt+1,nysg:nyng,0:1), & 599 600 w_m_l(nzb:nzt+1,nysg:nyng,0:1) ) 600 601 ENDIF 601 602 IF ( outflow_r ) THEN 602 ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &603 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &603 ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), & 604 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), & 604 605 w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) ) 605 606 ENDIF 606 607 IF ( outflow_l .OR. outflow_r ) THEN 607 ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng), &608 ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng), & 608 609 c_w(nzb:nzt+1,nysg:nyng) ) 609 610 ENDIF 610 611 IF ( outflow_s ) THEN 611 ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg), &612 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg), &612 ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg), & 613 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg), & 613 614 w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) ) 614 615 ENDIF 615 616 IF ( outflow_n ) THEN 616 ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &617 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &617 ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), & 618 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), & 618 619 w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) ) 619 620 ENDIF 620 621 IF ( outflow_s .OR. outflow_n ) THEN 621 ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg), &622 ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg), & 622 623 c_w(nzb:nzt+1,nxlg:nxrg) ) 623 624 ENDIF … … 682 683 !-- pres may need to be called outside the RK-substeps. Antti will 683 684 !-- check if this is really required. 684 ALLOCATE( weight_substep(1:intermediate_timestep_count_max), &685 ALLOCATE( weight_substep(1:intermediate_timestep_count_max), & 685 686 weight_pres(0:intermediate_timestep_count_max) ) 686 687 weight_substep = 1.0_wp … … 691 692 ! 692 693 !-- Initialize model variables 693 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &694 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. & 694 695 TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN 695 696 ! … … 819 820 CALL location_message( 'finished', .TRUE. ) 820 821 821 ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &822 ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) & 822 823 THEN 823 824 … … 825 826 ! 826 827 !-- Overwrite initial profiles in case of nudging 827 IF ( nudging ) THEN828 IF ( nudging ) THEN 828 829 pt_init = ptnudge(:,1) 829 830 u_init = unudge(:,1) … … 833 834 ENDIF 834 835 835 WRITE( message_string, * ) 'Initial profiles of u, v and ', &836 WRITE( message_string, * ) 'Initial profiles of u, v and ', & 836 837 'scalars from NUDGING_DATA are used.' 837 838 CALL message( 'init_3d_model', 'PA0370', 0, 0, 0, 6, 0 ) … … 936 937 CALL location_message( 'finished', .TRUE. ) 937 938 938 ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) &939 ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) & 939 940 THEN 940 941 … … 1018 1019 DO j=0,ny 1019 1020 DO i=0,nx 1020 id_random_array(j,i) = random_dummy + 1E6 * ( ensemble_member_nr - 1000 ) 1021 id_random_array(j,i) = random_dummy + 1E6 & 1022 * ( ensemble_member_nr - 1000 ) 1021 1023 random_dummy = random_dummy + 1 1022 1024 END DO … … 1047 1049 ! 1048 1050 !-- Initialize shf with data from external file LSF_DATA 1049 IF ( large_scale_forcing .AND. lsf_surf ) THEN1051 IF ( large_scale_forcing .AND. lsf_surf ) THEN 1050 1052 CALL ls_forcing_surf ( simulated_time ) 1051 1053 ENDIF … … 1068 1070 !-- Determine the near-surface water flux 1069 1071 IF ( humidity .OR. passive_scalar ) THEN 1070 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &1072 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 1071 1073 precipitation ) THEN 1072 1074 qrsws = 0.0_wp … … 1106 1108 IF ( humidity .OR. passive_scalar ) THEN 1107 1109 qswst = 0.0_wp 1108 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &1110 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 1109 1111 precipitation ) THEN 1110 1112 nrswst = 0.0_wp … … 1133 1135 z0 = roughness_length 1134 1136 z0h = z0h_factor * z0 1137 z0q = z0h_factor * z0 1135 1138 1136 1139 IF ( .NOT. constant_heatflux ) THEN … … 1145 1148 1146 1149 IF ( humidity .OR. passive_scalar ) THEN 1147 IF ( .NOT.constant_waterflux ) qsws = 0.0_wp1148 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &1150 IF ( .NOT. constant_waterflux ) qsws = 0.0_wp 1151 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 1149 1152 precipitation ) THEN 1150 1153 qrsws = 0.0_wp … … 1159 1162 !-- the reference state will be set (overwritten) in init_ocean) 1160 1163 IF ( use_single_reference_value ) THEN 1161 IF ( .NOT.humidity ) THEN1164 IF ( .NOT. humidity ) THEN 1162 1165 ref_state(:) = pt_reference 1163 1166 ELSE … … 1165 1168 ENDIF 1166 1169 ELSE 1167 IF ( .NOT.humidity ) THEN1170 IF ( .NOT. humidity ) THEN 1168 1171 ref_state(:) = pt_init(:) 1169 1172 ELSE … … 1216 1219 !-- If required, change the surface humidity/scalar at the start of the 3D 1217 1220 !-- run 1218 IF ( ( humidity .OR. passive_scalar ) .AND. &1221 IF ( ( humidity .OR. passive_scalar ) .AND. & 1219 1222 q_surface_initial_change /= 0.0_wp ) THEN 1220 1223 q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change … … 1229 1232 tq_m = 0.0_wp 1230 1233 q_p = q 1231 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &1234 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 1232 1235 precipitation ) THEN 1233 1236 tqr_m = 0.0_wp … … 1245 1248 CALL location_message( 'finished', .TRUE. ) 1246 1249 1247 ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. &1248 TRIM( initializing_actions ) == 'cyclic_fill' ) &1250 ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. & 1251 TRIM( initializing_actions ) == 'cyclic_fill' ) & 1249 1252 THEN 1250 1253 … … 1282 1285 ! 1283 1286 !-- Initialization of the turbulence recycling method 1284 IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. &1287 IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. & 1285 1288 turbulent_inflow ) THEN 1286 1289 ! … … 1346 1349 inflow_damping_height = hom_sum(nzb+6,pr_palm,0) 1347 1350 ELSE 1348 WRITE( message_string, * ) 'inflow_damping_height must be ', &1349 'explicitly specified because&the inversion height ', &1351 WRITE( message_string, * ) 'inflow_damping_height must be ', & 1352 'explicitly specified because&the inversion height ', & 1350 1353 'calculated by the prerun is zero.' 1351 1354 CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 ) … … 1382 1385 ! 1383 1386 !-- Inside buildings set velocities and TKE back to zero 1384 IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. &1387 IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. & 1385 1388 topography /= 'flat' ) THEN 1386 1389 ! … … 1415 1418 IF ( humidity .OR. passive_scalar ) THEN 1416 1419 q_p = q 1417 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &1420 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 1418 1421 precipitation ) THEN 1419 1422 qr_p = qr … … 1430 1433 IF ( humidity .OR. passive_scalar ) THEN 1431 1434 tq_m = 0.0_wp 1432 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &1435 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 1433 1436 precipitation ) THEN 1434 1437 tqr_m = 0.0_wp … … 1486 1489 DO j = nys, nyn 1487 1490 DO k = nzb_2d(j,nx)+1, nzt 1488 volume_flow_initial_l(1) = volume_flow_initial_l(1) + &1491 volume_flow_initial_l(1) = volume_flow_initial_l(1) + & 1489 1492 u_init(k) * dzw(k) 1490 1493 volume_flow_area_l(1) = volume_flow_area_l(1) + dzw(k) … … 1522 1525 DO j = nys, nyn 1523 1526 DO k = nzb_2d(j,nx)+1, nzt 1524 volume_flow_initial_l(1) = volume_flow_initial_l(1) + &1527 volume_flow_initial_l(1) = volume_flow_initial_l(1) + & 1525 1528 hom_sum(k,1,0) * dzw(k) 1526 1529 volume_flow_area_l(1) = volume_flow_area_l(1) + dzw(k) … … 1532 1535 DO i = nxl, nxr 1533 1536 DO k = nzb_2d(ny,i)+1, nzt 1534 volume_flow_initial_l(2) = volume_flow_initial_l(2) + &1537 volume_flow_initial_l(2) = volume_flow_initial_l(2) + & 1535 1538 hom_sum(k,2,0) * dzw(k) 1536 1539 volume_flow_area_l(2) = volume_flow_area_l(2) + dzw(k) … … 1568 1571 DO i = nxl, nxr 1569 1572 DO k = nzb_2d(ny,i)+1, nzt 1570 volume_flow_initial_l(2) = volume_flow_initial_l(2) + &1573 volume_flow_initial_l(2) = volume_flow_initial_l(2) + & 1571 1574 v(k,ny,i) * dzw(k) 1572 1575 volume_flow_area_l(2) = volume_flow_area_l(2) + dzw(k) … … 1605 1608 !-- Impose random perturbation on the horizontal velocity field and then 1606 1609 !-- remove the divergences from the velocity field at the initial stage 1607 IF ( create_disturbances .AND. disturbance_energy_limit /= 0.0_wp .AND. &1608 TRIM( initializing_actions ) /= 'read_restart_data' .AND. &1610 IF ( create_disturbances .AND. disturbance_energy_limit /= 0.0_wp .AND. & 1611 TRIM( initializing_actions ) /= 'read_restart_data' .AND. & 1609 1612 TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN 1610 1613 … … 1718 1721 rdf_sc = 0.0_wp 1719 1722 IF ( rayleigh_damping_factor /= 0.0_wp ) THEN 1720 IF ( .NOT.ocean ) THEN1723 IF ( .NOT. ocean ) THEN 1721 1724 DO k = nzb+1, nzt 1722 1725 IF ( zu(k) >= rayleigh_damping_height ) THEN 1723 rdf(k) = rayleigh_damping_factor * &1726 rdf(k) = rayleigh_damping_factor * & 1724 1727 ( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height ) & 1725 / ( zu(nzt) - rayleigh_damping_height ) )&1728 / ( zu(nzt) - rayleigh_damping_height ) ) & 1726 1729 )**2 1727 1730 ENDIF … … 1730 1733 DO k = nzt, nzb+1, -1 1731 1734 IF ( zu(k) <= rayleigh_damping_height ) THEN 1732 rdf(k) = rayleigh_damping_factor * &1735 rdf(k) = rayleigh_damping_factor * & 1733 1736 ( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) ) & 1734 / ( rayleigh_damping_height - zu(nzb+1)))&1737 / ( rayleigh_damping_height - zu(nzb+1) ) ) & 1735 1738 )**2 1736 1739 ENDIF … … 1774 1777 ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp * & 1775 1778 REAL( pt_damping_width - i * dx, KIND=wp ) / ( & 1776 REAL( pt_damping_width, KIND=wp ) 1779 REAL( pt_damping_width, KIND=wp ) ) ) )**2 1777 1780 ENDIF 1778 1781 ENDDO … … 1831 1834 1832 1835 IF ( dots_num > dots_max ) THEN 1833 WRITE( message_string, * ) 'number of time series quantities exceeds', &1834 ' its maximum of dots_max = ', dots_max, &1836 WRITE( message_string, * ) 'number of time series quantities exceeds', & 1837 ' its maximum of dots_max = ', dots_max, & 1835 1838 ' &Please increase dots_max in modules.f90.' 1836 1839 CALL message( 'init_3d_model', 'PA0194', 1, 2, 0, 6, 0 ) … … 1869 1872 !-- All xy-grid points 1870 1873 ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1 1871 mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) +&1872 zw(nzb_s_inner(j,i))1874 mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) & 1875 + zw(nzb_s_inner(j,i)) 1873 1876 ! 1874 1877 !-- xy-grid points above topography … … 1881 1884 ! 1882 1885 !-- All grid points of the total domain above topography 1883 ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) +&1884 ( nz - nzb_s_inner(j,i) + 2 )1886 ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) & 1887 + ( nz - nzb_s_inner(j,i) + 2 ) 1885 1888 ENDIF 1886 1889 ENDDO … … 1891 1894 #if defined( __parallel ) 1892 1895 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1893 CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM, &1896 CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM, & 1894 1897 comm2d, ierr ) 1895 1898 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1896 CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr, &1899 CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr, & 1897 1900 MPI_INTEGER, MPI_SUM, comm2d, ierr ) 1898 1901 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1899 CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0), &1902 CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0), & 1900 1903 (nz+2)*sr, MPI_INTEGER, MPI_SUM, comm2d, ierr ) 1901 1904 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1902 CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL, &1905 CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL, & 1903 1906 MPI_SUM, comm2d, ierr ) 1904 1907 ngp_3d_inner = INT( ngp_3d_inner_tmp, KIND = SELECTED_INT_KIND( 18 ) ) 1905 1908 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1906 CALL MPI_ALLREDUCE( mean_surface_level_height_l(0), &1907 mean_surface_level_height(0), sr, MPI_REAL, &1909 CALL MPI_ALLREDUCE( mean_surface_level_height_l(0), & 1910 mean_surface_level_height(0), sr, MPI_REAL, & 1908 1911 MPI_SUM, comm2d, ierr ) 1909 1912 mean_surface_level_height = mean_surface_level_height / REAL( ngp_2dh ) … … 1924 1927 !-- the respective subdomain lie below the surface topography 1925 1928 ngp_2dh_outer = MAX( 1, ngp_2dh_outer(:,:) ) 1926 ngp_3d_inner = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )), &1929 ngp_3d_inner = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )), & 1927 1930 ngp_3d_inner(:) ) 1928 1931 ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) ) 1929 1932 1930 DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l, &1933 DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l, & 1931 1934 ngp_3d_inner_l, ngp_3d_inner_tmp ) 1932 1935 -
palm/trunk/SOURCE/land_surface_model.f90
r1784 r1788 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! Bugfix: calculate lambda_surface based on temperature gradient between skin 22 ! layer and soil layer instead of Obukhov length 23 ! Changed: moved calculation of surface specific humidity to energy balance solver 24 ! New: water surfaces are available by using a fixed sea surface temperature. 25 ! The roughness lengths are calculated dynamically using the Charnock 26 ! parameterization. This involves the new roughness length for moisture z0q. 27 ! New: modified solution of the energy balance solver and soil model for 28 ! paved surfaces (i.e. asphalt concrete). 29 ! Syntax layout improved. 30 ! Changed: parameter dewfall removed. 21 31 ! 22 !23 32 ! Former revisions: 24 33 ! ----------------- … … 105 114 !> @todo Consider partial absorption of the net shortwave radiation by the 106 115 !> skin layer. 107 !> @todo Allow for water surfaces116 !> @todo Improve surface water parameterization 108 117 !> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0, 109 118 !> nzt_soil=3)). … … 118 127 119 128 USE arrays_3d, & 120 ONLY: hyp, ol, pt, pt_p, q, q_p, ql, qsws, shf, ts, us, vpt, z0, z0h 129 ONLY: hyp, ol, pt, pt_p, q, q_p, ql, qsws, shf, ts, us, vpt, z0, z0h, & 130 z0q 121 131 122 132 USE cloud_parameters, & … … 171 181 ! 172 182 !-- LSM variables 173 INTEGER(iwp) :: veg_type = 2, & !< vegetation type, 0: user-defined, 1-19: generic (see list) 174 soil_type = 3 !< soil type, 0: user-defined, 1-6: generic (see list) 183 INTEGER(iwp) :: veg_type = 2, & !< NAMELIST veg_type_2d 184 soil_type = 3 !< NAMELIST soil_type_2d 185 186 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: soil_type_2d, & !< soil type, 0: user-defined, 1-7: generic (see list) 187 veg_type_2d !< vegetation type, 0: user-defined, 1-19: generic (see list) 188 189 LOGICAL, DIMENSION(:,:), ALLOCATABLE :: water_surface, & !< flag parameter for water surfaces (classes 14+15) 190 pave_surface, & !< flag parameter for pavements (asphalt etc.) (class 20) 191 building_surface !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet) 175 192 176 193 LOGICAL :: conserve_water_content = .TRUE., & !< open or closed bottom surface for the soil model 177 dewfall = .TRUE., & !< allow/inhibit dewfall178 194 force_radiation_call_l = .FALSE., & !< flag parameter for unscheduled radiation model calls 179 195 land_surface = .FALSE. !< flag parameter indicating wheather the lsm is used … … 200 216 m_total = 0.0_wp, & !< weighted total water content of the soil (m3/m3) 201 217 n_vangenuchten = 9999999.9_wp, & !< NAMELIST n_vg 218 pave_depth = 9999999.9_wp, & !< depth of the pavement 219 pave_heat_capacity = 1.94E6_wp, & !< volumetric heat capacity of pavement (e.g. roads) 220 pave_heat_conductivity = 1.00_wp, & !< heat conductivity for pavements (e.g. roads) 202 221 q_s = 0.0_wp, & !< saturation specific humidity 203 222 residual_moisture = 9999999.9_wp, & !< NAMELIST m_res … … 210 229 wilting_point = 9999999.9_wp, & !< NAMELIST m_wilt 211 230 z0_eb = 9999999.9_wp, & !< NAMELIST z0 (lsm_par) 212 z0h_eb = 9999999.9_wp !< NAMELIST z0h (lsm_par) 231 z0h_eb = 9999999.9_wp, & !< NAMELIST z0h (lsm_par) 232 z0q_eb = 9999999.9_wp !< NAMELIST z0q (lsm_par) 213 233 214 234 REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: & … … 295 315 296 316 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 297 lambda_h, & !< heat conductivity of soil ( ?)317 lambda_h, & !< heat conductivity of soil (W/m/K) 298 318 lambda_w, & !< hydraulic diffusivity of soil (?) 299 gamma_w, & !< hydraulic conductivity of soil ( ?)319 gamma_w, & !< hydraulic conductivity of soil (W/m/K) 300 320 rho_c_total !< volumetric heat capacity of the actual soil matrix (?) 301 321 … … 327 347 ! 328 348 !-- Predefined Land surface classes (veg_type) 329 CHARACTER(26), DIMENSION(0: 19), PARAMETER :: veg_type_name = (/ &349 CHARACTER(26), DIMENSION(0:20), PARAMETER :: veg_type_name = (/ & 330 350 'user defined ', & ! 0 331 351 'crops, mixed farming ', & ! 1 … … 347 367 'deciduous shrubs ', & ! 17 348 368 'mixed forest/woodland ', & ! 18 349 'interrupted forest ' & ! 19 369 'interrupted forest ', & ! 19 370 'pavements/roads ' & ! 20 350 371 /) 351 372 … … 368 389 !-- Land surface parameters I 369 390 !-- r_canopy_min, lai, c_veg, g_d 370 REAL(wp), DIMENSION(0:3,1: 19), PARAMETER :: veg_pars = RESHAPE( (/ &391 REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: veg_pars = RESHAPE( (/ & 371 392 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 1 372 393 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp, & ! 2 … … 387 408 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp, & ! 17 388 409 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 18 389 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp & ! 19 390 /), (/ 4, 19 /) ) 391 392 ! 393 !-- Land surface parameters II z0, z0h 394 REAL(wp), DIMENSION(0:1,1:19), PARAMETER :: roughness_par = RESHAPE( (/ & 395 0.25_wp, 0.25E-2_wp, & ! 1 396 0.20_wp, 0.20E-2_wp, & ! 2 397 2.00_wp, 2.00_wp, & ! 3 398 2.00_wp, 2.00_wp, & ! 4 399 2.00_wp, 2.00_wp, & ! 5 400 2.00_wp, 2.00_wp, & ! 6 401 0.47_wp, 0.47E-2_wp, & ! 7 402 0.013_wp, 0.013E-2_wp, & ! 8 403 0.034_wp, 0.034E-2_wp, & ! 9 404 0.5_wp, 0.50E-2_wp, & ! 10 405 0.17_wp, 0.17E-2_wp, & ! 11 406 1.3E-3_wp, 1.3E-4_wp, & ! 12 407 0.83_wp, 0.83E-2_wp, & ! 13 408 0.00_wp, 0.00E-2_wp, & ! 14 409 0.00_wp, 0.00E-2_wp, & ! 15 410 0.10_wp, 0.10E-2_wp, & ! 16 411 0.25_wp, 0.25E-2_wp, & ! 17 412 2.00_wp, 2.00E-2_wp, & ! 18 413 1.10_wp, 1.10E-2_wp & ! 19 414 /), (/ 2, 19 /) ) 410 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp, & ! 19 411 0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp & ! 20 412 /), (/ 4, 20 /) ) 413 414 ! 415 !-- Land surface parameters II z0, z0h, z0q 416 REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: roughness_par = RESHAPE( (/ & 417 0.25_wp, 0.25E-2_wp, 0.25E-2_wp, & ! 1 418 0.20_wp, 0.20E-2_wp, 0.20E-2_wp, & ! 2 419 2.00_wp, 2.00_wp, 2.00_wp, & ! 3 420 2.00_wp, 2.00_wp, 2.00_wp, & ! 4 421 2.00_wp, 2.00_wp, 2.00_wp, & ! 5 422 2.00_wp, 2.00_wp, 2.00_wp, & ! 6 423 0.47_wp, 0.47E-2_wp, 0.47E-2_wp, & ! 7 424 0.013_wp, 0.013E-2_wp, 0.013E-2_wp, & ! 8 425 0.034_wp, 0.034E-2_wp, 0.034E-2_wp, & ! 9 426 0.5_wp, 0.50E-2_wp, 0.50E-2_wp, & ! 10 427 0.17_wp, 0.17E-2_wp, 0.17E-2_wp, & ! 11 428 1.3E-3_wp, 1.3E-4_wp, 1.3E-4_wp, & ! 12 429 0.83_wp, 0.83E-2_wp, 0.83E-2_wp, & ! 13 430 0.00_wp, 0.00_wp, 0.00_wp, & ! 14 431 0.00_wp, 0.00_wp, 0.00_wp, & ! 15 432 0.10_wp, 0.10E-2_wp, 0.10E-2_wp, & ! 16 433 0.25_wp, 0.25E-2_wp, 0.25E-2_wp, & ! 17 434 2.00_wp, 2.00E-2_wp, 2.00E-2_wp, & ! 18 435 1.10_wp, 1.10E-2_wp, 1.10E-2_wp, & ! 19 436 1.0E-4_wp, 1.0E-5_wp, 1.0E-5_wp & ! 20 437 /), (/ 3, 20 /) ) 415 438 416 439 ! 417 440 !-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in 418 REAL(wp), DIMENSION(0:2,1: 19), PARAMETER :: surface_pars = RESHAPE( (/ &441 REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: surface_pars = RESHAPE( (/ & 419 442 10.0_wp, 10.0_wp, 0.05_wp, & ! 1 420 443 10.0_wp, 10.0_wp, 0.05_wp, & ! 2 … … 430 453 58.0_wp, 58.0_wp, 0.00_wp, & ! 12 431 454 10.0_wp, 10.0_wp, 0.05_wp, & ! 13 432 1.0E 20_wp, 1.0E20_wp, 0.00_wp, & ! 14433 1.0E 20_wp, 1.0E20_wp, 0.00_wp, & ! 15455 1.0E10_wp, 1.0E10_wp, 0.00_wp, & ! 14 456 1.0E10_wp, 1.0E10_wp, 0.00_wp, & ! 15 434 457 10.0_wp, 10.0_wp, 0.05_wp, & ! 16 435 458 10.0_wp, 10.0_wp, 0.05_wp, & ! 17 436 459 20.0_wp, 15.0_wp, 0.03_wp, & ! 18 437 20.0_wp, 15.0_wp, 0.03_wp & ! 19 438 /), (/ 3, 19 /) ) 460 20.0_wp, 15.0_wp, 0.03_wp, & ! 19 461 0.0_wp, 0.0_wp, 0.00_wp & ! 20 462 /), (/ 3, 20 /) ) 439 463 440 464 ! 441 465 !-- Root distribution (sum = 1) level 1, level 2, level 3, level 4, 442 REAL(wp), DIMENSION(0:3,1: 19), PARAMETER :: root_distribution = RESHAPE( (/ &466 REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: root_distribution = RESHAPE( (/ & 443 467 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 1 444 468 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & ! 2 … … 459 483 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 17 460 484 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 18 461 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp & ! 19 462 /), (/ 4, 19 /) ) 485 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 19 486 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp & ! 20 487 /), (/ 4, 20 /) ) 463 488 464 489 ! … … 499 524 !-- Public parameters, constants and initial values 500 525 PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient, & 501 conserve_water_content, dewfall, field_capacity,&526 conserve_water_content, field_capacity, & 502 527 f_shortwave_incoming, hydraulic_conductivity, init_lsm, & 503 528 init_lsm_arrays, lambda_surface_stable, lambda_surface_unstable, & 504 529 land_surface, leaf_area_index, lsm_energy_balance, lsm_soil_model, & 505 530 lsm_swap_timelevel, l_vangenuchten, min_canopy_resistance, & 506 min_soil_resistance, n_vangenuchten, residual_moisture, rho_cp, & 531 min_soil_resistance, n_vangenuchten, pave_heat_capacity, & 532 pave_depth, pave_heat_conductivity, residual_moisture, rho_cp, & 507 533 rho_lv, root_fraction, saturation_moisture, skip_time_do_lsm, & 508 534 soil_moisture, soil_temperature, soil_type, soil_type_name, & 509 535 vegetation_coverage, veg_type, veg_type_name, wilting_point, z0_eb, & 510 z0h_eb 536 z0h_eb, z0q_eb 511 537 512 538 ! … … 586 612 !-- Allocate 2D vegetation model arrays 587 613 ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) ) 614 ALLOCATE ( building_surface(nysg:nyng,nxlg:nxrg) ) 588 615 ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) ) 589 616 ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) ) … … 601 628 ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) ) 602 629 ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) ) 630 ALLOCATE ( pave_surface(nysg:nyng,nxlg:nxrg) ) 603 631 ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) ) 604 632 ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) ) … … 613 641 ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) ) 614 642 ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) ) 643 ALLOCATE ( soil_type_2d(nysg:nyng,nxlg:nxrg) ) 644 ALLOCATE ( veg_type_2d(nysg:nyng,nxlg:nxrg) ) 645 ALLOCATE ( water_surface(nysg:nyng,nxlg:nxrg) ) 615 646 616 647 #if ! defined( __nopointer ) … … 650 681 ! 651 682 !-- If no cloud physics is used, rho_surface has not been calculated before 652 IF ( .NOT.cloud_physics ) THEN683 IF ( .NOT. cloud_physics ) THEN 653 684 rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn ) 654 685 ENDIF … … 762 793 ENDIF 763 794 795 ! 796 !-- Map values to the respective 2D arrays 764 797 alpha_vg = alpha_vangenuchten 765 798 l_vg = l_vangenuchten … … 796 829 797 830 ! 798 !-- Set artifical values for ts and us so that r_a has its initial value for799 !-- the first time step831 !-- Set artifical values for ts and us so that r_a has its initial value 832 !-- for the first time step 800 833 DO i = nxlg, nxrg 801 834 DO j = nysg, nyng … … 864 897 z0h_eb = roughness_par(1,veg_type) 865 898 ENDIF 866 z0h_factor = z0h_eb / z0_eb 899 IF ( z0q_eb == 9999999.9_wp ) THEN 900 z0q_eb = roughness_par(2,veg_type) 901 ENDIF 902 z0h_factor = z0h_eb / ( z0_eb + 1.0E-20_wp ) 867 903 868 904 IF ( ANY( root_fraction == 9999999.9_wp ) ) THEN … … 881 917 z0h_eb = z0_eb * z0h_factor 882 918 ENDIF 919 IF ( z0q_eb == 9999999.9_wp ) THEN 920 z0q_eb = z0_eb * z0h_factor 921 ENDIF 883 922 884 923 ENDIF 885 924 886 925 ! 887 !-- Initialize vegetation 926 !-- For surfaces covered with pavement, set depth of the pavement (with dry 927 !-- soil below). The depth must be greater than the first soil layer depth 928 IF ( veg_type == 20 ) THEN 929 IF ( pave_depth == 9999999.9_wp ) THEN 930 pave_depth = zs(nzb_soil) 931 ELSE 932 pave_depth = MAX( zs(nzb_soil), pave_depth ) 933 ENDIF 934 ENDIF 935 936 ! 937 !-- Map vegetation and soil types to 2D array to allow for heterogeneous 938 !-- surfaces via user interface see below 939 veg_type_2d = veg_type 940 soil_type_2d = soil_type 941 942 ! 943 !-- Map vegetation parameters to the respective 2D arrays 888 944 r_canopy_min = min_canopy_resistance 889 945 lai = leaf_area_index … … 895 951 z0 = z0_eb 896 952 z0h = z0h_eb 953 z0q = z0q_eb 897 954 898 955 ! … … 900 957 CALL user_init_land_surface 901 958 959 ! 960 !-- Set flag parameter if vegetation type was set to a water surface. Also 961 !-- set temperature to a constant value in all "soil" layers. 962 DO i = nxlg, nxrg 963 DO j = nysg, nyng 964 IF ( veg_type_2d(j,i) == 14 .OR. veg_type_2d(j,i) == 15 ) THEN 965 water_surface(j,i) = .TRUE. 966 ELSEIF ( veg_type_2d(j,i) == 20 ) THEN 967 pave_surface(j,i) = .TRUE. 968 m_soil(:,j,i) = 0.0_wp 969 ENDIF 970 971 ENDDO 972 ENDDO 973 974 ! 975 !-- Calculate new roughness lengths (for water surfaces only) 976 CALL calc_z0_water_surface 902 977 903 978 t_soil_p = t_soil … … 935 1010 INTEGER(iwp) :: k, ks !< running index 936 1011 937 REAL(wp) :: f1, & !< resistance correction term 1 1012 REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface 1013 f1, & !< resistance correction term 1 938 1014 f2, & !< resistance correction term 2 939 1015 f3, & !< resistance correction term 3 … … 965 1041 966 1042 ! 967 !-- Set lambda_surface according to stratification 968 IF ( ol(j,i) >= 0.0_wp ) THEN 969 lambda_surface = lambda_surface_s(j,i) 1043 !-- Set lambda_surface according to stratification between skin layer and soil 1044 IF ( .NOT. pave_surface(j,i) ) THEN 1045 1046 c_surface_tmp = c_surface 1047 1048 IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i)) THEN 1049 lambda_surface = lambda_surface_s(j,i) 1050 ELSE 1051 lambda_surface = lambda_surface_u(j,i) 1052 ENDIF 970 1053 ELSE 971 lambda_surface = lambda_surface_u(j,i) 1054 1055 c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp 1056 lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil) 1057 972 1058 ENDIF 973 1059 … … 1016 1102 ENDDO 1017 1103 1018 IF ( m_total > m_wilt(j,i) .AND.m_total < m_fc(j,i) ) THEN1104 IF ( m_total > m_wilt(j,i) .AND. m_total < m_fc(j,i) ) THEN 1019 1105 f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) ) 1020 1106 ELSEIF ( m_total >= m_fc(j,i) ) THEN … … 1050 1136 !-- Third step: calculate bare soil resistance r_soil. The Clapp & 1051 1137 !-- Hornberger parametrization does not consider c_veg. 1052 IF ( soil_type /= 7 ) THEN1138 IF ( soil_type_2d(j,i) /= 7 ) THEN 1053 1139 m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) * & 1054 1140 m_res(j,i) … … 1064 1150 1065 1151 ! 1066 !-- Calculate fraction of liquid water reservoir 1067 m_liq_eb_max = m_max_depth * lai(j,i) 1068 c_liq(j,i) = MIN(1.0_wp, m_liq_eb(j,i) / (m_liq_eb_max+1.0E-20_wp)) 1069 1152 !-- Calculate the maximum possible liquid water amount on plants and 1153 !-- bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is 1154 !-- assumed, while paved surfaces might hold up 1 mm of water. The 1155 !-- liquid water fraction for paved surfaces is calculated after 1156 !-- Noilhan & Planton (1989), while the ECMWF formulation is used for 1157 !-- vegetated surfaces and bare soils. 1158 IF ( pave_surface(j,i) ) THEN 1159 m_liq_eb_max = m_max_depth * 5.0_wp 1160 c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 ) 1161 ELSE 1162 m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i) & 1163 + (1.0_wp - c_veg(j,i)) ) 1164 c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max ) 1165 ENDIF 1070 1166 1071 1167 ! … … 1076 1172 !-- In case of dewfall, set evapotranspiration to zero 1077 1173 !-- All super-saturated water is then removed from the air 1078 IF ( humidity .AND.q_s <= qv1 ) THEN1174 IF ( humidity .AND. q_s <= qv1 ) THEN 1079 1175 r_canopy(j,i) = 0.0_wp 1080 1176 r_soil(j,i) = 0.0_wp … … 1083 1179 ! 1084 1180 !-- Calculate coefficients for the total evapotranspiration 1085 f_qsws_veg = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/ & 1086 (r_a(j,i) + r_canopy(j,i)) 1087 1088 f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) + & 1181 !-- In case of water surface, set vegetation and soil fluxes to zero. 1182 !-- For pavements, only evaporation of liquid water is possible. 1183 IF ( water_surface(j,i) ) THEN 1184 f_qsws_veg = 0.0_wp 1185 f_qsws_soil = 0.0_wp 1186 f_qsws_liq = rho_lv / r_a(j,i) 1187 ELSEIF ( pave_surface (j,i) ) THEN 1188 f_qsws_veg = 0.0_wp 1189 f_qsws_soil = 0.0_wp 1190 f_qsws_liq = rho_lv * c_liq(j,i) / r_a(j,i) 1191 ELSE 1192 f_qsws_veg = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/ & 1193 (r_a(j,i) + r_canopy(j,i)) 1194 f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) + & 1089 1195 r_soil(j,i)) 1090 f_qsws_liq = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i) 1091 1092 1196 f_qsws_liq = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i) 1197 ENDIF 1093 1198 ! 1094 1199 !-- If soil moisture is below wilting point, plants do no longer … … 1098 1203 ! ENDIF 1099 1204 1100 !1101 !-- If dewfall is deactivated, vegetation, soil and liquid water1102 !-- reservoir are not allowed to take up water from the super-saturated1103 !-- air.1104 IF ( humidity .AND. q_s <= qv1 .AND. .NOT. dewfall ) THEN1105 f_qsws_veg = 0.0_wp1106 f_qsws_soil = 0.0_wp1107 f_qsws_liq = 0.0_wp1108 ENDIF1109 1110 1205 f_shf = rho_cp / r_a(j,i) 1111 1206 f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq … … 1123 1218 rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i) 1124 1219 1220 ! 1221 !-- Calculate new skin temperature 1125 1222 IF ( humidity ) THEN 1126 1223 #if defined ( __rrtmg ) … … 1186 1283 !-- Implicit solution when the surface layer has no heat capacity, 1187 1284 !-- otherwise use RK3 scheme. 1188 t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface * & 1189 t_surface(j,i) ) / ( c_surface + coef_2 * dt_3d& 1190 * tsc(2) ) 1285 t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp * & 1286 t_surface(j,i) ) / ( c_surface_tmp + coef_2 & 1287 * dt_3d * tsc(2) ) 1288 1191 1289 ! 1192 1290 !-- Add RK3 term 1193 IF ( c_surface /= 0.0_wp ) THEN1291 IF ( c_surface_tmp /= 0.0_wp ) THEN 1194 1292 1195 1293 t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3) & … … 1211 1309 ENDIF 1212 1310 ENDIF 1213 1214 1311 ENDIF 1215 1312 … … 1223 1320 !-- often reached, when no oscillations would occur (causes immense 1224 1321 !-- computing time for the radiation code). 1225 IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp .AND.&1322 IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp .AND. & 1226 1323 unscheduled_radiation_calls ) THEN 1227 1324 force_radiation_call_l = .TRUE. … … 1249 1346 #endif 1250 1347 1251 1252 1348 ghf_eb(j,i) = lambda_surface * (t_surface_p(j,i) & 1253 1349 - t_soil(nzb_soil,j,i)) … … 1275 1371 * t_surface_p(j,i) ) 1276 1372 ENDIF 1373 1277 1374 ! 1278 1375 !-- Calculate the true surface resistance … … 1280 1377 r_s(j,i) = 1.0E10_wp 1281 1378 ELSE 1282 r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt &1379 r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt & 1283 1380 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) & 1284 1381 / qsws_eb(j,i) - r_a(j,i) … … 1297 1394 ! 1298 1395 !-- Add precipitation to liquid water reservoir, if possible. 1299 !-- Otherwise, add the water to bare soil fraction. 1396 !-- Otherwise, add the water to soil. In case of 1397 !-- pavements, the exceeding water amount is implicitely removed 1398 !-- as runoff as qsws_soil_eb is then not used in the soil model 1300 1399 IF ( m_liq_eb(j,i) /= m_liq_eb_max ) THEN 1301 1400 qsws_liq_eb(j,i) = qsws_liq_eb(j,i) & … … 1376 1475 !-- Make a logical OR for all processes. Force radiation call if at 1377 1476 !-- least one processor reached the threshold change in skin temperature 1378 IF ( unscheduled_radiation_calls .AND. intermediate_timestep_count&1477 IF ( unscheduled_radiation_calls .AND. intermediate_timestep_count & 1379 1478 == intermediate_timestep_count_max-1 ) THEN 1380 1479 #if defined( __parallel ) … … 1388 1487 ENDIF 1389 1488 1489 ! 1490 !-- Calculate surface specific humidity 1491 IF ( humidity ) THEN 1492 CALL calc_q_surface 1493 ENDIF 1494 1495 ! 1496 !-- Calculate new roughness lengths (for water surfaces only) 1497 CALL calc_z0_water_surface 1498 1390 1499 1391 1500 END SUBROUTINE lsm_energy_balance … … 1415 1524 DO i = nxlg, nxrg 1416 1525 DO j = nysg, nyng 1417 DO k = nzb_soil, nzt_soil 1418 ! 1419 !-- Calculate volumetric heat capacity of the soil, taking into 1420 !-- account water content 1421 rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) & 1422 + rho_c_water * m_soil(k,j,i)) 1423 1424 ! 1425 !-- Calculate soil heat conductivity at the center of the soil 1426 !-- layers 1427 lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) * & 1428 lambda_h_water ** m_soil(k,j,i) 1429 1430 ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i))) 1431 1432 lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) + & 1433 lambda_h_dry 1434 1435 ENDDO 1436 1437 ! 1438 !-- Calculate soil heat conductivity (lambda_h) at the _stag level 1439 !-- using linear interpolation 1440 DO k = nzb_soil, nzt_soil-1 1441 lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) * 0.5_wp 1442 ENDDO 1443 lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil) 1444 1445 ! 1446 !-- Prognostic equation for soil temperature t_soil 1447 tend(:) = 0.0_wp 1448 tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) * & 1449 ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i) & 1450 - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1) & 1451 + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil) 1452 1453 1454 1455 1456 DO k = nzb_soil+1, nzt_soil 1457 tend(k) = (1.0_wp/rho_c_total(k,j,i)) & 1458 * ( lambda_h(k,j,i) & 1459 * ( t_soil(k+1,j,i) - t_soil(k,j,i) ) & 1460 * ddz_soil(k+1) & 1461 - lambda_h(k-1,j,i) & 1462 * ( t_soil(k,j,i) - t_soil(k-1,j,i) ) & 1463 * ddz_soil(k) & 1464 ) * ddz_soil_stag(k) 1465 ENDDO 1466 1467 t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i) & 1468 + dt_3d * ( tsc(2) & 1469 * tend(:) + tsc(3) & 1470 * tt_soil_m(:,j,i) ) 1471 1472 ! 1473 !-- Calculate t_soil tendencies for the next Runge-Kutta step 1474 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1475 IF ( intermediate_timestep_count == 1 ) THEN 1476 DO k = nzb_soil, nzt_soil 1477 tt_soil_m(k,j,i) = tend(k) 1478 ENDDO 1479 ELSEIF ( intermediate_timestep_count < & 1480 intermediate_timestep_count_max ) THEN 1481 DO k = nzb_soil, nzt_soil 1482 tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp & 1483 * tt_soil_m(k,j,i) 1484 ENDDO 1485 ENDIF 1486 ENDIF 1487 1488 1489 DO k = nzb_soil, nzt_soil 1490 1491 ! 1492 !-- Calculate soil diffusivity at the center of the soil layers 1493 lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat & 1494 / m_sat(j,i) ) * ( MAX(m_soil(k,j,i), & 1495 m_wilt(j,i)) / m_sat(j,i) )**(b_ch + 2.0_wp) 1496 1497 ! 1498 !-- Parametrization of Van Genuchten 1499 IF ( soil_type /= 7 ) THEN 1500 ! 1501 !-- Calculate the hydraulic conductivity after Van Genuchten 1502 !-- (1980) 1503 h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) - & 1504 MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_vg(j,i) & 1505 / (n_vg(j,i)-1.0_wp)) - 1.0_wp & 1506 )**(1.0_wp/n_vg(j,i)) / alpha_vg(j,i) 1507 1508 1509 gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp + & 1510 (alpha_vg(j,i)*h_vg)**n_vg(j,i))**(1.0_wp & 1511 -1.0_wp/n_vg(j,i)) - (alpha_vg(j,i)*h_vg & 1512 )**(n_vg(j,i)-1.0_wp))**2 ) & 1513 / ( (1.0_wp + (alpha_vg(j,i)*h_vg & 1514 )**n_vg(j,i))**((1.0_wp - 1.0_wp/n_vg(j,i)) & 1515 *(l_vg(j,i) + 2.0_wp)) ) 1516 1517 ! 1518 !-- Parametrization of Clapp & Hornberger 1519 ELSE 1520 gamma_temp(k) = gamma_w_sat(j,i) * (m_soil(k,j,i) & 1521 / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp) 1522 ENDIF 1523 1524 ENDDO 1525 1526 1527 IF ( humidity ) THEN 1528 ! 1529 !-- Calculate soil diffusivity (lambda_w) at the _stag level 1530 !-- using linear interpolation. To do: replace this with 1531 !-- ECMWF-IFS Eq. 8.81 1526 1527 IF ( pave_surface(j,i) ) THEN 1528 rho_c_total(nzb_soil,j,i) = pave_heat_capacity 1529 lambda_temp(nzb_soil) = pave_heat_conductivity 1530 ENDIF 1531 1532 IF ( .NOT. water_surface(j,i) ) THEN 1533 DO k = nzb_soil, nzt_soil 1534 1535 1536 IF ( pave_surface(j,i) .AND. zs(k) <= pave_depth ) THEN 1537 1538 rho_c_total(k,j,i) = pave_heat_capacity 1539 lambda_temp(k) = pave_heat_conductivity 1540 1541 ELSE 1542 ! 1543 !-- Calculate volumetric heat capacity of the soil, taking 1544 !-- into account water content 1545 rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) & 1546 + rho_c_water * m_soil(k,j,i)) 1547 1548 ! 1549 !-- Calculate soil heat conductivity at the center of the soil 1550 !-- layers 1551 lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) * & 1552 lambda_h_water ** m_soil(k,j,i) 1553 1554 ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i))) 1555 1556 lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) + & 1557 lambda_h_dry 1558 ENDIF 1559 1560 ENDDO 1561 1562 ! 1563 !-- Calculate soil heat conductivity (lambda_h) at the _stag level 1564 !-- using linear interpolation. For pavement surface, the 1565 !-- true pavement depth is considered 1532 1566 DO k = nzb_soil, nzt_soil-1 1533 1534 lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) & 1535 * 0.5_wp 1536 gamma_w(k,j,i) = ( gamma_temp(k+1) + gamma_temp(k) ) & 1537 * 0.5_wp 1538 1567 IF ( pave_surface(j,i) .AND. zs(k) < pave_depth & 1568 .AND. zs(k+1) > pave_depth ) THEN 1569 lambda_h(k,j,i) = ( pave_depth - zs(k) ) / dz_soil(k+1) & 1570 * lambda_temp(k) & 1571 + ( 1.0_wp - ( pave_depth - zs(k) ) & 1572 / dz_soil(k+1) ) * lambda_temp(k+1) 1573 ELSE 1574 lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) & 1575 * 0.5_wp 1576 ENDIF 1539 1577 ENDDO 1540 1541 ! 1542 ! 1543 !-- In case of a closed bottom (= water content is conserved), set 1544 !-- hydraulic conductivity to zero to that no water will be lost 1545 !-- in the bottom layer. 1546 IF ( conserve_water_content ) THEN 1547 gamma_w(nzt_soil,j,i) = 0.0_wp 1548 ELSE 1549 gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil) 1550 ENDIF 1551 1552 !-- The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v)) 1553 !-- ensures the mass conservation for water. The transpiration of 1554 !-- plants equals the cumulative withdrawals by the roots in the 1555 !-- soil. The scheme takes into account the availability of water 1556 !-- in the soil layers as well as the root fraction in the 1557 !-- respective layer. Layer with moisture below wilting point will 1558 !-- not contribute, which reflects the preference of plants to 1559 !-- take water from moister layers. 1560 1561 ! 1562 !-- Calculate the root extraction (ECMWF 7.69, the sum of root_extr 1563 !-- = 1). The energy balance solver guarantees a positive 1564 !-- transpiration, so that there is no need for an additional check. 1565 DO k = nzb_soil, nzt_soil 1566 IF ( m_soil(k,j,i) > m_wilt(j,i) ) THEN 1567 m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i) 1568 ENDIF 1569 ENDDO 1570 1571 IF ( m_total > 0.0_wp ) THEN 1572 DO k = nzb_soil, nzt_soil 1573 IF ( m_soil(k,j,i) > m_wilt(j,i) ) THEN 1574 root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total 1575 ELSE 1576 root_extr(k) = 0.0_wp 1577 ENDIF 1578 ENDDO 1579 ENDIF 1580 1581 ! 1582 !-- Prognostic equation for soil water content m_soil 1578 lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil) 1579 1580 1581 1582 1583 ! 1584 !-- Prognostic equation for soil temperature t_soil 1583 1585 tend(:) = 0.0_wp 1584 tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * ( & 1585 m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )&1586 * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - (&1587 root_extr(nzb_soil) * qsws_veg_eb(j,i)&1588 + qsws_soil_eb(j,i) ) * drho_l_lv ) &1589 * ddz_soil_stag(nzb_soil) 1590 1591 DO k = nzb_soil+1, nzt_soil-11592 tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)&1593 - m_soil(k,j,i) ) * ddz_soil(k+1) - gamma_w(k,j,i)&1594 - lambda_w(k-1,j,i) * (m_soil(k,j,i) -&1595 m_soil(k-1,j,i)) * ddz_soil(k)&1596 + gamma_w(k-1,j,i) - (root_extr(k)&1597 * qsws_veg_eb(j,i) * drho_l_lv)&1598 ) * ddz_soil_stag(k)1586 1587 tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) * & 1588 ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i) & 1589 - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1) & 1590 + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil) 1591 1592 DO k = nzb_soil+1, nzt_soil 1593 tend(k) = (1.0_wp/rho_c_total(k,j,i)) & 1594 * ( lambda_h(k,j,i) & 1595 * ( t_soil(k+1,j,i) - t_soil(k,j,i) ) & 1596 * ddz_soil(k+1) & 1597 - lambda_h(k-1,j,i) & 1598 * ( t_soil(k,j,i) - t_soil(k-1,j,i) ) & 1599 * ddz_soil(k) & 1600 ) * ddz_soil_stag(k) 1599 1601 1600 1602 ENDDO 1601 tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i) & 1602 - lambda_w(nzt_soil-1,j,i) & 1603 * (m_soil(nzt_soil,j,i) & 1604 - m_soil(nzt_soil-1,j,i)) & 1605 * ddz_soil(nzt_soil) & 1606 + gamma_w(nzt_soil-1,j,i) - ( & 1607 root_extr(nzt_soil) & 1608 * qsws_veg_eb(j,i) * drho_l_lv ) & 1609 ) * ddz_soil_stag(nzt_soil) 1610 1611 m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)& 1612 + dt_3d * ( tsc(2) * tend(:) & 1613 + tsc(3) * tm_soil_m(:,j,i) ) 1614 1615 ! 1616 !-- Account for dry soils (find a better solution here!) 1617 DO k = nzb_soil, nzt_soil 1618 IF ( m_soil_p(k,j,i) < 0.0_wp ) m_soil_p(k,j,i) = 0.0_wp 1619 ENDDO 1620 1621 1622 1623 1624 ! 1625 !-- Calculate m_soil tendencies for the next Runge-Kutta step 1603 1604 t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)& 1605 + dt_3d * ( tsc(2) & 1606 * tend(nzb_soil:nzt_soil) & 1607 + tsc(3) & 1608 * tt_soil_m(:,j,i) ) 1609 1610 ! 1611 !-- Calculate t_soil tendencies for the next Runge-Kutta step 1626 1612 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1627 1613 IF ( intermediate_timestep_count == 1 ) THEN 1628 1614 DO k = nzb_soil, nzt_soil 1629 t m_soil_m(k,j,i) = tend(k)1615 tt_soil_m(k,j,i) = tend(k) 1630 1616 ENDDO 1631 1617 ELSEIF ( intermediate_timestep_count < & 1632 1618 intermediate_timestep_count_max ) THEN 1633 1619 DO k = nzb_soil, nzt_soil 1634 t m_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp &1635 * tm_soil_m(k,j,i)1620 tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp & 1621 * tt_soil_m(k,j,i) 1636 1622 ENDDO 1637 1623 ENDIF 1638 1624 ENDIF 1639 1625 1626 1627 DO k = nzb_soil, nzt_soil 1628 1629 ! 1630 !-- Calculate soil diffusivity at the center of the soil layers 1631 lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat & 1632 / m_sat(j,i) ) * ( MAX( m_soil(k,j,i), & 1633 m_wilt(j,i) ) / m_sat(j,i) )**( & 1634 b_ch + 2.0_wp ) 1635 1636 ! 1637 !-- Parametrization of Van Genuchten 1638 IF ( soil_type /= 7 ) THEN 1639 ! 1640 !-- Calculate the hydraulic conductivity after Van Genuchten 1641 !-- (1980) 1642 h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) - & 1643 MAX( m_soil(k,j,i), m_wilt(j,i) ) ) )**( & 1644 n_vg(j,i) / (n_vg(j,i) - 1.0_wp ) ) - 1.0_wp & 1645 )**( 1.0_wp / n_vg(j,i) ) / alpha_vg(j,i) 1646 1647 1648 gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp + & 1649 ( alpha_vg(j,i) * h_vg )**n_vg(j,i))**( & 1650 1.0_wp - 1.0_wp / n_vg(j,i) ) - ( & 1651 alpha_vg(j,i) * h_vg )**( n_vg(j,i) & 1652 - 1.0_wp) )**2 ) & 1653 / ( ( 1.0_wp + ( alpha_vg(j,i) * h_vg & 1654 )**n_vg(j,i) )**( ( 1.0_wp - 1.0_wp & 1655 / n_vg(j,i) ) *( l_vg(j,i) + 2.0_wp) ) ) 1656 1657 ! 1658 !-- Parametrization of Clapp & Hornberger 1659 ELSE 1660 gamma_temp(k) = gamma_w_sat(j,i) * ( m_soil(k,j,i) & 1661 / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp) 1662 ENDIF 1663 1664 ENDDO 1665 1666 ! 1667 !-- Prognostic equation for soil moisture content. Only performed, 1668 !-- when humidity is enabled in the atmosphere and the surface type 1669 !-- is not pavement (implies dry soil below). 1670 IF ( humidity .AND. .NOT. pave_surface(j,i) ) THEN 1671 ! 1672 !-- Calculate soil diffusivity (lambda_w) at the _stag level 1673 !-- using linear interpolation. To do: replace this with 1674 !-- ECMWF-IFS Eq. 8.81 1675 DO k = nzb_soil, nzt_soil-1 1676 1677 lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) & 1678 * 0.5_wp 1679 gamma_w(k,j,i) = ( gamma_temp(k+1) + gamma_temp(k) ) & 1680 * 0.5_wp 1681 1682 ENDDO 1683 1684 ! 1685 ! 1686 !-- In case of a closed bottom (= water content is conserved), set 1687 !-- hydraulic conductivity to zero to that no water will be lost 1688 !-- in the bottom layer. 1689 IF ( conserve_water_content ) THEN 1690 gamma_w(nzt_soil,j,i) = 0.0_wp 1691 ELSE 1692 gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil) 1693 ENDIF 1694 1695 !-- The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v)) 1696 !-- ensures the mass conservation for water. The transpiration of 1697 !-- plants equals the cumulative withdrawals by the roots in the 1698 !-- soil. The scheme takes into account the availability of water 1699 !-- in the soil layers as well as the root fraction in the 1700 !-- respective layer. Layer with moisture below wilting point will 1701 !-- not contribute, which reflects the preference of plants to 1702 !-- take water from moister layers. 1703 1704 ! 1705 !-- Calculate the root extraction (ECMWF 7.69, the sum of root_extr 1706 !-- = 1). The energy balance solver guarantees a positive 1707 !-- transpiration, so that there is no need for an additional check. 1708 DO k = nzb_soil, nzt_soil 1709 IF ( m_soil(k,j,i) > m_wilt(j,i) ) THEN 1710 m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i) 1711 ENDIF 1712 ENDDO 1713 1714 IF ( m_total > 0.0_wp ) THEN 1715 DO k = nzb_soil, nzt_soil 1716 IF ( m_soil(k,j,i) > m_wilt(j,i) ) THEN 1717 root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total 1718 ELSE 1719 root_extr(k) = 0.0_wp 1720 ENDIF 1721 ENDDO 1722 ENDIF 1723 1724 ! 1725 !-- Prognostic equation for soil water content m_soil. 1726 tend(:) = 0.0_wp 1727 1728 tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * ( & 1729 m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) ) & 1730 * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( & 1731 root_extr(nzb_soil) * qsws_veg_eb(j,i) & 1732 + qsws_soil_eb(j,i) ) * drho_l_lv ) & 1733 * ddz_soil_stag(nzb_soil) 1734 1735 DO k = nzb_soil+1, nzt_soil-1 1736 tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i) & 1737 - m_soil(k,j,i) ) * ddz_soil(k+1) & 1738 - gamma_w(k,j,i) & 1739 - lambda_w(k-1,j,i) * (m_soil(k,j,i) - & 1740 m_soil(k-1,j,i)) * ddz_soil(k) & 1741 + gamma_w(k-1,j,i) - (root_extr(k) & 1742 * qsws_veg_eb(j,i) * drho_l_lv) & 1743 ) * ddz_soil_stag(k) 1744 1745 ENDDO 1746 tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i) & 1747 - lambda_w(nzt_soil-1,j,i) & 1748 * (m_soil(nzt_soil,j,i) & 1749 - m_soil(nzt_soil-1,j,i)) & 1750 * ddz_soil(nzt_soil) & 1751 + gamma_w(nzt_soil-1,j,i) - ( & 1752 root_extr(nzt_soil) & 1753 * qsws_veg_eb(j,i) * drho_l_lv ) & 1754 ) * ddz_soil_stag(nzt_soil) 1755 1756 m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)& 1757 + dt_3d * ( tsc(2) * tend(:) & 1758 + tsc(3) * tm_soil_m(:,j,i) ) 1759 1760 ! 1761 !-- Account for dry soils (find a better solution here!) 1762 DO k = nzb_soil, nzt_soil 1763 IF ( m_soil_p(k,j,i) < 0.0_wp ) m_soil_p(k,j,i) = 0.0_wp 1764 ENDDO 1765 1766 ! 1767 !-- Calculate m_soil tendencies for the next Runge-Kutta step 1768 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1769 IF ( intermediate_timestep_count == 1 ) THEN 1770 DO k = nzb_soil, nzt_soil 1771 tm_soil_m(k,j,i) = tend(k) 1772 ENDDO 1773 ELSEIF ( intermediate_timestep_count < & 1774 intermediate_timestep_count_max ) THEN 1775 DO k = nzb_soil, nzt_soil 1776 tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp& 1777 * tm_soil_m(k,j,i) 1778 ENDDO 1779 ENDIF 1780 ENDIF 1781 1782 ENDIF 1783 1640 1784 ENDIF 1641 1785 1642 1786 ENDDO 1643 1787 ENDDO 1644 1645 !1646 !-- Calculate surface specific humidity1647 IF ( humidity ) THEN1648 CALL calc_q_surface1649 ENDIF1650 1651 1788 1652 1789 END SUBROUTINE lsm_soil_model … … 1656 1793 ! Description: 1657 1794 ! ------------ 1658 !> Calculation of specific humidity of the surface layer (surface) 1795 !> Calculation of roughness length for open water (lakes, ocean). The 1796 !> parameterization follows Charnock (1955). Two different implementations 1797 !> are available: as in ECMWF-IFS (Beljaars 1994) or as in FLake (Subin et al. 1798 !> 2012) 1799 !------------------------------------------------------------------------------! 1800 SUBROUTINE calc_z0_water_surface 1801 1802 USE control_parameters, & 1803 ONLY: g, kappa, molecular_viscosity 1804 1805 IMPLICIT NONE 1806 1807 INTEGER :: i !< running index 1808 INTEGER :: j !< running index 1809 1810 REAL(wp), PARAMETER :: alpha_ch = 0.018_wp !< Charnock constant (0.01-0.11). Use 0.01 for FLake and 0.018 for ECMWF 1811 ! REAL(wp), PARAMETER :: pr_number = 0.71_wp !< molecular Prandtl number in the Charnock parameterization (differs from prandtl_number) 1812 ! REAL(wp), PARAMETER :: sc_number = 0.66_wp !< molecular Schmidt number in the Charnock parameterization 1813 ! REAL(wp) :: re_0 !< near-surface roughness Reynolds number 1814 1815 1816 DO i = nxlg, nxrg 1817 DO j = nysg, nyng 1818 IF ( water_surface(j,i) ) THEN 1819 1820 ! 1821 !-- Disabled: FLake parameterization. Ideally, the Charnock 1822 !-- coefficient should depend on the water depth and the fetch 1823 !-- length 1824 ! re_0 = z0(j,i) * us(j,i) / molecular_viscosity 1825 ! 1826 ! z0(j,i) = MAX( 0.1_wp * molecular_viscosity / us(j,i), & 1827 ! alpha_ch * us(j,i) / g ) 1828 ! 1829 ! z0h(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 3.2_wp ) ) 1830 ! z0q(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 4.2_wp ) ) 1831 1832 ! 1833 !-- Set minimum roughness length for u* > 0.2 1834 ! IF ( us(j,i) > 0.2_wp ) THEN 1835 ! z0h(j,i) = MAX( 1.0E-5_wp, z0h(j,i) ) 1836 ! z0q(j,i) = MAX( 1.0E-5_wp, z0q(j,i) ) 1837 ! ENDIF 1838 1839 ! 1840 !-- ECMWF IFS model parameterization after Beljaars (1994). At low 1841 !-- wind speed, the sea surface becomes aerodynamically smooth and 1842 !-- the roughness scales with the viscosity. At high wind speed, the 1843 !-- Charnock relation is used. 1844 z0(j,i) = ( 0.11_wp * molecular_viscosity / us(j,i) ) & 1845 + ( alpha_ch * us(j,i)**2 / g ) 1846 1847 z0h(j,i) = 0.40_wp * molecular_viscosity / us(j,i) 1848 z0q(j,i) = 0.62_wp * molecular_viscosity / us(j,i) 1849 1850 ENDIF 1851 ENDDO 1852 ENDDO 1853 1854 END SUBROUTINE calc_z0_water_surface 1855 1856 1857 !------------------------------------------------------------------------------! 1858 ! Description: 1859 ! ------------ 1860 !> Calculation of specific humidity of the skin layer (surface). It is assumend 1861 !> that the skin is always saturated. 1659 1862 !------------------------------------------------------------------------------! 1660 1863 SUBROUTINE calc_q_surface … … 1665 1868 INTEGER :: j !< running index 1666 1869 INTEGER :: k !< running index 1870 1667 1871 REAL(wp) :: resistance !< aerodynamic and soil resistance term 1668 1872 … … 1701 1905 END SUBROUTINE calc_q_surface 1702 1906 1907 1703 1908 !------------------------------------------------------------------------------! 1704 1909 ! Description: -
palm/trunk/SOURCE/modules.f90
r1787 r1788 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Added roughness length for moisture (z0q) 22 22 ! 23 23 ! Former revisions: … … 384 384 saswsb, saswst, shf, tnudge, td_lsa_lpt, td_lsa_q, td_sub_lpt, & 385 385 td_sub_q, total_2d_a, total_2d_o, ts, tswst, ug_vert, unudge, us, & 386 usws, uswst, vnudge, vg_vert, vsws, vswst, wnudge, wsubs_vert, z0, z0h 386 usws, uswst, vnudge, vg_vert, vsws, vswst, wnudge, wsubs_vert, & 387 z0, & !< roughness length for momentum 388 z0h, & !< roughness length for heat 389 z0q !< roughness length for moisture 387 390 388 391 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & … … 432 435 USE kinds 433 436 434 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lwp_av, & !> Avg. liquid water path 435 precipitation_rate_av, & !> Avg. of precipitation rate 436 ol_av, & !> Avg. of Obukhov length 437 qsws_av, & !> Avg. of surface moisture flux 438 shf_av, & !> Avg. of surface heat flux 439 ts_av, & !> Avg. of characteristic temperature scale 440 us_av, & !> Avg. of friction velocity 441 z0_av, & !> Avg. of roughness length for momentum 442 z0h_av !> Avg. of roughness length for heat and moisture 437 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lwp_av, & !< Avg. liquid water path 438 precipitation_rate_av, & !< Avg. of precipitation rate 439 ol_av, & !< Avg. of Obukhov length 440 qsws_av, & !< Avg. of surface moisture flux 441 shf_av, & !< Avg. of surface heat flux 442 ts_av, & !< Avg. of characteristic temperature scale 443 us_av, & !< Avg. of friction velocity 444 z0_av, & !< Avg. of roughness length for momentum 445 z0h_av, & !< Avg. of roughness length for heat 446 z0q_av !< Avg. of roughness length for moisture 443 447 444 448 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & -
palm/trunk/SOURCE/package_parin.f90
r1787 r1788 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Parameter dewfall removed. 22 22 ! 23 23 ! Former revisions: … … 129 129 USE land_surface_model_mod, & 130 130 ONLY: alpha_vangenuchten, c_surface, canopy_resistance_coefficient, & 131 conserve_water_content, dewfall, field_capacity,&131 conserve_water_content, field_capacity, & 132 132 f_shortwave_incoming, hydraulic_conductivity, & 133 133 lambda_surface_stable, lambda_surface_unstable, leaf_area_index, & 134 134 land_surface, l_vangenuchten, min_canopy_resistance, & 135 min_soil_resistance, n_vangenuchten, residual_moisture, & 135 min_soil_resistance, n_vangenuchten, pave_heat_capacity, & 136 pave_depth, pave_heat_conductivity, residual_moisture, & 136 137 root_fraction, saturation_moisture, wilting_point, & 137 138 skip_time_do_lsm, soil_moisture, soil_temperature, soil_type, & 138 vegetation_coverage, veg_type, zs, z0_eb, z0h_eb 139 vegetation_coverage, veg_type, zs, z0_eb, z0h_eb, z0q_eb 139 140 140 141 USE kinds … … 205 206 NAMELIST /lsm_par/ alpha_vangenuchten, c_surface, & 206 207 canopy_resistance_coefficient, & 207 conserve_water_content, dewfall,&208 conserve_water_content, & 208 209 f_shortwave_incoming, field_capacity, & 209 210 hydraulic_conductivity, & … … 212 213 l_vangenuchten, min_canopy_resistance, & 213 214 min_soil_resistance, n_vangenuchten, & 215 pave_depth, pave_heat_capacity, & 216 pave_heat_conductivity, & 214 217 residual_moisture, root_fraction, & 215 218 saturation_moisture, skip_time_do_lsm, & 216 219 soil_moisture, soil_temperature, soil_type, & 217 220 vegetation_coverage, veg_type, wilting_point,& 218 zs, z0_eb, z0h_eb 221 zs, z0_eb, z0h_eb, z0q_eb 219 222 220 223 -
palm/trunk/SOURCE/radiation_model.f90
r1784 r1788 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Added new albedo class for pavements / roads. 22 22 ! 23 23 ! Former revisions: … … 112 112 113 113 #if defined ( __rrtmg ) 114 115 ! USE netcdf_interface, & 116 ! ONLY: nc_stat, netcdf_handle_error 117 114 118 USE parrrsw, & 115 119 ONLY: naerec, nbndsw … … 139 143 ! 140 144 !-- Predefined Land surface classes (albedo_type) after Briegleb (1992) 141 CHARACTER(37), DIMENSION(0:1 6), PARAMETER :: albedo_type_name = (/ &145 CHARACTER(37), DIMENSION(0:17), PARAMETER :: albedo_type_name = (/ & 142 146 'user defined ', & ! 0 143 147 'ocean ', & ! 1 … … 156 160 'land ice ', & ! 14 157 161 'sea ice ', & ! 15 158 'snow ' & ! 16 162 'snow ', & ! 16 163 'pavement/roads ' & ! 17 159 164 /) 160 165 … … 163 168 day_init = 172, & !< day of the year at model start (21/06) 164 169 dots_rad = 0 !< starting index for timeseries output 165 166 167 168 169 170 170 171 171 LOGICAL :: unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed … … 214 214 !-- Land surface albedos for solar zenith angle of 60° after Briegleb (1992) 215 215 !-- (shortwave, longwave, broadband): sw, lw, bb, 216 REAL(wp), DIMENSION(0:2,1:1 6), PARAMETER :: albedo_pars = RESHAPE( (/&216 REAL(wp), DIMENSION(0:2,1:17), PARAMETER :: albedo_pars = RESHAPE( (/& 217 217 0.06_wp, 0.06_wp, 0.06_wp, & ! 1 218 218 0.09_wp, 0.28_wp, 0.19_wp, & ! 2 … … 230 230 0.90_wp, 0.65_wp, 0.77_wp, & ! 14 231 231 0.90_wp, 0.65_wp, 0.77_wp, & ! 15 232 0.95_wp, 0.70_wp, 0.82_wp & ! 16 233 /), (/ 3, 16 /) ) 232 0.95_wp, 0.70_wp, 0.82_wp, & ! 16 233 0.08_wp, 0.08_wp, 0.08_wp & ! 17 234 /), (/ 3, 17 /) ) 234 235 235 236 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & … … 275 276 rrtm_iaer = 0, & !< aerosol option flag (0: no aerosol layers, for lw only: 6 (requires setting of rrtm_sw_ecaer), 10: one or more aerosol layers (not implemented) 276 277 rrtm_idrv = 1 !< longwave upward flux calculation option (0,1) 278 279 INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling 277 280 278 281 LOGICAL :: snd_exists = .FALSE. !< flag parameter to check whether a user-defined input files exists … … 1030 1033 rrtm_aldir(0,:,:) = aldif 1031 1034 rrtm_asdir(0,:,:) = asdif 1035 1036 ! 1037 !-- Asphalt 1038 ELSEIF ( albedo_type == 17 ) THEN 1039 rrtm_aldir(0,:,:) = aldif 1040 rrtm_asdir(0,:,:) = asdif 1032 1041 ! 1033 1042 !-- Land surfaces -
palm/trunk/SOURCE/read_3d_binary.f90
r1758 r1788 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Added z0q and z0q_av 22 22 ! 23 23 ! Former revisions: … … 97 97 rif_wall, shf, ts, tswst, u, u_m_l, u_m_n, u_m_r, u_m_s, us, & 98 98 usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s, vpt, vsws, vswst, & 99 w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h 99 w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h, z0q 100 100 101 101 USE averaging … … 354 354 !-- First compare the version numbers 355 355 READ ( 13 ) version_on_file 356 binary_version = '4. 2'356 binary_version = '4.3' 357 357 IF ( TRIM( version_on_file ) /= TRIM( binary_version ) ) THEN 358 358 WRITE( message_string, * ) 'version mismatch concerning data ', & … … 1358 1358 IF ( k == 1 ) READ ( 13 ) tmp_2d 1359 1359 z0h_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 1360 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 1361 1362 CASE ( 'z0q' ) 1363 IF ( k == 1 ) READ ( 13 ) tmp_2d 1364 z0q(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 1365 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 1366 1367 CASE ( 'z0q_av' ) 1368 IF ( .NOT. ALLOCATED( z0q_av ) ) THEN 1369 ALLOCATE( z0q_av(nysg:nyng,nxlg:nxrg) ) 1370 ENDIF 1371 IF ( k == 1 ) READ ( 13 ) tmp_2d 1372 z0q_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 1360 1373 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 1361 1374 -
palm/trunk/SOURCE/sum_up_3d_data.f90
r1694 r1788 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Added z0q and z0q_av 22 22 ! 23 23 ! Former revisions: … … 92 92 USE arrays_3d, & 93 93 ONLY: dzw, e, nr, ol, p, pt, q, qc, ql, ql_c, ql_v, qr, qsws, rho, sa,& 94 shf, ts, u, us, v, vpt, w, z0, z0h 94 shf, ts, u, us, v, vpt, w, z0, z0h, z0q 95 95 96 96 USE averaging, & … … 98 98 precipitation_rate_av, pt_av, q_av, qc_av, ql_av, ql_c_av, & 99 99 ql_v_av, ql_vp_av, qr_av, qsws_av, qv_av, rho_av, s_av, sa_av, & 100 shf_av, ts_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, z0h_av 100 shf_av, ts_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, z0h_av, & 101 z0q_av 101 102 102 103 USE cloud_parameters, & … … 494 495 ENDIF 495 496 z0h_av = 0.0_wp 497 498 CASE ( 'z0q*' ) 499 IF ( .NOT. ALLOCATED( z0q_av ) ) THEN 500 ALLOCATE( z0q_av(nysg:nyng,nxlg:nxrg) ) 501 ENDIF 502 z0q_av = 0.0_wp 496 503 497 504 CASE DEFAULT … … 1006 1013 ENDDO 1007 1014 1015 CASE ( 'z0q*' ) 1016 DO i = nxlg, nxrg 1017 DO j = nysg, nyng 1018 z0q_av(j,i) = z0q_av(j,i) + z0q(j,i) 1019 ENDDO 1020 ENDDO 1021 1008 1022 CASE DEFAULT 1009 1023 ! -
palm/trunk/SOURCE/surface_layer_fluxes.f90
r1758 r1788 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Added parameter z0q which replaces z0h in the similarity functions for 22 ! humidity. 23 ! Syntax layout improved. 22 24 ! 23 25 ! Former revisions: … … 117 119 USE arrays_3d, & 118 120 ONLY: e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, & 119 shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h 121 shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h, z0q 120 122 121 123 USE cloud_parameters, & … … 234 236 !-- Use either Newton iteration or a lookup table for the bulk Richardson 235 237 !-- number to calculate the Obukhov length 236 ELSEIF ( most_method == 'newton' .OR.most_method == 'lookup' ) THEN238 ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' ) THEN 237 239 238 240 CALL calc_uv_total … … 302 304 !-- Check for roughness heterogeneity. In that case terminate run and 303 305 !-- inform user 304 IF ( MINVAL( z0h ) /= MAXVAL( z0h ) .OR.&306 IF ( MINVAL( z0h ) /= MAXVAL( z0h ) .OR. & 305 307 MINVAL( z0 ) /= MAXVAL( z0 ) ) THEN 306 308 terminate_run_l = .TRUE. … … 347 349 ! 348 350 !-- Calculate stretching factor for the free convection range 349 DO WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )351 DO WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp ) 350 352 regr_old = regr 351 353 regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr ) & … … 533 535 !-- Ensure that the bulk Richardson number and the Obukhov 534 536 !-- lengtH have the same sign 535 IF ( rib(j,i) * ol(j,i) < 0.0_wp .OR.&537 IF ( rib(j,i) * ol(j,i) < 0.0_wp .OR. & 536 538 ABS( ol(j,i) ) == ol_max ) THEN 537 539 IF ( rib(j,i) > 0.0_wp ) ol(j,i) = 0.01_wp … … 638 640 639 641 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 640 !# WARNING: does not work on GPU so far because of DO WHILE construct642 !# WARNING: does not work on GPU so far because of DO WHILE construct 641 643 !!!!!!$acc kernels loop 642 644 DO i = nxl, nxr … … 656 658 l = l_bnd 657 659 IF ( rib_tab(l) - rib(j,i) > 0.0_wp ) THEN 658 DO WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp .AND.l > 0 )660 DO WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp .AND. l > 0 ) 659 661 l = l-1 660 662 ENDDO 661 663 ELSE 662 DO WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp&663 .AND. l < num_steps-1 )664 DO WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp & 665 .AND. l < num_steps-1 ) 664 666 l = l+1 665 667 ENDDO … … 669 671 ! 670 672 !-- Linear interpolation to find the correct value of z/L 671 ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) ) &672 / ( rib_tab(l) - rib_tab(l-1) ) &673 ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) ) & 674 / ( rib_tab(l) - rib_tab(l-1) ) & 673 675 * ( rib(j,i) - rib_tab(l-1) ) ) 674 676 … … 797 799 ! 798 800 !-- For a given surface temperature: 799 IF ( large_scale_forcing .AND.lsf_surf ) THEN801 IF ( large_scale_forcing .AND. lsf_surf ) THEN 800 802 !$OMP PARALLEL DO 801 803 !$acc kernels loop private( j, k ) … … 837 839 IF ( constant_waterflux ) THEN 838 840 ! 839 !-- For a given water flux in the Prandtl layer:841 !-- For a given water flux in the surface layer 840 842 !$OMP PARALLEL DO 841 843 !$acc kernels loop private( j ) … … 847 849 848 850 ELSE 849 coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND.&851 coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. & 850 852 run_coupled ) 851 853 852 IF ( large_scale_forcing .AND.lsf_surf ) THEN854 IF ( large_scale_forcing .AND. lsf_surf ) THEN 853 855 !$OMP PARALLEL DO 854 856 !$acc kernels loop private( j, k ) … … 862 864 863 865 !$OMP PARALLEL DO PRIVATE( e_s, k, z_mo ) 864 !$acc kernels loop independent present( nzb_s_inner, ol, pt, q, qs, qv1, zu, zw, z0 h) private( e_s, j, k, z_mo )866 !$acc kernels loop independent present( nzb_s_inner, ol, pt, q, qs, qv1, zu, zw, z0q ) private( e_s, j, k, z_mo ) 865 867 DO i = nxlg, nxrg 866 868 !$acc loop independent … … 882 884 IF ( cloud_physics ) THEN 883 885 qs(j,i) = kappa * ( qv1(j,i) - q(k,j,i) ) & 884 / ( LOG( z_mo / z0 h(j,i) ) &886 / ( LOG( z_mo / z0q(j,i) ) & 885 887 - psi_h( z_mo / ol(j,i) ) & 886 + psi_h( z0 h(j,i) / ol(j,i) ) )888 + psi_h( z0q(j,i) / ol(j,i) ) ) 887 889 888 890 ELSE 889 891 qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) & 890 / ( LOG( z_mo / z0 h(j,i) ) &892 / ( LOG( z_mo / z0q(j,i) ) & 891 893 - psi_h( z_mo / ol(j,i) ) & 892 + psi_h( z0 h(j,i) / ol(j,i) ) )894 + psi_h( z0q(j,i) / ol(j,i) ) ) 893 895 ENDIF 894 896 … … 901 903 ! 902 904 !-- If required compute qr* and nr* 903 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation ) THEN 905 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation ) & 906 THEN 904 907 905 908 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 906 !$acc kernels loop independent present( nr, nrs, nzb_s_inner, ol, qr, qrs, zu, zw, z0 h) private( j, k, z_mo )909 !$acc kernels loop independent present( nr, nrs, nzb_s_inner, ol, qr, qrs, zu, zw, z0q ) private( j, k, z_mo ) 907 910 DO i = nxlg, nxrg 908 911 !$acc loop independent … … 913 916 914 917 qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) & 915 / ( LOG( z_mo / z0 h(j,i) ) &918 / ( LOG( z_mo / z0q(j,i) ) & 916 919 - psi_h( z_mo / ol(j,i) ) & 917 + psi_h( z0 h(j,i) / ol(j,i) ) )920 + psi_h( z0q(j,i) / ol(j,i) ) ) 918 921 919 922 nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) & 920 / ( LOG( z_mo / z0 h(j,i) ) &923 / ( LOG( z_mo / z0q(j,i) ) & 921 924 - psi_h( z_mo / ol(j,i) ) & 922 + psi_h( z0 h(j,i) / ol(j,i) ) )925 + psi_h( z0q(j,i) / ol(j,i) ) ) 923 926 ENDDO 924 927 ENDDO … … 1003 1006 ! 1004 1007 !-- Compute the vertical kinematic heat flux 1005 IF ( .NOT. constant_heatflux .AND. ( simulated_time <=&1006 skip_time_do_lsm .OR. .NOT.land_surface ) ) THEN1008 IF ( .NOT. constant_heatflux .AND. ( simulated_time <= & 1009 skip_time_do_lsm .OR. .NOT. land_surface ) ) THEN 1007 1010 !$OMP PARALLEL DO 1008 1011 !$acc kernels loop independent present( shf, ts, us ) … … 1018 1021 ! 1019 1022 !-- Compute the vertical water/scalar flux 1020 IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar )&1021 .AND. ( simulated_time <= skip_time_do_lsm .OR. .NOT.&1022 land_surface ) ) THEN1023 IF ( .NOT. constant_waterflux .AND. ( humidity .OR. & 1024 passive_scalar ) .AND. ( simulated_time <= skip_time_do_lsm & 1025 .OR. .NOT. land_surface ) ) THEN 1023 1026 !$OMP PARALLEL DO 1024 1027 !$acc kernels loop independent present( qs, qsws, us ) … … 1089 1092 1090 1093 IF ( zeta < 0.0_wp ) THEN 1091 x = SQRT( SQRT( 1.0_wp - 16.0_wp * zeta ) )1094 x = SQRT( SQRT( 1.0_wp - 16.0_wp * zeta ) ) 1092 1095 psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2 & 1093 1096 * ( 1.0_wp + x**2 ) * 0.125_wp ) … … 1126 1129 1127 1130 IF ( zeta < 0.0_wp ) THEN 1128 x = SQRT( 1.0_wp - 16.0_wp * zeta )1131 x = SQRT( 1.0_wp - 16.0_wp * zeta ) 1129 1132 psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp ) 1130 1133 ELSE -
palm/trunk/SOURCE/write_3d_binary.f90
r1710 r1788 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Added z0q and z0q_av 22 22 ! 23 23 ! Former revisions: … … 90 90 rif_wall, shf, ts, tswst, u, u_m_l, u_m_n, u_m_r, u_m_s, us, & 91 91 usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s, vpt, vsws, vswst, & 92 w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h 92 w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h, z0q 93 93 94 94 USE averaging … … 142 142 ! 143 143 !-- Write arrays. 144 binary_version = '4. 2'144 binary_version = '4.3' 145 145 146 146 WRITE ( 14 ) binary_version … … 470 470 WRITE ( 14 ) 'z0h_av '; WRITE ( 14 ) z0h_av 471 471 ENDIF 472 WRITE ( 14 ) 'z0q '; WRITE ( 14 ) z0q 473 IF ( ALLOCATED( z0q_av ) ) THEN 474 WRITE ( 14 ) 'z0q_av '; WRITE ( 14 ) z0q_av 475 ENDIF 472 476 473 477 !
Note: See TracChangeset
for help on using the changeset viewer.
