Changeset 1551


Ignore:
Timestamp:
Mar 3, 2015 2:18:16 PM (7 years ago)
Author:
maronga
Message:

land surface model released

Location:
palm/trunk
Files:
26 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r1518 r1551  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# Bugfix: further adjustments for the land surface model and radiation model
    2323#
    2424# Former revisions:
     
    293293advec_w_pw.o: modules.o mod_kinds.o
    294294advec_w_up.o: modules.o mod_kinds.o
    295 average_3d_data.o: modules.o cpulog.o mod_kinds.o
     295average_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
     296                   radiation_model.o
    296297boundary_conds.o: modules.o mod_kinds.o
    297298buoyancy.o: modules.o mod_kinds.o
     
    318319data_output_spectra.o: modules.o cpulog.o mod_kinds.o
    319320data_output_tseries.o: modules.o cpulog.o mod_kinds.o
    320 data_output_2d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
    321 data_output_3d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
     321data_output_2d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o\
     322                  land_surface_model.o radiation_model.o
     323data_output_3d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o\
     324                  land_surface_model.o
    322325diffusion_e.o: modules.o mod_kinds.o
    323326diffusion_s.o: modules.o mod_kinds.o
     
    326329diffusion_w.o: modules.o mod_kinds.o wall_fluxes.o
    327330diffusivities.o: modules.o mod_kinds.o
    328 disturb_field.o: modules.o cpulog.o mod_kinds.o random_function.o random_generator_parallel.o
     331disturb_field.o: modules.o cpulog.o mod_kinds.o random_function.o\
     332                 random_generator_parallel.o
    329333disturb_heatflux.o: modules.o cpulog.o mod_kinds.o
    330334eqn_state_seawater.o: modules.o mod_kinds.o
     
    332336exchange_horiz_2d.o: modules.o cpulog.o mod_kinds.o
    333337fft_xy.o: cuda_fft_interfaces.o modules.o mod_kinds.o singleton.o temperton_fft.o
    334 flow_statistics.o: modules.o cpulog.o mod_kinds.o
     338flow_statistics.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
     339                   radiation_model.o
    335340global_min_max.o: modules.o mod_kinds.o
    336 header.o: modules.o cpulog.o mod_kinds.o plant_canopy_model.o subsidence.o
     341header.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
     342          plant_canopy_model.o radiation_model.o subsidence.o
    337343impact_of_latent_heat.o: modules.o mod_kinds.o
    338344inflow_turbulence.o: modules.o cpulog.o mod_kinds.o
    339345init_1d_model.o: modules.o mod_kinds.o
    340 init_3d_model.o: modules.o cpulog.o mod_kinds.o random_function.o advec_ws.o \
    341         land_surface_model.o ls_forcing.o lpm_init.o plant_canopy_model.o \
     346init_3d_model.o: modules.o cpulog.o mod_kinds.o random_function.o advec_ws.o\
     347        land_surface_model.o ls_forcing.o lpm_init.o plant_canopy_model.o\
    342348        radiation_model.o random_generator_parallel.o
    343349init_advec.o: modules.o mod_kinds.o
     
    395401mod_kinds.o: mod_kinds.f90
    396402mod_particle_attributes.o: mod_particle_attributes.f90 mod_kinds.o
    397 netcdf.o: modules.o mod_kinds.o
     403netcdf.o: modules.o mod_kinds.o land_surface_model.o
    398404nudging.o: modules.o cpulog.o mod_kinds.o
    399 package_parin.o: modules.o mod_kinds.o land_surface_model.o plant_canopy_model.o \
    400                 radiation_model.o
     405package_parin.o: modules.o mod_kinds.o land_surface_model.o\
     406                 plant_canopy_model.o radiation_model.o
    401407palm.o: modules.o cpulog.o ls_forcing.o mod_kinds.o nudging.o
    402408parin.o: modules.o cpulog.o mod_kinds.o progress_bar.o
     
    404410poisfft.o: modules.o cpulog.o fft_xy.o mod_kinds.o tridia_solver.o
    405411poismg.o: modules.o cpulog.o mod_kinds.o
    406 prandtl_fluxes.o: modules.o mod_kinds.o land_surface_model.o
     412prandtl_fluxes.o: modules.o mod_kinds.o
    407413pres.o: modules.o cpulog.o mod_kinds.o poisfft.o
    408414print_1d.o: modules.o cpulog.o mod_kinds.o
     
    419425random_gauss.o: mod_kinds.o random_function.o random_generator_parallel.o
    420426random_generator_parallel.o: mod_kinds.o
    421 read_3d_binary.o: modules.o cpulog.o mod_kinds.o random_function.o random_generator_parallel.o
     427read_3d_binary.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
     428                  radiation_model.o random_function.o\
     429                  random_generator_parallel.o
    422430read_var_list.o: modules.o mod_kinds.o plant_canopy_model.o
    423431run_control.o: modules.o cpulog.o mod_kinds.o
     
    426434sor.o: modules.o mod_kinds.o
    427435subsidence.o: modules.o mod_kinds.o
    428 sum_up_3d_data.o: modules.o cpulog.o mod_kinds.o
     436sum_up_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
     437                  radiation_model.o
    429438surface_coupler.o: modules.o cpulog.o mod_kinds.o
    430439swap_timelevel.o: modules.o cpulog.o mod_kinds.o land_surface_model.o
     
    468477user_statistics.o: modules.o mod_kinds.o user_module.o
    469478wall_fluxes.o: modules.o mod_kinds.o
    470 write_3d_binary.o: modules.o cpulog.o mod_kinds.o random_function.o random_generator_parallel.o
     479write_3d_binary.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
     480                   radiation_model.o random_function.o\
     481                   random_generator_parallel.o
    471482write_var_list.o: modules.o mod_kinds.o plant_canopy_model.o
  • palm/trunk/SOURCE/average_3d_data.f90

    r1323 r1551  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Added support for land surface and radiation model parameters.
    2323!
    2424! Former revisions:
     
    7575    USE kinds
    7676
     77    USE land_surface_model_mod,                                                &
     78        ONLY:  c_liq_av, c_soil_av, c_veg_av, ghf_eb_av, lai_av, m_liq_eb_av,  &
     79               m_soil_av, nzb_soil, nzt_soil, qsws_eb_av, qsws_liq_eb_av,      &
     80               qsws_soil_eb_av, qsws_veg_eb_av, shf_eb_av, t_soil_av
     81
     82    USE radiation_model_mod,                                                   &
     83        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av
    7784
    7885    IMPLICIT NONE
     
    98105       SELECT CASE ( TRIM( doav(ii) ) )
    99106
     107         CASE ( 'c_liq*' )
     108             DO  i = nxlg, nxrg
     109                DO  j = nysg, nyng
     110                   c_liq_av(j,i) = c_liq_av(j,i) / REAL( average_count_3d, KIND=wp )
     111                ENDDO
     112             ENDDO
     113
     114         CASE ( 'c_soil*' )
     115             DO  i = nxlg, nxrg
     116                DO  j = nysg, nyng
     117                   c_soil_av(j,i) = c_soil_av(j,i) / REAL( average_count_3d, KIND=wp )
     118                ENDDO
     119             ENDDO
     120
     121         CASE ( 'c_veg*' )
     122             DO  i = nxlg, nxrg
     123                DO  j = nysg, nyng
     124                   c_veg_av(j,i) = c_veg_av(j,i) / REAL( average_count_3d, KIND=wp )
     125                ENDDO
     126             ENDDO
     127
    100128          CASE ( 'e' )
    101129             DO  i = nxlg, nxrg
     
    107135             ENDDO
    108136
     137         CASE ( 'ghf_eb*' )
     138             DO  i = nxlg, nxrg
     139                DO  j = nysg, nyng
     140                   ghf_eb_av(j,i) = ghf_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
     141                ENDDO
     142             ENDDO
     143
    109144          CASE ( 'qsws*' )
    110145             DO  i = nxlg, nxrg
     
    114149             ENDDO
    115150
     151         CASE ( 'lai*' )
     152             DO  i = nxlg, nxrg
     153                DO  j = nysg, nyng
     154                   lai_av(j,i) = lai_av(j,i) / REAL( average_count_3d, KIND=wp )
     155                ENDDO
     156             ENDDO
     157
    116158          CASE ( 'lpt' )
    117159             DO  i = nxlg, nxrg
     
    127169                DO  j = nysg, nyng
    128170                   lwp_av(j,i) = lwp_av(j,i) / REAL( average_count_3d, KIND=wp )
     171                ENDDO
     172             ENDDO
     173
     174         CASE ( 'm_liq_eb*' )
     175             DO  i = nxlg, nxrg
     176                DO  j = nysg, nyng
     177                   m_liq_eb_av(j,i) = m_liq_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
     178                ENDDO
     179             ENDDO
     180
     181          CASE ( 'm_soil' )
     182             DO  i = nxlg, nxrg
     183                DO  j = nysg, nyng
     184                   DO  k = nzb_soil, nzt_soil
     185                      m_soil_av(k,j,i) = m_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     186                   ENDDO
    129187                ENDDO
    130188             ENDDO
     
    247305             ENDDO
    248306
     307         CASE ( 'qsws_eb*' )
     308             DO  i = nxlg, nxrg
     309                DO  j = nysg, nyng
     310                   qsws_eb_av(j,i) = qsws_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
     311                ENDDO
     312             ENDDO
     313
     314         CASE ( 'qsws_liq_eb*' )
     315             DO  i = nxlg, nxrg
     316                DO  j = nysg, nyng
     317                   qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
     318                ENDDO
     319             ENDDO
     320
     321         CASE ( 'qsws_soil_eb*' )
     322             DO  i = nxlg, nxrg
     323                DO  j = nysg, nyng
     324                   qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
     325                ENDDO
     326             ENDDO
     327
     328         CASE ( 'qsws_veg_eb*' )
     329             DO  i = nxlg, nxrg
     330                DO  j = nysg, nyng
     331                   qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
     332                ENDDO
     333             ENDDO
     334
    249335          CASE ( 'qv' )
    250336             DO  i = nxlg, nxrg
     
    256342             ENDDO
    257343
     344         CASE ( 'rad_sw_in*' )
     345             DO  i = nxlg, nxrg
     346                DO  j = nysg, nyng
     347                   rad_sw_in_av(j,i) = rad_sw_in_av(j,i) / REAL( average_count_3d, KIND=wp )
     348                ENDDO
     349             ENDDO
     350
     351         CASE ( 'rad_net*' )
     352             DO  i = nxlg, nxrg
     353                DO  j = nysg, nyng
     354                   rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, KIND=wp )
     355                ENDDO
     356             ENDDO
     357
    258358          CASE ( 'rho' )
    259359             DO  i = nxlg, nxrg
     
    294394                DO  j = nysg, nyng
    295395                   ts_av(j,i) = ts_av(j,i) / REAL( average_count_3d, KIND=wp )
     396                ENDDO
     397             ENDDO
     398
     399          CASE ( 't_soil' )
     400             DO  i = nxlg, nxrg
     401                DO  j = nysg, nyng
     402                   DO  k = nzb_soil, nzt_soil
     403                      t_soil_av(k,j,i) = t_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     404                   ENDDO
    296405                ENDDO
    297406             ENDDO
  • palm/trunk/SOURCE/check_open.f90

    r1469 r1551  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Removed redundant output for combine_plot_fields
    2323!
    2424! Former revisions:
     
    9999
    100100    USE indices,                                                               &
    101         ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz,   &
    102                nzb, nzt
     101        ONLY:  nbgp, nx, nxl, nxr, ny, nyn, nyng, nys, nysg, nz, nzb, nzt
    103102
    104103    USE kinds
     
    462461!
    463462!--          Specifications for combine_plot_fields
    464              WRITE ( 30 )  -nbgp,nx+nbgp,-nbgp,ny+nbgp, 0 ,nz_do3d
     463             WRITE ( 30 )  -nbgp,nx+nbgp,-nbgp,ny+nbgp
    465464             WRITE ( 30 )  0,nx+1,0,ny+1,0,nz_do3d
    466465#endif
  • palm/trunk/SOURCE/check_parameters.f90

    r1505 r1551  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Added various checks for land surface and radiation model. In the course of this
     23! action, the length of the variable var has to be increased
    2324!
    2425! Former revisions:
     
    278279
    279280    CHARACTER (LEN=1)   ::  sq                       !:
    280     CHARACTER (LEN=6)   ::  var                      !:
     281    CHARACTER (LEN=15)  ::  var                      !:
    281282    CHARACTER (LEN=7)   ::  unit                     !:
    282283    CHARACTER (LEN=8)   ::  date                     !:
     
    970971       IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' )  THEN
    971972          message_string = 'lsm requires setting of'//                         &
    972                            'bc_pt_b = "dirichlet" and '//                        &
     973                           'bc_pt_b = "dirichlet" and '//                      &
    973974                           'bc_q_b  = "dirichlet"'
    974975          CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
     
    988989             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    989990          ENDIF
     991 
     992          IF ( min_canopy_resistance == 9999999.9_wp)  THEN
     993             message_string = 'veg_type = 0 (user_defined)'//                  &
     994                              'requires setting of min_canopy_resistance'//    &
     995                              '/= 9999999.9'
     996             CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
     997          ENDIF
     998
     999          IF ( leaf_area_index == 9999999.9_wp)  THEN
     1000             message_string = 'veg_type = 0 (user_defined)'//                  &
     1001                              'requires setting of leaf_area_index'//          &
     1002                              '/= 9999999.9'
     1003             CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )
     1004          ENDIF
     1005
     1006          IF ( vegetation_coverage == 9999999.9_wp)  THEN
     1007             message_string = 'veg_type = 0 (user_defined)'//                  &
     1008                              'requires setting of vegetation_coverage'//      &
     1009                              '/= 9999999.9'
     1010             CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
     1011          ENDIF
     1012
     1013          IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
     1014             message_string = 'veg_type = 0 (user_defined)'//                  &
     1015                              'requires setting of'//                          &
     1016                              'canopy_resistance_coefficient /= 9999999.9'
     1017             CALL message( 'check_parameters', 'PA0418', 1, 2, 0, 6, 0 )
     1018          ENDIF
     1019
     1020          IF ( lambda_surface_stable == 9999999.9_wp)  THEN
     1021             message_string = 'veg_type = 0 (user_defined)'//                  &
     1022                              'requires setting of lambda_surface_stable'//    &
     1023                              '/= 9999999.9'
     1024             CALL message( 'check_parameters', 'PA0419', 1, 2, 0, 6, 0 )
     1025          ENDIF
     1026
     1027          IF ( lambda_surface_unstable == 9999999.9_wp)  THEN
     1028             message_string = 'veg_type = 0 (user_defined)'//                  &
     1029                              'requires setting of lambda_surface_unstable'//  &
     1030                              '/= 9999999.9'
     1031             CALL message( 'check_parameters', 'PA0420', 1, 2, 0, 6, 0 )
     1032          ENDIF
     1033
     1034          IF ( f_shortwave_incoming == 9999999.9_wp)  THEN
     1035             message_string = 'veg_type = 0 (user_defined)'//                  &
     1036                              'requires setting of f_shortwave_incoming'//     &
     1037                              '/= 9999999.9'
     1038             CALL message( 'check_parameters', 'PA0421', 1, 2, 0, 6, 0 )
     1039          ENDIF
     1040
     1041          IF ( z0_eb == 9999999.9_wp)  THEN
     1042             message_string = 'veg_type = 0 (user_defined)'//                  &
     1043                              'requires setting of z0_eb'//                   &
     1044                              '/= 9999999.9'
     1045             CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 )
     1046          ENDIF
     1047
     1048          IF ( z0h_eb == 9999999.9_wp)  THEN
     1049             message_string = 'veg_type = 0 (user_defined)'//                  &
     1050                              'requires setting of z0h_eb'//                  &
     1051                              '/= 9999999.9'
     1052             CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 )
     1053          ENDIF
     1054
     1055
     1056       ENDIF
     1057
     1058       IF ( soil_type == 0 )  THEN
     1059
     1060          IF ( alpha_vangenuchten == 9999999.9_wp)  THEN
     1061             message_string = 'soil_type = 0 (user_defined)'//                 &
     1062                              'requires setting of alpha_vangenuchten'//       &
     1063                              '/= 9999999.9'
     1064             CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 )
     1065          ENDIF
     1066
     1067          IF ( l_vangenuchten == 9999999.9_wp)  THEN
     1068             message_string = 'soil_type = 0 (user_defined)'//                 &
     1069                              'requires setting of l_vangenuchten'//           &
     1070                              '/= 9999999.9'
     1071             CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 )
     1072          ENDIF
     1073
     1074          IF ( n_vangenuchten == 9999999.9_wp)  THEN
     1075             message_string = 'soil_type = 0 (user_defined)'//                 &
     1076                              'requires setting of n_vangenuchten'//           &
     1077                              '/= 9999999.9'
     1078             CALL message( 'check_parameters', 'PA0424', 1, 2, 0, 6, 0 )
     1079          ENDIF
     1080
     1081          IF ( hydraulic_conductivity == 9999999.9_wp)  THEN
     1082             message_string = 'soil_type = 0 (user_defined)'//                 &
     1083                              'requires setting of hydraulic_conductivity'//   &
     1084                              '/= 9999999.9'
     1085             CALL message( 'check_parameters', 'PA0425', 1, 2, 0, 6, 0 )
     1086          ENDIF
     1087
     1088          IF ( saturation_moisture == 9999999.9_wp)  THEN
     1089             message_string = 'soil_type = 0 (user_defined)'//                 &
     1090                              'requires setting of saturation_moisture'//      &
     1091                              '/= 9999999.9'
     1092             CALL message( 'check_parameters', 'PA0426', 1, 2, 0, 6, 0 )
     1093          ENDIF
     1094
     1095          IF ( field_capacity == 9999999.9_wp)  THEN
     1096             message_string = 'soil_type = 0 (user_defined)'//                 &
     1097                              'requires setting of field_capacity'//           &
     1098                              '/= 9999999.9'
     1099             CALL message( 'check_parameters', 'PA0427', 1, 2, 0, 6, 0 )
     1100          ENDIF
     1101
     1102          IF ( wilting_point == 9999999.9_wp)  THEN
     1103             message_string = 'soil_type = 0 (user_defined)'//                 &
     1104                              'requires setting of wilting_point'//            &
     1105                              '/= 9999999.9'
     1106             CALL message( 'check_parameters', 'PA0428', 1, 2, 0, 6, 0 )
     1107          ENDIF
     1108
     1109          IF ( residual_moisture == 9999999.9_wp)  THEN
     1110             message_string = 'soil_type = 0 (user_defined)'//                 &
     1111                              'requires setting of residual_moisture'//        &
     1112                              '/= 9999999.9'
     1113             CALL message( 'check_parameters', 'PA0429', 1, 2, 0, 6, 0 )
     1114          ENDIF
     1115
    9901116       ENDIF
    9911117
     
    9961122       ENDIF
    9971123
    998 
    9991124    END IF
     1125
     1126    IF ( radiation )  THEN
     1127       IF ( radiation_scheme == 'constant' )  THEN
     1128          irad_scheme = 0
     1129       ELSEIF ( radiation_scheme == 'clear-sky' )  THEN
     1130          irad_scheme = 1
     1131       ELSEIF ( radiation_scheme == 'rrtm' )  THEN
     1132          irad_scheme = 2
     1133       ELSE
     1134          message_string = 'unknown radiation_scheme = '//                     &
     1135                           TRIM( radiation_scheme )
     1136          CALL message( 'check_parameters', 'PA0430', 1, 2, 0, 6, 0 )
     1137       ENDIF
     1138    ENDIF
    10001139
    10011140
     
    28623001             ENDIF
    28633002
     3003          CASE ( 't_soil', '#t_soil' )
     3004             IF ( .NOT. land_surface )  THEN
     3005                message_string = 'data_output_pr = ' //                        &
     3006                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     3007                                 'lemented for land_surface = .FALSE.'
     3008                CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
     3009             ELSE
     3010                dopr_index(i) = 89
     3011                dopr_unit(i)  = 'K'
     3012                hom(0:nzs-1,2,89,:)  = SPREAD( - zs, 2, statistic_regions+1 )
     3013                IF ( data_output_pr(i)(1:1) == '#' )  THEN
     3014                   dopr_initial_index(i) = 90
     3015                   hom(0:nzs-1,2,90,:)   = SPREAD( - zs, 2, statistic_regions+1 )
     3016                   data_output_pr(i)     = data_output_pr(i)(2:)
     3017                ENDIF
     3018             ENDIF
     3019
     3020          CASE ( 'm_soil', '#m_soil' )
     3021             IF ( .NOT. land_surface )  THEN
     3022                message_string = 'data_output_pr = ' //                        &
     3023                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     3024                                 'lemented for land_surface = .FALSE.'
     3025                CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
     3026             ELSE
     3027                dopr_index(i) = 91
     3028                dopr_unit(i)  = 'm3/m3'
     3029                hom(0:nzs-1,2,91,:)  = SPREAD( - zs, 2, statistic_regions+1 )
     3030                IF ( data_output_pr(i)(1:1) == '#' )  THEN
     3031                   dopr_initial_index(i) = 92
     3032                   hom(0:nzs-1,2,92,:)   = SPREAD( - zs, 2, statistic_regions+1 )
     3033                   data_output_pr(i)     = data_output_pr(i)(2:)
     3034                ENDIF
     3035             ENDIF
     3036
    28643037
    28653038          CASE DEFAULT
     
    29323105          ENDIF
    29333106       ENDIF
     3107
    29343108!
    29353109!--    Check for allowed value and set units
     
    29513125             ENDIF
    29523126             unit = 'K'
     3127
     3128          CASE ( 'm_soil' )
     3129             IF ( .NOT. land_surface )  THEN
     3130                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3131                         'land_surface = .TRUE.'
     3132                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     3133             ENDIF
     3134             unit = 'm3/m3'
    29533135
    29543136          CASE ( 'nr' )
     
    30763258             unit = 'psu'
    30773259
    3078           CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'qsws*', 'shf*', 'z0*', 'z0h*' )
     3260          CASE ( 't_soil' )
     3261             IF ( .NOT. land_surface )  THEN
     3262                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3263                         'land_surface = .TRUE.'
     3264                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     3265             ENDIF
     3266             unit = 'K'
     3267
     3268
     3269          CASE ( 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'lai*', 'lwp*',     &
     3270                 'm_liq_eb*', 'pra*', 'prr*', 'qsws*', 'qsws_eb*',             &
     3271                 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*',              &
     3272                 'rad_net*', 'rad_sw_in*', 'shf*', 'shf_eb*', 't*', 'u*',      &
     3273                 'z0*', 'z0h*' )
    30793274             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
    30803275                message_string = 'illegal value for data_output: "' //         &
     
    30833278                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
    30843279             ENDIF
     3280             IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT. land_surface )  THEN
     3281                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3282                                 'res land_surface = .TRUE.'
     3283                CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
     3284             ENDIF
     3285             IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT. land_surface )  THEN
     3286                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3287                                 'res land_surface = .TRUE.'
     3288                CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
     3289             ENDIF
     3290             IF ( TRIM( var ) == 'c_veg*'  .AND.  .NOT. land_surface )  THEN
     3291                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3292                                 'res land_surface = .TRUE.'
     3293                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
     3294             ENDIF
     3295             IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT. land_surface )  THEN
     3296                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3297                                 'res land_surface = .TRUE.'
     3298                CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
     3299             ENDIF
     3300             IF ( TRIM( var ) == 'lai*'  .AND.  .NOT. land_surface )  THEN
     3301                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3302                                 'res land_surface = .TRUE.'
     3303                CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
     3304             ENDIF
    30853305             IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. cloud_physics )  THEN
    30863306                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     
    30883308                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
    30893309             ENDIF
     3310             IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT. land_surface )  THEN
     3311                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3312                                 'res land_surface = .TRUE.'
     3313                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
     3314             ENDIF
    30903315             IF ( TRIM( var ) == 'pra*'  .AND.  .NOT. precipitation )  THEN
    30913316                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     
    31083333                CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
    31093334             ENDIF
    3110 
     3335             IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT. land_surface )  THEN
     3336                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3337                                 'res land_surface = .TRUE.'
     3338                CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
     3339             ENDIF
     3340             IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )  &
     3341             THEN
     3342                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3343                                 'res land_surface = .TRUE.'
     3344                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
     3345             ENDIF
     3346             IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT. land_surface ) &
     3347             THEN
     3348                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3349                                 'res land_surface = .TRUE.'
     3350                CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
     3351             ENDIF
     3352             IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )  &
     3353             THEN
     3354                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3355                                 'res land_surface = .TRUE.'
     3356                CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
     3357             ENDIF
     3358
     3359             IF ( TRIM( var ) == 'c_liq*' )  unit = 'none'
     3360             IF ( TRIM( var ) == 'c_soil*')  unit = 'none'
     3361             IF ( TRIM( var ) == 'c_veg*' )  unit = 'none'
     3362             IF ( TRIM( var ) == 'ghf_eb*')  unit = 'W/m2'
     3363             IF ( TRIM( var ) == 'lai*'   )  unit = 'none'
    31113364             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/kg*m'
    31123365             IF ( TRIM( var ) == 'pra*'   )  unit = 'mm'
    31133366             IF ( TRIM( var ) == 'prr*'   )  unit = 'mm/s'
    31143367             IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
     3368             IF ( TRIM( var ) == 'qsws_eb*'      ) unit = 'W/m2'
     3369             IF ( TRIM( var ) == 'qsws_liq_eb*'  ) unit = 'W/m2'
     3370             IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
     3371             IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
     3372             IF ( TRIM( var ) == 'rad_net*')       unit = 'W/m2'     
     3373             IF ( TRIM( var ) == 'rad_sw_in*')     unit = 'W/m2'   
    31153374             IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
     3375             IF ( TRIM( var ) == 'shf_eb*')  unit = 'W/m2'
    31163376             IF ( TRIM( var ) == 't*'     )  unit = 'K'
    31173377             IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
     
    31293389
    31303390          CASE DEFAULT
     3391
    31313392             CALL user_check_data_output( var, unit )
    31323393
     
    31373398                   CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
    31383399                ELSE
    3139                    message_string = 'illegal value for data_output =' //       &
     3400                   message_string = 'illegal value for data_output = "' //     &
    31403401                                    TRIM( data_output(i) ) // '"'
    31413402                   CALL message( 'check_parameters', 'PA0115', 1, 2, 0, 6, 0 )
  • palm/trunk/SOURCE/data_output_2d.f90

    r1360 r1551  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Added suppport for land surface model and radiation model output. In the course
     23! of this action, the limits for vertical loops have been changed (from nzb and
     24! nzt+1 to nzb_do and nzt_do, respectively in order to allow soil model output).
     25! Moreover, a new vertical grid zs was introduced.
    2326!
    2427! Former revisions:
     
    127130               
    128131    USE kinds
    129        
     132   
     133    USE land_surface_model_mod,                                                &
     134        ONLY:  c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb,            &
     135               ghf_eb_av, lai, lai_av, m_liq_eb, m_liq_eb_av, m_soil,          &
     136               m_soil_av, nzb_soil, nzt_soil, qsws_eb, qsws_eb_av,             &
     137               qsws_liq_eb, qsws_liq_eb_av, qsws_soil_eb, qsws_soil_eb_av,     &
     138               qsws_veg_eb, qsws_veg_eb_av, shf_eb, shf_eb_av, t_soil,         &
     139               t_soil_av, zs
     140   
    130141    USE netcdf_control
    131142
     
    135146   
    136147    USE pegrid
     148
     149    USE radiation_model_mod,                                                   &
     150        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av
    137151
    138152    IMPLICIT NONE
     
    157171    INTEGER(iwp) ::  n         !:
    158172    INTEGER(iwp) ::  ns        !:
     173    INTEGER(iwp) ::  nzb_do    !: lower limit of the data field (usually nzb)
     174    INTEGER(iwp) ::  nzt_do    !: upper limit of the data field (usually nzt+1)
    159175    INTEGER(iwp) ::  psi       !:
    160176    INTEGER(iwp) ::  s         !:
     
    343359
    344360       IF ( do2d_mode == mode )  THEN
     361
     362          nzb_do = nzb
     363          nzt_do = nzt+1
    345364!
    346365!--       Store the array chosen on the temporary array.
     
    356375                IF ( mode == 'xy' )  level_z = zu
    357376
     377             CASE ( 'c_liq*_xy' )        ! 2d-array
     378                IF ( av == 0 )  THEN
     379                   DO  i = nxlg, nxrg
     380                      DO  j = nysg, nyng
     381                         local_pf(i,j,nzb+1) = c_liq(j,i) * c_veg(j,i)
     382                      ENDDO
     383                   ENDDO
     384                ELSE
     385                   DO  i = nxlg, nxrg
     386                      DO  j = nysg, nyng
     387                         local_pf(i,j,nzb+1) = c_liq_av(j,i)
     388                      ENDDO
     389                   ENDDO
     390                ENDIF
     391                resorted = .TRUE.
     392                two_d = .TRUE.
     393                level_z(nzb+1) = zu(nzb+1)
     394
     395             CASE ( 'c_soil*_xy' )        ! 2d-array
     396                IF ( av == 0 )  THEN
     397                   DO  i = nxlg, nxrg
     398                      DO  j = nysg, nyng
     399                         local_pf(i,j,nzb+1) = 1.0_wp - c_veg(j,i)
     400                      ENDDO
     401                   ENDDO
     402                ELSE
     403                   DO  i = nxlg, nxrg
     404                      DO  j = nysg, nyng
     405                         local_pf(i,j,nzb+1) = c_soil_av(j,i)
     406                      ENDDO
     407                   ENDDO
     408                ENDIF
     409                resorted = .TRUE.
     410                two_d = .TRUE.
     411                level_z(nzb+1) = zu(nzb+1)
     412
     413             CASE ( 'c_veg*_xy' )        ! 2d-array
     414                IF ( av == 0 )  THEN
     415                   DO  i = nxlg, nxrg
     416                      DO  j = nysg, nyng
     417                         local_pf(i,j,nzb+1) = c_veg(j,i)
     418                      ENDDO
     419                   ENDDO
     420                ELSE
     421                   DO  i = nxlg, nxrg
     422                      DO  j = nysg, nyng
     423                         local_pf(i,j,nzb+1) = c_veg_av(j,i)
     424                      ENDDO
     425                   ENDDO
     426                ENDIF
     427                resorted = .TRUE.
     428                two_d = .TRUE.
     429                level_z(nzb+1) = zu(nzb+1)
     430
     431             CASE ( 'ghf_eb*_xy' )        ! 2d-array
     432                IF ( av == 0 )  THEN
     433                   DO  i = nxlg, nxrg
     434                      DO  j = nysg, nyng
     435                         local_pf(i,j,nzb+1) = ghf_eb(j,i)
     436                      ENDDO
     437                   ENDDO
     438                ELSE
     439                   DO  i = nxlg, nxrg
     440                      DO  j = nysg, nyng
     441                         local_pf(i,j,nzb+1) = ghf_eb_av(j,i)
     442                      ENDDO
     443                   ENDDO
     444                ENDIF
     445                resorted = .TRUE.
     446                two_d = .TRUE.
     447                level_z(nzb+1) = zu(nzb+1)
     448
     449             CASE ( 'lai*_xy' )        ! 2d-array
     450                IF ( av == 0 )  THEN
     451                   DO  i = nxlg, nxrg
     452                      DO  j = nysg, nyng
     453                         local_pf(i,j,nzb+1) = lai(j,i)
     454                      ENDDO
     455                   ENDDO
     456                ELSE
     457                   DO  i = nxlg, nxrg
     458                      DO  j = nysg, nyng
     459                         local_pf(i,j,nzb+1) = lai_av(j,i)
     460                      ENDDO
     461                   ENDDO
     462                ENDIF
     463                resorted = .TRUE.
     464                two_d = .TRUE.
     465                level_z(nzb+1) = zu(nzb+1)
     466
    358467             CASE ( 'lpt_xy', 'lpt_xz', 'lpt_yz' )
    359468                IF ( av == 0 )  THEN
     
    382491                two_d = .TRUE.
    383492                level_z(nzb+1) = zu(nzb+1)
     493
     494             CASE ( 'm_liq_eb*_xy' )        ! 2d-array
     495                IF ( av == 0 )  THEN
     496                   DO  i = nxlg, nxrg
     497                      DO  j = nysg, nyng
     498                         local_pf(i,j,nzb+1) = m_liq_eb(j,i)
     499                      ENDDO
     500                   ENDDO
     501                ELSE
     502                   DO  i = nxlg, nxrg
     503                      DO  j = nysg, nyng
     504                         local_pf(i,j,nzb+1) = m_liq_eb_av(j,i)
     505                      ENDDO
     506                   ENDDO
     507                ENDIF
     508                resorted = .TRUE.
     509                two_d = .TRUE.
     510                level_z(nzb+1) = zu(nzb+1)
     511
     512             CASE ( 'm_soil_xy', 'm_soil_xz', 'm_soil_yz' )
     513                nzb_do = nzb_soil
     514                nzt_do = nzt_soil
     515                IF ( av == 0 )  THEN
     516                   to_be_resorted => m_soil
     517                ELSE
     518                   to_be_resorted => m_soil_av
     519                ENDIF
     520                IF ( mode == 'xy' )  level_z = zs
    384521
    385522             CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )
     
    665802                level_z(nzb+1) = zu(nzb+1)
    666803
     804             CASE ( 'qsws_eb*_xy' )        ! 2d-array
     805                IF ( av == 0 ) THEN
     806                   DO  i = nxlg, nxrg
     807                      DO  j = nysg, nyng
     808                         local_pf(i,j,nzb+1) =  qsws_eb(j,i)
     809                      ENDDO
     810                   ENDDO
     811                ELSE
     812                   DO  i = nxlg, nxrg
     813                      DO  j = nysg, nyng
     814                         local_pf(i,j,nzb+1) =  qsws_eb_av(j,i)
     815                      ENDDO
     816                   ENDDO
     817                ENDIF
     818                resorted = .TRUE.
     819                two_d = .TRUE.
     820                level_z(nzb+1) = zu(nzb+1)
     821
     822             CASE ( 'qsws_liq_eb*_xy' )        ! 2d-array
     823                IF ( av == 0 ) THEN
     824                   DO  i = nxlg, nxrg
     825                      DO  j = nysg, nyng
     826                         local_pf(i,j,nzb+1) =  qsws_liq_eb(j,i)
     827                      ENDDO
     828                   ENDDO
     829                ELSE
     830                   DO  i = nxlg, nxrg
     831                      DO  j = nysg, nyng
     832                         local_pf(i,j,nzb+1) =  qsws_liq_eb_av(j,i)
     833                      ENDDO
     834                   ENDDO
     835                ENDIF
     836                resorted = .TRUE.
     837                two_d = .TRUE.
     838                level_z(nzb+1) = zu(nzb+1)
     839
     840             CASE ( 'qsws_soil_eb*_xy' )        ! 2d-array
     841                IF ( av == 0 ) THEN
     842                   DO  i = nxlg, nxrg
     843                      DO  j = nysg, nyng
     844                         local_pf(i,j,nzb+1) =  qsws_soil_eb(j,i)
     845                      ENDDO
     846                   ENDDO
     847                ELSE
     848                   DO  i = nxlg, nxrg
     849                      DO  j = nysg, nyng
     850                         local_pf(i,j,nzb+1) =  qsws_soil_eb_av(j,i)
     851                      ENDDO
     852                   ENDDO
     853                ENDIF
     854                resorted = .TRUE.
     855                two_d = .TRUE.
     856                level_z(nzb+1) = zu(nzb+1)
     857
     858             CASE ( 'qsws_veg_eb*_xy' )        ! 2d-array
     859                IF ( av == 0 ) THEN
     860                   DO  i = nxlg, nxrg
     861                      DO  j = nysg, nyng
     862                         local_pf(i,j,nzb+1) =  qsws_veg_eb(j,i)
     863                      ENDDO
     864                   ENDDO
     865                ELSE
     866                   DO  i = nxlg, nxrg
     867                      DO  j = nysg, nyng
     868                         local_pf(i,j,nzb+1) =  qsws_veg_eb_av(j,i)
     869                      ENDDO
     870                   ENDDO
     871                ENDIF
     872                resorted = .TRUE.
     873                two_d = .TRUE.
     874                level_z(nzb+1) = zu(nzb+1)
     875
    667876             CASE ( 'qv_xy', 'qv_xz', 'qv_yz' )
    668877                IF ( av == 0 )  THEN
     
    680889                IF ( mode == 'xy' )  level_z = zu
    681890
     891             CASE ( 'rad_net*_xy' )        ! 2d-array
     892                IF ( av == 0 ) THEN
     893                   DO  i = nxlg, nxrg
     894                      DO  j = nysg, nyng
     895                         local_pf(i,j,nzb+1) =  rad_net(j,i)
     896                      ENDDO
     897                   ENDDO
     898                ELSE
     899                   DO  i = nxlg, nxrg
     900                      DO  j = nysg, nyng
     901                         local_pf(i,j,nzb+1) =  rad_net_av(j,i)
     902                      ENDDO
     903                   ENDDO
     904                ENDIF
     905                resorted = .TRUE.
     906                two_d = .TRUE.
     907                level_z(nzb+1) = zu(nzb+1)
     908
     909             CASE ( 'rad_sw_in*_xy' )        ! 2d-array
     910                IF ( av == 0 ) THEN
     911                   DO  i = nxlg, nxrg
     912                      DO  j = nysg, nyng
     913                         local_pf(i,j,nzb+1) =  rad_sw_in(j,i)
     914                      ENDDO
     915                   ENDDO
     916                ELSE
     917                   DO  i = nxlg, nxrg
     918                      DO  j = nysg, nyng
     919                         local_pf(i,j,nzb+1) =  rad_sw_in_av(j,i)
     920                      ENDDO
     921                   ENDDO
     922                ENDIF
     923                resorted = .TRUE.
     924                two_d = .TRUE.
     925                level_z(nzb+1) = zu(nzb+1)
     926
    682927             CASE ( 'rho_xy', 'rho_xz', 'rho_yz' )
    683928                IF ( av == 0 )  THEN
     
    719964                level_z(nzb+1) = zu(nzb+1)
    720965
     966             CASE ( 'shf_eb*_xy' )        ! 2d-array
     967                IF ( av == 0 ) THEN
     968                   DO  i = nxlg, nxrg
     969                      DO  j = nysg, nyng
     970                         local_pf(i,j,nzb+1) =  shf_eb(j,i)
     971                      ENDDO
     972                   ENDDO
     973                ELSE
     974                   DO  i = nxlg, nxrg
     975                      DO  j = nysg, nyng
     976                         local_pf(i,j,nzb+1) =  shf_eb_av(j,i)
     977                      ENDDO
     978                   ENDDO
     979                ENDIF
     980                resorted = .TRUE.
     981                two_d = .TRUE.
     982                level_z(nzb+1) = zu(nzb+1)
     983
    721984             CASE ( 't*_xy' )        ! 2d-array
    722985                IF ( av == 0 )  THEN
     
    736999                two_d = .TRUE.
    7371000                level_z(nzb+1) = zu(nzb+1)
     1001
     1002             CASE ( 't_soil_xy', 't_soil_xz', 't_soil_yz' )
     1003                nzb_do = nzb_soil
     1004                nzt_do = nzt_soil
     1005                IF ( av == 0 )  THEN
     1006                   to_be_resorted => t_soil
     1007                ELSE
     1008                   to_be_resorted => t_soil_av
     1009                ENDIF
     1010                IF ( mode == 'xy' )  level_z = zs
    7381011
    7391012             CASE ( 'u_xy', 'u_xz', 'u_yz' )
     
    8391112!--             User defined quantity
    8401113                CALL user_data_output_2d( av, do2d(av,if), found, grid,        &
    841                                           local_pf, two_d )
     1114                                          local_pf, two_d, nzb_do, nzt_do )
    8421115                resorted = .TRUE.
    8431116
     
    8481121                ELSEIF ( grid == 'zu1' ) THEN
    8491122                   IF ( mode == 'xy' )  level_z(nzb+1) = zu(nzb+1)
     1123                ELSEIF ( grid == 'zs' ) THEN
     1124                   IF ( mode == 'xy' )  level_z = zs
    8501125                ENDIF
    8511126
     
    8631138             DO  i = nxlg, nxrg
    8641139                DO  j = nysg, nyng
    865                    DO  k = nzb, nzt+1
     1140                   DO  k = nzb_do, nzt_do
    8661141                      local_pf(i,j,k) = to_be_resorted(k,j,i)
    8671142                   ENDDO
     
    8741149!--       section mode chosen.
    8751150          is = 1
    876    loop1: DO  WHILE ( section(is,s) /= -9999  .OR.  two_d )
     1151   loop1: DO WHILE ( section(is,s) /= -9999  .OR.  two_d )
    8771152
    8781153             SELECT CASE ( mode )
     
    8851160                   ELSE
    8861161                      layer_xy = section(is,s)
     1162                   ENDIF
     1163
     1164!
     1165!--                Exit the loop for layers beyond the data output domain
     1166!--                (used for soil model)
     1167                   IF ( layer_xy .GT. nzt_do )  THEN
     1168                      EXIT loop1
    8871169                   ENDIF
    8881170
     
    9161198!
    9171199!--                   Carry out the averaging (all data are on the PE)
    918                       DO  k = nzb, nzt+1
     1200                      DO  k = nzb_do, nzt_do
    9191201                         DO  j = nysg, nyng
    9201202                            DO  i = nxlg, nxrg
     
    9241206                      ENDDO
    9251207
    926                       local_2d = local_2d / ( nzt -nzb + 2.0_wp)
     1208                      local_2d = local_2d / ( nzt_do - nzb_do + 1.0_wp)
    9271209
    9281210                   ELSE
     
    9671249                         DO  i = 0, io_blocks-1
    9681250                            IF ( i == io_group )  THEN
    969                                WRITE ( 21 )  nxlg, nxrg, nysg, nyng
     1251                               WRITE ( 21 )  nxlg, nxrg, nysg, nyng, nysg, nyng
    9701252                               WRITE ( 21 )  local_2d
    9711253                            ENDIF
     
    11031385                   IF ( section(is,s) == -1 )  THEN
    11041386
    1105                       ALLOCATE( local_2d_l(nxlg:nxrg,nzb:nzt+1) )
     1387                      ALLOCATE( local_2d_l(nxlg:nxrg,nzb_do:nzt_do) )
    11061388                      local_2d_l = 0.0_wp
    1107                       ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
     1389                      ngp = ( nxrg-nxlg + 1 ) * ( nzt_do-nzb_do + 1 )
    11081390!
    11091391!--                   First local averaging on the PE
    1110                       DO  k = nzb, nzt+1
     1392                      DO  k = nzb_do, nzt_do
    11111393                         DO  j = nys, nyn
    11121394                            DO  i = nxlg, nxrg
     
    11201402!--                   Now do the averaging over all PEs along y
    11211403                      IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1122                       CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb),                &
    1123                                           local_2d(nxlg,nzb), ngp, MPI_REAL,   &
     1404                      CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb_do),                &
     1405                                          local_2d(nxlg,nzb_do), ngp, MPI_REAL,   &
    11241406                                          MPI_SUM, comm1dy, ierr )
    11251407#else
     
    11361418                      IF ( section(is,s) >= nys  .AND.  section(is,s) <= nyn ) &
    11371419                      THEN
    1138                          local_2d = local_pf(:,section(is,s),nzb:nzt+1)
     1420                         local_2d = local_pf(:,section(is,s),nzb_do:nzt_do)
    11391421                      ENDIF
    11401422
     
    11571439!--                      output file afterwards to increase the performance.
    11581440                         DO  i = nxlg, nxrg
    1159                             DO  k = nzb, nzt+1
     1441                            DO  k = nzb_do, nzt_do
    11601442                               local_2d_sections_l(i,is,k) = local_2d(i,k)
    11611443                            ENDDO
     
    11841466                                      nys-1 == -1 ) )                          &
    11851467                               THEN
    1186                                   WRITE (22)  nxlg, nxrg, nzb, nzt+1
     1468                                  WRITE (22)  nxlg, nxrg, nzb_do, nzt_do, nzb, nzt+1
    11871469                                  WRITE (22)  local_2d
    11881470                               ELSE
    1189                                   WRITE (22)  -1, -1, -1, -1
     1471                                  WRITE (22)  -1, -1, -1, -1, -1, -1
    11901472                               ENDIF
    11911473                            ENDIF
     
    12031485                         CALL MPI_BARRIER( comm2d, ierr )
    12041486
    1205                          ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
     1487                         ngp = ( nxrg-nxlg + 1 ) * ( nzt_do-nzb_do + 1 )
    12061488                         IF ( myid == 0 )  THEN
    12071489!
     
    12111493                                 ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
    12121494                            THEN
    1213                                total_2d(nxlg:nxrg,nzb:nzt+1) = local_2d
     1495                               total_2d(nxlg:nxrg,nzb_do:nzt_do) = local_2d
    12141496                            ENDIF
    12151497!
     
    12401522!--                         Relocate the local array for the next loop increment
    12411523                            DEALLOCATE( local_2d )
    1242                             ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) )
     1524                            ALLOCATE( local_2d(nxlg:nxrg,nzb_do:nzt_do) )
    12431525
    12441526#if defined( __netcdf )
    12451527                            nc_stat = NF90_PUT_VAR( id_set_xz(av),          &
    12461528                                                 id_var_do2d(av,if),        &
    1247                                                  total_2d(0:nx+1,nzb:nzt+1),&
     1529                                                 total_2d(0:nx+1,nzb_do:nzt_do),&
    12481530                            start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
    1249                                              count = (/ nx+2, 1, nz+2, 1 /) )
     1531                                             count = (/ nx+2, 1, nzt_do-nzb_do+1, 1 /) )
    12501532                            CALL handle_netcdf_error( 'data_output_2d', 58 )
    12511533#endif
     
    12601542                            THEN
    12611543                               ind(1) = nxlg; ind(2) = nxrg
    1262                                ind(3) = nzb;   ind(4) = nzt+1
     1544                               ind(3) = nzb_do;   ind(4) = nzt_do
    12631545                            ELSE
    12641546                               ind(1) = -9999; ind(2) = -9999
     
    12701552!--                         If applicable, send data to PE0.
    12711553                            IF ( ind(1) /= -9999 )  THEN
    1272                                CALL MPI_SEND( local_2d(nxlg,nzb), ngp,         &
     1554                               CALL MPI_SEND( local_2d(nxlg,nzb_do), ngp,         &
    12731555                                              MPI_REAL, 0, 1, comm2d, ierr )
    12741556                            ENDIF
     
    12861568                   nc_stat = NF90_PUT_VAR( id_set_xz(av),                   &
    12871569                                           id_var_do2d(av,if),              &
    1288                                            local_2d(nxl:nxr+1,nzb:nzt+1),   &
     1570                                           local_2d(nxl:nxr+1,nzb_do:nzt_do),   &
    12891571                            start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
    1290                                            count = (/ nx+2, 1, nz+2, 1 /) )
     1572                                           count = (/ nx+2, 1, nzt_do-nzb_do+1, 1 /) )
    12911573                   CALL handle_netcdf_error( 'data_output_2d', 451 )
    12921574#endif
     
    13221604                   IF ( section(is,s) == -1 )  THEN
    13231605
    1324                       ALLOCATE( local_2d_l(nysg:nyng,nzb:nzt+1) )
     1606                      ALLOCATE( local_2d_l(nysg:nyng,nzb_do:nzt_do) )
    13251607                      local_2d_l = 0.0_wp
    1326                       ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 )
     1608                      ngp = ( nyng-nysg+1 ) * ( nzt_do-nzb_do+1 )
    13271609!
    13281610!--                   First local averaging on the PE
    1329                       DO  k = nzb, nzt+1
     1611                      DO  k = nzb_do, nzt_do
    13301612                         DO  j = nysg, nyng
    13311613                            DO  i = nxl, nxr
     
    13391621!--                   Now do the averaging over all PEs along x
    13401622                      IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1341                       CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb),                &
    1342                                           local_2d(nysg,nzb), ngp, MPI_REAL,   &
     1623                      CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb_do),                &
     1624                                          local_2d(nysg,nzb_do), ngp, MPI_REAL,   &
    13431625                                          MPI_SUM, comm1dx, ierr )
    13441626#else
     
    13551637                      IF ( section(is,s) >= nxl  .AND.  section(is,s) <= nxr ) &
    13561638                      THEN
    1357                          local_2d = local_pf(section(is,s),:,nzb:nzt+1)
     1639                         local_2d = local_pf(section(is,s),:,nzb_do:nzt_do)
    13581640                      ENDIF
    13591641
     
    13761658!--                      output file afterwards to increase the performance.
    13771659                         DO  j = nysg, nyng
    1378                             DO  k = nzb, nzt+1
     1660                            DO  k = nzb_do, nzt_do
    13791661                               local_2d_sections_l(is,j,k) = local_2d(j,k)
    13801662                            ENDDO
     
    14031685                                      nxl-1 == -1 ) )                          &
    14041686                               THEN
    1405                                   WRITE (23)  nysg, nyng, nzb, nzt+1
     1687                                  WRITE (23)  nysg, nyng, nzb_do, nzt_do, nzb, nzt+1
    14061688                                  WRITE (23)  local_2d
    14071689                               ELSE
    1408                                   WRITE (23)  -1, -1, -1, -1
     1690                                  WRITE (23)  -1, -1, -1, -1, -1, -1
    14091691                               ENDIF
    14101692                            ENDIF
     
    14221704                         CALL MPI_BARRIER( comm2d, ierr )
    14231705
    1424                          ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 )
     1706                         ngp = ( nyng-nysg+1 ) * ( nzt_do-nzb_do+1 )
    14251707                         IF ( myid == 0 )  THEN
    14261708!
     
    14301712                                 ( section(is,s) == -1  .AND.  nxl-1 == -1 ) ) &
    14311713                            THEN
    1432                                total_2d(nysg:nyng,nzb:nzt+1) = local_2d
     1714                               total_2d(nysg:nyng,nzb_do:nzt_do) = local_2d
    14331715                            ENDIF
    14341716!
     
    14591741!--                         Relocate the local array for the next loop increment
    14601742                            DEALLOCATE( local_2d )
    1461                             ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) )
     1743                            ALLOCATE( local_2d(nysg:nyng,nzb_do:nzt_do) )
    14621744
    14631745#if defined( __netcdf )
    14641746                            nc_stat = NF90_PUT_VAR( id_set_yz(av),          &
    14651747                                                 id_var_do2d(av,if),        &
    1466                                                  total_2d(0:ny+1,nzb:nzt+1),&
     1748                                                 total_2d(0:ny+1,nzb_do:nzt_do),&
    14671749                            start = (/ is, 1, 1, do2d_yz_time_count(av) /), &
    1468                                              count = (/ 1, ny+2, nz+2, 1 /) )
     1750                                             count = (/ 1, ny+2, nzt_do-nzb_do+1, 1 /) )
    14691751                            CALL handle_netcdf_error( 'data_output_2d', 61 )
    14701752#endif
     
    14791761                            THEN
    14801762                               ind(1) = nysg; ind(2) = nyng
    1481                                ind(3) = nzb;   ind(4) = nzt+1
     1763                               ind(3) = nzb_do;   ind(4) = nzt_do
    14821764                            ELSE
    14831765                               ind(1) = -9999; ind(2) = -9999
     
    14891771!--                         If applicable, send data to PE0.
    14901772                            IF ( ind(1) /= -9999 )  THEN
    1491                                CALL MPI_SEND( local_2d(nysg,nzb), ngp,         &
     1773                               CALL MPI_SEND( local_2d(nysg,nzb_do), ngp,         &
    14921774                                              MPI_REAL, 0, 1, comm2d, ierr )
    14931775                            ENDIF
     
    15051787                   nc_stat = NF90_PUT_VAR( id_set_yz(av),                   &
    15061788                                           id_var_do2d(av,if),              &
    1507                                            local_2d(nys:nyn+1,nzb:nzt+1),   &
     1789                                           local_2d(nys:nyn+1,nzb_do:nzt_do),   &
    15081790                            start = (/ is, 1, 1, do2d_xz_time_count(av) /), &
    1509                                            count = (/ 1, ny+2, nz+2, 1 /) )
     1791                                           count = (/ 1, ny+2, nzt_do-nzb_do+1, 1 /) )
    15101792                   CALL handle_netcdf_error( 'data_output_2d', 452 )
    15111793#endif
     
    15951877!
    15961878!--                      Distribute data over all PEs along y
    1597                          ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 ) * ns
     1879                         ngp = ( nxrg-nxlg+1 ) * ( nzt_do-nzb_do+1 ) * ns
    15981880                         IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
    1599                          CALL MPI_ALLREDUCE( local_2d_sections_l(nxlg,1,nzb),  &
    1600                                              local_2d_sections(nxlg,1,nzb),    &
     1881                         CALL MPI_ALLREDUCE( local_2d_sections_l(nxlg,1,nzb_do),  &
     1882                                             local_2d_sections(nxlg,1,nzb_do),    &
    16011883                                             ngp, MPI_REAL, MPI_SUM, comm1dy,  &
    16021884                                             ierr )
     
    16121894                                             id_var_do2d(av,if),               &
    16131895                                             local_2d_sections(nxl:nxr+1,1:ns, &
    1614                                                 nzb:nzt+1),                    &
     1896                                                nzb_do:nzt_do),                &
    16151897                                             start = (/ nxl+1, 1, 1,           &
    16161898                                                do2d_xz_time_count(av) /),     &
    1617                                              count = (/ nxr-nxl+2, ns, nzt+2,  &
     1899                                             count = (/ nxr-nxl+2, ns, nzt_do-nzb_do+1,  &
    16181900                                                        1 /) )
    16191901                      ELSE
     
    16211903                                             id_var_do2d(av,if),               &
    16221904                                             local_2d_sections(nxl:nxr,1:ns,   &
    1623                                                 nzb:nzt+1),                    &
     1905                                                nzb_do:nzt_do),                &
    16241906                                             start = (/ nxl+1, 1, 1,           &
    16251907                                                do2d_xz_time_count(av) /),     &
    1626                                              count = (/ nxr-nxl+1, ns, nzt+2,  &
     1908                                             count = (/ nxr-nxl+1, ns, nzt_do-nzb_do+1,  &
    16271909                                                1 /) )
    16281910                      ENDIF
     
    16471929                         ngp = ( nyng-nysg+1 ) * ( nzt-nzb + 2 ) * ns
    16481930                         IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
    1649                          CALL MPI_ALLREDUCE( local_2d_sections_l(1,nysg,nzb),  &
    1650                                              local_2d_sections(1,nysg,nzb),    &
     1931                         CALL MPI_ALLREDUCE( local_2d_sections_l(1,nysg,nzb_do),  &
     1932                                             local_2d_sections(1,nysg,nzb_do),    &
    16511933                                             ngp, MPI_REAL, MPI_SUM, comm1dx,  &
    16521934                                             ierr )
     
    16621944                                             id_var_do2d(av,if),               &
    16631945                                             local_2d_sections(1:ns,           &
    1664                                                 nys:nyn+1,nzb:nzt+1),          &
     1946                                                nys:nyn+1,nzb_do:nzt_do),      &
    16651947                                             start = (/ 1, nys+1, 1,           &
    16661948                                                do2d_yz_time_count(av) /),     &
    16671949                                             count = (/ ns, nyn-nys+2,         &
    1668                                                         nzt+2, 1 /) )
     1950                                                        nzt_do-nzb_do+1, 1 /) )
    16691951                      ELSE
    16701952                         nc_stat = NF90_PUT_VAR( id_set_yz(av),                &
    16711953                                             id_var_do2d(av,if),               &
    16721954                                             local_2d_sections(1:ns,nys:nyn,   &
    1673                                                 nzb:nzt+1),                    &
     1955                                                nzb_do:nzt_do),                &
    16741956                                             start = (/ 1, nys+1, 1,           &
    16751957                                                do2d_yz_time_count(av) /),     &
    16761958                                             count = (/ ns, nyn-nys+1,         &
    1677                                                         nzt+2, 1 /) )
     1959                                                        nzt_do-nzb_do+1, 1 /) )
    16781960                      ENDIF
    16791961
  • palm/trunk/SOURCE/data_output_3d.f90

    r1360 r1551  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Added suppport for land surface model and radiation model output. In the course
     23! of this action, the limits for vertical loops have been changed (from nzb and
     24! nzt+1 to nzb_do and nzt_do, respectively in order to allow soil model output).
     25! Moreover, a new vertical grid zs was introduced.
    2326!
    2427! Former revisions:
     
    112115    USE kinds
    113116   
     117    USE land_surface_model_mod,                                                &
     118        ONLY: m_soil, m_soil_av, nzb_soil, nzt_soil, t_soil, t_soil_av
     119
    114120    USE netcdf_control
    115121       
     
    130136    INTEGER(iwp) ::  k         !:
    131137    INTEGER(iwp) ::  n         !:
     138    INTEGER(iwp) ::  nzb_do    !: vertical lower limit for data output
     139    INTEGER(iwp) ::  nzt_do    !: vertical upper limit for data output
    132140    INTEGER(iwp) ::  pos       !:
    133141    INTEGER(iwp) ::  prec      !:
     
    183191
    184192!
    185 !-- Allocate a temporary array with the desired output dimensions.
    186     ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb:nz_do3d) )
    187 
    188 !
    189193!-- Update the netCDF time axis
    190194!-- In case of parallel output, this is only done by PE0 to increase the
     
    209213!--    Store the array chosen on the temporary array.
    210214       resorted = .FALSE.
     215       nzb_do = nzb
     216       nzt_do = nz_do3d
     217
     218!
     219!--    Allocate a temporary array with the desired output dimensions.
     220       ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) )
     221
    211222       SELECT CASE ( TRIM( do3d(av,if) ) )
    212223
     
    223234             ELSE
    224235                to_be_resorted => lpt_av
     236             ENDIF
     237
     238          CASE ( 'm_soil' )
     239             nzb_do = nzb_soil
     240             nzt_do = nzt_soil
     241!
     242!--          For soil model quantities, it is required to re-allocate local_pf
     243             DEALLOCATE ( local_pf )
     244             ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) )
     245
     246             IF ( av == 0 )  THEN
     247                to_be_resorted => m_soil
     248             ELSE
     249                to_be_resorted => m_soil_av
    225250             ENDIF
    226251
     
    251276                DO  i = nxlg, nxrg
    252277                   DO  j = nysg, nyng
    253                       DO  k = nzb, nz_do3d
     278                      DO  k = nzb_do, nzt_do
    254279                         local_pf(i,j,k) = tend(k,j,i)
    255280                      ENDDO
     
    267292                   DO  i = nxl, nxr
    268293                      DO  j = nys, nyn
    269                          DO  k = nzb, nz_do3d
     294                         DO  k = nzb_do, nzt_do
    270295                            number_of_particles = prt_count(k,j,i)
    271296                            IF (number_of_particles <= 0)  CYCLE
     
    296321                DO  i = nxlg, nxrg
    297322                   DO  j = nysg, nyng
    298                       DO  k = nzb, nz_do3d
     323                      DO  k = nzb_do, nzt_do
    299324                         local_pf(i,j,k) = tend(k,j,i)
    300325                      ENDDO
     
    336361                   DO  i = nxlg, nxrg
    337362                      DO  j = nysg, nyng
    338                          DO  k = nzb, nz_do3d
     363                         DO  k = nzb_do, nzt_do
    339364                            local_pf(i,j,k) = pt(k,j,i) + l_d_cp *             &
    340365                                                          pt_d_t(k) *          &
     
    389414                   DO  i = nxl, nxr
    390415                      DO  j = nys, nyn
    391                          DO  k = nzb, nz_do3d
     416                         DO  k = nzb_do, nzt_do
    392417                            number_of_particles = prt_count(k,j,i)
    393418                            IF (number_of_particles <= 0)  CYCLE
     
    409434                DO  i = nxlg, nxrg
    410435                   DO  j = nysg, nyng
    411                       DO  k = nzb, nz_do3d
     436                      DO  k = nzb_do, nzt_do
    412437                         local_pf(i,j,k) = tend(k,j,i)
    413438                      ENDDO
     
    431456                DO  i = nxlg, nxrg
    432457                   DO  j = nysg, nyng
    433                       DO  k = nzb, nz_do3d
     458                      DO  k = nzb_do, nzt_do
    434459                         local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
    435460                      ENDDO
     
    462487             ENDIF
    463488
     489          CASE ( 't_soil' )
     490             nzb_do = nzb_soil
     491             nzt_do = nzt_soil
     492!
     493!--          For soil model quantities, it is required to re-allocate local_pf
     494             DEALLOCATE ( local_pf )
     495             ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) )
     496
     497             IF ( av == 0 )  THEN
     498                to_be_resorted => t_soil
     499             ELSE
     500                to_be_resorted => t_soil_av
     501             ENDIF
     502
    464503          CASE ( 'u' )
    465504             IF ( av == 0 )  THEN
     
    494533!--          User defined quantity
    495534             CALL user_data_output_3d( av, do3d(av,if), found, local_pf,       &
    496                                        nz_do3d )
     535                                       nzb_do, nzt_do )
    497536             resorted = .TRUE.
    498537
     
    510549          DO  i = nxlg, nxrg
    511550             DO  j = nysg, nyng
    512                 DO  k = nzb, nz_do3d
     551                DO  k = nzb_do, nzt_do
    513552                   local_pf(i,j,k) = to_be_resorted(k,j,i)
    514553                ENDDO
     
    531570          DO  i = 0, io_blocks-1
    532571             IF ( i == io_group )  THEN
    533                 WRITE ( 30 )  nxlg, nxrg, nysg, nyng, nzb, nz_do3d
    534                 WRITE ( 30 )  local_pf
     572                WRITE ( 30 )  nxlg, nxrg, nysg, nyng, nzb_do, nzt_do
     573                WRITE ( 30 )  local_pf(:,:,nzb_do:nzt_do)
    535574             ENDIF
    536575#if defined( __parallel )
     
    547586          IF ( nxr == nx  .AND.  nyn /= ny )  THEN
    548587             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
    549                                local_pf(nxl:nxr+1,nys:nyn,nzb:nz_do3d),  &
    550                 start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /),  &
    551                 count = (/ nxr-nxl+2, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
     588                               local_pf(nxl:nxr+1,nys:nyn,nzb_do:nzt_do),    &
     589                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
     590                count = (/ nxr-nxl+2, nyn-nys+1, nzt_do-nzb_do+1, 1 /) )
    552591          ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
    553592             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
    554                                local_pf(nxl:nxr,nys:nyn+1,nzb:nz_do3d),  &
    555                 start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /),  &
    556                 count = (/ nxr-nxl+1, nyn-nys+2, nz_do3d-nzb+1, 1 /) )
     593                               local_pf(nxl:nxr,nys:nyn+1,nzb_do:nzt_do),    &
     594                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
     595                count = (/ nxr-nxl+1, nyn-nys+2, nzt_do-nzb_do+1, 1 /) )
    557596          ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
    558597             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
    559                              local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d),  &
    560                 start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /),  &
    561                 count = (/ nxr-nxl+2, nyn-nys+2, nz_do3d-nzb+1, 1 /) )
     598                             local_pf(nxl:nxr+1,nys:nyn+1,nzb_do:nzt_do  ),  &
     599                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
     600                count = (/ nxr-nxl+2, nyn-nys+2, nzt_do-nzb_do+1, 1 /) )
    562601          ELSE
    563602             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
    564                                  local_pf(nxl:nxr,nys:nyn,nzb:nz_do3d),  &
    565                 start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /),  &
    566                 count = (/ nxr-nxl+1, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
     603                                 local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do),    &
     604                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
     605                count = (/ nxr-nxl+1, nyn-nys+1, nzt_do-nzb_do+1, 1 /) )
    567606          ENDIF
    568607          CALL handle_netcdf_error( 'data_output_3d', 386 )
     
    572611#if defined( __netcdf )
    573612       nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),        &
    574                          local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d),      &
     613                         local_pf(nxl:nxr+1,nys:nyn+1,nzb_do:nzt_do),        &
    575614                         start = (/ 1, 1, 1, do3d_time_count(av) /),     &
    576                          count = (/ nx+2, ny+2, nz_do3d-nzb+1, 1 /) )
     615                         count = (/ nx+2, ny+2, nzt_do-nzb_do+1, 1 /) )
    577616       CALL handle_netcdf_error( 'data_output_3d', 446 )
    578617#endif
     
    581620       if = if + 1
    582621
     622!
     623!--    Deallocate temporary array
     624       DEALLOCATE ( local_pf )
     625
    583626    ENDDO
    584 
    585 !
    586 !-- Deallocate temporary array.
    587     DEALLOCATE( local_pf )
    588 
    589627
    590628    CALL cpu_log( log_point(14), 'data_output_3d', 'stop' )
  • palm/trunk/SOURCE/flow_statistics.f90

    r1499 r1551  
    2121! Current revisions:
    2222! -----------------
    23 !
     23! Added suppport for land surface model and radiation model output.
    2424!
    2525! Former revisions:
     
    132132       
    133133    USE cloud_parameters,                                                      &
    134         ONLY :  l_d_cp, prr, pt_d_t
     134        ONLY:   l_d_cp, prr, pt_d_t
    135135       
    136136    USE control_parameters,                                                    &
    137         ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
     137        ONLY:   average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
    138138                dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
    139139                large_scale_subsidence, max_pr_user, message_string, ocean,    &
     
    143143       
    144144    USE cpulog,                                                                &
    145         ONLY :  cpu_log, log_point
     145        ONLY:   cpu_log, log_point
    146146       
    147147    USE grid_variables,                                                        &
    148         ONLY :  ddx, ddy
     148        ONLY:   ddx, ddy
    149149       
    150150    USE indices,                                                               &
    151         ONLY :  ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
     151        ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
    152152                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner,        &
    153153                nzb_s_inner, nzt, nzt_diff
     
    155155    USE kinds
    156156   
     157    USE land_surface_model_mod,                                                &
     158        ONLY:   dots_soil, ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil,   &
     159                qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, shf_eb,       &
     160                t_soil
     161
    157162    USE pegrid
     163
     164    USE radiation_model_mod,                                                   &
     165        ONLY :  dots_rad, rad_net, rad_sw_in, radiation
    158166   
    159167    USE statistics
     
    703711                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
    704712                ENDIF
     713             ENDIF
     714
     715             IF ( land_surface )  THEN
     716                sums_l(nzb,93,tn) = sums_l(nzb,93,tn) + ghf_eb(j,i)
     717                sums_l(nzb,94,tn) = sums_l(nzb,94,tn) + shf_eb(j,i)
     718                sums_l(nzb,95,tn) = sums_l(nzb,95,tn) + qsws_eb(j,i)
     719                sums_l(nzb,96,tn) = sums_l(nzb,96,tn) + qsws_liq_eb(j,i)
     720                sums_l(nzb,97,tn) = sums_l(nzb,97,tn) + qsws_soil_eb(j,i)
     721                sums_l(nzb,98,tn) = sums_l(nzb,98,tn) + qsws_veg_eb(j,i)
     722             ENDIF
     723
     724             IF ( radiation )  THEN
     725                sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)  + rad_net(j,i)
     726                sums_l(nzb,100,tn) = sums_l(nzb,100,tn) + rad_sw_in(j,i)
    705727             ENDIF
    706728
     
    10741096       ENDIF
    10751097
     1098
     1099       IF ( land_surface )  THEN
     1100          !$OMP DO
     1101          DO  i = nxl, nxr
     1102             DO  j =  nys, nyn
     1103                DO  k = nzb_soil, nzt_soil
     1104                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i) * rmask(j,i,sr)
     1105                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i) * rmask(j,i,sr)
     1106                ENDDO
     1107             ENDDO
     1108          ENDDO
     1109       ENDIF
     1110       
     1111
    10761112!
    10771113!--    Calculate the user-defined profiles
     
    11331169          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
    11341170          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
    1135           sums(k,89:pr_palm-2)    = sums(k,89:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
     1171          sums(k,89:100)           = sums(k,89:100)     / ngp_2dh(sr)
     1172          sums(k,101:pr_palm-2)    = sums(k,101:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
    11361173       ENDDO
    11371174
     
    12481285       ENDIF
    12491286
     1287       IF ( land_surface )  THEN
     1288          hom(:,1,89,sr) = sums(:,89)              ! t_soil
     1289                                                   ! 90 is initial t_soil profile
     1290          hom(:,1,91,sr) = sums(:,91)              ! m_soil
     1291                                                   ! 92 is initial m_soil profile
     1292          hom(:,1,93,sr) = sums(:,93)              ! ghf_eb
     1293          hom(:,1,94,sr) = sums(:,94)              ! shf_eb
     1294          hom(:,1,95,sr) = sums(:,95)              ! qsws_eb
     1295          hom(:,1,96,sr) = sums(:,96)              ! qsws_liq_eb
     1296          hom(:,1,97,sr) = sums(:,97)              ! qsws_soil_eb
     1297          hom(:,1,98,sr) = sums(:,98)              ! qsws_veg_eb
     1298       ENDIF
     1299
     1300       IF ( radiation )  THEN
     1301          hom(:,1,99 ,sr) = sums(:,99)             ! rad_net
     1302          hom(:,1,100,sr) = sums(:,100)            ! rad_sw_in
     1303       ENDIF
     1304
    12501305       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
    12511306                                       ! upstream-parts u_x, u_y, u_z, v_x,
     
    13921447
    13931448       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
     1449
     1450!
     1451!--    Collect land surface model timeseries
     1452       IF ( land_surface )  THEN
     1453          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf_eb
     1454          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! shf_eb
     1455          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_eb
     1456          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_liq_eb
     1457          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! qsws_soil_eb
     1458          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! qsws_veg_eb
     1459       ENDIF
     1460!
     1461!--    Collect radiation model timeseries
     1462       IF ( radiation )  THEN
     1463          ts_value(dots_rad,sr)   = hom(nzb,1,99,sr)           ! rad_net
     1464          ts_value(dots_rad+1,sr) = hom(nzb,1,100,sr)          ! rad_sw_in
     1465       ENDIF
     1466
    13941467!
    13951468!--    Calculate additional statistics provided by the user interface
     
    28422915       ENDIF
    28432916
     2917
     2918       IF ( land_surface )  THEN
     2919          !$OMP DO
     2920          DO  i = nxl, nxr
     2921             DO  j =  nys, nyn
     2922                DO  k = nzb_soil, nzt_soil
     2923                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i) * rmask(j,i,sr)
     2924                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i) * rmask(j,i,sr)
     2925                ENDDO
     2926             ENDDO
     2927          ENDDO
     2928       ENDIF
     2929
     2930
    28442931!
    28452932!--    Calculate the user-defined profiles
     
    29032990          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
    29042991          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
    2905           sums(k,89:pr_palm-2)    = sums(k,89:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
     2992          sums(k,89:100)           = sums(k,89:100)     / ngp_2dh(sr)
     2993          sums(k,101:pr_palm-2)    = sums(k,101:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
    29062994       ENDDO
    29072995
     
    31623250
    31633251       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
     3252
     3253!
     3254!--    Collect land surface model timeseries
     3255       IF ( land_surface )  THEN
     3256          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf_eb
     3257          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! shf_eb
     3258          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_eb
     3259          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_liq_eb
     3260          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! qsws_soil_eb
     3261          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! qsws_veg_eb
     3262       ENDIF
     3263!
     3264!--    Collect radiation model timeseries
     3265       IF ( radiation )  THEN
     3266          ts_value(dots_rad,sr)   = hom(nzb,1,99,sr)           ! rad_net
     3267          ts_value(dots_rad+1,sr) = hom(nzb,1,100,sr)          ! rad_sw_in
     3268       ENDIF
     3269
    31643270!
    31653271!--    Calculate additional statistics provided by the user interface
  • palm/trunk/SOURCE/header.f90

    r1497 r1551  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Added informal output for land surface model and radiation model. Removed typo.
    2323!
    2424! Former revisions:
     
    170170! Description:
    171171! ------------
    172 ! Writing a header with all important informations about the actual run.
     172! Writing a header with all important information about the actual run.
    173173! This subroutine is called three times, two times at the beginning
    174174! (writing information on files RUN_CONTROL and HEADER) and one time at the
     
    200200       
    201201    USE kinds
    202    
     202   
     203    USE land_surface_model_mod,                                                &
     204        ONLY:  conserve_water_content, dewfall, land_surface, nzb_soil,        &
     205               nzt_soil, root_fraction, soil_moisture, soil_temperature,       &
     206               soil_type, soil_type_name, veg_type, veg_type_name, zs
     207 
    203208    USE model_1d,                                                              &
    204209        ONLY:  damp_level_ind_1d, dt_pr_1d, dt_run_control_1d, end_time_1d
     
    225230               lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, pch_index, &
    226231               plant_canopy
     232
     233    USE radiation_model_mod,                                                   &
     234        ONLY:  albedo, day_init, dt_radiation, lambda, net_radiation,          &
     235               radiation, radiation_scheme, time_utc_init
    227236   
    228237    USE spectrum,                                                              &
     
    263272    CHARACTER (LEN=86) ::  gradients           !:
    264273    CHARACTER (LEN=86) ::  leaf_area_density   !:
     274    CHARACTER (LEN=86) ::  roots               !:
    265275    CHARACTER (LEN=86) ::  slices              !:
    266276    CHARACTER (LEN=86) ::  temperatures        !:
     
    311321!
    312322!-- At the end of the run, output file (HEADER) will be rewritten with
    313 !-- new informations
     323!-- new information
    314324    IF ( io == 19 .AND. simulated_time_at_begin /= simulated_time ) REWIND( 19 )
    315325
     
    495505
    496506!
    497 !-- Runtime and timestep informations
     507!-- Runtime and timestep information
    498508    WRITE ( io, 200 )
    499509    IF ( .NOT. dt_fixed )  THEN
     
    850860
    851861
     862    IF ( land_surface )  THEN
     863
     864       temperatures = ''
     865       gradients    = '' ! use for humidity here
     866       coordinates  = '' ! use for height
     867       roots        = '' ! use for root fraction
     868       slices       = '' ! use for index
     869
     870       i = 1
     871       DO i = nzb_soil, nzt_soil
     872          WRITE (coor_chr,'(F10.2,7X)') soil_temperature(i)
     873          temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
     874
     875          WRITE (coor_chr,'(F10.2,7X)') soil_moisture(i)
     876          gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
     877
     878          WRITE (coor_chr,'(F10.2,7X)')  - zs(i)
     879          coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
     880
     881          WRITE (coor_chr,'(F10.2,7X)')  root_fraction(i)
     882          roots = TRIM( roots ) // ' '  // TRIM( coor_chr )
     883
     884          WRITE (coor_chr,'(I10,7X)')  i
     885          slices = TRIM( slices ) // ' '  // TRIM( coor_chr )
     886
     887
     888       ENDDO
     889
     890!
     891!--    Write land surface model header
     892       WRITE( io, 419 )
     893       IF ( conserve_water_content )  THEN
     894          WRITE( io, 440 )
     895       ELSE
     896          WRITE( io, 441 )
     897       ENDIF
     898
     899       IF ( dewfall )  THEN
     900          WRITE( io, 442 )
     901       ELSE
     902          WRITE( io, 443 )
     903       ENDIF
     904
     905       WRITE( io, 438 ) veg_type_name(veg_type), soil_type_name(soil_type)
     906       WRITE( io, 439 ) TRIM( coordinates ), TRIM( temperatures ),             &
     907                        TRIM( gradients ), TRIM( roots ), TRIM( slices )
     908
     909
     910    ENDIF
     911
     912    IF ( radiation )  THEN
     913!
     914!--    Write land surface model header
     915       WRITE( io, 444 )
     916
     917       IF ( radiation_scheme == "constant" )  THEN
     918          WRITE( io, 445 ) net_radiation
     919       ELSEIF ( radiation_scheme == "clear-sky" )  THEN
     920          WRITE( io, 446 )
     921       ELSE
     922          WRITE( io, 447 ) radiation_scheme
     923       ENDIF
     924
     925       WRITE( io, 448 ) albedo
     926       WRITE( io, 449 ) dt_radiation
     927
     928    ENDIF
     929
     930
    852931!
    853932!-- Boundary conditions
     
    877956
    878957    IF ( ibc_pt_b == 0 )  THEN
    879        runten = TRIM( runten ) // ' pt(0)   = pt_surface'
     958       IF ( land_surface )  THEN
     959          runten = TRIM( runten ) // ' pt(0)     = from soil model'
     960       ELSE
     961          runten = TRIM( runten ) // ' pt(0)     = pt_surface'
     962       ENDIF
    880963    ELSEIF ( ibc_pt_b == 1 )  THEN
    881        runten = TRIM( runten ) // ' pt(0)   = pt(1)'
     964       runten = TRIM( runten ) // ' pt(0)     = pt(1)'
    882965    ELSEIF ( ibc_pt_b == 2 )  THEN
    883        runten = TRIM( runten ) // ' pt(0) = from coupled model'
     966       runten = TRIM( runten ) // ' pt(0)     = from coupled model'
    884967    ENDIF
    885968    IF ( ibc_pt_t == 0 )  THEN
     
    9181001    IF ( humidity )  THEN
    9191002       IF ( ibc_q_b == 0 )  THEN
    920           runten = 'q(0)     = q_surface'
     1003          IF ( land_surface )  THEN
     1004             runten = 'q(0)     = from soil model'
     1005          ELSE
     1006             runten = 'q(0)     = q_surface'
     1007          ENDIF
     1008
    9211009       ELSE
    9221010          runten = 'q(0)     = q(1)'
     
    12251313             coordinates = '/'
    12261314!
    1227 !--          Building strings with index and coordinate informations of the
     1315!--          Building strings with index and coordinate information of the
    12281316!--          slices
    12291317             DO  WHILE ( section(i,1) /= -9999 )
     
    12711359             coordinates = '/'
    12721360!
    1273 !--          Building strings with index and coordinate informations of the
     1361!--          Building strings with index and coordinate information of the
    12741362!--          slices
    12751363             DO  WHILE ( section(i,2) /= -9999 )
     
    13131401             coordinates = '/'
    13141402!
    1315 !--          Building strings with index and coordinate informations of the
     1403!--          Building strings with index and coordinate information of the
    13161404!--          slices
    13171405             DO  WHILE ( section(i,3) /= -9999 )
     
    15711659!
    15721660!-- Geostrophic parameters
    1573     WRITE ( io, 410 )  omega, phi, f, fs
     1661    IF ( radiation .AND. radiation_scheme /= 'constant' )  THEN
     1662       WRITE ( io, 417 )  lambda
     1663    ENDIF
     1664    WRITE ( io, 410 )  phi, omega, f, fs
    15741665
    15751666!
    15761667!-- Other quantities
    15771668    WRITE ( io, 411 )  g
     1669    IF ( radiation .AND. radiation_scheme /= 'constant' )  THEN
     1670       WRITE ( io, 418 )  day_init, time_utc_init
     1671    ENDIF
     1672
    15781673    WRITE ( io, 412 )  TRIM( reference_state )
    15791674    IF ( use_single_reference_value )  THEN
     
    17321827
    17331828!
    1734 !-- User-defined informations
     1829!-- User-defined information
    17351830    CALL user_header( io )
    17361831
     
    18671962260 FORMAT (/' The model has a slope in x-direction. Inclination angle: ',F6.2,&
    18681963             ' degrees')
    1869 270 FORMAT (//' Topography informations:'/ &
    1870               ' -----------------------'// &
     1964270 FORMAT (//' Topography information:'/ &
     1965              ' ----------------------'// &
    18711966              1X,'Topography: ',A)
    18721967271 FORMAT (  ' Building size (x/y/z) in m: ',F5.1,' / ',F5.1,' / ',F5.1/ &
     
    19052000             ' -------------------'// &
    19062001             '                     p                    uv             ', &
    1907              '                   pt'// &
     2002             '                     pt'// &
    19082003             ' B. bound.: ',A/ &
    19092004             ' T. bound.: ',A)
     
    20472142400 FORMAT (//' Physical quantities:'/ &
    20482143              ' -------------------'/)
    2049 410 FORMAT ('    Angular velocity    :   omega = ',E9.3,' rad/s'/  &
    2050             '    Geograph. latitude  :   phi   = ',F4.1,' degr'/   &
    2051             '    Coriolis parameter  :   f     = ',F9.6,' 1/s'/    &
    2052             '                            f*    = ',F9.6,' 1/s')
    2053 411 FORMAT (/'    Gravity             :   g     = ',F4.1,' m/s**2')
     2144410 FORMAT ('    Geograph. latitude  :   phi    = ',F4.1,' degr'/   &
     2145            '    Angular velocity    :   omega  = ',E9.3,' rad/s'/  &
     2146            '    Coriolis parameter  :   f      = ',F9.6,' 1/s'/    &
     2147            '                            f*     = ',F9.6,' 1/s')
     2148411 FORMAT (/'    Gravity             :   g      = ',F4.1,' m/s**2')
    20542149412 FORMAT (/'    Reference state used in buoyancy terms: ',A)
    20552150413 FORMAT ('       Reference density in buoyancy terms: ',F8.3,' kg/m**3')
    20562151414 FORMAT ('       Reference temperature in buoyancy terms: ',F8.4,' K')
    2057 415 FORMAT (/'    Cloud physics parameters:'/ &
    2058              '    ------------------------'/)
    2059 416 FORMAT ('        Surface pressure   :   p_0   = ',F7.2,' hPa'/      &
    2060             '        Gas constant       :   R     = ',F5.1,' J/(kg K)'/ &
    2061             '        Density of air     :   rho_0 = ',F5.3,' kg/m**3'/  &
    2062             '        Specific heat cap. :   c_p   = ',F6.1,' J/(kg K)'/ &
    2063             '        Vapourization heat :   L_v   = ',E8.2,' J/kg')
     2152415 FORMAT (/' Cloud physics parameters:'/ &
     2153             ' ------------------------'/)
     2154416 FORMAT ('    Surface pressure   :   p_0   = ',F7.2,' hPa'/      &
     2155            '    Gas constant       :   R     = ',F5.1,' J/(kg K)'/ &
     2156            '    Density of air     :   rho_0 = ',F5.3,' kg/m**3'/  &
     2157            '    Specific heat cap. :   c_p   = ',F6.1,' J/(kg K)'/ &
     2158            '    Vapourization heat :   L_v   = ',E8.2,' J/kg')
     2159417 FORMAT ('    Geograph. longitude :   lambda = ',F4.1,' degr')
     2160418 FORMAT (/'    Day of the year at model start :   day_init      =     ',I3 &
     2161            /'    UTC time at model start        :   time_utc_init = ',F7.1' s')
     2162419 FORMAT (//' Land surface model information:'/ &
     2163              ' ------------------------------'/)
    20642164420 FORMAT (/'    Characteristic levels of the initial temperature profile:'// &
    20652165            '       Height:        ',A,'  m'/ &
     
    21202220                       '[0,1000] cm**2/s**3')
    21212221437 FORMAT ('    Droplet collision is switched off')
     2222438 FORMAT (' --> Land surface type  : ',A,/ &
     2223            ' --> Soil porosity type : ',A)
     2224439 FORMAT (/'    Initial soil temperature and moisture profile:'// &
     2225            '       Height:        ',A,'  m'/ &
     2226            '       Temperature:   ',A,'  K'/ &
     2227            '       Moisture:      ',A,'  m**3/m**3'/ &
     2228            '       Root fraction: ',A,'  '/ &
     2229            '       Gridpoint:     ',A)
     2230440 FORMAT (/' --> Dewfall is allowed (default)')
     2231441 FORMAT (' --> Dewfall is inhibited')
     2232442 FORMAT (' --> Soil bottom is closed (water content is conserved, default)')
     2233443 FORMAT (' --> Soil bottom is open (water content is not conserved)')
     2234444 FORMAT (//' Radiation model information:'/                                 &
     2235              ' ----------------------------'/)
     2236445 FORMAT (' --> Using constant net radiation: net_radiation = ', F6.2, '  W/m**2')
     2237446 FORMAT (' --> Simple radiation scheme for clear sky is used (no clouds,',  &
     2238                   ' default)')
     2239447 FORMAT (' --> Radiation scheme:', A)
     2240448 FORMAT (/'    Surface albedo: albedo = ', F5.3)
     2241449 FORMAT  ('    Timestep: dt_radiation = ', F5.2, '  s')
     2242
    21222243450 FORMAT (//' LES / Turbulence quantities:'/ &
    21232244              ' ---------------------------'/)
     
    22002321508 FORMAT ('    Ventilation effects on evaporation of rain drops')
    22012322509 FORMAT ('    Slope limiter used for sedimentation process')
    2202 510 FORMAT ('        Droplet density    :   N_c   = ',F6.1,' 1/cm**3')
    2203 511 FORMAT ('        Sedimentation Courant number:                  '/&
     2323510 FORMAT ('    Droplet density    :   N_c   = ',F6.1,' 1/cm**3')
     2324511 FORMAT ('    Sedimentation Courant number:                  '/&
    22042325            '                               C_s   = ',F3.1,'        ')
    22052326512 FORMAT (/' Date:                 ',A8,6X,'Run:       ',A20/      &
  • palm/trunk/SOURCE/init_3d_model.f90

    r1508 r1551  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Allocation of land surface arrays is now done in the subroutine init_lsm_arrays,
     23! which is part of land_surface_model.
    2324!
    2425! Former revisions:
     
    218219
    219220    USE land_surface_model_mod,                                                &
    220         ONLY:  init_lsm, land_surface
     221        ONLY:  init_lsm, init_lsm_arrays, land_surface
    221222 
    222223    USE ls_forcing_mod
     
    621622
    622623!
     624!-- Allocate land surface model arrays
     625    IF ( land_surface )  THEN
     626       CALL init_lsm_arrays
     627    ENDIF
     628
     629!
    623630!-- Allocate arrays containing the RK coefficient for calculation of
    624631!-- perturbation pressure and turbulent fluxes. At this point values are
  • palm/trunk/SOURCE/land_surface_model.f90

    r1514 r1551  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Flux calculation is now done in prandtl_fluxes. Added support for data output.
     23! Vertical indices have been replaced. Restart runs are now possible. Some
     24! variables have beem renamed. Bugfix in the prognostic equation for the surface
     25! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of
     26! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for
     27! the hydraulic conductivity. Bugfix for root fraction and extraction
     28! calculation
    2329!
    2430! Former revisions:
     
    4349! scheme implemented in the ECMWF IFS model, with modifications according to
    4450! H-TESSEL. The implementation is based on the formulation implemented in the
    45 ! DALES model.
     51! DALES and UCLA-LES models.
    4652!
    4753! To do list:
    4854! -----------
    49 ! - Add support for binary I/O support
    50 ! - Add support for lsm data output
    51 ! - Check for time step criterion
    52 ! - Check use with RK-2 and Euler time-stepping
    53 ! - Adaption for use with cloud physics (liquid water potential temperature)
    54 ! - Check reaction of plants at wilting point and at atmospheric saturation
    55 ! - Consider partial absorption of the net shortwave radiation by the skin layer
     55! - Check dewfall parametrization for fog simulations
     56! - Consider partial absorption of the net shortwave radiation by the surface layer
    5657! - Allow for water surfaces, check performance for bare soils
     58! - Invert indices (running from -3 to 0. Currently: nzb_soil=0, nzt_soil=3))
     59! - Implement surface runoff model (required when performing long-term LES with
     60!   considerable precipitation
     61! Notes:
     62! ------
     63! - No time step criterion is required as long as the soil layers do not become
     64!   too thin
    5765!------------------------------------------------------------------------------!
    5866     USE arrays_3d,                                                            &
     
    6068
    6169     USE cloud_parameters,                                                     &
    62          ONLY: cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v
     70         ONLY:  cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v
    6371
    6472     USE control_parameters,                                                   &
    65          ONLY: dt_3d, humidity, intermediate_timestep_count,                   &
    66                intermediate_timestep_count_max, precipitation, pt_surface,     &
    67                rho_surface, surface_pressure, timestep_scheme, tsc
     73         ONLY:  cloud_physics, dt_3d, humidity, intermediate_timestep_count,   &
     74                initializing_actions, intermediate_timestep_count_max,         &
     75                max_masks, precipitation, pt_surface, rho_surface,             &
     76                roughness_length, surface_pressure, timestep_scheme, tsc,      &
     77                z0h_factor
    6878
    6979     USE indices,                                                              &
    70          ONLY: nxlg, nxrg, nyng, nysg, nzb_s_inner
     80         ONLY:  nxlg, nxrg, nyng, nysg, nzb_s_inner
    7181
    7282     USE kinds
    7383
     84    USE netcdf_control,                                                        &
     85        ONLY:  dots_label, dots_num, dots_unit
     86
    7487     USE radiation_model_mod,                                                  &
    75          ONLY: Rn, SW_in, sigma_SB
    76 
     88         ONLY:  irad_scheme, rad_net, rad_sw_in, sigma_SB
     89
     90     USE statistics,                                                           &
     91         ONLY:  hom, statistic_regions
    7792
    7893    IMPLICIT NONE
     
    8095!
    8196!-- LSM model constants
    82     INTEGER(iwp), PARAMETER :: soil_layers = 4 !: number of soil layers (fixed for now)
     97    INTEGER(iwp), PARAMETER :: nzb_soil = 0, & !: bottom of the soil model (to be switched)
     98                               nzt_soil = 3, & !: top of the soil model (to be switched)
     99                               nzs = 4         !: number of soil layers (fixed for now)
     100
     101    INTEGER(iwp) :: dots_soil = 0  !: starting index for timeseries output
     102
     103    INTEGER(iwp), DIMENSION(0:1) :: id_dim_zs_xy, id_dim_zs_xz, id_dim_zs_yz,  &
     104                                    id_dim_zs_3d, id_var_zs_xy,                &
     105                                    id_var_zs_xz, id_var_zs_yz, id_var_zs_3d
     106                                   
     107    INTEGER(iwp), DIMENSION(1:max_masks,0:1) :: id_dim_zs_mask, id_var_zs_mask
    83108
    84109    REAL(wp), PARAMETER ::                     &
    85               b_CH               = 6.04_wp,    & ! Clapp & Hornberger exponent
    86               lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil
     110              b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
     111              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil   
    87112              lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix
    88113              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water
    89114              psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
    90               rhoC_soil          = 2.19E6_wp,  & ! volumetric heat capacity of soil
    91               rhoC_water         = 4.20E6_wp,  & ! volumetric heat capacity of water
     115              rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil
     116              rho_c_water        = 4.20E6_wp,  & ! volumetric heat capacity of water
    92117              m_max_depth        = 0.0002_wp     ! Maximum capacity of the water reservoir (m)
    93118
     
    99124
    100125    LOGICAL :: conserve_water_content = .TRUE., & !: open or closed bottom surface for the soil model
     126               dewfall = .TRUE.,                & !: allow/inhibit dewfall
    101127               land_surface = .FALSE.             !: flag parameter indicating wheather the lsm is used
    102128
    103129!   value 9999999.9_wp -> generic available or user-defined value must be set
    104130!   otherwise -> no generic variable and user setting is optional
    105     REAL(wp) :: alpha_VanGenuchten = 0.0_wp,            & !: NAMELIST alpha_VG
    106                 canopy_resistance_coefficient = 0.0_wp, & !: NAMELIST gD
    107                 C_skin   = 20000.0_wp,                  & !: Skin heat capacity
     131    REAL(wp) :: alpha_vangenuchten = 9999999.9_wp,      & !: NAMELIST alpha_vg
     132                canopy_resistance_coefficient = 9999999.9_wp, & !: NAMELIST g_d
     133                c_surface   = 20000.0_wp,               & !: Surface (skin) heat capacity
    108134                drho_l_lv,                              & !: (rho_l * l_v)**-1
    109135                exn,                                    & !: value of the Exner function
    110136                e_s = 0.0_wp,                           & !: saturation water vapour pressure
    111                 field_capacity = 0.0_wp,                & !: NAMELIST m_fc
    112                 f_shortwave_incoming = 9999999.9_wp,    & !: NAMELIST f_SW_in
    113                 hydraulic_conductivity = 0.0_wp,        & !: NAMELIST gamma_w_sat
    114                 Ke = 0.0_wp,                            & !: Kersten number
    115                 lambda_skin_stable = 9999999.9_wp,      & !: NAMELIST lambda_skin_s
    116                 lambda_skin_unstable = 9999999.9_wp,    & !: NAMELIST lambda_skin_u
    117                 leaf_area_index = 9999999.9_wp,         & !: NAMELIST LAI
    118                 l_VanGenuchten = 0.0_WP,                & !: NAMELIST l_VG
    119                 min_canopy_resistance = 110.0_wp,       & !: NAMELIST r_s_min
    120                 m_total = 0.0_wp,                       & !: weighed total water content of the soil (m3/m3)
    121                 n_VanGenuchten = 0.0_WP,                & !: NAMELIST n_VG
     137                field_capacity = 9999999.9_wp,          & !: NAMELIST m_fc
     138                f_shortwave_incoming = 9999999.9_wp,    & !: NAMELIST f_sw_in
     139                hydraulic_conductivity = 9999999.9_wp,  & !: NAMELIST gamma_w_sat
     140                ke = 0.0_wp,                            & !: Kersten number
     141                lambda_surface_stable = 9999999.9_wp,   & !: NAMELIST lambda_surface_s
     142                lambda_surface_unstable = 9999999.9_wp, & !: NAMELIST lambda_surface_u
     143                leaf_area_index = 9999999.9_wp,         & !: NAMELIST lai
     144                l_vangenuchten = 9999999.9_wp,          & !: NAMELIST l_vg
     145                min_canopy_resistance = 9999999.9_wp,   & !: NAMELIST r_canopy_min
     146                min_soil_resistance = 50.0_wp,          & !: NAMELIST r_soil_min
     147                m_total = 0.0_wp,                       & !: weighted total water content of the soil (m3/m3)
     148                n_vangenuchten = 9999999.9_wp,          & !: NAMELIST n_vg
    122149                q_s = 0.0_wp,                           & !: saturation specific humidity
    123                 residual_moisture = 0.0_wp,             & !: NAMELIST m_res
     150                residual_moisture = 9999999.9_wp,       & !: NAMELIST m_res
    124151                rho_cp,                                 & !: rho_surface * cp
    125152                rho_lv,                                 & !: rho * l_v
    126153                rd_d_rv,                                & !: r_d / r_v
    127                 saturation_moisture = 0.0_wp,           & !: NAMELIST m_sat
     154                saturation_moisture = 9999999.9_wp,     & !: NAMELIST m_sat
    128155                vegetation_coverage = 9999999.9_wp,     & !: NAMELIST c_veg
    129                 wilting_point = 0.0_wp                    !: NAMELIST m_wilt
    130 
    131     REAL(wp), DIMENSION(0:soil_layers-1) :: &
     156                wilting_point = 9999999.9_wp,           & !: NAMELIST m_wilt
     157                z0_eb  = 9999999.9_wp,                  & !: NAMELIST z0 (lsm_par)
     158                z0h_eb = 9999999.9_wp                    !: NAMELIST z0h (lsm_par)
     159
     160    REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: &
    132161              ddz_soil,                     & !: 1/dz_soil
    133162              ddz_soil_stag,                & !: 1/dz_soil_stag
     
    135164              dz_soil_stag,                 & !: soil grid spacing (edge-edge)
    136165              root_extr = 0.0_wp,           & !: root extraction
    137               root_fraction = (/0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp/), & !: distribution of root surface area to the individual soil layers
    138               soil_level = (/0.07_wp, 0.28_wp, 1.00_wp,  2.89_wp/),   & !: soil layer depths (m)
     166              root_fraction = (/9999999.9_wp, 9999999.9_wp,    &
     167                                9999999.9_wp, 9999999.9_wp /), & !: distribution of root surface area to the individual soil layers
     168              zs = (/0.07_wp, 0.28_wp, 1.00_wp,  2.89_wp/),    & !: soil layer depths (m)
    139169              soil_moisture = 0.0_wp          !: soil moisture content (m3/m3)
    140170
    141     REAL(wp), DIMENSION(0:soil_layers) ::   &
     171    REAL(wp), DIMENSION(nzb_soil:nzt_soil+1) ::   &
    142172              soil_temperature = 9999999.9_wp !: soil temperature (K)
    143173
    144174#if defined( __nopointer )
    145     REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: T_0,    & !: skin temperature (K)
    146                                                      T_0_p,  & !: progn. skin temperature (K)
    147                                                      m_liq,  & !: liquid water reservoir (m)
    148                                                      m_liq_p   !: progn. liquid water reservoir (m)
     175    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface,   & !: surface temperature (K)
     176                                                     t_surface_p, & !: progn. surface temperature (K)
     177                                                     m_liq_eb,    & !: liquid water reservoir (m)
     178                                                     m_liq_eb_av, & !: liquid water reservoir (m)
     179                                                     m_liq_eb_p     !: progn. liquid water reservoir (m)
    149180#else
    150     REAL(wp), DIMENSION(:,:), POINTER :: T_0,   &
    151                                          T_0_p, &
    152                                          m_liq, &
    153                                          m_liq_p
    154 
    155     REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: T_0_1, T_0_2,    &
    156                                                      m_liq_1, m_liq_2
     181    REAL(wp), DIMENSION(:,:), POINTER :: t_surface,      &
     182                                         t_surface_p,    &
     183                                         m_liq_eb,       &
     184                                         m_liq_eb_p
     185
     186    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface_1, t_surface_2, &
     187                                                     m_liq_eb_av,              &
     188                                                     m_liq_eb_1, m_liq_eb_2
    157189#endif
    158190
    159191!
    160192!-- Temporal tendencies for time stepping           
    161     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tT_0_m,  & !: skin temperature tendency (K)
    162                                              tm_liq_m   !: liquid water reservoir tendency (m)
     193    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_surface_m,  & !: surface temperature tendency (K)
     194                                             tm_liq_eb_m      !: liquid water reservoir tendency (m)
    163195
    164196!
    165197!-- Energy balance variables               
    166198    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::                                   &
    167               alpha_VG,      & !: coef. of Van Genuchten
    168               c_liq,         & !: liquid water coverage (of vegetated area)
    169               c_veg,         & !: vegetation coverage   
    170               f_SW_in,       & !: ?
    171               G,             & !: surface soil heat flux
    172               H,             & !: surface flux of sensible heat
    173               gamma_w_sat,   & !: hydraulic conductivity at saturation
    174               gD,            & !: coefficient for dependence of r_canopy on water vapour pressure deficit
    175               LAI,           & !: leaf area index
    176               LE,            & !: surface flux of latent heat (total)
    177               LE_veg,        & !: surface flux of latent heat (vegetation portion)
    178               LE_soil,       & !: surface flux of latent heat (soil portion)
    179               LE_liq,        & !: surface flux of latent heat (liquid water portion)
    180               lambda_h_sat,  & !: heat conductivity for dry soil
    181               lambda_skin_s, & !: coupling between skin and soil (depends on vegetation type)
    182               lambda_skin_u, & !: coupling between skin and soil (depends on vegetation type)
    183               l_VG,          & !: coef. of Van Genuchten
    184               m_fc,          & !: soil moisture at field capacity (m3/m3)
    185               m_res,         & !: residual soil moisture
    186               m_sat,         & !: saturation soil moisture (m3/m3)
    187               m_wilt,        & !: soil moisture at permanent wilting point (m3/m3)
    188               n_VG,          & !: coef. Van Genuchten 
    189               r_a,           & !: aerodynamic resistance
    190               r_canopy,      & !: canopy resistance
    191               r_soil,        & !: soil resitance
    192               r_soil_min,    & !: minimum soil resistance
    193               r_s,           & !: total surface resistance (combination of r_soil and r_canopy)         
    194               r_s_min          !: minimum canopy resistance
     199              alpha_vg,         & !: coef. of Van Genuchten
     200              c_liq,            & !: liquid water coverage (of vegetated area)
     201              c_liq_av,         & !: average of c_liq
     202              c_soil_av,        & !: average of c_soil
     203              c_veg,            & !: vegetation coverage
     204              c_veg_av,         & !: average of c_veg
     205              f_sw_in,          & !: fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
     206              ghf_eb,           & !: ground heat flux
     207              ghf_eb_av,        & !: average of ghf_eb
     208              gamma_w_sat,      & !: hydraulic conductivity at saturation
     209              g_d,              & !: coefficient for dependence of r_canopy on water vapour pressure deficit
     210              lai,              & !: leaf area index
     211              lai_av,           & !: average of lai
     212              lambda_h_sat,     & !: heat conductivity for dry soil
     213              lambda_surface_s,    & !: coupling between surface and soil (depends on vegetation type)
     214              lambda_surface_u,    & !: coupling between surface and soil (depends on vegetation type)
     215              l_vg,             & !: coef. of Van Genuchten
     216              m_fc,             & !: soil moisture at field capacity (m3/m3)
     217              m_res,            & !: residual soil moisture
     218              m_sat,            & !: saturation soil moisture (m3/m3)
     219              m_wilt,           & !: soil moisture at permanent wilting point (m3/m3)
     220              n_vg,             & !: coef. Van Genuchten 
     221              qsws_eb,          & !: surface flux of latent heat (total)
     222              qsws_eb_av,       & !: average of qsws_eb
     223              qsws_liq_eb,      & !: surface flux of latent heat (liquid water portion)
     224              qsws_liq_eb_av,   & !: average of qsws_liq_eb
     225              qsws_soil_eb,     & !: surface flux of latent heat (soil portion)
     226              qsws_soil_eb_av,  & !: average of qsws_soil_eb
     227              qsws_veg_eb,      & !: surface flux of latent heat (vegetation portion)
     228              qsws_veg_eb_av,   & !: average of qsws_veg_eb
     229              r_a,              & !: aerodynamic resistance
     230              r_canopy,         & !: canopy resistance
     231              r_soil,           & !: soil resitance
     232              r_soil_min,       & !: minimum soil resistance
     233              r_s,              & !: total surface resistance (combination of r_soil and r_canopy)         
     234              r_canopy_min,     & !: minimum canopy (stomatal) resistance
     235              shf_eb,           & !: surface flux of sensible heat
     236              shf_eb_av           !: average of shf_eb
    195237
    196238    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
     
    198240              lambda_w, &   !: hydraulic diffusivity of soil (?)
    199241              gamma_w,  &   !: hydraulic conductivity of soil (?)
    200               rhoC_total    !: volumetric heat capacity of the actual soil matrix (?)
     242              rho_c_total   !: volumetric heat capacity of the actual soil matrix (?)
    201243
    202244#if defined( __nopointer )
    203245    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
    204               T_soil,    & !: Soil temperature (K)
    205               T_soil_p,  & !: Prog. soil temperature (K)
     246              t_soil,    & !: Soil temperature (K)
     247              t_soil_av, & !: Average of t_soil
     248              t_soil_p,  & !: Prog. soil temperature (K)
    206249              m_soil,    & !: Soil moisture (m3/m3)
     250              m_soil_av, & !: Average of m_soil
    207251              m_soil_p     !: Prog. soil moisture (m3/m3)
    208252#else
    209253    REAL(wp), DIMENSION(:,:,:), POINTER ::                                     &
    210               T_soil, T_soil_p, &
     254              t_soil, t_soil_p, &
    211255              m_soil, m_soil_p   
    212256
    213257    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
    214               T_soil_1, T_soil_2, &
    215               m_soil_1, m_soil_2
     258              t_soil_av, t_soil_1, t_soil_2,                                  &
     259              m_soil_av, m_soil_1, m_soil_2
    216260
    217261
     
    220264
    221265    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
    222               tT_soil_m, & !: T_soil storage array
     266              tt_soil_m, & !: t_soil storage array
    223267              tm_soil_m, & !: m_soil storage array
    224268              root_fr      !: root fraction (sum=1)
    225269
    226 !
    227 !--  Land surface parameters according to the following classes (veg_type)
    228 !--  (0 user defined)
    229 !--  1 crops, mixed farming
    230 !--  2 short grass
    231 !--  3 evergreen needleleaf trees
    232 !--  4 deciduous needleleaf trees
    233 !--  5 evergreen broadleaf trees
    234 !--  6 deciduous broadleaf trees
    235 !--  7 tall grass
    236 !--  8 desert
    237 !--  9 tundra
    238 !-- 10 irrigated crops
    239 !-- 11 semidesert
    240 !-- 12 ice caps and glaciers
    241 !-- 13 bogs and marshes
    242 !-- 14 inland water
    243 !-- 15 ocean
    244 !-- 16 evergreen shrubs
    245 !-- 17 deciduous shrubs
    246 !-- 18 mixed forest/woodland
    247 !-- 19 interrupted forest
    248 
    249 !
    250 !-- Land surface parameters I     r_s_min,     LAI,   c_veg,      gD
     270
     271!
     272!-- Predefined Land surface classes (veg_type)
     273    CHARACTER(26), DIMENSION(0:19) :: veg_type_name = (/          &
     274                                   'user defined',                & ! 0
     275                                   'crops, mixed farming',        & !  1
     276                                   'short grass',                 & !  2
     277                                   'evergreen needleleaf trees',  & !  3
     278                                   'deciduous needleleaf trees',  & !  4
     279                                   'evergreen broadleaf trees' ,  & !  5
     280                                   'deciduous broadleaf trees',   & !  6
     281                                   'tall grass',                  & !  7
     282                                   'desert',                      & !  8
     283                                   'tundra',                      & !  9
     284                                   'irrigated crops',             & ! 10
     285                                   'semidesert',                  & ! 11
     286                                   'ice caps and glaciers' ,      & ! 12
     287                                   'bogs and marshes',            & ! 13
     288                                   'inland water',                & ! 14
     289                                   'ocean',                       & ! 15
     290                                   'evergreen shrubs',            & ! 16
     291                                   'deciduous shrubs',            & ! 17
     292                                   'mixed forest/woodland',       & ! 18
     293                                   'interrupted forest'           & ! 19
     294                                                       /)
     295
     296!
     297!-- Soil model classes (soil_type)
     298    CHARACTER(12), DIMENSION(0:7) :: soil_type_name = (/ &
     299                                   'user defined',        & ! 0
     300                                   'coarse',              & ! 1
     301                                   'medium',              & ! 2
     302                                   'medium-fine',         & ! 3
     303                                   'fine',                & ! 4
     304                                   'very fine' ,          & ! 5
     305                                   'organic',             & ! 6
     306                                   'loamy (CH)'           & ! 7
     307                                                        /)
     308!
     309!-- Land surface parameters according to the respective classes (veg_type)
     310
     311!
     312!-- Land surface parameters I
     313!--                          r_canopy_min,     lai,   c_veg,     g_d
    251314    REAL(wp), DIMENSION(0:3,1:19) :: veg_pars = RESHAPE( (/           &
    252315                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & !  1
     
    257320                                 240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp, & !  6
    258321                                 100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp, & !  7
    259                                  250.0_wp, 0.50_wp, 0.00_wp, 0.00_wp, & !  8
     322                                 250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp, & !  8
    260323                                  80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp, & !  9
    261324                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 10
     
    296359
    297360!
    298 !-- Land surface parameters III lambda_skin_s, lambda_skin_u, f_SW_in
    299     REAL(wp), DIMENSION(0:2,1:19) :: skin_pars = RESHAPE( (/           &
     361!-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
     362    REAL(wp), DIMENSION(0:2,1:19) :: surface_pars = RESHAPE( (/           &
    300363                                      10.0_wp,       10.0_wp, 0.05_wp, & !  1
    301364                                      10.0_wp,       10.0_wp, 0.05_wp, & !  2
     
    345408!
    346409!-- Soil parameters according to the following porosity classes (soil_type)
    347 !-- (0 user defined)
    348 !-- 1 coarse
    349 !-- 2 medium
    350 !-- 3 medium-fine
    351 !-- 4 fine
    352 !-- 5 very fine
    353 !-- 6 organic
    354 !
    355 !-- Soil parameters I           alpha_VG,      l_VG,    n_VG, gamma_w_sat
    356     REAL(wp), DIMENSION(0:3,1:6) :: soil_pars = RESHAPE( (/                &
     410
     411!
     412!-- Soil parameters I           alpha_vg,      l_vg,    n_vg, gamma_w_sat
     413    REAL(wp), DIMENSION(0:3,1:7) :: soil_pars = RESHAPE( (/                &
    357414                                 3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, & ! 1
    358415                                 3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, & ! 2
     
    360417                                 3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, & ! 4
    361418                                 2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, & ! 5
    362                                  1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp  & ! 6
    363                                  /), (/ 4, 6 /) )
     419                                 1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, & ! 6
     420                                 0.00_wp,  0.00_wp,  0.00_wp,  0.57E-6_wp  & ! 7
     421                                 /), (/ 4, 7 /) )
    364422
    365423!
    366424!-- Soil parameters II              m_sat,     m_fc,   m_wilt,    m_res 
    367     REAL(wp), DIMENSION(0:3,1:6) :: m_soil_pars = RESHAPE( (/            &
     425    REAL(wp), DIMENSION(0:3,1:7) :: m_soil_pars = RESHAPE( (/            &
    368426                                 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
    369427                                 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
     
    371429                                 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4
    372430                                 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5
    373                                  0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp  & ! 6
    374                                  /), (/ 4, 6 /) )
     431                                 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & ! 6
     432                                 0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp  & ! 7
     433                                 /), (/ 4, 7 /) )
    375434
    376435
     
    381440
    382441
    383     PUBLIC alpha_VanGenuchten, C_skin, canopy_resistance_coefficient,          &
    384            conserve_water_content,      field_capacity, f_shortwave_incoming,  &
    385            hydraulic_conductivity, init_lsm, lambda_skin_stable,               &
    386            lambda_skin_unstable, land_surface, leaf_area_index,                &
    387            lsm_energy_balance, lsm_soil_model, l_VanGenuchten,                 &
    388            min_canopy_resistance, n_VanGenuchten, residual_moisture,           &
    389            root_fraction, saturation_moisture, soil_level, soil_moisture,      &
    390            soil_temperature, soil_type, vegetation_coverage, veg_type,         &
    391            wilting_point
    392 
    393 #if defined( __nopointer )
    394     PUBLIC m_liq, m_liq_p, m_soil, m_soil_p, T_0, T_0_p, T_soil, T_soil_p
    395 #else
    396     PUBLIC m_liq, m_liq_1, m_liq_2, m_liq_p, m_soil, m_soil_1, m_soil_2,       &
    397            m_soil_p, T_0, T_0_1, T_0_2, T_0_p, T_soil, T_soil_1, T_soil_2,     &
    398            T_soil_p
    399 #endif
     442!
     443!-- Public parameters, constants and initial values
     444    PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient,       &
     445           conserve_water_content, dewfall, field_capacity,                    &
     446           f_shortwave_incoming, hydraulic_conductivity, init_lsm,             &
     447           init_lsm_arrays, lambda_surface_stable, lambda_surface_unstable,    &
     448           land_surface, leaf_area_index, lsm_energy_balance, lsm_soil_model,  &
     449           lsm_swap_timelevel, l_vangenuchten, min_canopy_resistance,          &
     450           min_soil_resistance, n_vangenuchten, residual_moisture, rho_cp,     &
     451           rho_lv, root_fraction, saturation_moisture, soil_moisture,          &
     452           soil_temperature, soil_type, soil_type_name, vegetation_coverage,   &
     453           veg_type, veg_type_name, wilting_point, z0_eb, z0h_eb
     454
     455!
     456!-- Public grid and NetCDF variables
     457    PUBLIC dots_soil, id_dim_zs_xy, id_dim_zs_xz, id_dim_zs_yz,                &
     458           id_dim_zs_3d, id_dim_zs_mask, id_var_zs_xy, id_var_zs_xz,           &
     459           id_var_zs_yz, id_var_zs_3d, id_var_zs_mask, nzb_soil, nzs, nzt_soil,&
     460           zs
     461
     462
     463!
     464!-- Public 2D output variables
     465    PUBLIC c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb, ghf_eb_av,     &
     466           lai, lai_av, qsws_eb, qsws_eb_av, qsws_liq_eb, qsws_liq_eb_av,      &
     467           qsws_soil_eb, qsws_soil_eb_av, qsws_veg_eb, qsws_veg_eb_av,         &
     468           shf_eb, shf_eb_av
     469
     470
     471!
     472!-- Public prognostic variables
     473    PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av
    400474
    401475
     
    412486    END INTERFACE lsm_soil_model
    413487
     488    INTERFACE lsm_swap_timelevel
     489       MODULE PROCEDURE lsm_swap_timelevel
     490    END INTERFACE lsm_swap_timelevel
    414491
    415492 CONTAINS
    416493
     494
     495!------------------------------------------------------------------------------!
     496! Description:
     497! ------------
     498!-- Allocate land surface model arrays and define pointers
     499!------------------------------------------------------------------------------!
     500    SUBROUTINE init_lsm_arrays
     501   
     502
     503       IMPLICIT NONE
     504
     505!
     506!--    Allocate surface and soil temperature / humidity
     507#if defined( __nopointer )
     508       ALLOCATE ( m_liq_eb(nysg:nyng,nxlg:nxrg) )
     509       ALLOCATE ( m_liq_eb_p(nysg:nyng,nxlg:nxrg) )
     510       ALLOCATE ( m_soil(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
     511       ALLOCATE ( m_soil_p(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
     512       ALLOCATE ( t_surface(nysg:nyng,nxlg:nxrg) )
     513       ALLOCATE ( t_surface_p(nysg:nyng,nxlg:nxrg) )
     514       ALLOCATE ( t_soil(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
     515       ALLOCATE ( t_soil_p(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
     516#else
     517       ALLOCATE ( m_liq_eb_1(nysg:nyng,nxlg:nxrg) )
     518       ALLOCATE ( m_liq_eb_2(nysg:nyng,nxlg:nxrg) )
     519       ALLOCATE ( m_soil_1(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
     520       ALLOCATE ( m_soil_2(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
     521       ALLOCATE ( t_surface_1(nysg:nyng,nxlg:nxrg) )
     522       ALLOCATE ( t_surface_2(nysg:nyng,nxlg:nxrg) )
     523       ALLOCATE ( t_soil_1(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
     524       ALLOCATE ( t_soil_2(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
     525#endif
     526
     527!
     528!--    Allocate intermediate timestep arrays
     529       ALLOCATE ( tm_liq_eb_m(nysg:nyng,nxlg:nxrg) )
     530       ALLOCATE ( tm_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
     531       ALLOCATE ( tt_surface_m(nysg:nyng,nxlg:nxrg) )
     532       ALLOCATE ( tt_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
     533
     534!
     535!--    Allocate 2D vegetation model arrays
     536       ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) )
     537       ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
     538       ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
     539       ALLOCATE ( f_sw_in(nysg:nyng,nxlg:nxrg) )
     540       ALLOCATE ( ghf_eb(nysg:nyng,nxlg:nxrg) )
     541       ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
     542       ALLOCATE ( g_d(nysg:nyng,nxlg:nxrg) )
     543       ALLOCATE ( lai(nysg:nyng,nxlg:nxrg) )
     544       ALLOCATE ( l_vg(nysg:nyng,nxlg:nxrg) )
     545       ALLOCATE ( lambda_h_sat(nysg:nyng,nxlg:nxrg) )
     546       ALLOCATE ( lambda_surface_u(nysg:nyng,nxlg:nxrg) )
     547       ALLOCATE ( lambda_surface_s(nysg:nyng,nxlg:nxrg) )
     548       ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
     549       ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
     550       ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
     551       ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
     552       ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) )
     553       ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) )
     554       ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) )
     555       ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) )
     556       ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) )
     557       ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
     558       ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
     559       ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
     560       ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
     561       ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
     562       ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) )
     563       ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) )
     564
     565#if ! defined( __nopointer )
     566!
     567!-- Initial assignment of the pointers
     568       t_soil    => t_soil_1;    t_soil_p    => t_soil_2
     569       t_surface => t_surface_1; t_surface_p => t_surface_2
     570       m_soil    => m_soil_1;    m_soil_p    => m_soil_2
     571       m_liq_eb  => m_liq_eb_1;  m_liq_eb_p  => m_liq_eb_2
     572#endif
     573
     574
     575    END SUBROUTINE init_lsm_arrays
    417576
    418577!------------------------------------------------------------------------------!
     
    432591
    433592!
     593!--    Calculate Exner function
     594       exn       = ( surface_pressure / 1000.0_wp )**0.286_wp
     595
     596
     597!
     598!--    If no cloud physics is used, rho_surface has not been calculated before
     599       IF ( .NOT. cloud_physics )  THEN
     600          rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
     601       ENDIF
     602
     603!
    434604!--    Calculate frequently used parameters
    435605       rho_cp    = cp * rho_surface
    436606       rd_d_rv   = r_d / r_v
    437607       rho_lv    = rho_surface * l_v
    438        drho_l_lv = 1.0 / (rho_l * l_v)
    439 
    440 !
    441 !--    Allocate skin and soil temperature / humidity
    442 #if defined( __nopointer )
    443        ALLOCATE ( T_0(nysg:nyng,nxlg:nxrg) )
    444        ALLOCATE ( T_0_p(nysg:nyng,nxlg:nxrg) )
    445 #else
    446        ALLOCATE ( T_0_1(nysg:nyng,nxlg:nxrg) )
    447        ALLOCATE ( T_0_2(nysg:nyng,nxlg:nxrg) )
    448 #endif
    449 
    450        ALLOCATE ( tT_0_m(nysg:nyng,nxlg:nxrg) )
    451 
    452 #if defined( __nopointer )
    453        ALLOCATE ( T_soil(0:soil_layers,nysg:nyng,nxlg:nxrg) )
    454        ALLOCATE ( T_soil_p(0:soil_layers,nysg:nyng,nxlg:nxrg) )
    455 #else
    456        ALLOCATE ( T_soil_1(0:soil_layers,nysg:nyng,nxlg:nxrg) )
    457        ALLOCATE ( T_soil_2(0:soil_layers,nysg:nyng,nxlg:nxrg) )
    458 #endif
    459 
    460        ALLOCATE ( tT_soil_m(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
    461 
    462 #if defined( __nopointer )
    463        ALLOCATE ( m_liq(nysg:nyng,nxlg:nxrg) )
    464        ALLOCATE ( m_liq_p(nysg:nyng,nxlg:nxrg) )
    465 #else
    466        ALLOCATE ( m_liq_1(nysg:nyng,nxlg:nxrg) )
    467        ALLOCATE ( m_liq_2(nysg:nyng,nxlg:nxrg) )
    468 #endif
    469 
    470        ALLOCATE ( tm_liq_m(nysg:nyng,nxlg:nxrg) )
    471 
    472 #if defined( __nopointer )
    473        ALLOCATE ( m_soil(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
    474        ALLOCATE ( m_soil_p(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
    475 #else
    476        ALLOCATE ( m_soil_1(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
    477        ALLOCATE ( m_soil_2(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
    478 #endif
    479 
    480        ALLOCATE ( tm_soil_m(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
    481 
    482 
    483 #if ! defined( __nopointer )
    484 !
    485 !--    Initial assignment of the pointers
    486        T_soil => T_soil_1; T_soil_p => T_soil_2
    487        T_0 => T_0_1; T_0_p => T_0_2
    488        m_soil => m_soil_1; m_soil_p => m_soil_2
    489        m_liq => m_liq_1; m_liq_p => m_liq_2
    490 #endif
    491 
    492        T_0    = 0.0_wp
    493        T_0_p  = 0.0_wp
    494        tT_0_m = 0.0_wp
    495 
    496        T_soil    = 0.0_wp
    497        T_soil_p  = 0.0_wp
    498        tT_soil_m = 0.0_wp
    499 
    500        m_liq    = 0.0_wp
    501        m_liq_p  = 0.0_wp
    502        tm_liq_m = 0.0_wp
    503 
    504        m_soil    = 0.0_wp
    505        m_soil_p  = 0.0_wp
    506        tm_soil_m = 0.0_wp
    507 
    508 !
    509 !--    Allocate 2D vegetation model arrays
    510        ALLOCATE ( alpha_VG(nysg:nyng,nxlg:nxrg) )
    511        ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
    512        ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
    513        ALLOCATE ( f_SW_in(nysg:nyng,nxlg:nxrg) )
    514        ALLOCATE ( G(nysg:nyng,nxlg:nxrg) )
    515        ALLOCATE ( H(nysg:nyng,nxlg:nxrg) )
    516        ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
    517        ALLOCATE ( gD(nysg:nyng,nxlg:nxrg) )
    518        ALLOCATE ( LAI(nysg:nyng,nxlg:nxrg) )
    519        ALLOCATE ( LE(nysg:nyng,nxlg:nxrg) )
    520        ALLOCATE ( LE_veg(nysg:nyng,nxlg:nxrg) )
    521        ALLOCATE ( LE_soil(nysg:nyng,nxlg:nxrg) )
    522        ALLOCATE ( LE_liq(nysg:nyng,nxlg:nxrg) )
    523        ALLOCATE ( l_VG(nysg:nyng,nxlg:nxrg) )
    524        ALLOCATE ( lambda_h_sat(nysg:nyng,nxlg:nxrg) )
    525        ALLOCATE ( lambda_skin_u(nysg:nyng,nxlg:nxrg) )
    526        ALLOCATE ( lambda_skin_s(nysg:nyng,nxlg:nxrg) )
    527        ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
    528        ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
    529        ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
    530        ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
    531        ALLOCATE ( n_VG(nysg:nyng,nxlg:nxrg) )
    532        ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
    533        ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
    534        ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
    535        ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
    536        ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
    537        ALLOCATE ( r_s_min(nysg:nyng,nxlg:nxrg) )
    538 
    539 !
    540 !--    Set initial and default values
    541        c_liq   = 0.0_wp
    542        c_veg   = 0.0_wp
    543        f_SW_in = 0.05_wp
    544        gD      = 0.0_wp
    545        LAI     = 0.0_wp
    546        lambda_skin_u = 10.0_wp
    547        lambda_skin_s = 10.0_wp
    548 
    549 
    550        G       = 0.0_wp
    551        H       = rho_cp * shf
     608       drho_l_lv = 1.0_wp / (rho_l * l_v)
     609
     610!
     611!--    Set inital values for prognostic quantities
     612       tt_surface_m = 0.0_wp
     613       tt_soil_m    = 0.0_wp
     614       tm_liq_eb_m  = 0.0_wp
     615       c_liq        = 0.0_wp
     616
     617       ghf_eb = 0.0_wp
     618       shf_eb = rho_cp * shf
    552619
    553620       IF ( humidity )  THEN
    554           LE = rho_l * l_v * qsws
     621          qsws_eb = rho_l * l_v * qsws
    555622       ELSE
    556           LE = 0.0_wp
     623          qsws_eb = 0.0_wp
    557624       ENDIF
    558625
    559        LE_veg  = 0.0_wp
    560        LE_soil = LE
    561        LE_liq  = 0.0_wp
     626       qsws_liq_eb  = 0.0_wp
     627       qsws_soil_eb = qsws_eb
     628       qsws_veg_eb  = 0.0_wp
    562629
    563630       r_a        = 50.0_wp
     631       r_s        = 50.0_wp
    564632       r_canopy   = 0.0_wp
    565633       r_soil     = 0.0_wp
    566        r_soil_min = 50.0_wp
    567        r_s        = 110.0_wp
    568        r_s_min    = min_canopy_resistance
    569634
    570635!
    571636!--    Allocate 3D soil model arrays
    572        ALLOCATE ( root_fr(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
    573        ALLOCATE ( lambda_h(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
    574        ALLOCATE ( rhoC_total(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
     637       ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
     638       ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
     639       ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
    575640
    576641       lambda_h = 0.0_wp
     
    578643!--    If required, allocate humidity-related variables for the soil model
    579644       IF ( humidity )  THEN
    580           ALLOCATE ( lambda_w(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
    581           ALLOCATE ( gamma_w(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )   
     645          ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
     646          ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )   
    582647
    583648          lambda_w = 0.0_wp
     
    588653!--    the center of the soil layers, whereas gradients/fluxes are defined
    589654!--    at the edges (_stag)
    590        dz_soil_stag(0) = soil_level(0)
    591 
    592        DO k = 1, soil_layers-1
    593           dz_soil_stag(k) = soil_level(k) - soil_level(k-1)
     655       dz_soil_stag(nzb_soil) = zs(nzb_soil)
     656
     657       DO k = nzb_soil+1, nzt_soil
     658          dz_soil_stag(k) = zs(k) - zs(k-1)
    594659       ENDDO
    595660
    596        DO k = 0, soil_layers-2
     661       DO k = nzb_soil, nzt_soil-1
    597662          dz_soil(k) = 0.5 * (dz_soil_stag(k+1) + dz_soil_stag(k))
    598663       ENDDO
    599        dz_soil(soil_layers-1) = dz_soil_stag(soil_layers-1)
    600 
    601        ddz_soil      = 1.0 / dz_soil
    602        ddz_soil_stag = 1.0 / dz_soil_stag
    603 !
    604 !--    Initialize soil
    605        IF ( soil_type .NE. 0 )  THEN   
    606           alpha_VG    = soil_pars(0,soil_type)
    607           l_VG        = soil_pars(1,soil_type)
    608           n_VG        = soil_pars(2,soil_type)   
    609           gamma_w_sat = soil_pars(3,soil_type)
    610           m_sat       = m_soil_pars(0,soil_type)
    611           m_fc        = m_soil_pars(1,soil_type)   
    612           m_wilt      = m_soil_pars(2,soil_type)
    613           m_res       = m_soil_pars(3,soil_type)
     664       dz_soil(nzt_soil) = dz_soil_stag(nzt_soil)
     665
     666       ddz_soil      = 1.0_wp / dz_soil
     667       ddz_soil_stag = 1.0_wp / dz_soil_stag
     668
     669!
     670!--    Initialize standard soil types. It is possible to overwrite each
     671!--    parameter by setting the respecticy NAMELIST variable to a
     672!--    value /= 9999999.9.
     673       IF ( soil_type /= 0 )  THEN 
     674 
     675          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
     676             alpha_vangenuchten = soil_pars(0,soil_type)
     677          ENDIF
     678
     679          IF ( l_vangenuchten == 9999999.9_wp )  THEN
     680             l_vangenuchten = soil_pars(1,soil_type)
     681          ENDIF
     682
     683          IF ( n_vangenuchten == 9999999.9_wp )  THEN
     684             n_vangenuchten = soil_pars(2,soil_type)           
     685          ENDIF
     686
     687          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
     688             hydraulic_conductivity = soil_pars(3,soil_type)           
     689          ENDIF
     690
     691          IF ( saturation_moisture == 9999999.9_wp )  THEN
     692             saturation_moisture = m_soil_pars(0,soil_type)           
     693          ENDIF
     694
     695          IF ( field_capacity == 9999999.9_wp )  THEN
     696             field_capacity = m_soil_pars(1,soil_type)           
     697          ENDIF
     698
     699          IF ( wilting_point == 9999999.9_wp )  THEN
     700             wilting_point = m_soil_pars(2,soil_type)           
     701          ENDIF
     702
     703          IF ( residual_moisture == 9999999.9_wp )  THEN
     704             residual_moisture = m_soil_pars(3,soil_type)       
     705          ENDIF
     706
     707       ENDIF   
     708
     709       alpha_vg      = alpha_vangenuchten
     710       l_vg          = l_vangenuchten
     711       n_vg          = n_vangenuchten
     712       gamma_w_sat   = hydraulic_conductivity
     713       m_sat         = saturation_moisture
     714       m_fc          = field_capacity
     715       m_wilt        = wilting_point
     716       m_res         = residual_moisture
     717       r_soil_min    = min_soil_resistance
     718
     719!
     720!--    Initial run actions
     721       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
     722
     723          t_soil    = 0.0_wp
     724          m_liq_eb  = 0.0_wp
     725          m_soil    = 0.0_wp
     726
     727!
     728!--       Map user settings of T and q for each soil layer
     729!--       (make sure that the soil moisture does not drop below the permanent
     730!--       wilting point) -> problems with devision by zero)
     731          DO k = nzb_soil, nzt_soil
     732             t_soil(k,:,:)    = soil_temperature(k)
     733             m_soil(k,:,:)    = MAX(soil_moisture(k),m_wilt(:,:))
     734             soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
     735          ENDDO
     736          t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1)
     737
     738!
     739!--       Calculate surface temperature
     740          t_surface = pt_surface * exn
     741
     742!
     743!--       Set artifical values for ts and us so that r_a has its initial value for
     744!--       the first time step
     745          DO  i = nxlg, nxrg
     746             DO  j = nysg, nyng
     747                k = nzb_s_inner(j,i)
     748
     749!
     750!--             Assure that r_a cannot be zero at model start
     751                IF ( pt(k+1,j,i) == pt(k,j,i) )  pt(k+1,j,i) = pt(k+1,j,i)     &
     752                                                 + 1.0E-10_wp
     753
     754                us(j,i) = 0.1_wp
     755                ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i)
     756                shf(j,i) = - us(j,i) * ts(j,i)
     757             ENDDO
     758          ENDDO
     759
     760!
     761!--    Actions for restart runs
    614762       ELSE
    615           alpha_VG    = alpha_VanGenuchten
    616           l_VG        = l_VanGenuchten
    617           n_VG        = n_VanGenuchten
    618           gamma_w_sat = hydraulic_conductivity
    619           m_sat       = saturation_moisture
    620           m_fc        = field_capacity
    621           m_wilt      = wilting_point
    622           m_res       = residual_moisture
    623        ENDIF   
    624 
    625 !
    626 !--    Map user settings of T and q for each soil layer
    627 !--    (make sure that the soil moisture does not drop below the permanent
    628 !--    wilting point) -> problems with devision by zero)
    629        DO k = 0, soil_layers-1
    630           T_soil(k,:,:)  = soil_temperature(k)
    631           m_soil(k,:,:)  = MAX(soil_moisture(k),m_wilt(:,:))
    632        ENDDO
    633        T_soil(soil_layers,:,:) = soil_temperature(soil_layers)
    634 
    635 
    636        exn = ( surface_pressure / 1000.0_wp )**0.286_wp
    637        T_0  = pt_surface * exn
    638 
    639        T_soil_p = T_soil
    640        m_soil_p = m_soil
     763
     764          DO  i = nxlg, nxrg
     765             DO  j = nysg, nyng
     766                k = nzb_s_inner(j,i)               
     767                t_surface(j,i) = pt(k,j,i) * exn
     768             ENDDO
     769          ENDDO
     770
     771       ENDIF
    641772
    642773!
     
    645776                           lambda_h_water ** m_sat(:,:)
    646777
     778
     779
     780
     781       DO k = nzb_soil, nzt_soil
     782          root_fr(k,:,:) = root_fraction(k)
     783       ENDDO
     784
     785       IF ( veg_type /= 0 )  THEN
     786          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
     787             min_canopy_resistance = veg_pars(0,veg_type)
     788          ENDIF
     789          IF ( leaf_area_index == 9999999.9_wp )  THEN
     790             leaf_area_index = veg_pars(1,veg_type)         
     791          ENDIF
     792          IF ( vegetation_coverage == 9999999.9_wp )  THEN
     793             vegetation_coverage = veg_pars(2,veg_type)     
     794          ENDIF
     795          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
     796              canopy_resistance_coefficient= veg_pars(3,veg_type)     
     797          ENDIF
     798          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
     799             lambda_surface_stable = surface_pars(0,veg_type)         
     800          ENDIF
     801          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
     802             lambda_surface_unstable = surface_pars(1,veg_type)       
     803          ENDIF
     804          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
     805             f_shortwave_incoming = surface_pars(2,veg_type)       
     806          ENDIF
     807          IF ( z0_eb == 9999999.9_wp )  THEN
     808             roughness_length = roughness_par(0,veg_type)
     809             z0_eb            = roughness_par(0,veg_type)
     810          ENDIF
     811          IF ( z0h_eb == 9999999.9_wp )  THEN
     812             z0h_eb = roughness_par(1,veg_type)
     813          ENDIF
     814          z0h_factor = z0h_eb / z0_eb
     815
     816          IF ( ANY( root_fraction == 9999999.9_wp ) )  THEN
     817             DO k = nzb_soil, nzt_soil
     818                root_fr(k,:,:) = root_distribution(k,veg_type)
     819                root_fraction(k) = root_distribution(k,veg_type)
     820             ENDDO
     821          ENDIF
     822
     823       ENDIF
     824
    647825!
    648826!--    Initialize vegetation
    649        IF ( veg_type .NE. 0 )  THEN
    650 
    651           r_s_min              = veg_pars(0,veg_type)
    652           LAI                  = veg_pars(1,veg_type)
    653           c_veg                = veg_pars(2,veg_type)
    654           gD                   = veg_pars(3,veg_type)
    655           lambda_skin_s        = skin_pars(0,veg_type)
    656           lambda_skin_u        = skin_pars(1,veg_type)
    657           f_SW_in              = skin_pars(2,veg_type)
    658           z0                   = roughness_par(0,veg_type)
    659           z0h                  = roughness_par(1,veg_type)
    660 
    661 
    662           DO k = 0, soil_layers-1
    663              root_fr(k,:,:) = root_distribution(k,veg_type)
    664           ENDDO
    665 
    666        ELSE
    667 
    668           DO k = 0, soil_layers-1
    669              root_fr(k,:,:) = root_fraction(k)
    670           ENDDO
    671 
    672        ENDIF
     827       r_canopy_min         = min_canopy_resistance
     828       lai                  = leaf_area_index
     829       c_veg                = vegetation_coverage
     830       g_d                  = canopy_resistance_coefficient
     831       lambda_surface_s     = lambda_surface_stable
     832       lambda_surface_u     = lambda_surface_unstable
     833       f_sw_in              = f_shortwave_incoming
     834       z0                   = z0_eb
     835       z0h                  = z0h_eb
    673836
    674837!
     
    676839       CALL user_init_land_surface
    677840
    678 !
    679 !--    Set artifical values for ts and us so that r_a has its initial value for
    680 !--    the first time step
    681        DO  i = nxlg, nxrg
    682           DO  j = nysg, nyng
    683              k = nzb_s_inner(j,i)
    684 !
    685 !--          Assure that r_a cannot be zero at model start
    686              IF ( pt(k+1,j,i) == pt(k,j,i) )  pt(k+1,j,i) = pt(k+1,j,i) + 1.0E-10_wp
    687 
    688              us(j,i) = 0.1_wp
    689              ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i)
    690              shf(j,i) = - us(j,i) * ts(j,i)
    691           ENDDO
    692        ENDDO
     841
     842       t_soil_p    = t_soil
     843       m_soil_p    = m_soil
     844       m_liq_eb_p  = m_liq_eb
     845
     846!--    Store initial profiles of t_soil and m_soil (assuming they are
     847!--    horizontally homogeneous on this PE)
     848       hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil(nzb_soil:nzt_soil,      &
     849                                                nysg,nxlg), 2,                 &
     850                                                statistic_regions+1 )
     851       hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil(nzb_soil:nzt_soil,      &
     852                                                nysg,nxlg), 2,                 &
     853                                                statistic_regions+1 )
    693854
    694855!
    695856!--    Calculate humidity at the surface
    696857       IF ( humidity )  THEN
    697           CALL calc_q0
     858          CALL calc_q_surface
    698859       ENDIF
     860
     861
     862
     863!
     864!--    Add timeseries for land surface model
     865       dots_label(dots_num+1) = "ghf_eb"
     866       dots_label(dots_num+2) = "shf_eb"
     867       dots_label(dots_num+3) = "qsws_eb"
     868       dots_label(dots_num+4) = "qsws_liq_eb"
     869       dots_label(dots_num+5) = "qsws_soil_eb"
     870       dots_label(dots_num+6) = "qsws_veg_eb"
     871       dots_unit(dots_num+1:dots_num+6) = "W/m2"
     872
     873       dots_soil = dots_num + 1
     874       dots_num  = dots_num + 6
     875
    699876
    700877       RETURN
     
    707884! Description:
    708885! ------------
    709 !
     886! Solver for the energy balance at the surface.
     887!
     888! Note: surface fluxes are calculated in the land surface model, but these are
     889! not used in the atmospheric code. The fluxes are calculated afterwards in
     890! prandtl_fluxes using the surface values of temperature and humidity as
     891! provided by the land surface model. In this way, the fluxes in the land
     892! surface model are not equal to the ones calculated in prandtl_fluxes
    710893!------------------------------------------------------------------------------!
    711894    SUBROUTINE lsm_energy_balance
     
    730913                   coef_1,      & !: coef. for prognostic equation
    731914                   coef_2,      & !: coef. for prognostic equation
    732                    f_LE,        & !: factor for LE
    733                    f_LE_veg,    & !: factor for LE_veg
    734                    f_LE_soil,   & !: factor for LE_soil
    735                    f_LE_liq,    & !: factor for LE_liq
    736                    f_H,         & !: factor for H
    737                    lambda_skin, & !: Current value of lambda_skin
    738                    m_liq_max      !: maxmimum value of the liquid water reservoir
     915                   f_qsws,      & !: factor for qsws_eb
     916                   f_qsws_veg,  & !: factor for qsws_veg_eb
     917                   f_qsws_soil, & !: factor for qsws_soil_eb
     918                   f_qsws_liq,  & !: factor for qsws_liq_eb
     919                   f_shf,       & !: factor for shf_eb
     920                   lambda_surface, & !: Current value of lambda_surface
     921                   m_liq_eb_max   !: maxmimum value of the liq. water reservoir
     922
    739923
    740924!
     
    748932
    749933!
    750 !--          Set lambda_skin according to stratification
     934!--          Set lambda_surface according to stratification
    751935             IF ( rif(j,i) >= 0.0_wp )  THEN
    752                 lambda_skin = lambda_skin_s(j,i)
     936                lambda_surface = lambda_surface_s(j,i)
    753937             ELSE
    754                 lambda_skin = lambda_skin_u(j,i)
     938                lambda_surface = lambda_surface_u(j,i)
    755939             ENDIF
    756940
     
    760944!--          time step is used here. Note that this formulation is the
    761945!--          equivalent to the ECMWF formulation using drag coefficients
    762              r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20)
     946             r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) +       &
     947                         1.0E-20_wp)
    763948
    764949!
     
    766951!--          f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
    767952 
    768 !--          f1: correction for incoming shortwave radiation
    769              f1 = MIN(1.0_wp, ( 0.004_wp * SW_in(j,i) + 0.05_wp ) /     &
    770                               (0.81_wp * (0.004_wp * SW_in(j,i) + 1.0_wp) ) )
    771 
    772 !
    773 !--          f2: correction for soil moisture f2=0 for very dry soil
     953!--          f1: correction for incoming shortwave radiation (stomata close at
     954!--          night)
     955             IF ( irad_scheme /= 0 )  THEN
     956                f1 = MIN(1.0_wp, ( 0.004_wp * rad_sw_in(j,i) + 0.05_wp ) /     &
     957                              (0.81_wp * (0.004_wp * rad_sw_in(j,i) + 1.0_wp) ))
     958             ELSE
     959                f1 = 1.0_wp
     960             ENDIF
     961
     962!
     963!--          f2: correction for soil moisture availability to plants (the
     964!--          integrated soil moisture must thus be considered here)
     965!--          f2 = 0 for very dry soils
    774966             m_total = 0.0_wp
    775              DO ks = 0, soil_layers-1
    776                  m_total = m_total + root_fr(ks,j,i) * m_soil(ks,j,i)
     967             DO ks = nzb_soil, nzt_soil
     968                 m_total = m_total + root_fr(ks,j,i)                           &
     969                           * MAX(m_soil(ks,j,i),m_wilt(j,i))
    777970             ENDDO
    778971
    779              IF (  m_total .GT. m_wilt(j,i) .AND. m_total .LE. m_fc(j,i) )  THEN
     972             IF ( m_total .GT. m_wilt(j,i) .AND. m_total .LT. m_fc(j,i) )  THEN
    780973                f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
     974             ELSEIF ( m_total .GE. m_fc(j,i) )  THEN
     975                f2 = 1.0_wp
    781976             ELSE
    782977                f2 = 1.0E-20_wp
     
    785980!
    786981!--          Calculate water vapour pressure at saturation
    787 !--          (T_0 should be replaced by liquid water temp?!)
    788              e_s = 0.01 * 610.78_wp * EXP( 17.269_wp * ( T_0(j,i) - 273.16_wp )&
    789                                            / ( T_0(j,i) - 35.86_wp ) )
     982             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
     983                           - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
    790984
    791985!
    792986!--          f3: correction for vapour pressure deficit
    793              IF ( gD(j,i) .NE. 0.0_wp )  THEN
     987             IF ( g_d(j,i) /= 0.0_wp )  THEN
    794988!
    795989!--             Calculate vapour pressure
    796                 e  = q_p(k+1,j,i) * surface_pressure / 0.622
    797                 f3 = EXP ( -gD(j,i) * (e_s - e) )
     990                e  = q_p(k+1,j,i) * surface_pressure / 0.622_wp
     991                f3 = EXP ( -g_d(j,i) * (e_s - e) )
    798992             ELSE
    799993                f3 = 1.0_wp
     
    801995
    802996!
     997!--          Calculate canopy resistance. In case that c_veg is 0 (bare soils),
     998!--          this calculation is obsolete, as r_canopy is not used below.
    803999!--          To do: check for very dry soil -> r_canopy goes to infinity
    804              r_canopy(j,i)  = r_s_min(j,i) / (LAI(j,i) * f1 * f2 * f3 + 1.0E-20)
    805 
    806 !
    807 !--          Third step: calculate bare soil resistance r_soil
    808              m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *        &
    809                      m_res(j,i)
    810 
    811              f2 = ( m_soil(0,j,i) - m_min ) / ( m_fc(j,i) - m_min )
     1000             r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3      &
     1001                                             + 1.0E-20_wp)
     1002
     1003!
     1004!--          Third step: calculate bare soil resistance r_soil. The Clapp &
     1005!--          Hornberger parametrization does not consider c_veg.
     1006             IF ( soil_type /= 7 )  THEN
     1007                m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *     &
     1008                        m_res(j,i)
     1009             ELSE
     1010                m_min = m_wilt(j,i)
     1011             ENDIF
     1012
     1013             f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
    8121014             f2 = MAX(f2,1.0E-20_wp)
     1015             f2 = MIN(f2,1.0_wp)
    8131016
    8141017             r_soil(j,i) = r_soil_min(j,i) / f2
     
    8161019!
    8171020!--          Calculate fraction of liquid water reservoir
    818              m_liq_max = m_max_depth * LAI(j,i)
    819              c_liq(j,i) = MIN(1.0_wp, m_liq(j,i)/m_liq_max)
    820 
     1021             m_liq_eb_max = m_max_depth * lai(j,i)
     1022             c_liq(j,i) = MIN(1.0_wp, m_liq_eb(j,i) / (m_liq_eb_max+1.0E-20_wp))
     1023             
     1024
     1025!
     1026!--          Calculate saturation specific humidity
    8211027             q_s = 0.622_wp * e_s / surface_pressure
    8221028
    8231029!
    824 !--          In case of dew fall, set resistances to zero.
    825 !--          To do: what does that physically reasoning behind this?
     1030!--          In case of dewfall, set evapotranspiration to zero
     1031!--          All super-saturated water is then removed from the air
     1032             IF ( humidity .AND. dewfall .AND. q_s .LE. q_p(k+1,j,i) )  THEN
     1033                r_canopy(j,i) = 0.0_wp
     1034                r_soil(j,i)   = 0.0_wp
     1035             ENDIF
     1036
     1037!
     1038!--          Calculate coefficients for the total evapotranspiration
     1039             f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/        &
     1040                           (r_a(j,i) + r_canopy(j,i))
     1041
     1042             f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +        &
     1043                                                          r_soil(j,i))
     1044             f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
     1045
     1046
     1047!
     1048!--          If soil moisture is below wilting point, plants do no longer
     1049!--          transpirate.
     1050!              IF ( m_soil(k,j,i) .LT. m_wilt(j,i) )  THEN
     1051!                 f_qsws_veg = 0.0_wp
     1052!              ENDIF
     1053
     1054!
     1055!--          If dewfall is deactivated, vegetation, soil and liquid water
     1056!--          reservoir are not allowed to take up water from the super-saturated
     1057!--          air.
    8261058             IF ( humidity )  THEN
    8271059                IF ( q_s .LE. q_p(k+1,j,i) )  THEN
    828                    r_canopy(j,i) = 0.0_wp
    829                    r_soil(j,i) = 0.0_wp
     1060                   IF ( .NOT. dewfall )  THEN
     1061                      f_qsws_veg  = 0.0_wp
     1062                      f_qsws_soil = 0.0_wp
     1063                      f_qsws_liq  = 0.0_wp
     1064                   ENDIF
    8301065                ENDIF
    8311066             ENDIF
    8321067
    833 
    834 !
    835 !--          Calculate coefficients for the total evapotranspiration
    836              f_LE_veg  = rho_lv * c_veg(j,i) * (1.0 - c_liq(j,i)) / (r_a(j,i)  &
    837                                                 + r_canopy(j,i))
    838              f_LE_soil = rho_lv * (1.0 - c_veg(j,i)) / (r_a(j,i) + r_soil(j,i))
    839              f_LE_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
    840 
    841 
    842 !
    843 !--          If soil moisture is below wilting point, plants do no longer
    844 !--          transpirate.
    845              IF ( m_soil(k,j,i) .LT. m_wilt(j,i) )  THEN
    846                 f_LE_veg = 0.0
     1068             f_shf  = rho_cp / r_a(j,i)
     1069             f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
     1070
     1071!
     1072!--          Calculate derivative of q_s for Taylor series expansion
     1073             e_s_dT = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -        &
     1074                              17.269_wp*(t_surface(j,i) - 273.16_wp)           &
     1075                              / (t_surface(j,i) - 35.86_wp)**2 )
     1076
     1077             dq_s_dT = 0.622_wp * e_s_dT / surface_pressure
     1078
     1079             T_1 = pt_p(k+1,j,i) * exn
     1080
     1081!
     1082!--          Add LW up so that it can be removed in prognostic equation
     1083             rad_net(j,i) = rad_net(j,i) + sigma_SB * t_surface(j,i) ** 4
     1084
     1085             IF ( humidity )  THEN
     1086
     1087!
     1088!--             Numerator of the prognostic equation
     1089                coef_1 = rad_net(j,i) + 3.0_wp * sigma_SB * t_surface(j,i) ** 4&
     1090                         + f_shf / exn * T_1 + f_qsws * ( q_p(k+1,j,i) - q_s   &
     1091                         + dq_s_dT * t_surface(j,i) ) + lambda_surface         &
     1092                         * t_soil(nzb_soil,j,i)
     1093
     1094!
     1095!--             Denominator of the prognostic equation
     1096                coef_2 = 4.0_wp * sigma_SB * t_surface(j,i) ** 3 + f_qsws      &
     1097                         * dq_s_dT + lambda_surface + f_shf / exn
     1098
     1099             ELSE
     1100
     1101!
     1102!--             Numerator of the prognostic equation
     1103                coef_1 = rad_net(j,i) + 3.0_wp * sigma_SB * t_surface(j,i) ** 4&
     1104                         + f_shf / exn * T_1 + lambda_surface                  &
     1105                         * t_soil(nzb_soil,j,i)
     1106
     1107!
     1108!--             Denominator of the prognostic equation
     1109                coef_2 = 4.0_wp * sigma_SB * t_surface(j,i) ** 3               &
     1110                         + lambda_surface + f_shf / exn
     1111
    8471112             ENDIF
    8481113
    849              f_H  = rho_cp / r_a(j,i)
    850              f_LE = f_LE_veg + f_LE_soil + f_LE_liq
    851        
    852 !
    853 !--          Calculate derivative of q_s for Taylor series expansion
    854              e_s_dT = e_s * ( 17.269_wp / (T_0(j,i) - 35.86_wp) -              &
    855                               17.269_wp*(T_0(j,i) - 273.16_wp) / (T_0(j,i)     &
    856                               - 35.86_wp)**2 )
    857 
    858              dq_s_dT = 0.622_wp * e_s_dT / surface_pressure
    859 
    860              T_1 = pt_p(k+1,j,i) * exn
    861 
    862 !
    863 !--          Add LW up so that it can be removed in prognostic equation
    864              Rn(j,i) = Rn(j,i) + sigma_SB * T_0(j,i) ** 4
    865 
    866              IF ( humidity )  THEN
    867 
    868 !
    869 !--             Numerator of the prognostic equation
    870                 coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H     &
    871                          / exn * T_1 + f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT   &
    872                          * T_0(j,i) ) + lambda_skin * T_soil(0,j,i)
    873 
    874 !
    875 !--             Denominator of the prognostic equation
    876                 coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3 + f_LE * dq_s_dT    &
    877                          + lambda_skin + f_H / exn
    878 
    879              ELSE
    880 
    881 !
    882 !--             Numerator of the prognostic equation
    883                 coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H /   &
    884                          exn * T_1 + lambda_skin * T_soil(0,j,i)
    885 
    886 !
    887 !--             Denominator of the prognostic equation
    888                 coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3                     &
    889                          + lambda_skin + f_H / exn
    890 
    891              ENDIF
    892 
    8931114             tend = 0.0_wp
    8941115
    8951116!
    896 !--          Implicit solution when the skin layer has no heat capacity,
     1117!--          Implicit solution when the surface layer has no heat capacity,
    8971118!--          otherwise use RK3 scheme.
    898              T_0_p(j,i) = ( coef_1 * dt_3d * tsc(2) + C_skin * T_0(j,i) ) /    &
    899                           ( C_skin + coef_2 * dt_3d * tsc(2) )
    900 
     1119             t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface *        &
     1120                                t_surface(j,i) ) / ( c_surface + coef_2 * dt_3d&
     1121                                * tsc(2) )
    9011122!
    9021123!--          Add RK3 term
    903              T_0_p(j,i) = T_0_p(j,i) + dt_3d * tsc(3) * tT_soil_m(0,j,i)
     1124             t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)              &
     1125                                * tt_surface_m(j,i)
    9041126
    9051127!
    9061128!--          Calculate true tendency
    907              tend = (T_0_p(j,i) - T_0(j,i) - tsc(3) * tT_0_m(j,i)) / (dt_3d    &
    908                       * tsc(2))
    909 
    910 !
    911 !--          Calculate T_0 tendencies for the next Runge-Kutta step
     1129             tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)        &
     1130                    * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
     1131!
     1132!--          Calculate t_surface tendencies for the next Runge-Kutta step
    9121133             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    9131134                IF ( intermediate_timestep_count == 1 )  THEN
    914                    tT_0_m(j,i) = tend
     1135                   tt_surface_m(j,i) = tend
    9151136                ELSEIF ( intermediate_timestep_count <                         &
    9161137                         intermediate_timestep_count_max )  THEN
    917                    tT_0_m(j,i) = -9.5625_wp * tend + 5.3125_wp * tT_0_m(j,i)
     1138                   tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp           &
     1139                                       * tt_surface_m(j,i)
    9181140                ENDIF
    9191141             ENDIF
    9201142
    921              pt_p(k,j,i) = T_0_p(j,i) / exn
     1143             pt_p(k,j,i) = t_surface_p(j,i) / exn
    9221144!
    9231145!--          Calculate fluxes
    924              Rn(j,i)        = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i)**4        &
    925                               - 4.0_wp * sigma_SB * T_0(j,i)**3 * T_0_p(j,i)
    926              G(j,i)         = lambda_skin * (T_0_p(j,i) - T_soil(0,j,i))
    927              H(j,i)         = - f_H  * ( pt_p(k+1,j,i) - pt_p(k,j,i) )
     1146             rad_net(j,i)   = rad_net(j,i) + 3.0_wp * sigma_SB                 &
     1147                              * t_surface(j,i)**4 - 4.0_wp * sigma_SB          &
     1148                              * t_surface(j,i)**3 * t_surface_p(j,i)
     1149             ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
     1150                              - t_soil(nzb_soil,j,i))
     1151             shf_eb(j,i)    = - f_shf  * ( pt_p(k+1,j,i) - pt_p(k,j,i) )
    9281152
    9291153             IF ( humidity )  THEN
    930                 LE(j,i)        = - f_LE      * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
    931                                    * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
    932 
    933                 LE_veg(j,i)    = - f_LE_veg  * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
    934                                    * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
    935                 LE_soil(j,i)   = - f_LE_soil * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
    936                                    * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
    937                 LE_liq(j,i)    = - f_LE_liq  * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
    938                                    * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
     1154                qsws_eb(j,i)  = - f_qsws    * ( q_p(k+1,j,i) - q_s + dq_s_dT   &
     1155                                * t_surface(j,i) - dq_s_dT * t_surface_p(j,i) )
     1156
     1157                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( q_p(k+1,j,i) - q_s       &
     1158                                    + dq_s_dT * t_surface(j,i) - dq_s_dT       &
     1159                                    * t_surface_p(j,i) )
     1160
     1161                qsws_soil_eb(j,i) = - f_qsws_soil * ( q_p(k+1,j,i) - q_s       &
     1162                                    + dq_s_dT * t_surface(j,i) - dq_s_dT       &
     1163                                    * t_surface_p(j,i) )
     1164
     1165                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( q_p(k+1,j,i) - q_s       &
     1166                                    + dq_s_dT * t_surface(j,i) - dq_s_dT       &
     1167                                    * t_surface_p(j,i) )
    9391168             ENDIF
    9401169
    941 !              IF ( i == 1 .AND. j == 1 )  THEN
    942 !                 PRINT*, "Rn", Rn(j,i)
    943 !                  PRINT*, "H", H(j,i)
    944 !                 PRINT*, "LE", LE(j,i)
    945 !                 PRINT*, "LE_liq", LE_liq(j,i)
    946 !                 PRINT*, "LE_veg", LE_veg(j,i)
    947 !                 PRINT*, "LE_soil", LE_soil(j,i)
    948 !                 PRINT*, "G", G(j,i)
    949 !              ENDIF
    950 
    9511170!
    9521171!--          Calculate the true surface resistance
    953              IF ( LE(j,i) .EQ. 0.0 )  THEN
    954                 r_s(j,i) = 1.0E10
     1172             IF ( qsws_eb(j,i) .EQ. 0.0_wp )  THEN
     1173                r_s(j,i) = 1.0E10_wp
    9551174             ELSE
    956                 r_s(j,i) = - rho_lv * ( q_p(k+1,j,i) - q_s + dq_s_dT * T_0(j,i)&
    957                            - dq_s_dT * T_0_p(j,i) ) / LE(j,i) - r_a(j,i)
     1175                r_s(j,i) = - rho_lv * ( q_p(k+1,j,i) - q_s + dq_s_dT           &
     1176                           * t_surface(j,i) - dq_s_dT * t_surface_p(j,i) )     &
     1177                           / qsws_eb(j,i) - r_a(j,i)
    9581178             ENDIF
    959 
    960 !
    961 !--          Calculate fluxes in the atmosphere
    962              shf(j,i) = H(j,i) / rho_cp
    9631179
    9641180!
     
    9671183             IF ( humidity )  THEN
    9681184!
    969 !--             If precipitation is activated, add rain water to LE_liq.
     1185!--             If precipitation is activated, add rain water to qsws_liq_eb.
    9701186!--             precipitation_rate is given in mm.
    9711187                IF ( precipitation )  THEN
    972                    LE_liq(j,i) = LE_liq(j,i) + precipitation_rate(j,i)         &
    973                                                * 0.001_wp * rho_l * l_v
     1188                   qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                         &
     1189                                       + precipitation_rate(j,i) * 0.001_wp    &
     1190                                       * rho_l * l_v
    9741191                ENDIF
    9751192!
     
    9771194                IF ( q_s .LE. q_p(k+1,j,i))  THEN
    9781195!
    979 !--                Check if reservoir is full (avoid values > m_liq_max)
    980 !--                In that case, LE_liq goes to LE_soil. In this case
    981 !--                LE_veg is zero anyway (because c_liq = 1), so that tend is
    982 !--                zero and no further check is needed
    983                    IF ( m_liq(j,i) .EQ. m_liq_max )  THEN
    984                       LE_soil(j,i) = LE_soil(j,i) + LE_liq(j,i)
    985                       LE_liq(j,i) = 0.0_wp
     1196!--                Check if reservoir is full (avoid values > m_liq_eb_max)
     1197!--                In that case, qsws_liq_eb goes to qsws_soil_eb. In this
     1198!--                case qsws_veg_eb is zero anyway (because c_liq = 1),       
     1199!--                so that tend is zero and no further check is needed
     1200                   IF ( m_liq_eb(j,i) .EQ. m_liq_eb_max )  THEN
     1201                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
     1202                                           + qsws_liq_eb(j,i)
     1203                      qsws_liq_eb(j,i)  = 0.0_wp
    9861204                   ENDIF
    9871205
    9881206!
    989 !--                In case LE_veg becomes negative (unphysical behavior), let
    990 !--                the water enter the liquid water reservoir as dew on the
     1207!--                In case qsws_veg_eb becomes negative (unphysical behavior),
     1208!--                let the water enter the liquid water reservoir as dew on the
    9911209!--                plant
    992                    IF ( LE_veg(j,i) .LT. 0.0_wp )  THEN
    993                       LE_liq(j,i) = LE_liq(j,i) + LE_veg(j,i)
    994                       LE_veg(j,i) = 0.0_wp
     1210                   IF ( qsws_veg_eb(j,i) .LT. 0.0_wp )  THEN
     1211                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
     1212                      qsws_veg_eb(j,i) = 0.0_wp
    9951213                   ENDIF
    9961214                ENDIF                   
    9971215 
    998                 tend = - LE_liq(j,i) * drho_l_lv
    999 
    1000                 m_liq_p(j,i) = m_liq(j,i) + dt_3d * ( tsc(2) * tend            &
    1001                                                    + tsc(3) * tm_liq_m(j,i) )
     1216                tend = - qsws_liq_eb(j,i) * drho_l_lv
     1217
     1218                m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend      &
     1219                                                   + tsc(3) * tm_liq_eb_m(j,i) )
    10021220
    10031221!
    10041222!--             Check if reservoir is overfull -> reduce to maximum
    10051223!--             (conservation of water is violated here)
    1006                 m_liq_p(j,i) = MIN(m_liq_p(j,i),m_liq_max)
     1224                m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
    10071225
    10081226!
    10091227!--             Check if reservoir is empty (avoid values < 0.0)
    10101228!--             (conservation of water is violated here)
    1011                 m_liq_p(j,i) = MAX(m_liq_p(j,i),0.0_wp)
    1012 
    1013 
    1014 !
    1015 !--             Calculate m_liq tendencies for the next Runge-Kutta step
     1229                m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
     1230
     1231
     1232!
     1233!--             Calculate m_liq_eb tendencies for the next Runge-Kutta step
    10161234                IF ( timestep_scheme(1:5) == 'runge' )  THEN
    10171235                   IF ( intermediate_timestep_count == 1 )  THEN
    1018                       tm_liq_m(j,i) = tend
     1236                      tm_liq_eb_m(j,i) = tend
    10191237                   ELSEIF ( intermediate_timestep_count <                      &
    10201238                            intermediate_timestep_count_max )  THEN
    1021                       tm_liq_m(j,i) = -9.5625_wp * tend + 5.3125_wp            &
    1022                                       * tm_liq_m(j,i)
     1239                      tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp         &
     1240                                      * tm_liq_eb_m(j,i)
    10231241                   ENDIF
    10241242                ENDIF
    10251243
    1026 !
    1027 !--             Calculate moisture flux in the atmosphere
    1028                 qsws(j,i) = LE(j,i) / rho_lv
    1029 
    10301244             ENDIF
    10311245
    10321246          ENDDO
    10331247       ENDDO
    1034 
    1035 
    10361248
    10371249    END SUBROUTINE lsm_energy_balance
     
    10421254! ------------
    10431255!
     1256! Soil model as part of the land surface model. The model predicts soil
     1257! temperature and water content.
    10441258!------------------------------------------------------------------------------!
    10451259    SUBROUTINE lsm_soil_model
     
    10521266       INTEGER(iwp) ::  k   !: running index
    10531267
    1054        REAL(wp)     :: h_VG !: Van Genuchten coef. h
    1055 
    1056        REAL(wp), DIMENSION(0:soil_layers-1) :: gamma_temp,  & !: temp. gamma
     1268       REAL(wp)     :: h_vg !: Van Genuchten coef. h
     1269
     1270       REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp,  & !: temp. gamma
    10571271                                               lambda_temp, & !: temp. lambda
    10581272                                               tend           !: tendency
     
    10601274       DO i = nxlg, nxrg   
    10611275          DO j = nysg, nyng
    1062              DO k = 0, soil_layers-1
     1276             DO k = nzb_soil, nzt_soil
    10631277!
    10641278!--             Calculate volumetric heat capacity of the soil, taking into
    10651279!--             account water content
    1066                 rhoC_total(k,j,i) = (rhoC_soil * (1.0 - m_sat(j,i))            &
    1067                                      + rhoC_water * m_soil(k,j,i))
     1280                rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i))       &
     1281                                     + rho_c_water * m_soil(k,j,i))
    10681282
    10691283!
    10701284!--             Calculate soil heat conductivity at the center of the soil
    10711285!--             layers
    1072                 Ke = 1.0 + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
    1073                 lambda_temp(k) = Ke * (lambda_h_sat(j,i) + lambda_h_dry) +     &
     1286                ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
     1287                lambda_temp(k) = ke * (lambda_h_sat(j,i) + lambda_h_dry) +     &
    10741288                                 lambda_h_dry
    10751289
     
    10791293!--          Calculate soil heat conductivity (lambda_h) at the _stag level
    10801294!--          using linear interpolation
    1081              DO k = 0, soil_layers-2
     1295             DO k = nzb_soil, nzt_soil-1
    10821296                 
    10831297                lambda_h(k,j,i) = lambda_temp(k) +                             &
     
    10861300
    10871301             ENDDO
    1088              lambda_h(soil_layers-1,j,i) = lambda_temp(soil_layers-1)
    1089 
    1090 !
    1091 !--          Prognostic equation for soil temperature T_soil
     1302             lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
     1303
     1304!
     1305!--          Prognostic equation for soil temperature t_soil
    10921306             tend(:) = 0.0_wp
    1093              tend(0) = (1.0/rhoC_total(0,j,i)) *                               &
    1094                        ( lambda_h(0,j,i) * ( T_soil(1,j,i) - T_soil(0,j,i) )   &
    1095                          * ddz_soil(0) + G(j,i) ) * ddz_soil_stag(0)
    1096 
    1097              DO  k = 1, soil_layers-1
    1098                 tend(k) = (1.0/rhoC_total(k,j,i))                              &
     1307             tend(0) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *                    &
     1308                       ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)     &
     1309                         - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil)         &
     1310                         + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
     1311
     1312             DO  k = 1, nzt_soil
     1313                tend(k) = (1.0_wp/rho_c_total(k,j,i))                          &
    10991314                          * (   lambda_h(k,j,i)                                &
    1100                               * ( T_soil(k+1,j,i) - T_soil(k,j,i) )            &
     1315                              * ( t_soil(k+1,j,i) - t_soil(k,j,i) )            &
    11011316                              * ddz_soil(k)                                    &
    11021317                              - lambda_h(k-1,j,i)                              &
    1103                               * ( T_soil(k,j,i) - T_soil(k-1,j,i) )            &
     1318                              * ( t_soil(k,j,i) - t_soil(k-1,j,i) )            &
    11041319                              * ddz_soil(k-1)                                  &
    11051320                            ) * ddz_soil_stag(k)
    11061321             ENDDO
    11071322
    1108              T_soil_p(0:soil_layers-1,j,i) = T_soil(0:soil_layers-1,j,i)       &
     1323             t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)   &
    11091324                                             + dt_3d * ( tsc(2)                &
    11101325                                             * tend(:) + tsc(3)                &
    1111                                              * tT_soil_m(:,j,i) )   
    1112 
    1113 !
    1114 !--          Calculate T_soil tendencies for the next Runge-Kutta step
     1326                                             * tt_soil_m(:,j,i) )   
     1327
     1328!
     1329!--          Calculate t_soil tendencies for the next Runge-Kutta step
    11151330             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    11161331                IF ( intermediate_timestep_count == 1 )  THEN
    1117                    DO  k = 0, soil_layers-1
    1118                       tT_soil_m(k,j,i) = tend(k)
     1332                   DO  k = nzb_soil, nzt_soil
     1333                      tt_soil_m(k,j,i) = tend(k)
    11191334                   ENDDO
    11201335                ELSEIF ( intermediate_timestep_count <                         &
    11211336                         intermediate_timestep_count_max )  THEN
    1122                    DO  k = 0, soil_layers-1
    1123                       tT_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp      &
    1124                                          * tT_soil_m(k,j,i)
     1337                   DO  k = nzb_soil, nzt_soil
     1338                      tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp      &
     1339                                         * tt_soil_m(k,j,i)
    11251340                   ENDDO
    11261341                ENDIF
     
    11281343
    11291344
    1130              DO k = 0, soil_layers-1
     1345             DO k = nzb_soil, nzt_soil
     1346
    11311347!
    11321348!--             Calculate soil diffusivity at the center of the soil layers
    1133                 lambda_temp(k) = (- b_CH * gamma_w_sat(j,i) * psi_sat          &
     1349                lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat          &
    11341350                                  / m_sat(j,i) ) * ( MAX(m_soil(k,j,i),        &
    1135                                   m_wilt(j,i)) / m_sat(j,i) )**(b_CH + 2.0_wp)
    1136 
    1137 !
    1138 !--             Calculate the hydraulic conductivity after Van Genuchten (1980)
    1139                 h_VG = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -          &
    1140                            MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_VG(j,i)      &
    1141                            / (n_VG(j,i)-1.0_wp)) - 1.0_wp                      &
    1142                        )**(1.0_wp/n_VG(j,i)) / alpha_VG(j,i)
    1143 
    1144                 gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +               &
    1145                                 (alpha_VG(j,i)*h_VG)**n_VG(j,i))**(1.0_wp      &
    1146                                 -1.0_wp/n_VG(j,i)) - (alpha_VG(j,i)*h_VG       &
    1147                                 )**(n_VG(j,i)-1.0_wp))**2 )                    &
    1148                                 / ( (1.0_wp + (alpha_VG(j,i)*h_VG)**n_VG(j,i)  &
    1149                                 )**((1.0_wp - 1.0_wp/n_VG(j,i))*(l_VG(j,i)     &
    1150                                 + 2.0)) )
     1351                                  m_wilt(j,i)) / m_sat(j,i) )**(b_ch + 2.0_wp)
     1352
     1353!
     1354!--             Parametrization of Van Genuchten
     1355                IF ( soil_type /= 7 )  THEN
     1356!
     1357!--                Calculate the hydraulic conductivity after Van Genuchten
     1358!--                (1980)
     1359                   h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -       &
     1360                              MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_vg(j,i)   &
     1361                              / (n_vg(j,i)-1.0_wp)) - 1.0_wp                   &
     1362                          )**(1.0_wp/n_vg(j,i)) / alpha_vg(j,i)
     1363
     1364                   gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +            &
     1365                                   (alpha_vg(j,i)*h_vg)**n_vg(j,i))**(1.0_wp   &
     1366                                   -1.0_wp/n_vg(j,i)) - (alpha_vg(j,i)*h_vg    &
     1367                                   )**(n_vg(j,i)-1.0_wp))**2 )                 &
     1368                                   / ( (1.0_wp + (alpha_vg(j,i)*h_vg           &
     1369                                   )**n_vg(j,i))**((1.0_wp - 1.0_wp/n_vg(j,i)) &
     1370                                   *(l_vg(j,i) + 2.0_wp)) )
     1371
     1372!
     1373!--             Parametrization of Clapp & Hornberger
     1374                ELSE
     1375                   gamma_temp(k) = gamma_w_sat(j,i) * (m_soil(k,j,i)           &
     1376                                   / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
     1377                ENDIF
    11511378
    11521379             ENDDO
     
    11561383!
    11571384!--             Calculate soil diffusivity (lambda_w) at the _stag level
    1158 !--             using linear interpolation
    1159                 DO k = 0, soil_layers-2
     1385!--             using linear interpolation. To do: replace this with
     1386!--             ECMWF-IFS Eq. 8.81
     1387                DO k = nzb_soil, nzt_soil-1
    11601388                     
    11611389                   lambda_w(k,j,i) = lambda_temp(k) +                          &
    11621390                                     ( lambda_temp(k+1) - lambda_temp(k) )     &
    1163                                      * 0.5 * dz_soil_stag(k) * ddz_soil(k+1)
     1391                                     * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1)
    11641392                   gamma_w(k,j,i)  = gamma_temp(k) +                           &
    11651393                                     ( gamma_temp(k+1) - gamma_temp(k) )       &
    1166                                      * 0.5 * dz_soil_stag(k) * ddz_soil(k+1)                 
     1394                                     * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1)                 
    11671395
    11681396                ENDDO
     
    11741402!--             in the bottom layer.
    11751403                IF ( conserve_water_content )  THEN
    1176                    gamma_w(soil_layers-1,j,i) = 0.0_wp
     1404                   gamma_w(nzt_soil,j,i) = 0.0_wp
    11771405                ELSE
    1178                    gamma_w(soil_layers-1,j,i) = lambda_temp(soil_layers-1)
     1406                   gamma_w(nzt_soil,j,i) = lambda_temp(nzt_soil)
    11791407                ENDIF     
    11801408
    1181 !--             The root extraction (= root_extr * LE_veg / (rho_l * l_v))
     1409!--             The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v))
    11821410!--             ensures the mass conservation for water. The transpiration of
    11831411!--             plants equals the cumulative withdrawals by the roots in the
    11841412!--             soil. The scheme takes into account the availability of water
    11851413!--             in the soil layers as well as the root fraction in the
    1186 !--             respective layer
    1187 
    1188 !
    1189 !--             Calculate the root extraction (ECMWF 7.69, with some
    1190 !--             modifications)
     1414!--             respective layer. Layer with moisture below wilting point will
     1415!--             not contribute, which reflects the preference of plants to
     1416!--             take water from moister layers.
     1417
     1418!
     1419!--             Calculate the root extraction (ECMWF 7.69, modified so that the
     1420!--             sum of root_extr = 1). The energy balance solver guarantees a
     1421!--             positive transpiration, so that there is no need for an
     1422!--             additional check.
     1423
    11911424                m_total = 0.0_wp
    1192                 DO k = 0, soil_layers-1
    1193                     m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i) *       &
    1194                               dz_soil_stag(k)
    1195 
     1425                DO k = nzb_soil, nzt_soil
     1426                    IF ( m_soil(k,j,i) .GT. m_wilt(j,i) )  THEN
     1427                       m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
     1428                    ENDIF
    11961429                ENDDO 
    11971430
    1198 !
    1199 !--             For conservation of mass, the sum of root_extr must be 1
    1200                 DO k = 0, soil_layers-1
    1201                    root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i)               &
    1202                                   * dz_soil_stag(k) / m_total
     1431                DO k = nzb_soil, nzt_soil
     1432                   IF ( m_soil(k,j,i) .GT. m_wilt(j,i) )  THEN
     1433                      root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total
     1434                   ELSE
     1435                      root_extr(k) = 0.0_wp
     1436                   ENDIF
    12031437                ENDDO
    1204 
    12051438
    12061439!
    12071440!--             Prognostic equation for soil water content m_soil
    12081441                tend(:) = 0.0_wp
    1209                 tend(0) = ( lambda_w(0,j,i) * ( m_soil(1,j,i) - m_soil(0,j,i) )&
    1210                             * ddz_soil(0) - gamma_w(0,j,i) - ( root_extr(0)    &
    1211                             * LE_veg(j,i) + LE_soil(j,i) ) * drho_l_lv         &
    1212                           ) * ddz_soil_stag(0)
    1213 
    1214                 DO  k = 1, soil_layers-2
     1442                tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (                  &
     1443                            m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
     1444                            * ddz_soil(nzb_soil) - gamma_w(nzb_soil,j,i) - (   &
     1445                            root_extr(nzb_soil) * qsws_veg_eb(j,i)             &
     1446                            + qsws_soil_eb(j,i) ) * drho_l_lv )                &
     1447                            * ddz_soil_stag(nzb_soil)
     1448
     1449                DO  k = nzb_soil+1, nzt_soil-1
    12151450                   tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)             &
    12161451                               - m_soil(k,j,i) ) * ddz_soil(k) - gamma_w(k,j,i)&
    12171452                               - lambda_w(k-1,j,i) * (m_soil(k,j,i) -          &
    12181453                               m_soil(k-1,j,i)) * ddz_soil(k-1)                &
    1219                                + gamma_w(k-1,j,i) - (root_extr(k) * LE_veg(j,i)&
    1220                                * drho_l_lv)                                    &
     1454                               + gamma_w(k-1,j,i) - (root_extr(k)              &
     1455                               * qsws_veg_eb(j,i) * drho_l_lv)                 &
    12211456                             ) * ddz_soil_stag(k)
    12221457
    12231458                ENDDO
    1224                 tend(soil_layers-1) = ( - gamma_w(soil_layers-1,j,i)           &
    1225                                         - lambda_w(soil_layers-2,j,i)          &
    1226                                         * (m_soil(soil_layers-1,j,i)           &
    1227                                         - m_soil(soil_layers-2,j,i))           &
    1228                                         * ddz_soil(soil_layers-2)              &
    1229                                         + gamma_w(soil_layers-2,j,i) - (       &
    1230                                           root_extr(soil_layers-1)             &
    1231                                         * LE_veg(j,i) * drho_l_lv      )       &
    1232                                       ) * ddz_soil_stag(soil_layers-1)             
    1233 
    1234                 m_soil_p(0:soil_layers-1,j,i) = m_soil(0:soil_layers-1,j,i)    &
     1459                tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                     &
     1460                                        - lambda_w(nzt_soil-1,j,i)             &
     1461                                        * (m_soil(nzt_soil,j,i)                &
     1462                                        - m_soil(nzt_soil-1,j,i))              &
     1463                                        * ddz_soil(nzt_soil-1)                 &
     1464                                        + gamma_w(nzt_soil-1,j,i) - (          &
     1465                                          root_extr(nzt_soil)                  &
     1466                                        * qsws_veg_eb(j,i) * drho_l_lv  )      &
     1467                                      ) * ddz_soil_stag(nzt_soil)             
     1468
     1469                m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
    12351470                                                + dt_3d * ( tsc(2) * tend(:)   &
    12361471                                                + tsc(3) * tm_soil_m(:,j,i) )   
     
    12441479                IF ( timestep_scheme(1:5) == 'runge' )  THEN
    12451480                   IF ( intermediate_timestep_count == 1 )  THEN
    1246                       DO  k = 0, soil_layers-1
     1481                      DO  k = nzb_soil, nzt_soil
    12471482                         tm_soil_m(k,j,i) = tend(k)
    12481483                      ENDDO
    12491484                   ELSEIF ( intermediate_timestep_count <                      &
    12501485                            intermediate_timestep_count_max )  THEN
    1251                       DO  k = 0, soil_layers-1
     1486                      DO  k = nzb_soil, nzt_soil
    12521487                         tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
    12531488                                     * tm_soil_m(k,j,i)
     
    12641499!--    Calculate surface specific humidity
    12651500       IF ( humidity )  THEN
    1266           CALL calc_q0
     1501          CALL calc_q_surface
    12671502       ENDIF
    12681503
     
    12741509! Description:
    12751510! ------------
    1276 !
     1511! Calculation of specific humidity of the surface layer (surface)
    12771512!------------------------------------------------------------------------------!
    1278     SUBROUTINE calc_q0
     1513    SUBROUTINE calc_q_surface
    12791514
    12801515       IMPLICIT NONE
     
    12881523          DO j = nysg, nyng
    12891524             k = nzb_s_inner(j,i)
    1290 !
    1291 !--          Temporary solution as long as T_0 is prescribed
    1292 
    1293              pt_p(k,j,i) = T_0(j,i) / exn
     1525
    12941526!
    12951527!--          Calculate water vapour pressure at saturation
    1296              e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( T_0(j,i) -         &
    1297                                               273.16_wp ) /  ( T_0(j,i) -      &
    1298                                               35.86_wp ) )
     1528             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
     1529                                              - 273.16_wp ) / ( t_surface(j,i) &
     1530                                              - 35.86_wp ) )
    12991531
    13001532!
     
    13131545       ENDDO
    13141546
    1315     END SUBROUTINE calc_q0
     1547    END SUBROUTINE calc_q_surface
     1548
     1549!------------------------------------------------------------------------------!
     1550! Description:
     1551! ------------
     1552! Swapping of timelevels
     1553!------------------------------------------------------------------------------!
     1554    SUBROUTINE lsm_swap_timelevel ( mod_count )
     1555
     1556       IMPLICIT NONE
     1557
     1558       INTEGER, INTENT(IN) :: mod_count
     1559
     1560#if defined( __nopointer )
     1561
     1562       t_surface    = t_surface_p
     1563       t_soil       = t_soil_p
     1564       IF ( humidity )  THEN
     1565          m_soil    = m_soil_p
     1566          m_liq_eb  = m_liq_eb_p
     1567       ENDIF
     1568
     1569#else
     1570   
     1571       SELECT CASE ( mod_count )
     1572
     1573          CASE ( 0 )
     1574
     1575             t_surface = t_surface_p
     1576             t_soil    = t_soil_p
     1577             IF ( humidity )  THEN
     1578                m_soil    = m_soil_p
     1579                m_liq_eb  = m_liq_eb_p
     1580             ENDIF
     1581
     1582
     1583          CASE ( 1 )
     1584
     1585             t_surface  => t_surface_1; t_surface_p  => t_surface_2
     1586             t_soil     => t_soil_1;    t_soil_p     => t_soil_2
     1587             IF ( humidity )  THEN
     1588                m_soil    => m_soil_1;   m_soil_p    => m_soil_2
     1589                m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_2
     1590             ENDIF
     1591
     1592       END SELECT
     1593#endif
     1594
     1595    END SUBROUTINE lsm_swap_timelevel
    13161596
    13171597
  • palm/trunk/SOURCE/modules.f90

    r1497 r1551  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Increased pr_palm to 120. Increased length of dots_unit and dots_label to 13
     23! digits. Increased length of domask, do2d, and do3d to 20 digits.
    2324!
    2425! Former revisions:
     
    568569    CHARACTER (LEN=20), DIMENSION(11)  ::  netcdf_precision = ' '
    569570
    570     CHARACTER (LEN=10), DIMENSION(max_masks,0:1,100) ::  domask = ' '
    571     CHARACTER (LEN=10), DIMENSION(0:1,100) ::  do2d = ' ', do3d = ' '
     571    CHARACTER (LEN=20), DIMENSION(max_masks,0:1,100) ::  domask = ' '
     572    CHARACTER (LEN=20), DIMENSION(0:1,100) ::  do2d = ' ', do3d = ' '
    572573
    573574    INTEGER(iwp) ::  abort_mode = 1, average_count_pr = 0, average_count_sp = 0, &
     
    11081109    INTEGER(iwp) ::  dots_num = 23
    11091110
    1110     CHARACTER (LEN=6), DIMENSION(dopr_norm_num) ::  dopr_norm_names =   &
    1111          (/ 'wpt0  ', 'ws2   ', 'tsw2  ', 'ws3   ', 'ws2tsw', 'wstsw2', &
     1111    CHARACTER (LEN=6), DIMENSION(dopr_norm_num) ::  dopr_norm_names =          &
     1112         (/ 'wpt0  ', 'ws2   ', 'tsw2  ', 'ws3   ', 'ws2tsw', 'wstsw2',        &
    11121113            'z_i   ' /)
    11131114
    1114     CHARACTER (LEN=6), DIMENSION(dopr_norm_num) ::  dopr_norm_longnames = &
    1115          (/ 'wpt0  ', 'w*2   ', 't*w2  ', 'w*3   ', 'w*2t*w', 'w*t*w2',   &
     1115    CHARACTER (LEN=6), DIMENSION(dopr_norm_num) ::  dopr_norm_longnames =      &
     1116         (/ 'wpt0  ', 'w*2   ', 't*w2  ', 'w*3   ', 'w*2t*w', 'w*t*w2',        &
    11161117            'z_i   ' /)
    11171118
    1118     CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_label = &
     1119    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_label =                   &
    11191120          (/ 'tnpt   ', 'x_     ', 'y_     ', 'z_     ', 'z_abs  ', 'u      ', &
    11201121             'v      ', 'w      ', 'u"     ', 'v"     ', 'w"     ', 'npt_up ', &
     
    11231124             'w*2    ', 'u"2    ', 'v"2    ', 'w"2    ', 'npt*2  ' /)
    11241125
    1125     CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_unit = &
     1126    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_unit =                    &
    11261127          (/ 'number ', 'm      ', 'm      ', 'm      ', 'm      ', 'm/s    ', &
    11271128             'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'number ', &
     
    11301131             'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'number2' /)
    11311132
    1132     CHARACTER (LEN=7), DIMENSION(dots_max) :: dots_label = &
    1133           (/ 'E      ', 'E*     ', 'dt     ', 'u*     ', 'th*    ', 'umax   ', &
    1134              'vmax   ', 'wmax   ', 'div_new', 'div_old', 'z_i_wpt', 'z_i_pt ', &
    1135              'w*     ', 'w"pt"0 ', 'w"pt"  ', 'wpt    ', 'pt(0)  ', 'pt(zp) ', &
    1136              'w"u"0  ', 'w"v"0  ', 'w"q"0  ', 'mo_L   ', 'q*     ',            &
     1133    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label =                    &
     1134          (/ 'E            ', 'E*           ', 'dt           ',                &
     1135             'u*           ', 'th*          ', 'umax         ',                &
     1136             'vmax         ', 'wmax         ', 'div_new      ',                &
     1137             'div_old      ', 'z_i_wpt      ', 'z_i_pt       ',                &
     1138             'w*           ', 'w"pt"0       ', 'w"pt"        ',                &
     1139             'wpt          ', 'pt(0)        ', 'pt(zp)       ',                &
     1140             'w"u"0        ', 'w"v"0        ', 'w"q"0        ',                &
     1141             'mo_L         ', 'q*           ',                                 &
    11371142             ( 'unknown', i9 = 1, dots_max-23 ) /)
    11381143
    1139     CHARACTER (LEN=7), DIMENSION(dots_max) :: dots_unit = &
    1140           (/ 'm2/s2  ', 'm2/s2  ', 's      ', 'm/s    ', 'K      ', 'm/s    ', &
    1141              'm/s    ', 'm/s    ', 's-1    ', 's-1    ', 'm      ', 'm      ', &
    1142              'm/s    ', 'K m/s  ', 'K m/s  ', 'K m/s  ', 'K      ', 'K      ', &
    1143              'm2/s2  ', 'm2/s2  ', 'kg m/s ', 'm      ', 'kg/kg  ',            &
     1144    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit =                     &
     1145          (/ 'm2/s2        ', 'm2/s2        ', 's            ',                &
     1146             'm/s          ', 'K            ', 'm/s          ',                &
     1147             'm/s          ', 'm/s          ', 's-1          ',                &
     1148             's-1          ', 'm            ', 'm            ',                &
     1149             'm/s          ', 'K m/s        ', 'K m/s        ',                &
     1150             'K m/s        ', 'K            ', 'K            ',                &
     1151             'm2/s2        ', 'm2/s2        ', 'kg m/s       ',                &
     1152             'm            ', 'kg/kg        ',                                 &
    11441153             ( 'unknown', i9 = 1, dots_max-23 ) /)
    11451154
     
    14221431
    14231432    CHARACTER (LEN=40) ::  region(0:9)
    1424     INTEGER(iwp) ::  pr_palm = 100, statistic_regions = 0
     1433    INTEGER(iwp) ::  pr_palm = 120, statistic_regions = 0
    14251434    INTEGER(iwp) ::  u_max_ijk(3) = -1, v_max_ijk(3) = -1, w_max_ijk(3) = -1
    14261435    LOGICAL ::  flow_statistics_called = .FALSE.
  • palm/trunk/SOURCE/netcdf.f90

    r1354 r1551  
    2323! Current revisions:
    2424! ------------------
    25 !
     25! Added support for land surface model and radiation model output. In the course
     26! of this action a new vertical grid zs (soil) was introduced.
    2627!
    2728! Former revisions:
     
    9798! In case of extend = .TRUE.:
    9899! Find out if dimensions and variables of an existing file match the values
    99 ! of the actual run. If so, get all necessary informations (ids, etc.) from
     100! of the actual run. If so, get all necessary information (ids, etc.) from
    100101! this file.
    101102!
     
    130131
    131132    USE kinds
     133
     134    USE land_surface_model_mod,                                                &
     135        ONLY: land_surface, nzb_soil, nzt_soil, id_dim_zs_xy, id_dim_zs_xz,    &
     136              id_dim_zs_yz, id_dim_zs_3d, id_dim_zs_mask, id_var_zs_xy,        &
     137              id_var_zs_xz, id_var_zs_yz ,id_var_zs_3d, id_var_zs_mask,        &
     138              nzs, zs
    132139
    133140    USE pegrid
     
    181188    INTEGER(iwp) ::  kk                                      !:
    182189    INTEGER(iwp) ::  ns                                      !:
     190    INTEGER(iwp) ::  ns_do                                   !: actual value of ns for soil model data
    183191    INTEGER(iwp) ::  ns_old                                  !:
    184192    INTEGER(iwp) ::  ntime_count                             !:
     
    439447!
    440448!--       In case of non-flat topography define 2d-arrays containing the height
    441 !--       informations
     449!--       information
    442450          IF ( TRIM( topography ) /= 'flat' )  THEN
    443451!
     
    478486
    479487          ENDIF             
    480 
     488 
     489          IF ( land_surface )  THEN
     490!
     491!--          Define vertical coordinate grid (zw grid)
     492             nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'zs_3d', &
     493                                     mask_size(mid,3), id_dim_zs_mask(mid,av) )
     494             CALL handle_netcdf_error( 'netcdf', 536 )
     495
     496             nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zs_3d', NF90_DOUBLE, &
     497                                     id_dim_zs_mask(mid,av), &
     498                                     id_var_zs_mask(mid,av) )
     499             CALL handle_netcdf_error( 'netcdf', 536 )
     500
     501             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_zs_mask(mid,av), &
     502                                  'units', 'meters' )
     503             CALL handle_netcdf_error( 'netcdf', 537 )
     504
     505          ENDIF
    481506
    482507!
     
    521546                   grid_y = 'y'
    522547                   grid_z = 'zw'
     548!
     549!--             soil grid
     550                CASE ( 't_soil', 'm_soil' )
     551
     552                   grid_x = 'x'
     553                   grid_y = 'y'
     554                   grid_z = 'zs'
    523555
    524556                CASE DEFAULT
     
    548580             ELSEIF ( grid_z == 'zw' )  THEN
    549581                id_z = id_dim_zw_mask(mid,av)
     582             ELSEIF ( grid_z == "zs" )  THEN
     583                id_z = id_dim_zs_mask(mid,av)
    550584             ENDIF
    551585
     
    692726
    693727          ENDIF
     728
     729          IF ( land_surface )  THEN
     730!
     731!--          Write zs data (vertical axes for soil model), use negative values
     732!--          to indicate soil depth
     733             ALLOCATE( netcdf_data(mask_size(mid,3)) )
     734
     735             netcdf_data = zs( mask_k_global(mid,:mask_size(mid,3)) )
     736
     737             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zs_mask(mid,av), &
     738                                     netcdf_data, start = (/ 1 /), &
     739                                     count = (/ mask_size(mid,3) /) )
     740             CALL handle_netcdf_error( 'netcdf', 538 )
     741
     742             DEALLOCATE( netcdf_data )
     743
     744          ENDIF
     745
    694746!
    695747!--       restore original parameter file_id (=formal parameter av) into av
     
    9821034!
    9831035!--       In case of non-flat topography define 2d-arrays containing the height
    984 !--       informations
     1036!--       information
    9851037          IF ( TRIM( topography ) /= 'flat' )  THEN
    9861038!
     
    10161068          ENDIF             
    10171069
     1070          IF ( land_surface )  THEN
     1071!
     1072!--          Define vertical coordinate grid (zs grid)
     1073             nc_stat = NF90_DEF_DIM( id_set_3d(av), 'zs_3d', nzt_soil-nzb_soil+1, &
     1074                                     id_dim_zs_3d(av) )
     1075             CALL handle_netcdf_error( 'netcdf', 70 )
     1076
     1077             nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zs_3d', NF90_DOUBLE, &
     1078                                     id_dim_zs_3d(av), id_var_zs_3d(av) )
     1079             CALL handle_netcdf_error( 'netcdf', 71 )
     1080
     1081             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zs_3d(av), 'units', &
     1082                                     'meters' )
     1083             CALL handle_netcdf_error( 'netcdf', 72 )
     1084
     1085          ENDIF
    10181086
    10191087!
     
    10581126                   grid_y = 'y'
    10591127                   grid_z = 'zw'
     1128!
     1129!--             soil grid
     1130                CASE ( 't_soil', 'm_soil' )
     1131
     1132                   grid_x = 'x'
     1133                   grid_y = 'y'
     1134                   grid_z = 'zs'
    10601135
    10611136                CASE DEFAULT
     
    10851160             ELSEIF ( grid_z == 'zw' )  THEN
    10861161                id_z = id_dim_zw_3d(av)
     1162             ELSEIF ( grid_z == 'zs' )  THEN
     1163                id_z = id_dim_zs_3d(av)
    10871164             ENDIF
    10881165
     
    12391316
    12401317             ENDIF
     1318
     1319             IF ( land_surface )  THEN
     1320!
     1321!--             Write zs grid
     1322                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zs_3d(av),  &
     1323                                        - zs(nzb_soil:nzt_soil), start = (/ 1 /), &
     1324                                        count = (/ nzt_soil-nzb_soil+1 /) )
     1325                CALL handle_netcdf_error( 'netcdf', 86 )
     1326             ENDIF
     1327
    12411328          ENDIF
    12421329
     
    14961583          CALL handle_netcdf_error( 'netcdf', 107 )
    14971584
     1585
     1586          IF ( land_surface )  THEN
     1587
     1588             ns_do = 0
     1589             DO WHILE ( section(ns_do+1,1) < nzs )
     1590                ns_do = ns_do + 1
     1591             ENDDO
     1592!
     1593!--          Define vertical coordinate grid (zs grid)
     1594             nc_stat = NF90_DEF_DIM( id_set_xy(av), 'zs_xy', ns_do, id_dim_zs_xy(av) )
     1595             CALL handle_netcdf_error( 'netcdf', 539 )
     1596
     1597             nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zs_xy', NF90_DOUBLE, &
     1598                                     id_dim_zs_xy(av), id_var_zs_xy(av) )
     1599             CALL handle_netcdf_error( 'netcdf', 540 )
     1600
     1601             nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zs_xy(av), 'units', &
     1602                                  'meters' )
     1603             CALL handle_netcdf_error( 'netcdf', 541 )
     1604
     1605          ENDIF
     1606
    14981607!
    14991608!--       Define a pseudo vertical coordinate grid for the surface variables
     
    15771686!
    15781687!--       In case of non-flat topography define 2d-arrays containing the height
    1579 !--       informations
     1688!--       information
    15801689          IF ( TRIM( topography ) /= 'flat' )  THEN
    15811690!
     
    16111720          ENDIF
    16121721
    1613 
    16141722!
    16151723!--       Define the variables
     
    16711779                         grid_y = 'y'
    16721780                         grid_z = 'zw'
     1781!
     1782!--                   soil grid
     1783                      CASE ( 't_soil_xy', 'm_soil_xy' )
     1784                         grid_x = 'x'
     1785                         grid_y = 'y'
     1786                         grid_z = 'zs'
    16731787
    16741788                      CASE DEFAULT
     
    16981812                   ELSEIF ( grid_z == 'zw' )  THEN
    16991813                      id_z = id_dim_zw_xy(av)
     1814                   ELSEIF ( grid_z == 'zs' )  THEN
     1815                      id_z = id_dim_zs_xy(av)
    17001816                   ENDIF
    17011817
     
    18131929
    18141930!
     1931!--             Write zs data
     1932             IF ( land_surface )  THEN
     1933                ns_do = 0
     1934                DO  i = 1, ns
     1935                   IF( section(i,1) == -1 )  THEN
     1936                      netcdf_data(i) = 1.0_wp  ! section averaged along z
     1937                      ns_do = ns_do + 1
     1938                   ELSEIF ( section(i,1) < nzs )  THEN
     1939                      netcdf_data(i) = - zs( section(i,1) )
     1940                      ns_do = ns_do + 1
     1941                   ENDIF
     1942                ENDDO
     1943
     1944                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zs_xy(av), &
     1945                                        netcdf_data(1:ns_do), start = (/ 1 /),    &
     1946                                        count = (/ ns_do /) )
     1947                CALL handle_netcdf_error( 'netcdf', 124 )
     1948
     1949             ENDIF
     1950
     1951!
    18151952!--          Write gridpoint number data
    18161953             netcdf_data(1:ns) = section(1:ns,1)
     
    18942031
    18952032             ENDIF
     2033
     2034
     2035
    18962036          ENDIF
    18972037
     
    22402380
    22412381!
    2242 !--       Define the two z-axes (zu and zw)
     2382!--       Define the three z-axes (zu, zw, and zs)
    22432383          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'zu', nz+2, id_dim_zu_xz(av) )
    22442384          CALL handle_netcdf_error( 'netcdf', 153 )
     
    22632403          CALL handle_netcdf_error( 'netcdf', 158 )
    22642404
     2405          IF ( land_surface )  THEN
     2406
     2407             nc_stat = NF90_DEF_DIM( id_set_xz(av), 'zs', nzs, id_dim_zs_xz(av) )
     2408             CALL handle_netcdf_error( 'netcdf', 542 )
     2409
     2410             nc_stat = NF90_DEF_VAR( id_set_xz(av), 'zs', NF90_DOUBLE, &
     2411                                     id_dim_zs_xz(av), id_var_zs_xz(av) )
     2412             CALL handle_netcdf_error( 'netcdf', 543 )
     2413
     2414             nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_zs_xz(av), 'units', &
     2415                                     'meters' )
     2416             CALL handle_netcdf_error( 'netcdf', 544 )
     2417
     2418          ENDIF
    22652419
    22662420!
     
    23082462                      grid_y = 'y'
    23092463                      grid_z = 'zw'
     2464
     2465!
     2466!--                soil grid
     2467                   CASE ( 't_soil_xz', 'm_soil_xz' )
     2468
     2469                      grid_x = 'x'
     2470                      grid_y = 'y'
     2471                      grid_z = 'zs'
    23102472
    23112473                   CASE DEFAULT
     
    23352497                ELSEIF ( grid_z == 'zw' )  THEN
    23362498                   id_z = id_dim_zw_xz(av)
     2499                ELSEIF ( grid_z == 'zs' )  THEN
     2500                   id_z = id_dim_zs_xz(av)
    23372501                ENDIF
    23382502
     
    25142678                                     count = (/ nz+2 /) )
    25152679             CALL handle_netcdf_error( 'netcdf', 167 )
     2680
     2681!
     2682!--          Write zs data
     2683             IF ( land_surface )  THEN
     2684                netcdf_data(0:nzs-1) = - zs(nzb_soil:nzt_soil)
     2685                nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zs_xz(av), &
     2686                                        netcdf_data(0:nzs), start = (/ 1 /),    &
     2687                                        count = (/ nzt_soil-nzb_soil+1 /) )
     2688               CALL handle_netcdf_error( 'netcdf', 548 )
     2689             ENDIF
     2690
    25162691
    25172692             DEALLOCATE( netcdf_data )
     
    29033078          CALL handle_netcdf_error( 'netcdf', 197 )
    29043079
     3080          IF ( land_surface )  THEN
     3081
     3082             nc_stat = NF90_DEF_DIM( id_set_yz(av), 'zs', nzs, id_dim_zs_yz(av) )
     3083             CALL handle_netcdf_error( 'netcdf', 545 )
     3084
     3085             nc_stat = NF90_DEF_VAR( id_set_yz(av), 'zs', NF90_DOUBLE, &
     3086                                     id_dim_zs_yz(av), id_var_zs_yz(av) )
     3087             CALL handle_netcdf_error( 'netcdf', 546 )
     3088
     3089             nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_zs_yz(av), 'units', &
     3090                                     'meters' )
     3091             CALL handle_netcdf_error( 'netcdf', 547 )
     3092
     3093          ENDIF
     3094
    29053095
    29063096!
     
    29483138                      grid_y = 'y'
    29493139                      grid_z = 'zw'
     3140!
     3141!--                soil grid
     3142                   CASE ( 't_soil_yz', 'm_soil_yz' )
     3143
     3144                      grid_x = 'x'
     3145                      grid_y = 'y'
     3146                      grid_z = 'zs'
    29503147
    29513148                   CASE DEFAULT
     
    29753172                ELSEIF ( grid_z == 'zw' )  THEN
    29763173                   id_z = id_dim_zw_yz(av)
     3174                ELSEIF ( grid_z == 'zs' )  THEN
     3175                   id_z = id_dim_zs_yz(av)
    29773176                ENDIF
    29783177
  • palm/trunk/SOURCE/package_parin.f90

    r1497 r1551  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Several changes in the liste for land surface model and radiation model
    2323!
    2424! Former revisions:
     
    100100
    101101    USE land_surface_model_mod,                                                &
    102         ONLY: alpha_VanGenuchten, C_skin, canopy_resistance_coefficient,       &
    103               conserve_water_content, f_shortwave_incoming,                    &
    104               hydraulic_conductivity, lambda_skin_stable, lambda_skin_unstable,&
    105               leaf_area_index, land_surface, l_VanGenuchten,                   &
    106               min_canopy_resistance, field_capacity, residual_moisture,        &
    107               saturation_moisture, wilting_point, n_VanGenuchten,              &
    108               root_fraction, soil_level, soil_moisture, soil_temperature,      &
    109               soil_type, vegetation_coverage, veg_type
     102        ONLY: alpha_vangenuchten, c_surface, canopy_resistance_coefficient,    &
     103              conserve_water_content, dewfall, f_shortwave_incoming,           &
     104              hydraulic_conductivity, lambda_surface_stable,                   &
     105              lambda_surface_unstable,                                         &
     106              leaf_area_index, land_surface, l_vangenuchten,                   &
     107              min_canopy_resistance, min_soil_resistance, field_capacity,      &
     108              residual_moisture, saturation_moisture, wilting_point,           &
     109              n_vangenuchten, root_fraction, soil_moisture, soil_temperature,  &
     110              soil_type, vegetation_coverage, veg_type, z0_eb, z0h_eb, zs
    110111
    111112    USE kinds
     
    134135
    135136    USE radiation_model_mod,                                                   &
    136         ONLY: albedo, day_init, dt_radiation, lambda, radiation,               &
    137               time_radiation, time_utc_init
     137        ONLY: albedo, day_init, dt_radiation, lambda, net_radiation, radiation,&
     138              radiation_scheme, time_radiation, time_utc_init
    138139               
    139140
     
    171172                                  vc_size_y, vc_size_z
    172173
    173     NAMELIST /lsm_par/            alpha_VanGenuchten, C_skin,                  &
     174    NAMELIST /lsm_par/            alpha_vangenuchten, c_surface,               &
    174175                                  canopy_resistance_coefficient,               &
    175                                   conserve_water_content, f_shortwave_incoming,&
    176                                   hydraulic_conductivity, lambda_skin_stable,  &
    177                                   lambda_skin_unstable, leaf_area_index,       &
    178                                   l_VanGenuchten, min_canopy_resistance,       &
    179                                   field_capacity, residual_moisture,           &
    180                                   saturation_moisture, wilting_point,          &
    181                                   n_VanGenuchten, root_fraction, soil_level,   &
     176                                  conserve_water_content, dewfall,             &
     177                                  f_shortwave_incoming,                        &
     178                                  hydraulic_conductivity,                      &
     179                                  lambda_surface_stable,                       &
     180                                  lambda_surface_unstable, leaf_area_index,    &
     181                                  l_vangenuchten, min_canopy_resistance,       &
     182                                  min_soil_resistance, field_capacity,         &
     183                                  residual_moisture, saturation_moisture,      &
     184                                  wilting_point, n_vangenuchten, root_fraction,&
    182185                                  soil_moisture, soil_temperature, soil_type,  &
    183                                   vegetation_coverage, veg_type
     186                                  vegetation_coverage, veg_type, zs, z0_eb,    &
     187                                  z0h_eb
    184188
    185189
     
    207211                                  write_particle_statistics
    208212
    209     NAMELIST /radiation_par/      lambda, albedo, day_init, dt_radiation,      &
    210                                   time_utc_init
     213    NAMELIST /radiation_par/      albedo, day_init, dt_radiation, lambda,      &
     214                                  net_radiation, radiation_scheme, time_utc_init
    211215
    212216    NAMELIST /spectra_par/        averaging_interval_sp, comp_spectra_level,   &
  • palm/trunk/SOURCE/prandtl_fluxes.f90

    r1497 r1551  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Removed land surface model part. The surface fluxes are now always calculated
     23! within prandtl_fluxes, based on the given surface temperature/humidity (which
     24! is either provided by the land surface model, by large scale forcing data, or
     25! directly prescribed by the user.
    2326!
    2427! Former revisions:
     
    9194
    9295    USE kinds
    93 
    94     USE land_surface_model_mod,                                                &
    95         ONLY: land_surface
    9696
    9797    IMPLICIT NONE
     
    482482!
    483483!-- Compute the vertical kinematic heat flux
    484     IF ( .NOT. constant_heatflux .AND. .NOT. land_surface )  THEN
     484    IF ( .NOT. constant_heatflux )  THEN
    485485       !$OMP PARALLEL DO
    486486       !$acc kernels loop independent
     
    495495!
    496496!-- Compute the vertical water/scalar flux
    497     IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar )       &
    498          .AND. .NOT. land_surface )  THEN
     497    IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar ) )  THEN
    499498       !$OMP PARALLEL DO
    500499       !$acc kernels loop independent
  • palm/trunk/SOURCE/radiation_model.f90

    r1497 r1551  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Added support for data output. Various variables have been renamed. Added
     23! interface for different radiation schemes (currently: clear-sky, constant, and
     24! RRTM (not yet implemented).
    2325!
    2426! Former revisions:
     
    4648    USE kinds
    4749
     50    USE netcdf_control,                                                        &
     51        ONLY:  dots_label, dots_num, dots_unit
     52
    4853
    4954    IMPLICIT NONE
    5055
     56    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtm'
     57
    5158    INTEGER(iwp) :: i, j, k
    5259
    5360
    54     INTEGER(iwp) :: day_init = 1 !: day of the year at model start
     61    INTEGER(iwp) :: day_init     = 172,  & !: day of the year at model start (21/06)
     62                    dots_rad     = 0,    & !: starting index for timeseries output
     63                    irad_scheme  = 0
    5564
    5665    LOGICAL :: radiation = .FALSE.  !: flag parameter indicating wheather the radiation model is used
     
    6170 
    6271    REAL(wp) :: albedo = 0.2_wp,             & !: NAMELIST alpha
    63                 dt_radiation = 9999999.9_wp, &
     72                dt_radiation = 0.0_wp,       & !: radiation model timestep
    6473                exn,                         & !: Exner function
    6574                lon = 0.0_wp,                & !: longitude in radians
     
    6978                decl_3,                      & !: declination coef. 3
    7079                time_utc,                    & !: current time in UTC
    71                 time_utc_init = 0.0_wp,      & !: UTC time at model start
     80                time_utc_init = 43200.0_wp,  & !: