Ignore:
Timestamp:
Aug 7, 2017 8:59:53 AM (7 years ago)
Author:
gronemeier
Message:

changes to the 1D model

File:
1 edited

Legend:

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

    r2334 r2337  
    2525! -----------------
    2626! $Id$
     27! revised calculation of mixing length
     28! removed rounding of time step
     29! corrected calculation of virtual potential temperature
     30!
     31! 2334 2017-08-04 11:57:04Z gronemeier
    2732! set c_m = 0.4 according to Detering and Etling (1985)
    2833!
     
    120125   
    121126    USE model_1d,                                                              &
    122         ONLY:  e1d, e1d_p, kh1d, km1d, l1d, l_black, qs1d, rif1d,              &
     127        ONLY:  e1d, e1d_p, kh1d, km1d, l1d, l1d_diss, l_black, qs1d, rif1d,    &
    123128               simulated_time_1d, te_e, te_em, te_u, te_um, te_v, te_vm, ts1d, &
    124129               u1d, u1d_p, us1d, usws1d, v1d, v1d_p, vsws1d, z01d, z0h1d
     
    141146    ALLOCATE( e1d(nzb:nzt+1),    e1d_p(nzb:nzt+1),                             &
    142147              kh1d(nzb:nzt+1),   km1d(nzb:nzt+1),                              &
    143               l_black(nzb:nzt+1), l1d(nzb:nzt+1),                              &
     148              l_black(nzb:nzt+1), l1d(nzb:nzt+1), l1d_diss(nzb:nzt+1),         &
    144149              rif1d(nzb:nzt+1),   te_e(nzb:nzt+1),                             &
    145150              te_em(nzb:nzt+1),  te_u(nzb:nzt+1),    te_um(nzb:nzt+1),         &
     
    183188       ENDIF
    184189    ENDIF
    185     l1d   = l_black
    186     u1d   = u_init
    187     u1d_p = u_init
    188     v1d   = v_init
    189     v1d_p = v_init
     190    l1d      = l_black
     191    l1d_diss = l_black
     192    u1d      = u_init
     193    u1d_p    = u_init
     194    v1d      = v_init
     195    v1d_p    = v_init
    190196
    191197!
     
    226232
    227233!
    228 !-- Integrate the 1D-model equations using the leap-frog scheme
     234!-- Integrate the 1D-model equations using the Runge-Kutta scheme
    229235    CALL time_integration_1d
    230236
     
    261267        ONLY:  current_timestep_number_1d, damp_level_ind_1d, dt_1d,           &
    262268               dt_pr_1d, dt_run_control_1d, e1d, e1d_p, end_time_1d,           &
    263                kh1d, km1d, l1d, l_black, qs1d, rif1d, simulated_time_1d,      &
     269               kh1d, km1d, l1d, l1d_diss, l_black, qs1d, rif1d, simulated_time_1d, &
    264270               stop_dt_1d, te_e, te_em, te_u, te_um, te_v, te_vm, time_pr_1d,  &
    265271               ts1d, time_run_control_1d, u1d, u1d_p, us1d, usws1d, v1d,       &
     
    340346                   pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
    341347                   flux = ( ( pt_init(k+1) - pt_init(k-1) ) +                  &
    342                             0.61_wp * pt_init(k) *                             &
    343                             ( q_init(k+1) - q_init(k-1) ) ) * dd2zu(k)
     348                            0.61_wp * ( pt_init(k+1) * q_init(k+1) -           &
     349                                        pt_init(k-1) * q_init(k-1)   )         &
     350                          ) * dd2zu(k)
    344351                ENDIF
    345352
    346353                IF ( dissipation_1d == 'detering' )  THEN
    347354!
    348 !--                According to Detering, c_e=0.064
    349                    dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d(k)
     355!--                According to Detering, c_e=c_m^3
     356                   dissipation = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    350357                ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
    351                    dissipation = ( 0.19_wp + 0.74_wp * l1d(k) / l_grid(k) )    &
    352                                  * e1d(k) * SQRT( e1d(k) ) / l1d(k)
     358                   dissipation = ( 0.19_wp                                     &
     359                                   + 0.74_wp * l1d_diss(k) /  l_grid(k)        &
     360                                 ) * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    353361                ENDIF
    354362
     
    380388                pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
    381389                flux = ( ( pt_init(k+1) - pt_init(k-1) ) +                     &
    382                          0.61_wp * pt_init(k) * ( q_init(k+1) - q_init(k-1) )  &
     390                         0.61_wp * ( pt_init(k+1) * q_init(k+1) -              &
     391                                     pt_init(k-1) * q_init(k-1)   )            &
    383392                       ) * dd2zu(k)
    384393             ENDIF
     
    387396!
    388397!--             According to Detering, c_e=0.064
    389                 dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d(k)
     398                dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    390399             ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
    391                 dissipation = ( 0.19_wp + 0.74_wp * l1d(k) / l_grid(k) )       &
    392                               * e1d(k) * SQRT( e1d(k) ) / l1d(k)
     400                dissipation = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l_grid(k) )  &
     401                              * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    393402             ENDIF
    394403
     
    513522                   b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /                 &
    514523                       zu(nzb+1) * z0h1d )
    515 !
    516 !--                In the borderline case the formula for stable stratification
    517 !--                must be applied, because otherwise a zero division would
    518 !--                occur in the argument of the logarithm.
    519                    IF ( a == 0.0_wp  .OR.  b == 0.0_wp )  THEN
    520                       ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /       &
    521                              ( LOG( zu(nzb+1) / z0h1d ) +                      &
    522                                5.0_wp * rif1d(nzb+1) *                         &
    523                                ( zu(nzb+1) - z0h1d ) / zu(nzb+1)               &
    524                              )
    525                    ELSE
    526                       ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /       &
    527                              LOG( (a-1.0_wp) / (a+1.0_wp) *                    &
    528                                   (b+1.0_wp) / (b-1.0_wp) )
    529                    ENDIF
     524
     525                   ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /          &
     526                          LOG( (a-1.0_wp) / (a+1.0_wp) *                       &
     527                               (b+1.0_wp) / (b-1.0_wp) )
    530528                ENDIF
    531529
     
    557555                   pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
    558556                   flux = ( ( pt_init(k+1) - pt_init(k-1) )                    &
    559                             + 0.61_wp * pt_init(k)                             &
    560                             * ( q_init(k+1) - q_init(k-1) )                    &
     557                            + 0.61_wp                                          &
     558                            * (   pt_init(k+1) * q_init(k+1)                   &
     559                                - pt_init(k-1) * q_init(k-1) )                 &
    561560                          ) * dd2zu(k)
    562561                ENDIF
     
    600599                   b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /  &
    601600                                                     zu(nzb+1) * z01d ) )
    602 !
    603 !--                In the borderline case the formula for stable stratification
    604 !--                must be applied, because otherwise a zero division would
    605 !--                occur in the argument of the logarithm.
    606                    IF ( a == 1.0_wp  .OR.  b == 1.0_wp )  THEN
    607                       us1d = kappa * uv_total / (                              &
    608                              LOG( zu(nzb+1) / z01d ) +                         &
    609                              5.0_wp * rif1d(nzb+1) * ( zu(nzb+1) - z01d ) /    &
    610                                                   zu(nzb+1) )
    611                    ELSE
    612                       us1d = kappa * uv_total / (                              &
    613                                  LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) /   &
    614                                       (1.0_wp+a) ) +                           &
    615                                  2.0_wp * ( ATAN( b ) - ATAN( a ) )            &
    616                                                 )
    617                    ENDIF
     601                   us1d = kappa * uv_total / (                                 &
     602                              LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) /      &
     603                                   (1.0_wp+a) ) +                              &
     604                              2.0_wp * ( ATAN( b ) - ATAN( a ) )               &
     605                                             )
    618606                ENDIF
    619607
     
    649637                      b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /              &
    650638                                         zu(nzb+1) * z0h1d )
    651 !
    652 !--                   In the borderline case the formula for stable stratification
    653 !--                   must be applied, because otherwise a zero division would
    654 !--                   occur in the argument of the logarithm.
    655                       IF ( a == 1.0_wp  .OR.  b == 1.0_wp )  THEN
    656                          qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /      &
    657                                 ( LOG( zu(nzb+1) / z0h1d ) +                   &
    658                                   5.0_wp * rif1d(nzb+1) *                      &
    659                                   ( zu(nzb+1) - z0h1d ) / zu(nzb+1)            &
    660                                 )
    661                       ELSE
    662                          qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /      &
    663                                 LOG( (a-1.0_wp) / (a+1.0_wp) *                 &
    664                                      (b+1.0_wp) / (b-1.0_wp) )
    665                       ENDIF
    666                    ENDIF               
     639                      qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
     640                             LOG( (a-1.0_wp) / (a+1.0_wp) *                    &
     641                                  (b+1.0_wp) / (b-1.0_wp) )
     642                   ENDIF
    667643                ELSE
    668644                   qs1d = 0.0_wp
    669                 ENDIF             
     645                ENDIF
    670646
    671647             ENDIF   !  constant_flux_layer
    672648
    673649!
    674 !--          Compute the diabatic mixing length
     650!--          Compute the diabatic mixing length. The unstable stratification
     651!--          must not be considered for l1d (km1d) as it is already considered
     652!--          in the dissipation of TKE via l1d_diss. Otherwise, km1d would be
     653!--          too large.
    675654             IF ( mixing_length_1d == 'blackadar' )  THEN
    676655                DO  k = nzb+1, nzt
    677656                   IF ( rif1d(k) >= 0.0_wp )  THEN
    678657                      l1d(k) = l_black(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
     658                      l1d_diss(k) = l1d(k)
    679659                   ELSE
    680                       l1d(k) = l_black(k) *                                    &
    681                                ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp
     660                      l1d(k) = l_black(k)
     661                      l1d_diss(k) = l_black(k) *                               &
     662                                    SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
    682663                   ENDIF
    683                    l1d(k) = l_black(k)
    684664                ENDDO
    685 
    686665             ELSEIF ( mixing_length_1d == 'as_in_3d_model' )  THEN
    687666                DO  k = nzb+1, nzt
     
    694673                   ENDIF
    695674                   l1d(k) = MIN( l_grid(k), l_stable )
     675                   l1d_diss(k) = l1d(k)
    696676                ENDDO
    697677             ENDIF
     
    700680!--          Compute the diffusion coefficients for momentum via the
    701681!--          corresponding Prandtl-layer relationship and according to
    702 !--          Prandtl-Kolmogorov, respectively. The unstable stratification is
    703 !--          computed via the adiabatic mixing length, for the unstability has
    704 !--          already been taken account of via the TKE (cf. also Diss.).
     682!--          Prandtl-Kolmogorov, respectively
    705683             IF ( constant_flux_layer )  THEN
    706684                IF ( rif1d(nzb+1) >= 0.0_wp )  THEN
     
    713691             ENDIF
    714692             DO  k = nzb_diff, nzt
    715                 km1d(k) = c_m * SQRT( e1d(k) )
    716 
    717                 IF ( rif1d(k) >= 0.0_wp )  THEN
    718                    km1d(k) = km1d(k) * l1d(k)
    719                 ELSE
    720                    km1d(k) = km1d(k) * l_black(k)
    721                 ENDIF
     693                km1d(k) = c_m * SQRT( e1d(k) ) * l1d(k)
    722694             ENDDO
    723695
     
    882854   
    883855    USE model_1d,                                                              &
    884         ONLY:  dt_1d, dt_max_1d, km1d, old_dt_1d, stop_dt_1d
     856        ONLY:  dt_1d, dt_max_1d, km1d, stop_dt_1d
    885857   
    886858    USE pegrid
     
    922894    ENDIF
    923895
    924 !
    925 !-- A more or less simple new time step value is obtained taking only the
    926 !-- first two significant digits
    927     div = 1000.0_wp
    928     DO  WHILE ( dt_1d < div )
    929        div = div / 10.0_wp
    930     ENDDO
    931     dt_1d = NINT( dt_1d * 100.0_wp / div ) * div / 100.0_wp
    932 
    933     old_dt_1d = dt_1d
    934 
    935 
    936896 END SUBROUTINE timestep_1d
    937897
Note: See TracChangeset for help on using the changeset viewer.