Changeset 2696 for palm/trunk/SOURCE/land_surface_model_mod.f90
- Timestamp:
- Dec 14, 2017 5:12:51 PM (6 years ago)
- Location:
- palm/trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk
-
palm/trunk/SOURCE
-
palm/trunk/SOURCE/land_surface_model_mod.f90
r2608 r2696 1 1 !> @file land_surface_model_mod.f90 2 2 !------------------------------------------------------------------------------! 3 ! This file is part of PALM.3 ! This file is part of the PALM model system. 4 4 ! 5 5 ! PALM is free software: you can redistribute it and/or modify it under the … … 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Bugfix: missing USE statement for calc_mean_profile 28 ! do not write surface temperatures onto pt array as this might cause 29 ! problems with nesting (MS) 30 ! Revised calculation of pt1 and qv1 (now done in surface_layer_fluxes). Bugfix 31 ! in calculation of surface density (cannot be done via an surface non-air 32 ! temperature) (BM) 33 ! Bugfix: g_d was NaN for non-vegetaed surface types (BM) 34 ! Bugfix initialization of c_veg and lai 35 ! Revise data output to enable _FillValues 36 ! Bugfix in calcultion of r_a and rad_net_l (MS) 37 ! Bugfix: rad_net is not updated in case of radiation_interaction and must thu 38 ! be calculated again from the radiative fluxes 39 ! Temporary fix for cases where no soil model is used on some PEs (BM) 40 ! Revised input and initialization of soil and surface paramters 41 ! pavement_depth is variable for each surface element 42 ! radiation quantities belong to surface type now 43 ! surface fractions initialized 44 ! Rename lsm_last_actions into lsm_write_restart_data (MS) 45 ! 46 ! 2608 2017-11-13 14:04:26Z schwenkel 27 47 ! Calculation of magnus equation in external module (diagnostic_quantities_mod). 28 48 ! Adjust calculation of vapor pressure and saturation mixing ratio that it is … … 265 285 !> DALES and UCLA-LES models. 266 286 !> 267 !> @todo Restart data for vertical natural land-surfaces268 287 !> @todo Extensive verification energy-balance solver for vertical surfaces, 269 288 !> e.g. parametrization of r_a … … 282 301 !> @note No time step criterion is required as long as the soil layers do not 283 302 !> become too thin. 303 !> @todo Attention, pavement_subpars_1/2 are hardcoded to 8 levels, in case 304 !> more levels are used this may cause an potential bug 305 !> @todo Routine calc_q_surface required? 284 306 !------------------------------------------------------------------------------! 285 307 MODULE land_surface_model_mod … … 287 309 USE arrays_3d, & 288 310 ONLY: hyp, pt, prr, q, q_p, ql, vpt, u, v, w 311 312 USE calc_mean_profile_mod, & 313 ONLY: calc_mean_profile 289 314 290 315 USE cloud_parameters, & … … 301 326 302 327 USE indices, & 303 ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb 328 ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb 329 330 USE netcdf_data_input_mod, & 331 ONLY : building_type_f, init_3d, input_pids_static, & 332 netcdf_data_input_interpolate, & 333 pavement_pars_f, pavement_subsurface_pars_f, pavement_type_f, & 334 root_area_density_lsm_f, soil_pars_f, soil_type_f, & 335 surface_fraction_f, vegetation_pars_f, vegetation_type_f, & 336 water_pars_f, water_type_f 304 337 305 338 USE kinds … … 308 341 309 342 USE radiation_model_mod, & 310 ONLY: albedo, albedo_type, emissivity, force_radiation_call, rad_net, & 311 rad_sw_in, rad_lw_out, rad_lw_out_change_0, radiation_scheme, & 312 unscheduled_radiation_calls 343 ONLY: albedo, albedo_type, emissivity, force_radiation_call, & 344 radiation_scheme, unscheduled_radiation_calls 313 345 314 346 USE statistics, & … … 316 348 317 349 USE surface_mod, & 318 ONLY : surf_lsm_h, surf_lsm_v, surf_type 350 ONLY : surf_lsm_h, surf_lsm_v, surf_type, surface_restore_elements 319 351 320 352 IMPLICIT NONE 321 353 322 354 TYPE surf_type_lsm 323 REAL(wp), DIMENSION(:), ALLOCATABLE :: var_1 D!< 1D prognostic variable324 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: var_2 D!< 2D prognostic variable355 REAL(wp), DIMENSION(:), ALLOCATABLE :: var_1d !< 1D prognostic variable 356 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: var_2d !< 2D prognostic variable 325 357 END TYPE surf_type_lsm 326 358 … … 355 387 ! 356 388 !-- LSM variables 357 CHARACTER(10) :: surface_type = ' vegetation'!< general classification. Allowed are:389 CHARACTER(10) :: surface_type = 'netcdf' !< general classification. Allowed are: 358 390 !< 'vegetation', 'pavement', ('building'), 359 391 !< 'water', and 'netcdf' … … 447 479 m_soil_v, & !< Soil moisture (m3/m3), vertical surface elements 448 480 m_soil_v_p !< Prog. soil moisture (m3/m3), vertical surface elements 481 449 482 #else 450 483 TYPE(surf_type_lsm), POINTER :: t_soil_h, & !< Soil temperature (K), horizontal surface elements … … 474 507 TYPE(surf_type_lsm), TARGET :: t_surface_h, & !< surface temperature (K), horizontal surface elements 475 508 t_surface_h_p, & !< progn. surface temperature (K), horizontal surface elements 476 m_liq_h, & !< liquid water reservoir (m), horizontal surface elements477 m_liq_h_p !< progn. liquid water reservoir (m), horizontal surface elements509 m_liq_h, & !< liquid water reservoir (m), horizontal surface elements 510 m_liq_h_p !< progn. liquid water reservoir (m), horizontal surface elements 478 511 479 512 TYPE(surf_type_lsm), DIMENSION(0:3), TARGET :: & … … 533 566 !-- Energy balance variables 534 567 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & 535 c_liq_av, & !< average of c_liq536 c_soil_av, & !< average of c_soil537 c_veg_av, & !< average of c_veg538 ghf_av, & !< average of ghf539 lai_av, & !< average of lai540 qsws_liq_av, & !< average of qsws_liq541 qsws_soil_av, & !< average of qsws_soil542 qsws_veg_av, & !< average of qsws_veg543 r_a_av, & !< average of r_a544 r_s_av !< average of r_s568 c_liq_av, & !< average of c_liq 569 c_soil_av, & !< average of c_soil 570 c_veg_av, & !< average of c_veg 571 ghf_av, & !< average of ghf 572 lai_av, & !< average of lai 573 qsws_liq_av, & !< average of qsws_liq 574 qsws_soil_av, & !< average of qsws_soil 575 qsws_veg_av, & !< average of qsws_veg 576 r_a_av, & !< average of r_a 577 r_s_av !< average of r_s 545 578 546 579 … … 548 581 !-- Predefined Land surface classes (vegetation_type) 549 582 CHARACTER(26), DIMENSION(0:18), PARAMETER :: vegetation_type_name = (/ & 550 'user defined ', & ! 0551 'bare soil ', & ! 1552 'crops, mixed farming ', & ! 2553 'short grass ', & ! 3554 'evergreen needleleaf trees', & ! 4555 'deciduous needleleaf trees', & ! 5556 'evergreen broadleaf trees ', & ! 6557 'deciduous broadleaf trees ', & ! 7558 'tall grass ', & ! 8559 'desert ', & ! 9560 'tundra ', & ! 10561 'irrigated crops ', & ! 11562 'semidesert ', & ! 12563 'ice caps and glaciers ', & ! 13564 'bogs and marshes ', & ! 14565 'evergreen shrubs ', & ! 15566 'deciduous shrubs ', & ! 16567 'mixed forest/woodland ', & ! 17568 'interrupted forest ' & ! 18583 'user defined ', & ! 0 584 'bare soil ', & ! 1 585 'crops, mixed farming ', & ! 2 586 'short grass ', & ! 3 587 'evergreen needleleaf trees', & ! 4 588 'deciduous needleleaf trees', & ! 5 589 'evergreen broadleaf trees ', & ! 6 590 'deciduous broadleaf trees ', & ! 7 591 'tall grass ', & ! 8 592 'desert ', & ! 9 593 'tundra ', & ! 10 594 'irrigated crops ', & ! 11 595 'semidesert ', & ! 12 596 'ice caps and glaciers ', & ! 13 597 'bogs and marshes ', & ! 14 598 'evergreen shrubs ', & ! 15 599 'deciduous shrubs ', & ! 16 600 'mixed forest/woodland ', & ! 17 601 'interrupted forest ' & ! 18 569 602 /) 570 603 … … 572 605 !-- Soil model classes (soil_type) 573 606 CHARACTER(12), DIMENSION(0:6), PARAMETER :: soil_type_name = (/ & 574 'user defined', & ! 0575 'coarse ', & ! 1576 'medium ', & ! 2577 'medium-fine ', & ! 3578 'fine ', & ! 4579 'very fine ', & ! 5580 'organic ' & ! 6607 'user defined', & ! 0 608 'coarse ', & ! 1 609 'medium ', & ! 2 610 'medium-fine ', & ! 3 611 'fine ', & ! 4 612 'very fine ', & ! 5 613 'organic ' & ! 6 581 614 /) 582 615 … … 600 633 'artifical turf (sports) ', & ! 14 601 634 'clay (sports) ' & ! 15 602 /) 603 604 605 606 607 635 /) 636 608 637 ! 609 638 !-- Water classes 610 639 CHARACTER(12), DIMENSION(0:5), PARAMETER :: water_type_name = (/ & 611 'user defined', & ! 0612 'lake ', & ! 1613 'river ', & ! 2614 'ocean ', & ! 3615 'pond ', & ! 4616 'fountain ' & ! 5640 'user defined', & ! 0 641 'lake ', & ! 1 642 'river ', & ! 2 643 'ocean ', & ! 3 644 'pond ', & ! 4 645 'fountain ' & ! 5 617 646 /) 618 619 ! 647 ! 648 !-- IDs for vegetation, pavement and water surfaces 649 INTEGER(iwp) :: ind_veg = 0 !< index for vegetation surfaces, used to assess surface-fraction, albedo, etc. 650 INTEGER(iwp) :: ind_pav = 1 !< index for pavement surfaces, used to assess surface-fraction, albedo, etc. 651 INTEGER(iwp) :: ind_wat = 2 !< index for vegetation surfaces, used to assess surface-fraction, albedo, etc. 652 653 ! 620 654 !-- Land surface parameters according to the respective classes (vegetation_type) 621 655 INTEGER(iwp) :: ind_v_rc_min = 0 !< index for r_canopy_min in vegetation_pars 656 INTEGER(iwp) :: ind_v_rc_lai = 1 !< index for LAI in vegetation_pars 657 INTEGER(iwp) :: ind_v_c_veg = 2 !< index for c_veg in vegetation_pars 658 INTEGER(iwp) :: ind_v_gd = 3 !< index for g_d in vegetation_pars 659 INTEGER(iwp) :: ind_v_z0 = 4 !< index for z0 in vegetation_pars 660 INTEGER(iwp) :: ind_v_z0qh = 5 !< index for z0h / z0q in vegetation_pars 661 INTEGER(iwp) :: ind_v_lambda_s = 6 !< index for lambda_s_s in vegetation_pars 662 INTEGER(iwp) :: ind_v_lambda_u = 7 !< index for lambda_s_u in vegetation_pars 663 INTEGER(iwp) :: ind_v_f_sw_in = 8 !< index for f_sw_in in vegetation_pars 664 INTEGER(iwp) :: ind_v_c_surf = 9 !< index for c_surface in vegetation_pars 665 INTEGER(iwp) :: ind_v_at = 10 !< index for albedo_type in vegetation_pars 666 INTEGER(iwp) :: ind_v_emis = 11 !< index for emissivity in vegetation_pars 667 668 INTEGER(iwp) :: ind_w_temp = 0 !< index for temperature in water_pars 669 INTEGER(iwp) :: ind_w_z0 = 1 !< index for z0 in water_pars 670 INTEGER(iwp) :: ind_w_z0h = 2 !< index for z0h in water_pars 671 INTEGER(iwp) :: ind_w_lambda_s = 3 !< index for lambda_s_s in water_pars 672 INTEGER(iwp) :: ind_w_lambda_u = 4 !< index for lambda_s_u in water_pars 673 INTEGER(iwp) :: ind_w_at = 5 !< index for albedo type in water_pars 674 INTEGER(iwp) :: ind_w_emis = 6 !< index for emissivity in water_pars 675 676 INTEGER(iwp) :: ind_p_z0 = 0 !< index for z0 in pavement_pars 677 INTEGER(iwp) :: ind_p_z0h = 1 !< index for z0h in pavement_pars 678 INTEGER(iwp) :: ind_p_at = 2 !< index for albedo type in pavement_pars 679 INTEGER(iwp) :: ind_p_emis = 3 !< index for emissivity in pavement_pars 680 INTEGER(iwp) :: ind_p_lambda_h = 0 !< index for lambda_h in pavement_subsurface_pars 681 INTEGER(iwp) :: ind_p_rho_c = 1 !< index for rho_c in pavement_pars 622 682 ! 623 683 !-- Land surface parameters … … 649 709 !-- level 1 - level 4 according to zs_ref 650 710 REAL(wp), DIMENSION(0:3,1:18), PARAMETER :: root_distribution = RESHAPE( (/ & 651 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 1652 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 2653 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & ! 3654 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp, & ! 4655 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp, & ! 5656 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp, & ! 6657 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp, & ! 7658 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp, & ! 8659 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 9660 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp, & ! 10661 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 11662 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp, & ! 12663 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 13664 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp, & ! 14665 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 15666 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 16667 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 17668 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp & ! 18711 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 1 712 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 2 713 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & ! 3 714 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp, & ! 4 715 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp, & ! 5 716 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp, & ! 6 717 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp, & ! 7 718 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp, & ! 8 719 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 9 720 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp, & ! 10 721 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 11 722 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp, & ! 12 723 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 13 724 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp, & ! 14 725 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 15 726 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 16 727 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 17 728 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp & ! 18 669 729 /), (/ 4, 18 /) ) 670 730 … … 682 742 1.30_wp, 0.400_wp, 1.20_wp, 0.93E-6_wp, 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp & ! 6 683 743 /), (/ 8, 6 /) ) 744 684 745 685 746 ! … … 747 808 ! 748 809 !-- TO BE FILLED 749 !-- Water parameters temperature, z0, z0h, albedo_type, emissivity,750 REAL(wp), DIMENSION(0: 4,1:5), PARAMETER :: water_pars = RESHAPE( (/ &751 283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 1752 283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 2753 283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 3754 283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 4755 283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp & ! 5756 /), (/ 5, 5 /) )810 !-- Water parameters temperature, z0, z0h, albedo_type, emissivity, 811 REAL(wp), DIMENSION(0:6,1:5), PARAMETER :: water_pars = RESHAPE( (/ & 812 283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 1 813 283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 2 814 283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 3 815 283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 4 816 283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp & ! 5 817 /), (/ 7, 5 /) ) 757 818 758 819 SAVE … … 764 825 ! 765 826 !-- Public functions 766 PUBLIC lsm_check_data_output, lsm_check_data_output_pr, & 827 PUBLIC lsm_boundary_condition, lsm_check_data_output, & 828 lsm_check_data_output_pr, & 767 829 lsm_check_parameters, lsm_define_netcdf_grid, lsm_3d_data_averaging,& 768 830 lsm_data_output_2d, lsm_data_output_3d, lsm_energy_balance, & 769 831 lsm_header, lsm_init, lsm_init_arrays, lsm_parin, lsm_soil_model, & 770 lsm_swap_timelevel, lsm_read_restart_data, lsm_ last_actions832 lsm_swap_timelevel, lsm_read_restart_data, lsm_write_restart_data 771 833 ! !vegetat 772 834 !-- Public parameters, constants and initial values … … 781 843 PUBLIC m_soil_h, t_soil_h 782 844 845 INTERFACE lsm_boundary_condition 846 MODULE PROCEDURE lsm_boundary_condition 847 END INTERFACE lsm_boundary_condition 783 848 784 849 INTERFACE lsm_check_data_output … … 842 907 END INTERFACE lsm_read_restart_data 843 908 844 INTERFACE lsm_ last_actions845 MODULE PROCEDURE lsm_ last_actions846 END INTERFACE lsm_ last_actions909 INTERFACE lsm_write_restart_data 910 MODULE PROCEDURE lsm_write_restart_data 911 END INTERFACE lsm_write_restart_data 847 912 848 913 CONTAINS 914 915 916 !------------------------------------------------------------------------------! 917 ! Description: 918 ! ------------ 919 !> Set internal Neumann boundary condition at outer soil grid points 920 !> for temperature and humidity. 921 !------------------------------------------------------------------------------! 922 SUBROUTINE lsm_boundary_condition 923 924 IMPLICIT NONE 925 926 INTEGER(iwp) :: i !< grid index x-direction 927 INTEGER(iwp) :: ioff !< offset index x-direction indicating location of soil grid point 928 INTEGER(iwp) :: j !< grid index y-direction 929 INTEGER(iwp) :: joff !< offset index x-direction indicating location of soil grid point 930 INTEGER(iwp) :: k !< grid index z-direction 931 INTEGER(iwp) :: koff !< offset index x-direction indicating location of soil grid point 932 INTEGER(iwp) :: l !< running index surface-orientation 933 INTEGER(iwp) :: m !< running index surface elements 934 935 koff = surf_lsm_h%koff 936 DO m = 1, surf_lsm_h%ns 937 i = surf_lsm_h%i(m) 938 j = surf_lsm_h%j(m) 939 k = surf_lsm_h%k(m) 940 pt(k+koff,j,i) = pt(k,j,i) 941 ENDDO 942 943 DO l = 0, 3 944 ioff = surf_lsm_v(l)%ioff 945 joff = surf_lsm_v(l)%joff 946 DO m = 1, surf_lsm_v(l)%ns 947 i = surf_lsm_v(l)%i(m) 948 j = surf_lsm_v(l)%j(m) 949 k = surf_lsm_v(l)%k(m) 950 pt(k,j+joff,i+ioff) = pt(k,j,i) 951 ENDDO 952 ENDDO 953 ! 954 !-- In case of humidity, set boundary conditions also for q and vpt. 955 IF ( humidity ) THEN 956 koff = surf_lsm_h%koff 957 DO m = 1, surf_lsm_h%ns 958 i = surf_lsm_h%i(m) 959 j = surf_lsm_h%j(m) 960 k = surf_lsm_h%k(m) 961 q(k+koff,j,i) = q(k,j,i) 962 vpt(k+koff,j,i) = vpt(k,j,i) 963 ENDDO 964 965 DO l = 0, 3 966 ioff = surf_lsm_v(l)%ioff 967 joff = surf_lsm_v(l)%joff 968 DO m = 1, surf_lsm_v(l)%ns 969 i = surf_lsm_v(l)%i(m) 970 j = surf_lsm_v(l)%j(m) 971 k = surf_lsm_v(l)%k(m) 972 q(k,j+joff,i+ioff) = q(k,j,i) 973 vpt(k,j+joff,i+ioff) = vpt(k,j,i) 974 ENDDO 975 ENDDO 976 ENDIF 977 978 END SUBROUTINE lsm_boundary_condition 849 979 850 980 !------------------------------------------------------------------------------! … … 975 1105 976 1106 END SUBROUTINE lsm_check_data_output 1107 977 1108 978 1109 … … 1060 1191 ONLY: bc_pt_b, bc_q_b, constant_flux_layer, message_string, & 1061 1192 most_method 1062 1063 USE radiation_model_mod, & 1064 ONLY: radiation 1065 1193 1066 1194 1067 1195 IMPLICIT NONE … … 1075 1203 TRIM( surface_type ) /= 'pavement' .AND. & 1076 1204 TRIM( surface_type ) /= 'water' .AND. & 1077 TRIM( surface_type ) /= 'netcdf' ) THEN 1078 message_string = 'unknown surface type surface_type = "' //&1205 TRIM( surface_type ) /= 'netcdf' ) THEN 1206 message_string = 'unknown surface type surface_type = "' // & 1079 1207 TRIM( surface_type ) // '"' 1080 1208 CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 ) 1081 1209 ENDIF 1082 1210 … … 1164 1292 ENDIF 1165 1293 ENDIF 1166 1294 1167 1295 IF ( vegetation_type == 1 ) THEN 1168 1296 IF ( vegetation_coverage /= 9999999.9_wp .AND. vegetation_coverage & … … 1260 1388 1261 1389 ENDIF 1390 1391 IF ( TRIM( surface_type ) == 'netcdf' ) THEN 1392 IF ( ANY( water_type_f%var /= water_type_f%fill ) .AND. & 1393 TRIM( most_method ) == 'lookup' ) THEN 1394 WRITE( message_string, * ) 'water-surfaces are not allowed in ' // & 1395 'combination with most_method = ', & 1396 TRIM( most_method ) 1397 CALL message( 'check_parameters', 'PA0999', 2, 2, 0, 6, 0 ) 1398 ENDIF 1399 ! 1400 !-- MS: Some problme here, after calling message everythings stucks at 1401 !-- MPI_FINALIZE call. 1402 IF ( ANY( pavement_type_f%var /= pavement_type_f%fill ) .AND. & 1403 ANY( dz_soil /= 9999999.9_wp ) ) THEN 1404 message_string = 'pavement-surfaces are not allowed in ' // & 1405 'combination with a non-default setting of dz_soil' 1406 CALL message( 'check_parameters', 'PA0999', 2, 2, 0, 6, 0 ) 1407 ENDIF 1408 ENDIF 1262 1409 1263 1410 ! 1264 1411 !-- Temporary message as long as NetCDF input is not available 1265 IF ( TRIM( surface_type ) == 'netcdf' ) THEN1266 message_string = 'surface_type = netcdf '// &1267 'is not supported at the moment'1268 1412 IF ( TRIM( surface_type ) == 'netcdf' .AND. .NOT. input_pids_static ) & 1413 THEN 1414 message_string = 'surface_type = netcdf requires static input file.' 1415 CALL message( 'check_parameters', 'PA0465', 1, 2, 0, 6, 0 ) 1269 1416 ENDIF 1270 1417 … … 1362 1509 DO k = nzb_soil, nzt_soil 1363 1510 IF ( soil_moisture(k) > saturation_moisture ) THEN 1364 message_string = 'soil_moisture must not exceed its saturation' // &1511 message_string = 'soil_moisture must not exceed its saturation' // & 1365 1512 ' value' 1366 1513 CALL message( 'check_parameters', 'PA0458', 1, 2, 0, 6, 0 ) … … 1378 1525 ALLOCATE ( dz_soil_center(nzb_soil:nzt_soil) ) 1379 1526 ALLOCATE ( zs(nzb_soil:nzt_soil+1) ) 1527 1380 1528 1381 1529 zs(nzb_soil) = 0.5_wp * dz_soil(nzb_soil) … … 1428 1576 INTEGER(iwp) :: k !< running index 1429 1577 INTEGER(iwp) :: k_off !< offset to determine index of surface element, seen from atmospheric grid point, for z 1430 INTEGER(iwp) :: k_rad !< index to access radiation array1431 1578 INTEGER(iwp) :: ks !< running index 1432 1579 INTEGER(iwp) :: l !< surface-facing index … … 1454 1601 lambda_soil, & !< Thermal conductivity of the uppermost soil layer (W/m2/K) 1455 1602 lambda_surface, & !< Current value of lambda_surface (W/m2/K) 1456 m_liq_max, & !< maxmimum value of the liq. water reservoir 1457 pt1, & !< potential temperature at first grid level 1458 qv1 !< specific humidity at first grid level 1603 m_liq_max !< maxmimum value of the liq. water reservoir 1459 1604 1460 1605 TYPE(surf_type_lsm), POINTER :: surf_t_surface … … 1472 1617 IF ( horizontal ) THEN 1473 1618 surf => surf_lsm_h 1619 1474 1620 surf_t_surface => t_surface_h 1475 1621 surf_t_surface_p => t_surface_h_p … … 1480 1626 surf_m_soil => m_soil_h 1481 1627 surf_t_soil => t_soil_h 1482 1483 k_off = -11484 j_off = 01485 i_off = 01486 1628 ELSE 1487 1629 surf => surf_lsm_v(l) 1630 1488 1631 surf_t_surface => t_surface_v(l) 1489 1632 surf_t_surface_p => t_surface_v_p(l) … … 1494 1637 surf_m_soil => m_soil_v(l) 1495 1638 surf_t_soil => t_soil_v(l) 1496 1497 k_off = 01498 IF ( l == 0 ) THEN1499 j_off = -11500 i_off = 01501 ELSEIF ( l == 1 ) THEN1502 j_off = 11503 i_off = 01504 ELSEIF ( l == 2 ) THEN1505 j_off = 01506 i_off = -11507 ELSEIF ( l == 3 ) THEN1508 j_off = 01509 i_off = 11510 ENDIF1511 1639 ENDIF 1640 1641 ! 1642 !-- Index offset of surface element point with respect to adjoining 1643 !-- atmospheric grid point 1644 k_off = surf%koff 1645 j_off = surf%joff 1646 i_off = surf%ioff 1512 1647 1513 1648 ! … … 1520 1655 j = surf%j(m) 1521 1656 k = surf%k(m) 1522 !1523 !-- Determine height index for radiation. Note, in clear-sky case radiation1524 !-- arrays have rank 0 in first dimensions, so index must be zero. In case1525 !-- of RRTMG radiation arrays have non-zero rank in first dimension, so that1526 !-- radiation can be obtained at respective height level1527 k_rad = MERGE( 0, k + k_off, radiation_scheme /= 'rrtmg' )1528 1657 1529 1658 ! … … 1589 1718 !-- time step is used here. Note that this formulation is the 1590 1719 !-- equivalent to the ECMWF formulation using drag coefficients 1591 IF ( cloud_physics ) THEN1592 pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)1593 qv1 = q(k,j,i) - ql(k,j,i)1594 ELSEIF ( cloud_droplets ) THEN1595 pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)1596 qv1 = q(k,j,i)1597 ELSE1598 pt1 = pt(k,j,i)1599 IF ( humidity ) THEN1600 qv1 = q(k,j,i)1601 ELSE1602 qv1 = 0.0_wp1603 ENDIF1604 ENDIF1720 ! IF ( cloud_physics ) THEN 1721 ! pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i) 1722 ! qv1 = q(k,j,i) - ql(k,j,i) 1723 ! ELSEIF ( cloud_droplets ) THEN 1724 ! pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i) 1725 ! qv1 = q(k,j,i) 1726 ! ELSE 1727 ! pt1 = pt(k,j,i) 1728 ! IF ( humidity ) THEN 1729 ! qv1 = q(k,j,i) 1730 ! ELSE 1731 ! qv1 = 0.0_wp 1732 ! ENDIF 1733 ! ENDIF 1605 1734 ! 1606 1735 !-- Calculate aerodynamical resistance. For horizontal and vertical … … 1608 1737 !-- can be obtain via parameterization of Mason (2000) / 1609 1738 !-- Krayenhoff and Voogt (2006). 1610 !-- To do: detailed investigation which approach is better! 1739 !-- To do: detailed investigation which approach gives more reliable results! 1740 !-- Please note, in case of very small friction velocity, e.g. in little 1741 !-- holes, the resistance can become negative. For this reason, limit r_a 1742 !-- to positive values. 1611 1743 IF ( horizontal .OR. .NOT. aero_resist_kray ) THEN 1612 surf%r_a(m) = ( pt1 - surf_lsm_h%pt_surface(m) ) /&1613 ( surf%ts(m) * surf%us(m) + 1.0E-20_wp)1744 surf%r_a(m) = ABS( ( surf%pt1(m) - surf%pt_surface(m) ) / & 1745 ( surf%ts(m) * surf%us(m) + 1.0E-20_wp ) ) 1614 1746 ELSE 1615 1747 surf%r_a(m) = 1.0_wp / ( 11.8_wp + 4.2_wp * & 1616 SQRT( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + & 1617 ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + & 1618 ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2 ) & 1748 SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + & 1749 ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + & 1750 ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2, & 1751 0.01_wp ) ) & 1619 1752 ) 1620 1753 ENDIF 1621 1754 ! 1622 1755 !-- Make sure that the resistance does not drop to zero for neutral 1623 !-- stratification 1624 IF ( ABS( surf%r_a(m)) < 1.0_wp ) surf%r_a(m) = 1.0_wp1756 !-- stratification. 1757 IF ( surf%r_a(m) < 1.0_wp ) surf%r_a(m) = 1.0_wp 1625 1758 ! 1626 1759 !-- Second step: calculate canopy resistance r_canopy … … 1629 1762 !-- f1: correction for incoming shortwave radiation (stomata close at 1630 1763 !-- night) 1631 f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k_rad,j+j_off,i+i_off)&1632 + 0.05_wp ) / (0.81_wp * (0.004_wp&1633 * rad_sw_in(k_rad,j+j_off,i+i_off)+ 1.0_wp)) )1764 f1 = MIN( 1.0_wp, ( 0.004_wp * surf%rad_sw_in(m) + 0.05_wp ) / & 1765 (0.81_wp * (0.004_wp * surf%rad_sw_in(m) & 1766 + 1.0_wp)) ) 1634 1767 ! 1635 1768 !-- f2: correction for soil moisture availability to plants (the … … 1665 1798 ! 1666 1799 !-- Calculate vapour pressure 1667 e = qv1 * surface_pressure / ( qv1+ 0.622_wp )1800 e = surf%qv1(m) * surface_pressure / ( surf%qv1(m) + 0.622_wp ) 1668 1801 f3 = EXP ( - surf%g_d(m) * (e_s - e) ) 1669 1802 ELSE … … 1682 1815 1683 1816 1684 f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min ) 1685 /( surf%m_fc(nzb_soil,m) - m_min )1817 f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min ) / & 1818 ( surf%m_fc(nzb_soil,m) - m_min ) 1686 1819 f2 = MAX( f2, 1.0E-20_wp ) 1687 1820 f2 = MIN( f2, 1.0_wp ) … … 1698 1831 IF ( surf%pavement_surface(m) ) THEN 1699 1832 m_liq_max = m_max_depth * 5.0_wp 1700 surf%c_liq(m) = MIN( 1.0_wp, ( surf_m_liq%var_1d(m) / m_liq_max)**0.67 )1833 surf%c_liq(m) = MIN( 1.0_wp, ( surf_m_liq%var_1d(m) / m_liq_max)**0.67 ) 1701 1834 ELSE 1702 1835 m_liq_max = m_max_depth * ( surf%c_veg(m) * surf%lai(m) & … … 1710 1843 !-- In case of dewfall, set evapotranspiration to zero 1711 1844 !-- All super-saturated water is then removed from the air 1712 IF ( humidity .AND. q_s <= qv1) THEN1845 IF ( humidity .AND. q_s <= surf%qv1(m) ) THEN 1713 1846 surf%r_canopy(m) = 0.0_wp 1714 1847 surf%r_soil(m) = 0.0_wp … … 1747 1880 dq_s_dt = 0.622_wp * e_s_dt / ( surface_pressure - e_s_dt ) 1748 1881 ! 1749 !-- Add LW up so that it can be removed in prognostic equation 1750 surf%rad_net_l(m) = rad_net(j,i) + rad_lw_out(nzb,j,i) 1882 !-- Calculate net radiation radiation without longwave outgoing flux because 1883 !-- it has a dependency on surface temperature and thus enters the prognostic 1884 !-- equations directly 1885 surf%rad_net_l(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m) & 1886 + surf%rad_lw_in(m) 1751 1887 ! 1752 1888 !-- Calculate new skin temperature … … 1754 1890 ! 1755 1891 !-- Numerator of the prognostic equation 1756 coef_1 = surf%rad_net_l(m) + rad_lw_out_change_0(j,i)&1757 * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)&1758 + f_shf * pt1 + f_qsws * ( qv1 - q_s&1892 coef_1 = surf%rad_net_l(m) + surf%rad_lw_out_change_0(m) & 1893 * surf_t_surface%var_1d(m) - surf%rad_lw_out(m) & 1894 + f_shf * surf%pt1(m) + f_qsws * ( surf%qv1(m) - q_s & 1759 1895 + dq_s_dt * surf_t_surface%var_1d(m) ) + lambda_surface & 1760 1896 * surf_t_soil%var_2d(nzb_soil,m) 1897 1761 1898 ! 1762 1899 !-- Denominator of the prognostic equation 1763 coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt&1900 coef_2 = surf%rad_lw_out_change_0(m) + f_qsws * dq_s_dt & 1764 1901 + lambda_surface + f_shf / exn 1765 1902 ELSE 1766 1903 ! 1767 1904 !-- Numerator of the prognostic equation 1768 coef_1 = surf%rad_net_l(m) + rad_lw_out_change_0(j,i)&1769 * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)&1770 + f_shf * pt1 + lambda_surface&1905 coef_1 = surf%rad_net_l(m) + surf%rad_lw_out_change_0(m) & 1906 * surf_t_surface%var_1d(m) - surf%rad_lw_out(m) & 1907 + f_shf * surf%pt1(m) + lambda_surface & 1771 1908 * surf_t_soil%var_2d(nzb_soil,m) 1772 1909 ! 1773 1910 !-- Denominator of the prognostic equation 1774 coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn1911 coef_2 = surf%rad_lw_out_change_0(m) + lambda_surface + f_shf / exn 1775 1912 1776 1913 ENDIF 1914 1777 1915 1778 1916 tend = 0.0_wp … … 1819 1957 !-- computing time for the radiation code). 1820 1958 IF ( ABS( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) ) & 1821 > 0.2_wp .AND. unscheduled_radiation_calls ) THEN 1959 > 0.2_wp .AND. & 1960 unscheduled_radiation_calls ) THEN 1822 1961 force_radiation_call_l = .TRUE. 1823 1962 ENDIF 1824 1963 1825 pt(k+k_off,j+j_off,i+i_off) = surf_t_surface_p%var_1d(m) / exn !is actually no air temperature 1964 1965 ! IF ( surf_t_surface_p%var_1d(m) < 270.0_wp .OR. surf_t_surface_p%var_1d(m) > 350.0_wp ) THEN 1966 ! PRINT*, "myid: ", myid 1967 ! PRINT*, "m: ", m 1968 ! PRINT*, "i,j,k: ", i, j, k 1969 ! PRINT*, "pt_p: ", surf_t_surface_p%var_1d(m) 1970 ! PRINT*, "f_shf: ", f_shf 1971 ! PRINT*, "f_qsws: ", f_qsws 1972 ! ENDIF 1973 1974 1975 ! pt(k+k_off,j+j_off,i+i_off) = surf_t_surface_p%var_1d(m) / exn !is actually no air temperature 1826 1976 surf%pt_surface(m) = surf_t_surface_p%var_1d(m) / exn 1827 1977 … … 1829 1979 !-- Calculate fluxes 1830 1980 surf%rad_net_l(m) = surf%rad_net_l(m) + & 1831 rad_lw_out_change_0(j,i)&1832 * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)&1833 - rad_lw_out_change_0(j,i) * surf_t_surface_p%var_1d(m)1834 1835 rad_net(j,i) = surf%rad_net_l(m)1836 rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) + rad_lw_out_change_0(j,i) *&1981 surf%rad_lw_out_change_0(m) & 1982 * surf_t_surface%var_1d(m) - surf%rad_lw_out(m) & 1983 - surf%rad_lw_out_change_0(m) * surf_t_surface_p%var_1d(m) 1984 1985 surf%rad_net(m) = surf%rad_net_l(m) 1986 surf%rad_lw_out(m) = surf%rad_lw_out(m) + surf%rad_lw_out_change_0(m) * & 1837 1987 ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) ) 1838 1988 … … 1840 1990 - surf_t_soil%var_2d(nzb_soil,m) ) 1841 1991 1842 surf%shf(m) = - f_shf * ( pt1 - surf%pt_surface(m) ) / cp 1843 1992 surf%shf(m) = - f_shf * ( surf%pt1(m) - surf%pt_surface(m) ) / cp 1844 1993 1845 1994 IF ( humidity ) THEN 1846 surf%qsws(m) = - f_qsws * ( qv1 - q_s + dq_s_dt&1995 surf%qsws(m) = - f_qsws * ( surf%qv1(m) - q_s + dq_s_dt & 1847 1996 * surf_t_surface%var_1d(m) - dq_s_dt * & 1848 1997 surf_t_surface_p%var_1d(m) ) 1849 1998 1850 surf%qsws_veg(m) = - f_qsws_veg * ( qv1 - q_s&1999 surf%qsws_veg(m) = - f_qsws_veg * ( surf%qv1(m) - q_s & 1851 2000 + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt & 1852 2001 * surf_t_surface_p%var_1d(m) ) 1853 2002 1854 surf%qsws_soil(m) = - f_qsws_soil * ( qv1 - q_s&2003 surf%qsws_soil(m) = - f_qsws_soil * ( surf%qv1(m) - q_s & 1855 2004 + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt & 1856 2005 * surf_t_surface_p%var_1d(m) ) 1857 2006 1858 surf%qsws_liq(m) = - f_qsws_liq * ( qv1 - q_s&2007 surf%qsws_liq(m) = - f_qsws_liq * ( surf%qv1(m) - q_s & 1859 2008 + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt & 1860 2009 * surf_t_surface_p%var_1d(m) ) … … 1866 2015 surf%r_s(m) = 1.0E10_wp 1867 2016 ELSE 1868 surf%r_s(m) = - rho_lv * ( qv1 - q_s + dq_s_dt&2017 surf%r_s(m) = - rho_lv * ( surf%qv1(m) - q_s + dq_s_dt & 1869 2018 * surf_t_surface%var_1d(m) - dq_s_dt * & 1870 2019 surf_t_surface_p%var_1d(m) ) / & … … 1954 2103 intermediate_timestep_count_max ) THEN 1955 2104 surf_tm_liq_m%var_1d(m) = -9.5625_wp * tend + & 1956 2105 5.3125_wp * surf_tm_liq_m%var_1d(m) 1957 2106 ENDIF 1958 2107 ENDIF … … 2015 2164 q_s = 0.622_wp * e_s / ( surface_pressure - e_s ) 2016 2165 2017 resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) )2166 resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) + 1E-5_wp ) 2018 2167 2019 2168 ! … … 2028 2177 q(k,j,i) 2029 2178 ENDIF 2030 2031 2179 ! 2032 2180 !-- Update virtual potential temperature … … 2098 2246 ENDIF 2099 2247 2100 WRITE( io, 4 ) TRIM( vegetation_type_name(vegetation_type) ), & 2101 TRIM (soil_type_name(soil_type) ) 2102 WRITE( io, 5 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ), & 2248 IF ( vegetation_type_f%from_file ) THEN 2249 WRITE( io, 5 ) 2250 ELSE 2251 WRITE( io, 4 ) TRIM( vegetation_type_name(vegetation_type) ), & 2252 TRIM (soil_type_name(soil_type) ) 2253 ENDIF 2254 WRITE( io, 6 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ), & 2103 2255 TRIM( m_soil_chr ), TRIM( roots_chr ), & 2104 2256 TRIM( vertical_index_chr ) … … 2111 2263 4 FORMAT (' --> Land surface type : ',A,/ & 2112 2264 ' --> Soil porosity type : ',A) 2113 5 FORMAT (/' Initial soil temperature and moisture profile:'// & 2265 5 FORMAT (' --> Land surface type : read from file' / & 2266 ' --> Soil porosity type : read from file' ) 2267 6 FORMAT (/' Initial soil temperature and moisture profile:'// & 2114 2268 ' Height: ',A,' m'/ & 2115 2269 ' Temperature: ',A,' K'/ & … … 2117 2271 ' Root fraction: ',A,' '/ & 2118 2272 ' Grid point: ',A) 2273 2119 2274 2120 2275 END SUBROUTINE lsm_header … … 2133 2288 IMPLICIT NONE 2134 2289 2135 INTEGER(iwp) :: i !< running index 2136 INTEGER(iwp) :: i_off !< index offset of surface element, seen from atmospheric grid point 2137 INTEGER(iwp) :: j !< running index 2138 INTEGER(iwp) :: j_off !< index offset of surface element, seen from atmospheric grid point 2139 INTEGER(iwp) :: k !< running index 2140 INTEGER(iwp) :: kn !< running index 2141 INTEGER(iwp) :: ko !< running index 2142 INTEGER(iwp) :: kroot !< running index 2143 INTEGER(iwp) :: kzs !< running index 2144 INTEGER(iwp) :: l !< running index surface facing 2145 INTEGER(iwp) :: m !< running index 2290 INTEGER(iwp) :: i !< running index 2291 INTEGER(iwp) :: i_off !< index offset of surface element, seen from atmospheric grid point 2292 INTEGER(iwp) :: j !< running index 2293 INTEGER(iwp) :: j_off !< index offset of surface element, seen from atmospheric grid point 2294 INTEGER(iwp) :: k !< running index 2295 INTEGER(iwp) :: kn !< running index 2296 INTEGER(iwp) :: ko !< running index 2297 INTEGER(iwp) :: kroot !< running index 2298 INTEGER(iwp) :: kzs !< running index 2299 INTEGER(iwp) :: l !< running index surface facing 2300 INTEGER(iwp) :: m !< running index 2301 INTEGER(iwp) :: st !< soil-type index 2146 2302 INTEGER(iwp) :: n_soil_layers_total !< temperature variable, stores the total number of soil layers + 4 2147 2148 2149 REAL(wp) :: pt1 !< potential temperature at first grid level 2303 INTEGER(iwp) :: n_surf !< number of surface types of given surface element 2150 2304 2151 2305 REAL(wp), DIMENSION(:), ALLOCATABLE :: bound, bound_root_fr !< temporary arrays for storing index bounds … … 2157 2311 !-- If no cloud physics is used, rho_surface has not been calculated before 2158 2312 IF ( .NOT. cloud_physics ) THEN 2159 rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn ) 2313 CALL calc_mean_profile( pt, 4 ) 2314 rho_surface = surface_pressure * 100.0_wp / ( r_d * hom(nzb+1,1,4,0) * exn ) 2160 2315 ENDIF 2161 2316 … … 2206 2361 surf_lsm_v(l)%r_soil = 0.0_wp 2207 2362 ENDDO 2208 2363 2364 ! 2365 !-- Set initial values for prognostic soil quantities 2366 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 2367 t_soil_h%var_2d = 0.0_wp 2368 m_soil_h%var_2d = 0.0_wp 2369 m_liq_h%var_1d = 0.0_wp 2370 2371 DO l = 0, 3 2372 t_soil_v(l)%var_2d = 0.0_wp 2373 m_soil_v(l)%var_2d = 0.0_wp 2374 m_liq_v(l)%var_1d = 0.0_wp 2375 ENDDO 2376 ENDIF 2209 2377 ! 2210 2378 !-- Allocate 3D soil model arrays … … 2261 2429 ENDIF 2262 2430 ENDDO 2263 2431 ! 2432 !-- Allocate albedo type and emissivity for vegetation, water and pavement 2433 !-- fraction. 2434 !-- Set default values at each surface element. 2435 ALLOCATE ( surf_lsm_h%albedo_type(0:2,1:surf_lsm_h%ns) ) 2436 ALLOCATE ( surf_lsm_h%emissivity(0:2,1:surf_lsm_h%ns) ) 2437 surf_lsm_h%albedo_type = 0 2438 surf_lsm_h%emissivity = emissivity 2439 DO l = 0, 3 2440 ALLOCATE ( surf_lsm_v(l)%albedo_type(0:2,1:surf_lsm_v(l)%ns) ) 2441 ALLOCATE ( surf_lsm_v(l)%emissivity(0:2,1:surf_lsm_v(l)%ns) ) 2442 surf_lsm_v(l)%albedo_type = 0 2443 surf_lsm_v(l)%emissivity = emissivity 2444 ENDDO 2445 ! 2446 !-- Allocate arrays for relative surface fraction. 2447 !-- 0 - vegetation fraction, 2 - water fraction, 1 - pavement fraction 2448 ALLOCATE( surf_lsm_h%frac(0:2,1:surf_lsm_h%ns) ) 2449 surf_lsm_h%frac = 0.0_wp 2450 DO l = 0, 3 2451 ALLOCATE( surf_lsm_v(l)%frac(0:2,1:surf_lsm_v(l)%ns) ) 2452 surf_lsm_v(l)%frac = 0.0_wp 2453 ENDDO 2454 ! 2455 !-- For vertical walls only - allocate special flag indicating if any building is on 2456 !-- top of any natural surfaces. Used for initialization only. 2457 DO l = 0, 3 2458 ALLOCATE( surf_lsm_v(l)%building_covered(1:surf_lsm_v(l)%ns) ) 2459 ENDDO 2264 2460 ! 2265 2461 !-- Set flag parameter for the prescribed surface type depending on user 2266 !-- input. 2267 DO m = 1, surf_lsm_h%ns 2268 2269 SELECT CASE ( TRIM( surface_type ) ) 2462 !-- input. Set surface fraction to 1 for the respective type. 2463 SELECT CASE ( TRIM( surface_type ) ) 2270 2464 2271 2465 CASE ( 'vegetation' ) 2272 2466 2273 surf_lsm_h%vegetation_surface(m) = .TRUE. 2467 surf_lsm_h%vegetation_surface = .TRUE. 2468 surf_lsm_h%frac(ind_veg,:) = 1.0_wp 2469 DO l = 0, 3 2470 surf_lsm_v(l)%vegetation_surface = .TRUE. 2471 surf_lsm_v(l)%frac(ind_veg,:) = 1.0_wp 2472 ENDDO 2274 2473 2275 2474 CASE ( 'water' ) 2276 2475 2277 surf_lsm_h%water_surface(m) = .TRUE. 2278 2279 CASE ( 'pavement' ) 2476 surf_lsm_h%water_surface = .TRUE. 2477 surf_lsm_h%frac(ind_wat,:) = 1.0_wp 2478 ! 2479 !-- Note, vertical water surface does not really make sense. 2480 DO l = 0, 3 2481 surf_lsm_v(l)%water_surface = .TRUE. 2482 surf_lsm_v(l)%frac(ind_wat,:) = 1.0_wp 2483 ENDDO 2484 2485 CASE ( 'pavement' ) 2280 2486 2281 surf_lsm_h%pavement_surface(m) = .TRUE. 2282 2283 END SELECT 2284 2285 ENDDO 2286 2287 DO l = 0, 3 2288 SELECT CASE ( TRIM( surface_type ) ) 2289 2290 CASE ( 'vegetation' ) 2291 2292 surf_lsm_v(l)%vegetation_surface = .TRUE. 2293 2294 CASE ( 'water' ) 2295 2296 surf_lsm_v(l)%water_surface = .TRUE. 2297 2298 CASE ( 'pavement' ) 2299 2300 surf_lsm_h%pavement_surface = .TRUE. 2301 2302 END SELECT 2303 2304 ENDDO 2305 ! 2306 !-- Initialize standard soil types. It is possible to overwrite each 2307 !-- parameter by setting the respecticy NAMELIST variable to a 2308 !-- value /= 9999999.9. 2487 surf_lsm_h%pavement_surface = .TRUE. 2488 surf_lsm_h%frac(ind_pav,:) = 1.0_wp 2489 DO l = 0, 3 2490 surf_lsm_v(l)%pavement_surface = .TRUE. 2491 surf_lsm_v(l)%frac(ind_pav,:) = 1.0_wp 2492 ENDDO 2493 2494 CASE ( 'netcdf' ) 2495 2496 DO m = 1, surf_lsm_h%ns 2497 i = surf_lsm_h%i(m) 2498 j = surf_lsm_h%j(m) 2499 IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill ) & 2500 surf_lsm_h%vegetation_surface(m) = .TRUE. 2501 IF ( pavement_type_f%var(j,i) /= pavement_type_f%fill ) & 2502 surf_lsm_h%pavement_surface(m) = .TRUE. 2503 IF ( water_type_f%var(j,i) /= water_type_f%fill ) & 2504 surf_lsm_h%water_surface(m) = .TRUE. 2505 ENDDO 2506 DO l = 0, 3 2507 DO m = 1, surf_lsm_v(l)%ns 2508 ! 2509 !-- Only for vertical surfaces. Check if natural walls at reference 2510 !-- grid point are covered by any building. This case, problems 2511 !-- with initialization will aris if index offsets are used. 2512 !-- In order to deal with this, set special flag. 2513 surf_lsm_v(l)%building_covered(m) = .FALSE. 2514 IF ( building_type_f%from_file ) THEN 2515 i = surf_lsm_v(l)%i(m) + surf_lsm_v(l)%ioff 2516 j = surf_lsm_v(l)%j(m) + surf_lsm_v(l)%joff 2517 IF ( building_type_f%var(j,i) /= 0 ) & 2518 surf_lsm_v(l)%building_covered(m) = .TRUE. 2519 ENDIF 2520 ! 2521 !-- Normally proceed with setting surface types. 2522 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 2523 surf_lsm_v(l)%building_covered(m) ) 2524 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 2525 surf_lsm_v(l)%building_covered(m) ) 2526 IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill ) & 2527 surf_lsm_v(l)%vegetation_surface(m) = .TRUE. 2528 IF ( pavement_type_f%var(j,i) /= pavement_type_f%fill ) & 2529 surf_lsm_v(l)%pavement_surface(m) = .TRUE. 2530 IF ( water_type_f%var(j,i) /= water_type_f%fill ) & 2531 surf_lsm_v(l)%water_surface(m) = .TRUE. 2532 ENDDO 2533 ENDDO 2534 2535 END SELECT 2536 ! 2537 !-- In case of netcdf input file, further initialize surface fractions. 2538 !-- At the moment only 1 surface is given at a location, so that the fraction 2539 !-- is either 0 or 1. This will be revised later. 2540 IF ( input_pids_static ) THEN 2541 DO m = 1, surf_lsm_h%ns 2542 i = surf_lsm_h%i(m) 2543 j = surf_lsm_h%j(m) 2544 ! 2545 !-- 0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction 2546 surf_lsm_h%frac(ind_veg,m) = surface_fraction_f%frac(ind_veg,j,i) 2547 surf_lsm_h%frac(ind_pav,m) = surface_fraction_f%frac(ind_pav,j,i) 2548 surf_lsm_h%frac(ind_wat,m) = surface_fraction_f%frac(ind_wat,j,i) 2549 2550 ENDDO 2551 DO l = 0, 3 2552 DO m = 1, surf_lsm_v(l)%ns 2553 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 2554 surf_lsm_v(l)%building_covered(m) ) 2555 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 2556 surf_lsm_v(l)%building_covered(m) ) 2557 ! 2558 !-- 0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction 2559 surf_lsm_v(l)%frac(ind_veg,m) = surface_fraction_f%frac(ind_veg,j,i) 2560 surf_lsm_v(l)%frac(ind_pav,m) = surface_fraction_f%frac(ind_pav,j,i) 2561 surf_lsm_v(l)%frac(ind_wat,m) = surface_fraction_f%frac(ind_wat,j,i) 2562 2563 ENDDO 2564 ENDDO 2565 ENDIF 2566 ! 2567 !-- Level 1, initialization of soil parameters. 2568 !-- It is possible to overwrite each parameter by setting the respecticy 2569 !-- NAMELIST variable to a value /= 9999999.9. 2309 2570 IF ( soil_type /= 0 ) THEN 2310 2571 … … 2339 2600 IF ( residual_moisture == 9999999.9_wp ) THEN 2340 2601 residual_moisture = soil_pars(7,soil_type) 2341 ENDIF2342 2343 ENDIF2344 2345 !2346 !-- Check whether parameters from the lookup tables are to be used2347 IF ( vegetation_type /= 0 ) THEN2348 2349 IF ( min_canopy_resistance == 9999999.9_wp ) THEN2350 min_canopy_resistance = vegetation_pars(0,vegetation_type)2351 ENDIF2352 2353 IF ( leaf_area_index == 9999999.9_wp ) THEN2354 leaf_area_index = vegetation_pars(1,vegetation_type)2355 ENDIF2356 2357 IF ( vegetation_coverage == 9999999.9_wp ) THEN2358 vegetation_coverage = vegetation_pars(2,vegetation_type)2359 ENDIF2360 2361 IF ( canopy_resistance_coefficient == 9999999.9_wp ) THEN2362 canopy_resistance_coefficient= vegetation_pars(3,vegetation_type)2363 ENDIF2364 2365 IF ( z0_vegetation == 9999999.9_wp ) THEN2366 z0_vegetation = vegetation_pars(4,vegetation_type)2367 ENDIF2368 2369 IF ( z0h_vegetation == 9999999.9_wp ) THEN2370 z0h_vegetation = vegetation_pars(5,vegetation_type)2371 ENDIF2372 2373 IF ( lambda_surface_stable == 9999999.9_wp ) THEN2374 lambda_surface_stable = vegetation_pars(6,vegetation_type)2375 ENDIF2376 2377 IF ( lambda_surface_unstable == 9999999.9_wp ) THEN2378 lambda_surface_unstable = vegetation_pars(7,vegetation_type)2379 ENDIF2380 2381 IF ( f_shortwave_incoming == 9999999.9_wp ) THEN2382 f_shortwave_incoming = vegetation_pars(8,vegetation_type)2383 ENDIF2384 2385 IF ( c_surface == 9999999.9_wp ) THEN2386 c_surface = vegetation_pars(9,vegetation_type)2387 ENDIF2388 2389 IF ( albedo_type == 9999999 .AND. albedo == 9999999.9_wp ) THEN2390 albedo_type = INT(vegetation_pars(10,vegetation_type))2391 ENDIF2392 2393 IF ( emissivity == 9999999.9_wp ) THEN2394 emissivity = vegetation_pars(11,vegetation_type)2395 ENDIF2396 2397 ENDIF2398 2399 IF ( water_type /= 0 ) THEN2400 2401 IF ( water_temperature == 9999999.9_wp ) THEN2402 water_temperature = water_pars(0,water_type)2403 ENDIF2404 2405 IF ( z0_water == 9999999.9_wp ) THEN2406 z0_water = water_pars(1,water_type)2407 ENDIF2408 2409 IF ( z0h_water == 9999999.9_wp ) THEN2410 z0h_water = water_pars(2,water_type)2411 ENDIF2412 2413 IF ( albedo_type == 9999999 .AND. albedo == 9999999.9_wp ) THEN2414 albedo_type = INT(water_pars(3,water_type))2415 ENDIF2416 2417 IF ( emissivity == 9999999.9_wp ) THEN2418 emissivity = water_pars(4,water_type)2419 ENDIF2420 2421 ENDIF2422 2423 IF ( pavement_type /= 0 ) THEN2424 2425 !2426 !-- When a pavement_type is used, overwrite a possible setting of2427 !-- the pavement depth as it is already defined by the pavement type2428 pavement_depth_level = 02429 2430 IF ( z0_pavement == 9999999.9_wp ) THEN2431 z0_pavement = pavement_pars(0,pavement_type)2432 ENDIF2433 2434 IF ( z0h_pavement == 9999999.9_wp ) THEN2435 z0h_pavement = pavement_pars(1,pavement_type)2436 ENDIF2437 2438 IF ( pavement_heat_conduct == 9999999.9_wp ) THEN2439 pavement_heat_conduct = pavement_subsurface_pars_1(0,pavement_type)2440 ENDIF2441 2442 IF ( pavement_heat_capacity == 9999999.9_wp ) THEN2443 pavement_heat_capacity = pavement_subsurface_pars_2(0,pavement_type)2444 ENDIF2445 2446 IF ( albedo_type == 9999999 .AND. albedo == 9999999.9_wp ) THEN2447 albedo_type = INT(pavement_pars(2,pavement_type))2448 ENDIF2449 2450 IF ( emissivity == 9999999.9_wp ) THEN2451 emissivity = pavement_pars(3,pavement_type)2452 ENDIF2453 2454 !2455 !-- If the depth level of the pavement is not set, determine it from2456 !-- lookup table.2457 IF ( pavement_depth_level == 0 ) THEN2458 DO k = nzb_soil, nzt_soil2459 IF ( pavement_subsurface_pars_1(k,pavement_type) == 9999999.9_wp &2460 .OR. pavement_subsurface_pars_2(k,pavement_type) &2461 == 9999999.9_wp) THEN2462 nzt_pavement = k-12463 EXIT2464 ENDIF2465 ENDDO2466 ELSE2467 nzt_pavement = pavement_depth_level2468 2602 ENDIF 2469 2603 … … 2494 2628 surf_lsm_v(l)%r_soil_min = min_soil_resistance 2495 2629 ENDDO 2496 2630 ! 2631 !-- Level 2, initialization of soil parameters via soil_type read from file. 2632 !-- Soil parameters are initialized for each (y,x)-grid point 2633 !-- individually using default paramter settings according to the given 2634 !-- soil type. 2635 IF ( soil_type_f%from_file ) THEN 2636 ! 2637 !-- Level of detail = 1, i.e. a homogeneous soil distribution along the 2638 !-- vertical dimension is assumed. 2639 IF ( soil_type_f%lod == 1 ) THEN 2640 ! 2641 !-- Horizontal surfaces 2642 DO m = 1, surf_lsm_h%ns 2643 i = surf_lsm_h%i(m) 2644 j = surf_lsm_h%j(m) 2645 2646 st = soil_type_f%var_2d(j,i) 2647 IF ( st /= soil_type_f%fill ) THEN 2648 surf_lsm_h%alpha_vg(:,m) = soil_pars(0,st) 2649 surf_lsm_h%l_vg(:,m) = soil_pars(1,st) 2650 surf_lsm_h%n_vg(:,m) = soil_pars(2,st) 2651 surf_lsm_h%gamma_w_sat(:,m) = soil_pars(3,st) 2652 surf_lsm_h%m_sat(:,m) = soil_pars(4,st) 2653 surf_lsm_h%m_fc(:,m) = soil_pars(5,st) 2654 surf_lsm_h%m_wilt(:,m) = soil_pars(6,st) 2655 surf_lsm_h%m_res(:,m) = soil_pars(7,st) 2656 ENDIF 2657 ENDDO 2658 ! 2659 !-- Vertical surfaces ( assumes the soil type given at respective (x,y) 2660 DO l = 0, 3 2661 DO m = 1, surf_lsm_v(l)%ns 2662 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 2663 surf_lsm_v(l)%building_covered(m) ) 2664 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 2665 surf_lsm_v(l)%building_covered(m) ) 2666 2667 st = soil_type_f%var_2d(j,i) 2668 IF ( st /= soil_type_f%fill ) THEN 2669 surf_lsm_v(l)%alpha_vg(:,m) = soil_pars(0,st) 2670 surf_lsm_v(l)%l_vg(:,m) = soil_pars(1,st) 2671 surf_lsm_v(l)%n_vg(:,m) = soil_pars(2,st) 2672 surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars(3,st) 2673 surf_lsm_v(l)%m_sat(:,m) = soil_pars(4,st) 2674 surf_lsm_v(l)%m_fc(:,m) = soil_pars(5,st) 2675 surf_lsm_v(l)%m_wilt(:,m) = soil_pars(6,st) 2676 surf_lsm_v(l)%m_res(:,m) = soil_pars(7,st) 2677 ENDIF 2678 ENDDO 2679 ENDDO 2680 ! 2681 !-- Level of detail = 2, i.e. soil type and thus the soil parameters 2682 !-- can be heterogeneous along the vertical dimension. 2683 ELSE 2684 ! 2685 !-- Horizontal surfaces 2686 DO m = 1, surf_lsm_h%ns 2687 i = surf_lsm_h%i(m) 2688 j = surf_lsm_h%j(m) 2689 2690 DO k = nzb_soil, nzt_soil 2691 st = soil_type_f%var_3d(k,j,i) 2692 IF ( st /= soil_type_f%fill ) THEN 2693 surf_lsm_h%alpha_vg(k,m) = soil_pars(0,st) 2694 surf_lsm_h%l_vg(k,m) = soil_pars(1,st) 2695 surf_lsm_h%n_vg(k,m) = soil_pars(2,st) 2696 surf_lsm_h%gamma_w_sat(k,m) = soil_pars(3,st) 2697 surf_lsm_h%m_sat(k,m) = soil_pars(4,st) 2698 surf_lsm_h%m_fc(k,m) = soil_pars(5,st) 2699 surf_lsm_h%m_wilt(k,m) = soil_pars(6,st) 2700 surf_lsm_h%m_res(k,m) = soil_pars(7,st) 2701 ENDIF 2702 ENDDO 2703 ENDDO 2704 ! 2705 !-- Vertical surfaces ( assumes the soil type given at respective (x,y) 2706 DO l = 0, 3 2707 DO m = 1, surf_lsm_v(l)%ns 2708 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 2709 surf_lsm_v(l)%building_covered(m) ) 2710 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 2711 surf_lsm_v(l)%building_covered(m) ) 2712 2713 DO k = nzb_soil, nzt_soil 2714 st = soil_type_f%var_3d(k,j,i) 2715 IF ( st /= soil_type_f%fill ) THEN 2716 surf_lsm_v(l)%alpha_vg(k,m) = soil_pars(0,st) 2717 surf_lsm_v(l)%l_vg(k,m) = soil_pars(1,st) 2718 surf_lsm_v(l)%n_vg(k,m) = soil_pars(2,st) 2719 surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars(3,st) 2720 surf_lsm_v(l)%m_sat(k,m) = soil_pars(4,st) 2721 surf_lsm_v(l)%m_fc(k,m) = soil_pars(5,st) 2722 surf_lsm_v(l)%m_wilt(k,m) = soil_pars(6,st) 2723 surf_lsm_v(l)%m_res(k,m) = soil_pars(7,st) 2724 ENDIF 2725 ENDDO 2726 ENDDO 2727 ENDDO 2728 ENDIF 2729 ENDIF 2730 ! 2731 !-- Level 3, initialization of single soil parameters at single z,x,y 2732 !-- position via soil_pars read from file. 2733 IF ( soil_pars_f%from_file ) THEN 2734 ! 2735 !-- Level of detail = 1, i.e. a homogeneous vertical distribution of soil 2736 !-- parameters is assumed. 2737 !-- Horizontal surfaces 2738 IF ( soil_pars_f%lod == 1 ) THEN 2739 ! 2740 !-- Horizontal surfaces 2741 DO m = 1, surf_lsm_h%ns 2742 i = surf_lsm_h%i(m) 2743 j = surf_lsm_h%j(m) 2744 2745 IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill ) & 2746 surf_lsm_h%alpha_vg(:,m) = soil_pars_f%pars_xy(0,j,i) 2747 IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill ) & 2748 surf_lsm_h%l_vg(:,m) = soil_pars_f%pars_xy(1,j,i) 2749 IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill ) & 2750 surf_lsm_h%n_vg(:,m) = soil_pars_f%pars_xy(2,j,i) 2751 IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill ) & 2752 surf_lsm_h%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i) 2753 IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill ) & 2754 surf_lsm_h%m_sat(:,m) = soil_pars_f%pars_xy(4,j,i) 2755 IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill ) & 2756 surf_lsm_h%m_fc(:,m) = soil_pars_f%pars_xy(5,j,i) 2757 IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill ) & 2758 surf_lsm_h%m_wilt(:,m) = soil_pars_f%pars_xy(6,j,i) 2759 IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill ) & 2760 surf_lsm_h%m_res(:,m) = soil_pars_f%pars_xy(7,j,i) 2761 2762 ENDDO 2763 ! 2764 !-- Vertical surfaces 2765 DO l = 0, 3 2766 DO m = 1, surf_lsm_v(l)%ns 2767 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 2768 surf_lsm_v(l)%building_covered(m) ) 2769 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 2770 surf_lsm_v(l)%building_covered(m) ) 2771 2772 IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill ) & 2773 surf_lsm_v(l)%alpha_vg(:,m) = soil_pars_f%pars_xy(0,j,i) 2774 IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill ) & 2775 surf_lsm_v(l)%l_vg(:,m) = soil_pars_f%pars_xy(1,j,i) 2776 IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill ) & 2777 surf_lsm_v(l)%n_vg(:,m) = soil_pars_f%pars_xy(2,j,i) 2778 IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill ) & 2779 surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i) 2780 IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill ) & 2781 surf_lsm_v(l)%m_sat(:,m) = soil_pars_f%pars_xy(4,j,i) 2782 IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill ) & 2783 surf_lsm_v(l)%m_fc(:,m) = soil_pars_f%pars_xy(5,j,i) 2784 IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill ) & 2785 surf_lsm_v(l)%m_wilt(:,m) = soil_pars_f%pars_xy(6,j,i) 2786 IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill ) & 2787 surf_lsm_v(l)%m_res(:,m) = soil_pars_f%pars_xy(7,j,i) 2788 2789 ENDDO 2790 ENDDO 2791 ! 2792 !-- Level of detail = 2, i.e. soil parameters can be set at each soil 2793 !-- layer individually. 2794 ELSE 2795 ! 2796 !-- Horizontal surfaces 2797 DO m = 1, surf_lsm_h%ns 2798 i = surf_lsm_h%i(m) 2799 j = surf_lsm_h%j(m) 2800 2801 DO k = nzb_soil, nzt_soil 2802 IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill ) & 2803 surf_lsm_h%alpha_vg(k,m) = soil_pars_f%pars_xyz(0,k,j,i) 2804 IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill ) & 2805 surf_lsm_h%l_vg(k,m) = soil_pars_f%pars_xyz(1,k,j,i) 2806 IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill ) & 2807 surf_lsm_h%n_vg(k,m) = soil_pars_f%pars_xyz(2,k,j,i) 2808 IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill ) & 2809 surf_lsm_h%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i) 2810 IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill ) & 2811 surf_lsm_h%m_sat(k,m) = soil_pars_f%pars_xyz(4,k,j,i) 2812 IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill ) & 2813 surf_lsm_h%m_fc(k,m) = soil_pars_f%pars_xyz(5,k,j,i) 2814 IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill ) & 2815 surf_lsm_h%m_wilt(k,m) = soil_pars_f%pars_xyz(6,k,j,i) 2816 IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill ) & 2817 surf_lsm_h%m_res(k,m) = soil_pars_f%pars_xyz(7,k,j,i) 2818 ENDDO 2819 2820 ENDDO 2821 ! 2822 !-- Vertical surfaces 2823 DO l = 0, 3 2824 DO m = 1, surf_lsm_v(l)%ns 2825 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 2826 surf_lsm_v(l)%building_covered(m) ) 2827 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 2828 surf_lsm_v(l)%building_covered(m) ) 2829 2830 DO k = nzb_soil, nzt_soil 2831 IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill ) & 2832 surf_lsm_v(l)%alpha_vg(k,m) = soil_pars_f%pars_xyz(0,k,j,i) 2833 IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill ) & 2834 surf_lsm_v(l)%l_vg(k,m) = soil_pars_f%pars_xyz(1,k,j,i) 2835 IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill ) & 2836 surf_lsm_v(l)%n_vg(k,m) = soil_pars_f%pars_xyz(2,k,j,i) 2837 IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill ) & 2838 surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i) 2839 IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill ) & 2840 surf_lsm_v(l)%m_sat(k,m) = soil_pars_f%pars_xyz(4,k,j,i) 2841 IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill ) & 2842 surf_lsm_v(l)%m_fc(k,m) = soil_pars_f%pars_xyz(5,k,j,i) 2843 IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill ) & 2844 surf_lsm_v(l)%m_wilt(k,m) = soil_pars_f%pars_xyz(6,k,j,i) 2845 IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill ) & 2846 surf_lsm_v(l)%m_res(k,m) = soil_pars_f%pars_xyz(7,k,j,i) 2847 ENDDO 2848 2849 ENDDO 2850 ENDDO 2851 2852 ENDIF 2853 ENDIF 2854 2855 ! 2856 !-- Level 1, initialization of vegetation parameters. A horizontally 2857 !-- homogeneous distribution is assumed here. 2858 IF ( vegetation_type /= 0 ) THEN 2859 2860 IF ( min_canopy_resistance == 9999999.9_wp ) THEN 2861 min_canopy_resistance = vegetation_pars(ind_v_rc_min,vegetation_type) 2862 ENDIF 2863 2864 IF ( leaf_area_index == 9999999.9_wp ) THEN 2865 leaf_area_index = vegetation_pars(ind_v_rc_lai,vegetation_type) 2866 ENDIF 2867 2868 IF ( vegetation_coverage == 9999999.9_wp ) THEN 2869 vegetation_coverage = vegetation_pars(ind_v_c_veg,vegetation_type) 2870 ENDIF 2871 2872 IF ( canopy_resistance_coefficient == 9999999.9_wp ) THEN 2873 canopy_resistance_coefficient= vegetation_pars(ind_v_gd,vegetation_type) 2874 ENDIF 2875 2876 IF ( z0_vegetation == 9999999.9_wp ) THEN 2877 z0_vegetation = vegetation_pars(ind_v_z0,vegetation_type) 2878 ENDIF 2879 2880 IF ( z0h_vegetation == 9999999.9_wp ) THEN 2881 z0h_vegetation = vegetation_pars(ind_v_z0qh,vegetation_type) 2882 ENDIF 2883 2884 IF ( lambda_surface_stable == 9999999.9_wp ) THEN 2885 lambda_surface_stable = vegetation_pars(ind_v_lambda_s,vegetation_type) 2886 ENDIF 2887 2888 IF ( lambda_surface_unstable == 9999999.9_wp ) THEN 2889 lambda_surface_unstable = vegetation_pars(ind_v_lambda_u,vegetation_type) 2890 ENDIF 2891 2892 IF ( f_shortwave_incoming == 9999999.9_wp ) THEN 2893 f_shortwave_incoming = vegetation_pars(ind_v_f_sw_in,vegetation_type) 2894 ENDIF 2895 2896 IF ( c_surface == 9999999.9_wp ) THEN 2897 c_surface = vegetation_pars(ind_v_c_surf,vegetation_type) 2898 ENDIF 2899 2900 IF ( albedo_type == 9999999 .AND. albedo == 9999999.9_wp ) THEN 2901 albedo_type = INT(vegetation_pars(ind_v_at,vegetation_type)) 2902 ENDIF 2903 2904 IF ( emissivity == 9999999.9_wp ) THEN 2905 emissivity = vegetation_pars(ind_v_emis,vegetation_type) 2906 ENDIF 2907 2908 ENDIF 2909 ! 2910 !-- Map values onto horizontal elemements 2911 DO m = 1, surf_lsm_h%ns 2912 IF ( surf_lsm_h%vegetation_surface(m) ) THEN 2913 surf_lsm_h%r_canopy_min(m) = min_canopy_resistance 2914 surf_lsm_h%lai(m) = leaf_area_index 2915 surf_lsm_h%c_veg(m) = vegetation_coverage 2916 surf_lsm_h%g_d(m) = canopy_resistance_coefficient 2917 surf_lsm_h%z0(m) = z0_vegetation 2918 surf_lsm_h%z0h(m) = z0h_vegetation 2919 surf_lsm_h%z0q(m) = z0h_vegetation 2920 surf_lsm_h%lambda_surface_s(m) = lambda_surface_stable 2921 surf_lsm_h%lambda_surface_u(m) = lambda_surface_unstable 2922 surf_lsm_h%f_sw_in(m) = f_shortwave_incoming 2923 surf_lsm_h%c_surface(m) = c_surface 2924 surf_lsm_h%albedo_type(ind_veg,m) = albedo_type 2925 surf_lsm_h%emissivity(ind_veg,m) = emissivity 2926 ELSE 2927 surf_lsm_h%lai(m) = 0.0_wp 2928 surf_lsm_h%c_veg(m) = 0.0_wp 2929 surf_lsm_h%g_d(m) = 0.0_wp 2930 ENDIF 2931 2932 ENDDO 2933 ! 2934 !-- Map values onto vertical elements, even though this does not make 2935 !-- much sense. 2936 DO l = 0, 3 2937 DO m = 1, surf_lsm_v(l)%ns 2938 IF ( surf_lsm_v(l)%vegetation_surface(m) ) THEN 2939 surf_lsm_v(l)%r_canopy_min(m) = min_canopy_resistance 2940 surf_lsm_v(l)%lai(m) = leaf_area_index 2941 surf_lsm_v(l)%c_veg(m) = vegetation_coverage 2942 surf_lsm_v(l)%g_d(m) = canopy_resistance_coefficient 2943 surf_lsm_v(l)%z0(m) = z0_vegetation 2944 surf_lsm_v(l)%z0h(m) = z0h_vegetation 2945 surf_lsm_v(l)%z0q(m) = z0h_vegetation 2946 surf_lsm_v(l)%lambda_surface_s(m) = lambda_surface_stable 2947 surf_lsm_v(l)%lambda_surface_u(m) = lambda_surface_unstable 2948 surf_lsm_v(l)%f_sw_in(m) = f_shortwave_incoming 2949 surf_lsm_v(l)%c_surface(m) = c_surface 2950 surf_lsm_v(l)%albedo_type(ind_veg,m) = albedo_type 2951 surf_lsm_v(l)%emissivity(ind_veg,m) = emissivity 2952 ELSE 2953 surf_lsm_v(l)%lai(m) = 0.0_wp 2954 surf_lsm_v(l)%c_veg(m) = 0.0_wp 2955 surf_lsm_v(l)%g_d(m) = 0.0_wp 2956 ENDIF 2957 ENDDO 2958 ENDDO 2959 2960 ! 2961 !-- Level 2, initialization of vegation parameters via vegetation_type read 2962 !-- from file. Vegetation parameters are initialized for each (y,x)-grid point 2963 !-- individually using default paramter settings according to the given 2964 !-- vegetation type. 2965 IF ( vegetation_type_f%from_file ) THEN 2966 ! 2967 !-- Horizontal surfaces 2968 DO m = 1, surf_lsm_h%ns 2969 i = surf_lsm_h%i(m) 2970 j = surf_lsm_h%j(m) 2971 2972 st = vegetation_type_f%var(j,i) 2973 IF ( st /= vegetation_type_f%fill .AND. st /= 0 ) THEN 2974 surf_lsm_h%r_canopy_min(m) = vegetation_pars(ind_v_rc_min,st) 2975 surf_lsm_h%lai(m) = vegetation_pars(ind_v_rc_lai,st) 2976 surf_lsm_h%c_veg(m) = vegetation_pars(ind_v_c_veg,st) 2977 surf_lsm_h%g_d(m) = vegetation_pars(ind_v_gd,st) 2978 surf_lsm_h%z0(m) = vegetation_pars(ind_v_z0,st) 2979 surf_lsm_h%z0h(m) = vegetation_pars(ind_v_z0qh,st) 2980 surf_lsm_h%z0q(m) = vegetation_pars(ind_v_z0qh,st) 2981 surf_lsm_h%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st) 2982 surf_lsm_h%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st) 2983 surf_lsm_h%f_sw_in(m) = vegetation_pars(ind_v_f_sw_in,st) 2984 surf_lsm_h%c_surface(m) = vegetation_pars(ind_v_c_surf,st) 2985 surf_lsm_h%albedo_type(ind_veg,m) = INT( vegetation_pars(ind_v_at,st) ) 2986 surf_lsm_h%emissivity(ind_veg,m) = vegetation_pars(ind_v_emis,st) 2987 ENDIF 2988 ENDDO 2989 ! 2990 !-- Vertical surfaces 2991 DO l = 0, 3 2992 DO m = 1, surf_lsm_v(l)%ns 2993 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 2994 surf_lsm_v(l)%building_covered(m) ) 2995 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 2996 surf_lsm_v(l)%building_covered(m) ) 2997 2998 st = vegetation_type_f%var(j,i) 2999 IF ( st /= vegetation_type_f%fill .AND. st /= 0 ) THEN 3000 surf_lsm_v(l)%r_canopy_min(m) = vegetation_pars(ind_v_rc_min,st) 3001 surf_lsm_v(l)%lai(m) = vegetation_pars(ind_v_rc_lai,st) 3002 surf_lsm_v(l)%c_veg(m) = vegetation_pars(ind_v_c_veg,st) 3003 surf_lsm_v(l)%g_d(m) = vegetation_pars(ind_v_gd,st) 3004 surf_lsm_v(l)%z0(m) = vegetation_pars(ind_v_z0,st) 3005 surf_lsm_v(l)%z0h(m) = vegetation_pars(ind_v_z0qh,st) 3006 surf_lsm_v(l)%z0q(m) = vegetation_pars(ind_v_z0qh,st) 3007 surf_lsm_v(l)%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st) 3008 surf_lsm_v(l)%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st) 3009 surf_lsm_v(l)%f_sw_in(m) = vegetation_pars(ind_v_f_sw_in,st) 3010 surf_lsm_v(l)%c_surface(m) = vegetation_pars(ind_v_c_surf,st) 3011 surf_lsm_v(l)%albedo_type(ind_veg,m) = INT( vegetation_pars(ind_v_at,st) ) 3012 surf_lsm_v(l)%emissivity(ind_veg,m) = vegetation_pars(ind_v_emis,st) 3013 ENDIF 3014 ENDDO 3015 ENDDO 3016 ENDIF 3017 ! 3018 !-- Level 3, initialization of vegation parameters at single (x,y) 3019 !-- position via vegetation_pars read from file. 3020 IF ( vegetation_pars_f%from_file ) THEN 3021 ! 3022 !-- Horizontal surfaces 3023 DO m = 1, surf_lsm_h%ns 3024 3025 i = surf_lsm_h%i(m) 3026 j = surf_lsm_h%j(m) 3027 ! 3028 !-- If surface element is not a vegetation surface and any value in 3029 !-- vegetation_pars is given, neglect this information and give an 3030 !-- informative message that this value will not be used. 3031 IF ( .NOT. surf_lsm_h%vegetation_surface(m) .AND. & 3032 ANY( vegetation_pars_f%pars_xy(:,j,i) /= & 3033 vegetation_pars_f%fill ) ) THEN 3034 WRITE( message_string, * ) & 3035 'surface element at grid point (j,i) = (', & 3036 j, i, ') is not a vegation surface, ', & 3037 'so that information given in ', & 3038 'vegetation_pars at this point is neglected.' 3039 CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 ) 3040 ELSE 3041 3042 IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /= & 3043 vegetation_pars_f%fill ) & 3044 surf_lsm_h%r_canopy_min(m) = & 3045 vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) 3046 IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /= & 3047 vegetation_pars_f%fill ) & 3048 surf_lsm_h%lai(m) = & 3049 vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) 3050 IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /= & 3051 vegetation_pars_f%fill ) & 3052 surf_lsm_h%c_veg(m) = & 3053 vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) 3054 IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /= & 3055 vegetation_pars_f%fill ) & 3056 surf_lsm_h%g_d(m) = & 3057 vegetation_pars_f%pars_xy(ind_v_gd,j,i) 3058 IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /= & 3059 vegetation_pars_f%fill ) & 3060 surf_lsm_h%z0(m) = & 3061 vegetation_pars_f%pars_xy(ind_v_z0,j,i) 3062 IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /= & 3063 vegetation_pars_f%fill ) THEN 3064 surf_lsm_h%z0h(m) = & 3065 vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) 3066 surf_lsm_h%z0q(m) = & 3067 vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) 3068 ENDIF 3069 IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /= & 3070 vegetation_pars_f%fill ) & 3071 surf_lsm_h%lambda_surface_s(m) = & 3072 vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) 3073 IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /= & 3074 vegetation_pars_f%fill ) & 3075 surf_lsm_h%lambda_surface_u(m) = & 3076 vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) 3077 IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /= & 3078 vegetation_pars_f%fill ) & 3079 surf_lsm_h%f_sw_in(m) = & 3080 vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) 3081 IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /= & 3082 vegetation_pars_f%fill ) & 3083 surf_lsm_h%c_surface(m) = & 3084 vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) 3085 IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /= & 3086 vegetation_pars_f%fill ) & 3087 surf_lsm_h%albedo_type(ind_veg,m) = & 3088 INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) ) 3089 IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /= & 3090 vegetation_pars_f%fill ) & 3091 surf_lsm_h%emissivity(ind_veg,m) = & 3092 vegetation_pars_f%pars_xy(ind_v_emis,j,i) 3093 ENDIF 3094 ENDDO 3095 ! 3096 !-- Vertical surfaces 3097 DO l = 0, 3 3098 DO m = 1, surf_lsm_v(l)%ns 3099 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 3100 surf_lsm_v(l)%building_covered(m) ) 3101 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 3102 surf_lsm_v(l)%building_covered(m) ) 3103 ! 3104 !-- If surface element is not a vegetation surface and any value in 3105 !-- vegetation_pars is given, neglect this information and give an 3106 !-- informative message that this value will not be used. 3107 IF ( .NOT. surf_lsm_v(l)%vegetation_surface(m) .AND. & 3108 ANY( vegetation_pars_f%pars_xy(:,j,i) /= & 3109 vegetation_pars_f%fill ) ) THEN 3110 WRITE( message_string, * ) & 3111 'surface element at grid point (j,i) = (', & 3112 j, i, ') is not a vegation surface, ', & 3113 'so that information given in ', & 3114 'vegetation_pars at this point is neglected.' 3115 CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 ) 3116 ELSE 3117 3118 IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /= & 3119 vegetation_pars_f%fill ) & 3120 surf_lsm_v(l)%r_canopy_min(m) = & 3121 vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) 3122 IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /= & 3123 vegetation_pars_f%fill ) & 3124 surf_lsm_v(l)%lai(m) = & 3125 vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) 3126 IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /= & 3127 vegetation_pars_f%fill ) & 3128 surf_lsm_v(l)%c_veg(m) = & 3129 vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) 3130 IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /= & 3131 vegetation_pars_f%fill ) & 3132 surf_lsm_v(l)%g_d(m) = & 3133 vegetation_pars_f%pars_xy(ind_v_gd,j,i) 3134 IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /= & 3135 vegetation_pars_f%fill ) & 3136 surf_lsm_v(l)%z0(m) = & 3137 vegetation_pars_f%pars_xy(ind_v_z0,j,i) 3138 IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /= & 3139 vegetation_pars_f%fill ) THEN 3140 surf_lsm_v(l)%z0h(m) = & 3141 vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) 3142 surf_lsm_v(l)%z0q(m) = & 3143 vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) 3144 ENDIF 3145 IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /= & 3146 vegetation_pars_f%fill ) & 3147 surf_lsm_v(l)%lambda_surface_s(m) = & 3148 vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) 3149 IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /= & 3150 vegetation_pars_f%fill ) & 3151 surf_lsm_v(l)%lambda_surface_u(m) = & 3152 vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) 3153 IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /= & 3154 vegetation_pars_f%fill ) & 3155 surf_lsm_v(l)%f_sw_in(m) = & 3156 vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) 3157 IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /= & 3158 vegetation_pars_f%fill ) & 3159 surf_lsm_v(l)%c_surface(m) = & 3160 vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) 3161 IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /= & 3162 vegetation_pars_f%fill ) & 3163 surf_lsm_v(l)%albedo_type(ind_veg,m) = & 3164 INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) ) 3165 IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /= & 3166 vegetation_pars_f%fill ) & 3167 surf_lsm_v(l)%emissivity(ind_veg,m) = & 3168 vegetation_pars_f%pars_xy(ind_v_emis,j,i) 3169 ENDIF 3170 3171 ENDDO 3172 ENDDO 3173 ENDIF 3174 3175 ! 3176 !-- Level 1, initialization of water parameters. A horizontally 3177 !-- homogeneous distribution is assumed here. 3178 IF ( water_type /= 0 ) THEN 3179 3180 IF ( water_temperature == 9999999.9_wp ) THEN 3181 water_temperature = water_pars(ind_w_temp,water_type) 3182 ENDIF 3183 3184 IF ( z0_water == 9999999.9_wp ) THEN 3185 z0_water = water_pars(ind_w_z0,water_type) 3186 ENDIF 3187 3188 IF ( z0h_water == 9999999.9_wp ) THEN 3189 z0h_water = water_pars(ind_w_z0h,water_type) 3190 ENDIF 3191 3192 IF ( albedo_type == 9999999 .AND. albedo == 9999999.9_wp ) THEN 3193 albedo_type = INT(water_pars(ind_w_at,water_type)) 3194 ENDIF 3195 3196 IF ( emissivity == 9999999.9_wp ) THEN 3197 emissivity = water_pars(ind_w_emis,water_type) 3198 ENDIF 3199 3200 ENDIF 3201 ! 3202 !-- Map values onto horizontal elemements 3203 DO m = 1, surf_lsm_h%ns 3204 IF ( surf_lsm_h%water_surface(m) ) THEN 3205 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) & 3206 t_soil_h%var_2d(:,m) = water_temperature 3207 surf_lsm_h%z0(m) = z0_water 3208 surf_lsm_h%z0h(m) = z0h_water 3209 surf_lsm_h%z0q(m) = z0h_water 3210 surf_lsm_h%lambda_surface_s(m) = 1.0E10_wp 3211 surf_lsm_h%lambda_surface_u(m) = 1.0E10_wp 3212 surf_lsm_h%c_surface(m) = 0.0_wp 3213 surf_lsm_h%albedo_type(ind_wat,m) = albedo_type 3214 surf_lsm_h%emissivity(ind_wat,m) = emissivity 3215 ENDIF 3216 ENDDO 3217 ! 3218 !-- Map values onto vertical elements, even though this does not make 3219 !-- much sense. 3220 DO l = 0, 3 3221 DO m = 1, surf_lsm_v(l)%ns 3222 IF ( surf_lsm_v(l)%water_surface(m) ) THEN 3223 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) & 3224 t_soil_v(l)%var_2d(:,m) = water_temperature 3225 surf_lsm_v(l)%z0(m) = z0_water 3226 surf_lsm_v(l)%z0h(m) = z0h_water 3227 surf_lsm_v(l)%z0q(m) = z0h_water 3228 surf_lsm_v(l)%lambda_surface_s(m) = 1.0E10_wp 3229 surf_lsm_v(l)%lambda_surface_u(m) = 1.0E10_wp 3230 surf_lsm_v(l)%c_surface(m) = 0.0_wp 3231 surf_lsm_v(l)%albedo_type(ind_wat,m) = albedo_type 3232 surf_lsm_v(l)%emissivity(ind_wat,m) = emissivity 3233 ENDIF 3234 ENDDO 3235 ENDDO 3236 ! 3237 ! 3238 !-- Level 2, initialization of water parameters via water_type read 3239 !-- from file. Water surfaces are initialized for each (y,x)-grid point 3240 !-- individually using default paramter settings according to the given 3241 !-- water type. 3242 !-- Note, parameter 3/4 of water_pars are albedo and emissivity, 3243 !-- whereas paramter 3/4 of water_pars_f are heat conductivities! 3244 IF ( water_type_f%from_file ) THEN 3245 ! 3246 !-- Horizontal surfaces 3247 DO m = 1, surf_lsm_h%ns 3248 i = surf_lsm_h%i(m) 3249 j = surf_lsm_h%j(m) 3250 3251 st = water_type_f%var(j,i) 3252 IF ( st /= water_type_f%fill .AND. st /= 0 ) THEN 3253 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) & 3254 t_soil_h%var_2d(:,m) = water_pars(ind_w_temp,st) 3255 surf_lsm_h%z0(m) = water_pars(ind_w_z0,st) 3256 surf_lsm_h%z0h(m) = water_pars(ind_w_z0h,st) 3257 surf_lsm_h%z0q(m) = water_pars(ind_w_z0h,st) 3258 surf_lsm_h%lambda_surface_s(m) = water_pars(ind_w_lambda_s,st) 3259 surf_lsm_h%lambda_surface_u(m) = water_pars(ind_w_lambda_u,st) 3260 surf_lsm_h%c_surface(m) = 0.0_wp 3261 surf_lsm_h%albedo_type(ind_wat,m) = INT( water_pars(ind_w_at,st) ) 3262 surf_lsm_h%emissivity(ind_wat,m) = water_pars(ind_w_emis,st) 3263 ENDIF 3264 ENDDO 3265 ! 3266 !-- Vertical surfaces 3267 DO l = 0, 3 3268 DO m = 1, surf_lsm_v(l)%ns 3269 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 3270 surf_lsm_v(l)%building_covered(m) ) 3271 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 3272 surf_lsm_v(l)%building_covered(m) ) 3273 3274 st = water_type_f%var(j,i) 3275 IF ( st /= water_type_f%fill .AND. st /= 0 ) THEN 3276 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) & 3277 t_soil_v(l)%var_2d(:,m) = water_pars(ind_w_temp,st) 3278 surf_lsm_v(l)%z0(m) = water_pars(ind_w_z0,st) 3279 surf_lsm_v(l)%z0h(m) = water_pars(ind_w_z0h,st) 3280 surf_lsm_v(l)%z0q(m) = water_pars(ind_w_z0h,st) 3281 surf_lsm_v(l)%lambda_surface_s(m) = & 3282 water_pars(ind_w_lambda_s,st) 3283 surf_lsm_v(l)%lambda_surface_u(m) = & 3284 water_pars(ind_w_lambda_u,st) 3285 surf_lsm_v(l)%c_surface(m) = 0.0_wp 3286 surf_lsm_v(l)%albedo_type(ind_wat,m) = & 3287 INT( water_pars(ind_w_at,st) ) 3288 surf_lsm_v(l)%emissivity(ind_wat,m) = & 3289 water_pars(ind_w_emis,st) 3290 ENDIF 3291 ENDDO 3292 ENDDO 3293 ENDIF 3294 3295 ! 3296 !-- Level 3, initialization of water parameters at single (x,y) 3297 !-- position via water_pars read from file. 3298 IF ( water_pars_f%from_file ) THEN 3299 ! 3300 !-- Horizontal surfaces 3301 DO m = 1, surf_lsm_h%ns 3302 i = surf_lsm_h%i(m) 3303 j = surf_lsm_h%j(m) 3304 ! 3305 !-- If surface element is not a water surface and any value in 3306 !-- water_pars is given, neglect this information and give an 3307 !-- informative message that this value will not be used. 3308 IF ( .NOT. surf_lsm_h%water_surface(m) .AND. & 3309 ANY( water_pars_f%pars_xy(:,j,i) /= water_pars_f%fill ) ) THEN 3310 WRITE( message_string, * ) & 3311 'surface element at grid point (j,i) = (', & 3312 j, i, ') is not a water surface, ', & 3313 'so that information given in ', & 3314 'water_pars at this point is neglected.' 3315 CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 ) 3316 ELSE 3317 IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /= & 3318 water_pars_f%fill .AND. & 3319 TRIM( initializing_actions ) /= 'read_restart_data' ) & 3320 t_soil_h%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i) 3321 3322 IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /= water_pars_f%fill ) & 3323 surf_lsm_h%z0(m) = water_pars_f%pars_xy(ind_w_z0,j,i) 3324 3325 IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /= water_pars_f%fill )& 3326 THEN 3327 surf_lsm_h%z0h(m) = water_pars_f%pars_xy(ind_w_z0h,j,i) 3328 surf_lsm_h%z0q(m) = water_pars_f%pars_xy(ind_w_z0h,j,i) 3329 ENDIF 3330 3331 IF ( water_pars_f%pars_xy(ind_w_lambda_s,j,i) /= & 3332 water_pars_f%fill ) & 3333 surf_lsm_h%lambda_surface_s(m) = & 3334 water_pars_f%pars_xy(ind_w_lambda_s,j,i) 3335 3336 IF ( water_pars_f%pars_xy(ind_w_lambda_u,j,i) /= & 3337 water_pars_f%fill ) & 3338 surf_lsm_h%lambda_surface_u(m) = & 3339 water_pars_f%pars_xy(ind_w_lambda_u,j,i) 3340 3341 IF ( water_pars_f%pars_xy(ind_w_at,j,i) /= & 3342 water_pars_f%fill ) & 3343 surf_lsm_h%albedo_type(ind_wat,m) = & 3344 INT( water_pars_f%pars_xy(ind_w_at,j,i) ) 3345 3346 IF ( water_pars_f%pars_xy(ind_w_emis,j,i) /= & 3347 water_pars_f%fill ) & 3348 surf_lsm_h%emissivity(ind_wat,m) = & 3349 water_pars_f%pars_xy(ind_w_emis,j,i) 3350 ENDIF 3351 ENDDO 3352 ! 3353 !-- Vertical surfaces 3354 DO l = 0, 3 3355 DO m = 1, surf_lsm_v(l)%ns 3356 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 3357 surf_lsm_v(l)%building_covered(m) ) 3358 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 3359 surf_lsm_v(l)%building_covered(m) ) 3360 ! 3361 !-- If surface element is not a water surface and any value in 3362 !-- water_pars is given, neglect this information and give an 3363 !-- informative message that this value will not be used. 3364 IF ( .NOT. surf_lsm_v(l)%water_surface(m) .AND. & 3365 ANY( water_pars_f%pars_xy(:,j,i) /= & 3366 water_pars_f%fill ) ) THEN 3367 WRITE( message_string, * ) & 3368 'surface element at grid point (j,i) = (', & 3369 j, i, ') is not a water surface, ', & 3370 'so that information given in ', & 3371 'water_pars at this point is neglected.' 3372 CALL message( 'land_surface_model_mod', 'PA0999', & 3373 0, 0, 0, 6, 0 ) 3374 ELSE 3375 3376 IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /= & 3377 water_pars_f%fill .AND. & 3378 TRIM( initializing_actions ) /= 'read_restart_data' ) & 3379 t_soil_v(l)%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i) 3380 3381 IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /= & 3382 water_pars_f%fill ) & 3383 surf_lsm_v(l)%z0(m) = water_pars_f%pars_xy(ind_w_z0,j,i) 3384 3385 IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /= & 3386 water_pars_f%fill ) THEN 3387 surf_lsm_v(l)%z0h(m) = water_pars_f%pars_xy(ind_w_z0h,j,i) 3388 surf_lsm_v(l)%z0q(m) = water_pars_f%pars_xy(ind_w_z0h,j,i) 3389 ENDIF 3390 3391 IF ( water_pars_f%pars_xy(ind_w_lambda_s,j,i) /= & 3392 water_pars_f%fill ) & 3393 surf_lsm_v(l)%lambda_surface_s(m) = & 3394 water_pars_f%pars_xy(ind_w_lambda_s,j,i) 3395 3396 IF ( water_pars_f%pars_xy(ind_w_lambda_u,j,i) /= & 3397 water_pars_f%fill ) & 3398 surf_lsm_v(l)%lambda_surface_u(m) = & 3399 water_pars_f%pars_xy(ind_w_lambda_u,j,i) 3400 3401 IF ( water_pars_f%pars_xy(ind_w_at,j,i) /= & 3402 water_pars_f%fill ) & 3403 surf_lsm_v(l)%albedo_type(ind_wat,m) = & 3404 INT( water_pars_f%pars_xy(ind_w_at,j,i) ) 3405 3406 IF ( water_pars_f%pars_xy(ind_w_emis,j,i) /= & 3407 water_pars_f%fill ) & 3408 surf_lsm_v(l)%emissivity(ind_wat,m) = & 3409 water_pars_f%pars_xy(ind_w_emis,j,i) 3410 ENDIF 3411 ENDDO 3412 ENDDO 3413 3414 ENDIF 3415 ! 3416 !-- Initialize pavement-type surfaces, level 1 3417 IF ( pavement_type /= 0 ) THEN 3418 3419 ! 3420 !-- When a pavement_type is used, overwrite a possible setting of 3421 !-- the pavement depth as it is already defined by the pavement type 3422 pavement_depth_level = 0 3423 3424 IF ( z0_pavement == 9999999.9_wp ) THEN 3425 z0_pavement = pavement_pars(ind_p_z0,pavement_type) 3426 ENDIF 3427 3428 IF ( z0h_pavement == 9999999.9_wp ) THEN 3429 z0h_pavement = pavement_pars(ind_p_z0h,pavement_type) 3430 ENDIF 3431 3432 IF ( pavement_heat_conduct == 9999999.9_wp ) THEN 3433 pavement_heat_conduct = pavement_subsurface_pars_1(0,pavement_type) 3434 ENDIF 3435 3436 IF ( pavement_heat_capacity == 9999999.9_wp ) THEN 3437 pavement_heat_capacity = pavement_subsurface_pars_2(0,pavement_type) 3438 ENDIF 3439 3440 IF ( albedo_type == 9999999 .AND. albedo == 9999999.9_wp ) THEN 3441 albedo_type = INT(pavement_pars(ind_p_at,pavement_type)) 3442 ENDIF 3443 3444 IF ( emissivity == 9999999.9_wp ) THEN 3445 emissivity = pavement_pars(ind_p_emis,pavement_type) 3446 ENDIF 3447 3448 ! 3449 !-- If the depth level of the pavement is not set, determine it from 3450 !-- lookup table. 3451 IF ( pavement_depth_level == 0 ) THEN 3452 DO k = nzb_soil, nzt_soil 3453 IF ( pavement_subsurface_pars_1(k,pavement_type) == 9999999.9_wp & 3454 .OR. pavement_subsurface_pars_2(k,pavement_type) == 9999999.9_wp)& 3455 THEN 3456 nzt_pavement = k-1 3457 EXIT 3458 ENDIF 3459 ENDDO 3460 ELSE 3461 nzt_pavement = pavement_depth_level 3462 ENDIF 3463 3464 ENDIF 3465 ! 3466 !-- Level 1 initialization of pavement type surfaces. Horizontally 3467 !-- homogeneous characteristics are assumed 3468 surf_lsm_h%nzt_pavement = pavement_depth_level 3469 DO m = 1, surf_lsm_h%ns 3470 IF ( surf_lsm_h%pavement_surface(m) ) THEN 3471 surf_lsm_h%nzt_pavement(m) = nzt_pavement 3472 surf_lsm_h%z0(m) = z0_pavement 3473 surf_lsm_h%z0h(m) = z0h_pavement 3474 surf_lsm_h%z0q(m) = z0h_pavement 3475 surf_lsm_h%lambda_surface_s(m) = pavement_heat_conduct & 3476 * ddz_soil(nzb_soil) & 3477 * 2.0_wp 3478 surf_lsm_h%lambda_surface_u(m) = pavement_heat_conduct & 3479 * ddz_soil(nzb_soil) & 3480 * 2.0_wp 3481 surf_lsm_h%c_surface(m) = pavement_heat_capacity & 3482 * dz_soil(nzb_soil) & 3483 * 0.25_wp 3484 3485 surf_lsm_h%albedo_type(ind_pav,m) = albedo_type 3486 surf_lsm_h%emissivity(ind_pav,m) = emissivity 3487 3488 IF ( pavement_type /= 0 ) THEN 3489 DO k = nzb_soil, surf_lsm_h%nzt_pavement(m) 3490 surf_lsm_h%lambda_h_def(k,m) = & 3491 pavement_subsurface_pars_1(k,pavement_type) 3492 surf_lsm_h%rho_c_total_def(k,m) = & 3493 pavement_subsurface_pars_2(k,pavement_type) 3494 ENDDO 3495 ELSE 3496 surf_lsm_v(l)%lambda_h_def(:,m) = pavement_heat_conduct 3497 surf_lsm_v(l)%rho_c_total_def(:,m) = pavement_heat_capacity 3498 ENDIF 3499 ENDIF 3500 ENDDO 3501 3502 DO l = 0, 3 3503 surf_lsm_v(l)%nzt_pavement = pavement_depth_level 3504 DO m = 1, surf_lsm_v(l)%ns 3505 IF ( surf_lsm_v(l)%pavement_surface(m) ) THEN 3506 surf_lsm_v(l)%nzt_pavement(m) = nzt_pavement 3507 surf_lsm_v(l)%z0(m) = z0_pavement 3508 surf_lsm_v(l)%z0h(m) = z0h_pavement 3509 surf_lsm_v(l)%z0q(m) = z0h_pavement 3510 surf_lsm_v(l)%lambda_surface_s(m) = pavement_heat_conduct & 3511 * ddz_soil(nzb_soil) & 3512 * 2.0_wp 3513 surf_lsm_v(l)%lambda_surface_u(m) = pavement_heat_conduct & 3514 * ddz_soil(nzb_soil) & 3515 * 2.0_wp 3516 surf_lsm_v(l)%c_surface(m) = pavement_heat_capacity & 3517 * dz_soil(nzb_soil) & 3518 * 0.25_wp 3519 3520 surf_lsm_v(l)%albedo_type(ind_pav,m) = albedo_type 3521 surf_lsm_v(l)%emissivity(ind_pav,m) = emissivity 3522 3523 IF ( pavement_type /= 0 ) THEN 3524 DO k = nzb_soil, surf_lsm_v(l)%nzt_pavement(m) 3525 surf_lsm_v(l)%lambda_h_def(k,m) = & 3526 pavement_subsurface_pars_1(k,pavement_type) 3527 surf_lsm_v(l)%rho_c_total_def(k,m) = & 3528 pavement_subsurface_pars_2(k,pavement_type) 3529 ENDDO 3530 ELSE 3531 surf_lsm_v(l)%lambda_h_def(:,m) = pavement_heat_conduct 3532 surf_lsm_v(l)%rho_c_total_def(:,m) = pavement_heat_capacity 3533 ENDIF 3534 ENDIF 3535 ENDDO 3536 ENDDO 3537 ! 3538 !-- Level 2 initialization of pavement type surfaces via pavement_type read 3539 !-- from file. Pavement surfaces are initialized for each (y,x)-grid point 3540 !-- individually. 3541 IF ( pavement_type_f%from_file ) THEN 3542 ! 3543 !-- Horizontal surfaces 3544 DO m = 1, surf_lsm_h%ns 3545 i = surf_lsm_h%i(m) 3546 j = surf_lsm_h%j(m) 3547 3548 st = pavement_type_f%var(j,i) 3549 IF ( st /= pavement_type_f%fill .AND. st /= 0 ) THEN 3550 ! 3551 !-- Determine deepmost index of pavement layer 3552 DO k = nzb_soil, nzt_soil 3553 IF ( pavement_subsurface_pars_1(k,st) == 9999999.9_wp & 3554 .OR. pavement_subsurface_pars_2(k,st) == 9999999.9_wp) & 3555 THEN 3556 surf_lsm_h%nzt_pavement(m) = k-1 3557 EXIT 3558 ENDIF 3559 ENDDO 3560 3561 surf_lsm_h%z0(m) = pavement_pars(ind_p_z0,st) 3562 surf_lsm_h%z0h(m) = pavement_pars(ind_p_z0h,st) 3563 surf_lsm_h%z0q(m) = pavement_pars(ind_p_z0h,st) 3564 3565 surf_lsm_h%lambda_surface_s(m) = & 3566 pavement_subsurface_pars_1(0,st) & 3567 * ddz_soil(nzb_soil) & 3568 * 2.0_wp 3569 surf_lsm_h%lambda_surface_u(m) = & 3570 pavement_subsurface_pars_1(0,st) & 3571 * ddz_soil(nzb_soil) & 3572 * 2.0_wp 3573 surf_lsm_h%c_surface(m) = & 3574 pavement_subsurface_pars_2(0,st)& 3575 * dz_soil(nzb_soil) & 3576 * 0.25_wp 3577 surf_lsm_h%albedo_type(ind_pav,m) = INT( pavement_pars(ind_p_at,st) ) 3578 surf_lsm_h%emissivity(ind_pav,m) = pavement_pars(ind_p_emis,st) 3579 3580 DO k = nzb_soil, surf_lsm_h%nzt_pavement(m) 3581 surf_lsm_h%lambda_h_def(k,m) = & 3582 pavement_subsurface_pars_1(k,pavement_type) 3583 surf_lsm_h%rho_c_total_def(k,m) = & 3584 pavement_subsurface_pars_2(k,pavement_type) 3585 ENDDO 3586 ENDIF 3587 ENDDO 3588 ! 3589 !-- Vertical surfaces 3590 DO l = 0, 3 3591 DO m = 1, surf_lsm_v(l)%ns 3592 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 3593 surf_lsm_v(l)%building_covered(m) ) 3594 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 3595 surf_lsm_v(l)%building_covered(m) ) 3596 3597 st = pavement_type_f%var(j,i) 3598 IF ( st /= pavement_type_f%fill .AND. st /= 0 ) THEN 3599 ! 3600 !-- Determine deepmost index of pavement layer 3601 DO k = nzb_soil, nzt_soil 3602 IF ( pavement_subsurface_pars_1(k,st) == 9999999.9_wp & 3603 .OR. pavement_subsurface_pars_2(k,st) == 9999999.9_wp) & 3604 THEN 3605 surf_lsm_v(l)%nzt_pavement(m) = k-1 3606 EXIT 3607 ENDIF 3608 ENDDO 3609 3610 surf_lsm_v(l)%z0(m) = pavement_pars(ind_p_z0,st) 3611 surf_lsm_v(l)%z0h(m) = pavement_pars(ind_p_z0h,st) 3612 surf_lsm_v(l)%z0q(m) = pavement_pars(ind_p_z0h,st) 3613 3614 surf_lsm_v(l)%lambda_surface_s(m) = & 3615 pavement_subsurface_pars_1(0,st) & 3616 * ddz_soil(nzb_soil) & 3617 * 2.0_wp 3618 surf_lsm_v(l)%lambda_surface_u(m) = & 3619 pavement_subsurface_pars_1(0,st) & 3620 * ddz_soil(nzb_soil) & 3621 * 2.0_wp 3622 3623 surf_lsm_v(l)%c_surface(m) = & 3624 pavement_subsurface_pars_2(0,st) & 3625 * dz_soil(nzb_soil) & 3626 * 0.25_wp 3627 surf_lsm_v(l)%albedo_type(ind_pav,m) = & 3628 INT( pavement_pars(ind_p_at,st) ) 3629 surf_lsm_v(l)%emissivity(ind_pav,m) = & 3630 pavement_pars(ind_p_emis,st) 3631 3632 DO k = nzb_soil, surf_lsm_h%nzt_pavement(m) 3633 surf_lsm_v(l)%lambda_h_def(k,m) = & 3634 pavement_subsurface_pars_1(k,pavement_type) 3635 surf_lsm_v(l)%rho_c_total_def(k,m) = & 3636 pavement_subsurface_pars_2(k,pavement_type) 3637 ENDDO 3638 ENDIF 3639 ENDDO 3640 ENDDO 3641 ENDIF 3642 ! 3643 !-- Level 3, initialization of pavement parameters at single (x,y) 3644 !-- position via pavement_pars read from file. 3645 IF ( pavement_pars_f%from_file ) THEN 3646 ! 3647 !-- Horizontal surfaces 3648 DO m = 1, surf_lsm_h%ns 3649 i = surf_lsm_h%i(m) 3650 j = surf_lsm_h%j(m) 3651 ! 3652 !-- If surface element is not a pavement surface and any value in 3653 !-- pavement_pars is given, neglect this information and give an 3654 !-- informative message that this value will not be used. 3655 IF ( .NOT. surf_lsm_h%pavement_surface(m) .AND. & 3656 ANY( pavement_pars_f%pars_xy(:,j,i) /= & 3657 pavement_pars_f%fill ) ) THEN 3658 WRITE( message_string, * ) & 3659 'surface element at grid point (j,i) = (', & 3660 j, i, ') is not a pavement surface, ', & 3661 'so that information given in ', & 3662 'pavement_pars at this point is neglected.' 3663 CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 ) 3664 ELSE 3665 IF ( pavement_pars_f%pars_xy(ind_p_z0,j,i) /= & 3666 pavement_pars_f%fill ) & 3667 surf_lsm_h%z0(m) = pavement_pars_f%pars_xy(ind_p_z0,j,i) 3668 IF ( pavement_pars_f%pars_xy(ind_p_z0h,j,i) /= & 3669 pavement_pars_f%fill ) THEN 3670 surf_lsm_h%z0h(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i) 3671 surf_lsm_h%z0q(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i) 3672 ENDIF 3673 IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i) & 3674 /= pavement_subsurface_pars_f%fill ) THEN 3675 surf_lsm_h%lambda_surface_s(m) = & 3676 pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)& 3677 * ddz_soil(nzb_soil) & 3678 * 2.0_wp 3679 surf_lsm_h%lambda_surface_u(m) = & 3680 pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)& 3681 * ddz_soil(nzb_soil) & 3682 * 2.0_wp 3683 ENDIF 3684 IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i) /= & 3685 pavement_subsurface_pars_f%fill ) THEN 3686 surf_lsm_h%c_surface(m) = & 3687 pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i) & 3688 * dz_soil(nzb_soil) & 3689 * 0.25_wp 3690 ENDIF 3691 IF ( pavement_pars_f%pars_xy(ind_p_at,j,i) /= & 3692 pavement_pars_f%fill ) & 3693 surf_lsm_h%albedo_type(ind_pav,m) = & 3694 INT( pavement_pars(ind_p_at,st) ) 3695 IF ( pavement_pars_f%pars_xy(ind_p_emis,j,i) /= & 3696 pavement_pars_f%fill ) & 3697 surf_lsm_h%emissivity(ind_pav,m) = & 3698 pavement_pars(ind_p_emis,st) 3699 ENDIF 3700 3701 ENDDO 3702 ! 3703 !-- Vertical surfaces 3704 DO l = 0, 3 3705 DO m = 1, surf_lsm_v(l)%ns 3706 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 3707 surf_lsm_v(l)%building_covered(m) ) 3708 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 3709 surf_lsm_v(l)%building_covered(m) ) 3710 ! 3711 !-- If surface element is not a pavement surface and any value in 3712 !-- pavement_pars is given, neglect this information and give an 3713 !-- informative message that this value will not be used. 3714 IF ( .NOT. surf_lsm_v(l)%pavement_surface(m) .AND. & 3715 ANY( pavement_pars_f%pars_xy(:,j,i) /= & 3716 pavement_pars_f%fill ) ) THEN 3717 WRITE( message_string, * ) & 3718 'surface element at grid point (j,i) = (', & 3719 j, i, ') is not a pavement surface, ', & 3720 'so that information given in ', & 3721 'pavement_pars at this point is neglected.' 3722 CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 ) 3723 ELSE 3724 3725 IF ( pavement_pars_f%pars_xy(ind_p_z0,j,i) /= & 3726 pavement_pars_f%fill ) & 3727 surf_lsm_v(l)%z0(m) = pavement_pars_f%pars_xy(ind_p_z0,j,i) 3728 IF ( pavement_pars_f%pars_xy(ind_p_z0h,j,i) /= & 3729 pavement_pars_f%fill ) THEN 3730 surf_lsm_v(l)%z0h(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i) 3731 surf_lsm_v(l)%z0q(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i) 3732 ENDIF 3733 IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)& 3734 /= pavement_subsurface_pars_f%fill ) THEN 3735 surf_lsm_v(l)%lambda_surface_s(m) = & 3736 pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)& 3737 * ddz_soil(nzb_soil) & 3738 * 2.0_wp 3739 surf_lsm_v(l)%lambda_surface_u(m) = & 3740 pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)& 3741 * ddz_soil(nzb_soil) & 3742 * 2.0_wp 3743 ENDIF 3744 IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i) & 3745 /= pavement_subsurface_pars_f%fill ) THEN 3746 surf_lsm_v(l)%c_surface(m) = & 3747 pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i)& 3748 * dz_soil(nzb_soil) & 3749 * 0.25_wp 3750 ENDIF 3751 IF ( pavement_pars_f%pars_xy(ind_p_at,j,i) /= & 3752 pavement_pars_f%fill ) & 3753 surf_lsm_v(l)%albedo_type(ind_pav,m) = & 3754 INT( pavement_pars(ind_p_at,st) ) 3755 3756 IF ( pavement_pars_f%pars_xy(ind_p_emis,j,i) /= & 3757 pavement_pars_f%fill ) & 3758 surf_lsm_v(l)%emissivity(ind_pav,m) = & 3759 pavement_pars(ind_p_emis,st) 3760 ENDIF 3761 ENDDO 3762 ENDDO 3763 ENDIF 3764 ! 3765 !-- Moreover, for grid points which are flagged with pavement-type 0 or whre 3766 !-- pavement_subsurface_pars_f is provided, soil heat conductivity and 3767 !-- capacity are initialized with parameters given in 3768 !-- pavement_subsurface_pars read from file. 3769 IF ( pavement_subsurface_pars_f%from_file ) THEN 3770 ! 3771 !-- Set pavement depth to nzt_soil. Please note, this is just a 3772 !-- workaround at the moment. 3773 DO m = 1, surf_lsm_h%ns 3774 IF ( surf_lsm_h%pavement_surface(m) ) THEN 3775 3776 i = surf_lsm_h%i(m) 3777 j = surf_lsm_h%j(m) 3778 3779 surf_lsm_h%nzt_pavement(m) = nzt_soil 3780 3781 DO k = nzb_soil, nzt_soil 3782 surf_lsm_h%lambda_h_def(k,m) = & 3783 pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,k,j,i) 3784 surf_lsm_h%rho_c_total_def(k,m) = & 3785 pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,k,j,i) 3786 ENDDO 3787 3788 ENDIF 3789 ENDDO 3790 DO l = 0, 3 3791 DO m = 1, surf_lsm_v(l)%ns 3792 IF ( surf_lsm_v(l)%pavement_surface(m) ) THEN 3793 3794 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 3795 surf_lsm_v(l)%building_covered(m) ) 3796 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 3797 surf_lsm_v(l)%building_covered(m) ) 3798 3799 surf_lsm_v(l)%nzt_pavement(m) = nzt_soil 3800 3801 DO k = nzb_soil, nzt_soil 3802 surf_lsm_v(l)%lambda_h_def(k,m) = & 3803 pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,k,j,i) 3804 surf_lsm_v(l)%rho_c_total_def(k,m) = & 3805 pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,k,j,i) 3806 ENDDO 3807 3808 ENDIF 3809 ENDDO 3810 ENDDO 3811 ENDIF 2497 3812 ! 2498 3813 !-- Initial run actions 2499 3814 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 2500 3815 ! 2501 !-- First, for horizontal surfaces 2502 ! 2503 !-- Set artifical values for ts and us so that r_a has its initial value 2504 !-- for the first time step. Only for interior core domain, not for ghost points 3816 !-- First, initialize soil temperature and moisture. 3817 !-- According to the initialization for surface and soil parameters, 3818 !-- initialize soil moisture and temperature via a level approach. This 3819 !-- is to assure that all surface elements are initialized, even if 3820 !-- data provided from input file contains fill values at some locations. 3821 !-- Level 1, initialization via profiles given in parameter file 2505 3822 DO m = 1, surf_lsm_h%ns 2506 2507 t_soil_h%var_2d(:,m) = 0.0_wp 2508 m_soil_h%var_2d(:,m) = 0.0_wp 2509 m_liq_h%var_1d(m) = 0.0_wp 2510 2511 !-- Map user settings of T and q for each soil layer 2512 !-- (make sure that the soil moisture does not drop below the permanent 2513 !-- wilting point) -> problems with devision by zero) 2514 IF ( surf_lsm_h%water_surface(m) ) THEN 2515 2516 surf_lsm_h%z0(m) = z0_water 2517 surf_lsm_h%z0h(m) = z0h_water 2518 surf_lsm_h%z0q(m) = z0h_water 2519 t_soil_h%var_2d(:,m) = water_temperature 2520 surf_lsm_h%lambda_surface_s(m) = 1.0E10_wp 2521 surf_lsm_h%lambda_surface_u(m) = 1.0E10_wp 2522 surf_lsm_h%c_surface(m) = 0.0_wp 2523 2524 ELSEIF ( surf_lsm_h%pavement_surface(m) ) THEN 2525 2526 surf_lsm_h%z0(m) = z0_pavement 2527 surf_lsm_h%z0h(m) = z0h_pavement 2528 surf_lsm_h%z0q(m) = z0h_pavement 3823 IF ( surf_lsm_h%vegetation_surface(m) .OR. & 3824 surf_lsm_h%pavement_surface(m) ) THEN 2529 3825 DO k = nzb_soil, nzt_soil 2530 t_soil_h%var_2d(k,m) 2531 m_soil_h%var_2d(k,m) 3826 t_soil_h%var_2d(k,m) = soil_temperature(k) 3827 m_soil_h%var_2d(k,m) = soil_moisture(k) 2532 3828 ENDDO 2533 3829 t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1) 2534 surf_lsm_h%lambda_surface_s(m) = pavement_heat_conduct & 2535 * ddz_soil(nzb_soil) & 2536 * 2.0_wp 2537 surf_lsm_h%lambda_surface_u(m) = surf_lsm_h%lambda_surface_s(m) 2538 surf_lsm_h%c_surface(m) = pavement_heat_capacity & 2539 * dz_soil(nzb_soil) & 2540 * 0.25_wp 2541 2542 IF ( pavement_type /= 0 ) THEN 3830 ENDIF 3831 ENDDO 3832 DO l = 0, 3 3833 DO m = 1, surf_lsm_v(l)%ns 3834 IF ( surf_lsm_v(l)%vegetation_surface(m) .OR. & 3835 surf_lsm_v(l)%pavement_surface(m) ) THEN 2543 3836 DO k = nzb_soil, nzt_soil 2544 surf_lsm_h%lambda_h_def(k,m) = pavement_subsurface_pars_1(k,pavement_type)2545 surf_lsm_h%rho_c_total_def(k,m) = pavement_subsurface_pars_2(k,pavement_type)3837 t_soil_v(l)%var_2d(k,m) = soil_temperature(k) 3838 m_soil_v(l)%var_2d(k,m) = soil_moisture(k) 2546 3839 ENDDO 2547 ELSE 2548 surf_lsm_h%lambda_h_def(:,m) = pavement_heat_conduct 2549 surf_lsm_h%rho_c_total_def(:,m) = pavement_heat_capacity 3840 t_soil_v(l)%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1) 2550 3841 ENDIF 2551 2552 ELSEIF ( surf_lsm_h%vegetation_surface(m) ) THEN 2553 2554 surf_lsm_h%z0(m) = z0_vegetation 2555 surf_lsm_h%z0h(m) = z0h_vegetation 2556 surf_lsm_h%z0q(m) = z0h_vegetation 2557 DO k = nzb_soil, nzt_soil 2558 t_soil_h%var_2d(k,m) = soil_temperature(k) 2559 m_soil_h%var_2d(k,m) = soil_moisture(k) 2560 ENDDO 2561 t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1) 2562 surf_lsm_h%lambda_surface_s(m) = lambda_surface_stable 2563 surf_lsm_h%lambda_surface_u(m) = lambda_surface_unstable 2564 surf_lsm_h%c_surface(m) = c_surface 2565 3842 ENDDO 3843 ENDDO 3844 ! 3845 !-- Level 2, if soil moisture and/or temperature are 3846 !-- provided from file, interpolate / extrapolate the provided data 3847 !-- onto the respective soil layers. Please note, both, zs as well as 3848 !-- init_3d%z_soil indicate a depth with positive values, so that no 3849 !-- distinction between atmosphere is required concerning interpolation. 3850 !-- Start with soil moisture 3851 IF ( init_3d%from_file_msoil ) THEN 3852 3853 IF ( init_3d%lod_msoil == 1 ) THEN 3854 DO m = 1, surf_lsm_h%ns 3855 3856 CALL netcdf_data_input_interpolate( & 3857 m_soil_h%var_2d(nzb_soil:nzt_soil,m), & 3858 init_3d%msoil_init(:), & 3859 zs(nzb_soil:nzt_soil), init_3d%z_soil, & 3860 nzb_soil, nzt_soil, & 3861 nzb_soil, init_3d%nzs-1 ) 3862 ENDDO 3863 DO l = 0, 3 3864 DO m = 1, surf_lsm_v(l)%ns 3865 3866 CALL netcdf_data_input_interpolate( & 3867 m_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),& 3868 init_3d%msoil_init(:), & 3869 zs(nzb_soil:nzt_soil), init_3d%z_soil, & 3870 nzb_soil, nzt_soil, & 3871 nzb_soil, init_3d%nzs-1 ) 3872 ENDDO 3873 ENDDO 3874 ELSE 3875 3876 DO m = 1, surf_lsm_h%ns 3877 i = surf_lsm_h%i(m) 3878 j = surf_lsm_h%j(m) 3879 3880 IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil ) & 3881 CALL netcdf_data_input_interpolate( & 3882 m_soil_h%var_2d(nzb_soil:nzt_soil,m), & 3883 init_3d%msoil(:,j,i), & 3884 zs(nzb_soil:nzt_soil), init_3d%z_soil, & 3885 nzb_soil, nzt_soil, & 3886 nzb_soil, init_3d%nzs-1 ) 3887 ENDDO 3888 DO l = 0, 3 3889 DO m = 1, surf_lsm_v(l)%ns 3890 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 3891 surf_lsm_v(l)%building_covered(m) ) 3892 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 3893 surf_lsm_v(l)%building_covered(m) ) 3894 3895 IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil ) & 3896 CALL netcdf_data_input_interpolate( & 3897 m_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),& 3898 init_3d%msoil(:,j,i), & 3899 zs(nzb_soil:nzt_soil), init_3d%z_soil, & 3900 nzb_soil, nzt_soil, & 3901 nzb_soil, init_3d%nzs-1 ) 3902 ENDDO 3903 ENDDO 2566 3904 ENDIF 2567 3905 3906 ENDIF 3907 ! 3908 !-- Soil temperature 3909 IF ( init_3d%from_file_tsoil ) THEN 3910 3911 IF ( init_3d%lod_tsoil == 1 ) THEN ! change to 1 if provided correctly by INIFOR 3912 DO m = 1, surf_lsm_h%ns 3913 3914 CALL netcdf_data_input_interpolate( & 3915 t_soil_h%var_2d(nzb_soil:nzt_soil,m), & 3916 init_3d%tsoil_init(:), & 3917 zs(nzb_soil:nzt_soil), init_3d%z_soil, & 3918 nzb_soil, nzt_soil, & 3919 nzb_soil, init_3d%nzs-1 ) 3920 t_soil_h%var_2d(nzt_soil+1,m) = t_soil_h%var_2d(nzt_soil,m) 3921 ENDDO 3922 DO l = 0, 3 3923 DO m = 1, surf_lsm_v(l)%ns 3924 3925 CALL netcdf_data_input_interpolate( & 3926 t_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),& 3927 init_3d%tsoil_init(:), & 3928 zs(nzb_soil:nzt_soil), init_3d%z_soil, & 3929 nzb_soil, nzt_soil, & 3930 nzb_soil, init_3d%nzs-1 ) 3931 t_soil_v(l)%var_2d(nzt_soil+1,m) = & 3932 t_soil_v(l)%var_2d(nzt_soil,m) 3933 ENDDO 3934 ENDDO 3935 ELSE 3936 3937 DO m = 1, surf_lsm_h%ns 3938 i = surf_lsm_h%i(m) 3939 j = surf_lsm_h%j(m) 3940 3941 IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil ) & 3942 CALL netcdf_data_input_interpolate( & 3943 t_soil_h%var_2d(nzb_soil:nzt_soil,m), & 3944 init_3d%tsoil(:,j,i), & 3945 zs(nzb_soil:nzt_soil), init_3d%z_soil, & 3946 nzb_soil, nzt_soil, & 3947 nzb_soil, init_3d%nzs-1 ) 3948 t_soil_h%var_2d(nzt_soil+1,m) = t_soil_h%var_2d(nzt_soil,m) 3949 ENDDO 3950 DO l = 0, 3 3951 DO m = 1, surf_lsm_v(l)%ns 3952 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 3953 surf_lsm_v(l)%building_covered(m) ) 3954 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 3955 surf_lsm_v(l)%building_covered(m) ) 3956 3957 IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil ) & 3958 CALL netcdf_data_input_interpolate( & 3959 t_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),& 3960 init_3d%tsoil(:,j,i), & 3961 zs(nzb_soil:nzt_soil), init_3d%z_soil, & 3962 nzb_soil, nzt_soil, & 3963 nzb_soil, init_3d%nzs-1 ) 3964 t_soil_v(l)%var_2d(nzt_soil+1,m) = & 3965 t_soil_v(l)%var_2d(nzt_soil,m) 3966 ENDDO 3967 ENDDO 3968 ENDIF 3969 ENDIF 3970 ! 3971 !-- Further initialization 3972 DO m = 1, surf_lsm_h%ns 3973 2568 3974 i = surf_lsm_h%i(m) 2569 3975 j = surf_lsm_h%j(m) … … 2575 3981 IF ( surf_lsm_h%lambda_surface_s(m) == 0.0_wp ) THEN 2576 3982 t_surface_h%var_1d(:) = t_soil_h%var_2d(nzb_soil,m) 2577 surf_lsm_h%pt_surface( :) = t_soil_h%var_2d(nzb_soil,m) / exn3983 surf_lsm_h%pt_surface(m) = t_soil_h%var_2d(nzb_soil,m) / exn 2578 3984 ELSE 2579 t_surface_h%var_1d( :) = pt_surface* exn2580 surf_lsm_h%pt_surface( :) = pt_surface3985 t_surface_h%var_1d(m) = pt(k-1,j,i) * exn 3986 surf_lsm_h%pt_surface(m) = pt(k-1,j,i) 2581 3987 ENDIF 2582 3988 2583 IF ( cloud_physics .OR. cloud_droplets )THEN2584 pt1= pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)3989 IF ( cloud_physics .OR. cloud_droplets ) THEN 3990 surf_lsm_h%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i) 2585 3991 ELSE 2586 pt1 = pt(k,j,i) 2587 ENDIF 3992 surf_lsm_h%pt1(m) = pt(k,j,i) 3993 ENDIF 3994 2588 3995 2589 3996 ! 2590 3997 !-- Assure that r_a cannot be zero at model start 2591 IF ( pt1 == surf_lsm_h%pt_surface(m) ) pt1 = pt1 + 1.0E-10_wp 3998 IF ( surf_lsm_h%pt1(m) == surf_lsm_h%pt_surface(m) ) & 3999 surf_lsm_h%pt1(m) = surf_lsm_h%pt1(m) + 1.0E-20_wp 2592 4000 2593 4001 surf_lsm_h%us(m) = 0.1_wp 2594 surf_lsm_h%ts(m) = ( pt1 - surf_lsm_h%pt_surface(m) )&4002 surf_lsm_h%ts(m) = ( surf_lsm_h%pt1(m) - surf_lsm_h%pt_surface(m) )& 2595 4003 / surf_lsm_h%r_a(m) 2596 4004 surf_lsm_h%shf(m) = - surf_lsm_h%us(m) * surf_lsm_h%ts(m) & 2597 4005 * rho_surface 2598 2599 4006 ENDDO 2600 2601 4007 ! 2602 4008 !-- Vertical surfaces 2603 4009 DO l = 0, 3 2604 4010 DO m = 1, surf_lsm_v(l)%ns 2605 2606 t_soil_v(l)%var_2d = 0.0_wp 2607 m_soil_v(l)%var_2d = 0.0_wp 2608 m_liq_v(l)%var_1d = 0.0_wp 2609 2610 2611 !-- Map user settings of T and q for each soil layer 2612 !-- (make sure that the soil moisture does not drop below the permanent 2613 !-- wilting point) -> problems with devision by zero) 2614 IF ( surf_lsm_v(l)%water_surface(m) ) THEN 2615 surf_lsm_v(l)%z0(m) = z0_water 2616 surf_lsm_v(l)%z0h(m) = z0h_water 2617 surf_lsm_v(l)%z0q(m) = z0h_water 2618 t_soil_v(l)%var_2d(:,m) = water_temperature 2619 surf_lsm_v(l)%lambda_surface_s(m) = 1.0E10_wp 2620 surf_lsm_v(l)%lambda_surface_u(m) = 1.0E10_wp 2621 surf_lsm_v(l)%c_surface(m) = 0.0_wp 2622 2623 ELSEIF ( surf_lsm_v(l)%pavement_surface(m) ) THEN 2624 surf_lsm_v(l)%z0(m) = z0_pavement 2625 surf_lsm_v(l)%z0h(m) = z0h_pavement 2626 surf_lsm_v(l)%z0q(m) = z0h_pavement 2627 DO k = nzb_soil, nzt_soil 2628 t_soil_v(l)%var_2d(k,m) = soil_temperature(k) 2629 m_soil_v(l)%var_2d(k,m) = soil_moisture(k) 2630 ENDDO 2631 t_soil_v(l)%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1) 2632 surf_lsm_v(l)%lambda_surface_s(m) = pavement_heat_conduct & 2633 * ddz_soil(nzb_soil) & 2634 * 2.0_wp 2635 surf_lsm_v(l)%lambda_surface_u(m) = surf_lsm_v(l)%lambda_surface_s(m) 2636 surf_lsm_v(l)%c_surface(m) = pavement_heat_capacity & 2637 * dz_soil(nzb_soil) & 2638 * 0.25_wp 2639 IF ( pavement_type /= 0 ) THEN 2640 DO k = nzb_soil, nzt_soil 2641 surf_lsm_v(l)%lambda_h_def(k,m) = pavement_subsurface_pars_1(k,pavement_type) 2642 surf_lsm_v(l)%rho_c_total_def(k,m) = pavement_subsurface_pars_2(k,pavement_type) 2643 ENDDO 2644 ELSE 2645 surf_lsm_v(l)%lambda_h_def(:,m) = pavement_heat_conduct 2646 surf_lsm_v(l)%rho_c_total_def(:,m) = pavement_heat_capacity 2647 ENDIF 2648 2649 2650 ELSEIF ( surf_lsm_v(l)%vegetation_surface(m) ) THEN 2651 2652 surf_lsm_v(l)%z0(m) = z0_vegetation 2653 surf_lsm_v(l)%z0h(m) = z0h_vegetation 2654 surf_lsm_v(l)%z0q(m) = z0h_vegetation 2655 DO k = nzb_soil, nzt_soil 2656 t_soil_v(l)%var_2d(k,m) = soil_temperature(k) 2657 m_soil_v(l)%var_2d(k,m) = soil_moisture(k) 2658 ENDDO 2659 t_soil_v(l)%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1) 2660 surf_lsm_v(l)%lambda_surface_s(m) = lambda_surface_stable 2661 surf_lsm_v(l)%lambda_surface_u(m) = lambda_surface_unstable 2662 surf_lsm_v(l)%c_surface(m) = c_surface 2663 ENDIF 2664 4011 i = surf_lsm_v(l)%i(m) 4012 j = surf_lsm_v(l)%j(m) 4013 k = surf_lsm_v(l)%k(m) 2665 4014 ! 2666 4015 !-- Calculate surface temperature. In case of bare soil, the surface … … 2668 4017 !-- layer 2669 4018 IF ( surf_lsm_v(l)%lambda_surface_s(m) == 0.0_wp ) THEN 2670 t_surface_v(l)%var_1d( :) = t_soil_v(l)%var_2d(nzb_soil,m)2671 surf_lsm_v(l)%pt_surface( :) = t_soil_v(l)%var_2d(nzb_soil,m) / exn4019 t_surface_v(l)%var_1d(m) = t_soil_v(l)%var_2d(nzb_soil,m) 4020 surf_lsm_v(l)%pt_surface(m) = t_soil_v(l)%var_2d(nzb_soil,m) / exn 2672 4021 ELSE 2673 t_surface_v(l)%var_1d(:) = pt_surface * exn 2674 surf_lsm_v(l)%pt_surface(:) = pt_surface 4022 j_off = surf_lsm_v(l)%joff 4023 i_off = surf_lsm_v(l)%ioff 4024 4025 t_surface_v(l)%var_1d(m) = pt(k,j+j_off,i+i_off) * exn 4026 surf_lsm_v(l)%pt_surface(m) = pt(k,j+j_off,i+i_off) 2675 4027 ENDIF 2676 4028 4029 4030 IF ( cloud_physics .OR. cloud_droplets ) THEN 4031 surf_lsm_v(l)%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i) 4032 ELSE 4033 surf_lsm_v(l)%pt1(m) = pt(k,j,i) 4034 ENDIF 4035 4036 ! 4037 !-- Assure that r_a cannot be zero at model start 4038 IF ( surf_lsm_v(l)%pt1(m) == surf_lsm_v(l)%pt_surface(m) ) & 4039 surf_lsm_v(l)%pt1(m) = surf_lsm_v(l)%pt1(m) + 1.0E-20_wp 2677 4040 ! 2678 4041 !-- Set artifical values for ts and us so that r_a has its initial value 2679 4042 !-- for the first time step. Only for interior core domain, not for ghost points 2680 2681 i = surf_lsm_v(l)%i(m)2682 j = surf_lsm_v(l)%j(m)2683 k = surf_lsm_v(l)%k(m)2684 2685 IF ( cloud_physics .OR. cloud_droplets ) THEN2686 pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)2687 ELSE2688 pt1 = pt(k,j,i)2689 ENDIF2690 2691 !2692 !-- Assure that r_a cannot be zero at model start2693 IF ( pt1 == surf_lsm_v(l)%pt_surface(m) ) pt1 = pt1 + 1.0E-10_wp2694 2695 4043 surf_lsm_v(l)%us(m) = 0.1_wp 2696 surf_lsm_v(l)%ts(m) = ( pt1 - surf_lsm_v(l)%pt_surface(m) ) / surf_lsm_v(l)%r_a(m) 2697 surf_lsm_v(l)%shf(m) = - surf_lsm_v(l)%us(m) * surf_lsm_v(l)%ts(m) * rho_surface 4044 surf_lsm_v(l)%ts(m) = ( surf_lsm_v(l)%pt1(m) - surf_lsm_v(l)%pt_surface(m) ) /& 4045 surf_lsm_v(l)%r_a(m) 4046 surf_lsm_v(l)%shf(m) = - surf_lsm_v(l)%us(m) * & 4047 surf_lsm_v(l)%ts(m) * rho_surface 4048 2698 4049 ENDDO 2699 4050 ENDDO … … 2715 4066 ! 2716 4067 !-- Set index offset of surface element, seen from atmospheric grid point 2717 IF ( l == 0 ) THEN 2718 j_off = -1 2719 i_off = 0 2720 ELSEIF ( l == 1 ) THEN 2721 j_off = 1 2722 i_off = 0 2723 ELSEIF ( l == 2 ) THEN 2724 j_off = 0 2725 i_off = -1 2726 ELSEIF ( l == 3 ) THEN 2727 j_off = 0 2728 i_off = 1 2729 ENDIF 4068 j_off = surf_lsm_v(l)%joff 4069 i_off = surf_lsm_v(l)%ioff 4070 2730 4071 DO m = 1, surf_lsm_v(l)%ns 2731 4072 i = surf_lsm_v(l)%i(m) … … 2738 4079 ENDIF 2739 4080 ! 2740 !-- Initialize root fraction2741 !-- Horizontal surfaces4081 !-- Level 1 initialization of root distribution - provided by the user via 4082 !-- via namelist. 2742 4083 DO m = 1, surf_lsm_h%ns 2743 i = surf_lsm_h%i(m)2744 j = surf_lsm_h%j(m)2745 2746 4084 DO k = nzb_soil, nzt_soil 2747 4085 surf_lsm_h%root_fr(k,m) = root_fraction(k) 2748 4086 ENDDO 2749 4087 ENDDO 2750 ! 2751 !-- Vertical surfaces 4088 2752 4089 DO l = 0, 3 2753 4090 DO m = 1, surf_lsm_v(l)%ns 2754 i = surf_lsm_v(l)%i(m)2755 j = surf_lsm_v(l)%j(m)2756 2757 4091 DO k = nzb_soil, nzt_soil 2758 4092 surf_lsm_v(l)%root_fr(k,m) = root_fraction(k) … … 2762 4096 2763 4097 ! 4098 !-- Level 2 initialization of root distribution. 2764 4099 !-- When no root distribution is given by the user, use look-up table to prescribe 2765 !-- the root fraction in the individual soil layers 4100 !-- the root fraction in the individual soil layers. 2766 4101 IF ( ALL( root_fraction == 9999999.9_wp ) ) THEN 2767 2768 4102 ! 2769 4103 !-- First, calculate the index bounds for integration … … 2820 4154 !-- Map calculated root fractions 2821 4155 DO m = 1, surf_lsm_h%ns 2822 i = surf_lsm_h%i(m)2823 j = surf_lsm_h%j(m)2824 4156 DO k = nzb_soil, nzt_soil 2825 4157 IF ( surf_lsm_h%pavement_surface(m) .AND. & 2826 k <= nzt_pavement) THEN4158 k <= surf_lsm_h%nzt_pavement(m) ) THEN 2827 4159 surf_lsm_h%root_fr(k,m) = 0.0_wp 2828 4160 ELSE 2829 4161 surf_lsm_h%root_fr(k,m) = root_fraction(k) 2830 4162 ENDIF 2831 ENDDO 2832 4163 4164 ENDDO 2833 4165 ! 2834 4166 !-- Normalize so that the sum = 1. Only relevant when the root … … 2840 4172 ENDDO 2841 4173 ENDIF 2842 2843 4174 ENDDO 2844 4175 DO l = 0, 3 2845 4176 DO m = 1, surf_lsm_v(l)%ns 2846 i = surf_lsm_v(l)%i(m)2847 j = surf_lsm_v(l)%j(m)2848 2849 4177 DO k = nzb_soil, nzt_soil 2850 4178 IF ( surf_lsm_v(l)%pavement_surface(m) .AND. & 2851 k <= nzt_pavement) THEN4179 k <= surf_lsm_h%nzt_pavement(m) ) THEN 2852 4180 surf_lsm_v(l)%root_fr(k,m) = 0.0_wp 2853 4181 ELSE … … 2864 4192 ENDDO 2865 4193 ENDIF 2866 2867 4194 ENDDO 2868 4195 ENDDO 2869 4196 ENDIF 2870 4197 ! 2871 !-- Map vegetation parameters to the respective 2D arrays 2872 surf_lsm_h%r_canopy_min = min_canopy_resistance 2873 surf_lsm_h%lai = leaf_area_index 2874 surf_lsm_h%c_veg = vegetation_coverage 2875 surf_lsm_h%g_d = canopy_resistance_coefficient 2876 surf_lsm_h%f_sw_in = f_shortwave_incoming 2877 2878 !-- Vertical surfaces 2879 DO l = 0, 3 2880 2881 ! 2882 !-- Map vegetation parameters to the respective 2D arrays 2883 surf_lsm_v(l)%r_canopy_min = min_canopy_resistance 2884 surf_lsm_v(l)%lai = leaf_area_index 2885 surf_lsm_v(l)%c_veg = vegetation_coverage 2886 surf_lsm_v(l)%g_d = canopy_resistance_coefficient 2887 surf_lsm_v(l)%f_sw_in = f_shortwave_incoming 2888 ENDDO 4198 !-- Level 3 initialization of root distribution. 4199 !-- Take value from file 4200 IF ( root_area_density_lsm_f%from_file ) THEN 4201 DO m = 1, surf_lsm_h%ns 4202 IF ( surf_lsm_h%vegetation_surface(m) ) THEN 4203 i = surf_lsm_h%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 4204 surf_lsm_v(l)%building_covered(m) ) 4205 j = surf_lsm_h%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 4206 surf_lsm_v(l)%building_covered(m) ) 4207 DO k = nzb_soil, nzt_soil 4208 surf_lsm_h%root_fr(k,m) = root_area_density_lsm_f%var(k,j,i) 4209 ENDDO 4210 4211 ENDIF 4212 ENDDO 4213 4214 DO l = 0, 3 4215 DO m = 1, surf_lsm_v(l)%ns 4216 IF ( surf_lsm_v(l)%vegetation_surface(m) ) THEN 4217 i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, & 4218 surf_lsm_v(l)%building_covered(m) ) 4219 j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, & 4220 surf_lsm_v(l)%building_covered(m) ) 4221 4222 DO k = nzb_soil, nzt_soil 4223 surf_lsm_v(l)%root_fr(k,m) = root_area_density_lsm_f%var(k,j,i) 4224 ENDDO 4225 4226 ENDIF 4227 ENDDO 4228 ENDDO 4229 4230 ENDIF 2889 4231 2890 4232 ! … … 2912 4254 !-- Store initial profiles of t_soil and m_soil (assuming they are 2913 4255 !-- horizontally homogeneous on this PE) 2914 hom(nzb_soil:nzt_soil,1,90,:) = SPREAD( t_soil_h%var_2d(nzb_soil:nzt_soil,1), & 2915 2, statistic_regions+1 ) 2916 hom(nzb_soil:nzt_soil,1,92,:) = SPREAD( m_soil_h%var_2d(nzb_soil:nzt_soil,1), & 2917 2, statistic_regions+1 ) 4256 !-- DEACTIVATED FOR NOW - leads to error when number of locations with 4257 !-- soil model is zero on a PE. 4258 ! hom(nzb_soil:nzt_soil,1,90,:) = SPREAD( t_soil_h%var_2d(nzb_soil:nzt_soil,1), & 4259 ! 2, statistic_regions+1 ) 4260 ! hom(nzb_soil:nzt_soil,1,92,:) = SPREAD( m_soil_h%var_2d(nzb_soil:nzt_soil,1), & 4261 ! 2, statistic_regions+1 ) 4262 4263 ! 4264 !-- Finally, make some consistency checks. 4265 !-- Ceck for eck for illegal combination of LAI and vegetation coverage. 4266 IF ( ANY( .NOT. surf_lsm_h%pavement_surface .AND. & 4267 surf_lsm_h%lai == 0.0_wp .AND. surf_lsm_h%c_veg == 1.0_wp ) & 4268 ) THEN 4269 message_string = 'For non-pavement surfaces the combination ' // & 4270 ' lai = 0.0 and c_veg = 1.0 is not allowed.' 4271 CALL message( 'lsm_read_restart_data', 'PA0999', 2, 2, 0, 6, 0 ) 4272 ENDIF 4273 4274 DO l = 0, 3 4275 IF ( ANY( .NOT. surf_lsm_v(l)%pavement_surface .AND. & 4276 surf_lsm_v(l)%lai == 0.0_wp .AND. & 4277 surf_lsm_v(l)%c_veg == 1.0_wp ) ) THEN 4278 message_string = 'For non-pavement surfaces the combination ' // & 4279 ' lai = 0.0 and c_veg = 1.0 is not allowed.' 4280 CALL message( 'lsm_read_restart_data', 'PA0999', 2, 2, 0, 6, 0 ) 4281 ENDIF 4282 ENDDO 4283 2918 4284 2919 4285 END SUBROUTINE lsm_init … … 2943 4309 ! 2944 4310 !-- Horizontal surfaces 2945 ALLOCATE ( m_liq_h%var_1d(1:surf_lsm_h%ns) )2946 4311 ALLOCATE ( m_liq_h_p%var_1d(1:surf_lsm_h%ns) ) 2947 4312 ALLOCATE ( t_surface_h%var_1d(1:surf_lsm_h%ns) ) 2948 4313 ALLOCATE ( t_surface_h_p%var_1d(1:surf_lsm_h%ns) ) 2949 ALLOCATE ( m_soil_h%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )2950 4314 ALLOCATE ( m_soil_h_p%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) ) 2951 ALLOCATE ( t_soil_h%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )2952 4315 ALLOCATE ( t_soil_h_p%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) ) 2953 4316 … … 2963 4326 ALLOCATE ( t_soil_v(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) ) 2964 4327 ALLOCATE ( t_soil_v_p(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) ) 4328 ENDDO 4329 ! 4330 !-- Allocate soil temperature and moisture. As these variables might be 4331 !-- already allocated in case of restarts, check this. 4332 IF ( .NOT. ALLOCATED( m_liq_h%var_1d ) ) & 4333 ALLOCATE ( m_liq_h%var_1d(1:surf_lsm_h%ns) ) 4334 IF ( .NOT. ALLOCATED( m_soil_h%var_2d ) ) & 4335 ALLOCATE ( m_soil_h%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) ) 4336 IF ( .NOT. ALLOCATED( t_soil_h%var_2d ) ) & 4337 ALLOCATE ( t_soil_h%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) ) 4338 4339 DO l = 0, 3 4340 IF ( .NOT. ALLOCATED( m_liq_v(l)%var_1d ) ) & 4341 ALLOCATE ( m_liq_v(l)%var_1d(1:surf_lsm_v(l)%ns) ) 4342 IF ( .NOT. ALLOCATED( m_soil_v(l)%var_2d ) ) & 4343 ALLOCATE ( m_soil_v(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) ) 4344 IF ( .NOT. ALLOCATED( t_soil_v(l)%var_2d ) ) & 4345 ALLOCATE ( t_soil_v(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) ) 2965 4346 ENDDO 2966 4347 #else … … 2988 4369 ENDDO 2989 4370 #endif 4371 ! 4372 !-- Allocate array for heat flux in W/m2, required for radiation? 4373 !-- Consider to remove this array 4374 ALLOCATE( surf_lsm_h%surfhf(1:surf_lsm_h%ns) ) 4375 DO l = 0, 3 4376 ALLOCATE( surf_lsm_v(l)%surfhf(1:surf_lsm_v(l)%ns) ) 4377 ENDDO 4378 2990 4379 2991 4380 ! … … 3004 4393 ALLOCATE ( tt_soil_v_m(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) ) 3005 4394 ENDDO 3006 !3007 !-- Allocate skin-surface temperature3008 ALLOCATE ( surf_lsm_h%pt_surface(1:surf_lsm_h%ns) )3009 DO l = 0, 33010 ALLOCATE ( surf_lsm_v(l)%pt_surface(1:surf_lsm_v(l)%ns) )3011 ENDDO3012 4395 3013 4396 ! 3014 4397 !-- Allocate 2D vegetation model arrays 3015 4398 !-- Horizontal surfaces 3016 ALLOCATE ( surf_lsm_h%building_surface(1:surf_lsm_h%ns) ) 3017 ALLOCATE ( surf_lsm_h%c_liq(1:surf_lsm_h%ns) ) 3018 ALLOCATE ( surf_lsm_h%c_surface(1:surf_lsm_h%ns) ) 3019 ALLOCATE ( surf_lsm_h%c_veg(1:surf_lsm_h%ns) ) 3020 ALLOCATE ( surf_lsm_h%f_sw_in(1:surf_lsm_h%ns) ) 3021 ALLOCATE ( surf_lsm_h%ghf(1:surf_lsm_h%ns) ) 3022 ALLOCATE ( surf_lsm_h%g_d(1:surf_lsm_h%ns) ) 3023 ALLOCATE ( surf_lsm_h%lai(1:surf_lsm_h%ns) ) 3024 ALLOCATE ( surf_lsm_h%lambda_surface_u(1:surf_lsm_h%ns) ) 3025 ALLOCATE ( surf_lsm_h%lambda_surface_s(1:surf_lsm_h%ns) ) 4399 ALLOCATE ( surf_lsm_h%building_surface(1:surf_lsm_h%ns) ) 4400 ALLOCATE ( surf_lsm_h%c_liq(1:surf_lsm_h%ns) ) 4401 ALLOCATE ( surf_lsm_h%c_surface(1:surf_lsm_h%ns) ) 4402 ALLOCATE ( surf_lsm_h%c_veg(1:surf_lsm_h%ns) ) 4403 ALLOCATE ( surf_lsm_h%f_sw_in(1:surf_lsm_h%ns) ) 4404 ALLOCATE ( surf_lsm_h%ghf(1:surf_lsm_h%ns) ) 4405 ALLOCATE ( surf_lsm_h%g_d(1:surf_lsm_h%ns) ) 4406 ALLOCATE ( surf_lsm_h%lai(1:surf_lsm_h%ns) ) 4407 ALLOCATE ( surf_lsm_h%lambda_surface_u(1:surf_lsm_h%ns) ) 4408 ALLOCATE ( surf_lsm_h%lambda_surface_s(1:surf_lsm_h%ns) ) 4409 ALLOCATE ( surf_lsm_h%nzt_pavement(1:surf_lsm_h%ns) ) 4410 ALLOCATE ( surf_lsm_h%pavement_surface(1:surf_lsm_h%ns) ) 4411 ALLOCATE ( surf_lsm_h%qsws_soil(1:surf_lsm_h%ns) ) 4412 ALLOCATE ( surf_lsm_h%qsws_liq(1:surf_lsm_h%ns) ) 4413 ALLOCATE ( surf_lsm_h%qsws_veg(1:surf_lsm_h%ns) ) 4414 ALLOCATE ( surf_lsm_h%rad_net_l(1:surf_lsm_h%ns) ) 4415 ALLOCATE ( surf_lsm_h%r_a(1:surf_lsm_h%ns) ) 4416 ALLOCATE ( surf_lsm_h%r_canopy(1:surf_lsm_h%ns) ) 4417 ALLOCATE ( surf_lsm_h%r_soil(1:surf_lsm_h%ns) ) 4418 ALLOCATE ( surf_lsm_h%r_soil_min(1:surf_lsm_h%ns) ) 4419 ALLOCATE ( surf_lsm_h%r_s(1:surf_lsm_h%ns) ) 4420 ALLOCATE ( surf_lsm_h%r_canopy_min(1:surf_lsm_h%ns) ) 3026 4421 ALLOCATE ( surf_lsm_h%vegetation_surface(1:surf_lsm_h%ns) ) 3027 ALLOCATE ( surf_lsm_h%pavement_surface(1:surf_lsm_h%ns) ) 3028 ALLOCATE ( surf_lsm_h%qsws_soil(1:surf_lsm_h%ns) ) 3029 ALLOCATE ( surf_lsm_h%qsws_liq(1:surf_lsm_h%ns) ) 3030 ALLOCATE ( surf_lsm_h%qsws_veg(1:surf_lsm_h%ns) ) 3031 ALLOCATE ( surf_lsm_h%rad_net_l(1:surf_lsm_h%ns) ) 3032 ALLOCATE ( surf_lsm_h%r_a(1:surf_lsm_h%ns) ) 3033 ALLOCATE ( surf_lsm_h%r_canopy(1:surf_lsm_h%ns) ) 3034 ALLOCATE ( surf_lsm_h%r_soil(1:surf_lsm_h%ns) ) 3035 ALLOCATE ( surf_lsm_h%r_soil_min(1:surf_lsm_h%ns) ) 3036 ALLOCATE ( surf_lsm_h%r_s(1:surf_lsm_h%ns) ) 3037 ALLOCATE ( surf_lsm_h%r_canopy_min(1:surf_lsm_h%ns) ) 3038 ALLOCATE ( surf_lsm_h%water_surface(1:surf_lsm_h%ns) ) 3039 3040 surf_lsm_h%water_surface = .FALSE. 3041 surf_lsm_h%pavement_surface = .FALSE. 4422 ALLOCATE ( surf_lsm_h%water_surface(1:surf_lsm_h%ns) ) 4423 4424 surf_lsm_h%water_surface = .FALSE. 4425 surf_lsm_h%pavement_surface = .FALSE. 3042 4426 surf_lsm_h%vegetation_surface = .FALSE. 3043 4427 ! 3044 4428 !-- Vertical surfaces 3045 4429 DO l = 0, 3 3046 ALLOCATE ( surf_lsm_v(l)%building_surface(1:surf_lsm_v(l)%ns) ) 3047 ALLOCATE ( surf_lsm_v(l)%c_liq(1:surf_lsm_v(l)%ns) ) 3048 ALLOCATE ( surf_lsm_v(l)%c_surface(1:surf_lsm_v(l)%ns) ) 3049 ALLOCATE ( surf_lsm_v(l)%c_veg(1:surf_lsm_v(l)%ns) ) 3050 ALLOCATE ( surf_lsm_v(l)%f_sw_in(1:surf_lsm_v(l)%ns) ) 3051 ALLOCATE ( surf_lsm_v(l)%ghf(1:surf_lsm_v(l)%ns) ) 3052 ALLOCATE ( surf_lsm_v(l)%g_d(1:surf_lsm_v(l)%ns) ) 3053 ALLOCATE ( surf_lsm_v(l)%lai(1:surf_lsm_v(l)%ns) ) 3054 ALLOCATE ( surf_lsm_v(l)%lambda_surface_u(1:surf_lsm_v(l)%ns) ) 3055 ALLOCATE ( surf_lsm_v(l)%lambda_surface_s(1:surf_lsm_v(l)%ns) ) 4430 ALLOCATE ( surf_lsm_v(l)%building_surface(1:surf_lsm_v(l)%ns) ) 4431 ALLOCATE ( surf_lsm_v(l)%c_liq(1:surf_lsm_v(l)%ns) ) 4432 ALLOCATE ( surf_lsm_v(l)%c_surface(1:surf_lsm_v(l)%ns) ) 4433 ALLOCATE ( surf_lsm_v(l)%c_veg(1:surf_lsm_v(l)%ns) ) 4434 ALLOCATE ( surf_lsm_v(l)%f_sw_in(1:surf_lsm_v(l)%ns) ) 4435 ALLOCATE ( surf_lsm_v(l)%ghf(1:surf_lsm_v(l)%ns) ) 4436 ALLOCATE ( surf_lsm_v(l)%g_d(1:surf_lsm_v(l)%ns) ) 4437 ALLOCATE ( surf_lsm_v(l)%lai(1:surf_lsm_v(l)%ns) ) 4438 ALLOCATE ( surf_lsm_v(l)%lambda_surface_u(1:surf_lsm_v(l)%ns) ) 4439 ALLOCATE ( surf_lsm_v(l)%lambda_surface_s(1:surf_lsm_v(l)%ns) ) 4440 ALLOCATE ( surf_lsm_v(l)%nzt_pavement(1:surf_lsm_v(l)%ns) ) 4441 ALLOCATE ( surf_lsm_v(l)%pavement_surface(1:surf_lsm_v(l)%ns) ) 4442 ALLOCATE ( surf_lsm_v(l)%qsws_soil(1:surf_lsm_v(l)%ns) ) 4443 ALLOCATE ( surf_lsm_v(l)%qsws_liq(1:surf_lsm_v(l)%ns) ) 4444 ALLOCATE ( surf_lsm_v(l)%qsws_veg(1:surf_lsm_v(l)%ns) ) 4445 ALLOCATE ( surf_lsm_v(l)%rad_net_l(1:surf_lsm_v(l)%ns) ) 4446 ALLOCATE ( surf_lsm_v(l)%r_a(1:surf_lsm_v(l)%ns) ) 4447 ALLOCATE ( surf_lsm_v(l)%r_canopy(1:surf_lsm_v(l)%ns) ) 4448 ALLOCATE ( surf_lsm_v(l)%r_soil(1:surf_lsm_v(l)%ns) ) 4449 ALLOCATE ( surf_lsm_v(l)%r_soil_min(1:surf_lsm_v(l)%ns) ) 4450 ALLOCATE ( surf_lsm_v(l)%r_s(1:surf_lsm_v(l)%ns) ) 4451 ALLOCATE ( surf_lsm_v(l)%r_canopy_min(1:surf_lsm_v(l)%ns) ) 3056 4452 ALLOCATE ( surf_lsm_v(l)%vegetation_surface(1:surf_lsm_v(l)%ns) ) 3057 ALLOCATE ( surf_lsm_v(l)%pavement_surface(1:surf_lsm_v(l)%ns) ) 3058 ALLOCATE ( surf_lsm_v(l)%qsws_soil(1:surf_lsm_v(l)%ns) ) 3059 ALLOCATE ( surf_lsm_v(l)%qsws_liq(1:surf_lsm_v(l)%ns) ) 3060 ALLOCATE ( surf_lsm_v(l)%qsws_veg(1:surf_lsm_v(l)%ns) ) 3061 ALLOCATE ( surf_lsm_v(l)%rad_net_l(1:surf_lsm_v(l)%ns) ) 3062 ALLOCATE ( surf_lsm_v(l)%r_a(1:surf_lsm_v(l)%ns) ) 3063 ALLOCATE ( surf_lsm_v(l)%r_canopy(1:surf_lsm_v(l)%ns) ) 3064 ALLOCATE ( surf_lsm_v(l)%r_soil(1:surf_lsm_v(l)%ns) ) 3065 ALLOCATE ( surf_lsm_v(l)%r_soil_min(1:surf_lsm_v(l)%ns) ) 3066 ALLOCATE ( surf_lsm_v(l)%r_s(1:surf_lsm_v(l)%ns) ) 3067 ALLOCATE ( surf_lsm_v(l)%r_canopy_min(1:surf_lsm_v(l)%ns) ) 3068 ALLOCATE ( surf_lsm_v(l)%water_surface(1:surf_lsm_v(l)%ns) ) 4453 ALLOCATE ( surf_lsm_v(l)%water_surface(1:surf_lsm_v(l)%ns) ) 3069 4454 3070 4455 surf_lsm_v(l)%water_surface = .FALSE. … … 3073 4458 3074 4459 ENDDO 4460 3075 4461 3076 4462 #if ! defined( __nopointer ) … … 3088 4474 m_soil_v => m_soil_v_1; m_soil_v_p => m_soil_v_2 3089 4475 m_liq_v => m_liq_v_1; m_liq_v_p => m_liq_v_2 4476 3090 4477 #endif 3091 4478 … … 3122 4509 residual_moisture, root_fraction, & 3123 4510 saturation_moisture, skip_time_do_lsm, & 3124 soil_moisture, soil_temperature, soil_type, & 4511 soil_moisture, soil_temperature, & 4512 soil_type, & 3125 4513 surface_type, & 3126 4514 vegetation_coverage, vegetation_type, & … … 3204 4592 IF ( horizontal ) THEN 3205 4593 surf => surf_lsm_h 4594 3206 4595 surf_m_soil => m_soil_h 3207 4596 surf_m_soil_p => m_soil_h_p … … 3212 4601 ELSE 3213 4602 surf => surf_lsm_v(l) 4603 3214 4604 surf_m_soil => m_soil_v(l) 3215 4605 surf_m_soil_p => m_soil_v_p(l) … … 3225 4615 DO k = nzb_soil, nzt_soil 3226 4616 3227 IF ( surf%pavement_surface(m) .AND. k <= nzt_pavement ) THEN 3228 4617 IF ( surf%pavement_surface(m) .AND. & 4618 k <= surf%nzt_pavement(m) ) THEN 4619 3229 4620 surf%rho_c_total(k,m) = surf%rho_c_total_def(k,m) 3230 4621 lambda_temp(k) = surf%lambda_h_def(k,m) … … 3266 4657 tend(:) = 0.0_wp 3267 4658 3268 tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) * & 3269 ( surf%lambda_h(nzb_soil,m) & 3270 * ( surf_t_soil%var_2d(nzb_soil+1,m) & 3271 - surf_t_soil%var_2d(nzb_soil,m) ) & 3272 * ddz_soil_center(nzb_soil) + surf%ghf(m) ) & 3273 * ddz_soil(nzb_soil) 4659 tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) * & 4660 ( surf%lambda_h(nzb_soil,m) * ( surf_t_soil%var_2d(nzb_soil+1,m) & 4661 - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil_center(nzb_soil) & 4662 + surf%ghf(m) ) * ddz_soil(nzb_soil) 3274 4663 3275 4664 DO k = nzb_soil+1, nzt_soil 3276 tend(k) = ( 1.0_wp / surf%rho_c_total(k,m) ) & 3277 * ( surf%lambda_h(k,m) & 3278 * ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) ) & 3279 * ddz_soil_center(k) - surf%lambda_h(k-1,m) & 3280 * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) ) & 3281 * ddz_soil_center(k-1) ) * ddz_soil(k) 4665 tend(k) = ( 1.0_wp / surf%rho_c_total(k,m) ) & 4666 * ( surf%lambda_h(k,m) & 4667 * ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) ) & 4668 * ddz_soil_center(k) & 4669 - surf%lambda_h(k-1,m) & 4670 * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) ) & 4671 * ddz_soil_center(k-1) & 4672 ) * ddz_soil(k) 3282 4673 3283 4674 ENDDO 3284 4675 3285 4676 surf_t_soil_p%var_2d(nzb_soil:nzt_soil,m) = & 3286 surf_t_soil%var_2d(nzb_soil:nzt_soil,m) & 3287 + dt_3d * ( tsc(2) & 3288 * tend(nzb_soil:nzt_soil) + tsc(3) & 3289 * surf_tt_soil_m%var_2d(nzb_soil:nzt_soil,m) ) 4677 surf_t_soil%var_2d(nzb_soil:nzt_soil,m) & 4678 + dt_3d * ( tsc(2) & 4679 * tend(nzb_soil:nzt_soil) & 4680 + tsc(3) & 4681 * surf_tt_soil_m%var_2d(nzb_soil:nzt_soil,m) ) 3290 4682 3291 4683 ! … … 3299 4691 intermediate_timestep_count_max ) THEN 3300 4692 DO k = nzb_soil, nzt_soil 3301 surf_tt_soil_m%var_2d(k,m) = -9.5625_wp * tend(k) 3302 + 5.3125_wp&3303 *surf_tt_soil_m%var_2d(k,m)4693 surf_tt_soil_m%var_2d(k,m) = -9.5625_wp * tend(k) + & 4694 5.3125_wp * & 4695 surf_tt_soil_m%var_2d(k,m) 3304 4696 ENDDO 3305 4697 ENDIF … … 3312 4704 !-- In order to prevent water tranport through paved surfaces, 3313 4705 !-- conductivity and diffusivity are set to zero 3314 IF ( surf%pavement_surface(m) .AND. k <= nzt_pavement ) THEN3315 4706 IF ( surf%pavement_surface(m) .AND. & 4707 k <= surf%nzt_pavement(m) ) THEN 3316 4708 lambda_temp(k) = 0.0_wp 3317 4709 gamma_temp(k) = 0.0_wp 3318 4710 3319 4711 ELSE 3320 4712 3321 4713 ! 3322 4714 !-- Calculate soil diffusivity at the center of the soil layers … … 3332 4724 h_vg = ( ( ( surf%m_res(k,m) - surf%m_sat(k,m) ) / & 3333 4725 ( surf%m_res(k,m) - & 3334 MAX( surf_m_soil%var_2d(k,m), surf%m_wilt(k,m))&4726 MAX( surf_m_soil%var_2d(k,m), surf%m_wilt(k,m) )& 3335 4727 ) & 3336 )**( surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp ) & 4728 )**( & 4729 surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp ) & 3337 4730 ) - 1.0_wp & 3338 4731 )**( 1.0_wp / surf%n_vg(k,m) ) / surf%alpha_vg(k,m) … … 3340 4733 gamma_temp(k) = surf%gamma_w_sat(k,m) * ( ( ( 1.0_wp + & 3341 4734 ( surf%alpha_vg(k,m) * h_vg )**surf%n_vg(k,m) & 3342 )**( 1.0_wp - 1.0_wp / surf%n_vg(k,m) ) & 3343 - ( surf%alpha_vg(k,m) * h_vg )**( surf%n_vg(k,m) & 3344 - 1.0_wp) )**2 ) / ( ( 1.0_wp + ( surf%alpha_vg(k,m) & 3345 * h_vg )**surf%n_vg(k,m) )**(( 1.0_wp - 1.0_wp & 3346 / surf%n_vg(k,m) ) * ( surf%l_vg(k,m) + 2.0_wp) )) 4735 )**( & 4736 1.0_wp - 1.0_wp / surf%n_vg(k,m)) - ( & 4737 surf%alpha_vg(k,m) * h_vg )**( surf%n_vg(k,m) & 4738 - 1.0_wp) )**2 ) & 4739 / ( ( 1.0_wp + ( surf%alpha_vg(k,m) * h_vg & 4740 )**surf%n_vg(k,m) )**( ( 1.0_wp - 1.0_wp & 4741 / surf%n_vg(k,m) ) * & 4742 ( surf%l_vg(k,m) + 2.0_wp) ) ) 4743 3347 4744 ENDIF 3348 4745 … … 3381 4778 !-- would cause accumulation of (non-existing) water in the lowest 3382 4779 !-- soil layer 3383 IF ( conserve_water_content .AND. surf_m_soil%var_2d(nzt_soil,m) /= 0.0_wp ) THEN 4780 IF ( conserve_water_content .AND. & 4781 surf_m_soil%var_2d(nzt_soil,m) /= 0.0_wp ) THEN 4782 3384 4783 surf%gamma_w(nzt_soil,m) = 0.0_wp 3385 4784 ELSE … … 3404 4803 IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(k,m) ) THEN 3405 4804 m_total = m_total + surf%root_fr(k,m) & 3406 * surf_m_soil%var_2d(k,m)4805 * surf_m_soil%var_2d(k,m) 3407 4806 ENDIF 3408 4807 ENDDO … … 3411 4810 IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(k,m) ) THEN 3412 4811 root_extr(k) = surf%root_fr(k,m) & 3413 4812 * surf_m_soil%var_2d(k,m) / m_total 3414 4813 ELSE 3415 4814 root_extr(k) = 0.0_wp … … 3428 4827 + surf%qsws_soil(m) ) * drho_l_lv ) & 3429 4828 * ddz_soil(nzb_soil) 4829 3430 4830 3431 4831 DO k = nzb_soil+1, nzt_soil-1 … … 3451 4851 surf_m_soil_p%var_2d(nzb_soil:nzt_soil,m) = & 3452 4852 surf_m_soil%var_2d(nzb_soil:nzt_soil,m) & 3453 + dt_3d * ( tsc(2) * tend(:) & 3454 + tsc(3) * surf_tm_soil_m%var_2d(:,m) ) 4853 + dt_3d * ( tsc(2) * tend(:) & 4854 + tsc(3) * surf_tm_soil_m%var_2d(:,m) ) 4855 3455 4856 ! 3456 4857 !-- Account for dry soils (find a better solution here!) … … 3473 4874 * surf_tm_soil_m%var_2d(k,m) 3474 4875 ENDDO 4876 3475 4877 ENDIF 3476 4878 ENDIF … … 3526 4928 m_liq_h => m_liq_h_1; m_liq_h_p => m_liq_h_2 3527 4929 ENDIF 4930 3528 4931 ! 3529 4932 !-- Vertical surfaces … … 3533 4936 m_soil_v => m_soil_v_1; m_soil_v_p => m_soil_v_2 3534 4937 m_liq_v => m_liq_v_1; m_liq_v_p => m_liq_v_2 4938 3535 4939 ENDIF 3536 4940 … … 3545 4949 m_soil_h => m_soil_h_2; m_soil_h_p => m_soil_h_1 3546 4950 m_liq_h => m_liq_h_2; m_liq_h_p => m_liq_h_1 4951 3547 4952 ENDIF 3548 4953 ! … … 3891 5296 ENDDO 3892 5297 ENDDO 5298 3893 5299 ! 3894 5300 !-- … … 4327 5733 !> Write restart data for land surface model 4328 5734 !------------------------------------------------------------------------------! 4329 SUBROUTINE lsm_ last_actions5735 SUBROUTINE lsm_write_restart_data 4330 5736 4331 5737 … … 4335 5741 4336 5742 IMPLICIT NONE 5743 5744 CHARACTER (LEN=1) :: dum !< dummy to create correct string for creating variable string 5745 INTEGER(iwp) :: l !< index variable for surface orientation 4337 5746 4338 5747 IF ( write_binary ) THEN … … 4352 5761 WRITE ( 14 ) 'lai_av '; WRITE ( 14 ) lai_av 4353 5762 ENDIF 4354 WRITE ( 14 ) 'm_liq '; WRITE ( 14 ) m_liq_h%var_1d4355 5763 IF ( ALLOCATED( m_liq_av ) ) THEN 4356 5764 WRITE ( 14 ) 'm_liq_av '; WRITE ( 14 ) m_liq_av 4357 5765 ENDIF 4358 WRITE ( 14 ) 'm_soil '; WRITE ( 14 ) m_soil_h%var_2d4359 5766 IF ( ALLOCATED( m_soil_av ) ) THEN 4360 5767 WRITE ( 14 ) 'm_soil_av '; WRITE ( 14 ) m_soil_av … … 4369 5776 WRITE ( 14 ) 'qsws_veg_av '; WRITE ( 14 ) qsws_veg_av 4370 5777 ENDIF 4371 WRITE ( 14 ) 't_soil '; WRITE ( 14 ) t_soil_h%var_2d4372 5778 IF ( ALLOCATED( t_soil_av ) ) THEN 4373 5779 WRITE ( 14 ) 't_soil_av '; WRITE ( 14 ) t_soil_av 4374 5780 ENDIF 4375 5781 5782 5783 WRITE ( 14 ) 'ns_h_on_file_lsm ' 5784 WRITE ( 14 ) surf_lsm_h%ns 5785 WRITE ( 14 ) 'ns_v_on_file_lsm ' 5786 WRITE ( 14 ) surf_lsm_v(0:3)%ns 5787 5788 WRITE ( 14 ) 'lsm_start_index_h ' 5789 WRITE ( 14 ) surf_lsm_h%start_index 5790 WRITE ( 14 ) 'lsm_end_index_h ' 5791 WRITE ( 14 ) surf_lsm_h%end_index 5792 WRITE ( 14 ) 't_soil_h ' 5793 WRITE ( 14 ) t_soil_h%var_2d 5794 5795 DO l = 0, 3 5796 WRITE ( 14 ) 'lsm_start_index_v ' 5797 WRITE ( 14 ) surf_lsm_v(l)%start_index 5798 WRITE ( 14 ) 'lsm_end_index_v ' 5799 WRITE ( 14 ) surf_lsm_v(l)%end_index 5800 WRITE( dum, '(I1)') l 5801 WRITE ( 14 ) 't_soil_v(' // dum // ') ' 5802 WRITE ( 14 ) t_soil_v(l)%var_2d 5803 ENDDO 5804 5805 WRITE ( 14 ) 'lsm_start_index_h ' 5806 WRITE ( 14 ) surf_lsm_h%start_index 5807 WRITE ( 14 ) 'lsm_end_index_h ' 5808 WRITE ( 14 ) surf_lsm_h%end_index 5809 WRITE ( 14 ) 'm_soil_h ' 5810 WRITE ( 14 ) m_soil_h%var_2d 5811 5812 DO l = 0, 3 5813 WRITE ( 14 ) 'lsm_start_index_v ' 5814 WRITE ( 14 ) surf_lsm_v(l)%start_index 5815 WRITE ( 14 ) 'lsm_end_index_v ' 5816 WRITE ( 14 ) surf_lsm_v(l)%end_index 5817 WRITE( dum, '(I1)') l 5818 WRITE ( 14 ) 'm_soil_v(' // dum // ') ' 5819 WRITE ( 14 ) m_soil_v(l)%var_2d 5820 ENDDO 5821 5822 5823 WRITE ( 14 ) 'lsm_start_index_h ' 5824 WRITE ( 14 ) surf_lsm_h%start_index 5825 WRITE ( 14 ) 'lsm_end_index_h ' 5826 WRITE ( 14 ) surf_lsm_h%end_index 5827 WRITE ( 14 ) 'm_liq_h ' 5828 WRITE ( 14 ) m_liq_h%var_1d 5829 5830 DO l = 0, 3 5831 WRITE ( 14 ) 'lsm_start_index_v ' 5832 WRITE ( 14 ) surf_lsm_v(l)%start_index 5833 WRITE ( 14 ) 'lsm_end_index_v ' 5834 WRITE ( 14 ) surf_lsm_v(l)%end_index 5835 WRITE( dum, '(I1)') l 5836 WRITE ( 14 ) 'm_liq_v(' // dum // ') ' 5837 WRITE ( 14 ) m_liq_v(l)%var_1d 5838 ENDDO 5839 4376 5840 WRITE ( 14 ) '*** end lsm *** ' 4377 5841 4378 5842 ENDIF 4379 5843 4380 END SUBROUTINE lsm_ last_actions5844 END SUBROUTINE lsm_write_restart_data 4381 5845 4382 5846 … … 4399 5863 CHARACTER (LEN=20) :: field_char !< 4400 5864 4401 INTEGER(iwp) :: i !< 4402 INTEGER(iwp) :: k !< 4403 INTEGER(iwp) :: nxlc !< 4404 INTEGER(iwp) :: nxlf !< 4405 INTEGER(iwp) :: nxl_on_file !< 4406 INTEGER(iwp) :: nxrc !< 4407 INTEGER(iwp) :: nxrf !< 4408 INTEGER(iwp) :: nxr_on_file !< 4409 INTEGER(iwp) :: nync !< 4410 INTEGER(iwp) :: nynf !< 4411 INTEGER(iwp) :: nyn_on_file !< 4412 INTEGER(iwp) :: nysc !< 4413 INTEGER(iwp) :: nysf !< 4414 INTEGER(iwp) :: nys_on_file !< 4415 INTEGER(iwp) :: overlap_count !< 5865 INTEGER(iwp) :: i !< 5866 INTEGER(iwp) :: k !< 5867 INTEGER(iwp) :: l !< running index surface orientation 5868 INTEGER(iwp) :: ns_h_on_file_lsm !< number of horizontal surface elements (natural type) on file 5869 INTEGER(iwp) :: nxlc !< 5870 INTEGER(iwp) :: nxlf !< 5871 INTEGER(iwp) :: nxl_on_file !< 5872 INTEGER(iwp) :: nxrc !< 5873 INTEGER(iwp) :: nxrf !< 5874 INTEGER(iwp) :: nxr_on_file !< 5875 INTEGER(iwp) :: nync !< 5876 INTEGER(iwp) :: nynf !< 5877 INTEGER(iwp) :: nyn_on_file !< 5878 INTEGER(iwp) :: nysc !< 5879 INTEGER(iwp) :: nysf !< 5880 INTEGER(iwp) :: nys_on_file !< 5881 INTEGER(iwp) :: overlap_count !< 5882 5883 INTEGER(iwp) :: ns_v_on_file_lsm(0:3) !< number of vertical surface elements (natural type) on file 4416 5884 4417 5885 INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxlfa !< … … 4422 5890 INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_ya !< 4423 5891 5892 INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: start_index_on_file 5893 INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) :: end_index_on_file 5894 4424 5895 REAL(wp), & 4425 5896 DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: & … … 4434 5905 tmp_3d2 !< 4435 5906 4436 REAL(wp), & 4437 DIMENSION(1:surf_lsm_h%ns) :: & 4438 tmp_walltype_1d !< 4439 4440 REAL(wp), & 4441 DIMENSION(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) :: & 4442 tmp_walltype_2d !< 4443 4444 REAL(wp), & 4445 DIMENSION(nzb_soil:nzt_soil,1:surf_lsm_h%ns) :: & 4446 tmp_walltype_2d2 !< 4447 5907 TYPE(surf_type_lsm) :: tmp_walltype_h_1d !< temporary 1D array containing the respective surface variable stored on file, horizontal surfaces 5908 TYPE(surf_type_lsm) :: tmp_walltype_h_2d !< temporary 2D array containing the respective surface variable stored on file, horizontal surfaces 5909 TYPE(surf_type_lsm) :: tmp_walltype_h_2d2 !< temporary 2D array containing the respective surface variable stored on file, horizontal surfaces 5910 5911 TYPE(surf_type_lsm), DIMENSION(0:3) :: tmp_walltype_v_1d !< temporary 1D array containing the respective surface variable stored on file, vertical surfaces 5912 TYPE(surf_type_lsm), DIMENSION(0:3) :: tmp_walltype_v_2d !< temporary 2D array containing the respective surface variable stored on file, vertical surfaces 5913 TYPE(surf_type_lsm), DIMENSION(0:3) :: tmp_walltype_v_2d2 !< temporary 2D array containing the respective surface variable stored on file, vertical surfaces 4448 5914 4449 5915 IF ( initializing_actions == 'read_restart_data' ) THEN 5916 READ ( 13 ) field_char 5917 5918 ! 5919 !-- At first, determine the number of surface elements (with certain orientation) on file 5920 IF ( TRIM( field_char ) /= 'ns_h_on_file_lsm' ) THEN 5921 ! 5922 !-- Add a proper error message 5923 ENDIF 5924 READ ( 13 ) ns_h_on_file_lsm 5925 5926 READ ( 13 ) field_char 5927 IF ( TRIM( field_char ) /= 'ns_v_on_file_lsm' ) THEN 5928 ! 5929 !-- Add a proper error message 5930 ENDIF 5931 READ ( 13 ) ns_v_on_file_lsm 5932 5933 ! 5934 !-- Allocate temporary arrays to store surface data 5935 ALLOCATE( tmp_walltype_h_1d%var_1d(1:ns_h_on_file_lsm) ) 5936 ALLOCATE( tmp_walltype_h_2d%var_2d(nzb_soil:nzt_soil+1,1:ns_h_on_file_lsm) ) 5937 ALLOCATE( tmp_walltype_h_2d2%var_2d(nzb_soil:nzt_soil,1:ns_h_on_file_lsm) ) 5938 5939 DO l = 0, 3 5940 ALLOCATE( tmp_walltype_v_1d(l)%var_1d(1:ns_v_on_file_lsm(l)) ) 5941 ALLOCATE( tmp_walltype_v_2d(l)%var_2d(nzb_soil:nzt_soil+1,1:ns_v_on_file_lsm(l)) ) 5942 ALLOCATE( tmp_walltype_v_2d2(l)%var_2d(nzb_soil:nzt_soil,1:ns_v_on_file_lsm(l)) ) 5943 ENDDO 5944 4450 5945 READ ( 13 ) field_char 4451 5946 … … 4463 5958 nync = nynfa(i,k) + offset_ya(i,k) 4464 5959 5960 4465 5961 SELECT CASE ( TRIM( field_char ) ) 4466 5962 … … 4495 5991 ENDIF 4496 5992 IF ( k == 1 ) READ ( 13 ) tmp_2d 4497 ghf_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &5993 ghf_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 4498 5994 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 4499 4500 CASE ( 'm_liq' )4501 IF ( k == 1 ) READ ( 13 ) tmp_walltype_1d !tmp_2d4502 m_liq_h%var_1d(1:surf_lsm_h%ns) = &4503 tmp_walltype_1d(1:surf_lsm_h%ns)4504 5995 4505 5996 CASE ( 'lai_av' ) … … 4518 6009 m_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 4519 6010 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 4520 4521 CASE ( 'm_soil' )4522 IF ( k == 1 ) READ ( 13 ) tmp_walltype_2d2(:,:)4523 m_soil_h%var_2d(:,1:surf_lsm_h%ns) = &4524 tmp_walltype_2d2(:,1:surf_lsm_h%ns)4525 6011 4526 6012 CASE ( 'm_soil_av' ) … … 4556 6042 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 4557 6043 4558 CASE ( 't_soil' )4559 IF ( k == 1 ) READ ( 13 ) tmp_walltype_2d(:,:)4560 t_soil_h%var_2d(:,1:surf_lsm_h%ns) = &4561 tmp_walltype_2d(:,1:surf_lsm_h%ns)4562 4563 6044 CASE ( 't_soil_av' ) 4564 6045 IF ( .NOT. ALLOCATED( t_soil_av ) ) THEN … … 4569 6050 tmp_3d2(:,nysf-nbgp:nynf+nbgp, & 4570 6051 nxlf-nbgp:nxrf+nbgp) 6052 6053 CASE ( 'lsm_start_index_h', 'lsm_start_index_v' ) 6054 IF ( k == 1 ) & 6055 READ ( 13 ) start_index_on_file 6056 6057 CASE ( 'lsm_end_index_h', 'lsm_end_index_v' ) 6058 IF ( k == 1 ) & 6059 READ ( 13 ) end_index_on_file 6060 6061 CASE ( 't_soil_h' ) 6062 6063 IF ( k == 1 ) THEN 6064 IF ( .NOT. ALLOCATED( t_soil_h%var_2d ) ) & 6065 ALLOCATE( t_soil_h%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) ) 6066 READ ( 13 ) tmp_walltype_h_2d%var_2d 6067 ENDIF 6068 CALL surface_restore_elements( & 6069 t_soil_h%var_2d, & 6070 tmp_walltype_h_2d%var_2d, & 6071 surf_lsm_h%start_index, & 6072 start_index_on_file, & 6073 end_index_on_file, & 6074 nxlc, nysc, & 6075 nxlf, nxrf, nysf, nynf, & 6076 nys_on_file, nyn_on_file, & 6077 nxl_on_file,nxr_on_file ) 6078 6079 CASE ( 't_soil_v(0)' ) 6080 6081 IF ( k == 1 ) THEN 6082 IF ( .NOT. ALLOCATED( t_soil_v(0)%var_2d ) ) & 6083 ALLOCATE( t_soil_v(0)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(0)%ns) ) 6084 READ ( 13 ) tmp_walltype_v_2d(0)%var_2d 6085 ENDIF 6086 CALL surface_restore_elements( & 6087 t_soil_v(0)%var_2d, & 6088 tmp_walltype_v_2d(0)%var_2d, & 6089 surf_lsm_v(0)%start_index, & 6090 start_index_on_file, & 6091 end_index_on_file, & 6092 nxlc, nysc, & 6093 nxlf, nxrf, nysf, nynf, & 6094 nys_on_file, nyn_on_file, & 6095 nxl_on_file,nxr_on_file ) 6096 6097 CASE ( 't_soil_v(1)' ) 6098 6099 IF ( k == 1 ) THEN 6100 IF ( .NOT. ALLOCATED( t_soil_v(1)%var_2d ) ) & 6101 ALLOCATE( t_soil_v(1)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(1)%ns) ) 6102 READ ( 13 ) tmp_walltype_v_2d(1)%var_2d 6103 ENDIF 6104 CALL surface_restore_elements( & 6105 t_soil_v(1)%var_2d, & 6106 tmp_walltype_v_2d(1)%var_2d, & 6107 surf_lsm_v(1)%start_index, & 6108 start_index_on_file, & 6109 end_index_on_file, & 6110 nxlc, nysc, & 6111 nxlf, nxrf, nysf, nynf, & 6112 nys_on_file, nyn_on_file, & 6113 nxl_on_file,nxr_on_file ) 6114 6115 CASE ( 't_soil_v(2)' ) 6116 6117 IF ( k == 1 ) THEN 6118 IF ( .NOT. ALLOCATED( t_soil_v(2)%var_2d ) ) & 6119 ALLOCATE( t_soil_v(2)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(2)%ns) ) 6120 READ ( 13 ) tmp_walltype_v_2d(2)%var_2d 6121 ENDIF 6122 CALL surface_restore_elements( & 6123 t_soil_v(2)%var_2d, & 6124 tmp_walltype_v_2d(2)%var_2d, & 6125 surf_lsm_v(2)%start_index, & 6126 start_index_on_file, & 6127 end_index_on_file, & 6128 nxlc, nysc, & 6129 nxlf, nxrf, nysf, nynf, & 6130 nys_on_file, nyn_on_file, & 6131 nxl_on_file,nxr_on_file ) 6132 6133 CASE ( 't_soil_v(3)' ) 6134 6135 IF ( k == 1 ) THEN 6136 IF ( .NOT. ALLOCATED( t_soil_v(3)%var_2d ) ) & 6137 ALLOCATE( t_soil_v(1)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(3)%ns) ) 6138 READ ( 13 ) tmp_walltype_v_2d(3)%var_2d 6139 ENDIF 6140 CALL surface_restore_elements( & 6141 t_soil_v(3)%var_2d, & 6142 tmp_walltype_v_2d(3)%var_2d, & 6143 surf_lsm_v(3)%start_index, & 6144 start_index_on_file, & 6145 end_index_on_file, & 6146 nxlc, nysc, & 6147 nxlf, nxrf, nysf, nynf, & 6148 nys_on_file, nyn_on_file, & 6149 nxl_on_file,nxr_on_file ) 6150 6151 CASE ( 'm_soil_h' ) 6152 6153 IF ( k == 1 ) THEN 6154 IF ( .NOT. ALLOCATED( m_soil_h%var_2d ) ) & 6155 ALLOCATE( m_soil_h%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) ) 6156 READ ( 13 ) tmp_walltype_h_2d2%var_2d 6157 ENDIF 6158 CALL surface_restore_elements( & 6159 m_soil_h%var_2d, & 6160 tmp_walltype_h_2d2%var_2d, & 6161 surf_lsm_h%start_index, & 6162 start_index_on_file, & 6163 end_index_on_file, & 6164 nxlc, nysc, & 6165 nxlf, nxrf, nysf, nynf, & 6166 nys_on_file, nyn_on_file, & 6167 nxl_on_file,nxr_on_file ) 6168 6169 CASE ( 'm_soil_v(0)' ) 6170 6171 IF ( k == 1 ) THEN 6172 IF ( .NOT. ALLOCATED( m_soil_v(0)%var_2d ) ) & 6173 ALLOCATE( m_soil_v(0)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(0)%ns) ) 6174 READ ( 13 ) tmp_walltype_v_2d2(0)%var_2d 6175 ENDIF 6176 CALL surface_restore_elements( & 6177 m_soil_v(0)%var_2d, & 6178 tmp_walltype_v_2d2(0)%var_2d, & 6179 surf_lsm_v(0)%start_index, & 6180 start_index_on_file, & 6181 end_index_on_file, & 6182 nxlc, nysc, & 6183 nxlf, nxrf, nysf, nynf, & 6184 nys_on_file, nyn_on_file, & 6185 nxl_on_file,nxr_on_file ) 6186 6187 CASE ( 'm_soil_v(1)' ) 6188 6189 IF ( k == 1 ) THEN 6190 IF ( .NOT. ALLOCATED( m_soil_v(1)%var_2d ) ) & 6191 ALLOCATE( m_soil_v(1)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(1)%ns) ) 6192 READ ( 13 ) tmp_walltype_v_2d2(1)%var_2d 6193 ENDIF 6194 CALL surface_restore_elements( & 6195 m_soil_v(1)%var_2d, & 6196 tmp_walltype_v_2d2(1)%var_2d, & 6197 surf_lsm_v(1)%start_index, & 6198 start_index_on_file, & 6199 end_index_on_file, & 6200 nxlc, nysc, & 6201 nxlf, nxrf, nysf, nynf, & 6202 nys_on_file, nyn_on_file, & 6203 nxl_on_file,nxr_on_file ) 6204 6205 6206 CASE ( 'm_soil_v(2)' ) 6207 6208 IF ( k == 1 ) THEN 6209 IF ( .NOT. ALLOCATED( m_soil_v(2)%var_2d ) ) & 6210 ALLOCATE( m_soil_v(2)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(2)%ns) ) 6211 READ ( 13 ) tmp_walltype_v_2d2(2)%var_2d 6212 ENDIF 6213 CALL surface_restore_elements( & 6214 m_soil_v(2)%var_2d, & 6215 tmp_walltype_v_2d2(2)%var_2d, & 6216 surf_lsm_v(2)%start_index, & 6217 start_index_on_file, & 6218 end_index_on_file, & 6219 nxlc, nysc, & 6220 nxlf, nxrf, nysf, nynf, & 6221 nys_on_file, nyn_on_file, & 6222 nxl_on_file,nxr_on_file ) 6223 6224 6225 CASE ( 'm_soil_v(3)' ) 6226 6227 IF ( k == 1 ) THEN 6228 IF ( .NOT. ALLOCATED( m_soil_v(3)%var_2d ) ) & 6229 ALLOCATE( m_soil_v(1)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(3)%ns) ) 6230 READ ( 13 ) tmp_walltype_v_2d2(3)%var_2d 6231 ENDIF 6232 CALL surface_restore_elements( & 6233 m_soil_v(3)%var_2d, & 6234 tmp_walltype_v_2d2(3)%var_2d, & 6235 surf_lsm_v(3)%start_index, & 6236 start_index_on_file, & 6237 end_index_on_file, & 6238 nxlc, nysc, & 6239 nxlf, nxrf, nysf, nynf, & 6240 nys_on_file, nyn_on_file, & 6241 nxl_on_file,nxr_on_file ) 6242 6243 6244 CASE ( 'm_liq_h' ) 6245 6246 IF ( k == 1 ) THEN 6247 IF ( .NOT. ALLOCATED( m_liq_h%var_1d ) ) & 6248 ALLOCATE( m_liq_h%var_1d(1:surf_lsm_h%ns) ) 6249 READ ( 13 ) tmp_walltype_h_1d%var_1d 6250 ENDIF 6251 CALL surface_restore_elements( & 6252 m_liq_h%var_1d, & 6253 tmp_walltype_h_1d%var_1d, & 6254 surf_lsm_h%start_index, & 6255 start_index_on_file, & 6256 end_index_on_file, & 6257 nxlc, nysc, & 6258 nxlf, nxrf, nysf, nynf, & 6259 nys_on_file, nyn_on_file, & 6260 nxl_on_file,nxr_on_file ) 6261 6262 6263 CASE ( 'm_liq_v(0)' ) 6264 6265 IF ( k == 1 ) THEN 6266 IF ( .NOT. ALLOCATED( m_liq_v(0)%var_1d ) ) & 6267 ALLOCATE( m_liq_v(0)%var_1d(1:surf_lsm_v(0)%ns) ) 6268 READ ( 13 ) tmp_walltype_v_1d(0)%var_1d 6269 ENDIF 6270 CALL surface_restore_elements( & 6271 m_liq_v(0)%var_1d, & 6272 tmp_walltype_v_1d(0)%var_1d, & 6273 surf_lsm_v(0)%start_index, & 6274 start_index_on_file, & 6275 end_index_on_file, & 6276 nxlc, nysc, & 6277 nxlf, nxrf, nysf, nynf, & 6278 nys_on_file, nyn_on_file, & 6279 nxl_on_file,nxr_on_file ) 6280 6281 6282 CASE ( 'm_liq_v(1)' ) 6283 6284 IF ( k == 1 ) THEN 6285 IF ( .NOT. ALLOCATED( m_liq_v(1)%var_1d ) ) & 6286 ALLOCATE( m_liq_v(1)%var_1d(1:surf_lsm_v(1)%ns) ) 6287 READ ( 13 ) tmp_walltype_v_1d(1)%var_1d 6288 ENDIF 6289 CALL surface_restore_elements( & 6290 m_liq_v(1)%var_1d, & 6291 tmp_walltype_v_1d(1)%var_1d, & 6292 surf_lsm_v(1)%start_index, & 6293 start_index_on_file, & 6294 end_index_on_file, & 6295 nxlc, nysc, & 6296 nxlf, nxrf, nysf, nynf, & 6297 nys_on_file, nyn_on_file, & 6298 nxl_on_file,nxr_on_file ) 6299 6300 6301 CASE ( 'm_liq_v(2)' ) 6302 6303 IF ( k == 1 ) THEN 6304 IF ( .NOT. ALLOCATED( m_liq_v(2)%var_1d ) ) & 6305 ALLOCATE( m_liq_v(2)%var_1d(1:surf_lsm_v(2)%ns) ) 6306 READ ( 13 ) tmp_walltype_v_1d(2)%var_1d 6307 ENDIF 6308 CALL surface_restore_elements( & 6309 m_liq_v(2)%var_1d, & 6310 tmp_walltype_v_1d(2)%var_1d, & 6311 surf_lsm_v(2)%start_index, & 6312 start_index_on_file, & 6313 end_index_on_file, & 6314 nxlc, nysc, & 6315 nxlf, nxrf, nysf, nynf, & 6316 nys_on_file, nyn_on_file, & 6317 nxl_on_file,nxr_on_file ) 6318 6319 CASE ( 'm_liq_v(3)' ) 6320 6321 IF ( k == 1 ) THEN 6322 IF ( .NOT. ALLOCATED( m_liq_v(3)%var_1d ) ) & 6323 ALLOCATE( m_liq_v(3)%var_1d(1:surf_lsm_v(3)%ns) ) 6324 READ ( 13 ) tmp_walltype_v_1d(3)%var_1d 6325 ENDIF 6326 CALL surface_restore_elements( & 6327 m_liq_v(3)%var_1d, & 6328 tmp_walltype_v_1d(3)%var_1d, & 6329 surf_lsm_v(3)%start_index, & 6330 start_index_on_file, & 6331 end_index_on_file, & 6332 nxlc, nysc, & 6333 nxlf, nxrf, nysf, nynf, & 6334 nys_on_file, nyn_on_file, & 6335 nxl_on_file,nxr_on_file ) 4571 6336 4572 6337 … … 4645 6410 !-- Charnock relation is used. 4646 6411 surf_lsm_h%z0(m) = ( 0.11_wp * molecular_viscosity / & 4647 surf_lsm_h%us(m) ) &6412 surf_lsm_h%us(m) ) & 4648 6413 + ( alpha_ch * surf_lsm_h%us(m)**2 / g ) 4649 6414
Note: See TracChangeset
for help on using the changeset viewer.
