Ignore:
Timestamp:
Dec 14, 2017 5:12:51 PM (4 years ago)
Author:
kanani
Message:

Merge of branch palm4u into trunk

Location:
palm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/model_1d_mod.f90

    r2339 r2696  
    11!> @file model_1d_mod.f90
    22!------------------------------------------------------------------------------!
    3 ! This file is part of PALM.
     3! This file is part of the PALM model system.
    44!
    55! PALM is free software: you can redistribute it and/or modify it under the
     
    2525! -----------------
    2626! $Id$
     27! implement TKE-e closure
     28! modification of dissipation production according to Detering and Etling
     29! reduced factor for timestep criterion to 0.125 and first dt to 1s (TG)
     30!
     31! 2339 2017-08-07 13:55:26Z gronemeier
    2732! corrected timestamp in header
    2833!
     
    131136               intermediate_timestep_count_max, kappa, km_constant,            &
    132137               message_string, mixing_length_1d, prandtl_number,               &
    133                roughness_length, run_description_header, simulated_time_chr, timestep_scheme, tsc, z0h_factor
     138               roughness_length, run_description_header, simulated_time_chr,   &
     139               timestep_scheme, tsc, z0h_factor
    134140
    135141    USE indices,                                                               &
     
    143149    IMPLICIT NONE
    144150
    145 
    146151    INTEGER(iwp) ::  current_timestep_number_1d = 0  !< current timestep number (1d-model)
    147     INTEGER(iwp) ::  damp_level_ind_1d               !< lower grid index of damping layer (1d-model) !#
     152    INTEGER(iwp) ::  damp_level_ind_1d               !< lower grid index of damping layer (1d-model)
    148153
    149154    LOGICAL ::  run_control_header_1d = .FALSE.  !< flag for output of run control header (1d-model)
    150155    LOGICAL ::  stop_dt_1d = .FALSE.             !< termination flag, used in case of too small timestep (1d-model)
    151156
     157    REAL(wp) ::  c_1 = 1.44_wp                 !< model constant
     158    REAL(wp) ::  c_2 = 1.92_wp                 !< model constant
     159    REAL(wp) ::  c_3 = 1.44_wp                 !< model constant
     160    REAL(wp) ::  c_h = 0.0015_wp               !< model constant according to Detering and Etling (1985)
    152161    REAL(wp) ::  c_m = 0.4_wp                  !< model constant, 0.4 according to Detering and Etling (1985)
    153     REAL(wp) ::  damp_level_1d = -1.0_wp       !< namelist parameter    !#
     162    REAL(wp) ::  c_mu = 0.09_wp                !< model constant
     163    REAL(wp) ::  damp_level_1d = -1.0_wp       !< namelist parameter
    154164    REAL(wp) ::  dt_1d = 60.0_wp               !< dynamic timestep (1d-model)
    155165    REAL(wp) ::  dt_max_1d = 300.0_wp          !< timestep limit (1d-model)
    156     REAL(wp) ::  dt_pr_1d = 9999999.9_wp       !< namelist parameter     !#
    157     REAL(wp) ::  dt_run_control_1d = 60.0_wp   !< namelist parameter     !#
    158     REAL(wp) ::  end_time_1d = 864000.0_wp     !< namelist parameter     !#
     166    REAL(wp) ::  dt_pr_1d = 9999999.9_wp       !< namelist parameter
     167    REAL(wp) ::  dt_run_control_1d = 60.0_wp   !< namelist parameter
     168    REAL(wp) ::  end_time_1d = 864000.0_wp     !< namelist parameter
    159169    REAL(wp) ::  qs1d                          !< characteristic humidity scale (1d-model)
    160170    REAL(wp) ::  simulated_time_1d = 0.0_wp    !< updated simulated time (1d-model)
     171    REAL(wp) ::  sig_diss = 1.3_wp             !< model constant
    161172    REAL(wp) ::  time_pr_1d = 0.0_wp           !< updated simulated time for profile output (1d-model)
    162173    REAL(wp) ::  time_run_control_1d = 0.0_wp  !< updated simulated time for run-control output (1d-model)
    163174    REAL(wp) ::  ts1d                          !< characteristic temperature scale (1d-model)
    164     REAL(wp) ::  us1d                          !< friction velocity (1d-model)                 !#
    165     REAL(wp) ::  usws1d                        !< u-component of the momentum flux (1d-model)  !#
    166     REAL(wp) ::  vsws1d                        !< v-component of the momentum flux (1d-model)  !#
     175    REAL(wp) ::  us1d                          !< friction velocity (1d-model)
     176    REAL(wp) ::  usws1d                        !< u-component of the momentum flux (1d-model)
     177    REAL(wp) ::  vsws1d                        !< v-component of the momentum flux (1d-model)
    167178    REAL(wp) ::  z01d                          !< roughness length for momentum (1d-model)
    168179    REAL(wp) ::  z0h1d                         !< roughness length for scalars (1d-model)
    169180
    170181
    171     REAL(wp), DIMENSION(:), ALLOCATABLE ::  e1d      !< tke (1d-model)                         !#
     182    REAL(wp), DIMENSION(:), ALLOCATABLE ::  diss1d   !< tke dissipation rate (1d-model)
     183    REAL(wp), DIMENSION(:), ALLOCATABLE ::  diss1d_p !< prognostic value of tke dissipation rate (1d-model)
     184    REAL(wp), DIMENSION(:), ALLOCATABLE ::  e1d      !< tke (1d-model)
    172185    REAL(wp), DIMENSION(:), ALLOCATABLE ::  e1d_p    !< prognostic value of tke (1d-model)
    173     REAL(wp), DIMENSION(:), ALLOCATABLE ::  kh1d     !< turbulent diffusion coefficient for heat (1d-model)    !#
    174     REAL(wp), DIMENSION(:), ALLOCATABLE ::  km1d     !< turbulent diffusion coefficient for momentum (1d-model)!#
     186    REAL(wp), DIMENSION(:), ALLOCATABLE ::  kh1d     !< turbulent diffusion coefficient for heat (1d-model)
     187    REAL(wp), DIMENSION(:), ALLOCATABLE ::  km1d     !< turbulent diffusion coefficient for momentum (1d-model)
    175188    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_black  !< mixing length Blackadar (1d-model)
    176     REAL(wp), DIMENSION(:), ALLOCATABLE ::  l1d      !< mixing length for turbulent diffusion coefficients (1d-model) !#
     189    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l1d      !< mixing length for turbulent diffusion coefficients (1d-model)
    177190    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l1d_diss !< mixing length for dissipation (1d-model)
    178     REAL(wp), DIMENSION(:), ALLOCATABLE ::  rif1d    !< Richardson flux number (1d-model)      !#
     191    REAL(wp), DIMENSION(:), ALLOCATABLE ::  rif1d    !< Richardson flux number (1d-model)
     192    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_diss  !< tendency of diss (1d-model)
     193    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_dissm !< weighted tendency of diss for previous sub-timestep (1d-model)
    179194    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_e     !< tendency of e (1d-model)
    180195    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_em    !< weighted tendency of e for previous sub-timestep (1d-model)
     
    183198    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_v     !< tendency of v (1d-model)
    184199    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_vm    !< weighted tendency of v for previous sub-timestep (1d-model)
    185     REAL(wp), DIMENSION(:), ALLOCATABLE ::  u1d      !< u-velocity component (1d-model)        !#
     200    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u1d      !< u-velocity component (1d-model)
    186201    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u1d_p    !< prognostic value of u-velocity component (1d-model)
    187     REAL(wp), DIMENSION(:), ALLOCATABLE ::  v1d      !< v-velocity component (1d-model)        !#
     202    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v1d      !< v-velocity component (1d-model)
    188203    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v1d_p    !< prognostic value of v-velocity component (1d-model)
    189204
     
    227242!
    228243!-- Public variables
    229     PUBLIC  damp_level_1d, damp_level_ind_1d, dt_pr_1d, dt_run_control_1d,     &
    230             e1d, end_time_1d, kh1d, km1d, l1d, rif1d, u1d, us1d, usws1d, v1d,  &
    231             vsws1d
     244    PUBLIC  damp_level_1d, damp_level_ind_1d, diss1d, dt_pr_1d,                &
     245            dt_run_control_1d, e1d, end_time_1d, kh1d, km1d, l1d, rif1d, u1d,  &
     246            us1d, usws1d, v1d, vsws1d
    232247
    233248
     
    238253    IMPLICIT NONE
    239254
    240     CHARACTER (LEN=9) ::  time_to_string  !<
     255    CHARACTER (LEN=9) ::  time_to_string  !< function to transform time from real to character string
     256
     257    INTEGER(iwp) ::  k  !< loop index
    241258   
    242     INTEGER(iwp) ::  k  !<
    243    
    244     REAL(wp) ::  lambda !<
     259    REAL(wp) ::  lambda !< maximum mixing length
    245260
    246261!
    247262!-- Allocate required 1D-arrays
    248     ALLOCATE( e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1),               &
     263    ALLOCATE( diss1d(nzb:nzt+1), diss1d_p(nzb:nzt+1),                          &
     264              e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1),               &
    249265              km1d(nzb:nzt+1), l_black(nzb:nzt+1), l1d(nzb:nzt+1),             &
    250               l1d_diss(nzb:nzt+1), rif1d(nzb:nzt+1), te_e(nzb:nzt+1),          &
     266              l1d_diss(nzb:nzt+1), rif1d(nzb:nzt+1), te_diss(nzb:nzt+1),       &
     267              te_dissm(nzb:nzt+1), te_e(nzb:nzt+1),                            &
    251268              te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1),             &
    252269              te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1),               &
     
    259276       kh1d = km_constant / prandtl_number
    260277    ELSE
     278       diss1d = 0.0_wp; diss1d_p = 0.0_wp
    261279       e1d = 0.0_wp; e1d_p = 0.0_wp
    262280       kh1d = 0.0_wp; km1d = 0.0_wp
     
    309327       us1d = 0.1_wp   ! without initial friction the flow would not change
    310328    ELSE
     329       diss1d(nzb+1) = 1.0_wp
    311330       e1d(nzb+1)  = 1.0_wp
    312331       km1d(nzb+1) = 1.0_wp
     
    323342!-- Tendencies must be preset in order to avoid runtime errors within the
    324343!-- first Runge-Kutta step
     344    te_dissm = 0.0_wp
    325345    te_em = 0.0_wp
    326346    te_um = 0.0_wp
     
    354374    IMPLICIT NONE
    355375
    356     CHARACTER (LEN=9) ::  time_to_string  !<
    357    
     376    CHARACTER (LEN=9) ::  time_to_string  !< function to transform time from real to character string
     377
    358378    INTEGER(iwp) ::  k  !< loop index
    359    
     379
    360380    REAL(wp) ::  a            !< auxiliary variable
    361381    REAL(wp) ::  b            !< auxiliary variable
    362     REAL(wp) ::  dissipation  !< dissipation of TKE
    363382    REAL(wp) ::  dpt_dz       !< vertical temperature gradient
    364383    REAL(wp) ::  flux         !< vertical temperature gradient
     
    372391!-- Determine the time step at the start of a 1D-simulation and
    373392!-- determine and printout quantities used for run control
    374     dt_1d = 10.0_wp
     393    dt_1d = 1.0_wp
    375394    CALL run_control_1d
    376395
     
    392411
    393412!
    394 !--       Compute all tendency terms. If a Prandtl-layer is simulated, k starts
    395 !--       at nzb+2.
     413!--       Compute all tendency terms. If a constant-flux layer is simulated,
     414!--       k starts at nzb+2.
    396415          DO  k = nzb_diff, nzt
    397416
     
    414433             DO  k = nzb_diff, nzt
    415434!
    416 !--             TKE
     435!--             TKE and dissipation rate
    417436                kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) )
    418437                kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) )
     
    428447                ENDIF
    429448
     449!
     450!--             Calculate dissipation rate if no prognostic equation is used for
     451!--             dissipation rate
    430452                IF ( dissipation_1d == 'detering' )  THEN
    431 !
    432 !--                According to Detering, c_e=c_m^3
    433                    dissipation = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
     453                   diss1d(k) = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    434454                ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
    435                    dissipation = ( 0.19_wp                                     &
    436                                    + 0.74_wp * l1d_diss(k) /  l_grid(k)        &
    437                                  ) * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    438                 ENDIF
    439 
     455                   diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) /  l_grid(k)  &
     456                               ) * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
     457                ENDIF
     458!
     459!--             TKE
    440460                te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
    441461                                    + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
     
    446466                                   - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
    447467                                                 ) * ddzw(k)                   &
    448                                    - dissipation
     468                                   - diss1d(k)
     469
     470                IF ( dissipation_1d == 'prognostic' )  THEN
     471!
     472!--                dissipation rate
     473                   te_diss(k) = km1d(k) *                                      &
     474                                  ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2  &
     475                                  + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2  &
     476                                  ) * c_1 * c_mu**0.75 / c_h * f / us1d        &
     477                                    * SQRT(e1d(k))                             &
     478                                  - g / pt_0 * kh1d(k) * flux * c_3            &
     479                                    * diss1d(k) / ( e1d(k) + 1.0E-20_wp )      &
     480                                  + ( kmzp * ( diss1d(k+1) - diss1d(k) )       &
     481                                           * ddzu(k+1)                         &
     482                                    - kmzm * ( diss1d(k) - diss1d(k-1) )       &
     483                                           * ddzu(k)                           &
     484                                    ) * ddzw(k) / sig_diss                     &
     485                                  - c_2 * diss1d(k)**2 / ( e1d(k) + 1.0E-20_wp )
     486
     487                ENDIF
     488
    449489             ENDDO
    450490          ENDIF
    451491
    452492!
    453 !--       Tendency terms at the top of the Prandtl-layer.
     493!--       Tendency terms at the top of the constant-flux layer.
    454494!--       Finite differences of the momentum fluxes are computed using half the
    455495!--       normal grid length (2.0*ddzw(k)) for the sake of enhanced accuracy
     
    470510             ENDIF
    471511
     512!
     513!--          Calculate dissipation rate if no prognostic equation is used for
     514!--          dissipation rate
    472515             IF ( dissipation_1d == 'detering' )  THEN
    473 !
    474 !--             According to Detering, c_e=0.064
    475                 dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
     516                diss1d(k) = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    476517             ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
    477                 dissipation = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l_grid(k) )  &
    478                               * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
     518                diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l_grid(k) )    &
     519                            * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    479520             ENDIF
    480521
     
    491532!
    492533!--          TKE
    493              te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2   &
    494                                  + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2   &
    495                                  )                                             &
    496                                  - g / pt_0 * kh1d(k) * flux                   &
    497                                  +           (                                 &
    498                                   kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)     &
    499                                 - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)       &
    500                                               ) * ddzw(k)                      &
    501                                 - dissipation
     534             IF ( .NOT. dissipation_1d == 'prognostic' )  THEN
     535                te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
     536                                    + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
     537                                    )                                          &
     538                                    - g / pt_0 * kh1d(k) * flux                &
     539                                    +           (                              &
     540                                     kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)  &
     541                                   - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
     542                                                 ) * ddzw(k)                   &
     543                                   - diss1d(k)
     544             ENDIF
     545
    502546          ENDIF
    503547
     
    513557          ENDDO
    514558          IF ( .NOT. constant_diffusion )  THEN
     559
    515560             DO  k = nzb+1, nzt
    516 
    517561                e1d_p(k) = e1d(k) + dt_1d * ( tsc(2) * te_e(k) + &
    518562                                              tsc(3) * te_em(k) )
    519 
    520563             ENDDO
     564
     565             IF ( dissipation_1d == 'prognostic' )  THEN
     566                DO  k = nzb_diff, nzt
     567                   diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + &
     568                                                 tsc(3) * te_dissm(k) )
     569                ENDDO
     570             ENDIF
    521571!
    522572!--          Eliminate negative TKE values, which can result from the
     
    540590                      te_em(k) = te_e(k)
    541591                   ENDDO
     592                   IF ( dissipation_1d == 'prognostic' )  THEN
     593                      DO k = nzb+1, nzt
     594                         te_dissm(k) = te_diss(k)
     595                      ENDDO
     596                   ENDIF
    542597                ENDIF
    543598
     
    554609                      te_em(k) = -9.5625_wp * te_e(k) + 5.3125_wp * te_em(k)
    555610                   ENDDO
     611                   IF ( dissipation_1d == 'prognostic' )  THEN
     612                      DO k = nzb+1, nzt
     613                         te_dissm(k) = -9.5625_wp * te_diss(k) + 5.3125_wp * te_dissm(k)
     614                      ENDDO
     615                   ENDIF
    556616                ENDIF
    557617
     
    562622!
    563623!--       Boundary conditions for the prognostic variables.
    564 !--       At the top boundary (nzt+1) u,v and e keep their initial values
    565 !--       (ug(nzt+1), vg(nzt+1), 0).
     624!--       At the top boundary (nzt+1) u, v, e, and diss keep their initial
     625!--       values (ug(nzt+1), vg(nzt+1), 0, 0).
    566626!--       At the bottom boundary, Dirichlet condition is used for u and v (0)
    567 !--       and Neumann condition for e (e(nzb)=e(nzb+1)).
     627!--       and Neumann condition for e and diss (e(nzb)=e(nzb+1)).
    568628          u1d_p(nzb) = 0.0_wp
    569629          v1d_p(nzb) = 0.0_wp
     
    575635          IF ( .NOT. constant_diffusion )  THEN
    576636             e1d  = e1d_p
     637             IF ( dissipation_1d == 'prognostic' )  THEN
     638                diss1d = diss1d_p
     639             ENDIF
    577640          ENDIF
    578641
     
    582645
    583646!
    584 !--          First compute the vertical fluxes in the Prandtl-layer
     647!--          First compute the vertical fluxes in the constant-flux layer
    585648             IF ( constant_flux_layer )  THEN
    586649!
     
    609672!
    610673!--          Compute the Richardson-flux numbers,
    611 !--          first at the top of the Prandtl-layer using u* of the previous
    612 !--          time step (+1E-30, if u* = 0), then in the remaining area. There
    613 !--          the rif-numbers of the previous time step are used.
     674!--          first at the top of the constant-flux layer using u* of the
     675!--          previous time step (+1E-30, if u* = 0), then in the remaining area.
     676!--          There the rif-numbers of the previous time step are used.
    614677
    615678             IF ( constant_flux_layer )  THEN
     
    689752
    690753!
    691 !--             Boundary condition for the turbulent kinetic energy at the top
    692 !--             of the Prandtl-layer. c_m = 0.4 according to Detering.
     754!--             Boundary condition for the turbulent kinetic energy and
     755!--             dissipation rate at the top of the constant-flux layer.
    693756!--             Additional Neumann condition de/dz = 0 at nzb is set to ensure
    694757!--             compatibility with the 3D model.
    695758                IF ( ibc_e_b == 2 )  THEN
    696759                   e1d(nzb+1) = ( us1d / c_m )**2
     760                ENDIF
     761                IF ( dissipation_1d == 'prognostic' )  THEN
     762                   e1d(nzb+1) = us1d**2 / SQRT( c_mu )
     763                   diss1d(nzb+1) = us1d**3 / ( kappa * zu(nzb+1) )
     764                   diss1d(nzb) = diss1d(nzb+1)
    697765                ENDIF
    698766                e1d(nzb) = e1d(nzb+1)
     
    767835                ENDIF
    768836             ENDIF
    769              DO  k = nzb_diff, nzt
    770                 km1d(k) = c_m * SQRT( e1d(k) ) * l1d(k)
    771              ENDDO
     837             IF ( dissipation_1d == 'prognostic' )  THEN
     838                DO  k = nzb_diff, nzt
     839                   km1d(k) = c_mu * e1d(k)**2 / ( diss1d(k) + 1.0E-30_wp )
     840                ENDDO
     841             ELSE
     842                DO  k = nzb_diff, nzt
     843                   km1d(k) = c_m * SQRT( e1d(k) ) * l1d(k)
     844                ENDDO
     845             ENDIF
    772846
    773847!
     
    918992!-- Compute the currently feasible time step according to the diffusion
    919993!-- criterion. At nzb+1 the half grid length is used.
    920     fac = 0.35_wp
     994    fac = 0.125  !0.35_wp                                                       !### changed from 0.35
    921995    dt_diff = dt_max_1d
    922996    DO  k = nzb+2, nzt
     
    9771051       DO  k = nzt+1, nzb, -1
    9781052          WRITE ( 17, 103)  k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), &
    979                             rif1d(k), km1d(k), kh1d(k), l1d(k), zu(k), k
     1053                            rif1d(k), km1d(k), kh1d(k), l1d(k), diss1d(k)
    9801054       ENDDO
    9811055       WRITE ( 17, 101 )
     
    9931067100 FORMAT ('# ',A/'#',10('-')/'# 1d-model profiles')
    9941068104 FORMAT (//'# Time: ',A)
    995 101 FORMAT ('#',79('-'))
    996 102 FORMAT ('#  k     zu      u      v     pt      e    rif    Km    Kh     ', &
    997             'l      zu      k')
    998 103 FORMAT (1X,I4,1X,F7.1,1X,F6.2,1X,F6.2,1X,F6.2,1X,F6.2,1X,F5.2,1X,F5.2, &
    999             1X,F5.2,1X,F6.2,1X,F7.1,2X,I4)
     1069101 FORMAT ('#',111('-'))
     1070102 FORMAT ('#  k     zu      u          v          pt         e          ',   &
     1071            'rif        Km         Kh         l          diss   ')
     1072103 FORMAT (1X,I4,1X,F7.1,9(1X,E10.3))
    10001073
    10011074
Note: See TracChangeset for help on using the changeset viewer.