Ignore:
Timestamp:
Sep 14, 2020 7:55:28 AM (7 weeks ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

File:
1 edited

Legend:

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

    r4586 r4677  
    11!> @file model_1d_mod.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
     
    2524! -----------------
    2625! $Id$
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4586 2020-07-01 16:16:43Z gronemeier
    2729! renamed Richardson flux number into gradient Richardson number
    2830!
     
    5355!>        the 1D model uses different turbulence closure approaches at least if
    5456!>        the 3D model is set to LES-mode.
    55 !------------------------------------------------------------------------------!
     57!--------------------------------------------------------------------------------------------------!
    5658 MODULE model_1d_mod
    5759
    58     USE arrays_3d,                                                             &
    59         ONLY:  dd2zu, ddzu, ddzw, dzu, dzw, pt_init, q_init, ug, u_init,       &
    60                vg, v_init, zu
    61 
    62     USE basic_constants_and_equations_mod,                                     &
     60    USE arrays_3d,                                                                                 &
     61        ONLY:  dd2zu, ddzu, ddzw, dzu, dzw, pt_init, q_init, ug, u_init, vg, v_init, zu
     62
     63    USE basic_constants_and_equations_mod,                                                         &
    6364        ONLY:  g, kappa, pi
    6465
    65     USE control_parameters,                                                    &
    66         ONLY:  constant_diffusion, constant_flux_layer, dissipation_1d, f,     &
    67                humidity, ibc_e_b, intermediate_timestep_count,                 &
    68                intermediate_timestep_count_max, km_constant,                   &
    69                message_string, mixing_length_1d, prandtl_number,               &
    70                roughness_length, run_description_header, simulated_time_chr,   &
    71                timestep_scheme, tsc, z0h_factor
    72 
    73     USE indices,                                                               &
     66    USE control_parameters,                                                                        &
     67        ONLY:  constant_diffusion, constant_flux_layer, dissipation_1d, f, humidity, ibc_e_b,      &
     68               intermediate_timestep_count, intermediate_timestep_count_max, km_constant,          &
     69               message_string, mixing_length_1d, prandtl_number, roughness_length,                 &
     70               run_description_header, simulated_time_chr, timestep_scheme, tsc, z0h_factor
     71
     72    USE indices,                                                                                   &
    7473        ONLY:  nzb, nzb_diff, nzt
    7574
    7675    USE kinds
    7776
    78     USE pegrid,                                                                &
     77    USE pegrid,                                                                                    &
    7978        ONLY:  myid
    8079
     
    176175!
    177176!-- Public variables
    178     PUBLIC  damp_level_1d, damp_level_ind_1d, diss1d, dt_pr_1d,                &
    179             dt_run_control_1d, e1d, end_time_1d, kh1d, km1d, l1d, ri1d, u1d,   &
    180             us1d, usws1d, v1d, vsws1d
     177    PUBLIC  damp_level_1d, damp_level_ind_1d, diss1d, dt_pr_1d, dt_run_control_1d, e1d,            &
     178            end_time_1d, kh1d, km1d, l1d, ri1d, u1d, us1d, usws1d, v1d, vsws1d
    181179
    182180
     
    185183 SUBROUTINE init_1d_model
    186184
    187     USE grid_variables,                                                        &
     185    USE grid_variables,                                                                            &
    188186        ONLY:  dx, dy
    189187
     
    196194!
    197195!-- Allocate required 1D-arrays
    198     ALLOCATE( diss1d(nzb:nzt+1), diss1d_p(nzb:nzt+1),                          &
    199               e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1),               &
    200               km1d(nzb:nzt+1), l1d(nzb:nzt+1), l1d_init(nzb:nzt+1),            &
    201               l1d_diss(nzb:nzt+1), ri1d(nzb:nzt+1), te_diss(nzb:nzt+1),        &
    202               te_dissm(nzb:nzt+1), te_e(nzb:nzt+1),                            &
    203               te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1),             &
    204               te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1),               &
     196    ALLOCATE( diss1d(nzb:nzt+1), diss1d_p(nzb:nzt+1),                                              &
     197              e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1),                                   &
     198              km1d(nzb:nzt+1), l1d(nzb:nzt+1), l1d_init(nzb:nzt+1),                                &
     199              l1d_diss(nzb:nzt+1), ri1d(nzb:nzt+1), te_diss(nzb:nzt+1),                            &
     200              te_dissm(nzb:nzt+1), te_e(nzb:nzt+1),                                                &
     201              te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1),                                 &
     202              te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1),                                   &
    205203              u1d_p(nzb:nzt+1),  v1d(nzb:nzt+1), v1d_p(nzb:nzt+1) )
    206204
     
    223221!--       Blackadar mixing length
    224222          IF ( f /= 0.0_wp )  THEN
    225              lambda = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /        &
    226                                ABS( f ) + 1E-10_wp
     223             lambda = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) / ABS( f ) + 1E-10_wp
    227224          ELSE
    228225             lambda = 30.0_wp
     
    258255
    259256!
    260 !-- Set initial horizontal velocities at the lowest grid levels to a very small
    261 !-- value in order to avoid too small time steps caused by the diffusion limit
    262 !-- in the initial phase of a run (at k=1, dz/2 occurs in the limiting formula!)
     257!-- Set initial horizontal velocities at the lowest grid levels to a very small value in order to
     258!-- avoid too small time steps caused by the diffusion limit in the initial phase of a run (at k=1,
     259!-- dz/2 occurs in the limiting formula!)
    263260    u1d(0:1)   = 0.1_wp
    264261    u1d_p(0:1) = 0.1_wp
     
    310307
    311308
    312 !------------------------------------------------------------------------------!
     309!--------------------------------------------------------------------------------------------------!
    313310! Description:
    314311! ------------
    315312!> Runge-Kutta time differencing scheme for the 1D-model.
    316 !------------------------------------------------------------------------------!
     313!--------------------------------------------------------------------------------------------------!
    317314
    318315 SUBROUTINE time_integration_1d
     
    335332
    336333!
    337 !-- Determine the time step at the start of a 1D-simulation and
    338 !-- determine and printout quantities used for run control
     334!-- Determine the time step at the start of a 1D-simulation and determine and printout quantities
     335!-- used for run control
    339336    dt_1d = 0.01_wp
    340337    CALL run_control_1d
     
    345342
    346343!
    347 !--    Depending on the timestep scheme, carry out one or more intermediate
    348 !--    timesteps
     344!--    Depending on the timestep scheme, carry out one or more intermediate timesteps
    349345
    350346       intermediate_timestep_count = 0
    351        DO  WHILE ( intermediate_timestep_count < &
    352                    intermediate_timestep_count_max )
     347       DO  WHILE ( intermediate_timestep_count < intermediate_timestep_count_max )
    353348
    354349          intermediate_timestep_count = intermediate_timestep_count + 1
     
    357352
    358353!
    359 !--       Compute all tendency terms. If a constant-flux layer is simulated,
    360 !--       k starts at nzb+2.
     354!--       Compute all tendency terms. If a constant-flux layer is simulated, k starts at nzb+2.
    361355          DO  k = nzb_diff, nzt
    362356
     
    365359!
    366360!--          u-component
    367              te_u(k) =  f * ( v1d(k) - vg(k) ) + ( &
    368                               kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) &
    369                             - kmzm * ( u1d(k) - u1d(k-1) ) * ddzu(k)   &
     361             te_u(k) =  f * ( v1d(k) - vg(k) ) + (                                                 &
     362                              kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1)                             &
     363                            - kmzm * ( u1d(k) - u1d(k-1) ) * ddzu(k)                               &
    370364                                                 ) * ddzw(k)
    371365!
    372366!--          v-component
    373              te_v(k) = -f * ( u1d(k) - ug(k) ) + (                     &
    374                               kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) &
    375                             - kmzm * ( v1d(k) - v1d(k-1) ) * ddzu(k)   &
     367             te_v(k) = -f * ( u1d(k) - ug(k) ) + (                                                 &
     368                              kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1)                             &
     369                            - kmzm * ( v1d(k) - v1d(k-1) ) * ddzu(k)                               &
    376370                                                 ) * ddzw(k)
    377371          ENDDO
     
    387381                ELSE
    388382                   pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
    389                    flux = ( ( pt_init(k+1) - pt_init(k-1) ) +                  &
    390                             0.61_wp * ( pt_init(k+1) * q_init(k+1) -           &
    391                                         pt_init(k-1) * q_init(k-1)   )         &
     383                   flux = ( ( pt_init(k+1) - pt_init(k-1) ) +                                      &
     384                            0.61_wp * ( pt_init(k+1) * q_init(k+1) -                               &
     385                                        pt_init(k-1) * q_init(k-1)   )                             &
    392386                          ) * dd2zu(k)
    393387                ENDIF
    394388
    395389!
    396 !--             Calculate dissipation rate if no prognostic equation is used for
    397 !--             dissipation rate
     390!--             Calculate dissipation rate if no prognostic equation is used for dissipation rate.
    398391                IF ( dissipation_1d == 'detering' )  THEN
    399392                   diss1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    400393                ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
    401                    diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) &
    402                                ) * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
     394                   diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) )                   &
     395                               * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    403396                ENDIF
    404397!
    405398!--             TKE
    406                 te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
    407                                     + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
    408                                     )                                          &
    409                                     - g / pt_0 * kh1d(k) * flux                &
    410                                     +            (                             &
    411                                      kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)  &
    412                                    - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
    413                                                  ) * ddzw(k) / sig_e           &
     399                te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2                    &
     400                                    + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2                    &
     401                                    )                                                              &
     402                                    - g / pt_0 * kh1d(k) * flux                                    &
     403                                    +            (                                                 &
     404                                     kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)                      &
     405                                   - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)                        &
     406                                                 ) * ddzw(k) / sig_e                               &
    414407                                   - diss1d(k)
    415408
     
    420413                      alpha_buoyancy = 1.0_wp - l1d(k) / lambda
    421414                   ELSE
    422                       alpha_buoyancy = 1.0_wp - ( 1.0_wp + ( c_2 - 1.0_wp )    &
    423                                                          / ( c_2 - c_1    ) )  &
     415                      alpha_buoyancy = 1.0_wp - ( 1.0_wp + ( c_2 - 1.0_wp )                        &
     416                                                         / ( c_2 - c_1    ) )                      &
    424417                                              * l1d(k) / lambda
    425418                   ENDIF
    426419                   c_3 = ( c_1 - c_2 ) * alpha_buoyancy + 1.0_wp
    427                    te_diss(k) = ( km1d(k) *                                    &
    428                                   ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2  &
    429                                   + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2  &
    430                                   ) * ( c_1 + (c_2 - c_1) * l1d(k) / lambda )  &
    431                                   - g / pt_0 * kh1d(k) * flux * c_3            &
    432                                   - c_2 * diss1d(k)                            &
    433                                 ) * diss1d(k) / ( e1d(k) + 1.0E-20_wp )        &
    434                                 + (   kmzp * ( diss1d(k+1) - diss1d(k) )       &
    435                                            * ddzu(k+1)                         &
    436                                     - kmzm * ( diss1d(k) - diss1d(k-1) )       &
    437                                            * ddzu(k)                           &
     420                   te_diss(k) = ( km1d(k) *                                                        &
     421                                  ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2                      &
     422                                  + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2                      &
     423                                  ) * ( c_1 + (c_2 - c_1) * l1d(k) / lambda )                      &
     424                                  - g / pt_0 * kh1d(k) * flux * c_3                                &
     425                                  - c_2 * diss1d(k)                                                &
     426                                ) * diss1d(k) / ( e1d(k) + 1.0E-20_wp )                            &
     427                                + (   kmzp * ( diss1d(k+1) - diss1d(k) )                           &
     428                                           * ddzu(k+1)                                             &
     429                                    - kmzm * ( diss1d(k) - diss1d(k-1) )                           &
     430                                           * ddzu(k)                                               &
    438431                                  ) * ddzw(k) / sig_diss
    439432
     
    445438!
    446439!--       Tendency terms at the top of the constant-flux layer.
    447 !--       Finite differences of the momentum fluxes are computed using half the
    448 !--       normal grid length (2.0*ddzw(k)) for the sake of enhanced accuracy
     440!--       Finite differences of the momentum fluxes are computed using half the normal grid length
     441!--       (2.0*ddzw(k)) for the sake of enhanced accuracy
    449442          IF ( constant_flux_layer )  THEN
    450443
     
    457450             ELSE
    458451                pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
    459                 flux = ( ( pt_init(k+1) - pt_init(k-1) ) +                     &
    460                          0.61_wp * ( pt_init(k+1) * q_init(k+1) -              &
    461                                      pt_init(k-1) * q_init(k-1)   )            &
     452                flux = ( ( pt_init(k+1) - pt_init(k-1) ) +                                         &
     453                         0.61_wp * ( pt_init(k+1) * q_init(k+1) -                                  &
     454                                     pt_init(k-1) * q_init(k-1)   )                                &
    462455                       ) * dd2zu(k)
    463456             ENDIF
    464457
    465458!
    466 !--          Calculate dissipation rate if no prognostic equation is used for
    467 !--          dissipation rate
     459!--          Calculate dissipation rate if no prognostic equation is used for dissipation rate.
    468460             IF ( dissipation_1d == 'detering' )  THEN
    469461                diss1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    470462             ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
    471                 diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) )  &
     463                diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) )                      &
    472464                            * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    473465             ENDIF
     
    475467!
    476468!--          u-component
    477              te_u(k) = f * ( v1d(k) - vg(k) ) + (                              &
    478                        kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) + usws1d       &
     469             te_u(k) = f * ( v1d(k) - vg(k) ) + (                                                  &
     470                       kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) + usws1d                           &
    479471                                                ) * 2.0_wp * ddzw(k)
    480472!
    481473!--          v-component
    482              te_v(k) = -f * ( u1d(k) - ug(k) ) + (                             &
    483                        kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) + vsws1d       &
     474             te_v(k) = -f * ( u1d(k) - ug(k) ) + (                                                 &
     475                       kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) + vsws1d                           &
    484476                                                 ) * 2.0_wp * ddzw(k)
    485477!
     
    490482                !>   while for u and v it is not?
    491483                !>   2018-04-23, gronemeier
    492                 te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
    493                                     + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
    494                                     )                                          &
    495                                     - g / pt_0 * kh1d(k) * flux                &
    496                                     +           (                              &
    497                                      kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)  &
    498                                    - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
    499                                                  ) * ddzw(k) / sig_e           &
     484                te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2                    &
     485                                    + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2                    &
     486                                    )                                                              &
     487                                    - g / pt_0 * kh1d(k) * flux                                    &
     488                                    +           (                                                  &
     489                                     kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)                      &
     490                                   - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)                        &
     491                                                 ) * ddzw(k) / sig_e                               &
    500492                                   - diss1d(k)
    501493             ENDIF
     
    507499          DO  k = nzb+1, nzt
    508500
    509              u1d_p(k) = u1d(k) + dt_1d * ( tsc(2) * te_u(k) + &
    510                                            tsc(3) * te_um(k) )
    511              v1d_p(k) = v1d(k) + dt_1d * ( tsc(2) * te_v(k) + &
    512                                            tsc(3) * te_vm(k) )
     501             u1d_p(k) = u1d(k) + dt_1d * ( tsc(2) * te_u(k) + tsc(3) * te_um(k) )
     502             v1d_p(k) = v1d(k) + dt_1d * ( tsc(2) * te_v(k) + tsc(3) * te_vm(k) )
    513503
    514504          ENDDO
     
    516506
    517507             DO  k = nzb+1, nzt
    518                 e1d_p(k) = e1d(k) + dt_1d * ( tsc(2) * te_e(k) + &
    519                                               tsc(3) * te_em(k) )
     508                e1d_p(k) = e1d(k) + dt_1d * ( tsc(2) * te_e(k) + tsc(3) * te_em(k) )
    520509             ENDDO
    521510
    522511!
    523 !--          Eliminate negative TKE values, which can result from the
    524 !--          integration due to numerical inaccuracies. In such cases the TKE
    525 !--          value is reduced to 10 percent of its old value.
     512!--          Eliminate negative TKE values, which can result from the integration due to numerical
     513!--          inaccuracies. In such cases the TKE value is reduced to 10 percent of its old value.
    526514             WHERE ( e1d_p < 0.0_wp )  e1d_p = 0.1_wp * e1d
    527515
    528516             IF ( dissipation_1d == 'prognostic' )  THEN
    529517                DO  k = nzb+1, nzt
    530                    diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + &
    531                                                        tsc(3) * te_dissm(k) )
     518                   diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + tsc(3) * te_dissm(k) )
    532519                ENDDO
    533520                WHERE ( diss1d_p < 0.0_wp )  diss1d_p = 0.1_wp * diss1d
     
    556543                ENDIF
    557544
    558              ELSEIF ( intermediate_timestep_count < &
    559                          intermediate_timestep_count_max )  THEN
     545             ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
    560546
    561547                DO  k = nzb+1, nzt
     
    570556                   IF ( dissipation_1d == 'prognostic' )  THEN
    571557                      DO k = nzb+1, nzt
    572                          te_dissm(k) = -9.5625_wp * te_diss(k)  &
    573                                      +  5.3125_wp * te_dissm(k)
     558                         te_dissm(k) = -9.5625_wp * te_diss(k) + 5.3125_wp * te_dissm(k)
    574559                      ENDDO
    575560                   ENDIF
     
    581566!
    582567!--       Boundary conditions for the prognostic variables.
    583 !--       At the top boundary (nzt+1) u, v, e, and diss keep their initial
    584 !--       values (ug(nzt+1), vg(nzt+1), 0, 0).
    585 !--       At the bottom boundary, Dirichlet condition is used for u and v (0)
    586 !--       and Neumann condition for e and diss (e(nzb)=e(nzb+1)).
     568!--       At the top boundary (nzt+1) u, v, e, and diss keep their initial values (ug(nzt+1),
     569!--       vg(nzt+1), 0, 0).
     570!--       At the bottom boundary, Dirichlet condition is used for u and v (0) and Neumann condition
     571!--       for e and diss (e(nzb)=e(nzb+1)).
    587572          u1d_p(nzb) = 0.0_wp
    588573          v1d_p(nzb) = 0.0_wp
     
    611596!
    612597!--                Stable stratification
    613                    ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /          &
    614                           ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) *  &
    615                                           ( zu(nzb+1) - z0h1d ) / zu(nzb+1)    &
     598                   ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /                              &
     599                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) *                      &
     600                                          ( zu(nzb+1) - z0h1d ) / zu(nzb+1)                        &
    616601                          )
    617602                ELSE
     
    619604!--                Unstable stratification
    620605                   a = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) )
    621                    b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) /                  &
    622                        zu(nzb+1) * z0h1d )
    623 
    624                    ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /          &
    625                           LOG( (a-1.0_wp) / (a+1.0_wp) *                       &
    626                                (b+1.0_wp) / (b-1.0_wp) )
     606                   b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / zu(nzb+1) * z0h1d )
     607
     608                   ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /                              &
     609                          LOG( (a-1.0_wp) / (a+1.0_wp) * (b+1.0_wp) / (b-1.0_wp) )
    627610                ENDIF
    628611
     
    634617!
    635618!--          Compute the gradient Richardson numbers,
    636 !--          first at the top of the constant-flux layer using u* of the
    637 !--          previous time step (+1E-30, if u* = 0), then in the remaining area.
     619!--          first at the top of the constant-flux layer using u* of the previous time step
     620!--          (+1E-30, if u* = 0), then in the remaining area.
    638621!--          There, the Ri numbers of the previous time step are used.
    639622
     
    646629                   flux = ts1d + 0.61_wp * pt_init(k) * qs1d
    647630                ENDIF
    648                 ri1d(nzb+1) = zu(nzb+1) * kappa * g * flux /                   &
    649                                ( pt_0 * ( us1d**2 + 1E-30_wp ) )
     631                ri1d(nzb+1) = zu(nzb+1) * kappa * g * flux / ( pt_0 * ( us1d**2 + 1E-30_wp ) )
    650632             ENDIF
    651633
     
    656638                ELSE
    657639                   pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
    658                    flux = ( ( pt_init(k+1) - pt_init(k-1) )                    &
    659                             + 0.61_wp                                          &
    660                             * (   pt_init(k+1) * q_init(k+1)                   &
    661                                 - pt_init(k-1) * q_init(k-1) )                 &
     640                   flux = ( ( pt_init(k+1) - pt_init(k-1) )                                        &
     641                            + 0.61_wp                                                              &
     642                            * (   pt_init(k+1) * q_init(k+1)                                       &
     643                                - pt_init(k-1) * q_init(k-1) )                                     &
    662644                          ) * dd2zu(k)
    663645                ENDIF
    664646                IF ( ri1d(k) >= 0.0_wp )  THEN
    665                    ri1d(k) = g / pt_0 * flux /                                 &
    666                               (  ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2     &
    667                                + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2     &
    668                                + 1E-30_wp                                      &
     647                   ri1d(k) = g / pt_0 * flux /                                                     &
     648                              (  ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2                         &
     649                               + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2                         &
     650                               + 1E-30_wp                                                          &
    669651                              )
    670652                ELSE
    671                    ri1d(k) = g / pt_0 * flux /                                 &
    672                               (  ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2     &
    673                                + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2     &
    674                                + 1E-30_wp                                      &
     653                   ri1d(k) = g / pt_0 * flux /                                                     &
     654                              (  ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2                         &
     655                               + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2                         &
     656                               + 1E-30_wp                                                          &
    675657                              ) * ( 1.0_wp - 16.0_wp * ri1d(k) )**0.25_wp
    676658                ENDIF
    677659             ENDDO
    678660!
    679 !--          Richardson numbers must remain restricted to a realistic value
    680 !--          range. It is exceeded excessively for very small velocities
    681 !--          (u,v --> 0).
     661!--          Richardson numbers must remain restricted to a realistic value range. It is exceeded
     662!--           excessively for very small velocities (u,v --> 0).
    682663             WHERE ( ri1d < -5.0_wp )  ri1d = -5.0_wp
    683664             WHERE ( ri1d > 1.0_wp )  ri1d = 1.0_wp
     
    691672!
    692673!--                Stable stratification
    693                    us1d = kappa * uv_total / (                                 &
    694                              LOG( zu(nzb+1) / z01d ) + 5.0_wp * ri1d(nzb+1) *  &
    695                                               ( zu(nzb+1) - z01d ) / zu(nzb+1) &
     674                   us1d = kappa * uv_total / ( LOG( zu(nzb+1) / z01d )                             &
     675                                               + 5.0_wp * ri1d(nzb+1) * ( zu(nzb+1) - z01d )       &
     676                                                 / zu(nzb+1)                                      &
    696677                                             )
    697678                ELSE
     
    699680!--                Unstable stratification
    700681                   a = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) ) )
    701                    b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) /   &
    702                                                      zu(nzb+1) * z01d ) )
    703                    us1d = kappa * uv_total / (                                 &
    704                               LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) /      &
    705                                    (1.0_wp+a) ) +                              &
    706                               2.0_wp * ( ATAN( b ) - ATAN( a ) )               &
     682                   b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / zu(nzb+1) * z01d ) )
     683                   us1d = kappa * uv_total / ( LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) /         &
     684                                                    (1.0_wp+a) ) +                                 &
     685                                               2.0_wp * ( ATAN( b ) - ATAN( a ) )                  &
    707686                                             )
    708687                ENDIF
     
    714693
    715694!
    716 !--             Boundary condition for the turbulent kinetic energy and
    717 !--             dissipation rate at the top of the constant-flux layer.
    718 !--             Additional Neumann condition de/dz = 0 at nzb is set to ensure
    719 !--             compatibility with the 3D model.
     695!--             Boundary condition for the turbulent kinetic energy and dissipation rate at the top
     696!--             of the constant-flux layer.
     697!--             Additional Neumann condition de/dz = 0 at nzb is set to ensure compatibility with
     698!--             the 3D model.
    720699                IF ( ibc_e_b == 2 )  THEN
    721700                   e1d(nzb+1) = ( us1d / c_0 )**2
     
    734713!
    735714!--                   Stable stratification
    736                       qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
    737                           ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) *  &
    738                                           ( zu(nzb+1) - z0h1d ) / zu(nzb+1)    &
    739                           )
     715                      qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /                             &
     716                             ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) *                   &
     717                                             ( zu(nzb+1) - z0h1d ) / zu(nzb+1)                     &
     718                             )
    740719                   ELSE
    741720!
    742721!--                   Unstable stratification
    743722                      a = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) )
    744                       b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) /               &
    745                                          zu(nzb+1) * z0h1d )
    746                       qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
    747                              LOG( (a-1.0_wp) / (a+1.0_wp) *                    &
    748                                   (b+1.0_wp) / (b-1.0_wp) )
     723                      b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / zu(nzb+1) * z0h1d )
     724                      qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /                             &
     725                             LOG( (a-1.0_wp) / (a+1.0_wp) * (b+1.0_wp) / (b-1.0_wp) )
    749726                   ENDIF
    750727                ELSE
     
    755732
    756733!
    757 !--          Compute the diabatic mixing length. The unstable stratification
    758 !--          must not be considered for l1d (km1d) as it is already considered
    759 !--          in the dissipation of TKE via l1d_diss. Otherwise, km1d would be
    760 !--          too large.
     734!--          Compute the diabatic mixing length. The unstable stratification must not be considered
     735!--          for l1d (km1d) as it is already considered in the dissipation of TKE via l1d_diss.
     736!--          Otherwise, km1d would be too large.
    761737             IF ( dissipation_1d /= 'prognostic' )  THEN
    762738                IF ( mixing_length_1d == 'blackadar' )  THEN
     
    767743                      ELSE
    768744                         l1d(k) = l1d_init(k)
    769                          l1d_diss(k) = l1d_init(k) *                           &
    770                                        SQRT( 1.0_wp - 16.0_wp * ri1d(k) )
     745                         l1d_diss(k) = l1d_init(k) * SQRT( 1.0_wp - 16.0_wp * ri1d(k) )
    771746                      ENDIF
    772747                   ENDDO
     
    775750                      dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
    776751                      IF ( dpt_dz > 0.0_wp )  THEN
    777                          l_stable = 0.76_wp * SQRT( e1d(k) )                   &
    778                                   / SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
     752                         l_stable = 0.76_wp * SQRT( e1d(k) )                                       &
     753                                    / SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
    779754                      ELSE
    780755                         l_stable = l1d_init(k)
     
    786761             ELSE
    787762                DO  k = nzb+1, nzt
    788                    l1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) )                   &
    789                           / ( diss1d(k) + 1.0E-30_wp )
     763                   l1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / ( diss1d(k) + 1.0E-30_wp )
    790764                ENDDO
    791765             ENDIF
    792766
    793767!
    794 !--          Compute the diffusion coefficients for momentum via the
    795 !--          corresponding Prandtl-layer relationship and according to
    796 !--          Prandtl-Kolmogorov, respectively
     768!--          Compute the diffusion coefficients for momentum via the corresponding Prandtl-layer
     769!--          relationship and according to Prandtl-Kolmogorov, respectively
    797770             IF ( constant_flux_layer )  THEN
    798771                IF ( ri1d(nzb+1) >= 0.0_wp )  THEN
    799                    km1d(nzb+1) = us1d * kappa * zu(nzb+1) /                    &
     772                   km1d(nzb+1) = us1d * kappa * zu(nzb+1) /                                        &
    800773                                 ( 1.0_wp + 5.0_wp * ri1d(nzb+1) )
    801774                ELSE
    802                    km1d(nzb+1) = us1d * kappa * zu(nzb+1) *                    &
     775                   km1d(nzb+1) = us1d * kappa * zu(nzb+1) *                                        &
    803776                                 ( 1.0_wp - 16.0_wp * ri1d(nzb+1) )**0.25_wp
    804777                ENDIF
     
    823796
    824797!
    825 !--          Compute the diffusion coefficient for heat via the relationship
    826 !--          kh = phim / phih * km
     798!--          Compute the diffusion coefficient for heat via the relationship kh = phim / phih * km
    827799             DO  k = nzb+1, nzt
    828800                IF ( ri1d(k) >= 0.0_wp )  THEN
     
    865837    ENDDO   ! time loop
    866838!
    867 !-- Set intermediate_timestep_count back to zero. This is required e.g. for
    868 !-- initial calls of calc_mean_profile.
     839!-- Set intermediate_timestep_count back to zero. This is required e.g. for initial calls of
     840!-- calc_mean_profile.
    869841    intermediate_timestep_count = 0
    870842
     
    872844
    873845
    874 !------------------------------------------------------------------------------!
     846!--------------------------------------------------------------------------------------------------!
    875847! Description:
    876848! ------------
    877849!> Compute and print out quantities for run control of the 1D model.
    878 !------------------------------------------------------------------------------!
     850!--------------------------------------------------------------------------------------------------!
    879851
    880852 SUBROUTINE run_control_1d
     
    922894       alpha = alpha / ( 2.0_wp * pi ) * 360.0_wp
    923895
    924        WRITE ( 15, 101 )  current_timestep_number_1d, simulated_time_chr, &
    925                           dt_1d, umax, vmax, us1d, alpha, energy
     896       WRITE ( 15, 101 )  current_timestep_number_1d, simulated_time_chr, dt_1d, umax, vmax, us1d, &
     897                          alpha, energy
    926898!
    927899!--    Write buffer contents to disc immediately
     
    932904!
    933905!-- formats
    934 100 FORMAT (///'1D run control output:'/ &
    935               &'------------------------------'// &
    936            &'ITER.   HH:MM:SS    DT      UMAX   VMAX    U*   ALPHA   ENERG.'/ &
    937            &'-------------------------------------------------------------')
     906100 FORMAT (///'1D run control output:'/                                                           &
     907               '------------------------------'//                                                  &
     908            'ITER.   HH:MM:SS    DT      UMAX   VMAX    U*   ALPHA   ENERG.'/                      &
     909            '-------------------------------------------------------------')
    938910101 FORMAT (I7,1X,A9,1X,F6.2,2X,F6.2,1X,F6.2,1X,F6.3,2X,F5.1,2X,F7.2)
    939911
     
    943915
    944916
    945 !------------------------------------------------------------------------------!
     917!--------------------------------------------------------------------------------------------------!
    946918! Description:
    947919! ------------
    948920!> Compute the time step w.r.t. the diffusion criterion
    949 !------------------------------------------------------------------------------!
     921!--------------------------------------------------------------------------------------------------!
    950922
    951923 SUBROUTINE timestep_1d
     
    965937
    966938!
    967 !-- Compute the currently feasible time step according to the diffusion
    968 !-- criterion. At nzb+1 the half grid length is used.
     939!-- Compute the currently feasible time step according to the diffusion criterion. At nzb+1 the half
     940!-- grid length is used.
    969941    fac = 0.125
    970942    dt_diff = dt_max_1d
     
    985957       stop_dt_1d = .TRUE.
    986958
    987        WRITE( message_string, * ) 'timestep has exceeded the lower limit&',    &
    988                                   'dt_1d = ',dt_1d,' s   simulation stopped!'
     959       WRITE( message_string, * ) 'timestep has exceeded the lower limit&', 'dt_1d = ',dt_1d,      &
     960                                  ' s   simulation stopped!'
    989961       CALL message( 'timestep_1d', 'PA0192', 1, 2, 0, 6, 0 )
    990962
     
    995967
    996968
    997 !------------------------------------------------------------------------------!
     969!--------------------------------------------------------------------------------------------------!
    998970! Description:
    999971! ------------
    1000972!> List output of profiles from the 1D-model
    1001 !------------------------------------------------------------------------------!
     973!--------------------------------------------------------------------------------------------------!
    1002974
    1003975 SUBROUTINE print_1d_model
     
    10291001       WRITE ( 17, 101 )
    10301002       DO  k = nzt+1, nzb, -1
    1031           WRITE ( 17, 103)  k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k),      &
    1032                             ri1d(k), km1d(k), kh1d(k), l1d(k), diss1d(k)
     1003          WRITE ( 17, 103)  k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), ri1d(k), km1d(k),        &
     1004                            kh1d(k), l1d(k), diss1d(k)
    10331005       ENDDO
    10341006       WRITE ( 17, 101 )
Note: See TracChangeset for help on using the changeset viewer.