Changeset 1817 for palm/trunk
- Timestamp:
- Apr 6, 2016 3:44:20 PM (9 years ago)
- Location:
- palm/trunk
- Files:
-
- 7 edited
- 1 moved
-
SCRIPTS/mbuild (modified) (3 diffs)
-
SOURCE/Makefile (modified) (13 diffs)
-
SOURCE/check_parameters.f90 (modified) (11 diffs)
-
SOURCE/header.f90 (modified) (5 diffs)
-
SOURCE/init_3d_model.f90 (modified) (5 diffs)
-
SOURCE/land_surface_model_mod.f90 (moved) (moved from palm/trunk/SOURCE/land_surface_model.f90) (14 diffs)
-
SOURCE/package_parin.f90 (modified) (4 diffs)
-
SOURCE/parin.f90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SCRIPTS/mbuild
r1805 r1817 22 22 # Current revisions: 23 23 # ------------------ 24 # 24 # Bugfix: Makefile_check removed (parameter file check) 25 25 # 26 26 # Former revisions: … … 427 427 cd $local_source_path 428 428 printf "\n\n *** tar of makefile and source files in $local_source_path" 429 430 cat Makefile_check|while read line431 do432 line=$(echo $line|grep SOURCES)433 if [[ $line == *"SOURCES"* ]]434 then435 line=$(echo $line|sed 's/SOURCES = //g')436 break437 fi438 done439 440 tar -cf ${mainprog}_sources_check.tar Makefile_check $line441 429 printf "\n" 442 430 fi … … 1481 1469 # FINAL WORK 1482 1470 rm -f hosts_found_in_config_file 1483 rm -f ${local_source_path}/${mainprog}_sources_check.tar -
palm/trunk/SOURCE/Makefile
r1809 r1817 20 20 # Current revisions: 21 21 # ------------------ 22 # 22 # Renamed land_surface_model to land_surface_model_mod. 23 23 # 24 24 # Former revisions: … … 259 259 init_grid.f90 init_masks.f90 init_ocean.f90 \ 260 260 init_pegrid.f90 init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \ 261 interaction_droplets_ptq.f90 land_surface_model .f90 \261 interaction_droplets_ptq.f90 land_surface_model_mod.f90 \ 262 262 local_stop.f90 local_system.f90 local_tremain.f90 \ 263 263 local_tremain_ini.f90 lpm.f90 lpm_advec.f90 lpm_boundary_conds.f90 \ … … 340 340 advec_w_pw.o: modules.o mod_kinds.o 341 341 advec_w_up.o: modules.o mod_kinds.o 342 average_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model .o\342 average_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model_mod.o\ 343 343 radiation_model.o 344 344 boundary_conds.o: modules.o mod_kinds.o … … 350 350 check_for_restart.o: modules.o mod_kinds.o pmc_interface.o 351 351 check_open.o: modules.o mod_kinds.o mod_particle_attributes.o netcdf_interface.o 352 check_parameters.o: modules.o mod_kinds.o land_surface_model .o \352 check_parameters.o: modules.o mod_kinds.o land_surface_model_mod.o \ 353 353 netcdf_interface.o plant_canopy_model.o pmc_interface.o radiation_model.o \ 354 354 subsidence.o … … 370 370 data_output_tseries.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o 371 371 data_output_2d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o\ 372 netcdf_interface.o land_surface_model .o radiation_model.o372 netcdf_interface.o land_surface_model_mod.o radiation_model.o 373 373 data_output_3d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o\ 374 netcdf_interface.o land_surface_model .o374 netcdf_interface.o land_surface_model_mod.o 375 375 diffusion_e.o: modules.o mod_kinds.o 376 376 diffusion_s.o: modules.o mod_kinds.o … … 386 386 exchange_horiz_2d.o: modules.o cpulog.o mod_kinds.o 387 387 fft_xy.o: cuda_fft_interfaces.o modules.o mod_kinds.o singleton.o temperton_fft.o 388 flow_statistics.o: modules.o cpulog.o mod_kinds.o land_surface_model .o\388 flow_statistics.o: modules.o cpulog.o mod_kinds.o land_surface_model_mod.o\ 389 389 netcdf_interface.o radiation_model.o 390 390 global_min_max.o: modules.o mod_kinds.o 391 header.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o land_surface_model .o\391 header.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o land_surface_model_mod.o\ 392 392 plant_canopy_model.o pmc_handle_communicator.o pmc_interface.o \ 393 393 radiation_model.o spectrum.o subsidence.o … … 395 395 inflow_turbulence.o: modules.o cpulog.o mod_kinds.o 396 396 init_1d_model.o: modules.o mod_kinds.o 397 init_3d_model.o: modules.o mod_kinds.o advec_ws.o cpulog.o land_surface_model .o \397 init_3d_model.o: modules.o mod_kinds.o advec_ws.o cpulog.o land_surface_model_mod.o \ 398 398 lpm_init.o ls_forcing.o netcdf_interface.o plant_canopy_model.o \ 399 399 radiation_model.o random_function.o random_generator_parallel.o \ … … 411 411 init_slope.o: modules.o mod_kinds.o 412 412 interaction_droplets_ptq.o: modules.o mod_kinds.o 413 land_surface_model .o: modules.o mod_kinds.o radiation_model.o413 land_surface_model_mod.o: modules.o mod_kinds.o radiation_model.o 414 414 local_stop.o: modules.o mod_kinds.o pmc_interface.o 415 415 local_tremain.o: modules.o cpulog.o mod_kinds.o … … 452 452 mod_particle_attributes.o: mod_particle_attributes.f90 mod_kinds.o 453 453 netcdf_interface.o: netcdf_interface.f90 modules.o mod_kinds.o \ 454 land_surface_model .o spectrum.o454 land_surface_model_mod.o spectrum.o 455 455 nudging.o: modules.o cpulog.o mod_kinds.o 456 package_parin.o: modules.o mod_kinds.o land_surface_model.o\456 package_parin.o: modules.o mod_kinds.o \ 457 457 plant_canopy_model.o radiation_model.o spectrum.o 458 458 palm.o: modules.o cpulog.o ls_forcing.o mod_kinds.o nudging.o\ 459 459 pmc_interface.o surface_layer_fluxes.o 460 parin.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o \460 parin.o: modules.o cpulog.o land_surface_model_mod.o mod_kinds.o netcdf_interface.o \ 461 461 pmc_interface.o progress_bar.o 462 462 plant_canopy_model.o: modules.o mod_kinds.o … … 487 487 random_gauss.o: mod_kinds.o random_function.o random_generator_parallel.o 488 488 random_generator_parallel.o: mod_kinds.o 489 read_3d_binary.o: modules.o cpulog.o mod_kinds.o land_surface_model .o\489 read_3d_binary.o: modules.o cpulog.o mod_kinds.o land_surface_model_mod.o\ 490 490 radiation_model.o random_function.o\ 491 491 random_generator_parallel.o … … 497 497 spectrum.o: spectrum.f90 modules.o mod_kinds.o cpulog.o fft_xy.o 498 498 subsidence.o: modules.o mod_kinds.o 499 sum_up_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model .o\499 sum_up_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model_mod.o\ 500 500 radiation_model.o 501 501 surface_coupler.o: modules.o cpulog.o mod_kinds.o 502 surface_layer_fluxes.o: modules.o mod_kinds.o land_surface_model .o503 swap_timelevel.o: modules.o cpulog.o mod_kinds.o land_surface_model .o\502 surface_layer_fluxes.o: modules.o mod_kinds.o land_surface_model_mod.o 503 swap_timelevel.o: modules.o cpulog.o mod_kinds.o land_surface_model_mod.o\ 504 504 pmc_interface.o 505 505 temperton_fft.o: modules.o mod_kinds.o 506 506 time_integration.o: modules.o advec_ws.o buoyancy.o calc_mean_profile.o \ 507 cpulog.o interaction_droplets_ptq.o land_surface_model .o \507 cpulog.o interaction_droplets_ptq.o land_surface_model_mod.o \ 508 508 ls_forcing.o mod_kinds.o nudging.o pmc_interface.o production_e.o \ 509 509 prognostic_equations.o progress_bar.o radiation_model.o \ … … 531 531 user_init_3d_model.o: modules.o mod_kinds.o user_module.o 532 532 user_init_grid.o: modules.o mod_kinds.o user_module.o 533 user_init_land_surface.o: modules.o mod_kinds.o user_module.o land_surface_model .o533 user_init_land_surface.o: modules.o mod_kinds.o user_module.o land_surface_model_mod.o 534 534 user_init_plant_canopy.o: modules.o mod_kinds.o user_module.o plant_canopy_model.o 535 535 user_init_radiation.o: modules.o mod_kinds.o user_module.o radiation_model.o … … 544 544 user_statistics.o: modules.o mod_kinds.o netcdf_interface.o user_module.o 545 545 wall_fluxes.o: modules.o mod_kinds.o 546 write_3d_binary.o: modules.o cpulog.o mod_kinds.o land_surface_model .o\546 write_3d_binary.o: modules.o cpulog.o mod_kinds.o land_surface_model_mod.o\ 547 547 radiation_model.o random_function.o\ 548 548 random_generator_parallel.o -
palm/trunk/SOURCE/check_parameters.f90
r1812 r1817 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Moved checks for land_surface model to the respective module 22 22 ! 23 23 ! Former revisions: … … 330 330 USE grid_variables 331 331 USE indices 332 USE land_surface_model_mod 332 USE land_surface_model_mod, & 333 ONLY: land_surface, lsm_check_data_output, lsm_check_data_output_pr, & 334 lsm_check_parameters, nzs, zs 333 335 USE kinds 334 336 USE model_1d … … 1035 1037 ENDIF 1036 1038 1037 1038 IF ( land_surface ) THEN 1039 1040 ! 1041 !-- Dirichlet boundary conditions are required as the surface fluxes are 1042 !-- calculated from the temperature/humidity gradients in the land surface 1043 !-- model 1044 IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' ) THEN 1045 message_string = 'lsm requires setting of'// & 1046 'bc_pt_b = "dirichlet" and '// & 1047 'bc_q_b = "dirichlet"' 1048 CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 ) 1049 ENDIF 1050 1051 IF ( .NOT. constant_flux_layer ) THEN 1052 message_string = 'lsm requires '// & 1053 'constant_flux_layer = .T.' 1054 CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 ) 1055 ENDIF 1056 1057 IF ( topography /= 'flat' ) THEN 1058 message_string = 'lsm cannot be used ' // & 1059 'in combination with topography /= "flat"' 1060 CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 ) 1061 ENDIF 1062 1063 IF ( ( veg_type == 14 .OR. veg_type == 15 ) .AND. & 1064 most_method == 'lookup' ) THEN 1065 WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ', & 1066 'allowed in combination with ', & 1067 'most_method = ', most_method 1068 CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 ) 1069 ENDIF 1070 1071 IF ( veg_type == 0 ) THEN 1072 IF ( SUM( root_fraction ) /= 1.0_wp ) THEN 1073 message_string = 'veg_type = 0 (user_defined)'// & 1074 'requires setting of root_fraction(0:3)'// & 1075 '/= 9999999.9 and SUM(root_fraction) = 1' 1076 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 1077 ENDIF 1078 1079 IF ( min_canopy_resistance == 9999999.9_wp ) THEN 1080 message_string = 'veg_type = 0 (user defined)'// & 1081 'requires setting of min_canopy_resistance'// & 1082 '/= 9999999.9' 1083 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 1084 ENDIF 1085 1086 IF ( leaf_area_index == 9999999.9_wp ) THEN 1087 message_string = 'veg_type = 0 (user_defined)'// & 1088 'requires setting of leaf_area_index'// & 1089 '/= 9999999.9' 1090 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 1091 ENDIF 1092 1093 IF ( vegetation_coverage == 9999999.9_wp ) THEN 1094 message_string = 'veg_type = 0 (user_defined)'// & 1095 'requires setting of vegetation_coverage'// & 1096 '/= 9999999.9' 1097 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 1098 ENDIF 1099 1100 IF ( canopy_resistance_coefficient == 9999999.9_wp) THEN 1101 message_string = 'veg_type = 0 (user_defined)'// & 1102 'requires setting of'// & 1103 'canopy_resistance_coefficient /= 9999999.9' 1104 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 1105 ENDIF 1106 1107 IF ( lambda_surface_stable == 9999999.9_wp ) THEN 1108 message_string = 'veg_type = 0 (user_defined)'// & 1109 'requires setting of lambda_surface_stable'// & 1110 '/= 9999999.9' 1111 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 1112 ENDIF 1113 1114 IF ( lambda_surface_unstable == 9999999.9_wp ) THEN 1115 message_string = 'veg_type = 0 (user_defined)'// & 1116 'requires setting of lambda_surface_unstable'// & 1117 '/= 9999999.9' 1118 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 1119 ENDIF 1120 1121 IF ( f_shortwave_incoming == 9999999.9_wp ) THEN 1122 message_string = 'veg_type = 0 (user_defined)'// & 1123 'requires setting of f_shortwave_incoming'// & 1124 '/= 9999999.9' 1125 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 1126 ENDIF 1127 1128 IF ( z0_eb == 9999999.9_wp ) THEN 1129 message_string = 'veg_type = 0 (user_defined)'// & 1130 'requires setting of z0_eb'// & 1131 '/= 9999999.9' 1132 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 1133 ENDIF 1134 1135 IF ( z0h_eb == 9999999.9_wp ) THEN 1136 message_string = 'veg_type = 0 (user_defined)'// & 1137 'requires setting of z0h_eb'// & 1138 '/= 9999999.9' 1139 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 1140 ENDIF 1141 1142 1143 ENDIF 1144 1145 IF ( soil_type == 0 ) THEN 1146 1147 IF ( alpha_vangenuchten == 9999999.9_wp ) THEN 1148 message_string = 'soil_type = 0 (user_defined)'// & 1149 'requires setting of alpha_vangenuchten'// & 1150 '/= 9999999.9' 1151 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1152 ENDIF 1153 1154 IF ( l_vangenuchten == 9999999.9_wp ) THEN 1155 message_string = 'soil_type = 0 (user_defined)'// & 1156 'requires setting of l_vangenuchten'// & 1157 '/= 9999999.9' 1158 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1159 ENDIF 1160 1161 IF ( n_vangenuchten == 9999999.9_wp ) THEN 1162 message_string = 'soil_type = 0 (user_defined)'// & 1163 'requires setting of n_vangenuchten'// & 1164 '/= 9999999.9' 1165 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1166 ENDIF 1167 1168 IF ( hydraulic_conductivity == 9999999.9_wp ) THEN 1169 message_string = 'soil_type = 0 (user_defined)'// & 1170 'requires setting of hydraulic_conductivity'// & 1171 '/= 9999999.9' 1172 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1173 ENDIF 1174 1175 IF ( saturation_moisture == 9999999.9_wp ) THEN 1176 message_string = 'soil_type = 0 (user_defined)'// & 1177 'requires setting of saturation_moisture'// & 1178 '/= 9999999.9' 1179 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1180 ENDIF 1181 1182 IF ( field_capacity == 9999999.9_wp ) THEN 1183 message_string = 'soil_type = 0 (user_defined)'// & 1184 'requires setting of field_capacity'// & 1185 '/= 9999999.9' 1186 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1187 ENDIF 1188 1189 IF ( wilting_point == 9999999.9_wp ) THEN 1190 message_string = 'soil_type = 0 (user_defined)'// & 1191 'requires setting of wilting_point'// & 1192 '/= 9999999.9' 1193 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1194 ENDIF 1195 1196 IF ( residual_moisture == 9999999.9_wp ) THEN 1197 message_string = 'soil_type = 0 (user_defined)'// & 1198 'requires setting of residual_moisture'// & 1199 '/= 9999999.9' 1200 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 1201 ENDIF 1202 1203 ENDIF 1204 1205 IF ( .NOT. radiation ) THEN 1206 message_string = 'lsm requires '// & 1207 'radiation = .T.' 1208 CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 ) 1209 ENDIF 1210 1211 END IF 1039 ! 1040 !-- When land surface model is used, peform addtional checks 1041 IF ( land_surface ) CALL lsm_check_parameters 1042 1212 1043 1213 1044 IF ( radiation ) THEN … … 3144 2975 ENDIF 3145 2976 3146 CASE ( 't_soil', '#t_soil' )3147 IF ( .NOT. land_surface ) THEN3148 message_string = 'data_output_pr = ' // &3149 TRIM( data_output_pr(i) ) // ' is not imp' // &3150 'lemented for land_surface = .FALSE.'3151 CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )3152 ELSE3153 dopr_index(i) = 893154 dopr_unit(i) = 'K'3155 hom(0:nzs-1,2,89,:) = SPREAD( - zs, 2, statistic_regions+1 )3156 IF ( data_output_pr(i)(1:1) == '#' ) THEN3157 dopr_initial_index(i) = 903158 hom(0:nzs-1,2,90,:) = SPREAD( - zs, 2, statistic_regions+1 )3159 data_output_pr(i) = data_output_pr(i)(2:)3160 ENDIF3161 ENDIF3162 3163 CASE ( 'm_soil', '#m_soil' )3164 IF ( .NOT. land_surface ) THEN3165 message_string = 'data_output_pr = ' // &3166 TRIM( data_output_pr(i) ) // ' is not imp' // &3167 'lemented for land_surface = .FALSE.'3168 CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )3169 ELSE3170 dopr_index(i) = 913171 dopr_unit(i) = 'm3/m3'3172 hom(0:nzs-1,2,91,:) = SPREAD( - zs, 2, statistic_regions+1 )3173 IF ( data_output_pr(i)(1:1) == '#' ) THEN3174 dopr_initial_index(i) = 923175 hom(0:nzs-1,2,92,:) = SPREAD( - zs, 2, statistic_regions+1 )3176 data_output_pr(i) = data_output_pr(i)(2:)3177 ENDIF3178 ENDIF3179 3180 2977 CASE ( 'rad_net' ) 3181 2978 IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )& … … 3306 3103 CASE DEFAULT 3307 3104 3308 CALL user_check_data_output_pr( data_output_pr(i), i, unit ) 3105 CALL lsm_check_data_output_pr( data_output_pr(i), i, unit, dopr_unit(i) ) 3106 3107 IF ( unit == 'illegal' ) THEN 3108 CALL user_check_data_output_pr( data_output_pr(i), i, unit ) 3109 ENDIF 3309 3110 3310 3111 IF ( unit == 'illegal' ) THEN … … 3394 3195 unit = 'K' 3395 3196 3396 CASE ( 'm_soil' )3397 IF ( .NOT. land_surface ) THEN3398 message_string = 'output of "' // TRIM( var ) // '" requi' // &3399 'land_surface = .TRUE.'3400 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )3401 ENDIF3402 unit = 'm3/m3'3403 3404 3197 CASE ( 'nr' ) 3405 3198 IF ( .NOT. cloud_physics ) THEN … … 3536 3329 unit = 'psu' 3537 3330 3538 CASE ( 't_soil' ) 3539 IF ( .NOT. land_surface ) THEN 3540 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3541 'land_surface = .TRUE.' 3542 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 3543 ENDIF 3544 unit = 'K' 3545 3546 3547 CASE ( 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'lai*', 'lwp*', & 3548 'm_liq_eb*', 'ol*', 'pra*', 'prr*', 'qsws*', 'qsws_eb*', & 3549 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*', 'rad_net*', & 3331 CASE ( 'p', 'pt', 'u', 'v', 'w' ) 3332 IF ( TRIM( var ) == 'p' ) unit = 'Pa' 3333 IF ( TRIM( var ) == 'pt' ) unit = 'K' 3334 IF ( TRIM( var ) == 'u' ) unit = 'm/s' 3335 IF ( TRIM( var ) == 'v' ) unit = 'm/s' 3336 IF ( TRIM( var ) == 'w' ) unit = 'm/s' 3337 CONTINUE 3338 3339 CASE ( 'lai*', 'lwp*', 'ol*', 'pra*', 'prr*', 'qsws*', 'rad_net*', & 3550 3340 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*', & 3551 'r_a*', 'r_s*', 'shf*', 'shf_eb*', 't*', 'u*', 'z0*', 'z0h*', & 3552 'z0q*' ) 3341 'shf*', 't*', 'u*', 'z0*', 'z0h*', 'z0q*' ) 3553 3342 IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN 3554 3343 message_string = 'illegal value for data_output: "' // & … … 3569 3358 ENDIF 3570 3359 ENDIF 3571 3572 IF ( TRIM( var ) == 'c_liq*' .AND. .NOT. land_surface ) THEN3573 message_string = 'output of "' // TRIM( var ) // '" requi' // &3574 'res land_surface = .TRUE.'3575 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )3576 ENDIF3577 IF ( TRIM( var ) == 'c_soil*' .AND. .NOT. land_surface ) THEN3578 message_string = 'output of "' // TRIM( var ) // '" requi' // &3579 'res land_surface = .TRUE.'3580 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )3581 ENDIF3582 IF ( TRIM( var ) == 'c_veg*' .AND. .NOT. land_surface ) THEN3583 message_string = 'output of "' // TRIM( var ) // '" requi' // &3584 'res land_surface = .TRUE.'3585 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )3586 ENDIF3587 IF ( TRIM( var ) == 'ghf_eb*' .AND. .NOT. land_surface ) THEN3588 message_string = 'output of "' // TRIM( var ) // '" requi' // &3589 'res land_surface = .TRUE.'3590 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )3591 ENDIF3592 3360 IF ( TRIM( var ) == 'lai*' .AND. .NOT. land_surface ) THEN 3593 3361 message_string = 'output of "' // TRIM( var ) // '" requi' // & … … 3600 3368 CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 ) 3601 3369 ENDIF 3602 IF ( TRIM( var ) == 'm_liq_eb*' .AND. .NOT. land_surface ) THEN3603 message_string = 'output of "' // TRIM( var ) // '" requi' // &3604 'res land_surface = .TRUE.'3605 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )3606 ENDIF3607 3370 IF ( TRIM( var ) == 'pra*' .AND. .NOT. precipitation ) THEN 3608 3371 message_string = 'output of "' // TRIM( var ) // '" requi' // & … … 3625 3388 CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 ) 3626 3389 ENDIF 3627 IF ( TRIM( var ) == 'qsws_eb*' .AND. .NOT. land_surface ) THEN 3628 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3629 'res land_surface = .TRUE.' 3630 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 3631 ENDIF 3632 IF ( TRIM( var ) == 'qsws_liq_eb*' .AND. .NOT. land_surface ) & 3633 THEN 3634 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3635 'res land_surface = .TRUE.' 3636 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 3637 ENDIF 3638 IF ( TRIM( var ) == 'qsws_soil_eb*' .AND. .NOT. land_surface ) & 3639 THEN 3640 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3641 'res land_surface = .TRUE.' 3642 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 3643 ENDIF 3644 IF ( TRIM( var ) == 'qsws_veg_eb*' .AND. .NOT. land_surface ) & 3645 THEN 3646 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3647 'res land_surface = .TRUE.' 3648 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 3649 ENDIF 3650 IF ( TRIM( var ) == 'r_a*' .AND. .NOT. land_surface ) & 3651 THEN 3652 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3653 'res land_surface = .TRUE.' 3654 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 3655 ENDIF 3656 IF ( TRIM( var ) == 'r_s*' .AND. .NOT. land_surface ) & 3657 THEN 3658 message_string = 'output of "' // TRIM( var ) // '" requi' // & 3659 'res land_surface = .TRUE.' 3660 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 3661 ENDIF 3662 3663 IF ( TRIM( var ) == 'c_liq*' ) unit = 'none' 3664 IF ( TRIM( var ) == 'c_soil*') unit = 'none' 3665 IF ( TRIM( var ) == 'c_veg*' ) unit = 'none' 3666 IF ( TRIM( var ) == 'ghf_eb*') unit = 'W/m2' 3390 3667 3391 IF ( TRIM( var ) == 'lai*' ) unit = 'none' 3668 3392 IF ( TRIM( var ) == 'lwp*' ) unit = 'kg/m2' 3669 IF ( TRIM( var ) == 'm_liq_eb*' ) unit = 'm'3670 3393 IF ( TRIM( var ) == 'ol*' ) unit = 'm' 3671 3394 IF ( TRIM( var ) == 'pra*' ) unit = 'mm' 3672 3395 IF ( TRIM( var ) == 'prr*' ) unit = 'mm/s' 3673 3396 IF ( TRIM( var ) == 'qsws*' ) unit = 'kgm/kgs' 3674 IF ( TRIM( var ) == 'qsws_eb*' ) unit = 'W/m2'3675 IF ( TRIM( var ) == 'qsws_liq_eb*' ) unit = 'W/m2'3676 IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'3677 IF ( TRIM( var ) == 'qsws_veg_eb*' ) unit = 'W/m2'3678 3397 IF ( TRIM( var ) == 'rad_net*' ) unit = 'W/m2' 3679 3398 IF ( TRIM( var ) == 'rrtm_aldif*' ) unit = '' … … 3681 3400 IF ( TRIM( var ) == 'rrtm_asdif*' ) unit = '' 3682 3401 IF ( TRIM( var ) == 'rrtm_asdir*' ) unit = '' 3683 IF ( TRIM( var ) == 'r_a*') unit = 's/m'3684 IF ( TRIM( var ) == 'r_s*') unit = 's/m'3685 3402 IF ( TRIM( var ) == 'shf*' ) unit = 'K*m/s' 3686 IF ( TRIM( var ) == 'shf_eb*') unit = 'W/m2'3687 3403 IF ( TRIM( var ) == 't*' ) unit = 'K' 3688 3404 IF ( TRIM( var ) == 'u*' ) unit = 'm/s' 3689 3405 IF ( TRIM( var ) == 'z0*' ) unit = 'm' 3690 IF ( TRIM( var ) == 'z0h*' ) unit = 'm' 3691 3692 3693 CASE ( 'p', 'pt', 'u', 'v', 'w' ) 3694 IF ( TRIM( var ) == 'p' ) unit = 'Pa' 3695 IF ( TRIM( var ) == 'pt' ) unit = 'K' 3696 IF ( TRIM( var ) == 'u' ) unit = 'm/s' 3697 IF ( TRIM( var ) == 'v' ) unit = 'm/s' 3698 IF ( TRIM( var ) == 'w' ) unit = 'm/s' 3699 CONTINUE 3700 3406 IF ( TRIM( var ) == 'z0h*' ) unit = 'm' 3407 3701 3408 CASE DEFAULT 3702 3409 3703 CALL user_check_data_output( var, unit ) 3410 CALL lsm_check_data_output ( var, unit, i, ilen, k ) 3411 3412 IF ( unit == 'illegal' ) THEN 3413 CALL user_check_data_output( var, unit ) 3414 ENDIF 3704 3415 3705 3416 IF ( unit == 'illegal' ) THEN -
palm/trunk/SOURCE/header.f90
r1809 r1817 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Moved land_surface_model header output to the respective module. 22 22 ! 23 23 ! Former revisions: … … 258 258 259 259 USE kinds 260 260 261 261 USE land_surface_model_mod, & 262 ONLY: conserve_water_content, land_surface, nzb_soil, & 263 nzt_soil, root_fraction, soil_moisture, soil_temperature, & 264 soil_type, soil_type_name, veg_type, veg_type_name, zs 262 ONLY: land_surface, lsm_header 265 263 266 264 USE model_1d, & … … 969 967 ENDIF 970 968 971 972 IF ( land_surface ) THEN 973 974 temperatures = '' 975 gradients = '' ! use for humidity here 976 coordinates = '' ! use for height 977 roots = '' ! use for root fraction 978 slices = '' ! use for index 979 980 i = 1 981 DO i = nzb_soil, nzt_soil 982 WRITE (coor_chr,'(F10.2,7X)') soil_temperature(i) 983 temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr ) 984 985 WRITE (coor_chr,'(F10.2,7X)') soil_moisture(i) 986 gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr ) 987 988 WRITE (coor_chr,'(F10.2,7X)') - zs(i) 989 coordinates = TRIM( coordinates ) // ' ' // TRIM( coor_chr ) 990 991 WRITE (coor_chr,'(F10.2,7X)') root_fraction(i) 992 roots = TRIM( roots ) // ' ' // TRIM( coor_chr ) 993 994 WRITE (coor_chr,'(I10,7X)') i 995 slices = TRIM( slices ) // ' ' // TRIM( coor_chr ) 996 997 998 ENDDO 999 1000 ! 1001 !-- Write land surface model header 1002 WRITE( io, 419 ) 1003 IF ( conserve_water_content ) THEN 1004 WRITE( io, 440 ) 1005 ELSE 1006 WRITE( io, 441 ) 1007 ENDIF 1008 1009 WRITE( io, 438 ) TRIM( veg_type_name(veg_type) ), & 1010 TRIM (soil_type_name(soil_type) ) 1011 WRITE( io, 439 ) TRIM( coordinates ), TRIM( temperatures ), & 1012 TRIM( gradients ), TRIM( roots ), TRIM( slices ) 1013 1014 1015 ENDIF 969 IF ( land_surface ) CALL lsm_header ( io ) 1016 970 1017 971 IF ( radiation ) THEN … … 2308 2262 418 FORMAT (/' Day of the year at model start : day_init = ',I3 & 2309 2263 /' UTC time at model start : time_utc_init = ',F7.1' s') 2310 419 FORMAT (//' Land surface model information:'/ &2311 ' ------------------------------'/)2312 2264 420 FORMAT (/' Characteristic levels of the initial temperature profile:'// & 2313 2265 ' Height: ',A,' m'/ & … … 2368 2320 '[0,1000] cm**2/s**3') 2369 2321 437 FORMAT (' Droplet collision is switched off') 2370 438 FORMAT (' --> Land surface type : ',A,/ &2371 ' --> Soil porosity type : ',A)2372 439 FORMAT (/' Initial soil temperature and moisture profile:'// &2373 ' Height: ',A,' m'/ &2374 ' Temperature: ',A,' K'/ &2375 ' Moisture: ',A,' m**3/m**3'/ &2376 ' Root fraction: ',A,' '/ &2377 ' Gridpoint: ',A)2378 440 FORMAT (' --> Soil bottom is closed (water content is conserved, default)')2379 441 FORMAT (' --> Soil bottom is open (water content is not conserved)')2380 2322 444 FORMAT (//' Radiation model information:'/ & 2381 2323 ' ----------------------------'/) -
palm/trunk/SOURCE/init_3d_model.f90
r1816 r1817 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Renamed lsm calls. 22 22 ! 23 23 ! Former revisions: … … 69 69 ! 70 70 ! 1551 2015-03-03 14:18:16Z maronga 71 ! Allocation of land surface arrays is now done in the subroutine init_lsm_arrays,71 ! Allocation of land surface arrays is now done in the subroutine lsm_init_arrays, 72 72 ! which is part of land_surface_model. 73 73 ! … … 266 266 267 267 USE land_surface_model_mod, & 268 ONLY: init_lsm, init_lsm_arrays, land_surface268 ONLY: lsm_init, lsm_init_arrays, land_surface 269 269 270 270 USE ls_forcing_mod … … 678 678 !-- Allocate land surface model arrays 679 679 IF ( land_surface ) THEN 680 CALL init_lsm_arrays680 CALL lsm_init_arrays 681 681 ENDIF 682 682 … … 1665 1665 IF ( land_surface ) THEN 1666 1666 CALL location_message( 'initializing land surface model', .FALSE. ) 1667 CALL init_lsm1667 CALL lsm_init 1668 1668 CALL location_message( 'finished', .TRUE. ) 1669 1669 ENDIF -
palm/trunk/SOURCE/land_surface_model_mod.f90
r1812 r1817 1 !> @file land_surface_model .f901 !> @file land_surface_model_mod.f90 2 2 !--------------------------------------------------------------------------------! 3 3 ! This file is part of PALM. … … 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Added interface for lsm_init_arrays. Added subroutines for check_parameters, 22 ! header, and parin. Renamed some subroutines. 22 23 ! 23 24 ! Former revisions: … … 139 140 ONLY: cloud_physics, dt_3d, humidity, intermediate_timestep_count, & 140 141 initializing_actions, intermediate_timestep_count_max, & 141 max_masks, precipitation, pt_surface, rho_surface,&142 r oughness_length, surface_pressure, timestep_scheme, tsc,&143 z0h_factor, time_since_reference_point142 max_masks, precipitation, pt_surface, & 143 rho_surface, roughness_length, surface_pressure, & 144 timestep_scheme, tsc, z0h_factor, time_since_reference_point 144 145 145 146 USE indices, & … … 523 524 PRIVATE 524 525 525 526 527 ! 528 !-- Public functions 529 PUBLIC lsm_check_data_output, lsm_check_data_output_pr, & 530 lsm_check_parameters, lsm_energy_balance, lsm_header, lsm_init, & 531 lsm_init_arrays, lsm_parin, lsm_soil_model 526 532 ! 527 533 !-- Public parameters, constants and initial values 528 534 PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient, & 529 535 conserve_water_content, field_capacity, & 530 f_shortwave_incoming, hydraulic_conductivity, init_lsm,&531 init_lsm_arrays, lambda_surface_stable, lambda_surface_unstable,&532 land_surface, leaf_area_index, lsm_energy_balance, lsm_soil_model,&533 lsm_swap_timelevel, l_vangenuchten, min_canopy_resistance,&534 min_ soil_resistance, n_vangenuchten, pave_heat_capacity,&535 pave_ depth, pave_heat_conductivity, residual_moisture, rho_cp,&536 r ho_lv, root_fraction, saturation_moisture, skip_time_do_lsm,&537 s oil_moisture, soil_temperature, soil_type, soil_type_name,&538 vegetation_coverage, veg_type, veg_type_name, wilting_point, z0_eb,&539 z0h_eb, z0q_eb536 f_shortwave_incoming, hydraulic_conductivity, & 537 lambda_surface_stable, lambda_surface_unstable, & 538 land_surface, leaf_area_index, & 539 lsm_swap_timelevel, l_vangenuchten, & 540 min_canopy_resistance, min_soil_resistance, n_vangenuchten, & 541 pave_heat_capacity, pave_depth, pave_heat_conductivity, & 542 residual_moisture, rho_cp, rho_lv, root_fraction, & 543 saturation_moisture, skip_time_do_lsm, soil_moisture, & 544 soil_temperature, soil_type, soil_type_name, vegetation_coverage, & 545 veg_type, veg_type_name, wilting_point, z0_eb, z0h_eb, z0q_eb 540 546 541 547 ! … … 554 560 PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av 555 561 556 INTERFACE init_lsm 557 MODULE PROCEDURE init_lsm 558 END INTERFACE init_lsm 559 562 563 INTERFACE lsm_check_data_output 564 MODULE PROCEDURE lsm_check_data_output 565 END INTERFACE lsm_check_data_output 566 567 INTERFACE lsm_check_data_output_pr 568 MODULE PROCEDURE lsm_check_data_output_pr 569 END INTERFACE lsm_check_data_output_pr 570 571 INTERFACE lsm_check_parameters 572 MODULE PROCEDURE lsm_check_parameters 573 END INTERFACE lsm_check_parameters 574 560 575 INTERFACE lsm_energy_balance 561 576 MODULE PROCEDURE lsm_energy_balance 562 577 END INTERFACE lsm_energy_balance 563 578 579 INTERFACE lsm_header 580 MODULE PROCEDURE lsm_header 581 END INTERFACE lsm_header 582 583 INTERFACE lsm_init 584 MODULE PROCEDURE lsm_init 585 END INTERFACE lsm_init 586 587 INTERFACE lsm_init_arrays 588 MODULE PROCEDURE lsm_init_arrays 589 END INTERFACE lsm_init_arrays 590 591 INTERFACE lsm_parin 592 MODULE PROCEDURE lsm_parin 593 END INTERFACE lsm_parin 594 564 595 INTERFACE lsm_soil_model 565 596 MODULE PROCEDURE lsm_soil_model … … 571 602 572 603 CONTAINS 604 605 !------------------------------------------------------------------------------! 606 ! Description: 607 ! ------------ 608 !> Check data output for land surface model 609 !------------------------------------------------------------------------------! 610 SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k ) 611 612 613 USE control_parameters, & 614 ONLY: data_output, message_string 615 616 IMPLICIT NONE 617 618 CHARACTER (LEN=*) :: unit !< 619 CHARACTER (LEN=*) :: var !< 620 621 INTEGER(iwp) :: i 622 INTEGER(iwp) :: ilen 623 INTEGER(iwp) :: k 624 625 SELECT CASE ( TRIM( var ) ) 626 627 CASE ( 'm_soil' ) 628 IF ( .NOT. land_surface ) THEN 629 message_string = 'output of "' // TRIM( var ) // '" requi' // & 630 'res land_surface = .TRUE.' 631 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 632 ENDIF 633 unit = 'm3/m3' 634 635 CASE ( 't_soil' ) 636 IF ( .NOT. land_surface ) THEN 637 message_string = 'output of "' // TRIM( var ) // '" requi' // & 638 'res land_surface = .TRUE.' 639 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 640 ENDIF 641 unit = 'K' 642 643 CASE ( 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', & 644 'm_liq_eb*', 'qsws_eb*', & 645 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*', & 646 'r_a*', 'r_s*', 'shf_eb*' ) 647 IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN 648 message_string = 'illegal value for data_output: "' // & 649 TRIM( var ) // '" & only 2d-horizontal ' // & 650 'cross sections are allowed for this value' 651 CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 ) 652 ENDIF 653 654 IF ( TRIM( var ) == 'c_liq*' .AND. .NOT. land_surface ) THEN 655 message_string = 'output of "' // TRIM( var ) // '" requi' // & 656 'res land_surface = .TRUE.' 657 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 658 ENDIF 659 IF ( TRIM( var ) == 'c_soil*' .AND. .NOT. land_surface ) THEN 660 message_string = 'output of "' // TRIM( var ) // '" requi' // & 661 'res land_surface = .TRUE.' 662 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 663 ENDIF 664 IF ( TRIM( var ) == 'c_veg*' .AND. .NOT. land_surface ) THEN 665 message_string = 'output of "' // TRIM( var ) // '" requi' // & 666 'res land_surface = .TRUE.' 667 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 668 ENDIF 669 IF ( TRIM( var ) == 'ghf_eb*' .AND. .NOT. land_surface ) THEN 670 message_string = 'output of "' // TRIM( var ) // '" requi' // & 671 'res land_surface = .TRUE.' 672 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 673 ENDIF 674 IF ( TRIM( var ) == 'm_liq_eb*' .AND. .NOT. land_surface ) THEN 675 message_string = 'output of "' // TRIM( var ) // '" requi' // & 676 'res land_surface = .TRUE.' 677 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 678 ENDIF 679 IF ( TRIM( var ) == 'qsws_eb*' .AND. .NOT. land_surface ) THEN 680 message_string = 'output of "' // TRIM( var ) // '" requi' // & 681 'res land_surface = .TRUE.' 682 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 683 ENDIF 684 IF ( TRIM( var ) == 'qsws_liq_eb*' .AND. .NOT. land_surface ) & 685 THEN 686 message_string = 'output of "' // TRIM( var ) // '" requi' // & 687 'res land_surface = .TRUE.' 688 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 689 ENDIF 690 IF ( TRIM( var ) == 'qsws_soil_eb*' .AND. .NOT. land_surface ) & 691 THEN 692 message_string = 'output of "' // TRIM( var ) // '" requi' // & 693 'res land_surface = .TRUE.' 694 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 695 ENDIF 696 IF ( TRIM( var ) == 'qsws_veg_eb*' .AND. .NOT. land_surface ) & 697 THEN 698 message_string = 'output of "' // TRIM( var ) // '" requi' // & 699 'res land_surface = .TRUE.' 700 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 701 ENDIF 702 IF ( TRIM( var ) == 'r_a*' .AND. .NOT. land_surface ) & 703 THEN 704 message_string = 'output of "' // TRIM( var ) // '" requi' // & 705 'res land_surface = .TRUE.' 706 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 707 ENDIF 708 IF ( TRIM( var ) == 'r_s*' .AND. .NOT. land_surface ) & 709 THEN 710 message_string = 'output of "' // TRIM( var ) // '" requi' // & 711 'res land_surface = .TRUE.' 712 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 ) 713 ENDIF 714 715 IF ( TRIM( var ) == 'c_liq*' ) unit = 'none' 716 IF ( TRIM( var ) == 'c_soil*') unit = 'none' 717 IF ( TRIM( var ) == 'c_veg*' ) unit = 'none' 718 IF ( TRIM( var ) == 'ghf_eb*') unit = 'W/m2' 719 IF ( TRIM( var ) == 'm_liq_eb*' ) unit = 'm' 720 IF ( TRIM( var ) == 'qsws_eb*' ) unit = 'W/m2' 721 IF ( TRIM( var ) == 'qsws_liq_eb*' ) unit = 'W/m2' 722 IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2' 723 IF ( TRIM( var ) == 'qsws_veg_eb*' ) unit = 'W/m2' 724 IF ( TRIM( var ) == 'r_a*') unit = 's/m' 725 IF ( TRIM( var ) == 'r_s*') unit = 's/m' 726 IF ( TRIM( var ) == 'shf_eb*') unit = 'W/m2' 727 728 CASE DEFAULT 729 unit = 'illegal' 730 731 END SELECT 732 733 734 END SUBROUTINE lsm_check_data_output 735 736 SUBROUTINE lsm_check_data_output_pr( variable, var_count, unit, dopr_unit ) 737 738 USE control_parameters, & 739 ONLY: data_output_pr, message_string 740 741 USE indices 742 743 USE profil_parameter 744 745 USE statistics 746 747 IMPLICIT NONE 748 749 CHARACTER (LEN=*) :: unit !< 750 CHARACTER (LEN=*) :: variable !< 751 CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit 752 753 INTEGER(iwp) :: user_pr_index !< 754 INTEGER(iwp) :: var_count !< 755 756 SELECT CASE ( TRIM( variable ) ) 757 758 CASE ( 't_soil', '#t_soil' ) 759 IF ( .NOT. land_surface ) THEN 760 message_string = 'data_output_pr = ' // & 761 TRIM( data_output_pr(var_count) ) // ' is' // & 762 'not implemented for land_surface = .FALSE.' 763 CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 ) 764 ELSE 765 dopr_index(var_count) = 89 766 dopr_unit = 'K' 767 hom(0:nzs-1,2,89,:) = SPREAD( - zs, 2, statistic_regions+1 ) 768 IF ( data_output_pr(var_count)(1:1) == '#' ) THEN 769 dopr_initial_index(var_count) = 90 770 hom(0:nzs-1,2,90,:) = SPREAD( - zs, 2, statistic_regions+1 ) 771 data_output_pr(var_count) = data_output_pr(var_count)(2:) 772 ENDIF 773 unit = dopr_unit 774 ENDIF 775 776 CASE ( 'm_soil', '#m_soil' ) 777 IF ( .NOT. land_surface ) THEN 778 message_string = 'data_output_pr = ' // & 779 TRIM( data_output_pr(var_count) ) // ' is' // & 780 ' not implemented for land_surface = .FALSE.' 781 CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 ) 782 ELSE 783 dopr_index(var_count) = 91 784 dopr_unit = 'm3/m3' 785 hom(0:nzs-1,2,91,:) = SPREAD( - zs, 2, statistic_regions+1 ) 786 IF ( data_output_pr(var_count)(1:1) == '#' ) THEN 787 dopr_initial_index(var_count) = 92 788 hom(0:nzs-1,2,92,:) = SPREAD( - zs, 2, statistic_regions+1 ) 789 data_output_pr(var_count) = data_output_pr(var_count)(2:) 790 ENDIF 791 unit = dopr_unit 792 ENDIF 793 794 795 CASE DEFAULT 796 unit = 'illegal' 797 798 END SELECT 799 800 801 END SUBROUTINE lsm_check_data_output_pr 802 803 804 !------------------------------------------------------------------------------! 805 ! Description: 806 ! ------------ 807 !> Check parameters routine for land surface model 808 !------------------------------------------------------------------------------! 809 SUBROUTINE lsm_check_parameters 810 811 USE control_parameters, & 812 ONLY: bc_pt_b, bc_q_b, constant_flux_layer, message_string, & 813 most_method, topography 814 815 USE radiation_model_mod, & 816 ONLY: radiation 817 818 819 IMPLICIT NONE 820 821 822 ! 823 !-- Dirichlet boundary conditions are required as the surface fluxes are 824 !-- calculated from the temperature/humidity gradients in the land surface 825 !-- model 826 IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' ) THEN 827 message_string = 'lsm requires setting of'// & 828 'bc_pt_b = "dirichlet" and '// & 829 'bc_q_b = "dirichlet"' 830 CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 ) 831 ENDIF 832 833 IF ( .NOT. constant_flux_layer ) THEN 834 message_string = 'lsm requires '// & 835 'constant_flux_layer = .T.' 836 CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 ) 837 ENDIF 838 839 IF ( topography /= 'flat' ) THEN 840 message_string = 'lsm cannot be used ' // & 841 'in combination with topography /= "flat"' 842 CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 ) 843 ENDIF 844 845 IF ( ( veg_type == 14 .OR. veg_type == 15 ) .AND. & 846 most_method == 'lookup' ) THEN 847 WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ', & 848 'allowed in combination with ', & 849 'most_method = ', most_method 850 CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 ) 851 ENDIF 852 853 IF ( veg_type == 0 ) THEN 854 IF ( SUM( root_fraction ) /= 1.0_wp ) THEN 855 message_string = 'veg_type = 0 (user_defined)'// & 856 'requires setting of root_fraction(0:3)'// & 857 '/= 9999999.9 and SUM(root_fraction) = 1' 858 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 859 ENDIF 860 861 IF ( min_canopy_resistance == 9999999.9_wp ) THEN 862 message_string = 'veg_type = 0 (user defined)'// & 863 'requires setting of min_canopy_resistance'// & 864 '/= 9999999.9' 865 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 866 ENDIF 867 868 IF ( leaf_area_index == 9999999.9_wp ) THEN 869 message_string = 'veg_type = 0 (user_defined)'// & 870 'requires setting of leaf_area_index'// & 871 '/= 9999999.9' 872 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 873 ENDIF 874 875 IF ( vegetation_coverage == 9999999.9_wp ) THEN 876 message_string = 'veg_type = 0 (user_defined)'// & 877 'requires setting of vegetation_coverage'// & 878 '/= 9999999.9' 879 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 880 ENDIF 881 882 IF ( canopy_resistance_coefficient == 9999999.9_wp) THEN 883 message_string = 'veg_type = 0 (user_defined)'// & 884 'requires setting of'// & 885 'canopy_resistance_coefficient /= 9999999.9' 886 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 887 ENDIF 888 889 IF ( lambda_surface_stable == 9999999.9_wp ) THEN 890 message_string = 'veg_type = 0 (user_defined)'// & 891 'requires setting of lambda_surface_stable'// & 892 '/= 9999999.9' 893 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 894 ENDIF 895 896 IF ( lambda_surface_unstable == 9999999.9_wp ) THEN 897 message_string = 'veg_type = 0 (user_defined)'// & 898 'requires setting of lambda_surface_unstable'// & 899 '/= 9999999.9' 900 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 901 ENDIF 902 903 IF ( f_shortwave_incoming == 9999999.9_wp ) THEN 904 message_string = 'veg_type = 0 (user_defined)'// & 905 'requires setting of f_shortwave_incoming'// & 906 '/= 9999999.9' 907 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 908 ENDIF 909 910 IF ( z0_eb == 9999999.9_wp ) THEN 911 message_string = 'veg_type = 0 (user_defined)'// & 912 'requires setting of z0_eb'// & 913 '/= 9999999.9' 914 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 915 ENDIF 916 917 IF ( z0h_eb == 9999999.9_wp ) THEN 918 message_string = 'veg_type = 0 (user_defined)'// & 919 'requires setting of z0h_eb'// & 920 '/= 9999999.9' 921 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 ) 922 ENDIF 923 924 925 ENDIF 926 927 IF ( soil_type == 0 ) THEN 928 929 IF ( alpha_vangenuchten == 9999999.9_wp ) THEN 930 message_string = 'soil_type = 0 (user_defined)'// & 931 'requires setting of alpha_vangenuchten'// & 932 '/= 9999999.9' 933 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 934 ENDIF 935 936 IF ( l_vangenuchten == 9999999.9_wp ) THEN 937 message_string = 'soil_type = 0 (user_defined)'// & 938 'requires setting of l_vangenuchten'// & 939 '/= 9999999.9' 940 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 941 ENDIF 942 943 IF ( n_vangenuchten == 9999999.9_wp ) THEN 944 message_string = 'soil_type = 0 (user_defined)'// & 945 'requires setting of n_vangenuchten'// & 946 '/= 9999999.9' 947 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 948 ENDIF 949 950 IF ( hydraulic_conductivity == 9999999.9_wp ) THEN 951 message_string = 'soil_type = 0 (user_defined)'// & 952 'requires setting of hydraulic_conductivity'// & 953 '/= 9999999.9' 954 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 955 ENDIF 956 957 IF ( saturation_moisture == 9999999.9_wp ) THEN 958 message_string = 'soil_type = 0 (user_defined)'// & 959 'requires setting of saturation_moisture'// & 960 '/= 9999999.9' 961 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 962 ENDIF 963 964 IF ( field_capacity == 9999999.9_wp ) THEN 965 message_string = 'soil_type = 0 (user_defined)'// & 966 'requires setting of field_capacity'// & 967 '/= 9999999.9' 968 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 969 ENDIF 970 971 IF ( wilting_point == 9999999.9_wp ) THEN 972 message_string = 'soil_type = 0 (user_defined)'// & 973 'requires setting of wilting_point'// & 974 '/= 9999999.9' 975 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 976 ENDIF 977 978 IF ( residual_moisture == 9999999.9_wp ) THEN 979 message_string = 'soil_type = 0 (user_defined)'// & 980 'requires setting of residual_moisture'// & 981 '/= 9999999.9' 982 CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 ) 983 ENDIF 984 985 ENDIF 986 987 IF ( .NOT. radiation ) THEN 988 message_string = 'lsm requires '// & 989 'radiation = .T.' 990 CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 ) 991 ENDIF 992 993 994 END SUBROUTINE lsm_check_parameters 995 996 !------------------------------------------------------------------------------! 997 ! Description: 998 ! ------------ 999 !> Solver for the energy balance at the surface. 1000 !------------------------------------------------------------------------------! 1001 SUBROUTINE lsm_energy_balance 1002 1003 1004 IMPLICIT NONE 1005 1006 INTEGER(iwp) :: i !< running index 1007 INTEGER(iwp) :: j !< running index 1008 INTEGER(iwp) :: k, ks !< running index 1009 1010 REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface 1011 f1, & !< resistance correction term 1 1012 f2, & !< resistance correction term 2 1013 f3, & !< resistance correction term 3 1014 m_min, & !< minimum soil moisture 1015 e, & !< water vapour pressure 1016 e_s, & !< water vapour saturation pressure 1017 e_s_dt, & !< derivate of e_s with respect to T 1018 tend, & !< tendency 1019 dq_s_dt, & !< derivate of q_s with respect to T 1020 coef_1, & !< coef. for prognostic equation 1021 coef_2, & !< coef. for prognostic equation 1022 f_qsws, & !< factor for qsws_eb 1023 f_qsws_veg, & !< factor for qsws_veg_eb 1024 f_qsws_soil, & !< factor for qsws_soil_eb 1025 f_qsws_liq, & !< factor for qsws_liq_eb 1026 f_shf, & !< factor for shf_eb 1027 lambda_surface, & !< Current value of lambda_surface 1028 m_liq_eb_max, & !< maxmimum value of the liq. water reservoir 1029 pt1, & !< potential temperature at first grid level 1030 qv1 !< specific humidity at first grid level 1031 1032 ! 1033 !-- Calculate the exner function for the current time step 1034 exn = ( surface_pressure / 1000.0_wp )**0.286_wp 1035 1036 DO i = nxlg, nxrg 1037 DO j = nysg, nyng 1038 k = nzb_s_inner(j,i) 1039 1040 ! 1041 !-- Set lambda_surface according to stratification between skin layer and soil 1042 IF ( .NOT. pave_surface(j,i) ) THEN 1043 1044 c_surface_tmp = c_surface 1045 1046 IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i)) THEN 1047 lambda_surface = lambda_surface_s(j,i) 1048 ELSE 1049 lambda_surface = lambda_surface_u(j,i) 1050 ENDIF 1051 ELSE 1052 1053 c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp 1054 lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil) 1055 1056 ENDIF 1057 1058 ! 1059 !-- First step: calculate aerodyamic resistance. As pt, us, ts 1060 !-- are not available for the prognostic time step, data from the last 1061 !-- time step is used here. Note that this formulation is the 1062 !-- equivalent to the ECMWF formulation using drag coefficients 1063 IF ( cloud_physics ) THEN 1064 pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i) 1065 qv1 = q(k+1,j,i) - ql(k+1,j,i) 1066 ELSE 1067 pt1 = pt(k+1,j,i) 1068 qv1 = q(k+1,j,i) 1069 ENDIF 1070 1071 r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp) 1072 1073 ! 1074 !-- Make sure that the resistance does not drop to zero 1075 IF ( ABS(r_a(j,i)) < 1.0E-10_wp ) r_a(j,i) = 1.0E-10_wp 1076 1077 ! 1078 !-- Second step: calculate canopy resistance r_canopy 1079 !-- f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation 1080 1081 !-- f1: correction for incoming shortwave radiation (stomata close at 1082 !-- night) 1083 IF ( radiation_scheme /= 'constant' ) THEN 1084 f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) / & 1085 (0.81_wp * (0.004_wp * rad_sw_in(k,j,i) & 1086 + 1.0_wp)) ) 1087 ELSE 1088 f1 = 1.0_wp 1089 ENDIF 1090 1091 1092 ! 1093 !-- f2: correction for soil moisture availability to plants (the 1094 !-- integrated soil moisture must thus be considered here) 1095 !-- f2 = 0 for very dry soils 1096 m_total = 0.0_wp 1097 DO ks = nzb_soil, nzt_soil 1098 m_total = m_total + root_fr(ks,j,i) & 1099 * MAX(m_soil(ks,j,i),m_wilt(j,i)) 1100 ENDDO 1101 1102 IF ( m_total > m_wilt(j,i) .AND. m_total < m_fc(j,i) ) THEN 1103 f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) ) 1104 ELSEIF ( m_total >= m_fc(j,i) ) THEN 1105 f2 = 1.0_wp 1106 ELSE 1107 f2 = 1.0E-20_wp 1108 ENDIF 1109 1110 ! 1111 !-- Calculate water vapour pressure at saturation 1112 e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i) & 1113 - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) ) 1114 1115 ! 1116 !-- f3: correction for vapour pressure deficit 1117 IF ( g_d(j,i) /= 0.0_wp ) THEN 1118 ! 1119 !-- Calculate vapour pressure 1120 e = qv1 * surface_pressure / 0.622_wp 1121 f3 = EXP ( -g_d(j,i) * (e_s - e) ) 1122 ELSE 1123 f3 = 1.0_wp 1124 ENDIF 1125 1126 ! 1127 !-- Calculate canopy resistance. In case that c_veg is 0 (bare soils), 1128 !-- this calculation is obsolete, as r_canopy is not used below. 1129 !-- To do: check for very dry soil -> r_canopy goes to infinity 1130 r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3 & 1131 + 1.0E-20_wp) 1132 1133 ! 1134 !-- Third step: calculate bare soil resistance r_soil. The Clapp & 1135 !-- Hornberger parametrization does not consider c_veg. 1136 IF ( soil_type_2d(j,i) /= 7 ) THEN 1137 m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) * & 1138 m_res(j,i) 1139 ELSE 1140 m_min = m_wilt(j,i) 1141 ENDIF 1142 1143 f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min ) 1144 f2 = MAX(f2,1.0E-20_wp) 1145 f2 = MIN(f2,1.0_wp) 1146 1147 r_soil(j,i) = r_soil_min(j,i) / f2 1148 1149 ! 1150 !-- Calculate the maximum possible liquid water amount on plants and 1151 !-- bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is 1152 !-- assumed, while paved surfaces might hold up 1 mm of water. The 1153 !-- liquid water fraction for paved surfaces is calculated after 1154 !-- Noilhan & Planton (1989), while the ECMWF formulation is used for 1155 !-- vegetated surfaces and bare soils. 1156 IF ( pave_surface(j,i) ) THEN 1157 m_liq_eb_max = m_max_depth * 5.0_wp 1158 c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 ) 1159 ELSE 1160 m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i) & 1161 + (1.0_wp - c_veg(j,i)) ) 1162 c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max ) 1163 ENDIF 1164 1165 ! 1166 !-- Calculate saturation specific humidity 1167 q_s = 0.622_wp * e_s / surface_pressure 1168 1169 ! 1170 !-- In case of dewfall, set evapotranspiration to zero 1171 !-- All super-saturated water is then removed from the air 1172 IF ( humidity .AND. q_s <= qv1 ) THEN 1173 r_canopy(j,i) = 0.0_wp 1174 r_soil(j,i) = 0.0_wp 1175 ENDIF 1176 1177 ! 1178 !-- Calculate coefficients for the total evapotranspiration 1179 !-- In case of water surface, set vegetation and soil fluxes to zero. 1180 !-- For pavements, only evaporation of liquid water is possible. 1181 IF ( water_surface(j,i) ) THEN 1182 f_qsws_veg = 0.0_wp 1183 f_qsws_soil = 0.0_wp 1184 f_qsws_liq = rho_lv / r_a(j,i) 1185 ELSEIF ( pave_surface (j,i) ) THEN 1186 f_qsws_veg = 0.0_wp 1187 f_qsws_soil = 0.0_wp 1188 f_qsws_liq = rho_lv * c_liq(j,i) / r_a(j,i) 1189 ELSE 1190 f_qsws_veg = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/ & 1191 (r_a(j,i) + r_canopy(j,i)) 1192 f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) + & 1193 r_soil(j,i)) 1194 f_qsws_liq = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i) 1195 ENDIF 1196 ! 1197 !-- If soil moisture is below wilting point, plants do no longer 1198 !-- transpirate. 1199 ! IF ( m_soil(k,j,i) < m_wilt(j,i) ) THEN 1200 ! f_qsws_veg = 0.0_wp 1201 ! ENDIF 1202 1203 f_shf = rho_cp / r_a(j,i) 1204 f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq 1205 1206 ! 1207 !-- Calculate derivative of q_s for Taylor series expansion 1208 e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) - & 1209 17.269_wp*(t_surface(j,i) - 273.16_wp) & 1210 / (t_surface(j,i) - 35.86_wp)**2 ) 1211 1212 dq_s_dt = 0.622_wp * e_s_dt / surface_pressure 1213 1214 ! 1215 !-- Add LW up so that it can be removed in prognostic equation 1216 rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i) 1217 1218 ! 1219 !-- Calculate new skin temperature 1220 IF ( humidity ) THEN 1221 #if defined ( __rrtmg ) 1222 ! 1223 !-- Numerator of the prognostic equation 1224 coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1225 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1226 + f_shf * pt1 + f_qsws * ( qv1 - q_s & 1227 + dq_s_dt * t_surface(j,i) ) + lambda_surface & 1228 * t_soil(nzb_soil,j,i) 1229 1230 ! 1231 !-- Denominator of the prognostic equation 1232 coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt & 1233 + lambda_surface + f_shf / exn 1234 #else 1235 1236 ! 1237 !-- Numerator of the prognostic equation 1238 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1239 * t_surface(j,i) ** 4 & 1240 + f_shf * pt1 + f_qsws * ( qv1 & 1241 - q_s + dq_s_dt * t_surface(j,i) ) & 1242 + lambda_surface * t_soil(nzb_soil,j,i) 1243 1244 ! 1245 !-- Denominator of the prognostic equation 1246 coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws & 1247 * dq_s_dt + lambda_surface + f_shf / exn 1248 1249 #endif 1250 ELSE 1251 1252 #if defined ( __rrtmg ) 1253 ! 1254 !-- Numerator of the prognostic equation 1255 coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1256 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1257 + f_shf * pt1 + lambda_surface & 1258 * t_soil(nzb_soil,j,i) 1259 1260 ! 1261 !-- Denominator of the prognostic equation 1262 coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn 1263 #else 1264 1265 ! 1266 !-- Numerator of the prognostic equation 1267 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1268 * t_surface(j,i) ** 4 + f_shf * pt1 & 1269 + lambda_surface * t_soil(nzb_soil,j,i) 1270 1271 ! 1272 !-- Denominator of the prognostic equation 1273 coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 & 1274 + lambda_surface + f_shf / exn 1275 #endif 1276 ENDIF 1277 1278 tend = 0.0_wp 1279 1280 ! 1281 !-- Implicit solution when the surface layer has no heat capacity, 1282 !-- otherwise use RK3 scheme. 1283 t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp * & 1284 t_surface(j,i) ) / ( c_surface_tmp + coef_2 & 1285 * dt_3d * tsc(2) ) 1286 1287 ! 1288 !-- Add RK3 term 1289 IF ( c_surface_tmp /= 0.0_wp ) THEN 1290 1291 t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3) & 1292 * tt_surface_m(j,i) 1293 1294 ! 1295 !-- Calculate true tendency 1296 tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3) & 1297 * tt_surface_m(j,i)) / (dt_3d * tsc(2)) 1298 ! 1299 !-- Calculate t_surface tendencies for the next Runge-Kutta step 1300 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1301 IF ( intermediate_timestep_count == 1 ) THEN 1302 tt_surface_m(j,i) = tend 1303 ELSEIF ( intermediate_timestep_count < & 1304 intermediate_timestep_count_max ) THEN 1305 tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp & 1306 * tt_surface_m(j,i) 1307 ENDIF 1308 ENDIF 1309 ENDIF 1310 1311 ! 1312 !-- In case of fast changes in the skin temperature, it is possible to 1313 !-- update the radiative fluxes independently from the prescribed 1314 !-- radiation call frequency. This effectively prevents oscillations, 1315 !-- especially when setting skip_time_do_radiation /= 0. The threshold 1316 !-- value of 0.2 used here is just a first guess. This method should be 1317 !-- revised in the future as tests have shown that the threshold is 1318 !-- often reached, when no oscillations would occur (causes immense 1319 !-- computing time for the radiation code). 1320 IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp .AND. & 1321 unscheduled_radiation_calls ) THEN 1322 force_radiation_call_l = .TRUE. 1323 ENDIF 1324 1325 pt(k,j,i) = t_surface_p(j,i) / exn 1326 1327 ! 1328 !-- Calculate fluxes 1329 #if defined ( __rrtmg ) 1330 rad_net_l(j,i) = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1331 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1332 - rad_lw_out_change_0(j,i) * t_surface_p(j,i) 1333 1334 IF ( rrtm_idrv == 1 ) THEN 1335 rad_net(j,i) = rad_net_l(j,i) 1336 rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) & 1337 + rad_lw_out_change_0(j,i) & 1338 * ( t_surface_p(j,i) - t_surface(j,i) ) 1339 ENDIF 1340 #else 1341 rad_net_l(j,i) = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1342 * t_surface(j,i)**4 - 4.0_wp * sigma_sb & 1343 * t_surface(j,i)**3 * t_surface_p(j,i) 1344 #endif 1345 1346 ghf_eb(j,i) = lambda_surface * (t_surface_p(j,i) & 1347 - t_soil(nzb_soil,j,i)) 1348 1349 shf_eb(j,i) = - f_shf * ( pt1 - pt(k,j,i) ) 1350 1351 shf(j,i) = shf_eb(j,i) / rho_cp 1352 1353 IF ( humidity ) THEN 1354 qsws_eb(j,i) = - f_qsws * ( qv1 - q_s + dq_s_dt & 1355 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) 1356 1357 qsws(j,i) = qsws_eb(j,i) / rho_lv 1358 1359 qsws_veg_eb(j,i) = - f_qsws_veg * ( qv1 - q_s & 1360 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1361 * t_surface_p(j,i) ) 1362 1363 qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s & 1364 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1365 * t_surface_p(j,i) ) 1366 1367 qsws_liq_eb(j,i) = - f_qsws_liq * ( qv1 - q_s & 1368 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1369 * t_surface_p(j,i) ) 1370 ENDIF 1371 1372 ! 1373 !-- Calculate the true surface resistance 1374 IF ( qsws_eb(j,i) == 0.0_wp ) THEN 1375 r_s(j,i) = 1.0E10_wp 1376 ELSE 1377 r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt & 1378 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) & 1379 / qsws_eb(j,i) - r_a(j,i) 1380 ENDIF 1381 1382 ! 1383 !-- Calculate change in liquid water reservoir due to dew fall or 1384 !-- evaporation of liquid water 1385 IF ( humidity ) THEN 1386 ! 1387 !-- If precipitation is activated, add rain water to qsws_liq_eb 1388 !-- and qsws_soil_eb according the the vegetation coverage. 1389 !-- precipitation_rate is given in mm. 1390 IF ( precipitation ) THEN 1391 1392 ! 1393 !-- Add precipitation to liquid water reservoir, if possible. 1394 !-- Otherwise, add the water to soil. In case of 1395 !-- pavements, the exceeding water amount is implicitely removed 1396 !-- as runoff as qsws_soil_eb is then not used in the soil model 1397 IF ( m_liq_eb(j,i) /= m_liq_eb_max ) THEN 1398 qsws_liq_eb(j,i) = qsws_liq_eb(j,i) & 1399 + c_veg(j,i) * prr(k,j,i) * hyrho(k) & 1400 * 0.001_wp * rho_l * l_v 1401 ELSE 1402 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) & 1403 + c_veg(j,i) * prr(k,j,i) * hyrho(k) & 1404 * 0.001_wp * rho_l * l_v 1405 ENDIF 1406 1407 !-- Add precipitation to bare soil according to the bare soil 1408 !-- coverage. 1409 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) * (1.0_wp & 1410 - c_veg(j,i)) * prr(k,j,i) * hyrho(k) & 1411 * 0.001_wp * rho_l * l_v 1412 ENDIF 1413 1414 ! 1415 !-- If the air is saturated, check the reservoir water level 1416 IF ( qsws_eb(j,i) < 0.0_wp ) THEN 1417 1418 ! 1419 !-- Check if reservoir is full (avoid values > m_liq_eb_max) 1420 !-- In that case, qsws_liq_eb goes to qsws_soil_eb. In this 1421 !-- case qsws_veg_eb is zero anyway (because c_liq = 1), 1422 !-- so that tend is zero and no further check is needed 1423 IF ( m_liq_eb(j,i) == m_liq_eb_max ) THEN 1424 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) & 1425 + qsws_liq_eb(j,i) 1426 qsws_liq_eb(j,i) = 0.0_wp 1427 ENDIF 1428 1429 ! 1430 !-- In case qsws_veg_eb becomes negative (unphysical behavior), 1431 !-- let the water enter the liquid water reservoir as dew on the 1432 !-- plant 1433 IF ( qsws_veg_eb(j,i) < 0.0_wp ) THEN 1434 qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i) 1435 qsws_veg_eb(j,i) = 0.0_wp 1436 ENDIF 1437 ENDIF 1438 1439 tend = - qsws_liq_eb(j,i) * drho_l_lv 1440 1441 m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend & 1442 + tsc(3) * tm_liq_eb_m(j,i) ) 1443 1444 ! 1445 !-- Check if reservoir is overfull -> reduce to maximum 1446 !-- (conservation of water is violated here) 1447 m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max) 1448 1449 ! 1450 !-- Check if reservoir is empty (avoid values < 0.0) 1451 !-- (conservation of water is violated here) 1452 m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp) 1453 1454 1455 ! 1456 !-- Calculate m_liq_eb tendencies for the next Runge-Kutta step 1457 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1458 IF ( intermediate_timestep_count == 1 ) THEN 1459 tm_liq_eb_m(j,i) = tend 1460 ELSEIF ( intermediate_timestep_count < & 1461 intermediate_timestep_count_max ) THEN 1462 tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp & 1463 * tm_liq_eb_m(j,i) 1464 ENDIF 1465 ENDIF 1466 1467 ENDIF 1468 1469 ENDDO 1470 ENDDO 1471 1472 ! 1473 !-- Make a logical OR for all processes. Force radiation call if at 1474 !-- least one processor reached the threshold change in skin temperature 1475 IF ( unscheduled_radiation_calls .AND. intermediate_timestep_count & 1476 == intermediate_timestep_count_max-1 ) THEN 1477 #if defined( __parallel ) 1478 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1479 CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call, & 1480 1, MPI_LOGICAL, MPI_LOR, comm2d, ierr ) 1481 #else 1482 force_radiation_call = force_radiation_call_l 1483 #endif 1484 force_radiation_call_l = .FALSE. 1485 ENDIF 1486 1487 ! 1488 !-- Calculate surface specific humidity 1489 IF ( humidity ) THEN 1490 CALL calc_q_surface 1491 ENDIF 1492 1493 ! 1494 !-- Calculate new roughness lengths (for water surfaces only) 1495 CALL calc_z0_water_surface 1496 1497 1498 END SUBROUTINE lsm_energy_balance 1499 1500 !------------------------------------------------------------------------------! 1501 ! Description: 1502 ! ------------ 1503 !> Header output for land surface model 1504 !------------------------------------------------------------------------------! 1505 SUBROUTINE lsm_header ( io ) 1506 1507 1508 IMPLICIT NONE 1509 1510 CHARACTER (LEN=86) :: t_soil_chr !< String for soil temperature profile 1511 CHARACTER (LEN=86) :: roots_chr !< String for root profile 1512 CHARACTER (LEN=86) :: vertical_index_chr !< String for the vertical index 1513 CHARACTER (LEN=86) :: m_soil_chr !< String for soil moisture 1514 CHARACTER (LEN=86) :: soil_depth_chr !< String for soil depth 1515 CHARACTER (LEN=10) :: coor_chr !< Temporary string 1516 1517 INTEGER(iwp) :: i !< Loop index over soil layers 1518 1519 INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file 1520 1521 t_soil_chr = '' 1522 m_soil_chr = '' 1523 soil_depth_chr = '' 1524 roots_chr = '' 1525 vertical_index_chr = '' 1526 1527 i = 1 1528 DO i = nzb_soil, nzt_soil 1529 WRITE (coor_chr,'(F10.2,7X)') soil_temperature(i) 1530 t_soil_chr = TRIM( t_soil_chr ) // ' ' // TRIM( coor_chr ) 1531 1532 WRITE (coor_chr,'(F10.2,7X)') soil_moisture(i) 1533 m_soil_chr = TRIM( m_soil_chr ) // ' ' // TRIM( coor_chr ) 1534 1535 WRITE (coor_chr,'(F10.2,7X)') - zs(i) 1536 soil_depth_chr = TRIM( soil_depth_chr ) // ' ' // TRIM( coor_chr ) 1537 1538 WRITE (coor_chr,'(F10.2,7X)') root_fraction(i) 1539 roots_chr = TRIM( roots_chr ) // ' ' // TRIM( coor_chr ) 1540 1541 WRITE (coor_chr,'(I10,7X)') i 1542 vertical_index_chr = TRIM( vertical_index_chr ) // ' ' // & 1543 TRIM( coor_chr ) 1544 ENDDO 1545 1546 ! 1547 !-- Write land surface model header 1548 WRITE( io, 1 ) 1549 IF ( conserve_water_content ) THEN 1550 WRITE( io, 2 ) 1551 ELSE 1552 WRITE( io, 3 ) 1553 ENDIF 1554 1555 WRITE( io, 4 ) TRIM( veg_type_name(veg_type) ), & 1556 TRIM (soil_type_name(soil_type) ) 1557 WRITE( io, 5 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ), & 1558 TRIM( m_soil_chr ), TRIM( roots_chr ), & 1559 TRIM( vertical_index_chr ) 1560 1561 1 FORMAT (//' Land surface model information:'/ & 1562 ' ------------------------------'/) 1563 2 FORMAT (' --> Soil bottom is closed (water content is conserved, default)') 1564 3 FORMAT (' --> Soil bottom is open (water content is not conserved)') 1565 4 FORMAT (' --> Land surface type : ',A,/ & 1566 ' --> Soil porosity type : ',A) 1567 5 FORMAT (/' Initial soil temperature and moisture profile:'// & 1568 ' Height: ',A,' m'/ & 1569 ' Temperature: ',A,' K'/ & 1570 ' Moisture: ',A,' m**3/m**3'/ & 1571 ' Root fraction: ',A,' '/ & 1572 ' Gridpoint: ',A) 1573 1574 1575 1576 END SUBROUTINE lsm_header 1577 1578 1579 !------------------------------------------------------------------------------! 1580 ! Description: 1581 ! ------------ 1582 !> Initialization of the land surface model 1583 !------------------------------------------------------------------------------! 1584 SUBROUTINE lsm_init 1585 1586 1587 IMPLICIT NONE 1588 1589 INTEGER(iwp) :: i !< running index 1590 INTEGER(iwp) :: j !< running index 1591 INTEGER(iwp) :: k !< running index 1592 1593 REAL(wp) :: pt1 !< potential temperature at first grid level 1594 1595 1596 ! 1597 !-- Calculate Exner function 1598 exn = ( surface_pressure / 1000.0_wp )**0.286_wp 1599 1600 1601 ! 1602 !-- If no cloud physics is used, rho_surface has not been calculated before 1603 IF ( .NOT. cloud_physics ) THEN 1604 rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn ) 1605 ENDIF 1606 1607 ! 1608 !-- Calculate frequently used parameters 1609 rho_cp = cp * rho_surface 1610 rd_d_rv = r_d / r_v 1611 rho_lv = rho_surface * l_v 1612 drho_l_lv = 1.0_wp / (rho_l * l_v) 1613 1614 ! 1615 !-- Set inital values for prognostic quantities 1616 tt_surface_m = 0.0_wp 1617 tt_soil_m = 0.0_wp 1618 tm_soil_m = 0.0_wp 1619 tm_liq_eb_m = 0.0_wp 1620 c_liq = 0.0_wp 1621 1622 ghf_eb = 0.0_wp 1623 shf_eb = rho_cp * shf 1624 1625 IF ( humidity ) THEN 1626 qsws_eb = rho_l * l_v * qsws 1627 ELSE 1628 qsws_eb = 0.0_wp 1629 ENDIF 1630 1631 qsws_liq_eb = 0.0_wp 1632 qsws_soil_eb = 0.0_wp 1633 qsws_veg_eb = 0.0_wp 1634 1635 r_a = 50.0_wp 1636 r_s = 50.0_wp 1637 r_canopy = 0.0_wp 1638 r_soil = 0.0_wp 1639 1640 ! 1641 !-- Allocate 3D soil model arrays 1642 ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 1643 ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 1644 ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 1645 1646 lambda_h = 0.0_wp 1647 ! 1648 !-- If required, allocate humidity-related variables for the soil model 1649 IF ( humidity ) THEN 1650 ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 1651 ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 1652 1653 lambda_w = 0.0_wp 1654 ENDIF 1655 1656 ! 1657 !-- Calculate grid spacings. Temperature and moisture are defined at 1658 !-- the edges of the soil layers (_stag), whereas gradients/fluxes are defined 1659 !-- at the centers 1660 dz_soil(nzb_soil) = zs(nzb_soil) 1661 1662 DO k = nzb_soil+1, nzt_soil 1663 dz_soil(k) = zs(k) - zs(k-1) 1664 ENDDO 1665 dz_soil(nzt_soil+1) = dz_soil(nzt_soil) 1666 1667 DO k = nzb_soil, nzt_soil-1 1668 dz_soil_stag(k) = 0.5_wp * (dz_soil(k+1) + dz_soil(k)) 1669 ENDDO 1670 dz_soil_stag(nzt_soil) = dz_soil(nzt_soil) 1671 1672 ddz_soil = 1.0_wp / dz_soil 1673 ddz_soil_stag = 1.0_wp / dz_soil_stag 1674 1675 ! 1676 !-- Initialize standard soil types. It is possible to overwrite each 1677 !-- parameter by setting the respecticy NAMELIST variable to a 1678 !-- value /= 9999999.9. 1679 IF ( soil_type /= 0 ) THEN 1680 1681 IF ( alpha_vangenuchten == 9999999.9_wp ) THEN 1682 alpha_vangenuchten = soil_pars(0,soil_type) 1683 ENDIF 1684 1685 IF ( l_vangenuchten == 9999999.9_wp ) THEN 1686 l_vangenuchten = soil_pars(1,soil_type) 1687 ENDIF 1688 1689 IF ( n_vangenuchten == 9999999.9_wp ) THEN 1690 n_vangenuchten = soil_pars(2,soil_type) 1691 ENDIF 1692 1693 IF ( hydraulic_conductivity == 9999999.9_wp ) THEN 1694 hydraulic_conductivity = soil_pars(3,soil_type) 1695 ENDIF 1696 1697 IF ( saturation_moisture == 9999999.9_wp ) THEN 1698 saturation_moisture = m_soil_pars(0,soil_type) 1699 ENDIF 1700 1701 IF ( field_capacity == 9999999.9_wp ) THEN 1702 field_capacity = m_soil_pars(1,soil_type) 1703 ENDIF 1704 1705 IF ( wilting_point == 9999999.9_wp ) THEN 1706 wilting_point = m_soil_pars(2,soil_type) 1707 ENDIF 1708 1709 IF ( residual_moisture == 9999999.9_wp ) THEN 1710 residual_moisture = m_soil_pars(3,soil_type) 1711 ENDIF 1712 1713 ENDIF 1714 1715 ! 1716 !-- Map values to the respective 2D arrays 1717 alpha_vg = alpha_vangenuchten 1718 l_vg = l_vangenuchten 1719 n_vg = n_vangenuchten 1720 gamma_w_sat = hydraulic_conductivity 1721 m_sat = saturation_moisture 1722 m_fc = field_capacity 1723 m_wilt = wilting_point 1724 m_res = residual_moisture 1725 r_soil_min = min_soil_resistance 1726 1727 ! 1728 !-- Initial run actions 1729 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 1730 1731 t_soil = 0.0_wp 1732 m_liq_eb = 0.0_wp 1733 m_soil = 0.0_wp 1734 1735 ! 1736 !-- Map user settings of T and q for each soil layer 1737 !-- (make sure that the soil moisture does not drop below the permanent 1738 !-- wilting point) -> problems with devision by zero) 1739 DO k = nzb_soil, nzt_soil 1740 t_soil(k,:,:) = soil_temperature(k) 1741 m_soil(k,:,:) = MAX(soil_moisture(k),m_wilt(:,:)) 1742 soil_moisture(k) = MAX(soil_moisture(k),wilting_point) 1743 ENDDO 1744 t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1) 1745 1746 ! 1747 !-- Calculate surface temperature 1748 t_surface = pt_surface * exn 1749 1750 ! 1751 !-- Set artifical values for ts and us so that r_a has its initial value 1752 !-- for the first time step 1753 DO i = nxlg, nxrg 1754 DO j = nysg, nyng 1755 k = nzb_s_inner(j,i) 1756 1757 IF ( cloud_physics ) THEN 1758 pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i) 1759 ELSE 1760 pt1 = pt(k+1,j,i) 1761 ENDIF 1762 1763 ! 1764 !-- Assure that r_a cannot be zero at model start 1765 IF ( pt1 == pt(k,j,i) ) pt1 = pt1 + 1.0E-10_wp 1766 1767 us(j,i) = 0.1_wp 1768 ts(j,i) = (pt1 - pt(k,j,i)) / r_a(j,i) 1769 shf(j,i) = - us(j,i) * ts(j,i) 1770 ENDDO 1771 ENDDO 1772 1773 ! 1774 !-- Actions for restart runs 1775 ELSE 1776 1777 DO i = nxlg, nxrg 1778 DO j = nysg, nyng 1779 k = nzb_s_inner(j,i) 1780 t_surface(j,i) = pt(k,j,i) * exn 1781 ENDDO 1782 ENDDO 1783 1784 ENDIF 1785 1786 DO k = nzb_soil, nzt_soil 1787 root_fr(k,:,:) = root_fraction(k) 1788 ENDDO 1789 1790 IF ( veg_type /= 0 ) THEN 1791 IF ( min_canopy_resistance == 9999999.9_wp ) THEN 1792 min_canopy_resistance = veg_pars(0,veg_type) 1793 ENDIF 1794 IF ( leaf_area_index == 9999999.9_wp ) THEN 1795 leaf_area_index = veg_pars(1,veg_type) 1796 ENDIF 1797 IF ( vegetation_coverage == 9999999.9_wp ) THEN 1798 vegetation_coverage = veg_pars(2,veg_type) 1799 ENDIF 1800 IF ( canopy_resistance_coefficient == 9999999.9_wp ) THEN 1801 canopy_resistance_coefficient= veg_pars(3,veg_type) 1802 ENDIF 1803 IF ( lambda_surface_stable == 9999999.9_wp ) THEN 1804 lambda_surface_stable = surface_pars(0,veg_type) 1805 ENDIF 1806 IF ( lambda_surface_unstable == 9999999.9_wp ) THEN 1807 lambda_surface_unstable = surface_pars(1,veg_type) 1808 ENDIF 1809 IF ( f_shortwave_incoming == 9999999.9_wp ) THEN 1810 f_shortwave_incoming = surface_pars(2,veg_type) 1811 ENDIF 1812 IF ( z0_eb == 9999999.9_wp ) THEN 1813 roughness_length = roughness_par(0,veg_type) 1814 z0_eb = roughness_par(0,veg_type) 1815 ENDIF 1816 IF ( z0h_eb == 9999999.9_wp ) THEN 1817 z0h_eb = roughness_par(1,veg_type) 1818 ENDIF 1819 IF ( z0q_eb == 9999999.9_wp ) THEN 1820 z0q_eb = roughness_par(2,veg_type) 1821 ENDIF 1822 z0h_factor = z0h_eb / ( z0_eb + 1.0E-20_wp ) 1823 1824 IF ( ANY( root_fraction == 9999999.9_wp ) ) THEN 1825 DO k = nzb_soil, nzt_soil 1826 root_fr(k,:,:) = root_distribution(k,veg_type) 1827 root_fraction(k) = root_distribution(k,veg_type) 1828 ENDDO 1829 ENDIF 1830 1831 ELSE 1832 1833 IF ( z0_eb == 9999999.9_wp ) THEN 1834 z0_eb = roughness_length 1835 ENDIF 1836 IF ( z0h_eb == 9999999.9_wp ) THEN 1837 z0h_eb = z0_eb * z0h_factor 1838 ENDIF 1839 IF ( z0q_eb == 9999999.9_wp ) THEN 1840 z0q_eb = z0_eb * z0h_factor 1841 ENDIF 1842 1843 ENDIF 1844 1845 ! 1846 !-- For surfaces covered with pavement, set depth of the pavement (with dry 1847 !-- soil below). The depth must be greater than the first soil layer depth 1848 IF ( veg_type == 20 ) THEN 1849 IF ( pave_depth == 9999999.9_wp ) THEN 1850 pave_depth = zs(nzb_soil) 1851 ELSE 1852 pave_depth = MAX( zs(nzb_soil), pave_depth ) 1853 ENDIF 1854 ENDIF 1855 1856 ! 1857 !-- Map vegetation and soil types to 2D array to allow for heterogeneous 1858 !-- surfaces via user interface see below 1859 veg_type_2d = veg_type 1860 soil_type_2d = soil_type 1861 1862 ! 1863 !-- Map vegetation parameters to the respective 2D arrays 1864 r_canopy_min = min_canopy_resistance 1865 lai = leaf_area_index 1866 c_veg = vegetation_coverage 1867 g_d = canopy_resistance_coefficient 1868 lambda_surface_s = lambda_surface_stable 1869 lambda_surface_u = lambda_surface_unstable 1870 f_sw_in = f_shortwave_incoming 1871 z0 = z0_eb 1872 z0h = z0h_eb 1873 z0q = z0q_eb 1874 1875 ! 1876 !-- Possibly do user-defined actions (e.g. define heterogeneous land surface) 1877 CALL user_init_land_surface 1878 1879 ! 1880 !-- Set flag parameter if vegetation type was set to a water surface. Also 1881 !-- set temperature to a constant value in all "soil" layers. 1882 DO i = nxlg, nxrg 1883 DO j = nysg, nyng 1884 IF ( veg_type_2d(j,i) == 14 .OR. veg_type_2d(j,i) == 15 ) THEN 1885 water_surface(j,i) = .TRUE. 1886 ELSEIF ( veg_type_2d(j,i) == 20 ) THEN 1887 pave_surface(j,i) = .TRUE. 1888 m_soil(:,j,i) = 0.0_wp 1889 ENDIF 1890 1891 ENDDO 1892 ENDDO 1893 1894 ! 1895 !-- Calculate new roughness lengths (for water surfaces only) 1896 CALL calc_z0_water_surface 1897 1898 t_soil_p = t_soil 1899 m_soil_p = m_soil 1900 m_liq_eb_p = m_liq_eb 1901 t_surface_p = t_surface 1902 1903 1904 1905 !-- Store initial profiles of t_soil and m_soil (assuming they are 1906 !-- horizontally homogeneous on this PE) 1907 hom(nzb_soil:nzt_soil,1,90,:) = SPREAD( t_soil(nzb_soil:nzt_soil, & 1908 nysg,nxlg), 2, & 1909 statistic_regions+1 ) 1910 hom(nzb_soil:nzt_soil,1,92,:) = SPREAD( m_soil(nzb_soil:nzt_soil, & 1911 nysg,nxlg), 2, & 1912 statistic_regions+1 ) 1913 1914 END SUBROUTINE lsm_init 573 1915 574 1916 … … 578 1920 !> Allocate land surface model arrays and define pointers 579 1921 !------------------------------------------------------------------------------! 580 SUBROUTINE init_lsm_arrays1922 SUBROUTINE lsm_init_arrays 581 1923 582 1924 … … 658 2000 659 2001 660 END SUBROUTINE init_lsm_arrays 2002 END SUBROUTINE lsm_init_arrays 2003 661 2004 662 2005 !------------------------------------------------------------------------------! 663 2006 ! Description: 664 2007 ! ------------ 665 !> Initialization of theland surface model2008 !> Parin for &lsmpar for land surface model 666 2009 !------------------------------------------------------------------------------! 667 SUBROUTINE init_lsm668 2010 SUBROUTINE lsm_parin 2011 669 2012 670 2013 IMPLICIT NONE 671 2014 672 INTEGER(iwp) :: i !< running index 673 INTEGER(iwp) :: j !< running index 674 INTEGER(iwp) :: k !< running index 675 676 REAL(wp) :: pt1 !< potential temperature at first grid level 677 678 679 ! 680 !-- Calculate Exner function 681 exn = ( surface_pressure / 1000.0_wp )**0.286_wp 682 683 684 ! 685 !-- If no cloud physics is used, rho_surface has not been calculated before 686 IF ( .NOT. cloud_physics ) THEN 687 rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn ) 688 ENDIF 689 690 ! 691 !-- Calculate frequently used parameters 692 rho_cp = cp * rho_surface 693 rd_d_rv = r_d / r_v 694 rho_lv = rho_surface * l_v 695 drho_l_lv = 1.0_wp / (rho_l * l_v) 696 697 ! 698 !-- Set inital values for prognostic quantities 699 tt_surface_m = 0.0_wp 700 tt_soil_m = 0.0_wp 701 tm_soil_m = 0.0_wp 702 tm_liq_eb_m = 0.0_wp 703 c_liq = 0.0_wp 704 705 ghf_eb = 0.0_wp 706 shf_eb = rho_cp * shf 707 708 IF ( humidity ) THEN 709 qsws_eb = rho_l * l_v * qsws 710 ELSE 711 qsws_eb = 0.0_wp 712 ENDIF 713 714 qsws_liq_eb = 0.0_wp 715 qsws_soil_eb = 0.0_wp 716 qsws_veg_eb = 0.0_wp 717 718 r_a = 50.0_wp 719 r_s = 50.0_wp 720 r_canopy = 0.0_wp 721 r_soil = 0.0_wp 722 723 ! 724 !-- Allocate 3D soil model arrays 725 ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 726 ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 727 ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 728 729 lambda_h = 0.0_wp 730 ! 731 !-- If required, allocate humidity-related variables for the soil model 732 IF ( humidity ) THEN 733 ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 734 ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 735 736 lambda_w = 0.0_wp 737 ENDIF 738 739 ! 740 !-- Calculate grid spacings. Temperature and moisture are defined at 741 !-- the edges of the soil layers (_stag), whereas gradients/fluxes are defined 742 !-- at the centers 743 dz_soil(nzb_soil) = zs(nzb_soil) 744 745 DO k = nzb_soil+1, nzt_soil 746 dz_soil(k) = zs(k) - zs(k-1) 2015 CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file 2016 2017 NAMELIST /lsm_par/ alpha_vangenuchten, c_surface, & 2018 canopy_resistance_coefficient, & 2019 conserve_water_content, & 2020 f_shortwave_incoming, field_capacity, & 2021 hydraulic_conductivity, & 2022 lambda_surface_stable, & 2023 lambda_surface_unstable, leaf_area_index, & 2024 l_vangenuchten, min_canopy_resistance, & 2025 min_soil_resistance, n_vangenuchten, & 2026 pave_depth, pave_heat_capacity, & 2027 pave_heat_conductivity, & 2028 residual_moisture, root_fraction, & 2029 saturation_moisture, skip_time_do_lsm, & 2030 soil_moisture, soil_temperature, soil_type, & 2031 vegetation_coverage, veg_type, wilting_point,& 2032 zs, z0_eb, z0h_eb, z0q_eb 2033 2034 line = ' ' 2035 2036 ! 2037 !-- Try to find land surface model package 2038 REWIND ( 11 ) 2039 line = ' ' 2040 DO WHILE ( INDEX( line, '&lsm_par' ) == 0 ) 2041 READ ( 11, '(A)', END=10 ) line 747 2042 ENDDO 748 dz_soil(nzt_soil+1) = dz_soil(nzt_soil) 749 750 DO k = nzb_soil, nzt_soil-1 751 dz_soil_stag(k) = 0.5_wp * (dz_soil(k+1) + dz_soil(k)) 752 ENDDO 753 dz_soil_stag(nzt_soil) = dz_soil(nzt_soil) 754 755 ddz_soil = 1.0_wp / dz_soil 756 ddz_soil_stag = 1.0_wp / dz_soil_stag 757 758 ! 759 !-- Initialize standard soil types. It is possible to overwrite each 760 !-- parameter by setting the respecticy NAMELIST variable to a 761 !-- value /= 9999999.9. 762 IF ( soil_type /= 0 ) THEN 763 764 IF ( alpha_vangenuchten == 9999999.9_wp ) THEN 765 alpha_vangenuchten = soil_pars(0,soil_type) 766 ENDIF 767 768 IF ( l_vangenuchten == 9999999.9_wp ) THEN 769 l_vangenuchten = soil_pars(1,soil_type) 770 ENDIF 771 772 IF ( n_vangenuchten == 9999999.9_wp ) THEN 773 n_vangenuchten = soil_pars(2,soil_type) 774 ENDIF 775 776 IF ( hydraulic_conductivity == 9999999.9_wp ) THEN 777 hydraulic_conductivity = soil_pars(3,soil_type) 778 ENDIF 779 780 IF ( saturation_moisture == 9999999.9_wp ) THEN 781 saturation_moisture = m_soil_pars(0,soil_type) 782 ENDIF 783 784 IF ( field_capacity == 9999999.9_wp ) THEN 785 field_capacity = m_soil_pars(1,soil_type) 786 ENDIF 787 788 IF ( wilting_point == 9999999.9_wp ) THEN 789 wilting_point = m_soil_pars(2,soil_type) 790 ENDIF 791 792 IF ( residual_moisture == 9999999.9_wp ) THEN 793 residual_moisture = m_soil_pars(3,soil_type) 794 ENDIF 795 796 ENDIF 797 798 ! 799 !-- Map values to the respective 2D arrays 800 alpha_vg = alpha_vangenuchten 801 l_vg = l_vangenuchten 802 n_vg = n_vangenuchten 803 gamma_w_sat = hydraulic_conductivity 804 m_sat = saturation_moisture 805 m_fc = field_capacity 806 m_wilt = wilting_point 807 m_res = residual_moisture 808 r_soil_min = min_soil_resistance 809 810 ! 811 !-- Initial run actions 812 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 813 814 t_soil = 0.0_wp 815 m_liq_eb = 0.0_wp 816 m_soil = 0.0_wp 817 818 ! 819 !-- Map user settings of T and q for each soil layer 820 !-- (make sure that the soil moisture does not drop below the permanent 821 !-- wilting point) -> problems with devision by zero) 822 DO k = nzb_soil, nzt_soil 823 t_soil(k,:,:) = soil_temperature(k) 824 m_soil(k,:,:) = MAX(soil_moisture(k),m_wilt(:,:)) 825 soil_moisture(k) = MAX(soil_moisture(k),wilting_point) 826 ENDDO 827 t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1) 828 829 ! 830 !-- Calculate surface temperature 831 t_surface = pt_surface * exn 832 833 ! 834 !-- Set artifical values for ts and us so that r_a has its initial value 835 !-- for the first time step 836 DO i = nxlg, nxrg 837 DO j = nysg, nyng 838 k = nzb_s_inner(j,i) 839 840 IF ( cloud_physics ) THEN 841 pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i) 842 ELSE 843 pt1 = pt(k+1,j,i) 844 ENDIF 845 846 ! 847 !-- Assure that r_a cannot be zero at model start 848 IF ( pt1 == pt(k,j,i) ) pt1 = pt1 + 1.0E-10_wp 849 850 us(j,i) = 0.1_wp 851 ts(j,i) = (pt1 - pt(k,j,i)) / r_a(j,i) 852 shf(j,i) = - us(j,i) * ts(j,i) 853 ENDDO 854 ENDDO 855 856 ! 857 !-- Actions for restart runs 858 ELSE 859 860 DO i = nxlg, nxrg 861 DO j = nysg, nyng 862 k = nzb_s_inner(j,i) 863 t_surface(j,i) = pt(k,j,i) * exn 864 ENDDO 865 ENDDO 866 867 ENDIF 868 869 DO k = nzb_soil, nzt_soil 870 root_fr(k,:,:) = root_fraction(k) 871 ENDDO 872 873 IF ( veg_type /= 0 ) THEN 874 IF ( min_canopy_resistance == 9999999.9_wp ) THEN 875 min_canopy_resistance = veg_pars(0,veg_type) 876 ENDIF 877 IF ( leaf_area_index == 9999999.9_wp ) THEN 878 leaf_area_index = veg_pars(1,veg_type) 879 ENDIF 880 IF ( vegetation_coverage == 9999999.9_wp ) THEN 881 vegetation_coverage = veg_pars(2,veg_type) 882 ENDIF 883 IF ( canopy_resistance_coefficient == 9999999.9_wp ) THEN 884 canopy_resistance_coefficient= veg_pars(3,veg_type) 885 ENDIF 886 IF ( lambda_surface_stable == 9999999.9_wp ) THEN 887 lambda_surface_stable = surface_pars(0,veg_type) 888 ENDIF 889 IF ( lambda_surface_unstable == 9999999.9_wp ) THEN 890 lambda_surface_unstable = surface_pars(1,veg_type) 891 ENDIF 892 IF ( f_shortwave_incoming == 9999999.9_wp ) THEN 893 f_shortwave_incoming = surface_pars(2,veg_type) 894 ENDIF 895 IF ( z0_eb == 9999999.9_wp ) THEN 896 roughness_length = roughness_par(0,veg_type) 897 z0_eb = roughness_par(0,veg_type) 898 ENDIF 899 IF ( z0h_eb == 9999999.9_wp ) THEN 900 z0h_eb = roughness_par(1,veg_type) 901 ENDIF 902 IF ( z0q_eb == 9999999.9_wp ) THEN 903 z0q_eb = roughness_par(2,veg_type) 904 ENDIF 905 z0h_factor = z0h_eb / ( z0_eb + 1.0E-20_wp ) 906 907 IF ( ANY( root_fraction == 9999999.9_wp ) ) THEN 908 DO k = nzb_soil, nzt_soil 909 root_fr(k,:,:) = root_distribution(k,veg_type) 910 root_fraction(k) = root_distribution(k,veg_type) 911 ENDDO 912 ENDIF 913 914 ELSE 915 916 IF ( z0_eb == 9999999.9_wp ) THEN 917 z0_eb = roughness_length 918 ENDIF 919 IF ( z0h_eb == 9999999.9_wp ) THEN 920 z0h_eb = z0_eb * z0h_factor 921 ENDIF 922 IF ( z0q_eb == 9999999.9_wp ) THEN 923 z0q_eb = z0_eb * z0h_factor 924 ENDIF 925 926 ENDIF 927 928 ! 929 !-- For surfaces covered with pavement, set depth of the pavement (with dry 930 !-- soil below). The depth must be greater than the first soil layer depth 931 IF ( veg_type == 20 ) THEN 932 IF ( pave_depth == 9999999.9_wp ) THEN 933 pave_depth = zs(nzb_soil) 934 ELSE 935 pave_depth = MAX( zs(nzb_soil), pave_depth ) 936 ENDIF 937 ENDIF 938 939 ! 940 !-- Map vegetation and soil types to 2D array to allow for heterogeneous 941 !-- surfaces via user interface see below 942 veg_type_2d = veg_type 943 soil_type_2d = soil_type 944 945 ! 946 !-- Map vegetation parameters to the respective 2D arrays 947 r_canopy_min = min_canopy_resistance 948 lai = leaf_area_index 949 c_veg = vegetation_coverage 950 g_d = canopy_resistance_coefficient 951 lambda_surface_s = lambda_surface_stable 952 lambda_surface_u = lambda_surface_unstable 953 f_sw_in = f_shortwave_incoming 954 z0 = z0_eb 955 z0h = z0h_eb 956 z0q = z0q_eb 957 958 ! 959 !-- Possibly do user-defined actions (e.g. define heterogeneous land surface) 960 CALL user_init_land_surface 961 962 ! 963 !-- Set flag parameter if vegetation type was set to a water surface. Also 964 !-- set temperature to a constant value in all "soil" layers. 965 DO i = nxlg, nxrg 966 DO j = nysg, nyng 967 IF ( veg_type_2d(j,i) == 14 .OR. veg_type_2d(j,i) == 15 ) THEN 968 water_surface(j,i) = .TRUE. 969 ELSEIF ( veg_type_2d(j,i) == 20 ) THEN 970 pave_surface(j,i) = .TRUE. 971 m_soil(:,j,i) = 0.0_wp 972 ENDIF 973 974 ENDDO 975 ENDDO 976 977 ! 978 !-- Calculate new roughness lengths (for water surfaces only) 979 CALL calc_z0_water_surface 980 981 t_soil_p = t_soil 982 m_soil_p = m_soil 983 m_liq_eb_p = m_liq_eb 984 t_surface_p = t_surface 985 986 987 988 !-- Store initial profiles of t_soil and m_soil (assuming they are 989 !-- horizontally homogeneous on this PE) 990 hom(nzb_soil:nzt_soil,1,90,:) = SPREAD( t_soil(nzb_soil:nzt_soil, & 991 nysg,nxlg), 2, & 992 statistic_regions+1 ) 993 hom(nzb_soil:nzt_soil,1,92,:) = SPREAD( m_soil(nzb_soil:nzt_soil, & 994 nysg,nxlg), 2, & 995 statistic_regions+1 ) 996 997 END SUBROUTINE init_lsm 998 999 1000 1001 !------------------------------------------------------------------------------! 1002 ! Description: 1003 ! ------------ 1004 !> Solver for the energy balance at the surface. 1005 !------------------------------------------------------------------------------! 1006 SUBROUTINE lsm_energy_balance 1007 1008 1009 IMPLICIT NONE 1010 1011 INTEGER(iwp) :: i !< running index 1012 INTEGER(iwp) :: j !< running index 1013 INTEGER(iwp) :: k, ks !< running index 1014 1015 REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface 1016 f1, & !< resistance correction term 1 1017 f2, & !< resistance correction term 2 1018 f3, & !< resistance correction term 3 1019 m_min, & !< minimum soil moisture 1020 e, & !< water vapour pressure 1021 e_s, & !< water vapour saturation pressure 1022 e_s_dt, & !< derivate of e_s with respect to T 1023 tend, & !< tendency 1024 dq_s_dt, & !< derivate of q_s with respect to T 1025 coef_1, & !< coef. for prognostic equation 1026 coef_2, & !< coef. for prognostic equation 1027 f_qsws, & !< factor for qsws_eb 1028 f_qsws_veg, & !< factor for qsws_veg_eb 1029 f_qsws_soil, & !< factor for qsws_soil_eb 1030 f_qsws_liq, & !< factor for qsws_liq_eb 1031 f_shf, & !< factor for shf_eb 1032 lambda_surface, & !< Current value of lambda_surface 1033 m_liq_eb_max, & !< maxmimum value of the liq. water reservoir 1034 pt1, & !< potential temperature at first grid level 1035 qv1 !< specific humidity at first grid level 1036 1037 ! 1038 !-- Calculate the exner function for the current time step 1039 exn = ( surface_pressure / 1000.0_wp )**0.286_wp 1040 1041 DO i = nxlg, nxrg 1042 DO j = nysg, nyng 1043 k = nzb_s_inner(j,i) 1044 1045 ! 1046 !-- Set lambda_surface according to stratification between skin layer and soil 1047 IF ( .NOT. pave_surface(j,i) ) THEN 1048 1049 c_surface_tmp = c_surface 1050 1051 IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i)) THEN 1052 lambda_surface = lambda_surface_s(j,i) 1053 ELSE 1054 lambda_surface = lambda_surface_u(j,i) 1055 ENDIF 1056 ELSE 1057 1058 c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp 1059 lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil) 1060 1061 ENDIF 1062 1063 ! 1064 !-- First step: calculate aerodyamic resistance. As pt, us, ts 1065 !-- are not available for the prognostic time step, data from the last 1066 !-- time step is used here. Note that this formulation is the 1067 !-- equivalent to the ECMWF formulation using drag coefficients 1068 IF ( cloud_physics ) THEN 1069 pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i) 1070 qv1 = q(k+1,j,i) - ql(k+1,j,i) 1071 ELSE 1072 pt1 = pt(k+1,j,i) 1073 qv1 = q(k+1,j,i) 1074 ENDIF 1075 1076 r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp) 1077 1078 ! 1079 !-- Make sure that the resistance does not drop to zero 1080 IF ( ABS(r_a(j,i)) < 1.0E-10_wp ) r_a(j,i) = 1.0E-10_wp 1081 1082 ! 1083 !-- Second step: calculate canopy resistance r_canopy 1084 !-- f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation 1085 1086 !-- f1: correction for incoming shortwave radiation (stomata close at 1087 !-- night) 1088 IF ( radiation_scheme /= 'constant' ) THEN 1089 f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) / & 1090 (0.81_wp * (0.004_wp * rad_sw_in(k,j,i) & 1091 + 1.0_wp)) ) 1092 ELSE 1093 f1 = 1.0_wp 1094 ENDIF 1095 1096 1097 ! 1098 !-- f2: correction for soil moisture availability to plants (the 1099 !-- integrated soil moisture must thus be considered here) 1100 !-- f2 = 0 for very dry soils 1101 m_total = 0.0_wp 1102 DO ks = nzb_soil, nzt_soil 1103 m_total = m_total + root_fr(ks,j,i) & 1104 * MAX(m_soil(ks,j,i),m_wilt(j,i)) 1105 ENDDO 1106 1107 IF ( m_total > m_wilt(j,i) .AND. m_total < m_fc(j,i) ) THEN 1108 f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) ) 1109 ELSEIF ( m_total >= m_fc(j,i) ) THEN 1110 f2 = 1.0_wp 1111 ELSE 1112 f2 = 1.0E-20_wp 1113 ENDIF 1114 1115 ! 1116 !-- Calculate water vapour pressure at saturation 1117 e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i) & 1118 - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) ) 1119 1120 ! 1121 !-- f3: correction for vapour pressure deficit 1122 IF ( g_d(j,i) /= 0.0_wp ) THEN 1123 ! 1124 !-- Calculate vapour pressure 1125 e = qv1 * surface_pressure / 0.622_wp 1126 f3 = EXP ( -g_d(j,i) * (e_s - e) ) 1127 ELSE 1128 f3 = 1.0_wp 1129 ENDIF 1130 1131 ! 1132 !-- Calculate canopy resistance. In case that c_veg is 0 (bare soils), 1133 !-- this calculation is obsolete, as r_canopy is not used below. 1134 !-- To do: check for very dry soil -> r_canopy goes to infinity 1135 r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3 & 1136 + 1.0E-20_wp) 1137 1138 ! 1139 !-- Third step: calculate bare soil resistance r_soil. The Clapp & 1140 !-- Hornberger parametrization does not consider c_veg. 1141 IF ( soil_type_2d(j,i) /= 7 ) THEN 1142 m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) * & 1143 m_res(j,i) 1144 ELSE 1145 m_min = m_wilt(j,i) 1146 ENDIF 1147 1148 f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min ) 1149 f2 = MAX(f2,1.0E-20_wp) 1150 f2 = MIN(f2,1.0_wp) 1151 1152 r_soil(j,i) = r_soil_min(j,i) / f2 1153 1154 ! 1155 !-- Calculate the maximum possible liquid water amount on plants and 1156 !-- bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is 1157 !-- assumed, while paved surfaces might hold up 1 mm of water. The 1158 !-- liquid water fraction for paved surfaces is calculated after 1159 !-- Noilhan & Planton (1989), while the ECMWF formulation is used for 1160 !-- vegetated surfaces and bare soils. 1161 IF ( pave_surface(j,i) ) THEN 1162 m_liq_eb_max = m_max_depth * 5.0_wp 1163 c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 ) 1164 ELSE 1165 m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i) & 1166 + (1.0_wp - c_veg(j,i)) ) 1167 c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max ) 1168 ENDIF 1169 1170 ! 1171 !-- Calculate saturation specific humidity 1172 q_s = 0.622_wp * e_s / surface_pressure 1173 1174 ! 1175 !-- In case of dewfall, set evapotranspiration to zero 1176 !-- All super-saturated water is then removed from the air 1177 IF ( humidity .AND. q_s <= qv1 ) THEN 1178 r_canopy(j,i) = 0.0_wp 1179 r_soil(j,i) = 0.0_wp 1180 ENDIF 1181 1182 ! 1183 !-- Calculate coefficients for the total evapotranspiration 1184 !-- In case of water surface, set vegetation and soil fluxes to zero. 1185 !-- For pavements, only evaporation of liquid water is possible. 1186 IF ( water_surface(j,i) ) THEN 1187 f_qsws_veg = 0.0_wp 1188 f_qsws_soil = 0.0_wp 1189 f_qsws_liq = rho_lv / r_a(j,i) 1190 ELSEIF ( pave_surface (j,i) ) THEN 1191 f_qsws_veg = 0.0_wp 1192 f_qsws_soil = 0.0_wp 1193 f_qsws_liq = rho_lv * c_liq(j,i) / r_a(j,i) 1194 ELSE 1195 f_qsws_veg = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/ & 1196 (r_a(j,i) + r_canopy(j,i)) 1197 f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) + & 1198 r_soil(j,i)) 1199 f_qsws_liq = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i) 1200 ENDIF 1201 ! 1202 !-- If soil moisture is below wilting point, plants do no longer 1203 !-- transpirate. 1204 ! IF ( m_soil(k,j,i) < m_wilt(j,i) ) THEN 1205 ! f_qsws_veg = 0.0_wp 1206 ! ENDIF 1207 1208 f_shf = rho_cp / r_a(j,i) 1209 f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq 1210 1211 ! 1212 !-- Calculate derivative of q_s for Taylor series expansion 1213 e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) - & 1214 17.269_wp*(t_surface(j,i) - 273.16_wp) & 1215 / (t_surface(j,i) - 35.86_wp)**2 ) 1216 1217 dq_s_dt = 0.622_wp * e_s_dt / surface_pressure 1218 1219 ! 1220 !-- Add LW up so that it can be removed in prognostic equation 1221 rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i) 1222 1223 ! 1224 !-- Calculate new skin temperature 1225 IF ( humidity ) THEN 1226 #if defined ( __rrtmg ) 1227 ! 1228 !-- Numerator of the prognostic equation 1229 coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1230 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1231 + f_shf * pt1 + f_qsws * ( qv1 - q_s & 1232 + dq_s_dt * t_surface(j,i) ) + lambda_surface & 1233 * t_soil(nzb_soil,j,i) 1234 1235 ! 1236 !-- Denominator of the prognostic equation 1237 coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt & 1238 + lambda_surface + f_shf / exn 1239 #else 1240 1241 ! 1242 !-- Numerator of the prognostic equation 1243 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1244 * t_surface(j,i) ** 4 & 1245 + f_shf * pt1 + f_qsws * ( qv1 & 1246 - q_s + dq_s_dt * t_surface(j,i) ) & 1247 + lambda_surface * t_soil(nzb_soil,j,i) 1248 1249 ! 1250 !-- Denominator of the prognostic equation 1251 coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws & 1252 * dq_s_dt + lambda_surface + f_shf / exn 1253 1254 #endif 1255 ELSE 1256 1257 #if defined ( __rrtmg ) 1258 ! 1259 !-- Numerator of the prognostic equation 1260 coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1261 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1262 + f_shf * pt1 + lambda_surface & 1263 * t_soil(nzb_soil,j,i) 1264 1265 ! 1266 !-- Denominator of the prognostic equation 1267 coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn 1268 #else 1269 1270 ! 1271 !-- Numerator of the prognostic equation 1272 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1273 * t_surface(j,i) ** 4 + f_shf * pt1 & 1274 + lambda_surface * t_soil(nzb_soil,j,i) 1275 1276 ! 1277 !-- Denominator of the prognostic equation 1278 coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 & 1279 + lambda_surface + f_shf / exn 1280 #endif 1281 ENDIF 1282 1283 tend = 0.0_wp 1284 1285 ! 1286 !-- Implicit solution when the surface layer has no heat capacity, 1287 !-- otherwise use RK3 scheme. 1288 t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp * & 1289 t_surface(j,i) ) / ( c_surface_tmp + coef_2 & 1290 * dt_3d * tsc(2) ) 1291 1292 ! 1293 !-- Add RK3 term 1294 IF ( c_surface_tmp /= 0.0_wp ) THEN 1295 1296 t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3) & 1297 * tt_surface_m(j,i) 1298 1299 ! 1300 !-- Calculate true tendency 1301 tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3) & 1302 * tt_surface_m(j,i)) / (dt_3d * tsc(2)) 1303 ! 1304 !-- Calculate t_surface tendencies for the next Runge-Kutta step 1305 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1306 IF ( intermediate_timestep_count == 1 ) THEN 1307 tt_surface_m(j,i) = tend 1308 ELSEIF ( intermediate_timestep_count < & 1309 intermediate_timestep_count_max ) THEN 1310 tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp & 1311 * tt_surface_m(j,i) 1312 ENDIF 1313 ENDIF 1314 ENDIF 1315 1316 ! 1317 !-- In case of fast changes in the skin temperature, it is possible to 1318 !-- update the radiative fluxes independently from the prescribed 1319 !-- radiation call frequency. This effectively prevents oscillations, 1320 !-- especially when setting skip_time_do_radiation /= 0. The threshold 1321 !-- value of 0.2 used here is just a first guess. This method should be 1322 !-- revised in the future as tests have shown that the threshold is 1323 !-- often reached, when no oscillations would occur (causes immense 1324 !-- computing time for the radiation code). 1325 IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp .AND. & 1326 unscheduled_radiation_calls ) THEN 1327 force_radiation_call_l = .TRUE. 1328 ENDIF 1329 1330 pt(k,j,i) = t_surface_p(j,i) / exn 1331 1332 ! 1333 !-- Calculate fluxes 1334 #if defined ( __rrtmg ) 1335 rad_net_l(j,i) = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1336 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1337 - rad_lw_out_change_0(j,i) * t_surface_p(j,i) 1338 1339 IF ( rrtm_idrv == 1 ) THEN 1340 rad_net(j,i) = rad_net_l(j,i) 1341 rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) & 1342 + rad_lw_out_change_0(j,i) & 1343 * ( t_surface_p(j,i) - t_surface(j,i) ) 1344 ENDIF 1345 #else 1346 rad_net_l(j,i) = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1347 * t_surface(j,i)**4 - 4.0_wp * sigma_sb & 1348 * t_surface(j,i)**3 * t_surface_p(j,i) 1349 #endif 1350 1351 ghf_eb(j,i) = lambda_surface * (t_surface_p(j,i) & 1352 - t_soil(nzb_soil,j,i)) 1353 1354 shf_eb(j,i) = - f_shf * ( pt1 - pt(k,j,i) ) 1355 1356 shf(j,i) = shf_eb(j,i) / rho_cp 1357 1358 IF ( humidity ) THEN 1359 qsws_eb(j,i) = - f_qsws * ( qv1 - q_s + dq_s_dt & 1360 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) 1361 1362 qsws(j,i) = qsws_eb(j,i) / rho_lv 1363 1364 qsws_veg_eb(j,i) = - f_qsws_veg * ( qv1 - q_s & 1365 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1366 * t_surface_p(j,i) ) 1367 1368 qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s & 1369 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1370 * t_surface_p(j,i) ) 1371 1372 qsws_liq_eb(j,i) = - f_qsws_liq * ( qv1 - q_s & 1373 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1374 * t_surface_p(j,i) ) 1375 ENDIF 1376 1377 ! 1378 !-- Calculate the true surface resistance 1379 IF ( qsws_eb(j,i) == 0.0_wp ) THEN 1380 r_s(j,i) = 1.0E10_wp 1381 ELSE 1382 r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt & 1383 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) & 1384 / qsws_eb(j,i) - r_a(j,i) 1385 ENDIF 1386 1387 ! 1388 !-- Calculate change in liquid water reservoir due to dew fall or 1389 !-- evaporation of liquid water 1390 IF ( humidity ) THEN 1391 ! 1392 !-- If precipitation is activated, add rain water to qsws_liq_eb 1393 !-- and qsws_soil_eb according the the vegetation coverage. 1394 !-- precipitation_rate is given in mm. 1395 IF ( precipitation ) THEN 1396 1397 ! 1398 !-- Add precipitation to liquid water reservoir, if possible. 1399 !-- Otherwise, add the water to soil. In case of 1400 !-- pavements, the exceeding water amount is implicitely removed 1401 !-- as runoff as qsws_soil_eb is then not used in the soil model 1402 IF ( m_liq_eb(j,i) /= m_liq_eb_max ) THEN 1403 qsws_liq_eb(j,i) = qsws_liq_eb(j,i) & 1404 + c_veg(j,i) * prr(k,j,i) * hyrho(k) & 1405 * 0.001_wp * rho_l * l_v 1406 ELSE 1407 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) & 1408 + c_veg(j,i) * prr(k,j,i) * hyrho(k) & 1409 * 0.001_wp * rho_l * l_v 1410 ENDIF 1411 1412 !-- Add precipitation to bare soil according to the bare soil 1413 !-- coverage. 1414 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) * (1.0_wp & 1415 - c_veg(j,i)) * prr(k,j,i) * hyrho(k) & 1416 * 0.001_wp * rho_l * l_v 1417 ENDIF 1418 1419 ! 1420 !-- If the air is saturated, check the reservoir water level 1421 IF ( qsws_eb(j,i) < 0.0_wp ) THEN 1422 1423 ! 1424 !-- Check if reservoir is full (avoid values > m_liq_eb_max) 1425 !-- In that case, qsws_liq_eb goes to qsws_soil_eb. In this 1426 !-- case qsws_veg_eb is zero anyway (because c_liq = 1), 1427 !-- so that tend is zero and no further check is needed 1428 IF ( m_liq_eb(j,i) == m_liq_eb_max ) THEN 1429 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) & 1430 + qsws_liq_eb(j,i) 1431 qsws_liq_eb(j,i) = 0.0_wp 1432 ENDIF 1433 1434 ! 1435 !-- In case qsws_veg_eb becomes negative (unphysical behavior), 1436 !-- let the water enter the liquid water reservoir as dew on the 1437 !-- plant 1438 IF ( qsws_veg_eb(j,i) < 0.0_wp ) THEN 1439 qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i) 1440 qsws_veg_eb(j,i) = 0.0_wp 1441 ENDIF 1442 ENDIF 1443 1444 tend = - qsws_liq_eb(j,i) * drho_l_lv 1445 1446 m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend & 1447 + tsc(3) * tm_liq_eb_m(j,i) ) 1448 1449 ! 1450 !-- Check if reservoir is overfull -> reduce to maximum 1451 !-- (conservation of water is violated here) 1452 m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max) 1453 1454 ! 1455 !-- Check if reservoir is empty (avoid values < 0.0) 1456 !-- (conservation of water is violated here) 1457 m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp) 1458 1459 1460 ! 1461 !-- Calculate m_liq_eb tendencies for the next Runge-Kutta step 1462 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1463 IF ( intermediate_timestep_count == 1 ) THEN 1464 tm_liq_eb_m(j,i) = tend 1465 ELSEIF ( intermediate_timestep_count < & 1466 intermediate_timestep_count_max ) THEN 1467 tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp & 1468 * tm_liq_eb_m(j,i) 1469 ENDIF 1470 ENDIF 1471 1472 ENDIF 1473 1474 ENDDO 1475 ENDDO 1476 1477 ! 1478 !-- Make a logical OR for all processes. Force radiation call if at 1479 !-- least one processor reached the threshold change in skin temperature 1480 IF ( unscheduled_radiation_calls .AND. intermediate_timestep_count & 1481 == intermediate_timestep_count_max-1 ) THEN 1482 #if defined( __parallel ) 1483 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1484 CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call, & 1485 1, MPI_LOGICAL, MPI_LOR, comm2d, ierr ) 1486 #else 1487 force_radiation_call = force_radiation_call_l 1488 #endif 1489 force_radiation_call_l = .FALSE. 1490 ENDIF 1491 1492 ! 1493 !-- Calculate surface specific humidity 1494 IF ( humidity ) THEN 1495 CALL calc_q_surface 1496 ENDIF 1497 1498 ! 1499 !-- Calculate new roughness lengths (for water surfaces only) 1500 CALL calc_z0_water_surface 1501 1502 1503 END SUBROUTINE lsm_energy_balance 2043 BACKSPACE ( 11 ) 2044 2045 ! 2046 !-- Read user-defined namelist 2047 READ ( 11, lsm_par ) 2048 2049 ! 2050 !-- Set flag that indicates that the land surface model is switched on 2051 land_surface = .TRUE. 2052 2053 10 CONTINUE 2054 2055 2056 END SUBROUTINE lsm_parin 1504 2057 1505 2058 … … 1555 2108 lambda_h_water ** m_soil(k,j,i) 1556 2109 1557 ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i))) 2110 ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) & 2111 / m_sat(j,i))) 1558 2112 1559 2113 lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) + & … … 1687 2241 ! 1688 2242 ! 1689 !-- In case of a closed bottom (= water content is conserved), set1690 !-- hydraulic conductivity to zero to that no water will be lost1691 !-- in the bottom layer.2243 !-- In case of a closed bottom (= water content is conserved), 2244 !-- set hydraulic conductivity to zero to that no water will be 2245 !-- lost in the bottom layer. 1692 2246 IF ( conserve_water_content ) THEN 1693 2247 gamma_w(nzt_soil,j,i) = 0.0_wp … … 1696 2250 ENDIF 1697 2251 1698 !-- The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v)) 1699 !-- ensures the mass conservation for water. The transpiration of 1700 !-- plants equals the cumulative withdrawals by the roots in the 1701 !-- soil. The scheme takes into account the availability of water 1702 !-- in the soil layers as well as the root fraction in the 1703 !-- respective layer. Layer with moisture below wilting point will 1704 !-- not contribute, which reflects the preference of plants to 1705 !-- take water from moister layers. 1706 1707 ! 1708 !-- Calculate the root extraction (ECMWF 7.69, the sum of root_extr 1709 !-- = 1). The energy balance solver guarantees a positive 1710 !-- transpiration, so that there is no need for an additional check. 2252 !-- The root extraction (= root_extr * qsws_veg_eb / (rho_l 2253 !-- * l_v)) ensures the mass conservation for water. The 2254 !-- transpiration of plants equals the cumulative withdrawals by 2255 !-- the roots in the soil. The scheme takes into account the 2256 !-- availability of water in the soil layers as well as the root 2257 !-- fraction in the respective layer. Layer with moisture below 2258 !-- wilting point will not contribute, which reflects the 2259 !-- preference of plants to take water from moister layers. 2260 2261 ! 2262 !-- Calculate the root extraction (ECMWF 7.69, the sum of 2263 !-- root_extr = 1). The energy balance solver guarantees a 2264 !-- positive transpiration, so that there is no need for an 2265 !-- additional check. 1711 2266 DO k = nzb_soil, nzt_soil 1712 2267 IF ( m_soil(k,j,i) > m_wilt(j,i) ) THEN … … 1718 2273 DO k = nzb_soil, nzt_soil 1719 2274 IF ( m_soil(k,j,i) > m_wilt(j,i) ) THEN 1720 root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total 2275 root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) & 2276 / m_total 1721 2277 ELSE 1722 2278 root_extr(k) = 0.0_wp … … 1792 2348 END SUBROUTINE lsm_soil_model 1793 2349 2350 2351 !------------------------------------------------------------------------------! 2352 ! Description: 2353 ! ------------ 2354 !> Swapping of timelevels 2355 !------------------------------------------------------------------------------! 2356 SUBROUTINE lsm_swap_timelevel ( mod_count ) 2357 2358 IMPLICIT NONE 2359 2360 INTEGER, INTENT(IN) :: mod_count 2361 2362 #if defined( __nopointer ) 2363 2364 t_surface = t_surface_p 2365 t_soil = t_soil_p 2366 IF ( humidity ) THEN 2367 m_soil = m_soil_p 2368 m_liq_eb = m_liq_eb_p 2369 ENDIF 2370 2371 #else 2372 2373 SELECT CASE ( mod_count ) 2374 2375 CASE ( 0 ) 2376 2377 t_surface => t_surface_1; t_surface_p => t_surface_2 2378 t_soil => t_soil_1; t_soil_p => t_soil_2 2379 IF ( humidity ) THEN 2380 m_soil => m_soil_1; m_soil_p => m_soil_2 2381 m_liq_eb => m_liq_eb_1; m_liq_eb_p => m_liq_eb_2 2382 ENDIF 2383 2384 2385 CASE ( 1 ) 2386 2387 t_surface => t_surface_2; t_surface_p => t_surface_1 2388 t_soil => t_soil_2; t_soil_p => t_soil_1 2389 IF ( humidity ) THEN 2390 m_soil => m_soil_2; m_soil_p => m_soil_1 2391 m_liq_eb => m_liq_eb_2; m_liq_eb_p => m_liq_eb_1 2392 ENDIF 2393 2394 END SELECT 2395 #endif 2396 2397 END SUBROUTINE lsm_swap_timelevel 2398 1794 2399 1795 2400 !------------------------------------------------------------------------------! … … 1909 2514 1910 2515 1911 !------------------------------------------------------------------------------!1912 ! Description:1913 ! ------------1914 !> Swapping of timelevels1915 !------------------------------------------------------------------------------!1916 SUBROUTINE lsm_swap_timelevel ( mod_count )1917 1918 IMPLICIT NONE1919 1920 INTEGER, INTENT(IN) :: mod_count1921 1922 #if defined( __nopointer )1923 1924 t_surface = t_surface_p1925 t_soil = t_soil_p1926 IF ( humidity ) THEN1927 m_soil = m_soil_p1928 m_liq_eb = m_liq_eb_p1929 ENDIF1930 1931 #else1932 1933 SELECT CASE ( mod_count )1934 1935 CASE ( 0 )1936 1937 t_surface => t_surface_1; t_surface_p => t_surface_21938 t_soil => t_soil_1; t_soil_p => t_soil_21939 IF ( humidity ) THEN1940 m_soil => m_soil_1; m_soil_p => m_soil_21941 m_liq_eb => m_liq_eb_1; m_liq_eb_p => m_liq_eb_21942 ENDIF1943 1944 1945 CASE ( 1 )1946 1947 t_surface => t_surface_2; t_surface_p => t_surface_11948 t_soil => t_soil_2; t_soil_p => t_soil_11949 IF ( humidity ) THEN1950 m_soil => m_soil_2; m_soil_p => m_soil_11951 m_liq_eb => m_liq_eb_2; m_liq_eb_p => m_liq_eb_11952 ENDIF1953 1954 END SELECT1955 #endif1956 1957 END SUBROUTINE lsm_swap_timelevel1958 1959 1960 2516 END MODULE land_surface_model_mod -
palm/trunk/SOURCE/package_parin.f90
r1789 r1817 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Reading of &lsmpar moved to land_surface_model_mod. 22 22 ! 23 23 ! Former revisions: … … 129 129 vc_alpha, vc_gradient_normals, vc_mode, vc_size_x, vc_size_y, & 130 130 vc_size_z 131 132 USE land_surface_model_mod, &133 ONLY: alpha_vangenuchten, c_surface, canopy_resistance_coefficient, &134 conserve_water_content, field_capacity, &135 f_shortwave_incoming, hydraulic_conductivity, &136 lambda_surface_stable, lambda_surface_unstable, leaf_area_index, &137 land_surface, l_vangenuchten, min_canopy_resistance, &138 min_soil_resistance, n_vangenuchten, pave_heat_capacity, &139 pave_depth, pave_heat_conductivity, residual_moisture, &140 root_fraction, saturation_moisture, wilting_point, &141 skip_time_do_lsm, soil_moisture, soil_temperature, soil_type, &142 vegetation_coverage, veg_type, zs, z0_eb, z0h_eb, z0q_eb143 131 144 132 USE kinds … … 206 194 vc_gradient_normals, vc_mode, vc_size_x, & 207 195 vc_size_y, vc_size_z 208 209 NAMELIST /lsm_par/ alpha_vangenuchten, c_surface, &210 canopy_resistance_coefficient, &211 conserve_water_content, &212 f_shortwave_incoming, field_capacity, &213 hydraulic_conductivity, &214 lambda_surface_stable, &215 lambda_surface_unstable, leaf_area_index, &216 l_vangenuchten, min_canopy_resistance, &217 min_soil_resistance, n_vangenuchten, &218 pave_depth, pave_heat_capacity, &219 pave_heat_conductivity, &220 residual_moisture, root_fraction, &221 saturation_moisture, skip_time_do_lsm, &222 soil_moisture, soil_temperature, soil_type, &223 vegetation_coverage, veg_type, wilting_point,&224 zs, z0_eb, z0h_eb, z0q_eb225 226 227 196 228 197 NAMELIST /particles_par/ alloc_factor, bc_par_b, bc_par_lr, & … … 345 314 40 CONTINUE 346 315 347 !348 !-- Try to find land surface model package349 REWIND ( 11 )350 line = ' '351 DO WHILE ( INDEX( line, '&lsm_par' ) == 0 )352 READ ( 11, '(A)', END=50 ) line353 ENDDO354 BACKSPACE ( 11 )355 356 !357 !-- Read user-defined namelist358 READ ( 11, lsm_par )359 360 !361 !-- Set flag that indicates that the land surface model is switched on362 land_surface = .TRUE.363 364 50 CONTINUE365 316 366 317 ! -
palm/trunk/SOURCE/parin.f90
r1805 r1817 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Added call to land surface model for reading of &lsmpar 22 22 ! 23 23 ! Former revisions: … … 208 208 USE kinds 209 209 210 USE land_surface_model_mod, & 211 ONLY: lsm_parin 212 210 213 USE model_1d, & 211 214 ONLY: damp_level_1d, dt_pr_1d, dt_run_control_1d, end_time_1d … … 450 453 ENDIF 451 454 455 456 ! 457 !-- Check if land surface model is used and read &lsmpar if required 458 CALL lsm_parin 459 460 452 461 ! 453 462 !-- Read control parameters for optionally used model software packages
Note: See TracChangeset
for help on using the changeset viewer.
