Changeset 1976 for palm/trunk/SOURCE/land_surface_model_mod.f90
- Timestamp:
- Jul 27, 2016 1:28:04 PM (8 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/land_surface_model_mod.f90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/land_surface_model_mod.f90
r1973 r1976 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Parts of the code have been reformatted. Use of radiation model output is 22 ! generalized and simplified. Added more output quantities due to modularization 22 23 ! 23 24 ! Former revisions: … … 178 179 179 180 USE radiation_model_mod, & 180 ONLY: force_radiation_call, radiation_scheme, rad_net, rad_sw_in, & 181 rad_lw_out, rad_lw_out_change_0, sigma_sb, & 182 unscheduled_radiation_calls 181 ONLY: force_radiation_call, rad_net, rad_sw_in, rad_lw_out, & 182 rad_lw_out_change_0, unscheduled_radiation_calls 183 183 184 #if defined ( __rrtmg )185 USE radiation_model_mod, &186 ONLY: rrtm_idrv187 #endif188 189 184 USE statistics, & 190 185 ONLY: hom, statistic_regions … … 636 631 END INTERFACE lsm_read_restart_data 637 632 638 639 633 INTERFACE lsm_last_actions 640 634 MODULE PROCEDURE lsm_last_actions … … 648 642 !> Check data output for land surface model 649 643 !------------------------------------------------------------------------------! 650 644 SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k ) 651 645 652 646 653 USE control_parameters, & 654 ONLY: data_output, message_string 655 656 IMPLICIT NONE 657 658 CHARACTER (LEN=*) :: unit !< 659 CHARACTER (LEN=*) :: var !< 660 661 INTEGER(iwp) :: i 662 INTEGER(iwp) :: ilen 663 INTEGER(iwp) :: k 664 665 SELECT CASE ( TRIM( var ) ) 666 667 CASE ( 'm_soil' ) 668 IF ( .NOT. land_surface ) THEN 669 message_string = 'output of "' // TRIM( var ) // '" requi' // & 670 'res land_surface = .TRUE.' 671 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 672 ENDIF 673 unit = 'm3/m3' 647 USE control_parameters, & 648 ONLY: data_output, message_string 649 650 IMPLICIT NONE 651 652 CHARACTER (LEN=*) :: unit !< 653 CHARACTER (LEN=*) :: var !< 654 655 INTEGER(iwp) :: i 656 INTEGER(iwp) :: ilen 657 INTEGER(iwp) :: k 658 659 SELECT CASE ( TRIM( var ) ) 660 661 CASE ( 'm_soil' ) 662 IF ( .NOT. land_surface ) THEN 663 message_string = 'output of "' // TRIM( var ) // '" requi' // & 664 'res land_surface = .TRUE.' 665 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 666 ENDIF 667 unit = 'm3/m3' 668 669 CASE ( 't_soil' ) 670 IF ( .NOT. land_surface ) THEN 671 message_string = 'output of "' // TRIM( var ) // '" requi' // & 672 'res land_surface = .TRUE.' 673 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 674 ENDIF 675 unit = 'K' 674 676 675 CASE ( 't_soil' ) 676 IF ( .NOT. land_surface ) THEN 677 message_string = 'output of "' // TRIM( var ) // '" requi' // & 678 'res land_surface = .TRUE.' 679 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 680 ENDIF 681 unit = 'K' 677 CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'm_liq_eb*',& 678 'qsws_eb*', 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*', & 679 'r_a*', 'r_s*', 'shf_eb*' ) 680 IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN 681 message_string = 'illegal value for data_output: "' // & 682 TRIM( var ) // '" & only 2d-horizontal ' // & 683 'cross sections are allowed for this value' 684 CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 ) 685 ENDIF 686 IF ( TRIM( var ) == 'lai*' .AND. .NOT. land_surface ) THEN 687 message_string = 'output of "' // TRIM( var ) // '" requi' // & 688 'res land_surface = .TRUE.' 689 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 690 ENDIF 691 IF ( TRIM( var ) == 'c_liq*' .AND. .NOT. land_surface ) THEN 692 message_string = 'output of "' // TRIM( var ) // '" requi' // & 693 'res land_surface = .TRUE.' 694 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 695 ENDIF 696 IF ( TRIM( var ) == 'c_soil*' .AND. .NOT. land_surface ) THEN 697 message_string = 'output of "' // TRIM( var ) // '" requi' // & 698 'res land_surface = .TRUE.' 699 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 700 ENDIF 701 IF ( TRIM( var ) == 'c_veg*' .AND. .NOT. land_surface ) THEN 702 message_string = 'output of "' // TRIM( var ) // '" requi' // & 703 'res land_surface = .TRUE.' 704 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 705 ENDIF 706 IF ( TRIM( var ) == 'ghf_eb*' .AND. .NOT. land_surface ) THEN 707 message_string = 'output of "' // TRIM( var ) // '" requi' // & 708 'res land_surface = .TRUE.' 709 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 710 ENDIF 711 IF ( TRIM( var ) == 'm_liq_eb*' .AND. .NOT. land_surface ) THEN 712 message_string = 'output of "' // TRIM( var ) // '" requi' // & 713 'res land_surface = .TRUE.' 714 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 715 ENDIF 716 IF ( TRIM( var ) == 'qsws_eb*' .AND. .NOT. land_surface ) THEN 717 message_string = 'output of "' // TRIM( var ) // '" requi' // & 718 'res land_surface = .TRUE.' 719 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 720 ENDIF 721 IF ( TRIM( var ) == 'qsws_liq_eb*' .AND. .NOT. land_surface ) & 722 THEN 723 message_string = 'output of "' // TRIM( var ) // '" requi' // & 724 'res land_surface = .TRUE.' 725 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 726 ENDIF 727 IF ( TRIM( var ) == 'qsws_soil_eb*' .AND. .NOT. land_surface ) & 728 THEN 729 message_string = 'output of "' // TRIM( var ) // '" requi' // & 730 'res land_surface = .TRUE.' 731 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 732 ENDIF 733 IF ( TRIM( var ) == 'qsws_veg_eb*' .AND. .NOT. land_surface ) & 734 THEN 735 message_string = 'output of "' // TRIM( var ) // '" requi' // & 736 'res land_surface = .TRUE.' 737 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 738 ENDIF 739 IF ( TRIM( var ) == 'r_a*' .AND. .NOT. land_surface ) & 740 THEN 741 message_string = 'output of "' // TRIM( var ) // '" requi' // & 742 'res land_surface = .TRUE.' 743 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 744 ENDIF 745 IF ( TRIM( var ) == 'r_s*' .AND. .NOT. land_surface ) & 746 THEN 747 message_string = 'output of "' // TRIM( var ) // '" requi' // & 748 'res land_surface = .TRUE.' 749 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 750 ENDIF 751 752 IF ( TRIM( var ) == 'lai*' ) unit = 'none' 753 IF ( TRIM( var ) == 'c_liq*' ) unit = 'none' 754 IF ( TRIM( var ) == 'c_soil*') unit = 'none' 755 IF ( TRIM( var ) == 'c_veg*' ) unit = 'none' 756 IF ( TRIM( var ) == 'ghf_eb*') unit = 'W/m2' 757 IF ( TRIM( var ) == 'm_liq_eb*' ) unit = 'm' 758 IF ( TRIM( var ) == 'qsws_eb*' ) unit = 'W/m2' 759 IF ( TRIM( var ) == 'qsws_liq_eb*' ) unit = 'W/m2' 760 IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2' 761 IF ( TRIM( var ) == 'qsws_veg_eb*' ) unit = 'W/m2' 762 IF ( TRIM( var ) == 'r_a*') unit = 's/m' 763 IF ( TRIM( var ) == 'r_s*') unit = 's/m' 764 IF ( TRIM( var ) == 'shf_eb*') unit = 'W/m2' 682 765 683 CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'm_liq_eb*',& 684 'qsws_eb*', 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*', & 685 'r_a*', 'r_s*', 'shf_eb*' ) 686 IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN 687 message_string = 'illegal value for data_output: "' // & 688 TRIM( var ) // '" & only 2d-horizontal ' // & 689 'cross sections are allowed for this value' 690 CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 ) 691 ENDIF 692 IF ( TRIM( var ) == 'lai*' .AND. .NOT. land_surface ) THEN 693 message_string = 'output of "' // TRIM( var ) // '" requi' // & 694 'res land_surface = .TRUE.' 695 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 696 ENDIF 697 IF ( TRIM( var ) == 'c_liq*' .AND. .NOT. land_surface ) THEN 698 message_string = 'output of "' // TRIM( var ) // '" requi' // & 699 'res land_surface = .TRUE.' 700 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 701 ENDIF 702 IF ( TRIM( var ) == 'c_soil*' .AND. .NOT. land_surface ) THEN 703 message_string = 'output of "' // TRIM( var ) // '" requi' // & 704 'res land_surface = .TRUE.' 705 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 706 ENDIF 707 IF ( TRIM( var ) == 'c_veg*' .AND. .NOT. land_surface ) THEN 708 message_string = 'output of "' // TRIM( var ) // '" requi' // & 709 'res land_surface = .TRUE.' 710 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 711 ENDIF 712 IF ( TRIM( var ) == 'ghf_eb*' .AND. .NOT. land_surface ) THEN 713 message_string = 'output of "' // TRIM( var ) // '" requi' // & 714 'res land_surface = .TRUE.' 715 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 716 ENDIF 717 IF ( TRIM( var ) == 'm_liq_eb*' .AND. .NOT. land_surface ) THEN 718 message_string = 'output of "' // TRIM( var ) // '" requi' // & 719 'res land_surface = .TRUE.' 720 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 721 ENDIF 722 IF ( TRIM( var ) == 'qsws_eb*' .AND. .NOT. land_surface ) THEN 723 message_string = 'output of "' // TRIM( var ) // '" requi' // & 724 'res land_surface = .TRUE.' 725 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 726 ENDIF 727 IF ( TRIM( var ) == 'qsws_liq_eb*' .AND. .NOT. land_surface ) & 728 THEN 729 message_string = 'output of "' // TRIM( var ) // '" requi' // & 730 'res land_surface = .TRUE.' 731 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 732 ENDIF 733 IF ( TRIM( var ) == 'qsws_soil_eb*' .AND. .NOT. land_surface ) & 734 THEN 735 message_string = 'output of "' // TRIM( var ) // '" requi' // & 736 'res land_surface = .TRUE.' 737 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 738 ENDIF 739 IF ( TRIM( var ) == 'qsws_veg_eb*' .AND. .NOT. land_surface ) & 740 THEN 741 message_string = 'output of "' // TRIM( var ) // '" requi' // & 742 'res land_surface = .TRUE.' 743 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 744 ENDIF 745 IF ( TRIM( var ) == 'r_a*' .AND. .NOT. land_surface ) & 746 THEN 747 message_string = 'output of "' // TRIM( var ) // '" requi' // & 748 'res land_surface = .TRUE.' 749 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 750 ENDIF 751 IF ( TRIM( var ) == 'r_s*' .AND. .NOT. land_surface ) & 752 THEN 753 message_string = 'output of "' // TRIM( var ) // '" requi' // & 754 'res land_surface = .TRUE.' 755 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 756 ENDIF 757 758 IF ( TRIM( var ) == 'lai*' ) unit = 'none' 759 IF ( TRIM( var ) == 'c_liq*' ) unit = 'none' 760 IF ( TRIM( var ) == 'c_soil*') unit = 'none' 761 IF ( TRIM( var ) == 'c_veg*' ) unit = 'none' 762 IF ( TRIM( var ) == 'ghf_eb*') unit = 'W/m2' 763 IF ( TRIM( var ) == 'm_liq_eb*' ) unit = 'm' 764 IF ( TRIM( var ) == 'qsws_eb*' ) unit = 'W/m2' 765 IF ( TRIM( var ) == 'qsws_liq_eb*' ) unit = 'W/m2' 766 IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2' 767 IF ( TRIM( var ) == 'qsws_veg_eb*' ) unit = 'W/m2' 768 IF ( TRIM( var ) == 'r_a*') unit = 's/m' 769 IF ( TRIM( var ) == 'r_s*') unit = 's/m' 770 IF ( TRIM( var ) == 'shf_eb*') unit = 'W/m2' 771 772 CASE DEFAULT 773 unit = 'illegal' 774 775 END SELECT 776 777 778 END SUBROUTINE lsm_check_data_output 779 766 CASE DEFAULT 767 unit = 'illegal' 768 769 END SELECT 770 771 772 END SUBROUTINE lsm_check_data_output 773 774 775 !------------------------------------------------------------------------------! 776 ! Description: 777 ! ------------ 778 !> Check data output of profiles for land surface model 779 !------------------------------------------------------------------------------! 780 780 SUBROUTINE lsm_check_data_output_pr( variable, var_count, unit, dopr_unit ) 781 781 782 783 ONLY:data_output_pr, message_string784 785 786 787 788 789 790 791 782 USE control_parameters, & 783 ONLY: data_output_pr, message_string 784 785 USE indices 786 787 USE profil_parameter 788 789 USE statistics 790 791 IMPLICIT NONE 792 792 793 794 795 793 CHARACTER (LEN=*) :: unit !< 794 CHARACTER (LEN=*) :: variable !< 795 CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit 796 796 797 798 799 800 797 INTEGER(iwp) :: user_pr_index !< 798 INTEGER(iwp) :: var_count !< 799 800 SELECT CASE ( TRIM( variable ) ) 801 801 802 CASE ( 't_soil', '#t_soil' ) 803 IF ( .NOT. land_surface ) THEN 804 message_string = 'data_output_pr = ' // & 805 TRIM( data_output_pr(var_count) ) // ' is' // & 806 'not implemented for land_surface = .FALSE.' 807 CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 ) 808 ELSE 809 dopr_index(var_count) = 89 810 dopr_unit = 'K' 811 hom(0:nzs-1,2,89,:) = SPREAD( - zs, 2, statistic_regions+1 ) 812 IF ( data_output_pr(var_count)(1:1) == '#' ) THEN 813 dopr_initial_index(var_count) = 90 814 hom(0:nzs-1,2,90,:) = SPREAD( - zs, 2, statistic_regions+1 ) 815 data_output_pr(var_count) = data_output_pr(var_count)(2:) 816 ENDIF 817 unit = dopr_unit 802 CASE ( 't_soil', '#t_soil' ) 803 IF ( .NOT. land_surface ) THEN 804 message_string = 'data_output_pr = ' // & 805 TRIM( data_output_pr(var_count) ) // ' is' // & 806 'not implemented for land_surface = .FALSE.' 807 CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 ) 808 ELSE 809 dopr_index(var_count) = 89 810 dopr_unit = 'K' 811 hom(0:nzs-1,2,89,:) = SPREAD( - zs, 2, statistic_regions+1 ) 812 IF ( data_output_pr(var_count)(1:1) == '#' ) THEN 813 dopr_initial_index(var_count) = 90 814 hom(0:nzs-1,2,90,:) = SPREAD( - zs, 2, statistic_regions+1 ) 815 data_output_pr(var_count) = data_output_pr(var_count)(2:) 818 816 ENDIF 819 820 CASE ( 'm_soil', '#m_soil' )821 IF ( .NOT. land_surface ) THEN 822 message_string = 'data_output_pr = ' // &823 TRIM( data_output_pr(var_count) ) // ' is' // &824 ' not implemented for land_surface = .FALSE.'825 CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )826 ELSE827 dopr_index(var_count) = 91828 dopr_unit = 'm3/m3'829 hom(0:nzs-1,2,91,:) = SPREAD( - zs, 2, statistic_regions+1 )830 IF ( data_output_pr(var_count)(1:1) == '#' ) THEN831 dopr_initial_index(var_count) = 92832 hom(0:nzs-1,2,92,:) = SPREAD( - zs, 2, statistic_regions+1 )833 data_output_pr(var_count) = data_output_pr(var_count)(2:)834 ENDIF835 unit = dopr_unit817 unit = dopr_unit 818 ENDIF 819 820 CASE ( 'm_soil', '#m_soil' ) 821 IF ( .NOT. land_surface ) THEN 822 message_string = 'data_output_pr = ' // & 823 TRIM( data_output_pr(var_count) ) // ' is' // & 824 ' not implemented for land_surface = .FALSE.' 825 CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 ) 826 ELSE 827 dopr_index(var_count) = 91 828 dopr_unit = 'm3/m3' 829 hom(0:nzs-1,2,91,:) = SPREAD( - zs, 2, statistic_regions+1 ) 830 IF ( data_output_pr(var_count)(1:1) == '#' ) THEN 831 dopr_initial_index(var_count) = 92 832 hom(0:nzs-1,2,92,:) = SPREAD( - zs, 2, statistic_regions+1 ) 833 data_output_pr(var_count) = data_output_pr(var_count)(2:) 836 834 ENDIF 837 838 839 CASE DEFAULT 840 unit = 'illegal' 841 842 END SELECT 843 844 845 END SUBROUTINE lsm_check_data_output_pr 835 unit = dopr_unit 836 ENDIF 837 838 839 CASE DEFAULT 840 unit = 'illegal' 841 842 END SELECT 843 844 845 END SUBROUTINE lsm_check_data_output_pr 846 846 847 847 … … 851 851 !> Check parameters routine for land surface model 852 852 !------------------------------------------------------------------------------! 853 854 855 USE control_parameters,&856 ONLY: bc_pt_b, bc_q_b, constant_flux_layer, message_string,&857 853 SUBROUTINE lsm_check_parameters 854 855 USE control_parameters, & 856 ONLY: bc_pt_b, bc_q_b, constant_flux_layer, message_string, & 857 most_method, topography 858 858 859 USE radiation_model_mod,&860 ONLY:radiation859 USE radiation_model_mod, & 860 ONLY: radiation 861 861 862 862 863 863 IMPLICIT NONE 864 864 865 865 866 866 ! 867 !-- Dirichlet boundary conditions are required as the surface fluxes are 868 !-- calculated from the temperature/humidity gradients in the land surface 869 !-- model 870 IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' ) THEN 871 message_string = 'lsm requires setting of'// & 872 'bc_pt_b = "dirichlet" and '// & 873 'bc_q_b = "dirichlet"' 874 CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 ) 867 !-- Dirichlet boundary conditions are required as the surface fluxes are 868 !-- calculated from the temperature/humidity gradients in the land surface 869 !-- model 870 IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' ) THEN 871 message_string = 'lsm requires setting of'// & 872 'bc_pt_b = "dirichlet" and '// & 873 'bc_q_b = "dirichlet"' 874 CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 ) 875 ENDIF 876 877 IF ( .NOT. constant_flux_layer ) THEN 878 message_string = 'lsm requires '// & 879 'constant_flux_layer = .T.' 880 CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 ) 881 ENDIF 882 883 IF ( topography /= 'flat' ) THEN 884 message_string = 'lsm cannot be used ' // & 885 'in combination with topography /= "flat"' 886 CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 ) 887 ENDIF 888 889 IF ( ( veg_type == 14 .OR. veg_type == 15 ) .AND. & 890 most_method == 'lookup' ) THEN 891 WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ', & 892 'allowed in combination with ', & 893 'most_method = ', most_method 894 CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 ) 895 ENDIF 896 897 IF ( veg_type == 0 ) THEN 898 IF ( SUM( root_fraction ) /= 1.0_wp ) THEN 899 message_string = 'veg_type = 0 (user_defined)'// & 900 'requires setting of root_fraction(0:3)'// & 901 '/= 9999999.9 and SUM(root_fraction) = 1' 902 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 875 903 ENDIF 876 877 IF ( .NOT. constant_flux_layer ) THEN 878 message_string = 'lsm requires '// & 879 'constant_flux_layer = .T.' 880 CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 ) 904 905 IF ( min_canopy_resistance == 9999999.9_wp ) THEN 906 message_string = 'veg_type = 0 (user defined)'// & 907 'requires setting of min_canopy_resistance'// & 908 '/= 9999999.9' 909 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 881 910 ENDIF 882 911 883 IF ( topography /= 'flat' ) THEN 884 message_string = 'lsm cannot be used ' // & 885 'in combination with topography /= "flat"' 886 CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 ) 912 IF ( leaf_area_index == 9999999.9_wp ) THEN 913 message_string = 'veg_type = 0 (user_defined)'// & 914 'requires setting of leaf_area_index'// & 915 '/= 9999999.9' 916 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 887 917 ENDIF 888 918 889 IF ( ( veg_type == 14 .OR. veg_type == 15 ) .AND. & 890 most_method == 'lookup' ) THEN 891 WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ', & 892 'allowed in combination with ', & 893 'most_method = ', most_method 894 CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 ) 919 IF ( vegetation_coverage == 9999999.9_wp ) THEN 920 message_string = 'veg_type = 0 (user_defined)'// & 921 'requires setting of vegetation_coverage'// & 922 '/= 9999999.9' 923 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 895 924 ENDIF 896 925 897 IF ( veg_type == 0 ) THEN 898 IF ( SUM( root_fraction ) /= 1.0_wp ) THEN 899 message_string = 'veg_type = 0 (user_defined)'// & 900 'requires setting of root_fraction(0:3)'// & 901 '/= 9999999.9 and SUM(root_fraction) = 1' 902 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 903 ENDIF 904 905 IF ( min_canopy_resistance == 9999999.9_wp ) THEN 906 message_string = 'veg_type = 0 (user defined)'// & 907 'requires setting of min_canopy_resistance'// & 908 '/= 9999999.9' 909 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 910 ENDIF 911 912 IF ( leaf_area_index == 9999999.9_wp ) THEN 913 message_string = 'veg_type = 0 (user_defined)'// & 914 'requires setting of leaf_area_index'// & 915 '/= 9999999.9' 916 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 917 ENDIF 918 919 IF ( vegetation_coverage == 9999999.9_wp ) THEN 920 message_string = 'veg_type = 0 (user_defined)'// & 921 'requires setting of vegetation_coverage'// & 922 '/= 9999999.9' 923 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 924 ENDIF 925 926 IF ( canopy_resistance_coefficient == 9999999.9_wp) THEN 927 message_string = 'veg_type = 0 (user_defined)'// & 928 'requires setting of'// & 929 'canopy_resistance_coefficient /= 9999999.9' 930 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 931 ENDIF 932 933 IF ( lambda_surface_stable == 9999999.9_wp ) THEN 934 message_string = 'veg_type = 0 (user_defined)'// & 935 'requires setting of lambda_surface_stable'// & 936 '/= 9999999.9' 937 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 938 ENDIF 939 940 IF ( lambda_surface_unstable == 9999999.9_wp ) THEN 941 message_string = 'veg_type = 0 (user_defined)'// & 942 'requires setting of lambda_surface_unstable'// & 943 '/= 9999999.9' 944 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 945 ENDIF 946 947 IF ( f_shortwave_incoming == 9999999.9_wp ) THEN 948 message_string = 'veg_type = 0 (user_defined)'// & 949 'requires setting of f_shortwave_incoming'// & 950 '/= 9999999.9' 951 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 952 ENDIF 953 954 IF ( z0_eb == 9999999.9_wp ) THEN 955 message_string = 'veg_type = 0 (user_defined)'// & 956 'requires setting of z0_eb'// & 957 '/= 9999999.9' 958 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 959 ENDIF 960 961 IF ( z0h_eb == 9999999.9_wp ) THEN 962 message_string = 'veg_type = 0 (user_defined)'// & 963 'requires setting of z0h_eb'// & 964 '/= 9999999.9' 965 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 966 ENDIF 967 968 926 IF ( canopy_resistance_coefficient == 9999999.9_wp) THEN 927 message_string = 'veg_type = 0 (user_defined)'// & 928 'requires setting of'// & 929 'canopy_resistance_coefficient /= 9999999.9' 930 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 969 931 ENDIF 970 932 971 IF ( soil_type == 0 ) THEN 972 973 IF ( alpha_vangenuchten == 9999999.9_wp ) THEN 974 message_string = 'soil_type = 0 (user_defined)'// & 975 'requires setting of alpha_vangenuchten'// & 976 '/= 9999999.9' 977 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 978 ENDIF 979 980 IF ( l_vangenuchten == 9999999.9_wp ) THEN 981 message_string = 'soil_type = 0 (user_defined)'// & 982 'requires setting of l_vangenuchten'// & 983 '/= 9999999.9' 984 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 985 ENDIF 986 987 IF ( n_vangenuchten == 9999999.9_wp ) THEN 988 message_string = 'soil_type = 0 (user_defined)'// & 989 'requires setting of n_vangenuchten'// & 990 '/= 9999999.9' 991 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 992 ENDIF 993 994 IF ( hydraulic_conductivity == 9999999.9_wp ) THEN 995 message_string = 'soil_type = 0 (user_defined)'// & 996 'requires setting of hydraulic_conductivity'// & 997 '/= 9999999.9' 998 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 999 ENDIF 1000 1001 IF ( saturation_moisture == 9999999.9_wp ) THEN 1002 message_string = 'soil_type = 0 (user_defined)'// & 1003 'requires setting of saturation_moisture'// & 1004 '/= 9999999.9' 1005 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1006 ENDIF 1007 1008 IF ( field_capacity == 9999999.9_wp ) THEN 1009 message_string = 'soil_type = 0 (user_defined)'// & 1010 'requires setting of field_capacity'// & 1011 '/= 9999999.9' 1012 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1013 ENDIF 1014 1015 IF ( wilting_point == 9999999.9_wp ) THEN 1016 message_string = 'soil_type = 0 (user_defined)'// & 1017 'requires setting of wilting_point'// & 1018 '/= 9999999.9' 1019 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1020 ENDIF 1021 1022 IF ( residual_moisture == 9999999.9_wp ) THEN 1023 message_string = 'soil_type = 0 (user_defined)'// & 1024 'requires setting of residual_moisture'// & 1025 '/= 9999999.9' 1026 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1027 ENDIF 1028 933 IF ( lambda_surface_stable == 9999999.9_wp ) THEN 934 message_string = 'veg_type = 0 (user_defined)'// & 935 'requires setting of lambda_surface_stable'// & 936 '/= 9999999.9' 937 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 1029 938 ENDIF 1030 939 1031 IF ( .NOT. radiation ) THEN 1032 message_string = 'lsm requires '// & 1033 'radiation = .T.' 1034 CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 ) 940 IF ( lambda_surface_unstable == 9999999.9_wp ) THEN 941 message_string = 'veg_type = 0 (user_defined)'// & 942 'requires setting of lambda_surface_unstable'// & 943 '/= 9999999.9' 944 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 1035 945 ENDIF 946 947 IF ( f_shortwave_incoming == 9999999.9_wp ) THEN 948 message_string = 'veg_type = 0 (user_defined)'// & 949 'requires setting of f_shortwave_incoming'// & 950 '/= 9999999.9' 951 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 952 ENDIF 953 954 IF ( z0_eb == 9999999.9_wp ) THEN 955 message_string = 'veg_type = 0 (user_defined)'// & 956 'requires setting of z0_eb'// & 957 '/= 9999999.9' 958 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 959 ENDIF 960 961 IF ( z0h_eb == 9999999.9_wp ) THEN 962 message_string = 'veg_type = 0 (user_defined)'// & 963 'requires setting of z0h_eb'// & 964 '/= 9999999.9' 965 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 966 ENDIF 967 968 969 ENDIF 970 971 IF ( soil_type == 0 ) THEN 972 973 IF ( alpha_vangenuchten == 9999999.9_wp ) THEN 974 message_string = 'soil_type = 0 (user_defined)'// & 975 'requires setting of alpha_vangenuchten'// & 976 '/= 9999999.9' 977 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 978 ENDIF 979 980 IF ( l_vangenuchten == 9999999.9_wp ) THEN 981 message_string = 'soil_type = 0 (user_defined)'// & 982 'requires setting of l_vangenuchten'// & 983 '/= 9999999.9' 984 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 985 ENDIF 986 987 IF ( n_vangenuchten == 9999999.9_wp ) THEN 988 message_string = 'soil_type = 0 (user_defined)'// & 989 'requires setting of n_vangenuchten'// & 990 '/= 9999999.9' 991 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 992 ENDIF 993 994 IF ( hydraulic_conductivity == 9999999.9_wp ) THEN 995 message_string = 'soil_type = 0 (user_defined)'// & 996 'requires setting of hydraulic_conductivity'// & 997 '/= 9999999.9' 998 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 999 ENDIF 1000 1001 IF ( saturation_moisture == 9999999.9_wp ) THEN 1002 message_string = 'soil_type = 0 (user_defined)'// & 1003 'requires setting of saturation_moisture'// & 1004 '/= 9999999.9' 1005 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1006 ENDIF 1007 1008 IF ( field_capacity == 9999999.9_wp ) THEN 1009 message_string = 'soil_type = 0 (user_defined)'// & 1010 'requires setting of field_capacity'// & 1011 '/= 9999999.9' 1012 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1013 ENDIF 1014 1015 IF ( wilting_point == 9999999.9_wp ) THEN 1016 message_string = 'soil_type = 0 (user_defined)'// & 1017 'requires setting of wilting_point'// & 1018 '/= 9999999.9' 1019 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1020 ENDIF 1021 1022 IF ( residual_moisture == 9999999.9_wp ) THEN 1023 message_string = 'soil_type = 0 (user_defined)'// & 1024 'requires setting of residual_moisture'// & 1025 '/= 9999999.9' 1026 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1027 ENDIF 1028 1029 ENDIF 1030 1031 IF ( .NOT. radiation ) THEN 1032 message_string = 'lsm requires '// & 1033 'radiation = .T.' 1034 CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 ) 1035 ENDIF 1036 1036 1037 1037 1038 1038 END SUBROUTINE lsm_check_parameters 1039 1039 1040 1040 !------------------------------------------------------------------------------! … … 1043 1043 !> Solver for the energy balance at the surface. 1044 1044 !------------------------------------------------------------------------------! 1045 SUBROUTINE lsm_energy_balance 1046 1047 1048 IMPLICIT NONE 1049 1050 INTEGER(iwp) :: i !< running index 1051 INTEGER(iwp) :: j !< running index 1052 INTEGER(iwp) :: k, ks !< running index 1053 1054 REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface 1055 f1, & !< resistance correction term 1 1056 f2, & !< resistance correction term 2 1057 f3, & !< resistance correction term 3 1058 m_min, & !< minimum soil moisture 1059 e, & !< water vapour pressure 1060 e_s, & !< water vapour saturation pressure 1061 e_s_dt, & !< derivate of e_s with respect to T 1062 tend, & !< tendency 1063 dq_s_dt, & !< derivate of q_s with respect to T 1064 coef_1, & !< coef. for prognostic equation 1065 coef_2, & !< coef. for prognostic equation 1066 f_qsws, & !< factor for qsws_eb 1067 f_qsws_veg, & !< factor for qsws_veg_eb 1068 f_qsws_soil, & !< factor for qsws_soil_eb 1069 f_qsws_liq, & !< factor for qsws_liq_eb 1070 f_shf, & !< factor for shf_eb 1071 lambda_surface, & !< Current value of lambda_surface 1072 m_liq_eb_max, & !< maxmimum value of the liq. water reservoir 1073 pt1, & !< potential temperature at first grid level 1074 qv1 !< specific humidity at first grid level 1075 1076 ! 1077 !-- Calculate the exner function for the current time step 1078 exn = ( surface_pressure / 1000.0_wp )**0.286_wp 1079 1080 DO i = nxlg, nxrg 1081 DO j = nysg, nyng 1082 k = nzb_s_inner(j,i) 1083 1084 ! 1085 !-- Set lambda_surface according to stratification between skin layer and soil 1086 IF ( .NOT. pave_surface(j,i) ) THEN 1087 1088 c_surface_tmp = c_surface 1089 1090 IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i)) THEN 1091 lambda_surface = lambda_surface_s(j,i) 1092 ELSE 1093 lambda_surface = lambda_surface_u(j,i) 1094 ENDIF 1045 SUBROUTINE lsm_energy_balance 1046 1047 1048 IMPLICIT NONE 1049 1050 INTEGER(iwp) :: i !< running index 1051 INTEGER(iwp) :: j !< running index 1052 INTEGER(iwp) :: k, ks !< running index 1053 1054 REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface 1055 f1, & !< resistance correction term 1 1056 f2, & !< resistance correction term 2 1057 f3, & !< resistance correction term 3 1058 m_min, & !< minimum soil moisture 1059 e, & !< water vapour pressure 1060 e_s, & !< water vapour saturation pressure 1061 e_s_dt, & !< derivate of e_s with respect to T 1062 tend, & !< tendency 1063 dq_s_dt, & !< derivate of q_s with respect to T 1064 coef_1, & !< coef. for prognostic equation 1065 coef_2, & !< coef. for prognostic equation 1066 f_qsws, & !< factor for qsws_eb 1067 f_qsws_veg, & !< factor for qsws_veg_eb 1068 f_qsws_soil, & !< factor for qsws_soil_eb 1069 f_qsws_liq, & !< factor for qsws_liq_eb 1070 f_shf, & !< factor for shf_eb 1071 lambda_surface, & !< Current value of lambda_surface 1072 m_liq_eb_max, & !< maxmimum value of the liq. water reservoir 1073 pt1, & !< potential temperature at first grid level 1074 qv1 !< specific humidity at first grid level 1075 1076 ! 1077 !-- Calculate the exner function for the current time step 1078 exn = ( surface_pressure / 1000.0_wp )**0.286_wp 1079 1080 DO i = nxlg, nxrg 1081 DO j = nysg, nyng 1082 k = nzb_s_inner(j,i) 1083 1084 ! 1085 !-- Set lambda_surface according to stratification between skin layer and soil 1086 IF ( .NOT. pave_surface(j,i) ) THEN 1087 1088 c_surface_tmp = c_surface 1089 1090 IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i)) THEN 1091 lambda_surface = lambda_surface_s(j,i) 1095 1092 ELSE 1096 1097 c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp 1098 lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil) 1099 1093 lambda_surface = lambda_surface_u(j,i) 1100 1094 ENDIF 1101 1102 ! 1103 !-- First step: calculate aerodyamic resistance. As pt, us, ts 1104 !-- are not available for the prognostic time step, data from the last 1105 !-- time step is used here. Note that this formulation is the 1106 !-- equivalent to the ECMWF formulation using drag coefficients 1107 IF ( cloud_physics ) THEN 1108 pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i) 1109 qv1 = q(k+1,j,i) - ql(k+1,j,i) 1110 ELSE 1111 pt1 = pt(k+1,j,i) 1112 qv1 = q(k+1,j,i) 1113 ENDIF 1114 1115 r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp) 1116 1117 ! 1118 !-- Make sure that the resistance does not drop to zero for neutral 1119 !-- stratification 1120 IF ( ABS(r_a(j,i)) < 1.0_wp ) r_a(j,i) = 1.0_wp 1121 1122 ! 1123 !-- Second step: calculate canopy resistance r_canopy 1124 !-- f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation 1095 ELSE 1096 1097 c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp 1098 lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil) 1099 1100 ENDIF 1101 1102 ! 1103 !-- First step: calculate aerodyamic resistance. As pt, us, ts 1104 !-- are not available for the prognostic time step, data from the last 1105 !-- time step is used here. Note that this formulation is the 1106 !-- equivalent to the ECMWF formulation using drag coefficients 1107 IF ( cloud_physics ) THEN 1108 pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i) 1109 qv1 = q(k+1,j,i) - ql(k+1,j,i) 1110 ELSE 1111 pt1 = pt(k+1,j,i) 1112 qv1 = q(k+1,j,i) 1113 ENDIF 1114 1115 r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp) 1116 1117 ! 1118 !-- Make sure that the resistance does not drop to zero for neutral 1119 !-- stratification 1120 IF ( ABS(r_a(j,i)) < 1.0_wp ) r_a(j,i) = 1.0_wp 1121 1122 ! 1123 !-- Second step: calculate canopy resistance r_canopy 1124 !-- f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation 1125 1125 1126 !-- f1: correction for incoming shortwave radiation (stomata close at 1127 !-- night) 1128 IF ( radiation_scheme /= 'constant' ) THEN 1129 f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) / & 1130 (0.81_wp * (0.004_wp * rad_sw_in(k,j,i) & 1131 + 1.0_wp)) ) 1132 ELSE 1133 f1 = 1.0_wp 1134 ENDIF 1135 1136 1137 ! 1138 !-- f2: correction for soil moisture availability to plants (the 1139 !-- integrated soil moisture must thus be considered here) 1140 !-- f2 = 0 for very dry soils 1141 m_total = 0.0_wp 1142 DO ks = nzb_soil, nzt_soil 1143 m_total = m_total + root_fr(ks,j,i) & 1144 * MAX(m_soil(ks,j,i),m_wilt(j,i)) 1145 ENDDO 1146 1147 IF ( m_total > m_wilt(j,i) .AND. m_total < m_fc(j,i) ) THEN 1148 f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) ) 1149 ELSEIF ( m_total >= m_fc(j,i) ) THEN 1150 f2 = 1.0_wp 1151 ELSE 1152 f2 = 1.0E-20_wp 1153 ENDIF 1154 1155 ! 1156 !-- Calculate water vapour pressure at saturation 1157 e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i) & 1158 - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) ) 1159 1160 ! 1161 !-- f3: correction for vapour pressure deficit 1162 IF ( g_d(j,i) /= 0.0_wp ) THEN 1163 ! 1164 !-- Calculate vapour pressure 1165 e = qv1 * surface_pressure / 0.622_wp 1166 f3 = EXP ( -g_d(j,i) * (e_s - e) ) 1167 ELSE 1168 f3 = 1.0_wp 1169 ENDIF 1170 1171 ! 1172 !-- Calculate canopy resistance. In case that c_veg is 0 (bare soils), 1173 !-- this calculation is obsolete, as r_canopy is not used below. 1174 !-- To do: check for very dry soil -> r_canopy goes to infinity 1175 r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3 & 1176 + 1.0E-20_wp) 1177 1178 ! 1179 !-- Third step: calculate bare soil resistance r_soil. The Clapp & 1180 !-- Hornberger parametrization does not consider c_veg. 1181 IF ( soil_type_2d(j,i) /= 7 ) THEN 1182 m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) * & 1183 m_res(j,i) 1184 ELSE 1185 m_min = m_wilt(j,i) 1186 ENDIF 1187 1188 f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min ) 1189 f2 = MAX(f2,1.0E-20_wp) 1190 f2 = MIN(f2,1.0_wp) 1191 1192 r_soil(j,i) = r_soil_min(j,i) / f2 1193 1194 ! 1195 !-- Calculate the maximum possible liquid water amount on plants and 1196 !-- bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is 1197 !-- assumed, while paved surfaces might hold up 1 mm of water. The 1198 !-- liquid water fraction for paved surfaces is calculated after 1199 !-- Noilhan & Planton (1989), while the ECMWF formulation is used for 1200 !-- vegetated surfaces and bare soils. 1201 IF ( pave_surface(j,i) ) THEN 1202 m_liq_eb_max = m_max_depth * 5.0_wp 1203 c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 ) 1204 ELSE 1205 m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i) & 1206 + (1.0_wp - c_veg(j,i)) ) 1207 c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max ) 1208 ENDIF 1209 1210 ! 1211 !-- Calculate saturation specific humidity 1212 q_s = 0.622_wp * e_s / surface_pressure 1213 1214 ! 1215 !-- In case of dewfall, set evapotranspiration to zero 1216 !-- All super-saturated water is then removed from the air 1217 IF ( humidity .AND. q_s <= qv1 ) THEN 1218 r_canopy(j,i) = 0.0_wp 1219 r_soil(j,i) = 0.0_wp 1220 ENDIF 1221 1222 ! 1223 !-- Calculate coefficients for the total evapotranspiration 1224 !-- In case of water surface, set vegetation and soil fluxes to zero. 1225 !-- For pavements, only evaporation of liquid water is possible. 1226 IF ( water_surface(j,i) ) THEN 1227 f_qsws_veg = 0.0_wp 1228 f_qsws_soil = 0.0_wp 1229 f_qsws_liq = rho_lv / r_a(j,i) 1230 ELSEIF ( pave_surface (j,i) ) THEN 1231 f_qsws_veg = 0.0_wp 1232 f_qsws_soil = 0.0_wp 1233 f_qsws_liq = rho_lv * c_liq(j,i) / r_a(j,i) 1234 ELSE 1235 f_qsws_veg = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/ & 1236 (r_a(j,i) + r_canopy(j,i)) 1237 f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) + & 1238 r_soil(j,i)) 1239 f_qsws_liq = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i) 1240 ENDIF 1241 ! 1242 !-- If soil moisture is below wilting point, plants do no longer 1243 !-- transpirate. 1244 ! IF ( m_soil(k,j,i) < m_wilt(j,i) ) THEN 1245 ! f_qsws_veg = 0.0_wp 1246 ! ENDIF 1247 1248 f_shf = rho_cp / r_a(j,i) 1249 f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq 1250 1251 ! 1252 !-- Calculate derivative of q_s for Taylor series expansion 1253 e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) - & 1254 17.269_wp*(t_surface(j,i) - 273.16_wp) & 1255 / (t_surface(j,i) - 35.86_wp)**2 ) 1256 1257 dq_s_dt = 0.622_wp * e_s_dt / surface_pressure 1258 1259 ! 1260 !-- Add LW up so that it can be removed in prognostic equation 1261 rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i) 1262 1263 ! 1264 !-- Calculate new skin temperature 1265 IF ( humidity ) THEN 1266 #if defined ( __rrtmg ) 1267 ! 1268 !-- Numerator of the prognostic equation 1269 coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1270 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1271 + f_shf * pt1 + f_qsws * ( qv1 - q_s & 1272 + dq_s_dt * t_surface(j,i) ) + lambda_surface & 1273 * t_soil(nzb_soil,j,i) 1274 1275 ! 1276 !-- Denominator of the prognostic equation 1277 coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt & 1278 + lambda_surface + f_shf / exn 1279 #else 1280 1281 ! 1282 !-- Numerator of the prognostic equation 1283 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1284 * t_surface(j,i) ** 4 & 1285 + f_shf * pt1 + f_qsws * ( qv1 & 1286 - q_s + dq_s_dt * t_surface(j,i) ) & 1287 + lambda_surface * t_soil(nzb_soil,j,i) 1288 1289 ! 1290 !-- Denominator of the prognostic equation 1291 coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws & 1292 * dq_s_dt + lambda_surface + f_shf / exn 1293 1294 #endif 1295 ELSE 1296 1297 #if defined ( __rrtmg ) 1298 ! 1299 !-- Numerator of the prognostic equation 1300 coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1301 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1302 + f_shf * pt1 + lambda_surface & 1303 * t_soil(nzb_soil,j,i) 1304 1305 ! 1306 !-- Denominator of the prognostic equation 1307 coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn 1308 #else 1309 1310 ! 1311 !-- Numerator of the prognostic equation 1312 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1313 * t_surface(j,i) ** 4 + f_shf * pt1 & 1314 + lambda_surface * t_soil(nzb_soil,j,i) 1315 1316 ! 1317 !-- Denominator of the prognostic equation 1318 coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 & 1319 + lambda_surface + f_shf / exn 1320 #endif 1321 ENDIF 1322 1323 tend = 0.0_wp 1324 1325 ! 1326 !-- Implicit solution when the surface layer has no heat capacity, 1327 !-- otherwise use RK3 scheme. 1328 t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp * & 1329 t_surface(j,i) ) / ( c_surface_tmp + coef_2 & 1126 !-- f1: correction for incoming shortwave radiation (stomata close at 1127 !-- night) 1128 f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) / & 1129 (0.81_wp * (0.004_wp * rad_sw_in(k,j,i) & 1130 + 1.0_wp)) ) 1131 1132 1133 1134 ! 1135 !-- f2: correction for soil moisture availability to plants (the 1136 !-- integrated soil moisture must thus be considered here) 1137 !-- f2 = 0 for very dry soils 1138 m_total = 0.0_wp 1139 DO ks = nzb_soil, nzt_soil 1140 m_total = m_total + root_fr(ks,j,i) & 1141 * MAX(m_soil(ks,j,i),m_wilt(j,i)) 1142 ENDDO 1143 1144 IF ( m_total > m_wilt(j,i) .AND. m_total < m_fc(j,i) ) THEN 1145 f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) ) 1146 ELSEIF ( m_total >= m_fc(j,i) ) THEN 1147 f2 = 1.0_wp 1148 ELSE 1149 f2 = 1.0E-20_wp 1150 ENDIF 1151 1152 ! 1153 !-- Calculate water vapour pressure at saturation 1154 e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i) & 1155 - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) ) 1156 1157 ! 1158 !-- f3: correction for vapour pressure deficit 1159 IF ( g_d(j,i) /= 0.0_wp ) THEN 1160 ! 1161 !-- Calculate vapour pressure 1162 e = qv1 * surface_pressure / 0.622_wp 1163 f3 = EXP ( -g_d(j,i) * (e_s - e) ) 1164 ELSE 1165 f3 = 1.0_wp 1166 ENDIF 1167 1168 ! 1169 !-- Calculate canopy resistance. In case that c_veg is 0 (bare soils), 1170 !-- this calculation is obsolete, as r_canopy is not used below. 1171 !-- To do: check for very dry soil -> r_canopy goes to infinity 1172 r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3 & 1173 + 1.0E-20_wp) 1174 1175 ! 1176 !-- Third step: calculate bare soil resistance r_soil. The Clapp & 1177 !-- Hornberger parametrization does not consider c_veg. 1178 IF ( soil_type_2d(j,i) /= 7 ) THEN 1179 m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) * & 1180 m_res(j,i) 1181 ELSE 1182 m_min = m_wilt(j,i) 1183 ENDIF 1184 1185 f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min ) 1186 f2 = MAX(f2,1.0E-20_wp) 1187 f2 = MIN(f2,1.0_wp) 1188 1189 r_soil(j,i) = r_soil_min(j,i) / f2 1190 1191 ! 1192 !-- Calculate the maximum possible liquid water amount on plants and 1193 !-- bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is 1194 !-- assumed, while paved surfaces might hold up 1 mm of water. The 1195 !-- liquid water fraction for paved surfaces is calculated after 1196 !-- Noilhan & Planton (1989), while the ECMWF formulation is used for 1197 !-- vegetated surfaces and bare soils. 1198 IF ( pave_surface(j,i) ) THEN 1199 m_liq_eb_max = m_max_depth * 5.0_wp 1200 c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 ) 1201 ELSE 1202 m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i) & 1203 + (1.0_wp - c_veg(j,i)) ) 1204 c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max ) 1205 ENDIF 1206 1207 ! 1208 !-- Calculate saturation specific humidity 1209 q_s = 0.622_wp * e_s / surface_pressure 1210 1211 ! 1212 !-- In case of dewfall, set evapotranspiration to zero 1213 !-- All super-saturated water is then removed from the air 1214 IF ( humidity .AND. q_s <= qv1 ) THEN 1215 r_canopy(j,i) = 0.0_wp 1216 r_soil(j,i) = 0.0_wp 1217 ENDIF 1218 1219 ! 1220 !-- Calculate coefficients for the total evapotranspiration 1221 !-- In case of water surface, set vegetation and soil fluxes to zero. 1222 !-- For pavements, only evaporation of liquid water is possible. 1223 IF ( water_surface(j,i) ) THEN 1224 f_qsws_veg = 0.0_wp 1225 f_qsws_soil = 0.0_wp 1226 f_qsws_liq = rho_lv / r_a(j,i) 1227 ELSEIF ( pave_surface (j,i) ) THEN 1228 f_qsws_veg = 0.0_wp 1229 f_qsws_soil = 0.0_wp 1230 f_qsws_liq = rho_lv * c_liq(j,i) / r_a(j,i) 1231 ELSE 1232 f_qsws_veg = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/ & 1233 (r_a(j,i) + r_canopy(j,i)) 1234 f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) + & 1235 r_soil(j,i)) 1236 f_qsws_liq = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i) 1237 ENDIF 1238 ! 1239 !-- If soil moisture is below wilting point, plants do no longer 1240 !-- transpirate. 1241 ! IF ( m_soil(k,j,i) < m_wilt(j,i) ) THEN 1242 ! f_qsws_veg = 0.0_wp 1243 ! ENDIF 1244 1245 f_shf = rho_cp / r_a(j,i) 1246 f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq 1247 1248 ! 1249 !-- Calculate derivative of q_s for Taylor series expansion 1250 e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) - & 1251 17.269_wp*(t_surface(j,i) - 273.16_wp) & 1252 / (t_surface(j,i) - 35.86_wp)**2 ) 1253 1254 dq_s_dt = 0.622_wp * e_s_dt / surface_pressure 1255 1256 ! 1257 !-- Add LW up so that it can be removed in prognostic equation 1258 rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i) 1259 1260 ! 1261 !-- Calculate new skin temperature 1262 IF ( humidity ) THEN 1263 1264 ! 1265 !-- Numerator of the prognostic equation 1266 coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1267 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1268 + f_shf * pt1 + f_qsws * ( qv1 - q_s & 1269 + dq_s_dt * t_surface(j,i) ) + lambda_surface & 1270 * t_soil(nzb_soil,j,i) 1271 1272 ! 1273 !-- Denominator of the prognostic equation 1274 coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt & 1275 + lambda_surface + f_shf / exn 1276 ELSE 1277 1278 ! 1279 !-- Numerator of the prognostic equation 1280 coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1281 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1282 + f_shf * pt1 + lambda_surface & 1283 * t_soil(nzb_soil,j,i) 1284 1285 ! 1286 !-- Denominator of the prognostic equation 1287 coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn 1288 1289 ENDIF 1290 1291 tend = 0.0_wp 1292 1293 ! 1294 !-- Implicit solution when the surface layer has no heat capacity, 1295 !-- otherwise use RK3 scheme. 1296 t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp * & 1297 t_surface(j,i) ) / ( c_surface_tmp + coef_2 & 1330 1298 * dt_3d * tsc(2) ) 1331 1299 1332 1300 ! 1333 !-- Add RK3 term 1334 IF ( c_surface_tmp /= 0.0_wp ) THEN 1335 1336 t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3) & 1337 * tt_surface_m(j,i) 1338 1339 ! 1340 !-- Calculate true tendency 1341 tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3) & 1342 * tt_surface_m(j,i)) / (dt_3d * tsc(2)) 1343 ! 1344 !-- Calculate t_surface tendencies for the next Runge-Kutta step 1345 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1346 IF ( intermediate_timestep_count == 1 ) THEN 1347 tt_surface_m(j,i) = tend 1348 ELSEIF ( intermediate_timestep_count < & 1349 intermediate_timestep_count_max ) THEN 1350 tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp & 1351 * tt_surface_m(j,i) 1352 ENDIF 1301 !-- Add RK3 term 1302 IF ( c_surface_tmp /= 0.0_wp ) THEN 1303 1304 t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3) & 1305 * tt_surface_m(j,i) 1306 1307 ! 1308 !-- Calculate true tendency 1309 tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3) & 1310 * tt_surface_m(j,i)) / (dt_3d * tsc(2)) 1311 ! 1312 !-- Calculate t_surface tendencies for the next Runge-Kutta step 1313 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1314 IF ( intermediate_timestep_count == 1 ) THEN 1315 tt_surface_m(j,i) = tend 1316 ELSEIF ( intermediate_timestep_count < & 1317 intermediate_timestep_count_max ) THEN 1318 tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp & 1319 * tt_surface_m(j,i) 1353 1320 ENDIF 1354 1321 ENDIF 1355 1356 ! 1357 !-- In case of fast changes in the skin temperature, it is possible to 1358 !-- update the radiative fluxes independently from the prescribed 1359 !-- radiation call frequency. This effectively prevents oscillations, 1360 !-- especially when setting skip_time_do_radiation /= 0. The threshold 1361 !-- value of 0.2 used here is just a first guess. This method should be 1362 !-- revised in the future as tests have shown that the threshold is 1363 !-- often reached, when no oscillations would occur (causes immense 1364 !-- computing time for the radiation code). 1365 IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp .AND. & 1366 unscheduled_radiation_calls ) THEN 1367 force_radiation_call_l = .TRUE. 1322 ENDIF 1323 1324 ! 1325 !-- In case of fast changes in the skin temperature, it is possible to 1326 !-- update the radiative fluxes independently from the prescribed 1327 !-- radiation call frequency. This effectively prevents oscillations, 1328 !-- especially when setting skip_time_do_radiation /= 0. The threshold 1329 !-- value of 0.2 used here is just a first guess. This method should be 1330 !-- revised in the future as tests have shown that the threshold is 1331 !-- often reached, when no oscillations would occur (causes immense 1332 !-- computing time for the radiation code). 1333 IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp .AND. & 1334 unscheduled_radiation_calls ) THEN 1335 force_radiation_call_l = .TRUE. 1336 ENDIF 1337 1338 pt(k,j,i) = t_surface_p(j,i) / exn 1339 1340 ! 1341 !-- Calculate fluxes 1342 rad_net_l(j,i) = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1343 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1344 - rad_lw_out_change_0(j,i) * t_surface_p(j,i) 1345 1346 rad_net(j,i) = rad_net_l(j,i) 1347 rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) + rad_lw_out_change_0(j,i) & 1348 * ( t_surface_p(j,i) - t_surface(j,i) ) 1349 1350 ghf_eb(j,i) = lambda_surface * (t_surface_p(j,i) & 1351 - t_soil(nzb_soil,j,i)) 1352 1353 shf_eb(j,i) = - f_shf * ( pt1 - pt(k,j,i) ) 1354 1355 shf(j,i) = shf_eb(j,i) / rho_cp 1356 1357 IF ( humidity ) THEN 1358 qsws_eb(j,i) = - f_qsws * ( qv1 - q_s + dq_s_dt & 1359 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) 1360 1361 qsws(j,i) = qsws_eb(j,i) / rho_lv 1362 1363 qsws_veg_eb(j,i) = - f_qsws_veg * ( qv1 - q_s & 1364 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1365 * t_surface_p(j,i) ) 1366 1367 qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s & 1368 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1369 * t_surface_p(j,i) ) 1370 1371 qsws_liq_eb(j,i) = - f_qsws_liq * ( qv1 - q_s & 1372 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1373 * t_surface_p(j,i) ) 1374 ENDIF 1375 1376 ! 1377 !-- Calculate the true surface resistance 1378 IF ( qsws_eb(j,i) == 0.0_wp ) THEN 1379 r_s(j,i) = 1.0E10_wp 1380 ELSE 1381 r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt & 1382 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) & 1383 / qsws_eb(j,i) - r_a(j,i) 1384 ENDIF 1385 1386 ! 1387 !-- Calculate change in liquid water reservoir due to dew fall or 1388 !-- evaporation of liquid water 1389 IF ( humidity ) THEN 1390 ! 1391 !-- If precipitation is activated, add rain water to qsws_liq_eb 1392 !-- and qsws_soil_eb according the the vegetation coverage. 1393 !-- precipitation_rate is given in mm. 1394 IF ( precipitation ) THEN 1395 1396 ! 1397 !-- Add precipitation to liquid water reservoir, if possible. 1398 !-- Otherwise, add the water to soil. In case of 1399 !-- pavements, the exceeding water amount is implicitely removed 1400 !-- as runoff as qsws_soil_eb is then not used in the soil model 1401 IF ( m_liq_eb(j,i) /= m_liq_eb_max ) THEN 1402 qsws_liq_eb(j,i) = qsws_liq_eb(j,i) & 1403 + c_veg(j,i) * prr(k,j,i) * hyrho(k) & 1404 * 0.001_wp * rho_l * l_v 1405 ELSE 1406 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) & 1407 + c_veg(j,i) * prr(k,j,i) * hyrho(k) & 1408 * 0.001_wp * rho_l * l_v 1409 ENDIF 1410 1411 !-- Add precipitation to bare soil according to the bare soil 1412 !-- coverage. 1413 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) + (1.0_wp & 1414 - c_veg(j,i)) * prr(k,j,i) * hyrho(k) & 1415 * 0.001_wp * rho_l * l_v 1368 1416 ENDIF 1369 1417 1370 pt(k,j,i) = t_surface_p(j,i) / exn 1371 1372 ! 1373 !-- Calculate fluxes 1374 #if defined ( __rrtmg ) 1375 rad_net_l(j,i) = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1376 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1377 - rad_lw_out_change_0(j,i) * t_surface_p(j,i) 1378 1379 IF ( rrtm_idrv == 1 ) THEN 1380 rad_net(j,i) = rad_net_l(j,i) 1381 rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) & 1382 + rad_lw_out_change_0(j,i) & 1383 * ( t_surface_p(j,i) - t_surface(j,i) ) 1418 ! 1419 !-- If the air is saturated, check the reservoir water level 1420 IF ( qsws_eb(j,i) < 0.0_wp ) THEN 1421 1422 ! 1423 !-- Check if reservoir is full (avoid values > m_liq_eb_max) 1424 !-- In that case, qsws_liq_eb goes to qsws_soil_eb. In this 1425 !-- case qsws_veg_eb is zero anyway (because c_liq = 1), 1426 !-- so that tend is zero and no further check is needed 1427 IF ( m_liq_eb(j,i) == m_liq_eb_max ) THEN 1428 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) & 1429 + qsws_liq_eb(j,i) 1430 1431 qsws_liq_eb(j,i) = 0.0_wp 1432 ENDIF 1433 1434 ! 1435 !-- In case qsws_veg_eb becomes negative (unphysical behavior), 1436 !-- let the water enter the liquid water reservoir as dew on the 1437 !-- plant 1438 IF ( qsws_veg_eb(j,i) < 0.0_wp ) THEN 1439 qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i) 1440 qsws_veg_eb(j,i) = 0.0_wp 1441 ENDIF 1442 ENDIF 1443 1444 tend = - qsws_liq_eb(j,i) * drho_l_lv 1445 m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend & 1446 + tsc(3) * tm_liq_eb_m(j,i) ) 1447 1448 ! 1449 !-- Check if reservoir is overfull -> reduce to maximum 1450 !-- (conservation of water is violated here) 1451 m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max) 1452 1453 ! 1454 !-- Check if reservoir is empty (avoid values < 0.0) 1455 !-- (conservation of water is violated here) 1456 m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp) 1457 1458 1459 ! 1460 !-- Calculate m_liq_eb tendencies for the next Runge-Kutta step 1461 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1462 IF ( intermediate_timestep_count == 1 ) THEN 1463 tm_liq_eb_m(j,i) = tend 1464 ELSEIF ( intermediate_timestep_count < & 1465 intermediate_timestep_count_max ) THEN 1466 tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp & 1467 * tm_liq_eb_m(j,i) 1468 ENDIF 1384 1469 ENDIF 1470 1471 ENDIF 1472 1473 ENDDO 1474 ENDDO 1475 1476 ! 1477 !-- Make a logical OR for all processes. Force radiation call if at 1478 !-- least one processor reached the threshold change in skin temperature 1479 IF ( unscheduled_radiation_calls .AND. intermediate_timestep_count & 1480 == intermediate_timestep_count_max-1 ) THEN 1481 #if defined( __parallel ) 1482 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1483 CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call, & 1484 1, MPI_LOGICAL, MPI_LOR, comm2d, ierr ) 1385 1485 #else 1386 rad_net_l(j,i) = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1387 * t_surface(j,i)**4 - 4.0_wp * sigma_sb & 1388 * t_surface(j,i)**3 * t_surface_p(j,i) 1486 force_radiation_call = force_radiation_call_l 1389 1487 #endif 1390 1391 ghf_eb(j,i) = lambda_surface * (t_surface_p(j,i) & 1392 - t_soil(nzb_soil,j,i)) 1393 1394 shf_eb(j,i) = - f_shf * ( pt1 - pt(k,j,i) ) 1395 1396 shf(j,i) = shf_eb(j,i) / rho_cp 1397 1398 IF ( humidity ) THEN 1399 qsws_eb(j,i) = - f_qsws * ( qv1 - q_s + dq_s_dt & 1400 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) 1401 1402 qsws(j,i) = qsws_eb(j,i) / rho_lv 1403 1404 qsws_veg_eb(j,i) = - f_qsws_veg * ( qv1 - q_s & 1405 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1406 * t_surface_p(j,i) ) 1407 1408 qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s & 1409 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1410 * t_surface_p(j,i) ) 1411 1412 qsws_liq_eb(j,i) = - f_qsws_liq * ( qv1 - q_s & 1413 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1414 * t_surface_p(j,i) ) 1415 ENDIF 1416 1417 ! 1418 !-- Calculate the true surface resistance 1419 IF ( qsws_eb(j,i) == 0.0_wp ) THEN 1420 r_s(j,i) = 1.0E10_wp 1421 ELSE 1422 r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt & 1423 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) & 1424 / qsws_eb(j,i) - r_a(j,i) 1425 ENDIF 1426 1427 ! 1428 !-- Calculate change in liquid water reservoir due to dew fall or 1429 !-- evaporation of liquid water 1430 IF ( humidity ) THEN 1431 ! 1432 !-- If precipitation is activated, add rain water to qsws_liq_eb 1433 !-- and qsws_soil_eb according the the vegetation coverage. 1434 !-- precipitation_rate is given in mm. 1435 IF ( precipitation ) THEN 1436 1437 ! 1438 !-- Add precipitation to liquid water reservoir, if possible. 1439 !-- Otherwise, add the water to soil. In case of 1440 !-- pavements, the exceeding water amount is implicitely removed 1441 !-- as runoff as qsws_soil_eb is then not used in the soil model 1442 IF ( m_liq_eb(j,i) /= m_liq_eb_max ) THEN 1443 qsws_liq_eb(j,i) = qsws_liq_eb(j,i) & 1444 + c_veg(j,i) * prr(k,j,i) * hyrho(k) & 1445 * 0.001_wp * rho_l * l_v 1446 ELSE 1447 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) & 1448 + c_veg(j,i) * prr(k,j,i) * hyrho(k) & 1449 * 0.001_wp * rho_l * l_v 1450 ENDIF 1451 1452 !-- Add precipitation to bare soil according to the bare soil 1453 !-- coverage. 1454 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) + (1.0_wp & 1455 - c_veg(j,i)) * prr(k,j,i) * hyrho(k) & 1456 * 0.001_wp * rho_l * l_v 1457 ENDIF 1458 1459 ! 1460 !-- If the air is saturated, check the reservoir water level 1461 IF ( qsws_eb(j,i) < 0.0_wp ) THEN 1462 1463 ! 1464 !-- Check if reservoir is full (avoid values > m_liq_eb_max) 1465 !-- In that case, qsws_liq_eb goes to qsws_soil_eb. In this 1466 !-- case qsws_veg_eb is zero anyway (because c_liq = 1), 1467 !-- so that tend is zero and no further check is needed 1468 IF ( m_liq_eb(j,i) == m_liq_eb_max ) THEN 1469 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) & 1470 + qsws_liq_eb(j,i) 1471 1472 qsws_liq_eb(j,i) = 0.0_wp 1473 ENDIF 1474 1475 ! 1476 !-- In case qsws_veg_eb becomes negative (unphysical behavior), 1477 !-- let the water enter the liquid water reservoir as dew on the 1478 !-- plant 1479 IF ( qsws_veg_eb(j,i) < 0.0_wp ) THEN 1480 qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i) 1481 qsws_veg_eb(j,i) = 0.0_wp 1482 ENDIF 1483 ENDIF 1484 1485 tend = - qsws_liq_eb(j,i) * drho_l_lv 1486 1487 m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend & 1488 + tsc(3) * tm_liq_eb_m(j,i) ) 1489 1490 ! 1491 !-- Check if reservoir is overfull -> reduce to maximum 1492 !-- (conservation of water is violated here) 1493 m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max) 1494 1495 ! 1496 !-- Check if reservoir is empty (avoid values < 0.0) 1497 !-- (conservation of water is violated here) 1498 m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp) 1499 1500 1501 ! 1502 !-- Calculate m_liq_eb tendencies for the next Runge-Kutta step 1503 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1504 IF ( intermediate_timestep_count == 1 ) THEN 1505 tm_liq_eb_m(j,i) = tend 1506 ELSEIF ( intermediate_timestep_count < & 1507 intermediate_timestep_count_max ) THEN 1508 tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp & 1509 * tm_liq_eb_m(j,i) 1510 ENDIF 1511 ENDIF 1512 1513 ENDIF 1514 1515 ENDDO 1516 ENDDO 1517 1518 ! 1519 !-- Make a logical OR for all processes. Force radiation call if at 1520 !-- least one processor reached the threshold change in skin temperature 1521 IF ( unscheduled_radiation_calls .AND. intermediate_timestep_count & 1522 == intermediate_timestep_count_max-1 ) THEN 1523 #if defined( __parallel ) 1524 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1525 CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call, & 1526 1, MPI_LOGICAL, MPI_LOR, comm2d, ierr ) 1527 #else 1528 force_radiation_call = force_radiation_call_l 1529 #endif 1530 force_radiation_call_l = .FALSE. 1531 ENDIF 1532 1533 ! 1534 !-- Calculate surface specific humidity 1535 IF ( humidity ) THEN 1536 CALL calc_q_surface 1537 ENDIF 1538 1539 ! 1540 !-- Calculate new roughness lengths (for water surfaces only) 1541 CALL calc_z0_water_surface 1542 1543 1544 END SUBROUTINE lsm_energy_balance 1488 force_radiation_call_l = .FALSE. 1489 ENDIF 1490 1491 ! 1492 !-- Calculate surface specific humidity 1493 IF ( humidity ) THEN 1494 CALL calc_q_surface 1495 ENDIF 1496 1497 ! 1498 !-- Calculate new roughness lengths (for water surfaces only) 1499 CALL calc_z0_water_surface 1500 1501 1502 END SUBROUTINE lsm_energy_balance 1503 1545 1504 1546 1505 !------------------------------------------------------------------------------! … … 2451 2410 ! Description: 2452 2411 ! ------------ 2453 !> S oubroutine for averaging 3D data2412 !> Subroutine for averaging 3D data 2454 2413 !------------------------------------------------------------------------------! 2455 2414 SUBROUTINE lsm_3d_data_averaging( mode, variable ) … … 2805 2764 ! Description: 2806 2765 ! ------------ 2807 !> S oubroutine definesappropriate grid for netcdf variables.2766 !> Subroutine defining appropriate grid for netcdf variables. 2808 2767 !> It is called out from subroutine netcdf. 2809 2768 !------------------------------------------------------------------------------! 2810 2769 SUBROUTINE lsm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z ) 2811 2770 2812 2813 2814 2815 2816 2817 2818 2819 2820 2821 2822 ! 2823 !-- 2824 2825 2826 CASE ( 'm_soil', 't_soil' )2827 found = .TRUE.2828 2829 2830 2831 2832 2833 2834 2835 2836 2837 2838 2839 2771 IMPLICIT NONE 2772 2773 CHARACTER (LEN=*), INTENT(IN) :: var !< 2774 LOGICAL, INTENT(OUT) :: found !< 2775 CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< 2776 CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< 2777 CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< 2778 2779 found = .TRUE. 2780 2781 ! 2782 !-- Check for the grid 2783 SELECT CASE ( TRIM( var ) ) 2784 2785 CASE ( 'm_soil', 't_soil', 'm_soil_xy', 't_soil_xy', 'm_soil_xz', & 2786 't_soil_xz', 'm_soil_yz', 't_soil_yz' ) 2787 grid_x = 'x' 2788 grid_y = 'y' 2789 grid_z = 'zs' 2790 2791 CASE DEFAULT 2792 found = .FALSE. 2793 grid_x = 'none' 2794 grid_y = 'none' 2795 grid_z = 'none' 2796 END SELECT 2797 2798 END SUBROUTINE lsm_define_netcdf_grid 2840 2799 2841 2800 !------------------------------------------------------------------------------! … … 2843 2802 ! Description: 2844 2803 ! ------------ 2845 !> S oubroutine defines3D output variables2804 !> Subroutine defining 3D output variables 2846 2805 !------------------------------------------------------------------------------! 2847 2806 SUBROUTINE lsm_data_output_2d( av, variable, found, grid, mode, local_pf, & 2848 2807 two_d, nzb_do, nzt_do ) 2849 2808 2850 2851 2809 USE indices 2852 2810 2853 2811 USE kinds 2854 2812 2855 USE user2856 2813 2857 2814 IMPLICIT NONE … … 3174 3131 ! Description: 3175 3132 ! ------------ 3176 !> S oubroutine defines3D output variables3133 !> Subroutine defining 3D output variables 3177 3134 !------------------------------------------------------------------------------! 3178 3135 SUBROUTINE lsm_data_output_3d( av, variable, found, local_pf )
Note: See TracChangeset
for help on using the changeset viewer.
