Changeset 1788


Ignore:
Timestamp:
Mar 10, 2016 11:01:04 AM (6 years ago)
Author:
maronga
Message:

added support for water and paved surfaced in land surface model / minor changes

Location:
palm/trunk/SOURCE
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/check_parameters.f90

    r1787 r1788  
    1919! Current revisions:
    2020! -----------------
    21 !
    22 !
     21! Added check for use of most_method = 'lookup' in combination with water
     22! surface presribed in the land surface model. Added output of z0q.
     23! Syntax layout improved.
     24!
    2325! Former revisions:
    2426! -----------------
     
    397399
    398400       IF ( dt_coupling == 9999999.9_wp )  THEN
    399           message_string = 'dt_coupling is not set but required for coup' // &
     401          message_string = 'dt_coupling is not set but required for coup' //   &
    400402                           'ling mode "' //  TRIM( coupling_mode ) // '"'
    401403          CALL message( 'check_parameters', 'PA0003', 1, 2, 0, 6, 0 )
     
    419421#if ! defined( __check )
    420422       IF ( myid == 0 ) THEN
    421           CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, &
     423          CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter,  &
    422424                         ierr )
    423           CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, &
     425          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter,       &
    424426                         status, ierr )
    425427       ENDIF
     
    436438          IF ( myid == 0  ) THEN
    437439             CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
    438              CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, &
     440             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter,    &
    439441                            status, ierr )
    440442          ENDIF   
     
    451453          CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &
    452454                         ierr )
    453           CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, &
     455          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter,       &
    454456                         status, ierr )
    455457       ENDIF
     
    464466#if ! defined( __check )
    465467       IF ( myid == 0 ) THEN
    466           CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, &
     468          CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter,   &
    467469                         ierr )
    468           CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, &
     470          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter,       &
    469471                         status, ierr )
    470472       ENDIF   
     
    483485          CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, &
    484486                         14, comm_inter, ierr )
    485           CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, &
     487          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter,       &
    486488                         status, ierr )   
    487489       ENDIF
     
    499501       IF ( myid == 0 ) THEN
    500502          CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
    501           CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, &
     503          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter,       &
    502504                                                             status, ierr )
    503505       ENDIF
     
    508510
    509511          IF ( dx < remote ) THEN
    510              WRITE( message_string, * ) 'coupling mode "', &
    511                    TRIM( coupling_mode ),                  &
     512             WRITE( message_string, * ) 'coupling mode "',                     &
     513                   TRIM( coupling_mode ),                                      &
    512514           '": dx in Atmosphere is not equal to or not larger then dx in ocean'
    513515             CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
     
    515517
    516518          IF ( (nx_a+1)*dx /= (nx_o+1)*remote )  THEN
    517              WRITE( message_string, * ) 'coupling mode "', &
    518                     TRIM( coupling_mode ), &
     519             WRITE( message_string, * ) 'coupling mode "',                     &
     520                    TRIM( coupling_mode ),                                     &
    519521             '": Domain size in x-direction is not equal in ocean and atmosphere'
    520522             CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
     
    526528       IF ( myid == 0) THEN
    527529          CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )
    528           CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, &
     530          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter,       &
    529531                         status, ierr )
    530532       ENDIF
     
    534536
    535537          IF ( dy < remote )  THEN
    536              WRITE( message_string, * ) 'coupling mode "', &
    537                     TRIM( coupling_mode ), &
     538             WRITE( message_string, * ) 'coupling mode "',                     &
     539                    TRIM( coupling_mode ),                                     &
    538540                 '": dy in Atmosphere is not equal to or not larger then dy in ocean'
    539541             CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
     
    541543
    542544          IF ( (ny_a+1)*dy /= (ny_o+1)*remote )  THEN
    543              WRITE( message_string, * ) 'coupling mode "', &
    544                    TRIM( coupling_mode ), &
     545             WRITE( message_string, * ) 'coupling mode "',                     &
     546                   TRIM( coupling_mode ),                                      &
    545547             '": Domain size in y-direction is not equal in ocean and atmosphere'
    546548             CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )
     
    548550
    549551          IF ( MOD(nx_o+1,nx_a+1) /= 0 )  THEN
    550              WRITE( message_string, * ) 'coupling mode "', &
    551                    TRIM( coupling_mode ), &
     552             WRITE( message_string, * ) 'coupling mode "',                     &
     553                   TRIM( coupling_mode ),                                      &
    552554             '": nx+1 in ocean is not divisible without remainder with nx+1 in', &
    553555             ' atmosphere'
     
    556558
    557559          IF ( MOD(ny_o+1,ny_a+1) /= 0 )  THEN
    558              WRITE( message_string, * ) 'coupling mode "', &
    559                    TRIM( coupling_mode ), &
     560             WRITE( message_string, * ) 'coupling mode "',                     &
     561                   TRIM( coupling_mode ),                                      &
    560562             '": ny+1 in ocean is not divisible without remainder with ny+1 in', &
    561563             ' atmosphere'
     
    565567       ENDIF
    566568#else
    567        WRITE( message_string, * ) 'coupling requires PALM to be called with', &
     569       WRITE( message_string, * ) 'coupling requires PALM to be called with',  &
    568570            ' ''mrun -K parallel'''
    569571       CALL message( 'check_parameters', 'PA0141', 1, 2, 0, 6, 0 )
     
    575577!-- Exchange via intercommunicator
    576578    IF ( coupling_mode == 'atmosphere_to_ocean' .AND. myid == 0 )  THEN
    577        CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter, &
     579       CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter,     &
    578580                      ierr )
    579581    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' .AND. myid == 0)  THEN
    580        CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19, &
     582       CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19,          &
    581583                      comm_inter, status, ierr )
    582584    ENDIF
     
    633635
    634636       CASE DEFAULT
    635           message_string = 'illegal value given for loop_optimization: "' // &
     637          message_string = 'illegal value given for loop_optimization: "' //   &
    636638                           TRIM( loop_optimization ) // '"'
    637639          CALL message( 'check_parameters', 'PA0013', 1, 2, 0, 6, 0 )
     
    659661    IF ( topography /= 'flat' )  THEN
    660662       action = ' '
    661        IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme'     &
     663       IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme'      &
    662664      .AND. scalar_advec /= 'ws-scheme-mono' )  THEN
    663665          WRITE( action, '(A,A)' )  'scalar_advec = ', scalar_advec
    664666       ENDIF
    665        IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' ) &
     667       IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' )&
    666668       THEN
    667669          WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
     
    686688       ENDIF
    687689       IF ( action /= ' ' )  THEN
    688           message_string = 'a non-flat topography does not allow ' // &
     690          message_string = 'a non-flat topography does not allow ' //          &
    689691                           TRIM( action )
    690692          CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
     
    695697!--    is applicable. If this is not possible, abort.
    696698       IF ( TRIM( topography_grid_convention ) == ' ' )  THEN
    697           IF ( TRIM( topography ) /= 'single_building' .AND.  &
    698                TRIM( topography ) /= 'single_street_canyon' .AND.  &
     699          IF ( TRIM( topography ) /= 'single_building' .AND.                   &
     700               TRIM( topography ) /= 'single_street_canyon' .AND.              &
    699701               TRIM( topography ) /= 'read_from_file' )  THEN
    700702!--          The default value is not applicable here, because it is only valid
    701703!--          for the two standard cases 'single_building' and 'read_from_file'
    702704!--          defined in init_grid.
    703              WRITE( message_string, * )  &
    704                   'The value for "topography_grid_convention" ',  &
    705                   'is not set. Its default value is & only valid for ',  &
    706                   '"topography" = ''single_building'', ',  &
    707                   '''single_street_canyon'' & or ''read_from_file''.',  &
     705             WRITE( message_string, * )                                        &
     706                  'The value for "topography_grid_convention" ',               &
     707                  'is not set. Its default value is & only valid for ',        &
     708                  '"topography" = ''single_building'', ',                      &
     709                  '''single_street_canyon'' & or ''read_from_file''.',         &
    708710                  ' & Choose ''cell_edge'' or ''cell_center''.'
    709711             CALL message( 'user_check_parameters', 'PA0239', 1, 2, 0, 6, 0 )
     
    711713!--          The default value is applicable here.
    712714!--          Set convention according to topography.
    713              IF ( TRIM( topography ) == 'single_building' .OR.  &
     715             IF ( TRIM( topography ) == 'single_building' .OR.                 &
    714716                  TRIM( topography ) == 'single_street_canyon' )  THEN
    715717                topography_grid_convention = 'cell_edge'
     
    718720             ENDIF
    719721          ENDIF
    720        ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND.  &
     722       ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND.        &
    721723                TRIM( topography_grid_convention ) /= 'cell_center' )  THEN
    722           WRITE( message_string, * )  &
    723                'The value for "topography_grid_convention" is ', &
     724          WRITE( message_string, * )                                           &
     725               'The value for "topography_grid_convention" is ',               &
    724726               'not recognized. & Choose ''cell_edge'' or ''cell_center''.'
    725727          CALL message( 'user_check_parameters', 'PA0240', 1, 2, 0, 6, 0 )
     
    738740       ENDIF
    739741
    740     ELSEIF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.  &
     742    ELSEIF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.                       &
    741743             TRIM( coupling_char ) == '_O' )  THEN
    742744
     
    744746!--    Check whether an (uncoupled) atmospheric run has been declared as an
    745747!--    ocean run (this setting is done via mrun-option -y)
    746 
    747        message_string = 'ocean = .F. does not allow coupling_char = "' // &
     748       message_string = 'ocean = .F. does not allow coupling_char = "' //      &
    748749                        TRIM( coupling_char ) // '" set by mrun-option "-y"'
    749750       CALL message( 'check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
     
    777778          gamma_mg = 1
    778779       ELSE
    779           message_string = 'unknown multigrid cycle: cycle_mg = "' // &
     780          message_string = 'unknown multigrid cycle: cycle_mg = "' //          &
    780781                           TRIM( cycle_mg ) // '"'
    781782          CALL message( 'check_parameters', 'PA0020', 1, 2, 0, 6, 0 )
     
    783784    ENDIF
    784785
    785     IF ( fft_method /= 'singleton-algorithm'  .AND.  &
    786          fft_method /= 'temperton-algorithm'  .AND.  &
    787          fft_method /= 'fftw'                 .AND.  &
     786    IF ( fft_method /= 'singleton-algorithm'  .AND.                            &
     787         fft_method /= 'temperton-algorithm'  .AND.                            &
     788         fft_method /= 'fftw'                 .AND.                            &
    788789         fft_method /= 'system-specific' )  THEN
    789        message_string = 'unknown fft-algorithm: fft_method = "' // &
     790       message_string = 'unknown fft-algorithm: fft_method = "' //             &
    790791                        TRIM( fft_method ) // '"'
    791792       CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
    792793    ENDIF
    793794   
    794     IF( momentum_advec == 'ws-scheme' .AND. &
     795    IF( momentum_advec == 'ws-scheme' .AND.                                    &
    795796        .NOT. call_psolver_at_all_substeps  ) THEN
    796         message_string = 'psolver must be called at each RK3 substep when "'//&
     797        message_string = 'psolver must be called at each RK3 substep when "'// &
    797798                      TRIM(momentum_advec) // ' "is used for momentum_advec'
    798799        CALL message( 'check_parameters', 'PA0344', 1, 2, 0, 6, 0 )
     
    802803    IF ( momentum_advec /= 'pw-scheme'  .AND.  momentum_advec /= 'ws-scheme' ) &
    803804    THEN
    804        message_string = 'unknown advection scheme: momentum_advec = "' // &
     805       message_string = 'unknown advection scheme: momentum_advec = "' //      &
    805806                        TRIM( momentum_advec ) // '"'
    806807       CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
     
    811812                   timestep_scheme == 'runge-kutta-2' ) )                      &
    812813    THEN
    813        message_string = 'momentum_advec or scalar_advec = "' &
     814       message_string = 'momentum_advec or scalar_advec = "'                   &
    814815         // TRIM( momentum_advec ) // '" is not allowed with timestep_scheme = "' // &
    815816         TRIM( timestep_scheme ) // '"'
     
    819820         scalar_advec /= 'ws-scheme-mono' .AND. scalar_advec /= 'bc-scheme' )  &
    820821    THEN
    821        message_string = 'unknown advection scheme: scalar_advec = "' // &
     822       message_string = 'unknown advection scheme: scalar_advec = "' //        &
    822823                        TRIM( scalar_advec ) // '"'
    823824       CALL message( 'check_parameters', 'PA0024', 1, 2, 0, 6, 0 )
     
    825826    IF ( scalar_advec == 'bc-scheme'  .AND.  loop_optimization == 'cache' )    &
    826827    THEN
    827        message_string = 'advection_scheme scalar_advec = "' &
     828       message_string = 'advection_scheme scalar_advec = "'                    &
    828829         // TRIM( scalar_advec ) // '" not implemented for & loop_optimization = "' // &
    829830         TRIM( loop_optimization ) // '"'
     
    851852!-- Set LOGICAL switches to enhance performance
    852853    IF ( momentum_advec == 'ws-scheme' )       ws_scheme_mom = .TRUE.
    853     IF ( scalar_advec   == 'ws-scheme' .OR.                                   &
     854    IF ( scalar_advec   == 'ws-scheme' .OR.                                    &
    854855         scalar_advec   == 'ws-scheme-mono' )  ws_scheme_sca = .TRUE.
    855856    IF ( scalar_advec   == 'ws-scheme-mono' )  monotonic_adjustment = .TRUE.
     
    908909    IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
    909910
    910     IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
     911    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
    911912         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
    912913!
    913914!--    No restart run: several initialising actions are possible
    914915       action = initializing_actions
    915        DO WHILE ( TRIM( action ) /= '' )
     916       DO  WHILE ( TRIM( action ) /= '' )
    916917          position = INDEX( action, ' ' )
    917918          SELECT CASE ( action(1:position-1) )
    918919
    919              CASE ( 'set_constant_profiles', 'set_1d-model_profiles', &
     920             CASE ( 'set_constant_profiles', 'set_1d-model_profiles',          &
    920921                    'by_user', 'initialize_vortex',     'initialize_ptanom' )
    921922                action = action(position+1:)
    922923
    923924             CASE DEFAULT
    924                 message_string = 'initializing_action = "' // &
     925                message_string = 'initializing_action = "' //                  &
    925926                                 TRIM( action ) // '" unkown or not allowed'
    926927                CALL message( 'check_parameters', 'PA0030', 1, 2, 0, 6, 0 )
     
    930931    ENDIF
    931932
    932     IF ( TRIM( initializing_actions ) == 'initialize_vortex' .AND. &
     933    IF ( TRIM( initializing_actions ) == 'initialize_vortex'  .AND.            &
    933934         conserve_volume_flow ) THEN
    934          message_string = 'initializing_actions = "initialize_vortex"' // &
     935         message_string = 'initializing_actions = "initialize_vortex"' //      &
    935936                        ' ist not allowed with conserve_volume_flow = .T.'
    936937       CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
     
    938939
    939940
    940     IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND. &
     941    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.    &
    941942         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
    942        message_string = 'initializing_actions = "set_constant_profiles"' // &
    943                         ' and "set_1d-model_profiles" are not allowed ' //  &
     943       message_string = 'initializing_actions = "set_constant_profiles"' //    &
     944                        ' and "set_1d-model_profiles" are not allowed ' //     &
    944945                        'simultaneously'
    945946       CALL message( 'check_parameters', 'PA0031', 1, 2, 0, 6, 0 )
    946947    ENDIF
    947948
    948     IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND. &
     949    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.    &
    949950         INDEX( initializing_actions, 'by_user' ) /= 0 )  THEN
    950        message_string = 'initializing_actions = "set_constant_profiles"' // &
     951       message_string = 'initializing_actions = "set_constant_profiles"' //    &
    951952                        ' and "by_user" are not allowed simultaneously'
    952953       CALL message( 'check_parameters', 'PA0032', 1, 2, 0, 6, 0 )
    953954    ENDIF
    954955
    955     IF ( INDEX( initializing_actions, 'by_user' ) /= 0  .AND. &
     956    IF ( INDEX( initializing_actions, 'by_user' ) /= 0  .AND.                  &
    956957         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
    957        message_string = 'initializing_actions = "by_user" and ' // &
     958       message_string = 'initializing_actions = "by_user" and ' //             &
    958959                        '"set_1d-model_profiles" are not allowed simultaneously'
    959960       CALL message( 'check_parameters', 'PA0033', 1, 2, 0, 6, 0 )
    960961    ENDIF
    961962
    962     IF ( cloud_physics  .AND.  .NOT. humidity )  THEN
    963        WRITE( message_string, * ) 'cloud_physics = ', cloud_physics, ' is ', &
     963    IF ( cloud_physics  .AND.  .NOT.  humidity )  THEN
     964       WRITE( message_string, * ) 'cloud_physics = ', cloud_physics, ' is ',   &
    964965              'not allowed with humidity = ', humidity
    965966       CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 )
     
    967968
    968969    IF ( precipitation  .AND.  .NOT.  cloud_physics )  THEN
    969        WRITE( message_string, * ) 'precipitation = ', precipitation, ' is ', &
     970       WRITE( message_string, * ) 'precipitation = ', precipitation, ' is ',   &
    970971              'not allowed with cloud_physics = ', cloud_physics
    971972       CALL message( 'check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
     
    973974
    974975    IF ( humidity  .AND.  sloping_surface )  THEN
    975        message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' // &
     976       message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' //   &
    976977                        'are not allowed simultaneously'
    977978       CALL message( 'check_parameters', 'PA0036', 1, 2, 0, 6, 0 )
     
    979980
    980981    IF ( passive_scalar  .AND.  humidity )  THEN
    981        message_string = 'humidity = .TRUE. and passive_scalar = .TRUE. ' // &
     982       message_string = 'humidity = .TRUE. and passive_scalar = .TRUE. ' //    &
    982983                        'is not allowed simultaneously'
    983984       CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )
     
    10321033!--    calculated from the temperature/humidity gradients in the land surface
    10331034!--    model
    1034        IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' )  THEN
     1035       IF ( bc_pt_b == 'neumann'  .OR. bc_q_b == 'neumann' )  THEN
    10351036          message_string = 'lsm requires setting of'//                         &
    10361037                           'bc_pt_b = "dirichlet" and '//                      &
     
    10391040       ENDIF
    10401041
    1041        IF ( .NOT. constant_flux_layer )  THEN
     1042       IF (  .NOT. constant_flux_layer )  THEN
    10421043          message_string = 'lsm requires '//                                   &
    10431044                           'constant_flux_layer = .T.'
     
    10461047
    10471048       IF ( topography /= 'flat' )  THEN
    1048           message_string = 'lsm cannot be used ' //  &
     1049          message_string = 'lsm cannot be used ' //                            &
    10491050                           'in combination with  topography /= "flat"'
    10501051          CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
    10511052       ENDIF
    10521053
     1054       IF ( ( veg_type == 14  .OR.  veg_type == 15 ) .AND.                       &
     1055              most_method == 'lookup' )  THEN
     1056           WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ',     &
     1057                                      'allowed in combination with ',          &
     1058                                      'most_method = ', most_method
     1059          CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
     1060       ENDIF
     1061
    10531062       IF ( veg_type == 0 )  THEN
    1054           IF ( SUM(root_fraction) /= 1.0_wp)  THEN
     1063          IF ( SUM( root_fraction ) /= 1.0_wp )  THEN
    10551064             message_string = 'veg_type = 0 (user_defined)'//                  &
    10561065                              'requires setting of root_fraction(0:3)'//       &
     
    10591068          ENDIF
    10601069 
    1061           IF ( min_canopy_resistance == 9999999.9_wp)  THEN
     1070          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
    10621071             message_string = 'veg_type = 0 (user defined)'//                  &
    10631072                              'requires setting of min_canopy_resistance'//    &
     
    10661075          ENDIF
    10671076
    1068           IF ( leaf_area_index == 9999999.9_wp)  THEN
     1077          IF ( leaf_area_index == 9999999.9_wp )  THEN
    10691078             message_string = 'veg_type = 0 (user_defined)'//                  &
    10701079                              'requires setting of leaf_area_index'//          &
     
    10731082          ENDIF
    10741083
    1075           IF ( vegetation_coverage == 9999999.9_wp)  THEN
     1084          IF ( vegetation_coverage == 9999999.9_wp )  THEN
    10761085             message_string = 'veg_type = 0 (user_defined)'//                  &
    10771086                              'requires setting of vegetation_coverage'//      &
     
    10871096          ENDIF
    10881097
    1089           IF ( lambda_surface_stable == 9999999.9_wp)  THEN
     1098          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
    10901099             message_string = 'veg_type = 0 (user_defined)'//                  &
    10911100                              'requires setting of lambda_surface_stable'//    &
     
    10941103          ENDIF
    10951104
    1096           IF ( lambda_surface_unstable == 9999999.9_wp)  THEN
     1105          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
    10971106             message_string = 'veg_type = 0 (user_defined)'//                  &
    10981107                              'requires setting of lambda_surface_unstable'//  &
     
    11011110          ENDIF
    11021111
    1103           IF ( f_shortwave_incoming == 9999999.9_wp)  THEN
     1112          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
    11041113             message_string = 'veg_type = 0 (user_defined)'//                  &
    11051114                              'requires setting of f_shortwave_incoming'//     &
     
    11081117          ENDIF
    11091118
    1110           IF ( z0_eb == 9999999.9_wp)  THEN
     1119          IF ( z0_eb == 9999999.9_wp )  THEN
    11111120             message_string = 'veg_type = 0 (user_defined)'//                  &
    11121121                              'requires setting of z0_eb'//                   &
     
    11151124          ENDIF
    11161125
    1117           IF ( z0h_eb == 9999999.9_wp)  THEN
     1126          IF ( z0h_eb == 9999999.9_wp )  THEN
    11181127             message_string = 'veg_type = 0 (user_defined)'//                  &
    11191128                              'requires setting of z0h_eb'//                  &
     
    11271136       IF ( soil_type == 0 )  THEN
    11281137
    1129           IF ( alpha_vangenuchten == 9999999.9_wp)  THEN
     1138          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
    11301139             message_string = 'soil_type = 0 (user_defined)'//                 &
    11311140                              'requires setting of alpha_vangenuchten'//       &
     
    11341143          ENDIF
    11351144
    1136           IF ( l_vangenuchten == 9999999.9_wp)  THEN
     1145          IF ( l_vangenuchten == 9999999.9_wp )  THEN
    11371146             message_string = 'soil_type = 0 (user_defined)'//                 &
    11381147                              'requires setting of l_vangenuchten'//           &
     
    11411150          ENDIF
    11421151
    1143           IF ( n_vangenuchten == 9999999.9_wp)  THEN
     1152          IF ( n_vangenuchten == 9999999.9_wp )  THEN
    11441153             message_string = 'soil_type = 0 (user_defined)'//                 &
    11451154                              'requires setting of n_vangenuchten'//           &
     
    11481157          ENDIF
    11491158
    1150           IF ( hydraulic_conductivity == 9999999.9_wp)  THEN
     1159          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
    11511160             message_string = 'soil_type = 0 (user_defined)'//                 &
    11521161                              'requires setting of hydraulic_conductivity'//   &
     
    11551164          ENDIF
    11561165
    1157           IF ( saturation_moisture == 9999999.9_wp)  THEN
     1166          IF ( saturation_moisture == 9999999.9_wp )  THEN
    11581167             message_string = 'soil_type = 0 (user_defined)'//                 &
    11591168                              'requires setting of saturation_moisture'//      &
     
    11621171          ENDIF
    11631172
    1164           IF ( field_capacity == 9999999.9_wp)  THEN
     1173          IF ( field_capacity == 9999999.9_wp )  THEN
    11651174             message_string = 'soil_type = 0 (user_defined)'//                 &
    11661175                              'requires setting of field_capacity'//           &
     
    11691178          ENDIF
    11701179
    1171           IF ( wilting_point == 9999999.9_wp)  THEN
     1180          IF ( wilting_point == 9999999.9_wp )  THEN
    11721181             message_string = 'soil_type = 0 (user_defined)'//                 &
    11731182                              'requires setting of wilting_point'//            &
     
    11761185          ENDIF
    11771186
    1178           IF ( residual_moisture == 9999999.9_wp)  THEN
     1187          IF ( residual_moisture == 9999999.9_wp )  THEN
    11791188             message_string = 'soil_type = 0 (user_defined)'//                 &
    11801189                              'requires setting of residual_moisture'//        &
     
    11851194       ENDIF
    11861195
    1187        IF ( .NOT. radiation )  THEN
     1196       IF (  .NOT. radiation )  THEN
    11881197          message_string = 'lsm requires '//                                   &
    11891198                           'radiation = .T.'
     
    11941203
    11951204    IF ( radiation )  THEN
    1196        IF ( radiation_scheme /= 'constant'  .AND.                              &
    1197             radiation_scheme /= 'clear-sky' .AND.                              &
     1205       IF ( radiation_scheme /= 'constant'   .AND.                             &
     1206            radiation_scheme /= 'clear-sky'  .AND.                             &
    11981207            radiation_scheme /= 'rrtmg' )  THEN
    11991208          message_string = 'unknown radiation_scheme = '//                     &
     
    12151224
    12161225       ENDIF
    1217        IF ( albedo_type == 0 .AND. albedo == 9999999.9_wp .AND.                &
     1226       IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.             &
    12181227            radiation_scheme == 'clear-sky')  THEN
    12191228          message_string = 'radiation_scheme = "clear-sky" in combination' //  &
     
    12221231          CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
    12231232       ENDIF
    1224        IF ( albedo_type == 0 .AND. radiation_scheme == 'rrtmg' .AND.           &
     1233       IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.        &
    12251234          (    albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
    12261235          .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp&
     
    12351244       ENDIF
    12361245       IF ( topography /= 'flat' )  THEN
    1237           message_string = 'radiation scheme cannot be used ' //  &
     1246          message_string = 'radiation scheme cannot be used ' //               &
    12381247                           'in combination with  topography /= "flat"'
    12391248          CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
     
    12441253                 loop_optimization == 'vector' )                               &
    12451254         .AND.  cloud_physics  .AND.  icloud_scheme == 0 )  THEN
    1246        message_string = 'cloud_scheme = seifert_beheng requires ' // &
     1255       message_string = 'cloud_scheme = seifert_beheng requires ' //           &
    12471256                        'loop_optimization = "cache" or "vector"'
    12481257       CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
     
    12701279       gradient = 0.0_wp
    12711280
    1272        IF ( .NOT. ocean )  THEN
     1281       IF (  .NOT. ocean )  THEN
    12731282
    12741283          ug_vertical_gradient_level_ind(1) = 0
    12751284          ug(0) = ug_surface
    12761285          DO  k = 1, nzt+1
    1277              IF ( i < 11 ) THEN
    1278                 IF ( ug_vertical_gradient_level(i) < zu(k)  .AND. &
     1286             IF ( i < 11 )  THEN
     1287                IF ( ug_vertical_gradient_level(i) < zu(k)  .AND.              &
    12791288                     ug_vertical_gradient_level(i) >= 0.0_wp )  THEN
    12801289                   gradient = ug_vertical_gradient(i) / 100.0_wp
     
    12991308          ug(nzt+1) = ug_surface
    13001309          DO  k = nzt, nzb, -1
    1301              IF ( i < 11 ) THEN
    1302                 IF ( ug_vertical_gradient_level(i) > zu(k)  .AND. &
     1310             IF ( i < 11 )  THEN
     1311                IF ( ug_vertical_gradient_level(i) > zu(k)  .AND.              &
    13031312                     ug_vertical_gradient_level(i) <= 0.0_wp )  THEN
    13041313                   gradient = ug_vertical_gradient(i) / 100.0_wp
     
    13341343       gradient = 0.0_wp
    13351344
    1336        IF ( .NOT. ocean )  THEN
     1345       IF (  .NOT. ocean )  THEN
    13371346
    13381347          vg_vertical_gradient_level_ind(1) = 0
    13391348          vg(0) = vg_surface
    13401349          DO  k = 1, nzt+1
    1341              IF ( i < 11 ) THEN
    1342                 IF ( vg_vertical_gradient_level(i) < zu(k)  .AND. &
     1350             IF ( i < 11 )  THEN
     1351                IF ( vg_vertical_gradient_level(i) < zu(k)  .AND.              &
    13431352                     vg_vertical_gradient_level(i) >= 0.0_wp )  THEN
    13441353                   gradient = vg_vertical_gradient(i) / 100.0_wp
     
    13631372          vg(nzt+1) = vg_surface
    13641373          DO  k = nzt, nzb, -1
    1365              IF ( i < 11 ) THEN
    1366                 IF ( vg_vertical_gradient_level(i) > zu(k)  .AND. &
     1374             IF ( i < 11 )  THEN
     1375                IF ( vg_vertical_gradient_level(i) > zu(k)  .AND.              &
    13671376                     vg_vertical_gradient_level(i) <= 0.0_wp )  THEN
    13681377                   gradient = vg_vertical_gradient(i) / 100.0_wp
     
    14151424
    14161425             IF ( kk < 100 )  THEN
    1417                 DO WHILE ( uv_heights(kk+1) <= zu(k) )
     1426                DO  WHILE ( uv_heights(kk+1) <= zu(k) )
    14181427                   kk = kk + 1
    14191428                   IF ( kk == 100 )  EXIT
     
    14211430             ENDIF
    14221431
    1423              IF ( kk < 100 .AND. uv_heights(kk+1) /= 9999999.9_wp )  THEN
     1432             IF ( kk < 100  .AND. uv_heights(kk+1) /= 9999999.9_wp )  THEN
    14241433                u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
    14251434                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
     
    14441453!
    14451454!--    Compute initial temperature profile using the given temperature gradients
    1446        IF ( .NOT. neutral )  THEN
     1455       IF (  .NOT. neutral )  THEN
    14471456
    14481457          i = 1
    14491458          gradient = 0.0_wp
    14501459
    1451           IF ( .NOT. ocean )  THEN
     1460          IF (  .NOT. ocean )  THEN
    14521461
    14531462             pt_vertical_gradient_level_ind(1) = 0
    14541463             DO  k = 1, nzt+1
    1455                 IF ( i < 11 ) THEN
    1456                    IF ( pt_vertical_gradient_level(i) < zu(k)  .AND. &
     1464                IF ( i < 11 )  THEN
     1465                   IF ( pt_vertical_gradient_level(i) < zu(k)  .AND.           &
    14571466                        pt_vertical_gradient_level(i) >= 0.0_wp )  THEN
    14581467                      gradient = pt_vertical_gradient(i) / 100.0_wp
     
    14761485             pt_vertical_gradient_level_ind(1) = nzt+1
    14771486             DO  k = nzt, 0, -1
    1478                 IF ( i < 11 ) THEN
    1479                    IF ( pt_vertical_gradient_level(i) > zu(k)  .AND. &
     1487                IF ( i < 11 )  THEN
     1488                   IF ( pt_vertical_gradient_level(i) > zu(k)  .AND.           &
    14801489                        pt_vertical_gradient_level(i) <= 0.0_wp )  THEN
    14811490                      gradient = pt_vertical_gradient(i) / 100.0_wp
     
    15331542          q_vertical_gradient_level_ind(1) = 0
    15341543          DO  k = 1, nzt+1
    1535              IF ( i < 11 ) THEN
    1536                 IF ( q_vertical_gradient_level(i) < zu(k)  .AND. &
     1544             IF ( i < 11 )  THEN
     1545                IF ( q_vertical_gradient_level(i) < zu(k)  .AND.               &
    15371546                     q_vertical_gradient_level(i) >= 0.0_wp )  THEN
    15381547                   gradient = q_vertical_gradient(i) / 100.0_wp
     
    15791588          sa_vertical_gradient_level_ind(1) = nzt+1
    15801589          DO  k = nzt, 0, -1
    1581              IF ( i < 11 ) THEN
    1582                 IF ( sa_vertical_gradient_level(i) > zu(k)  .AND. &
     1590             IF ( i < 11 )  THEN
     1591                IF ( sa_vertical_gradient_level(i) > zu(k)  .AND.              &
    15831592                     sa_vertical_gradient_level(i) <= 0.0_wp )  THEN
    15841593                   gradient = sa_vertical_gradient(i) / 100.0_wp
     
    16061615!
    16071616!-- Check if the control parameter use_subsidence_tendencies is used correctly
    1608     IF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_subsidence )  THEN
    1609        message_string = 'The usage of use_subsidence_tendencies ' // &
     1617    IF ( use_subsidence_tendencies  .AND.  .NOT.  large_scale_subsidence )  THEN
     1618       message_string = 'The usage of use_subsidence_tendencies ' //           &
    16101619                            'requires large_scale_subsidence = .T..'
    16111620       CALL message( 'check_parameters', 'PA0396', 1, 2, 0, 6, 0 )
    16121621    ELSEIF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_forcing )  THEN
    1613        message_string = 'The usage of use_subsidence_tendencies ' // &
     1622       message_string = 'The usage of use_subsidence_tendencies ' //           &
    16141623                            'requires large_scale_forcing = .T..'
    16151624       CALL message( 'check_parameters', 'PA0397', 1, 2, 0, 6, 0 )
     
    16191628!-- Initialize large scale subsidence if required
    16201629    If ( large_scale_subsidence )  THEN
    1621        IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp .AND. &
    1622                                      .NOT. large_scale_forcing )  THEN
     1630       IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp  .AND.            &
     1631                                     .NOT.  large_scale_forcing )  THEN
    16231632          CALL init_w_subsidence
    16241633       ENDIF
     
    16271636!--    are read in from file LSF_DATA
    16281637
    1629        IF ( subs_vertical_gradient_level(1) == -9999999.9_wp .AND. &
    1630                                      .NOT. large_scale_forcing )  THEN
    1631           message_string = 'There is no default large scale vertical ' // &
    1632                            'velocity profile set. Specify the subsidence ' // &
    1633                            'velocity profile via subs_vertical_gradient and ' // &
    1634                            'subs_vertical_gradient_level.'
     1638       IF ( subs_vertical_gradient_level(1) == -9999999.9_wp  .AND.            &
     1639            .NOT. large_scale_forcing )  THEN
     1640          message_string = 'There is no default large scale vertical ' //      &
     1641                           'velocity profile set. Specify the subsidence ' //  &
     1642                           'velocity profile via subs_vertical_gradient ' //  &
     1643                           'and subs_vertical_gradient_level.'
    16351644          CALL message( 'check_parameters', 'PA0380', 1, 2, 0, 6, 0 )
    16361645       ENDIF
    16371646    ELSE
    16381647        IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp )  THEN
    1639            message_string = 'Enable usage of large scale subsidence by ' // &
     1648           message_string = 'Enable usage of large scale subsidence by ' //    &
    16401649                            'setting large_scale_subsidence = .T..'
    16411650          CALL message( 'check_parameters', 'PA0381', 1, 2, 0, 6, 0 )
     
    16591668       vpt_reference = pt_reference * ( 1.0_wp + 0.61_wp * q_surface )
    16601669    ELSE
    1661        message_string = 'illegal value for reference_state: "' // &
     1670       message_string = 'illegal value for reference_state: "' //              &
    16621671                        TRIM( reference_state ) // '"'
    16631672       CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 )
     
    16861695    IF ( alpha_surface /= 0.0_wp )  THEN
    16871696       IF ( ABS( alpha_surface ) > 90.0_wp )  THEN
    1688           WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface, &
     1697          WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface,   &
    16891698                                     ' ) must be < 90.0'
    16901699          CALL message( 'check_parameters', 'PA0043', 1, 2, 0, 6, 0 )
     
    17161725          ENDIF
    17171726       ELSE
    1718           WRITE( message_string, * ) 'cfl_factor = ', cfl_factor, &
     1727          WRITE( message_string, * ) 'cfl_factor = ', cfl_factor,              &
    17191728                 ' out of range & 0.0 < cfl_factor <= 1.0 is required'
    17201729          CALL message( 'check_parameters', 'PA0045', 1, 2, 0, 6, 0 )
     
    17401749!-- Set wind speed in the Galilei-transformed system
    17411750    IF ( galilei_transformation )  THEN
    1742        IF ( use_ug_for_galilei_tr .AND.                     &
    1743             ug_vertical_gradient_level(1) == 0.0_wp  .AND.  &
    1744             ug_vertical_gradient(1) == 0.0_wp  .AND.        &
    1745             vg_vertical_gradient_level(1) == 0.0_wp  .AND.  &
     1751       IF ( use_ug_for_galilei_tr                    .AND.                     &
     1752            ug_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
     1753            ug_vertical_gradient(1) == 0.0_wp        .AND.                     &
     1754            vg_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
    17461755            vg_vertical_gradient(1) == 0.0_wp )  THEN
    17471756          u_gtrans = ug_surface * 0.6_wp
    17481757          v_gtrans = vg_surface * 0.6_wp
    1749        ELSEIF ( use_ug_for_galilei_tr  .AND.                     &
    1750                 ( ug_vertical_gradient_level(1) /= 0.0_wp  .OR.  &
     1758       ELSEIF ( use_ug_for_galilei_tr  .AND.                                   &
     1759                ( ug_vertical_gradient_level(1) /= 0.0_wp  .OR.                &
    17511760                ug_vertical_gradient(1) /= 0.0_wp ) )  THEN
    1752           message_string = 'baroclinity (ug) not allowed simultaneously' // &
     1761          message_string = 'baroclinity (ug) not allowed simultaneously' //    &
    17531762                           ' with galilei transformation'
    17541763          CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 )
    1755        ELSEIF ( use_ug_for_galilei_tr  .AND.                     &
    1756                 ( vg_vertical_gradient_level(1) /= 0.0_wp  .OR.  &
     1764       ELSEIF ( use_ug_for_galilei_tr  .AND.                                   &
     1765                ( vg_vertical_gradient_level(1) /= 0.0_wp  .OR.                &
    17571766                vg_vertical_gradient(1) /= 0.0_wp ) )  THEN
    1758           message_string = 'baroclinity (vg) not allowed simultaneously' // &
     1767          message_string = 'baroclinity (vg) not allowed simultaneously' //    &
    17591768                           ' with galilei transformation'
    17601769          CALL message( 'check_parameters', 'PA0047', 1, 2, 0, 6, 0 )
    17611770       ELSE
    1762           message_string = 'variable translation speed used for galilei-' // &
    1763              'transformation, which may cause & instabilities in stably ' // &
     1771          message_string = 'variable translation speed used for galilei-' //   &
     1772             'transformation, which may cause & instabilities in stably ' //   &
    17641773             'stratified regions'
    17651774          CALL message( 'check_parameters', 'PA0048', 0, 1, 0, 6, 0 )
     
    17821791    IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
    17831792       IF ( psolver(1:9) /= 'multigrid' )  THEN
    1784           message_string = 'non-cyclic lateral boundaries do not allow ' // &
     1793          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
    17851794                           'psolver = "' // TRIM( psolver ) // '"'
    17861795          CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
    17871796       ENDIF
    1788        IF ( momentum_advec /= 'pw-scheme' .AND.                               &
    1789           ( momentum_advec /= 'ws-scheme' .AND.                               &
    1790             momentum_advec /= 'ws-scheme-mono' )                              &
     1797       IF ( momentum_advec /= 'pw-scheme'  .AND.                               &
     1798          ( momentum_advec /= 'ws-scheme'  .AND.                               &
     1799            momentum_advec /= 'ws-scheme-mono' )                               &
    17911800          )  THEN
    17921801
    1793           message_string = 'non-cyclic lateral boundaries do not allow ' // &
     1802          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
    17941803                           'momentum_advec = "' // TRIM( momentum_advec ) // '"'
    17951804          CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
    17961805       ENDIF
    1797        IF ( scalar_advec /= 'pw-scheme' .AND.                                  &
    1798           ( scalar_advec /= 'ws-scheme' .AND.                                  &
     1806       IF ( scalar_advec /= 'pw-scheme'  .AND.                                 &
     1807          ( scalar_advec /= 'ws-scheme'  .AND.                                 &
    17991808            scalar_advec /= 'ws-scheme-mono' )                                 &
    18001809          )  THEN
    1801           message_string = 'non-cyclic lateral boundaries do not allow ' // &
     1810          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
    18021811                           'scalar_advec = "' // TRIM( scalar_advec ) // '"'
    18031812          CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 )
    18041813       ENDIF
    18051814       IF ( galilei_transformation )  THEN
    1806           message_string = 'non-cyclic lateral boundaries do not allow ' // &
     1815          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
    18071816                           'galilei_transformation = .T.'
    18081817          CALL message( 'check_parameters', 'PA0054', 1, 2, 0, 6, 0 )
     
    18191828          bc_e_b = 'neumann'
    18201829          ibc_e_b = 1
    1821           message_string = 'boundary condition bc_e_b changed to "' // &
     1830          message_string = 'boundary condition bc_e_b changed to "' //         &
    18221831                           TRIM( bc_e_b ) // '"'
    18231832          CALL message( 'check_parameters', 'PA0057', 0, 1, 0, 6, 0 )
    18241833       ENDIF
    18251834    ELSE
    1826        message_string = 'unknown boundary condition: bc_e_b = "' // &
     1835       message_string = 'unknown boundary condition: bc_e_b = "' //            &
    18271836                        TRIM( bc_e_b ) // '"'
    18281837       CALL message( 'check_parameters', 'PA0058', 1, 2, 0, 6, 0 )
     
    18361845       ibc_p_b = 1
    18371846    ELSE
    1838        message_string = 'unknown boundary condition: bc_p_b = "' // &
     1847       message_string = 'unknown boundary condition: bc_p_b = "' //            &
    18391848                        TRIM( bc_p_b ) // '"'
    18401849       CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 )
     
    18471856       ibc_p_t = 1
    18481857    ELSE
    1849        message_string = 'unknown boundary condition: bc_p_t = "' // &
     1858       message_string = 'unknown boundary condition: bc_p_t = "' //            &
    18501859                        TRIM( bc_p_t ) // '"'
    18511860       CALL message( 'check_parameters', 'PA0061', 1, 2, 0, 6, 0 )
     
    18621871          ibc_pt_b = 1
    18631872       ELSE
    1864           message_string = 'unknown boundary condition: bc_pt_b = "' // &
     1873          message_string = 'unknown boundary condition: bc_pt_b = "' //        &
    18651874                           TRIM( bc_pt_b ) // '"'
    18661875          CALL message( 'check_parameters', 'PA0062', 1, 2, 0, 6, 0 )
     
    18771886       ibc_pt_t = 3
    18781887    ELSE
    1879        message_string = 'unknown boundary condition: bc_pt_t = "' // &
     1888       message_string = 'unknown boundary condition: bc_pt_t = "' //           &
    18801889                        TRIM( bc_pt_t ) // '"'
    18811890       CALL message( 'check_parameters', 'PA0063', 1, 2, 0, 6, 0 )
     
    18891898          ELSEIF ( ibc_pt_b == 1 )  THEN
    18901899             constant_heatflux = .TRUE.
    1891              IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND.    &
    1892                   .NOT. land_surface )  THEN
     1900             IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.   &
     1901                  .NOT.  land_surface )  THEN
    18931902                surface_heatflux = shf_surf(1)
    18941903             ELSE
     
    18991908    ELSE
    19001909        constant_heatflux = .TRUE.
    1901         IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND.         &
    1902              large_scale_forcing .AND. .NOT. land_surface )  THEN
     1910        IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.        &
     1911             large_scale_forcing  .AND.  .NOT. land_surface )  THEN
    19031912           surface_heatflux = shf_surf(1)
    19041913        ENDIF
     
    19091918    IF ( neutral )  THEN
    19101919
    1911        IF ( surface_heatflux /= 0.0_wp  .AND. surface_heatflux /= 9999999.9_wp ) &
     1920       IF ( surface_heatflux /= 0.0_wp  .AND.                                  &
     1921            surface_heatflux /= 9999999.9_wp )  THEN
     1922          message_string = 'heatflux must not be set for pure neutral flow'
     1923          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
     1924       ENDIF
     1925
     1926       IF ( top_heatflux /= 0.0_wp  .AND.  top_heatflux /= 9999999.9_wp )      &
    19121927       THEN
    19131928          message_string = 'heatflux must not be set for pure neutral flow'
     
    19151930       ENDIF
    19161931
    1917        IF ( top_heatflux /= 0.0_wp  .AND.  top_heatflux /= 9999999.9_wp ) &
    1918        THEN
    1919           message_string = 'heatflux must not be set for pure neutral flow'
    1920           CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
    1921        ENDIF
    1922 
    1923     ENDIF
    1924 
    1925     IF ( top_momentumflux_u /= 9999999.9_wp  .AND.  &
     1932    ENDIF
     1933
     1934    IF ( top_momentumflux_u /= 9999999.9_wp  .AND.                             &
    19261935         top_momentumflux_v /= 9999999.9_wp )  THEN
    19271936       constant_top_momentumflux = .TRUE.
    1928     ELSEIF (  .NOT. ( top_momentumflux_u == 9999999.9_wp  .AND.  &
     1937    ELSEIF (  .NOT. ( top_momentumflux_u == 9999999.9_wp  .AND.                &
    19291938           top_momentumflux_v == 9999999.9_wp ) )  THEN
    1930        message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' // &
     1939       message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' //  &
    19311940                        'must be set'
    19321941       CALL message( 'check_parameters', 'PA0064', 1, 2, 0, 6, 0 )
     
    19371946!-- temperature. In this case specification of a constant heat flux is
    19381947!-- forbidden.
    1939     IF ( ibc_pt_b == 0  .AND.   constant_heatflux  .AND. &
     1948    IF ( ibc_pt_b == 0  .AND.  constant_heatflux  .AND.                        &
    19401949         surface_heatflux /= 0.0_wp )  THEN
    19411950       message_string = 'boundary_condition: bc_pt_b = "' // TRIM( bc_pt_b ) //&
     
    19441953    ENDIF
    19451954    IF ( constant_heatflux  .AND.  pt_surface_initial_change /= 0.0_wp )  THEN
    1946        WRITE ( message_string, * )  'constant_heatflux = .TRUE. is not allo', &
    1947                'wed with pt_surface_initial_change (/=0) = ', &
     1955       WRITE ( message_string, * )  'constant_heatflux = .TRUE. is not allo',  &
     1956               'wed with pt_surface_initial_change (/=0) = ',                  &
    19481957               pt_surface_initial_change
    19491958       CALL message( 'check_parameters', 'PA0066', 1, 2, 0, 6, 0 )
     
    19541963!-- temperature. In this case specification of a constant heat flux is
    19551964!-- forbidden.
    1956     IF ( ibc_pt_t == 0  .AND.   constant_top_heatflux  .AND. &
     1965    IF ( ibc_pt_t == 0  .AND.  constant_top_heatflux  .AND.                    &
    19571966         top_heatflux /= 0.0_wp )  THEN
    19581967       message_string = 'boundary_condition: bc_pt_t = "' // TRIM( bc_pt_t ) //&
     
    19691978          ibc_sa_t = 1
    19701979       ELSE
    1971           message_string = 'unknown boundary condition: bc_sa_t = "' // &
     1980          message_string = 'unknown boundary condition: bc_sa_t = "' //        &
    19721981                           TRIM( bc_sa_t ) // '"'
    19731982          CALL message( 'check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
     
    19751984
    19761985       IF ( top_salinityflux == 9999999.9_wp )  constant_top_salinityflux = .FALSE.
    1977        IF ( ibc_sa_t == 1  .AND.   top_salinityflux == 9999999.9_wp )  THEN
    1978           message_string = 'boundary condition: bc_sa_t = "' // &
    1979                            TRIM( bc_sa_t ) // '" requires to set ' // &
     1986       IF ( ibc_sa_t == 1  .AND.  top_salinityflux == 9999999.9_wp )  THEN
     1987          message_string = 'boundary condition: bc_sa_t = "' //                &
     1988                           TRIM( bc_sa_t ) // '" requires to set ' //          &
    19801989                           'top_salinityflux'
    19811990          CALL message( 'check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
     
    19861995!--    salinity. In this case specification of a constant salinity flux is
    19871996!--    forbidden.
    1988        IF ( ibc_sa_t == 0  .AND.   constant_top_salinityflux  .AND. &
     1997       IF ( ibc_sa_t == 0  .AND.  constant_top_salinityflux  .AND.            &
    19891998            top_salinityflux /= 0.0_wp )  THEN
    1990           message_string = 'boundary condition: bc_sa_t = "' // &
    1991                            TRIM( bc_sa_t ) // '" is not allowed with ' // &
     1999          message_string = 'boundary condition: bc_sa_t = "' //                &
     2000                           TRIM( bc_sa_t ) // '" is not allowed with ' //      &
    19922001                           'constant_top_salinityflux = .TRUE.'
    19932002          CALL message( 'check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
     
    20102019          ibc_q_b = 1
    20112020       ELSE
    2012           message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // &
     2021          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
    20132022                           '_b ="' // TRIM( bc_q_b ) // '"'
    20142023          CALL message( 'check_parameters', 'PA0071', 1, 2, 0, 6, 0 )
     
    20212030          ibc_q_t = 3
    20222031       ELSE
    2023           message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // &
     2032          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
    20242033                           '_t ="' // TRIM( bc_q_t ) // '"'
    20252034          CALL message( 'check_parameters', 'PA0072', 1, 2, 0, 6, 0 )
     
    20402049       ELSE
    20412050          constant_waterflux = .TRUE.
    2042           IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
    2043                  large_scale_forcing ) THEN
     2051          IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.      &
     2052               large_scale_forcing ) THEN
    20442053             surface_waterflux = qsws_surf(1)
    20452054          ENDIF
     
    20582067       IF ( constant_waterflux  .AND.  q_surface_initial_change /= 0.0_wp )  THEN
    20592068          WRITE( message_string, * )  'a prescribed surface flux is not allo', &
    2060                  'wed with ', sq, '_surface_initial_change (/=0) = ', &
     2069                 'wed with ', sq, '_surface_initial_change (/=0) = ',          &
    20612070                 q_surface_initial_change
    20622071          CALL message( 'check_parameters', 'PA0074', 1, 2, 0, 6, 0 )
     
    20712080       ibc_uv_b = 1
    20722081       IF ( constant_flux_layer )  THEN
    2073           message_string = 'boundary condition: bc_uv_b = "' // &
     2082          message_string = 'boundary condition: bc_uv_b = "' //                &
    20742083               TRIM( bc_uv_b ) // '" is not allowed with constant_flux_layer'  &
    20752084               // ' = .TRUE.'
     
    20772086       ENDIF
    20782087    ELSE
    2079        message_string = 'unknown boundary condition: bc_uv_b = "' // &
     2088       message_string = 'unknown boundary condition: bc_uv_b = "' //           &
    20802089                        TRIM( bc_uv_b ) // '"'
    20812090       CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )
     
    21062115          ibc_uv_t = 3
    21072116       ELSE
    2108           message_string = 'unknown boundary condition: bc_uv_t = "' // &
     2117          message_string = 'unknown boundary condition: bc_uv_t = "' //        &
    21092118                           TRIM( bc_uv_t ) // '"'
    21102119          CALL message( 'check_parameters', 'PA0077', 1, 2, 0, 6, 0 )
     
    21172126       rayleigh_damping_factor = 0.0_wp
    21182127    ELSE
    2119        IF ( rayleigh_damping_factor < 0.0_wp .OR. rayleigh_damping_factor > 1.0_wp ) &
    2120        THEN
    2121           WRITE( message_string, * )  'rayleigh_damping_factor = ', &
     2128       IF ( rayleigh_damping_factor < 0.0_wp  .OR.                            &
     2129            rayleigh_damping_factor > 1.0_wp )  THEN
     2130          WRITE( message_string, * )  'rayleigh_damping_factor = ',            &
    21222131                              rayleigh_damping_factor, ' out of range [0.0,1.0]'
    21232132          CALL message( 'check_parameters', 'PA0078', 1, 2, 0, 6, 0 )
     
    21262135
    21272136    IF ( rayleigh_damping_height == -1.0_wp )  THEN
    2128        IF ( .NOT. ocean )  THEN
     2137       IF (  .NOT. ocean )  THEN
    21292138          rayleigh_damping_height = 0.66666666666_wp * zu(nzt)
    21302139       ELSE
     
    21322141       ENDIF
    21332142    ELSE
    2134        IF ( .NOT. ocean )  THEN
    2135           IF ( rayleigh_damping_height < 0.0_wp  .OR. &
     2143       IF (  .NOT. ocean )  THEN
     2144          IF ( rayleigh_damping_height < 0.0_wp  .OR.                          &
    21362145               rayleigh_damping_height > zu(nzt) )  THEN
    2137              WRITE( message_string, * )  'rayleigh_damping_height = ', &
     2146             WRITE( message_string, * )  'rayleigh_damping_height = ',         &
    21382147                   rayleigh_damping_height, ' out of range [0.0,', zu(nzt), ']'
    21392148             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
    21402149          ENDIF
    21412150       ELSE
    2142           IF ( rayleigh_damping_height > 0.0_wp  .OR. &
     2151          IF ( rayleigh_damping_height > 0.0_wp  .OR.                          &
    21432152               rayleigh_damping_height < zu(nzb) )  THEN
    2144              WRITE( message_string, * )  'rayleigh_damping_height = ', &
     2153             WRITE( message_string, * )  'rayleigh_damping_height = ',         &
    21452154                   rayleigh_damping_height, ' out of range [0.0,', zu(nzb), ']'
    21462155             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
     
    21542163!-- be opened (cf. check_open)
    21552164    IF ( statistic_regions > 9  .OR.  statistic_regions < 0 )  THEN
    2156        WRITE ( message_string, * ) 'number of statistic_regions = ', &
     2165       WRITE ( message_string, * ) 'number of statistic_regions = ',           &
    21572166                   statistic_regions+1, ' but only 10 regions are allowed'
    21582167       CALL message( 'check_parameters', 'PA0082', 1, 2, 0, 6, 0 )
    21592168    ENDIF
    2160     IF ( normalizing_region > statistic_regions  .OR. &
     2169    IF ( normalizing_region > statistic_regions  .OR.                          &
    21612170         normalizing_region < 0)  THEN
    2162        WRITE ( message_string, * ) 'normalizing_region = ', &
     2171       WRITE ( message_string, * ) 'normalizing_region = ',                    &
    21632172                normalizing_region, ' must be >= 0 and <= ',statistic_regions, &
    21642173                ' (value of statistic_regions)'
     
    21842193!
    21852194!-- Set the default skip time intervals for data output, if necessary
    2186     IF ( skip_time_dopr    == 9999999.9_wp ) &
     2195    IF ( skip_time_dopr    == 9999999.9_wp )                                   &
    21872196                                       skip_time_dopr    = skip_time_data_output
    2188     IF ( skip_time_dosp    == 9999999.9_wp ) &
     2197    IF ( skip_time_dosp    == 9999999.9_wp )                                   &
    21892198                                       skip_time_dosp    = skip_time_data_output
    2190     IF ( skip_time_do2d_xy == 9999999.9_wp ) &
     2199    IF ( skip_time_do2d_xy == 9999999.9_wp )                                   &
    21912200                                       skip_time_do2d_xy = skip_time_data_output
    2192     IF ( skip_time_do2d_xz == 9999999.9_wp ) &
     2201    IF ( skip_time_do2d_xz == 9999999.9_wp )                                   &
    21932202                                       skip_time_do2d_xz = skip_time_data_output
    2194     IF ( skip_time_do2d_yz == 9999999.9_wp ) &
     2203    IF ( skip_time_do2d_yz == 9999999.9_wp )                                   &
    21952204                                       skip_time_do2d_yz = skip_time_data_output
    2196     IF ( skip_time_do3d    == 9999999.9_wp ) &
     2205    IF ( skip_time_do3d    == 9999999.9_wp )                                   &
    21972206                                       skip_time_do3d    = skip_time_data_output
    2198     IF ( skip_time_data_output_av == 9999999.9_wp ) &
     2207    IF ( skip_time_data_output_av == 9999999.9_wp )                            &
    21992208                                skip_time_data_output_av = skip_time_data_output
    22002209    DO  mid = 1, max_masks
    2201        IF ( skip_time_domask(mid) == 9999999.9_wp ) &
     2210       IF ( skip_time_domask(mid) == 9999999.9_wp )                            &
    22022211                                skip_time_domask(mid)    = skip_time_data_output
    22032212    ENDDO
     
    22072216!-- spectra)
    22082217    IF ( averaging_interval > dt_data_output_av )  THEN
    2209        WRITE( message_string, * )  'averaging_interval = ', &
     2218       WRITE( message_string, * )  'averaging_interval = ',                    &
    22102219             averaging_interval, ' must be <= dt_data_output = ', dt_data_output
    22112220       CALL message( 'check_parameters', 'PA0085', 1, 2, 0, 6, 0 )
     
    22172226
    22182227    IF ( averaging_interval_pr > dt_dopr )  THEN
    2219        WRITE( message_string, * )  'averaging_interval_pr = ', &
     2228       WRITE( message_string, * )  'averaging_interval_pr = ',                 &
    22202229             averaging_interval_pr, ' must be <= dt_dopr = ', dt_dopr
    22212230       CALL message( 'check_parameters', 'PA0086', 1, 2, 0, 6, 0 )
     
    22272236
    22282237    IF ( averaging_interval_sp > dt_dosp )  THEN
    2229        WRITE( message_string, * )  'averaging_interval_sp = ', &
     2238       WRITE( message_string, * )  'averaging_interval_sp = ',                 &
    22302239             averaging_interval_sp, ' must be <= dt_dosp = ', dt_dosp
    22312240       CALL message( 'check_parameters', 'PA0087', 1, 2, 0, 6, 0 )
     
    22522261!-- Check the sample rate for averaging (first for 3d-data, then for profiles)
    22532262    IF ( dt_averaging_input > averaging_interval )  THEN
    2254        WRITE( message_string, * )  'dt_averaging_input = ', &
    2255                 dt_averaging_input, ' must be <= averaging_interval = ', &
     2263       WRITE( message_string, * )  'dt_averaging_input = ',                    &
     2264                dt_averaging_input, ' must be <= averaging_interval = ',       &
    22562265                averaging_interval
    22572266       CALL message( 'check_parameters', 'PA0088', 1, 2, 0, 6, 0 )
     
    22592268
    22602269    IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
    2261        WRITE( message_string, * )  'dt_averaging_input_pr = ', &
     2270       WRITE( message_string, * )  'dt_averaging_input_pr = ',                 &
    22622271                dt_averaging_input_pr, ' must be <= averaging_interval_pr = ', &
    22632272                averaging_interval_pr
     
    22722281       ELSE
    22732282          IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
    2274              WRITE( message_string, * )  'precipitation_amount_interval = ', &
    2275                  precipitation_amount_interval, ' must not be larger than ', &
     2283             WRITE( message_string, * )  'precipitation_amount_interval = ',   &
     2284                 precipitation_amount_interval, ' must not be larger than ',   &
    22762285                 'dt_do2d_xy = ', dt_do2d_xy
    22772286             CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
     
    25322541          CASE ( 's', '#s' )
    25332542             IF ( .NOT. passive_scalar )  THEN
    2534                 message_string = 'data_output_pr = ' // &
     2543                message_string = 'data_output_pr = ' //                        &
    25352544                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
    25362545                                 'lemented for passive_scalar = .FALSE.'
     
    25732582          CASE ( 'lpt', '#lpt' )
    25742583             IF ( .NOT. cloud_physics ) THEN
    2575                 message_string = 'data_output_pr = ' // &
     2584                message_string = 'data_output_pr = ' //                        &
    25762585                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
    25772586                                 'lemented for cloud_physics = .FALSE.'
     
    26162625
    26172626          CASE ( 'w"q"' )
    2618              IF ( .NOT. humidity )  THEN
    2619                 message_string = 'data_output_pr = ' // &
     2627             IF (  .NOT. humidity )  THEN
     2628                message_string = 'data_output_pr = ' //                        &
    26202629                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
    26212630                                 'lemented for humidity = .FALSE.'
     
    26282637
    26292638          CASE ( 'w*q*' )
    2630              IF ( .NOT. humidity )  THEN
    2631                 message_string = 'data_output_pr = ' // &
     2639             IF (  .NOT. humidity )  THEN
     2640                message_string = 'data_output_pr = ' //                        &
    26322641                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
    26332642                                 'lemented for humidity = .FALSE.'
     
    26402649
    26412650          CASE ( 'wq' )
    2642              IF ( .NOT. humidity )  THEN
    2643                 message_string = 'data_output_pr = ' // &
     2651             IF (  .NOT. humidity )  THEN
     2652                message_string = 'data_output_pr = ' //                        &
    26442653                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
    26452654                                 'lemented for humidity = .FALSE.'
     
    26522661
    26532662          CASE ( 'w"s"' )
    2654              IF ( .NOT. passive_scalar ) THEN
    2655                 message_string = 'data_output_pr = ' // &
     2663             IF (  .NOT.  passive_scalar ) THEN
     2664                message_string = 'data_output_pr = ' //                        &
    26562665                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
    26572666                                 'lemented for passive_scalar = .FALSE.'
     
    26642673
    26652674          CASE ( 'w*s*' )
    2666              IF ( .NOT. passive_scalar ) THEN
    2667                 message_string = 'data_output_pr = ' // &
     2675             IF (  .NOT.  passive_scalar ) THEN
     2676                message_string = 'data_output_pr = ' //                        &
    26682677                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
    26692678                                 'lemented for passive_scalar = .FALSE.'
     
    26762685
    26772686          CASE ( 'ws' )
    2678              IF ( .NOT. passive_scalar ) THEN
    2679                 message_string = 'data_output_pr = ' // &
     2687             IF (  .NOT.  passive_scalar ) THEN
     2688                message_string = 'data_output_pr = ' //                        &
    26802689                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
    26812690                                 'lemented for passive_scalar = .FALSE.'
     
    26882697
    26892698          CASE ( 'w"qv"' )
    2690              IF ( humidity  .AND.  .NOT. cloud_physics ) &
    2691              THEN
     2699             IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
    26922700                dopr_index(i) = 48
    26932701                dopr_unit(i)  = 'kg/kg m/s'
    26942702                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
    2695              ELSEIF( humidity .AND. cloud_physics ) THEN
     2703             ELSEIF ( humidity  .AND.  cloud_physics ) THEN
    26962704                dopr_index(i) = 51
    26972705                dopr_unit(i)  = 'kg/kg m/s'
    26982706                hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
    26992707             ELSE
    2700                 message_string = 'data_output_pr = ' // &
     2708                message_string = 'data_output_pr = ' //                        &
    27012709                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
    27022710                                 'lemented for cloud_physics = .FALSE. an&' // &
     
    27062714
    27072715          CASE ( 'w*qv*' )
    2708              IF ( humidity  .AND.  .NOT. cloud_physics ) &
     2716             IF ( humidity  .AND.  .NOT. cloud_physics )                       &
    27092717             THEN
    27102718                dopr_index(i) = 49
     
    27162724                hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
    27172725             ELSE
    2718                 message_string = 'data_output_pr = ' // &
     2726                message_string = 'data_output_pr = ' //                        &
    27192727                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
    27202728                                 'lemented for cloud_physics = .FALSE. an&' // &
     
    27242732
    27252733          CASE ( 'wqv' )
    2726              IF ( humidity  .AND.  .NOT. cloud_physics ) &
    2727              THEN
     2734             IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
    27282735                dopr_index(i) = 50
    27292736                dopr_unit(i)  = 'kg/kg m/s'
    27302737                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
    2731              ELSEIF( humidity .AND. cloud_physics ) THEN
     2738             ELSEIF ( humidity  .AND.  cloud_physics ) THEN
    27322739                dopr_index(i) = 53
    27332740                dopr_unit(i)  = 'kg/kg m/s'
     
    27422749
    27432750          CASE ( 'ql' )
    2744              IF ( .NOT. cloud_physics  .AND.  .NOT. cloud_droplets )  THEN
     2751             IF (  .NOT.  cloud_physics  .AND.  .NOT. cloud_droplets )  THEN
    27452752                message_string = 'data_output_pr = ' //                        &
    27462753                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    28012808
    28022809          CASE ( 'rho' )
    2803              IF ( .NOT. ocean ) THEN
     2810             IF (  .NOT. ocean ) THEN
    28042811                message_string = 'data_output_pr = ' //                        &
    28052812                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    28192826
    28202827          CASE ( 'w"sa"' )
    2821              IF ( .NOT. ocean ) THEN
     2828             IF (  .NOT. ocean ) THEN
    28222829                message_string = 'data_output_pr = ' //                        &
    28232830                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    28312838
    28322839          CASE ( 'w*sa*' )
    2833              IF ( .NOT. ocean ) THEN
     2840             IF (  .NOT. ocean ) THEN
    28342841                message_string = 'data_output_pr = ' //                        &
    28352842                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    28432850
    28442851          CASE ( 'wsa' )
    2845              IF ( .NOT. ocean ) THEN
     2852             IF (  .NOT. ocean ) THEN
    28462853                message_string = 'data_output_pr = ' //                        &
    28472854                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    28652872
    28662873          CASE ( 'q*2' )
    2867              IF ( .NOT. humidity )  THEN
     2874             IF (  .NOT. humidity )  THEN
    28682875                message_string = 'data_output_pr = ' //                        &
    28692876                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    28772884
    28782885          CASE ( 'prho' )
    2879              IF ( .NOT. ocean ) THEN
     2886             IF (  .NOT. ocean ) THEN
    28802887                message_string = 'data_output_pr = ' //                        &
    28812888                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    28942901
    28952902          CASE ( 'nr' )
    2896              IF ( .NOT. cloud_physics )  THEN
     2903             IF (  .NOT. cloud_physics )  THEN
    28972904                message_string = 'data_output_pr = ' //                        &
    28982905                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    29042911                                 'lemented for cloud_scheme /= seifert_beheng'
    29052912                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
    2906              ELSEIF ( .NOT. precipitation )  THEN
     2913             ELSEIF (  .NOT. precipitation )  THEN
    29072914                message_string = 'data_output_pr = ' //                        &
    29082915                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    29162923
    29172924          CASE ( 'qr' )
    2918              IF ( .NOT. cloud_physics )  THEN
     2925             IF (  .NOT. cloud_physics )  THEN
    29192926                message_string = 'data_output_pr = ' //                        &
    29202927                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    29262933                                 'lemented for cloud_scheme /= seifert_beheng'
    29272934                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
    2928              ELSEIF ( .NOT. precipitation )  THEN
     2935             ELSEIF (  .NOT. precipitation )  THEN
    29292936                message_string = 'data_output_pr = ' //                        &
    29302937                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    29382945
    29392946          CASE ( 'qc' )
    2940              IF ( .NOT. cloud_physics )  THEN
     2947             IF (  .NOT. cloud_physics )  THEN
    29412948                message_string = 'data_output_pr = ' //                        &
    29422949                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    29552962
    29562963          CASE ( 'prr' )
    2957              IF ( .NOT. cloud_physics )  THEN
     2964             IF (  .NOT. cloud_physics )  THEN
    29582965                message_string = 'data_output_pr = ' //                        &
    29592966                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    29652972                                 'lemented for cloud_scheme /= seifert_beheng'
    29662973                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
    2967              ELSEIF ( .NOT. precipitation )  THEN
     2974             ELSEIF (  .NOT. precipitation )  THEN
    29682975                message_string = 'data_output_pr = ' //                        &
    29692976                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    29882995
    29892996          CASE ( 'w_subs' )
    2990              IF ( .NOT. large_scale_subsidence )  THEN
     2997             IF (  .NOT. large_scale_subsidence )  THEN
    29912998                message_string = 'data_output_pr = ' //                        &
    29922999                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    30003007
    30013008          CASE ( 'td_lsa_lpt' )
    3002              IF ( .NOT. large_scale_forcing )  THEN
     3009             IF (  .NOT. large_scale_forcing )  THEN
    30033010                message_string = 'data_output_pr = ' //                        &
    30043011                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    30123019
    30133020          CASE ( 'td_lsa_q' )
    3014              IF ( .NOT. large_scale_forcing )  THEN
     3021             IF (  .NOT. large_scale_forcing )  THEN
    30153022                message_string = 'data_output_pr = ' //                        &
    30163023                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    30243031
    30253032          CASE ( 'td_sub_lpt' )
    3026              IF ( .NOT. large_scale_forcing )  THEN
     3033             IF (  .NOT. large_scale_forcing )  THEN
    30273034                message_string = 'data_output_pr = ' //                        &
    30283035                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    30363043
    30373044          CASE ( 'td_sub_q' )
    3038              IF ( .NOT. large_scale_forcing )  THEN
     3045             IF (  .NOT. large_scale_forcing )  THEN
    30393046                message_string = 'data_output_pr = ' //                        &
    30403047                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    30483055
    30493056          CASE ( 'td_nud_lpt' )
    3050              IF ( .NOT. nudging )  THEN
     3057             IF (  .NOT. nudging )  THEN
    30513058                message_string = 'data_output_pr = ' //                        &
    30523059                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    30603067
    30613068          CASE ( 'td_nud_q' )
    3062              IF ( .NOT. nudging )  THEN
     3069             IF (  .NOT. nudging )  THEN
    30633070                message_string = 'data_output_pr = ' //                        &
    30643071                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    30723079
    30733080          CASE ( 'td_nud_u' )
    3074              IF ( .NOT. nudging )  THEN
     3081             IF (  .NOT. nudging )  THEN
    30753082                message_string = 'data_output_pr = ' //                        &
    30763083                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    30843091
    30853092          CASE ( 'td_nud_v' )
    3086              IF ( .NOT. nudging )  THEN
     3093             IF (  .NOT. nudging )  THEN
    30873094                message_string = 'data_output_pr = ' //                        &
    30883095                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    30963103
    30973104          CASE ( 't_soil', '#t_soil' )
    3098              IF ( .NOT. land_surface )  THEN
     3105             IF (  .NOT. land_surface )  THEN
    30993106                message_string = 'data_output_pr = ' //                        &
    31003107                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    31133120
    31143121          CASE ( 'm_soil', '#m_soil' )
    3115              IF ( .NOT. land_surface )  THEN
     3122             IF (  .NOT. land_surface )  THEN
    31163123                message_string = 'data_output_pr = ' //                        &
    31173124                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     
    31303137
    31313138          CASE ( 'rad_net' )
    3132              IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
     3139             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
     3140             THEN
    31333141                message_string = 'data_output_pr = ' //                        &
    31343142                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
     
    31433151
    31443152          CASE ( 'rad_lw_in' )
    3145              IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
     3153             IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
     3154             THEN
    31463155                message_string = 'data_output_pr = ' //                        &
    31473156                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
     
    31563165
    31573166          CASE ( 'rad_lw_out' )
    3158              IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
     3167             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
     3168             THEN
    31593169                message_string = 'data_output_pr = ' //                        &
    31603170                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
     
    31693179
    31703180          CASE ( 'rad_sw_in' )
    3171              IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
     3181             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
     3182             THEN
    31723183                message_string = 'data_output_pr = ' //                        &
    31733184                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
     
    31823193
    31833194          CASE ( 'rad_sw_out')
    3184              IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
     3195             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' ) &
     3196             THEN
    31853197                message_string = 'data_output_pr = ' //                        &
    31863198                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
     
    31953207
    31963208          CASE ( 'rad_lw_cs_hr' )
    3197              IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' )  THEN
     3209             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
     3210             THEN
    31983211                message_string = 'data_output_pr = ' //                        &
    31993212                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
     
    32083221
    32093222          CASE ( 'rad_lw_hr' )
    3210              IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' )  THEN
     3223             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
     3224             THEN
    32113225                message_string = 'data_output_pr = ' //                        &
    32123226                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
     
    32213235
    32223236          CASE ( 'rad_sw_cs_hr' )
    3223              IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' )  THEN
     3237             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
     3238             THEN
    32243239                message_string = 'data_output_pr = ' //                        &
    32253240                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
     
    32343249
    32353250          CASE ( 'rad_sw_hr' )
    3236              IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' )  THEN
     3251             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
     3252             THEN
    32373253                message_string = 'data_output_pr = ' //                        &
    32383254                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
     
    33293345
    33303346          CASE ( 'lpt' )
    3331              IF ( .NOT. cloud_physics )  THEN
     3347             IF (  .NOT. cloud_physics )  THEN
    33323348                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    33333349                         'res cloud_physics = .TRUE.'
     
    33373353
    33383354          CASE ( 'm_soil' )
    3339              IF ( .NOT. land_surface )  THEN
     3355             IF (  .NOT. land_surface )  THEN
    33403356                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    33413357                         'land_surface = .TRUE.'
     
    33453361
    33463362          CASE ( 'nr' )
    3347              IF ( .NOT. cloud_physics )  THEN
     3363             IF (  .NOT. cloud_physics )  THEN
    33483364                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    33493365                         'res cloud_physics = .TRUE.'
     
    33573373
    33583374          CASE ( 'pc', 'pr' )
    3359              IF ( .NOT. particle_advection )  THEN
     3375             IF (  .NOT. particle_advection )  THEN
    33603376                message_string = 'output of "' // TRIM( var ) // '" requir' // &
    33613377                   'es a "particles_par"-NAMELIST in the parameter file (PARIN)'
     
    33663382
    33673383          CASE ( 'prr' )
    3368              IF ( .NOT. cloud_physics )  THEN
     3384             IF (  .NOT. cloud_physics )  THEN
    33693385                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    33703386                         'res cloud_physics = .TRUE.'
     
    33743390                         'res cloud_scheme = seifert_beheng'
    33753391                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
    3376              ELSEIF ( .NOT. precipitation )  THEN
     3392             ELSEIF (  .NOT. precipitation )  THEN
    33773393                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    33783394                                 'res precipitation = .TRUE.'
     
    33823398
    33833399          CASE ( 'q', 'vpt' )
    3384              IF ( .NOT. humidity )  THEN
     3400             IF (  .NOT. humidity )  THEN
    33853401                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    33863402                                 'res humidity = .TRUE.'
     
    33913407
    33923408          CASE ( 'qc' )
    3393              IF ( .NOT. cloud_physics )  THEN
     3409             IF (  .NOT. cloud_physics )  THEN
    33943410                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    33953411                         'res cloud_physics = .TRUE.'
     
    34033419
    34043420          CASE ( 'ql' )
    3405              IF ( .NOT. ( cloud_physics  .OR.  cloud_droplets ) )  THEN
     3421             IF ( .NOT.  ( cloud_physics  .OR.  cloud_droplets ) )  THEN
    34063422                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    34073423                         'res cloud_physics = .TRUE. or cloud_droplets = .TRUE.'
     
    34113427
    34123428          CASE ( 'ql_c', 'ql_v', 'ql_vp' )
    3413              IF ( .NOT. cloud_droplets )  THEN
     3429             IF (  .NOT. cloud_droplets )  THEN
    34143430                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    34153431                                 'res cloud_droplets = .TRUE.'
     
    34213437
    34223438          CASE ( 'qr' )
    3423              IF ( .NOT. cloud_physics )  THEN
     3439             IF (  .NOT. cloud_physics )  THEN
    34243440                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    34253441                         'res cloud_physics = .TRUE.'
     
    34293445                         'res cloud_scheme = seifert_beheng'
    34303446                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
    3431              ELSEIF ( .NOT. precipitation )  THEN
     3447             ELSEIF (  .NOT. precipitation )  THEN
    34323448                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    34333449                                 'res precipitation = .TRUE.'
     
    34373453
    34383454          CASE ( 'qv' )
    3439              IF ( .NOT. cloud_physics )  THEN
     3455             IF (  .NOT. cloud_physics )  THEN
    34403456                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    34413457                                 'res cloud_physics = .TRUE.'
     
    34463462          CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_lw_cs_hr', 'rad_lw_hr',       &
    34473463                 'rad_sw_in', 'rad_sw_out', 'rad_sw_cs_hr', 'rad_sw_hr' )
    3448              IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' )  THEN
     3464             IF (  .NOT.  radiation  .OR. radiation_scheme /= 'rrtmg' )  THEN
    34493465                message_string = '"output of "' // TRIM( var ) // '" requi' // &
    34503466                                 'res radiation = .TRUE. and ' //              &
     
    34553471
    34563472          CASE ( 'rho' )
    3457              IF ( .NOT. ocean )  THEN
     3473             IF (  .NOT. ocean )  THEN
    34583474                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    34593475                                 'res ocean = .TRUE.'
     
    34633479
    34643480          CASE ( 's' )
    3465              IF ( .NOT. passive_scalar )  THEN
     3481             IF (  .NOT. passive_scalar )  THEN
    34663482                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    34673483                                 'res passive_scalar = .TRUE.'
     
    34713487
    34723488          CASE ( 'sa' )
    3473              IF ( .NOT. ocean )  THEN
     3489             IF (  .NOT. ocean )  THEN
    34743490                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    34753491                                 'res ocean = .TRUE.'
     
    34793495
    34803496          CASE ( 't_soil' )
    3481              IF ( .NOT. land_surface )  THEN
     3497             IF (  .NOT. land_surface )  THEN
    34823498                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    34833499                         'land_surface = .TRUE.'
     
    34913507                 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*', 'rad_net*',  &
    34923508                 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*',   &
    3493                  'r_a*', 'r_s*', 'shf*', 'shf_eb*', 't*', 'u*', 'z0*', 'z0h*' )
     3509                 'r_a*', 'r_s*', 'shf*', 'shf_eb*', 't*', 'u*', 'z0*', 'z0h*', &
     3510                 'z0q*' )
    34943511             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
    34953512                message_string = 'illegal value for data_output: "' //         &
     
    34983515                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
    34993516             ENDIF
    3500              IF ( .NOT. radiation .OR. radiation_scheme /= "rrtmg" )  THEN
    3501                 IF ( TRIM( var ) == 'rrtm_aldif*' .OR.                         &
    3502                      TRIM( var ) == 'rrtm_aldir*' .OR.                         &
    3503                      TRIM( var ) == 'rrtm_asdif*' .OR.                         &
     3517             IF (  .NOT.  radiation  .OR. radiation_scheme /= "rrtmg" )  THEN
     3518                IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
     3519                     TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
     3520                     TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
    35043521                     TRIM( var ) == 'rrtm_asdir*'      )                       &
    35053522                THEN
     
    35113528             ENDIF
    35123529
    3513              IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT. land_surface )  THEN
     3530             IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT.  land_surface )  THEN
    35143531                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    35153532                                 'res land_surface = .TRUE.'
    35163533                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    35173534             ENDIF
    3518              IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT. land_surface )  THEN
     3535             IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT.  land_surface )  THEN
    35193536                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    35203537                                 'res land_surface = .TRUE.'
     
    35263543                CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    35273544             ENDIF
    3528              IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT. land_surface )  THEN
     3545             IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT.  land_surface )  THEN
    35293546                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    35303547                                 'res land_surface = .TRUE.'
    35313548                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    35323549             ENDIF
    3533              IF ( TRIM( var ) == 'lai*'  .AND.  .NOT. land_surface )  THEN
     3550             IF ( TRIM( var ) == 'lai*'  .AND.  .NOT.  land_surface )  THEN
    35343551                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    35353552                                 'res land_surface = .TRUE.'
     
    35413558                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
    35423559             ENDIF
    3543              IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT. land_surface )  THEN
     3560             IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT.  land_surface )  THEN
    35443561                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    35453562                                 'res land_surface = .TRUE.'
     
    35563573                CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
    35573574             ENDIF
    3558              IF ( TRIM( var ) == 'prr*'  .AND.  .NOT. precipitation )  THEN
     3575             IF ( TRIM( var ) == 'prr*'  .AND.  .NOT.  precipitation )  THEN
    35593576                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    35603577                                 'res precipitation = .TRUE.'
    35613578                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
    35623579             ENDIF
    3563              IF ( TRIM( var ) == 'qsws*'  .AND.  .NOT. humidity )  THEN
     3580             IF ( TRIM( var ) == 'qsws*'  .AND.  .NOT.  humidity )  THEN
    35643581                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    35653582                                 'res humidity = .TRUE.'
    35663583                CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
    35673584             ENDIF
    3568              IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT. land_surface )  THEN
     3585             IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT.  land_surface )  THEN
    35693586                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    35703587                                 'res land_surface = .TRUE.'
    35713588                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    35723589             ENDIF
    3573              IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )  &
     3590             IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )   &
    35743591             THEN
    35753592                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     
    35773594                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    35783595             ENDIF
    3579              IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT. land_surface ) &
     3596             IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT.  land_surface ) &
    35803597             THEN
    35813598                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     
    35833600                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    35843601             ENDIF
    3585              IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )  &
     3602             IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )   &
    35863603             THEN
    35873604                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     
    35893606                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    35903607             ENDIF
    3591              IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT. land_surface ) &
     3608             IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT.  land_surface ) &
    35923609             THEN
    35933610                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     
    35953612                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    35963613             ENDIF
    3597              IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT. land_surface ) &
     3614             IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT.  land_surface ) &
    35983615             THEN
    35993616                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     
    37573774!-- Check mask conditions
    37583775    DO mid = 1, max_masks
    3759        IF ( data_output_masks(mid,1) /= ' ' .OR.                               &
     3776       IF ( data_output_masks(mid,1) /= ' '  .OR.                              &
    37603777            data_output_masks_user(mid,1) /= ' ' ) THEN
    37613778          masks = masks + 1
     
    37633780    ENDDO
    37643781   
    3765     IF ( masks < 0 .OR. masks > max_masks )  THEN
     3782    IF ( masks < 0  .OR. masks > max_masks )  THEN
    37663783       WRITE( message_string, * )  'illegal value: masks must be >= 0 and ',   &
    37673784            '<= ', max_masks, ' (=max_masks)'
     
    39083925!
    39093926!-- Check random generator
    3910     IF ( (random_generator /= 'system-specific'     .AND.                      &
    3911           random_generator /= 'random-parallel'   ) .AND.                      &
     3927    IF ( (random_generator /= 'system-specific'      .AND.                     &
     3928          random_generator /= 'random-parallel'   )  .AND.                     &
    39123929          random_generator /= 'numerical-recipes' )  THEN
    39133930       message_string = 'unknown random generator: random_generator = "' //    &
     
    39193936!-- Determine upper and lower hight level indices for random perturbations
    39203937    IF ( disturbance_level_b == -9999999.9_wp )  THEN
    3921        IF ( ocean ) THEN
     3938       IF ( ocean )  THEN
    39223939          disturbance_level_b     = zu((nzt*2)/3)
    39233940          disturbance_level_ind_b = ( nzt * 2 ) / 3
     
    41824199!
    41834200!-- Check pressure gradient conditions
    4184     IF ( dp_external .AND. conserve_volume_flow )  THEN
     4201    IF ( dp_external  .AND. conserve_volume_flow )  THEN
    41854202       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo', &
    41864203            'w are .TRUE. but one of them must be .FALSE.'
     
    41884205    ENDIF
    41894206    IF ( dp_external )  THEN
    4190        IF ( dp_level_b < zu(nzb) .OR. dp_level_b > zu(nzt) )  THEN
     4207       IF ( dp_level_b < zu(nzb)  .OR. dp_level_b > zu(nzt) )  THEN
    41914208          WRITE( message_string, * )  'dp_level_b = ', dp_level_b, ' is out ', &
    41924209               ' of range'
     
    41994216       ENDIF
    42004217    ENDIF
    4201     IF ( ANY( dpdxy /= 0.0_wp ) .AND. .NOT. dp_external )  THEN
     4218    IF ( ANY( dpdxy /= 0.0_wp )  .AND.  .NOT. dp_external )  THEN
    42024219       WRITE( message_string, * )  'dpdxy is nonzero but dp_external is ',     &
    42034220            '.FALSE., i.e. the external pressure gradient & will not be applied'
     
    42304247       ENDIF
    42314248    ENDIF
    4232     IF ( ( u_bulk /= 0.0_wp .OR. v_bulk /= 0.0_wp ) .AND.                      &
    4233          ( .NOT. conserve_volume_flow .OR.                                     &
     4249    IF ( ( u_bulk /= 0.0_wp  .OR.  v_bulk /= 0.0_wp )  .AND.                   &
     4250         ( .NOT. conserve_volume_flow  .OR.                                    &
    42344251         TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) )  THEN
    42354252       WRITE( message_string, * )  'nonzero bulk velocity requires ',          &
     
    42704287!
    42714288!-- Check nudging and large scale forcing from external file
    4272     IF ( nudging .AND. ( .NOT. large_scale_forcing ) )  THEN
     4289    IF ( nudging  .AND.  (  .NOT. large_scale_forcing ) )  THEN
    42734290       message_string = 'Nudging requires large_scale_forcing = .T.. &'//      &
    42744291                        'Surface fluxes and geostrophic wind should be &'//    &
     
    42774294    ENDIF
    42784295
    4279     IF ( large_scale_forcing .AND. ( bc_lr /= 'cyclic'  .OR.                   &
     4296    IF ( large_scale_forcing  .AND.  ( bc_lr /= 'cyclic'  .OR.                 &
    42804297                                    bc_ns /= 'cyclic' ) )  THEN
    42814298       message_string = 'Non-cyclic lateral boundaries do not allow for &' //  &
     
    42844301     ENDIF
    42854302
    4286     IF ( large_scale_forcing .AND. ( .NOT. humidity ) )  THEN
     4303    IF ( large_scale_forcing  .AND.  (  .NOT. humidity ) )  THEN
    42874304       message_string = 'The usage of large scale forcing from external &'//   &
    42884305                        'file LSF_DATA requires humidity = .T..'
     
    42904307     ENDIF
    42914308
    4292     IF ( large_scale_forcing .AND. topography /= 'flat' )  THEN
     4309    IF ( large_scale_forcing  .AND. topography /= 'flat' )  THEN
    42934310       message_string = 'The usage of large scale forcing from external &'//   &
    42944311                        'file LSF_DATA is not implemented for non-flat topography'
     
    42964313    ENDIF
    42974314
    4298     IF ( large_scale_forcing .AND.  ocean  )  THEN
     4315    IF ( large_scale_forcing  .AND.  ocean  )  THEN
    42994316       message_string = 'The usage of large scale forcing from external &'//   &
    43004317                        'file LSF_DATA is not implemented for ocean runs'
     
    43204337!-- Check for valid setting of most_method
    43214338    IF ( TRIM( most_method ) /= 'circular'  .AND.                              &
    4322          TRIM( most_method ) /= 'newton'  .AND.                                &
     4339         TRIM( most_method ) /= 'newton'    .AND.                              &
    43234340         TRIM( most_method ) /= 'lookup' )  THEN
    4324        message_string = 'most_method = "' // TRIM( most_method ) //      &
     4341       message_string = 'most_method = "' // TRIM( most_method ) //            &
    43254342                        '" is unknown'
    43264343       CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )
     
    43524369       IF ( dt_do == 0.0_wp )  THEN
    43534370          IF ( dt_fixed )  THEN
    4354              WRITE( message_string, '(A,F9.4,A)' )  'Output at every '  //  &
    4355                     'timestep is desired (' // dt_do_name // ' = 0.0).&'//  &
    4356                     'Setting the output interval to the fixed timestep '//  &
     4371             WRITE( message_string, '(A,F9.4,A)' )  'Output at every '  //     &
     4372                    'timestep is desired (' // dt_do_name // ' = 0.0).&'//     &
     4373                    'Setting the output interval to the fixed timestep '//     &
    43574374                    'dt = ', dt, 's.'
    43584375             CALL message( 'check_parameters', 'PA0060', 0, 0, 0, 6, 0 )
    43594376             dt_do = dt
    43604377          ELSE
    4361              message_string = dt_do_name // ' = 0.0 while using a ' //      &
    4362                               'variable timestep and parallel netCDF4 ' //  &
     4378             message_string = dt_do_name // ' = 0.0 while using a ' //         &
     4379                              'variable timestep and parallel netCDF4 ' //     &
    43634380                              'is not allowed.'
    43644381             CALL message( 'check_parameters', 'PA0081', 1, 2, 0, 6, 0 )
  • palm/trunk/SOURCE/data_output_2d.f90

    r1784 r1788  
    1919! Current revisions:
    2020! -----------------
     21! Added output of z0q
    2122!
    22 !
    2323! Former revisions:
    2424! -----------------
     
    129129    USE arrays_3d,                                                             &
    130130        ONLY:  dzw, e, nr, ol, p, pt, q, qc, ql, ql_c, ql_v, ql_vp, qr, qsws,  &
    131                rho, sa, shf, tend, ts, u, us, v, vpt, w, z0, z0h, zu, zw
     131               rho, sa, shf, tend, ts, u, us, v, vpt, w, z0, z0h, z0q, zu, zw
    132132       
    133133    USE averaging
     
    12471247                      DO  j = nysg, nyng
    12481248                         local_pf(i,j,nzb+1) =  z0h_av(j,i)
     1249                      ENDDO
     1250                   ENDDO
     1251                ENDIF
     1252                resorted = .TRUE.
     1253                two_d = .TRUE.
     1254                level_z(nzb+1) = zu(nzb+1)
     1255
     1256             CASE ( 'z0q*_xy' )        ! 2d-array
     1257                IF ( av == 0 ) THEN
     1258                   DO  i = nxlg, nxrg
     1259                      DO  j = nysg, nyng
     1260                         local_pf(i,j,nzb+1) =  z0q(j,i)
     1261                      ENDDO
     1262                   ENDDO
     1263                ELSE
     1264                   DO  i = nxlg, nxrg
     1265                      DO  j = nysg, nyng
     1266                         local_pf(i,j,nzb+1) =  z0q_av(j,i)
    12491267                      ENDDO
    12501268                   ENDDO
  • palm/trunk/SOURCE/header.f90

    r1787 r1788  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Parameter dewfall removed
    2222!
    2323! Former revisions:
     
    248248   
    249249    USE land_surface_model_mod,                                                &
    250         ONLY:  conserve_water_content, dewfall, land_surface, nzb_soil,        &
     250        ONLY:  conserve_water_content, land_surface, nzb_soil,                 &
    251251               nzt_soil, root_fraction, soil_moisture, soil_temperature,       &
    252252               soil_type, soil_type_name, veg_type, veg_type_name, zs
     
    964964       ELSE
    965965          WRITE( io, 441 )
    966        ENDIF
    967 
    968        IF ( dewfall )  THEN
    969           WRITE( io, 442 )
    970        ELSE
    971           WRITE( io, 443 )
    972966       ENDIF
    973967
     
    23412335            '       Root fraction: ',A,'  '/ &
    23422336            '       Gridpoint:     ',A)
    2343 440 FORMAT (/' --> Dewfall is allowed (default)')
    2344 441 FORMAT (' --> Dewfall is inhibited')
    2345 442 FORMAT (' --> Soil bottom is closed (water content is conserved, default)')
    2346 443 FORMAT (' --> Soil bottom is open (water content is not conserved)')
     2337440 FORMAT (' --> Soil bottom is closed (water content is conserved, default)')
     2338441 FORMAT (' --> Soil bottom is open (water content is not conserved)')
    23472339444 FORMAT (//' Radiation model information:'/                                 &
    23482340              ' ----------------------------'/)
  • palm/trunk/SOURCE/init_3d_model.f90

    r1784 r1788  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Added z0q.
     22! Syntax layout improved.
    2223!
    2324! Former revisions:
     
    320321!
    321322!-- Allocate arrays
    322     ALLOCATE( mean_surface_level_height(0:statistic_regions),               &
    323               mean_surface_level_height_l(0:statistic_regions),             &
    324               ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), &
    325               ngp_3d(0:statistic_regions),                                  &
    326               ngp_3d_inner(0:statistic_regions),                            &
    327               ngp_3d_inner_l(0:statistic_regions),                          &
    328               ngp_3d_inner_tmp(0:statistic_regions),                        &
    329               sums_divnew_l(0:statistic_regions),                           &
     323    ALLOCATE( mean_surface_level_height(0:statistic_regions),                  &
     324              mean_surface_level_height_l(0:statistic_regions),                &
     325              ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions),    &
     326              ngp_3d(0:statistic_regions),                                     &
     327              ngp_3d_inner(0:statistic_regions),                               &
     328              ngp_3d_inner_l(0:statistic_regions),                             &
     329              ngp_3d_inner_tmp(0:statistic_regions),                           &
     330              sums_divnew_l(0:statistic_regions),                              &
    330331              sums_divold_l(0:statistic_regions) )
    331332    ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt) )
    332     ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                 &
    333               ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),               &
    334               ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),               &
    335               ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),             &
    336               rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),               &
    337               sums(nzb:nzt+1,pr_palm+max_pr_user),                          &
    338               sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1),   &
    339               sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), &
    340               sums_up_fraction_l(10,3,0:statistic_regions),                 &
    341               sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                &
     333    ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                    &
     334              ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),                  &
     335              ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),                  &
     336              ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),                &
     337              rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),                  &
     338              sums(nzb:nzt+1,pr_palm+max_pr_user),                             &
     339              sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1),      &
     340              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1),    &
     341              sums_up_fraction_l(10,3,0:statistic_regions),                    &
     342              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                   &
    342343              ts_value(dots_max,0:statistic_regions) )
    343344    ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) )
    344345
    345     ALLOCATE( ol(nysg:nyng,nxlg:nxrg), shf(nysg:nyng,nxlg:nxrg),      &
    346               ts(nysg:nyng,nxlg:nxrg), tswst(nysg:nyng,nxlg:nxrg),    &
    347               us(nysg:nyng,nxlg:nxrg), usws(nysg:nyng,nxlg:nxrg),     &
    348               uswst(nysg:nyng,nxlg:nxrg), vsws(nysg:nyng,nxlg:nxrg),  &
    349               vswst(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg),    &
    350               z0h(nysg:nyng,nxlg:nxrg) )
    351 
    352     ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr),         &
    353               kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
    354               km(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
    355               p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),     &
     346    ALLOCATE( ol(nysg:nyng,nxlg:nxrg), shf(nysg:nyng,nxlg:nxrg),               &
     347              ts(nysg:nyng,nxlg:nxrg), tswst(nysg:nyng,nxlg:nxrg),             &
     348              us(nysg:nyng,nxlg:nxrg), usws(nysg:nyng,nxlg:nxrg),              &
     349              uswst(nysg:nyng,nxlg:nxrg), vsws(nysg:nyng,nxlg:nxrg),           &
     350              vswst(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg),             &
     351              z0h(nysg:nyng,nxlg:nxrg), z0q(nysg:nyng,nxlg:nxrg) )
     352
     353    ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr),                                    &
     354              kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                               &
     355              km(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                               &
     356              p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
    356357              tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    357358
    358359#if defined( __nopointer )
    359     ALLOCATE( e(nzb:nzt+1,nysg:nyng,nxlg:nxrg),     &
    360               e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    361               pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
    362               pt_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    363               u(nzb:nzt+1,nysg:nyng,nxlg:nxrg),     &
    364               u_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    365               v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),     &
    366               v_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    367               w(nzb:nzt+1,nysg:nyng,nxlg:nxrg),     &
    368               w_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    369               te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    370               tpt_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
    371               tu_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    372               tv_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     360    ALLOCATE( e(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
     361              e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     362              pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                               &
     363              pt_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
     364              u(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
     365              u_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     366              v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
     367              v_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     368              w(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
     369              w_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     370              te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
     371              tpt_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                            &
     372              tu_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
     373              tv_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
    373374              tw_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    374375#else
    375     ALLOCATE( e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    376               e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    377               e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    378               pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    379               pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    380               u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    381               u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    382               u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    383               v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    384               v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    385               v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    386               w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    387               w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
     376    ALLOCATE( e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     377              e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     378              e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     379              pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
     380              pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
     381              u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     382              u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     383              u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     384              v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     385              v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     386              v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     387              w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     388              w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
    388389              w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    389     IF ( .NOT. neutral )  THEN
     390    IF (  .NOT. neutral )  THEN
    390391       ALLOCATE( pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    391392    ENDIF
     
    412413    ENDIF
    413414
    414     IF ( humidity  .OR.  passive_scalar ) THEN
     415    IF ( humidity  .OR.  passive_scalar )  THEN
    415416!
    416417!--    2D-humidity/scalar arrays
    417        ALLOCATE ( qs(nysg:nyng,nxlg:nxrg),   &
    418                   qsws(nysg:nyng,nxlg:nxrg), &
     418       ALLOCATE ( qs(nysg:nyng,nxlg:nxrg),                                     &
     419                  qsws(nysg:nyng,nxlg:nxrg),                                   &
    419420                  qswst(nysg:nyng,nxlg:nxrg) )
    420421
     
    422423!--    3D-humidity/scalar arrays
    423424#if defined( __nopointer )
    424        ALLOCATE( q(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    425                  q_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     425       ALLOCATE( q(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
     426                 q_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
    426427                 tq_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    427428#else
    428        ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
    429                  q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     429       ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
     430                 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
    430431                 q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    431432#endif
     
    440441#endif
    441442
    442           IF ( cloud_physics ) THEN
     443          IF ( cloud_physics )  THEN
    443444
    444445!
     
    451452!
    452453!--          Precipitation amount and rate (only needed if output is switched)
    453              ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg), &
     454             ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg),              &
    454455                       precipitation_rate(nysg:nyng,nxlg:nxrg) )
    455456
     
    457458!
    458459!--             1D-arrays
    459                 ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1), &
     460                ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1),                 &
    460461                           q_1d(nzb:nzt+1), qc_1d(nzb:nzt+1) )
    461462!
     
    474475!
    475476!--                2D-rain water content and rain drop concentration arrays
    476                    ALLOCATE ( qrs(nysg:nyng,nxlg:nxrg),                 &
    477                               qrsws(nysg:nyng,nxlg:nxrg),               &
    478                               qrswst(nysg:nyng,nxlg:nxrg),              &
    479                               nrs(nysg:nyng,nxlg:nxrg),                 &
    480                               nrsws(nysg:nyng,nxlg:nxrg),               &
     477                   ALLOCATE ( qrs(nysg:nyng,nxlg:nxrg),                        &
     478                              qrsws(nysg:nyng,nxlg:nxrg),                      &
     479                              qrswst(nysg:nyng,nxlg:nxrg),                     &
     480                              nrs(nysg:nyng,nxlg:nxrg),                        &
     481                              nrsws(nysg:nyng,nxlg:nxrg),                      &
    481482                              nrswst(nysg:nyng,nxlg:nxrg) )
    482483!
    483484!--                3D-rain water content, rain drop concentration arrays
    484485#if defined( __nopointer )
    485                    ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),         &
    486                              nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
    487                              qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),         &
    488                              qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
    489                              tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),      &
     486                   ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
     487                             nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     488                             qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
     489                             qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     490                             tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),             &
    490491                             tqr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    491492#else
    492                    ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
    493                              nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
    494                              nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
    495                              qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
    496                              qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
     493                   ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     494                             nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     495                             nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     496                             qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     497                             qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
    497498                             qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    498499#endif
     
    509510!--          Liquid water content, change in liquid water content
    510511#if defined( __nopointer )
    511              ALLOCATE ( ql(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     512             ALLOCATE ( ql(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                     &
    512513                        ql_c(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    513514#else
    514              ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     515             ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                   &
    515516                        ql_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    516517#endif
    517518!
    518519!--          Real volume of particles (with weighting), volume of particles
    519              ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     520             ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                   &
    520521                        ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    521522          ENDIF
     
    526527
    527528    IF ( ocean )  THEN
    528        ALLOCATE( saswsb(nysg:nyng,nxlg:nxrg), &
     529       ALLOCATE( saswsb(nysg:nyng,nxlg:nxrg),                                  &
    529530                 saswst(nysg:nyng,nxlg:nxrg) )
    530531#if defined( __nopointer )
    531        ALLOCATE( prho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    532                  rho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
    533                  sa(nzb:nzt+1,nysg:nyng,nxlg:nxrg),     &
    534                  sa_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
     532       ALLOCATE( prho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
     533                 rho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
     534                 sa(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                            &
     535                 sa_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
    535536                 tsa_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    536537#else
    537        ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
    538                  rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    539                  sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
    540                  sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
     538       ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                        &
     539                 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                         &
     540                 sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
     541                 sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
    541542                 sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    542543       prho => prho_1
     
    553554!-- 3D-array for storing the dissipation, needed for calculating the sgs
    554555!-- particle velocities
    555     IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence  .OR.  &
     556    IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence  .OR.      &
    556557         num_acc_per_node > 0 )  THEN
    557558       ALLOCATE( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     
    559560
    560561    IF ( dt_dosp /= 9999999.9_wp )  THEN
    561        ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ), &
     562       ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ),                             &
    562563                 spectrum_y( 1:ny/2, 1:10, 1:10 ) )
    563564       spectrum_x = 0.0_wp
     
    595596!-- are needed for radiation boundary conditions
    596597    IF ( outflow_l )  THEN
    597        ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2), &
    598                  v_m_l(nzb:nzt+1,nysg:nyng,0:1), &
     598       ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2),                               &
     599                 v_m_l(nzb:nzt+1,nysg:nyng,0:1),                               &
    599600                 w_m_l(nzb:nzt+1,nysg:nyng,0:1) )
    600601    ENDIF
    601602    IF ( outflow_r )  THEN
    602        ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &
    603                  v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &
     603       ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                           &
     604                 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                           &
    604605                 w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) )
    605606    ENDIF
    606607    IF ( outflow_l  .OR.  outflow_r )  THEN
    607        ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng), &
     608       ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng),           &
    608609                 c_w(nzb:nzt+1,nysg:nyng) )
    609610    ENDIF
    610611    IF ( outflow_s )  THEN
    611        ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg), &
    612                  v_m_s(nzb:nzt+1,1:2,nxlg:nxrg), &
     612       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg),                               &
     613                 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg),                               &
    613614                 w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) )
    614615    ENDIF
    615616    IF ( outflow_n )  THEN
    616        ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &
    617                  v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &
     617       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                           &
     618                 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                           &
    618619                 w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) )
    619620    ENDIF
    620621    IF ( outflow_s  .OR.  outflow_n )  THEN
    621        ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg), &
     622       ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg),           &
    622623                 c_w(nzb:nzt+1,nxlg:nxrg) )
    623624    ENDIF
     
    682683!--        pres may need to be called  outside the RK-substeps. Antti will
    683684!--        check if this is really required.
    684     ALLOCATE( weight_substep(1:intermediate_timestep_count_max), &
     685    ALLOCATE( weight_substep(1:intermediate_timestep_count_max),               &
    685686              weight_pres(0:intermediate_timestep_count_max) )
    686687    weight_substep = 1.0_wp
     
    691692!
    692693!-- Initialize model variables
    693     IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
     694    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
    694695         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
    695696!
     
    819820          CALL location_message( 'finished', .TRUE. )
    820821
    821        ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
     822       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
    822823       THEN
    823824
     
    825826!
    826827!--       Overwrite initial profiles in case of nudging
    827           IF ( nudging ) THEN
     828          IF ( nudging )  THEN
    828829             pt_init = ptnudge(:,1)
    829830             u_init  = unudge(:,1)
     
    833834             ENDIF
    834835
    835              WRITE( message_string, * ) 'Initial profiles of u, v and ', &
     836             WRITE( message_string, * ) 'Initial profiles of u, v and ',       &
    836837                 'scalars from NUDGING_DATA are used.'
    837838             CALL message( 'init_3d_model', 'PA0370', 0, 0, 0, 6, 0 )
     
    936937          CALL location_message( 'finished', .TRUE. )
    937938
    938        ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) &
     939       ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 )                  &
    939940       THEN
    940941
     
    10181019          DO j=0,ny
    10191020             DO i=0,nx
    1020                 id_random_array(j,i) = random_dummy + 1E6 * ( ensemble_member_nr - 1000 )
     1021                id_random_array(j,i) = random_dummy + 1E6                      &
     1022                                       * ( ensemble_member_nr - 1000 )
    10211023                random_dummy = random_dummy + 1
    10221024             END DO
     
    10471049!
    10481050!--             Initialize shf with data from external file LSF_DATA
    1049                 IF ( large_scale_forcing .AND. lsf_surf ) THEN
     1051                IF ( large_scale_forcing .AND. lsf_surf )  THEN
    10501052                   CALL ls_forcing_surf ( simulated_time )
    10511053                ENDIF
     
    10681070!--       Determine the near-surface water flux
    10691071          IF ( humidity  .OR.  passive_scalar )  THEN
    1070              IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
     1072             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
    10711073                  precipitation )  THEN
    10721074                qrsws = 0.0_wp
     
    11061108             IF ( humidity  .OR.  passive_scalar )  THEN
    11071109                qswst = 0.0_wp
    1108                 IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
     1110                IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.           &
    11091111                     precipitation ) THEN
    11101112                   nrswst = 0.0_wp
     
    11331135          z0 = roughness_length
    11341136          z0h = z0h_factor * z0
     1137          z0q = z0h_factor * z0
    11351138
    11361139          IF ( .NOT. constant_heatflux )  THEN
     
    11451148
    11461149          IF ( humidity  .OR.  passive_scalar )  THEN
    1147              IF ( .NOT. constant_waterflux )  qsws   = 0.0_wp
    1148              IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
     1150             IF (  .NOT. constant_waterflux )  qsws   = 0.0_wp
     1151             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
    11491152                  precipitation )  THEN
    11501153                qrsws = 0.0_wp
     
    11591162!--    the reference state will be set (overwritten) in init_ocean)
    11601163       IF ( use_single_reference_value )  THEN
    1161           IF ( .NOT. humidity )  THEN
     1164          IF (  .NOT. humidity )  THEN
    11621165             ref_state(:) = pt_reference
    11631166          ELSE
     
    11651168          ENDIF
    11661169       ELSE
    1167           IF ( .NOT. humidity )  THEN
     1170          IF (  .NOT. humidity )  THEN
    11681171             ref_state(:) = pt_init(:)
    11691172          ELSE
     
    12161219!--    If required, change the surface humidity/scalar at the start of the 3D
    12171220!--    run
    1218        IF ( ( humidity .OR. passive_scalar ) .AND. &
     1221       IF ( ( humidity .OR. passive_scalar ) .AND.                             &
    12191222            q_surface_initial_change /= 0.0_wp )  THEN
    12201223          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
     
    12291232          tq_m = 0.0_wp
    12301233          q_p = q
    1231           IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
     1234          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
    12321235               precipitation )  THEN
    12331236             tqr_m = 0.0_wp
     
    12451248       CALL location_message( 'finished', .TRUE. )
    12461249
    1247     ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.    &
    1248          TRIM( initializing_actions ) == 'cyclic_fill' )  &
     1250    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.         &
     1251         TRIM( initializing_actions ) == 'cyclic_fill' )                       &
    12491252    THEN
    12501253
     
    12821285!
    12831286!--    Initialization of the turbulence recycling method
    1284        IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.  &
     1287       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.               &
    12851288            turbulent_inflow )  THEN
    12861289!
     
    13461349                inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
    13471350             ELSE
    1348                 WRITE( message_string, * ) 'inflow_damping_height must be ',&
    1349                      'explicitly specified because&the inversion height ', &
     1351                WRITE( message_string, * ) 'inflow_damping_height must be ',   &
     1352                     'explicitly specified because&the inversion height ',     &
    13501353                     'calculated by the prerun is zero.'
    13511354                CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
     
    13821385!
    13831386!--    Inside buildings set velocities and TKE back to zero
    1384        IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.  &
     1387       IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.                &
    13851388            topography /= 'flat' )  THEN
    13861389!
     
    14151418       IF ( humidity  .OR.  passive_scalar )  THEN
    14161419          q_p = q
    1417           IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
     1420          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
    14181421               precipitation )  THEN
    14191422             qr_p = qr
     
    14301433       IF ( humidity  .OR.  passive_scalar )  THEN
    14311434          tq_m = 0.0_wp
    1432           IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
     1435          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
    14331436               precipitation )  THEN
    14341437             tqr_m = 0.0_wp
     
    14861489             DO  j = nys, nyn
    14871490                DO  k = nzb_2d(j,nx)+1, nzt
    1488                    volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
     1491                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
    14891492                                              u_init(k) * dzw(k)
    14901493                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
     
    15221525             DO  j = nys, nyn
    15231526                DO  k = nzb_2d(j,nx)+1, nzt
    1524                    volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
     1527                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
    15251528                                              hom_sum(k,1,0) * dzw(k)
    15261529                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
     
    15321535             DO  i = nxl, nxr
    15331536                DO  k = nzb_2d(ny,i)+1, nzt 
    1534                    volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
     1537                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
    15351538                                              hom_sum(k,2,0) * dzw(k)
    15361539                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
     
    15681571             DO  i = nxl, nxr
    15691572                DO  k = nzb_2d(ny,i)+1, nzt 
    1570                    volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
     1573                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
    15711574                                              v(k,ny,i) * dzw(k)
    15721575                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
     
    16051608!-- Impose random perturbation on the horizontal velocity field and then
    16061609!-- remove the divergences from the velocity field at the initial stage
    1607     IF ( create_disturbances  .AND.  disturbance_energy_limit /= 0.0_wp  .AND.&
    1608          TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
     1610    IF ( create_disturbances  .AND.  disturbance_energy_limit /= 0.0_wp  .AND. &
     1611         TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
    16091612         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
    16101613
     
    17181721    rdf_sc = 0.0_wp
    17191722    IF ( rayleigh_damping_factor /= 0.0_wp )  THEN
    1720        IF ( .NOT. ocean )  THEN
     1723       IF (  .NOT. ocean )  THEN
    17211724          DO  k = nzb+1, nzt
    17221725             IF ( zu(k) >= rayleigh_damping_height )  THEN
    1723                 rdf(k) = rayleigh_damping_factor * &
     1726                rdf(k) = rayleigh_damping_factor *                             &
    17241727                      ( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height ) &
    1725                                          / ( zu(nzt) - rayleigh_damping_height ) )&
     1728                             / ( zu(nzt) - rayleigh_damping_height ) )         &
    17261729                      )**2
    17271730             ENDIF
     
    17301733          DO  k = nzt, nzb+1, -1
    17311734             IF ( zu(k) <= rayleigh_damping_height )  THEN
    1732                 rdf(k) = rayleigh_damping_factor * &
     1735                rdf(k) = rayleigh_damping_factor *                             &
    17331736                      ( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) ) &
    1734                                          / ( rayleigh_damping_height - zu(nzb+1)))&
     1737                             / ( rayleigh_damping_height - zu(nzb+1) ) )       &
    17351738                      )**2
    17361739             ENDIF
     
    17741777             ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp *              &
    17751778                            REAL( pt_damping_width - i * dx, KIND=wp ) / (     &
    1776                             REAL( pt_damping_width, KIND=wp )            ) ) )**2
     1779                            REAL( pt_damping_width, KIND=wp ) ) ) )**2
    17771780          ENDIF
    17781781       ENDDO
     
    18311834
    18321835    IF ( dots_num > dots_max )  THEN
    1833        WRITE( message_string, * ) 'number of time series quantities exceeds', &
    1834                                   ' its maximum of dots_max = ', dots_max,    &
     1836       WRITE( message_string, * ) 'number of time series quantities exceeds',  &
     1837                                  ' its maximum of dots_max = ', dots_max,     &
    18351838                                  ' &Please increase dots_max in modules.f90.'
    18361839       CALL message( 'init_3d_model', 'PA0194', 1, 2, 0, 6, 0 )   
     
    18691872!--             All xy-grid points
    18701873                ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
    1871                 mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) + &
    1872                                                   zw(nzb_s_inner(j,i))
     1874                mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) &
     1875                                                  + zw(nzb_s_inner(j,i))
    18731876!
    18741877!--             xy-grid points above topography
     
    18811884!
    18821885!--             All grid points of the total domain above topography
    1883                 ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) + &
    1884                                      ( nz - nzb_s_inner(j,i) + 2 )
     1886                ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr)                       &
     1887                                     + ( nz - nzb_s_inner(j,i) + 2 )
    18851888             ENDIF
    18861889          ENDDO
     
    18911894#if defined( __parallel )
    18921895    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1893     CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM,   &
     1896    CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM,    &
    18941897                        comm2d, ierr )
    18951898    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1896     CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr,  &
     1899    CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr,   &
    18971900                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
    18981901    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1899     CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0),         &
     1902    CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0),          &
    19001903                        (nz+2)*sr, MPI_INTEGER, MPI_SUM, comm2d, ierr )
    19011904    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1902     CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL, &
     1905    CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL,  &
    19031906                        MPI_SUM, comm2d, ierr )
    19041907    ngp_3d_inner = INT( ngp_3d_inner_tmp, KIND = SELECTED_INT_KIND( 18 ) )
    19051908    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1906     CALL MPI_ALLREDUCE( mean_surface_level_height_l(0),                       &
    1907                         mean_surface_level_height(0), sr, MPI_REAL,           &
     1909    CALL MPI_ALLREDUCE( mean_surface_level_height_l(0),                        &
     1910                        mean_surface_level_height(0), sr, MPI_REAL,            &
    19081911                        MPI_SUM, comm2d, ierr )
    19091912    mean_surface_level_height = mean_surface_level_height / REAL( ngp_2dh )
     
    19241927!-- the respective subdomain lie below the surface topography
    19251928    ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   )
    1926     ngp_3d_inner    = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )),            &
     1929    ngp_3d_inner    = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )),             &
    19271930                           ngp_3d_inner(:) )
    19281931    ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) )
    19291932
    1930     DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l,      &
     1933    DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l,       &
    19311934                ngp_3d_inner_l, ngp_3d_inner_tmp )
    19321935
  • palm/trunk/SOURCE/land_surface_model.f90

    r1784 r1788  
    1919! Current revisions:
    2020! -----------------
     21! Bugfix: calculate lambda_surface based on temperature gradient between skin
     22! layer and soil layer instead of Obukhov length
     23! Changed: moved calculation of surface specific humidity to energy balance solver
     24! New: water surfaces are available by using a fixed sea surface temperature.
     25! The roughness lengths are calculated dynamically using the Charnock
     26! parameterization. This involves the new roughness length for moisture z0q.
     27! New: modified solution of the energy balance solver and soil model for
     28! paved surfaces (i.e. asphalt concrete).
     29! Syntax layout improved.
     30! Changed: parameter dewfall removed.
    2131!
    22 !
    2332! Former revisions:
    2433! -----------------
     
    105114!> @todo Consider partial absorption of the net shortwave radiation by the
    106115!>       skin layer.
    107 !> @todo Allow for water surfaces
     116!> @todo Improve surface water parameterization
    108117!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0,
    109118!>       nzt_soil=3)).
     
    118127 
    119128    USE arrays_3d,                                                             &
    120         ONLY:  hyp, ol, pt, pt_p, q, q_p, ql, qsws, shf, ts, us, vpt, z0, z0h
     129        ONLY:  hyp, ol, pt, pt_p, q, q_p, ql, qsws, shf, ts, us, vpt, z0, z0h, &
     130               z0q
    121131
    122132    USE cloud_parameters,                                                      &
     
    171181!
    172182!-- LSM variables
    173     INTEGER(iwp) :: veg_type  = 2, & !< vegetation type, 0: user-defined, 1-19: generic (see list)
    174                     soil_type = 3    !< soil type, 0: user-defined, 1-6: generic (see list)
     183    INTEGER(iwp) :: veg_type  = 2, & !< NAMELIST veg_type_2d
     184                    soil_type = 3    !< NAMELIST soil_type_2d
     185
     186    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  soil_type_2d, &  !< soil type, 0: user-defined, 1-7: generic (see list)
     187                                                  veg_type_2d      !< vegetation type, 0: user-defined, 1-19: generic (see list)
     188
     189    LOGICAL, DIMENSION(:,:), ALLOCATABLE :: water_surface,     & !< flag parameter for water surfaces (classes 14+15)
     190                                            pave_surface,      & !< flag parameter for pavements (asphalt etc.) (class 20)
     191                                            building_surface     !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet)
    175192
    176193    LOGICAL :: conserve_water_content = .TRUE.,  & !< open or closed bottom surface for the soil model
    177                dewfall = .TRUE.,                 & !< allow/inhibit dewfall
    178194               force_radiation_call_l = .FALSE., & !< flag parameter for unscheduled radiation model calls
    179195               land_surface = .FALSE.              !< flag parameter indicating wheather the lsm is used
     
    200216                m_total = 0.0_wp,                       & !< weighted total water content of the soil (m3/m3)
    201217                n_vangenuchten = 9999999.9_wp,          & !< NAMELIST n_vg
     218                pave_depth = 9999999.9_wp,              & !< depth of the pavement
     219                pave_heat_capacity = 1.94E6_wp,         & !< volumetric heat capacity of pavement (e.g. roads)
     220                pave_heat_conductivity = 1.00_wp,       & !< heat conductivity for pavements (e.g. roads)
    202221                q_s = 0.0_wp,                           & !< saturation specific humidity
    203222                residual_moisture = 9999999.9_wp,       & !< NAMELIST m_res
     
    210229                wilting_point = 9999999.9_wp,           & !< NAMELIST m_wilt
    211230                z0_eb  = 9999999.9_wp,                  & !< NAMELIST z0 (lsm_par)
    212                 z0h_eb = 9999999.9_wp                     !< NAMELIST z0h (lsm_par)
     231                z0h_eb = 9999999.9_wp,                  & !< NAMELIST z0h (lsm_par)
     232                z0q_eb = 9999999.9_wp                     !< NAMELIST z0q (lsm_par)
    213233
    214234    REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: &
     
    295315
    296316    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
    297               lambda_h, &   !< heat conductivity of soil (?)                           
     317              lambda_h, &   !< heat conductivity of soil (W/m/K)                           
    298318              lambda_w, &   !< hydraulic diffusivity of soil (?)
    299               gamma_w,  &   !< hydraulic conductivity of soil (?)
     319              gamma_w,  &   !< hydraulic conductivity of soil (W/m/K)
    300320              rho_c_total   !< volumetric heat capacity of the actual soil matrix (?)
    301321
     
    327347!
    328348!-- Predefined Land surface classes (veg_type)
    329     CHARACTER(26), DIMENSION(0:19), PARAMETER :: veg_type_name = (/ &
     349    CHARACTER(26), DIMENSION(0:20), PARAMETER :: veg_type_name = (/ &
    330350                                   'user defined              ',    & ! 0
    331351                                   'crops, mixed farming      ',    & !  1
     
    347367                                   'deciduous shrubs          ',    & ! 17
    348368                                   'mixed forest/woodland     ',    & ! 18
    349                                    'interrupted forest        '     & ! 19
     369                                   'interrupted forest        ',    & ! 19
     370                                   'pavements/roads           '     & ! 20
    350371                                                                 /)
    351372
     
    368389!-- Land surface parameters I
    369390!--                          r_canopy_min,     lai,   c_veg,     g_d
    370     REAL(wp), DIMENSION(0:3,1:19), PARAMETER :: veg_pars = RESHAPE( (/ &
     391    REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: veg_pars = RESHAPE( (/ &
    371392                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & !  1
    372393                                 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,  & !  2
     
    387408                                 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,  & ! 17
    388409                                 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & ! 18
    389                                  175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp   & ! 19
    390                                  /), (/ 4, 19 /) )
    391 
    392 !
    393 !-- Land surface parameters II          z0,         z0h
    394     REAL(wp), DIMENSION(0:1,1:19), PARAMETER :: roughness_par = RESHAPE( (/ &
    395                                    0.25_wp,  0.25E-2_wp,                    & !  1
    396                                    0.20_wp,  0.20E-2_wp,                    & !  2
    397                                    2.00_wp,     2.00_wp,                    & !  3
    398                                    2.00_wp,     2.00_wp,                    & !  4
    399                                    2.00_wp,     2.00_wp,                    & !  5
    400                                    2.00_wp,     2.00_wp,                    & !  6
    401                                    0.47_wp,  0.47E-2_wp,                    & !  7
    402                                   0.013_wp, 0.013E-2_wp,                    & !  8
    403                                   0.034_wp, 0.034E-2_wp,                    & !  9
    404                                     0.5_wp,  0.50E-2_wp,                    & ! 10
    405                                    0.17_wp,  0.17E-2_wp,                    & ! 11
    406                                  1.3E-3_wp,   1.3E-4_wp,                    & ! 12
    407                                    0.83_wp,  0.83E-2_wp,                    & ! 13
    408                                    0.00_wp,  0.00E-2_wp,                    & ! 14
    409                                    0.00_wp,  0.00E-2_wp,                    & ! 15
    410                                    0.10_wp,  0.10E-2_wp,                    & ! 16
    411                                    0.25_wp,  0.25E-2_wp,                    & ! 17
    412                                    2.00_wp,  2.00E-2_wp,                    & ! 18
    413                                    1.10_wp,  1.10E-2_wp                     & ! 19
    414                                  /), (/ 2, 19 /) )
     410                                 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp,  & ! 19
     411                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp   & ! 20
     412                                 /), (/ 4, 20 /) )
     413
     414!
     415!-- Land surface parameters II          z0,         z0h,         z0q
     416    REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: roughness_par = RESHAPE( (/ &
     417                                   0.25_wp,  0.25E-2_wp,  0.25E-2_wp,       & !  1
     418                                   0.20_wp,  0.20E-2_wp,  0.20E-2_wp,       & !  2
     419                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  3
     420                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  4
     421                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  5
     422                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  6
     423                                   0.47_wp,  0.47E-2_wp,  0.47E-2_wp,       & !  7
     424                                  0.013_wp, 0.013E-2_wp, 0.013E-2_wp,       & !  8
     425                                  0.034_wp, 0.034E-2_wp, 0.034E-2_wp,       & !  9
     426                                    0.5_wp,  0.50E-2_wp,  0.50E-2_wp,       & ! 10
     427                                   0.17_wp,  0.17E-2_wp,  0.17E-2_wp,       & ! 11
     428                                 1.3E-3_wp,   1.3E-4_wp,   1.3E-4_wp,       & ! 12
     429                                   0.83_wp,  0.83E-2_wp,  0.83E-2_wp,       & ! 13
     430                                   0.00_wp,     0.00_wp,     0.00_wp,       & ! 14
     431                                   0.00_wp,     0.00_wp,     0.00_wp,       & ! 15
     432                                   0.10_wp,  0.10E-2_wp,  0.10E-2_wp,       & ! 16
     433                                   0.25_wp,  0.25E-2_wp,  0.25E-2_wp,       & ! 17
     434                                   2.00_wp,  2.00E-2_wp,  2.00E-2_wp,       & ! 18
     435                                   1.10_wp,  1.10E-2_wp,  1.10E-2_wp,       & ! 19
     436                                 1.0E-4_wp,   1.0E-5_wp,   1.0E-5_wp        & ! 20
     437                                 /), (/ 3, 20 /) )
    415438
    416439!
    417440!-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
    418     REAL(wp), DIMENSION(0:2,1:19), PARAMETER :: surface_pars = RESHAPE( (/ &
     441    REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: surface_pars = RESHAPE( (/ &
    419442                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  1
    420443                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  2
     
    430453                                      58.0_wp,       58.0_wp, 0.00_wp,     & ! 12
    431454                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 13
    432                                     1.0E20_wp,     1.0E20_wp, 0.00_wp,     & ! 14
    433                                     1.0E20_wp,     1.0E20_wp, 0.00_wp,     & ! 15
     455                                    1.0E10_wp,     1.0E10_wp, 0.00_wp,     & ! 14
     456                                    1.0E10_wp,     1.0E10_wp, 0.00_wp,     & ! 15
    434457                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 16
    435458                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 17
    436459                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 18
    437                                       20.0_wp,       15.0_wp, 0.03_wp      & ! 19
    438                                       /), (/ 3, 19 /) )
     460                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 19
     461                                       0.0_wp,        0.0_wp, 0.00_wp      & ! 20
     462                                      /), (/ 3, 20 /) )
    439463
    440464!
    441465!-- Root distribution (sum = 1)  level 1, level 2, level 3, level 4,
    442     REAL(wp), DIMENSION(0:3,1:19), PARAMETER :: root_distribution = RESHAPE( (/ &
     466    REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: root_distribution = RESHAPE( (/ &
    443467                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  1
    444468                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  2
     
    459483                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 17
    460484                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 18
    461                                  0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp             & ! 19
    462                                  /), (/ 4, 19 /) )
     485                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 19
     486                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp             & ! 20
     487                                 /), (/ 4, 20 /) )
    463488
    464489!
     
    499524!-- Public parameters, constants and initial values
    500525    PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient,       &
    501            conserve_water_content, dewfall, field_capacity,                    &
     526           conserve_water_content, field_capacity,                             &
    502527           f_shortwave_incoming, hydraulic_conductivity, init_lsm,             &
    503528           init_lsm_arrays, lambda_surface_stable, lambda_surface_unstable,    &
    504529           land_surface, leaf_area_index, lsm_energy_balance, lsm_soil_model,  &
    505530           lsm_swap_timelevel, l_vangenuchten, min_canopy_resistance,          &
    506            min_soil_resistance, n_vangenuchten, residual_moisture, rho_cp,     &
     531           min_soil_resistance, n_vangenuchten, pave_heat_capacity,            &
     532           pave_depth, pave_heat_conductivity, residual_moisture, rho_cp,      &
    507533           rho_lv, root_fraction, saturation_moisture, skip_time_do_lsm,       &
    508534           soil_moisture, soil_temperature, soil_type, soil_type_name,         &
    509535           vegetation_coverage, veg_type, veg_type_name, wilting_point, z0_eb, &
    510            z0h_eb
     536           z0h_eb, z0q_eb
    511537
    512538!
     
    586612!--    Allocate 2D vegetation model arrays
    587613       ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) )
     614       ALLOCATE ( building_surface(nysg:nyng,nxlg:nxrg) )
    588615       ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
    589616       ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
     
    601628       ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
    602629       ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) )
     630       ALLOCATE ( pave_surface(nysg:nyng,nxlg:nxrg) )
    603631       ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) )
    604632       ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) )
     
    613641       ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) )
    614642       ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) )
     643       ALLOCATE ( soil_type_2d(nysg:nyng,nxlg:nxrg) )
     644       ALLOCATE ( veg_type_2d(nysg:nyng,nxlg:nxrg) )
     645       ALLOCATE ( water_surface(nysg:nyng,nxlg:nxrg) )
    615646
    616647#if ! defined( __nopointer )
     
    650681!
    651682!--    If no cloud physics is used, rho_surface has not been calculated before
    652        IF ( .NOT. cloud_physics )  THEN
     683       IF (  .NOT. cloud_physics )  THEN
    653684          rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
    654685       ENDIF
     
    762793       ENDIF   
    763794
     795!
     796!--    Map values to the respective 2D arrays
    764797       alpha_vg      = alpha_vangenuchten
    765798       l_vg          = l_vangenuchten
     
    796829
    797830!
    798 !--       Set artifical values for ts and us so that r_a has its initial value for
    799 !--       the first time step
     831!--       Set artifical values for ts and us so that r_a has its initial value
     832!--       for the first time step
    800833          DO  i = nxlg, nxrg
    801834             DO  j = nysg, nyng
     
    864897             z0h_eb = roughness_par(1,veg_type)
    865898          ENDIF
    866           z0h_factor = z0h_eb / z0_eb
     899          IF ( z0q_eb == 9999999.9_wp )  THEN
     900             z0q_eb = roughness_par(2,veg_type)
     901          ENDIF
     902          z0h_factor = z0h_eb / ( z0_eb + 1.0E-20_wp )
    867903
    868904          IF ( ANY( root_fraction == 9999999.9_wp ) )  THEN
     
    881917             z0h_eb = z0_eb * z0h_factor
    882918          ENDIF
     919          IF ( z0q_eb == 9999999.9_wp )  THEN
     920             z0q_eb = z0_eb * z0h_factor
     921          ENDIF
    883922
    884923       ENDIF
    885924
    886925!
    887 !--    Initialize vegetation
     926!--    For surfaces covered with pavement, set depth of the pavement (with dry
     927!--    soil below). The depth must be greater than the first soil layer depth
     928       IF ( veg_type == 20 )  THEN
     929          IF ( pave_depth == 9999999.9_wp )  THEN
     930             pave_depth = zs(nzb_soil) 
     931          ELSE
     932             pave_depth = MAX( zs(nzb_soil), pave_depth )
     933          ENDIF
     934       ENDIF
     935
     936!
     937!--    Map vegetation and soil types to 2D array to allow for heterogeneous
     938!--    surfaces via user interface see below
     939       veg_type_2d = veg_type
     940       soil_type_2d = soil_type
     941
     942!
     943!--    Map vegetation parameters to the respective 2D arrays
    888944       r_canopy_min         = min_canopy_resistance
    889945       lai                  = leaf_area_index
     
    895951       z0                   = z0_eb
    896952       z0h                  = z0h_eb
     953       z0q                  = z0q_eb
    897954
    898955!
     
    900957       CALL user_init_land_surface
    901958
     959!
     960!--    Set flag parameter if vegetation type was set to a water surface. Also
     961!--    set temperature to a constant value in all "soil" layers.
     962       DO  i = nxlg, nxrg
     963          DO  j = nysg, nyng
     964             IF ( veg_type_2d(j,i) == 14  .OR.  veg_type_2d(j,i) == 15 )  THEN
     965                water_surface(j,i) = .TRUE.
     966             ELSEIF ( veg_type_2d(j,i) == 20 )  THEN
     967                pave_surface(j,i) = .TRUE.
     968                m_soil(:,j,i) = 0.0_wp
     969             ENDIF
     970
     971          ENDDO
     972       ENDDO
     973
     974!
     975!--    Calculate new roughness lengths (for water surfaces only)
     976       CALL calc_z0_water_surface
    902977
    903978       t_soil_p    = t_soil
     
    9351010       INTEGER(iwp) ::  k, ks     !< running index
    9361011
    937        REAL(wp) :: f1,          & !< resistance correction term 1
     1012       REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
     1013                   f1,          & !< resistance correction term 1
    9381014                   f2,          & !< resistance correction term 2
    9391015                   f3,          & !< resistance correction term 3
     
    9651041
    9661042!
    967 !--          Set lambda_surface according to stratification
    968              IF ( ol(j,i) >= 0.0_wp )  THEN
    969                 lambda_surface = lambda_surface_s(j,i)
     1043!--          Set lambda_surface according to stratification between skin layer and soil
     1044             IF (  .NOT.  pave_surface(j,i) )  THEN
     1045
     1046                c_surface_tmp = c_surface
     1047
     1048                IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i))  THEN
     1049                   lambda_surface = lambda_surface_s(j,i)
     1050                ELSE
     1051                   lambda_surface = lambda_surface_u(j,i)
     1052                ENDIF
    9701053             ELSE
    971                 lambda_surface = lambda_surface_u(j,i)
     1054
     1055                c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
     1056                lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
     1057
    9721058             ENDIF
    9731059
     
    10161102             ENDDO
    10171103
    1018              IF ( m_total > m_wilt(j,i) .AND. m_total < m_fc(j,i) )  THEN
     1104             IF ( m_total > m_wilt(j,i)  .AND. m_total < m_fc(j,i) )  THEN
    10191105                f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
    10201106             ELSEIF ( m_total >= m_fc(j,i) )  THEN
     
    10501136!--          Third step: calculate bare soil resistance r_soil. The Clapp &
    10511137!--          Hornberger parametrization does not consider c_veg.
    1052              IF ( soil_type /= 7 )  THEN
     1138             IF ( soil_type_2d(j,i) /= 7 )  THEN
    10531139                m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *     &
    10541140                        m_res(j,i)
     
    10641150
    10651151!
    1066 !--          Calculate fraction of liquid water reservoir
    1067              m_liq_eb_max = m_max_depth * lai(j,i)
    1068              c_liq(j,i) = MIN(1.0_wp, m_liq_eb(j,i) / (m_liq_eb_max+1.0E-20_wp))
    1069              
     1152!--          Calculate the maximum possible liquid water amount on plants and
     1153!--          bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
     1154!--          assumed, while paved surfaces might hold up 1 mm of water. The
     1155!--          liquid water fraction for paved surfaces is calculated after
     1156!--          Noilhan & Planton (1989), while the ECMWF formulation is used for
     1157!--          vegetated surfaces and bare soils.
     1158             IF ( pave_surface(j,i) )  THEN
     1159                m_liq_eb_max = m_max_depth * 5.0_wp
     1160                c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 )
     1161             ELSE
     1162                m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i)           &
     1163                               + (1.0_wp - c_veg(j,i)) )
     1164                c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max )
     1165             ENDIF
    10701166
    10711167!
     
    10761172!--          In case of dewfall, set evapotranspiration to zero
    10771173!--          All super-saturated water is then removed from the air
    1078              IF ( humidity .AND. q_s <= qv1 )  THEN
     1174             IF ( humidity  .AND. q_s <= qv1 )  THEN
    10791175                r_canopy(j,i) = 0.0_wp
    10801176                r_soil(j,i)   = 0.0_wp
     
    10831179!
    10841180!--          Calculate coefficients for the total evapotranspiration
    1085              f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/        &
    1086                            (r_a(j,i) + r_canopy(j,i))
    1087 
    1088              f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +        &
     1181!--          In case of water surface, set vegetation and soil fluxes to zero.
     1182!--          For pavements, only evaporation of liquid water is possible.
     1183             IF ( water_surface(j,i) )  THEN
     1184                f_qsws_veg  = 0.0_wp
     1185                f_qsws_soil = 0.0_wp
     1186                f_qsws_liq  = rho_lv / r_a(j,i)
     1187             ELSEIF ( pave_surface (j,i) )  THEN
     1188                f_qsws_veg  = 0.0_wp
     1189                f_qsws_soil = 0.0_wp
     1190                f_qsws_liq  = rho_lv * c_liq(j,i) / r_a(j,i)
     1191             ELSE
     1192                f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/     &
     1193                              (r_a(j,i) + r_canopy(j,i))
     1194                f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +     &
    10891195                                                          r_soil(j,i))
    1090              f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
    1091 
    1092 
     1196                f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
     1197             ENDIF
    10931198!
    10941199!--          If soil moisture is below wilting point, plants do no longer
     
    10981203!              ENDIF
    10991204
    1100 !
    1101 !--          If dewfall is deactivated, vegetation, soil and liquid water
    1102 !--          reservoir are not allowed to take up water from the super-saturated
    1103 !--          air.
    1104              IF ( humidity .AND. q_s <= qv1 .AND. .NOT. dewfall )  THEN
    1105                 f_qsws_veg  = 0.0_wp
    1106                 f_qsws_soil = 0.0_wp
    1107                 f_qsws_liq  = 0.0_wp
    1108              ENDIF
    1109 
    11101205             f_shf  = rho_cp / r_a(j,i)
    11111206             f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
     
    11231218             rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
    11241219
     1220!
     1221!--          Calculate new skin temperature
    11251222             IF ( humidity )  THEN
    11261223#if defined ( __rrtmg )
     
    11861283!--          Implicit solution when the surface layer has no heat capacity,
    11871284!--          otherwise use RK3 scheme.
    1188              t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface *        &
    1189                                 t_surface(j,i) ) / ( c_surface + coef_2 * dt_3d&
    1190                                 * tsc(2) )
     1285             t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *    &
     1286                                t_surface(j,i) ) / ( c_surface_tmp + coef_2    &
     1287                                * dt_3d * tsc(2) )
     1288
    11911289!
    11921290!--          Add RK3 term
    1193              IF ( c_surface /= 0.0_wp )  THEN
     1291             IF ( c_surface_tmp /= 0.0_wp )  THEN
    11941292
    11951293                t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)           &
     
    12111309                   ENDIF
    12121310                ENDIF
    1213 
    12141311             ENDIF
    12151312
     
    12231320!--          often reached, when no oscillations would occur (causes immense
    12241321!--          computing time for the radiation code).
    1225              IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp .AND.      &
     1322             IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp  .AND.     &
    12261323                  unscheduled_radiation_calls )  THEN
    12271324                force_radiation_call_l = .TRUE.
     
    12491346#endif
    12501347
    1251 
    12521348             ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
    12531349                              - t_soil(nzb_soil,j,i))
     
    12751371                                    * t_surface_p(j,i) )
    12761372             ENDIF
     1373
    12771374!
    12781375!--          Calculate the true surface resistance
     
    12801377                r_s(j,i) = 1.0E10_wp
    12811378             ELSE
    1282                 r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                   &
     1379                r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                    &
    12831380                           * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )     &
    12841381                           / qsws_eb(j,i) - r_a(j,i)
     
    12971394!
    12981395!--                Add precipitation to liquid water reservoir, if possible.
    1299 !--                Otherwise, add the water to bare soil fraction.
     1396!--                Otherwise, add the water to soil. In case of
     1397!--                pavements, the exceeding water amount is implicitely removed
     1398!--                as runoff as qsws_soil_eb is then not used in the soil model
    13001399                   IF ( m_liq_eb(j,i) /= m_liq_eb_max )  THEN
    13011400                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                      &
     
    13761475!--    Make a logical OR for all processes. Force radiation call if at
    13771476!--    least one processor reached the threshold change in skin temperature
    1378        IF ( unscheduled_radiation_calls .AND. intermediate_timestep_count      &
     1477       IF ( unscheduled_radiation_calls  .AND.  intermediate_timestep_count    &
    13791478            == intermediate_timestep_count_max-1 )  THEN
    13801479#if defined( __parallel )
     
    13881487       ENDIF
    13891488
     1489!
     1490!--    Calculate surface specific humidity
     1491       IF ( humidity )  THEN
     1492          CALL calc_q_surface
     1493       ENDIF
     1494
     1495!
     1496!--    Calculate new roughness lengths (for water surfaces only)
     1497       CALL calc_z0_water_surface
     1498
    13901499
    13911500    END SUBROUTINE lsm_energy_balance
     
    14151524       DO  i = nxlg, nxrg
    14161525          DO  j = nysg, nyng
    1417              DO  k = nzb_soil, nzt_soil
    1418 !
    1419 !--             Calculate volumetric heat capacity of the soil, taking into
    1420 !--             account water content
    1421                 rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i))       &
    1422                                      + rho_c_water * m_soil(k,j,i))
    1423 
    1424 !
    1425 !--             Calculate soil heat conductivity at the center of the soil
    1426 !--             layers
    1427                 lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) *          &
    1428                                lambda_h_water ** m_soil(k,j,i)
    1429 
    1430                 ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
    1431 
    1432                 lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +          &
    1433                                  lambda_h_dry
    1434 
    1435              ENDDO
    1436 
    1437 !
    1438 !--          Calculate soil heat conductivity (lambda_h) at the _stag level
    1439 !--          using linear interpolation
    1440              DO  k = nzb_soil, nzt_soil-1
    1441                 lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) * 0.5_wp
    1442              ENDDO
    1443              lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
    1444 
    1445 !
    1446 !--          Prognostic equation for soil temperature t_soil
    1447              tend(:) = 0.0_wp
    1448              tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *             &
    1449                        ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)     &
    1450                          - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1)       &
    1451                          + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
    1452 
    1453 
    1454 
    1455 
    1456              DO  k = nzb_soil+1, nzt_soil
    1457                 tend(k) = (1.0_wp/rho_c_total(k,j,i))                          &
    1458                           * (   lambda_h(k,j,i)                                &
    1459                               * ( t_soil(k+1,j,i) - t_soil(k,j,i) )            &
    1460                               * ddz_soil(k+1)                                  &
    1461                               - lambda_h(k-1,j,i)                              &
    1462                               * ( t_soil(k,j,i) - t_soil(k-1,j,i) )            &
    1463                               * ddz_soil(k)                                    &
    1464                             ) * ddz_soil_stag(k)
    1465              ENDDO
    1466 
    1467              t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)   &
    1468                                              + dt_3d * ( tsc(2)                &
    1469                                              * tend(:) + tsc(3)                &
    1470                                              * tt_soil_m(:,j,i) )   
    1471 
    1472 !
    1473 !--          Calculate t_soil tendencies for the next Runge-Kutta step
    1474              IF ( timestep_scheme(1:5) == 'runge' )  THEN
    1475                 IF ( intermediate_timestep_count == 1 )  THEN
    1476                    DO  k = nzb_soil, nzt_soil
    1477                       tt_soil_m(k,j,i) = tend(k)
    1478                    ENDDO
    1479                 ELSEIF ( intermediate_timestep_count <                         &
    1480                          intermediate_timestep_count_max )  THEN
    1481                    DO  k = nzb_soil, nzt_soil
    1482                       tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp      &
    1483                                          * tt_soil_m(k,j,i)
    1484                    ENDDO
    1485                 ENDIF
    1486              ENDIF
    1487 
    1488 
    1489              DO  k = nzb_soil, nzt_soil
    1490 
    1491 !
    1492 !--             Calculate soil diffusivity at the center of the soil layers
    1493                 lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat          &
    1494                                   / m_sat(j,i) ) * ( MAX(m_soil(k,j,i),        &
    1495                                   m_wilt(j,i)) / m_sat(j,i) )**(b_ch + 2.0_wp)
    1496 
    1497 !
    1498 !--             Parametrization of Van Genuchten
    1499                 IF ( soil_type /= 7 )  THEN
    1500 !
    1501 !--                Calculate the hydraulic conductivity after Van Genuchten
    1502 !--                (1980)
    1503                    h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -       &
    1504                               MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_vg(j,i)   &
    1505                               / (n_vg(j,i)-1.0_wp)) - 1.0_wp                   &
    1506                           )**(1.0_wp/n_vg(j,i)) / alpha_vg(j,i)
    1507 
    1508 
    1509                    gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +            &
    1510                                    (alpha_vg(j,i)*h_vg)**n_vg(j,i))**(1.0_wp   &
    1511                                    -1.0_wp/n_vg(j,i)) - (alpha_vg(j,i)*h_vg    &
    1512                                    )**(n_vg(j,i)-1.0_wp))**2 )                 &
    1513                                    / ( (1.0_wp + (alpha_vg(j,i)*h_vg           &
    1514                                    )**n_vg(j,i))**((1.0_wp - 1.0_wp/n_vg(j,i)) &
    1515                                    *(l_vg(j,i) + 2.0_wp)) )
    1516 
    1517 !
    1518 !--             Parametrization of Clapp & Hornberger
    1519                 ELSE
    1520                    gamma_temp(k) = gamma_w_sat(j,i) * (m_soil(k,j,i)           &
    1521                                    / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
    1522                 ENDIF
    1523 
    1524              ENDDO
    1525 
    1526 
    1527              IF ( humidity )  THEN
    1528 !
    1529 !--             Calculate soil diffusivity (lambda_w) at the _stag level
    1530 !--             using linear interpolation. To do: replace this with
    1531 !--             ECMWF-IFS Eq. 8.81
     1526
     1527             IF ( pave_surface(j,i) )  THEN
     1528                rho_c_total(nzb_soil,j,i) = pave_heat_capacity
     1529                lambda_temp(nzb_soil)     = pave_heat_conductivity
     1530             ENDIF
     1531
     1532             IF (  .NOT.  water_surface(j,i) )  THEN
     1533                DO  k = nzb_soil, nzt_soil
     1534
     1535
     1536                   IF ( pave_surface(j,i)  .AND.  zs(k) <= pave_depth )  THEN
     1537                   
     1538                      rho_c_total(k,j,i) = pave_heat_capacity
     1539                      lambda_temp(k)     = pave_heat_conductivity   
     1540
     1541                   ELSE           
     1542!
     1543!--                   Calculate volumetric heat capacity of the soil, taking
     1544!--                   into account water content
     1545                      rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) &
     1546                                           + rho_c_water * m_soil(k,j,i))
     1547
     1548!
     1549!--                   Calculate soil heat conductivity at the center of the soil
     1550!--                   layers
     1551                      lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) *    &
     1552                                     lambda_h_water ** m_soil(k,j,i)
     1553
     1554                      ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
     1555
     1556                      lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +    &
     1557                                       lambda_h_dry
     1558                   ENDIF
     1559
     1560                ENDDO
     1561
     1562!
     1563!--             Calculate soil heat conductivity (lambda_h) at the _stag level
     1564!--             using linear interpolation. For pavement surface, the
     1565!--             true pavement depth is considered
    15321566                DO  k = nzb_soil, nzt_soil-1
    1533                      
    1534                    lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )     &
    1535                                      * 0.5_wp
    1536                    gamma_w(k,j,i)  = ( gamma_temp(k+1) + gamma_temp(k) )       &
    1537                                      * 0.5_wp
    1538 
     1567                   IF ( pave_surface(j,i)  .AND.  zs(k)   < pave_depth         &
     1568                                           .AND.  zs(k+1) > pave_depth )  THEN
     1569                      lambda_h(k,j,i) = ( pave_depth - zs(k) ) / dz_soil(k+1)  &
     1570                                        * lambda_temp(k)                       &
     1571                                        + ( 1.0_wp - ( pave_depth - zs(k) )    &
     1572                                        / dz_soil(k+1) ) * lambda_temp(k+1)
     1573                   ELSE
     1574                      lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
     1575                                        * 0.5_wp
     1576                   ENDIF
    15391577                ENDDO
    1540 
    1541 !
    1542 !
    1543 !--             In case of a closed bottom (= water content is conserved), set
    1544 !--             hydraulic conductivity to zero to that no water will be lost
    1545 !--             in the bottom layer.
    1546                 IF ( conserve_water_content )  THEN
    1547                    gamma_w(nzt_soil,j,i) = 0.0_wp
    1548                 ELSE
    1549                    gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil)
    1550                 ENDIF     
    1551 
    1552 !--             The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v))
    1553 !--             ensures the mass conservation for water. The transpiration of
    1554 !--             plants equals the cumulative withdrawals by the roots in the
    1555 !--             soil. The scheme takes into account the availability of water
    1556 !--             in the soil layers as well as the root fraction in the
    1557 !--             respective layer. Layer with moisture below wilting point will
    1558 !--             not contribute, which reflects the preference of plants to
    1559 !--             take water from moister layers.
    1560 
    1561 !
    1562 !--             Calculate the root extraction (ECMWF 7.69, the sum of root_extr
    1563 !--             = 1). The energy balance solver guarantees a positive
    1564 !--             transpiration, so that there is no need for an additional check.
    1565                 DO  k = nzb_soil, nzt_soil
    1566                     IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
    1567                        m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
    1568                     ENDIF
    1569                 ENDDO 
    1570 
    1571                 IF ( m_total > 0.0_wp )  THEN
    1572                    DO  k = nzb_soil, nzt_soil
    1573                       IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
    1574                          root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total
    1575                       ELSE
    1576                          root_extr(k) = 0.0_wp
    1577                       ENDIF
    1578                    ENDDO
    1579                 ENDIF
    1580 
    1581 !
    1582 !--             Prognostic equation for soil water content m_soil
     1578                lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
     1579
     1580
     1581
     1582
     1583!
     1584!--             Prognostic equation for soil temperature t_soil
    15831585                tend(:) = 0.0_wp
    1584                 tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (                  &
    1585                             m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
    1586                             * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( &
    1587                             root_extr(nzb_soil) * qsws_veg_eb(j,i)             &
    1588                             + qsws_soil_eb(j,i) ) * drho_l_lv )                &
    1589                             * ddz_soil_stag(nzb_soil)
    1590 
    1591                 DO  k = nzb_soil+1, nzt_soil-1
    1592                    tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)             &
    1593                              - m_soil(k,j,i) ) * ddz_soil(k+1) - gamma_w(k,j,i)&
    1594                              - lambda_w(k-1,j,i) * (m_soil(k,j,i) -            &
    1595                              m_soil(k-1,j,i)) * ddz_soil(k)                    &
    1596                              + gamma_w(k-1,j,i) - (root_extr(k)                &
    1597                              * qsws_veg_eb(j,i) * drho_l_lv)                   &
    1598                              ) * ddz_soil_stag(k)
     1586
     1587                tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *          &
     1588                          ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i) &
     1589                            - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1)    &
     1590                            + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
     1591
     1592                DO  k = nzb_soil+1, nzt_soil
     1593                   tend(k) = (1.0_wp/rho_c_total(k,j,i))                       &
     1594                             * (   lambda_h(k,j,i)                             &
     1595                                 * ( t_soil(k+1,j,i) - t_soil(k,j,i) )         &
     1596                                 * ddz_soil(k+1)                               &
     1597                                 - lambda_h(k-1,j,i)                           &
     1598                                 * ( t_soil(k,j,i) - t_soil(k-1,j,i) )         &
     1599                                 * ddz_soil(k)                                 &
     1600                               ) * ddz_soil_stag(k)
    15991601
    16001602                ENDDO
    1601                 tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                     &
    1602                                         - lambda_w(nzt_soil-1,j,i)             &
    1603                                         * (m_soil(nzt_soil,j,i)                &
    1604                                         - m_soil(nzt_soil-1,j,i))              &
    1605                                         * ddz_soil(nzt_soil)                   &
    1606                                         + gamma_w(nzt_soil-1,j,i) - (          &
    1607                                           root_extr(nzt_soil)                  &
    1608                                         * qsws_veg_eb(j,i) * drho_l_lv  )      &
    1609                                       ) * ddz_soil_stag(nzt_soil)             
    1610 
    1611                 m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
    1612                                                 + dt_3d * ( tsc(2) * tend(:)   &
    1613                                                 + tsc(3) * tm_soil_m(:,j,i) )   
    1614 
    1615 !
    1616 !--             Account for dry soils (find a better solution here!)
    1617                 DO  k = nzb_soil, nzt_soil
    1618                    IF ( m_soil_p(k,j,i) < 0.0_wp )  m_soil_p(k,j,i) = 0.0_wp
    1619                 ENDDO
    1620 
    1621 
    1622 
    1623 
    1624 !
    1625 !--             Calculate m_soil tendencies for the next Runge-Kutta step
     1603
     1604                t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)&
     1605                                                  + dt_3d * ( tsc(2)           &
     1606                                                  * tend(nzb_soil:nzt_soil)    &
     1607                                                  + tsc(3)                     &
     1608                                                  * tt_soil_m(:,j,i) )   
     1609
     1610!
     1611!--             Calculate t_soil tendencies for the next Runge-Kutta step
    16261612                IF ( timestep_scheme(1:5) == 'runge' )  THEN
    16271613                   IF ( intermediate_timestep_count == 1 )  THEN
    16281614                      DO  k = nzb_soil, nzt_soil
    1629                          tm_soil_m(k,j,i) = tend(k)
     1615                         tt_soil_m(k,j,i) = tend(k)
    16301616                      ENDDO
    16311617                   ELSEIF ( intermediate_timestep_count <                      &
    16321618                            intermediate_timestep_count_max )  THEN
    16331619                      DO  k = nzb_soil, nzt_soil
    1634                          tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
    1635                                      * tm_soil_m(k,j,i)
     1620                         tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
     1621                                         * tt_soil_m(k,j,i)
    16361622                      ENDDO
    16371623                   ENDIF
    16381624                ENDIF
    16391625
     1626
     1627                DO  k = nzb_soil, nzt_soil
     1628
     1629!
     1630!--                Calculate soil diffusivity at the center of the soil layers
     1631                   lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat       &
     1632                                     / m_sat(j,i) ) * ( MAX( m_soil(k,j,i),    &
     1633                                     m_wilt(j,i) ) / m_sat(j,i) )**(           &
     1634                                     b_ch + 2.0_wp )
     1635
     1636!
     1637!--                Parametrization of Van Genuchten
     1638                   IF ( soil_type /= 7 )  THEN
     1639!
     1640!--                   Calculate the hydraulic conductivity after Van Genuchten
     1641!--                   (1980)
     1642                      h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -    &
     1643                                 MAX( m_soil(k,j,i), m_wilt(j,i) ) ) )**(      &
     1644                                 n_vg(j,i) / (n_vg(j,i) - 1.0_wp ) ) - 1.0_wp  &
     1645                             )**( 1.0_wp / n_vg(j,i) ) / alpha_vg(j,i)
     1646
     1647
     1648                      gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +         &
     1649                                      ( alpha_vg(j,i) * h_vg )**n_vg(j,i))**(  &
     1650                                      1.0_wp - 1.0_wp / n_vg(j,i) ) - (        &
     1651                                      alpha_vg(j,i) * h_vg )**( n_vg(j,i)      &
     1652                                      - 1.0_wp) )**2 )                         &
     1653                                      / ( ( 1.0_wp + ( alpha_vg(j,i) * h_vg    &
     1654                                      )**n_vg(j,i) )**( ( 1.0_wp  - 1.0_wp     &
     1655                                      / n_vg(j,i) ) *( l_vg(j,i) + 2.0_wp) ) )
     1656
     1657!
     1658!--                Parametrization of Clapp & Hornberger
     1659                   ELSE
     1660                      gamma_temp(k) = gamma_w_sat(j,i) * ( m_soil(k,j,i)       &
     1661                                      / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
     1662                   ENDIF
     1663
     1664                ENDDO
     1665
     1666!
     1667!--             Prognostic equation for soil moisture content. Only performed,
     1668!--             when humidity is enabled in the atmosphere and the surface type
     1669!--             is not pavement (implies dry soil below).
     1670                IF ( humidity  .AND.  .NOT.  pave_surface(j,i) )  THEN
     1671!
     1672!--                Calculate soil diffusivity (lambda_w) at the _stag level
     1673!--                using linear interpolation. To do: replace this with
     1674!--                ECMWF-IFS Eq. 8.81
     1675                   DO  k = nzb_soil, nzt_soil-1
     1676                     
     1677                      lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
     1678                                        * 0.5_wp
     1679                      gamma_w(k,j,i)  = ( gamma_temp(k+1) + gamma_temp(k) )    &
     1680                                        * 0.5_wp
     1681
     1682                   ENDDO
     1683
     1684!
     1685!
     1686!--                In case of a closed bottom (= water content is conserved), set
     1687!--                hydraulic conductivity to zero to that no water will be lost
     1688!--                in the bottom layer.
     1689                   IF ( conserve_water_content )  THEN
     1690                      gamma_w(nzt_soil,j,i) = 0.0_wp
     1691                   ELSE
     1692                      gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil)
     1693                   ENDIF     
     1694
     1695!--                The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v))
     1696!--                ensures the mass conservation for water. The transpiration of
     1697!--                plants equals the cumulative withdrawals by the roots in the
     1698!--                soil. The scheme takes into account the availability of water
     1699!--                in the soil layers as well as the root fraction in the
     1700!--                respective layer. Layer with moisture below wilting point will
     1701!--                not contribute, which reflects the preference of plants to
     1702!--                take water from moister layers.
     1703
     1704!
     1705!--                Calculate the root extraction (ECMWF 7.69, the sum of root_extr
     1706!--                = 1). The energy balance solver guarantees a positive
     1707!--                transpiration, so that there is no need for an additional check.
     1708                   DO  k = nzb_soil, nzt_soil
     1709                       IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
     1710                          m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
     1711                       ENDIF
     1712                   ENDDO 
     1713
     1714                   IF ( m_total > 0.0_wp )  THEN
     1715                      DO  k = nzb_soil, nzt_soil
     1716                         IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
     1717                            root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total
     1718                         ELSE
     1719                            root_extr(k) = 0.0_wp
     1720                         ENDIF
     1721                      ENDDO
     1722                   ENDIF
     1723
     1724!
     1725!--                Prognostic equation for soil water content m_soil.
     1726                   tend(:) = 0.0_wp
     1727
     1728                   tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (               &
     1729                            m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
     1730                            * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( &
     1731                               root_extr(nzb_soil) * qsws_veg_eb(j,i)          &
     1732                               + qsws_soil_eb(j,i) ) * drho_l_lv )             &
     1733                               * ddz_soil_stag(nzb_soil)
     1734
     1735                   DO  k = nzb_soil+1, nzt_soil-1
     1736                      tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)          &
     1737                                - m_soil(k,j,i) ) * ddz_soil(k+1)              &
     1738                                - gamma_w(k,j,i)                               &
     1739                                - lambda_w(k-1,j,i) * (m_soil(k,j,i) -         &
     1740                                m_soil(k-1,j,i)) * ddz_soil(k)                 &
     1741                                + gamma_w(k-1,j,i) - (root_extr(k)             &
     1742                                * qsws_veg_eb(j,i) * drho_l_lv)                &
     1743                                ) * ddz_soil_stag(k)
     1744
     1745                   ENDDO
     1746                   tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                  &
     1747                                           - lambda_w(nzt_soil-1,j,i)          &
     1748                                           * (m_soil(nzt_soil,j,i)             &
     1749                                           - m_soil(nzt_soil-1,j,i))           &
     1750                                           * ddz_soil(nzt_soil)                &
     1751                                           + gamma_w(nzt_soil-1,j,i) - (       &
     1752                                             root_extr(nzt_soil)               &
     1753                                           * qsws_veg_eb(j,i) * drho_l_lv  )   &
     1754                                     ) * ddz_soil_stag(nzt_soil)             
     1755
     1756                   m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
     1757                                                   + dt_3d * ( tsc(2) * tend(:)   &
     1758                                                   + tsc(3) * tm_soil_m(:,j,i) )   
     1759   
     1760!
     1761!--                Account for dry soils (find a better solution here!)
     1762                   DO  k = nzb_soil, nzt_soil
     1763                      IF ( m_soil_p(k,j,i) < 0.0_wp )  m_soil_p(k,j,i) = 0.0_wp
     1764                   ENDDO
     1765
     1766!
     1767!--                Calculate m_soil tendencies for the next Runge-Kutta step
     1768                   IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1769                      IF ( intermediate_timestep_count == 1 )  THEN
     1770                         DO  k = nzb_soil, nzt_soil
     1771                            tm_soil_m(k,j,i) = tend(k)
     1772                         ENDDO
     1773                      ELSEIF ( intermediate_timestep_count <                   &
     1774                               intermediate_timestep_count_max )  THEN
     1775                         DO  k = nzb_soil, nzt_soil
     1776                            tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp&
     1777                                     * tm_soil_m(k,j,i)
     1778                         ENDDO
     1779                      ENDIF
     1780                   ENDIF
     1781
     1782                ENDIF
     1783
    16401784             ENDIF
    16411785
    16421786          ENDDO
    16431787       ENDDO
    1644 
    1645 !
    1646 !--    Calculate surface specific humidity
    1647        IF ( humidity )  THEN
    1648           CALL calc_q_surface
    1649        ENDIF
    1650 
    16511788
    16521789    END SUBROUTINE lsm_soil_model
     
    16561793! Description:
    16571794! ------------
    1658 !> Calculation of specific humidity of the surface layer (surface)
     1795!> Calculation of roughness length for open water (lakes, ocean). The
     1796!> parameterization follows Charnock (1955). Two different implementations
     1797!> are available: as in ECMWF-IFS (Beljaars 1994) or as in FLake (Subin et al.
     1798!> 2012)
     1799!------------------------------------------------------------------------------!
     1800    SUBROUTINE calc_z0_water_surface
     1801
     1802       USE control_parameters,                                                 &
     1803           ONLY: g, kappa, molecular_viscosity
     1804
     1805       IMPLICIT NONE
     1806
     1807       INTEGER :: i  !< running index
     1808       INTEGER :: j  !< running index
     1809
     1810       REAL(wp), PARAMETER :: alpha_ch  = 0.018_wp !< Charnock constant (0.01-0.11). Use 0.01 for FLake and 0.018 for ECMWF
     1811!       REAL(wp), PARAMETER :: pr_number = 0.71_wp !< molecular Prandtl number in the Charnock parameterization (differs from prandtl_number)
     1812!       REAL(wp), PARAMETER :: sc_number = 0.66_wp !< molecular Schmidt number in the Charnock parameterization
     1813!       REAL(wp) :: re_0 !< near-surface roughness Reynolds number
     1814
     1815
     1816       DO  i = nxlg, nxrg   
     1817          DO  j = nysg, nyng
     1818             IF ( water_surface(j,i) )  THEN
     1819
     1820!
     1821!--             Disabled: FLake parameterization. Ideally, the Charnock
     1822!--             coefficient should depend on the water depth and the fetch
     1823!--             length
     1824!                re_0 = z0(j,i) * us(j,i) / molecular_viscosity
     1825!       
     1826!                z0(j,i) = MAX( 0.1_wp * molecular_viscosity / us(j,i),            &
     1827!                              alpha_ch * us(j,i) / g )
     1828!
     1829!                z0h(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 3.2_wp ) )
     1830!                z0q(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 4.2_wp ) )
     1831
     1832!
     1833!--              Set minimum roughness length for u* > 0.2
     1834!                IF ( us(j,i) > 0.2_wp )  THEN
     1835!                   z0h(j,i) = MAX( 1.0E-5_wp, z0h(j,i) )
     1836!                   z0q(j,i) = MAX( 1.0E-5_wp, z0q(j,i) )
     1837!                ENDIF
     1838
     1839!
     1840!--             ECMWF IFS model parameterization after Beljaars (1994). At low
     1841!--             wind speed, the sea surface becomes aerodynamically smooth and
     1842!--             the roughness scales with the viscosity. At high wind speed, the
     1843!--             Charnock relation is used.
     1844                z0(j,i) =   ( 0.11_wp * molecular_viscosity / us(j,i) )        &
     1845                          + ( alpha_ch * us(j,i)**2 / g )
     1846
     1847                z0h(j,i) = 0.40_wp * molecular_viscosity / us(j,i)
     1848                z0q(j,i) = 0.62_wp * molecular_viscosity / us(j,i)
     1849
     1850             ENDIF
     1851          ENDDO
     1852       ENDDO
     1853
     1854    END SUBROUTINE calc_z0_water_surface
     1855
     1856
     1857!------------------------------------------------------------------------------!
     1858! Description:
     1859! ------------
     1860!> Calculation of specific humidity of the skin layer (surface). It is assumend
     1861!> that the skin is always saturated.
    16591862!------------------------------------------------------------------------------!
    16601863    SUBROUTINE calc_q_surface
     
    16651868       INTEGER :: j              !< running index
    16661869       INTEGER :: k              !< running index
     1870
    16671871       REAL(wp) :: resistance    !< aerodynamic and soil resistance term
    16681872
     
    17011905    END SUBROUTINE calc_q_surface
    17021906
     1907
    17031908!------------------------------------------------------------------------------!
    17041909! Description:
  • palm/trunk/SOURCE/modules.f90

    r1787 r1788  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Added roughness length for moisture (z0q)
    2222!
    2323! Former revisions:
     
    384384          saswsb, saswst, shf, tnudge, td_lsa_lpt, td_lsa_q, td_sub_lpt,       &
    385385          td_sub_q, total_2d_a, total_2d_o, ts, tswst, ug_vert, unudge, us,    &
    386           usws, uswst, vnudge, vg_vert, vsws, vswst, wnudge, wsubs_vert, z0, z0h
     386          usws, uswst, vnudge, vg_vert, vsws, vswst, wnudge, wsubs_vert,       &
     387          z0,                                                                  & !< roughness length for momentum
     388          z0h,                                                                 & !< roughness length for heat
     389          z0q                                                                    !< roughness length for moisture
    387390
    388391    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
     
    432435    USE kinds
    433436
    434     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lwp_av,                & !> Avg. liquid water path
    435                                               precipitation_rate_av, & !> Avg. of precipitation rate
    436                                               ol_av,                 & !> Avg. of Obukhov length
    437                                               qsws_av,               & !> Avg. of surface moisture flux
    438                                               shf_av,                & !> Avg. of surface heat flux
    439                                               ts_av,                 & !> Avg. of characteristic temperature scale
    440                                               us_av,                 & !> Avg. of friction velocity
    441                                               z0_av,                 & !> Avg. of roughness length for momentum
    442                                               z0h_av                   !> Avg. of roughness length for heat and moisture
     437    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lwp_av,                & !< Avg. liquid water path
     438                                              precipitation_rate_av, & !< Avg. of precipitation rate
     439                                              ol_av,                 & !< Avg. of Obukhov length
     440                                              qsws_av,               & !< Avg. of surface moisture flux
     441                                              shf_av,                & !< Avg. of surface heat flux
     442                                              ts_av,                 & !< Avg. of characteristic temperature scale
     443                                              us_av,                 & !< Avg. of friction velocity
     444                                              z0_av,                 & !< Avg. of roughness length for momentum
     445                                              z0h_av,                & !< Avg. of roughness length for heat
     446                                              z0q_av                   !< Avg. of roughness length for moisture
    443447
    444448    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
  • palm/trunk/SOURCE/package_parin.f90

    r1787 r1788  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Parameter dewfall removed.
    2222!
    2323! Former revisions:
     
    129129    USE land_surface_model_mod,                                                &
    130130        ONLY: alpha_vangenuchten, c_surface, canopy_resistance_coefficient,    &
    131               conserve_water_content, dewfall, field_capacity,                 &
     131              conserve_water_content, field_capacity,                          &
    132132              f_shortwave_incoming, hydraulic_conductivity,                    &
    133133              lambda_surface_stable, lambda_surface_unstable, leaf_area_index, &
    134134              land_surface, l_vangenuchten, min_canopy_resistance,             &
    135               min_soil_resistance, n_vangenuchten, residual_moisture,          &
     135              min_soil_resistance, n_vangenuchten, pave_heat_capacity,         &
     136              pave_depth, pave_heat_conductivity, residual_moisture,           &
    136137              root_fraction, saturation_moisture, wilting_point,               &
    137138              skip_time_do_lsm, soil_moisture, soil_temperature, soil_type,    &
    138               vegetation_coverage, veg_type, zs, z0_eb, z0h_eb
     139              vegetation_coverage, veg_type, zs, z0_eb, z0h_eb, z0q_eb
    139140
    140141    USE kinds
     
    205206    NAMELIST /lsm_par/            alpha_vangenuchten, c_surface,               &
    206207                                  canopy_resistance_coefficient,               &
    207                                   conserve_water_content, dewfall,             &
     208                                  conserve_water_content,                      &
    208209                                  f_shortwave_incoming, field_capacity,        &
    209210                                  hydraulic_conductivity,                      &
     
    212213                                  l_vangenuchten, min_canopy_resistance,       &
    213214                                  min_soil_resistance, n_vangenuchten,         &
     215                                  pave_depth, pave_heat_capacity,              &
     216                                  pave_heat_conductivity,                      &
    214217                                  residual_moisture, root_fraction,            &
    215218                                  saturation_moisture, skip_time_do_lsm,       &
    216219                                  soil_moisture, soil_temperature, soil_type,  &
    217220                                  vegetation_coverage, veg_type, wilting_point,&
    218                                   zs, z0_eb, z0h_eb
     221                                  zs, z0_eb, z0h_eb, z0q_eb
    219222
    220223
  • palm/trunk/SOURCE/radiation_model.f90

    r1784 r1788  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Added new albedo class for pavements / roads.
    2222!
    2323! Former revisions:
     
    112112
    113113#if defined ( __rrtmg )
     114
     115!     USE netcdf_interface,                                                      &
     116!         ONLY:  nc_stat, netcdf_handle_error
     117
    114118    USE parrrsw,                                                               &
    115119        ONLY:  naerec, nbndsw
     
    139143!
    140144!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
    141     CHARACTER(37), DIMENSION(0:16), PARAMETER :: albedo_type_name = (/      &
     145    CHARACTER(37), DIMENSION(0:17), PARAMETER :: albedo_type_name = (/      &
    142146                                   'user defined                         ', & !  0
    143147                                   'ocean                                ', & !  1
     
    156160                                   'land ice                             ', & ! 14
    157161                                   'sea ice                              ', & ! 15
    158                                    'snow                                 '  & ! 16
     162                                   'snow                                 ', & ! 16
     163                                   'pavement/roads                       '  & ! 17
    159164                                                         /)
    160165
     
    163168                    day_init     = 172,  & !< day of the year at model start (21/06)
    164169                    dots_rad     = 0       !< starting index for timeseries output
    165 
    166 
    167 
    168 
    169 
    170170
    171171    LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
     
    214214!-- Land surface albedos for solar zenith angle of 60° after Briegleb (1992)     
    215215!-- (shortwave, longwave, broadband):   sw,      lw,      bb,
    216     REAL(wp), DIMENSION(0:2,1:16), PARAMETER :: albedo_pars = RESHAPE( (/&
     216    REAL(wp), DIMENSION(0:2,1:17), PARAMETER :: albedo_pars = RESHAPE( (/&
    217217                                   0.06_wp, 0.06_wp, 0.06_wp,            & !  1
    218218                                   0.09_wp, 0.28_wp, 0.19_wp,            & !  2
     
    230230                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 14
    231231                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 15
    232                                    0.95_wp, 0.70_wp, 0.82_wp             & ! 16
    233                                  /), (/ 3, 16 /) )
     232                                   0.95_wp, 0.70_wp, 0.82_wp,            & ! 16
     233                                   0.08_wp, 0.08_wp, 0.08_wp             & ! 17
     234                                 /), (/ 3, 17 /) )
    234235
    235236    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
     
    275276                    rrtm_iaer = 0,     & !< aerosol option flag (0: no aerosol layers, for lw only: 6 (requires setting of rrtm_sw_ecaer), 10: one or more aerosol layers (not implemented)
    276277                    rrtm_idrv = 1        !< longwave upward flux calculation option (0,1)
     278
     279    INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
    277280
    278281    LOGICAL :: snd_exists = .FALSE.      !< flag parameter to check whether a user-defined input files exists
     
    10301033                 rrtm_aldir(0,:,:) = aldif
    10311034                 rrtm_asdir(0,:,:) = asdif
     1035
     1036!
     1037!--        Asphalt
     1038           ELSEIF ( albedo_type == 17 )  THEN
     1039                 rrtm_aldir(0,:,:) = aldif
     1040                 rrtm_asdir(0,:,:) = asdif
    10321041!
    10331042!--        Land surfaces
  • palm/trunk/SOURCE/read_3d_binary.f90

    r1758 r1788  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Added z0q and z0q_av
    2222!
    2323! Former revisions:
     
    9797               rif_wall, shf, ts, tswst, u, u_m_l, u_m_n, u_m_r, u_m_s, us,    &
    9898               usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s, vpt, vsws, vswst,   &
    99                w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h
     99               w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h, z0q
    100100
    101101    USE averaging
     
    354354!--    First compare the version numbers
    355355       READ ( 13 )  version_on_file
    356        binary_version = '4.2'
     356       binary_version = '4.3'
    357357       IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
    358358          WRITE( message_string, * ) 'version mismatch concerning data ',      &
     
    13581358                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    13591359                   z0h_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     1360                       tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     1361
     1362                CASE ( 'z0q' )
     1363                   IF ( k == 1 )  READ ( 13 )  tmp_2d
     1364                   z0q(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     1365                     tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     1366
     1367                CASE ( 'z0q_av' )
     1368                   IF ( .NOT. ALLOCATED( z0q_av ) )  THEN
     1369                      ALLOCATE( z0q_av(nysg:nyng,nxlg:nxrg) )
     1370                   ENDIF
     1371                   IF ( k == 1 )  READ ( 13 )  tmp_2d
     1372                   z0q_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
    13601373                       tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    13611374
  • palm/trunk/SOURCE/sum_up_3d_data.f90

    r1694 r1788  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Added z0q and z0q_av
    2222!
    2323! Former revisions:
     
    9292    USE arrays_3d,                                                             &
    9393        ONLY:  dzw, e, nr, ol, p, pt, q, qc, ql, ql_c, ql_v, qr, qsws, rho, sa,&
    94                shf, ts, u, us, v, vpt, w, z0, z0h
     94               shf, ts, u, us, v, vpt, w, z0, z0h, z0q
    9595
    9696    USE averaging,                                                             &
     
    9898               precipitation_rate_av, pt_av, q_av, qc_av, ql_av, ql_c_av,      &
    9999               ql_v_av, ql_vp_av, qr_av, qsws_av, qv_av, rho_av, s_av, sa_av,  &
    100                shf_av, ts_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, z0h_av
     100               shf_av, ts_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, z0h_av,  &
     101               z0q_av
    101102
    102103    USE cloud_parameters,                                                      &
     
    494495                ENDIF
    495496                z0h_av = 0.0_wp
     497
     498             CASE ( 'z0q*' )
     499                IF ( .NOT. ALLOCATED( z0q_av ) )  THEN
     500                   ALLOCATE( z0q_av(nysg:nyng,nxlg:nxrg) )
     501                ENDIF
     502                z0q_av = 0.0_wp
    496503
    497504             CASE DEFAULT
     
    10061013             ENDDO
    10071014
     1015          CASE ( 'z0q*' )
     1016             DO  i = nxlg, nxrg
     1017                DO  j = nysg, nyng
     1018                   z0q_av(j,i) = z0q_av(j,i) + z0q(j,i)
     1019                ENDDO
     1020             ENDDO
     1021
    10081022          CASE DEFAULT
    10091023!
  • palm/trunk/SOURCE/surface_layer_fluxes.f90

    r1758 r1788  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Added parameter z0q which replaces z0h in the similarity functions for
     22! humidity.
     23! Syntax layout improved.
    2224!
    2325! Former revisions:
     
    117119    USE arrays_3d,                                                             &
    118120        ONLY:  e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, &
    119                shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h
     121               shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h, z0q
    120122
    121123    USE cloud_parameters,                                                      &
     
    234236!--    Use either Newton iteration or a lookup table for the bulk Richardson
    235237!--    number to calculate the Obukhov length
    236        ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' )  THEN
     238       ELSEIF ( most_method == 'newton'  .OR. most_method == 'lookup' )  THEN
    237239
    238240          CALL calc_uv_total
     
    302304!--       Check for roughness heterogeneity. In that case terminate run and
    303305!--       inform user
    304           IF ( MINVAL( z0h ) /= MAXVAL( z0h ) .OR.                             &
     306          IF ( MINVAL( z0h ) /= MAXVAL( z0h )  .OR.                            &
    305307               MINVAL( z0  ) /= MAXVAL( z0  ) )  THEN
    306308             terminate_run_l = .TRUE.
     
    347349!
    348350!--       Calculate stretching factor for the free convection range
    349           DO WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
     351          DO  WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
    350352             regr_old = regr
    351353             regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr )   &
     
    533535!--             Ensure that the bulk Richardson number and the Obukhov
    534536!--             lengtH have the same sign
    535                 IF ( rib(j,i) * ol(j,i) < 0.0_wp .OR.                        &
     537                IF ( rib(j,i) * ol(j,i) < 0.0_wp  .OR.                         &
    536538                     ABS( ol(j,i) ) == ol_max )  THEN
    537539                   IF ( rib(j,i) > 0.0_wp ) ol(j,i) =  0.01_wp
     
    638640
    639641          !$OMP PARALLEL DO PRIVATE( k, z_mo )
    640           !# WARNING: does not work on GPU so far because of DO WHILE construct
     642          !# WARNING: does not work on GPU so far because of DO  WHILE construct
    641643          !!!!!!$acc kernels loop
    642644          DO  i = nxl, nxr
     
    656658                l = l_bnd
    657659                IF ( rib_tab(l) - rib(j,i) > 0.0_wp )  THEN
    658                    DO WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp .AND. l > 0 )
     660                   DO  WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp  .AND. l > 0 )
    659661                      l = l-1
    660662                   ENDDO
    661663                ELSE
    662                    DO WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp                &
    663                               .AND. l < num_steps-1 )
     664                   DO  WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp                  &
     665                              .AND.  l < num_steps-1 )
    664666                      l = l+1
    665667                   ENDDO
     
    669671!
    670672!--             Linear interpolation to find the correct value of z/L
    671                 ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) )       &
    672                             / (  rib_tab(l) - rib_tab(l-1) )                &
     673                ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) )          &
     674                            / (  rib_tab(l) - rib_tab(l-1) )                   &
    673675                            * ( rib(j,i) - rib_tab(l-1) ) )
    674676
     
    797799!
    798800!--       For a given surface temperature:
    799           IF ( large_scale_forcing .AND. lsf_surf )  THEN
     801          IF ( large_scale_forcing  .AND. lsf_surf )  THEN
    800802             !$OMP PARALLEL DO
    801803             !$acc kernels loop private( j, k )
     
    837839          IF ( constant_waterflux )  THEN
    838840!
    839 !--          For a given water flux in the Prandtl layer:
     841!--          For a given water flux in the surface layer
    840842             !$OMP PARALLEL DO
    841843             !$acc kernels loop private( j )
     
    847849
    848850          ELSE
    849              coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND.      &
     851             coupled_run = ( coupling_mode == 'atmosphere_to_ocean'  .AND.     &
    850852                             run_coupled )
    851853
    852              IF ( large_scale_forcing .AND. lsf_surf )  THEN
     854             IF ( large_scale_forcing  .AND. lsf_surf )  THEN
    853855                !$OMP PARALLEL DO
    854856                !$acc kernels loop private( j, k )
     
    862864
    863865             !$OMP PARALLEL DO PRIVATE( e_s, k, z_mo )
    864              !$acc kernels loop independent present( nzb_s_inner, ol, pt, q, qs, qv1, zu, zw, z0h ) private( e_s, j, k, z_mo )
     866             !$acc kernels loop independent present( nzb_s_inner, ol, pt, q, qs, qv1, zu, zw, z0q ) private( e_s, j, k, z_mo )
    865867             DO  i = nxlg, nxrg
    866868                !$acc loop independent
     
    882884                   IF ( cloud_physics )  THEN
    883885                      qs(j,i) = kappa * ( qv1(j,i) - q(k,j,i) )                &
    884                                         / ( LOG( z_mo / z0h(j,i) )             &
     886                                        / ( LOG( z_mo / z0q(j,i) )             &
    885887                                            - psi_h( z_mo / ol(j,i) )          &
    886                                             + psi_h( z0h(j,i) / ol(j,i) ) )
     888                                            + psi_h( z0q(j,i) / ol(j,i) ) )
    887889
    888890                   ELSE
    889891                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) )              &
    890                                         / ( LOG( z_mo / z0h(j,i) )             &
     892                                        / ( LOG( z_mo / z0q(j,i) )             &
    891893                                            - psi_h( z_mo / ol(j,i) )          &
    892                                             + psi_h( z0h(j,i) / ol(j,i) ) )
     894                                            + psi_h( z0q(j,i) / ol(j,i) ) )
    893895                   ENDIF
    894896
     
    901903!
    902904!--    If required compute qr* and nr*
    903        IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation )  THEN
     905       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )   &
     906       THEN
    904907
    905908          !$OMP PARALLEL DO PRIVATE( k, z_mo )
    906           !$acc kernels loop independent present( nr, nrs, nzb_s_inner, ol, qr, qrs, zu, zw, z0h ) private( j, k, z_mo )
     909          !$acc kernels loop independent present( nr, nrs, nzb_s_inner, ol, qr, qrs, zu, zw, z0q ) private( j, k, z_mo )
    907910          DO  i = nxlg, nxrg
    908911             !$acc loop independent
     
    913916
    914917                qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) )              &
    915                                  / ( LOG( z_mo / z0h(j,i) )                 &
     918                                 / ( LOG( z_mo / z0q(j,i) )                 &
    916919                                     - psi_h( z_mo / ol(j,i) )              &
    917                                   &nb