Changeset 1817


Ignore:
Timestamp:
Apr 6, 2016 3:44:20 PM (6 years ago)
Author:
maronga
Message:

further modularization of land surface model

Location:
palm/trunk
Files:
7 edited
1 moved

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/mbuild

    r1805 r1817  
    2222# Current revisions:
    2323# ------------------
    24 #
     24# Bugfix: Makefile_check removed (parameter file check)
    2525#
    2626# Former revisions:
     
    427427    cd  $local_source_path
    428428    printf "\n\n  *** tar of makefile and source files in $local_source_path" 
    429  
    430     cat Makefile_check|while read line
    431     do
    432        line=$(echo $line|grep SOURCES)
    433        if [[ $line == *"SOURCES"* ]]
    434        then
    435           line=$(echo $line|sed 's/SOURCES = //g')
    436           break
    437        fi
    438     done
    439 
    440     tar -cf  ${mainprog}_sources_check.tar  Makefile_check  $line
    441429    printf "\n"
    442430 fi
     
    14811469    # FINAL WORK
    14821470 rm -f  hosts_found_in_config_file
    1483  rm -f ${local_source_path}/${mainprog}_sources_check.tar
  • palm/trunk/SOURCE/Makefile

    r1809 r1817  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# Renamed land_surface_model to land_surface_model_mod.
    2323#
    2424# Former revisions:
     
    259259        init_grid.f90 init_masks.f90 init_ocean.f90 \
    260260        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 \
    262262        local_stop.f90 local_system.f90 local_tremain.f90 \
    263263        local_tremain_ini.f90 lpm.f90 lpm_advec.f90 lpm_boundary_conds.f90 \
     
    340340advec_w_pw.o: modules.o mod_kinds.o
    341341advec_w_up.o: modules.o mod_kinds.o
    342 average_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
     342average_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model_mod.o\
    343343                   radiation_model.o
    344344boundary_conds.o: modules.o mod_kinds.o
     
    350350check_for_restart.o: modules.o mod_kinds.o pmc_interface.o
    351351check_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 \
     352check_parameters.o: modules.o mod_kinds.o land_surface_model_mod.o \
    353353        netcdf_interface.o plant_canopy_model.o pmc_interface.o radiation_model.o \
    354354        subsidence.o
     
    370370data_output_tseries.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o
    371371data_output_2d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o\
    372    netcdf_interface.o land_surface_model.o radiation_model.o
     372   netcdf_interface.o land_surface_model_mod.o radiation_model.o
    373373data_output_3d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o\
    374    netcdf_interface.o land_surface_model.o
     374   netcdf_interface.o land_surface_model_mod.o
    375375diffusion_e.o: modules.o mod_kinds.o
    376376diffusion_s.o: modules.o mod_kinds.o
     
    386386exchange_horiz_2d.o: modules.o cpulog.o mod_kinds.o
    387387fft_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\
     388flow_statistics.o: modules.o cpulog.o mod_kinds.o land_surface_model_mod.o\
    389389   netcdf_interface.o radiation_model.o
    390390global_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\
     391header.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o land_surface_model_mod.o\
    392392   plant_canopy_model.o pmc_handle_communicator.o pmc_interface.o \
    393393   radiation_model.o spectrum.o subsidence.o
     
    395395inflow_turbulence.o: modules.o cpulog.o mod_kinds.o
    396396init_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 \
     397init_3d_model.o: modules.o mod_kinds.o advec_ws.o cpulog.o land_surface_model_mod.o \
    398398   lpm_init.o ls_forcing.o netcdf_interface.o plant_canopy_model.o \
    399399   radiation_model.o random_function.o random_generator_parallel.o \
     
    411411init_slope.o: modules.o mod_kinds.o
    412412interaction_droplets_ptq.o: modules.o mod_kinds.o
    413 land_surface_model.o: modules.o mod_kinds.o radiation_model.o
     413land_surface_model_mod.o: modules.o mod_kinds.o radiation_model.o
    414414local_stop.o: modules.o mod_kinds.o pmc_interface.o
    415415local_tremain.o: modules.o cpulog.o mod_kinds.o
     
    452452mod_particle_attributes.o: mod_particle_attributes.f90 mod_kinds.o
    453453netcdf_interface.o: netcdf_interface.f90 modules.o mod_kinds.o \
    454    land_surface_model.o spectrum.o
     454   land_surface_model_mod.o spectrum.o
    455455nudging.o: modules.o cpulog.o mod_kinds.o
    456 package_parin.o: modules.o mod_kinds.o land_surface_model.o \
     456package_parin.o: modules.o mod_kinds.o \
    457457   plant_canopy_model.o radiation_model.o spectrum.o
    458458palm.o: modules.o cpulog.o ls_forcing.o mod_kinds.o nudging.o\
    459459        pmc_interface.o surface_layer_fluxes.o
    460 parin.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o \
     460parin.o: modules.o cpulog.o land_surface_model_mod.o mod_kinds.o netcdf_interface.o \
    461461   pmc_interface.o progress_bar.o
    462462plant_canopy_model.o: modules.o mod_kinds.o
     
    487487random_gauss.o: mod_kinds.o random_function.o random_generator_parallel.o
    488488random_generator_parallel.o: mod_kinds.o
    489 read_3d_binary.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
     489read_3d_binary.o: modules.o cpulog.o mod_kinds.o land_surface_model_mod.o\
    490490                  radiation_model.o random_function.o\
    491491                  random_generator_parallel.o
     
    497497spectrum.o: spectrum.f90 modules.o mod_kinds.o cpulog.o fft_xy.o
    498498subsidence.o: modules.o mod_kinds.o
    499 sum_up_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
     499sum_up_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model_mod.o\
    500500                  radiation_model.o
    501501surface_coupler.o: modules.o cpulog.o mod_kinds.o
    502 surface_layer_fluxes.o: modules.o mod_kinds.o land_surface_model.o
    503 swap_timelevel.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
     502surface_layer_fluxes.o: modules.o mod_kinds.o land_surface_model_mod.o
     503swap_timelevel.o: modules.o cpulog.o mod_kinds.o land_surface_model_mod.o\
    504504   pmc_interface.o
    505505temperton_fft.o: modules.o mod_kinds.o
    506506time_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 \
    508508        ls_forcing.o mod_kinds.o nudging.o pmc_interface.o production_e.o \
    509509        prognostic_equations.o progress_bar.o radiation_model.o \
     
    531531user_init_3d_model.o: modules.o mod_kinds.o user_module.o
    532532user_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.o
     533user_init_land_surface.o: modules.o mod_kinds.o user_module.o land_surface_model_mod.o
    534534user_init_plant_canopy.o: modules.o mod_kinds.o user_module.o plant_canopy_model.o
    535535user_init_radiation.o: modules.o mod_kinds.o user_module.o radiation_model.o
     
    544544user_statistics.o: modules.o mod_kinds.o netcdf_interface.o user_module.o
    545545wall_fluxes.o: modules.o mod_kinds.o
    546 write_3d_binary.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
     546write_3d_binary.o: modules.o cpulog.o mod_kinds.o land_surface_model_mod.o\
    547547                   radiation_model.o random_function.o\
    548548                   random_generator_parallel.o
  • palm/trunk/SOURCE/check_parameters.f90

    r1812 r1817  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Moved checks for land_surface model to the respective module
    2222!
    2323! Former revisions:
     
    330330    USE grid_variables
    331331    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
    333335    USE kinds
    334336    USE model_1d
     
    10351037    ENDIF
    10361038
    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
    12121043
    12131044    IF ( radiation )  THEN
     
    31442975             ENDIF
    31452976
    3146           CASE ( 't_soil', '#t_soil' )
    3147              IF (  .NOT.  land_surface )  THEN
    3148                 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              ELSE
    3153                 dopr_index(i) = 89
    3154                 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) == '#' )  THEN
    3157                    dopr_initial_index(i) = 90
    3158                    hom(0:nzs-1,2,90,:)   = SPREAD( - zs, 2, statistic_regions+1 )
    3159                    data_output_pr(i)     = data_output_pr(i)(2:)
    3160                 ENDIF
    3161              ENDIF
    3162 
    3163           CASE ( 'm_soil', '#m_soil' )
    3164              IF (  .NOT.  land_surface )  THEN
    3165                 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              ELSE
    3170                 dopr_index(i) = 91
    3171                 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) == '#' )  THEN
    3174                    dopr_initial_index(i) = 92
    3175                    hom(0:nzs-1,2,92,:)   = SPREAD( - zs, 2, statistic_regions+1 )
    3176                    data_output_pr(i)     = data_output_pr(i)(2:)
    3177                 ENDIF
    3178              ENDIF
    3179 
    31802977          CASE ( 'rad_net' )
    31812978             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
     
    33063103          CASE DEFAULT
    33073104
    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
    33093110
    33103111             IF ( unit == 'illegal' )  THEN
     
    33943195             unit = 'K'
    33953196
    3396           CASE ( 'm_soil' )
    3397              IF (  .NOT.  land_surface )  THEN
    3398                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    3399                          'land_surface = .TRUE.'
    3400                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    3401              ENDIF
    3402              unit = 'm3/m3'
    3403 
    34043197          CASE ( 'nr' )
    34053198             IF (  .NOT.  cloud_physics )  THEN
     
    35363329             unit = 'psu'
    35373330
    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*',   &
    35503340                 '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*' )
    35533342             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
    35543343                message_string = 'illegal value for data_output: "' //         &
     
    35693358                ENDIF
    35703359             ENDIF
    3571 
    3572              IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT.  land_surface )  THEN
    3573                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    3574                                  'res land_surface = .TRUE.'
    3575                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    3576              ENDIF
    3577              IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT.  land_surface )  THEN
    3578                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    3579                                  'res land_surface = .TRUE.'
    3580                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    3581              ENDIF
    3582              IF ( TRIM( var ) == 'c_veg*'  .AND.  .NOT. land_surface )  THEN
    3583                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    3584                                  'res land_surface = .TRUE.'
    3585                 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    3586              ENDIF
    3587              IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT.  land_surface )  THEN
    3588                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    3589                                  'res land_surface = .TRUE.'
    3590                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    3591              ENDIF
    35923360             IF ( TRIM( var ) == 'lai*'  .AND.  .NOT.  land_surface )  THEN
    35933361                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     
    36003368                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
    36013369             ENDIF
    3602              IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT.  land_surface )  THEN
    3603                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    3604                                  'res land_surface = .TRUE.'
    3605                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    3606              ENDIF
    36073370             IF ( TRIM( var ) == 'pra*'  .AND.  .NOT. precipitation )  THEN
    36083371                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     
    36253388                CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
    36263389             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
    36673391             IF ( TRIM( var ) == 'lai*'   )  unit = 'none'
    36683392             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/m2'
    3669              IF ( TRIM( var ) == 'm_liq_eb*'     )  unit = 'm'
    36703393             IF ( TRIM( var ) == 'ol*'   )   unit = 'm'
    36713394             IF ( TRIM( var ) == 'pra*'   )  unit = 'mm'
    36723395             IF ( TRIM( var ) == 'prr*'   )  unit = 'mm/s'
    36733396             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'
    36783397             IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'   
    36793398             IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''   
     
    36813400             IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = ''
    36823401             IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = ''
    3683              IF ( TRIM( var ) == 'r_a*')     unit = 's/m'     
    3684              IF ( TRIM( var ) == 'r_s*')     unit = 's/m'
    36853402             IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
    3686              IF ( TRIM( var ) == 'shf_eb*')  unit = 'W/m2'
    36873403             IF ( TRIM( var ) == 't*'     )  unit = 'K'
    36883404             IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
    36893405             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             
    37013408          CASE DEFAULT
    37023409
    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
    37043415
    37053416             IF ( unit == 'illegal' )  THEN
  • palm/trunk/SOURCE/header.f90

    r1809 r1817  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Moved land_surface_model header output to the respective module.
    2222!
    2323! Former revisions:
     
    258258       
    259259    USE kinds
    260   
     260 
    261261    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
    265263 
    266264    USE model_1d,                                                              &
     
    969967    ENDIF
    970968
    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 )
    1016970
    1017971    IF ( radiation )  THEN
     
    23082262418 FORMAT (/'    Day of the year at model start :   day_init      =     ',I3 &
    23092263            /'    UTC time at model start        :   time_utc_init = ',F7.1' s')
    2310 419 FORMAT (//' Land surface model information:'/ &
    2311               ' ------------------------------'/)
    23122264420 FORMAT (/'    Characteristic levels of the initial temperature profile:'// &
    23132265            '       Height:        ',A,'  m'/ &
     
    23682320                       '[0,1000] cm**2/s**3')
    23692321437 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)')
    23802322444 FORMAT (//' Radiation model information:'/                                 &
    23812323              ' ----------------------------'/)
  • palm/trunk/SOURCE/init_3d_model.f90

    r1816 r1817  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Renamed lsm calls.
    2222!
    2323! Former revisions:
     
    6969!
    7070! 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,
    7272! which is part of land_surface_model.
    7373!
     
    266266
    267267    USE land_surface_model_mod,                                                &
    268         ONLY:  init_lsm, init_lsm_arrays, land_surface
     268        ONLY:  lsm_init, lsm_init_arrays, land_surface
    269269 
    270270    USE ls_forcing_mod
     
    678678!-- Allocate land surface model arrays
    679679    IF ( land_surface )  THEN
    680        CALL init_lsm_arrays
     680       CALL lsm_init_arrays
    681681    ENDIF
    682682
     
    16651665    IF ( land_surface )  THEN
    16661666       CALL location_message( 'initializing land surface model', .FALSE. )
    1667        CALL init_lsm
     1667       CALL lsm_init
    16681668       CALL location_message( 'finished', .TRUE. )
    16691669    ENDIF
  • palm/trunk/SOURCE/land_surface_model_mod.f90

    r1812 r1817  
    1 !> @file land_surface_model.f90
     1!> @file land_surface_model_mod.f90
    22!--------------------------------------------------------------------------------!
    33! This file is part of PALM.
     
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Added interface for lsm_init_arrays. Added subroutines for check_parameters,
     22! header, and parin. Renamed some subroutines.
    2223!
    2324! Former revisions:
     
    139140        ONLY:  cloud_physics, dt_3d, humidity, intermediate_timestep_count,    &
    140141               initializing_actions, intermediate_timestep_count_max,          &
    141                max_masks, precipitation, pt_surface, rho_surface,              &
    142                roughness_length, surface_pressure, timestep_scheme, tsc,       &
    143                z0h_factor, time_since_reference_point
     142               max_masks, precipitation, pt_surface,                           &
     143               rho_surface, roughness_length, surface_pressure,                &
     144               timestep_scheme, tsc, z0h_factor, time_since_reference_point
    144145
    145146    USE indices,                                                               &
     
    523524    PRIVATE
    524525
    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
    526532!
    527533!-- Public parameters, constants and initial values
    528534    PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient,       &
    529535           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            rho_lv, root_fraction, saturation_moisture, skip_time_do_lsm,       &
    537            soil_moisture, soil_temperature, soil_type, soil_type_name,         &
    538            vegetation_coverage, veg_type, veg_type_name, wilting_point, z0_eb, &
    539            z0h_eb, z0q_eb
     536           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
    540546
    541547!
     
    554560    PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av
    555561
    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   
    560575    INTERFACE lsm_energy_balance
    561576       MODULE PROCEDURE lsm_energy_balance
    562577    END INTERFACE lsm_energy_balance
    563578
     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   
    564595    INTERFACE lsm_soil_model
    565596       MODULE PROCEDURE lsm_soil_model
     
    571602
    572603 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
     15611   FORMAT (//' Land surface model information:'/                              &
     1562              ' ------------------------------'/)
     15632   FORMAT (' --> Soil bottom is closed (water content is conserved, default)')
     15643   FORMAT (' --> Soil bottom is open (water content is not conserved)')         
     15654   FORMAT (' --> Land surface type  : ',A,/                                   &
     1566            ' --> Soil porosity type : ',A)
     15675   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
    5731915
    5741916
     
    5781920!> Allocate land surface model arrays and define pointers
    5791921!------------------------------------------------------------------------------!
    580     SUBROUTINE init_lsm_arrays
     1922    SUBROUTINE lsm_init_arrays
    5811923   
    5821924
     
    6582000
    6592001
    660     END SUBROUTINE init_lsm_arrays
     2002    END SUBROUTINE lsm_init_arrays
     2003
    6612004
    6622005!------------------------------------------------------------------------------!
    6632006! Description:
    6642007! ------------
    665 !> Initialization of the land surface model
     2008!> Parin for &lsmpar for land surface model
    6662009!------------------------------------------------------------------------------!
    667     SUBROUTINE init_lsm
    668    
     2010    SUBROUTINE lsm_parin
     2011
    6692012
    6702013       IMPLICIT NONE
    6712014
    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
    7472042       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
    15042057
    15052058
     
    15552108                                     lambda_h_water ** m_soil(k,j,i)
    15562109
    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)))
    15582112
    15592113                      lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +    &
     
    16872241!
    16882242!
    1689 !--                In case of a closed bottom (= water content is conserved), set
    1690 !--                hydraulic conductivity to zero to that no water will be lost
    1691 !--                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.
    16922246                   IF ( conserve_water_content )  THEN
    16932247                      gamma_w(nzt_soil,j,i) = 0.0_wp
     
    16962250                   ENDIF     
    16972251
    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.
    17112266                   DO  k = nzb_soil, nzt_soil
    17122267                       IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
     
    17182273                      DO  k = nzb_soil, nzt_soil
    17192274                         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
    17212277                         ELSE
    17222278                            root_extr(k) = 0.0_wp
     
    17922348    END SUBROUTINE lsm_soil_model
    17932349
     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
    17942399
    17952400!------------------------------------------------------------------------------!
     
    19092514
    19102515
    1911 !------------------------------------------------------------------------------!
    1912 ! Description:
    1913 ! ------------
    1914 !> Swapping of timelevels
    1915 !------------------------------------------------------------------------------!
    1916     SUBROUTINE lsm_swap_timelevel ( mod_count )
    1917 
    1918        IMPLICIT NONE
    1919 
    1920        INTEGER, INTENT(IN) :: mod_count
    1921 
    1922 #if defined( __nopointer )
    1923 
    1924        t_surface    = t_surface_p
    1925        t_soil       = t_soil_p
    1926        IF ( humidity )  THEN
    1927           m_soil    = m_soil_p
    1928           m_liq_eb  = m_liq_eb_p
    1929        ENDIF
    1930 
    1931 #else
    1932    
    1933        SELECT CASE ( mod_count )
    1934 
    1935           CASE ( 0 )
    1936 
    1937              t_surface  => t_surface_1; t_surface_p  => t_surface_2
    1938              t_soil     => t_soil_1;    t_soil_p     => t_soil_2
    1939              IF ( humidity )  THEN
    1940                 m_soil    => m_soil_1;   m_soil_p    => m_soil_2
    1941                 m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_2
    1942              ENDIF
    1943 
    1944 
    1945           CASE ( 1 )
    1946 
    1947              t_surface  => t_surface_2; t_surface_p  => t_surface_1
    1948              t_soil     => t_soil_2;    t_soil_p     => t_soil_1
    1949              IF ( humidity )  THEN
    1950                 m_soil    => m_soil_2;   m_soil_p    => m_soil_1
    1951                 m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
    1952              ENDIF
    1953 
    1954        END SELECT
    1955 #endif
    1956 
    1957     END SUBROUTINE lsm_swap_timelevel
    1958 
    1959 
    19602516 END MODULE land_surface_model_mod
  • palm/trunk/SOURCE/package_parin.f90

    r1789 r1817  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Reading of &lsmpar moved to land_surface_model_mod.
    2222!
    2323! Former revisions:
     
    129129               vc_alpha, vc_gradient_normals, vc_mode, vc_size_x, vc_size_y,   &
    130130               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_eb
    143131
    144132    USE kinds
     
    206194                                  vc_gradient_normals, vc_mode, vc_size_x,     &
    207195                                  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_eb
    225 
    226 
    227196
    228197    NAMELIST /particles_par/      alloc_factor, bc_par_b, bc_par_lr,           &
     
    345314 40 CONTINUE
    346315
    347 !
    348 !-- Try to find land surface model package
    349     REWIND ( 11 )
    350     line = ' '
    351     DO   WHILE ( INDEX( line, '&lsm_par' ) == 0 )
    352        READ ( 11, '(A)', END=50 )  line
    353     ENDDO
    354     BACKSPACE ( 11 )
    355 
    356 !
    357 !-- Read user-defined namelist
    358     READ ( 11, lsm_par )
    359 
    360 !
    361 !-- Set flag that indicates that the land surface model is switched on
    362     land_surface = .TRUE.
    363 
    364  50 CONTINUE
    365316
    366317!
  • palm/trunk/SOURCE/parin.f90

    r1805 r1817  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Added call to land surface model for reading of &lsmpar
    2222!
    2323! Former revisions:
     
    208208    USE kinds
    209209
     210    USE land_surface_model_mod,                                                &
     211        ONLY: lsm_parin
     212   
    210213    USE model_1d,                                                              &
    211214        ONLY:  damp_level_1d, dt_pr_1d, dt_run_control_1d, end_time_1d
     
    450453          ENDIF
    451454
     455         
     456!
     457!--       Check if land surface model is used and read &lsmpar if required
     458          CALL lsm_parin
     459         
     460         
    452461!
    453462!--       Read control parameters for optionally used model software packages
Note: See TracChangeset for help on using the changeset viewer.